Calculate MSD from unwraped XYZ file (Python/Fortran)

Unwraped xyz file can be extracted from ase trajectory using the following python script:

import ase.io
# save traj in .xyz --- the .traj is a non human readable database file
out_traj = ase.io.read('test.traj', ':')
for at in out_traj:
#    at.wrap()    #! This function wraps the atoms into primitive box for each frame
    if 'momenta' in at.arrays: del at.arrays['momenta']
ase.io.write('test.xyz', out_traj, 'xyz')

If you are analysing the XDATCAR (VASP) trajectory, use this python script to convert it into unwraped xyz file.

Usage: Python filename.py

You should have XDATCAR file in the working path.

#lipai@mail.ustc.edu.cn
#convert XDATCAR to unwraped xyz file
import numpy as np
from copy import deepcopy

xdatcar = open('XDATCAR', 'r')
xyz = open('XDATCAR.xyz', 'w')

system = xdatcar.readline()
scale = float(xdatcar.readline().rstrip('\n'))
print(scale)

#get lattice vectors
a1 = np.array([ float(s)*scale for s in xdatcar.readline().rstrip('\n').split() ])
a2 = np.array([ float(s)*scale for s in xdatcar.readline().rstrip('\n').split() ])
a3 = np.array([ float(s)*scale for s in xdatcar.readline().rstrip('\n').split() ])

comment='Lattice=\"'+str(a1[0])+' '+str(a1[1])+' '+str(a1[2])
comment=comment+str(a2[0])+' '+str(a2[1])+' '+str(a2[2])
comment=comment+str(a3[0])+' '+str(a3[1])+' '+str(a3[2])+'\"'

#Read xdatcar
element_names = xdatcar.readline().rstrip('\n').split()

element_dict = {}
element_numbers = xdatcar.readline().rstrip('\n').split()

Natom = 0
Ntype = len(element_names)
Nname=[]
for t in range(Ntype):
    Natom += int(element_numbers[t])
    for i in range(int(element_numbers[t])):
        Nname.append(element_names[t])

print(Ntype,Natom)

f_prev=np.zeros([Natom,3])
f_next=np.zeros([Natom,3])
while True:
    line = xdatcar.readline()
    if len(line) == 0:
        break
    xyz.write(str(Natom) + "\n"+comment+"\n")

    for atom in range(Natom):
        p = xdatcar.readline().rstrip('\n').split()
        f_next[atom,:] = np.array([ float(s) for s in p ])
        f_next[atom,:] -= np.around(f_next[atom,:]-f_prev[atom,:])

        c_coords=f_next[atom,0]*a1+f_next[atom,1]*a2+f_next[atom,2]*a3
        xyz.write(Nname[atom]+" "+str(c_coords[0])+" "+str(c_coords[1])+" "+str(c_coords[2])+"\n")

    f_prev=deepcopy(f_next)

xdatcar.close()
xyz.close()

The code for calculating the MSD from xyz file with full sampling is as follows (Recommended):

Usage: ./a.out filename.xyz Nind

which calculates the MSD of the Nind-th atom. There fore Nind should be a positive integer and not larger than the total number of atoms. If one set Nind to be 0, it calculates MSD of all atoms and save msd results in the folder msd_all.

!lipai@mail.ustc.edu.cn
program xyz2msd
implicit none
integer :: ArgNum,i,j,k
character(len=20) :: filename,str1
integer :: Nind,AtomNum,FrameNum,Nline

real,dimension(3,3)::Lattice

real,dimension(:,:,:),allocatable :: xyz
real,dimension(:),allocatable   :: msd
real,dimension(:,:),allocatable :: msd_i
character(2) :: temp
integer :: stat1

ArgNum=command_argument_count()
if(ArgNum/=2) then
   write(*,*) "Wrong input argument!"
   stop
end if
call get_command_argument(1,filename)
filename=trim(adjustl(filename))
call get_command_argument(2,str1)
read(str1,*) Nind

open(unit=111,file=filename,status='old')

read(111,*) AtomNum
Nline=1
do
  read(111,*,iostat=stat1)
  Nline=Nline+1
  if(stat1/=0) exit
end do
FrameNum=Nline/(AtomNum+2)
write(*,*) AtomNum,FrameNum
allocate(xyz(FrameNum,AtomNum,3))

rewind(111)
do i=1,FrameNum
  read(111,*)
  read(111,*)
  do j=1,AtomNum
    read(111,*) temp, xyz(i,j,:)
  end do
end do
close(111)

allocate(msd(FrameNum),msd_i(FrameNum,3))
if(Nind==0) then
   call system("mkdir msd_all")
   do i=1,AtomNum
       call get_msd_ind(i,pre="msd_all/")
   end do
else
   if(Nind>AtomNum) then
       write(*,*) "Error! Nind>AtomNum"
       stop
   end if
   call get_msd_ind(Nind)
end if

contains
subroutine get_msd_ind(ind,pre)
implicit none
integer,intent(in) :: ind
character(*),optional,intent(in)  :: pre
character(20) :: outname
integer :: del_t,n,m,Fd
write(outname,*) ind
if(present(pre)) then
    outname=pre//trim(adjustl(outname))
else
    outname=trim(adjustl(outname))
end if
write(*,*) "Dealing with "//outname//".msd"
write(*,*) ind
msd=0
msd_i=0
do del_t=1,FrameNum-1
  Fd=FrameNum-del_t
  if(mod(del_t,100)==0) write(*,*) del_t
  do n=1,Fd
    do m=1,3
        msd_i(del_t,m)=msd_i(del_t,m)+((xyz(n+del_t,ind,m)-xyz(n,ind,m))**2)/Fd
    end do
  end do
  do m=1,3
      msd(del_t)=msd(del_t)+msd_i(del_t,m)
  end do
end do
open(99,file=trim(adjustl(outname))//".msd",status="replace")
open(98,file=trim(adjustl(outname))//".msd_i",status="replace")
do del_t=1,framenum-1
    write(99,*) msd(del_t)
    write(98,*) msd_i(del_t,:)
end do
close(99)
close(98)
end subroutine

end program

You may also try the following code which takes lease sampling. Its algorithm is the same with thatprovided by Muratahan Aykol, whose script can also be found in this website.

!lipai@mail.ustc.edu.cn
program xyz2msd
implicit none
integer :: ArgNum,i,j,k
character(len=20) :: filename,str1
integer :: Nind,AtomNum,FrameNum

real,dimension(3,3)::Lattice

real,dimension(:,:,:),allocatable :: xyz
real,dimension(:),allocatable   :: msd
real,dimension(:,:),allocatable :: msd_i
character(2) :: temp
integer :: stat1

ArgNum=command_argument_count()
if(ArgNum/=2) then
   write(*,*) "Wrong input argument!"
   stop
end if
call get_command_argument(1,filename)
filename=trim(adjustl(filename))
call get_command_argument(2,str1)
read(str1,*) Nind

open(unit=111,file=filename,status='old')

read(111,*) AtomNum
do
  read(111,*,iostat=stat1)
  Nline=Nline+1
  if(stat1/=0) exit
end do
FrameNum=Nline/(AtomNum+2)
write(*,*) AtomNum,FrameNum
allocate(xyz(FrameNum,AtomNum,3))

rewind(111)
do i=1,FrameNum
  read(111,*)
  read(111,*)
  do j=1,AtomNum
    read(111,*) temp, xyz(i,j,:)
  end do
end do
close(111)

allocate(msd(FrameNum),msd_i(FrameNum,3))
if(Nind==0) then
   call system("mkdir msd_all")
   do i=1,AtomNum
       call get_msd_ind(i,pre="msd_all/")
   end do
else
   if(Nind>AtomNum) then
       write(*,*) "Error! Nind>AtomNum"
       stop
   end if
   call get_msd_ind(Nind)
end if

contains
subroutine get_msd_ind(ind,pre)
implicit none
integer,intent(in) :: ind
character(*),optional,intent(in)  :: pre
character(20) :: outname
integer :: del_t,n,m,Fd
write(outname,*) ind
if(present(pre)) then
    outname=pre//trim(adjustl(outname))
else
    outname=trim(adjustl(outname))
end if
write(*,*) "Dealing with "//outname//".msd"
write(*,*) ind
msd=0
msd_i=0
do del_t=1,FrameNum-1
  Fd=FrameNum-del_t
  if(mod(del_t,100)==0) write(*,*) del_t
  do m=1,3
      msd_i(del_t,m)=(xyz(del_t+1,ind,m)-xyz(1,ind,m))**2
  end do
  do m=1,3
      msd(del_t)=msd(del_t)+msd_i(del_t,m)
  end do
end do
open(99,file=trim(adjustl(outname))//".msd",status="replace")
open(98,file=trim(adjustl(outname))//".msd_i",status="replace")
do del_t=1,framenum-1
    write(99,*) msd(del_t)
    write(98,*) msd_i(del_t,:)
end do
close(99)
close(98)
end subroutine

end program