|
|
本帖最后由 Lacrimosa 于 2024-11-15 15:22 编辑
如有疏漏请务必指出，我将及时更新，感谢各位。
理论篇
1.自由能计算方法:
Thermodynamic Integration
沿等温线积分：
沿等压线积分:
沿等容线积分:
Hamiltonian Integration
只要在两个状态(例如：LJ流体到水, 或者理想气体到LJ流体)间找到任意一条可逆的路径(注意要避免相变), 对该路径进行积分即可求出两个状态间的自由能变化. 此处通过定义一个耦合参数λ, 将两个状态的哈密顿量关联起来。
上式中q为粒子的配分函数。
2.水的自由能计算
液态水自由能的绝对值可以通过Hamiltonian integration进行计算。 其中始末状态分别选为LJ流体以及水即可。目前主流的水模型中均只有氧原子有LJ参数, 氢原子为点电荷。 因此只需要逐渐关闭水分子的静电作用, 即减小原子电荷, 便可实现从LJ流体到水的可逆变化。 该过程的自由能如下:
其中, 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(λ).
操作篇
以下内容基于文章：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设置，为了更接近原理这里没有用到。
跑平衡过程
- Pressure=("743" "4280") #bar
- Temperature=("225" "443") #K
- N_site=4 #3 for 3 site model, 4 for 4 site model.
- water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...
- create_system=yes
- N=360
- density=1.05 #g/cm^3
- box_size=1.9
- conf=conf
- top=topo
- state=liquid #solid or liquid
- Minimization=yes
- Equilibration=yes
- #--Hamiltonian Integration---
- #
- Integration=no
- num_points=14 #number of point for integration
- Analyze=no #compute integral
- #----------------------------
- #--Thermodynamic Integration---
- #
- TI=no
- #------------------------------
Hamitlonian Integration
- Pressure=("743" "4280") #bar
- Temperature=("225" "443") #K
- N_site=4 #3 for 3 site model, 4 for 4 site model.
- water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...
- create_system=no
- N=360
- density=1.05 #g/cm^3
- box_size=1.9
- conf=conf
- top=topo
- state=liquid #solid or liquid
- Minimization=no
- Equilibration=no
- #--Hamiltonian Integration---
- #
- Integration=yes
- num_points=14 #number of point for integration
- Analyze=no #compute integral
- #----------------------------
- #--Thermodynamic Integration---
- #
- TI=no
- #------------------------------
分析数据：结果见路径下的results.dat
- Pressure=("743" "4280") #bar
- Temperature=("225" "443") #K
- N_site=4 #3 for 3 site model, 4 for 4 site model.
- water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...
- create_system=no
- N=360
- density=1.05 #g/cm^3
- box_size=1.9
- conf=conf
- top=topo
- state=liquid #solid or liquid
- Minimization=no
- Equilibration=no
- #--Hamiltonian Integration---
- #
- Integration=no
- num_points=14 #number of point for integration
- Analyze=yes #compute integral
- #----------------------------
- #--Thermodynamic Integration---
- #
- TI=no
- #------------------------------
5.Thermodynamic Integration：from 225K to 443K
从200K至450K每隔50K对#1中跑平衡的结构进行NVT模拟5ns，计算每个温度下的平均Potential，拟合U/NkT^2随温度的变化，带入等容线方程进行积分，最终可获得一条自由能差（以200K下水的自由能为0点）随温度的变化曲线，将225K和443K下的自由能带入，即为水的自由能随温度的变化曲线。
脚本设置如下：
跑平衡：
- Pressure=("none" "none" "none" "none" "none" "none") #bar
- Temperature=("200" "250" "300" "350" "400" "450") #K
- N_site=4 #3 for 3 site model, 4 for 4 site model.
- water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...
- create_system=yes
- N=360
- density=1.05 #g/cm^3
- box_size=1.9
- conf=conf
- top=topo
- state=liquid #solid or liquid
- Minimization=yes
- Equilibration=yes
- #--Hamiltonian Integration---
- #
- Integration=no
- num_points=14 #number of point for integration
- Analyze=yes #compute integral
- #----------------------------
- #--Thermodynamic Integration---
- #
- TI=no
- #------------------------------
- Pressure=("none" "none" "none" "none" "none" "none") #bar
- Temperature=("200" "250" "300" "350" "400" "450") #K
- N_site=4 #3 for 3 site model, 4 for 4 site model.
- water_model=tip4p #spce/tip4p/tip4p2005/tip4pice...
- create_system=no
- N=360
- density=1.05 #g/cm^3
- box_size=1.9
- conf=conf
- top=topo
- state=liquid #solid or liquid
- Minimization=no
- Equilibration=no
- #--Hamiltonian Integration---
- #
- Integration=no
- num_points=14 #number of point for integration
- Analyze=yes #compute integral
- #----------------------------
- #--Thermodynamic Integration---
- #
- TI=yes
- #------------------------------
数据分析见图：分析数据的jupyter-notebook见附件
TIP4P2005和TIP4PIce的拓扑文件参考了卢天老师的写法，在此表示感谢。
simulation.sh
(24.25 KB, 下载次数 Times of downloads: 2)
Analyze.ipynb
(112.17 KB, 下载次数 Times of downloads: 3)
|
评分 Rate
-
查看全部评分 View all ratings
|
发表于 Post on 2024-11-14 23:21:25
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author