手把手教你用Gromacs完成溶菌酶在水中的动力学模拟

2022年8月9日 浏览：1909 评论：1
随着中国科研的发展，有限的超算中心和自建计算集群的方式已无法满足超算用户计算业务和创新的需要，特别对降低超算总体拥有成本，无需排队获取资源，弹性规模，软件开箱即用和安全运维等核心需求，云上超算是唯一解决之道。

今天，将为大家介绍 如何在北鲲云超算平台上快速完成动力学模拟计算，模拟溶菌酶（lysozyme）在水相中的结构变化 。北鲲云超算平台是前段时 间小伙伴安利的，经过一段时间测试，效果还不错。 只需要一台上网本就可以在短时间内完成动力学模拟所有内容，而且还预装了Gromacs在内的300+应用软件，实现了开箱即用，非常的灵活便捷。接下来，就正式进入今天的主题。
手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图1
水中的溶菌酶
 
01

 结构处理


溶菌酶又称胞壁质酶或N-乙酰胞壁质聚糖水解酶，是一种 能水解细菌中黏 多糖的碱性酶。溶菌酶还可与带负电荷的病毒蛋白直接结合，与DNA、RNA、脱辅基蛋白形成复合体，使病毒失活。这里我们从PDB数据库下载溶菌酶蛋白质三维结构，编号为1AKI，或则利用PyMOL直接从命令行获取三维结构，然后利用PyMOL检查结构，是否缺 失，并且删除除开蛋白质外的杂原子和水分子，并保存为protein.pdb文件。 
fetch 1AKI


02

 将输入文件上传到北鲲云超算平台


打开 北鲲云超算 平台，注册登录后在/home/cloudam/jobs下新建文件夹Lysozyme，将准备好的protein.pdb文件直接拖拽到该文件中，即可完成上传。上传文件和下载文件都可以利用sftp非常方便完成。

除了命令行，北鲲云还有可视化和图形界面作业提交方式。可视化无需代码，比较适合无计算机专业背景的同学，图形界面适 合于单节点小规模计算任务。这两种作业提交方式大家就自行体验吧（北鲲云针对新用户有赠送200元体验金的活动，大家可以关注一下）。
手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图2 北鲲云超算平 台作业提交入口
手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图3
北鲲云超算平台操作终端窗口


03

 生成拓扑文件


在确保蛋白质没有杂原子且结构完整后，将其上传至计算平台，然后我们利用Gromacs开始进行动力学模拟， 在此之前 需要选择所使用的Gromacs版本，在北鲲云超算平台上可以使用的Gromacs如下： 
  

   GROMACS/2019.6-fosscuda GROMACS/2020-foss-2019b 
GROMACS/2021-gromacs-cpu-new
GROMACS/gromacs-5.1.5-gpu
GROMACS/2019.6-intel-2019b    
GROMACS/2020-fosscuda-2019b   
GROMACS/4.5.5-openmpi-1.3.1   
GROMACS/2019.6-intelmpi-cuda
GROMACS/2021-fosscuda-2019b
GROMACS/gromacs-5.1.2

这里我们使用Gromacs2020版进行本次计算，也可自由选择其他版本进行模拟（如果没有你需要的版本，也可找技术支持安装），加载命令如下： 
module add Gromacs/2020-fosscuda-2019b

可以通过通过键入gmx或者gmx_mpi命令查看是否加载成功，如果出现类似如下信息则表明加载成功，可以进行后续模拟： 
  

   [cloudam@master jobs]$ gmx_mpi
   
                        
   :-) GROMACS - gmx_mpi, 
   2021 (-:
   

   
                            GROMACS is written 
   by:
   
     Andrey Alekseenko Emile Apol Rossen Apostolov
   
         Paul Bauer Herman J.C. Berendsen Par Bjelkmar
   
       Christian Blau Viacheslav Bolnykh Kevin Boyd
   
     Aldert van Buuren Rudi van Drunen Anton Feenstra
   
    Gilles Gouaillardet Alan Gray Gerrit Groenhof
   
       Anca Hamuraru Vincent Hindriksen M. Eric Irrgang
   
      Aleksei Iupinov Christoph Junghans Joe Jordan
   
    Dimitrios Karkoulis Peter Kasson Jiri Kraus
   
      Carsten Kutzner Per Larsson Justin A. Lemkul
   
       Viveca Lindahl Magnus Lundborg Erik Marklund
   
        Pascal Merz Pieter Meulenhoff Teemu Murtola
   
        Szilard Pall Sander Pronk Roland Schulz
   
       Michael Shirts Alexey Shvetsov Alfons Sijbers
   
       Peter Tieleman Jon Vincent Teemu Virolainen
   
     Christian Wennberg Maarten Wolf Artem Zhmurov
   
                           
   and the project 
   leaders:
   
        Mark Abraham, Berk Hess, Erik Lindahl, 
   and David van der Spoel

pdb2gmx是我们用到的第一个GROMACS模块, 它的作用的是产生三个文件:分子拓扑文件topol.top、位置限制文件posre.itp、后处理结构文件protein_processed.gro,具体命令如下： 
gmx pdb2gmx -f protein.pdb -o protein_processed.gro -water tip3p -ignh -merge all
此时会出现如下选择力场的提示： 
  

   Select the 
   Force 
   Field:
   

   

   From 
   '/public/software/.local/easybuild/software/GROMACS/2021-gromacs-cpu-new/share/gromacs/top':
   
 
   1: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 
   24, 
   1999
   -2012, 
   2003)
   
 
   2: AMBER94 
   force 
   field (Cornell et al., JACS 
   117, 
   5179
   -5197, 
   1995)
   
 
   3: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 
   29, 
   461
   -469, 
   1996)
   
 
   4: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 
   21, 
   1049
   -1074, 
   2000)
   
 
   5: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 
   65, 
   712
   -725, 
   2006)
   
 
   6: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 
   78, 
   1950
   -58, 
   2010)
   
 
   7: AMBERGS 
   force 
   field (Garcia & Sanbonmatsu, PNAS 
   99, 
   2782
   -2787, 
   2002)
   
 
   8: CHARMM27 all-atom 
   force 
   field (CHARM22 plus CMAP 
   for proteins)
   
 
   9: GROMOS96 
   43a1 
   force 
   field
   

   10: GROMOS96 
   43a2 
   force 
   field (improved alkane dihedrals)
   

   11: GROMOS96 
   45a3 
   force 
   field (Schuler JCC 
   2001 
   22 
   1205)
   

   12: GROMOS96 
   53a5 
   force 
   field (JCC 
   2004 vol 
   25 pag 
   1656)
   

   13: GROMOS96 
   53a6 
   force 
   field (JCC 
   2004 vol 
   25 pag 
   1656)
   

   14: GROMOS96 
   54a7 
   force 
   field (Eur. Biophys. J. (
   2011), 
   40,, 
   843
   -856, DOI: 
   10.1007/s00249
   -011
   -0700
   -9)
   

   15: OPLS-AA/L all-atom 
   force 
   field (
   2001 aminoacid dihedrals)

这里选择力场AMBER99来处理蛋白质，键入4回车，检查生成文件应该包含如下内容： 
posre.itp  protein.pdb  protein_processed.gro  topol.top

04

定义单位盒子并填充溶剂 


使用editconf来定义盒子，命令如下 ： 
gmx editconf -fprotein_processed.gro -o pro_newbox.gro -c -d 1.0 -bt cubic其中-c表示将蛋白质放在盒子中心，-d 1.0表示蛋白质距离盒子边缘至少1埃距离，-bt表示盒子形状为立方体。

然后使用solvate来填充水分子，命令如下： 
gmx solvate -cppro_newbox.gro -cs spc216.gro -o pro_solv.gro -p topol.top
-cp指定输入文件，-cs指定溶剂模型文件，这里是通用的三位点溶剂模型spc216.gro

05

添加离子 


完成上述步骤后，已经获得了带有电荷的溶液体系，topol.top文件记录的电荷信息。因此需要添加对应的离子来进行体系中和，至于需要添加多少离子，gromacs就比较方便了，是不需要像Amber一样进行计算的。添加离子需要两步命令，第一步需要一个ion.mdp参数文件，这里我们将所有的mdp参数文件放在MDP文件夹下，命令如下： 
gmx grompp -f../MDP/ions.mdp -c pro_solv.gro -p topol.top -o ions.tpr
上述命令会获得 一个tpr结尾的文件，该文件它提供了我们体系的原子级别的描述。

 第二步进行利用Na离子和Cl离子替换溶剂部分水分子，使得体系电荷被中和。命令如下: 
gmx grompp -f ../MDP/minim.mdp -c pro_solv_ions.gro -p topol.top -o em.tpr

溶剂SOL的组编号为13,因此选择13然后回车，如下所示：
  

   Reading file ions.tpr, VERSION 2021 (single precision)
   
Reading file ions.tpr, VERSION 2021 (single precision)
   
Will try to add 0 NA ions and 8 CL ions.
   

   Select a continuous 
   group 
   of solvent molecules
   

   Group     
   0 ( 
   System) has 
   33892 elements
   

   Group     
   1 ( Protein) has 
   1960 elements
   

   Group     
   2 ( Protein-H) has 
   1001 elements
   

   Group     
   3 ( C-alpha) has 
   129 elements
   

   Group     
   4 ( Backbone) has 
   387 elements
   

   Group     
   5 ( MainChain) has 
   517 elements
   

   Group     
   6 ( MainChain+Cb) has 
   634 elements
   

   Group     
   7 ( MainChain+H) has 
   646 elements
   

   Group     
   8 ( SideChain) has 
   1314 elements
   

   Group     
   9 ( SideChain-H) has 
   484 elements
   

   Group    
   10 ( Prot-Masses) has 
   1960 elements
   

   Group    
   11 ( non-Protein) has 
   31932 elements
   

   Group    
   12 ( Water) has 
   31932 elements
   

   Group    
   13 ( SOL) has 
   31932 elements
   

   Group    
   14 ( non-Water) has 
   1960 elements
   

   Select a 
   group: 
   13

参数文件内容：

微信截图_20220809170450.png

手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图5手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图6手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图7

手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图8

手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图9
手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图10
VMD检查体系中蛋白质和离子分布

06

体系能量最小化 


命令如下，本命令同样需要mdp参数文件，名为minim.mdp， 能量最小化过程与添加离子过程相似，利用grompp将 拓扑和模拟参数写入一个二进制的输入文件中(.tpr), 然后使用GROMACS MD引擎的mdrun模块来进行能量最小化 
gmx grompp -f ../MDP/minim.mdp -c pro_solv_ions.gro -p topol.top -o em.tpr

MDP参数文件内容：

2.png


然后键入如下命令进行能量最小化，这一步应当提交到计算节点进行，不要在管理节点直接运行。 
gmx mdrun -v -deffnm em
能量最小化开始，屏幕中出现运行进度提示： 
  

   Steepest Descent
   s:
   
   Tolerance (Fmax) = 
   1.00000
   e+
   03
   
   Number of steps = 
   50000
   
Step= 
   0, Dmax= 
   1.0
   e-
   02 
   nm, Epot= -
   2.94025
   e+
   05 Fmax= 
   1.53483
   e+
   05, atom= 
   19880
   
Step= 
   1, Dmax= 
   1.0
   e-
   02 
   nm, Epot= -
   3.12614
   e+
   05 Fmax= 
   6.09036
   e+
   04, atom= 
   19880
   
Step= 
   2, Dmax= 
   1.2
   e-
   02 
   nm, Epot= -
   3.35310
   e+
   05 Fmax= 
   2.40700
   e+
   04, atom= 
   19880
   
Step= 
   3, Dmax= 
   1.4
   e-
   02 
   nm, Epot= -
   3.61996
   e+
   05 Fmax= 
   1.00755
   e+
   04, atom= 
   24188
   
Step= 
   4, Dmax= 
   1.7
   e-
   02 
   nm, Epot= -
   3.87330
   e+
   05 Fmax= 
   4.66750
   e+
   03, atom= 
   983
   
Step= 
   5, Dmax= 
   2.1
   e-
   02 
   nm, Epot= -
   4.08800
   e+
   05 Fmax= 
   1.30036
   e+
   04, atom= 
   547
   
Step= 
   6, Dmax= 
   2.5
   e-
   02 
   nm, Epot= -
   4.13336
   e+
   05 Fmax= 
   2.07112
   e+
   04, atom= 
   547
   
Step= 
   7, Dmax= 
   3.0
   e-
   02 
   nm, Epot= -
   4.17185
   e+
   05 Fmax= 
   1.90699
   e+
   04, atom= 
   547
   
Step= 
   8, Dmax= 
   3.6
   e-
   02 
   nm, Epot= -
   4.18940
   e+
   05 Fmax= 
   2.96782
   e+
   04, atom= 
   547
   
Step= 
   9, Dmax= 
   4.3
   e-
   02 
   nm, Epot= -
   4.21615
   e+
   05 Fmax= 
   2.83858
   e+
   04, atom= 
   547
   
Step= 
   11, Dmax= 
   2.6
   e-
   02 
   nm, Epot= -
   4.26495
   e+
   05 Fmax= 
   6.89298
   e+
   03, atom= 
   547
   
Step= 
   12, Dmax= 
   3.1
   e-
   02 
   nm, Epot= -
   4.27186
   e+
   05 Fmax= 
   3.75623
   e+
   04, atom= 
   1785
   
Step= 
   13, Dmax= 
   3.7
   e-
   02 
   nm, Epot= -
   4.33871
   e+
   05 Fmax= 
   1.81256
   e+
   04, atom= 
   1785
   
Step= 
   15, Dmax= 
   2.2
   e-
   02 
   nm, Epot= -
   4.35813
   e+
   05 Fmax= 
   1.54546
   e+
   04, atom= 
   1785
   
Step= 
   16, Dmax= 
   2.7
   e-
   02 
   nm, Epot= -
   4.37056
   e+
   05 Fmax= 
   2.45185
   e+
   04, atom= 
   1785

利用如下命令检查是否有效的能量最小化： 
echo 10 0 | gmx energy -f em.edr -o potential.xvg
利用xmgrace打 开potential.xvg文件，结果如如下所 示：
手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图12
体系势能随时间变化曲线

07

NVT平衡


NVT系综(粒子数, 体积和温度都是恒定的)系综也被称为等温等容系综或正则系综。NVT平衡过程的需要的时间与体系的构成有关, 但在NVT系综中, 体系的温度应达到预期值并基本保持不变. 如果温度仍然没有稳定, 那就需要更多的时间。通常情况下, 50 ps到100 ps就足够了, 因此在本例中我们进行100 ps的NVT平衡. 根据机器的不同, 运行可能需要一段时间(在双核的MacBook上刚刚超过一小时，而在北鲲云上几分钟即可完成)。命令如下： 
  

   gmx grompp -
   f ../MDP/nvt.mdp -
   c 
   em.gro -r 
   em.gro -
   p topol.top -
   o nvt.tpr
  

  

   gmx mdrun -deffnm nvt

MDP参数文件内容如下： 
  

   title                   = OPLS Lysozyme NVT equilibration
   

   define                  = -DPOSRES ; position restrain the protein
   

   ; Run parameters
   

   integrator              = md ; leap-frog integrator
   

   nsteps                  = 
   50000     ; 
   2 * 
   50000 = 
   100 ps
   

   dt                      = 
   0.002     ; 
   2 fs
   

   ; Output control
   

   nstxout                 = 
   500       ; save coordinates every 
   1.0 ps
   

   nstvout                 = 
   500       ; save velocities every 
   1.0 ps
   

   nstenergy               = 
   500       ; save energies every 
   1.0 ps
   

   nstlog                  = 
   500       ; update log file every 
   1.0 ps
   

   ; Bond parameters
   

   continuation            = 
   no        ; first dynamics run
   

   constraint_algorithm    = lincs ; holonomic constraints
   

   constraints             = h-bonds ; bonds involving H are constrained
   

   lincs_iter              = 
   1         ; accuracy of LINCS
   

   lincs_order             = 
   4         ; also related to accuracy
   

   ; Nonbonded settings
   

   cutoff-scheme           = Verlet ; Buffered nei
   ghor searching
   

   ns_type                 = grid ; search nei
   ghoring grid cells
   

   nstlist                 = 
   10        ; 
   20 fs, largely irrelevant with Verlet
   

   rcoulomb                = 
   1.0       ; short-range electrostatic cut
   off (in nm)
   

   rvdw                    = 
   1.0       ; short-range van der Waals cut
   off (in nm)
   

   DispCorr                = EnerPres ; account for cut-
   off vdW scheme
   

   ; Electrostatics
   

   coulombtype             = PME ; Particle Mesh Ewald for long-range electrostatics
   

   pme_order               = 
   4         ; cubic interpolation
   

   fourierspacing          = 
   0.16      ; grid spacing for FFT
   

   ; Temperature coupling is on
   

   tcoupl                  = V-rescale ; modified Berendsen thermostat
   

   tc-grps                 = Protein Non-Protein ; two coupling groups - more accurate
   

   tau_t                   = 
   0.1     
   0.1           ; time constant, in ps
   

   ref_t                   = 
   300     
   300           ; reference temperature, 
   one for each group, in K
   

   ; Pressure coupling is off
   

   pcoupl                  = 
   no        ; 
   no pressure coupling in NVT
   

   ; Periodic boundary conditions
   

   pbc                     = xyz ; 
   3-D PBC
   

   ; Velocity generation
   

   gen_vel                 = 
   yes       ; assign velocities from Maxwell distribution
   

   gen_temp                = 
   300       ; temperature for Maxwell distribution
   

   gen_seed                = -
   1        ; generate a random seed

同样通过使用energy来检查平衡是否完成，命令如下： 
echo 16 0 |gmx energy -f nvt.edr -o temperature.xvg
结果如下所示，温度基本保持在300K，如果波动太厉害可以考虑增加nstes步数来进一步获取平衡体系。
手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图13
体系温度随时间变化曲线

08

NPT平衡


前一步的NVT平衡稳定了体系的温度. 在采集数据之前, 我们还需要稳定体系的压力(因此还包括密度)。压力平衡是在NPT系综下进行的, 其中粒子数, 压力和温度都保持不变。这个系综也被称为等温等压系综, 最接近实验条件。100 ps NPT平衡的命令如下： 
  

   gmx grompp -f ../MDP/npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
   

   
gmx mdrun -deffnm npt

MDP参数内容： 
  

   title                   = OPLS Lysozyme NPT equilibration
   

   define                  = -DPOSRES ; position restrain the protein
   

   ; Run parameters
   

   integrator              = md ; leap-frog integrator
   

   nsteps                  = 
   50000     ; 
   2 * 
   50000 = 
   100 ps
   

   dt                      = 
   0.002     ; 
   2 fs
   

   ; Output control
   

   nstxout                 = 
   500       ; save coordinates every 
   1.0 ps
   

   nstvout                 = 
   500       ; save velocities every 
   1.0 ps
   

   nstenergy               = 
   500       ; save energies every 
   1.0 ps
   

   nstlog                  = 
   500       ; update log file every 
   1.0 ps
   

   ; Bond parameters
   

   continuation            = 
   yes       ; Restarting after NVT
   

   constraint_algorithm    = lincs ; holonomic constraints
   

   constraints             = h-bonds ; bonds involving H are constrained
   

   lincs_iter              = 
   1         ; accuracy of LINCS
   

   lincs_order             = 
   4         ; also related to accuracy
   

   ; Nonbonded settings
   

   cutoff-scheme           = Verlet ; Buffered nei
   ghor searching
   

   ns_type                 = grid ; search nei
   ghoring grid cells
   

   nstlist                 = 
   10        ; 
   20 fs, largely irrelevant with Verlet scheme
   

   rcoulomb                = 
   1.0       ; short-range electrostatic cut
   off (in nm)
   

   rvdw                    = 
   1.0       ; short-range van der Waals cut
   off (in nm)
   

   DispCorr                = EnerPres ; account for cut-
   off vdW scheme
   

   ; Electrostatics
   

   coulombtype             = PME ; Particle Mesh Ewald for long-range electrostatics
   

   pme_order               = 
   4         ; cubic interpolation
   

   fourierspacing          = 
   0.16      ; grid spacing for FFT
   

   ; Temperature coupling is on
   

   tcoupl                  = V-rescale ; modified Berendsen thermostat
   

   tc-grps                 = Protein Non-Protein ; two coupling groups - more accurate
   

   tau_t                   = 
   0.1     
   0.1           ; time constant, in ps
   

   ref_t                   = 
   300     
   300           ; reference temperature, 
   one for each group, in K
   

   ; Pressure coupling is on
   

   pcoupl                  = Parrinello-Rahman ; Pressure coupling 
   on in NPT
   

   pcoupltype              = isotropic ; uniform scaling of box vectors
   

   tau_p                   = 
   2.0                   ; time constant, in ps
   

   ref_p                   = 
   1.0                   ; reference pressure, in bar
   

   compressibility         = 
   4.5e-
   5                ; isothermal compressibility of water, bar^-
   1
   

   refcoord_scaling        = com
   

   ; Periodic boundary conditions
   

   pbc                     = xyz ; 
   3-D PBC
   

   ; Velocity generation
   

   gen_vel                 = 
   no        ; Velocity generation is 
   off

检查平衡结果命令： 
echo 18 0| gmx energy -f npt.edr -o pressure.xvg
 
手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图14
体系压力随时间变化曲线


09

动力学模拟成品模拟


随着两个平衡阶段的完成, 体系已经在需要的温度和压强下平衡好了。我们现在可以放开位置限制并进行成品MD以收集数据了， 这个过程跟前面的类似。运行grompp时, 我们还要用到检查点文件(在这种情况下,其中包含了压力耦合信息)。 我们要进行一个20 ns的MD模拟, 命令如下: 
  

   gmx grompp -f ../MDP/md.mdp -c npt.gro -t npt.cpt -p topol.top -o md.tpr
   

   
gmx mdrun -v -deffnm md

MDP参数文件内容为： 
  

   title                   = OPLS Lysozyme NPT equilibration
   

   ; Run parameters
   

   integrator              = md ; leap-frog integrator
   

   nsteps                  = 
   10000000    ; 
   2 * 
   10000000 = 
   20000 ps (
   20 ns)
   

   dt                      = 
   0.002     ; 
   2 fs
   

   ; Output control
   

   nstxout                 = 
   0         ; suppress bulky .trr file by specifying
   

   nstvout                 = 
   0         ; 
   0 for output frequency of nstxout,
   

   nstfout                 = 
   0         ; nstvout, and nstfout
   

   nstenergy               = 
   5000      ; save energies every 
   10.0 ps
   

   nstlog                  = 
   5000      ; update log file every 
   10.0 ps
   

   nstxout-compressed      = 
   5000      ; save compressed coordinates every 
   10.0 ps
   

   compressed-x-grps       = System ; save the whole system
   

   ; Bond parameters
   

   continuation            = 
   yes       ; Restarting after NPT
   

   constraint_algorithm    = lincs ; holonomic constraints
   

   constraints             = h-bonds ; bonds involving H are constrained
   

   lincs_iter              = 
   1         ; accuracy of LINCS
   

   lincs_order             = 
   4         ; also related to accuracy
   

   ; Neighorsearching
   

   cutoff-scheme           = Verlet ; Buffered nei
   ghor searching
   

   ns_type                 = grid ; search nei
   ghoring grid cells
   

   nstlist                 = 
   10        ; 
   20 fs, largely irrelevant with Verlet scheme
   

   rcoulomb                = 
   1.0       ; short-range electrostatic cut
   off (in nm)
   

   rvdw                    = 
   1.0       ; short-range van der Waals cut
   off (in nm)
   

   ; Electrostatics
   

   coulombtype             = PME ; Particle Mesh Ewald for long-range electrostatics
   

   pme_order               = 
   4         ; cubic interpolation
   

   fourierspacing          = 
   0.16      ; grid spacing for FFT
   

   ; Temperature coupling is on
   

   tcoupl                  = V-rescale ; modified Berendsen thermostat
   

   tc-grps                 = Protein Non-Protein ; two coupling groups - more accurate
   

   tau_t                   = 
   0.1     
   0.1           ; time constant, in ps
   

   ref_t                   = 
   300     
   300           ; reference temperature, 
   one for each group, in K
   

   ; Pressure coupling is on
   

   pcoupl                  = Parrinello-Rahman ; Pressure coupling 
   on in NPT
   

   pcoupltype              = isotropic ; uniform scaling of box vectors
   

   tau_p                   = 
   2.0                   ; time constant, in ps
   

   ref_p                   = 
   1.0                   ; reference pressure, in bar
   

   compressibility         = 
   4.5e-
   5                ; isothermal compressibility of water, bar^-
   1
   

   ; Periodic boundary conditions
   

   pbc                     = xyz ; 
   3-D PBC
   

   ; Dispersion correction
   

   DispCorr                = EnerPres ; account for cut-
   off vdW scheme
   

   ; Velocity generation
   

   gen_vel                 = 
   no        ; Velocity generation is 
   off

最后模拟耗时较长，而在北鲲云超算平台上供用户选择的GPU型号非常多，比如这里教程使用的是V100显卡，计算20ns的上述体系只需要

手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图15
北鲲云计算平台可供选择GPU型号示例

动力学模拟运行完毕后提示运行20ns的动力学模拟使用了 3个小时20分钟，速度可达每天143.41 ns，可谓速度飞快~ 
  

   starting mdrun 
   'HYDROLASE in water'
   

   10000000 steps, 
   20000.0 ps.
   
step 
   80: timed 
   with pme grid 
   44 
   44 
   44, coulomb cutoff 
   1.000: 
   140.1 M-cycles
   
step 
   160: timed 
   with pme grid 
   40 
   40 
   40, coulomb cutoff 
   1.091: 
   135.9 M-cycles
   
step 
   240: timed 
   with pme grid 
   36 
   36 
   36, coulomb cutoff 
   1.212: 
   147.5 M-cycles
   
step 
   320: timed 
   with pme grid 
   32 
   32 
   32, coulomb cutoff 
   1.363: 
   175.3 M-cycles
   
step 
   400: timed 
   with pme grid 
   36 
   36 
   36, coulomb cutoff 
   1.212: 
   145.9 M-cycles
   
step 
   480: timed 
   with pme grid 
   40 
   40 
   40, coulomb cutoff 
   1.091: 
   138.7 M-cycles
   
step 
   560: timed 
   with pme grid 
   42 
   42 
   42, coulomb cutoff 
   1.039: 
   141.3 M-cycles
   
step 
   640: timed 
   with pme grid 
   44 
   44 
   44, coulomb cutoff 
   1.000: 
   143.6 M-cycles
   
              optimal pme grid 
   40 
   40 
   40, coulomb cutoff 
   1.091
   
step 
   9999900, remaining wall clock time: 
   0 s
   
Writing final coordinates.
   
step 
   10000000, remaining wall clock time: 
   0 s
   

   
NOTE: The GPU has >
   25% less load than the CPU. This imbalance causes
   
      performance loss.
   

   
               Core t (s) Wall t (s) (%)
   
       Time: 
   322366.259    
   12049.380     
   2675.4 
  

  

   3h20:
   49
   
                 (ns/day) (hour/ns)
   
Performance: 
             143.410        
   0.167
   

   
gcq
   #450: "Above all, don't fear difficult moments. The best comes from them." (Rita Levi-Montalcini)

利用Chimera的MD movie工具或者其他可视化软件对动力学模拟成品进行动画展示，模拟部分结果如下所示，绿色原子表示添加的离子。
Chimera展示动力学模拟动画

利用rms工 具计算体系中蛋白质骨架的RMSD随时间波动情况，命令如下 
echo 4 4| gmx rms -f md.xtc -s md.tpr -o rmsd

利用xmgrace查看结果 
xmgrace -nxy rmsd.xvg

手把手教你用Gromacs完成溶菌酶在水中的动力学模拟的图16
动力学分析RMSD随模拟时间变化曲线


10

北鲲云超算平台作业提交


上述教程非常详细的展示了如何用Gromacs进行动力学模拟，值得注意的是，上述教程中的命令可以在单机完成，也可以上述所有命令写成作业脚本进行提交。提交脚本命令： 
sbatch -N 1 -p g-v100-1 -c 12 md-gromacs.sh
其中，-N为节点的数量，这里输入的是 1。-p为选择的PARTITION，这里使用的是V100卡（g-v100-1 ）。

md-gromacs.sh脚本的内容涵盖上述教程中的所有命令，根据北鲲云超算平台的指南需要在脚本开头加上导入gromacs模块，如果申请了GPU需要将GPU模块也导入(1-6行)，具体脚本内容如下： 
  

   module 
   add GROMACS/
   2020-fosscuda-
   2019
   b
   
export GMX_GPU_DD_COMMS=true
   
export GMX_GPU_PME_PP_COMMS=true
   
export GMX_FORCE_UPDATE_DEFAULT_GPU=true
   
ntomp=
   "$SLURM_CPUS_PER_TASK"
   
export OMP_NUM_THREADS=$ntomp
   

   

   echo 
   4 | gmx pdb2gmx -
   f protein.pdb -
   o protein_processed.gro -water tip3p -ignh -merge 
   all
   

   
gmx editconf -
   f protein_processed.gro -
   o pro_newbox.gro -
   c -d 
   1.0 -bt cubic
   

   
gmx solvate -
   cp pro_newbox.gro -
   cs spc216.gro -
   o pro_solv.gro -
   p topol.top
   

   
gmx grompp -
   f ../MDP/ions.mdp -
   c pro_solv.gro -
   p topol.top -
   o ions.tpr
   

   

   echo 
   13| gmx genion -s ions.tpr -
   o pro_solv_ions.gro -
   p topol.top -pname NA -nname CL -neutral
   
          
   
gmx grompp -
   f ../MDP/minim.mdp -
   c pro_solv_ions.gro -
   p topol.top -
   o 
   em.tpr
   
          
   
gmx mdrun -v -deffnm 
   em
   

   

   echo 
   10 
   0 | gmx energy -
   f 
   em.edr -
   o potential.xvg
   
 
   
gmx grompp -
   f ../MDP/nvt.mdp -
   c 
   em.gro -r 
   em.gro -
   p topol.top -
   o nvt.tpr
   

   
gmx mdrun -deffnm nvt
   
          
   

   

   echo 
   16 
   0 |gmx energy -
   f nvt.edr -
   o temperature.xvg
   

   
gmx grompp -
   f ../MDP/npt.mdp -
   c nvt.gro -r nvt.gro -t nvt.cpt -
   p topol.top -
   o npt.tpr
   

   
gmx mdrun -deffnm npt
   

   

   echo 
   18 
   0| gmx energy -
   f npt.edr -
   o pressure.xvg
   
 
   
gmx grompp -
   f ../MDP/md.mdp -
   c npt.gro -t npt.cpt -
   p topol.top -
   o md.tpr
   

   
gmx mdrun -v -deffnm md
   

   

   echo 
   4 
   4| gmx rms -
   f md.xtc -s md.tpr -
   o rmsd.xvg


以上就是本次教程的所有内容，而所有操作只需要可以登录北鲲云超算平台在线操作即可，无需自己配备高性能的计算机，和为繁琐的工具安装浪费时间了。

经过最近这段时间的测试， 以下是北鲲云超算平台吸引我的几点优势：

  • 高性能计算资源，极大的节省计算成本

  • 支持海量的计算工具，且开箱即用

  • 计算任务进度实时追踪提示的贴心服务

  •  详细耐心的技术咨询

登录后免费查看全文
立即登录
土石混合体多相多物理场耦合数值仿真 CINNO Research | 2022年上半年紫光展锐中国智能机SoC出货量同比增加38%
默认 最新
当前暂无评论，小编等你评论哦！
点赞 1 评论 1 收藏
关注