openmmslicer.smc module

This is the main module in openmmslicer, where one can find the SequentialSampler, which makes use of all other modules.

class openmmslicer.smc.SequentialSampler(coordinates, structure, integrator, moves, platform=None, platform_properties=None, npt=True, checkpoint=None, md_config=None, alch_config=None)

Bases: object

A complete sequential Monte Carlo sampler which can enhance the sampling of certain degrees of freedom.

Parameters

• coordinates (str) – Path to a file which contains all coordinates.

• structure (parmed.Structure) – An object containing all structural information.

• integrator (openmm.Integrator) – The integrator which should be used during the sampling step.

• moves (openmmslicer.moves.MoveList or [openmmslicer.moves.Move]) – All moves which must be applied at lambda = 0.

• platform (str) – The platform which should be used for simulation.

• platform_properties (dict) – Additional platform properties.

• npt (bool) – Whether to add a barostat at 1 atm.

• checkpoint (str) – A path to a pickled checkpoint file to load SequentialSampler from. If this is None, SequentialSampler is initialised normally.

• md_config (dict) – Additional parameters passed to generateSystem().

• alch_config (dict) – Additional parameters passed to generateAlchSystem().

coordinates

Path to a file which contains all coordinates.

Type

str

moves

All moves which must be applied at lambda = 0.

Type

[Move]

structure

An object containing all structural information.

Type

parmed.Structure

system

The initial unmodified system.

Type

openmm.System

alch_system

The modified alchemical system.

Type

openmm.System

platform

The currently used platform.

Type

str

platform_properties

The currently used platform properties.

Type

dict

integrator

The integrator used during the sampling step.

Type

openmm.integrator

simulation

A simulation object of the current alchemical system.

Type

openmm.Simulation

lambda_

The current lambda value.

Type

float

total_sampling_steps

A counter keeping track of the number of times the integrator was called.

Type

int

lambda_history

A list containing all past lambda values.

Type

list

deltaE_history

A list containing all past deltaE values.

Type

list

weight_history

A list containing all past weights used for resampling.

Type

list

reporter_history

A list containing paths to all past trajectory files.

Type

list

current_states

A list containing all current states.

Type

[int] or [openmm.State]

logZ

The current estimate of the dimensionless free energy difference.

Type

float

classmethod addBarostat(system, temperature=Quantity(value=298, unit=kelvin), pressure=Quantity(value=1, unit=atmosphere), frequency=25, **kwargs)

Adds a MonteCarloBarostat to the MD system.

Parameters

• system (openmm.System) – The OpenMM System object corresponding to the reference system.

• temperature (openmm.unit.Quantity) – Temperature (Kelvin) to be simulated at.

• pressure (openmm.unit.Quantity) – Pressure (atm) for Barostat for NPT simulations.

• frequency (int, default=25) – Frequency at which Monte Carlo pressure changes should be attempted (in time steps)

Returns

system – The OpenMM System with the MonteCarloBarostat attached.

Return type

openmm.System

property alchemical_atoms

The absolute indices of all alchemical atoms.

Type

[int]

calculateStateEnergies(lambda_, states=None, extra_conformers=None, *args, **kwargs)

Calculates the reduced potential energies of all states for a given lambda value.

Parameters

• lambda (float) – The desired lambda value.

• states ([int] or [openmm.State] or None) – Which states need to be used. If None, self.current_states are used. Otherwise, these could be in any format supported by setState().

• extra_conformers (list) – Extra conformers relevant for lambda = 0. These could be either a list of openmm.State or a list of transformations which can be generated dynamically.

• args – Positional arguments to be passed to setState().

• kwargs – Keyword arguments to be passed to setState().

equilibrate(equilibration_steps=100000, restrain_backbone=True, force_constant=Quantity(value=5.0, unit=kilocalorie / angstrom ** 2 * mole), output_interval=100000, reporter_filename=None)

Equilibrates the system at lambda = 0.

Parameters

• equilibration_steps (int) – The number of equilibration steps.

• restrain_backbone (bool) – Whether to restrain all atoms with the following names: ‘CA’, ‘C’, ‘N’.

• force_constant (openmm.unit.Quantity) – The magnitude of the restraint force constant.

• interval (output) – How often to output to a trajectory file.

• reporter_filename (str) – The path to the trajectory file.

static generateAlchSystem(system, atom_indices, softcore_alpha=0.5, softcore_a=1, softcore_b=1, softcore_c=6, softcore_beta=0.5, softcore_d=1, softcore_e=1, softcore_f=2, alchemical_torsions=True, annihilate_electrostatics=True, annihilate_sterics=True, disable_alchemical_dispersion_correction=True, alchemical_pme_treatment='direct-space', suppress_warnings=True, **kwargs)

Returns the OpenMM System for alchemical perturbations. This function calls openmmtools.alchemy.AbsoluteAlchemicalFactory and openmmtools.alchemy.AlchemicalRegion to generate the System for the SMC simulation.

Parameters

• system (openmm.System) – The OpenMM System object corresponding to the reference system.

• atom_indices (list of int) – Atom indicies of the move or designated for which the nonbonded forces (both sterics and electrostatics components) have to be alchemically modified.

• annihilate_electrostatics (bool, optional) – If True, electrostatics should be annihilated, rather than decoupled (default is True).

• annihilate_sterics (bool, optional) – If True, sterics (Lennard-Jones or Halgren potential) will be annihilated, rather than decoupled (default is False).

• softcore_alpha (float, optional) – Alchemical softcore parameter for Lennard-Jones (default is 0.5).

• softcore_a (float, optional) – Parameters modifying softcore Lennard-Jones form. Introduced in Eq. 13 of Ref. [TTPham-JChemPhys135-2011] (default is 1).

• softcore_b (float, optional) – Parameters modifying softcore Lennard-Jones form. Introduced in Eq. 13 of Ref. [TTPham-JChemPhys135-2011] (default is 1).

• softcore_c (float, optional) – Parameters modifying softcore Lennard-Jones form. Introduced in Eq. 13 of Ref. [TTPham-JChemPhys135-2011] (default is 1).

• softcore_beta (float, optional) – Alchemical softcore parameter for electrostatics. Set this to zero to recover standard electrostatic scaling (default is 0.0).

• softcore_d (float, optional) – Parameters modifying softcore electrostatics form (default is 1).

• softcore_e (float, optional) – Parameters modifying softcore electrostatics form (default is 1).

• softcore_f (float, optional) – Parameters modifying softcore electrostatics form (default is 1).

• disable_alchemical_dispersion_correction (bool, optional, default=True) – If True, the long-range dispersion correction will not be included for the alchemical region to avoid the need to recompute the correction (a CPU operation that takes ~ 0.5 s) every time ‘lambda_sterics’ is changed. If using nonequilibrium protocols, it is recommended that this be set to True since this can lead to enormous (100x) slowdowns if the correction must be recomputed every time step.

• alchemical_pme_treatment (str, optional, default = 'direct-space') – Controls how alchemical region electrostatics are treated when PME is used. Options are ‘direct-space’, ‘coulomb’, ‘exact’. - ‘direct-space’ only models the direct space contribution - ‘coulomb’ includes switched Coulomb interaction - ‘exact’ includes also the reciprocal space contribution, but it’s only possible to annihilate the charges and the softcore parameters controlling the electrostatics are deactivated. Also, with this method, modifying the global variable lambda_electrostatics is not sufficient to control the charges. The recommended way to change them is through the AlchemicalState class.

Returns

alch_system – System to be used for the SMC simulation.

Return type

alchemical_system

References

TTPham-JChemPhys135-2011(1,2,3)

20. 20. Pham and M. R. Shirts, J. Chem. Phys 135, 034114 (2011). http://dx.doi.org/10.1063/1.3607597

generateAlchemicalRegion()

Makes sure that all rotated dihedrals are also made alchemical.

classmethod generateSimFromStruct(structure, system, integrator, platform=None, properties=None, **kwargs)

Generate the OpenMM Simulation objects from a given parmed.Structure()

Parameters

• structure (parmed.Structure) – ParmEd Structure object of the entire system to be simulated.

• system (openmm.System) – The OpenMM System object corresponding to the reference system.

• integrator (openmm.Integrator) – The OpenMM Integrator object for the simulation.

• platform (str, default = None) – Valid choices: ‘Auto’, ‘OpenCL’, ‘CUDA’ If None is specified, the fastest available platform will be used.

• properties (dict) – Additional platform properties.

Returns

simulation – The generated OpenMM Simulation from the parmed.Structure, openmm.System, amd the integrator.

Return type

openmm.Simulation

static generateSystem(structure, **kwargs)

Construct an OpenMM System representing the topology described by the prmtop file. This function is just a wrapper for parmed Structure.createSystem().

Parameters

• structure (parmed.Structure()) – The parmed.Structure of the molecular system to be simulated

• nonbondedMethod (cutoff method) – This is the cutoff method. It can be either the NoCutoff, CutoffNonPeriodic, CutoffPeriodic, PME, or Ewald objects from the simtk.openmm.app namespace

• nonbondedCutoff (float or distance Quantity) – The nonbonded cutoff must be either a floating point number (interpreted as nanometers) or a Quantity with attached units. This is ignored if nonbondedMethod is NoCutoff.

• switchDistance (float or distance Quantity) – The distance at which the switching function is turned on for van der Waals interactions. This is ignored when no cutoff is used, and no switch is used if switchDistance is 0, negative, or greater than the cutoff

• constraints (None, app.HBonds, app.HAngles, or app.AllBonds) – Which type of constraints to add to the system (e.g., SHAKE). None means no bonds are constrained. HBonds means bonds with hydrogen are constrained

• rigidWater (bool=True) – If True, water is kept rigid regardless of the value of constraints. A value of False is ignored if constraints is not None.

• implicitSolvent (None, app.HCT, app.OBC1, app.OBC2, app.GBn, app.GBn2) – The Generalized Born implicit solvent model to use.

• implicitSolventKappa (float or 1/distance Quantity = None) – This is the Debye kappa property related to modeling saltwater conditions in GB. It should have units of 1/distance (1/nanometers is assumed if no units present). A value of None means that kappa will be calculated from implicitSolventSaltConc (below)

• implicitSolventSaltConc (float or amount/volume Quantity=0 moles/liter) – If implicitSolventKappa is None, the kappa will be computed from the salt concentration. It should have units compatible with mol/L

• temperature (float or temperature Quantity = 298.15 kelvin) – This is only used to compute kappa from implicitSolventSaltConc

• soluteDielectric (float=1.0) – The dielectric constant of the protein interior used in GB

• solventDielectric (float=78.5) – The dielectric constant of the water used in GB

• useSASA (bool=False) – If True, use the ACE non-polar solvation model. Otherwise, use no SASA-based nonpolar solvation model.

• removeCMMotion (bool=True) – If True, the center-of-mass motion will be removed periodically during the simulation. If False, it will not.

• hydrogenMass (float or mass quantity = None) – If not None, hydrogen masses will be changed to this mass and the difference subtracted from the attached heavy atom (hydrogen mass repartitioning)

• ewaldErrorTolerance (float=0.0005) – When using PME or Ewald, the Ewald parameters will be calculated from this value

• flexibleConstraints (bool=True) – If False, the energies and forces from the constrained degrees of freedom will NOT be computed. If True, they will (but those degrees of freedom will still be constrained).

• verbose (bool=False) – If True, the progress of this subroutine will be printed to stdout

• splitDihedrals (bool=False) – If True, the dihedrals will be split into two forces – proper and impropers. This is primarily useful for debugging torsion parameter assignments.

Returns

system – System formatted according to the prmtop file.

Return type

openmm.System

Notes

This function calls prune_empty_terms if any Topology lists have changed.

property kT

The current temperature multiplied by the gas constant.

Type

openmm.unit.Quantity

property moves

The list of moves which must be performed at lambda = 0.

property n_walkers

The number of current walkers.

Type

int

run(final_decorrelation_step=True, *args, **kwargs)

Performs a complete sequential Monte Carlo run until lambda = 1.

Parameters

• final_decorrelation_step (bool) – Whether to decorrelate the final resampled walkers for another number of default_decorrelation_steps.

• args – Positional arguments to be passed to runSingleIteration().

• kwargs – Keyword arguments to be passed to runSingleIteration().

runSingleIteration(sampling_metric=<class 'openmmslicer.sampling_metrics.EnergyCorrelation'>, resampling_method=<class 'openmmslicer.resampling_methods.SystematicResampler'>, resampling_metric=<class 'openmmslicer.resampling_metrics.WorstCaseSampleSize'>, target_metric_value=None, target_metric_value_initial=None, target_metric_tol=None, maximum_metric_evaluations=20, default_dlambda=0.1, minimum_dlambda=0.01, default_decorrelation_steps=500, maximum_decorrelation_steps=5000, reporter_filename=None, n_walkers=1000, n_conformers_per_walker=100, equilibration_options=None, dynamically_generate_conformers=True, keep_walkers_in_memory=False, write_checkpoint=None, write_instant_checkpoint=None, load_instant_checkpoint=None)

Performs a single iteration of sequential Monte Carlo.

Parameters

• sampling_metric (class) – A sampling metric with callable methods as described in openmmslicer.sampling_metrics. This metric is used to adaptively determine the optimal sampling time. None removes adaptive sampling.

• resampling_method (class) – A resampling method with callable methods as described in openmmslicer.resampling_methods.

• resampling_metric (class) – A resampling metric with callable methods as described in openmmslicer.resampling_metrics. This metric is used to adaptively determine the optimal next lambda. None removes adaptive resampling.

• target_metric_value (float) – The threshold for the resampling metric. None uses the default value given by the class.

• target_metric_value_initial (float) – Same as target_metric_value, but only valid for lambda = 0.

• target_metric_tol (float) – The relative tolerance for the resampling metric. None uses the default value given by the class.

• maximum_metric_evaluations (int) – The maximum number of energy evaluations to determine the resampling metric.

• default_dlambda (float) – Determines the next lambda value. Only used if resampling_metric is None.

• minimum_dlambda (float) – The minimum allowed change in lambda.

• default_decorrelation_steps (int) – The default number of decorrelation steps. If sampling metric is None, these denote the true number of decorrelation steps.

• maximum_decorrelation_steps (int) – The maximum number of decorrelation steps. Only used with adaptive sampling.

• reporter_filename (str) – A template containing the reporter filename. Must contain curly brackets, so that python’s format() can be called.

• n_walkers (int) – The number of walkers to be resampled.

• n_conformers_per_walker (int) – How many conformers to generate for each walker using self.moves. Only used at lambda = 0.

• equilibration_options (dict) – Parameters to be passed to equilibrate(). Only used at lambda = 0.

• dynamically_generate_conformers (bool) – Whether to store the extra conformers as states or as transformations. The former is much faster, but also extremely memory-intensive. Only set to False if you are certain that you have enough memory.

• keep_walkers_in_memory (bool) – Whether to keep the walkers as states or load them dynamically from the hard drive. The former is much faster during the energy evaluation step but is more memory-intensive. Only set to True if you are certain that you have enough memory.

• write_checkpoint (str) – Describes a path to a checkpoint file written after completing the whole SMC iteration. None means no such file will be written.

• write_instant_checkpoint (str) – Describes a path to a checkpoint file written immediately after the decorrelation steps for each walker. None means no such file will be written.

• load_instant_checkpoint (str) – Describes a path to a checkpoint file written with write_instant_checkpoint. In order to use this, the same arguments must be passed to this function as when the instant checkpoint file was written. This option also typically requires that the SequentialSampler was instantiated with a particular checkpoint. None means that this function proceeds as normal.

saveCheckpoint(filename='checkpoint.pickle', *args, **kwargs)

Saves a pickled checkpoint file.

Parameters

• filename (str) – Path to the pickle filename to save the checkpoint to.

• args – Positional arguments to be passed to pickle.dump().

• kwargs – Keyword arguments to be passed to pickle.dump().

setState(state, context, reporter_filename=None, transform=None)

Sets a given state to the current context.

Parameters

• state (int, openmm.State) – Either sets the state from a State object or from a frame of a trajectory file given by reporter_filename.

• context (openmm.Context) – The context to which the state needs to be applied.

• reporter_filename (str) – The path to the trajectory file containing the relevant frame, if applicable.

• transform – Optionally generate a transform dynamically from a format, specific to the underlying moves.

property temperature

The temperature of the current integrator.

Type

openmm.unit.Quantity