gromacs.cbook – Gromacs Cook Book

The cbook (cook book) module contains short recipes for tasks that are often repeated. In the simplest case this is just one of the gromacs tools with a certain set of default command line options.

By abstracting and collecting these invocations here, errors can be reduced and the code snippets can also serve as canonical examples for how to do simple things.

Miscellaneous canned Gromacs commands

Simple commands with new default options so that they solve a specific problem (see also Manipulating trajectories and structures):

gromacs.cbook.rmsd_backbone([s="md.tpr", f="md.xtc"[, ...]])

Computes the RMSD of the “Backbone” atoms after fitting to the “Backbone” (including both translation and rotation).

Manipulating trajectories and structures

Standard invocations for manipulating trajectories.

gromacs.cbook.trj_compact([s="md.tpr", f="md.xtc", o="compact.xtc"[, ...]])

Writes an output trajectory or frame with a compact representation of the system centered on the protein. It centers on the group “Protein” and outputs the whole “System” group.

gromacs.cbook.trj_xyfitted([s="md.tpr", f="md.xtc"[, ...]])

Writes a trajectory centered and fitted to the protein in the XY-plane only.

This is useful for membrane proteins. The system must be oriented so that the membrane is in the XY plane. The protein backbone is used for the least square fit, centering is done for the whole protein., but this can be changed with the input = ('backbone', 'protein','system') keyword.

Note

Gromacs 4.x only

gromacs.cbook.trj_fitandcenter(xy=False, **kwargs)

Center everything and make a compact representation (pass 1) and fit the system to a reference (pass 2).

Keywords :

s

input structure file (tpr file required to make molecule whole); if a list or tuple is provided then s[0] is used for pass 1 (should be a tpr) and s[1] is used for the fitting step (can be a pdb of the whole system)

If a second structure is supplied then it is assumed that the fitted trajectory should not be centered.

f

input trajectory

o

output trajectory

input

A list with three groups. The default is

[‘backbone’, ‘protein’,’system’]

The fit command uses all three (1st for least square fit, 2nd for centering, 3rd for output), the centered/make-whole stage use 2nd for centering and 3rd for output.

input1

If input1 is supplied then input is used exclusively for the fitting stage (pass 2) and input1 for the centering (pass 1).

n

Index file used for pass 1 and pass 2.

n1

If n1 is supplied then index n1 is only used for pass 1 (centering) and n for pass 2 (fitting).

xy : boolean

If True then only do a rot+trans fit in the xy plane (good for membrane simulations); default is False.

kwargs

All other arguments are passed to Trjconv.

Note that here we first center the protein and create a compact box, using -pbc mol -ur compact -center -boxcenter tric and write an intermediate xtc. Then in a second pass we perform a rotation+translation fit (or restricted to the xy plane if xy = True is set) on the intermediate xtc to produce the final trajectory. Doing it in this order has the disadvantage that the solvent box is rotating around the protein but the opposite order (with center/compact second) produces strange artifacts where columns of solvent appear cut out from the box—it probably means that after rotation the information for the periodic boundaries is not correct any more.

Most kwargs are passed to both invocations of gromacs.tools.Trjconv so it does not really make sense to use eg skip; in this case do things manually.

By default the input to the fit command is (‘backbone’, ‘protein’,’system’); the compact command always uses the second and third group for its purposes or if this fails, prompts the user.

Both steps cannot performed in one pass; this is a known limitation of trjconv. An intermediate temporary XTC files is generated which should be automatically cleaned up unless bad things happened.

The function tries to honour the input/output formats. For instance, if you want trr output you need to supply a trr file as input and explicitly give the output file also a trr suffix.

Note

For big trajectories it can take a very long time and consume a large amount of temporary diskspace.

We follow the g_spatial documentation in preparing the trajectories:

trjconv -s a.tpr -f a.xtc -o b.xtc -center -boxcenter tric -ur compact -pbc mol
trjconv -s a.tpr -f b.xtc -o c.xtc -fit rot+trans

gromacs.cbook.cat(prefix='md', dirname='.', partsdir='parts', fulldir='full', resolve_multi='pass')

Concatenate all parts of a simulation.

The xtc, trr, and edr files in dirname such as prefix.xtc, prefix.part0002.xtc, prefix.part0003.xtc, ... are

1. moved to the partsdir (under dirname)
2. concatenated with the Gromacs tools to yield prefix.xtc, prefix.trr, prefix.edr, prefix.gro (or prefix.md) in dirname
3. Store these trajectories in fulldir

Note

Trajectory files are never deleted by this function to avoid data loss in case of bugs. You will have to clean up yourself by deleting dirname/partsdir.

Symlinks for the trajectories are not handled well and break the function. Use hard links instead.

Warning

If an exception occurs when running this function then make doubly and triply sure where your files are before running this function again; otherwise you might overwrite data. Possibly you will need to manually move the files from partsdir back into the working directory dirname; this should onlu overwrite generated files so far but check carefully!

Keywords :

prefix

deffnm of the trajectories [md]

*resolve_multi”

how to deal with multiple “final” gro or pdb files: normally there should only be one but in case of restarting from the checkpoint of a finished simulation one can end up with multiple identical ones.

• “pass” : do nothing and log a warning

• “guess” : take prefix.pdb or prefix.gro if it exists, otherwise the one of

  prefix.partNNNN.gro|pdb with the highes NNNN

dirname

change to dirname and assume all tarjectories are located there [.]

partsdir

directory where to store the input files (they are moved out of the way); partsdir must be manually deleted [parts]

fulldir

directory where to store the final results [full]

class gromacs.cbook.Frames(structure, trj, maxframes=None, format='pdb', **kwargs)

A iterator that transparently provides frames from a trajectory.

The iterator chops a trajectory into individual frames for analysis tools that only work on separate structures such as gro or pdb files. Instead of turning the whole trajectory immediately into pdb files (and potentially filling the disk), the iterator can be instructed to only provide a fixed number of frames and compute more frames when needed.

Note

Setting a limit on the number of frames on disk can lead to longish waiting times because trjconv must re-seek to the middle of the trajectory and the only way it can do this at the moment is by reading frames sequentially. This might still be preferrable to filling up a disk, though.

Warning

The maxframes option is not implemented yet; use the dt option or similar to keep the number of frames manageable.

Set up the Frames iterator.

Arguments :

structure

name of a structure file (tpr, pdb, ...)

trj

name of the trajectory (xtc, trr, ...)

format

output format for the frames, eg “pdb” or “gro” [pdb]

maxframes : int

maximum number of frames that are extracted to disk at one time; set to None to extract the whole trajectory at once. [None]

kwargs

All other arguments are passed to class:~gromacs.tools.Trjconv; the only options that cannot be changed are sep and the output file name o.

all_frames

Unordered list of all frames currently held on disk.

cleanup()

Clean up all temporary frames (which can be HUGE).

delete_frames()

Delete all frames.

extract()

Extract frames from the trajectory to the temporary directory.

framenumber = None

Holds the current frame number of the currently extracted batch of frames. Increases when iterating.

totalframes = None

Total number of frames read so far; only important when maxframes > 0 is used.

class gromacs.cbook.Transformer(s='topol.tpr', f='traj.xtc', n=None, force=None, dirname='.', outdir=None)

Class to handle transformations of trajectories.

1. Center, compact, and fit to reference structure in tpr (optionally, only center in the xy plane): center_fit()
2. Write compact xtc and tpr with water removed: strip_water()
3. Write compact xtc and tpr only with protein: keep_protein_only()

Set up Transformer with structure and trajectory.

Supply n = tpr, f = xtc (and n = ndx) relative to dirname.

Keywords :

s

tpr file (or similar); note that this should not contain position restraints if it is to be used with a reduced system (see strip_water())

f

trajectory (xtc, trr, ...)

n

index file (it is typically safe to leave this as None; in cases where a trajectory needs to be centered on non-standard groups this should contain those groups)

force

Set the default behaviour for handling existing files:

• True: overwrite existing trajectories
• False: throw a IOError exception
• None: skip existing and log a warning [default]

dirname

directory in which all operations are performed, relative paths are interpreted relative to dirname [.]

outdir

directory under which output files are placed; by default the same directory where the input files live

center_fit(**kwargs)

Write compact xtc that is fitted to the tpr reference structure.

See gromacs.cbook.trj_fitandcenter() for details and description of kwargs. The most important ones are listed here but in most cases the defaults should work.

Keywords :

s

Input structure (typically the default tpr file but can be set to some other file with a different conformation for fitting)

n

Alternative index file.

o

Name of the output trajectory.

xy : Boolean

If True then only fit in xy-plane (useful for a membrane normal to z). The default is False.

force

• True: overwrite existing trajectories
• False: throw a IOError exception
• None: skip existing and log a warning [default]

Returns :

dictionary with keys tpr, xtc, which are the names of the the new files

fit(xy=False, **kwargs)

Write xtc that is fitted to the tpr reference structure.

Runs gromacs.tools.trjconv with appropriate arguments for fitting. The most important kwargs are listed here but in most cases the defaults should work.

Note that the default settings do not include centering or periodic boundary treatment as this often does not work well with fitting. It is better to do this as a separate step (see center_fit() or gromacs.cbook.trj_fitandcenter())

Keywords :

s

Input structure (typically the default tpr file but can be set to some other file with a different conformation for fitting)

n

Alternative index file.

o

Name of the output trajectory. A default name is created. If e.g. dt = 100 is one of the kwargs then the default name includes “_dt100ps”.

xy : boolean

If True then only do a rot+trans fit in the xy plane (good for membrane simulations); default is False.

force

True: overwrite existing trajectories False: throw a IOError exception None: skip existing and log a warning [default]

fitgroup

index group to fit on [“backbone”]

kwargs

kwargs are passed to trj_xyfitted()

Returns :

dictionary with keys tpr, xtc, which are the names of the the new files

keep_protein_only(os=None, o=None, on=None, compact=False, groupname='proteinonly', **kwargs)

Write xtc and tpr only containing the protein.

Keywords :

os

Name of the output tpr file; by default use the original but insert “proteinonly” before suffix.

o

Name of the output trajectory; by default use the original name but insert “proteinonly” before suffix.

on

Name of a new index file.

compact

True: write a compact and centered trajectory False: use trajectory as it is [False]

groupname

Name of the protein-only group.

keepalso

List of literal make_ndx selections of additional groups that should be kept, e.g. [‘resname DRUG’, ‘atom 6789’].

force : Boolean

• True: overwrite existing trajectories
• False: throw a IOError exception
• None: skip existing and log a warning [default]

kwargs

are passed on to gromacs.cbook.trj_compact() (unless the values have to be set to certain values such as s, f, n, o keywords). The input keyword is always mangled: Only the first entry (the group to centre the trajectory on) is kept, and as a second group (the output group) groupname is used.

Returns :

dictionary with keys tpr, xtc, ndx which are the names of the the new files

Warning

The input tpr file should not have any position restraints; otherwise Gromacs will throw a hissy-fit and say

Software inconsistency error: Position restraint coordinates are missing

(This appears to be a bug in Gromacs 4.x.)

outfile(p)

Path for an output file.

If outdir is set then the path is outdir/basename(p) else just p

rp(*args)

Return canonical path to file under dirname with components args

If args form an absolute path then just return it as the absolute path.

strip_fit(**kwargs)

Strip water and fit to the remaining system.

First runs strip_water() and then fit(); see there for arguments.

• strip_input is used for strip_water() (but is only useful in special cases, e.g. when there is no Protein group defined. Then set strip_input = ['Other'].
• input is passed on to fit() and can contain the [center_group, fit_group, output_group]
• By default fit = “rot+trans” (and fit is passed to fit(), together with the xy = False keyword)

strip_water(os=None, o=None, on=None, compact=False, resn='SOL', groupname='notwater', **kwargs)

Write xtc and tpr with water (by resname) removed.

Keywords :

os

Name of the output tpr file; by default use the original but insert “nowater” before suffix.

o

Name of the output trajectory; by default use the original name but insert “nowater” before suffix.

on

Name of a new index file (without water).

compact

True: write a compact and centered trajectory False: use trajectory as it is [False]

resn

Residue name of the water molecules; all these residues are excluded.

groupname

Name of the group that is generated by subtracting all waters from the system.

force : Boolean

• True: overwrite existing trajectories
• False: throw a IOError exception
• None: skip existing and log a warning [default]

kwargs

are passed on to gromacs.cbook.trj_compact() (unless the values have to be set to certain values such as s, f, n, o keywords). The input keyword is always mangled: Only the first entry (the group to centre the trajectory on) is kept, and as a second group (the output group) groupname is used.

Returns :

dictionary with keys tpr, xtc, ndx which are the names of the the new files

Warning

The input tpr file should not have any position restraints; otherwise Gromacs will throw a hissy-fit and say

Software inconsistency error: Position restraint coordinates are missing

(This appears to be a bug in Gromacs 4.x.)

gromacs.cbook.get_volume(f)

Return the volume in nm^3 of structure file f.

(Uses gromacs.editconf(); error handling is not good)

Processing output

There are cases when a script has to to do different things depending on the output from a Gromacs tool.

For instance, a common case is to check the total charge after grompping a tpr file. The grompp_qtot function does just that.

gromacs.cbook.grompp_qtot(*args, **kwargs)

Run gromacs.grompp and return the total charge of the system.

Arguments :The arguments are the ones one would pass to gromacs.grompp(). Returns :The total charge as reported

Some things to keep in mind:

• The stdout output of grompp is only shown when an error occurs. For debugging, look at the log file or screen output and try running the normal gromacs.grompp() command and analyze the output if the debugging messages are not sufficient.

• Check that qtot is correct. Because the function is based on pattern matching of the informative output of grompp it can break when the output format changes. This version recognizes lines like

  '  System has non-zero total charge: -4.000001e+00'
  

  using the regular expression System has non-zero total charge: *(?P<qtot>[-+]?d*.d+([eE][-+]d+)?).

gromacs.cbook.get_volume(f)

Return the volume in nm^3 of structure file f.

(Uses gromacs.editconf(); error handling is not good)

gromacs.cbook.parse_ndxlist(output)

Parse output from make_ndx to build list of index groups:

groups = parse_ndxlist(output)

output should be the standard output from make_ndx, e.g.:

rc,output,junk = gromacs.make_ndx(..., input=('', 'q'), stdout=False, stderr=True)

(or simply use

rc,output,junk = cbook.make_ndx_captured(...)

which presets input, stdout and stderr; of course input can be overriden.)

Returns :

The function returns a list of dicts (groups) with fields

name

name of the groups

nr

number of the group (starts at 0)

natoms

number of atoms in the group

Working with topologies and mdp files

gromacs.cbook.create_portable_topology(topol, struct, **kwargs)

Create a processed topology.

The processed (or portable) topology file does not contain any #include statements and hence can be easily copied around. It also makes it possible to re-grompp without having any special itp files available.

Arguments :

topol

topology file

struct

coordinat (structure) file

Keywords :

processed

name of the new topology file; if not set then it is named like topol but with pp_ prepended

includes

path or list of paths of directories in which itp files are searched for

grompp_kwargs*

other options for grompp such as maxwarn=2 can also be supplied

Returns :

full path to the processed topology

gromacs.cbook.edit_mdp(mdp, new_mdp=None, extend_parameters=None, **substitutions)

Change values in a Gromacs mdp file.

Parameters and values are supplied as substitutions, eg nsteps=1000.

By default the template mdp file is overwritten in place.

If a parameter does not exist in the template then it cannot be substituted and the parameter/value pair is returned. The user has to check the returned list in order to make sure that everything worked as expected. At the moment it is not possible to automatically append the new values to the mdp file because of ambiguities when having to replace dashes in parameter names with underscores (see the notes below on dashes/underscores).

If a parameter is set to the value None then it will be ignored.

Arguments :

mdp : filename

filename of input (and output filename of new_mdp=None)

new_mdp : filename

filename of alternative output mdp file [None]

extend_parameters : string or list of strings

single parameter or list of parameters for which the new values should be appended to the existing value in the mdp file. This makes mostly sense for a single parameter, namely ‘include’, which is set as the default. Set to [] to disable. [‘include’]

substitutions

parameter=value pairs, where parameter is defined by the Gromacs mdp file; dashes in parameter names have to be replaced by underscores. If a value is a list-like object then the items are written as a sequence, joined with spaces, e.g.

ref_t=[310,310,310] --->  ref_t = 310 310 310

Returns :

Dict of parameters that have not been substituted.

Example

edit_mdp('md.mdp', new_mdp='long_md.mdp', nsteps=100000, nstxtcout=1000, lincs_iter=2)

Note

• Dashes in Gromacs mdp parameters have to be replaced by an underscore when supplied as python keyword arguments (a limitation of python). For example the MDP syntax is lincs-iter = 4 but the corresponding keyword would be lincs_iter = 4.
• If the keyword is set as a dict key, eg mdp_params['lincs-iter']=4 then one does not have to substitute.
• Parameters aa_bb and aa-bb are considered the same (although this should not be a problem in practice because there are no mdp parameters that only differ by a underscore).
• This code is more compact in Perl as one can use s/// operators: s/^(\s*${key}\s*=\s*).*/$1${val}/

See also

One can also load the mdp file with gromacs.formats.MDP, edit the object (a dict), and save it again.

gromacs.cbook.add_mdp_includes(topology=None, kwargs=None)

Set the mdp include key in the kwargs dict.

1. Add the directory containing topology.
2. Add all directories appearing under the key includes
3. Generate a string of the form “-Idir1 -Idir2 ...” that is stored under the key include (the corresponding mdp parameter)

By default, the directories . and .. are also added to the include string for the mdp; when fed into gromacs.cbook.edit_mdp() it will result in a line such as

include = -I. -I.. -I../topology_dir ....

Note that the user can always override the behaviour by setting the include keyword herself; in this case this function does nothing.

If no kwargs were supplied then a dict is generated with the single include entry.

Arguments :

topology : top filename

Topology file; the name of the enclosing directory is added to the include path (if supplied) [None]

kwargs : dict

Optional dictionary of mdp keywords; will be modified in place. If it contains the includes keyword with either a single string or a list of strings then these paths will be added to the include statement.

Returns :

kwargs with the include keyword added if it did not exist previously; if the keyword already existed, nothing happens.

Note

The kwargs dict is modified in place. This function is a bit of a hack. It might be removed once all setup functions become methods in a nice class.

gromacs.cbook.grompp_qtot(*args, **kwargs)

Run gromacs.grompp and return the total charge of the system.

Arguments :The arguments are the ones one would pass to gromacs.grompp(). Returns :The total charge as reported

Some things to keep in mind:

• The stdout output of grompp is only shown when an error occurs. For debugging, look at the log file or screen output and try running the normal gromacs.grompp() command and analyze the output if the debugging messages are not sufficient.

• Check that qtot is correct. Because the function is based on pattern matching of the informative output of grompp it can break when the output format changes. This version recognizes lines like

  '  System has non-zero total charge: -4.000001e+00'
  

  using the regular expression System has non-zero total charge: *(?P<qtot>[-+]?d*.d+([eE][-+]d+)?).

Working with index files

Manipulation of index files (ndx) can be cumbersome because the make_ndx program is not very sophisticated (yet) compared to full-fledged atom selection expression as available in Charmm, VMD, or MDAnalysis. Some tools help in building and interpreting index files.

See also

The gromacs.formats.NDX class can solve a number of index problems in a cleaner way than the classes and functions here.

class gromacs.cbook.IndexBuilder(struct=None, selections=None, names=None, name_all=None, ndx=None, out_ndx='selection.ndx', offset=0)

Build an index file with specified groups and the combined group.

This is not a full blown selection parser a la Charmm, VMD or MDAnalysis but a very quick hack.

Example

How to use the IndexBuilder:

G = gromacs.cbook.IndexBuilder('md_posres.pdb',
              ['S312:OG','T313:OG1','A38:O','A309:O','@a62549 & r NA'],
              offset=-9, out_ndx='selection.ndx')
groupname, ndx = G.combine()
del G

The residue numbers are given with their canonical resids from the sequence or pdb. offset=-9 says that one calculates Gromacs topology resids by subtracting 9 from the canonical resid.

The combined selection is OR ed by default and written to selection.ndx. One can also add all the groups in the initial ndx file (or the make_ndx default groups) to the output (see the defaultgroups keyword for IndexBuilder.combine()).

Generating an index file always requires calling combine() even if there is only a single group.

Deleting the class removes all temporary files associated with it (see IndexBuilder.indexfiles).

Raises :

If an empty group is detected (which does not always work) then a gromacs.BadParameterWarning is issued.

Bugs :

If make_ndx crashes with an unexpected error then this is fairly hard to diagnose. For instance, in certain cases it segmentation faults when a tpr is provided as a struct file and the resulting error messages becomes

GromacsError: [Errno -11] Gromacs tool failed
Command invocation: make_ndx -o /tmp/tmp_Na1__NK7cT3.ndx -f md_posres.tpr

In this case run the command invocation manually to see what the problem could be.

See also

In some cases it might be more straightforward to use gromacs.formats.NDX.

Build a index group from the selection arguments.

If selections and a structure file are supplied then the individual selections are constructed with separate calls to gromacs.make_ndx(). Use IndexBuilder.combine() to combine them into a joint selection or IndexBuilder.write() to simply write out the individual named selections (useful with names).

Arguments :

struct : filename

Structure file (tpr, pdb, ...)

selections : list

The list must contain strings or tuples, which must be be one of the following constructs:

“<1-letter aa code><resid>[:<atom name]”

Selects the CA of the residue or the specified atom name.

example: "S312:OA" or "A22" (equivalent to "A22:CA")

(“<1-letter aa code><resid>”, “<1-letter aa code><resid>, [“<atom name>”])

Selects a range of residues. If only two residue identifiers are provided then all atoms are selected. With an optional third atom identifier, only this atom anme is selected for each residue in the range. [EXPERIMENTAL]

“@<make_ndx selection>”

The @ letter introduces a verbatim make_ndx command. It will apply the given selection without any further processing or checks.

example: "@a 6234 - 6238" or '@"SOL"' (note the quoting) or "@r SER & r 312 & t OA".

names : list

Strings to name the selections; if not supplied or if individuals are None then a default name is created. When simply using IndexBuilder.write() then these should be supplied.

name_all : string

Name of the group that is generated by IndexBuilder.combine().

offset : int, dict

This number is added to the resids in the first selection scheme; this allows names to be the same as in a crystal structure. If offset is a dict then it is used to directly look up the resids.

ndx : filename or list of filenames

Optional input index file(s).

out_ndx : filename

Output index file.

combine(name_all=None, out_ndx=None, operation='|', defaultgroups=False)

Combine individual groups into a single one and write output.

Keywords :

name_all : string

Name of the combined group, None generates a name. [None]

out_ndx : filename

Name of the output file that will contain the individual groups and the combined group. If None then default from the class constructor is used. [None]

operation : character

Logical operation that is used to generate the combined group from the individual groups: “|” (OR) or “&” (AND); if set to False then no combined group is created and only the individual groups are written. [“|”]

defaultgroups : bool

True: append everything to the default groups produced by make_ndx (or rather, the groups provided in the ndx file on initialization — if this was None then these are truly default groups); False: only use the generated groups

Returns :

(combinedgroup_name, output_ndx), a tuple showing the actual group name and the name of the file; useful when all names are autogenerated.

Warning

The order of the atom numbers in the combined group is not guaranteed to be the same as the selections on input because make_ndx sorts them ascending. Thus you should be careful when using these index files for calculations of angles and dihedrals. Use gromacs.formats.NDX in these cases.

See also

IndexBuilder.write().

gmx_resid(resid)

Returns resid in the Gromacs index by transforming with offset.

gromacs.cbook.parse_ndxlist(output)

Parse output from make_ndx to build list of index groups:

groups = parse_ndxlist(output)

output should be the standard output from make_ndx, e.g.:

rc,output,junk = gromacs.make_ndx(..., input=('', 'q'), stdout=False, stderr=True)

(or simply use

rc,output,junk = cbook.make_ndx_captured(...)

which presets input, stdout and stderr; of course input can be overriden.)

Returns :

The function returns a list of dicts (groups) with fields

name

name of the groups

nr

number of the group (starts at 0)

natoms

number of atoms in the group

gromacs.cbook.get_ndx_groups(ndx, **kwargs)

Return a list of index groups in the index file ndx.

Arguments :

• ndx is a Gromacs index file.
• kwargs are passed to make_ndx_captured().

Returns :

list of groups as supplied by parse_ndxlist()

Alternatively, load the index file with gromacs.formats.NDX for full control.

gromacs.cbook.make_ndx_captured(**kwargs)

make_ndx that captures all output

Standard make_ndx() command with the input and output pre-set in such a way that it can be conveniently used for parse_ndxlist().

Example::

ndx_groups = parse_ndxlist(make_ndx_captured(n=ndx)[0])

Note that the convenient get_ndx_groups() function does exactly that and can probably used in most cases.

Arguments :keywords are passed on to make_ndx() Returns :(returncode, output, None)

File editing functions

It is often rather useful to be able to change parts of a template file. For specialized cases the two following functions are useful:

gromacs.cbook.edit_mdp(mdp, new_mdp=None, extend_parameters=None, **substitutions)

Change values in a Gromacs mdp file.

Parameters and values are supplied as substitutions, eg nsteps=1000.

By default the template mdp file is overwritten in place.

If a parameter does not exist in the template then it cannot be substituted and the parameter/value pair is returned. The user has to check the returned list in order to make sure that everything worked as expected. At the moment it is not possible to automatically append the new values to the mdp file because of ambiguities when having to replace dashes in parameter names with underscores (see the notes below on dashes/underscores).

If a parameter is set to the value None then it will be ignored.

Arguments :

mdp : filename

filename of input (and output filename of new_mdp=None)

new_mdp : filename

filename of alternative output mdp file [None]

extend_parameters : string or list of strings

single parameter or list of parameters for which the new values should be appended to the existing value in the mdp file. This makes mostly sense for a single parameter, namely ‘include’, which is set as the default. Set to [] to disable. [‘include’]

substitutions

parameter=value pairs, where parameter is defined by the Gromacs mdp file; dashes in parameter names have to be replaced by underscores. If a value is a list-like object then the items are written as a sequence, joined with spaces, e.g.

ref_t=[310,310,310] --->  ref_t = 310 310 310

Returns :

Dict of parameters that have not been substituted.

Example

edit_mdp('md.mdp', new_mdp='long_md.mdp', nsteps=100000, nstxtcout=1000, lincs_iter=2)

Note

• Dashes in Gromacs mdp parameters have to be replaced by an underscore when supplied as python keyword arguments (a limitation of python). For example the MDP syntax is lincs-iter = 4 but the corresponding keyword would be lincs_iter = 4.
• If the keyword is set as a dict key, eg mdp_params['lincs-iter']=4 then one does not have to substitute.
• Parameters aa_bb and aa-bb are considered the same (although this should not be a problem in practice because there are no mdp parameters that only differ by a underscore).
• This code is more compact in Perl as one can use s/// operators: s/^(\s*${key}\s*=\s*).*/$1${val}/

See also

One can also load the mdp file with gromacs.formats.MDP, edit the object (a dict), and save it again.

gromacs.cbook.edit_txt(filename, substitutions, newname=None)

Primitive text file stream editor.

This function can be used to edit free-form text files such as the topology file. By default it does an in-place edit of filename. If newname is supplied then the edited file is written to newname.

Arguments :

filename

input text file

substitutions

substitution commands (see below for format)

newname

output filename; if None then filename is changed in place [None]

substitutions is a list of triplets; the first two elements are regular expression strings, the last is the substitution value. It mimics sed search and replace. The rules for substitutions:

substitutions        ::=  "[" search_replace_tuple, ... "]"
search_replace_tuple ::=  "(" line_match_RE "," search_RE "," replacement ")"
line_match_RE        ::=  regular expression that selects the line (uses match)
search_RE            ::=  regular expression that is searched in the line
replacement          ::=  replacement string for search_RE

Running edit_txt() does pretty much what a simple

sed /line_match_RE/s/search_RE/replacement/

with repeated substitution commands does.

Special replacement values: - None: the rule is ignored - False: the line is deleted (even if other rules match)

Note

• No sanity checks are performed and the substitutions must be supplied exactly as shown.
• All substitutions are applied to a line; thus the order of the substitution commands may matter when one substitution generates a match for a subsequent rule.
• If replacement is set to None then the whole expression is ignored and whatever is in the template is used. To unset values you must provided an empty string or similar.
• Delete a matching line if replacement=``False``.