Problem： The first residue has type 'Protein', while residue KAL10 is of type 'Other'. in pdb2gmx
Installation
cd gromacs-2019.1
mkdir build
cd build
cmake .. -DGMX_BUILD_OWN_FFTW=ON -DCMAKE_INSTALL_PREFIX=/home/Geng/Program/GROMACS
or
cmake .. -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DGMX_MPI=on -DGMX_BUILD_MDRUN_ONLY=ON -DCMAKE_INSTALL_PREFIX=/home/Geng/Program/GROMACS
or
cmake .. -DGMX_OPENMP=ON -DCMAKE_INSTALL_PREFIX=/home/Geng/Program/GROMACS
make -j 4
make check
make install
source /home/Geng/Program/GROMACS/bin/GMXRC
NOTE
Compiling with parallelization options
2. MPI support
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.
For example, depending on your actual MPI library, use
cmake -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DGMX_MPI=on.
GROMACS的安装方法
First release: 2019-Jan-2 Last update: 2020-Jan-6
从GROMACS 2020开始，要求gcc版本>=5，而CentOS 7.x的gcc版本是4.8.5，因此2020版没法直接在CentOS 7.x下装。要么升级gcc版本（有一定风险，方法自行google），要么用老一点的GROMACS，要么用CentOS >=8.0版。
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的时候出现乱七八糟的提示安装不了，把操作系统重启一下往往就好了。
然后运行make -j install开始编译，过一会儿编译完毕后，就出现了/sob/fftw338目录。然后可以把FFTW的解压目录和压缩包删掉了。
此时程序就被编译和安装到了/sob/gmx2018.4目录下。修改用户目录下的.bashrc文件（比如运行gedit ~/.bashrc），在末尾加入source /sob/gmx2018.4/bin/GMXRC，然后保存。
实际上，FFTW库可以不必手动安装，因为可以在安装GROMACS时自动下载并安装FFTW库。但由于国内链接FFTW官网服务器往往比较慢，自动下载FFTW库可能中途卡住或者过程巨慢，因此还是建议手动方式安装FFTW库。如果你确实打算自动安装FFTW库的话，将上文第2节直接忽略掉，也不必设export CMAKE_PREFIX_PATH=/sob/fftw338，把第3节的cmake3那一步额外加上-DGMX_BUILD_OWN_FFTW=ON选项即可，这样编译GROMACS时就会自动在FFTW官网上下载FFTW包并自动编译之。
一般计算只需要按照前述编译单精度版本就够了，但如果模拟刚开始就崩溃，有时候用双精度版本可解决，但计算比单精度版慢将近一倍、trr/edr等文件体积大一倍。另外，做正则振动分析时在能量极小化和对角化Hessian矩阵的时候一般也需要用双精度版以确保数值精度。注意，编译双精度版本时不支持GPU加速。
要编译双精度版本的话，先按照前文方式编译一遍单精度版本，毕竟这之后在研究中肯定也得用。然后再重复一遍编译过程，但是在编译FFTW时去掉--enable-float，并且在使用cmake3命令时额外加上-DGMX_DOUBLE=ON选项。双精度版本的GROMACS可执行文件是gmx_d，而单精度是gmx，因此单精度和双精度的可执行文件可以同时存在于同一目录，互不冲突。
顺带一提，笔者在答疑时经常看到有人明明用的是单机并行，却非要装个MPI版GROMACS，这需要批评。因为这不仅需要多做一步，而且比起用默认方式基于thread-MPI和OpenMP的并行方式效率还更低，因此单机并行装MPI版完全是自取其辱。
important link to gromacs
If you have prmtop and inpcrd（prmcrd/mdrest） files for amber, just convert them to gromacs formate with the following script
./use-SOL
use-SOL
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 nrexclCV9 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.078000Or; 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
[ moleculetype ]; Name nrexclCV9 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
submit script
6 trajactory files (trr)这姑且成为轨道文件，读取trr文件对我来说走了很多弯路，曾经试图尝试使用guassview和vmd1.8.4等读图工具读取。后来跑到smth得知 读取trr在gromacs中有一种自带程序 ngmx。
gro file to pdb file
open trr file with pymol
title = Protein-ligand complex MD simulation
; Run parameters
integrator = md ; #指定积分算法；md:蛙跳牛顿积分算法，用于平衡动力学积分
nsteps = 500000 ; #积分或能量最小化步数
dt = 0.002 ; #积分步长
; Output control
nstxout = 0 ; #坐标保存到轨迹文件的频率
nstvout = 0 ; #速度保存到轨迹文件的频率
nstenergy = 5000 ; #能量保存到轨迹文件的频率，必须是nstcalcenergy的倍数
nstlog = 5000 ; # log文件更新频率
nstxout-compressed = 5000 ; #
compressed-x-grps = System
energygrps = Protein JZ4 #保存能量的组
; Bond parameters
continuation = yes ; #初试构象不约束，不复位，用于精确的继续计算或重计算
constraint_algorithm = lincs ; #约束算法 lincs ：不能用于角度约束
constraints = all-bonds ; #键约束，all-bonds:所有键约束
lincs_iter = 1 ; #迭代次数，用于LINCS约束精度，默认1
lincs_order = 4 ; #约束偶合矩阵阶次，用于LINCS约束精度，默认4
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid ; #邻近列表搜索方法
nstlist = 10 ; #邻近列表更新频率
rcoulomb = 1.4 ; #短程库伦截断
rvdw = 1.4 ; #短程范德华力截断
; Electrostatics
coulombtype = PME ; #库伦计算方式
pme_order = 4 ; #PME插值，默认4表示3次插值
fourierspacing = 0.16 ; #FFT傅里叶变换格点间距，默认。。。，与PME同时使用
; Temperature coupling
tcoupl = V-rescale ; #指定热耦合方法
tc-grps = Protein_JZ4 Water_and_ions ; #热偶合组
ref_t = 300 300 ; #参考温度--恒温值，个数对应组
; Pressure coupling
pcoupl = Parrinello-Rahman ; #指定压力耦合方式 ；Parrinello-Rahman：
pcoupltype = isotropic ; #isotropic:盒子各向同性
ref_p = 1.0 ; #参考压力---恒压值，一般为1 bar
compressibility = 4.5e-5 ; #水可压缩性，1 bar300K时为4.5e-5 bar-1
; Periodic boundary conditions
pbc = xyz ; #周期性边界条件；xyz:使用周期性边界条件
; Dispersion correction
DispCorr = EnerPres ; #色散校正
; Velocity generation
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)
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), 会使用对立方盒子进行角截断的方法生成截角八面体, 为转换这种截角八面体文件, 可使用以下命令:
gmx editconf -f in -rotate 0 45 35.264 -bt o -box veclen -o out
其中veclen是立方盒子大小乘以sqrt(3)/2.
在使用pdb2gmx创建了模拟分子体系之后, 可以使用editconf为你的分子创建一个模拟盒子, 也可以认为是使用editconf将分子放进一个盒子中. 这样, 你就可以往盒子里面添加水分子, 离子或者其他溶剂等等了.
-princ这个选项可以用来对齐分子, 比如使分子沿X轴对齐. 例如, 你想将分子中的两个残基沿Y轴对齐, 那么就在索引文件中将这俩个残基标记一下, 然后使用-princ, 根据提示就能对齐分子了.
你需要先了解一下用GROMACS做分子动力学模拟的一般过程:
本示例教程将引导GROMACS新用户进行一次模拟, 模拟的体系是水盒子中的蛋白质(溶菌酶, lysozyme)及离子. 我们会解释每一步中用到的输入文件以及得到的输出文件. 这些输入文件中所用的设置都很典型, 因此, 你在进行模拟时一般也可以采用这些设置.
本教程假定你正在使用5.x版本的GROMACS.
GROMACS基础
随着5.0版本GROMACS的发布, 所有的工具程序现在都成为了gmx程序中的一个模块. 而在以前版本的GROMACS中, 这些工具都是可以单独使用的, 并且有自己的命令. 在5.0版本中, 仍然可以通过符号链接使用这些命令, 但将来的版本会废弃这种作法. 因此, 你最好还是习惯这种新的使用方式. 要得到GROMACS某一模块的帮助信息, 你可以使用下面命令中的任何一个
gmx help (module)
gmx (module) -h
使用时只要将其中的(module)替换为要查询的命令的实际名称即可. 相关信息会输出到终端, 其中包括可用的算法, 选项, 需要的文件格式, 已知的缺陷和限制, 等等. 对GROMACS的新用户来说, 查看常用命令的帮助信息是了解每个命令功能的有效方式.
溶菌酶教程
你会看到下面的网页
网页上包含了很多与鸡蛋白溶菌酶有关的信息, 可以帮助你更好地了解这个蛋白质. 我们需要的是这个蛋白质的晶体结构文件, 因此点击右边的Download Files, 然后下载PDB格式的晶体结构文件.
得到结构文件之后, 你可以使用VMD, Chimera, PyMOL等可视化程序来查看一下这个蛋白质的结构. 如果还不熟悉这些程序的使用, 你可以参考网上的教程以及示例视频. 在查看轨迹方面, VMD功能非常强大, 而使用PyMOL你可以轻松地得到高质量的蛋白质结构图片. 本教程首页的蛋白质结构图片就是利用PyMOL得到的. 如果你不是很喜欢VMD和PyMOL程序的操作模式, 你或许可以试试Jmol. 此外, 还有非常多的分子结构可视化程序. 如果你只是用于查看结构, 选择你喜欢的那个即可. 上面提到的这几个程序都支持脚本, 利用脚本, 你可以更好地操控分子, 并能进行一些处理.
下面是VMD中各种显示模式下的蛋白质三维结构图
查看过这个蛋白质分子之后, 你要去掉晶体结构中的结晶水. 请使用普通的文本编辑器, 如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模块, 它的作用的是产生三个文件:
使用如下命令执行pdb2gmx程序:
pdb2gmx将处理结构, 输出一些相关信息后, 提示你选择一个力场:
力场包含了将要写入到拓扑文件中的信息. 这个选择非常重要! 你必须仔细了解每个力场, 并决定哪个力场最适用于你的体系. 在本教程中, 我们选用全原子OPLS力场, 因此在命令提示行中输入15, 然后回车.
pdb2gmx程序还接受很多其他选项, 下面列出常用的几个:
最后的注意事项: 许多用户认为.gro文件是必需的. 事实并非如此. GROMACS可处理很多不同的文件类型, .gro不过是一些程序输出坐标文件时所用的默认格式. 这种格式非常紧凑, 但精度有限. 如果你更愿意使用其他格式, 如PDB格式, 为你的输出文件指定合适的文件名称, 并使用.pdb作为扩展名即可. 使用pdb2gmx程序的目的在于生成与力场兼容的拓扑文件, 输出的结构不过是此目的的副产品, 只是为了方便用户. 输出结构的格式可以任意选择(参看GROMACS手册上对不同格式的说明).
此行调用了OPLS-AA力场的参数, 它位于文件的开头, 这意味着下面的所有参数都来自OPLS-AA力场. 下一重要行是[ moleculetype ], 后面是
“Protein_A”定义了分子名称, 这是因为这个蛋白质在PDB文件中被标定为A链. 对键合近邻的排除数为3. 关于排除的更多信息可从GROMACS手册上找到. 对此信息的讨论超出了本教程的范围.
下一节定义了蛋白质中的[ atoms ], 信息按列给出:
这些信息的解释如下:
下面几节包括[ bonds ], [ pairs ], [ angles ]和[ dihedrals ]. 其中一些无需解释(键, 键角, 二面角). 这些节中的参数和函数类型可参看GROMACS手册的第5章. 特殊的1–4相互作用包含于“pairs”(GROMACS手册5.3.4节).
正如你看到的, 通过使用值为1000 kJ mol-1 nm-2的力常数(kpr), 也可以对水分子进行位置限制.
接下来包含了离子的参数:
最后是体系级别的定义. [ system ]指令给出了体系的名称, 在模拟中此名称将被写入到输出文件中. [ molecules ]指令列出了体系中的所有分子.
[ molecule ]指令的几个关键注意点:
任何时候, 如果不能满足这些明确的要求, 运行grompp程序(后面介绍)时都会产生致命错误, 如mismatched names, molecules not being found或其他各种错误.
现在我们已经检查了拓扑文件的内容, 可以继续构建体系了.
现在你已经熟悉了GROMACS的拓扑文件, 可以继续创建体系了. 在本例中, 我们将要模拟一个简单的水溶液体系. 我们也可以模拟处于其他不同溶剂中的蛋白质或其他分子, 只要涉及到的物种有合适的力场参数.
定义一个模拟用的盒子并添加溶剂要分两步完成:
现在你需要决定使用哪种晶胞. 对于本教程的目的而言, 我们将使用一个简单的立方盒子作为晶胞. 当你对周期性的边界条件与盒子类型有了更多了解后, 我强烈推荐你使用菱形十二面体晶胞, 因为在周期性距离相同的情况下, 它的体积大约只有立方体晶胞的71%, 因此可以减少需要加入的溶剂水分子的数目.
让我们使用editconf来定义盒子:
上面的命令将蛋白质置于盒子的中心(-c), 并且它到盒子边缘的距离至少为1.0 nm(-d 1.0). 盒子类型是立方体(-bt cubic). 到盒子边缘的距离是一个重要参数. 因为我们要使用周期性边界条件, 必须满足最小映象约定, 也就是说, 一个蛋白质永远不能“看到”它自身的周期性映象(不能与其自身有相互作用), 否则计算的力就会含有虚假的部分. 指定溶质与盒子之间的距离为1.0 nm意味着, 蛋白质分子的任意两个周期性映象之间的距离至少是2.0 nm. 对于模拟中常用的任何截断方案, 这个距离都足够了.
现在我们已经定义好了模拟盒子, 可以用溶剂(水)填充它了. 使用solvate模块添加溶剂:
solvate记录了增加的水分子数目, 并将其写入拓扑文件中以显示它所做的更改. 注意, 如果你使用其他的(非水)溶剂, solvate不会在拓扑文件中写入这些信息! 它自动记录更新水分子的功能是直接写在源代码中的.
GROMACS中添加离子的工具是genion. genion的功能是读取拓扑信息, 然后将体系中的一些水分子替换为指定的离子. genion需要的输入文件称为运行输入文件, 扩展名为.tpr. 这个文件可使用GROMACS的grompp(GROMacs Pre-Processor)模块产生, 而且后面我们运行模拟时也会用它. grompp的功能是处理坐标文件和拓扑(它描述了分子)以产生原子级别的输入文件(.tpr). .tpr文件包含了体系中所有原子的所有参数.
为了用grompp产生.tpr文件, 我们还需要一个扩展名为.mdp(molecular dynamics parameter)的输入文件. grompp会将坐标和拓扑信息与.mdp文件中设定的参数组合起来生成.tpr`文件.
实际上, 在这个步骤中所用的.mdp文件中可使用任何合理的参数. 我通常会使用能量最小化的参数设置, 因为它非常简单而且不涉及任何复杂的参数组合. 请注意 本教程中所用的文件可能 只 适用于OPLS-AA力场. 其他力场的参数设置, 特别是非键参数设置可能很不一样.
使用下面的命令来产生.tpr文件
现在我们得到了一个二进制的.tpr文件, 它提供了我们体系的原子级别的描述. 将此文件用于genion:
出现提示后, 选择13 “SOL”. 这意味这我们将用离子替换一些溶剂分子. 你肯定不想用离子去替换蛋白质的某一部分.
在genion命令中, 我们以结构/状态文件(-s)作为输入, 以.gro文件作为输出(-o), 以拓扑文件(-p)来反映去除的水分子和增加的离子, 并且定义了阳离子和阴离子的名称(分别使用-pname和-nname), 告诉genion只需要添加必要的离子来中和蛋白质所带的净电荷, 添加的阴离子数目为8(-nn 8). 对于genion, 除了简单地中和体系所带的净电荷以外, 你也可以同时指定-neutral和-conc选项来添加指定浓度的离子. 关于如何使用这些选项, 请参考genion的说明.
在以前版本的Gromacs中, 使用-pname和-nname指定的离子名称由力场决定, 但从4.5版本开始就完全标准化了. 指定的离子名称始终是大写的元素符号, 与[ moleculetype ]中的名称一致, 并会写入拓扑文件. 残基名称或原子名称可能会带有电荷符号(+/-), 也可能不带, 取决于力场. 不要在genion命令中使用原子名称或残基名称, 否则在下面的步骤中会导致错误.
[ molecules ]指令现在看起来应该这样:
使用分子结构可视化软件查看一下现在的体系
现在, 我们已经添加了溶剂分子和离子, 得到了一个电中性的体系. 在开始动力学模拟之前, 我们必须保证体系的结构正常, 原子之间的距离不会过近, 几何构型合理. 对结构进行弛豫可以达到这些要求, 这个过程称为能量最小化(EM, energy minimization).
能量最小化过程与添加离子过程差不多. 我们要再次使用grompp将结构, 拓扑和模拟参数写入一个二进制的输入文件中(.tpr), 但这次我们不需要将.tpr文件传递给genion, 而是使用GROMACS MD引擎的mdrun模块来进行能量最小化.
使用grompp处理这个参数文件, 以便得到二进制的输入文件:
现在我们可以调用mdrun来进行能量最小化了:
gmx mdrun -v -deffnm em
现在我们的体系已经处于能量最小点了, 可以用它进行真正的动力学模拟了.
EM可保证我们的初始结构在几何构型和溶剂分子取向等方面都合理. 为了开始真正的动力学模拟, 我们必须对蛋白质周围的溶剂和离子进行平衡. 如果我们在这时就尝试进行非限制的动力学模拟, 体系可能会崩溃. 原因在于我们基本上只是优化了溶剂分子自身, 而没有考虑溶质. 我们需要将体系置于设定的模拟温度下, 以确定溶质(蛋白质)的合理取向. 达到正确的温度(基于动能)之后, 我们要对体系施加压力直到它达到合适的密度.
我们将使用grompp和mdrun, 像在能量最小化过程中做的一样
除注释外, 所用参数的完整解释可以在GROMACS手册中找到. 注意.mdp文件中下面的这几个参数:
让我们来分析温度变化情况, 再次使用energy模块:
提示时输入15 0来选择体系温度并退出. 得到的结果应该和下面的差不多:
从上图可以清楚地看出, 体系的温度很快就达到了目标温度(300 K), 并在平衡过程中后面的时间内保持稳定. 对于这个体系, 更短的平衡时间(50 ps)也足够了.
前一步的NVT平衡稳定了体系的温度. 在采集数据之前, 我们还需要稳定体系的压力(因此还包括密度). 压力平衡是在NPT系综下进行的, 其中粒子数, 压力和温度都保持不变. 这个系综也被称为等温等压系综, 最接近实验条件.
该文件与NVT平衡时所用的参数文件没有太大不同. 注意添加的压力耦合部分, 其中使用了Parrinello-Rahman控压器.
其他几项改动如下:
我们使用grompp和mdrun, 像在NVT平衡所做的那样. 注意, 我们现在要使用-t选项以包括NVT平衡过程中的产生的检查点文件. 这个文件包含了继续模拟所需要的所有状态变量. 为使用NVT过程中得到的速度我们必须包含这个文件. 坐标文件(-c)是NVT模拟的最终输出文件.
gmx mdrun -deffnm npt
让我们来分析压力变化情况, 再次使用energy模块:
提示时输入16 0来选择体系压力并退出. 结果应与下图类似:
在100 ps的平衡过程中压力值涨落很大, 这并不意外. 图中的红线为数据的移动平均值. 在整个平衡过程中, 压力的平均值为1.05 bar.
让我们再来看看密度, 使用energy模块并在提示时输入22 0
跟压力一样, 红线是密度的移动平均值. 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中的一样.
依次运行下面的命令:
gmx mdrun -deffnm md_0_1
在GPU上运行GROMACS
假定你有一个可用的GPU, 要利用它可使用下面的mdrun命令:
gmx mdrun -deffnm md_0_1 -nb gpu
如果可用的GPU卡超过一个, 或者想利用GROMACS支持的杂合并行方案对计算进行划分, 请参考GROMACS手册以及网络上的资料. 这些技术细节超出了本教程的范围.
现在已经完成了对蛋白质的模拟, 我们应该来分析一下我们的体系. 哪些类型的数据才是重要的呢? 这是在模拟前就要思考的一个重要问题, 所以你应该对自己的体系需要采集哪些数据类型有自己的想法. 在本教程中, 我们只介绍一些基本工具.
第一个模块是trjconv, 这是一个后处理工具, 用于处理坐标, 修正周期性或手动调整轨迹(时间单位, 帧频率等). 在本教程中, 我们要使用trjconv来处理体系中的任何周期性. 蛋白质在单元晶胞中扩散, 可能看起来会在盒子两边之间进行“跳跃”. 我们使用下面的命令来处理这种情况:
选择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%完美的.
我们也可以将初始构型与模拟后的构型进行比较, 这样可以更直观地看出二者的区别.
去掉水分子可以看得更清楚一些
如果将两者进行最小二乘叠合更容易看出区别
蛋白质的回旋半径Rg可衡量其密实度. 如果蛋白质的折叠很稳定, 其Rg将保持一个相对稳定的值. 如果蛋白质去折叠, 它的Rg将随时间变化. 我们来分析一下模拟的溶菌酶的回旋半径:
可以看到, Rg值基本不变, 这预示着在温度为300 K时, 1 ns的时间内蛋白质很稳定, 处于紧密(折叠)的形式. 这一结果并非意外, 但说明了GROMACS具有先进的分析功能.
现在你已经用GROMACS完成了一个分子动力学模拟过程, 并分析了一些结果, 本教程并不全面. 你还可以用GROMACS完成更多类型的模拟(自由能计算, 非平衡MD, 简正模式分析, 等等). 你应该阅读一些文献以及GROMACS手册, 试着调整这里提供的.mdp文件中的参数来以便使模拟更有效, 更精确.
除非指定了选项, 否则不会使用可选文件. 非可选文件与此不同, 当不指定选项时, 会使用默认的文件名.
所有GROMACS程序接受的文件选项都可以忽略扩展名或文件名. 在这种情况下, 将会使用默认的文件名. 可使用多种输入文件类型, 如通用结构格式时, 将会搜索目录下具有指定名称或默认名称的每种类型的文件. 若没有发现这样的文件, 或使用输出文件时, 将会使用第一种文件类型.
除mdrun, trjcat和eneconv外, 所有GROMACS程序都会坚持命令行选项是否合理. 如果不合理, 程序将会中断执行.
所有GROMACS程序都有4个隐藏的选项:
枚举选项(enum)应使用选项说明中列出的变量之一, 变量可以缩略. 将会选择使用第一个与列表中最短变量匹配的变量.
向量选项可以使用1或3个参数. 当只提供一个参数时, 其他两个参数使用相同的值.
所有GROMACS程序都可以读取压缩的或gzip压缩的文件. 当读取压缩的.xtc, .trr和.trj文件时, 可以会有问题, 但这些文件无论如何也不能很好地进行压缩.
大多数GROMACS程序都可以处理原子数比输入运行文件或结构文件中的原子数少的轨迹, 但轨迹只能包含输入运行或结构文件中的头n个原子.