- GROMACS流程（使用GMX-5.0）
gmx help (module) 或 gmx (module) -h
学了AMBER之后来理解GROMACS更舒服了，而且中文讲解更清晰，AMBER在核酸领域比GROMACS更强，蛋白质差不多，材料领域用LAMMPS
本来想学gaussian的，但是gaussian还没弄到win版，先学gromacs吧
pdb文件——gmx=tleap 产生gro和top两个文件——gro处理后用于溶剂化+坐标文件mdp写复合文件tpr
预平衡（以毒素肽为例）
不是必需的，看初始结构是否平衡。包括：- 能量最小化：使得MD时快速达到平衡，也可以补充需要的离子然后min
- 蛋白质位置限制性模拟：限制重原子维持蛋白质结构，放开蛋白质分子的位置限制进行模拟
NVT预平衡：先做NVT减小盒子内压力，防NPT时体系崩溃，再做NPT
综上：min-NVT-NPT处理pdb
- 上网下载：wget http://www.rcsb.org/pdb/files/1OMB.pdb
需要处理的方面包括但不限于:
· 氢原子
· C端氧原子
· 结晶水
· 缺失原子/残基/侧链
· 二硫键
· 带电残基对于1OMB：去氢加上末尾C端结束的氧原子及其类型OXT使用win上 deepview转换，并去掉SPDBV开头行，直接到END就好
Build->Add C-terminal Oxygen (OXT)
File->Save->Current Layer注意，优秀的ftp系统能够处理Win-Linux的文本转换，但最好
to_unix filename filenameRedHat系统还可以用：dos2unix
fws.pdb文件中已经不含氢原子, 删除了SEQRES, SHEET, SSBOND, CONECT信息, 并且添加了OXT原子.
gmx获取拓扑—-pdb2gmx module
gmx-5.x
gmx pdb2gmx -f fws.pdb -o fws.gro -p fws.top -i fws.itp -ff amber99sb-ildn -water tip3p -ignh
忽略氢，力场amber99sb-ILDN file：pdb，输出gromacs文件，输出top文件，水盒子
产生结构文件gro，拓扑文件top，位置限制文件itp
可以不用ff，会问你用哪个力场-ignh: 忽略氢原子,因本pdb文件来自NMR，含氢。统一使用GROMOS力场的氢原子命名规则, 以免名称不一致.
-ff: 指定使用的力场. 用Amber99SB-ILDN力场. 也可以使用其他力场, 根据体系选择.
-f: 指定预处理的蛋白质结构文件
-o: 指定新生成的GROMACS文件名(也可以是其它的文件类型, 如pdb)
-p: 指定新生成的拓扑文件名, 包含所有原子及原子间的相互作用参数
-water: 指定水模型. 使用TIP3P水模型. 也可以使用如SPC/E.gmx pdb2gmx [-f [<.gro/.g96/…>]] [-o [<.gro/.g96/…>]] [-p [<.top>]]
[-i [<.itp>]] [-n [<.ndx>]] [-q [<.gro/.g96/...>]] [-chainsep <enum>] [-merge <enum>] [-ff <string>] [-water <enum>] [-[no]inter] [-[no]ss] [-[no]ter] [-[no]lys] [-[no]arg] [-[no]asp] [-[no]glu] [-[no]gln] [-[no]his] [-angle <real>] [-dist <real>] [-[no]una] [-[no]ignh] [-[no]missing] [-[no]v] [-posrefc <real>] [-vsite <enum>] [-[no]heavyh] [-[no]deuterate] [-[no]chargegrp] [-[no]cmap] [-[no]renum] [-[no]rtpres] -o [<.gro/.g96/...>] (conf.gro) Structure file: gro g96 pdb brk ent esp -p [<.top>] (topol.top) Topology file -i [<.itp>] (posre.itp) Include file for topology建模拟溶剂盒子—-editconf module(或许可以删除)
gmx-5.x
gmx editconf -f fws.gro -o fws-PB.gro -bt dodecahedron -d 1.2
-f: 输入蛋白结构
-o: 输出带模拟周期性溶剂盒子信息的结构文件 perodic box PB
-bt: 默认使用长方盒子, 可以使用-bt选项改变盒子类型, 如八面体, 十二面体等，菱形十二面体盒子接近球形, 计算效率最高
-d: 蛋白与模拟盒子在X, Y, Z方向上的最小距离，不能小于0.9 nm
生成一个 xxx-PBC.gro文件,得到了置于周期性菱形盒子中的蛋白质分子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]-bt
(triclinic) Box type for -box and -d: triclinic, cubic, dodecahedron, octahedron-box
(0 0 0) Box vector lengths (a,b,c)-angles
(90 90 90) Angles between the box vectors (bc,ac,ab)-d
(0) Distance between the solute and the box 可用于转换GROMACS文件(*.gro)和pdb文件(*.pdb) ''' editconf -f file.gro -o file.pdb '''（添加离子？）gmx genion -s
真空能量最小化 (vac min)准备—-grompp module
gmx grompp -f fws-VM.mdp -c fws-PB.gro -p fws.top -o fws-VM.tpr -maxwarn 1
参数文件mdp，GROMACS预处理器grompp(即gromacs pre-processor的缩写)进行处理。注：shift/cutoff的单位为nm, UA=united atom, AA=all-atom, DispCorr只能和周期性边界条件一起使用.
vac-min-fws.mdp
; 传递给预处理器的一些定义
define = -DFLEXIBLE ; 使用柔性水模型而非刚性模型, 这样最陡下降法可进一步最小化能量; 模拟类型, 结束控制, 输出控制参数
integrator = steep ; 指定使用最陡下降法进行能量最小化. 若设为cg则使用共轭梯度法
emtol = 500.0 ; 若力的最大值小于此值则认为能量最小化收敛(单位kJ mol^-1^ nm^-1^)
emstep = 0.01 ; 初始步长(nm)
nsteps = 1000 ; 在能量最小化中, 指定最大迭代次数
nstenergy = 10 ; 能量写出频率
energygrps = System ; 要写出的能量组; 近邻列表, 相互作用计算参数
nstlist = 1 ; 更新近邻列表的频率. 1表示每步都更新
ns_type = grid ; 近邻列表确定方法(simple或grid)
coulombtype = PME ; 计算长程静电的方法. PME为粒子网格Ewald方法, 还可以使用cut-off
rlist = 1.0 ; 短程力近邻列表的截断值
rcoulomb = 1.0 ; 长程库仑力的截断值
vdwtype = cut-off ; 计算范德华作用的方法
rvdw = 1.0 ; 范德华距离截断值
constraints = none ; 设置模型中使用的约束
pbc = xyz ; 3维周期性边界条件这里有一个初学者的bug，似乎与top文件电荷未处理有关，论坛指出加一个-maxwarn 1就好了，产生mdout.mdp和运行输入文件fws-VM.tpr（动力学输入文件，包括压缩的初始结构，拓扑，模拟力场信息）
gmx grompp [-f [<.mdp>]] [-c [<.gro/.g96/…>]] [-r [<.gro/.g96/…>]]
[-rb [<.gro/.g96/...>]] [-n [<.ndx>]] [-p [<.top>]] [-t [<.trr/.cpt/...>]] [-e [<.edr>]] [-ref [<.trr/.cpt/...>]] [-po [<.mdp>]] [-pp [<.top>]] [-o [<.tpr>]] [-imd [<.gro>]] [-[no]v] [-time <real>] [-[no]rmvsbds] [-maxwarn <int>] [-[no]zero] [-[no]renum]Options to specify input files:
-f [<.mdp>] (grompp.mdp)grompp input file with MD parameters-c [<.gro/.g96/…>] (conf.gro)
Structure file: gro g96 pdb brk ent esp tpr-r [<.gro/.g96/…>] (restraint.gro) (Opt.)
Structure file: gro g96 pdb brk ent esp tpr-rb [<.gro/.g96/…>] (restraint.gro) (Opt.)
Structure file: gro g96 pdb brk ent esp tpr-n [<.ndx>] (index.ndx) (Opt.)
Index file-p [<.top>] (topol.top)
Topology file-t [<.trr/.cpt/…>] (traj.trr) (Opt.)
Full precision trajectory: trr cpt tng-e [<.edr>] (ener.edr) (Opt.)
Energy fileOptions to specify input/output files:
-ref [<.trr/.cpt/…>] (rotref.trr) (Opt.)Full precision trajectory: trr cpt tngOptions to specify output files:
-po [<.mdp>] (mdout.mdp)grompp input file with MD parameters-pp [<.top>] (processed.top) (Opt.)
Topology file-o [<.tpr>] (topol.tpr)
Portable xdr run input file-imd [<.gro>] (imdgroup.gro) (Opt.)
Coordinate file in Gromos-87 formatOther options:
-[no]v (no)Be loud and noisy-time
(-1) Take frame at or first after this time.-[no]rmvsbds (yes)
Remove constant bonded interactions with virtual sites-maxwarn
(0) Number of allowed warnings during input processing. Not for normal use and may generate unstable systems-[no]zero (no)
Set parameters for bonded interactions without defaults to zero instead of generating an error-[no]renum (yes)
Renumber atomtypes and minimize number of atomtypes跑个动力学能量最小化—-mdrun(值得研究参数)
gmx mdrun -v -deffnm fws-VM
- 用tpr文件做md能量最小化
- -deffnm指定默认的文件名称, -v显示模拟过程中的信息
- 注：我们得到了日志文件fws-VM.log, 全精度轨迹文件fws-VM.trr, 能量文件fws-VM.edr, 能量最小化后的结构文件fws-VM.gro.
- 说明: 如果最小化不收敛, 可以将nsteps的数值增大再重新提交一次. 要重新提交作业的话, 你需要重新运行grompp.
对能量最小化的讨论 http://mdbbs.org/thread-22150-1-1.html
gmx mdrun [-s [<.tpr>]] [-cpi [<.cpt>]] [-table [<.xvg>]] [-tablep [<.xvg>]] [-tableb [<.xvg> [...]]] [-rerun [<.xtc/.trr/...>]] [-ei [<.edi>]] [-multidir [<dir> [...]]] [-awh [<.xvg>]] [-membed [<.dat>]] [-mp [<.top>]] [-mn [<.ndx>]] [-o [<.trr/.cpt/...>]] [-x [<.xtc/.tng>]] [-cpo [<.cpt>]] [-c [<.gro/.g96/...>]] [-e [<.edr>]] [-g [<.log>]] [-dhdl [<.xvg>]] [-field [<.xvg>]] [-tpi [<.xvg>]] [-tpid [<.xvg>]] [-eo [<.xvg>]] [-devout [<.xvg>]] [-runav [<.xvg>]] [-px [<.xvg>]] [-pf [<.xvg>]] [-ro [<.xvg>]] [-ra [<.log>]] [-rs [<.log>]] [-rt [<.log>]] [-mtx [<.mtx>]] [-if [<.xvg>]] [-swap [<.xvg>]] [-deffnm <string>] [-xvg <enum>] [-dd <vector>] [-ddorder <enum>] [-npme <int>] [-nt <int>] [-ntmpi <int>] [-ntomp <int>] [-ntomp_pme <int>] [-pin <enum>] [-pinoffset <int>] [-pinstride <int>] [-gpu_id <string>] [-gputasks <string>] [-[no]ddcheck] [-rdd <real>] [-rcon <real>] [-dlb <enum>] [-dds <real>] [-gcom <int>] [-nb <enum>] [-nstlist <int>] [-[no]tunepme] [-pme <enum>] [-pmefft <enum>] [-bonded <enum>] [-[no]v] [-pforce <real>] [-[no]reprod] [-cpt <real>] [-[no]cpnum] [-[no]append] [-nsteps <int>] [-maxh <real>] [-replex <int>] [-nex <int>] [-reseed <int>]-deffnm 统一所有文件名称
-x 可以将所有文件压缩进xtc或者是tng文件
-mtx 导出hessian matrix-s [<.tpr>] (topol.tpr) Portable xdr run input file -cpi [<.cpt>] (state.cpt) (Opt.) Checkpoint file -table [<.xvg>] (table.xvg) (Opt.) xvgr/xmgr file -tablep [<.xvg>] (tablep.xvg) (Opt.) xvgr/xmgr file -tableb [<.xvg> [...]] (table.xvg) (Opt.) xvgr/xmgr file -rerun [<.xtc/.trr/...>] (rerun.xtc) (Opt.) Trajectory: xtc trr cpt gro g96 pdb tng -ei [<.edi>] (sam.edi) (Opt.) ED sampling input -multidir [<dir> [...]] (rundir) (Opt.) Run directory -awh [<.xvg>] (awhinit.xvg) (Opt.) xvgr/xmgr file -membed [<.dat>] (membed.dat) (Opt.) Generic data file -mp [<.top>] (membed.top) (Opt.) Topology file -mn [<.ndx>] (membed.ndx) (Opt.) Index file Options to specify output files: -o [<.trr/.cpt/...>] (traj.trr) Full precision trajectory: trr cpt tng -x [<.xtc/.tng>] (traj_comp.xtc) (Opt.) Compressed trajectory (tng format or portable xdr format) -cpo [<.cpt>] (state.cpt) (Opt.) Checkpoint file -c [<.gro/.g96/...>] (confout.gro) Structure file: gro g96 pdb brk ent esp -e [<.edr>] (ener.edr) Energy file -g [<.log>] (md.log) Log file -dhdl [<.xvg>] (dhdl.xvg) (Opt.) xvgr/xmgr file -field [<.xvg>] (field.xvg) (Opt.) xvgr/xmgr file -tpi [<.xvg>] (tpi.xvg) (Opt.) xvgr/xmgr file -tpid [<.xvg>] (tpidist.xvg) (Opt.) xvgr/xmgr file -eo [<.xvg>] (edsam.xvg) (Opt.) xvgr/xmgr file -devout [<.xvg>] (deviatie.xvg) (Opt.) xvgr/xmgr file -runav [<.xvg>] (runaver.xvg) (Opt.) xvgr/xmgr file -px [<.xvg>] (pullx.xvg) (Opt.) xvgr/xmgr file -pf [<.xvg>] (pullf.xvg) (Opt.) xvgr/xmgr file -ro [<.xvg>] (rotation.xvg) (Opt.) xvgr/xmgr file -ra [<.log>] (rotangles.log) (Opt.) Log file -rs [<.log>] (rotslabs.log) (Opt.) Log file -rt [<.log>] (rottorque.log) (Opt.) Log file -mtx [<.mtx>] (nm.mtx) (Opt.) Hessian matrix -if [<.xvg>] (imdforces.xvg) (Opt.) xvgr/xmgr file -swap [<.xvg>] (swapions.xvg) (Opt.) xvgr/xmgr file Other options: -deffnm <string> Set the default filename for all file options -xvg <enum> (xmgrace) xvg plot formatting: xmgrace, xmgr, none -dd <vector> (0 0 0) Domain decomposition grid, 0 is optimize -ddorder <enum> (interleave) DD rank order: interleave, pp_pme, cartesian -npme <int> (-1) Number of separate ranks to be used for PME, -1 is guess -nt <int> (0) Total number of threads to start (0 is guess) -ntmpi <int> (0) Number of thread-MPI ranks to start (0 is guess) -ntomp <int> (0) Number of OpenMP threads per MPI rank to start (0 is guess) -ntomp_pme <int> (0) Number of OpenMP threads per MPI rank to start (0 is -ntomp) -pin <enum> (auto) Whether mdrun should try to set thread affinities: auto, on, off -pinoffset <int> (0) The lowest logical core number to which mdrun should pin the first thread -pinstride <int> (0) Pinning distance in logical cores for threads, use 0 to minimize the number of threads per physical core -gpu_id <string> List of unique GPU device IDs available to use -gputasks <string> List of GPU device IDs, mapping each PP task on each node to a device -[no]ddcheck (yes) Check for all bonded interactions with DD -rdd <real> (0) The maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates -rcon <real> (0) Maximum distance for P-LINCS (nm), 0 is estimate -dlb <enum> (auto) Dynamic load balancing (with DD): auto, no, yes -dds <real> (0.8) Fraction in (0,1) by whose reciprocal the initial DD cell size will be increased in order to provide a margin in which dynamic load balancing can act while preserving the minimum cell size. -gcom <int> (-1) Global communication frequency -nb <enum> (auto) Calculate non-bonded interactions on: auto, cpu, gpu -nstlist <int> (0) Set nstlist when using a Verlet buffer tolerance (0 is guess) -[no]tunepme (yes) Optimize PME load between PP/PME ranks or GPU/CPU (only with the Verlet cut-off scheme) -pme <enum> (auto) Perform PME calculations on: auto, cpu, gpu -pmefft <enum> (auto) Perform PME FFT calculations on: auto, cpu, gpu -bonded <enum> (auto) Perform bonded calculations on: auto, cpu, gpu -[no]v (no) Be loud and noisy -pforce <real> (-1) Print all forces larger than this (kJ/mol nm) -[no]reprod (no) Try to avoid optimizations that affect binary reproducibility -cpt <real> (15) Checkpoint interval (minutes) -[no]cpnum (no) Keep and number checkpoint files -[no]append (yes) Append to previous output files when continuing from checkpoint instead of adding the simulation part number to all file names -nsteps <int> (-2) Run this number of steps (-1 means infinite, -2 means use mdp option, smaller is invalid) -maxh <real> (-1) Terminate after 0.99 times this time (hours) -replex <int> (0) Attempt replica exchange periodically with this period (steps) -nex <int> (0) Number of random exchanges to carry out each exchange interval (N^3 is one suggestion). -nex zero or not specified gives neighbor replica exchange. -reseed <int> (-1) Seed for replica exchange, -1 is generate a seed溶剂化再最小化
gmx-5.x
gmx solvate -cp fws-VM.gro -cs spc216.gro -p fws.top -o fws-ion.gro-cp指定需要填充水分子的体系, 带模拟蛋白盒子, -cs指定使用SPC水模型进行填充, spc216是GROMACS统一的三位点水分子结构, -p修改体系的拓扑文件, 加入相应水分子的物理参数 -o指定填充水分子后的输出文件.
最小化参数文件
fws-VM-sol.mdp define = -DFLEXIBLE integrator = steep emtol = 250.0 nsteps = 5000 nstenergy = 10 energygrps = System nstlist = 1 ns_type = grid coulombtype = PME rlist = 1.0 rcoulomb = 1.0 rvdw = 1.0 constraints = none pbc = xyz准备命令同理
gmx grompp -f fws-VM-sol.mdp -c fws-ion.gro -p fws.top -o ion.tpr -maxwarn 1基于tpr文件加离子—-genion
gmx genion -s ion.tpr -o fws-b4em.gro -neutral -conc 0.15 -p fws.top
- -neutral选项保证体系总的净电荷为零, 体系呈电中性, -conc选项设定需要的离子浓度(这里为0.15 M). -g选项指定输出日志文件的名称. -norandom选项会依据静电势来放置离子而不是随机放置(默认), 但我们在这里使用随机放置方法.
- 运行这个命令时, 会提示选择一个连续的溶剂分子组, 选择13 (SOL), 回车, 程序会告知你有两个溶剂分子被氯离子代替. fws.top中会包含NA和CL离子. 注意, 这里出现的分子顺序必须与坐标文件中的一致.
说明: genion默认使用的盐为NaCl. 如果你需要使用不同的阳离子和阴离子, 可使用-pname(阳离子)和-nname(阴离子)分别指定阴阳离子的名称(根据相应力场的ions.itp文件中离子的设定), 还可以使用-pn和-nn分别指定添加的阴离子数目.
gmx grompp -f fws-VM-sol.mdp -c fws-b4em.gro -p fws.top -o fws-VM-sol.tpr
gmx mdrun -v -deffnm fws-VM-sol
位置限制
NVT
gmx grompp -f nvt-pr-md.mdp -c fws-VM-sol.gro -p addoxt.top -o nvt-pr.tpr -r em-sol.gro -nt 8
报错：
Cannot find position restraint file restraint.gro (option -r).From GROMACS-2018, you need to specify the position restraint coordinate files explicitly to avoid mistakes, although you can still use the same file as you specify for the -c option.版本迭代问题，在最新的教程下看到确实是加入了-r em-sol.gro
报错：
Fatal error:There is no domain decomposition for 32 ranks that is compatible with the given box and a minimum cell size of 0.95325 nm Change the number of ranks or mdrun option -rcon or -dds or your LINCS settings与core的thread问题有关 命令后加个 -nt 1
https://jerkwin.github.io/GMX/GMXtut-10/#5-%E9%94%99%E8%AF%AF%E4%BF%A1%E6%81%AF 对问题的探讨
或者不想那么多，改一个 -nt 2试试，-nt 1 真的太慢了感觉。
改成 -nt 8 线程明显快多了，不知道最高可以多少线程
On 1 MPI rank, each using 8 OpenMP threads- 并行运算
需要把盒子划分成小单元并行计算，使用PME方法计算计算静电相互作用时, 还可以指定一些节点单独用于PME计算, 这些节点被称为PME节点. GROMACS通常会自动进行这些设置, 但当模拟体系具有很高或很低的空间不均匀性时(例如, 模拟真空中的板块), 可能需要自己对这些设置进行调整以达到更好的性能.
- 并行运算
位置限制性预平衡模拟，使溶剂和离子弛豫运动，然后等压升温。走一个NVT和一个NPT
说明: 当模拟运行时间较长时, 你可以使用后台运行方法, 如
gmx-4.x: nohup mdrun -deffnm npt-pr &
gmx-5.x: nohup gmx mdrun -deffnm npt-pr &
为检查模拟运行进度, 可使用tail命令查看日志文件, 如
tail -n 25 npt-pr.log
或者
tail -f npt-pr.log持续输出
这个教程挺好的
http://www.mdtutorials.com/gmx/lysozyme/index.htmlhttp://www.mdtutorials.com/gmx/lysozyme/index.html
Screen
screen -S 文件名
比nohup好用 screen ctrl+a +d退出 ctrl+r继续 screen -ls有哪些窗口 screen -r 数字恢复
用top查看正在进行的进程
screen命令详解，https://www.cnblogs.com/mchina/archive/2013/01/30/2880680.html
调mdp文件integrator为nm，而不是steep或者md
Hessian十进制化
近来用GROMACS做简正模式分析（NMA）,做完之后，想做后续工作需要用到Hessian矩阵和特征值和特征向量。但是GMX输出的特征值是十进制格式，输出的Hessian(.mtx)和特征向量(.trr)都是二进制文件，无法直接读出。
结果花了好多天的力气，曾试图改GMX底层代码，看了好多天他的c文件，结果头大。无奈又翻了几遍手册，发现居然有专门读取GMX二进制文件的命令~ 具体操作如下
gromacs 4.5版本 gmxdump -f xxx.trr -om xxx.mdp >>xxx.txt
gromacs 5+版本 gmx dump -f xxx.trr -om xxx.mdp >>xxx.txt
这里如果不加上>>后边的这段，仅仅是写上xxx.mdp，是不会输出任何文件。我也不清楚为啥，只有将它转化为txt才有输出
双精度文本不支持gpu加速
Gromacs有Python包gromacs wrapper，Fl studio也有python包
数字特征分析：比如主要分布在哪里，（或者可以切片然后整体分析，值得考虑）
某个区和某个区配对排序，或者不同正则运动方式配对排序
排序后列矩阵，互为音高音强
每个音乐的特点
数字特征
音乐
所以能不能用python脚本处理输入的蛋白/分子结构然后输出个中国风音乐或者是电音？
音阶了解一下
进一步的问题是如何写这个音乐了，如何模块化实现呢？
不过可以先不在python框架下，而是先把流程实现了。
顺便可以整点web的东西？
（弯曲振动为音高，伸缩振动为音强）
eigenfreq.xvg
eigenval.xvg
eigenvec.trr
