在lammps模拟过程中的常用势函数设置

文章目录

1 lj/cut
- 1.1 lj/cut在in文件中使用方法
- 1.2 lj/cut在data文件中使用方法
- 1.3 lj/cut参数查询方法
- 1.4 lj/cut参数单位转换方法
- 1.5 lj/cut不同原子之间的参数
- 1.6 lj/cut参数设置案例
2 lj/cut/coul
- 2.1 lj/cut/coul/cut：短程作用
- 2.2 lj/cut/coul/long：长程作用
3 morse
4 tersoff
- 4.1 基本语法
- 4.2 tersoff势设置案例
5 EAM
- 5.1 eam
- 5.2 eam/alloy
- 5.3 eam/fs
6.MEAM/C
- 6.1 library.meam文件格式
- 6.2 meam文件格式
- 6.3 meam/c设置方法
7 CVFF
- 7.1 建模并转换为data文件
- 7.2 in文件的写法
8 PCFF
- 8.1 建模并转换为data文件
- 8.2 in文件的写法
9 reaxff
- 9.1 reaxff设置方法
- 9.2 reaxff文件下载
10 TIP4P水分子势
- 10.1 TIP4P介绍
- 10.2 TIP4P设置案例
- 10.3 代码注释
11 混合力场设置
- 11.1 eam/lj混合势
- 11.2 tersoff/eam/lj混合势
- 11.3 常见的错误
- 11.4 hybrid/overlay设置叠加力场
- 8.5 hybrid/scaled设置不同权重的多种力场
12 高熵合金势函数设置方法
- 12.1 下载专用势函数
- 12.2 使用混合势
- 12.3 拟合势参数
- 12.4 官方自带EAM势拟合方法
13 OPLS
- 13.1 OPLS参数设置
14 born
- 14.1 born势设置
15 Buckingham势
16 comb3

在lammps模拟中，势函数的使用必不可少。如何选择合适的势函数？如何设置势函数？这对初学者来说较为困难，本文介绍在lammps模拟中常用的势函数如何使用，希望能够帮助初学者加深对势函数的了解。

1 lj/cut

lj/cut力场公式：
E=4ϵ[(σr)12−(σr)6]r<rcE=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]\quad\quad r<r_{c}E=4ϵ[(rσ)12−(rσ)6]r<rc
在lammps使用lj/cut只需要设置3个参数：ϵ\epsilonϵ、σ\sigmaσ、rcr_{c}rc
lj/cut有两种写法，第一种是写到in文件中，另一种是写到data文件中，两种写法稍有不同。

1.1 lj/cut在in文件中使用方法

在in文件中使用lj/cut，需要明确的指定相互作用的原子类型，如：

pair_style   lj/cut 10  #截断半径常用10-15，最常用10，12，全局
pair_coeff   * * 0.02 3.12
pair_coeff   1 1 0.03 3.22 8.5  #局部截断半径

上面两句指定全部原子之间的lj/cut力场，截断半径为10，ϵ\epsilonϵ=0.02、σ\sigmaσ=3.12，1和1原子截断半径为8.5，ϵ\epsilonϵ=0.03、σ\sigmaσ=3.22。

1.2 lj/cut在data文件中使用方法

在data文件中，只能指定同种原子之间的lj/cut参数，每一行设置一种原子，不能设置不同原子之间的力场参数。
例如：

Pair Coeffs # lj/cut
1    0.038     2.44
2    0.019     3.01

上述命令分别设置了原子1之间的受力、原子2之间的受力，没有设置原子和原子之间的受力。

1.3 lj/cut参数查询方法

大多数原子lj/cut参数都可以从文献中查到，推荐一个LJ参数查询网站：https://lammpstube.com/mdpotentials/。

1.4 lj/cut参数单位转换方法

如果是在metal单位下进行Si的lammps模拟，需要把kcal/mol单位转换为eV，推荐另一个网站：https://www.colby.edu/chemistry/PChem/Hartree.html，进入网站后在相应的位置输入能量值，可以自动转换为其他单位。

1.5 lj/cut不同原子之间的参数

对于1和2之间的参数，lammps提供了三种拟合公式：
1.mix geometric (默认公式)
ϵij=ϵiϵjσij=σiσj\epsilon_{ij}=\sqrt{\epsilon_{i}\epsilon_{j}}\\ \sigma_{ij}=\sqrt{\sigma_{i}\sigma_{j}}ϵij=ϵiϵjσij=σiσj
2.mix arithmetic
ϵij=ϵiϵjσij=12(σi+σj)\epsilon_{ij}=\sqrt{\epsilon_{i}\epsilon_{j}}\\ \sigma_{ij}=\frac{1}{2}\left(\sigma_{i}+\sigma_{j}\right)ϵij=ϵiϵjσij=21(σi+σj)
3.mix sixthpower
ϵij=2ϵiϵjσi3σj3σi6+σj6σij=(12(σi6+σj6))16\epsilon_{ij}=\frac{2\sqrt{\epsilon_{i}\epsilon_{j}}\sigma_{i}^{3}\sigma_{j}^{3}}{\sigma_{i}^{6}+\sigma_{j}^{6}}\\\sigma_{ij}=\left(\frac{1}{2}\left(\sigma_{i}^{6}+\sigma_{j}^{6}\right)\right)^{\frac{1}{6}}ϵij=σi6+σj62ϵiϵjσi3σj3σij=(21(σi6+σj6))61
这个拟合过程由lammps自动完成，不需要人为干预，pair_modify命令可以选择拟合方式。

1.6 lj/cut参数设置案例

例如：查文献得到Cr-Cr和Fe-Fe的参数为：（ϵ\epsilonϵ、σ\sigmaσ）

Cr-Cr: 0.502 2.336
Fe-Fe: 0.527 2.321

Cr-Fe的lj/cut参数就可以使用默认的公式进行拟合，经过拟合后的Cr-Fe参数为：

Cr-Fe: 0.514 2.3285

lj/cut力场参数比较简单，模拟结果虽然不如专用的力场精确，但参数获取比较方便，但找不到专用力场参数时，使用lj/cut也是一种比较好的选择。

2 lj/cut/coul

普通的lj势只考虑了原子之间的吸引力与排斥力，没有考虑原子之间的电荷作用。lj/cut/coul力场在普通lj势的基础上增加库仑力，其中常用的为lj/cut/coul/cut和lj/cut/cout/long两种力场。

2.1 lj/cut/coul/cut：短程作用

lj/cut/coul/cut力场公式：
E=4ϵ[(σr)12−(σr)6]r<rcE=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]\quad\quad r<r_{c}E=4ϵ[(rσ)12−(rσ)6]r<rc
E=Cqiqjϵrr<rcE=\frac{Cq_{i}q_{j}}{\epsilon r}\quad\quad r<r_{c}E=ϵrCqiqjr<rc
lj/cut/coul/cut的截断半径有两个：lj势的截断半径与库仑力的截断半径，若只写一个截断半径，默认cutoff2 = cutoff。

# 只设置一个截断半径
pair_style   lj/cut/coul/cut 10.0   #全局截断半径
pair_coeff   * * 100.0 3.0# 设置两个截断半径
pair_style   lj/cut/coul/cut 10.0 8.0
pair_coeff   1 1 100.0 3.5 9.0 9.0      #原子1-1之间使用截断半径(局部)9.0 9.0

2.2 lj/cut/coul/long：长程作用

lj/cut/coul/long和lj/cut/coul/cut相比，不仅要计算邻近原子之间的库仑力，还要计算邻近原子之外的原子对中心原子的库仑力。

pair_style   lj/cut/coul/long 10.0      #全局截断半径
pair_coeff   * * 100.0 3.0
kspace_style   pppm 1.0e-5
#计算长程库伦力时，需要配合kspace_style pppm或kspace_style ewald使用

长程库仑力的计算需要特定算法，kspace_style命令。但是这个算法经常提示错误，最常见的错误，“Out of range atoms - cannot compute PPPM”。出现原因分为两种情况：
1.原子移动速度相对较低，在一个时间步内移动距离超过1/2的neighbor bin，同时中心原子的邻居列表更新不及时。解决方法：提高邻居列表的更新频率，减小时间步长。
2.原子移动速度相对较高，当原子重叠度较高或者力场参数设置不对，原子之间的受力过大，导致原子移动速度过快，在一个时间步长内跑出盒子。解决方法：重新建模，若用MS建模，可以在MS跑平衡，获得一个近似平衡稳定的体系。

3 morse

morse力场与lj力场相似，用来描述没有键相互连接的原子之间的作用力。
morse力场公式：E=D0[e−2α(r−r0)−2e−α(r−r0)]r<rcE=D_{0}\left[e^{-2\alpha\left(r-r_{0}\right)}-2e^{-\alpha\left(r-r_{0}\right)}\right]\quad\quad\quad\quad r<r_{c}E=D0[e−2α(r−r0)−2e−α(r−r0)]r<rc
在使用morse力场时需要设置4个参数，分别为：D0D_{0}D0、α\alphaα、r0r_{0}r0、rcr_{c}rc。
morse力场使用方法：

pair_style   morse cutoff
pair_coeff   * * D0 a r0 [cutoff]

例如:

pair_style   morse 2.5          #全局截断半径
pair_coeff   * * 100.0 2.0 1.5
pair_coeff   1 1 100.0 2.0 1.5 3.0
#1-1原子之间使用局部截断半径3.0

morse势参数查询网址：http://apot.mgedata.cn/InteratomicPotentials/index。

4 tersoff

tersoff势是一种非键接(non-bond)势，在SiC、GaAs、GaN等体系中用的较多。
tersoff势参数保存在一个文本文件中，通常以“.tersoff”为后缀名，因此，在lammps中，不需要设置tersoff势的具体参数，仅需指定对应的原子类型即可。
下面介绍不同情况下，tersoff势设置方法。

4.1 基本语法

假设体系中只包含Si、C两种原子，对应的原子类型分别为：Si(type 1)、Si(type 2)、C(type 3)。
tersoff势写法：

pair_style    tersoff
pair_coeff    * * SiC.tersoff Si Si C

第一行pair_style指定势函数类型为tersoff。
第二行pair_coeff映射原子类型，pair_coeff命令后面必须为“∗*∗”，不能写具体的原子类型(如1 1)。
“∗*∗”后面为tersoff势文件名称，最后一部分为原子列表。
“Si Si C”表示前两种原子类型为Si，第3种原子为C，lammps在积分运算时会自动根据这个映射关系到“SiC.tersoff”文件找出原子间的参数。
这部分顺序必须与data文件或者体系模型中的原子类型相对应，否则会计算出错。
如体系中C原子类型为2，Si原子类型分别为1和3，则代码要改为：

pair_coeff    * * SiC.tersoff Si C Si

4.2 tersoff势设置案例

SiC切削Si基体，原子类型：

1 28.0855 #Si
2 28.0855 #Si
3 28.0855 #Si
4 28.0855 #Si
5 12.0107 #C

势函数写法：

pair_style   tersoff
pair_coeff   * * SiC.tersoff Si Si Si Si C

5 EAM

金属原子之间没有键连接，因此，在lammps模拟中，金属体系的势函数类型为pair_style，而不是bond_style。
模拟金属体系时，可以用lj势描述金属原子之间的受力，不过更精确的是嵌入原子势(EAM)，eam势函数公式为：
Ei=Fα(∑j≠iρβ(rij))+12∑j≠iϕαβ(rij)E_{i}=F_{\alpha}\left(\sum_{j\ne i}\rho_{\beta}\left(r_{ij}\right)\right)+\frac{1}{2}\sum_{j\ne i}\phi_{\alpha\beta}\left(r_{ij}\right)Ei=Fα⎝⎛j=i∑ρβ(rij)⎠⎞+21j=i∑ϕαβ(rij)
eam势由两部分组成，在原子对势(pair)的基础上添加了电子云密度相关项，比单纯的对势精确度更高。
eam势函数写在一个以“.eam”为后缀的文件中，lammps自带的势函数包含一部分eam势文件，也可以到网站下载eam文件，下载eam文件后，保存到in文件所在的文件夹。
常用的势函数下载网站有：
http://www.ctcms.nist.gov/potentials
https://openkim.org

5.1 eam

当模拟体系中只包含一种金属原子时，势函数的设置比较简单。
单原子eam势文件名后不需要进行原子类型映射，不用写原子类型列表。

pair_style   eam
pair_coeff   1 1 Cu_u6.eam

5.2 eam/alloy

对于合金体系，对于的eam势函数为eam/alloy或者eam/fs，写法稍有不同。
pair_style指明eam合金势类型，pair_coeff映射原子类型。合金体系必须在pair_coeff语句中势文件名后面把所有的原子类型全部列出，顺序和in文件中原子类型要保持一致。
例如，体系中包含Ni、Al两种原子，Ni的原子类型分别为1、3、4，Al原子类型为2，eam/alloy类型写法：

pair_style   eam/alloy
pair_coeff   * * NiAlH.eam.alloy Ni Al Ni Ni

5.3 eam/fs

eam/fs类型，Ni原子类型为1、2、3，Al原子类型为4：

pair_style   eam/fs
pair_coeff   * * NiAlH_jea.eam.fs Ni Ni Ni Al

6.MEAM/C

合金体系的势函数除了eam势，还有meam势。在新版本的lammps中，meam势类型已经改为meam/c。
和普通的势文件不同，meam/c势有两个函数势文件：library.meam和**.meam**，表示不同的势函数名称。

6.1 library.meam文件格式

library.meam类似于参数库，存储了多种元素的meam参数，每一行元素占3行，共19个通用参数。
例如，下面一段代码是ZrCu合金的library.meam文件：

#elt lat z   iel atwt
#       alpha       b0  b1  b2  b3      alat    esub    asub
#       t0  t1  t2  t3  rhozero ibar'Zr'  'hcp'  12 1   91.22004.4501908328 2.450   1.000   3.000   2.000   3.2000000000    6.360   0.6801.00   6.300   -3.300  -10.00  1.000   3'Cu' 'fcc'  12 1   63.54605.1548300830 3.830   2.200   6.000   2.200   3.6133156519    3.540   0.9401.00   2.720   3.040   1.950   1.000   3

library.meam前三行是注释部分，说明了各行参数的定义，后面分别是Zr和Cu对应的参数。

6.2 meam文件格式

第二个meam文件存储合金元素专用的参数，描述合金原子之间的相互作用。
如ZrCu.meam参数如下：

rc = 5.0
delr = 0.1
augt1 = 0
erose_form = 2
ialloy = 2zbl(1,1) = 0
nn2(1,1) = 1
rho0(1) = 1.000
Ec(1,1) = 6.360
re(1,1) = 3.2000
alpha(1,1) = 4.45019083
repuls(1,1) = 0.00
attrac(1,1) = 0.00
Cmin(1,1,1) = 1.00
Cmax(1,1,1) = 1.44...

可以看出，在ZrCu.meam中没有出现Zr和Cu元素名称，而是用序号1、2表示。
在这里，1和2为索引号，并不是in文件中原子类型编号。

6.3 meam/c设置方法

把library.meam和ZrCu.meam文件保存到in文件同一文件夹。
假设in文件中只有两种原子，原子类型1为Zr，2为Cu，势函数设置为：

pair_style   meam/c
pair_coeff   * * library.meam Zr Cu ZrCu.meam Zr Cu

假设原子类型1为Cu，2为Zr，势函数设置为：

pair_style   meam/c
pair_coeff   * * library.meam Zr Cu ZrCu.meam Cu Zr

假设有四种原子类型，1和2为Cu，3和4为Zr，势函数设置为：

pair_style   meam/c
pair_coeff   * * library.meam Zr Cu ZrCu.meam Cu Cu Zr Zr

7 CVFF

cvff势是由pair、bond、angle、dihedral、improper等势组成，在ms中直接设置cvff势即可，但是在lammps中，需要分别设置以上各部分势。
在lammps中，cvff没有势文件，只要设置对应的势类型和参数即可。一般情况下，cvff势不需要自己找参数。
最简单的方式是在ms中建立模型，设置cvff势后，导出为car文件。 使用免费的msi2lmp转换工具，把car文件转换为lammps可识别的data文件。转换完成后，data文件内自带cvff势参数。
下面以Cu和聚乙烯复合物为例，介绍cvff势的具体使用方法。

7.1 建模并转换为data文件

ms建模完成后导出文件为layer.car，使用msi2lmp转换为data文件：

msi2lmp layer -class I -frc cvff -i >layer.data

转换之后得到 layer.data。
data文件中的pair_coeff只需写出同种原子之间的势参数，不同原子之间的势参数自动计算，之后的bond、angle同理。

7.2 in文件的写法

既然data文件已经自动设置cvff势，在in文件中只需写明势的类型即可，势的类型就是data中各种势“#”后面的名称。
cvff势默认的设置语句要放到read_data命令的前面。 cvff默认的pair势有长程库仑力，因此需要设置kspace_style。
in文件cvff势具体设置为：

pair_style      lj/cut/coul/long   10 12  #截断半径10 12
bond_style      harmonic     #键相互作用
angle_style     harmonic     #角相互作用
dihedral_style  harmonic   #二面角相互作用
kspace_style    pppm  1e-4   #库仑势算法，(pppm)适用于周期性边界条件ppp
read_data       layer.data

#非周期性边界条件
pair_style      lj/cut/coul/cut   10 12  #截断半径10 12
bond_style      harmonic     #键相互作用
angle_style     harmonic     #角相互作用
dihedral_style  harmonic   #二面角相互作用
read_data       layer.data

读取文件后，如果不需要替换参数，直接就可以进行弛豫计算。

8 PCFF

cvff、pcff是ms文件转换为lammps data文件最常用的两种势。相比于cvff势，pcff势参数更多，但是在设置方式上和cvff势过程是一样的。
以沥青材料为例，介绍pcff势设置方法。

8.1 建模并转换为data文件

在ms中使用AC模块建立沥青模型，使用forcite模块设置pcff力场，导出为asphalt.car。
使用msi2lmp转换为data文件：

msi2lmp asphalt -class II -frc pcff -i >asphalt.data

转换之后得到 asphalt.data。
pcff在lammps中对应的势函数多为class2类型。
pcff势对应的lj势为9-6势，如果需要修正参数的话，需要注意参数的转换。

8.2 in文件的写法

data文件中势函数部分已经把全部的参数列出，在in文件中只需写明势的类型即可，势的类型是data中各种势"#"后面的名称。
pcff势的设置语句要放到read_data命令的前面。 pcff默认的pair势有长程库仑力，因此需要设置kapace_style。
in文件pcff势具体设置为：

pair_style      lj/class2/coul/long 10 12 #截断半径10-12
bond_style      class2  #键相互作用
angle_style     class2  #角相互作用
dihedral_style  class2  #二面角相互作用
improper_style  class2   #非正常二面角相互作用
kspace_style    pppm 1e-4
read_data       asphalt.data

读取文件之后，如果不需要替换参数，直接就可以进行弛豫计算。

9 reaxff

在大多数的lammps模拟中，两个原子之间只要有键(bond)连接，这个键就默认永久存在。在模拟过程中，这个键不允许断开，同时也不会有新键生成。
但当键连接的两个原子距离过大时，超过势函数允许范围时，就会产生"bond missing"错误，迫使模拟中止。但reaxff反应势允许键的生成与断裂，可以用于模拟聚合物的交联或者裂解。

9.1 reaxff设置方法

1.设置原子类型为charge

units        real                  #单位选取real类型
atom_style   charge

2.设置pair_style

pair_style   reaxff lmp_control

reax/c为反应势类型说明，lmp_control为输出控制文件，用于控制输出轨迹文件trj，若不需要，可以设置为NULL。
3.设置pair_coeff

pair_coeff  * * ffield.reax.cho H C O
# 力场文件后列出原子类型即可
# 势函数文件 ffield.reax.cho 存放在potential文件夹下

4.设置电荷平衡

fix   2 all qeq/reaxff 1 0.0 10.0 1e-6 param.qeq
# 每隔1步进行一次电荷平衡qeq，锥度半径0.0-10.0，平衡电荷精度为1e-6
# param.qeq存放所需原子的参数，自左向右分别为原子类型、电负性(eV)、自库伦势、轨道指数，若该数值存在于势文件中，可将param.qeq替换为reaxff# 如果不需要进行电荷平衡，可以在pair_style中设置
pair_style reaxff controlfile checkqeq no
# checkqeq no 不对电荷平衡(qeq)进行fix检查

9.2 reaxff文件下载

scm网站是lammps官方推荐的一个反应势下载网站，https://www.scm.com/doc/ReaxFF/Included_Forcefields.html。
打开之后，点击参考文献，跳转到论文页面，在该界面找到"Support Information"，下载后的势文件多为pdf文件，之后可以复制文本转换为txt文件。

10 TIP4P水分子势

在lammps模拟中，水的类型比较多，有TIP3P、SPC、TIP4P等类型。不同的类型对应的水分子结构和力场参数均有所差别，在设置时需要注意区分。

10.1 TIP4P介绍

TIP4P水在正常水分子的基础上增加了一个虚原子，具体可查阅lammps官方手册(https://docs.lammps.org/Howto_tip4p.html)，虚原子不是实体原子，因此在水分子建模中不需要建出虚原子，仅在pair势中设置即可。
这两种势对应的公式一样的，由两部分组成：LJ势和库仑势。
LJ势对应公式：
E=4ϵ[(σr)12−(σr)6]r<rcE=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]\quad\quad r<r_{c}E=4ϵ[(rσ)12−(rσ)6]r<rc
库仑势对应公式：
E=Cqiqjϵrr<rcE=\frac{Cq_{i}q_{j}}{\epsilon r}\quad\quad r<r_{c}E=ϵrCqiqjr<rc

10.2 TIP4P设置案例

units               real                #real单位体系
atom_style          full                #full原子类型
boundary            p p p             #周期性边界条件pair_style          lj/cut/tip4p/long   1 2 1 1 0.1546 12.0
#pair_style 参数设置：氧原子类型，氢原子类型，键类型，角类型，虚原子距离，截断半径
bond_style          harmonic     #键类型
angle_style         harmonic     #角类型
kspace_style        pppm/tip4p 1.0e-4    #长程库仑力算法
pair_modify         shift yes mix arithmetic
read_data           H20.data#设置力场参数
pair_coeff          1 1 0.1852 3.1589   #O-O原子
pair_coeff          2 2 0 0  #H-H原子
bond_coeff          1 1000 0.9572  #H-O键
angle_coeff         1 1000 104.52  #H-O-H角
fix                 1 all shake 0.0001 50 0 b 1 a 1
#通过SHAKE算法对键角进行约束，迭代精度差0.0001，迭代步数50，不输出运行结果，对类型1的键和类型2的键约束

10.3 代码注释

第一行中pair_style 明确类型lj/cut/tip4p/long。
1 2 1 1表示：O原子类型、H原子类型、O-H键类型、H-O-H角类型，这些与data文件中类型相对应。
0.15表示虚原子与水分子中O原子的距离
12.0表示LJ势中的截断半径(cutoff1)
10.0表示库仑力中的截断半径(cutoff2)，如果该参数省略，则默认等于cutoff1。
TIP4P势在设置势需要注意以下几点
(1)原子排列顺序
对应系统中的每个TIP4P水分子，O和2个H原子的原子ID必须是连续的，并且O原子在前，顺序为O H H。
例如，如果TIP4P水分子中一个O原子的原子ID为500，则另外两个H原子的ID必须为501和502。
(2)设置kspace_style
```
kspace_style    pppm/tip4p  1.0e-4
```
(3)固定键角
在模拟过程中需要固定水分子的键和角，可使用以下命令：
```
fix      1 all shake 0.0001 20 10 b 1 a 1
```
在lammps模拟中，水分子力场设置相对复杂，在使用过程中一定要注意参数设置是否正确，否则可能影响计算结果。

11 混合力场设置

在分子动力学模拟中，势函数描述了原子之间的作用力，不同的原子需要设置不同的势函数类型。
对于复合材料来说，因为具有多种原子类型，如果没有专用的势函数，还可以设置混合势函数来解决问题。在lammps中，使用hybrid可以设置两种以上的势函数。

11.1 eam/lj混合势

以Al/Graphene 复合物烧结过程为例，详解混合势函数的设置方法。模型为8个Al原子组成的纳米球，中间插入一层石墨烯，温度由300K升温到1000K结束。
在本例中，需要设置三种势函数

pair_style      hybrid  eam/fs airebo 3.0 lj/cut 10
pair_coeff      * * eam/fs Al_nm.eam.fs Al NULL   #Al原子之间，* *须注明所有的原子列表，NULL占位符
pair_coeff      * * airebo CH.aireo NULL C        #C原子之间
pair_coeff      1 2 lj/cut 0.023 2.852            #C和AL原子之间

11.2 tersoff/eam/lj混合势

如果体系中还包含使用其他势的原子，则需要使用混合势写法。
以SiC和Cu体系为例，体系中原子类型分别为：Si(tyoe 1)、C(type 2)、Cu(type 3)，对应的写法：

pair_style        hybrid tersoff eam lj/cut 10.0
pair_coeff        * * tersoff SiC.tersoff Si C NULL
pair_coeff        3 3 eam Cu.eam
pair_coeff        1 3 lj/cut 0.0034 3.12
pair_coeff        2 3 lj/cut 0.053 2.98

pair_style命令需要指出设置的三种势函数类型
第2行表示前2种原子(Si、C)使用tersoff势，第三种原子位置为NULL表示该原子(Cu)不使用SiC势。
第3行设置Cu为eam势
第4、5行设置Si-Cu、C-Cu为lj/cut势。

11.3 常见的错误

比较常见的一个错误是：All pair coeffs are not set，错误提示为“不是所有的原子都设置了势参数”，表示有的原子可能没有设置势参数。
出现这种错误主要有两种原因：
1.原子数量多，的确少写了某种原子的势参数。
2.设置了全部的原子势函数，但是设置方法出错。
下面一个简单的Cu-Al模型为例，给出这种错误的解决方法。
Cu原子类型为1，Al原子类型为2，Cu和Al均使用meam势，Cu-Al之间使用LJ势。
势函数设置如下：

pair_style    hybrid meam/c lj/cut 10
pair_coeff    * * meam/c library.Cu.meam Cu Cu.meam Cu NULL  #Cu
pair_coeff    * * meam/c library-Al.meam Al Al.meam NULL Al  #Al
pair_coeff    1 2 lj/cut 0.4 2.47  #Cu-Al

如果按照以上代码设置，就会提示“All pair coeffs are not set”。
出现这种错位，主要是因为同一in文件中使用了两个meam势，第二个Al的meam势会覆盖掉Cu的meam势，导致Cu的meam势参数丢失。
解决方法：
如果在同一in文件中使用多个同种类型的势，**为防止覆盖，需要对相同类型的势进行编号区分，**所以正确写法：

pair_style    hybrid meam/c meam/c lj/cut 10
pair_coeff    * * meam/c 1 library.Cu.meam Cu Cu.meam Cu NULL  #Cu
pair_coeff    * * meam/c 2 library-Al.meam Al Al.meam NULL Al  #Al
pair_coeff    1 2 lj/cut 0.4 2.47  #Cu-Al

11.4 hybrid/overlay设置叠加力场

hybrid为在模拟体系中使用多种力场，但对于体系中的原子来说，只设置了一种力场。而hybrid/overlay则为体系中的原子同时设置多种力场，实现力场叠加。

pair_style   hybrid/overlay lj/cut 2.5 coul/long 2.0
pair_coeff   * * lj/cut 1.0 1.0
pair_coeff   * * coul/long
# coul/long力场并不会覆盖之前的lj/cut力场
# 原子之间受力同时受到lj/cut和coul/long两种力场控制

相当于以下代码：

pair_style   lj/cut/coul/long 2.5 2.0
pair_coeff   * * 1.0 1.0

hybrid/overlay将多种力场叠加在一起，但在力的计算中，这些力场之间的权重相同。

8.5 hybrid/scaled设置不同权重的多种力场

pair_style  hybrid/scaled 0.3 tesoff 0.7 sw
pair_coeff  * * tersoff Si.tersoff
pair_coeff  * * sw.Si.sw Si
# Si受到tersoff和sw力场控制，tersoff权重0.3，sw权重0.7

12 高熵合金势函数设置方法

高熵合金包含的原子数较多，势函数的设置相对复杂。下面介绍三种高熵合金势函数设置方法。

12.1 下载专用势函数

下面网站包含了大多数原子的势函数：https://www.ctcms.nist.gov/potentials/

12.2 使用混合势

如果上面的网站找不到专用的势函数，可以下载多个势函数，使用混合命令组合到一起。
假设在FeCMnSi中加入Ti，组成一种新的合金FeCMnSiTi合金，但是并不能找到这种合金的势函数。
但是能找到FeCMnSi(FeCMnSi.eam.alloy)的势和Ti(Ti.eam.fs)的势，可以使用hybrid命令合并到一起。 Ti原子与FeCMnSi原子之间使用LJ势。

12.3 拟合势参数

lammps官方安装文件中自带一个拟合程序，该程序可自动生成所需要的合金势，代码收录在官方Ubuntu版lammps源代码tools/eam_database目录内。
可以拟合的合金原子包括：Cu，Ag，Au，Ni，Pd，Pt，Al，Pb，Fe，Mo，Ta，W，Mg，Co，Ti，Zr。

12.4 官方自带EAM势拟合方法

1.编写creat.f文件
在linux系统中进入tool/eam_database目录，运行以下代码进行编译：

gfortran create.f
# 若没有权限可以
# sudo gfortran create.f

运行结束后，生成a.out程序。
2.编写EAM.input文件

&funccard
atomtype='Al'
&end
&funccard
atomtype=...
&end
&funccard
...
&end
# 依据给出的例子添加或删除合金元素

3.运行a.out程序

a.out < EAM.input

生成对应的合金势文件，文件名：合金元素.set。

13 OPLS

13.1 OPLS参数设置

OPLS力场不需要力场文件，但需要设置力场参数。OPLS包含了键(bond)、角(angle)、非正常二面角(dihedral或torison)以及非键接势(non-bonded)。在lammps设置中，这些势需要单独设置。
OPLS对应公式：
Bondstretching:Ebond=∑bondsKr(r−req)2Bond stretching:\quad \quad E_{bond}=\sum_{bonds}K_{r}\left(r-r_{eq}\right)^{2}Bondstretching:Ebond=bonds∑Kr(r−req)2
Anglebending:Eangle=∑anglesKθ(θ−θeq)2Angle bending:\quad \quad E_{angle}=\sum_{angles}K_{\theta}\left(\theta-\theta_{eq}\right)^{2}Anglebending:Eangle=angles∑Kθ(θ−θeq)2
Torsion:E(ϕ)=V12[1+cos(ϕ+f1)]+V22[1−cos(2ϕ+f2)]+V32[1+cos(3ϕ+f3)]Torsion:\quad E\left(\phi \right) =\frac{V_{1}}{2}\left[1+cos\left(\phi+f1\right)\right]+\frac{V_{2}}{2}\left[1-cos\left(2\phi+f2\right)\right]+\frac{V_{3}}{2}\left[1+cos\left(3\phi+f3\right)\right]Torsion:E(ϕ)=2V1[1+cos(ϕ+f1)]+2V2[1−cos(2ϕ+f2)]+2V3[1+cos(3ϕ+f3)]
Non−bonded:Eab=∑iona∑jonb[qiqje2rij+4εij(σij12rij12−σij6rij6)]fijNon-bonded: \quad E_{ab}=\sum^{ona}_{i}\sum^{onb}_{j}\left[\frac{q_{i}q_{j}e^{2}}{r_{ij}}+4\varepsilon_{ij}\left(\frac{\sigma^{12}_{ij}}{r^{12}_{ij}}-\frac{\sigma^{6}_{ij}}{r^{6}_{ij}}\right)\right]f_{ij}Non−bonded:Eab=i∑onaj∑onb[rijqiqje2+4εij(rij12σij12−rij6σij6)]fij
fij=0.5ifi,jare1,4;otherwise,fij=1.0\quad \quad \quad \quad f_{ij}=0.5\quad if\quad i,j\quad are \quad1,4;\quad otherwise,\quad f_{ij}=1.0fij=0.5ifi,jare1,4;otherwise,fij=1.0
(1)bond

bond_style    harmonic
bond_coeff    5 80.0 1.2
# bond势为谐振势 harmonic，或者势弹簧式

(2)angle

angle_style   harmonic
angle_coeff   1 300.0 107.0
# angle是谐振势harmonic

(3)dihedral或torsion

dihedral_style    opls                     # opls势
dihedral_coeff    1 1.740 -0.157 0.279 0.00

(4)non-bonded

pair_style       lj/cut/coul/long 10.0 8.0         #lj势+库仑势
pair_coeff       1 1 100.0 3.5 9.0 9.0
# 设置势的权重
special_bonds    lj/coul 0 0 0.5
# 1-2原子对 1-3原子对的权重为0，1-4原子对的权重为0.5

14 born

14.1 born势设置

born势对应公式：
E=Aexp(σ−rρ)−Cr6+Dr8r<rcE=Aexp\left(\frac{\sigma -r}{\rho}\right)-\frac{C}{r^{6}}+\frac{D}{r^{8}}\quad \quad r<r_{c}E=Aexp(ρσ−r)−r6C+r8Dr<rc
&enap;&enap;&enap;&enap;born势不需要势文件，仅需要设置势类型和势参数。

# 不计算库仑力
pair_style    born [cutoff]
# 计算库仑力
pair_style    born/coul/long [cutoff]

15 Buckingham势

Buckingham势基本公式：
E=Ae−rρ−Cr6r<rcE=Ae^{-\frac{r}{\rho}}-\frac{C}{r^{6}}\quad \quad r<r_{c}E=Ae−ρr−r6Cr<rc
pair_style buck命令包括多个类型，pair_style buck/coul/cut(包括长程库仑力)，pair_style buck/coul/long(包括长程库仑力)，pair_style buck/coul/msm(包括长程MSM库仑力)。
举例：

pair_style   buck/coul/long 10.0
pair_coeff   1 1 100.0 1.5 200.0
# 或者
pair_style   buck/coul/cut 10.0 8.0
pair_coeff   * * 100.0 1.5 200.0
# SiO2设置
pair_style   buck/coul/long 10
pair_coeff   1 1 0 0.1 0       #si-si
pair_coeff   2 2 32025.73059 0.3623188 4035.5875      #O-O
pair_coeff   1 2 415158.1815 0.2052124 3079.453555     #si-O
kspace_style     ewald 1e-4        #长程库仑力算法

16 comb3

lammps支持comb3势函数，并且在potentials文件夹中提供了comb3的势函数文件：ffield.comb3，但是其支持的原子类型有限，只支持O Cu N C H Ti Zn Z这几种原子。

pair_style   comb polar_off
pair_coeff   * * ffield.comb3 Cu C C

因为体系中还含有C原子，需要加载lib.comb3文件，将ffield.comb3 lib.comb3两个文件复制到模拟文件夹中。
pair_style 确定势函数类型，polar_off表示不包含原子极化。
pair_coeff确定势参数，* *表示所有原子类型，ffield.comb3为势函数文件名，Cu C C映射原子类型。