Example 1: Construct and equilibrate a box of methanol

0. Prerequisites

Install gmxbatch.

Create a directory and put the following files in it:

Change to that directory.

1. Initialization

In the first part of a simulation script with gmxbatch is the initialization. This is probably the most important part as all parameters are decided here.

Let us first import the required packages

>>> import gmxbatch
>>> import matplotlib.pyplot as plt  # for plotting

Of course we will need a GMX engine:

>>> engine = gmxbatch.MDEngine('gmx')  # use the default 'gmx' executable on the $PATH

Define the force field we want to use. We choose the CHARMM36m force field, because it has parameters for methanol through CGenFF. You should already have downloaded and extracted it to your working directory. and extract it in your working directory. You should have a directory called ‘charmm36-mar2019.ff’.

The force field object in gmxbatch comprises three aspects: the atom and interaction definitions (through the forcefield.itp topology include file), default MDP options, and a list of known molecule types. Molecule types for common solvents and ions are usually found in the force field directory in corresponding .itp files. When initializing a force field object, two required parameters are the ITP file of the force field itself and a search path for ITP files containing molecule types.

At present gmxbatch can handle Amber, CHARMM and GROMOS force fields, but it should be straightforward to extend the code with others. We use here the just downloaded CHARMM36m force field, and add the force field directory and the current directory as search paths for molecule type ITPs.

>>> ff = gmxbatch.CHARMM(itp='charmm36-mar2019.ff/forcefield.itp', moltypespath=['charmm36-mar2019.ff', '.'])

Now we define the environment, complete with a thermostat set at 300 K with 1 ps coupling time, and an isotropic Parrinello-Rahman barostat at 1 bar with 5 ps coupling time.

>>> env = gmxbatch.Environment(
...           thermostat=gmxbatch.Thermostat(
...               groups=['System'],  # only one coupling group: only solvent
...               ref_temperature=300,  # K
...               tau=1,  # ps
...               algorithm='V-rescale'  # Velocity rescaling algorithm for the _production_ run
...           ),
...           barostat=gmxbatch.Barostat(
...               ref_pressure=1,  # bar
...               tau=5,  # ps
...               couplingtype='isotropic',  # isotropic system
...               compressibility=4.5e-5,  # 1/bar; compressibility of water
...               algorithm='Parrinello-Rahman'  # used for the _production_ run only
...           ))

Now assemble the system. Take the initial state from meoh.gro in your working directory

>>> conf = gmxbatch.Coordinates('meoh.gro')

Assemble the system

>>> system = gmxbatch.System(
...              name='Methanol',
...              forcefield=ff,
...              conf=conf,
...              moleculetypes=[ff.moleculetype('MEOH', 1, 'meoh.itp')],
...              indexgroups=None)

Finally, create the Simulation instance:

>>> sim = gmxbatch.Simulation(
...           engine=engine,
...           system=system,
...           environment=env)

Sort the atoms in the coordinate set to match the topology:

>>> sim.sortAtoms()

Now the initialization part is done.

2. Simulation

The strength of gmxbatch is in its simplicity when coming to this part. Many simulation tasks involving more calls to gmx subprograms are represented as one function call, as you will see below.

a. Adjust box size

First we adjust the box size. We put our single methanol molecule in a cubic box with 0.5 nm minimal distance from the molecule in the center.

>>> sim.rebox('cubic', 0.5)

Under the hood, a new file is generated, and sim.conf is updated to show the re-boxed structure.

b. Energy minimization of a single methanol molecule

Executing MD runs is really cheap:

>>> results = sim.em(nsteps=10000)

This also updates the underlying system. Additionally, each MD run returns a Results object with the trajectory, the final state of the run and a number of energy terms. We will see them later.

c. Replicate the molecules

Make a small box and do the replication:

>>> sim.rebox('cubic', 0.1)
>>> print(sim.system.moleculetypes[0])
MoleculeType: MEOH
   kind: Solvent
   count: 1
   itp: ../moleculetypes/charmm/meoh.itp
   atoms: 6
   posres macro: POSRES_MEOH
>>> print(len(sim.system.conf))
6
>>> sim.repeat(10, 10, 10, rot=True)  # 10×10×10 molecules, each is randomly oriented
>>> print(sim.system.moleculetypes[0])
MoleculeType: MEOH
   kind: Solvent
   count: 1000
   itp: ../moleculetypes/charmm/meoh.itp
   atoms: 6
   posres macro: POSRES_MEOH
>>> print(len(sim.system.conf))
6000

As you see, both the topology and the coordinate set has been updated.

d. Energy minimization

Do a short energy minimization on the box.

>>> emresults = sim.em(nsteps=50000)
>>> f=plt.figure()  # create a new Matplotlib figure
>>> emresults.energy['potential'].plot(figure=fig)  # plot the Potential energy vs. "time"

e. Equilibration in the NVT ensemble

Perform equilibration in the NVT ensemble:

>>> nvtresults = sim.nvt(100, deffnm='nvt')  # 100 ps
>>> nvtresults.energy['potential'].plot(figure=plt.figure())  # plot the potential energy vs. the time
>>> nvtresults.energy['temperature'].plot(figure=plt.figure())  # see if the temperature has been stabilized or not

f. Compress the system

The methanol molecules are still very much apart. Apply a large pressure (e.g. 1000 bar) to compress them.

>>> env.barostat.ref_pressure = 1000  # bar
>>> npt1000results = sim.npt(200, deffnm='npt1000bar')  # 0.2 ns
>>> npt1000results.energy['potential'].plot(figure=plt.figure())  # plot the potential energy vs. the time
>>> npt1000results.energy['temperature'].plot(figure=plt.figure())  # plot the temperature vs. the time
>>> npt1000results.energy['volume'].plot(figure=plt.figure())  # plot the volume vs. the time

g. Equilibration in the NpT ensemble

Now we set the reference pressure back to 1 bar and re-equilibrate the system at normal atmospheric pressure.

>>> env.barostat.ref_pressure = 1
>>> nptresults = sim.npt(200, deffnm='npt')  # 0.2 ns
>>> nptresults.energy['potential'].plot(figure=plt.figure())  # plot the potential energy vs. the time
>>> nptresults.energy['temperature'].plot(figure=plt.figure())  # plot the temperature vs. the time
>>> nptresults.energy['volume'].plot(figure=plt.figure())  # plot the volume vs. the time

h. Production MD

Do a short production run on the equilibrated system

>>> mdresults = sim.md(1000, deffnm='md')  # 1 ns
>>> mdresults.energy['potential'].plot(figure=plt.figure())  # plot the potential energy vs. the time
>>> mdresults.energy['temperature'].plot(figure=plt.figure())  # plot the temperature vs. the time
>>> mdresults.energy['volume'].plot(figure=plt.figure())  # plot the volume vs. the time
>>> sim.conf.write('meoh1000final.gro')  # write the resulting state to an output file