Example 1: Construct and equilibrate a box of methanol ====================================================== 0. Prerequisites ~~~~~~~~~~~~~~~~ Install `gmxbatch`. Create a directory and put the following files in it: - From the source distribution of `gmxbatch` (at https://gitlab.com/awacha/gmxbatch ): - `meoh.gro `_ from the `test/meoh` directory - `meoh.itp `_ from the `test/moleculetypes/charmm` directory - The `CHARMM36m force field bundle `_ from Alexander MacKerell's site. Download and extract 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