前言
本教程以甲烷在水中的溶剂化自由能和配体与受体蛋白质结合自由能为入门和进阶例子对使用分子动力学模拟方法来计算自由能的常见方法给出了示例，教程不会过多讲解各种计算自由能方法原理和分析原理，关于自由能计算原理和方法可以参考自由能专题1：原理及常见方法分析原理可以常考自由能专题2：计算与分析指南。该教程全部都经过本人的亲自的模拟和分析，确保无误，从实际操作出发，一步一步手动实现自由能的计算与分析，有助于理解自由能计算原理和结果分析。本教程主要采用GROMACS (5.1.4)分子动力学软件包来进行操作。以下案例的全部模拟配置文件均可在问末的链接中获得。
1.常见自由能计算方法和工具链
分子模拟计算自由能的方法众多，这里简单整理出一个表格，来展示本教程将要运用的自由能计算方法和后处理方法及相关的工具链。如下：
| 计算方法 | 后处理方法 | 自动实现分析工具/脚本 |
|---|---|---|
| PMF | PMF(表达式) | 图表数据处理工具 |
| Umbrella Sampling(US:求解PMF) | WHAM, | gmx wham（gromacs自带） ,PLUMED,Colvars |
| BAR | gmx bar(自带)，alchemical_analysis | |
| MBAR(0宽度的WHAM) | alchemical_analysis | |
| FEP(Zwanzig) | DEXP | alchemical_analysis（下同） |
| IEXP | ||
| GINS | ||
| GDEL | ||
| BAR | ||
| UBAR | ||
| RBAR | ||
| MBAR | ||
| TI(Kirkwood) | TI | alchemical_analysis |
| TI-CUBIC | ||
| ＭＭ/PBSA | gmx_mm/pbsa（脚本） | |
| g_mm/pbsa（脚本） | ||
| Jarzynski Equality(Jarzynski) | gmx analyze（gromacs自带） | |
| alchemical_analysis | ||
| Metadynamics（MtD） | PLUMED | |
| Colvars | ||
3.自由能计算实例
3.1 PMF：计算甲烷-甲烷之间的平均力势
在本教程中, 我将向您展示如何创建一个包含多个OPLS甲烷的系统, 将其置于TIP4PEW水盒子中, 并根据模拟所得信息计算甲烷-甲烷之间的平均力势.参考Barnett GROMACS教程3: 水中的多个甲烷,既然主题是计算自由能了，那么默认实操者能独立进行简单的分子模模拟任务，并熟知模拟流程，在这里关于模拟所需的三种重要基本配置文件（结构文件(.gro/.pdb), 拓扑文件(.itp,.top)和参数文件(.mdp).）如何创建就不再赘述了。这里会直接给出，只做简单的展示。（下同）
设置
甲烷的拓扑文件如下：
; Topology for methane in TIP4P
#include "oplsaa.ff/forcefield.itp"
#include "oplsaa.ff/tip4pew.itp"
[ moleculetype ]
; Name nrexcl
Methane 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
1 opls_138 1 ALAB CB 1 -0.240 12.011
2 opls_140 1 ALAB HB1 2 0.060 1.008
3 opls_140 1 ALAB HB2 3 0.060 1.008
4 opls_140 1 ALAB HB3 4 0.060 1.008
5 opls_140 1 ALAB HB4 5 0.060 1.008
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 5 1
[ angles ]
; ai aj ak funct c0 c1 c2 c3
2 1 3 1
2 1 4 1
2 1 5 1
3 1 4 1
3 1 5 1
4 1 5 1
; water topology
[ system ]
; Name
Methane in water
[ molecules ]
; Compound #mols
Methane 10
1.首先, 将10个甲烷插入3.2 nm x 3.2 nm x 3.2 nm的盒子中:
$ gmx insert-molecules -ci methane.pdb -o box.gro -nmol 10 -box 3.2 3.2 3.2
2.溶剂化向盒子里加入1000个水分子
$ gmx solvate -cs tip4p -cp box.gro -o conf.gro -p topol.top -maxsol 1000
gmx solvate会自动更新topol.top
3.参数文件
以上过程的所有的文件均可以在1得到：关于PMF配置文件如下：
configura文件：模拟使用的gromacs5.1.4版本
├── methane.pdb :甲烷结构文件
├── box.gro :步骤1生成的文件
├── conf.gro :步骤2生成的文件
├── min1.mdp :第一次能量最小化参数文件
├── min2.mdp :第二次能量最小化
├── eql.mdp :NVT 预平衡
├── eql2.mdp :NPT 预平衡
├── prd.mdp :成品模拟，这里注意你可以改一下时长，我自己仅跑了5ns用于展示，建议跑100ns,rdf更平滑
├── topol.top :拓扑文件，包含10个甲烷和1000个水分子
└── run_pmf.bsh :自动化执行分子模拟流程的脚本
模拟
1.使用run_pmf.bsh脚本进行自动化模拟：
#1.给脚本赋予可执行的权限
$ chmod +x run_pmf.sh
#2.在当前目录下（configura）执行./run_pmf.bash脚本
$ ./run_pmf.sh
分析
1.让我们计算一下所谓的径向分布函数. 首先, 我们需要创建一个索引文件(index.ndx):
#以甲烷中的C原子代表整个分子
$ gmx make_ndx -f conf.gro
> a C
> q
2.计算甲烷-甲烷之间径向分布函数
$ gmx rdf -f prd.xtc -n index.ndx
然后选择 Group 6 “C” 两次，即选择两次6，Ctrl+d 计算生成rdf.xvg 文件
3.安装绘图工具gnuplot,在linux的ubuntu 系统安装：
#在终端输入以下命令自动安装：
$ sudo apt-get install gnuplot
4.绘制RDF(g®)图：
rdf曲线平滑依赖于你保存了多少帧, 运行了多长时间的模拟(原教程运行了100 ns, 保存了25k帧), 如果只是简单地绘制RDF, 得到的结果应该如下所示:（这里放原教程的图）
请注意, 它没有收敛到1. 我们多算了一个, 所以我们应该对其进行修正. 下面是gnuplot命令:
#1.在终端启动gnuplot
$ gnuplot
#2.进行修正:(我们有10个甲烷, 但计算RDF时只用了9个, 所以我们将g的值乘以N/(N-1), 即10/9. 它应该如下所示)
#输入>后面的代码即可
gnuplot> plot 'rdf.xvg' u 1:($2*10/9) w l
修正后的图像：
5.绘制甲烷-甲烷的平均力势（PMF）图：
现在, 要从甲烷-甲烷RDF中得到甲烷-甲烷之间的平均力势, 我们要计算：$ w=−kTln(g)w=−kTln(g).$
因此, 要用gnuplot绘图, 请执行以下操作:
#1.在终端启动gnuplot
$ gnuplot
#2.绘制w(r)
gnuplot> plot 'rdf.xvg' u 1:(-8.314e-3*298.15*log($2*10/9)) w l
得到的图如下：
总结
注意到我们的抽样存在一个小的空白区. 这意味着, 在模拟中甲烷没有在那个距离有过相互作用. 我们将在后面的教程中使用伞形采样来解决这个问题。
3.2 Umbrella Sampling(US)：使用伞形采样计算甲烷-甲烷的PMF
在本教程中, 我们将使用窗口采样, 有时也称为伞状采样, 来计算甲烷-甲烷的平均力势PMF. 在教程1中, 我们通过模拟几个甲烷分子计算径向分布函数的方法直接获得了PMF. 像这样对几个溶质分子进行采样的方法并不总是可行的.
本教程参考gromacs官方伞形采样(GMX 5.0)和Barnett GROMACS教程5: 使用伞形采样计算甲烷-甲烷的平均力势PMF
模拟所需的甲烷分子拓扑文件如上，同这两个参考教程不同之处是，本教程使用更清晰文件组织来进行模拟，有利于初学者的学习。
设置
1.创建盒子
#将两个甲烷分子插入盒子中
$ gmx insert-molecules -ci methane.pdb -o box.gro -nmol 2 -box 3.2 3.2 3.2
2.溶剂化
$ gmx solvate -cs tip4p -cp box.gro -o conf.gro -p topol.top
3.创建索引文件：
我们还需要创建一个索引文件, 其中包含我们感兴趣的两个组, 它们之间受伞形势限制. 使用gmx make_ndx创建一个索引文件, 其中包含每个甲基碳的组, 并将其分别命名为CA和CB.
$ gmx make_ndx -f conf.gro
> splitres 2 #本人操作是时CH4残基在组2，总之选择CH4残基所在组即可，将残基开成两个CH4（6和7）
> 6 & a C #单独取出两个甲烷分子中的C原子（下同） 生成 8
> 7 & a C #生成 9
> name 8 CA #重命名 8 （下同)
> name 9 CA
>q #保存退出 ，生成index.ndx索引文件
使用gmx make_ndx的其他方法同样可以达到目的. 关键是, 需要将每个甲烷的碳原子放在单独的索引组中, 并分别将它们命名为CA和CB.
4.参数文件
我们在很大程度上重复使用前一个教程（PMF）中的参数文件, 但是我们将添加一个与质心(COM)牵引有关的部分. 牵引代码可以将甲烷维持在指定距离。可能有几种不同的设置方法, 但对于这个系统, 我们将手动指定两个甲烷间的距离.我们将在一般模拟时所需要的所有mdp运行控制文件（能量最小化，npt和nvt预平衡，成品模拟mdp文件）（如1教程 PMF）中添加以下代码：
;{
US采样代码：质心牵引
;===============================================================================
pull = Yes ; Yes/No 使用牵引代码
pull-ngroups = 2 ; 牵引组数目,这里是CA和CB两个
pull-ncoords = 1 ; 只沿一个坐标方向进行牵引
pull-group1-name = CA ; 牵引组名称, 来自top或ndx；.对我们来说,
; 这是其中一个甲烷的碳, 尽管我们也可以选择整个甲烷.
;如果那样做的话, 就会沿着两个甲烷分子的质心方向进行牵引
pull-group2-name = CB ;另一个甲烷的碳原子
pull-coord1-type = umbrella
; umbrella:简谐势
; constraint:刚性约束
; constant-force:恒力, 无须pull-init和pull-rate
pull-coord1-geometry = distance
; distance:沿质心间的矢量, 用pull-dim选择分量
; direction:沿pull-vec方向
; direction-periodic:同direction, 但距离可超过盒长一半. 使用时盒子在牵引方向不能变 化
; cylinder:圆柱
pull-coord1-groups = 1 2 ; 对此坐标牵引的两个组(定义见后面). 实际上你可以定义
;更多的牵引坐标, 也可以在不同的分子组之间进行牵引, 但我们在这里不使用这些选项
pull-coord1-rate = 0.0 ; 速率(nm/ps)
;我们不需要两个组沿坐标方向运动, 因此设置为零
pull-coord1-k = 5000.0 ; 力常数(kJ/mol-nm^2)
pull-coord1-init = WINDOW ; 零时刻的基准距离,这是个变量
;两个组的初始间距. 这是关键字, 我们将使用bash脚本将其替换为每个窗口的相应值
pull-coord1-start = no ; No:不修改pull-init; Yes:将初始构型的质心距离添加到pull-init
;我们手动指定每个窗口的距离, 因此计算时不需要考虑质心距离
;}==============================================================================
5.生成构型文件
为了运行伞形采样, 我们必须沿反应坐标 ζ 生成一系列构型. 其中的一些构型将作为伞形采样窗口的初始构型, 并进行独立的模拟. 下图解释了这些原则. 图的上部显示了我们将要运行的牵引模拟, 目的在于产生一系列沿反应坐标的构型, 模拟完成后会抽取这些构型(图片上部和中部之间的虚线箭头). 图的中部对应着每一采样窗口内的独立模拟, 其中使用伞形偏离势将自由甲烷的质心限制在窗口中. 图的底部展示了构型直方图的理想结果, 相邻窗口存在重叠, 这样以后就可以从这些模拟得到连续的能量函数.
对本教程的例子, 反应坐标为单一方向. 为产生这些构型, 我们必须牵引甲烷CA, 使其远离甲烷CB.我们模拟了26个窗口, 距离从0.05 nm到大约1.3 nm. 请注意, 我为每次模拟都添加了牵引力输出文件选项-pf和牵引距离输出文件选项-px. 这是因为使用-deffnm选项时GROMACS会覆盖这两个文件. 此外, 我使用-n指定了索引文件, 因为gmx grompp需要获取我们使用牵引参数指定的组。
对于其他体系, 保存构型的频率可能需要更高或者更低, 只要足够就好. 关键是沿反应坐标获得足够的构型, 以使伞形采样窗口的间距正常, 窗口间距以甲烷CA和CB的质心(用C来代替了)距离表示, 其中CB为参考组.
由于我们需要为26个窗口进行独立的模拟，那么就需要27份模拟所需要的配置文件，但是我们将使用脚本来简单地更新参数文件模板中的适当值（pull-coord1-init）, 因此实际上只需要为模拟过程的每个部分提供一个模板。以上过程所需及生成的文件从1可以获得，其配置文件结构如下：
**configura 文件下的目录结构**
├── collect_dat.sh :批量收集分析所需要数据的脚本
├── batche.sh :批量执行job0.sh -->job26.sh 的脚本
├── write_mdp.sh :批量生成mdp文件的脚本，主要修改mdp 中的pull-coord1-init=WINDOWS 的值
├── write_sh.sh :批量生成job文件的脚本,主要修改job文件中 WINDOW=Wnum的值，
├── job.sh :为每一窗口独立运行整个模拟过程（EM->NVT->NPT->PRD），即WINDOW0，
├── MDP :总运行参数mdp文件夹
│ ├── EM :能量最小化参数mdp文件夹
│ │ ├── 1em.mdp :第1次能量最小化参数mdp模板文件
│ │ ├── 2em.mdp :第2次能量最小化参数mdp模板文件
│ │ └── write_mdp.sh : 用法：./write_mdp.sh 1em 即可后接mdp文件前缀 ./write_mdp.sh 2em
│ ├── NVT :等温等容预平衡参数mdp文件夹
│ │ ├── nvt.mdp :等温等容预平衡参数mdp模板文件
│ │ ├── write_mdp.sh : 用法：./write_mdp.sh nvt 即可后接mdp文件前缀
│ ├── NPT :等温等压预平衡参数mdp文件夹
│ │ ├── npt.mdp :等温等压预平衡参数mdp模板文件
│ │ └── write_mdp.sh : 用法：./write_mdp.sh npt 即可后接mdp文件前缀
│ └── PRD :成品模拟参数mdp文件夹
│ ├── prd.mdp :成品模拟参数模板mdp文件
│ └── write_mdp.sh : 用法：./write_mdp.sh prd 即可后接mdp文件前缀
├── Methane :结构文件和拓扑文件夹
│ ├── box.gro :10个甲烷盒子结构文件
│ ├── conf.gro :模拟体系结构文件10个甲烷+水
│ ├── index.ndx :甲烷中所有原子的索引文件
│ ├── methane.pdb :甲烷结构文件
│ └── topol.top :体系top文件
└── REDERME.md
5.1.使用write_mdp.sh脚本在文件夹EM/NVT/NPT/PRD下分别生成0到26个.mdp文件
脚本如下:
#!/bin/bash
#对于模拟, 我们将使用shell脚本替换.mdp文件中的WINDOW关键字
# 打开一个通用的.mdp文件，并替换以下值
#26个窗口, 距离从0.05 nm到大约1.3 nm.既WINDOW有27个值，增量为0.05
for ((i=0;i<27;i++))
do
# 小数点前补0
x=$(echo "0.05*$(($i+1))" | bc |awk '{printf "%.2f",$0}');
sed 's/WINDOW/'$x'/g' $1.mdp >$1$i.mdp
done
使用示例：例如在./configura/MDP/EM下批量生成1em.mdp和2em.mdp各27份
#给脚本赋予可执行权限
$chmod +x write_mdp.sh
#批量生成27份1em.mdp文件
$ ./write_mdp.sh 1em
#批量生成27份2em.mdp文件
$ ./write_mdp.sh 2em
接下来以上述方式生成整个模拟所需要的mdp文件(同上)
模拟
1.对于每一个独立窗口我们可以用脚本（job.sh）进行自动化的模拟整个过程如下：
#!/bin/bash
# 设置一些环境变量，即你所有模拟配置文件所在的根目录（里面包括MDP，Methane,.sh脚本，在命令行使用pwd即可获得）
#"/home/eric/free_energy/04.US/configura" 即为我所有模拟配置文件所在根目录
FREE_ENERGY=/home/eric/free_energy/04.US/configura
echo "Free energy home directory set to $FREE_ENERGY"
MDP=$FREE_ENERGY/MDP
echo ".mdp files are stored in $MDP"
#Wnum 为窗口数，将会被write_sh.sh进行批量替换成27个不同的值
WINDOW=Wnum
# 将为WINDOW的每个值和工作流程中的每个步骤创建一个新的目录，以实现最大限度的组织。
mkdir Window_$WINDOW
cd Window_$WINDOW
#################################
# 能量最小化第1步:STEEP #
#################################
echo "Starting minimization for Window = $WINDOW..."
mkdir EM_1
cd EM_1
# Iterative calls to grompp and mdrun to run the simulations
gmx grompp -f $MDP/EM/1em$WINDOW.mdp -c $FREE_ENERGY/Methane/conf.gro -p $FREE_ENERGY/Methane/topol.top -n $FREE_ENERGY/Methane/index.ndx -o min$WINDOW.tpr -maxwarn 1
gmx mdrun -nt 2 -v -deffnm min$WINDOW -pf pullf-min$WINDOW -px pullx-min$WINDOW
sleep 10
#################################
# 能量最小化第 2步: STEEP #
#################################
cd ../
mkdir EM_2
cd EM_2
#我们在这里使用了 -maxwarn 1，因为 grompp 错误地抱怨使用了一个普通的截止值。
# 这是一个小问题，将在未来的Gromacs版本中得到修正。
gmx grompp -f $MDP/EM/2em$WINDOW.mdp -c ../EM_1/min$WINDOW.gro -p $FREE_ENERGY/Methane/topol.top -n $FREE_ENERGY/Methane/index.ndx -o min$WINDOW.tpr -maxwarn 1
gmx mdrun -nt 2 -v -deffnm min$WINDOW -pf pullf-min$WINDOW -px pullx-min$WINDOW
echo "Minimization complete."
sleep 10
#####################
# NVT 等容平衡 #
#####################
echo "Starting constant volume equilibration..."
cd ../
mkdir NVT
cd NVT
gmx grompp -f $MDP/NVT/nvt$WINDOW.mdp -c ../EM_2/min$WINDOW.gro -p $FREE_ENERGY/Methane/topol.top -n $FREE_ENERGY/Methane/index.ndx -o nvt$WINDOW.tpr -maxwarn 1
gmx mdrun -nt 2 -v -deffnm nvt$WINDOW -pf pullf-nvt$WINDOW -px pullx-nvt$WINDOW
echo "Constant volume equilibration complete."
sleep 10
#####################
# NPT 等压平衡 #
#####################
echo "Starting constant pressure equilibration..."
cd ../
mkdir NPT
cd NPT
gmx grompp -f $MDP/NPT/npt$WINDOW.mdp -c ../NVT/nvt$WINDOW.gro -p $FREE_ENERGY/Methane/topol.top -n $FREE_ENERGY/Methane/index.ndx -t ../NVT/nvt$WINDOW.cpt -o npt$WINDOW.tpr -maxwarn 1
gmx mdrun -nt 2 -v -deffnm npt$WINDOW -pf pullf-npt$WINDOW -px pullx-npt$WINDOW
echo "Constant pressure equilibration complete."
sleep 10
#################
# 成品模拟 MD #
#################
echo "Starting production MD simulation..."
cd ../
mkdir Production_MD
cd Production_MD
gmx grompp -f $MDP/PRD/prd$WINDOW.mdp -c ../NPT/npt$WINDOW.gro -p $FREE_ENERGY/Methane/topol.top -n $FREE_ENERGY/Methane/index.ndx -t ../NPT/npt$WINDOW.cpt -o prd$WINDOW.tpr -maxwarn 1
gmx mdrun -nt 2 -v -deffnm prd$WINDOW -pf pullf-prd$WINDOW -px pullx-prd$WINDOW
echo "Production MD complete."
# End
echo "Ending. Job completed for Window= $WINDOW"
2.对于以上脚本你仅需要修改：FREE_ENERGY=/home/eric/free_energy/04.US/configura 中的“/home/eric/free_energy/04.US”部分。这一部分为你配置文件configura所在根目录，每个人都是不一样的。
3.以上只是一个独立窗口的完整的模拟脚本，而我们需要模拟27次，即需要下列脚本（write_sh.sh ）批量生成27个job.sh脚本：
#!/bin/bash
#对于模拟, 我们将使用shell脚本替换.job文件中的Wnum关键字
# 打开一个通用的job.sh文件，并替换Wnum的值
#26个窗口, 即Wnum有26个值，从0到26
for ((i=0;i<27;i++))
do
sed 's/Wnum/'$i'/g' $1.sh >$1$i.sh
done
- 准备工作全部做完后如下开始进行27个窗口独立的模拟，使用一个简单脚本（ batche.sh）批量自动化这一过程如下：
#!/bin/bash
#批量执行job脚本,这里0时已经被测式执行来。从1开始
#给新生成的26个Job脚本赋予可执行的权限
chmod +x job*.sh
for num in {
0..26}
do
./job$num.sh
done
5.开始模拟
#1.赋予脚本执行权限
$ chmod +x batche.sh
#2.开始全部的模拟（共27,这里只是用单机电脑来跑的，如果使用集群可以27个任务并行运行，出结果更快）本机4核8线程差不多三天完成。
$ ./batche.sh
分析
1.正常模拟结束将会产生27个Window_0 到Window_27文件夹，该文件夹内容如下，以Window_0为例：
├── Window_0
│ ├── EM_1
│ │ ├── min0.gro
│ │ └── ....
│ ├── EM_2
│ │ ├── min0.gro
│ │ └── ....
│ ├── NVT
│ │ ├──nvt0.gro
│ │ └── .....
│ ├── NPT
│ │ ├── npt0.gro
│ │ └── .....
│ └── PRD
│ ├── prd0.gro
│ ├── prd.tpr ;分析所需文件
│ ├── pullf-prd0.xvg ;分析所需文件
│ └── ....
2.收集26个窗口，即27份独立的模拟成品模拟文件，使用一个简单脚本(collect_dat.sh)将27Window文件夹下的数据收集起来如下：
#!/bin/bash
#创建保存分析需要用的prd0.tpr 到prd26.tpr 和 pullf-prd0.xvg 到pullf-prd26.xvg文件的文件夹 analysis_us
mkdir analysis_us
cd analysis_us
for num in {
0..26}
do
cp ../Window_$num/Production_MD/prd$num.tpr ./
cp ../Window_$num/Production_MD/pullf-prd$num.xvg ./
done
3.我们将使用gmx wham来获取PMF. 该程序接受一个包含.tpr文件列表的文件和另一个包含.xvg文件列表的文件, 后者中包含了输出的牵引力.要创建这两个文件列表的文件, 在analysis_us文件夹下执行以下命令：
$ ls prd.*.tpr > tpr.dat
$ ls pullf-prd.*.xvg > pullf.dat
4.使用GROMACS实现的加权直方图分析(WHAM)来重建PMF：
$ gmx wham -it tpr.dat -f pullf.dat
运行gmx wham之后, 平均力势的计算结果保存在名为profile.xvg的文件中. 绘图应该如下所示:
5.图像修正
我们希望在较远的距离处相互作用为零. 但是, 由于我们使用了三维的偏置势, 需要进行校正. 想象一下将其中一个甲烷作为参考点. 另一个甲烷可以在该点周围距离 r处进行采样, 覆盖了半径为 rr的某个球体的表面. 这为我们的采样增加了额外的构型空间, 降低了熵. 需要消除这些额外熵对PMF的贡献. 回想一下, 等温等压系综中的吉布斯自由能是 −kTln(W), 其中 W 为分配函数. 对我们的情况, 甲烷绕着球体表面采样时, w 与该球体的表面积成正比. 因此, 需要增加 2kTln® 的校正. 此外, 我们需要将图形向上移动, 使其尾部变为零. 对我的结果需要增加的数值大约为7, 但你需要增加的数值可能不同. 要在gnuplot中绘制此图, 请在gnuplot终端中执行以下操作:
#1.在终端启动gnuplot
$ gnuplot
#2.进行修正:输入>后面的代码即可
gnuplot> plot 'profile.xvg' u 1:($2+2*8.314e-3*298.15*log($1)+77) w l
您的PMF现在应该如下所示:
总结
与教程1中的PMF相比, 可以看到它们几乎相同:
区别在于, 使用教程3的直接方法, 采样不可能接近窗口采样. 两个甲烷不会自然地相互靠近, 这就是为什么我们必须增加伞形势以维持它们之间距离的原因。
另一个输出是histo.xvg, 它有助于确定窗口之间是否有足够的重叠. 以下是我模拟的每个直方图的绘图:
很显然, 我们的窗口已经足够重叠了. 否则的话, 我们可能必须选择较小的窗口, 或选择缺少模拟的特定点进行模拟.
3.3 FEP:计算甲烷在水中的溶剂化自由能
该教程在网上三个类似的教程，笔者在这里先简单例举出来，然后再说明这些教程的异同，和参考这些教程哪些内容
| 参考教程链接 参考内容 | ||
|---|---|---|
| 教程1：gromacs 官方教程6： 自由能计算流程，分析方法(bar) | ||
| 教程2：译：自由能计算: 水中的甲烷(GMX 5.0) 与教程1其他都相同，算升级版，多分析了两组数据 | ||
| 教程3：译：Barnett GROMACS教程4: 甲烷的溶剂化自由能 使用python脚本对自由能计算方法自动化分析 | ||
教程3和教程（1和2）不同之处有四点：
一是前者使用水模型同3.1和3.2一样是TIP4P,而后者是TIP3P；
二是前者在描述自由能的变化时考虑了范德华和静电相互作用，而后者只考虑范德华相互作用；
三是前者使用分析自由能的方法不同，前者是TI,后者是使用了gromacs自带的BAR后处理方法;
四是前者只进行了15次独立模拟，后者进行21次独立模拟。
根据FEP和TI计算原理可知，使用分子模拟计算自由能时都是通过多个 λ \lambda λ值控制不同相互作用类型的进度来处理自由能计算，因此这样就可以分别控制库仑、LJ和约束变换等自由能计算，所以无论使用FEP和TI那一种计算方法，在gromacs它们使用的都是同一套自由能计算代码，只是对自由能的后处理分析不同而已。
本教程主要基于教程（1和2）模拟设置进行模拟，，将指导用户计算一个简单的自由能变化, 中性甲烷与一盒子水之间范德华相互作用的去耦合(即移除)，并使用gromacs自带工具和python分析工具（alchemical_analysis）中基于FEP计算自由能的8种后处理方法来简要分析这一过程中自由能的变化，会涉及相关工具的使用和安装。
体系从状态A转变到状态B过程的自由能变化 Δ G A B \Delta G_{AB} ΔGAB是耦合参数 λ \lambda λ的函数, 它代表了状态A和B之间转变的程度, 也就是哈密顿量的微扰程度以及体系已经变化的程度. 利用不同 λ \lambda λ 对应的模拟可以绘制 ∂ H / ∂ λ \partial H/\partial \lambda ∂H/∂λ曲线, 并由此得到 Δ G A B \Delta G_{AB} ΔGAB自由能计算的第一步是确定使用多少个点描述从状态A( λ \lambda λ= 0)到状态B( λ = 1 \lambda = 1 λ=1)的转变, 其目的在于收集足够多的数据, 完全充分地对相空间进行取样, 并得到合理的 ∂ H / ∂ λ \partial H/\partial \lambda ∂H/∂λ 曲线:
如上图所示，状态A是在水中完全相互作用的甲烷. 状态B将是一个甲烷被隐去, 表示它与水的所有范德华相互作用和库仑相互作用（本例只考虑范德华力）都被关闭; 甲烷可以与自身相互作用, 水可以与自身相互作用, 但甲烷和水不会相互作用. 就好像我们将甲烷从水中取出并置于真空中, 而在其他地方我们仍然拥有甲烷所在的水. 我们总共有21个状态。
设置
1.准备甲烷的拓扑文件:注意此模拟使用的TIP3P水模型
; Topology for methane in TIP3P
#include "oplsaa.ff/forcefield.itp"
[ moleculetype ]
; Name nrexcl
Methane 3
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
1 opls_138 1 ALAB CB 1 0.000 12.011
2 opls_140 1 ALAB HB1 2 0.000 1.008
3 opls_140 1 ALAB HB2 3 0.000 1.008
4 opls_140 1 ALAB HB3 4 0.000 1.008
5 opls_140 1 ALAB HB4 5 0.000 1.008
[ bonds ]
; ai aj funct c0 c1 c2 c3
1 2 1
1 3 1
1 4 1
1 5 1
[ angles ]
; ai aj ak funct c0 c1 c2 c3
2 1 3 1
2 1 4 1
2 1 5 1
3 1 4 1
3 1 5 1
4 1 5 1
; water topology
#include "oplsaa.ff/tip3p.itp"
[ system ]
; Name
Methane in water
[ molecules ]
; Compound #mols
Methane 1
SOL 596
你会注意到, 所有的电荷都设置为零. 这种设置有一个实际的原因, 正常情况下, 在关闭范德华相互作用之前, 会先关闭溶质与水之间的电荷相互作用. 如果仅仅保留电荷相互作用而关闭了Lennard-Jones项, 正负电荷相互接近的距离可能无限小, 导致体系非常不稳定, 并可能崩溃. 本教程的步骤基本假定在此之前已经正确地关闭了电荷项, 我们将只关闭溶质与溶剂分子之间的范德华相互作用.这一点下面在设置自由能计算参数（couple-lambda0 ）会提到。
2.参数文件：
同教程3.2.US类似，只是我们添加的是一个计算自由能部分而前者是伞形拉伸的部分, 以便慢慢关闭甲烷与水的相互作用。
另外, 我们还需要为每个 λ \lambda λ状态的准备独立模拟所必备的基本配置文件. 我们需要21个最小化, 21个预平衡,21个成品模拟等的参数文件. 但是我们将使用脚本来简单地更新参数文件模板中的适当值（init_lambda_state ）, 因此实际上只需要为模拟过程的每个部分提供一个模板。
以em_steep.mdp模板为例：（前部分为一般模拟所学控制参数，请重点关注自由能计算部分控制参数！！！）
; ;{
预处理, 使用 C++ 语法
;===============================================================================
include = ; 含引用文件的目录, 拓扑文件中可引用其中的文件, 可多项
; 例: -I/home/joe/joy -I/home/tom/tim -IC:/GMX/top
define = -DFLEXIBLE
; 预定义, 默认无, 可多项, 区分大小写
; -DPOSRES: 使用位置限制文件进行位置限制动力学模拟
; -DFLEXIBLE: 启用柔性水, steep效果更好, cg, l-bfgs或简正分析须开启
;}==============================================================================
;{
能量最小化运行控制
;===============================================================================
integrator = l-bfgs ;能量最小化方法, steep:最陡下降; cg:共轭梯度; lbfgs:二次收敛
nsteps = 5000 ; 最大积分步数, 默认0, -1:无限制
emtol = 100 ; 力的最大容差(kJ/mol-nm), 对壳层值应小于1
emstep = 0.01 ; 初始步长(nm)
nbfgscorr = 10 ; L-BFGS最小化的校正步数
;}==============================================================================
;{
壳层动力学
;===============================================================================
niter = 20 ; 优化壳层位置和柔性约束的最大迭代次数
;}==============================================================================
;{
输出控制
;===============================================================================
nstlog = 1 ; 日志文件输出频率
nstenergy = 1 ; 能量文件输出频率
;}==============================================================================
;{
邻区搜索
;===============================================================================
cutoff-scheme = verlet ; 截断方式, Verlet:粒子截断; Group:电荷组
nstlist = 1 ; 邻区列表更新频率, 0:真空模拟; -1:自动
ns_type = grid ; 邻区搜索算法, Grid:较快; Simple:仅与Group联用
rlist = 1.2 ; 邻区列表截断距离(nm)
rcoulomb = 1.2 ; 静电截断半径, 不超过最小盒子边长一半
rvdw = 1.2 ; 范德华截断半径
;}==============================================================================
;{
周期性
;===============================================================================
pbc = XYZ ; 周期性边界条件, XYZ; XY; No:忽略盒子, 截断与nstlist置零
;}==============================================================================
;{
静电与范德华
;===============
最低0.47元/天 解锁文章
3079
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
