4. `mdpow.equil` — Setting up and running equilibrium MD¶

The mdpow.equil module facilitates the setup of equilibrium molecular dynamics simulations of a compound molecule in a simulation box of water or other solvent such as octanol.

It requires as input

the itp file for the compound

a coordinate (structure) file (in pdb or gro format)

By default it uses the OPLS/AA forcefield and the TIP4P water model.

Warning

Other forcefields than OPLS/AA are currently not officially supported; it is not hard to do but requires tedious changes to a few paths in template scripts.

class mdpow.equil.Simulation(molecule=None, forcefield: Union[mdpow.forcefields.Forcefield, str] = 'OPLS-AA', **kwargs)[source]¶

Simple MD simulation of a single compound molecule in water.

Typical use

S = Simulation(molecule='DRUG')
S.topology(itp='drug.itp')
S.solvate(struct='DRUG-H.pdb')
S.energy_minimize()
S.MD_relaxed()
S.MD()

Note

The OPLS/AA force field and the TIP4P water molecule is the default; changing this is possible but will require provision of customized itp, mdp and template top files at various stages.

Set up Simulation instance.

The molecule of the compound molecule should be supplied. Existing files (which have been generated in previous runs) can also be supplied.

Keywords:

molecule: Identifier for the compound molecule. This is the same as the entry in the [ molecule ] section of the itp file. [“DRUG”]
filename: If provided and molecule is None then load the instance from the pickle file filename, which was generated with save().
dirname: base directory; all other directories are created under it
forcefield: A Forcefield, or the string name of a packaged forcefield: ‘OPLS-AA’, ‘CHARMM’ or ‘AMBER’.
solvent: ‘water’ or ‘octanol’ or ‘cyclohexane’ or ‘wetoctanol’ or ‘toluene’
solventmodel: None chooses the default (e.g, the mdpow.forcefields.Forcefield default water model for solvent == "water". Other options are the models defined in mdpow.forcefields.GROMACS_WATER_MODELS. At the moment, there are no alternative parameterizations included for other solvents.
mdp: dict with keys corresponding to the stages energy_minimize, MD_restrained, MD_relaxed, MD_NPT and values mdp file names (if no entry then the package defaults are used)
distance: minimum distance between solute and closest box face
kwargs: advanced keywords for short-circuiting; see mdpow.equil.Simulation.filekeys.

Changed in version 0.9.0: forcefield may now be either a Forcefield or the string name of a builtin forcefield.

MD(**kwargs)¶

Short NPT MD simulation.

See documentation of gromacs.setup.MD() for details such as runtime or specific queuing system options. The following keywords can not be changed: top, mdp, ndx, mainselection.

Note

If the system crashes (with LINCS errors), try initial equilibration with timestep dt = 0.0001 ps (0.1 fs instead of 2 fs) and runtime = 5 ps as done in MD_relaxed()

Keywords:

struct: starting conformation; by default, the struct is the last frame from the position restraints run, or, if this file cannot be found (e.g. because Simulation.MD_restrained() was not run) it falls back to the relaxed and then the solvated system.
mdp: MDP run parameter file for Gromacs
runtime: total run time in ps
qscript: list of queuing system scripts to prepare; available values are in gromacs.config.templates or you can provide your own filename(s) in the current directory (see gromacs.qsub for the format of the templates)
qname: name of the job as shown in the queuing system
startdir: advanced uses: path of the directory on a remote system, which will be hard-coded into the queuing system script(s); see gromacs.setup.MD() and gromacs.manager.Manager

MD_NPT(**kwargs)[source]¶

Short NPT MD simulation.

See documentation of gromacs.setup.MD() for details such as runtime or specific queuing system options. The following keywords can not be changed: top, mdp, ndx, mainselection.

Note

If the system crashes (with LINCS errors), try initial equilibration with timestep dt = 0.0001 ps (0.1 fs instead of 2 fs) and runtime = 5 ps as done in MD_relaxed()

Keywords:

struct: starting conformation; by default, the struct is the last frame from the position restraints run, or, if this file cannot be found (e.g. because Simulation.MD_restrained() was not run) it falls back to the relaxed and then the solvated system.
mdp: MDP run parameter file for Gromacs
runtime: total run time in ps
qscript: list of queuing system scripts to prepare; available values are in gromacs.config.templates or you can provide your own filename(s) in the current directory (see gromacs.qsub for the format of the templates)
qname: name of the job as shown in the queuing system
startdir: advanced uses: path of the directory on a remote system, which will be hard-coded into the queuing system script(s); see gromacs.setup.MD() and gromacs.manager.Manager

MD_relaxed(**kwargs)[source]¶

Short MD simulation with timestep = 0.1 fs to relax strain.

Energy minimization does not always remove all problems and LINCS constraint errors occur. A very short runtime = 5 ps MD with very short integration time step dt tends to solve these problems.

Keywords:

struct: starting coordinates (typically guessed)
mdp: MDP run parameter file for Gromacs
qscript: list of queuing system submission scripts; probably a good idea to always include the default “local.sh” even if you have your own [“local.sh”]
qname: name of the job as shown in the queuing system
startdir: advanced uses: path of the directory on a remote system, which will be hard-coded into the queuing system script(s); see gromacs.setup.MD() and gromacs.manager.Manager

MD_restrained(**kwargs)[source]¶

Short MD simulation with position restraints on compound.

See documentation of gromacs.setup.MD_restrained() for details. The following keywords can not be changed: top, mdp, ndx, mainselection

Note

Position restraints are activated with -DPOSRES directives for gromacs.grompp(). Hence this will only work if the compound itp file does indeed contain a [ posres ] section that is protected by a #ifdef POSRES clause.

Keywords:

struct: starting coordinates (leave empty for inspired guess of file name)
mdp: MDP run parameter file for Gromacs
qscript: list of queuing system submission scripts; probably a good idea to always include the default “local.sh” even if you have your own [“local.sh”]
qname: name of the job as shown in the queuing system
startdir: advanced uses: path of the directory on a remote system, which will be hard-coded into the queuing system script(s); see gromacs.setup.MD() and gromacs.manager.Manager

coordinate_structures = ('solvated', 'energy_minimized', 'MD_relaxed', 'MD_restrained', 'MD_NPT')¶: Coordinate files of the full system in increasing order of advancement of the protocol; the later the better. The values are keys into Simulation.files.

energy_minimize(**kwargs)[source]¶

Energy minimize the solvated structure on the local machine.

kwargs are passed to gromacs.setup.energ_minimize() but if solvate() step has been carried out previously all the defaults should just work.

filekeys = ('topology', 'processed_topology', 'structure', 'solvated', 'ndx', 'energy_minimized', 'MD_relaxed', 'MD_restrained', 'MD_NPT')¶: Keyword arguments to pre-set some file names; they are keys in Simulation.files.

get_last_checkpoint()[source]¶: Returns the checkpoint of the most advanced step in the protocol. Relies on md.gro being present from previous simulation, assumes that checkpoint is then present.

get_last_structure()[source]¶: Returns the coordinates of the most advanced step in the protocol.

load(filename=None)[source]¶: Re-instantiate class from pickled file.

make_paths_relative(prefix='.')[source]¶: Hack to be able to copy directories around: prune basedir from paths.

Warning

This is not guaranteed to work for all paths. In particular, check mdpow.equil.Simulation.dirs.includes and adjust manually if necessary.

mdp_defaults = {'MD_NPT': 'NPT_opls.mdp', 'MD_relaxed': 'NPT_opls.mdp', 'MD_restrained': 'NPT_opls.mdp', 'energy_minimize': 'em_opls.mdp'}¶: Default Gromacs MDP run parameter files for the different stages. (All are part of the package and are found with mdpow.config.get_template().)

processed_topology(**kwargs)[source]¶: Create a portable topology file from the topology and the solvated system.

protocols = ('MD_NPT', 'MD_NPT_run', 'MD_relaxed', 'MD_relaxed_run', 'MD_restrained', 'MD_restrained_run', 'energy_minimize', 'solvate', 'topology')¶: Check list of all methods that can be run as an independent protocol; see also Simulation.get_protocol() and restart.Journal

save(filename=None)[source]¶

Save instance to a pickle file.

The default filename is the name of the file that was last loaded from or saved to.

solvate(struct=None, **kwargs)[source]¶

Solvate structure struct in a box of solvent.

The solvent is determined with the solvent keyword to the constructor.

Keywords:

struct

pdb or gro coordinate file (if not supplied, the value is used that was supplied to the constructor of Simulation)

distance

minimum distance between solute and the closes box face; the default depends on the solvent but can be set explicitly here, too.

bt

any box type understood by gromacs.editconf() (-bt):

“triclinic” is a triclinic box,
“cubic” is a rectangular box with all sides equal;
“dodecahedron” represents a rhombic dodecahedron;
“octahedron” is a truncated octahedron.

The default is “dodecahedron”.

kwargs

All other arguments are passed on to gromacs.setup.solvate(), but set to sensible default values. top and water are always fixed.

topology(itp='drug.itp', prm=None, **kwargs)[source]¶

Generate a topology for compound molecule.

Keywords:	itp Gromacs itp file; will be copied to topology dir and included in topology prm Gromacs prm file; if given, will be copied to topology dir and included in topology dirname name of the topology directory [“top”] kwargs see source for top_template, topol

class mdpow.equil.WaterSimulation(molecule=None, forcefield: Union[mdpow.forcefields.Forcefield, str] = 'OPLS-AA', **kwargs)[source]¶

Equilibrium MD of a solute in a box of water.

Set up Simulation instance.

The molecule of the compound molecule should be supplied. Existing files (which have been generated in previous runs) can also be supplied.

Keywords:

molecule: Identifier for the compound molecule. This is the same as the entry in the [ molecule ] section of the itp file. [“DRUG”]
filename: If provided and molecule is None then load the instance from the pickle file filename, which was generated with save().
dirname: base directory; all other directories are created under it
forcefield: A Forcefield, or the string name of a packaged forcefield: ‘OPLS-AA’, ‘CHARMM’ or ‘AMBER’.
solvent: ‘water’ or ‘octanol’ or ‘cyclohexane’ or ‘wetoctanol’ or ‘toluene’
solventmodel: None chooses the default (e.g, the mdpow.forcefields.Forcefield default water model for solvent == "water". Other options are the models defined in mdpow.forcefields.GROMACS_WATER_MODELS. At the moment, there are no alternative parameterizations included for other solvents.
mdp: dict with keys corresponding to the stages energy_minimize, MD_restrained, MD_relaxed, MD_NPT and values mdp file names (if no entry then the package defaults are used)
distance: minimum distance between solute and closest box face
kwargs: advanced keywords for short-circuiting; see mdpow.equil.Simulation.filekeys.

Changed in version 0.9.0: forcefield may now be either a Forcefield or the string name of a builtin forcefield.

class mdpow.equil.OctanolSimulation(molecule=None, forcefield: Union[mdpow.forcefields.Forcefield, str] = 'OPLS-AA', **kwargs)[source]¶

Equilibrium MD of a solute in a box of octanol.

Set up Simulation instance.

The molecule of the compound molecule should be supplied. Existing files (which have been generated in previous runs) can also be supplied.

Keywords:

molecule: Identifier for the compound molecule. This is the same as the entry in the [ molecule ] section of the itp file. [“DRUG”]
filename: If provided and molecule is None then load the instance from the pickle file filename, which was generated with save().
dirname: base directory; all other directories are created under it
forcefield: A Forcefield, or the string name of a packaged forcefield: ‘OPLS-AA’, ‘CHARMM’ or ‘AMBER’.
solvent: ‘water’ or ‘octanol’ or ‘cyclohexane’ or ‘wetoctanol’ or ‘toluene’
solventmodel: None chooses the default (e.g, the mdpow.forcefields.Forcefield default water model for solvent == "water". Other options are the models defined in mdpow.forcefields.GROMACS_WATER_MODELS. At the moment, there are no alternative parameterizations included for other solvents.
mdp: dict with keys corresponding to the stages energy_minimize, MD_restrained, MD_relaxed, MD_NPT and values mdp file names (if no entry then the package defaults are used)
distance: minimum distance between solute and closest box face
kwargs: advanced keywords for short-circuiting; see mdpow.equil.Simulation.filekeys.

Changed in version 0.9.0: forcefield may now be either a Forcefield or the string name of a builtin forcefield.

mdpow.equil.DIST = {'cyclohexane': 1.5, 'octanol': 1.5, 'toluene': 1.5, 'water': 1.0, 'wetoctanol': 1.5}¶: minimum distance between solute and box surface (in nm)

4. mdpow.equil — Setting up and running equilibrium MD¶

4. `mdpow.equil` — Setting up and running equilibrium MD¶