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

阅读笔记。

原始参考文献

标题: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)分别计算了一下平衡态下石墨烯和甘氨酸之间的相互作用能的大小,发现这两种方式基本一致。下面也简化了一个计算能量的小脚本以供参考:

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
#!/bin/bash

# 定义组
g1="GRA"
g2="UNL"
plus="GRA_UNL"
begin=8000
end=10000

# 计算单组和两组和的静电与范德华相互作用
for i in ${g1} ${g2} ${plus}
do
echo -e "${i}" |gmx trjconv -f md.trr -s md.tpr -n index.ndx -b ${begin} -e ${end} -o inter_${i}.xtc
echo -e "${i}" |gmx convert-tpr -s md.tpr -n index.ndx -o new_${i}.tpr
gmx mdrun -rerun inter_${i}.xtc -s new_${i}.tpr -e inter_${i}.edr -g inter_${i}.log
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 > inter_f${i}.txt
awk 'NR==1 {print $3}' inter_f${i}.txt > LJ_SR_${i}.dat
awk 'NR==2 {print $3}' inter_f${i}.txt > Disper.-corr._${i}.dat
awk 'NR==3 {print $3}' inter_f${i}.txt > Coulomb_SR_${i}.dat
awk 'NR==4 {print $3}' inter_f${i}.txt > Coul.-recip._${i}.dat
done

# 相减得到相互作用能
VDW_energy=$(echo "(`cat LJ_SR_${plus}.dat`)+(`cat Disper.-corr._${plus}.dat`)-(`cat LJ_SR_${g1}.dat`)-(`cat Disper.-corr._${g1}.dat`)-(`cat LJ_SR_${g2}.dat`)-(`cat Disper.-corr._${g2}.dat`)" | bc)
Coulomb_energy=$(echo "(`cat Coulomb_SR_${plus}.dat`)+(`cat Coul.-recip._${plus}.dat`)-(`cat Coulomb_SR_${g1}.dat`)-(`cat Coul.-recip._${g1}.dat`)-(`cat Coulomb_SR_${g2}.dat`)-(`cat Coul.-recip._${g2}.dat`)" | bc)
Total_energy=$(echo "${VDW_energy}+${Coulomb_energy}" |bc)

# 清除中间文件
rm -rf inter_* Coul* LJ_SR* Disper* 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