Gromacs 伞形采样

GROMACS教程：伞形采样|Jerkwin参考教程如下，准备计算蛋白质结合的自由能。

GROMACS教程：伞形采样|Jerkwin

新建文件夹 opt/work5 做为计算文件夹

第一步: 准备拓扑文件

下载2BEG文件，放在work5文件夹，重命名为input.pdb

pymol里看起来是这样的，Pymol可以官网免费下载教育版。

先改文件夹权限，方便后面不会permission denied

sudo chmod 777 -R /opt/work5

使用gmx pdb2gmx处理结构:

gmx pdb2gmx -f input.pdb -ignh -ter -o complex.gro

选水模型选1

选N端选2

下一个选 COO-

会选很多次，选错会报错。

现在的文件

在topol_Protein_chain_B.itp文件的末尾加上下面几行:

#ifdef POSRES_B
#include "posre_Protein_chain_B.itp"
#endif

看一下complex文件

第二步: 定义单元晶胞

gmx editconf -f complex.gro -o newbox.gro -center 3.280 2.181 2.4775 -box 6.560 4.362 12

第三步: 添加溶剂和离子

gmx solvate -cp newbox.gro -cs spc216.gro -o solv.gro -p topol.top

加完溶液的样子

添加离子使用ions.mdp文件

; ions.mdp - used as input into grompp to generate ions.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep     ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0    ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01      ; Energy step size
nsteps      = 50000     ; Maximum number of (minimization) steps to perform; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist     = 1         ; Frequency to update the neighbor list and long range forces
ns_type     = grid      ; Method to determine neighbor list (simple, grid)
rlist       = 1.4       ; Cut-off for making neighbor list (short range forces)
coulombtype = PME       ; Treatment of long range electrostatic interactions
rcoulomb    = 1.4       ; Short-range electrostatic cut-off
rvdw        = 1.4       ; Short-range Van der Waals cut-off
pbc         = xyz       ; Periodic Boundary Conditions

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral -conc 0.1

第一步有一个warning 按提示忽略，新增的文件：

第四步: 能量最小化与平衡

gmx grompp -f minim.mdp -c solv_ions.gro -p topol.top -o em.tpr
gmx mdrun -v -deffnm em
gmx grompp -f npt.mdp -c em.gro -p topol.top -o npt.tpr
gmx mdrun -deffnm npt

2021年11月20日16:13:36

2021年11月20日16:51:31

这一步38min

把gro文件拖入pymol，再把trr文件拖入，可以看到能量最小化过程。

第五步: 生成构型

gmx make_ndx -f npt.gro
(>为make_ndx的提示符)
> r 1-27
> name 19 Chain_A
> r 28-54
> name 20 Chain_B
> q

gmx grompp -f md_pull.mdp -c npt.gro -p topol.top -n index.ndx -t npt.cpt -o pull.tpr
gmx mdrun -s pull.tpr

最后一部需要时间较久，用北鲲云计算平台，只需要提交pull.tpr文件

7min

把xtc 和gro文件导入pymol可以看到牵引流程

gmx trjconv -s pull.tpr -f traj.xtc -o conf.gro -sep

选group 0

产生了501个gro文件

使用插件计算每个文件质心距离，插件保存在Desktop\软件\插件

perl distances.pl

生成文件

第六步: 伞形采样模拟

0	-0.012
254	-0.107
279	-0.01
306	0.109
354	0.214
411	0.308
417	0.41
457	0.492

选择这些序列的文件进行采样。大概按间隔0.1选取，这里得到的参数和原教程差别很大，推荐的脚本也不是很好用，选择手动输入，

gmx grompp -f npt_umbrella.mdp -c conf0.gro -p topol.top -n index.ndx -o npt0.tpr -r conf0.gro -maxwarn 5

gmx grompp -f npt_umbrella.mdp -c conf254.gro -p topol.top -n index.ndx -o npt1.tpr -r conf254.gro -maxwarn 5

gmx grompp -f npt_umbrella.mdp -c conf306.gro -p topol.top -n index.ndx -o npt2.tpr -r conf306.gro -maxwarn 5

gmx grompp -f npt_umbrella.mdp -c conf354.gro -p topol.top -n index.ndx -o npt3.tpr -r conf354.gro -maxwarn 5

gmx grompp -f npt_umbrella.mdp -c conf411.gro -p topol.top -n index.ndx -o npt4.tpr -r conf411.gro -maxwarn 5

gmx grompp -f npt_umbrella.mdp -c conf417.gro -p topol.top -n index.ndx -o npt5.tpr -r conf417.gro -maxwarn 5

gmx grompp -f npt_umbrella.mdp -c conf457.gro -p topol.top -n index.ndx -o npt6.tpr -r conf457.gro -maxwarn 5

生成七个npt文件，传入云计算。

很久很久以后

npt模拟后的gro：直接从前面的扭曲成后面这样。有点怪

第六个：不知道为啥pymol导出的gif有黑条。

对现有的npt1、3、4、5、6线进行下一步

gmx grompp -f md_umbrella.mdp -c npt1.gro -p topol.top -n index.ndx -o umbrella1.tpr -r npt1.gro -maxwarn 3

gmx grompp -f md_umbrella.mdp -c npt3.gro -p topol.top -n index.ndx -o umbrella3.tpr -r npt3.gro -maxwarn 3

gmx grompp -f md_umbrella.mdp -c npt4.gro -p topol.top -n index.ndx -o umbrella4.tpr -r npt4.gro -maxwarn 3

gmx grompp -f md_umbrella.mdp -c npt5.gro -p topol.top -n index.ndx -o umbrella5.tpr -r npt5.gro -maxwarn 3

gmx grompp -f md_umbrella.mdp -c npt6.gro -p topol.top -n index.ndx -o umbrella6.tpr -r npt6.gro -maxwarn 3

产生了很多tpr文件

0和2 也好了

gmx grompp -f md_umbrella.mdp -c npt0.gro -p topol.top -n index.ndx -o umbrella0.tpr -r npt0.gro -maxwarn 3

gmx grompp -f md_umbrella.mdp -c npt2.gro -p topol.top -n index.ndx -o umbrella2.tpr -r npt2.gro -maxwarn 3

北鲲云计算步骤：

点Gromacs 然后点右边的提交作业

然后选可视化模板，点本地上传，传一个tpr文件就可以了，然后点下一步，选择合适的核数，然后点两次提交就可以等着下载out文件了

结束后点击作业后面的查看（小眼睛图标），然后下载需要的文件

一次只能上传一个tpr文件，优选16核的海牛。

把输出文件复制到工作文件夹

第七步: 数据分析

新建文件pullf-files.dat和pullx-files.dat

文件内容：

umbrella0_pullf.xvg
umbrella1_pullf.xvg
umbrella2_pullf.xvg
umbrella3_pullf.xvg
umbrella4_pullf.xvg
umbrella5_pullf.xvg
umbrella6_pullf.xvg
第二个类似

运行命令：

gmx wham -it tpr-files.dat -if pullf-files.dat -o -hist -unit kCal

报错

Fatal error:
. Should be tpr, xvg, or pdo.0.tpr

手输dat文件的问题，使用下面的命令生成dat文件

ls umbrella*_pullx.xvg > pullx.dat

告别手输，都用这个命令就成功了

最后的图，和教程差别蛮大

能量曲线也没有收敛，应该是取的区间太少。刚开始的distance也不一样

重新按照英文教程走一遍，数据基本正确，

Umbrella Sampling (mdtutorials.com)http://www.mdtutorials.com/gmx/umbrella/07_analysis.html

取样很少，但看得出来收敛的趋势，而且实践过程几乎没有报warning。