Modules

pmx is a essentially a collection of classes and functions that allow the manipulation of molecular structure and topology files. It has been thought to be mainly used in conjunction with the Gromacs simulation package. In fact, pmx makes use of some Gromacs functions to read/write trajectory data and to allow fast neighborsearching from a python script. With time, pmx has also become a library for alchemical free energy calculations. Here below are the main Modules that the user might want to use to facilitate the setup and analysis of alchemical free energy calculations in Gromacs.

pmx.alchemy

This module contains functions that can be used to setup the hybrid structures and topologies needed for alchemical free energy calculations.

class pmx.alchemy.AbsRestraints(protein=None, ligand=None, complex=None, ligname=None, pro_ids=None, lig_ids=None, kbond=4184, kangle=41.84, kdihedral=41.84, T=298.15, seed=None)

Bases: object

Identifies a set of restraints as defined by Boresch between the protein/host and ligand/guest molecule. 5

(Pro)  p3 -- p2          l2   (Lig)
              \         /  \
               p1 --- l1    l3

The process of choosing atoms involved in the restraints is stochastic. Thus, every time you call this function, you might get a different set of restraints. If you would like to have a deterministic behaviour provide a seed.

If you have a structure of the protein-ligand complex and want to choose restraints based on this, you can do so by providing a Model of the complex as well as the residue name of the ligand (ligname).

If you know which atoms to use for the restraints, you can provide their atom indices as they are in the input Model objects with the arguments pro_ids and lig_ids. The expected order of atoms is l1, l2, l3 as in the scheme above for the ligand, and p1, p2, p3 for the protein.

Parameters

  • protein (Model, optional) – model of the protein/host molecule.

  • ligand (Model, optional) – model of the ligand/guest molecule.

  • complex (Model, optional) – model of the protein-ligand or host-guest complex.

  • ligname (str, optional) – name of ligand/host residue. This is required if a complex is provided.

  • pro_ids (list, optional) – list with three indeces.

  • lig_ids (list, optional) – list with three indeces.

  • kbond (float, optional) – force constant to use for the distance restraint. Default is 4184 kJ/(mol * nm^2).

  • kangle (float, optional) – force constant to use for the angle restraints. Default is 41.84 kJ/(mol * rad^2).

  • kdihedral (float, optional) – force constant to use for the dihedral restraints. Default is 41.84 kJ/(mol * rad^2).

  • T (float, optional) – temperature in Kelvin. Default is 298.15 K.

  • seed (bool, optional) – random seed.

lig_atoms

list of three ligand/guest Atoms selected.

Type

list

pro_atoms

list of three protein/host Atoms selected.

Type

list

dist

distance l1-p1 in nm.

Type

float

angle1

angle l2-l1-p1 in degrees.

Type

float

angle2

angle l1-p1-p2 in degrees.

Type

float

dihedral1

dihedral l3-l2-l1-p1 in degrees.

Type

float

dihedral2

dihedral l2-l1-p1-p2 in degrees.

Type

float

dihedral3

dihedral l1-p1-p2-p3 in degrees.

Type

float

dg

free energy of restraining the ligand while decoupled. This is calculate with equation 32 in Boresch et al. 5

Type

float

calc_dg(release=False)

Calculates the free energy contribution of the restraints given the distance, angles, and dihedrals and their force constants as defined in the instance.

Parameters

release (bool, optional) – whether to calculate the free energy of adding (False) or releasing (True) the restraints.

make_ii(switch_on=True, ligand_first=True)

Returns the data needed to add an intermolecular interactions section in a topology file to be used as restraints.

Parameters

  • switch_on (bool, optional) – whether you want to switch the restraints on (restraints present in state B) or off (restraints present in state A). Default is True (switching on).

  • ligand_first (bool, optional) – whether in the input gro/pdb file of the protein-ligand complex you are preparing the ligand atoms will come before (True) or after (False) the protein atoms. This is needed in order to adjust the atoms indices correctly. If you provided a Model of the complex (self.is_complex is True), this argument is disregarded as indices as provided will be used.

Returns

ii – dictionary with “bonds”, “angles”, and “dihedrals” keys with intermolecular_interactions information that can be passed to a Topology object.

Return type

dict

select_restraints()

Automated restraints selection. For the ligand/guest, three heavy atoms close to its center of mass are selected. For the protein, the backbone N-CA-C atoms closest to the selected ligand atoms are picked. If the protein/host is another small molecule, then three atoms close to the ligand atoms are selected.

write_summary(fname='restraints.info')

Write restraints information, including free energy of restraining, to file.

Parameters

fname (str, optional) – name of output file. Default is “restraints.info”.

pmx.alchemy.gen_hybrid_top(topol, recursive=True, verbose=False, scaleDih=1.0)

Fills the bstate of a topology file containing pmx hybrid residues. This can be either a top or itp file. If the file contains other itp files via include statements, the function can iterate through them if the recursive flag is set to True.

It returns a Topology instance with b states filled for the input topol instance, and optionally a list of additional Topology instances for the included itp files. If there are no itp files, or if recursive is set to False, then an empty list is returned.

Parameters

  • topol (Topology) – input Topology instance

  • recursive (bool, optional) – whether to fill b states also for all itp files included in the topology file. Default is True.

  • verbose (bool) – whether to print out info. Default is False.

  • scaleDih (float) – scaling for dihedrals with a dummy.

Returns

top, itps – output Topology instance with b states assigned (top), and list of Topology instances for the additional itp files (itps) included in the input topology.

Return type

obj, list of obj

pmx.alchemy.mutate(m, mut_resid, mut_resname, ff, mut_chain=None, refB=None, inplace=False, verbose=False)

Creates an hybrid structure file.

Parameters

  • m (Model) – the model to be mutated. See pmx.model.Model.

  • mut_resid (int) – the ID of the residue to be mutated.

  • mut_resname (str) – the name of the target residue.

  • ff (str) – the forcefield to use.

  • mut_chain (str, optional) – the chain ID for the residue you want to mutate. This is needed if you have multiple chains and are providing a Model that did not have its residues renumbered.

  • refB (str, optional) – reference structure file of the B state in PDB or GRO format. If it is provided, the dummy atoms will be built based on the position of the atoms in this structure. This option is available only for Protein mutations.

  • inplace (bool, optional) – whether to modify the input Model in place. Default is False.

  • verbose (bool, optional) – whether to print info. Default is False.

Returns

m2 – model with mutated residues. If argument ‘inplace’ is True, nothing is returned.

Return type

Model

pmx.alchemy.write_split_top(pmxtop, outfile='pmxtop.top', scale_mass=False, verbose=False)

Write three topology files to be used for three separate free energy calculations: charges off, vdw on, changes on. This can be useful when one wants to avoid using a soft-core for the electrostatic interactions.

Parameters

  • pmxtop (Topology) – Topology object with B states defined.

  • outfile (str) – name of output topology file.

  • scale_mass (bool) – whether to scale the masses of morphing atoms

  • verbose (bool) – whether to print information or not

Returns

Return type

None

pmx.analysis

pmx.analysis.integrate_dgdl(fn, ndata=- 1, lambda0=0, invert_values=False, sigmoid=0.0)

Integrates the data in a dgdl.xvg file.

Parameters

  • fn (str) – the inpur dgdl.xvg file from Gromacs.

  • ndata (int, optional) – number of datapoints in file. If -1, then ??? default is -1.

  • lambda0 ([0,1]) – whether the simulations started from lambda 0 or 1. Default is 0.

  • invert_values (bool) – whether to invert the sign of the returned work value.

Returns

  • integr (float) – result of the integration performed using Simpson’s rule.

  • ndata (int) – number of data points in the input file.

pmx.analysis.ks_norm_test(data, alpha=0.05, refks=None)

Performs a Kolmogorov-Smirnov test of normality.

Parameters

  • data (array_like) – a one-dimensional array of values. This is the distribution tested for normality.

  • alpha (float) – significance level of the statistics. Default if 0.05.

  • refks –

    ???

Returns

  • Q (float)

  • lam0 (float)

  • check (float)

  • bOk (bool)

pmx.analysis.plot_work_dist(wf, wr, fname='Wdist.png', nbins=20, dG=None, dGerr=None, units='kJ/mol', dpi=300)

Plots forward and reverse work distributions. Optionally, it adds the estimate of the free energy change and its uncertainty on the plot.

Parameters

  • wf (list) – list of forward work values.

  • wr (list) – list of reverse work values.

  • fname (str, optional) – filename of the saved image. Default is ‘Wdist.png’.

  • nbins (int, optional) – number of bins to use for the histogram. Default is 20.

  • dG (float, optional) – free energy estimate.

  • dGerr (float, optional) – uncertainty of the free energy estimate.

  • units (str, optional) – the units of dG and dGerr. Default is ‘kJ/mol’.

  • dpi (int) – resolution of the saved image file.

Returns

Return type

None

pmx.analysis.read_dgdl_files(lst, lambda0=0, invert_values=False, verbose=True, sigmoid=0.0)

Takes a list of dgdl.xvg files and returns the integrated work values.

Parameters

  • lst (list) – list containing the paths to the dgdl.xvg files.

  • lambda0 ([0,1]) – whether the simulations started from lambda 0 or 1. Default is 0.

  • invert_values (bool) – whether to invert the sign of the returned work value.

Returns

  • w (array) – array of work values.

  • lst (list) – sorted list of input dgdl.avg files corresponding to the work values in w.

pmx.estimators

class pmx.estimators.BAR(wf, wr, T, nboots=0, nblocks=1)

Bases: object

Bennett acceptance ratio (BAR) estimator. 3

The standard error is calculated via bootstrap when nboots>0, and by separating the work values into groups when nblocks>1.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • T (float, optional) – temperature in Kelvin. Default is 298.15 K.

  • nboots (int, optional) – how many bootstrap samples to draw for estimating the standard error. Default is zero (do not estimate the error).

  • nblocks (int, optional) – how many blocks to divide the input work values into for the estimation of the standard error. Default is one (do not estimate the error).

Examples

>>> estimate = BAR(wf, wr, T=300, nboots=1000, nblocks=10)
>>> dg = estimate.dg
>>> dg_err1 = estimate.err
>>> dg_err2 = estimate.err_boot
>>> dg_err3 = estimate.err_blocks
>>> dg_convergece = estimate.conv
dg

the free energy estimate.

Type

float

err

analytical estimate of the standard error.

Type

float

err_boot

standard error of the free energy estimate calculated via bootstrap.

Type

float

err_blocks

standard error of the free energy estimate calculated by separating the input work values into groups/blocks.

Type

float

conv

convergence of the free energy estimate. 4

Type

float

conv_err_boot

standard error of the convergence estimate calculated via bootstrap.

Type

float

static calc_conv(dg, wf, wr, T)

Evaluates BAR convergence as described by Hahn & Then. 4

Returns a value between -1 and 1: the closer this value to zero the better the BAR convergence.

Parameters

  • dg (float) – the BAR free energy estimate

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • T (float) – temperature

Returns

conv – convergence estimate.

Return type

float

static calc_conv_err_boot(dg, wf, wr, nboots, T)

Calculates the error in the convergence measure by bootstrapping.

Parameters

  • dg (float) – the BAR free energy estimate

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • T (float) – temperature

  • nboots (int) – number of bootstrap samples.

Returns

sderr – the standard error of the convergence measure.

Return type

float

static calc_dg(wf, wr, T)

Estimates and returns the free energy difference.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • T (float) – temperature

Returns

dg – the BAR free energy estimate.

Return type

float

static calc_err(dg, wf, wr, T)

Calculates the analytical error estimate.

Parameters

  • dg (float) – the BAR free energy estimate

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • T (float) – temperature

Returns

sderr – the standard error of the estimate.

Return type

float

static calc_err_blocks(wf, wr, nblocks, T)

Calculates the standard error based on a number of blocks the work values are divided into. It is useful when you run independent equilibrium simulations, so that you can then use their respective work values to compute the standard error based on the repeats.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • T (float) – temperature

  • nblocks (int) – number of blocks to divide the data into. This can be for instance the number of independent equilibrium simulations you ran.

Returns

sderr – the standard error of the estimate.

Return type

float

static calc_err_boot(wf, wr, nboots, T)

Calculates the error by bootstrapping.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • T (float) – temperature

  • nboots (int) – number of bootstrap samples.

Returns

sderr – the standard error of the estimate.

Return type

float

class pmx.estimators.Crooks(wf, wr, nboots=0, nblocks=1)

Bases: object

Crooks Gaussian Intersection (CGI) estimator. 2

The forward and reverse work values are fitted to Gaussian functions and their intersection is taken as the free energy estimate. In some cases, when the two Gaussians are very close to each other, the intersection cannot be taken and the average of the two Gaussian means is taken as the free energy estimate insted.

The standard error is calculated via bootstrap when nboots>0, and by separating the work values into groups when nblocks>1. Bootstrap errors are calculating with two different bootstrap procedures: a parametric one in which samples are drawn from Gaussian distributions fitted to the original work distributions, and a nonparametric one in which the work values themself are drawn at random with replacement.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • nboots (int, optional) – how many bootstrap samples to draw for estimating the standard error. Default is zero (do not estimate the error).

  • nblocks (int, optional) – how many blocks to divide the input work values into for the estimation of the standard error. Default is one (do not estimate the error).

Examples

>>> estimate = Crooks(wf, wr, nboots=1000, nblocks=10)
>>> dg = estimate.dg
>>> dg_err1 = estimate.err_boot1
>>> dg_err2 = estimate.err_boot2
>>> dg_err3 = estimate.err_blocks
dg

the free energy estimate.

Type

float

err_boot1

standard error of the free energy estimate calculated via parametric bootstrap.

Type

float

err_boot2

standard error of the free energy estimate calculated via nonparametric bootstrap.

Type

float

err_blocks

standard error of the free energy estimate calculated by separating the input work values into groups/blocks.

Type

float

inters_bool

whether the interection could be taken. If False, the free energy estimate is the average of the two Gaussian means.

Type

bool

Af

height of the forward Gaussian.

Type

float

mf

mean of the forward Gaussian.

Type

float

devf

standard deviation of the forward Gaussian.

Type

float

Ar

height of the reverse Gaussian.

Type

float

mr

mean of the reverse Gaussian.

Type

float

devr

standard deviation of the reverse Gaussian.

Type

float

static calc_dg(wf, wr)

Calculates the free energy difference using the Crooks Gaussian Intersection method. It finds the intersection of two Gaussian functions. If the intersection cannot be computed, the average of the two Gaussian locations is returned.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

Returns

  • float – location of the intersection.

  • bool – whether the intersection could be calculated. If the intersection was calculated as expected a True value is returned. If the Gaussians are too close to each other, the intersection cannot be calculated and a False value is returned; in this case, the first float value retured is the average of the Gaussian means.

static calc_err_blocks(wf, wr, nblocks)

Calculates the standard error based on a number of blocks the work values are divided into. It is useful when you run independent equilibrium simulations, so that you can then use their respective work values to compute the standard error based on the repeats.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • nblocks (int) – number of blocks to divide the data into. This can be for instance the number of independent equilibrium simulations you ran.

Returns

sderr – the standard error of the estimate.

Return type

float

static calc_err_boot1(m1, s1, n1, m2, s2, n2, nboots)

Calculates the standard error of the Crooks Gaussian Intersection via parametric bootstrap. Given the parameters of the forward and reverse Gaussian distributions, multiple (nboots) bootstrap samples are built by random sampling from these two Gaussian distributions. The CGI free energy is then calculated for each bootstrap sample (forward and reverse Gaussians). The standard error of the estimate is returned as the standard deviation of the bootstrapped free energies.

Parameters

  • m1 (float) – mean of the forward Gaussian.

  • s1 (float) – standard deviation of the forward Gaussian.

  • n1 (int) – number of work values to which the first Gaussian was fit.

  • m2 (float) – mean of the reverse Gaussian.

  • s2 (float) – standard deviation of the reverse Gaussian.

  • n2 (int) – number of work values to which the second Gaussian was fit.

  • nboots (int) – number of bootstrap samples to use for the error estimate. Parametric bootstrap is used where work values are resampled from two Gaussians.

Returns

sderr – the standard error of the estimate.

Return type

float

static calc_err_boot2(wf, wr, nboots)

Calculates the standard error of the Crooks Gaussian Intersection via non-parametric bootstrap. The work values are resampled randomly with replacement multiple (nboots) times, and the CGI free energy recalculated for each bootstrap samples. The standard error of the estimate is returned as the standard deviation of the bootstrapped free energies.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • nboots (int) – number of bootstrap samples to use for the error estimate.

Returns

sderr – the standard error of the estimate.

Return type

float

class pmx.estimators.Jarz(wf, wr, T=298.15, nboots=0, nblocks=1)

Bases: object

Jarzynski estimator. 1

Both the forward and reverse estimates, as well as the mean of the two are calculated.

The standard error is calculated via bootstrap when nboots>0, and by separating the work values into groups when nblocks>1.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • T (float, optional) – temperature in Kelvin. Default is 298.15 K.

  • nboots (int, optional) – how many bootstrap samples to draw for estimating the standard error. Default is zero (do not estimate the error).

  • nblocks (int, optional) – how many blocks to divide the input work values into for the estimation of the standard error. Default is one (do not estimate the error).

Examples

>>> estimate = Jarz(wf, wr, T=300, nboots=1000, nblocks=10)
>>> dg_forward = estimate.dg_for
>>> dg_reverse = estimate.dg_rev
>>> dg_mean = estimate.dg_mean
>>> dg_forward_err1 = estimate.err_boot_for
>>> dg_forward_err2 = estimate.err_blocks_for
dg_for

the forward free energy estimate.

Type

float

dg_rev

the reverse free energy estimate.

Type

float

dg_mean

the mean of the forward and reverse free energy estimates.

Type

float

err_boot_for

standard error of the forward free energy estimate calculated via bootstrap.

Type

float

err_boot_rev

standard error of the reverse free energy estimate calculated via bootstrap.

Type

float

err_blocks_for

standard error of the forward free energy estimate calculated by separating the input work values into groups/blocks.

Type

float

err_blocks_rev

standard error of the reverse free energy estimate calculated by separating the input work values into groups/blocks.

Type

float

static calc_dg(w, T, bReverse=False)

Calculates the free energy difference using Jarzynski’s equality. 1

Parameters

  • w (array_like) – array of work values.

  • T (float) – temperature in Kelvin.

  • bReverse (bool, optional) – whether the work values provided are for the reverse transition. Default if False. If they are for the reverse transition, set it to True.

Returns

dg – estimate of the free energy difference.

Return type

float

static calc_err_blocks(w, T, nblocks, bReverse=False)

Calculates the standard error based on a number of blocks the work values are divided into. It is useful when you run independent equilibrium simulations, so that you can then use their respective work values to compute the standard error based on the repeats.

Parameters

  • w (array_like) – array of work values.

  • T (float) – temperature.

  • bReverse (bool, optional) – whether the work values provided are for the reverse transition. Default if False. If they are for the reverse transition, set it to True.

  • nblocks (int) – number of blocks to divide the data into. This can be for instance the number of independent equilibrium simulations you ran.

Returns

sderr – the standard error of the estimate.

Return type

float

static calc_err_boot(w, T, nboots, bReverse=False)

Calculates the standard error via bootstrap. The work values are resampled randomly with replacement multiple (nboots) times, and the Jarzinski free energy recalculated for each bootstrap samples. The standard error of the estimate is returned as the standard d eviation of the bootstrapped free energies.

Parameters

  • w (array_like) – work values.

  • T (float) – temperature.

  • bReverse (bool, optional) – whether the work values provided are for the reverse transition. Default if False. If they are for the reverse transition, set it to True.

  • nboots (int) – number of bootstrap samples to use for the error estimate.

Returns

sderr – the standard error of the estimate.

Return type

float

class pmx.estimators.JarzGauss(wf, wr, T=298.15, nboots=0, nblocks=1)

Bases: object

Jarzynski estimator using a Gaussian approximation. 6

Both the forward and reverse estimates.

The standard error is calculated using the analytical expression derived by Hummer (2001). 6 When nboots>0, the error is also calculated by bootstrap. When nblocks>1, it is calculated by separating the work values into blocks.

Parameters

  • wf (array_like) – array of forward work values.

  • wr (array_like) – array of reverse work values.

  • T (float, optional) – temperature in Kelvin. Default is 298.15 K.

  • nboots (int, optional) – how many bootstrap samples to draw for estimating the standard error. Default is zero (do not estimate the error).

  • nblocks (int, optional) – how many blocks to divide the input work values into for the estimation of the standard error. Default is one (do not estimate the error).

Examples

>>> estimate = JarzGauss(wf, wr, T=300, nboots=1000, nblocks=10)
>>> dg_forward = estimate.dg_for
>>> dg_reverse = estimate.dg_rev
>>> dg_forward_err1 = estimate.err_for
>>> dg_forward_err2 = estimate.err_boot_for
>>> dg_forward_err3 = estimate.err_blocks_for
dg_for

the forward free energy estimate.

Type

float

dg_rev

the reverse free energy estimate.

Type

float

err_for

analytical estimate of the standard error of the forward free energy estimate.

Type

float

err_rev

analytical estimate of the standard error of the reverse free energy estimate.

Type

float

err_boot_for

standard error of the forward free energy estimate calculated via bootstrap.

Type

float

err_boot_rev

standard error of the reverse free energy estimate calculated via bootstrap.

Type

float

err_blocks_for

standard error of the forward free energy estimate calculated by separating the input work values into groups/blocks.

Type

float

err_blocks_rev

standard error of the reverse free energy estimate calculated by separating the input work values into groups/blocks.

Type

float

static calc_dg(w, T, bReverse=False)

Calculates the free energy difference using the Jarzynski estimator with a Gaussian approximation. 6

Parameters

  • w (array_like) – array of work values.

  • T (float) – temperature in Kelvin.

  • bReverse (bool, optional) – whether the work values provided are for the reverse transition. Default if False. If they are for the reverse transition, set it to True.

Returns

dg – estimate of the free energy difference.

Return type

float

static calc_err(w, T, bReverse=False)

Calculates the standard error via an analytic expression. The expression is derived by Hummer, 2001, JChemPhys. 6

Parameters

  • w (array_like) – work values.

  • T (float) – temperature.

  • bReverse (bool, optional) – whether the work values provided are for the reverse transition. Default if False. If they are for the reverse transition, set it to True.

Returns

err – standard deviation of the estimator.

Return type

float

static calc_err_blocks(w, T, nblocks, bReverse=False)

Calculates the standard error based on a number of blocks the work values are divided into. It is useful when you run independent equilibrium simulations, so that you can then use their respective work values to compute the standard error based on the repeats.

Parameters

  • w (array_like) – array of work values.

  • T (float) – temperature.

  • nblocks (int) – number of blocks to divide the data into. This can be for instance the number of independent equilibrium simulations you ran.

  • bReverse (bool, optional) – whether the work values provided are for the reverse transition. Default if False. If they are for the reverse transition, set it to True.

Returns

sderr – the standard error of the estimate.

Return type

float

static calc_err_boot(w, T, nboots, bReverse=False)

Calculates the standard error via bootstrap. The work values are resampled randomly with replacement multiple (nboots) times, and the Gaussian approximation for Jarzinski free energy is recalculated for each bootstrap samples. The standard error of the estimate is returned as the standard d eviation of the bootstrapped free energies.

Parameters

  • w (array_like) – work values.

  • T (float) – temperature.

  • nboots (int) – number of bootstrap samples to use for the error estimate.

  • bReverse (bool, optional) – whether the work values provided are for the reverse transition. Default if False. If they are for the reverse transition, set it to True.

Returns

err – standard error of the mean.

Return type

float

pmx.gmx

This module contains wrapper functions for some of the most used Gromacs tools.

pmx.gmx.editconf(f, o='editconf.gro', bt='cubic', d=1.2, other_flags='')

Simple gmx editconf wrapper.

Parameters

  • f (str) – input structure file

  • o (str, optional) – name of output structure file. Default is “solvate.gro”

  • bt (str) – box type:triclinic, cubic, dodecahedron, or octahedron

  • d (float) – distance between the solute and the box (nm)

  • other_flags (str, optional) – additional flags to pass as you would type them in the shell

Returns

Return type

None

pmx.gmx.genion(s, p, o='genion.gro', np=0, nn=0, conc=0.15, neutral=True, other_flags='')

Simple gmx genion wrapper. By default, group ‘SOL’ will be replaced by ions.

Parameters

  • s (str) – input tpr file

  • p (str) – input topology file.

  • o (str, optional) – name of output structure file. Default is “genion.gro”

  • np (int, optional) – number of positive ions. Default is 0.

  • nn (int, optional.) – number of negative ions. Default is 0.

  • conc (float, optional) – specify salt concentration (mol/liter). Default is 0.15 M.

  • neutral (bool) – whether to add enough ions to neutralise the system. These ions are added on top of those specified with -np/-nn or -conc. Default is True.

  • other_flags (str, optional) – additional flags to pass as you would type them in the shell

Returns

Return type

None

pmx.gmx.get_gmx()

Gets path to gmx executable, and throws error if not found.

pmx.gmx.grompp(f, c, p, o='grompp.tpr', maxwarn=0, other_flags='')

Simple gmx grompp wrapper.

Parameters

  • f (str) – input mdp file

  • c (str) – input structure file

  • p (str) – input topology file. Default is “topol.top”

  • o (str, optional) – output tpr file. Default is “grompp.tpr”

  • maxwarn (int, optional) – number of allowed warnings. Default is 0.

  • other_flags (str, optional) – additional flags to pass as you would type them in the shell

Returns

Return type

None

pmx.gmx.mdrun(s, deffnm='md', verbose=False, other_flags='')

Simple gmx mdrun wrapper.

Parameters

  • s (str) – input tpr file

  • deffnm (str, optional) – set the default filename for all file options. Default is ‘md’.

  • verbose (bool, optional) – whether to activate verbose flag in Gromacs mdrun. Default is False.

  • other_flags (str, optional) – additional flags to pass as you would type them in the shell.

Returns

Return type

None

pmx.gmx.pdb2gmx(f, o='pdb2gmx.gro', p='topol.top', ff='amber99sb-star-ildn-mut', water='tip3p', other_flags='')

Simple gmx pdb2gmx wrapper.

Parameters

  • f (str) – input structure file

  • o (str, optional) – output structure file. Default is “pdb2gmx.gro”.

  • p (str, optional) – name of topology file. Default is “topol.top”.

  • ff (str, optional) – forcefield. Default is “amber99sb-star-ildn-mut”.

  • water (str, optional) – water model. Default is “tip3p”.

  • other_flags (str, optional) – additional flags to pass as you would type them in the shell

Returns

Return type

None

pmx.gmx.set_gmxlib()

Sets the environment variable GMXLIB to the default/expected pmx location.

pmx.gmx.solvate(cp, cs='spc216.gro', p='topol.top', o='solvate.gro', other_flags='')

Simple gmx solvate wrapper.

Parameters

  • cp (str) – input structure file

  • cs (str, optional) – structure file of the solvent. Default is “spc216.gro”.

  • p (str, optional) – name of topology file. Default is “topol.top”

  • o (str, optional) – name of output structure file. Default is “solvate.gro”

  • other_flags (str, optional) – additional flags to pass as you would type them in the shell

Returns

Return type

None

pmx.gmx.trjconv(f, s, o='trjconv.xtc', ur='compact', pbc='none', fit='none', out_grp='System', fit_grp='C-alpha', sep=False, other_flags='')

Simple gmx trjconv wrapper.

Parameters

  • f (str) – input structure of trajectory file

  • s (str) – input tpr file

  • o (str, optional) – output trajectory/structure file. Default is “trjconv.xtc”

  • ur (str, optional) – unit-cell representation: rect, tric, compact. Default is ‘compact’.

  • pbc (str, optional) – PBC treatment: none, mol, res, atom, nojump, cluster, whole. Default is ‘none’.

  • fit (str, optional) – fit molecule to ref structure in the structure file: none, rot+trans, rotxy+transxy, translation, transxy, progressive. Default is ‘none’.

  • out_grp (str, optional) – output group. Defauls is ‘System’.

  • fit_grp (str, optional) – group to use for the fitting if ‘fit’ is not none. Default is ‘C-alpha’.

  • sep (bool, optional) – write each frame to a separate .gro, .g96 or .pdb file. Default is False.

  • other_flags (str, optional) – additional flags to pass as you would type them in the shell

Returns

Return type

None

pmx.gmx.write_mdp(mdp, fout='mdpfile.mdp', nsteps=10000, cutoff=1.0, T=300)

Writes a few standard mdp files.

With the argument mdp you can choose from a few standard predefined mdp file: “enmin” (energy minimisation); “npt” (simulation in NPT ensemble); “npt-restr” (simulation in NPT ensemble with -DPOSRES defined).

Parameters

  • mdp (str) – what type of mdp file to load and write. Options available are: ‘enmin’, ‘npt-restr’, ‘npt’.

  • fout (str) – filename of the mdp file to be written.

  • nsteps (int) – number of steps.

  • cutoff (float) – short-range cutoff (nm) for vdw and coulomb interactions.

  • T (float) – temperature in Kelvin.

Returns

Return type

None

pmx.builder

This module comtains some functions to build protein structures from scratch.

Most important: - build_chain

Usage:
>>> from pmx.builder import *
>>> ch = build_chain('AAAAAA')    # builds an ordinary polyalanine
>>> ch.write('polyala.pdb')       # store file
>>> ch = build_chain('FGHRTCV',hydrogens = False) # build chain without H
>>> ch = build_chain('FGHRTCV',ss='HHHHHHH') build as helix
>>> ch = build_chain('FGHRTCV',dihedrals = ((phi1,psi1,omega1),(...)))
build chain with defined dihedral angles
pmx.builder.add_bp(m, strand=None, bRNA=False)
pmx.builder.attach_aminoacid(mol, resname, hydrogens=True, omega=180.0, phi=- 139.0, psi=135.0, chirality='L')
pmx.builder.attach_group(atom, mol)
pmx.builder.build_chain(sequence, dihedrals=None, hydrogens=True, ss=None, chain_id='X', bCyclic=False, bCap=False, chirality=None)
pmx.builder.build_dna_strand(seq)
pmx.builder.build_rna_strand(seq)
pmx.builder.get_fragments()
pmx.builder.make_3ter(r)
pmx.builder.make_5ter(r)
pmx.builder.make_residue(key, hydrogens=True, chirality='L')

returns a molecule object with default geometry

pmx.builder.make_stereoisomer_residue(m)
pmx.builder.read_pdb_with_connect(f)
pmx.builder.set_omega(mol1, mol2, phi)
pmx.builder.set_phi(mol1, mol2, phi)
pmx.builder.set_psi(mol, phi)
pmx.builder.write_pdb_with_connect(mol, f, n=1)

pmx.utils

Various utility functions and classes

exception pmx.utils.MissingTopolParamError(s, atoms)

Bases: Exception

Raised when certain parameters for terms defined in the topology cannot ba found.

exception pmx.utils.UnknownResidueError(s)

Bases: Exception

Class for unknown residue Exceptions.

pmx.utils.create_folder(path, fname=False)
pmx.utils.data2gauss(data)

Takes a one dimensional array and fits a Gaussian.

Parameters

data (list) – 1D array of values

Returns

  • float – mean of the distribution.

  • float – standard deviation of the distribution.

  • float – height of the curve’s peak.

pmx.utils.doLog(fp, s, commandName='ligandHybridTop')
pmx.utils.ff_selection(gmxlib=None)
pmx.utils.ffopen(filename, mode='r', backup=True)
pmx.utils.gauss_func(A, mean, dev, x)

Given the parameters of a Gaussian and a range of the x-values, returns the y-values of the Gaussian function

pmx.utils.get_ff_path(ff, verbose=False)

Get path of force field.

Parameters

ff (str) – force field name

Returns

ff_path – absolute path to the force field

Return type

str

pmx.utils.get_mtp_file(residue, ff)

Returns the correct mutres mtp file depending on the nature of the residue.

Parameters

  • residue (Molecule instance) – the residue that will be mutated

  • ff (str) – the force field to use

Returns

mtp_file – path to the mtp file

Return type

str

pmx.utils.get_pmxdata(fname)

Returns the absolute path for the file with ‘fname’ inside the pmx data folder.

exception pmx.utils.importError(s)

Bases: Exception

Exceptions class for importing a module.

pmx.utils.initialise_logger(logname, logfile, level='INFO')
pmx.utils.killBackups(arg, dirname, fname)
pmx.utils.list2file(lst, fname)

Writes a list of strings to a file.

pmx.utils.listDirs(dir='./')

returns a list of directrories in directory dir

pmx.utils.listFiles(dir='./', ext=None, abs=True, backups=False)

returns a list of files in directory dir, optionally only certain file types

exception pmx.utils.mtpError(s)

Bases: Exception

Exceptions class for …

pmx.utils.multiple_replace(text, word_dic, isfile=False)

Takes a text and replace words that match a key in a dictionary with the associated value, returning the modified text.

Parameters

  • text (str) – input text

  • word_dic (dict) – dictionary where keys strings are replaced by their value strings

  • isfile (bool, optional) – if isfile is True, the input ‘text’ will be considered a textfile, which will be opened, the patterns replaced, and written again.

Returns

newtext – output text with strings replaced

Return type

str

Examples

>>> a = 'hello, John'
>>> multiple_replace(a, {'hello':'goodbye', 'John':'Mark'})
'goodbye, Mark'
pmx.utils.natural_sort(l)
pmx.utils.removeBackups(dir, check=True)
pmx.utils.remove_netmount(string)
pmx.utils.show_ff(gmxlib=None)

Prints the list of forcefields available in GMXLIB.

Parameters

gmxlib (str, optional) – Path to force field library. If not set explicitly, it is taken from the enviroment variable GMXLIB.

pmx.utils.which(program)

Equivalent of UNIX which command. Returns the path of the executable if it is found, otherwise returns None.

Parameters

program (str) – name of executable you are looking for

Returns

path – string of path if found, otherwise None.

Return type

str

pmx.atomselection

This module contains the Atomselection class. It contains methods that are inherited to the Molecule, Chain and Model classes. Take a look there for details…

class pmx.atomselection.Atomselection(**kwargs)

Bases: object

Basic class to handle sets of atoms. Atoms are stored in a list <atoms>

a2nm()
atomlistFromTop(topDic)

return a list of atom objects found in topDic

com(vector_only=False)

move atoms to center of mass or return vector only

coords()
fetch_atoms(key, how='byname', wildcard=False, inv=False)

Find atoms given their name, element, or id.

Parameters

  • key (str|int) – key used to find the atoms.

  • how (str, optional) – the way atoms are identified. Options are “byname”, “byelem”, or “byid”. Default is “byname”.

  • wildcard (bool, optional) – whether to loosely match the atom name. In this case, if the key provided is within an atom name, the atom will be selected. E.g. the key=’H1’ will select atoms named, for instance, “H12”, “H13”, and “HH1”.

  • inv (bool, optional) – whether to invert the selection. Default is False.

Returns

atoms – list of Atom objects found.

Return type

list

get_b13()
get_b14()
get_by_id(lst)
get_long_name()
get_order()
get_symbol()
make_mol2_bondlist()
max_crd()
nm2a()
random_rotation()
rename_atoms_to_gmx()

Renames atoms to comply with Gromacs syntax. If the name starts with a digit, the digit is moved at the end of the name. E.g. “1CA” becomes “CA1”.

renumber_atoms(start=1)

Renumber all atoms starting from a chosen number (default is 1). The original atom IDs are stored in the attribute orig_id

Parameters

start (int, optional) – integer from which to start the indexing of the atoms

search_neighbors(cutoff=8.0, build_bonds=True)
translate(vec)
write(fn, title='', nr=1)
writeGRO(filename, title='')
writePDB(fname, title='', nr=1, bPDBTER=False, bAssignChainIDs=False, resnrlist=[])

Writes Atomselection to file in PDB format.

Parameters

  • fname (str) – filename

  • title (str, optional) – title of the PDB file. If not given, it is taken from the instance title

  • nr (int, optional) – what is nr? it seems to choose whether to overwrite or append to a file?

  • bPDBTER (bool, optional) – whether to separate chains with TER entries.

  • bAssignChainIDs (bool, optional) – whether to write the chains IDs to file. Default is False.

  • resnrlist (list, optional) – list of residues to write to file. If empty list provided, all residues are written to file, otherwise only the ones in the list will be written.

pmx.ffparser

docs to be added

class pmx.ffparser.ATPParser(fn='ffamber99sb.atp')

Bases: object

parse()
class pmx.ffparser.BondedParser(filename=None, version='old')

Bases: object

assign_params(model)
get_angle_param(type1, type2, type3)
get_bond_param(type1, type2)
get_dihedral_param(type1, type2, type3, type4, func)
parse(filename)
class pmx.ffparser.NBParser(filename=None, version='new', ff='amber')

Bases: object

assign_params(model)
parse(filename, version)
class pmx.ffparser.RTPParser(filename=None)

Bases: object

add_entry(name, rtp_dic)
assign_dihedral_params(model, directives)
assign_params(model)
next()
parse(filename)
write(out_file)

pmx.geometry

This file contains the Rotation class. A Rotation instance is built with two vectors.

Examples

>>> v1 = [1,2,3]
>>> v2 = [2,3,4]
>>> r = Rotation(v1,v2)  # create rotation object around v2-v1
>>> v3 = [4,5,6]
>>> v3 = r.apply(v3)     # rotate v3 around v2-v1
class pmx.geometry.Rotation(v1, v2)

Bases: object

apply(v, phi)
pmx.geometry.apply_fit_R(atoms, R)
pmx.geometry.bb_super(mol1, mol2, use_orig_mc_coords=True)

superpose mol2 on mol1

pmx.geometry.calc_fit_R(cs1, cs2, m)
pmx.geometry.center_vector(v)
pmx.geometry.fit(model1, model2, atom_names=[])
pmx.geometry.fit_atoms(fit_atoms1, fit_atoms2, rot_atoms2)
pmx.geometry.fit_by_ndx(ref, model, ndx1, ndx2)
pmx.geometry.nuc_super(mol1, mol2, name1=None, name2=None)

superpose mol2 on mol1

pmx.geometry.planarity(atom_list)
pmx.geometry.translate_by_ndx(struct, ndx)
pmx.geometry.vec_ang(v1, v2)

pmx.library

Library for useful and less useful things needed by the pmx packages.

pmx.library.pmx_data_file(filename, verbose=False)

pmx.mutdb

Functions to read the mutation database

pmx.mutdb.read_mtp(filename='ffoplsaa.mtp')
pmx.mutdb.read_mtp_entry(entry, filename='ffamber99sb.mtp', version='old')
pmx.mutdb.read_mutpdb(filename='mutations_oplsaa.pdb')
pmx.mutdb.read_new_mtp_entry(entry, filename='mutres.mtp')

pmx.parser

Some functions to write parsers for different data formats.

Usage:

>>> lst = open(infile,'r').readlines()     # read file
>>> lst = kickOutComments(lst,'#')       # remove comments with # or ;
>>> lst = parseList('ifs',lst)            # parse each line and return
[integer, float, string] triples

>>> subl = readSection(lst,'[ begin ]','[ end ]') # get all lines between
[ begin ] and [ end ]
exception pmx.parser.ParserError(s)

Bases: Exception

pmx.parser.kickOutComments(lines, comment='#')
pmx.parser.parseList(format_string, lst, ignore_missing=False)
pmx.parser.readSection(lines, begin, end)
pmx.parser.read_and_format(filename, format_string, comment='#', ignore_missing=False)
format_stringstr

e.g. ‘is’ = [int, str]; ‘iis’ = [int, int, str] etc.

pmx.parser.read_fasta(fn)
pmx.parser.read_xvg(fn, style='xy')

pmx.rotamer

This file contains stuff to deal with the Dunbrack rotamer library

pmx.rotamer.check_overlaps(model, mol, nb_list)
pmx.rotamer.get_rotamers(bbdep, resname, phi, psi, residue=False, hydrogens=True, full=False)
pmx.rotamer.load_bbdep()
pmx.rotamer.make_bbdep(min_val=0.01)
pmx.rotamer.mini_nb(model, mol, cutoff)
pmx.rotamer.mutate(residue, new_aa, bbdep)
pmx.rotamer.real_resname(r)
pmx.rotamer.select_best_rotamer(model, rotamers)

pmx.xdrfile

class pmx.xdrfile.Frame(n, mode, x=None, box=None, units=None, v=None, f=None)

Bases: object

get_box(mode='std')
get_natoms()
get_prec()
get_step()
get_time()
update(atom_sel)
update_atoms(atom_sel)
update_box(box)
class pmx.xdrfile.XDRFile(fn, mode='Auto', ft='Auto', atomNum=False)

Bases: object

close()

Close the xdr file. Call this when you are done reading/writing to avoid a memory leak.

exdr3DX = 7
exdrCLOSE = 8
exdrDOUBLE = 3
exdrENDOFFILE = 11
exdrFLOAT = 5
exdrHEADER = 1
exdrINT = 4
exdrMAGIC = 9
exdrNOMEM = 10
exdrNR = 12
exdrOK = 0
exdrSTRING = 2
exdrUINT = 6
write_xtc_frame(step=0, time=0.0, prec=1000.0, lam=0.0, box=False, x=False, units='A', bTrr=False)

Notes

1(1,2)

Jarzynski C (1997) Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach. Phys Rev E 56:5018–5035

2

Goette M, Grubmüller H (2009) Accuracy and convergence of free energy differences calculated from nonequilibrium switching processes. J Comput Chem 30(3):447–456

3

Shirts MR, Bair E, Hooker G, Pande VS (2003) Equilibrium free energies from non-equilibrium measurements using maximum-likelihood methods. Phys Rev Lett 91 (14):140601

4(1,2)

Hahn AM, Then H (2010) Measuring the convergence of Monte Carlo free-energy calculations. Phys Rev E 81(4):041117

5(1,2)

Boresch S, Tettinger F, Leitgeb M, Karplus M (2003) Absolute Binding Free Energies:  A Quantitative Approach for Their Calculation. J Phys Chem B 107(35):9535-9551

6(1,2,3,4)

Hummer, G (2001) Fast-growth thermodynamic integration: Error and efficiency analysis J. Chem. Phys. 114