计算化学公社

标题: 【教程】计算固液界面的自由能-以水为例 [打印本页]
作者
Author: Lacrimosa    时间: 2024-11-15 23:20
标题: 【教程】计算固液界面的自由能-以水为例
本帖最后由 Lacrimosa 于 2024-11-17 00:54 编辑


Mold Integration Method

理论部分

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

在这里势阱选为一个U型的平底势u_{pw}：
(, 下载次数 Times of downloads: 1)

根据Hamiltonian integration方法可将始末态分别选为：
(, 下载次数 Times of downloads: 0)

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

随后通过对λ进行积分即可获得两个状态间的自由能ΔG^m，其中<U_{pm}>=被占据的势阱数(N_{fw})×势阱深度(ε) 。
(, 下载次数 Times of downloads: 0)
通过以上步骤求得的自由能变化ΔG^m实际上包含了两部分： (相界面形成带来的自由能变化)ΔG+(势阱被填充带来的自由能变化)ΔG_w=2γA - N_w*ε ，N_w为总势阱数。

最终表面张力可写为：
(, 下载次数 Times of downloads: 0)
上式中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

  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

  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

  18. ;
  19. ;
  20. #include “./tip4pice.itp”

  22. [ moleculetype ]
  23. ; Name          nrexcl
  24. MOLD         1

  26. [ atoms ]
  27. ;   nr       type      resnr residue  atom   cgnr     charge       mass
  28.     1        MOLD    1     MOLD     MOLD      2         

  30. [ system ]
  31. Mold Integration

  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的时候即认为水分子进入势阱)。最后对数据进行积分，带入理论部分中的公式计算即可。


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


(, 下载次数 Times of downloads: 3)










作者
Author: Lacrimosa    时间: 7 day ago
本帖最后由 Lacrimosa 于 2024-11-17 00:53 编辑

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



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3