两种状态下的扩散系数计算公式（扩散系数的计算公式） - 原点资讯

之前有朋友在公众号后台留言问VASP怎么做AIMD。笔者这里用VASP官网做Si不同温度下扩散系数的例子举例。扩散系数的计算理论作者之前在介绍用MD程序LAMMPS计算MSD的例子中已经讲过，可以参考一下。

1. 静力学结构优化(略)

2. 不同温度下的AIMD计算

3. 数据处理

第2步自动修改并提交不同温度下计算AIMDINCAR的shell脚本如下：

for i in 2000 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800docat >INCAR <<!SYSTEM = Si# electronic degrees LREAL = A # real space projectionPREC = Normal # chose Low only after testsEDIFF = 1E-5 # do not use default (too large drift)ISMEAR = -1 ; SIGMA = 0.130 # Fermi smearing: 1500 K 0.086 10-3ALGO = Very Fast # recommended for MD (fall back ALGO = Fast)MAXMIX = 40 # reuse mixer from one MD step to nextISYM = 0 # no symmetry NELMIN = 4 # minimum 4 steps per time step, avoid breaking after 2 steps# MD (do little writing to save disc space)IBRION = 0 # main molecular dynamics tagNSW = 400 # number of MD stepsPOTIM = 3 # time step of MDNWRITE = 0 # controls outputNBLOCK = 10 # after ten steps pair correlation function is written outLCHARG = .FALSE. # no charge density written outLWAVE = .FALSE. # no wave function coefficients written outTEBEG = $i # starting temperature for MDTEEND = $i # end temperature for MD# canonic (Nose) MD with XDATCAR updated every 10 stepsMDALGO = 2 ä switch to select thermostatSMASS = 3 # Nose massISIF = 2 # this tag selects the ensemble in combination with the thermostat!mpirun -np 2 /path/to/your/vasp/executablecp XDATCAR XDATCAR.$icp OUTCAR OUTCAR.$icp PCDAT PCDAT.$icp CONTCAR CONTCAR.$icp POSCAR POSCAR.$icp OSZICAR OSZICAR.$icp CONTCAR POSCARdone

计算完成后会生成一系列XDATCAR_温度文件，改文件为VASP计算过程中的轨迹文件，用它就可以处理MSD并拟合扩散系数。处理的python程序diffusion_coefficient.py 如下：

#!/usr/bin/python
import sysimport reimport math
#setting grid for histogram
potim = 3 #timestep from INCAR filereadfile = open(sys.argv[1],"r") #input XDATCAR file in format XDATCAR.TEMPtemp=re.sub("XDATCAR.",,sys.argv[1]) #extracts temperature from input file namez=0 #counternatoms=0 #number of atoms in XDATCAR fileposion = #atom positions in Cartesian coordinatesconfcount = 0 #number of structures in XDATCAR filedirect= #number of time steps for each structure in XDATCAR filea= #lattice parameter in 1st dimensionb= #lattice parameter in 2nd dimensionc= #lattice parameter in 3rd dimension#read in XDATCAR fileline=readfile.readlinewhile (line): z=z 1 line.strip line=re.sub('^',' ',line) y=line.split if (z==2): scale=float(y[0]) if (z==3): a.append(float(y[0])) a.append(float(y[1])) a.append(float(y[2])) a_len=(a[0]*a[0] a[1]*a[1] a[2]*a[2])**0.5 if (z==4): b.append(float(y[0])) b.append(float(y[1])) b.append(float(y[2])) b_len=(b[0]*b[0] b[1]*b[1] b[2]*b[2])**0.5 if (z==5): c.append(float(y[0])) c.append(float(y[1])) c.append(float(y[2])) c_len=(c[0]*c[0] c[1]*c[1] c[2]*c[2])**0.5 if (z==7): natoms=int(y[0]) if (y[0]=="Direct"): direct.append(int(y[2])) posion.append([]) for i in range(0,natoms): line=readfile.readline line.strip line=re.sub('^',' ',line) f=line.split cartpos_x=a[0]*float(f[0]) a[1]*float(f[1]) a[2]*float(f[2]) cartpos_y=b[0]*float(f[0]) b[1]*float(f[1]) b[2]*float(f[2]) cartpos_z=c[0]*float(f[0]) c[1]*float(f[1]) c[2]*float(f[2]) #positions of ions for each structure are obtained here posion[confcount].append([cartpos_x,cartpos_y,cartpos_z]) confcount=confcount 1 line=readfile.readlinereadfile.close
#calculate diffusion coefficient#skip first 10 configurations corresponding to 300 fsd=0.0for i in range(10,confcount): for j in range(0,natoms): x_diff=posion[i][j][0]-posion[0][j][0] #if length is larger than 0.5 (in crystallographic coordinates) then we have to shift atom #due to periodic image to obtain the shortest distance. if (abs(x_diff)>(0.5*a_len)): if (x_diff<0): x_diff=x_diff a_len elif (x_diff>0): x_diff=x_diff-a_len y_diff=posion[i][j][1]-posion[0][j][1] if (abs(y_diff)>(0.5*b_len)): if (y_diff<0): y_diff=y_diff b_len elif (y_diff>0): y_diff=y_diff-b_len z_diff=posion[i][j][2]-posion[0][j][2] if (abs(z_diff)>(0.5*c_len)): if (z_diff<0): z_diff=z_diff c_len elif (x_diff>0): z_diff=z_diff-c_len d=d x_diff**2.0 y_diff**2.0 z_diff**2.0
#print diffusion coefficient (in Ang^2/ps) vs temperature (in K)d=d/(confcount-1-10)/natoms/6.0time=(direct[confcount-1]-direct[10])*potim/10**3.0 #conversion to psprint temp,d/time

运行方式为：

python diffusion_coefficient.py XDATCAR_温度值

批量处理并绘制扩散系数与温度关系的shell脚本

#!/bin/bash
if test -f "diff_coeff.dat"; then rm diff_coeff.datfi
touch diff_coeff.dat
for i in 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000; do diffusion_coefficient.py XDATCAR.$i >> diff_coeff.datdone
gnuplot -e "set terminal jpeg; set key left; set xlabel 'temperature (K)'; set ylabel 'D (Ang^2/ps)'; set style data lines; plot 'diff_coeff.dat' " > diff_coeff.jpg

运行结果如下

两种状态下的扩散系数计算公式,扩散系数的计算公式(1)