Jing's Notes
0%

【GROMACS】4.GROMACS中的命令

发表于 更新于 分类于 ，
Valine： 本文字数： 15k 阅读时长 ≈ 14 分钟
Gromacs的命令 `gmx xxx`

总体

常用的参数

  • -f：结构文件（gro、pdb、tpr等)或轨迹文件(tr、xtc等），或xvg文件，或edr文件。
  • -s、-c：结构文件 （gro、pdb、tpr等)。具体视程序而定
  • -n：索引文件(ndx）
  • -b、-e：设定轨迹中被考虑的部分是从哪里开始（begin)到哪里结束(end)
  • -dt：轨迹每多长时间被读入一次

dump

用来把GROMACS的各种私有的二进制文件转化成可读形式输出，便于检查和提取数据。

  • -s: .tpr
  • -f: .xtc/.trr
  • -e: .edr
  • -cp:.cpt

例子：从info.out中将看到md.cpt中包含的原子坐标/速度、盒子尺寸、当前模拟时间等各种信息。

1
gmx dump-cp mid.cpt > info.out

从tpr中可以获取原子间相互作用参数信息，从edr中可以获取每一帧各种能量项，从trr中可以获取每一帧坐标信息、速度、受力信息(若储存了的话)。

editconf

1
2
3
4
5
6
7
8
9
10
gmx editconf [-f [<.gro/.g96/...>]] [-n [<.ndx>]] [-bf [<.dat>]]
             [-o [<.gro/.g96/...>]] [-mead [<.pqr>]] [-[no]w]
             [-[no]ndef] [-bt <enum>] [-box <vector>]
             [-angles <vector>] [-d <real>] [-[no]c]
             [-center <vector>] [-aligncenter <vector>]
             [-align <vector>] [-translate <vector>]
             [-rotate <vector>] [-[no]princ] [-scale <vector>]
             [-density <real>] [-[no]pbc] [-resnr <int>] [-[no]grasp]
             [-rvdw <real>] [-[no]sig56] [-[no]vdwread] [-[no]atom]
             [-[no]legend] [-label <string>] [-[no]conect]

用来设置模拟体系的盒子，也可以用来平移和旋转结构，或者转换文件格式。

  • -f：输入结构文件
  • -n：输入索引文件（可选）
  • -o：输出结构文件（可选）

设置盒子，会自动把体系放在盒子中央，除非使用了-noc

  • -bt：盒子类型，triclinic（三斜。默认)，cubic（立方）dodecahedron（梭形十二面体），octahedron（截角八面体）
  • -box <vec>：对triclinic盒子自定义盒子三个边长(a,b,c)，对其它类型盒子只需定义一个边长值
  • -angles α β γ：使用-box时同时自定义盒子向量 (bc,ac,ab) 之间的夹角
  • -d <real>：将分子边界向X/Y/Z方向延展此值来定义盒子

修改结构：

  • -translate <vector>：对体系在X/Y/Z方向进行平移的距离
  • -rotate <vector>：绕X, Y和Z轴的旋转角度, 单位为度
  • -scale <vector>：将原子坐标在三个方向乘以此数值(可达到比如调节密度的目的）
  • -[no]princ：使某个组的最长的主轴向X轴对齐。
  • -mead foo.pqr：如果读入的是.tpr文件，则可以输出.pqr文件，格式类似pdb，但同时记录了原子电荷和原子半径信息

energy

用于从能量文件(.edr)中提取动能、势能的各种成分、温度、盒子尺寸、压力、维里值、体积、密度、能量组之间的非键作用能等数据，输出随时间的变化，并进行统计分析得到平均值、RMSD、漂移量，并可以计算一些体系的热力学性质。

  • -f:.edr文件

  • -o:.xvg文件(默认为energy.xvg)

  • -dp:在xvg文件中以更高精度输出数据

    运行比如gmx energy -f md.edr，然后将屏幕上感兴趣的量的序号依次输入;每输入一个按一次回车，最后以空行或0结束，则被选中的量随时间变化就会被提取并一起记录到energy.xvg中。使用xmgrace -nxy energy.xvg即可作到一张图上(若只有一套Y数据则不需要-nxy)

  • -fluct_props:计算与性质波动有关的量，诸如热容、等温可压缩系数、热膨胀系数

  • -nmol <int>:体系中的分子数，默认为1。对于计算热力学量时应当设为实际值，否则计算出的结果可能错误
    只有选择全了相应条目的Energy terms needed下面的项，才会在提取数据后将相应的量也算出来。

  • -fee：选项会计算体系与理想气体状态时的自由能差值

    $\Delta A = A(N,V,T) - A_{idealgas}(N,V,T) = kT \ln(<\exp(U_{pot}/kT)>)$

    $\Delta G = G(N,p,T) - G_{idealgas}(N,p,T) = kT \ln()$

    在界面里应当选择Potential 。在输出信息中的 -kT ln列下方的数值即是要取的值

  • -f2：提供另一个edr文件，用于计算两种参数下体系的自由能差。两个edr文件里的帧数需相同，应来自同一体系。

    $\Delta G = -kTln(_A),$

    界面里应选择Total Energy。

  • -fetemp <real> :估算自由能时用的参考温度，默认为300 K

genion

将结构中指定的组(一般是溶剂组)的部分分子随机替换为单原子离子。

  • -s：.tpr文件
  • -n：索引文件(可选)
  • -p：拓扑文件(会对之进行恰当修改)
  • -o：结构文件

用法：默认为加入Na+、Cl-，若自设离子可用-pname [阳离子名] -pq [阳离子电荷数] -nname [阴离子名] -nq [阴离子电荷数]，离子名均为大写，如NA、CL、CA

  • -np 阳离子数 -nn 阴离子数
  • -conc 盐浓度(mol/L)
  • -neutral 加入足够离子使体系电荷中性化

genion如果要求加入的离子数非常多的话，可能会出现No more replaceable solvent的提示导致运行失败：

  • 可以被genion替换的溶剂分子是距离离子超过rmin的溶剂分子。当被替换的溶剂分子已经很多，没有其它满足条件的可被替换的溶剂分子时，就会出现上述错误。
  • 解决方法是加上-rmin [离子间最小距离]，将离子间最小距离设成小于默认的0.6 nm的值。但也不能设得太小，否则可能导致带电相同的离子间出现位置太近，由于静电互斥过强导致模拟一开始崩溃。

genion命令无法加入多原子离子。可以利用packmol构建这种体系，也可以用inisert-molecules命令通过-replace选项对原有分子进行替换。

例子：

  • 产生含有5个Ga2+和10个NO3-的3×3×3nm水盒子

    1
    2
    3
    gmx solvate -box 3 3 3 -o waterbox.gro
    gmx insert-molecules -f waterbox.gro -ci NO3-.pdb -nmol 10 -replace -o salt.gro
    gmx insert-molecules -f salt.gro -ci Ca2+.pdb -nmol 5 -replace -o salt.gro

    每次都选择water组

grompp

grompp(GROMACS preprocessor)，读取运行参数文件、拓扑文件和结构文件，将信息处理、整合，产生mdrun的输入文件.tpr。

  • -f:.mdp文件
  • -c:结构文件
  • -p：拓扑文件
  • -t(非必需)：.trr轨迹文件(默认读取其最后一帧信息)或.cpt文件
    • 精确延续之前的体系状态应当用-t提供.cpt文件
  • -r：位置限制文件(使用位置限制时才需要)
  • -n(非必须)：索引文件
  • -o：.tpr文件
  • -pp选项可以输出一个processed.top,是预处理器处理top文件后的结果，所有include的文件都会被展开，对ga_1之类参数的引用也直接被实际的参数所替换。
  • -maxwarn [数目]：可以设允许最多出现多少个Warmning而不终止运行。默认只要出现Wamning则grompp就会终止。
    • grompp给出的所有Warning都应当留意，往往暗示出输入文件里的错误，若轻易无视之可能模拟崩溃或者得到无意义的结果。

grompp还会自动产生mdout.mdp文件，是参数完整的.mdp文件，会体现没设定的参数的默认值，并包含许多注释。

grompp是按照top文件里的[ molecules ]记录的分子顺序展开成的原子顺序构建tpr文件的，从结构文件中只是读取坐标。grompp会按顺序检查top和结构文件里的原子名的对应关系，如果不符合将会给出警告并使用top里的原子名。

常见Warning

关于净电荷的常见Warning

  • You are using Ewald electrostatics in a system withnet charge. This can lead to severe artifacts, suchas ions moving into regions with low dielectric, dueto the uniform background charge. We suggest toneutralize your system with counter ions, possibly incombination with a physiological salt concentration.
    • 这说明的当前用的是Ewald、PME之类考虑无穷远静电相互作用的方法，但是体系的净电荷不为0，此时模拟可能明显不合理，可能需要修正。
    • 此时应当检查grompp提示的System has non-zerototal charge后面显示的净电荷是多少，有三种情况：
      • 数值非常接近于0：如0.0013。可以用-maxwarn无视此警告继续产生tpr
      • 数值是不为0的整数：如3.0、-1.0。说明需要用诸如gmx genion加抗衡离子使得整个体系电中性化
      • 数值是一个偏离整数明显的数：如0.167。这必定是因为有的moleculetype里[ atoms ]定义的所有原子电荷加不为整数所致。哪怕每一个分子的净电荷偏离0仅有诸如0.002的程度，但分子数很多时也会造成整体净电荷偏离0非常明显。需要自行检查体系中涉及的各个moleculetype的净电荷(手动加和其中的原子电荷)，对于净电荷轻微偏离0的分子，可手动调整其中某一个或多个原子电荷以令分子的净电荷恰为整数。

关于原子顺序的常见Warning

  • 结构文件里的原子顺序必须和拓扑文件里的[ molecules ]及各类分子的[ atoms ]所对应的原子顺序完全一致，否则模拟时参数会错乱。而拓扑文件和结构文件中的原子名、残基名可以相同也可以不同，不同时在用grompp时会出现warning，可以无视，此时mdrun最终产生的gro文件里的原子和残基名将与拓扑文件里的一致。

insert-molecules

将指定的分子插入体系空挡中，或构建一个填充指定分子的盒子。与solvate命令关键不同之处在于此命令中被填入的分子只需提供单个分子的结构文件。

  • -f：被填充的结构文件
  • -ci：要填入的分子的结构文件
  • -o：填充后的结构文件，默认为out.gro
  • -box <vector>：构建新盒子并往里填充
    • 如果-f和-box一起用，则是把-f文件里的盒子改成-box的尺寸后再往空隙里填。

随机插入

随机插入，此为默认情况。和solvate命令一样原子间距离肯定不会小于它们之间范德华半径之和(且可以用-scale调节半径的倍数)

  • -nmol <int>：最多填入的分子数

  • -try <int>：尝试插入-nmol乘以-try次。因此此值越大可能填充进去的分子越多。默认为10

  • -rot：默认为xyz，即填充时分子在三个方向可随意旋转。也可以为z(只允许绕Z轴旋转)、none(不能旋转)

例子

  • 将formamide.pdb里的分子填入3×4×4 nm的空盒子，输出solv.gro。尝试最多填入1000个，如果填不下则能填多少填多少，尽可能填满

    1
    gmx insert-molecules -box 3 4 4 -ci formamide.pdb -o solv.gro -nmol 1000

  • 将ethanol.gro里的分子填入protein.pdb当中溶质与盒子间的空隙，尝试最多填入200个

    1
    gmx insert-molecules -f protein.pdb -ci ethariol.gro -o solv.gro -nmol 200

  • 把solv.gro的盒子(假设原先为2×2×2nm) 改为2×2×5nm，往多出来的2×2×3nm空间内填充CH4.pdb，尝试最多填入1000个(此法可以构建两相界面体系)

    1
    gmx insert-molecules -f solv.gro -ci CH4.pdb -box 2 2 5 -o biphase.gro -nmol 1000

指定位置插入

在一批指定的位置插入分子，需自己提供一个.dat后缀的文件，每一行指定个分子插入的X，Y，Z坐标（相对于-ci文件里的坐标原点而言），用-ip选项接上此文件的路径。

1
gmx insert-molecules -box 1 2 2 -ci formamide.pdb -o new.gro -ip pos.dat

pos.dat的内容类似于：

1
2
3
4
0.25 0.5.0.5
8.25 0.5 1.0
0.25 1.8 0.5
0.25 1.0 1.0
  • 通常-ci指定的结构文件里的分子的中心位置应当处于坐标原点，以令.dat设的恰是各个分子中心出现的位置

mdrun

是GROMACS最关键的命令，用于做各种计算模拟任务。

  • -s:.tpr文件
  • -cpi:.cpt文件(可选，用于精确续算目的。若提供则完整的状态信息从这里读，否则用.tpr里的信息)
    • 运行时每隔一定时间(用-cpt来设，默认15分钟)往.cpt里写一次完整的状态数据用于续算。前一次的.cpt文件会被改名为_prev结尾
  • -rerun:轨迹文件(可选，计算已有轨迹的能量才用)
  • -o:输出的.trr文件
  • -x:输出的.xtc文件(可选)
  • -cpo:输出的.cpt文件(可选)
  • -c：输出的结构文件
  • -e:输出的.edr文件
  • -g：输出的log文件
  • -deffnm：设定对所有选项都用的文件名
  • -v:在终端不断输出当前已跑步数以及预计的完成时间

例子

  • 常规计算：将运行md1.tpr,输出md1.trr/xtc/log/cpt/edr/gro。gro文件应轨迹最后一帧结构而且是分子保留完整的状态。运行进度预计完成时间不断输出到屏幕上。四个线程将被利用。

    1
    gmx mdrun -v -deffnm md1 -nt 4 -pin on
    • 如果要在远程服务器上运行，用比如：nohup gmx mdrun -deffnm md1 nt 36 -pin on &，之后若用exit命令断开与服务器的连接，任务在服务器上会一直跑完。运行进度可以用tail -f md1.log来实时监控。

  • 跑完意外中断的任务md1，轨迹、能量、日志文件会接续之前文件断点处继续写入，最终得到的这些文件和一次性完整运行的相同。

    1
    gmx mdrun -v -deffnm md1 -cpi md1.cpt

  • md1已完整跑完，延续之前的模拟参数，再跑额外的5ns模拟，新任务叫md2

    1
    2
    gmx convert-tpr -s md1.tpr -extend 5000 -o md2.tpr
    gmx mdrun -v -deffnm md2 -cpi md1.cpt -noappend

    最终会得到md2.part0002.trr/xtc/log/edr/gro和md2.cpt

  • 同上，但是直接在之前md1的输出文件上续写

    1
    2
    gmx convert-tpr -s md1.tpr -extend 5000 -o md1.tpr
    gmx mdrun -v -deffnm md1 -cpi md1.cpt

    注:convert-tpr还可以结合-until xxx将模拟终止时间延长到xxx

  • heat任务已跑完，想延续其最后的状态结合其它的mdp设定跑prod任务

    1
    2
    gmx grompp -f prod.mdp -c heat.gro -t heat.cpt -p NANA.top -o prod.tpr
    gmx mdrun -v -deffnm prod
    • 这样产生的轨迹的初始时间由mdp里的tinit定义，如果没定义则默认从0开始计。
    • 注：多次延长模拟时间时，若为了让每一个任务都有独立的从而易于管理的文件，可以每一段延长的模拟都设置对应的mdp文件(如prod2.mdp、prod3.mdp…)并用此法产生相应的tpr文件。

  • 用UR.tpr里的参数和设定计算config.pdb结构的能量(即单点能)

    1
    gmx mdrun -deffnm UR -rerun config.pdb

  • 用SSR.tpr里的参数和设定计算traj.trr里每一帧的能量并写入SSR.edr中

    1
    gmx mdrun -deffnm SSR -rerun traj.trr
    • 使用-rerun计算某个结构文件或轨迹文件里各个结构的能量时，提供的结构/轨迹文件里的原子顺序必须与.tpr文件精确一致
    • 如果用-reprod，可以避免一些自动优化，使得在相同软硬件环境下(包括GROMACS版本)对相同tpr每次跑的结果精确相同。

并行机制

  • MPI：是最普的并行化技术。GROMACS4.5版以前的并行无论节点内还是节点间一律通过MPI实现
  • thread-MPI:是MPI的轻量级子集，专门用于节点内并行，从4.5版起开始被GROMACS内置
  • OpenMP:是最方便、流行的实现节点内并行的技术，从4.6版开始被GROMACS支持。日前版本在节点内并行时thread-MPI与OpenMP会自动结合使用，每个thread-MPI线程下属会有多个OpenMP线程，总并行线程数是二者的乘积。另外GROMACS的一些个别分析命令如hbond也支持OpenMP并行

基于MPI并行需要在编译时提供MPI库(如OpenMPI)，而基于thread-MPI和OpenMP在节点内并行在编译时不需要额外的库。因此，如果只需要节点内并行则无需提供MPI库并专门编译MPI版，而且用MPI版来并行计算速度还会更慢。

参数：

  • -nt：设定用几个线程并行。默认用所有CPU逻辑核数来并行。也可以用-ntmpi和-ntomp分别明确指定其中thread-MPI线程数和下属的OpenMP线程数，二者之积对应总线程数(-nt值)。
    • 当节点内n核并行时，top命令看到的gmx进程的CPU占用率原理上应接近n*100%
  • -pin on:将线程和CPU核心绑定，避免系统随意调度导致性能损失。
  • -pinoffset [id]:-pin on时绑定的CPU核心从第id个核开始算，第一个核心id为0。此选项默认为0。如果同时跑多个任务且用-pin on，总应当结合-pinoffset避免有的核心同时做两个任务
  • -pinstride [int]:-pin on时绑定的核的序号的间隔。通常设为1，使得被绑定的核的序号挨着

对于支持超线程的CPU一般建议设成-nt [CPU物理核数] -pin on。没有必要关闭超线程。
对于跑非常长、很耗时的轨迹，可以加上比如-nsteps 20000只跑20000步，测试一下以下情况的差异以寻找性能最佳的组合用来跑当前任务:

  • 用和不用-nt [物理核心数]
  • 用和不用-pin on
  • -ntmpi与-ntomp的不同组合(二者之积应对应实际调用的总核数)

报错：

  • Fatal error：Your choice of number of MPI ranks and amount of resourcesresults in using 36 OpenMP threads per rank, which is mostlikely inefficient. The optimum is usually between 1 and 6threads per rank. If you want to run with this setup,specify the -ntomp option. But we suggest to change thenumber of MPI ranks (option -ntmpi).

    • 写上-ntomp 6使得每个thread-MPI线程下属的OpenMP线程在程序建议的范围。OpenMP线程数太高的话性能会打折扣，将hread-MPI和OpenMP恰当搭配并行时性能才最好。

  • 模拟小体系时可能有以下提示 NOTE: Parallelization is limited by the smallnumber of atoms, only starting 2 thread-MPIranks, You can use the -nt and/or -ntmpi option to optimize the number of threads

    • 说明体系原子数太少，仅能利用2个thread-MP1线程进行计算，不用管。
    • 如果实际模拟时发现CPU利用率很低，可尝试用-ntmpi 2结合-ntomp [物理核心数/2]，尽可能充分靠OpenMP来并行。

GPU加速相关

当mdrun检测到有GPU可以利用时GPU版会自动利用GPU计算，当有些功能不支持GPU导致报错时可以明确指定用CPU算。

  • -bonded、-nb、-pme [auto/cpu/gpu]:指定成键、非键作用的实空间部分、PME算的非键作用的倒易空间部分用CPU还是GPU算。默认为auto
  • GPU版mdrun可以通过-nb cpu -pme cpu -bonded cpu强行要求只在CPU上算(2018版没-bonded选项)
  • -gpu_id [GPUid]:有多块GPU且想用来分别跑多个任务时，指定任务在哪个GPU上跑，避免冲突
  • -update gpu:令约束以及坐标更新也放到GPU上算，从而进一步加快速度，这称为GPU-resident模式，但功能上有限制。
    • 从2020版开始支持，从2023版开始默认开启。2023版以前默认为--update auto，总是在CPU上做。用export GMX_FORCE_UPDATE _DEFAULT_GPU=true设置环境变量可以总是强行在GPU上做。
    • 注:为了尽可能避免通信，GPU上做update需要满足updategroups条件。此时如果用了约束，几个原子若共同与某个原子有约束关系，那么在[ atoms ]中这儿个原子应当连续地排在那个原子的后头。例如用constraints = hbonds时，氢原子应当连续地排在它相连的重原子的后头。

例子

  • 16核32线程机子上用两个GPU分别跑两个任务，且各用不同的8个CPU物理核心的最佳的运行方式：任务1:由0-7号CPU核心连同0号GPU运行；任务2由8-15号CPU核心连同1号GPU运行

    1
    2
    gmx mdrun [命令] -nt 8 -pin on -pinstride 1 -gpu_id 0
    gmx mdrun [命令] -nt 8 -pin on -pinoffset 8 -pinstride 1 -gpu_id 1

pdb2gmx

pdb2gmx是GROMACS中模拟生物大分子体系的关键性工具。它载入生物大分子结构文件，从相关力场库文件中读取所需信息，给体系恰当加氢和重排原子顺序，最后产生生物分子的拓扑文件。

  • -f：结构文件，一般用pdb（输入）
    • pdb2gmx也适合用来产生其它类型的由单体聚合而成的体系，如聚合物，但需要自已恰当修改rtp文件后pdb2gmx才能处理。
    • pdb2gmx不是普性的拓扑文件产生工具，对rtp文件里没有定义的小分子，不可能用pdb2gmx产生拓扑文件
  • -o：整理后的结构文件，默认为conf.gro
  • -p：拓扑文件，默认为topol.top
  • -i：位置限制文件，默认为posre.itp
  • -ff <string>：设定所用力场，默认是交互式选择
  • -water <enum>：设定所用水模型，默认是交互式选择
  • -heavyh:把氢原子质量设为4(并从相连重原子上扣除相应质量)，使氢的运动减缓，从而能用更大步长
  • -ignh：忽略输入文件里所有氢。NMR测定的pdb里虽然有所有的氢，某些程序(如分子对接程序)产生的结构文件里也往往有氢，但原子名大概率不标准，应当用此选项让pdb2gmx按标准规则加氢
    • 通常建议始终带上-ignh
  • -ter：交互式选择如何处理末端氨基和羧基，默认是搞成带电状态。对AMBER力场不适用(因为它是直接按照r2b文件对残基处在氮端和碳端的情况进行设定)
  • -merge:默认是no，代表每个链作为单独的[ moleculetype ]。也可以用all都弄到个[ moleculetype ]里，或interactive由用户下面则

以下如果都要交互式选择，则用-inter选项：

  • -lys、-arg、-asp、-glu、-gin、-his：交互式选择这些残基的质子化态。默认是除了HIS都设为带电状态，HIS自动根据形成氢键最优方式判断是否两个氮都质子化(带电)或者中性时氢在哪个氮上
  • -ss:交互式方式选择如何设定二硫键

pdb2gmx利用的文件

（在力场文件夹中，pdb2gmx按照此顺序处理: r2b→hdb→ rtp→tdb）

  • aminoacids.r2b：

    • 含有将结构文件里的残基名转化到building block 名的对应关系。同一个残基在不间力场里的名字有时不同, 所以需要先做这个转换再套用 rtp

    • 例子：G54A7力场的 r2b文件：

      1
      2
      3
      4
      5
      6
      7
      ;GMX   Force-field
      CYS    CYSH
      HISD   HISA
      HISE   HISB
      LYS    LYSH
      LYSN   LYS
      HEM    HEME

    • 例子：AMBER力场比较特殊，第3 、 4 列说明如果残 基在氮端和碳端时分别改名成什么。第5列是对于单个两亲氨基酸的悄况，当前未定义：

      1
      2
      3
      4
      5
      ; rtp residue to rtp building block table

      ;      main  N-ter C-ter 2-ter
      ALA    ALA   NALA  CALA  -
      ARG    ARG   NARG  CARG  -

  • aminoacids.hdb：hydrogen database 。设定氢原子 如何连接 到重原子上,用于自动补氢

    • 例子：

      1
      2
      3
      4
      5
      6
      7
      ALA     1       
      1	1	H	N	-C	CA
      ARG     4  →下面要定义的项数     
      1	1	H	N	-C	CA
      1	1	HE	NE	CD	CZ
      2	3	HH1	NH1	CZ	NE
      2	3	HH2	NH2	CZ	NE
      • 第一列：这一项要加氢的数目
      • 第二列：加氢的规则
      • 第三列：加入的氢原子名。此项加两个氢，因此名称将为HH21和HH22
      • 第4-6列：加氢时候用于定位几何位置的原子

  • aminoacids.rtp,：residue topology parameter。记录残基(以 及水、个别小分子、离子)的拓扑信息, pdb2gmx载入结构文 件时会依照此文件构建拓扑信息。

    • 例子：AMBER力场中 rtp文件的一部分

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      [ bondedtypes ] ;设定成键项的函数类型以 及产生拓扑信息的规则
      ; bonds  angles  dihedrals  impropers all_dihedrals nrexcl HH14 RemoveDih
           1       1          9          4        1         3      1     0
      ...
      [ GLY ] ; HAx atoms assigned new ff03 atom type
       [ atoms ] 
       ; 原子名 原子类型		  原子电荷     序号
           N    N           -0.374282    1
           H    H            0.253981    2
      ...
       [ bonds ]
           N     H
           N    CA
      ...
           C     O
          -C     N  ;原子名带负号、正号分别是指上一个、下一个残基的原子
       [ impropers ]
          -C    CA     N     H
          CA    +N     C     O

    • pdb里原了顺序与 rtp里可以 不同, pdb2gmx会按照原子 名来匹配。

    • pdb2gmx会根据 rtp里设定的成键关系,自动产生所有该 有的[ angles ] 、[ dihedrals ] 、 [ pairs ] 项。
      • rtp里若直接定义了[ angles ]等项,会稷盖自动产生的。 GROMOS 力场的 rtp直接就定义完整了所有成键项。

  • aminoacids.n.tdb和aminoacids.c.tdb：termini database, 是氨基酸的氮端和碳端的末端数据库,设定对末端残基如何删除、 添加、替换原子和连键

    • 注：Amber力场tdb文件为空,是因为它在 rtp 里就对末端残基直接定义了不同类 型,不需要再按照tdb的规则对末端残基在普通残基基础上做额外处理。

    • 例子：此例说明把 N原子的原子类型改为 NL, CA原子的原子类型改为CH1。 按照4号规则,以 N 、 CA、 C为参考 加入3个H类型氢原子,分别名为 H1 、 H2 、 H3 。原先残基中名为 H的原子 则被删除。然后设定新加入的原子 与其它原子间的成键关系。

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      ; G54A7力场的氮端tdb文件的一部分 aminoacids.n.tdb
      [ None ]

      [ NH3+ ]
      [ replace ]
      N		NL	14.0067	0.129
      CA		CH1	13.019	0.127	0
      [ add ]
      3	4	H	N	CA	C	  ;H代表新加入的氢的名称前缀
      	H	1.008	0.248     ;新加入氢原子类型、质量、电荷
      [ delete ]
      H
      [ bonds ]
      N	H1	gb_2
      N	H2	gb_2
      N	H3	gb_2
      [ angles ]
      H1	N	H2	ga_10
      H2	N	H3	ga_10
      ...
      [ dihedrals ]
      H1	N	CA	C	gd_39
      ...

    • [ none ] 代表对于非末端残基 什么都不做修改

    • [ NH3+ ] 代表对于氮端的且质子化的 氨基酸要做哪些处理

    • 例子：G54A7力场的碳端tdb文件的一部分 aminoacids.c.tdb。此例说明把C原子的原子名重置为C。按照8号规则，以C、CA、N做为参考位置加入2个OM类型氧原子，分别名为O1和O2。但如果原本结构中就已经有O和OXT原子了，则直接将之改成名为O1和O2的OM类型原子。

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      [ COO- ]
      [ replace ] ;此处代表对于碳端的且羧基要设成带电状态的氨基酸要做哪些处理。
      C	C	C	12.011	0.27
      O	O1	OM	15.9994	-0.635
      OXT	O2	OM	15.9994	-0.635
      [ add ]
      2	8	O	C	CA	N
      	OM	15.9994	-0.635
      [ bonds ]			;然后设定新加入的原子与其它原子间的成键关系。
      C	O1	gb_6
      C	O2	gb_6
      [ angles ]
      O1	C	O2	ga_38
      CA	C	O1	ga_22
      CA	C	O2	ga_22
      [ dihedrals ]
      N	CA	C	O2	gd_45
      N       CA      C       O2      gd_42
      [ impropers ]
      C	CA	O2	O1	gi_1

  • aminoacids.vsd：virtual site database文件,用于设定症拟点, 一般用不着

注意

常见提示：WARNING: Residue 18 named GLU of a molecule in the input file was mapped to an entry in the topology database, but the atom CG used in i!,at entry Is not found in the input file. Perhaps your atom and/or residue naming needs to befixed.

  • 这说明 pdb2gmx在将输入的结构文件里的 第 18号残基 GLU 与 rtp文件里的 [GLU] 字段 匹配时,发现结构文件里这个残基少了名为 CG 的原子,显然需要补全结构文件里的这 个缺失的原子后再用 pdb2gmx 。但也有可 能这个原子实际存在,只是结构文件里的原 子名不对,此时需要修改原子名。
  • 注意这里的残基号是从 1 开始排的,和 pdb 文件里记录的残基号未必相同。

solvate

可以生成一个充满指定分子的盒子，也可以在-cp指定的结构文件中加溶剂分子把盒子空隙填满。
默认填入的溶剂是share/gromacs/top/spc216.gro，这是对216个SPC水组成的盒子经过动力学模拟后己跑到充分平衡状态的文件。如果要填入其它溶剂分子则用-cs[文件名]来指定，文件也应当是经过NPT动力学模拟后的单一纯溶剂体系结构文件。GROMACS会先把整个模拟盒子用此文件里的体系通过平移复制填湾，然后扣掉和溶质有重叠的溶剂分子。

  • -maxsol [数目]可以设定只填入多少溶剂分子，填充会很不均匀。
  • -shell [厚度]可以填充指定厚度的溶剂层。
  • -p [top文件名]还会同时对指定的拓扑文件进行修改，将其中的溶剂分子数设为当前值。
  • -scale <real>：当溶剂与溶质之间距离小于它们的范德华半径和，则溶剂会被去掉。范德华半径是share/gromacs/top/vdwradii.dat的半径值乘以0.57（默认值）。对水中的蛋白质, 使用默认值0.57可以得到接近1000 g/l的密度值。

例子

  • 生成2×5×4nm的水盒子体系

    1
    gmx solvate -box 2 5 4 -o test.gro

  • 在protein.gro的盒子空档中填充水分子

    1
    gmx solvate -cp protein.gro -o solv.gro

  • 在protein.gro中的溶质附近0.5nm范围内填满水分子(前提是.gro里的盒子己足够大)

    1
    gmx solvate -cp protein.gro -o solv.gro -shell 0.5

trjconv

此命令用于轨迹文件的处理，主要用处：

  • 转换轨格式(xtc/trr/gro/pdb之间转换)
  • 减少帧数
  • 只保留指定组的原子
  • 截取部分轨迹
  • 修改轨迹中对周期性的描述
  • 将轨迹向参考结构叠合

合并轨迹用trjcat

参数：

  • -f:被处理的轨迹文件
  • -s:.tpr/gro/pdb文件(可选，输出结构文件及-fit时必须)
  • -n:.ndx文件(可选)
  • -o:新产生的轨迹文件
  • -skip <int>:每几帧输出一次，默认为1
  • -dump <time>：输出离指定时间最近的帧
  • -pbc <选项>:修改PBC处理方式
    • 默认值为none，即不对轨迹的PBC描述做任何修
    • atom:只要原子越过盒子边界就返回另一边。这是mdrun直接产生出的轨迹对应的情况。这会导致边界处分子/残基被截断
    • mol/res:分子/残基的质心如果跨越了盒子，则分子/残基就挪到盒子另一头，这保证了分子残基的完整性。如果因为用了-trans和-center使得有的分子/残基处于盒子外，也可以由此得到修正
      • 注:mdrun直接产生的gro文件对应的是mol状态
    • nojump:使得原子即便跨越了盒子边界也照样运动(相对于第1帧而言)，而不会被挪到另一边。这保证了轨迹的连续性。
      • GROMACS自带的分析工具基本都能自动、恰当地在分析轨迹时考虑PBC问题。但第三方的分析程序则可能没法恰当处理PBC轨迹，如果不先搞成nojump的形式则会得到错误的结果。
    • cluster:选定一个组(通过index文件定义)，使得这个组保持完整。主要用于令轨迹中分子团簇(包括诸如蛋白质+配体复合物)保持完整。使用前应先使用mol保证分子的完整性
  • -ur <选项>:修改盒子描述方式
  • -center:平移体系使得指定组的几何中心处在盒子中央。盒子中心可以由-boxcenter具体定义
  • -trans <vector>:对所有坐标按照矢量平移
  • -fit <选项>:对轨迹中某个组根据参考结构叠合
  • -step:把轨迹每一帧输出为单独的.gro/pdb文件

读的文件里有速度信息就会被读/写。但受力不会被读/写，除非写了-force

例子

  • 将all.trr中3 ns到5 ns的轨迹每隔100ps一次提取出来到3_5.xtc

    1
    gmx trjcon -f all.trr -b 3000 -e 5000 -dt 100 -0 3_5.xtc

  • 从all.trr中提取出grp.ndx里定义的某个组到sub.trr

    1
    gmx trjconv -f all.trr -n grp.ndx -o sub.trr
    • 如果用-s提供了结构文件，即便没用-n也会提示选择要输出的组

  • 提取md.trr中最接近3000ps的结构和速度到3000ps.gro

    1
    gmx trjconv -f md.trr -s md.tpr -o 3000ps.gro -dump 3000

  • 将轨迹按照index.ndx里的某个组进行叠合输出到new.xtc，参考结构是ref.gro

    1
    gmx trjconv -f sol.trr -s ref.gro -n index.ndx -fit rot+trans -o new.xtc
    • 程序会让选择两次组，第一次选择对哪个组做叠合，第二次选择输出文件里包含哪个组的信息

  • 把md.xtc中从50ns开始到末尾的轨迹输出到cen.xtc，并对每一帧结构进行平移使得选定的某个组的质心恒处在盒子中心

    1
    gmx trjconv -f md.xtc -b 50000 -o cen.xtc -n index.ndx -center

  • 令sol.trr中的分子保持完整(必须提供.tpr)

    1
    gmx trjconv -f sol.trr -s sol.tpr -pbc mol -o new.xtc

    • 若再用以下命令，则分子运动轨迹会保持连续，而不被盒子边界所截断

      1
      gmx trjconv -f new.xtc -pbc nojump -o yuri.xtc
有用可戳(●ˇ∀ˇ●)