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.
Simple commands with new default options so that they solve a specific problem (see also Manipulating trajectories and structures):
Computes the RMSD of the “Backbone” atoms after fitting to the “Backbone” (including both translation and rotation).
Standard invocations for manipulating trajectories.
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.
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
Center everything and make a compact representation (pass 1) and fit the system to a reference (pass 2).
| Keywords : |
|
|---|
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
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
- moved to the partsdir (under dirname)
- concatenated with the Gromacs tools to yield prefix.xtc, prefix.trr, prefix.edr, prefix.gro (or prefix.md) in dirname
- 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 : |
|
|---|
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 : |
|
|---|
Unordered list of all frames currently held on disk.
Clean up all temporary frames (which can be HUGE).
Delete all frames.
Extract frames from the trajectory to the temporary directory.
Holds the current frame number of the currently extracted batch of frames. Increases when iterating.
Total number of frames read so far; only important when maxframes > 0 is used.
Class to handle transformations of trajectories.
Set up Transformer with structure and trajectory.
Supply n = tpr, f = xtc (and n = ndx) relative to dirname.
| Keywords : |
|
|---|
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 : |
|
|---|---|
| Returns : | dictionary with keys tpr, xtc, which are the names of the the new files |
Write xtc that is fitted to the tpr reference structure.
See gromacs.cbook.trj_xyfitted() for details and description of kwargs. The most important ones are listed here but in most cases the defaults should work.
| Keywords : |
|
|---|---|
| Returns : | dictionary with keys tpr, xtc, which are the names of the the new files |
Write xtc and tpr only containing the protein.
| Keywords : |
|
|---|---|
| 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.)
Return canonical path to file under dirname with components args
If args form an absolute path then just return it as the absolute path.
Write xtc and tpr with water (by resname) removed.
| Keywords : |
|
|---|---|
| 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.)
Return the volume in nm^3 of structure file f.
(Uses gromacs.editconf(); error handling is not good)
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.
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:
Return the volume in nm^3 of structure file f.
(Uses gromacs.editconf(); error handling is not good)
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
|
|---|
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 : |
|
|---|---|
| Keywords : |
|
| Returns : | full path to the processed topology |
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 : |
|
|---|---|
| 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
See also
One can also load the mdp file with gromacs.formats.MDP, edit the object (a dict), and save it again.
Set the mdp include key in the kwargs dict.
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 : |
|
|---|---|
| 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.
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:
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.
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 GThe 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.
| Arguments : |
|
|---|
Combine individual groups into a single one and write output.
| Keywords : |
|
|---|---|
| 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.
Returns resid in the Gromacs index by transforming with offset.
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
|
|---|
Return a list of index groups in the index file ndx.
| Arguments : |
|
|---|---|
| Returns : | list of groups as supplied by parse_ndxlist() |
Alternatively, load the index file with gromacs.formats.NDX for full control.
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().
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) |
It is often rather useful to be able to change parts of a template file. For specialized cases the two following functions are useful:
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 : |
|
|---|---|
| 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
See also
One can also load the mdp file with gromacs.formats.MDP, edit the object (a dict), and save it again.
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 : |
|
|---|
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