【教程】计算固液界面的自由能-以水为例

Lacrimosa · 发表于 Post on 2024-11-15 23:20:20

本帖最后由 Lacrimosa 于 2024-11-17 00:54 编辑

Mold Integration Method

理论部分

在固液相平衡条件下, 液体中形成一层固体所产生的自由能变化ΔG = 2γA，其中γ为固液界面自由能，2A为固液界面的面积（固体晶面上下两侧均与液体接触）。因此若能算出ΔG，便可知γ。由于均相成核十分困难, 因此需要使用一个势阱诱导固相的形成; 为了使粒子按晶格排列, 首先势阱需要位于晶体的平衡位置；同时这个势阱需要足够窄, 仅够一个粒子进入其中; 通过逐渐开启势阱与粒子的相互作用, 可获得一条可逆的路径（纯液体，粒子间的相互作用->固液混合的slab，粒子间的相互作用+粒子填入势阱的相互作用，如图）.

在这里势阱选为一个U型的平底势u_{pw}：

根据Hamiltonian integration方法可将始末态分别选为：

其中N为粒子数，N_w为势阱的个数，U_{pp}为粒子间的势函数，λ为放缩系数，U_{pm}为全部势阱与全部粒子相互作用的势函数，其可以写为每个势阱与每个粒子的相互作用之和。

随后通过对λ进行积分即可获得两个状态间的自由能ΔG^m，其中<U_{pm}>=被占据的势阱数(N_{fw})×势阱深度(ε) 。

通过以上步骤求得的自由能变化ΔG^m实际上包含了两部分： (相界面形成带来的自由能变化)ΔG+(势阱被填充带来的自由能变化)ΔG_w=2γA - N_w*ε ，N_w为总势阱数。

最终表面张力可写为：

上式中p_x表示半各向异性控压条件（固液共存体系需要该设置）

模拟部分

模拟数据全部来源于此文： J. Phys. Chem. C 2016, 120, 8068−8075
该计算方法目前已经整合进LAMMPS: doi:10.21105/joss.06083
在GROMACS中进行计算需要用到列表势（关于列表势的介绍见：【教程】如何在gmx中使用列表势 http://bbs.keinsci.com/forum.php ... 35740&fromuid=10773）

1. 准备初始结构：

准备一个冰Ih的晶体：可通过genice程序获得，二楼也提供了一个本文中所用的结构

将所考察的晶面转向z轴：对于晶胞的{001}，{010}，{100}面，可以使用vmd实现旋转，例如：在tkconsole中输入：[atomselect top all] move [transaxis x 90] 可实现绕x轴旋转90度，y/z轴以此类推。转动到合适位置后再使用pbc set {a b c alpha beta gamma}调整盒子三个方向的参数即可。若是{111}这种晶面，则需要重新定义晶胞，一个可行的方案是使用M$里的Redefine Lattice实现。

设置势阱的位置：在VMD中打开转动好的晶体结构，选择其中任意一层（位于z轴中央比较好，后面作图的时候看着方便），将其resname设置为MOLD（在tkconsole中输入：[atomselect top “z>18 and z<22 and name OW”] set resname MOLD），保存为conf.pdb。使用gromacs做一个该结构的index.ndx留作备用。在后续的模拟中，resname为MOLD的原子位置将用于放置势阱。

2. 找到熔点下MOLD的平衡位置：
对#1中准备好的结构进行能量最小化，在熔点下跑1ns NPT平衡模拟，找到该晶面在熔点下的表面积A。最后再做一次能量最小化，使所有原子都达到该密度下的平衡位置。最后提取出该结构的MOLD坐标（使用#1中的index.ndx），保存为mold.pdb

3.构建水-势阱体系：
创建一个和#2中盒子尺寸同等大小的水盒子，保持水分子数量与#2中一致。进行能量最小化后，将结构与上一步的mold.pdb合并。

4.创建列表势：
见附件python脚本，根据文章中的数据，r_w=0.083 nm alpha=0.0017

相关的mdp设置见： http://bbs.keinsci.com/forum.php ... 35740&fromuid=10773
此外还需使用半各向同性parrinello-rahman控压，cutoff均设为1.4 nm

top文件设置可以参考这里:
lambda为理论部分公式中的λ

[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 1 yes 0.5 0.8333
[ atomtypes ]
; tip4pice C6, C12
H_ice 1 1.008 0.0000 A 0.00000e+00 0.00000e+00
O_ice 8 16.00 0.0000 A 0.00000e+00 0.00000e+00
MW 0 0.0000 0.0000 D 0.00000e+00 0.00000e+00
MOLD 6 12.00 0.0000 A 0.0000 0.0000
[ nonbond_params ]
; i j func C(=4 pi sig^6) A(=4 i sig^12)
O_ice O_ice 1 3.55885e-03 3.58949e-06
MOLD O_ice 1 <font color="#ff0000"><b>lambda</b></font> 0.000
MOLD MOLD 1 0.0000 0.0000
;
;
#include “./tip4pice.itp”
[ moleculetype ]
; Name nrexcl
MOLD 1
[ atoms ]
; nr type resnr residue atom cgnr charge mass
1 MOLD 1 MOLD MOLD 2
[ system ]
Mold Integration
[ molecules ]
SOL 1792
MOLD 128

复制代码

5. Mold Integration：
在不同lambda（见#4 top设置）下进行1.5ns NPT模拟, lambda选取范围为0-8RT

6.结果分析：
#5中获得了一系列lambda下的轨迹，取后1ns，统计填入势阱内的分子数(当势阱到氧原子距离小于r_w的时候即认为水分子进入势阱)。最后对数据进行积分，带入理论部分中的公式计算即可。

模拟所用到的脚本以及数据对比将在后续发布

table.py (1.31 KB, 下载次数 Times of downloads: 3)

Lacrimosa · 发表于 Post on 7 day ago

本帖最后由 Lacrimosa 于 2024-11-17 00:53 编辑

本文所用到的结构

conf.zip (84.73 KB, 下载次数 Times of downloads: 3)

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[GROMACS] 【教程】计算固液界面的自由能-以水为例

评分 Rate