Common Workflows

As may be expected, most of the workflows used to estimate the physical properties within the framework make use of very similar workflows. This page aims to document the built-in ‘template’ workflows from which the more complex physical property estimation workflows are constructed.

Direct Simulation

Properties being estimated using the direct simulation calculation layer typically base their workflows off of the generate_simulation_protocols() template.

Note

This template currently assumes that a liquid phase property is being computed.

The workflow produced by this template proceeds as follows:

  1. 1000 molecules are inserted into a simulation box with an approximate density of 0.95 g / mL using packmol (BuildCoordinatesPackmol).

  2. the system is parameterized using either the OpenFF toolkit, TLeap or LigParGen depending on the force field being employed (BuildSmirnoffSystem, BuildTLeapSystem or BuildLigParGenSystem).

  3. an energy minimization is performed using the default OpenMM energy minimizer (OpenMMEnergyMinimisation).

  4. the system is equilibrated by running a short NPT simulation for 100000 steps using a timestep of 2 fs and using the OpenMM simulation engine (OpenMMSimulation).

  5. while the uncertainty in the average observable is greater than the requested tolerance (if specified):

    5a) a longer NPT production simulation is run for 1000000 steps with a timestep of 2 fs and using the OpenMM simulation protocol (OpenMMSimulation) with its default Langevin integrator and Monte Carlo barostat.

    5b) the correlated samples are removed from the simulation outputs and the average value of the observable of interest and its uncertainty are computed by bootstrapping with replacement for 250 iterations (AverageObservable). See [1] for details of the decorrelation procedure.

    5c) steps 5a) and 5b) are repeated until the uncertainty condition (if applicable) is met.

The decorrelated simulation outputs are then made available ready to be cached by a storage backend (DecorrelateObservables, DecorrelateTrajectory).

MBAR Reweighting

Properties being estimated using the MBAR reweighting calculation layer typically base their workflows off of the generate_reweighting_protocols() template.

The workflow produced by this template proceeds as follows:

  1. for each stored simulation data:

    1a) the cached data is retrieved from disk (UnpackStoredSimulationData)

  2. the cached data from is concatenated together to form a single trajectory of configurations and observables (ConcatenateTrajectories, ConcatenateStatistics).

  3. for each stored simulation data:

    3a) the system is parameterized using the force field parameters which were used when originally generating the cached data i.e. one of the reference states (BuildSmirnoffSystem, BuildTLeapSystem or BuildLigParGenSystem).

    3b) the reduced potential of each configuration in the concatenated trajectory is evaluated using the parameterized system (OpenMMEvaluateEnergies).

  4. the system is parameterized using the force field parameters with which the property of interest should be calculated using i.e. of the target state (BuildSmirnoffSystem, BuildTLeapSystem or BuildLigParGenSystem) and the reduced potential of each configuration in the concatenated trajectory is evaluated using the parameterized system (OpenMMEvaluateEnergies).

    4a) (optional) if the observable of interest is a function of the force field parameters it is recomputed using the target state parameters. These recomputed values then replace the original concatenated observables loaded from the cached data.

  5. the reference potentials, target potentials and the joined observables are sub-sampled to only retain equilibrated, uncorrelated samples (AverageObservable, DecorrelateObservables, DecorrelateTrajectory). See [1] for details of the decorrelation procedure.

  6. the MBAR method is employed to compute the average value of the observable of interest and its uncertainty at the target state, taking the reference state reduced potentials as input. See [2] for the theory behind this approach. An exception is raised if there are not enough effective samples to reweight (ReweightObservable).

In more specialised cases the generate_base_reweighting_protocols() template (which generate_reweighting_protocols() is built off of) is instead used due to its greater flexibility.

References

1(1,2)

John D Chodera. A simple method for automated equilibration detection in molecular simulations. Journal of chemical theory and computation, 12(4):1799–1805, 2016.

2

Richard A Messerly, S Mostafa Razavi, and Michael R Shirts. Configuration-sampling-based surrogate models for rapid parameterization of non-bonded interactions. Journal of Chemical Theory and Computation, 14(6):3144–3162, 2018.