使用GULP(General Utility Lattice Program)拟合经典势能模型的方法

GULP(General Utility Lattice Program)

链接地址：http://gulp.curtin.edu.au/gulp/

这是一个采用经典势能模型的计算平台，可以在经典势场模型的基础上计算很多物性，包括结构优化，声子计算，热输运，缺陷以及分子动力学等相关计算。

在这些物性计算的基础上，就可以进一步开展经典势能模型的参数拟合，因为经典分子动力学计算的准确程度完全取决于原子间作用模型的准确性和参数的合理性。

GULP使用需要制备一个.gin文件，内部包含Keywords和Options，所有关键词和选项参数可参考链接 http://gulp.curtin.edu.au/gulp/help/help_52_txt.html

势能拟合关键词只需要fit，拟合变量为相应势能模型的参数，拟合目标物理量分为两类：

第1类是structure相关的derivatives,即晶格参数如cell或vector内自由度以及晶胞内原子位置的自由度，由相应flag控制是否拟合，拟合时认为给出的cell或vector或原子相应位置处derivatives为0，可理解为给出结构处各原子在平衡位置，晶胞参数在相应取值处应力为0，相应的权重利用default_weight grad 1.0设置，即默认每个flag控制的权重为1；

第2类是除了structure相关derivatives的物理量，即在每组结构下observable内设置，可以是给出结构的能量，弹性系数，声子谱等。

因此，在开始势能拟合前，先要了解gulp内相应势能模型的表达式和参数设置，进行势能拟合可在.gin文本文件内使用fit关键词，然后利用选项

default_weight

改变势能拟合权重的默认值：

angle angle observable 1.0
bond bond length observable 1.0
cell_length cell parameter weight in relax fit 1000.0
cell_angle cell parameter weight in relax fit 1000.0
coord Cartesian coordinate weight (0-2-D) 1000.0
dipole dipole moment 10.0
elastic elastic constant tensor component 0.01
energy energy (absolute or reaction) 1.0
frac fractional coordinate in relax fit (3-D) 10000.0
freq vibrational frequency or phonon mode 1.0
grad gradient in a conventional fit 1.0
modulus bulk or Young's modulus 1.0
stress stress value in conventional fit 1.0

通常的拟合值在结构后的observable内设置，即每个独立结构（包含cell或vector以及原子坐标cood或frac），都设置observable，可以拟合的物理量包括：

elastic hfdlc sdlc energy bulk_modulus shear_modulus weight gradients hfrefractive_index srefractive_index piezoelectric frequency gradient potential entropy bornq monopoleq cv stress qreaxff fbond fangle reaction young poisson coordno sqomega volume fenergy

如果需要考虑晶格常数以及原子坐标等，可以额外添加一个晶胞结构，并打开vector或者cell后相应flags即可，注意可更改default_weight中晶格默认权重1000。如果要拟合晶胞，则关键词fit后不能添加noflag。

因此GULP中提供了势能模型参数的拟合，这里以sw势能模型为例，简单做一个介绍：

1.采用不同组能量进行拟合

关键词用fit，为了提高效率，可以再加入noflag，这样就不计算每组结构的力，仅对能量进行拟合

关键词后开始放入多组用于拟合的晶体结构和能量，每组包括描述结构的vectors或者cell，如：

vectors
7.2836 0.000000 0.000000
0.000000 9.0228 0.000000
0.000000 0.000000 20
0 0 0 0 0 0      #晶格参数的变化flag，用于晶格优化
Cartesian   #使用时默认单位angstrom
S2 core 0.899778 1.01243 1.62439 0 1 0 0 0 0    #最后3个为每个原子坐标变化的自由度flag，即在想yxz方向是否可动
S2 core 4.56721 1.01518 1.63963 0 1 0 0 0 0
Ge1 core 2.73844 2.6256 1.39913 0 1 0 0 0 0
Ge1 core 6.38759 2.65297 1.39224 0 1 0 0 0 0
.....

observable
energy
-41564.246628349894 1.0
end
shift
-41551.228686
注意：每组的能量都是放置在observable....end之间，即表示拟合用的观测量，这里是能量，用energy，由于通常采用第一性原理计算，所以拟合时能量的准确数值没有必要完全相同，因此在第一组结构后还加入shift，表示在这个量上减去偏移量，再进行拟合，这个shift也可以作为一个拟合变量，需要在最后一组结构后加入：

variables
shift
end
表明shift也是一个拟合变量

最后根据势能模型，给出势能初始参数以及是否变化的flag，比如对sw，包括2体势能sw2和sw3的3体描述角度的势能

sw2
Ge1 S2    2.090   0.350  18.371   0.000   2.926 1 0 1
S1  Ge2   2.090   0.350  18.371   0.000   2.926 1 0 1
.....

拟合时要求数据量大于变量数目，测试用4800组第一性的计算结果拟合，并且也可以在每组结构中加入space空间群的描述，就可以拟合中考虑对称性

2.采用能量和力进行拟合

同上，只是关键词不需要加入noflag

注意，在每组原子结构中

Cartesian
S2 core 0.899778 1.01243 1.62439 0 1 0 0 0 0
......

observables
energy
-41564.246628349894 1.0
gradients
1 0.006756431721086923 -0.0066830013827219575 0.003624256995566136 0.03125
....

end

并且在observables中加入gradients，对应每个原子序号以及x,y,z的受力，拟合用的权重

当然还可以加入其它一些力学性质等构建拟合的损失函数，详见官网给出的教程中observables参数

3.采用声子谱计算数据进行拟合

关键词还是fit

给出晶体结构和原子坐标描述后，在observables中加入frequency，比如：

observable
frequency 72
1 0.0000 1 1.0
2 0.0000 1 1.0
3 0.0000 1 1.0
......

end

frequency后的数字表示拟合用声子频率的数据数目，即之后数据的行数

每一行包括：声子是哪一支(晶胞内4个原子，则有4*3=12支声子) 频率拟合值(cm-1)  k点的序号  权重

其中k点序号是对应于之后kpoints参数里的k点排列位置，如最后给出kpoints的信息：

kpoints
0.000000  0.000000  0.000000
kpoints
-0.500000  0.000000  0.000000

.......