Skip to content

Latest commit

 

History

History
120 lines (107 loc) · 5.67 KB

gromacs.md

File metadata and controls

120 lines (107 loc) · 5.67 KB
Raw

Running MD with GROMACS

DP/MM Simulation

This part gives a simple tutorial on how to run a DP/MM simulation for methane in water, which means using DP for methane and TIP3P for water. All relevant files can be found in examples/methane.

Topology Preparation

Similar to QM/MM simulation, the internal interactions (including bond, angle, dihedrals, LJ, Columb) of the region descibed by a neural network potential (NNP) have to be turned off. In GROMACS, bonded interactions can be turned off by modifying [ bonds ], [ angles ], [ dihedrals ] and [ pairs ] sections. And LJ and Columb interactions must be turned off by [ exclusions ] section.

For example, if one wants to simulate ethane in water, using DeepPotential for methane and TIP3P for water, the topology of methane should be like the following (as presented in examples/methane/methane.itp):

[ atomtypes ]
;name btype  mass  charge ptype    sigma  epsilon
  c3    c3   0.0     0.0     A 0.339771 0.451035
  hc    hc   0.0     0.0     A 0.260018 0.087027

[ moleculetype ]
;name            nrexcl
 methane          3

[ atoms ]
; nr type  resnr residue atom  cgnr  charge   mass
  1   c3      1     MOL   C1     1 -0.1068 12.010
  2   hc      1     MOL   H1     2  0.0267  1.008
  3   hc      1     MOL   H2     3  0.0267  1.008
  4   hc      1     MOL   H3     4  0.0267  1.008
  5   hc      1     MOL   H4     5  0.0267  1.008

[ bonds ]
; i  j  func  b0  kb
 1  2     5        
 1  3     5        
 1  4     5        
 1  5     5        

[ exclusions ]
; ai  aj1  aj2  aj3  aj4
  1    2    3    4    5
  2    1    3    4    5
  3    1    2    4    5
  4    1    2    3    5
  5    1    2    3    4

For comparsion, the original topology file genearted by acpype will be:

; methane_GMX.itp created by acpype (v: 2021-02-05T22:15:50CET) on Wed Sep  8 01:21:53 2021

[ atomtypes ]
;name   bond_type     mass     charge   ptype   sigma         epsilon       Amb
 c3       c3          0.00000  0.00000   A     3.39771e-01   4.51035e-01 ; 1.91  0.1078
 hc       hc          0.00000  0.00000   A     2.60018e-01   8.70272e-02 ; 1.46  0.0208

[ moleculetype ]
;name            nrexcl
 methane          3

[ atoms ]
;   nr  type  resi  res  atom  cgnr     charge      mass       ; qtot   bond_type
     1   c3     1   MOL    C1    1    -0.106800     12.01000 ; qtot -0.107
     2   hc     1   MOL    H1    2     0.026700      1.00800 ; qtot -0.080
     3   hc     1   MOL    H2    3     0.026700      1.00800 ; qtot -0.053
     4   hc     1   MOL    H3    4     0.026700      1.00800 ; qtot -0.027
     5   hc     1   MOL    H4    5     0.026700      1.00800 ; qtot 0.000

[ bonds ]
;   ai     aj funct   r             k
     1      2   1    1.0970e-01    3.1455e+05 ;     C1 - H1    
     1      3   1    1.0970e-01    3.1455e+05 ;     C1 - H2    
     1      4   1    1.0970e-01    3.1455e+05 ;     C1 - H3    
     1      5   1    1.0970e-01    3.1455e+05 ;     C1 - H4    

[ angles ]
;   ai     aj     ak    funct   theta         cth
     2      1      3      1    1.0758e+02    3.2635e+02 ;     H1 - C1     - H2    
     2      1      4      1    1.0758e+02    3.2635e+02 ;     H1 - C1     - H3    
     2      1      5      1    1.0758e+02    3.2635e+02 ;     H1 - C1     - H4    
     3      1      4      1    1.0758e+02    3.2635e+02 ;     H2 - C1     - H3    
     3      1      5      1    1.0758e+02    3.2635e+02 ;     H2 - C1     - H4    
     4      1      5      1    1.0758e+02    3.2635e+02 ;     H3 - C1     - H4

DeepMD Settings

Before running simulation, we need to tell GROMACS to use DeepPotential by setting environment variable GMX_DEEPMD_INPUT_JSON:

export GMX_DEEPMD_INPUT_JSON=input.json

Then, in your working directories, we have to write input.json file:

{
    "graph_file": "/path/to/graph.pb",
    "type_file": "type.raw",
    "index_file": "index.raw",
    "lambda": 1.0,
    "pbc": false
}

Here is an explanation for these settings:

  • graph_file : The graph file (with suffix .pb) generated by dp freeze command
  • type_file : File to specify DP atom types (in space-sepreated format). Here, type.raw looks like
1 0 0 0 0
  • index_file : File containing indices of DP atoms (in space-seperated format), which should be in consistent with indices' order in .gro file but starting from zero. Here, index.raw looks like
0 1 2 3 4
  • lambda: Optional, default 1.0. Used in alchemical calculations.
  • pbc: Optional, default true. If true, the GROMACS peroidic condition is passed to DeepMD.

Run Simulation

Finally, you can run GROMACS using gmx mdrun as usual.

All-atom DP Simulation

This part gives an example on how to run a simulation with all atoms described by a DeepPotential with Gromacs, taking water as an example. Instead of using [ exclusions ] to turn off the non-bonded energies, we can simply do this by setting LJ parameters (i.e. epsilon and sigma) and partial charges to 0, as shown in examples/water/gmx/water.top:

[ atomtypes ]
; name      at.num  mass     charge ptype  sigma      epsilon
HW           1       1.008   0.0000  A   0.00000e+00  0.00000e+00
OW           8      16.00    0.0000  A   0.00000e+00  0.00000e+00

As mentioned in the above section, input.json and relevant files (index.raw, type.raw) should also be created. Then, we can start the simulation under NVT ensemble and plot the radial distribution function (RDF) by gmx rdf command. We can see that the RDF given by Gromacs+DP matches perfectly with Lammps+DP, which further provides an evidence on the validity of our simulation. rdf

However, we still recommend you run all-atom DP simulation using LAMMPS since it is more stable and efficient.