Cellulose nanofibril MD tutorial

Introduction

原教程连接为https://wiki.aalto.fi/display/IMM/Cellulose+nanofibril+MD+tutorial，在此对这篇教程进行翻译，并添加必要的解释和笔记。这篇教程是以纤维素纳米纤维为例进行的指导，但实际上这个方法可以被应用到所有需要PBC的纤维素结构模拟上，比如5条链并排的单层纤维素平面，只需要调整一些参数即可，在后文对定位置会用“😇”标注出来。

小标题是译者自己起的。

Cellulose Structure Building

首先，我们通过 Cellulose-builder 搭建一个纤维素纳米纤维的结构（生成结构的部分在此略去，可以参考这篇教程）。

最终结构如下图所示（图片源自原教程）：纤维二糖、纤维素纳米纤维的结构参数见标注。😇请注意，如果你要用其他结构搭建 top，那么请合理计算出你的纤维素结构的参数 (长，宽，高)（这里的“长”指的是纤维素单链的长度，“宽”和“高”是不言自明的；“长”是必须合理计算出来的，但是“宽”和“高”不必要（算出来也没什么坏处）），这些参数要用于后续设置 box 大小，以便于生成具有PBC条件的纤维素链😇

Modify the Force Field Parameter

在开始模拟之前，我们必须要有纤维素的基本结构单元葡萄糖（虽然理论上是纤维二糖，但是这里力场的最小单元可以是葡萄糖）的键合和非键参数，在本教程中是用的是Amber14 力场，你可以在这里下载它。😇这些参数来自于这篇文献 DOI: 10.3390/
nano12193379，是 Vaibhav Modi 等人基于 GLYCAM06 力场进行量化计算优化后的参数（毕竟 GLYCAM06 只有糖单体的力场，而且用 tleap 生成这么大的体系属实有些麻烦）；如果你想用其他立场的话，也可以去自己找参数，然后添加到相应文件里面（当然如果作者提供了全套的就更好了，就像上面 Amber14 那个）；这里译者再提供一个 CHARMM 的力场参数 DOI: 10.1021/ct900242e（在 Supporting 里），但是由于原作者力场编写格式是 CHARMM 的，所以需要把1-4连接的葡萄糖当成一个非标准残基添加到用于 Gromacs 的 CHARMM36c 里，具体教程请移步此处😇

下载完之后，请解压，会有一个力场 amber14sb_OL15.ff 和两个文件 residuetypes.dat、specbond.dat 。接下来你可以选择：①只在你希望的文件夹里用，那就把解压出来的文件放到那个文件夹里就行。②希望在哪都能用这个力场，把解压的所有文件都放到<path-to-force-field-directory>这个文件夹里，然后：

echo export GMXLIB=”<path-to-force-field-directory>” >> ~/.bashrc
# or
export GMXLIB=”<path-to-force-field-directory>”

这里我们先看一下这个力场与普通的Amber力场有何区别。首先是多了三个与糖相关的力场参数文件：

New-glc.rtp: Residue topology with partial charges and bond descriptions 
New-glc.r2b: Residue name translation database 
New-glc.arn: Atom renaming database

另外 residuetypes.dat、specbond.dat 本来是出现在 Gromacs 安装目录下的文件，解压目录下的 residuetypes.dat 多了下面这些，用以表示残基归属（方便pdb2gmx时自动命名）：

RGB	Cellulose
NGB	Cellulose
IBG	Cellulose

specbond.dat 多了：

RBG     O4      1       IBG     C1      1       1.475   RBG     IBG
IBG		O4		1		IBG		C1		1		1.475   IBG		IBG
IBG     O4      1       NBG     C1      1       1.475   IBG     NBG

此外还多了atomname2type.n2t，这个在后面会详细解释。

Build the System

Preliminary and Secondary Topology Files

在生成 top 之先修改一下 Cellulose-builder 生成的 crystal.pdb ，因为它里面的残基命名是 BGLC 而不是力场里的 IBG（有时候他也会命名成BGLCX，记得自己看一下 pdb）：

sed -i "s/BGLC/IBG /" crystal.pdb 

随后用pdb2gmx生成一个初步的 topology files，该 top [atom]字段将用于最终的top：

gmx_mpi pdb2gmx -f crystal.pdb -o pdb2gmx.gro -p topol-pdb2gmx.top -water tip3p

请忽略这里的 long bond warning，因为pdb2gmx不支持周期性共价键的生成，这也是为什么我们只用到 topol-pdb2gmx.top 的[atom]；后面将用x2top生成剩下的参数。

接着设置box大小：

gmx_mpi editconf -f pdb2gmx.gro -o mfc-20mono.gro -box 7 8 10.38

😇如果 Cellulose-builder 给你生成的体系的糖链没有沿着你想要的方向，或者不在你想要的平面的话，使用下面这个指令进行摆正😇

gmx editconf -f pdb2gmx.gro -o mfc-20mono.gro -princ
# then
gmx editconf -f mfc-20mono.gro  -o mfc-20mono.gro -box 8 10.38 6 -center 4 5.19 1

现在使用x2top而不是pdb2gmx生成一个具有周期性键连参数的拓扑。

gmx x2top -f mfc-20mono.gro -o topol-x2top.top -ff select -pbc -name IBG -noparam -alldih

Final Topology Files

对 topol-x2top.top 中的一些参数进行修改：

sed -i '/\[ dihedrals \]/,/\[ system \]/s/1 *$/9/' topol-x2top.top
# x2top不会自动匹配到相应的angles或者dihedrals函数类型，所以需要修改这一部分以确保和topol-pdb2gmx.top的内容相同

😇如果你用的是先前提到的修改版 CHARMM36c 体系参数，那么需要把angles部分的设置修改一下（见下）。当然，无论你用的是什么力场，最重要的是要检查一下 topol-pdb2gmx.top 和 topol-x2top.top 让他们在 angles或者dihedrals 的函数类型上一致😇

sed -i '/\[ angles \]/,/\[ system \]/s/1 *$/5/' topol-x2top.top

随后将 topol-pdb2gmx.top 与 topol-x2top.top 整合到一起。

sed -n '/Include forcefield parameters/,/\[ bonds \]/p' topol-pdb2gmx.top | sed '/bond/d' >& newtopol.top && sed -n '/\[ bonds \]/,/\[ system \]/p' topol-x2top.top | sed '/\[ system \]/d' >> newtopol.top && sed -n '/Include Position restraint file/,//p' topol-pdb2gmx.top >> newtopol.top

😇至此我们就搭建好了最基础的模拟体系，如果你还希望加点其他东西（比如蛋白或者其他分子的），就像常规的蛋白质-配体体系搭建方法一样，整合 .gro 文件，在 .top 中添加 .itp 文件的引入，可以参考这里😇

Solvate

接下来的操作与常规的Gromacs模拟类似，首先溶剂化体系：

gmx_mpi solvate -cp mfc-20mono.gro -cs spc216.gro -p newtopol.top -o box-sol.gro

😇如果需要的话，接着添加离子中和体系😇

gmx grompp -f ions.mdp -c box-sol.gro  -p newtopol.top -o ions.tpr

gmx genion -s ions.tpr -o solv_ions.gro -p newtopol.top -pname NA -nname CL -neutral -conc 0.15

记得检查一下 newtopol.top 的 [ molecules ] 是不是正确写入了。

Molecular Dynamic Stimulation

Energy minimization

在开始动力学之前，使用能量最小化来松弛系统并消除任何空间冲突。模拟周期性/无限（periodic/infinite）纤维素链参数所需的一个关键参数是在mdp文件中添加“periodic_molecules”。使用此mdp参数文件生成tpr文件：

gmx_mpi grompp -f em.mdp -p newtopol.top -c box-sol.gro -o em

gmx_mpi mdrun -deffnm em

😇你可以用gmx_mpi energy来检查一下这个过程，这种操作比较基础在这里就忽略了😇

Equilibration

我们将对溶剂化的纤维进行两步平衡，首先用NVT系综进行500ps温度平衡，然后进行500ps的恒温恒压NPT平衡。短时间足以平衡小型/中型系统，但我们仍需要检查输出以重新确定平衡时间。NVT所需 .mdp 在此：

gmx_mpi grompp -f nvt.mdp -p newtopol.top -c em.gro -r em.gro -o nvt 

gmx_mpi mdrun -deffnm nvt -dlb yes	# -dlb也可以设置为 auto

在将系统的温度提高到300K左右的稳定范围后，我们现在将使用温浴、压浴相结合，在与实验条件非常相似的环境中对系统进行建模。在第二步平衡中，我们将使用 V-rescale 恒温器、弱耦合 Berendsen 压浴和并对重原子进行位置约束，以运行 500ps 的NPT系综。注意，我们采取了半各向同性压力耦合（pcoupltype = semiisotropic），以分离原纤维沿着链轴（Z轴）和其他两个轴感受压力的方式。压力耦合在X和Y方向上是各向同性的。NPT 所需 .mdp 文件在此：

# .mdp的压浴设置
; Pressure coupling 
pcoupl                  = Berendsen     ;Parrinello-Rahman for production run (NPT) 
pcoupltype              = semiisotropic ; uniform scaling of box vectors in x-y and independent z 
tau_p                   = 2.0           ; time constant, in ps 
ref_p                   = 1.0 1.0       ; reference pressure, in bar 
compressibility         = 4.5e-5 4.5e-5 ; isothermal compressibility of water, bar^-1 
refcoord_scaling        = com

执行 NPT 平衡：

gmx_mpi grompp -f npt.mdp -p newtopol.top -c nvt.gro -r nvt.gro -t nvt.cpt -o npt 

gmx_mpi mdrun -deffnm npt -dlb yes 

正如预期的那样，压力值将大幅波动，但平均值应更接近1巴的参考压力值。Gromacs手册和教程中详细解释了波动背后的原因。

Production run

我们将使用与 NPT 相同的设置进行MD，但不施加任何位置限制，同时使用 Parrionello Rahman 进行更紧密的压力耦合。您可以将 npt.mdp 复制到 md.mdp 并进行以下更改，也可以在此处下载参数文件。

pcoupl  = Parrinello-Rahman 

进行MD模拟：

gmx_mpi grompp -f md.mdp -p newtopol.top -c npt.gro -o topol 

gmx_mpi mdrun -s topol -dlb yes 

随后你就可以对模拟结果进行你想做的分析了。

Happy simulating!