vasp格式计算msd

分子动力学 > 分析方法 > msd

发布时间 : 2024-05-25 09:03

字数:6.4k

1. MSD 包括所有起源和一半一半和单个起源
2. NPT的MSD计算（因为每一步的晶格常数都会变，所以不能用上面的
3. pymatgen计算MSD
4. 续算的分子动力学XDATCAR处理
5. 计算单个原子的总MSD判断跳跃次数

1. MSD 包括所有起源和一半一半和单个起源
2. NPT的MSD计算（因为每一步的晶格常数都会变，所以不能用上面的
3. pymatgen计算MSD
4. 续算的分子动力学XDATCAR处理
5. 计算单个原子的总MSD判断跳跃次数

1. MSD 包括所有起源和一半一半和单个起源

==去质心漂移对氢原子扩散重要==

重要的是TMSD，不管统计多少，TMSD总体是不变的，但是MSD需要除以一个数字，就是离子数目，来作为每个离子平均的TMSD，这时候，在系数项乘数时，需要乘前面的离子数目。所以叫做均方根位移，这里的均是对有效的假定的离子数目做平均

!November 12,2015, by Peihao Huang --<<Molecular dynamics simulation_elementary methods>>-Wiley
program main
implicit none
integer::i,j,Natoms,step,atom,Origin,n,Totalstep,Nkindsofatom,Nstep,Nreject,nna  !整型
real(kind=8)
::factor,a(3,3),t,timestep,DeltaX,DeltaY,DeltaZ,DeltaRx,DeltaRy,DeltaRz,DeltaR,Xcenter0
real(kind=8)   ::Ycenter0,Zcenter0,mh,mo,m3,m4,msum,dif,logdif,sigma
!实型
real(kind=8),allocatable ::X(:,:),Y(:,:),Z(:,:),MSD(:),x0(:),y0(:),z0(:),Xcenter(:),Ycenter(:),Zcenter(:)
!实型双精度可分配数组
integer,allocatable ::natom(:)
character(len=6),allocatable ::temp(:)
!长度为6的字符串数组
character(len=10) ::temp1
mo=39.10
mh=88.91
m3=28.09
m4=16.00
msum=30*mo+6*mh+24*m3+72*m4

!These parameters should be input by hand 
Totalstep=180000
Nstep=170000
Nreject=10000
timestep=0.002
Nkindsofatom=4
allocate(natom(Nkindsofatom),temp(Nkindsofatom))
!natom为四维整型数组，temp为四维字符型数组
!Read the original structure,only the fractional coordination is needed for the box is fixed
open(unit=30,file='POSCAR')
read(30,*)
read(30,*)
read(30,*)
read(30,*)
read(30,*) ! for earlier versions of VASP, you should delete on line of this
read(30,*)
read(30,*) natom
!将原子个数存在七维实型数组里
Natoms=0
do i=1,Nkindsofatom
  Natoms=Natoms+natom(i)
  write(*,*)natom(i)
end do
!计算总的原子个数
write(*,*)Natoms
allocate(x0(Natoms),y0(Natoms),z0(Natoms))
!分配一个natoms维的数组
!read(30,*)
read(30,*)
do atom=1,Natoms
  read(30,*)x0(atom),y0(atom),z0(atom)
end do
close(30)
!将所有的初始坐标存储
!Read the structural information of the trajactory
open(unit=10,file='XDATCAR')
read(10,*)temp1
read(10,*)factor
read(10,*)a(1,1),a(1,2),a(1,3)
read(10,*)a(2,1),a(2,2),a(2,3)
read(10,*)a(3,1),a(3,2),a(3,3)
do i=1,3
  do j=1,3
    a(i,j)=a(i,j)*factor
  end do
end do
!读入晶格常数
read(10,*)temp
read(10,*)
allocate(X(Totalstep,Natoms),Y(Totalstep,Natoms),Z(Totalstep,Natoms))
!分配二维数组
do step=1,Totalstep
  read(10,*)
  do atom=1,Natoms
    read(10,*)X(step,atom),Y(step,atom),Z(step,atom)
  end do
end do
close(10)
!将所有原子坐标读入
!Unfold the atoms in the box with PBC
do atom=1,Natoms
  if(X(1,atom)-x0(atom)<-0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)+1.0
    end do
  else if(X(1,atom)-x0(atom)>0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)-1.0
    end do
  end if
  if(Y(1,atom)-y0(atom)<-0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)+1.0
    end do
  else if(Y(1,atom)-y0(atom)>0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)-1.0
    end do
  end if
  if(Z(1,atom)-z0(atom)<-0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)+1.0
    end do
  else if(Z(1,atom)-z0(atom)>0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)-1.0
    end do
  end if
end do

do step=2,Totalstep
  do atom=1,Natoms
    if(X(step,atom)-X(step-1,atom)<-0.5) then
      do n=step,Totalstep
        X(n,atom)=X(n,atom)+1.0
      end do
!从第二步开始计算，第二部比第一步多或少
    else if(X(step,atom)-X(step-1,atom)>0.5) then
      do n=step,Totalstep
        X(n,atom)=X(n,atom)-1.0
      end do
    end if
    if(Y(step,atom)-Y(step-1,atom)<-0.5) then
      do n=step,Totalstep
        Y(n,atom)=Y(n,atom)+1.0
      end do
    else if(Y(step,atom)-Y(step-1,atom)>0.5) then
      do n=step,Totalstep
        Y(n,atom)=Y(n,atom)-1.0
      end do
    end if
    if(Z(step,atom)-Z(step-1,atom)<-0.5) then
      do n=step,Totalstep
        Z(n,atom)=Z(n,atom)+1.0
      end do
    else if(Z(step,atom)-Z(step-1,atom)>0.5) then
      do n=step,Totalstep
        Z(n,atom)=Z(n,atom)-1.0
      end do
    end if
  end do
end do
!考虑周期性
!在第一帧上全部加减
!Remove the drift of center of mass
Xcenter0=0
Ycenter0=0
Zcenter0=0
do atom=1,natom(1)
  Xcenter0=Xcenter0+x0(atom)*mo/msum
  Ycenter0=Ycenter0+y0(atom)*mo/msum
  zcenter0=Zcenter0+z0(atom)*mo/msum
end do
!在natom中提取第一个原子的数量，并且计算质心
do atom=natom(1)+1,natom(1)+natom(2)
  Xcenter0=Xcenter0+x0(atom)*mh/msum
  Ycenter0=Ycenter0+y0(atom)*mh/msum
  zcenter0=Zcenter0+z0(atom)*mh/msum
end do
!提取第二个原子质量并求质心
do atom=natom(1)+natom(2)+1,natom(1)+natom(2)+natom(3)
  Xcenter0=Xcenter0+x0(atom)*m3/msum
  Ycenter0=Ycenter0+y0(atom)*m3/msum
  zcenter0=Zcenter0+z0(atom)*m3/msum
end do
do atom=natom(1)+natom(2)+natom(3)+1,natom(1)+natom(2)+natom(3)+natom(4)
   Xcenter0=Xcenter0+x0(atom)*m4/msum
   Ycenter0=Ycenter0+y0(atom)*m4/msum
   zcenter0=Zcenter0+z0(atom)*m4/msum
 end do
 !计算质心
deallocate(x0,y0,z0)

allocate(Xcenter(Totalstep),Ycenter(Totalstep),Zcenter(Totalstep))
open(unit=400,file='Mc-drift.dat')
do step=1,Totalstep
  Xcenter(step)=0
  Ycenter(step)=0
  Zcenter(step)=0
  do atom=1,natom(1)
    Xcenter(step)=Xcenter(step)+X(step,atom)*mo/msum
    Ycenter(step)=Ycenter(step)+Y(step,atom)*mo/msum
    Zcenter(step)=Zcenter(step)+Z(step,atom)*mo/msum
  end do
  do atom=natom(1)+1,natom(1)+natom(2)
    Xcenter(step)=Xcenter(step)+X(step,atom)*mh/msum
    Ycenter(step)=Ycenter(step)+Y(step,atom)*mh/msum
    Zcenter(step)=Zcenter(step)+Z(step,atom)*mh/msum
  end do
  do atom=natom(1)+natom(2)+1,natom(1)+natom(2)+natom(3)
    Xcenter(step)=Xcenter(step)+X(step,atom)*m3/msum
    Ycenter(step)=Ycenter(step)+Y(step,atom)*m3/msum
    Zcenter(step)=Zcenter(step)+Z(step,atom)*m3/msum
  end do

  do atom=natom(1)+natom(2)+natom(3)+1,natom(1)+natom(2)+natom(3)+natom(4)
    Xcenter(step)=Xcenter(step)+X(step,atom)*m4/msum
    Ycenter(step)=Ycenter(step)+Y(step,atom)*m4/msum
    Zcenter(step)=Zcenter(step)+Z(step,atom)*m4/msum
  end do
  Xcenter(step)=Xcenter(step)-Xcenter0
  Ycenter(step)=Ycenter(step)-Ycenter0
  Zcenter(step)=zcenter(step)-Zcenter0
  t=step*timestep
  write(400,"(4(F9.6(2x)))")t,Xcenter(step),Ycenter(step),Zcenter(step)
  !计算质心运动
  do atom=1,Natoms
    X(step,atom)=X(step,atom)-Xcenter(step)
    Y(step,atom)=Y(step,atom)-Ycenter(step)
    Z(step,atom)=Z(step,atom)-Zcenter(step)
  end do
  !去除质心运动
end do
deallocate(Xcenter,Ycenter,Zcenter)

!Reject some steps before equilibrium
do step=1,Nstep
  do atom=1,Natoms
    X(step,atom)=X(step+Nreject,atom)
    Y(step,atom)=Y(step+Nreject,atom)
    Z(step,atom)=Z(step+Nreject,atom)
  end do
end do
!去除部分轨迹
!Multile time origins
allocate(MSD(Nstep-1))
open(unit=20,file='msd-zong-t.dat')
do step=1,Nstep-1
  MSD(step)=0.0
  do Origin=1,Nstep-step
    do atom=1,natom(1)
      DeltaX=X(Origin+step,atom)-X(Origin,atom)
      DeltaY=Y(Origin+step,atom)-Y(Origin,atom)
      DeltaZ=Z(Origin+step,atom)-Z(Origin,atom)
      DeltaRx=DeltaX*a(1,1)+DeltaY*a(2,1)+DeltaZ*a(3,1)
      DeltaRy=DeltaX*a(1,2)+DeltaY*a(2,2)+DeltaZ*a(3,2)
      DeltaRz=DeltaX*a(1,3)+DeltaY*a(2,3)+DeltaZ*a(3,3)
      DeltaR=DeltaRx**2+DeltaRy**2+DeltaRz**2
      MSD(step)=MSD(step)+DeltaR
    end do
  end do
  MSD(step)=MSD(step)/REAL(natom(1))/REAL(Nstep-step)
  t=step*timestep
  dif=MSD(step)/(t*6)*0.0001
  logdif=log10(dif)
  sigma=30*1.6*1.6*dif/1.38/400/1796.86
!  sigma=sigma/(10**(-9))
  write(20,*)t,MSD(step)
end do
close(20)


!一半一半

open(unit=50,file='msd-zong-2-t.dat')
do step=1,Nstep/2
  MSD(step)=0.0
  do Origin=1,Nstep/2
    do atom=1,natom(1)
      DeltaX=X(Origin+step,atom)-X(Origin,atom)
      DeltaY=Y(Origin+step,atom)-Y(Origin,atom)
      DeltaZ=Z(Origin+step,atom)-Z(Origin,atom)
      DeltaRx=DeltaX*a(1,1)+DeltaY*a(2,1)+DeltaZ*a(3,1)
      DeltaRy=DeltaX*a(1,2)+DeltaY*a(2,2)+DeltaZ*a(3,2)
      DeltaRz=DeltaX*a(1,3)+DeltaY*a(2,3)+DeltaZ*a(3,3)
      DeltaR=DeltaRx**2+DeltaRy**2+DeltaRz**2
      MSD(step)=MSD(step)+DeltaR
    end do
  end do
  MSD(step)=MSD(step)*2/REAL(natom(1))/REAL(Nstep)
  t=step*timestep
  dif=MSD(step)/(t*6)*0.0001
  logdif=log10(dif)
  sigma=30*1.6*1.6*dif/1.38/400/1796.86
!  sigma=sigma/(10**(-9))
  write(50,*)t,MSD(step)
end do
close(50)



!只和第一帧比较
open(unit=60,file='msd-zong-t.dat')
do step=0,Nstep-1
  MSD(step)=0.0
  do atom=1,natom(1)
    DeltaX=X(1+step,atom)-X(1,atom)
    DeltaY=Y(1+step,atom)-Y(1,atom)
    DeltaZ=Z(1+step,atom)-Z(1,atom)
    DeltaRx=DeltaX*a(1,1)+DeltaY*a(2,1)+DeltaZ*a(3,1)
    DeltaRy=DeltaX*a(1,2)+DeltaY*a(2,2)+DeltaZ*a(3,2)
    DeltaRz=DeltaX*a(1,3)+DeltaY*a(2,3)+DeltaZ*a(3,3)
    DeltaR=DeltaRx**2+DeltaRy**2+DeltaRz**2
    MSD(step)=MSD(step)+DeltaR
  end do
  MSD(step)=MSD(step)/REAL(natom(1))
  t=step*timestep
  dif=MSD(step)/(t*6)*0.0001
  logdif=log10(dif)
  sigma=30*1.6*1.6*dif/1.38/400/1796.86
!  sigma=sigma/(10**(-9))
  write(60,*)t,MSD(step)
end do
close(60)





deallocate(MSD)
deallocate(X,Y,Z)
write(*,*)"success!"
end program

同一段轨迹的不同算法会得出不同的结果（真实VS统计结果）

平均算法会低估前面时间的msd

2. NPT的MSD计算（因为每一步的晶格常数都会变，所以不能用上面的

program main
implicit none
integer::i,j,Natoms,step,atom,Origin,n,Totalstep,Nkindsofatom,Nstep,Nreject,nna
!整型
real(kind=8)  ::factor,a(3,3),t,timestep,DeltaX,DeltaY,DeltaZ,DeltaRx,DeltaRy,DeltaRz,DeltaR,Xcenter0
real(kind=8)   ::Ycenter0,Zcenter0,mh,mo,m3,m4,msum,dif,logdif,sigma
!实型
real(kind=8),allocatable  ::X(:,:),Y(:,:),Z(:,:),MSD(:),x0(:),y0(:),z0(:),Xcenter(:),Ycenter(:),Zcenter(:),MM(:,:)
real(kind=8),allocatable  ::NN(:,:),QQ(:,:)
!实型双精度可分配数组
integer,allocatable ::natom(:)
character(len=6),allocatable ::temp(:)
!长度为6的字符串数组
character(len=10) ::temp1
mo=40.08
mh=1.008
m3=28.09
m4=16.00
msum=256*mo+1472*mh

!These parameters should be input by hand 
Totalstep=451
Nstep=400
Nreject=0
timestep=0.002
Nkindsofatom=2
allocate(natom(Nkindsofatom),temp(Nkindsofatom))
!natom为四维整型数组，temp为四维字符型数组
!Read the original structure,only the fractional coordination is needed
!for the box is fixed
open(unit=30,file='POSCAR')
read(30,*)
read(30,*)
read(30,*)a(1,1),a(1,2),a(1,3)
read(30,*)a(2,1),a(2,2),a(2,3)
read(30,*)a(3,1),a(3,2),a(3,3)! for earlier versions of VASP, you should delete on line ofthis
read(30,*)
read(30,*) natom
!将原子个数存在七维实型数组里
Natoms=0
do i=1,Nkindsofatom
  Natoms=Natoms+natom(i)
  write(*,*)natom(i)
end do
!计算总的原子个数
write(*,*)Natoms
allocate(x0(Natoms),y0(Natoms),z0(Natoms))
!分配一个natoms维的数组
!read(30,*)
read(30,*)
do atom=1,Natoms
  read(30,*)x0(atom),y0(atom),z0(atom)
end do
close(30)
!将所有的初始坐标存储
!Read the structural information of the trajactory
open(unit=10,file='XDATCAR')
allocate(X(Totalstep,Natoms),Y(Totalstep,Natoms),Z(Totalstep,Natoms))
allocate(MM(Totalstep,3),NN(Totalstep,3),QQ(Totalstep,3))
do step=1,Totalstep
  read(10,*)
  read(10,*)
  read(10,*) MM(step,1),MM(step,2),MM(step,3)
  read(10,*) NN(step,1),NN(step,2),NN(step,3)
  read(10,*) QQ(step,1),QQ(step,2),QQ(step,3)
  read(10,*)
  read(10,*)
  read(10,*)
  do atom=1,Natoms
    read(10,*)X(step,atom),Y(step,atom),Z(step,atom)
  end do
end do
close(10)
!将所有原子坐标读入
!Unfold the atoms in the box with PBC
do atom=1,Natoms
  if(X(1,atom)-x0(atom)<-0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)+1.0
    end do
  else if(X(1,atom)-x0(atom)>0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)-1.0
    end do
  end if
  if(Y(1,atom)-y0(atom)<-0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)+1.0
    end do
  else if(Y(1,atom)-y0(atom)>0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)-1.0
    end do
  end if
  if(Z(1,atom)-z0(atom)<-0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)+1.0
    end do
  else if(Z(1,atom)-z0(atom)>0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)-1.0
    end do
  end if
end do

do step=2,Totalstep
  do atom=1,Natoms
    if(X(step,atom)-X(step-1,atom)<-0.5) then
      do n=step,Totalstep
        X(n,atom)=X(n,atom)+1.0
      end do
!从第二步开始计算，第二部比第一步多或少
    else if(X(step,atom)-X(step-1,atom)>0.5) then
      do n=step,Totalstep
        X(n,atom)=X(n,atom)-1.0
      end do
    end if
    if(Y(step,atom)-Y(step-1,atom)<-0.5) then
      do n=step,Totalstep
        Y(n,atom)=Y(n,atom)+1.0
      end do
    else if(Y(step,atom)-Y(step-1,atom)>0.5) then
      do n=step,Totalstep
        Y(n,atom)=Y(n,atom)-1.0
      end do
    end if
    if(Z(step,atom)-Z(step-1,atom)<-0.5) then
      do n=step,Totalstep
        Z(n,atom)=Z(n,atom)+1.0
      end do
    else if(Z(step,atom)-Z(step-1,atom)>0.5) then
      do n=step,Totalstep
        Z(n,atom)=Z(n,atom)-1.0
      end do
    end if
  end do
end do
open(unit=100,file='XDATCAR-3')
  do step=1,Nstep
   write(100,*) step,"sucess"
  do atom =1,Natoms
 write(100,"(F11.8,F15.8,F15.8)") X(step,atom),Y(step,atom),Z(step,atom)
  end do
  end do
close(100)


do step=1,Totalstep
   do atom=1,Natoms
     X(step,atom)=MM(step,1)*X(step,atom)+NN(step,1)*Y(step,atom)+QQ(step,1)*Z(step,atom)
     Y(step,atom)=MM(step,2)*X(step,atom)+NN(step,2)*Y(step,atom)+QQ(step,2)*Z(step,atom)
     Z(step,atom)=MM(step,3)*X(step,atom)+NN(step,3)*Y(step,atom)+QQ(step,3)*Z(step,atom)
   end do
end do

do atom=1,Natoms
  x0(atom)=a(1,1)*x0(atom)
  y0(atom)=a(2,2)*y0(atom)
  z0(atom)=a(3,3)*z0(atom)
end do



open(unit=100,file='XDATCAR-2')
  do step=1,Nstep
   write(100,*) step ,"success"
  do atom =1,Natoms
 write(100,"(F11.8,F15.8,F15.8)") X(step,atom),Y(step,atom),Z(step,atom)
  end do
  end do
close(100)

!考虑周期性
!在第一帧上全部加减
!Remove the drift of center of mass
Xcenter0=0
Ycenter0=0
Zcenter0=0
do atom=1,natom(1)
  Xcenter0=Xcenter0+x0(atom)*mo/msum
  Ycenter0=Ycenter0+y0(atom)*mo/msum
  zcenter0=Zcenter0+z0(atom)*mo/msum
end do
!在natom中提取第一个原子的数量，并且计算质心
do atom=natom(1)+1,natom(2)
  Xcenter0=Xcenter0+x0(atom)*mh/msum
  Ycenter0=Ycenter0+y0(atom)*mh/msum
  zcenter0=Zcenter0+z0(atom)*mh/msum
end do
!提取第二个原子质量并求质心
 !计算质心
deallocate(x0,y0,z0)

allocate(Xcenter(Totalstep),Ycenter(Totalstep),Zcenter(Totalstep))
open(unit=400,file='Mc-drift.dat')
do step=1,Totalstep
  Xcenter(step)=0
  Ycenter(step)=0
  Zcenter(step)=0
  do atom=1,natom(1)
    Xcenter(step)=Xcenter(step)+X(step,atom)*mo/msum
    Ycenter(step)=Ycenter(step)+Y(step,atom)*mo/msum
    Zcenter(step)=Zcenter(step)+Z(step,atom)*mo/msum
  end do
  do atom=natom(1)+1,natom(1)+natom(2)
    Xcenter(step)=Xcenter(step)+X(step,atom)*mh/msum
    Ycenter(step)=Ycenter(step)+Y(step,atom)*mh/msum
    Zcenter(step)=Zcenter(step)+Z(step,atom)*mh/msum
  end do

  Xcenter(step)=Xcenter(step)-Xcenter0
  Ycenter(step)=Ycenter(step)-Ycenter0
  Zcenter(step)=zcenter(step)-Zcenter0
  t=step*timestep
  write(400,"(4(F9.6(2x)))")t,Xcenter(step),Ycenter(step),Zcenter(step)
  !计算质心运动
  do atom=1,Natoms
    X(step,atom)=X(step,atom)-Xcenter(step)
    Y(step,atom)=Y(step,atom)-Ycenter(step)
    Z(step,atom)=Z(step,atom)-Zcenter(step)
  end do
  !去除质心运动
end do
deallocate(Xcenter,Ycenter,Zcenter)

!Reject some steps before equilibrium
do step=1,Nstep
  do atom=1,Natoms
    X(step,atom)=X(step+Nreject,atom)
    Y(step,atom)=Y(step+Nreject,atom)
    Z(step,atom)=Z(step+Nreject,atom)
  end do
end do
!去除部分轨迹
!Multile time origins
allocate(MSD(Nstep-1))
open(unit=20,file='msd-zong-t-2.dat')
do step=1,Nstep-1
  MSD(step)=0.0
  do Origin=1,Nstep-step
    do atom=natom(1)+1,natom(1)+natom(2)
      DeltaX=X(Origin+step,atom)-X(Origin,atom)
      DeltaY=Y(Origin+step,atom)-Y(Origin,atom)
      DeltaZ=Z(Origin+step,atom)-Z(Origin,atom)
      DeltaR=DeltaX**2+DeltaY**2+DeltaZ**2
      MSD(step)=MSD(step)+DeltaR
    end do
  end do
  MSD(step)=MSD(step)/REAL(natom(2))/REAL(Nstep-step)
  t=step*timestep
  dif=MSD(step)/(t*6)*0.0001
  logdif=log10(dif)
  sigma=30*1.6*1.6*dif/1.38/400/1796.86
!  sigma=sigma/(10**(-9))
  write(20,*)t,MSD(step)
end do
close(20)
!
!
!!一半一半
!
!open(unit=50,file='msd-zong-2-t.dat')
!do step=1,Nstep/2
!  MSD(step)=0.0
!  do Origin=1,Nstep/2
!    do atom=1,natom(1)
!      DeltaX=X(Origin+step,atom)-X(Origin,atom)
!      DeltaY=Y(Origin+step,atom)-Y(Origin,atom)
!      DeltaZ=Z(Origin+step,atom)-Z(Origin,atom)
!      DeltaRx=DeltaX*a(1,1)+DeltaY*a(2,1)+DeltaZ*a(3,1)
!      DeltaRy=DeltaX*a(1,2)+DeltaY*a(2,2)+DeltaZ*a(3,2)
!      DeltaRz=DeltaX*a(1,3)+DeltaY*a(2,3)+DeltaZ*a(3,3)
!      DeltaR=DeltaRx**2+DeltaRy**2+DeltaRz**2
!      MSD(step)=MSD(step)+DeltaR
!    end do
!  end do
!  MSD(step)=MSD(step)*2/REAL(natom(1))/REAL(Nstep)
!  t=step*timestep
!  dif=MSD(step)/(t*6)*0.0001
!  logdif=log10(dif)
!  sigma=30*1.6*1.6*dif/1.38/400/1796.86
!!  sigma=sigma/(10**(-9))
!  write(50,*)t,MSD(step)
!end do
!close(50)



!只和第一帧比较
open(unit=60,file='msd-zong-t.dat')
do step=0,Nstep-1
  MSD(step)=0.0
  do atom=natom(1),natom(1)+natom(2)
    DeltaX=X(1+step,atom)-X(1,atom)
    DeltaY=Y(1+step,atom)-Y(1,atom)
    DeltaZ=Z(1+step,atom)-Z(1,atom)
    DeltaR=DeltaX**2+DeltaY**2+DeltaZ**2
    MSD(step)=MSD(step)+DeltaR
  end do
  MSD(step)=MSD(step)/REAL(natom(2))
  t=step*timestep
  dif=MSD(step)/(t*6)*0.0001
  logdif=log10(dif)
  sigma=30*1.6*1.6*dif/1.38/400/1796.86
!  sigma=sigma/(10**(-9))
  write(60,*)t,MSD(step)
end do
close(60)





deallocate(MSD)
deallocate(X,Y,Z)
write(*,*)"success!"
end program

3. pymatgen计算MSD

import os
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen.core import Structure
from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer
import numpy as np
import pickle
from pymatgen.analysis.diffusion_analyzer import get_conversion_factor
traj=Trajectory.from_file('XDATCAR')
diff=DiffusionAnalyzer.from_structures(traj,'H',600,1,1,smoothed='constant',avg_nsteps=1)
#转换系数
conv_factor = get_conversion_factor(diff.structure, diff.specie, diff.temperature)
with open('conv_factor.dat','w') as f:
    f.write('conv_factor\t%f' %(conv_factor))
#get summary
summary = diff.get_summary_dict()
import json
json.dump(summary,open('summary.json','w'))
#msd数据
diff.export_msdt('msddata.dat')
#mscd数据
#mscd = diff.mscd
#with open('mscd.dat','w') as f:
#    for row in mscd :
#       f.write(' '.join([str(row)]))
#       f.write('\n')

4. 续算的分子动力学XDATCAR处理

!November 12,2015, by Peihao Huang --<<Molecular dynamics simulation_elementary methods>>-Wiley
program main
implicit none
integer::i,j,Natoms,step,atom,Origin,n,Totalstep,Nkindsofatom,Nstep,Nreject  !整型
real(kind=8)             ::factor,a(3,3),t,timestep,DeltaX,DeltaY,DeltaZ,DeltaRx,DeltaRy,DeltaRz,DeltaR,Xcenter0
real(kind=8)             ::Ycenter0,Zcenter0,mh,mo,m3,m4,msum,dif,logdif,sigma
!实型
real(kind=8),allocatable ::X(:,:),Y(:,:),Z(:,:),MSD(:),x0(:),y0(:),z0(:),Xcenter(:),Ycenter(:),Zcenter(:)
!实型双精度可分配数组
integer,allocatable ::natom(:)
character(len=6),allocatable ::temp(:)
!长度为6的字符串数组
character(len=7) ::temp1
mo=6.94
mh=69.72
m3=78.96
m4=16
msum=8*mo+8*mh+16*m3+48*m4

!These parameters should be input by hand 
Totalstep=187554
Nstep=187554
Nreject=0
timestep=0.002
Nkindsofatom=4

allocate(natom(Nkindsofatom),temp(Nkindsofatom))
!natom为四维整型数组，temp为四维字符型数组
!Read the original structure,only the fractional coordination is needed for the box is fixed
open(unit=30,file='POSCAR')
read(30,*)
read(30,*)
read(30,*)
read(30,*)
read(30,*) ! for earlier versions of VASP, you should delete on line of this
read(30,*)
read(30,*) natom
!将原子个数存在四维实型数组里
Natoms=0
do i=1,Nkindsofatom
  Natoms=Natoms+natom(i)
end do
!计算总的原子个数
write(*,*)Natoms
allocate(x0(Natoms),y0(Natoms),z0(Natoms))
!分配一个natoms维的数组
!read(30,*)
read(30,*)
do atom=1,Natoms
  read(30,*)x0(atom),y0(atom),z0(atom)
end do
close(30)
!将所有的初始坐标存储
!Read the structural information of the trajactory
open(unit=10,file='XDATCAR')
read(10,*)temp1
read(10,*)factor
read(10,*)a(1,1),a(1,2),a(1,3)
read(10,*)a(2,1),a(2,2),a(2,3)
read(10,*)a(3,1),a(3,2),a(3,3)
do i=1,3
  do j=1,3
    a(i,j)=a(i,j)*factor
  end do
end do
!读入晶格常数
read(10,*)temp
read(10,*)temp
allocate(X(Totalstep,Natoms),Y(Totalstep,Natoms),Z(Totalstep,Natoms))
!分配二维数组
do step=1,Totalstep
  read(10,*)
  do atom=1,Natoms
    read(10,*)X(step,atom),Y(step,atom),Z(step,atom)
  end do
end do
close(10)
!将所有原子坐标读入
!Unfold the atoms in the box with PBC
do atom=1,Natoms
  if(X(1,atom)-x0(atom)<-0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)+1.0
    end do
  else if(X(1,atom)-x0(atom)>0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)-1.0
    end do
  end if
  if(Y(1,atom)-y0(atom)<-0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)+1.0
    end do
  else if(Y(1,atom)-y0(atom)>0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)-1.0
    end do
  end if
  if(Z(1,atom)-z0(atom)<-0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)+1.0
    end do
  else if(Z(1,atom)-z0(atom)>0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)-1.0
    end do
  end if
end do

open(unit=100,file='XDATCAR-1')
  do step=1,Nstep
   write(100,*) step
  do atom =1,Natoms
 write(100,"(F11.8,F15.8,F15.8)") X(step,atom),Y(step,atom),Z(step,atom)
  end do
  end do
close(100)


do atom=1,Natoms
  if(X(1,atom)-x0(atom)<-0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)+1.0
    end do
  else if(X(1,atom)-x0(atom)>0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)-1.0
    end do
  end if
  if(Y(1,atom)-y0(atom)<-0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)+1.0
    end do
  else if(Y(1,atom)-y0(atom)>0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)-1.0
    end do
  end if
  if(Z(1,atom)-z0(atom)<-0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)+1.0
    end do
  else if(Z(1,atom)-z0(atom)>0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)-1.0
    end do
  end if
end do

open(unit=100,file='XDATCAR-2')
  do step=1,Nstep
   write(100,*) step
  do atom =1,Natoms
 write(100,"(F11.8,F15.8,F15.8)") X(step,atom),Y(step,atom),Z(step,atom)
  end do
  end do
close(100)

do atom=1,Natoms
  if(X(1,atom)-x0(atom)<-0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)+1.0
    end do
  else if(X(1,atom)-x0(atom)>0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)-1.0
    end do
  end if
  if(Y(1,atom)-y0(atom)<-0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)+1.0
    end do
  else if(Y(1,atom)-y0(atom)>0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)-1.0
    end do
  end if
  if(Z(1,atom)-z0(atom)<-0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)+1.0
    end do
  else if(Z(1,atom)-z0(atom)>0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)-1.0
    end do
  end if
end do

open(unit=100,file='XDATCAR-3')
  do step=1,Nstep
   write(100,*) step
  do atom =1,Natoms
 write(100,"(F11.8,F15.8,F15.8)") X(step,atom),Y(step,atom),Z(step,atom)
  end do
  end do
close(100)

do atom=1,Natoms
  if(X(1,atom)-x0(atom)<-0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)+1.0
    end do
  else if(X(1,atom)-x0(atom)>0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)-1.0
    end do
  end if
  if(Y(1,atom)-y0(atom)<-0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)+1.0
    end do
  else if(Y(1,atom)-y0(atom)>0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)-1.0
    end do
  end if
  if(Z(1,atom)-z0(atom)<-0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)+1.0
    end do
  else if(Z(1,atom)-z0(atom)>0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)-1.0
    end do
  end if
end do

open(unit=100,file='XDATCAR-4')
  do step=1,Nstep
   write(100,*) step
  do atom =1,Natoms
 write(100,"(F11.8,F15.8,F15.8)") X(step,atom),Y(step,atom),Z(step,atom)
  end do
  end do
close(100)

deallocate(X,Y,Z)
write(*,*)"success!"
end program

5. 计算单个原子的总MSD判断跳跃次数

楼下 /work/home/liz/workspace/1-system-MD/1-Na5YSi4O12/9a-smass=2potim=1/600/4-msd-for-every-atom/msd-atom.f90

!November 12,2015, by Peihao Huang --<<Molecular dynamics simulation_elementary methods>>-Wiley
program main
implicit none
integer::i,j,Natoms,step,atom,Origin,n,Totalstep,Nkindsofatom,Nstep,Nreject,nna  !整型
real(kind=8）::factor,a(3,3),t,timestep,DeltaX,DeltaY,DeltaZ,DeltaRx,DeltaRy,DeltaRz,DeltaR,Xcenter0
real(kind=8) ::Ycenter0,Zcenter0,mh,mo,m3,m4,msum,dif,logdif,sigma
!实型
real(kind=8),allocatable::X(:,:),Y(:,:),Z(:,:),MSD(:),x0(:),y0(:),z0(:),Xcenter(:),Ycenter(:),Zcenter(:)
!实型双精度可分配数组
integer,allocatable ::natom(:)
character(len=6),allocatable ::temp(:)
!长度为6的字符串数组
character(len=10) ::temp1
character(len=20) :: filename

!These parameters should be input by hand 
Totalstep=270000
Nstep=270000
Nreject=0
timestep=0.002
Nkindsofatom=4

allocate(natom(Nkindsofatom),temp(Nkindsofatom))
!natom为四维整型数组，temp为四维字符型数组
!Read the original structure,only the fractional coordination is needed for the box is fixed
open(unit=30,file='POSCAR')
read(30,*)
read(30,*)
read(30,*)
read(30,*)
read(30,*) ! for earlier versions of VASP, you should delete on line of this
read(30,*)
read(30,*) natom
!将原子个数存在七维实型数组里
Natoms=0
do i=1,Nkindsofatom
  Natoms=Natoms+natom(i)
end do
!计算总的原子个数
write(*,*)Natoms
allocate(x0(Natoms),y0(Natoms),z0(Natoms))
!分配一个natoms维的数组
!read(30,*)
read(30,*)
do atom=1,Natoms
  read(30,*)x0(atom),y0(atom),z0(atom)
end do
close(30)
!将所有的初始坐标存储
!Read the structural information of the trajactory
open(unit=10,file='XDATCAR')
read(10,*)temp1
read(10,*)factor
read(10,*)a(1,1),a(1,2),a(1,3)
read(10,*)a(2,1),a(2,2),a(2,3)
read(10,*)a(3,1),a(3,2),a(3,3)
do i=1,3
  do j=1,3
    a(i,j)=a(i,j)*factor
  end do
end do
!读入晶格常数
read(10,*)temp
read(10,*)
allocate(X(Totalstep,Natoms),Y(Totalstep,Natoms),Z(Totalstep,Natoms))
!分配二维数组
do step=1,Totalstep
  read(10,*)
  do atom=1,Natoms
    read(10,*)X(step,atom),Y(step,atom),Z(step,atom)
  end do
end do
close(10)
!将所有原子坐标读入
!Unfold the atoms in the box with PBC
do atom=1,Natoms
  if(X(1,atom)-x0(atom)<-0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)+1.0
    end do
  else if(X(1,atom)-x0(atom)>0.5) then
    do step=1,Totalstep
      X(step,atom)=X(step,atom)-1.0
    end do
  end if
  if(Y(1,atom)-y0(atom)<-0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)+1.0
    end do
  else if(Y(1,atom)-y0(atom)>0.5) then
    do step=1,Totalstep
      Y(step,atom)=Y(step,atom)-1.0
    end do
  end if
  if(Z(1,atom)-z0(atom)<-0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)+1.0
    end do
  else if(Z(1,atom)-z0(atom)>0.5) then
    do step=1,Totalstep
      Z(step,atom)=Z(step,atom)-1.0
    end do
  end if
end do

do step=2,Totalstep
  do atom=1,Natoms
    if(X(step,atom)-X(step-1,atom)<-0.5) then
      do n=step,Totalstep
        X(n,atom)=X(n,atom)+1.0
      end do
!从第二步开始计算，第二部比第一步多或少
    else if(X(step,atom)-X(step-1,atom)>0.5) then
      do n=step,Totalstep
        X(n,atom)=X(n,atom)-1.0
      end do
    end if
    if(Y(step,atom)-Y(step-1,atom)<-0.5) then
      do n=step,Totalstep
        Y(n,atom)=Y(n,atom)+1.0
      end do
    else if(Y(step,atom)-Y(step-1,atom)>0.5) then
      do n=step,Totalstep
        Y(n,atom)=Y(n,atom)-1.0
      end do
    end if
    if(Z(step,atom)-Z(step-1,atom)<-0.5) then
      do n=step,Totalstep
        Z(n,atom)=Z(n,atom)+1.0
      end do
    else if(Z(step,atom)-Z(step-1,atom)>0.5) then
      do n=step,Totalstep
        Z(n,atom)=Z(n,atom)-1.0
      end do
    end if
  end do
end do
!考虑周期性
!在第一帧上全部加减
!Remove the drift of center of mass

!Reject some steps before equilibrium
do step=1,Nstep
  do atom=1,Natoms
    X(step,atom)=X(step+Nreject,atom)
    Y(step,atom)=Y(step+Nreject,atom)
    Z(step,atom)=Z(step+Nreject,atom)
  end do
end do
!去除部分轨迹
!Multile time origins
allocate(MSD(Nstep-1))
!open(unit=20,file='msd-zong-t.dat')
!do step=1,Nstep-1
!  MSD(step)=0.0
!  atom=1
!  DeltaX=X(step,atom)-X(1,atom)
!  DeltaY=Y(step,atom)-Y(1,atom)
!  DeltaZ=Z(step,atom)-Z(1,atom)
!  DeltaRx=DeltaX*a(1,1)+DeltaY*a(2,1)+DeltaZ*a(3,1)
!  DeltaRy=DeltaX*a(2,2)+DeltaY*a(2,2)+DeltaZ*a(3,2)
!  DeltaRz=DeltaX*a(1,3)+DeltaY*a(2,3)+DeltaZ*a(3,3)
!  DeltaR=DeltaRx**2+DeltaRy**2+DeltaRz**2
!  MSD(step)=MSD(step)+DeltaR
!    
!  t=step*timestep
!  dif=MSD(step)/(t*6)*0.0001
!  logdif=log10(dif)
!  sigma=30*1.6*1.6*dif/1.38/400/1796.86
!!  sigma=sigma/(10**(-9))
!  write(20,*)t,MSD(step)
!end do
!close(20)
do atom=1, 30
  write(filename, "(A,I0,A)") "msd-", atom, ".dat"
  open(unit=20, file=trim(filename))
  do step=1, Nstep-1
    MSD(step) = 0.0
    DeltaX = X(step,atom) - X(1,atom)
    DeltaY = Y(step,atom) - Y(1,atom)
    DeltaZ = Z(step,atom) - Z(1,atom)
    DeltaRx = DeltaX * a(1,1) + DeltaY * a(2,1) + DeltaZ * a(3,1)
    DeltaRy = DeltaX * a(1,2) + DeltaY * a(2,2) + DeltaZ * a(3,2)
    DeltaRz = DeltaX * a(1,3) + DeltaY * a(2,3) + DeltaZ * a(3,3)
    DeltaR = DeltaRx**2 + DeltaRy**2 + DeltaRz**2
    MSD(step) = MSD(step) + DeltaR

    t = step * timestep
    dif = MSD(step) / (t * 6) * 0.0001
    logdif = log10(dif)
    sigma = 30 * 1.6 * 1.6 * dif / 1.38 / 400 / 1796.86
    !  sigma = sigma / (10**(-9))
    write(20, *) t, MSD(step)
  end do
  close(20)
end do

deallocate(MSD)
deallocate(X,Y,Z)
write(*,*)"success!"
end program

转载请注明来源有问题可通过github提交issue