随着中国科研的发展，有限的超算中心和自建计算集群的方式已无法满足超算用户计算业务和创新的需要，特别对降低超算总体拥有成本，无需排队获取资源，弹性规模，软件开箱即用和安全运维等核心需求，云上超算是唯一解决之道。
今天，将为大家介绍如何在北鲲云超算平台上快速完成动力学模拟计算，模拟溶菌酶（lysozyme）在水相中的结构变化。北鲲云超算平台是前段时间小伙伴安利的，经过一段时间测试，效果还不错。只需要一台上网本就可以在短时间内完成动力学模拟所有内容，而且还预装了Gromacs在内的300+应用软件，实现了开箱即用，非常的灵活便捷。接下来，就正式进入今天的主题。
水中的溶菌酶
01 结构处理
溶菌酶又称胞壁质酶或N-乙酰胞壁质聚糖水解酶，是一种能水解细菌中黏多糖的碱性酶。溶菌酶还可与带负电荷的病毒蛋白直接结合，与DNA、RNA、脱辅基蛋白形成复合体，使病毒失活。这里我们从PDB数据库下载溶菌酶蛋白质三维结构，编号为1AKI，或则利用PyMOL直接从命令行获取三维结构，然后利用PyMOL检查结构，是否缺失，并且删除除开蛋白质外的杂原子和水分子，并保存为protein.pdb文件。
fetch 1AKI
02 将输入文件上传到北鲲云超算平台
打开北鲲云超算平台，注册登录后在/home/cloudam/jobs下新建文件夹Lysozyme，将准备好的protein.pdb文件直接拖拽到该文件中，即可完成上传。上传文件和下载文件都可以利用sftp非常方便完成。
除了命令行，北鲲云还有可视化和图形界面作业提交方式。可视化无需代码，比较适合无计算机专业背景的同学，图形界面适合于单节点小规模计算任务。这两种作业提交方式大家就自行体验吧（北鲲云针对新用户有赠送200元体验金的活动，大家可以关注一下）。
北鲲云超算平台作业提交入口
北鲲云超算平台操作终端窗口
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
参数文件内容：
; 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 ; Minimization 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
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = cutoff ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
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参数文件内容：
; minim.mdp - used as input into grompp to generate em.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 ; Minimization 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
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.0 ; Short-range electrostatic cut-off
rvdw = 1.0 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions
然后键入如下命令进行能量最小化，这一步应当提交到计算节点进行，不要在管理节点直接运行。
gmx mdrun -v -deffnm em
能量最小化开始，屏幕中出现运行进度提示：
Steepest Descents:
Tolerance (Fmax) = 1.00000e+03
Number of steps = 50000
Step= 0, Dmax= 1.0e-02 nm, Epot= -2.94025e+05 Fmax= 1.53483e+05, atom= 19880
Step= 1, Dmax= 1.0e-02 nm, Epot= -3.12614e+05 Fmax= 6.09036e+04, atom= 19880
Step= 2, Dmax= 1.2e-02 nm, Epot= -3.35310e+05 Fmax= 2.40700e+04, atom= 19880
Step= 3, Dmax= 1.4e-02 nm, Epot= -3.61996e+05 Fmax= 1.00755e+04, atom= 24188
Step= 4, Dmax= 1.7e-02 nm, Epot= -3.87330e+05 Fmax= 4.66750e+03, atom= 983
Step= 5, Dmax= 2.1e-02 nm, Epot= -4.08800e+05 Fmax= 1.30036e+04, atom= 547
Step= 6, Dmax= 2.5e-02 nm, Epot= -4.13336e+05 Fmax= 2.07112e+04, atom= 547
Step= 7, Dmax= 3.0e-02 nm, Epot= -4.17185e+05 Fmax= 1.90699e+04, atom= 547
Step= 8, Dmax= 3.6e-02 nm, Epot= -4.18940e+05 Fmax= 2.96782e+04, atom= 547
Step= 9, Dmax= 4.3e-02 nm, Epot= -4.21615e+05 Fmax= 2.83858e+04, atom= 547
Step= 11, Dmax= 2.6e-02 nm, Epot= -4.26495e+05 Fmax= 6.89298e+03, atom= 547
Step= 12, Dmax= 3.1e-02 nm, Epot= -4.27186e+05 Fmax= 3.75623e+04, atom= 1785
Step= 13, Dmax= 3.7e-02 nm, Epot= -4.33871e+05 Fmax= 1.81256e+04, atom= 1785
Step= 15, Dmax= 2.2e-02 nm, Epot= -4.35813e+05 Fmax= 1.54546e+04, atom= 1785
Step= 16, Dmax= 2.7e-02 nm, Epot= -4.37056e+05 Fmax= 2.45185e+04, atom= 1785
利用如下命令检查是否有效的能量最小化：
echo 10 0 | gmx energy -f em.edr -o potential.xvg
利用xmgrace打开potential.xvg文件，结果如如下所示：
体系势能随时间变化曲线
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 neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (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步数来进一步获取平衡体系。
体系温度随时间变化曲线
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 neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (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
体系压力随时间变化曲线
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
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (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的上述体系只需要
北鲲云计算平台可供选择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.43h20: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工具或者其他可视化软件对动力学模拟成品进行动画展示，模拟部分结果如下所示，绿色原子表示添加的离子。
，时长00:11
Chimera展示动力学模拟动画
利用rms工具计算体系中蛋白质骨架的RMSD随时间波动情况，命令如下
echo 4 4| gmx rms -f md.xtc -s md.tpr -o rmsd
利用xmgrace查看结果
xmgrace -nxy rmsd.xvg
动力学分析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-2019b
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
以上就是本次教程的所有内容，而所有操作只需要可以登录北鲲云超算平台在线操作即可，无需自己配备高性能的计算机，和为繁琐的工具安装浪费时间了。
经过最近这段时间的测试，以下是北鲲云超算平台吸引我的几点优势：
高性能计算资源，极大的节省计算成本
支持海量的计算工具，且开箱即用
计算任务进度实时追踪提示的贴心服务
详细耐心的技术咨询