Gromacs - 课题组新闻 - 汕头大学医学院东庚课题组

Problem： The first residue has type 'Protein', while residue KAL10 is of type 'Other'. in pdb2gmx

Check the residuetypes.dat file in the GROMACS library directory path/to/gromacs/share/gromacs/top/

copy residuetypes.dat into current directory, then add "KAl protein" into residuetypes.data file.

cmake .. -DGMX_BUILD_OWN_FFTW=ON -DCMAKE_INSTALL_PREFIX=/home/Geng/Program/GROMACS

^{cmake .. -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DGMX_MPI=on -DGMX_BUILD_MDRUN_ONLY=ON -DCMAKE_INSTALL_PREFIX=/home/Geng/Program/GROMACS}

GROMACS Tutorials

http://www.strodel.info/index_files/lecture/html/tutorial.html

For maximum performance you will need to examine how you will use GROMACS and what hardware you plan to run on. Often OpenMP parallelism is an advantage for GROMACS, but support for this is generally built into your compiler and detected automatically.

GROMACS can run in parallel on multiple cores of a single workstation using its built-in thread-MPI. No user action is required in order to enable this.

If you wish to run in parallel on multiple machines across a network, you will need to have an MPI library installed that supports the MPI 1.3 standard, and wrapper compilers that will compile code using that library.

The GROMACS team recommends OpenMPI version 1.6 (or higher), MPICH version 1.4.1 (or higher), or your hardware vendor’s MPI installation. The most recent version of either of these is likely to be the best. More specialized networks might depend on accelerations only available in the vendor’s library. LAM-MPI might work, but since it has been deprecated for years, it is not supported.

Typical installation

As above, and with further details below, but you should consider using the following CMake options with the appropriate value instead of xxx :

GROMACS的安装方法

文/Sobereva@北京科音

First release: 2019-Jan-2 Last update: 2020-Jan-6

本文对最流行的分子动力学GROMACS程序在Linux下的安装方法进行详细说明。当新出的GROMACS版本的安装方法和本文所述出现较大差异时，本文也会做相应的更新。PS：之前笔者也写过老版本GROMACS安装方法，见《Gromacs 5.1.1与4.6.7编译方法》（http://sobereva.com/247）和《Gromacs 4.0.4、4.5.5版安装方法》（http://sobereva.com/29），但这俩文章对于目前版本来说已经没有意义了。

下面介绍的是GROMACS 2018.4版的安装，对GROMACS 2018及之后的其它版本经测试也完全适用。本文是针对计算化学工作者最常用的CentOS 7.6操作系统写的，对于其它Linux系统的用户，安装方法可能与本文有异。本文的CentOS 7.6是按照《在VMware 15中安装CentOS 7.6的完整过程视频演示》（http://sobereva.com/454）演示的方式完全新装的。本文假定读者是root用户，程序将被安装到/sob目录下。如果你是普通用户，请随机应变，恰当设置路径。本文所示的安装过程中主机全程都能访问Internet。本文编译用的C++编译器是操作系统自带的gcc，虽然用Intel C++编译器也可以，但编译出的程序的速度没有显著差别。

从GROMACS 2020开始，要求gcc版本>=5，而CentOS 7.x的gcc版本是4.8.5，因此2020版没法直接在CentOS 7.x下装。要么升级gcc版本（有一定风险，方法自行google），要么用老一点的GROMACS，要么用CentOS >=8.0版。

下面安装的是纯CPU运算、单精度、能单机并行但不能跨节点并行的版本。如果需要GPU加速或跨节点或双精度运算，看文末的附注。本文用make的时候使用了-j选项以通过并行编译降低编译耗时，但个别情况下可能导致编译出错，在虚拟机下编译还有卡住的可能，届时请去掉-j再试。下文的安装过程有全程视频演示，初学者不熟悉Linux的话请严格效仿者安装：https://www.bilibili.com/video/av39749252/。

1 安装cmake 3.x

GROMACS 2018需要系统里有cmake 3.x才能编译。CentOS 7.6自带的cmake版本太老，因此需要先装cmake 3.x。

首先运行以下命令，添加EPEL源
yum install epel-release

然后在终端里输入yum install cmake3即可下载和安装cmake包，遇到提示的时候都输入y。之后输入cmake3 /V命令，如果显示出了3.x的版本号就说明没问题了。

注1：如果用yum的时候出现乱七八糟的提示安装不了，把操作系统重启一下往往就好了。

注2：如果你没有管理员权限或者不方便用yum安装cmake，可以去http://www.cmake.org/cmake/resources/software.html下载cmake包，然后解压，进入其目录，运行
./bootstrap --prefix=/sob/cmake3
make -j
make install
就新产生了/sob/cmake3目录。然后删掉cmake解压的目录。之后在~/.bashrc里加入export PATH=$PATH:/sob/cmake3/bin。重新进入终端后，cmake命令就可以用了。

2 安装FFTW库

GROMACS 2018依赖于快速傅立叶变换库FFTW 3.3.8，可以在http://www.fftw.org/fftw-3.3.8.tar.gz下载。将其压缩包解压，进入此目录后运行
./configure --prefix=/sob/fftw338 --enable-sse2 --enable-avx --enable-float --enable-shared
以上语句代表FFTW将被安装到/sob/fftw338目录。如果你的CPU相对较新，支持AVX2指令集，可再加上--enable-avx2选项以获得更好性能。

然后运行make -j install开始编译，过一会儿编译完毕后，就出现了/sob/fftw338目录。然后可以把FFTW的解压目录和压缩包删掉了。

3 安装GROMACS

下载GROMACS 2018.4压缩包，地址为ftp://ftp.gromacs.org/pub/gromacs/gromacs-2018.4.tar.gz。然后将之拷到/sob目录下解压。进入解压后的目录，在终端里依次运行
mkdir build
cd build
export CMAKE_PREFIX_PATH=/sob/fftw338
cmake3 .. -DCMAKE_INSTALL_PREFIX=/sob/gmx2018.4
make install -j
Intel四核机子下，不到10分钟就能编译完毕。

此时程序就被编译和安装到了/sob/gmx2018.4目录下。修改用户目录下的.bashrc文件（比如运行gedit ~/.bashrc），在末尾加入source /sob/gmx2018.4/bin/GMXRC，然后保存。

之后关闭终端窗口，再次打开终端，输入gmx -version，看看是否输出了GROMACS的相关信息，是的话就说明安装成功了。之后可以把GROMACS压缩包和解压出来的目录删掉。

注1：在安装GROMACS过程中自动安装FFTW库

实际上，FFTW库可以不必手动安装，因为可以在安装GROMACS时自动下载并安装FFTW库。但由于国内链接FFTW官网服务器往往比较慢，自动下载FFTW库可能中途卡住或者过程巨慢，因此还是建议手动方式安装FFTW库。如果你确实打算自动安装FFTW库的话，将上文第2节直接忽略掉，也不必设export CMAKE_PREFIX_PATH=/sob/fftw338，把第3节的cmake3那一步额外加上-DGMX_BUILD_OWN_FFTW=ON选项即可，这样编译GROMACS时就会自动在FFTW官网上下载FFTW包并自动编译之。

注2：编译支持CUDA GPU加速的版本

GROMACS目前支持对nVidia的GPU通过CUDA方式的加速，也支持以OpenCL方式对其它厂商的GPU实现GPU加速。对于用CUDA方式加速，先去https://developer.nvidia.com/cuda-downloads下载CUDA toolkit并安装到默认路径，之后编译GROMACS方法同前，区别仅是cmake3这一步额外加上-DGMX_GPU=ON -DCUDA_TOOLKIT_ROOT_DIR=/usr/local/cuda选项（以实际CUDA tookit安装路径为准）
之后运行gmx mdrun运算时，会自动使用机子里的GPU进行加速。如果又不想使用GPU加速了，那还得按照上文方式编译一个只支持CPU运算的版本并放到不同的路径，并且把.bashrc改成source那个版本的GMXRC。

注3：编译双精度版本

一般计算只需要按照前述编译单精度版本就够了，但如果模拟刚开始就崩溃，有时候用双精度版本可解决，但计算比单精度版慢将近一倍、trr/edr等文件体积大一倍。另外，做正则振动分析时在能量极小化和对角化Hessian矩阵的时候一般也需要用双精度版以确保数值精度。注意，编译双精度版本时不支持GPU加速。

要编译双精度版本的话，先按照前文方式编译一遍单精度版本，毕竟这之后在研究中肯定也得用。然后再重复一遍编译过程，但是在编译FFTW时去掉--enable-float，并且在使用cmake3命令时额外加上-DGMX_DOUBLE=ON选项。双精度版本的GROMACS可执行文件是gmx_d，而单精度是gmx，因此单精度和双精度的可执行文件可以同时存在于同一目录，互不冲突。

注4：编译GROMACS的MPI版本

GROMACS跨节点并行计算需要MPI库，支持OpenMPI>=1.6、MPICH>=1.4.1。在编译这种GROMACS之前首先要安装MPI库，这里用OpenMPI。去http://www.open-mpi.org下载OpenMPI最新版本，解压并进入此目录后运行以下命令，就会编译并安装OpenMPI到/sob/openmpi目录：
./configure --prefix=/sob/openmpi
make all install -j
之后在用户目录下的.bashrc末尾加入以下两行
export PATH=$PATH:/sob/openmpi/bin
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/sob/openmpi/lib
然后重新进入终端使以上语句生效。之后编译GROMACS的方法同前，但在cmake3那一步额外加上-DGMX_MPI=ON选项。编译出来的可执行文件是gmx_mpi，比单机版本的可执行文件多了_mpi后缀。运行时候使用比如这样的命令：mpirun -np 16 gmx_mpi mdrun。
注：对于root用户，OpenMPI要求每次执行mpirun命令都得带着-allow-run-as-root选项才行，这很烦人，但可以通过在编译OpenMPI之前修改OpenMPI的源代码来避免，见《root用户在用openmpi并行计算时避免加--allow-run-as-root的方法》（http://sobereva.com/409）。

顺带一提，笔者在答疑时经常看到有人明明用的是单机并行，却非要装个MPI版GROMACS，这需要批评。因为这不仅需要多做一步，而且比起用默认方式基于thread-MPI和OpenMP的并行方式效率还更低，因此单机并行装MPI版完全是自取其辱。

注5：关于AVX512指令集

编译GROMACS时，如果你的CPU相对较新，支持AVX512指令集的话，程序可能会自动基于AVX512指令集进行编译。但如果你的机子上的gcc版本太老，不支持AVX512的话就会报错。如果根据屏幕上的提示确认是由于gcc不支持AVX512指令集而无法进行编译，可以按照http://bbs.keinsci.com/thread-12687-1-1.html中的做法更新gcc，或者索性不用AVX512指令集，即cmake3那一步加上-DGMX_SIMD=AVX2_256来使用AVX2指令集。

If you have prmtop and inpcrd（prmcrd/mdrest） files for amber, just convert them to gromacs formate with the following script

Jason Swails' parmed program can make this an easy task. The firstthing you need to do is create your prmtop and inpcrd files using tleap.After that, write a script (let's call it convert.py):

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

做genrestr、index
    gmx genrestr -f ligand.gro -o posre_ligand.itp -fc 1000 1000 1000

    gmx make_ndx -f em.gro -o index.ndx
    >1|13                (protein|ligand)
    >q
    注意，index文件中添加的离子需要没有添加到溶剂当中，需要手动添加，也就是将离子的原子信息直接添加至溶剂的原子信息后。若不修改，则会报错，无法识别到那几个离子。

two ways to do restraints depends on the format of .top file for ligand, i.e. whether the atom number of ligand is from 1 or following protein. e.g.

[ moleculetype ]
; Name nrexcl
CV9 3

[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
; residue 1 CV9 rtp CV9 q 0.0
1 c3 1 CV9 C37 1 -0.08510000 12.010000 ; qtot -0.085100
3 hc 1 CV9 H371 2 0.04170000 1.008000 ; qtot -0.043400
3 hc 1 CV9 H372 3 0.04170000 1.008000 ; qtot -0.001700
4 hc 1 CV9 H373 4 0.04170000 1.008000 ; qtot 0.040000
5 c3 1 CV9 C10 5 -0.02750000 12.010000 ; qtot 0.012500
6 hc 1 CV9 H10 6 0.06770000 1.008000 ; qtot 0.080200
7 c2 1 CV9 C9 7 -0.15820000 12.010000 ; qtot -0.078000

Or

; residue 1899 GLN rtp CGLN q -1.0
1683 N 1899 GLN N 1683 -0.3821 14.01 ; qtot -7.382
1684 H 1899 GLN H 1684 0.2681 1.008 ; qtot -7.114
1685 CT 1899 GLN CA 1685 -0.2248 12.01 ; qtot -7.339
1686 H1 1899 GLN HA 1686 0.1232 1.008 ; qtot -7.216
1687 CT 1899 GLN CB 1687 -0.0664 12.01 ; qtot -7.282
1688 HC 1899 GLN HB1 1688 0.0452 1.008 ; qtot -7.237
1689 HC 1899 GLN HB2 1689 0.0452 1.008 ; qtot -7.192
1690 CT 1899 GLN CG 1690 -0.021 12.01 ; qtot -7.213
1691 HC 1899 GLN HG1 1691 0.0203 1.008 ; qtot -7.192
1692 HC 1899 GLN HG2 1692 0.0203 1.008 ; qtot -7.172
1693 C 1899 GLN CD 1693 0.7093 12.01 ; qtot -6.463
1694 O 1899 GLN OE1 1694 -0.6098 16 ; qtot -7.072
1695 N 1899 GLN NE2 1695 -0.9574 14.01 ; qtot -8.03
1696 H 1899 GLN HE21 1696 0.4304 1.008 ; qtot -7.599
1697 H 1899 GLN HE22 1697 0.4304 1.008 ; qtot -7.169
1698 C 1899 GLN C 1698 0.7775 12.01 ; qtot -6.392
1699 O2 1899 GLN OC1 1699 -0.8042 16 ; qtot -7.196
1700 O2 1899 GLN OC2 1700 -0.8042 16 ; qtot -8
; LIGAND
1701 nc 1 LIG N1 1701 -0.32110 14.010000 dnc 0.00000 14.010000
1702 cd 1 LIG C2 1702 0.33780 12.010000 dcd 0.00000 12.010000
1703 cd 1 LIG C3 1703 0.37350 12.010000 dcd 0.00000 12.010000
1704 nc 1 LIG N4 1704 -0.42800 14.010000 dnc 0.00000 14.010000
1705 cd 1 LIG C5 1705 0.01760 12.010000 dcd 0.00000 12.010000
1706 cc 1 LIG C6 1706 -0.33950 12.010000 dcc 0.00000 12.010000
1707 cc 1 LIG C7 1707 0.44880 12.010000 dcc 0.00000 12.010000
1708 nd 1 LIG N8 1708 -0.56910 14.010000 dnd 0.00000 14.010000
1709 na 1 LIG N9 1709 0.24150 14.010000 dna 0.00000 14.010000
1710 c3 1 LIG C10 1710 -0.24340 12.010000 dc3 0.00000 12.010000

if the former case, you can put the restraints in the [moleculetype], and renumber the atome number in the .itp file from 1 to * same as them in the .top file. then run gmx grompp -f enmin.41.mdp -c solv_ions.gro -p gromacs.top -r solv_ions.gro -o em.tpr

if the latter case, you can put the restraints in the any postion in the .top file, and the index.ndx file must be included in the tpr file with gmx grompp -f enmin.41.mdp -c solv_ions.gro -p gromacs.top -r solv_ions.gro -n index.ndx -o em.tpr

#for gromacs
source /home/Geng/Program/GROMACS/bin/GMXRC

gmx grompp -f enmin.00.mdp -c solv_ions.gro -p gromacs.top -o em.tpr
gmx mdrun -deffnm em -ntmpi 1 -ntomp 16
gmx grompp -f nvt.00.mdp -c em.gro -p gromacs.top -o nvt.tpr
gmx mdrun -deffnm nvt -ntmpi 1 -ntomp 16
gmx grompp -f npt.00.mdp -c nvt.gro -p gromacs.top -o npt.tpr
gmx mdrun -deffnm npt -ntmpi 1 -ntomp 16
gmx grompp -f prod.00.mdp -c npt.gro -p gromacs.top -o prod.tpr -maxwarn 1
gmx mdrun -deffnm prod -ntmpi 1 -ntomp 16

常用的gromacs文件有以下类型：
1.molecular topology files
初学分子模拟的我暂时称之为分子拓扑文件（不知道别人都是怎么翻译的），top file是pdb2gmx自动讲PDB格式文件转化生成的

2.molecular structure files
这种结构文件包括pdb gro两种，pdb2gmx将pdb文件转化为top文件的同时也将pdb转化为了gro。gro和pdb都是结构文件，它们的不同在于格式，gro file还保存了velocities当然，如果你用不到v时，也可直接使用pdb格式的结构文件
另外，genbox程序用于生成研究对象周围箱中的溶剂分子，首先需要定义合适的盒子的尺寸大小。genbox的output file是gro。genbox同时还讲原有的top文件修改为假如溶剂分子后的top文件

3.molecular dynamics parameter files (mdp)
mdp files其实是一个参数文件，包括了要做分子动力学模拟所需的例如time-step number of steps tempratures pressure等等的参速

5.run input file (tpr)
将上面提到的四种文件组合生成tpr文件，也就是
top + gro + mdp + ndx = tpr
tpr file包括了要run分子动力学模拟所需要的所有信息
grompp程序处理所有的input文件生成tpr

6 trajactory files (trr)
这姑且成为轨道文件，读取trr文件对我来说走了很多弯路，曾经试图尝试使用guassview和vmd1.8.4等读图工具读取。后来跑到smth得知读取trr在gromacs中有一种自带程序 ngmx。

# topology from PDB file, trajectory from DCD file
load C:\Users\Administrator.SC-201811140916\Desktop\sampletrajectory.pdb
load_traj C:\Users\Administrator.SC-201811140916\Desktop\sampletrajectory.dcd
# gromacs trajectory, using "mytraj" as object name
load sampletrajectory.gro, mytraj
load_traj sampletrajectory.xtc, mytraj
# desmond trajectory
load sample-out.cms, mytraj
load_traj sample_trj/clickme.dtr, mytraj
# playing through states, memory optimized (but eventually slower)
set defer_builds_mode, 3
mplay

integrator = md ; #指定积分算法；md:蛙跳牛顿积分算法，用于平衡动力学积分

nstenergy = 5000 ; #能量保存到轨迹文件的频率，必须是nstcalcenergy的倍数

continuation = yes ; #初试构象不约束，不复位，用于精确的继续计算或重计算

fourierspacing = 0.16 ; #FFT傅里叶变换格点间距，默认。。。，与PME同时使用

gen_vel = no ; #速度生成；no:不生成速度。输入文件没有速度，则为0；

If you run MD simulations with a small ligand using Charmm force field, you should first get the parameters.

1. Get the structure with MOL2 format (You can use chimera to convert a PDB to mol2)

python(2 or 3) cgenff_charmm2gmx.py LIGAND ligand_fix.mol2 ligand.str charmm36-mar2019.ff

gmx editconf: 编辑模拟盒子以及转换和操控结构文件(翻译: 严立京)

gmx editconf的主要功能是对体系结构进行编辑, 也可以将通用结构格式保存或转换为.gro, .g96或.pdb等其他格式.

在分子动力学模拟中, 通常会给体系添加一个周期性的模拟盒子. gmx editconf有许多控制盒子的选项.

利用选项-box, -d和-angles可以对盒子进行修改. 除非明确使用了-noc选项, -box和-d都可以使体系在盒子内居中.

选项-bt设定盒子类型: triclinic为三斜盒子, cubic为所有边长都相等的长方体盒子(即立方体盒子), dodecahedron代表菱形十二面体盒子(等边十二面体), octahedron为截角八面体盒子(即将两个底面重合的四面体切去方向相反的两头, 同时保证所有的边长相等). 后面两种盒子是三斜盒子的特殊情况. 截角八面体三个盒向量的长度是两个相对六边形之间的最短距离. 相对于具有周期性映象距离的立方盒子, 具有相同周期距离的菱形十二面体盒子的体积是立方盒子的71%, 而截角八面体盒子的体积是立方盒子的77%.

对一般的三斜盒子, -box的参数是三个实数, 为长方体的边长. 对于立方盒子, 菱形十二面体盒子或者截面八面体盒子, 选项-box只需要提供一个数值, 即盒子边长.

-d选项指定体系中的原子到盒子编边界的最小距离. 使用-d选项时, 对三斜盒子会使用体系在x, y和z方向的大小, 对立方盒子, 菱形十二面体盒子或截角八面体盒子, 盒子的大小被设定为体系直径(原子间的最大距离)加上两倍的指定距离.

选项-angles只能与选项-box和三斜盒子一起使用才有意义, 而且不能和选项-d一起使用.

当使用-n或-ndef时, 可以指定一个索引文件, 并选择其中的一个组来计算大小和几何中心, 否则会使用整个体系的大小和几何中心.

-rotate选项可以对坐标和速度进行旋转. 如-rotate 0 30 0表示将体系绕Y轴沿顺时针方向旋转30度.

-princ选项将体系(或体系某一部分)的主轴与坐标轴平齐, 并且最长的轴沿x轴方向. 这可以减小盒子的体积, 特别当分子为长条形时. 但是注意分子在纳秒的时间尺度内可能发生明显的旋转, 所以使用时要小心.

缩放会在任何其他操作之前进行. 可以对盒子和坐标进行缩放以得到一定的密度(选项-density). 注意如果输入是.gro文件的话, 密度可能不够精确. -scale选项的一个特性是, 当某一维度的缩放因子为-1时, 可以得到体系相对于一个平面的镜面映象. 当三个维度的缩放因子都是-1时, 可以获得体系相对于坐标原点的对称映象.

组的选择是在其他所有操作都完成之后进行的. 在程序输出时, 可以只输出体系中的某一个组, 或者某一个部分, 还可以建立划分更细致的索引文件, 以便进行更加细致的选择.

可以粗略地去除体系的周期性. 当去除周期性时, 输入文件最底部的盒向量必须保证正确, 这非常重要, 因为gmx editconf去除周期性的算法十分简单, 只是将原子坐标直接减去盒子边长.

当输出.pdb文件时, 可以使用-bf选项添加B因子. B因子可以从文件中读取, 格式如下: 第一行声明文件中所含B因子数值的个数, 从第二行开始, 每行声明一个索引号, 后面跟着B因子. 默认情况下, B因子将附加到每个残基上, 每个残基一个数值, 除非索引大于残基数目或者设定了-atom选项. 显然, 可以添加任何类型的数值数据而不仅仅是B因子. -legend选项将生成一列CA原子, 其B因子的范围为所用数据的最小值到最大值, 可以有效地作为查看的图例, 便于可视化软件显示.

使用-mead选项时可以生成一个特殊的.pdb文件(.pqr), 它可用于MEAD静电程序(泊松玻尔兹曼方程求解器). 使用这个选项的前提条件是输入文件必须为运行输入文件(如tpr), 因为这样的文件中才包含了力场参数. 输出文件中的B因子段为原子的范德华半径而占有率段为原子的电荷.

-grasp选项的作用与上一选项类似, 只不过互换了电荷与半径的位置, 电荷位于B因子段, 而半径位于占有率段.

选项-align可以将特定组的主轴与给定的向量平齐, -aligncenter选项指定可选的旋转中心.

最后, 使用选项-label, gmx editconf可以为.pdb文件添加一个链标识符. 如果一个文件中不同残基属于不同肽链, 那么这个选项可以为残基指定肽链归属, 这样不但有利于可视化, 在使用一些程序如Rasmol进行分析时也很有帮助, 在建立模拟体系时也十分方便.

对一些软件包(如GROMOS), 会使用对立方盒子进行角截断的方法生成截角八面体, 为转换这种截角八面体文件, 可使用以下命令:

已知问题

补充说明

输入/输出文件选项
选项	默认值	类型	说明
`-f [<.gro/.g96/...>]`	conf.gro	输入	结构文件: gro g96 pdb brk ent esp tpr tpb tpa
`-n [<.ndx>]`	index.ndx	输入, 可选	索引文件
`-o [<.gro/.g96/...>]`	out.gro	输出, 可选	结构文件: gro g96 pdb brk ent esp. 放进盒子里面的分子坐标文件.
`-mead [<.pqr>]`	mead.pqr	输出, 可选	用于MEAD的坐标文件
`-bf [<.dat>]`	bfact.dat	输入, 可选	通用数据文件

控制选项
选项	默认值	说明
`-nice <int>`	0	设置优先级
`-[no]w`	no	程序结束自动打开输出的.xvg, .xpm, .eps和.pdb文件
`-[no]ndef`	no	从默认索引组选择输出
`-bt <enum>`	triclinic	用于`-box`和`-d`的盒子类型: triclinic, cubic, dodecahedron, octahedron
`-box <vector>`	0 0 0	盒向量长度 (a,b,c). 即自定义的盒子大小
`-angles <vector>`	90 90 90	盒向量之间的角度 (bc,ac,ab)
`-d <real>`	0	溶质分子与盒子之间的距离
`-[no]c`	no	使分子在盒子内居中(`-box`和`-d`选项暗含此选项)
`-center <vector>`	0 0 0	几何中心的坐标
`-aligncenter <vector>`	0 0 0	平齐的旋转中心
`-align <vector>`	0 0 0	与目标向量平齐
`-translate <vector>`	0 0 0	平移
`-rotate <vector>`	0 0 0	绕X, Y和Z轴的旋转角度, 单位为度
`-[no]princ`	no	使分子取向沿其主轴
`-scale <vector>`	1 1 1	缩放因子
`-density <real>`	1000	通过缩放使输出盒子的密度(g/L)为指定值
`-[no]pbc`	no	移除周期性(使分子保持完整)
`-resnr <int>`	-1	从resnr开始重新对残基进行编号
`-[no]grasp`	no	在B因子段存储原子电荷, 在占有率段存储原子半径.
`-rvdw <real>`	0.12	如果在数据库中找不到范德华半径或者拓扑文件中不存在参数, 将使用默认的范德华半径(单位nm). 可用于处理缺少力场参数的原子.
`-[no]sig56`	no	使用rmin/2(范德华势能的最小点对应距离的一半)而不是σ/2(范德华半径的一半)
`-[no]vdwread`	no	从`vdwradii.dat`文件中读取范德华半径, 而不是根据力场计算半径.
`-[no]atom`	no	强制为每个原子附加B因子
`-[no]legend`	no	创建B因子图例
`-label <string>`	A	为所有残基添加链标识符, 以便指定其肽链归属
`-[no]conect`	no	当写入的时候将CONECT记录添加到`.pdb`文件中. 只有当拓扑文件存在时才可以.

在使用pdb2gmx创建了模拟分子体系之后, 可以使用editconf为你的分子创建一个模拟盒子, 也可以认为是使用editconf将分子放进一个盒子中. 这样, 你就可以往盒子里面添加水分子, 离子或者其他溶剂等等了.

-princ这个选项可以用来对齐分子, 比如使分子沿X轴对齐. 例如, 你想将分子中的两个残基沿Y轴对齐, 那么就在索引文件中将这俩个残基标记一下, 然后使用-princ, 根据提示就能对齐分子了.

三、周期性边界条件设定
      在添加溶剂环境之前，我们需要先设定模拟体型的大小：
     gmx editconf -f complex.gro -o newbox.gro -c -bt cubic -d 1.0
      注意在这里一定要加上 -c ，让你的蛋白处于盒子的中心位置，否则后面的模拟过程蛋白会跑出盒子，极其头痛，盒子的形状有几个可以选，dodecahedron用的比较多，但我这里用的是cubic，蛋白与盒子边界的最小距离设定为1.0 nm。

四、溶剂条件设定
       盒子设定好之后，就可以添加溶剂了，溶剂模型有好几种，若无特殊要求，一般使用水溶剂 spc216，同时此处需要将溶剂写进拓扑文件里：
      gmx solvate -cp newbox.gro -cs spc216.gro -p topol.top -o solv.gro。
      注意，在加完盒子、溶剂之后，一定要查看是否合理，当然你也可以自信的不去查看，但倘若此处的盒子不合理，后面的所有工作都将白费，这里有两种方法：将gro文件转换为pdb，用Pymol查看；其二，直接用VMD打开gro文件查看盒子及水分子的合理情况。-cs 可选择tip3p.gro。

五、添加离子至体系平衡
      到这里，我们已经拥有了溶剂化的蛋白质，但是体系并没有呈电中性，因此需要添加一些正离子或负离子，我们设定的是NA+和CL-，至电中性：
     gmx grompp -f em.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
     >SOL
     注意，在进行这一步之前需先将力场参数进行修改，首先将包含此力场（Amber03）的整个文件夹拷贝至当前目录，文件夹在gromacs的安装目录下（gromacs---share---amber03ff文件夹），然后修改文件夹中的forcefield.itp文件：ligand.itp中的【atomtypes】部分替换掉
forcefield.itp中的【atomtypes】部分。更为重要的是ligand.itp中【defult】、【atomtypes】、【system】、【molecules】需要删掉，否则会与立场文件及拓扑文件相冲突，会有很多报错，本人就在这些问题上兜兜转转了好几天。

  六、溶剂化体系的优化
      这个时候，我们已经拥有了一个相对较为平衡的体系，但是由于添加了溶剂、离子，容易产生原子重叠以及同种电荷相近等，还需要对溶剂化体系进行优化：
    gmx grompp -f em_real.mdp -c solv_ions.gro -p topol.top -o em.tpr
    gmx mdrun -v -deffnm em
    在开始执行优化时，-v 是为了看运行此进程需要花费的时间，当然你要想关掉窗口不关进程，也可以使用nohup：nohup gmx mdrun  -deffnm em&。
    注意，文中提到的所有.mdp文件，gromacs官网的手册中都有提供，可自行下载，记得修改文件中小分子及受体蛋白的名称，执行的时间长短也可以根据自身需要进行修改。

七、做genrestr、index
    gmx genrestr -f ligand.gro -o posre_ligand.itp -fc 1000 1000 1000

    gmx make_ndx -f em.gro -o index.ndx
    >1|13                (protein|ligand)
    >q
    注意，index文件中添加的离子需要没有添加到溶剂当中，需要手动添加，也就是将离子的原子信息直接添加至溶剂的原子信息后。若不修改，则会报错，无法识别到那几个离子。

八、做温度耦合
    现在体系已经相对松弛，现在需要引入温度因素，让系统达到一个新的松弛状态：
    gmx grompp -f nvt.mdp -c em.gro -p topol.top -n index.ndx -o nvt.tpr
    gmx mdrun -deffnm nvt
    注意，这个地方体系易崩溃，造成崩溃的原因有4：初始态是否合理；能量最小化是否收敛；盒子大小问题；温度耦合系数是否合理，步长是否过长。

九、做压力耦合
    gmx grompp -f npt.mdp -c nvt.gro -p topol.top -n index.ndx -o npt.tpr
    gmx mdrun -deffnm npt

十、正式模拟
    到这一步的时候也基本修成正果了，历经了九九八十一难，终于可以开始正式模拟:
gmx grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -n index.ndx -o md.tpr
gmx mdrun -dennm md

概述

本示例教程将引导GROMACS新用户进行一次模拟, 模拟的体系是水盒子中的蛋白质(溶菌酶, lysozyme)及离子. 我们会解释每一步中用到的输入文件以及得到的输出文件. 这些输入文件中所用的设置都很典型, 因此, 你在进行模拟时一般也可以采用这些设置.

第一步: 准备拓扑

随着5.0版本GROMACS的发布, 所有的工具程序现在都成为了gmx程序中的一个模块. 而在以前版本的GROMACS中, 这些工具都是可以单独使用的, 并且有自己的命令. 在5.0版本中, 仍然可以通过符号链接使用这些命令, 但将来的版本会废弃这种作法. 因此, 你最好还是习惯这种新的使用方式. 要得到GROMACS某一模块的帮助信息, 你可以使用下面命令中的任何一个

使用时只要将其中的(module)替换为要查询的命令的实际名称即可. 相关信息会输出到终端, 其中包括可用的算法, 选项, 需要的文件格式, 已知的缺陷和限制, 等等. 对GROMACS的新用户来说, 查看常用命令的帮助信息是了解每个命令功能的有效方式.

首先, 我们必须下载所用蛋白质的结构文件. 在本教程中, 我们将使用鸡蛋白溶菌酶(PDB代码为1AKI). 打开RCSB网站, 输入鸡蛋白溶菌酶的PDB编号, 然后点击Go,

网页上包含了很多与鸡蛋白溶菌酶有关的信息, 可以帮助你更好地了解这个蛋白质. 我们需要的是这个蛋白质的晶体结构文件, 因此点击右边的Download Files, 然后下载PDB格式的晶体结构文件.

得到结构文件之后, 你可以使用VMD, Chimera, PyMOL等可视化程序来查看一下这个蛋白质的结构. 如果还不熟悉这些程序的使用, 你可以参考网上的教程以及示例视频. 在查看轨迹方面, VMD功能非常强大, 而使用PyMOL你可以轻松地得到高质量的蛋白质结构图片. 本教程首页的蛋白质结构图片就是利用PyMOL得到的. 如果你不是很喜欢VMD和PyMOL程序的操作模式, 你或许可以试试Jmol. 此外, 还有非常多的分子结构可视化程序. 如果你只是用于查看结构, 选择你喜欢的那个即可. 上面提到的这几个程序都支持脚本, 利用脚本, 你可以更好地操控分子, 并能进行一些处理.

Points点模式

Line线模式Bonds键模式CPK模式VDW模式Tube管线模式Ribbons飘带模式NewCartoons新卡通模式

查看过这个蛋白质分子之后, 你要去掉晶体结构中的结晶水. 请使用普通的文本编辑器, 如vi/vim, emacs(Linux/Mac)或者记事本程序(Windows). 不要使用文字处理软件(例如Windows下的Word)! 它们非常笨重, 不适合快速地查看编辑纯文本文件. 如果你不喜欢vi/vim或emacs的操作模式, 可以试试gVIM或是其他软件. 对Windows的用户, 系综自带的记事本程序不是很好用, 你可以选择一些和它类似, 但功能更强大的程序, 如Notepad2, Notepad++, EmEditor, EditPlus, UltrEdit, 等等. 选择非常广泛, 请选择一个你喜欢的, 并尽可能地将熟悉它的各种功能, 这样在处理各种文件时才能得心应手.

去除结构中的结晶水时, 直接删掉PDB文件中与结晶水相关的行(残基为“HOH”的行)即可. 注意, 并不是任何时候都需要进行这个过程(比如对与活性位点结合的水分子就不能去除). 对本教程而言, 我们不需要结晶水, 所以可以将它们都去掉.

始终要注意检查.pdb文件中MISSING注释下面列出的项, 这些项代表了在晶体结构文件中缺失的那些原子或残基. 在模拟中, 末端区域的缺失可能并不会引起问题, 但缺失原子的不完整内部序列或任何氨基酸残基都将导致pdb2gmx程序运行失败. 必须使用其他软件对这些缺失的原子/残基进行建模并补充完整. 还要注意pdb2gmx不是万能的, 它无法为任意分子生成拓扑文件, 而只能用于力场中已经定义的残基(在*.rtp文件中, 一般包括蛋白质, 核酸和 非常有限 的辅酶因子, 如NAD(H)和ATP).

现在结构中已经没有结晶水了, 我们也确认了没有缺少任何需要的原子. PDB文件中应该只包含蛋白质原子, 这样就可以将其用作pdb2gmx的输入. pdb2gmx是我们用到的第一个GROMACS模块, 它的作用的是产生三个文件:

拓扑文件(默认为topol.top)包含了定义模拟中所用分子的所有必需信息, 包括非键参数(原子类型和电荷)和键合参数(键长, 键角和二面角). 当得到拓扑文件后, 我们会更详细地对它进行说明.

pdb2gmx将处理结构, 输出一些相关信息后, 提示你选择一个力场:

力场包含了将要写入到拓扑文件中的信息. 这个选择非常重要! 你必须仔细了解每个力场, 并决定哪个力场最适用于你的体系. 在本教程中, 我们选用全原子OPLS力场, 因此在命令提示行中输入15, 然后回车.

现在你已经得到了三个新的文件, 1AKI_processed.gro, topol.top和posre.itp. 1AKI_processed.gro是GROMACS格式的结构文件, 包含了力场中定义的所有原子(即, 已经将氢原子加到蛋白质中的氨基酸上了). topol.top文件是体系的拓扑文件(稍后会解释). posre.itp文件包含了用于限制重原子位置的信息(后面解释).

最后的注意事项: 许多用户认为.gro文件是必需的. 事实并非如此. GROMACS可处理很多不同的文件类型, .gro不过是一些程序输出坐标文件时所用的默认格式. 这种格式非常紧凑, 但精度有限. 如果你更愿意使用其他格式, 如PDB格式, 为你的输出文件指定合适的文件名称, 并使用.pdb作为扩展名即可. 使用pdb2gmx程序的目的在于生成与力场兼容的拓扑文件, 输出的结构不过是此目的的副产品, 只是为了方便用户. 输出结构的格式可以任意选择(参看GROMACS手册上对不同格式的说明).

第二步: 检查拓扑文件

让我们来看一下输出的拓扑文件(topol.top)中有些什么. 使用普通文本编辑器来检查其内容. 在几行注释(前面有分号;标注)之后, 你可以看到下面的语句:

此行调用了OPLS-AA力场的参数, 它位于文件的开头, 这意味着下面的所有参数都来自OPLS-AA力场. 下一重要行是[ moleculetype ], 后面是

“Protein_A”定义了分子名称, 这是因为这个蛋白质在PDB文件中被标定为A链. 对键合近邻的排除数为3. 关于排除的更多信息可从GROMACS手册上找到. 对此信息的讨论超出了本教程的范围.

下面几节包括[ bonds ], [ pairs ], [ angles ]和[ dihedrals ]. 其中一些无需解释(键, 键角, 二面角). 这些节中的参数和函数类型可参看GROMACS手册的第5章. 特殊的1–4相互作用包含于“pairs”(GROMACS手册5.3.4节).

从位置限制开始, 文件的其余部分涉及一些有用的/必须的拓扑的定义. pdb2gmx命令生成的“posre.itp”文件定义了平衡时用于维持原子位置的力常数(后面会详细解释).

至此“Protein_A”分子类型的定义结束. 拓扑文件的其余部分用于定义其他分子并提供体系级别的说明. 下一分子类型(默认)是溶剂, 在本例中为SPC/E模型的水分子. 水的其他典型模型包括SPC, TIP3P和TIP4P. 通过在pdb2gmx命令中使用“-water spce”选项我们选择了SPC/E水模型. Water Models对许多不同的水模型的进行了很好的总结. 但是要注意GROMACS并没有包含所有的水模型.

正如你看到的, 通过使用值为1000 kJ mol^-1 nm^-2的力常数(k_pr), 也可以对水分子进行位置限制.

最后是体系级别的定义. [ system ]指令给出了体系的名称, 在模拟中此名称将被写入到输出文件中. [ molecules ]指令列出了体系中的所有分子.

任何时候, 如果不能满足这些明确的要求, 运行grompp程序(后面介绍)时都会产生致命错误, 如mismatched names, molecules not being found或其他各种错误.

第三步: 定义单位盒子并填充溶剂

现在你已经熟悉了GROMACS的拓扑文件, 可以继续创建体系了. 在本例中, 我们将要模拟一个简单的水溶液体系. 我们也可以模拟处于其他不同溶剂中的蛋白质或其他分子, 只要涉及到的物种有合适的力场参数.

现在你需要决定使用哪种晶胞. 对于本教程的目的而言, 我们将使用一个简单的立方盒子作为晶胞. 当你对周期性的边界条件与盒子类型有了更多了解后, 我强烈推荐你使用菱形十二面体晶胞, 因为在周期性距离相同的情况下, 它的体积大约只有立方体晶胞的71%, 因此可以减少需要加入的溶剂水分子的数目.

上面的命令将蛋白质置于盒子的中心(-c), 并且它到盒子边缘的距离至少为1.0 nm(-d 1.0). 盒子类型是立方体(-bt cubic). 到盒子边缘的距离是一个重要参数. 因为我们要使用周期性边界条件, 必须满足最小映象约定, 也就是说, 一个蛋白质永远不能“看到”它自身的周期性映象(不能与其自身有相互作用), 否则计算的力就会含有虚假的部分. 指定溶质与盒子之间的距离为1.0 nm意味着, 蛋白质分子的任意两个周期性映象之间的距离至少是2.0 nm. 对于模拟中常用的任何截断方案, 这个距离都足够了.

现在我们已经定义好了模拟盒子, 可以用溶剂(水)填充它了. 使用solvate模块添加溶剂:

使用的蛋白质构型文件(-cp)来自前面editconf步骤中的输出文件, 而溶剂的构型文件(-cs)来自标准安装的GROMACS. 我们使用的spc216.gro是通用的已平衡的三位点溶剂模型. 你也可以使用spc216.gro作为SPC, SPC/E或TIP3P水模型的溶剂构型, 因为它们都是三位点的水模型. 输出文件的名称为1AKI_solv.gro, 并且我们为solvate模块指定了拓扑文件的名称(topol.top), 这样它就能修改拓扑文件. 注意topol.top中[ molecules ]的变化:

solvate记录了增加的水分子数目, 并将其写入拓扑文件中以显示它所做的更改. 注意, 如果你使用其他的(非水)溶剂, solvate不会在拓扑文件中写入这些信息! 它自动记录更新水分子的功能是直接写在源代码中的.

第四步: 添加离子

现在我们已经有了一个带电荷的溶液体系. pdb2gmx程序的输出文件显示, 我们所用的蛋白质带有+8e的净电荷(根据它的氨基酸残基计算得到). 如果你忽略了pdb2gmx输出的这个信息, 查看一下topol.top文件中[ atoms ]指令的的最后一行, 它应该含有“qtot 8.”这一信息. 由于生命体系中不存在净电荷, 所以我们必须往我们体系中添加离子, 以保证总电荷为零.

GROMACS中添加离子的工具是genion. genion的功能是读取拓扑信息, 然后将体系中的一些水分子替换为指定的离子. genion需要的输入文件称为运行输入文件, 扩展名为.tpr. 这个文件可使用GROMACS的grompp(GROMacs Pre-Processor)模块产生, 而且后面我们运行模拟时也会用它. grompp的功能是处理坐标文件和拓扑(它描述了分子)以产生原子级别的输入文件(.tpr). .tpr文件包含了体系中所有原子的所有参数.

为了用grompp产生.tpr文件, 我们还需要一个扩展名为.mdp(molecular dynamics parameter)的输入文件. grompp会将坐标和拓扑信息与.mdp文件中设定的参数组合起来生成.tpr`文件.

.mdp文件通常用于运行能量最小化(EM)或分子动力学模拟(MD), 但在这里我们只是简单地用它来生成系统的原子描述. 一个.mdp文件的示例(后面我们将使用它)如下:

实际上, 在这个步骤中所用的.mdp文件中可使用任何合理的参数. 我通常会使用能量最小化的参数设置, 因为它非常简单而且不涉及任何复杂的参数组合. 请注意 本教程中所用的文件可能只适用于OPLS-AA力场. 其他力场的参数设置, 特别是非键参数设置可能很不一样.

现在我们得到了一个二进制的.tpr文件, 它提供了我们体系的原子级别的描述. 将此文件用于genion:

gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -nn 8

出现提示后, 选择13 “SOL”. 这意味这我们将用离子替换一些溶剂分子. 你肯定不想用离子去替换蛋白质的某一部分.

在genion命令中, 我们以结构/状态文件(-s)作为输入, 以.gro文件作为输出(-o), 以拓扑文件(-p)来反映去除的水分子和增加的离子, 并且定义了阳离子和阴离子的名称(分别使用-pname和-nname), 告诉genion只需要添加必要的离子来中和蛋白质所带的净电荷, 添加的阴离子数目为8(-nn 8). 对于genion, 除了简单地中和体系所带的净电荷以外, 你也可以同时指定-neutral和-conc选项来添加指定浓度的离子. 关于如何使用这些选项, 请参考genion的说明.

在以前版本的Gromacs中, 使用-pname和-nname指定的离子名称由力场决定, 但从4.5版本开始就完全标准化了. 指定的离子名称始终是大写的元素符号, 与[ moleculetype ]中的名称一致, 并会写入拓扑文件. 残基名称或原子名称可能会带有电荷符号(+/-), 也可能不带, 取决于力场. 不要在genion命令中使用原子名称或残基名称, 否则在下面的步骤中会导致错误.

第五步: 能量最小化

现在, 我们已经添加了溶剂分子和离子, 得到了一个电中性的体系. 在开始动力学模拟之前, 我们必须保证体系的结构正常, 原子之间的距离不会过近, 几何构型合理. 对结构进行弛豫可以达到这些要求, 这个过程称为能量最小化(EM, energy minimization).

能量最小化过程与添加离子过程差不多. 我们要再次使用grompp将结构, 拓扑和模拟参数写入一个二进制的输入文件中(.tpr), 但这次我们不需要将.tpr文件传递给genion, 而是使用GROMACS MD引擎的mdrun模块来进行能量最小化.

确保在运行genbox和genion时你已经更新了topol.top文件, 否则你会得到一堆错误信息(“number of coordinates in coordinate file does not match topology”, 坐标文件中的坐标与拓扑不匹配, 等等).

之所以使用-v选项是因为我们没什么耐心, 急于看到运行结果: 它使mdrun输出更多信息, 这样就会在屏幕上输出每步运行的情况. -deffnm选项定义了输入文件和输出文件的名称. 因此, 如果你没有对grompp输出的em.tpr进行命名, 你必须使用mdrun的-s选项明确指定它的名称. 就我们而言, 我们将得到以下文件:

有两个重要的指标来决定能量最小化是否成功. 第一个是势能(在能量最小化过程的最后输出, 即使你未使用-v选项). E_pot应当是负值, 根据体系大小和水分子的多少, 大约在10⁵–10⁶的数量级(对水中的单个蛋白质而言). 第二个重要的指标是力的最大值F_max. 我们在minim.mdp中设置的目标是emtol=1000.0, 这表示F_max的目标值不能大于1000 kJ mol^-1 nm^-1. 能量最小化完成后, 你有可能得到一个合理的E_pot, 但F_max>emtol. 如果是这样, 用于模拟时你的体系可能不够稳定. 思考一下为什么会这样, 可能需要更改一下能量最小化的参数设置(integrator, emstep等), 再试试重新进行能量最小化过程.

让我们做一些分析. em.edr文件中包含了GROMACS在能量最小化过程中记录的所有能量项. 你可以使用GROMACS的energy模块来分析任何一个.edr文件:

提示时, 输入10 0来选择势能Potential(10), 并用零(0)来结束输入. 屏幕上会显示E_pot的平均值, 得到的能量值会写入potential.xvg文件. 要利用这些数据绘图, 你可以试试Xmgrace绘图工具. 得到的结果应该和下面的差不多, 从中可以看到E_pot收敛得很好, 而且稳定.

现在我们的体系已经处于能量最小点了, 可以用它进行真正的动力学模拟了.

第六步: NVT平衡

EM可保证我们的初始结构在几何构型和溶剂分子取向等方面都合理. 为了开始真正的动力学模拟, 我们必须对蛋白质周围的溶剂和离子进行平衡. 如果我们在这时就尝试进行非限制的动力学模拟, 体系可能会崩溃. 原因在于我们基本上只是优化了溶剂分子自身, 而没有考虑溶质. 我们需要将体系置于设定的模拟温度下, 以确定溶质(蛋白质)的合理取向. 达到正确的温度(基于动能)之后, 我们要对体系施加压力直到它达到合适的密度.

还记得好久以前我们用pdb2gmx生成的posre.itp文件么? 现在它要派上用场了. posre.itp文件的目的在于对蛋白质中的重原子(非氢原子)施加位置限制(position restraining)力. 这些原子不会移动, 除非增加非常大的能量. 位置限制的用途在于, 我们可以平衡蛋白质周围的溶剂分子, 而不引起蛋白质结构的变化.

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

除注释外, 所用参数的完整解释可以在GROMACS手册中找到. 注意.mdp文件中下面的这几个参数:

提示时输入15 0来选择体系温度并退出. 得到的结果应该和下面的差不多:

从上图可以清楚地看出, 体系的温度很快就达到了目标温度(300 K), 并在平衡过程中后面的时间内保持稳定. 对于这个体系, 更短的平衡时间(50 ps)也足够了.

第七步: NPT平衡

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

该文件与NVT平衡时所用的参数文件没有太大不同. 注意添加的压力耦合部分, 其中使用了Parrinello-Rahman控压器.

我们使用grompp和mdrun, 像在NVT平衡所做的那样. 注意, 我们现在要使用-t选项以包括NVT平衡过程中的产生的检查点文件. 这个文件包含了继续模拟所需要的所有状态变量. 为使用NVT过程中得到的速度我们必须包含这个文件. 坐标文件(-c)是NVT模拟的最终输出文件.

在100 ps的平衡过程中压力值涨落很大, 这并不意外. 图中的红线为数据的移动平均值. 在整个平衡过程中, 压力的平均值为1.05 bar.

跟压力一样, 红线是密度的移动平均值. 100 ps过程中密度的平均值为998.3 kg m^-3, 比较接近实验值1000 kg m^-3与SPC/E水模型的值1008 kg m^-3. SPC/E水模型的参数给出的密度值接近水的实验值. 在整个过程中密度值都很稳定, 意味着体系的压力和密度下都平衡得很好.

请注意, 经常有人问我为什么他得到的密度值与我的结果不同. 与压力有关的性质收敛很慢, 因此你运行NPT平衡的时间必须比这里指定的稍长一些.

【李继存注】经常有人问上面两个图中的红线怎么得到. 如果你要使用累积平均值来画, 那可能需要一小段代码来完成. 但如果你只是像图中一样, 使用移动平均值来简单地平滑一下, 就很简单了.

在Xmgrace中, 依次点击菜单 Data -> Transformations -> Running averages..., 在弹出的对话框中设定Length of average, Accept即可. Length of average的具体数值要根据具体数据的特点来设, 越大, 得到的平均线越平滑. 自己试几次就知道了.

在Origin中, 依次点击菜单 分析 -> 平滑 -> 相邻平均 或 FFT滤波器, 设定平滑的点数即可. 具体数值的设置原则, 和Xmgrace中的一样.

第八步: 成品MD

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

PME负载估计可指示应该使用多少处理器进行PME计算, 多少处理器进行PP计算. 详细情况请参考GROMACS 4的相关论文和GROMACS手册. 对立方盒子, 最佳设置的PME负载为0.25(3:1 PP:PME, 我们很幸运!), 对十二面体盒子, 最佳PME负载为0.33(2:1 PP:PME). 当执行mdrun的时候, 程序会自动分配用于PP和PME计算的处理器数目. 因此, 确保计算时使用了合适的节点数(-np X选项中的值), 这样性能最好. 对本教程中的这个体系, 在24个CPU上(18 PP, 6 PME)我得到的计算速度大约是14 ns/day.

自4.0版本开始, GROMACS运行MD模拟时可以使用GPU加速器. 非键相互作用使用GPU进行计算, 而键合与PME相互作用则使用标准的CPU硬件进行计算. 当安装GROMACS(参考www.gromacs.org上的安装指导)的时候, 会自动检测存在的GPU硬件设备. 使用GPU加速的最低要求为CUDA库和SDK, 以及具有2.0计算能力的GPU卡. 这里列出了一些更常见的卡及其配置. 要使用GPU, 上面.mdp文件唯一要做的修改是添加下面一行以确保使用Verlet截断方案(GPU不支持旧的组方案):

如果可用的GPU卡超过一个, 或者想利用GROMACS支持的杂合并行方案对计算进行划分, 请参考GROMACS手册以及网络上的资料. 这些技术细节超出了本教程的范围.

第九步: 分析

现在已经完成了对蛋白质的模拟, 我们应该来分析一下我们的体系. 哪些类型的数据才是重要的呢? 这是在模拟前就要思考的一个重要问题, 所以你应该对自己的体系需要采集哪些数据类型有自己的想法. 在本教程中, 我们只介绍一些基本工具.

第一个模块是trjconv, 这是一个后处理工具, 用于处理坐标, 修正周期性或手动调整轨迹(时间单位, 帧频率等). 在本教程中, 我们要使用trjconv来处理体系中的任何周期性. 蛋白质在单元晶胞中扩散, 可能看起来会在盒子两边之间进行“跳跃”. 我们使用下面的命令来处理这种情况:

gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -ur compact

选择0("System")用于输出. 我们要基于这个“修正”后的轨迹进行分析. 先来看看结构稳定性. GROMACS内置的rms模块可用于计算RMSD, 使用下面的命令来运行这个工具:

计算最小二乘拟合RMSD和组RMSD时, 都选择4("Backbone"). -tu选项设定输出结果的时间单位为ns, 即便轨迹文件以ps为单位输出. 这是为了使输出文件更加清晰(尤其当模拟时间很长时, 100 ns比起1e+05 ps更美观). 输出显示了MD模拟前后溶菌酶结构的RMSD:

如果我们要计算相对于晶体结构的RMSD值, 可以使用下面的命令:

上面两个图都显示出RMSD大约是0.1 nm(1Å), 这表示蛋白质的结构非常稳定. 两图之间的微小差异意味着, 当t=0 ns时的蛋白质的结构与晶体结构稍有不同. 这是预期结果, 因为它已经进行了能量最小化, 而且如我们前面讨论的, 位置限制并不是100%完美的.

我们也可以将初始构型与模拟后的构型进行比较, 这样可以更直观地看出二者的区别.

蛋白质的回旋半径R_g可衡量其密实度. 如果蛋白质的折叠很稳定, 其R_g将保持一个相对稳定的值. 如果蛋白质去折叠, 它的R_g将随时间变化. 我们来分析一下模拟的溶菌酶的回旋半径:

可以看到, R_g值基本不变, 这预示着在温度为300 K时, 1 ns的时间内蛋白质很稳定, 处于紧密(折叠)的形式. 这一结果并非意外, 但说明了GROMACS具有先进的分析功能.

总结

现在你已经用GROMACS完成了一个分子动力学模拟过程, 并分析了一些结果, 本教程并不全面. 你还可以用GROMACS完成更多类型的模拟(自由能计算, 非平衡MD, 简正模式分析, 等等). 你应该阅读一些文献以及GROMACS手册, 试着调整这里提供的.mdp文件中的参数来以便使模拟更有效, 更精确.

如果你对改进这个教程有些建议, 如果你发现了错误, 或者你觉得有些地方不够清楚, 请给我发邮件jalemkul@vt.edu, 不要客气. 请注意: 这不是邀请你因为GROMACS的问题而给我发邮件. 我并不是作为一个私人家教或个人客服在为自己打广告. 那是GROMACS用户列表的事. 我可能会在那里帮助你, 但那只是作为对整个社区的服务, 而不只针对最终用户.

GROMACS程序选项概述

GROMACS各类程序(功能分类)

分析轨迹

创建拓扑与坐标文件

运行模拟

查看轨迹

处理能量

选项	类型	默认	说明
`-[no]hidden`	bool	`yes`	[隐藏] 打印隐藏选项
`-[no]quiet`	bool	`no`	[隐藏] 不打印帮助信息. 打印尽可能少的信息, 忽略所有参数输出
`-man`	enum	`html`	[隐藏] 将程序手册打印到指定文件并退出. 支持的文件格式: no, html, tex, nroff, java, ascii或completion
`-[no]debug`	bool	`no`	[隐藏] 将调试信息写入文件

gmx enemat - 从能量文件中提取能量矩阵(翻译: 赵丙春) 原始文档
gmx energy - 将能量写入xvg文件并显示平均值(翻译: 姚闯) 原始文档
gmx mdrun - 利用-rerun选项(重新)计算轨迹中每帧的能量(翻译: 王浩博) 原始文档

GROMACS Tutorials