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

Lacrimosa

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为理论部分公式中的λ

1. [ defaults ]
   
2. ; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
   
3. 1               1               yes             0.5     0.8333
   
4. 
   
5. [ atomtypes ]
   
6. ; tip4pice                           C6, C12
   
7.   H_ice     1       1.008   0.0000  A   0.00000e+00  0.00000e+00
   
8.   O_ice     8       16.00   0.0000  A   0.00000e+00  0.00000e+00
   
9.   MW     0      0.0000   0.0000  D   0.00000e+00  0.00000e+00
   
10.   MOLD     6       12.00   0.0000  A   0.0000 0.0000
    
11. 
    
12. [ nonbond_params ]
    
13. ; i      j       func    C(=4 pi sig^6) A(=4 i sig^12)
    
14.   O_ice     O_ice       1       3.55885e-03  3.58949e-06
    
15.   MOLD   O_ice       1      <font color="#ff0000"><b>lambda</b></font> 0.000
    
16.   MOLD   MOLD     1       0.0000  0.0000
    
17. 
    
18. ;
    
19. ;
    
20. #include “./tip4pice.itp”
    
21. 
    
22. [ moleculetype ]
    
23. ; Name          nrexcl
    
24. MOLD         1
    
25. 
    
26. [ atoms ]
    
27. ;   nr       type      resnr residue  atom   cgnr     charge       mass
    
28.     1        MOLD    1     MOLD     MOLD      2         
    
29. 
    
30. [ system ]
    
31. Mold Integration
    
32. 
    
33. [ molecules ]
    
34. SOL         1792
    
35. 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)

6 day ago 上传 Uploaded

点击下载
Click to download

Lacrimosa

本文所用到的结构
conf.zip (84.73 KB, 下载次数 Times of downloads: 3)

6 day ago 上传 Uploaded

点击下载
Click to download