计算化学公社

标题: 【教程】使用GROMACS计算液体的自由能--以水为例 [打印本页]
作者
Author: Lacrimosa    时间: 2024-11-14 23:21
标题: 【教程】使用GROMACS计算液体的自由能--以水为例
本帖最后由 Lacrimosa 于 2025-6-19 15:54 编辑


如有疏漏请务必指出，我将及时更新，感谢各位。


理论篇
1.自由能计算方法:

Thermodynamic Integration
沿等温线积分：
(, 下载次数 Times of downloads: 81)
沿等压线积分:
(, 下载次数 Times of downloads: 75)
沿等容线积分:
(, 下载次数 Times of downloads: 79)

Hamiltonian Integration
只要在两个状态(例如：LJ流体到水, 或者理想气体到LJ流体)间找到任意一条可逆的路径(注意要避免相变), 对该路径进行积分即可求出两个状态间的自由能变化. 此处通过定义一个耦合参数λ, 将两个状态的哈密顿量关联起来。
 (, 下载次数 Times of downloads: 77)
上式中q为粒子的配分函数。

2.水的自由能计算
液态水自由能的绝对值可以通过Hamiltonian integration进行计算。 其中始末状态分别选为LJ流体以及水即可。目前主流的水模型中均只有氧原子有LJ参数, 氢原子为点电荷。 因此只需要逐渐关闭水分子的静电作用, 即减小原子电荷, 便可实现从LJ流体到水的可逆变化。 该过程的自由能如下:
 (, 下载次数 Times of downloads: 80)
其中, U_{Coul}为体系的静电势能, 任意温度压力下, LJ流体的自由能可以通过状态方程算出(见Johnson J K, Zollweg J A and Gubbins K E 1993 Mol. Phys. 78-591. The Lennard-Jones equation of state revisited.) 以上公式的含义为: 首先在不同的λ下进行MD模拟产生模拟轨迹, 随后对这些轨迹计算λ=1下的静电势能, 最后求积分即可. 在MD模拟中, 逐步关闭静电作用是通过减小原子电荷实现的, 即给电荷乘一个缩放因子, 由库仑定律可知, 电荷的缩放因子应为sqrt(λ).
 (, 下载次数 Times of downloads: 76)


操作篇
以下内容基于文章：J.Phys.:Condens.Matter 20(2008)153101
1.建模跑平衡:
创建一个360个水分子的水盒子（TIP4P水模型），盒子尺寸为2.172nm (密度为1.05g/cm^3) , 在225K下进行2ns的NVT模拟.


2.Hamitlonian Integration：from LJ to liquid water
使用Gauss-Legendre quadrature formula生成14个λ以及相应的权重（10个也可以），将TIP4P的itp文件中的原子电荷乘以sqrt(λ)，随后在每个λ下分别进行能量最小化, 以及5ns的NVT模拟. 随后使用TIP4P原始的itp制作tpr文件，对产生相中的11段轨迹进行rerun, 计算后4ns的平均静电势能,将平均值乘以相应λ下的权重后求和，最终获得225K下水与LJ的自由能差.


3.使用LJ的状态方程计算225K，1.05g/cm^3下的自由能，与#2中结果相加，即为该条件下水的自由能


4.重复在443K同样密度下重复#1-#3，最终获得443K下水的自由能

脚本设置：对附件中simulation.sh中的setting部分进行如下修改即可。注意GROMACS中带有调节λ的mdp设置，为了更接近原理这里没有用到。
跑平衡过程
  1. Pressure=("743" "4280") #bar
  2. Temperature=("225" "443") #K

  4. N_site=4 #3 for 3 site model, 4 for 4 site model.
  5. water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...

  7. create_system=yes
  8. N=360
  9. density=1.05 #g/cm^3 如果要做npt模拟则不需要指定密度，直接写no就可以。
  10. box_size=1.9
  11. conf=conf
  12. top=topo
  13. state=liquid #solid or liquid

  15. Minimization=yes
  16. Equilibration=yes

  18. #--Hamiltonian Integration---
  19. #
  20. Integration=no
  21. num_points=14 #number of point for integration
  22. Analyze=no #compute integral
  23. #----------------------------

  25. #--Thermodynamic Integration---
  26. #
  27. TI=no
  28. #------------------------------
复制代码



Hamitlonian Integration
  1. Pressure=("743" "4280") #bar
  2. Temperature=("225" "443") #K

  4. N_site=4 #3 for 3 site model, 4 for 4 site model.
  5. water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...

  7. create_system=no
  8. N=360
  9. density=1.05 #g/cm^3
  10. box_size=1.9
  11. conf=conf
  12. top=topo
  13. state=liquid #solid or liquid

  15. Minimization=no
  16. Equilibration=no

  18. #--Hamiltonian Integration---
  19. #
  20. Integration=yes
  21. num_points=14 #number of point for integration
  22. Analyze=no #compute integral
  23. #----------------------------

  25. #--Thermodynamic Integration---
  26. #
  27. TI=no
  28. #------------------------------
复制代码




分析数据：结果见路径下的results.dat
  1. Pressure=("743" "4280")    #bar
  2. Temperature=("225" "443")  #K

  4. N_site=4 #3 for 3 site model, 4 for 4 site model.
  5. water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...

  7. create_system=no
  8. N=360
  9. density=1.05 #g/cm^3
  10. box_size=1.9
  11. conf=conf
  12. top=topo
  13. state=liquid #solid or liquid

  15. Minimization=no
  16. Equilibration=no

  18. #--Hamiltonian Integration---
  19. #
  20. Integration=no
  21. num_points=14 #number of point for integration
  22. Analyze=yes #compute integral
  23. #----------------------------

  25. #--Thermodynamic Integration---
  26. #
  27. TI=no
  28. #------------------------------

复制代码




5.Thermodynamic Integration：from 225K to 443K
从200K至450K每隔50K对#1中跑平衡的结构进行NVT模拟5ns，计算每个温度下的平均Potential，拟合U/NkT^2随温度的变化，带入等容线方程进行积分，最终可获得一条自由能差（以200K下水的自由能为0点）随温度的变化曲线，将225K和443K下的自由能带入，即为水的自由能随温度的变化曲线。


脚本设置如下：
跑平衡：
  1. Pressure=("none" "none" "none" "none" "none" "none") #bar
  2. Temperature=("200" "250" "300" "350" "400" "450") #K

  4. N_site=4 #3 for 3 site model, 4 for 4 site model.
  5. water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...

  7. create_system=yes
  8. N=360
  9. density=1.05 #g/cm^3
  10. box_size=1.9
  11. conf=conf
  12. top=topo
  13. state=liquid #solid or liquid

  15. Minimization=yes
  16. Equilibration=yes

  18. #--Hamiltonian Integration---
  19. #
  20. Integration=no
  21. num_points=14 #number of point for integration
  22. Analyze=yes #compute integral
  23. #----------------------------

  25. #--Thermodynamic Integration---
  26. #
  27. TI=no
  28. #------------------------------
复制代码
分析数据：见results.dat
  1. Pressure=("none" "none" "none" "none" "none" "none")    #bar
  2. Temperature=("200" "250" "300" "350" "400" "450")  #K

  4. N_site=4 #3 for 3 site model, 4 for 4 site model.
  5. water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...

  7. create_system=no
  8. N=360
  9. density=1.05 #g/cm^3
  10. box_size=1.9
  11. conf=conf
  12. top=topo
  13. state=liquid #solid or liquid

  15. Minimization=no
  16. Equilibration=no

  18. #--Hamiltonian Integration---
  19. #
  20. Integration=no
  21. num_points=14 #number of point for integration
  22. Analyze=yes #compute integral
  23. #----------------------------

  25. #--Thermodynamic Integration---
  26. #
  27. TI=yes
  28. #------------------------------
复制代码


数据分析见图：分析数据的jupyter-notebook见附件 (, 下载次数 Times of downloads: 71)



TIP4P2005和TIP4PIce的拓扑文件参考了卢天老师的写法，在此表示感谢。

2025/6/19更新了该文件，修复了一些bug。
(, 下载次数 Times of downloads: 15)


(, 下载次数 Times of downloads: 9)
作者
Author: Lacrimosa    时间: 2024-11-14 23:36
本帖最后由 Lacrimosa 于 2025-7-22 20:46 编辑

注意事项：
1. 在高压下计算的时候需要注意关闭电荷的时候LJ流体可能会结晶，导致积分路径不连续，过程不可逆。 比较好的解决方式是在低压下计算绝对值，然后沿等温线积分到高压区。

作者
Author: popelrain    时间: 2025-5-14 21:38
非常感谢分享~



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