SW势能是分子动力学模型中相对较为简单的3体经典势能模型,其计算量小,能够描述结构中2体的键和3体的键角等信息。
通常是在gulp中实现sw势能模型中各参数的拟合,但是gulp中sw2和sw3的公式与lammps中公式有区别,可以分别查阅gulp和lammps的手册,很容易实现两者的相互转换,
但是本文要指出的是:将gulp势能参数直接按照上述公式对应转换为lammp势能文件时,可能导致lammps计算中能量实际上为gulp中的两倍!!,虽然表面上不影响分动的系综计算,但是在具体应用中,如EMD计算热导率性质,通过热流自关联函数的时间积分计算热导率,此时系统中的热流实际上为真实的2倍,导致最终热导率数值会比理论值大4倍!!!,如gulp中通过alamode实现RTA下热导率数值为5.72,但是用lammps的EMD得到数值为12.02(显然不正确,因为RTA近似下热导率仅仅考虑了2阶和3阶力常数,EMD分动计算包括所有非谐效应,所以EMD的热导率数值必然要小于RTA的结果),实际上lammps的计算数值是放大了4倍的结果,应为12.02/3=3.0(这个结果可以接受);以上出现4倍问题的原因是(仅针对sw模型):
gulp中的sw势能模型是通过sw2和sw3选项参数实现的,gulp作为分子静力学计算软件,其在原胞或者super命令建立的超胞内,通过寻找所有独立的2体对和3体对,然后从sw2和sw3中找到其对应的能量,最后求和得到总能量;
但是lammps作为分子动力学软件,是针对胞内每个原子,然后通过cutoff建立近邻环境的原子列表,因此能量的计算首先是针对某个原子,然后从近邻环境的列表循环中给出针对某个原子的所有2体对,3体对,其中3体对是针对近邻环境列表中第2个原子和第3个原子的2重循环,这点从lammps的sw势能文件中所有项目的的罗列可以看出,因此实际上一个3体对在上述循环中出现了2次,同样2体对的能量对2体中每个原子都计入了一遍,因此导致最终胞的总能是gulp计算的2倍!,所以在sw势能参数转换时,除了公式中相应参数一一对应转换以外,还要在lammps的sw势能文件中每个条目的总强度系数上要除2,这样分动计算中势能模型才和gulp完全一致,热导率等涉及热流计算,以及杨氏模量等力学性能都会受上述影响。
具体转换的python脚本可以和我联系索取,qq1796634232,也欢迎和我大家交流。