一条通往开源GROMACS分子动力学模拟仿真捷径（附视频教程）

;; Note the strange order of atoms to make it faster in gromacs.;[ moleculetype ]; molnamenrexclSOL2 [ atoms ]; idat typeres nr residu nameat namecg nrcharge1       O        1       SOL      O     1       0.002       H        1       SOL      H     1       0.58973       H        1       SOL      H     1       0.58974       M        1       SOL      M     1      -1.1794 [constraints];i j funct doh  dhh1       2       1       0.095721       3       1       0.095722       3       1       0.15139 [ exclusions ]1234213431244123 ; The position of the virtual site is computed as follows:;;O;        ;    D;  ;HH;; const = distance (OD) / [ cos (angle(DOH)) * distance (OH) ];  0.01577 nm/ [ cos (52.26 deg)* 0.09572 nm] ; Vsite pos x4 = x1 + a*(x2-x1) + b*(x3-x1) [ virtual_sites3 ]; Vsite fromfunctab412310.13458335      0.13458335

#include "oplsaa.ff/forcefield.itp"#include "oplsaa.ff/tip4p-ice.itp"

[ atoms ]; idat typeres nr residu nameat namecg nrcharge1       opls_113        1       SOL      O     1       0.002       opls_114        1       SOL      H     1       0.58973       opls_114        1       SOL      H     1       0.58974       opls_115        1       SOL      M     1      -1.1794

tolerance 2.0output icebox.pdbfiletype pdb structure ice.pdb number 1 inside box 0. 0. 0. 28. 24. 23.end structure structure ice.pdb number 1 inside box 0. 0. 46. 28. 23. 69.end structure structure water.pdb number 600 inside box 0. 0. 23. 28. 23. 46.end structure

CRYST1   28   23   69  90.00  90.00  90.00 P 1           1

gmx grompp -f em.mdp -c icebox.pdb -p icebox.top -o em.tpr -maxwarn 2gmx mdrun -v -deffnm em

cat em.gro |awk '{print $0;if($2=="OW"){a=$1;b=$3;c=$4;d=$5;e=$6;}if($2=="HW2")printf("%8s     MW %4d%8.3f%8.3f%8.3f\n",a,b,c,d,e)}' > newicebox.gro

gmx grompp -f em.mdp -c newicebox.gro -p icebox4.top -o em.tpr -maxwarn 2gmx mdrun -v -deffnm em

; Run parametersintegrator  = md        ; leap-frog integratornsteps      = 100000000  ; 1 * 100000000 = 10000 ps (10 ns)dt          = 0.001     ; 1 fs; Output controlnstxout     = 4000       ; save coordinates every 0.4 psnstvout     = 4000       ; save velocities every 0.4 psnstenergy   = 4000       ; save energies every 0.4 psnstlog      = 4000       ; update log file every 0.4 psnstxtcout   = 40000      ; xtc compressed trajectory output every 40 ps; Bond parameterscontinuation            = yes       ; Restarting after NVTconstraint_algorithm    = lincs     ; holonomic constraintsconstraints             = h-bonds   ; only h-bonds constrainedlincs_iter              = 1         ; accuracy of LINCSlincs_order             = 4         ; also related to accuracy; Neighborsearchingns_type     = grid      ; search neighboring grid celsnstlist     = 5         ; 5 fsrlist       = 1.0       ; short-range neighborlist cutoff (in nm)rcoulomb    = 1.0       ; short-range electrostatic cutoff (in nm)rvdw        = 1.0       ; short-range van der Waals cutoff (in nm); Electrostaticscoulombtype     = PME       ; Particle Mesh Ewald for long-range electrostaticspme_order       = 4         ; cubic interpolationfourierspacing  = 0.16      ; grid spacing for FFT; Temperature coupling is ontcoupl      = v-rescale                       ; More accurate thermostattc-grps     = ICE        SOL      ; two coupling groups - more accuratetau_t       = 0.5           0.5                 ; time constant, in psref_t       = 260           260                 ; reference temperature, one for each group, in K; Pressure coupling is onpcoupl      = Parrinello-Rahman     ; Pressure coupling on in NPTpcoupltype  = semiisotropic         ; uniform scaling of x-y box vectors, independent ztau_p       = 5.0                   ; time constant, in psref_p       = 450.0   450.0         ; reference pressure, x-y, z (in bar)compressibility = 4.5e-5    4.5e-5  ; isothermal compressibility, bar^-1; Periodic boundary conditionspbc         = xyz       ; 3-D PBC; Dispersion correctionDispCorr    = EnerPres  ; account for cut-off vdW scheme; Velocity generationgen_vel     = no        ; Velocity generation is off; COM motion removal; These options remove motion of the ice relative to the solvent/ionsnstcomm         = 1comm-mode       = Linearcomm-grps       = ICE SOL; Scale COM of reference coordinatesrefcoord_scaling = com