msd.F90

Go to the documentation of this file.

! This file is part of the 'atomes' software.
!
! 'atomes' is free software: you can redistribute it and/or modify it under the terms
! of the GNU Affero General Public License as published by the Free Software Foundation,
! either version 3 of the License, or (at your option) any later version.
!
! 'atomes' is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
! without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
! See the GNU General Public License for more details.
!
! You should have received a copy of the GNU Affero General Public License along with 'atomes'.
! If not, see <https://www.gnu.org/licenses/>
!
! Copyright (C) 2022-2024 by CNRS and University of Strasbourg
!
 
INTEGER (KIND=c_int) FUNCTION msd (DLT, NDTS) bind (C,NAME='msd_')
 
!
! Mean Square Displacement
!
 
USE parameters
 
#ifdef OPENMP
!$ USE OMP_LIB
#endif
IMPLICIT NONE
 
INTEGER (KIND=c_int), INTENT(IN) :: ndts
real(kind=c_double), INTENT(IN) :: dlt
DOUBLE PRECISION :: masstot
DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: msdtab
#ifdef OPENMP
INTEGER :: numth
#endif
 
INTERFACE
  LOGICAL FUNCTION allocmsd()
  END FUNCTION
END INTERFACE
 
! Calcul du déplacement carré moyen
if(.not.transpo()) then
  msd=0
  goto 001
endif
 
if (.not. allocmsd()) then
  msd=0
  goto 001
endif
 
masstot=0.0d0
do j=1, nsp
  masstot=masstot+nbspbs(j)*mass(j)
enddo
 
#ifdef OPENMP
numth = omp_get_max_threads()
if (ns.lt.numth) numth=ns
 
! OpemMP on MD steps
!$OMP PARALLEL NUM_THREADS(NUMTH) DEFAULT (NONE) &
!$OMP& PRIVATE(RCm, RCm2, R2Cor, Dij, Rij, i, j, k, l, m, n, o, p, Vij) &
!$OMP& SHARED(NUMTH, NS, NA, D2i, D2dir, DRIFT, LOT, MASS, MASSTOT, NDTS, DLT, NFULLPOS)
!$OMP DO SCHEDULE(STATIC,NS-1/NUMTH)
#endif
do j=1, ns-1
 
  do m=1, 3
    rcm(m)=0.0d0
  enddo
 
  do i=1, na
 
    k=lot(i)
    do m=1, 3
      rcm(m)=rcm(m)+mass(k)*nfullpos(i,m,j)
    enddo
  enddo
 
  do m=1, 3
    rcm(m)=rcm(m)/masstot
  enddo
 
  do k=j+1, ns
 
    do m=1, 3
      rcm2(m)=0.0d0
    enddo
 
    do i=1, na
      n=lot(i)
      do m=1, 3
        rcm2(m)=rcm2(m)+mass(n)*nfullpos(i,m,k)
      enddo
    enddo
 
    do m=1, 3
      rcm2(m)=rcm2(m)/masstot
    enddo
 
    do i=1, na
 
      o=lot(i)
      dij=0.0d0
      do m=1,3
        rij(m)=nfullpos(i,m,k)-nfullpos(i,m,j)
        r2cor(m) = rcm2(m) - rcm(m)
        vij = (rij(m)-r2cor(m))**2
#ifdef OPENMP
        !$OMP ATOMIC
#endif
        d2dir(o,m,k-j)=d2dir(o,m,k-j)+vij
        dij=dij+vij
      enddo
#ifdef OPENMP
      !$OMP ATOMIC
#endif
      d2i(o,k-j)=d2i(o,k-j)+dij
 
#ifdef OPENMP
      !$OMP CRITICAL
#endif
      p=4
      do m=1, 2
      do n=m+1, 3
        d2dir(o,p,k-j)=d2dir(o,m,k-j)+d2dir(o,n,k-j)
        p=p+1
      enddo
      enddo
#ifdef OPENMP
      !$OMP END CRITICAL
#endif
    enddo
 
    if (k .eq. j+1) then
 
      do i=1, na
        o=lot(i)
        do m=1,3
          drift(m,k)=drift(m,k)+1e5*(nfullpos(i,m,k)-nfullpos(i,m,j))*mass(o)/(ndts*dlt)
        enddo
      enddo
 
      do m=1,3
        drift(m,k)=drift(m,k)/masstot
      enddo
 
    endif
 
  enddo
enddo
#ifdef OPENMP
!$OMP END DO NOWAIT
!$OMP END PARALLEL
#endif
 
do k=1, ns-1
 
  l=k+1
  do m=1,3
    rcm(m)=0.0d0
    rcm2(m)=0.0d0
  enddo
 
  do i=1, na
    o=lot(i)
    do m=1, 3
        rcm(m)=rcm(m)+ mass(o)*nfullpos(i,m,k)
        rcm2(m)=rcm2(m)+ mass(o)*nfullpos(i,m,l)
    enddo
  enddo
 
  do m=1, 3
    rcm(m)=rcm(m)/masstot
    rcm2(m)=rcm2(m)/masstot
  enddo
 
  do i=1, na
 
    o=lot(i)
    dij=0.0d0
    do m=1,3
      r2cor(m) = rcm2(m) - rcm(m)
      rij(m)=nfullpos(i,m,l)-nfullpos(i,m,k)
      d2dirnac(o,m,k)=(rij(m)-r2cor(m))**2
      dij=dij+(rij(m)-r2cor(m))**2
      cor(m,k)=cor(m,k)+r2cor(m)
    enddo
 
    p=4
    do m=1, 2
    do n=m+1, 3
      d2dirnac(o,p,k)=d2dirnac(o,m,k)+d2dirnac(o,n,k)
      p=p+1
    enddo
    enddo
    d2inac(o,k)=d2inac(o,k)+dij
 
  enddo
enddo
 
do k=1, ns-1
  do i=1, 6
    do l=1, nsp
      d2dir(l,i,k)=d2dir(l,i,k)/(ns-k)/nbspbs(l)
      d2dirnac(l,i,k)=d2dirnac(l,i,k)/nbspbs(l)
    enddo
  enddo
enddo
 
do k=1, ns-1
  do i=1, nsp
    d2i(i,k)=d2i(i,k)/(ns-k)/nbspbs(i)
    d2inac(i,k)=d2inac(i,k)/nbspbs(i)
  enddo
enddo
 
! Evaluation of the diffusion constant D
! A minimum of dynamic time is needed for a 'physical' meaning of D
! so if the simulation time is big enough we compute D
! Here arbitrary we choose 30ps as time limit (30000 fs)
! We also decide to evaluate only if there is more than 1000 MD steps
if (ns*ndts*dlt .gt. 30000 .and. ns.gt.1000) then
! then take car to cut the first 1 ps of the simulation
! not to be anymore in the ballistic regime
! afterwards take care to cut the last 1 ps of the calculation
! usually the correlations are unclear in this zone of the MSD
  z=ns*ndts*dlt
  z=z-1.0
  k=int(z/(ndts*dlt))
  l=int(1.0/(ndts*dlt))
  do i=1, nsp
    dcte(i)=(d2i(i,k)-d2i(i,l))/((k-l)*ndts*dlt)
  enddo
 
else
 
  do i=1, nsp
    dcte(i)=0.d0
  enddo
 
endif
 
allocate(msdtab(ns-1), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: MSD"//char(0), "Table: MSDTAB"//char(0))
  msd=0
  goto 001
endif
 
k=0
do i=1, nsp
  do j=1, ns-1
    msdtab(j)=d2i(i,j)
  enddo
  call save_curve (ns-1, msdtab, k, idmsd)
  k=k+1
  do j=1, ns-1
    msdtab(j)=d2inac(i,j)
  enddo
  call save_curve (ns-1, msdtab, k, idmsd)
  k=k+1
enddo
 
do i=1, nsp
  do l=1, 6
    do j=1, ns-1
      msdtab(j)=d2dir(i,l,j)
    enddo
    call save_curve (ns-1, msdtab, k, idmsd)
    k=k+1
  enddo
enddo
 
do i=1, nsp
  do l=1, 6
    do j=1, ns-1
      msdtab(j)=d2dirnac(i,l,j)
    enddo
    call save_curve (ns-1, msdtab, k, idmsd)
    k=k+1
  enddo
enddo
 
do i=1, 3
  do j=1, ns-1
    msdtab(j)=cor(i,j)
  enddo
  call save_curve (ns-1, msdtab, k, idmsd)
  k=k+1
enddo
do i=1, 3
  do j=1, ns-1
    msdtab(j)=drift(i,j)
  enddo
  call save_curve (ns-1, msdtab, k, idmsd)
  k=k+1
enddo
 
msd=1
 
001 continue
 
if (allocated(nfullpos)) deallocate(nfullpos)
 
call deallocmsd
 
CONTAINS
 
LOGICAL FUNCTION transpo()
 
USE parameters
 
if (allocated(nfullpos)) deallocate(nfullpos)
allocate(nfullpos(na,3,ns), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: TRANSPO"//char(0), "Table: NFULLPOS"//char(0))
  transpo=.false.
  goto 001
endif
 
if (allocated(poa)) deallocate(poa)
allocate(poa(na,3), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: TRANSPO"//char(0), "Table: POA"//char(0))
  transpo=.false.
  goto 001
endif
 
if (allocated(pob)) deallocate(pob)
allocate(pob(na,3), stat=err)
if (err .ne. 0) then
  call show_error ("Impossible to allocate memory"//char(0), &
                   "Function: TRANSPO"//char(0), "Table: POB"//char(0))
  transpo=.false.
  goto 001
endif
 
do noa=1, ns
do nob=1, na
do noc=1, 3
  nfullpos(nob,noc,noa) = fullpos(nob,noc,noa)
enddo
enddo
enddo
 
do noa=1, na
  do nob=1, 3
    poa(noa,nob)=nfullpos(noa,nob,1)
  enddo
enddo
 
do noc=2, ns
  do noa=1, na
    do nob=1, 3
      pob(noa,nob)=nfullpos(noa,nob,noc)
    enddo
    if (ncells .gt. 1) then
      call calcrij(noa, noa, -1, -1, noc-1)
    else
      call calcrij(noa, noa, -1, -1, 1)
    endif
    do nob=1, 3
      poa(noa,nob)=poa(noa,nob)+rij(nob)
      nfullpos(noa,nob,noc)=poa(noa,nob)
    enddo
  enddo
enddo
 
transpo = .true.
 
001 continue
 
if (allocated(poa)) deallocate(poa)
if (allocated(pob)) deallocate(pob)
 
END FUNCTION
 
END FUNCTION