文章目录
- 1.列举常用的lammps自带函数
- 2.如何进行应力-应变的计算及绘制应力云图
- 3.如何计算均方位移MSD
- 4.如何计算径向分布函数g(r )以及fix ave/time命令
- 5.在lammps流体模拟中如何计算温度
- 6.如何计算单个原子的体积、位移和总能量
- 7.如何计算原子组的相互作用
- 8.如何绘制温度云图、速度云图
1.列举常用的lammps自带函数
1.xcm() 计算原子组中心坐标
xcm(group_ID,x|y|z)
#返回tool原子组的重心x坐标，存入dx变量中
variable dx equal xcm(tool,x) xcm()返回原子组group_ID重心的某一方向坐标，如需返回xyz三个方向，则调用3次即可。
2.fcm() 计算原子组受力
fcm(group_ID,x|y|z)
#返回tool原子组在y方向受力，存入fyy中
variable fyy equal fcm(tool,y) fcm()返回原子组group_ID在某一方向的受力。
3.bound() 返回原子组的边界
bound(group_ID,xmin|xmax|ymin|ymax|zmin|zmax)
#返回tool原子组在x方向最小坐标，存入变量lox中
variable lox equal bound(tool,xmin) bound()返回一个原子组的边界，可以通过设置不同的参数来求在xyz三个方向上的最大坐标和最小坐标。
4.random() ：随机数函数
random(lo,hi,seed) random()生成一个介于(lo,hi)之间的随机数，seed为随机数种子。
5.vdisplace() 位置更新函数
vdisplace(value0,velocity)
# 返回值：
# value=value0+velocity*(timestep-startstep)*dt vdisplace()根据设定的初始位置和移动速度，返回某时刻新位置。value0为初始位置，velocity为速度。
在纳米压痕模拟中，纳米压球的下压可以设置如下：
variable y equal vdisplace(50,-0.5)
fix 2 all indent 10 sphere 0 v_y 0 20 units box
# indent 在模拟框内插入一个压头，力常数为10，压头为球状，球半径20，球心(0,y,0)2.如何进行应力-应变的计算及绘制应力云图
1.应变的计算
variable tmp equal "lz" #tmp变量等于z方向的长度
variable L0 equal ${tmp} #z方向初值储存在L0中
variable strain equal "(lz-v_L0)/v_L0" #计算应变(z方向)2.应力的计算
# 在metal单位下，除以10000转换为GPa
variable stressx equal "-pxx/10000" # x方向应力
variable stressy equal "-pyy/10000" # y方向
variable stressz equal "-pzz/10000" # z方向
# 应力应变保存为文件
fix def3 all print 100 "${strain} ${stressz}" file stress_strain.dat screen no
# 每100步保存应力应变，不输出屏幕# 计算某个原子组的应力
compute s mobile stress/atom NULL #计算mobile原子组的单原子的应力
compute 2 mobile reduce sum c_s[1] c_s[2] c_s[3] #compute reduce sum命令实现单原子量的求和计算
variable stressx equal c_2[1]/(v_Vol) #compute stress/atom计算的量带有体积项，需要除以体积
variable stressx equal c_2[2]/(v_Vol)
variable stressx equal c_2[3]/(v_Vol)3.应力云图绘制
compute 1 all stress/atom NULL #计算所有原子的单原子受力
compute v all voronoi/atom #计算单原子体积
variable stressx atom c_1[1]/c_v[1]/10000 #单原子x方向应力
variable stressy atom c_1[2]/c_v[1]/10000 #y方向
variable stressz atom c_1[3]/c_v[1]/10000 #z方向
dump 1 all custom 100 all.xyz id type x y z v_stressx v_stressy v_stressz #每100步保存所有原子的坐标，xyz方向应力计算后的原子应力通过dump保存到轨迹文件中，导入ovito即可绘制应力云图。
3.如何计算均方位移MSD
在lammps扩散模拟中，大部分需要计算均方位移MSD，对MSD曲线拟合斜率可得扩算系数。
compute 1 all msd com yes
#计算所有原子的均方位移，com yes计算每个原子位移之前，减去原子组质心的漂移影响
variable msdx equal c_1[1] # x方向msd
variable msdy equal c_1[2] # y方向msd
variable msdz equal c_1[3] # z方向msd
variable msd equal c_1[4] #平均msd
variable istep equal step #时间步长
fix 2 all print 1 "${istep} ${msdx} ${msdy} ${msdz} ${msd}" file msd.dat screen no
# 每1步保存所有原子的均方位移数据，不输出屏幕4.如何计算径向分布函数g(r )以及fix ave/time命令
lammps计算径向分布函数g(r )，径向分布函数(Radial distribution function)是指给定某个粒子的坐标，其他粒子在空间的分布几率。
1.compute rdf：
compute ID group-ID rdf Nbin itype1 jtype1 itype1 jtype2 ... keyword/value ...
# Nbin 分片数目，一般选择100-500之间
# itype1 表示中心原子 jtype1 表示分布原子
# itype和jtype需成对设置，表示计算第i种原子周围出现第j种原子的概率
# 如果不设置itype和jtype，则表示全部原子周围出现其他原子的概率
compute myRDF all rdf 100
fix 1 all ave/time 100 1 100 c_myRDF[*] file tmp.rdf mode vector
# 100 1 100 每100步保存一次数据，“100 1”表示以100步为基准，每隔100步采1个样，共采1个样以c_myRDF[1]为横坐标，c_myRDF[2]为纵坐标，即可绘制g(r )曲线。
2.平均值输出命令fix ave/time：
# fix ave/time 按照设定的步数对输出变量求平均值，并将平均值保存为文件
fix ID group-ID ave/time Nevery Nrepeat Nfreq value1 value2 ... keyword args ...以输出摩擦过程中的摩擦球受力为例：
compute reduce sum命令实现单原子量的求和计算
compute fx0 ball reduce sum fx #原子在x方向受力之和
compute fy0 ball reduce sum fy #y方向
compute fz0 ball reduce sum fz #z方向
variable step0 equal step #时间步长
fix 1 all ave/time 10 5 100 v_step0 ,c_fx0 c_fy0 c_fz0 file friction.txt
# 10 5 100 每100步保存一次平均值数据，“10 5”表示以100步为基准向前每隔10步采1个样，共采5个样，即计算60 70 80 90 100这5步的平均值5.在lammps流体模拟中如何计算温度
在分子动力学模拟中，体系的温度与原子的速度有关。温度可由公式EK=(dim/2)kNT求得，EK为原子组总动能(1/2)mv2，dim为模拟纬度，N为原子总数，k为玻尔兹曼常数，T为温度。
但是在lammps流体模拟中，流体原子按照一定速度流动，所以在流体的温度计算中一般去掉流体的流动速度。
(1)compute temp/partial
compute ID group-ID temp/partial xflag yflag zflag
# ID为compute命令的ID，group-ID为需计算温度的原子组ID
# xflag yflag zflag为是否计算该方向温度，1表示包含，0表示不包含
compute myTemp flow temp/partial 1 0 1
#计算flow原子组温度时，不计算y方向速度，储存在计算myTemp中(2)compute temp/com
先计算原子组的速度，然后在计算温度时，减去质心速度。
compute myTemp mobile temp/com
# 计算mobile原子组在减去质心速度后的温度，储存在计算myTemp中
fix 1 flow nvt temp 300 300 0.1
fix_modify 1 temp mytemp
# 在nvt系综中使用新定义的计算温度myTemp6.如何计算单个原子的体积、位移和总能量
1.compute voronoi/atom计算单个原子体积
根据voronoi算法，将单个原子所占据的空间划分为一个多边形，称为泰森多边形，多边形的体积即该原子的体积。
compute 1 all voronoi/atom
# 输出两个计算结果，第一个为单原子的体积，用c_1[1]，第二个为相邻原子数目，即该原子多面体的面数应用案例：
# 计算并输出单原子应力
compute 1 all stress/atom NULL #计算单原子的应力张量，输出单位：压力*体积，NULL动能不包含在应力中
compute 2 all voronoi/atom #计算单原子的体积
variable stressx atom c_1[1]/c_2[1] # X方向的应力
dump 1 all custom 1000 dump.xyz id type x y z v_stressx #每1000步保存原子x方向应力2.compute displace/atom计算原子位移
原子位移可以用来反映材料的变形程度，单原子位移可以在dump命令中输出。
compute 1 all displace/atom
# 计算结果为4个值，分别是原子在xyz三个方向的位移和平均位移应用案例：
# 输出所有原子的平均位移
compute 1 all displace/atom #计算单原子位移
dump 1 all custom 100 cu.xyz id type x y z c_1[4] #每100步保存所有原子平均位移
# 仅输出位移大于0.6的原子
variable Dhop equal 0.6
variable check atom "c_dsp[4]>v_Dhop" #原子位移大于0.6返回1
compute dsp all diaplace/atom refresh check #refresh 与dump_modify refresh生成增量转储文件
dump 1 all custom 100 tmp.dump id type x y z #每100步保存原子坐标
dump_modify 1 append yes thresh c_dsp[4] > ${Dhop} refresh c_dsp delay 100
#append追写文件，thresh 只有位移大于0.6才会输出，refresh在每次转储结束后更新，delay延迟特定时间后输出3.计算单原子总能量
原子总能量=动能+势能
compute Ke_atom all ke/atom #计算单原子动能
compute Pe_atom all pe/atom #计算单原子势能
variable E_total atom c_Ke_atom+c_Pe_atom #atom原子矢量
dump 1 all custom 100 npt.xyz id type x y z v_E_total #每100步保存原子坐标和总能量将原子能量值保存在轨迹文件中，可以用ovito显示总能量云图。
7.如何计算原子组的相互作用
在lammps模拟中，常常需要计算原子组之间的相互作用，如摩擦球与基体之间，流体与壁面之间。而compute group/group命令可以计算原子组之间的作用能和作用力。
compute ID group-ID1 group/group group-ID2
# 返回四个结果，c_1为原子组之间的作用能，c_1[1] c_1[2] c_1[3]为原子组在x y z方向上的作用力8.如何绘制温度云图、速度云图
1.温度云图
compute KE all ke/atom #计算所有原子的单原子动能
variable KB equal 8.625e-5
variable TEMP atom c_KE/1.5/${KB} #将动能准换为温度
dump 1 all custom 1000 Cu.xyz id type x y z v_TEMP #每1000步保存原子轨迹和温度通过dump命令将原子温度保存到文件中，导入ovito绘制温度云图。
2.速度云图
虽然可以根据单原子的参数绘制云图，但是由于原子之间相对离散，绘制得到的图像过渡不平滑。
根据有限元的思想，将整个体系划分为若干个小单元，一个单元可包含多个原子，计算原子的平均值，以此绘图，得到比较平滑的图像。
compute chunk/atom将整个体系划分为若干个单元(chunk)，使用fix ave/chunk对每个单元内的参数进行平均计算并保存。
compute 1 all chunk/atom bin/2d x lower 2 y lower 2 units box #bin/2d从两个维度划分块，lower从最低点划分，分块厚度为2，使用盒子单位
fix 1 all ave/chunk 1 200 200 1 vx file velocity.txt # 1 200 200每200步保存一次数据，“”“1 200”每1步采一个样，共采集100个样，将vx x方向速度保存为文件模拟完成后，使用origin绘制等高线即可得到速度云图。
