平均力势----以石墨烯-甘氨酸体系为例

阅读笔记。

原始参考文献

标题:Adsorption of amino acids on graphene: assessment of current force fields
DOI:10.1039/C8SM02621A

流程

整体而言内容比较丰富,但是归纳起来也就一句话:研究石墨烯与各种氨基酸的相互作用

设置

初始体系构建有一定小窍门,对整个体系进行分段构建比较实用。

初始体系

  • 依次通过石墨烯构建,氨基酸插入,溶剂化,设置真空区域就能完成初始体系的构建。对于石墨烯片而言,我前面的博文也提及过要注意周期性完整的问题,具体看链接
  • 首先通过博文提及的石墨烯在线网站创建2.45950 x 2.13000 nm^2石墨烯单元(注意这里的长和宽尺寸,即盒子长和宽对于构建无限大的石墨片尤为重要)。
  • 然后直接把甘氨酸(经过文献中所说的封端处理)的gro文件坐标复制到石墨烯的gro文件的后面(盒子尺寸之前)并更改整个文件最前面的原子数使之对应。同时改变盒子z方向长度为8.0,即最后的第三个数字。
  • 读入vmd软件,通过快捷键6,鼠标点击并拖动甘氨酸残基,不要使它与石墨烯的结构重叠并且具有一定距离。保存为complex.gro
  • 通过gmx solvate -cp complex.gro -cs spc216.gro -o sol.gro命令加满水分子。
  • vmd读入上面的sol.gro文件,然后点击Extensions -> TKConsole调出TK终端,然后按照该博文所讲的方法选取盒子z方向的上下一定区域(这里我分别选择z<10z>70),并将其移除,然后将剩余部分保存为lef.pdb,作为gmx模拟的初始结构文件。例如下图:

拓扑文件

采用的amber99sb-ildn力场,拓扑文件:甘氨酸通过acpype,石墨烯通过gmx x2top -pbc,所有石墨烯原子均视为不带电,并进行fix(.mdp中的冻结参数)。然后进行拓扑文件整合,前面博文已提及过,不在详述。

正式模拟

分两步走:

  • 第一步是正常模拟,使体系达到平衡状态,得到最终石墨烯与氨基酸的位置结构文件。
  • 第二步是伞形采样。首先通过拉伸代码得到一系列构象,然后选取一部分构象进行伞形采样模拟即可(前文有具体的,此处略)。
  • 注意:根据文献描述,拉伸步骤分两步骤完成,当反应坐标为0.4 nm<=$\xi$<=0.8 nm时,采用8000 $kJ mol^{-1} nm^{-2}$力常数,窗口间距0.05 nm;0.9 nm<=$\xi$<=2.0 nm时,采用4000 $kJ mol^{-1} nm^{-2}$力常数,窗口间距0.1 nm;并进行五次独立采样统计误差。

结果

文献中的结果分析很丰富,我也没完全弄明白。
进行了其中一个甘氨酸的模拟,绘制的PMF曲线的结果对比如下:

看上去相似,与文献中的自由能相差很小,文献中通过计算也表明体系的大小并不影响自由能的计算。估计还是模拟时间窗口不够导致的。无非就是计算资源。
另外,文献中用的W$\xi [k_BT]$表示的自由能,在温度为298 K时与kJ/mol的转换关系为:

1 $k_BT$ = 2.476 $kJ mol^{-1}$

补充

补充一点静电和范德华相互作用能的计算步骤。详细的方法在这里提及过了,也考虑了使用GROMACS计算分子间相互作用一文中的介绍。这里我采用了第一种(用定义)和第三种(用能量组,使用cut-off)分别计算了一下平衡态下石墨烯和甘氨酸之间的相互作用能的大小,发现这两种方式基本一致。下面也简化了一个计算能量的小脚本以供参考(屏幕最终红色字体输出为轨迹平均能量,生成的t_vs_energy.dat中包括了能量随时间的变化。测试了2019.6没问题):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
#!/bin/bash

# 定义组,不要使用组名代替组编号,因为可能ndx中的组名有重复导致计算失败
xtc="md.xtc" # trajectory filename
tpr="md.tpr" # tpr filename
ndx="index.ndx" # ndx filename
g1=1 # protein group number
g2=13 # ligand group number
plus=22 # complex group number
begin=30000 # starting time (ps)
end=50000 # end time (ps)
DispCorr="N" # 'N' is not turn on dispersion correction(See Your MDP file), 'Y' is turn on


# 计算单组和两组和的静电与范德华相互作用
for i in ${g1} ${g2} ${plus}
do
echo -e "${i}" |gmx trjconv -f $xtc -s $tpr -n $ndx -b ${begin} -e ${end} -o inter_${i}.xtc
echo -e "${i}" |gmx convert-tpr -s $tpr -n $ndx -o new_${i}.tpr
gmx mdrun -rerun inter_${i}.xtc -s new_${i}.tpr -e inter_${i}.edr -g inter_${i}.log

if [[ ${DispCorr} == "N" ]]; then
echo -e "LJ-(SR) \n Coulomb-(SR) \n Coul.-recip." |gmx energy -f inter_${i}.edr -o inter_${i}.xvg |tail -n 3 | awk '{print $3}' > energy_${i}.dat
elif [[ ${DispCorr} == "Y" ]]; then
echo -e "LJ-(SR) \n Disper.-corr. \n Coulomb-(SR) \n Coul.-recip." |gmx energy -f inter_${i}.edr -o inter_${i}.xvg |tail -n 4 | awk '{print $3}' > energy_${i}.dat
fi
done

# 相减得到平均相互作用能
paste energy_${g1}.dat energy_${g2}.dat energy_${plus}.dat > all.dat
if [[ ${DispCorr} == "N" ]]; then
VDW_energy=`awk 'NR==1 {print $3-$2-$1}' all.dat`
Coulomb_energy=`awk 'NR>=2 {C_plus+=$3; C_g1+=$1; C_g2+=$2} END{print (C_plus-C_g1-C_g2)/2.0}' all.dat`
elif [[ ${DispCorr} == "Y" ]]; then
VDW_energy=`awk 'NR<=2 {V_plus+=$3; V_g1+=$1; V_g2+=$2} END {print V_plus-V_g1-V_g2}' all.dat`
Coulomb_energy=`awk 'NR>=3 {C_plus+=$3; C_g1+=$1; C_g2+=$2} END{print (C_plus-C_g1-C_g2)/2.0}' all.dat`
fi
Total_energy=`awk 'BEGIN {VDW_energy = "'${VDW_energy}'"+0; Coulomb_energy = "'${Coulomb_energy}'"+0; print VDW_energy+Coulomb_energy}'`

# 得到随时间模拟时间变化的相互作用能
if [[ ${DispCorr} == "N" ]]; then
for i in ${g1} ${g2} ${plus}
do
awk 'BEGIN{RS="\r?\n"} /^[^@#]/{print}' inter_${i}.xvg > temp_${i}.dat
done
paste temp_${g1}.dat temp_${g2}.dat temp_${plus}.dat > temp_all.dat
awk 'BEGIN{RS="\r?\n"; print "#Time\tVDW\tElec\tTotal"} {print $1, $10-$6-$2, ($11+$12-$7-$8-$3-$4)/2.0, $10-$6-$2+($11+$12-$7-$8-$3-$4)/2.0}' temp_all.dat > t_vs_energy.dat
elif [[ ${DispCorr} == "Y" ]]; then
for i in ${g1} ${g2} ${plus}
do
awk 'BEGIN{RS="\r?\n"} /^[^@#]/{print}' inter_${i}.xvg > temp_${i}.dat
done
paste temp_${g1}.dat temp_${g2}.dat temp_${plus}.dat > temp_all.dat
awk 'BEGIN{RS="\r?\n"; print "#Time\tVDW\tElec\tTotal"} {print $1, $12+$13-$7-$8-$2-$3, ($14+$15-$9-$10-$4-$5)/2.0, $12+$13-$7-$8-$2-$3+($14+$15-$9-$10-$4-$5)/2.0}' temp_all.dat > t_vs_energy.dat
fi

# 清除中间文件
rm -rf inter_* all.dat temp_*.dat energy_*.dat new_*.tpr *#

# 输出
echo -e "\e[1;31mVDW_energy is ${VDW_energy} kJ/mol\e[0m"
echo -e "\e[1;31mCoulomb_energy is ${Coulomb_energy} kJ/mol\e[0m"
echo -e "\e[1;31mTotal interaction energy is ${Total_energy} kJ/mol\e[0m"

其他资料

  1. http://cmb.bio.uni-goettingen.de/pub/Hub_JCTC2010.pdf
  2. http://sci-hub.tw/https://pubs.acs.org/doi/pdf/10.1021/jp401584u?rand=mc36q5kw
  3. http://stuchebrukhov.ucdavis.edu/utilities/ConversionTables.html