【教程】计算原子晶体的自由能-以LJ晶体为例

Lacrimosa

Einstein Crystal Method

理论部分

通过Hamiltonian integration，可以获得任意物质相对于某个参考态的自由能差，如果这个参考态的自由能具有解析形式，那么便可获得任意物质的自由能。对于晶体而言，这个参考态可以选为Einstein Crystal（一种没有粒子间相互作用的晶体，其粒子在平衡位置附近的振动可以用谐振函数表示，因此其自由能具有解析形式）。由Einstein Crystal到任意晶体的路径即为：Einstein Crystal--#1-->任意晶体+谐振势--#2-->任意晶体。可以看出#1是一个开启粒子间相互作用的过程，#2是一个关闭谐振势的过程，这两步可以均可以通过Hamiltonian integration实现。在实际操作中，通常还需要固定体系的质心，因此上述路径将变为：Einstein Crystal---->Einstein Crystal（固定质心）--#1-->任意晶体（固定质心）+谐振势--#2-->任意晶体。

公式推导详情见：J. Phys.: Condens. Matter 20 (2008) 153101
其中Einstein Crystal的自由能表达式为：

上式中q_r, q_v, q_e分别为粒子的转动，振动和电子配分函数，对于两相共存体系，可以假设其不影响两相的相对自由能，因此可以设为1。p_i为粒子i的动量，d1...dN中的数字指每个粒子相对于其平衡位置的坐标。
U_{Ein-id}为Einstein Crystal的势函数，其表达式如下：



#0 对于原子晶体，转动项（U_{Ein,or}）为0。最终Einstein Crystal的平动自由能（以NkT为单位）可以写为：

其中最后的V/N是考虑周期性边界条件后引入的项。LAMBDA为thermal de Broglie wavelength


#1 由Hamiltonian integration可知，Einstein Crystal（固定质心）--#1-->任意晶体（固定质心）+谐振势的自由能变化为：

注：这一步自由能的公式与文献中不同，在实际操作中需要很硬的弹簧才能应用文献中的公式


#2 任意晶体（固定质心）+谐振势--#2-->任意晶体的自由能变化为：



对以上三项求和，即为任意晶体的自由能。



模拟部分


以下内容基于文章：THE JOURNAL OF CHEMICAL PHYSICS 137,146101(2012)
计算LJ晶体的自由能
1.解析项的计算
在所需的温度压力下对晶体结构进行NpT平衡模拟，随后能量最小化获得该条件下的平衡位置坐标。将温度(T)，体积(V)，粒子数（N），LAMBDA和弹簧力常数（LAMBDA_E）代入#0中的公式。


2.Einstein Crystal-->Einstein Crystal + 粒子间势函数
这一步对应于理论部分的#1。首先以上一步中获得的能量最小化结构为平衡位置，使用gmx genrestr生成一个posre.itp用于position restraint(即Einstein Crystal中的弹簧)，调整限制粒子坐标的力常数到一个较大的数值。使用Gauss-Legendre quadrature formula生成14个λ以及相应的权重，对itp文件LJ参数中的epsilon进行放缩。随后在每个λ下分别进行5ns的NVT模拟（模拟中使用position restraint）. 再制作一个没有position restraint的LJ的tpr，使用该tpr对不同λ下产生的轨迹进行rerun，计算后4ns的势能,将平均值乘以相应λ下的权重后求和，即为ΔA_1。


3.Einstein Crystal + 粒子间势函数-->LJ晶体
这一步对应于理论部分的#2。使用上一步中生成的λ对posre.itp中的弹簧力常数进行放缩，完全开启粒子间的LJ势函数，在position resitraint下进行5ns的NVT模拟。随后使用gmx energy提取后4ns轨迹中position restraint项的弹簧势能，除以λ即可获得#2公式中的U_{Spring}。最后将平均值乘以相应λ下的权重后求和，取负值即为ΔA_2。




在THE JOURNAL OF CHEMICAL PHYSICS 137,146101(2012)中作者提供了GROMACS和LAMMPS的输入文件，
本文中所用到的脚本以及数据对比将在后续发布

Lacrimosa

2L留着写冰的例子