|
| 1 | +# Running MD with GROMACS |
| 2 | +## DP/MM Simulation |
| 3 | +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`. |
| 4 | +### Topology Preparation |
| 5 | +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. |
| 6 | + |
| 7 | +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`): |
| 8 | +``` |
| 9 | +[ atomtypes ] |
| 10 | +;name btype mass charge ptype sigma epsilon |
| 11 | + c3 c3 0.0 0.0 A 0.339771 0.451035 |
| 12 | + hc hc 0.0 0.0 A 0.260018 0.087027 |
| 13 | +
|
| 14 | +[ moleculetype ] |
| 15 | +;name nrexcl |
| 16 | + methane 3 |
| 17 | +
|
| 18 | +[ atoms ] |
| 19 | +; nr type resnr residue atom cgnr charge mass |
| 20 | + 1 c3 1 MOL C1 1 -0.1068 12.010 |
| 21 | + 2 hc 1 MOL H1 2 0.0267 1.008 |
| 22 | + 3 hc 1 MOL H2 3 0.0267 1.008 |
| 23 | + 4 hc 1 MOL H3 4 0.0267 1.008 |
| 24 | + 5 hc 1 MOL H4 5 0.0267 1.008 |
| 25 | +
|
| 26 | +[ bonds ] |
| 27 | +; i j func b0 kb |
| 28 | + 1 2 5 |
| 29 | + 1 3 5 |
| 30 | + 1 4 5 |
| 31 | + 1 5 5 |
| 32 | +
|
| 33 | +[ exclusions ] |
| 34 | +; ai aj1 aj2 aj3 aj4 |
| 35 | + 1 2 3 4 5 |
| 36 | + 2 1 3 4 5 |
| 37 | + 3 1 2 4 5 |
| 38 | + 4 1 2 3 5 |
| 39 | + 5 1 2 3 4 |
| 40 | +``` |
| 41 | +For comparsion, the original topology file genearted by `acpype` will be: |
| 42 | +``` |
| 43 | +; methane_GMX.itp created by acpype (v: 2021-02-05T22:15:50CET) on Wed Sep 8 01:21:53 2021 |
| 44 | +
|
| 45 | +[ atomtypes ] |
| 46 | +;name bond_type mass charge ptype sigma epsilon Amb |
| 47 | + c3 c3 0.00000 0.00000 A 3.39771e-01 4.51035e-01 ; 1.91 0.1078 |
| 48 | + hc hc 0.00000 0.00000 A 2.60018e-01 8.70272e-02 ; 1.46 0.0208 |
| 49 | +
|
| 50 | +[ moleculetype ] |
| 51 | +;name nrexcl |
| 52 | + methane 3 |
| 53 | +
|
| 54 | +[ atoms ] |
| 55 | +; nr type resi res atom cgnr charge mass ; qtot bond_type |
| 56 | + 1 c3 1 MOL C1 1 -0.106800 12.01000 ; qtot -0.107 |
| 57 | + 2 hc 1 MOL H1 2 0.026700 1.00800 ; qtot -0.080 |
| 58 | + 3 hc 1 MOL H2 3 0.026700 1.00800 ; qtot -0.053 |
| 59 | + 4 hc 1 MOL H3 4 0.026700 1.00800 ; qtot -0.027 |
| 60 | + 5 hc 1 MOL H4 5 0.026700 1.00800 ; qtot 0.000 |
| 61 | +
|
| 62 | +[ bonds ] |
| 63 | +; ai aj funct r k |
| 64 | + 1 2 1 1.0970e-01 3.1455e+05 ; C1 - H1 |
| 65 | + 1 3 1 1.0970e-01 3.1455e+05 ; C1 - H2 |
| 66 | + 1 4 1 1.0970e-01 3.1455e+05 ; C1 - H3 |
| 67 | + 1 5 1 1.0970e-01 3.1455e+05 ; C1 - H4 |
| 68 | +
|
| 69 | +[ angles ] |
| 70 | +; ai aj ak funct theta cth |
| 71 | + 2 1 3 1 1.0758e+02 3.2635e+02 ; H1 - C1 - H2 |
| 72 | + 2 1 4 1 1.0758e+02 3.2635e+02 ; H1 - C1 - H3 |
| 73 | + 2 1 5 1 1.0758e+02 3.2635e+02 ; H1 - C1 - H4 |
| 74 | + 3 1 4 1 1.0758e+02 3.2635e+02 ; H2 - C1 - H3 |
| 75 | + 3 1 5 1 1.0758e+02 3.2635e+02 ; H2 - C1 - H4 |
| 76 | + 4 1 5 1 1.0758e+02 3.2635e+02 ; H3 - C1 - H4 |
| 77 | +``` |
| 78 | +### DeepMD Settings |
| 79 | +Before running simulation, we need to tell GROMACS to use DeepPotential by setting environment variable `GMX_DEEPMD_INPUT_JSON`: |
| 80 | +```bash |
| 81 | +export GMX_DEEPMD_INPUT_JSON=input.json |
| 82 | +``` |
| 83 | +Then, in your working directories, we have to write `input.json` file: |
| 84 | +```json |
| 85 | +{ |
| 86 | + "graph_file": "/path/to/graph.pb", |
| 87 | + "type_file": "type.raw", |
| 88 | + "index_file": "index.raw", |
| 89 | + "lambda": 1.0, |
| 90 | + "pbc": false |
| 91 | +} |
| 92 | +``` |
| 93 | +Here is an explanation for these settings: |
| 94 | ++ `graph_file` : The graph file (with suffix .pb) generated by `dp freeze` command |
| 95 | ++ `type_file` : File to specify DP atom types (in space-sepreated format). Here, `type.raw` looks like |
| 96 | +``` |
| 97 | +1 0 0 0 0 |
| 98 | +``` |
| 99 | ++ `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 |
| 100 | +``` |
| 101 | +0 1 2 3 4 |
| 102 | +``` |
| 103 | ++ `lambda`: Optional, default 1.0. Used in alchemical calculations. |
| 104 | ++ `pbc`: Optional, default true. If true, the GROMACS peroidic condition is passed to DeepMD. |
| 105 | + |
| 106 | +### Run Simulation |
| 107 | +Finally, you can run GROMACS using `gmx mdrun` as usual. |
| 108 | + |
| 109 | +## All-atom DP Simulation |
| 110 | +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`: |
| 111 | +``` |
| 112 | +[ atomtypes ] |
| 113 | +; name at.num mass charge ptype sigma epsilon |
| 114 | +HW 1 1.008 0.0000 A 0.00000e+00 0.00000e+00 |
| 115 | +OW 8 16.00 0.0000 A 0.00000e+00 0.00000e+00 |
| 116 | +``` |
| 117 | +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. |
| 118 | + |
| 119 | + |
| 120 | +However, we still recommend you run all-atom DP simulation using LAMMPS since it is more stable and efficient. |
0 commit comments