Molly API

Main module for Molly.jl - molecular simulation in Julia.

See https://juliamolsim.github.io/Molly.jl/stable for documentation.

ANIPotential(path; ensemble_idx=nothing)

The ANI potential implemented as an AtomsCalculators.jl calculator.

See Smith et al. 2017. Load an ANI-2x neural network potential from an HDF5 file exported by test/torchani_reference.py, such as the one provided by Molly.ani2x_data_dir. Requires Lux and HDF5 to be loaded.

The system's atoms_data is required, since the element of each atom is read from atoms_data[i].element. The supported elements are H, C, N, O, S, F and Cl.

By default all ensemble members are loaded and their energies averaged. Pass ensemble_idx (one-indexed, 1:8 for ANI-2x) to load only a single member.

Coordinates without units are treated as nm following the Molly convention, and converted internally to the Å the ANI parameters use. Periodic systems must use a neighbor finder with the eligible pairs set.

The ANI-2x weights are Float32, so the energy/force paths run in Float32 internally regardless of the system's coordinate type. Not compatible with virial calculation.

ASECalculator(; <keyword arguments>)

A Python ASE calculator.

This calculator is only available when PythonCall is imported. It is the user's responsibility to have the required Python packages installed. This includes ASE and any packages providing the calculator.

Contrary to the rest of Molly, unitless quantities are assumed to have ASE units: Å for length, eV for energy, u for mass, and Å sqrt(u/eV) for time. Unitful quantities will be converted as appropriate.

The virial is calculated only if the ASE calculator supports stress calculation. Otherwise the virial contribution is ignored with a warning.

Arguments

• ase_calc: the ASE calculator created with PythonCall.
• atoms: the atoms, or atom equivalents, in the system.
• coords: the coordinates of the atoms in the system. Typically a vector of SVectors of 2 or 3 dimensions.
• boundary: the bounding box in which the simulation takes place.
• elements=nothing: vector of atom elements as a string, either elements or atoms_data (which contains element data) must be provided.
• atoms_data=nothing: other data associated with the atoms.
• velocities=nothing: the velocities of the atoms in the system, only required if the velocities contribute to the potential energy or forces.

AWHSimulation(awh_state::AWHState{T};
              <keyword arguments>)

Prepare an Accelerated Weight Histogram (AWH) simulation, which is run with simulate!.

This struct stores the parameters controlling the AWH updates. It can optionally update a PMFDeconvolution object on the fly to obtain an unbiased PMF, as described in Lindahl et al. 2014. Deconvolution is recommended for simulations where the λ coordinate represents a biased reaction coordinate (e.g. an umbrella potential). It is typically not required for standard alchemical transformations.

Arguments

• awh_state::AWHState{T}: The AWHState containing the λ windows and accumulators.
• n_md_steps::Int=10: Number of MD integration steps to perform between AWH coordinate sampling and reweighting.
• update_freq::Int=1: Number of samples to collect before applying an update to the AWH bias.
• well_tempered_factor::Real=10.0: If finite, the AWH target distribution (ρ) is dynamically scaled to favor low-energy λ windows. Smaller values accentuate this behavior. Use Inf to disable.
• coverage_threshold::Real=1.0: Proportion of λ windows that must be visited to double the fictitious sample size and advance the initial stage.
• significant_weight::Real=0.1: The fractional weight threshold (relative to an ideally uniform distribution) required for a λ window to be considered "visited". This filters out sampling noise.
• log_freq::Int=100: Number of AWH iterations between logging statistics.
• loggers=(): Loggers attached to the active system used for AWH dynamics.
• pmf=nothing: Optional PMFDeconvolution object created from awh_state.
• initial_step::Int=0: Absolute MD step for a new or resumed simulation.

AWHState(thermo_states::AbstractArray{<:ThermoState};
         <keyword arguments>)

State of an Accelerated Weight Histogram (AWH) simulation, see Lindahl et al. 2014.

Maintains the physical state of the system across multiple λ windows, as well as the accumulated statistical weights, free energy estimates, and target distribution parameters.

Arguments

• thermo_states::AbstractArray{<:ThermoState}: Iterable containing the ThermoState structs defining each λ window of the AWH simulation.
• first_state::Int=1: The index of the initial active λ state.
• n_bias::Int=100: Fictitious effective sampling size used during the initial stage. Smaller values yield more aggressive free energy updates.
• ρ::Union{Nothing, AbstractArray{T}}=nothing: Target distribution array along λ. If nothing, a uniform distribution is used.
• reuse_neighbors::Bool=true: Whether to reuse the active system's neighbor list when calculating energies for adjacent λ windows. Generally improves performance.

Loggers are attached by AWHSimulation, not by AWHState. History-dependent general interactions such as Metadynamics are not supported.

Base type for boundary types.

Parameters are dimension, float type and coordinate type. The built-in boundary types are CubicBoundary, TriclinicBoundary and RectangularBoundary.

Base type for non-bonded interaction cutoffs.

The parameter P is the number of cutoff points: 0 for NoCutoff, 1 for cutoffs with a dist_cutoff field such as DistanceCutoff, and 2 for cutoffs that also have a dist_activation field such as CubicSplineCutoff. Custom cutoffs should subtype this.

AbstractMetadynamicsMemory

Abstract type for the different ways a Metadynamics can store deposited Gaussian hills.

Subtypes must have the fields cvs, the tuple of CVs, and height, the default hill height, and implement potential_energy(memory, cv_value; kwargs...), Molly.bias_gradient(memory, cv_value; kwargs...) and Molly.add_hill!(memory, cv_value, height), where cv_value is a scalar for one CV and a tuple for several. The keyword arguments include the strictness of the simulation. Subtypes that support periodic CVs take the periods from Molly.cv_period in their constructor. Optionally, they implement Molly.can_deposit(memory, cv_value; kwargs...) to skip hills (e.g. outside a grid), and bias_points_cache with bias_at_points for a cheaper offset of a MetadynamicsLogger. Built-in subtypes are ListHills and GridHills.

AbstractTempering

Controls how the height of each newly deposited hill is scaled before being added to a Metadynamics's memory.

Subtypes must implement Molly.tempering_height, the height to actually deposit, and Molly.free_energy_scale, the prefactor of the free energy estimate. A MetadynamicsLogger takes its temperature from WellTempered, or from its temperature keyword for other temperings.

Built-in subtypes are Molly.NoTempering (the default) and WellTempered.

ActiveThermoState(space, first_state=1; loggers=())

Mutable active state for simulations over the thermodynamic states of a state space, such as the state_space field of a TSSState or an AWHState.

The active state owns a simulation system and integrator corresponding to the current thermodynamic state index, starting at first_state, and is moved between states without reallocating the active wrapper. loggers are attached to the active system. For example, explicit active states can be given as replica_active_states to TSSSimulation.

AndersenThermostat(temperature, coupling_const) <: AbstractThermostat

The Andersen thermostat for controlling temperature.

The velocity of each atom is randomly changed each time step with probability 1 - exp(-dt / coupling_const), which is approximately dt / coupling_const, to a velocity drawn from the Maxwell-Boltzmann distribution. See Andersen 1980.

Note that when used with constraints this can lead to velocities that violate the constraints.

AngleConstraint(i, j, k, angle_ijk, dist_ij, dist_jk)

Constraint between three atoms that maintains a fixed angle and two bond lengths.

Given to a constraint algorithm such as SHAKE_RATTLE or LINCS. Atoms i and k should be connected to central atom j with fixed bond distances given by dist_ij and dist_jk, forming the angle angle_ijk in radians. Internally, an AngleConstraint is converted into 3 distance constraints. None of the atoms in this constraint should be constrained with atoms not part of this constraint.

For example, a water molecule can be defined as AngleConstraint(1, 2, 3, deg2rad(104.5), 0.9572u"Å", 0.9572u"Å") where atom 2 is oxygen and atoms 1/3 are hydrogen.

Linear molecules like CO2 cannot be constrained.

AshbaughHatch(; cutoff=NoCutoff(), use_neighbors=true, shortcut=Molly.LJZeroShortcut(),
              σ_mixing=Molly.LorentzMixing(), ϵ_mixing=Molly.LorentzMixing(),
              λ_mixing=Molly.LorentzMixing(), weight_special=1)

The Ashbaugh-Hatch potential ($V_{\text{AH}}$) is a modified Lennard-Jones ($V_{\text{LJ}}$) 6-12 interaction between two atoms.

See Ashbaugh and Hatch 2008. The potential energy is defined as

\[V_{\text{LJ}}(r_{ij}) = 4\varepsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]\]

\[V_{\text{AH}}(r_{ij}) = \begin{cases} V_{\text{LJ}}(r_{ij}) + \varepsilon_{ij}(1-\lambda_{ij}), & r_{ij} \leq 2^{1/6}\sigma_{ij} \\ \lambda_{ij} V_{\text{LJ}}(r_{ij}), & r_{ij} > 2^{1/6}\sigma_{ij} \end{cases}\]

and the force on each atom by

\[\vec{F}_{\text{AH}} = \begin{cases} \vec{F}_{\text{LJ}}(r_{ij}), & r_{ij} \leq 2^{1/6}\sigma_{ij} \\ \lambda_{ij} \vec{F}_{\text{LJ}}(r_{ij}), & r_{ij} > 2^{1/6}\sigma_{ij} \end{cases}\]

where

\[\vec{F}_{\text{LJ}}(r_{ij}) = \frac{24\varepsilon_{ij}}{r_{ij}^2} \left[2\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} -\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right] \vec{r}_{ij}\]

$\lambda_{ij}$ is the λ of the atoms mixed with λ_mixing, which here is a per-atom hydropathy parameter rather than an alchemical coupling. If $\lambda_{ij}$ is 1 this gives the standard LennardJones potential.

Atom(; <keyword arguments>)

An atom and its associated information.

Properties unused in the simulation or in analysis can be left with their default values. The types used should be bits types if the GPU is going to be used.

Arguments

• index::Int32=1: the index of the atom in the system. This only needs to be set if it is used in the interactions. The order of atoms is determined by their order in the atom vector.
• atom_type::T=1: the type of the atom. This only needs to be set if it is used in the interactions.
• mass::M=1.0u"g/mol": the mass of the atom.
• charge::C=0.0: the charge of the atom, used for electrostatic interactions.
• σ::S=0.0u"nm": the Lennard-Jones finite distance at which the inter-particle potential is zero.
• ϵ::E=0.0u"kJ * mol^-1": the Lennard-Jones depth of the potential well.
• λ::L=1.0: the λ of the atom. In an alchemical system the alchemical atoms carry the global λ, which the λ interactions combine with the alchemical roles to scale the interactions, and other atoms keep 1.
• alch_role::Int32=Molly.EnvRole: role of the atom in an alchemical transformation. EnvRole, the default, is an atom that is not alchemical. The other roles are set by AbsoluteFESystem and RelativeFESystem, see Alchemical roles.

AtomData(; atom_type="?", atom_name="?", res_number=1, res_name="???",
         chain_id="A", element="?", hetero_atom=false, charge_type="")

Data associated with an atom.

charge_type labels the atoms that share a charge parameter in extract_parameters and inject_gradients, see the charge_types argument of System. If empty, the atom type is used.

Storing this separately allows the Atom types to be bits types and hence work on the GPU.

AverageObservableLogger(observable::Function, T::DataType, n_steps::Integer;
                        n_blocks::Integer=1024)

A logger that periodically records observations of a system and keeps a running empirical average.

While GeneralObservableLogger holds a full record of observations, AverageObservableLogger does not. In addition, calling values(logger::AverageObservableLogger; std::Bool=true) returns two values: the current running average, and an estimate of the standard deviation for this average based on the block averaging method described in Flyvbjerg and Petersen 1989.

Arguments

• observable::Function: the observable whose mean is recorded, must support the method observable(s::System, neighbors, step_n, buffers; n_threads::Integer).
• T::DataType: the type returned by observable.
• n_steps::Integer: number of simulation steps between observations.
• n_blocks::Integer=1024: the number of blocks used in the block averaging method, should be an even number.

BerendsenBarostat(pressure, coupling_const;
                  coupling_type=:isotropic,
                  compressibility=4.6e-5u"bar^-1",
                  max_scale_frac=0.1, n_steps=30) <: AbstractBarostat

The Berendsen barostat for controlling pressure.

See Berendsen et al. 1984. The scaling factor for the box every n_steps steps is

\[\mu_{ij} = \delta_{ij} - \frac{\Delta t}{3\, \tau_p} \kappa_{ij} \left(P_{0ij} - P_{ij}(t) \right).\]

with the fractional change limited to max_scale_frac.

The scaling factor $\mu$ is a matrix and $\delta$ is a Kronecker delta, allowing non-isotropic pressure control. Available options for coupling_type are :isotropic, :semiisotropic and :anisotropic.

This barostat should be used with caution as it is known not to properly sample isobaric ensembles and therefore can lead to simulation artifacts.

BerendsenThermostat(temperature, coupling_const) <: AbstractThermostat

The Berendsen thermostat for controlling temperature.

See Berendsen et al. 1984. The scaling factor for the velocities each step is

\[\lambda^2 = 1 + \frac{\delta t}{\tau} \left( \frac{T_0}{T} - 1 \right)\]

This thermostat should be used with caution as it can lead to simulation artifacts.

BiasPotential(cv_type, bias_type)

A potential to bias a simulation along a collective variable (CV), implemented as an AtomsCalculators.jl calculator.

The cv_type could for example be CalcDist and the bias_type could be LinearBias.

Forces resulting from the bias potential are evaluated in two steps, specifically by (1) calculating the gradient of the bias potential with respect to the value of the CV using Molly.bias_gradient, and (2) calculating the gradient of the CV with respect to the atomic coordinates using Molly.cv_gradient.

Gradients can be calculated with either automatic differentiation or explicitly defined gradient functions. Enzyme should be imported in the first case.

Molecules are made whole before calculating the CV if the correction of cv_type is :pbc. The virial contribution is only calculated for CV types that define it, which includes the built-in CV types when has_virial is true.

Buckingham(; cutoff=NoCutoff(), use_neighbors=true,
           shortcut=Molly.BuckinghamZeroShortcut(), A_mixing=Molly.GeometricMixing(),
           B_mixing=Molly.InverseMixing(), C_mixing=Molly.GeometricMixing(),
           weight_special=1)

The Buckingham interaction between two atoms.

The potential energy is defined as

\[V(r_{ij}) = A_{ij} \exp(-B_{ij} r_{ij}) - \frac{C_{ij}}{r_{ij}^6}\]

and the force on each atom by

\[\vec{F}_i = \left( A_{ij} B_{ij} \exp(-B_{ij} r_{ij}) - 6 \frac{C_{ij}}{r_{ij}^7} \right) \frac{\vec{r}_{ij}}{r_{ij}}\]

By default, using Molly.GeometricMixing for A and C and Molly.InverseMixing for B, the parameters are derived from the atom parameters according to

\[\begin{aligned} A_{ij} &= (A_{ii} A_{jj})^{1/2} \\ B_{ij} &= \frac{2}{\frac{1}{B_{ii}} + \frac{1}{B_{jj}}} \\ C_{ij} &= (C_{ii} C_{jj})^{1/2} \end{aligned}\]

so atoms that use this interaction should have fields A, B and C available. The default shortcut is Molly.BuckinghamZeroShortcut.

BuckinghamZeroShortcut()

Skip a pairwise interaction when the A or C of either atom is zero.

Give this when setting up a pairwise interaction such as Buckingham.

CMAPTorsion(index, size)

Torsional correction map (CMAP) for sets of five atoms, for example protein ϕ and ψ backbone torsion angles.

See Mackerell et al. 2004. The CMAP data is stored in the data field of the associated InteractionList5Atoms.

Only compatible with 3D systems.

CMAPTorsionλ(; index, size, λ_mixing=Molly.MinimumMixing(),
             scheduler=Molly.DefaultLambdaScheduler())

Torsional correction map (CMAP) for sets of five atoms in an alchemical system, the λ version of CMAPTorsion.

AbsoluteFESystem and RelativeFESystem convert each CMAPTorsion to this type. The potential energy is defined as

\[V(\phi, \psi) = w V_{\mathrm{CMAP}}(\phi, \psi)\]

where $V_{\mathrm{CMAP}}$ is the potential energy of the CMAPTorsion with the same index and size, and the CMAP data is stored as for CMAPTorsion. The weight $w$ depends on the alchemical roles of the atoms and on $\lambda$, the λ of the atoms mixed with λ_mixing. $w = 1 - \lambda$ if the interaction involves a core atom with the parameters of state A, $w = \lambda$ if it involves a core atom with the parameters of state B and $w = 1$ if it involves only unique and environment atoms.

scheduler is a Molly.LambdaScheduler, which has to use dual topology. Single topology is not supported for this interaction, so RelativeFESystem throws an error if the scheduler has dual=false.

Only compatible with 3D systems.

CRescaleBarostat(pressure, coupling_const;
                 coupling_type=:isotropic,
                 compressibility=4.6e-5u"bar^-1",
                 max_scale_frac=0.1, n_steps=30) <: AbstractBarostat

The stochastic cell rescale barostat.

See Bernetti and Bussi 2020 and Del Tatto et al. 2022. In brief, this is an extension of the BerendsenBarostat that includes a stochastic term to the scaling matrix. This allows proper sampling of isobaric ensembles.

\[\mu = \rm{exp}\left[ \frac{-\kappa_T \cdot \Delta t}{\tau_P} \cdot (P(t) - P_0) + \sqrt{\frac{2 \cdot k_BT \cdot \kappa_T \cdot \Delta t}{V(t) \cdot \tau_P}} \cdot dW \right]\]

where $\kappa_T$ is the isothermal compressibility, $\tau_P$ is the barostat coupling constant and $\rm{dW}$ represents a Wiener process. The fractional change is limited to max_scale_frac.

The scaling factor $\mu$ is a matrix, allowing non-isotropic pressure control. Available options for coupling_type are :isotropic, :semiisotropic and :anisotropic. Generally, n_steps multiplied by the simulator dt should fit at least 25 times inside $\tau_P$.

CalcAngle(atom_inds::AbstractVector{Int}=[], correction=:pbc, has_virial::Bool=true)

A collective variable that calculates the angle defined by three atoms.

Given as an argument to BiasPotential. The angle is defined by the angle between the vectors (j -> i) and (j -> k), where the indices are given by atom_inds.

Fields

• atom_inds::AbstractVector{Int}: The indices of the three atoms (i, j, k) defining the angle.
• correction::Symbol: The method used to handle periodic boundary conditions. Defaults to :pbc.
• has_virial::Bool: Whether the virial contribution should be calculated for this collective variable. Defaults to true.

CalcCMDist(calc_type=:closest)

Bias the distance between the centers of mass of two groups of atoms.

Given as an argument to CalcDist. By default, distances are calculated between the closest periodic images. Setting calc_type=:raw means that distances are calculated ignoring PBCs.

Should generally be used with molecule unwrapping since it assumes that the atoms within each group are in the same periodic box. If distances are evaluated using the minimum image convention on an unwrapped system, raw coordinates must be within a distance of 1.5x the box length of each other to ensure correct results.

CalcDist(atom_inds_1, atom_inds_2, dist_type=CalcMinDist(), correction=:pbc,
         has_virial=true)

Bias the distance between two atoms or groups of atoms.

Given as an argument to BiasPotential.

Arguments

• atom_inds_1: indices of the atom(s) in the first group.
• atom_inds_2: indices of the atom(s) in the second group.
• dist_type=CalcMinDist(): type of distance to calculate, one of CalcMinDist, CalcMaxDist, CalcCMDist or CalcSingleDist.
• correction=:pbc: the correction to be applied to the molecules. :pbc keeps molecules whole, :wrap wraps all atoms inside the simulation box. If using multiple atoms in a group, they should generally be in the same molecule and :pbc should be used.
• has_virial=true: whether the virial contribution should be calculated for this collective variable.

CalcMaxDist(calc_type=:closest)

Bias the maximum distance between two groups of atoms.

Given as an argument to CalcDist. By default, distances are calculated between the closest periodic images. Setting calc_type=:raw means that distances are calculated ignoring PBCs.

If distances are evaluated using the minimum image convention on an unwrapped system, raw coordinates must be within a distance of 1.5x the box length of each other to ensure correct results.

CalcMinDist(calc_type=:closest)

Bias the minimum distance between two groups of atoms.

Given as an argument to CalcDist. By default, distances are calculated between the closest periodic images. Setting calc_type=:raw means that distances are calculated ignoring PBCs.

If distances are evaluated using the minimum image convention on an unwrapped system, raw coordinates must be within a distance of 1.5x the box length of each other to ensure correct results.

CalcRMSD(ref_coords, atom_inds=[], ref_atom_inds=[], correction=:pbc, has_virial=true)

Bias the root-mean-square deviation (RMSD) between the coordinates of a group of atoms and a set of reference coordinates.

Given as an argument to BiasPotential. The two sets of coordinates are superimposed using the Kabsch algorithm, as for rmsd.

Arguments

• ref_coords: reference coordinates.
• atom_inds=[]: indices of the atoms in the group, [] uses all atoms.
• ref_atom_inds=[]: indices of the reference coordinates to use, [] uses all coordinates.
• correction=:pbc: the correction to be applied to the molecules. :pbc keeps molecules whole, :wrap wraps all atoms inside the simulation box. Generally atoms in a group should be in the same molecule and :pbc should be used.

CalcRg(atom_inds=[], correction=:pbc, has_virial=true)

Bias the radius of gyration of a group of atoms.

Given as an argument to BiasPotential. The radius of gyration is calculated as for radius_gyration.

Arguments

• atom_inds=[]: indices of the atoms in the group, [] uses all atoms.
• correction=:pbc: the correction to be applied to the molecules. :pbc keeps molecules whole, :wrap wraps all atoms inside the simulation box. Generally atoms in a group should be in the same molecule and :pbc should be used.
• has_virial=true: whether the virial contribution should be calculated for this collective variable.

CalcSingleDist(calc_type=:closest)

Bias the distance between two atoms.

Given as an argument to CalcDist. By default, distances are calculated between the closest periodic images. Setting calc_type=:raw means that distances are calculated ignoring PBCs.

If distances are evaluated using the minimum image convention on an unwrapped system, raw coordinates must be within a distance of 1.5x the box length of each other to ensure correct results.

CalcTorsion(atom_inds::AbstractVector{Int}=[], correction=:pbc, has_virial::Bool=true;
            gradient_singularity_tol=1e-6)

A collective variable that calculates the torsion angle (dihedral) defined by four atoms.

Given as an argument to BiasPotential. The angle is the angle between the planes formed by atoms (i, j, k) and (j, k, l), where the indices are given by atom_inds, as for torsion_angle. The torsion gradient is regularized near collinear geometries using gradient_singularity_tol, a dimensionless relative tolerance applied to the bond-vector norms.

Fields

• atom_inds::AbstractVector{Int}: The indices of the four atoms (i, j, k, l) defining the torsion.
• correction::Symbol: The method used to handle periodic boundary conditions. Defaults to :pbc.
• has_virial::Bool: Whether the virial contribution should be calculated for this collective variable. Defaults to true.
• gradient_singularity_tol::Float64: Relative tolerance used to cap torsion gradients near collinear geometries.

CellListMapNeighborFinder(; n_atoms, dist_cutoff, boundary, excluded_pairs=(),
                            special_pairs=(), auto_update=true, n_steps=10,
                            x0=nothing, number_of_batches=(0, 0), strictness=:warn)
CellListMapNeighborFinder(; eligible, dist_cutoff, boundary, special=nothing,
                            auto_update=true, n_steps=10, x0=nothing,
                            number_of_batches=(0, 0), strictness=:warn)

Find close atoms by distance using a cell list algorithm from CellListMap.jl.

This is the recommended neighbor finder on CPU. dist_cutoff is the neighbor search distance, which should be the interaction cutoff distance plus a buffer distance since the list is reused until it is rebuilt. auto_update and n_steps determine when the list is rebuilt, as for DistanceNeighborFinder. The pairs of atoms that can interact are given by n_atoms, excluded_pairs and special_pairs, or by the eligible and special matrices, as for DistanceNeighborFinder. The pairs are stored on the CPU. x0 are optional initial coordinates that improve the first approximation of the cell list structure. The number of dimensions dims is inferred from the boundary or x0, or assumed to be 3 otherwise.

The boundary parameter is required, and must be a Molly.AbstractBoundary. Infinite boundaries are only accepted in all dimensions.

CellListMap.jl chooses how many batches to split the work over from the number of Julia threads, so the n_threads argument to find_neighbors only selects between a serial (n_threads=1) and a parallel run. number_of_batches can be given as a tuple of the number of batches used to build the cell lists and to map over the pairs, which caps the number of tasks used, with (0, 0), the default, meaning to use the CellListMap.jl heuristics.

CosineAngle(; k, θ0)

A cosine bond angle between three atoms.

θ0 is in radians. The potential energy is defined as

\[V(\theta) = k(1 + \cos(\theta - \theta_0))\]

CosineAngleλ(; k, θ0, λ_mixing=Molly.MinimumMixing(),
             scheduler=Molly.DefaultLambdaScheduler())

A cosine bond angle between three atoms in an alchemical system, the λ version of CosineAngle.

AbsoluteFESystem and RelativeFESystem convert each CosineAngle to this type. The potential energy is defined as

\[V(\theta) = w k(1 + \cos(\theta - \theta_0))\]

where θ0 is as for CosineAngle, and the weight $w$ and the parameters depend on the alchemical roles of the atoms and on $\lambda$, the λ of the atoms mixed with λ_mixing. scheduler is a Molly.LambdaScheduler, which selects dual or single topology.

In dual topology k and θ0 are single values. $w = 1$ if the angle involves an atom unique to one state or only environment atoms, otherwise $w = 1 - \lambda$ if it involves a core atom with the parameters of state A and $w = \lambda$ if it involves a core atom with the parameters of state B.

In single topology $w = 1$ and k and θ0 are tuples of the values in states A and B. Each parameter $p$ is $(1 - \lambda) p^A + \lambda p^B$ if the angle involves a core or unique atom and $p^A$ otherwise.

Coulomb(; cutoff=NoCutoff(), use_neighbors=true, weight_special=1,
        coulomb_const=Molly.coulomb_const)

The Coulomb electrostatic interaction between two atoms.

The potential energy is defined as

\[V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}}\]

where coulomb_const is $1 / 4 \pi \varepsilon_0$. The interaction between special pairs is scaled by weight_special.

CoulombEwald(; dist_cutoff, error_tol=0.0005, ϵr=1, use_neighbors=true,
             weight_special=1, coulomb_const=Molly.coulomb_const, approximate_erfc=true)

The short range Ewald electrostatic interaction between two atoms.

Should be used alongside the Ewald or PME general interaction, which provide the long-range term, and the EwaldExclusion specific interaction, which provides the exclusions for bonded atoms. dist_cutoff, error_tol, ϵr and coulomb_const should match these interactions.

The potential energy is defined as

\[V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 \varepsilon_r} \frac{\mathrm{erfc}(\alpha r_{ij})}{r_{ij}}, \quad \alpha = \frac{\sqrt{-\ln(2 \epsilon_{\mathrm{tol}})}}{r_c}\]

for $r_{ij} \le r_c$, where $r_c$ is dist_cutoff, the cutoff distance for short range interactions, $\epsilon_{\mathrm{tol}}$ is error_tol and $\varepsilon_r$ is ϵr, the relative permittivity. Special pairs use the undamped Coulomb interaction scaled by weight_special. approximate_erfc determines whether to use a fast approximation to the erfc function.

CoulombEwaldScaled(; dist_cutoff, error_tol=0.0005, use_neighbors=true,
                   λ_mixing=Molly.MinimumMixing(),
                   scheduler=Molly.DefaultLambdaScheduler(), weight_special=1,
                   coulomb_const=Molly.coulomb_const, approximate_erfc=true)

The short-range Ewald electrostatic interaction between two atoms scaled by λ without a soft core, the λ version of CoulombEwald.

AbsoluteFESystem and RelativeFESystem convert each CoulombEwald to this type with coulomb_softcore=:scaled. The potential energy is that of CoulombEwald multiplied by $\lambda$, the electrostatic coupling of the pair, which scheduler, a Molly.LambdaScheduler, gives from the λ of the atoms mixed with λ_mixing and their alchemical roles. Should be used alongside Ewald or PME configured with the same electrostatic scheduler.

CoulombReactionField(; dist_cutoff, solvent_dielectric=78.3, solute_dielectric=1.0,
                     use_neighbors=true, weight_special=1,
                     coulomb_const=Molly.coulomb_const)

The Coulomb electrostatic interaction modified using the reaction field approximation between two atoms.

See Tironi et al. 1995. The potential energy is defined as

\[V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 \varepsilon_\mathrm{r}} \left( \frac{1}{r_{ij}} + k_\mathrm{rf} r_{ij}^2 - c_\mathrm{rf} \right)\]

where

\[k_\mathrm{rf} = \frac{1}{r_\mathrm{c}^3} \frac{\varepsilon_\mathrm{rf} - \varepsilon_\mathrm{r}}{2\varepsilon_\mathrm{rf} + \varepsilon_\mathrm{r}}, \quad c_\mathrm{rf} = \frac{1}{r_\mathrm{c}} + k_\mathrm{rf} r_\mathrm{c}^2 = \frac{1}{r_\mathrm{c}} \frac{3\varepsilon_\mathrm{rf}}{2\varepsilon_\mathrm{rf} + \varepsilon_\mathrm{r}}\]

solvent_dielectric corresponds to $\varepsilon_\mathrm{rf}$, the relative permittivity of the continuum beyond the cutoff, and solute_dielectric to $\varepsilon_\mathrm{r}$, the relative permittivity within the cutoff. coulomb_const is divided by solute_dielectric, so it should not include $\varepsilon_\mathrm{r}$ itself. Setting solvent_dielectric=Inf gives conducting boundary conditions ($k_\mathrm{rf} = 1/(2r_\mathrm{c}^3)$, $c_\mathrm{rf} = 3/(2r_\mathrm{c})$). solvent_dielectric must not be less than solute_dielectric. solute_dielectric defaults to one, which gives the reaction field of OpenMM.

Special (1-4) pairs use the plain Coulomb interaction without the reaction field, scaled by weight_special and, like all pairs, divided by solute_dielectric.

SetupCoulombReactionField provides reaction field parameters when setting up a system from a file.

CoulombReactionFieldExclusion(; self=false)

The reaction field correction for an excluded atom pair, or for an atom with itself if self is true.

The dielectric continuum beyond the cutoff responds to all charges within the cutoff, including those of excluded pairs, so the reaction field is applied to excluded pairs and to each atom with itself, without the direct Coulomb term. For an excluded pair within the cutoff the potential energy is

\[V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 \varepsilon_\mathrm{r}} \left( k_\mathrm{rf} r_{ij}^2 - c_\mathrm{rf} \right)\]

and for each atom it is

\[V_\mathrm{self} = -\frac{q_i^2}{4 \pi \varepsilon_0 \varepsilon_\mathrm{r}} \frac{c_\mathrm{rf}}{2}\]

which gives no force. $k_\mathrm{rf}$ and $c_\mathrm{rf}$ are as for CoulombReactionField. Excluded pairs beyond the cutoff contribute nothing, so the energy of an excluded pair jumps by $q_i q_j / (4 \pi \varepsilon_0 \varepsilon_\mathrm{r} r_\mathrm{c})$ at the cutoff. This matches GROMACS, whereas OpenMM does not apply the correction.

Since the properties of the interaction are the same for all pairs, they are stored once in the data field of the InteractionList2Atoms with Molly.CoulombReactionFieldExclusionData. Molly.reaction_field_exclusions sets up the list of excluded pairs and atoms.

CoulombReactionFieldExclusionData(; dist_cutoff, solvent_dielectric=78.3,
                                  solute_dielectric=1.0,
                                  coulomb_const=Molly.coulomb_const,
                                  scheduler=Molly.DefaultLambdaScheduler())
CoulombReactionFieldExclusionData(inter; scheduler)

The parameters of the CoulombReactionFieldExclusion interactions, stored in the data field of their InteractionList2Atoms.

The parameters should match the reaction field pairwise interaction, inter, from which they can be taken. scheduler, a Molly.LambdaScheduler, gives the charges of alchemical atoms as for EwaldExclusion and defaults to that of inter, if any.

CoulombReactionFieldScaled(; dist_cutoff, solvent_dielectric=78.3,
                           solute_dielectric=1.0, use_neighbors=true,
                           λ_mixing=Molly.MinimumMixing(),
                           scheduler=Molly.DefaultLambdaScheduler(), weight_special=1,
                           coulomb_const=Molly.coulomb_const)

The reaction field Coulomb interaction between two atoms scaled by λ without a soft core, the λ version of CoulombReactionField.

AbsoluteFESystem and RelativeFESystem convert each CoulombReactionField to this type with coulomb_softcore=:scaled. The potential energy is that of CoulombReactionField multiplied by $\lambda$, the electrostatic coupling of the pair, which scheduler, a Molly.LambdaScheduler, gives from the λ of the atoms mixed with λ_mixing and their alchemical roles. In single topology the charges are interpolated between the end states.

CoulombScaled(; cutoff=NoCutoff(), use_neighbors=true, λ_mixing=Molly.MinimumMixing(),
              scheduler=Molly.DefaultLambdaScheduler(), weight_special=1,
              coulomb_const=Molly.coulomb_const)

The Coulomb electrostatic interaction between two atoms scaled by λ without a soft core, the λ version of Coulomb.

AbsoluteFESystem and RelativeFESystem convert each Coulomb to this type with coulomb_softcore=:scaled. The potential energy is defined as

\[V(r_{ij}) = \lambda \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}}\]

where $\lambda$ is the electrostatic coupling of the pair, which scheduler, a Molly.LambdaScheduler, gives from the λ of the atoms mixed with λ_mixing and their alchemical roles. In single topology the charges are interpolated between the end states.

CoulombSoftCoreBeutler(; cutoff=NoCutoff(), α=0.3, use_neighbors=true,
                       σ_mixing=Molly.LorentzMixing(), ϵ_mixing=Molly.GeometricMixing(),
                       λ_mixing=Molly.MinimumMixing(),
                       scheduler=Molly.DefaultLambdaScheduler(), weight_special=1,
                       coulomb_const=Molly.coulomb_const)

The Coulomb electrostatic interaction between two atoms with a soft core, used for the appearing and disappearing of atoms.

See Beutler et al. 1994. AbsoluteFESystem and RelativeFESystem convert each Coulomb to this type with coulomb_softcore=:beutler. The potential energy is defined as

\[V(r_{ij}) = \lambda \frac{1}{4\pi\varepsilon_0} \frac{q_iq_j}{r_Q^{1/6}}\]

and the force on each atom by

\[\vec{F}_i = \lambda \frac{1}{4\pi\varepsilon_0} \frac{q_iq_jr_{ij}^5}{r_Q^{7/6}}\frac{\vec{r}_{ij}}{r_{ij}}\]

where

\[r_Q = \alpha(1-\lambda)\sigma_{ij}^6 + r_{ij}^6\]

and $\sigma_{ij}$ is mixed with σ_mixing. $\lambda$ is the electrostatic coupling of the pair, which scheduler, a Molly.LambdaScheduler, gives from the λ of the atoms mixed with λ_mixing and their alchemical roles. If $\lambda$ is 1 this gives the standard Coulomb potential, and if it is 0 the interaction is turned off. $\alpha$ determines the strength of softening the function.

CoulombSoftCoreBeutlerEwald(; dist_cutoff, error_tol=0.0005, α=0.3, use_neighbors=true,
                            σ_mixing=Molly.LorentzMixing(),
                            ϵ_mixing=Molly.GeometricMixing(),
                            λ_mixing=Molly.MinimumMixing(),
                            scheduler=Molly.DefaultLambdaScheduler(), weight_special=1,
                            coulomb_const=Molly.coulomb_const, approximate_erfc=true)

The short-range Ewald electrostatic interaction between two atoms with the Beutler soft-core distance modification.

AbsoluteFESystem and RelativeFESystem convert each CoulombEwald to this type with coulomb_softcore=:beutler. The soft core and $\lambda$, the electrostatic coupling of the pair, are as for CoulombSoftCoreBeutler. Should be used alongside Ewald or PME configured with the same electrostatic scheduler. Special pairs are excluded from the reciprocal-space calculation, so their interaction uses the undamped soft-core form.

CoulombSoftCoreBeutlerReactionField(; dist_cutoff, solvent_dielectric=78.3,
                                    solute_dielectric=1.0, α=0.3, use_neighbors=true,
                                    σ_mixing=Molly.LorentzMixing(),
                                    ϵ_mixing=Molly.GeometricMixing(),
                                    λ_mixing=Molly.MinimumMixing(),
                                    scheduler=Molly.DefaultLambdaScheduler(),
                                    weight_special=1, coulomb_const=Molly.coulomb_const)

The Coulomb electrostatic interaction with Beutler soft core and reaction field correction, used for the appearing and disappearing of atoms in alchemical simulations.

AbsoluteFESystem and RelativeFESystem convert each CoulombReactionField to this type with coulomb_softcore=:beutler. The soft core and $\lambda$, the electrostatic coupling of the pair, are as for CoulombSoftCoreBeutler. At $\lambda = 1$ this reduces to CoulombReactionField. At $\lambda = 0$ the interaction is zero. Special (1-4) pairs disable the reaction field correction and use plain soft-core Coulomb scaled by weight_special and divided by solute_dielectric.

CoulombSoftCoreGapsys(; cutoff=NoCutoff(), α=0.3, σQ=1.0u"nm", use_neighbors=true,
                      λ_mixing=Molly.MinimumMixing(),
                      scheduler=Molly.DefaultLambdaScheduler(), weight_special=1,
                      coulomb_const=Molly.coulomb_const)

The Coulomb electrostatic interaction between two atoms with a soft core, used for the appearing and disappearing of atoms.

See Gapsys et al. 2012. AbsoluteFESystem and RelativeFESystem convert each Coulomb to this type with coulomb_softcore=:gapsys. The potential energy is defined as

\[V(r_{ij}) = \left\{ \begin{array}{cl} \lambda \frac{1}{4\pi\varepsilon_0} \frac{q_iq_j}{r_{ij}}, & \text{if} & r_{ij} \ge r_Q \\ \lambda \frac{1}{4\pi\varepsilon_0} (\frac{q_iq_j}{r_Q^3}r_{ij}^2-\frac{3q_iq_j}{r_Q^2}r_{ij}+\frac{3q_iq_j}{r_Q}), & \text{if} & r_{ij} \lt r_Q \\ \end{array} \right.\]

and the force on each atom by

\[\vec{F}_i = \left\{ \begin{array}{cl} \lambda \frac{1}{4\pi\varepsilon_0} \frac{q_iq_j}{r_{ij}^2}\frac{\vec{r}_{ij}}{r_{ij}}, & \text{if} & r_{ij} \ge r_Q \\ \lambda \frac{1}{4\pi\varepsilon_0} (\frac{-2q_iq_j}{r_Q^3}r_{ij}+\frac{3q_iq_j}{r_Q^2})\frac{\vec{r}_{ij}}{r_{ij}}, & \text{if} & r_{ij} \lt r_Q \\ \end{array} \right.\]

where

\[r_Q = \alpha(1-\lambda)^{1/6}(1+\sigma_Q|q_i q_j|)\]

and $\sigma_Q$ is σQ. $\lambda$ is the electrostatic coupling of the pair, which scheduler, a Molly.LambdaScheduler, gives from the λ of the atoms mixed with λ_mixing and their alchemical roles. If $\lambda$ is 1 this gives the standard Coulomb potential, and if it is 0 the interaction is turned off. $\alpha$ determines the strength of softening the function.

CoulombSoftCoreGapsysEwald(; dist_cutoff, error_tol=0.0005, α=0.3, σQ=1.0u"nm",
                           use_neighbors=true, λ_mixing=Molly.MinimumMixing(),
                           scheduler=Molly.DefaultLambdaScheduler(), weight_special=1,
                           coulomb_const=Molly.coulomb_const, approximate_erfc=true)

The short-range Ewald electrostatic interaction between two atoms with the Gapsys soft-core polynomial near the origin.

AbsoluteFESystem and RelativeFESystem convert each CoulombEwald to this type with coulomb_softcore=:gapsys. The soft core and $\lambda$, the electrostatic coupling of the pair, are as for CoulombSoftCoreGapsys. Should be used alongside Ewald or PME configured with the same electrostatic scheduler. Special pairs are excluded from the reciprocal-space calculation, so their interaction uses the undamped soft-core form.

CoulombSoftCoreGapsysReactionField(; dist_cutoff, solvent_dielectric=78.3,
                                   solute_dielectric=1.0, α=0.3, σQ=1.0u"nm",
                                   use_neighbors=true, λ_mixing=Molly.MinimumMixing(),
                                   scheduler=Molly.DefaultLambdaScheduler(),
                                   weight_special=1, coulomb_const=Molly.coulomb_const)

The Coulomb electrostatic interaction with Gapsys soft core and reaction field correction, used for the appearing and disappearing of atoms in alchemical simulations.

AbsoluteFESystem and RelativeFESystem convert each CoulombReactionField to this type with coulomb_softcore=:gapsys. The soft core and $\lambda$, the electrostatic coupling of the pair, are as for CoulombSoftCoreGapsys. At $\lambda = 1$ this reduces to CoulombReactionField. At $\lambda = 0$ the interaction is zero. Special (1-4) pairs disable the reaction field correction and use the standard Gapsys polynomial scaled by weight_special and divided by solute_dielectric.

CubicBoundary(x, y, z)
CubicBoundary(x)

Cubic 3D bounding box defined by three side lengths.

If one length is given then all three sides will have that length. Setting one or more values to Inf gives no boundary in that dimension.

CubicSplineCutoff(dist_activation, dist_cutoff)

Cutoff that interpolates between the true potential at an activation distance and zero at a cutoff distance using a cubic Hermite spline.

\[\begin{aligned} V_c(r) &= \begin{cases} V(r), r \le r_a \\ (2t^3 - 3t^2 + 1) V(r_a) + (t^3 - 2t^2 + t) (r_c - r_a) V'(r_a), r_a < r \le r_c \\ 0, r > r_c \end{cases} \\ F_c(r) &= \begin{cases} F(r), r \le r_a \\ \frac{-(6t^2 - 6t) V(r_a)}{r_c - r_a} - (3t^2 - 4t + 1) V'(r_a), r_a < r \le r_c \\ 0, r > r_c \end{cases} \\ t &= \frac{r - r_a}{r_c - r_a} \end{aligned}\]

DExpDispersionCorrection(atoms, cutoff, α, β, σ_mixing=Molly.LorentzMixing(),
                         ϵ_mixing=Molly.GeometricMixing())

The long-range dispersion correction for the DoubleExponential interaction, the analogue of LJDispersionCorrection.

Use it with the same α, β, mixing rules and cutoff as the pairwise interaction. cutoff can be a DistanceCutoff (or a distance), a ShiftedPotentialCutoff or a PolynomialCutoff. Assuming a uniform density beyond the cutoff, the potential energy for a DistanceCutoff is

\[E = \frac{2 \pi N^2}{V} \left< \int_{r_c}^\infty r^2 V_{ij}(r) \, dr \right> = \frac{2 \pi N^2}{V} \left< \frac{\varepsilon_{ij} r_m^3}{\alpha - \beta} \left[ \beta G_\alpha\left(\frac{r_c}{r_m}\right) - \alpha G_\beta\left(\frac{r_c}{r_m}\right) \right] \right>\]

where

\[G_k(x) = e^{k(1 - x)} \left( \frac{x^2}{k} + \frac{2x}{k^2} + \frac{2}{k^3} \right)\]

$V_{ij}$ is the pair potential, $V$ the volume and $r_c$ the cutoff distance. The averages run over all pairs, including excluded pairs and self pairs. The other cutoffs also include the switched-out or shifted part of the potential, as for LJDispersionCorrection. The forces are zero, with an isotropic virial contribution.

The number of atoms is assumed not to change after setup. inject_gradients recomputes the correction from injected σ, ϵ, α and β, which share their keys with the pairwise interaction, including its gradients. Only compatible with 3D systems.

DExpDispersionCorrectionλ(atoms, cutoff, α, β, λ_mix, σ_mix, ϵ_mix; scheduler)

The long-range dispersion correction for the DoubleExponentialSoftCore interaction, the λ version of DExpDispersionCorrection.

AbsoluteFESystem and RelativeFESystem convert each DExpDispersionCorrection to this type. For a DistanceCutoff the potential energy is

\[E = \frac{2 \pi N^2}{V} \left< \lambda_{ij} \int_{r_c}^\infty r^2 V_{ij}(r; \alpha_s, \beta_s) \, dr \right>\]

where $V_{ij}(r; \alpha_s, \beta_s)$ is the pair potential with the soft-core shape parameters of DoubleExponentialSoftCore. Unlike for Lennard-Jones, the soft core changes the potential at every distance, so the coupling enters through $\alpha_s$ and $\beta_s$ as well as $\lambda_{ij}$. The averages, $\lambda_{ij}$ and the effective number of atoms $N$ are as for Molly.LJDispersionCorrectionλ, and the other cutoffs as for DExpDispersionCorrection. The forces are zero, with an isotropic virial contribution.

The factors are computed from the λ of atoms at setup; Molly.set_lambda! recomputes them. The number of atoms is otherwise assumed not to change after setup, and inject_gradients takes up injected σ, ϵ, α and β as for DExpDispersionCorrection. Only compatible with 3D systems.

DPDInteraction(; a=25.0, γ=4.5, σ=3.0, r_c=1.0, dt=0.01, use_neighbors=true,
               key=rand(UInt64))

The dissipative particle dynamics (DPD) interaction between two particles.

Combines conservative, dissipative, and random pairwise forces as described in Groot and Warren 1997.

The total pairwise force is

\[\vec{F}_{ij} = \vec{F}_{ij}^C + \vec{F}_{ij}^D + \vec{F}_{ij}^R\]

where the conservative force is

\[\vec{F}_{ij}^C = a \left(1 - \frac{r_{ij}}{r_c}\right) \hat{r}_{ij}\]

the dissipative force is

\[\vec{F}_{ij}^D = -\gamma \left(1 - \frac{r_{ij}}{r_c}\right)^2 (\hat{r}_{ij} \cdot \vec{v}_{ij}) \hat{r}_{ij}\]

and the random force is

\[\vec{F}_{ij}^R = \sigma \left(1 - \frac{r_{ij}}{r_c}\right) \xi_{ij} \Delta t^{-1/2} \hat{r}_{ij}\]

All forces are zero for $r_{ij} \geq r_c$. The weight functions satisfy $w^D(r) = [w^R(r)]^2$ with $w^R(r) = 1 - r/r_c$.

The fluctuation-dissipation relation $\sigma^2 = 2 \gamma k_B T$ must be satisfied by the user's choice of parameters to correctly thermostat at temperature $T$.

The conservative potential energy is

\[V(r_{ij}) = \frac{a}{2} r_c \left(1 - \frac{r_{ij}}{r_c}\right)^2\]

Pairwise random numbers are generated from a hash of the key, the particle indices, and the step number, ensuring momentum conservation.

Should be used with a simulator like DPDVelocityVerlet. When using a neighbor list, set the dist_cutoff of the neighbor finder to at least 1.5 * r_c to provide a skin distance that accounts for particle movement between neighbor list rebuilds.

Arguments

• a: the conservative force strength parameter.
• γ: the dissipative force strength parameter.
• σ: the random force strength parameter (must satisfy σ² = 2γkBT).
• r_c: the cutoff distance beyond which all forces are zero.
• dt: the simulation timestep (needed for random force scaling as Δt⁻¹ᐟ²).
• use_neighbors::Bool=true: whether to use the neighbor list.
• key::UInt64=rand(UInt64): the seed for the random forces.

DPDVelocityVerlet(; <keyword arguments>)

The modified velocity Verlet integrator for dissipative particle dynamics.

Implements the Groot-Warren modified velocity-Verlet (MVV) algorithm from Groot and Warren 1997.

Because DPD dissipative forces depend on particle velocities, a velocity prediction step is used before recomputing forces at the new positions. The λ parameter controls this prediction: v_predicted = v(t) + λ * dt * a(t). A value of 0.65 is commonly used.

Should be used with DPDInteraction as the pairwise interaction. Does not use the provided rng during simulation.

Arguments

• dt::T: the time step of the simulation.
• λ::L=0.65: the velocity prediction parameter (typically 0.5–0.65).
• coupling::C=nothing: the coupling which applies during the simulation.
• remove_CM_motion=1: remove the center of mass motion every this number of steps, set to false or 0 to not remove center of mass motion.

DefaultLambdaScheduler(; dual=true, LJindividual=false, LJspecial=false,
                       Cindividual=false, Cspecial=false, intraLJ=false,
                       Tscaled=true)

Lambda scheduler that turns off the electrostatics of an atom before its sterics.

Atoms that only exist in the first state lose their charges between global_λ = 0 and 0.5 and their Lennard-Jones interactions between 0.5 and 1. Atoms that only exist in the second state gain their Lennard-Jones interactions between 0 and 0.5 and their charges between 0.5 and 1, so an atom is never charged without its Lennard-Jones core. Core atoms, present in both states, are interpolated linearly.

A Molly.LambdaScheduler, see there for the keyword arguments.

DefaultSchedule()

Schedule of a Molly.DefaultLambdaScheduler, passed to Molly.LambdaScheduler.

DisplacementsLogger(n_steps, coords_start; n_steps_update::Integer=10)

Log the displacements of atoms in a system throughout a simulation, useful for calculating properties like mean square displacement in periodic systems.

The displacement accumulates as an atom moves over the periodic boundary. Displacements are updated every n_steps_update steps and a copy is saved every n_steps steps. coords_start are the initial reference positions and should match the coordinate type in the system.

It is assumed that a particle does not cross half the box size in n_steps_update steps. By default n_steps_update is set to 10 to mitigate this assumption, but it can be set to a higher value to reduce cost. n_steps must be a multiple of n_steps_update.

DistanceConstraint(i, j, dist)

Constraint between two atoms that maintains a fixed distance between the atoms.

Given to a constraint algorithm such as SHAKE_RATTLE or LINCS.

DistanceCutoff(dist_cutoff)

Cutoff that sets the potential and force to be zero past a specified cutoff distance.

\[\begin{aligned} V_c(r) &= \begin{cases} V(r), r \le r_c \\ 0, r > r_c \end{cases} \\ F_c(r) &= \begin{cases} F(r), r \le r_c \\ 0, r > r_c \end{cases} \end{aligned}\]

DistanceNeighborFinder(; n_atoms, dist_cutoff, excluded_pairs=(), special_pairs=(),
                       auto_update=true, n_steps=10, array_type=nothing,
                       strictness=:warn)
DistanceNeighborFinder(; eligible, dist_cutoff, special=nothing, auto_update=true,
                       n_steps=10, array_type=nothing, strictness=:warn)

Find close atoms by distance.

This is the recommended neighbor finder on non-NVIDIA GPUs when the box is too small for GPUCellListNeighborFinder. It checks every pair of atoms, so the time it takes grows with the square of the number of atoms.

dist_cutoff is the neighbor search distance, which should be the interaction cutoff distance plus a buffer distance since the list is reused until it is rebuilt. With auto_update=true, the default, the list is rebuilt when the atoms have moved far enough since the last rebuild that a pair within the cutoff of the interactions that use the list could be missing from it, and n_steps is ignored. With auto_update=false the list is rebuilt every n_steps steps, so the buffer distance should be larger than the distance an atom can move in that time, or every step with n_steps=1.

The pairs of atoms that can interact are described by n_atoms, the number of atoms, along with excluded_pairs, the pairs excluded from the pairwise interactions such as bonded atoms, and special_pairs, the pairs that use the special interaction path such as 1-4 atoms. These are iterables of (i, j) pairs or sparse matrices. They are stored as SparsePairMatrixs, so the memory used is proportional to the number of pairs. array_type is the array type used to store them, Array by default, and should match the System that the neighbor finder is used with, for example CuArray for a system on a NVIDIA GPU.

Alternatively, eligible can be given as a symmetric n_atoms x n_atoms Boolean matrix that is true for the pairs that can interact, along with special as a matrix of the same size that is true for the special pairs, or nothing if no pairs are special. If array_type is given the matrices are moved to it. Dense matrices take memory proportional to n_atoms^2, so they are only suitable for small systems. strictness determines the behavior when a dense matrix is given for more than 50,000 atoms, options are :warn to emit a warning, :nowarn to suppress it or :error to error.

DoubleExponential(; cutoff=NoCutoff(), use_neighbors=true, α, β,
                  shortcut=Molly.LJZeroShortcut(), σ_mixing=Molly.LorentzMixing(),
                  ϵ_mixing=Molly.GeometricMixing(), weight_special=1)

The double exponential interaction between two atoms.

See Wu and Brooks 2019. The potential energy is defined as

\[V(r_{ij}) = \varepsilon_{ij} \left[ \frac{\beta e^\alpha}{\alpha - \beta} e^{-\alpha r_{ij} / r_m} - \frac{\alpha e^\beta}{\alpha - \beta} e^{-\beta r_{ij} / r_m} \right]\]

and the force on each atom by

\[\vec{F}_i = \frac{\varepsilon_{ij}}{r_m} \left[ \frac{\alpha\beta e^\alpha}{\alpha - \beta} e^{-\alpha r_{ij} / r_m} - \frac{\alpha\beta e^\beta}{\alpha - \beta} e^{-\beta r_{ij} / r_m} \right] \frac{\vec{r}_{ij}}{r_{ij}}\]

where

\[r_m = 2^{1/6}\sigma_{ij}\]

and $\sigma_{ij}$ and $\varepsilon_{ij}$ are mixed with σ_mixing and ϵ_mixing.

DoubleExponentialSoftCore(; cutoff=NoCutoff(), use_neighbors=true, α, β,
                          shortcut=Molly.LJZeroShortcut(),
                          σ_mixing=Molly.LorentzMixing(),
                          ϵ_mixing=Molly.GeometricMixing(),
                          λ_mixing=Molly.MinimumMixing(),
                          scheduler=Molly.DefaultLambdaScheduler(), weight_special=1)

The double exponential interaction between two atoms with a soft core, used for the appearing and disappearing of atoms.

See Horton et al. 2023. The potential energy is defined as

\[V(r_{ij}) = \lambda\varepsilon_{ij} \left[ \frac{\beta_s e^{\alpha_s}}{\alpha_s - \beta_s} e^{-\alpha_s r_{ij} / r_m} - \frac{\alpha_s e^{\beta_s}}{\alpha_s - \beta_s} e^{-\beta_s r_{ij} / r_m} \right]\]

and the force on each atom by

\[\vec{F}_i = \frac{\lambda\varepsilon_{ij}}{r_m} \left[ \frac{\alpha_s\beta_s e^{\alpha_s}}{\alpha_s - \beta_s} e^{-\alpha_s r_{ij} / r_m} - \frac{\alpha_s\beta_s e^{\beta_s}}{\alpha_s - \beta_s} e^{-\beta_s r_{ij} / r_m} \right] \frac{\vec{r}_{ij}}{r_{ij}}\]

where

\[r_m = 2^{1/6}\sigma_{ij}\]

\[\alpha_s = 1.1 + \lambda(\alpha - 1.1)\]

\[\beta_s = 1 + \lambda(\beta - 1)\]

$\lambda$ is the steric coupling of the pair, which scheduler, a Molly.LambdaScheduler, gives from the λ of the atoms mixed with λ_mixing and their alchemical roles. If $\lambda$ is 1 this gives the standard DoubleExponential potential, and if it is 0 the interaction is turned off. $\alpha$ and $\beta$ determine the shape of the double exponential potential.

EleScaledLambdaScheduler(; dual=true, LJindividual=false, LJspecial=false,
                         Cindividual=false, Cspecial=false, intraLJ=false,
                         Tscaled=true)

Lambda scheduler with the stages of Molly.DefaultLambdaScheduler and a non-linear electrostatic coupling.

The charges of atoms that only exist in the first state are coupled by 1 - (2 global_λ)^2 between global_λ = 0 and 0.5, and those of atoms that only exist in the second state by sqrt(2 (global_λ - 0.5)) between 0.5 and 1. The charges therefore change slowly close to full coupling and quickly close to decoupling.

A Molly.LambdaScheduler, see there for the keyword arguments.

EleScaledSchedule()

Schedule of a Molly.EleScaledLambdaScheduler, passed to Molly.LambdaScheduler.

ElectronegativityEqualizationMethod(coupling, type_ids, χ, η; charges=nothing, tolerance=1e-6, max_iterations=200)

Electronegativity equalization method (EEM) for calculating partial charges on atoms (Mortier et al. 1986 and Rappé and Goddard 1991).

type_ids contains a one-based type index for each atom. χ and η contain the electronegativity and hardness of each type in eV, with charges expressed in elementary-charge units.

coupling supplies any required off-diagonal Coulomb coefficients through Molly.eem_hij.

χ and η must have the same floating point precision, which also sets the overall solver precision. All input arrays must be on the same device already, and charges, if provided, must use the same precision also. charges is the initial per-atom charge vector and is used as mutable storage for the solver. Its sum also sets the total charge that is preserved during charge equilibration.

tolerance controls the convergence of the projected preconditioned residual, scaled by the norm of the electronegativities vector. If the solver fails to converge within max_iterations, it raises an error and restores the charges from before the solve.

Include this calculator in a system's general_inters interactions block to equilibrate charges before any pairwise and specific interactions are evaluated. Energy and force evaluations write the converged charges from the solve into sys.atoms. EEM tallies the charge equilibration self-energy into the global energy, but the corresponding electrostatic pair energy and forces are supplied by another interaction. For ReaxFF, use ElectronegativityEqualizationMethod(reaxff, atoms_data; ...).

ElectronegativityEqualizationMethod(reaxff::ReaxFF, atoms_data; charges=nothing, tolerance=1e-6,
    eem_float_type=Float64, max_iterations=200)

Construct ElectronegativityEqualizationMethod charge equilibration using a ReaxFF force field's electronegativities, hardnesses, and tapered, shielded Coulomb interaction.

atoms_data contains an entry for each atom with an element field matching the force field. Arrays are placed on the same device as reaxff, with solver precision set by eem_float_type. charges gives the initial numerical charges in elementary-charge units and defaults to zeros; their sum sets the conserved total charge. tolerance and max_iterations control the EEM solver.

System(coord_file_path, reaxff; ...) calls this constructor automatically.

EnsembleSystem(coordinate_file, trajectory_file, force_field; <keyword arguments>)
EnsembleSystem(system, trajectory_file)

An object allowing data to be read from a trajectory or ensemble associated with a System.

The keyword arguments are the same as System setup from a file. In the case of passing a System directly, a copy of the system is made.

Energies calculated from the frames are affected by the precision of the coordinates in the trajectory file, which is reduced in many formats such as PDB, XTC and DCD.

Ewald(dist_cutoff; error_tol=0.0005, ϵr=1.0, coulomb_const=Molly.coulomb_const,
      scheduler=Molly.DefaultLambdaScheduler())

Ewald summation for long range electrostatics implemented as an AtomsCalculators.jl calculator.

Should be used alongside the CoulombEwald pairwise interaction, which provides the short range term, and the EwaldExclusion specific interaction, which provides the exclusions for bonded atoms. dist_cutoff, error_tol, ϵr and coulomb_const should match these interactions. SetupEwald provides Ewald parameters when setting up a system from a file.

dist_cutoff is the cutoff distance for short range interactions, ϵr is the relative permittivity and coulomb_const is the Coulomb constant. The charge of each atom is scaled by its electrostatic coupling, which scheduler, a Molly.LambdaScheduler, gives from the λ and alchemical role of the atom; this has no effect for atoms that are not alchemical. This algorithm is O(N^2) and in general PME should be used instead. Only compatible with 3D systems and CubicBoundary. Not compatible with infinite boundaries. Runs on the CPU, even for GPU systems.

EwaldExclusion()

Exclusions for bonded interactions for long range electrostatics.

Should be used alongside the Ewald or PME general interaction, which provide the long-range term, and the CoulombEwald pairwise interaction, which provides the short range term.

Since the properties of the interaction are the same for all pairs, they are stored once in the data field of the InteractionList2Atoms with EwaldExclusionData. dist_cutoff, error_tol, ϵr and coulomb_const for EwaldExclusionData should match the above interactions.

Only compatible with 3D systems.

ExceptionList(dict)
ExceptionList(keys, values)

Define a set of atom pairs that have a specified mixing value.

For example, giving Dict((1, 2) => 0.3u"nm") means that pairs of atom type 1 and 2 have a σ value of 0.3 nm. Provide when setting up a Molly.MixingException.

FENEBond(; k, r0, σ, ϵ)

A finitely extensible non-linear elastic (FENE) bond between two atoms, see Kremer and Grest 1990.

The potential energy is defined as

\[V(r) = -\frac{1}{2} k r^2_0 \ln \left( 1 - \left( \frac{r}{r_0} \right) ^2 \right) + V_{\text{WCA}}(r)\]

where the WCA contribution is given by

\[V_{\text{WCA}}(r) = \begin{cases} 4\varepsilon \left[ \left( \frac{\sigma}{r} \right) ^{12} - \left( \frac{\sigma}{r} \right) ^6 \right] + \varepsilon & r < 2^{1/6}\sigma\\ 0 & r \geq 2^{1/6}\sigma\\ \end{cases}\]

FIREMinimizer(; <keyword arguments>)

Energy minimization using the fast inertial relaxation engine (FIRE).

Molecular dynamics is run with an adaptive time step, with the velocity mixed towards the force direction when the system moves downhill and reset to zero when it moves uphill, as described in Bitzek et al. 2006. Internal velocities starting from zero are used, so the velocities of the system are not modified.

The potential energy is not required by the algorithm, so it is only calculated when log_stream is set or the system has loggers, in which case loggers requiring the potential energy calculate it themselves.

Arguments

• dt::D=0.001u"ps": the initial time step.
• dt_max::DM=0.01u"ps": the maximum time step.
• max_steps::Integer=1_000: the maximum number of steps.
• tol::F=1_000.0u"kJ * mol^-1 * nm^-1": the maximum force below which to finish minimization.
• alpha_start::Float64=0.1: the initial velocity mixing parameter, must be between 0 and 1.
• f_alpha::Float64=0.99: the factor to multiply the velocity mixing parameter by when moving downhill.
• f_inc::Float64=1.1: the factor to multiply the time step by when moving downhill.
• f_dec::Float64=0.5: the factor to multiply the time step by when moving uphill.
• n_min::Integer=5: the number of consecutive downhill steps required before the time step is increased.
• constraint_bond_constant::K=500_000.0u"kJ * mol^-1 * nm^-2": the force constant for the harmonic bonds that are used instead of constraints during minimisation. Set to nothing to not use harmonic bonds and ignore constraints. Unused if the system does not have constraints.
• log_stream::L=devnull: stream to print minimization progress to.

FenderHalseyMixing()

Mix ϵ parameters according to the Fender-Halsey rule.

See Fender and Halsey 1962.

\[\varepsilon_{ij} = \frac{2 \varepsilon_i \varepsilon_j}{\varepsilon_i + \varepsilon_j}\]

Give this when setting up a pairwise interaction such as LennardJones.

FlatBottomSquareBias(k, r_fb, cv_target)

A flat-bottomed square (harmonic) bias on a collective variable (CV) towards a target value.

The bias is zero when the value of the collective variable does not deviate from cv_target by more than r_fb, and is square (harmonic) outside this range.

The potential energy is defined as

\[V(\boldsymbol{s}) = \frac{1}{2} k (|\boldsymbol{s} - \boldsymbol{s}_t| - r_{fb})^2 H\]

where $s$ and $s_t$ are the system and target CV values respectively, and

\[H = \left\{ \begin{array}{cl} 0 & \text{if} & |\boldsymbol{s} - \boldsymbol{s}_t| < r_{fb} \\ 1 & \text{if} & |\boldsymbol{s} - \boldsymbol{s}_t| \geq r_{fb} \\ \end{array} \right.\]

Arguments

• k: The energy constant for the bias. Must be compliant with the System energy units.
• r_fb: Width of flat-bottom potential well. Inside this region the bias potential is always 0.
• cv_target: The target value of the collective variable.

ForceFieldXMLError(msg)

An error when reading a force field XML file.

GPUCellListNeighborFinder(; n_atoms, dist_cutoff, excluded_pairs=(), special_pairs=(),
                          auto_update=true, n_steps=10, force_float32=false, array_type,
                          max_neighbors=nothing, output=:molly_pairs, ragged=true,
                          strictness=:warn)
GPUCellListNeighborFinder(; eligible, dist_cutoff, special=nothing, auto_update=true,
                          n_steps=10, force_float32=false, array_type=nothing,
                          max_neighbors=nothing, output=:molly_pairs, ragged=true,
                          strictness=:warn)

Neighbor finder for GPU systems that bins atoms into a uniform grid of cells and searches the 3x3x3 stencil of cells around each atom.

Returns a GPUCellListNeighborList. It runs on any GPU backend and is the recommended neighbor finder for non-CUDA GPUs, where GPUNeighborFinder is not available. Three-dimensional CubicBoundary and TriclinicBoundary systems are supported, provided the grid has at least three cells along every box axis, i.e. opposite box faces are at least three times the internally padded search cutoff apart. For a TriclinicBoundary, that distance is smaller than the length of the corresponding basis vector.

dist_cutoff is the neighbor search distance, which should be the interaction cutoff distance plus a buffer distance since the list is reused until it is rebuilt. With auto_update=true, the default, the list is rebuilt when the atoms have moved far enough since the last rebuild that a pair within the cutoff of the interactions could be missing, and n_steps is ignored. With auto_update=false the list is rebuilt every n_steps steps, so the buffer distance should be larger than the distance an atom can move in that time.

Neighbor finding uses the system's floating-point type by default. Set force_float32=true to use Float32 for neighbor finding.

The neighbor search distance is increased internally by a small margin, based on both neighbor finding precision and box size, to guard against missing pairs due to rounding. This is particularly relevant when neighbor finding is done in Float32 when the system uses Float64. The returned neighbor list may therefore include pairs that are slightly beyond dist_cutoff, so interactions should also use their own physical distance cutoffs. The other neighbor finders do not currently apply this same margin, so different finders can have different behaviours for pairs extremely close to dist_cutoff.

The pairs of atoms that can interact are given by n_atoms, excluded_pairs and special_pairs, or by the eligible and special matrices, as for DistanceNeighborFinder. The exceptions are stored on the device as SparsePairMatrixs, so the memory they use is proportional to the number of excluded and special pairs. array_type is the array type used to store them, for example CuArray, and is required unless eligible is given on the GPU, in which case its array type is used. A SparsePairMatrix given as eligible is used as it is, including one that lists the pairs that can interact rather than the excluded ones. Like GPUNeighborFinder, dense matrices are converted once at construction and are not retained. A missing special matrix means that no pairs are special. strictness determines the behavior when a dense matrix is given for more than 50,000 atoms, options are :warn to emit a warning, :nowarn to suppress it or :error to error.

Output modes

• :molly_pairs: half-pair list with excluded pairs removed and the special flag set for special pairs. This is the only mode that can be used for the pairwise interactions of a System, since it is the only one that applies exclusions.
• :geometric_pairs: half-pair list of every pair within dist_cutoff, with all special flags false. The exceptions are ignored, so pairs that a bonded topology excludes are still returned and using this mode for a system with bonded exclusions gives wrong forces.
• :ragged: per-atom padded neighbor matrix only, read with ragged_neighbors, and no flat pair list. The exceptions are ignored. This mode is for querying neighbors directly, it cannot be used for the pairwise interactions of a System but can be useful for machine learning potentials.

ragged=false finds the pairs of the :molly_pairs and :geometric_pairs modes without storing the per-atom neighbor matrix, by counting the pairs in one pass of the search and writing them in a second. max_neighbors is the initial per-atom capacity of the neighbor matrix. When it is nothing the capacity is estimated from the global number density and dist_cutoff with a safety factor of 1.5, rounded up to a multiple of 32. An explicit positive integer overrides this estimate. Either way the capacity is grown automatically, and kept for later calls, if an atom turns out to have more neighbors than it can hold. The pair list is sized from the number of pairs found, with a little room to grow.

The device buffers behind the returned list are owned by that list and are reused, and hence overwritten, when it is passed back in as current_neighbors. A list that is passed to find_neighbors should therefore be treated as consumed, since any other reference to it, for example one kept by a simulator to compare against a trial list, would see the new neighbors with a stale pair count.

GPUCellListNeighborList(ragged_counts, ragged_neighbors, n, list, state)

The result of find_neighbors with a GPUCellListNeighborFinder, containing a padded per-atom neighbor matrix, unless the finder uses ragged=false, and, unless the finder uses output=:ragged, a flat half-pair list whose first n entries of list are valid.

Use neighbor_pairs to get the pairs and ragged_neighbors to get the per-atom matrix rather than reading the fields, since the layout is not part of the interface. state holds the device buffers behind the other fields and is internal; they are reused if this list is passed back to find_neighbors as current_neighbors, see GPUCellListNeighborFinder for what that means for holding on to a list.

GPUNeighborFinder(; n_atoms, dist_cutoff, excluded_pairs=(), special_pairs=(),
                  auto_update=true, n_steps=10, array_type, initialized=false,
                  strictness=:warn)
GPUNeighborFinder(; eligible, dist_cutoff, special=nothing, auto_update=true,
                  n_steps=10, array_type=nothing, initialized=false, strictness=:warn)

Neighbor finder for CUDA systems that uses tiled pairwise kernels.

See Eastman and Pande 2010. GPUNeighborFinder does not materialize a conventional per-atom neighbor list. Instead, the CUDA pairwise force and energy paths reorder atoms on a Morton curve and build a compact list of interacting 32x32 tiles directly on the device. The exclusions and special pairs are read from per-atom sparse lists inside the tiles that contain them, so the memory used grows linearly with the number of atoms. For that reason, find_neighbors returns nothing for this neighbor finder. With auto_update=true, the default, the force and energy kernels check each time they run whether the atoms have moved far enough since the tile list was built that a pair within the cutoff of the pairwise interactions could be missing, and refresh the tile list if so. With auto_update=false, find_neighbors should still be called at the same points as for the other neighbor finders, since the tile list is rebuilt at the next force or energy evaluation after find_neighbors is called on a step that is a multiple of n_steps or with force_recompute=true.

This is the recommended neighbor finder for CuArray systems. On other GPUs use GPUCellListNeighborFinder. The exceptions are stored in the eligible and special fields as SparsePairMatrixs on the device. Updating the exception pairs should reset initialized so the cached tile data is rebuilt on the next GPU force or energy evaluation.

Keyword arguments

• n_atoms: the number of atoms in the system.
• dist_cutoff: the neighbor search distance used by the pairwise kernels. This should be the interaction cutoff distance plus a buffer distance, since the tile list is reused until it is refreshed.
• excluded_pairs: an iterable of (i, j) pairs, or a sparse matrix, of the pairs that are excluded from the pairwise interactions, for example bonded atoms.
• special_pairs: an iterable of (i, j) pairs, or a sparse matrix, of the pairs that use the special interaction path, for example 1-4 atoms.
• auto_update=true: whether to refresh the tile list only when the atoms have moved far enough since the last refresh that a pair within the cutoff of the pairwise interactions could be missing, in which case n_steps is ignored. The Morton order is recomputed at the same time.
• n_steps=10: with auto_update=false, the number of steps between Morton reorder and tile-list refresh passes. A forced rebuild, for example after a barostat changes the coordinates, refreshes the tile list without changing the Morton order. The buffer distance should then be larger than the distance an atom can move in that time.
• array_type: the array type used to store the excluded and special pairs on the GPU, for example CuArray. It is required unless eligible is given on the GPU, in which case its array type is used.
• eligible, special: eligible and special matrices can be given instead of n_atoms, excluded_pairs and special_pairs, see DistanceNeighborFinder. They are converted once at construction into sparse matrices and are not retained.
• initialized: whether the current sparse-exception preprocessing state can be reused.
• strictness=:warn: determines behavior when a dense eligible or special matrix is given for more than 50,000 atoms, options are :warn to emit warnings, :nowarn to suppress warnings or :error to error.

GROMACSABFESchedule()

Schedule of a Molly.GROMACSLambdaABFEScheduler, passed to Molly.LambdaScheduler.

GROMACSLambdaABFEScheduler(; dual=true, LJindividual=false, LJspecial=false,
                           Cindividual=false, Cspecial=false, intraLJ=false,
                           Tscaled=true)

Lambda scheduler for absolute free energy calculations, following GROMACS.

The charges of the alchemical atoms are scaled between global_λ = 0 and 0.5 and their Lennard-Jones interactions between 0.5 and 1, so decoupled atoms lose their charges first. With AbsoluteFESystem, PME is replaced by the two grid scheme of GROMACS, where the grids of the two end states are mixed with the electrostatic coupling. Inserted atoms follow the same stages and would gain their charges before their Lennard-Jones interactions, so use Molly.DefaultLambdaScheduler for insertions.

A Molly.LambdaScheduler, see there for the keyword arguments.

GROMACSLambdaRBFEScheduler(; dual=true, LJindividual=false, LJspecial=false,
                           Cindividual=false, Cspecial=false, intraLJ=false,
                           Tscaled=true)

Lambda scheduler for relative free energy calculations, following GROMACS.

The electrostatics and sterics are scaled linearly, as in Molly.LinearLambdaScheduler. With RelativeFESystem, PME is replaced by the two grid scheme of GROMACS, where the grids of the two end states are mixed with the electrostatic coupling.

A Molly.LambdaScheduler, see there for the keyword arguments.

GROMACSRBFESchedule()

Schedule of a Molly.GROMACSLambdaRBFEScheduler, passed to Molly.LambdaScheduler.

GeneralObservableLogger(observable::Function, T, n_steps)

A logger which holds a record of regularly sampled observations of a system.

observable should return an object of type T and support the method observable(s::System, neighbors, step_n, buffers; n_threads::Integer)::T.

GeometricMixing()

Mix parameters according to their geometric mean.

For example:

\[\varepsilon_{ij} = \sqrt{\varepsilon_i \varepsilon_j}\]

Give this when setting up a pairwise interaction such as LennardJones.

Gravity(; cutoff=NoCutoff(), G=Unitful.G, use_neighbors=true)

The gravitational interaction between two atoms.

The potential energy is defined as

\[V(r_{ij}) = -\frac{G m_i m_j}{r_{ij}}\]

GridHills(cvs, height, sigma, grid_min, grid_max, n_points; cutoff=6,
          out_of_grid=:global_strictness)

Metadynamics memory for the CVs in the tuple cvs, accumulating deposited hills onto a discretized grid spanning [grid_min, grid_max], evaluated by cubic Hermite interpolation, so the forces are continuous. O(1) per evaluation regardless of the number of deposited hills; see ListHills for an exact but O(n_hills) alternative.

For a single CV, sigma, grid_min and grid_max are scalars and n_points an Integer, giving a 1D grid. For several CVs, pass same-length tuples (one entry per CV in cvs) and n_points as an Integer or a matching tuple, giving a dense N-dimensional grid. Memory and the cost per deposit scale as 2^N times the grid size and the hill footprint, since the exact derivatives of the hills are stored at the grid points, so this is only practical for a handful of CVs (2-3); with 3 CVs on a fine grid, use a coarser grid or a larger n_steps. Use ListHills for more CVs.

Arguments

• cvs: Tuple of CV descriptors, e.g. CalcDist; a single CV may be passed without the tuple.
• height: Default hill height, used unless Molly.add_hill! is given an explicit height.
• sigma: Gaussian width (standard deviation), scalar or one per dimension.
• grid_min, grid_max: Extent of the grid, scalar or one per dimension.
• n_points: Number of grid points, at least 2 in every dimension; a single Integer or one per dimension.
• cutoff=6: Standard deviations beyond which a deposited hill's contribution is ignored.
• out_of_grid::Symbol=:global_strictness: What a CV value outside [grid_min, grid_max] does: :error throws, :warn warns once per grid and CV and :nowarn is silent. Without an error, the bias is held at its edge value along a CV outside the grid, so it exerts no force along that CV, and no hills are deposited. With several CVs, the others still feel the bias at the edge. :global_strictness follows the strictness of the simulation.

The grid of a periodic CV such as a torsion (see Molly.cv_period) must span exactly one period, e.g. -π to π, with grid_max the same point as grid_min: its n_points points are then spaced by period / n_points, hills and interpolation wrap around, and CV values outside [grid_min, grid_max] are wrapped into the grid.

HarmonicAngle(; k, θ0)

A harmonic bond angle between three atoms.

θ0 is in radians. The second atom is the middle atom. The potential energy is defined as

\[V(\theta) = \frac{1}{2} k (\theta - \theta_0)^2\]

HarmonicAngleλ(; k, θ0, λ_mixing=Molly.MinimumMixing(),
               scheduler=Molly.DefaultLambdaScheduler())

A harmonic bond angle between three atoms in an alchemical system, the λ version of HarmonicAngle.

AbsoluteFESystem and RelativeFESystem convert each HarmonicAngle to this type. The potential energy is defined as

\[V(\theta) = \frac{w}{2} k (\theta - \theta_0)^2\]

where θ0 and the atom order are as for HarmonicAngle, and the weight $w$ and the parameters depend on the alchemical roles of the atoms and on $\lambda$, the λ of the atoms mixed with λ_mixing. scheduler is a Molly.LambdaScheduler, which selects dual or single topology.

In dual topology k and θ0 are single values. $w = 1$ if the angle involves an atom unique to one state or only environment atoms, otherwise $w = 1 - \lambda$ if it involves a core atom with the parameters of state A and $w = \lambda$ if it involves a core atom with the parameters of state B.

In single topology $w = 1$ and k and θ0 are tuples of the values in states A and B. Each parameter $p$ is $(1 - \lambda) p^A + \lambda p^B$ if the angle involves a core or unique atom and $p^A$ otherwise.

HarmonicBond(; k, r0)

A harmonic bond between two atoms.

The potential energy is defined as

\[V(r) = \frac{1}{2} k (r - r_0)^2\]

HarmonicBondλ(; k, r0, λ_mixing=Molly.MinimumMixing(),
              scheduler=Molly.DefaultLambdaScheduler())

A harmonic bond between two atoms in an alchemical system, the λ version of HarmonicBond.

AbsoluteFESystem and RelativeFESystem convert each HarmonicBond to this type. The potential energy is defined as

\[V(r) = \frac{w}{2} k (r - r_0)^2\]

where the weight $w$ and the parameters depend on the alchemical roles of the atoms and on $\lambda$, the λ of the atoms mixed with λ_mixing. scheduler is a Molly.LambdaScheduler, which selects dual or single topology.

In dual topology k and r0 are single values. $w = 1$ if the bond involves an atom unique to one state or only environment atoms, otherwise $w = 1 - \lambda$ if it involves a core atom with the parameters of state A and $w = \lambda$ if it involves a core atom with the parameters of state B.

In single topology $w = 1$ and k and r0 are tuples of the values in states A and B. Each parameter $p$ is $(1 - \lambda) p^A + \lambda p^B$ if the bond involves a core or unique atom and $p^A$ otherwise.

HarmonicPositionRestraint(; k, x0)

A harmonic position restraint on an atom to coordinate x0.

The potential energy is defined as

\[V(\boldsymbol{x}) = \frac{1}{2} k |\boldsymbol{x} - \boldsymbol{x}_0|^2\]

Does not contribute to the virial. x0 does not change when the box size changes, so be careful when using this alongside a barostat.

HarmonicTorsion(; k, θ0)

A harmonic torsion angle between four atoms, often used for improper torsions.

θ0 is in radians. The potential energy is defined as

\[V(\theta) = k (\theta - \theta_0)^2\]

where θ is the angle between the planes defined by atoms (i, j, k) and (j, k, l).

Only compatible with 3D systems.

HarmonicTorsionλ(; k, θ0, λ_mixing=Molly.MinimumMixing(),
                 scheduler=Molly.DefaultLambdaScheduler())

A harmonic torsion angle between four atoms in an alchemical system, the λ version of HarmonicTorsion.

AbsoluteFESystem and RelativeFESystem convert each HarmonicTorsion to this type. The potential energy is defined as

\[V(\theta) = w k (\theta - \theta_0)^2\]

where θ and θ0 are as for HarmonicTorsion, and the weight $w$ and the parameters depend on the alchemical roles of the atoms and on $\lambda$, the λ of the atoms mixed with λ_mixing. scheduler is a Molly.LambdaScheduler, which selects dual or single topology.

In dual topology k and θ0 are single values. $w = 1 - \lambda$ if the torsion involves a core atom with the parameters of state A, $w = \lambda$ if it involves a core atom with the parameters of state B and $w = 1$ if it involves only unique and environment atoms.

In single topology $w = 1$ and k and θ0 are tuples of the values in states A and B. Each parameter $p$ is $(1 - \lambda) p^A + \lambda p^B$ if the torsion involves a core atom and $p^A$ otherwise.

Only compatible with 3D systems.

ImmediateThermostat(temperature) <: AbstractThermostat

The immediate velocity rescaling thermostat for controlling temperature.

Velocities are immediately rescaled to match a target temperature (isokinetic ensemble). The scaling factor for the velocities each step is

\[\lambda = \sqrt{\frac{T_0}{T}}\]

This thermostat should be used with caution as it can lead to simulation artifacts.

ImplicitSolventGBN2(atoms, atoms_data, bonds=nothing;
                    solvent_dielectric=Molly.gb_solvent_dielectric,
                    solute_dielectric=Molly.gb_solute_dielectric, kappa=0.0u"nm^-1",
                    offset=Molly.gbn2_offset, dist_cutoff=0.0u"nm",
                    probe_radius=Molly.gb_probe_radius, sa_factor=Molly.gb_sa_factor,
                    use_ACE=true, neck_scale=Molly.gbn2_neck_scale,
                    neck_cut=Molly.gbn2_neck_cut,
                    element_to_radius=Molly.mbondi2_element_to_radius,
                    element_to_screen=Molly.gbn2_element_to_screen,
                    element_to_screen_nucleic=Molly.gbn2_element_to_screen_nucleic,
                    atom_params=Molly.gbn2_atom_params,
                    atom_params_nucleic=Molly.gbn2_atom_params_nucleic,
                    data_d0=Molly.gbn2_data_d0, data_m0=Molly.gbn2_data_m0,
                    coulomb_const=Molly.coulomb_const, n_threads=Threads.nthreads())

GBn2 solvation model implemented as an AtomsCalculators.jl calculator.

See Nguyen et al. 2013. Should be used along with a Coulomb interaction such as Coulomb. SetupImplicitSolventGBN2 provides parameters when setting up a system from a file.

atoms_data gives the element and residue of each atom and bonds the bonded pairs, including any constrained bonds. solvent_dielectric and solute_dielectric are the dielectric constants of the solvent and the solute. kappa is the Debye-Hückel screening parameter, the inverse of the Debye length, and is zero when there is no salt screening. offset is the dielectric offset subtracted from the intrinsic radii. dist_cutoff is the cutoff distance for the Born radius and solvation sums, with zero meaning no cutoff. probe_radius is the radius of the solvent probe and sa_factor the surface tension used by the solvent accessible surface area term, which use_ACE switches on and off. neck_scale scales the neck correction that GBn2 adds to the Born radii and neck_cut is the separation beyond which that correction is zero. element_to_radius maps an element to its intrinsic radius, and element_to_screen and atom_params map an element to its screening parameter and to its α, β and γ parameters, with element_to_screen_nucleic and atom_params_nucleic used instead for atoms in nucleic acid residues. The key "-" gives the value used for elements not in the dictionary. data_d0 and data_m0 are the tabulated neck function parameters, interpolated over the radii of each pair of atoms. coulomb_const is the Coulomb constant, which should match that of the Coulomb interaction. n_threads determines how many per-thread force buffers are allocated.

Only compatible with 3D systems. Not compatible with virial calculation.

ImplicitSolventOBC(atoms, atoms_data, bonds=nothing;
                   solvent_dielectric=Molly.gb_solvent_dielectric,
                   solute_dielectric=Molly.gb_solute_dielectric, kappa=0.0u"nm^-1",
                   offset=Molly.obc_offset, dist_cutoff=0.0u"nm",
                   probe_radius=Molly.gb_probe_radius, sa_factor=Molly.gb_sa_factor,
                   use_ACE=true, use_OBC2=false,
                   element_to_radius=Molly.mbondi2_element_to_radius,
                   element_to_screen=Molly.obc_element_to_screen,
                   coulomb_const=Molly.coulomb_const, n_threads=Threads.nthreads())

Onufriev-Bashford-Case GBSA model implemented as an AtomsCalculators.jl calculator.

See Onufriev et al. 2004. Should be used along with a Coulomb interaction such as Coulomb. SetupImplicitSolventOBC provides parameters when setting up a system from a file.

atoms_data gives the element and residue of each atom and bonds the bonded pairs, including any constrained bonds. solvent_dielectric and solute_dielectric are the dielectric constants of the solvent and the solute. kappa is the Debye-Hückel screening parameter, the inverse of the Debye length, and is zero when there is no salt screening. offset is the dielectric offset subtracted from the intrinsic radii. dist_cutoff is the cutoff distance for the Born radius and solvation sums, with zero meaning no cutoff. probe_radius is the radius of the solvent probe and sa_factor the surface tension used by the solvent accessible surface area term, which use_ACE switches on and off. use_OBC2 determines whether to use parameter set I (false, the default) or II (true). element_to_radius and element_to_screen map an element to its intrinsic radius and its screening parameter, the key "-" giving the value used for elements not in the dictionary. coulomb_const is the Coulomb constant, which should match that of the Coulomb interaction. n_threads determines how many per-thread force buffers are allocated.

Only compatible with 3D systems. Not compatible with virial calculation.

InteractionList1Atoms(is, inters)
InteractionList1Atoms(is, inters, types, data=nothing)
InteractionList1Atoms(inter_type)

A list of specific interactions that involve one atom such as position restraints.

Note that the virial contribution from these interactions is treated as zero, since only one coordinate is available and the virial in periodic space requires displacements.

InteractionList2Atoms(is, js, inters)
InteractionList2Atoms(is, js, inters, types, data=nothing)
InteractionList2Atoms(inter_type)

A list of specific interactions that involve two atoms such as bond potentials.

InteractionList3Atoms(is, js, ks, inters)
InteractionList3Atoms(is, js, ks, inters, types, data=nothing)
InteractionList3Atoms(inter_type)

A list of specific interactions that involve three atoms such as bond angle potentials.

InteractionList4Atoms(is, js, ks, ls, inters)
InteractionList4Atoms(is, js, ks, ls, inters, types, data=nothing)
InteractionList4Atoms(inter_type)

A list of specific interactions that involve four atoms such as torsion potentials.

InteractionList5Atoms(is, js, ks, ls, ms, inters)
InteractionList5Atoms(is, js, ks, ls, ms, inters, types, data=nothing)
InteractionList5Atoms(inter_type)

A list of specific interactions that involve five atoms such as CMAP torsion potentials.

InverseMixing()

Mix B parameters according to their harmonic mean.

\[B_{ij} = \frac{2}{\frac{1}{B_i} + \frac{1}{B_j}}\]

Give this when setting up a pairwise interaction such as Buckingham.

LBFGSMinimizer(; <keyword arguments>)

Energy minimization using the limited-memory BFGS (L-BFGS) algorithm.

A low rank approximation to the inverse Hessian is built from the coordinate and gradient differences of the last n_history steps and applied with the two-loop recursion of Nocedal 1980, giving a quasi-Newton search direction along which a backtracking line search enforcing the Armijo sufficient decrease condition is run. The history is reset when the search direction is not a descent direction or when the line search fails, and minimization stops early if the line search also fails along the steepest descent direction.

Arguments

• step_size::D=0.01u"nm": the maximum displacement of any atom in a trial step.
• max_steps::Integer=1_000: the maximum number of steps.
• tol::F=1_000.0u"kJ * mol^-1 * nm^-1": the maximum force below which to finish minimization.
• n_history::Integer=10: the number of coordinate and gradient difference pairs to store.
• max_line_search_steps::Integer=10: the maximum number of potential energy evaluations in the backtracking line search of each step.
• c1::Float64=1e-4: the Armijo sufficient decrease parameter.
• backtrack_factor::Float64=0.5: the factor to multiply the step length by on each line search backtrack.
• constraint_bond_constant::K=500_000.0u"kJ * mol^-1 * nm^-2": the force constant for the harmonic bonds that are used instead of constraints during minimisation. Set to nothing to not use harmonic bonds and ignore constraints. Unused if the system does not have constraints.
• log_stream::L=devnull: stream to print minimization progress to.

LINCS(; masses, dist_tolerance=1e-8u"nm", vel_tolerance=1e-8u"nm^2 * ps^-1",
      dist_constraints=nothing, angle_constraints=nothing, n_rec=4, n_iter=1,
      iter_vel_correction=false, gpu_block_size=128)

Constrain bond distances during a simulation using the LINCS (LINear Constraint Solver) algorithm.

LINCS is a non-iterative constraint algorithm that uses matrix expansion to approximate the inverse of the constraint coupling matrix. It is typically faster than SHAKE_RATTLE for large systems but is approximate for ring topologies. Either or both of dist_constraints and angle_constraints must be given. The masses and constraints should always be on the CPU, even for a GPU system; the System constructor moves the constraint data to the device of the system. SetupLINCS provides LINCS parameters when setting up a system from a file.

Velocity constraints are applied by projecting out the velocity components along the constraints, see iter_vel_correction. See Hess et al. 1997 for the original LINCS paper.

Not compatible with gradient calculation using Enzyme.

Arguments

• masses: vector of atom masses.
• dist_tolerance=1e-8u"nm": the tolerance for checking position constraints, should have the same units as the coordinates.
• vel_tolerance=1e-8u"nm^2 * ps^-1": the tolerance for checking velocity constraints, should have the same units as the velocities times the coordinates.
• dist_constraints: a vector of DistanceConstraint objects.
• angle_constraints: a vector of AngleConstraint objects. Each angle constraint is converted into 3 distance constraints internally. LINCS requires that angle constraints are isolated: none of their atoms may participate in distance constraints or in other angle constraints.
• n_rec=4: order of the matrix expansion for coupling matrix inversion. Higher values improve accuracy for coupled constraints at the cost of performance.
• n_iter=1: number of outer correction iterations for rotational lengthening. Higher values improve accuracy for strongly perturbed bonds.
• iter_vel_correction=false: whether to use iterative velocity constraint solving. When false, velocity correction uses one LINCS projection as in GROMACS. When true, additional iterative velocity corrections are performed.
• gpu_block_size=128: the number of threads per block to use for GPU calculations.

LJDispersionCorrection(atoms, cutoff, σ_mixing=Molly.LorentzMixing(),
                       ϵ_mixing=Molly.GeometricMixing())

The long-range dispersion correction for the LennardJones interaction.

See Shirts et al. 2007. Approximately represents the part of the potential that the cutoff leaves out, assuming a uniform density. Should be used alongside the LennardJones pairwise interaction when the long-range correction to the potential energy is required. cutoff is the cutoff of the LennardJones interaction, which can be a DistanceCutoff, a ShiftedPotentialCutoff, a PolynomialCutoff or a distance, which is treated as a DistanceCutoff. For a DistanceCutoff the potential energy is defined as

\[E = \frac{8 \pi N^2}{V} \left( \frac{\left< \epsilon_{ij} \sigma_{ij}^{12} \right>}{9 r_c^9} - \frac{\left< \epsilon_{ij} \sigma_{ij}^{6} \right>}{3 r_c^3} \right)\]

where the averages run over all pairs of the $N$ atoms, including excluded pairs and the pair of each atom with itself, $V$ is the volume and $r_c$ is the cutoff distance. This is the correction of OpenMM. For a PolynomialCutoff the part of the potential removed by the switching function is also included, as in OpenMM. For a ShiftedPotentialCutoff the shift of the pairs within the cutoff is also included, as in GROMACS. GROMACS leaves excluded pairs out of the averages and only corrects the $r^{-6}$ term, so its correction differs slightly. The forces are zero, with an isotropic virial contribution.

The number of atoms is assumed not to change after setup, and changing the σ and ϵ of the atoms directly has no effect; inject_gradients recomputes the correction from injected values, including its gradients. Only compatible with 3D systems.

LJDispersionCorrectionλ(atoms, cutoff, λ_mix, σ_mix, ϵ_mix; scheduler)

The long-range dispersion correction for the LennardJones interaction, scaled by λ for alchemical transformations, the λ version of LJDispersionCorrection.

AbsoluteFESystem and RelativeFESystem convert each LJDispersionCorrection to this type. Approximately represents contributions from beyond the cutoff distance. For a DistanceCutoff the potential energy is defined as

\[E = \frac{8 \pi N^2}{V} \left( \frac{\left< \lambda_{ij} \epsilon_{ij} \sigma_{ij}^{12} \right>}{9 r_c^9} - \frac{\left< \lambda_{ij} \epsilon_{ij} \sigma_{ij}^{6} \right>}{3 r_c^3} \right)\]

and the other cutoffs are treated as in LJDispersionCorrection. The averages run over all pairs, $\lambda_{ij}$ is the steric coupling of the pair that scheduler gives from the λ of the atoms mixed with λ_mix and their alchemical roles, and $N$ is the effective number of atoms, since dual topology duplicates an atom as a real atom and a virtual one whose couplings sum to 1.

$\lambda_{ij}$ is applied exactly once per pair, so the correction is as linear in the coupling as the explicit pairwise term. Beyond the cutoff a soft core has converged to the plain Lennard-Jones potential, so no soft-core factor enters here.

The factors are computed once from the λ of atoms, so changing the λ of the atoms directly has no effect; Molly.set_lambda! recomputes them. The number of atoms is assumed not to change after setup, and inject_gradients recomputes the correction from injected σ and ϵ values, including its gradients. Only compatible with 3D systems.

LJZeroShortcut()

Skip a pairwise interaction when the σ or ϵ of either atom is zero.

Give this when setting up a pairwise interaction such as LennardJones.

LambdaScheduler(schedule=Molly.DefaultSchedule(); dual=true, LJindividual=false,
                LJspecial=false, Cindividual=false, Cspecial=false, intraLJ=false,
                Tscaled=true)

Turns global_λ into the couplings of the steric, electrostatic and bonded interactions of each alchemical role.

schedule selects how the couplings are staged along global_λ. Each schedule has an alias that is normally used instead: Molly.DefaultLambdaScheduler, Molly.LinearLambdaScheduler, Molly.GROMACSLambdaABFEScheduler, Molly.GROMACSLambdaRBFEScheduler, Molly.OpenFEScheduler, Molly.NAMDLambdaScheduler, Molly.QuartersLambdaScheduler and Molly.EleScaledLambdaScheduler, e.g. Molly.OpenFEScheduler(dual=false). Schedulers are passed to AbsoluteFESystem and RelativeFESystem.

Arguments

• dual=true: whether to use dual topology (energy scaling) or single topology (parameter scaling). false by default for Molly.OpenFEScheduler.
• LJindividual=false: in single topology, whether the Lennard-Jones parameters are scaled for each atom before mixing (true) or mixed for each state and then interpolated (false).
• LJspecial=false: the same choice for the Lennard-Jones 1-4 interactions.
• Cindividual=false: the same choice for the Coulomb interactions. true by default for Molly.OpenFEScheduler.
• Cspecial=false: the same choice for the Coulomb 1-4 interactions.
• intraLJ=false: whether the Lennard-Jones interactions between alchemical atoms are kept on (true) or scaled with the atoms (false).
• Tscaled=true: whether a torsion that spans core and unique atoms follows the core atom and is scaled (true) or follows the unique atom and stays on at every global_λ (false). false by default for Molly.OpenFEScheduler, which does not scale these torsions, so that Molly reproduces it exactly. Torsions of unique atoms only are never scaled either way.

Coulomb interactions between alchemical atoms are always scaled with the atoms, as the PME mesh scales the charge of each atom.

Langevin(; <keyword arguments>)

The Langevin integrator, based on the Langevin Middle Integrator in OpenMM.

See Zhang et al. 2019. In the notation of LangevinSplitting this is the BAOA splitting rather than the leapfrog Verlet integrator, but the velocities are stored as in leapfrog, offset by half a time step behind the positions.

Arguments

• dt::S: the time step of the simulation.
• temperature::K: the equilibrium temperature of the simulation.
• friction::F: the friction coefficient of the simulation, e.g. 1.0u"ps^-1".
• coupling::C=nothing: the coupling which applies during the simulation.
• remove_CM_motion=1: remove the center of mass motion every this number of steps, set to false or 0 to not remove center of mass motion.

LangevinSplitting(; <keyword arguments>)

The Langevin simulator using a general splitting scheme.

This consists of a succession of A, B and O steps, corresponding respectively to updates in position, velocity for the potential part, and velocity for the thermal fluctuation-dissipation part. The Langevin and VelocityVerlet simulators without coupling correspond to the BAOA and BAB schemes respectively. For more information on the sampling properties of splitting schemes, see Fass et al. 2018.

Not compatible with constraints, will print a warning and continue without applying constraints. Not compatible with virial calculation.

Whenever A steps have moved the atoms since the neighbors were last found, the neighbors are found on the schedule of the neighbor finder before the forces are next calculated and at the end of the step.

Arguments

• dt::S: the time step of the simulation.
• temperature::K: the equilibrium temperature of the simulation.
• friction::F: the friction coefficient. If units are used, it should have a dimensionality of mass per time. Note this is a different friction from Langevin, which has dimensionality of inverse time.
• splitting::W: the splitting specifier. Should be a string consisting of the characters A, B and O. Strings with no Os reduce to deterministic symplectic schemes.
• remove_CM_motion=1: remove the center of mass motion every this number of steps, set to false or 0 to not remove center of mass motion.

LennardJones(; cutoff=NoCutoff(), use_neighbors=true, shortcut=Molly.LJZeroShortcut(),
             σ_mixing=Molly.LorentzMixing(), ϵ_mixing=Molly.GeometricMixing(),
             weight_special=1)

The Lennard-Jones 6-12 interaction between two atoms.

The potential energy is defined as

\[V(r_{ij}) = 4\varepsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]\]

and the force on each atom by

\[\begin{aligned} \vec{F}_i &= 24\varepsilon_{ij} \left(2\frac{\sigma_{ij}^{12}}{r_{ij}^{13}} - \frac{\sigma_{ij}^6}{r_{ij}^{7}}\right) \frac{\vec{r}_{ij}}{r_{ij}} \\ &= \frac{24\varepsilon_{ij}}{r_{ij}^2} \left[2\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} -\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right] \vec{r}_{ij} \end{aligned}\]

Should be used alongside the LJDispersionCorrection general interaction when the long-range correction to the potential energy is required.

LennardJones14(σ14_mixed, ϵ14_mixed, weight_14)

The Lennard-Jones 6-12 interaction of a 1-4 pair with its own σ and ϵ, as a specific interaction between the two atoms.

Force fields such as CHARMM give the 1-4 pairs σ and ϵ of their own instead of scaling the parameters of the pair, so these pairs are kept in a list of this interaction and left out of the pairwise Lennard-Jones, which has weight_special = 0. σ14_mixed and ϵ14_mixed are already mixed for the pair and weight_14 is the 1-4 scaling factor of the force field.

The potential energy is defined as

\[V(r_{ij}) = 4 w_{14} \varepsilon^{14}_{ij} \left[\left(\frac{\sigma^{14}_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma^{14}_{ij}}{r_{ij}}\right)^{6}\right]\]

with $w_{14}$ the 1-4 weight.

The λ versions for alchemical systems are Molly.LennardJones14Scaled, LennardJones14SoftCoreBeutler and LennardJones14SoftCoreGapsys.

LennardJones14Scaled(; σ14_mixed, ϵ14_mixed, weight_14, λ_mixing=Molly.MinimumMixing(),
                     scheduler=Molly.DefaultLambdaScheduler())

The 1-4 Lennard-Jones interaction of a pair of atoms scaled by λ without a soft core, the λ version of LennardJones14, which AbsoluteFESystem and RelativeFESystem convert to this type with lj_softcore=:scaled.

Force fields that give the 1-4 pairs their own σ and ϵ, such as CHARMM, keep them in a LennardJones14 list rather than in the pairwise Lennard-Jones. The energy and the force are those of Molly.LennardJonesScaled with the σ and ϵ of the pair, weighted by weight_14, and the pair is coupled like any other Lennard-Jones pair of the alchemical system.

LennardJones14SoftCoreBeutler(; σ14_mixed, ϵ14_mixed, weight_14, α=0.5,
                              λ_mixing=Molly.MinimumMixing(),
                              scheduler=Molly.DefaultLambdaScheduler())

The 1-4 Lennard-Jones interaction of a pair of atoms with a Beutler soft core, the λ version of LennardJones14, which AbsoluteFESystem and RelativeFESystem convert to this type with lj_softcore=:beutler.

Force fields that give the 1-4 pairs their own σ and ϵ, such as CHARMM, keep them in a LennardJones14 list rather than in the pairwise Lennard-Jones. The energy and the force are those of LennardJonesSoftCoreBeutler with the σ and ϵ of the pair, weighted by weight_14, and the pair is coupled like any other Lennard-Jones pair of the alchemical system.

LennardJones14SoftCoreGapsys(; σ14_mixed, ϵ14_mixed, weight_14, α=0.85,
                             λ_mixing=Molly.MinimumMixing(),
                             scheduler=Molly.DefaultLambdaScheduler())

The 1-4 Lennard-Jones interaction of a pair of atoms with a Gapsys soft core, the λ version of LennardJones14, which AbsoluteFESystem and RelativeFESystem convert to this type with lj_softcore=:gapsys.

Force fields that give the 1-4 pairs their own σ and ϵ, such as CHARMM, keep them in a LennardJones14 list rather than in the pairwise Lennard-Jones. The energy and the force are those of LennardJonesSoftCoreGapsys with the σ and ϵ of the pair, weighted by weight_14, and the pair is coupled like any other Lennard-Jones pair of the alchemical system.

LennardJonesScaled(; cutoff=NoCutoff(), use_neighbors=true,
                   shortcut=Molly.LJZeroShortcut(), σ_mixing=Molly.LorentzMixing(),
                   ϵ_mixing=Molly.GeometricMixing(), λ_mixing=Molly.MinimumMixing(),
                   scheduler=Molly.DefaultLambdaScheduler(), weight_special=1)

The Lennard-Jones 6-12 interaction between two atoms scaled by λ without a soft core, the λ version of LennardJones.

AbsoluteFESystem and RelativeFESystem convert each LennardJones to this type with lj_softcore=:scaled. The potential energy is defined as

\[V(r_{ij}) = 4\lambda\varepsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right]\]

and the force on each atom by

\[\vec{F}_i = \frac{24\lambda\varepsilon_{ij}}{r_{ij}^2} \left[2\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} -\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right] \vec{r}_{ij}\]

where $\lambda$ is the steric coupling of the pair, which scheduler, a Molly.LambdaScheduler, gives from the λ of the atoms mixed with λ_mixing and their alchemical roles.

If $\lambda$ is 1 this gives the standard LennardJones potential, and if it is 0 the interaction is turned off. As there is no soft core the energy diverges when an atom that is being inserted or deleted overlaps another one, so this is for a transformation where that does not happen, for example one where only parameters change. In single topology $\lambda$ is one and $\sigma$ and $\varepsilon$ are interpolated between the end states instead.

LennardJonesSoftCoreBeutler(; cutoff=NoCutoff(), α=0.3, use_neighbors=true,
                            shortcut=Molly.LJZeroShortcut(),
                            σ_mixing=Molly.LorentzMixing(),
                            ϵ_mixing=Molly.GeometricMixing(),
                            λ_mixing=Molly.MinimumMixing(),
                            scheduler=Molly.DefaultLambdaScheduler(), weight_special=1)

The Lennard-Jones 6-12 interaction between two atoms with a soft core, used for the appearing and disappearing of atoms.

See Beutler et al. 1994. AbsoluteFESystem and RelativeFESystem convert each LennardJones to this type with lj_softcore=:beutler. The potential energy is defined as

\[V(r_{ij}) = \lambda \left(\frac{C^{(12)}}{r_{LJ}^{12}} - \frac{C^{(6)}}{r_{LJ}^{6}}\right)\]

and the force on each atom by

\[\vec{F}_i = \lambda \left(\left(\frac{12C^{(12)}}{r_{LJ}^{13}} - \frac{6C^{(6)}}{r_{LJ}^7}\right)\left(\frac{r_{ij}}{r_{LJ}}\right)^5\right) \frac{\vec{r}_{ij}}{r_{ij}}\]

where

\[r_{LJ} = \left(\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}}+r_{ij}^6\right)^{1/6}, \quad C^{(12)} = 4\varepsilon_{ij}\sigma_{ij}^{12}, \quad C^{(6)} = 4\varepsilon_{ij}\sigma_{ij}^{6}\]

and $\sigma_{ij}$ and $\varepsilon_{ij}$ are mixed with σ_mixing and ϵ_mixing. $\lambda$ is the steric coupling of the pair, which scheduler, a Molly.LambdaScheduler, gives from the λ of the atoms mixed with λ_mixing and their alchemical roles. If $\lambda$ is 1 this gives the standard LennardJones potential, and if it is 0 the interaction is turned off. $\alpha$ determines the strength of softening the function.

LennardJonesSoftCoreGapsys(; cutoff=NoCutoff(), α=0.85, use_neighbors=true,
                           shortcut=Molly.LJZeroShortcut(),
                           σ_mixing=Molly.LorentzMixing(),
                           ϵ_mixing=Molly.GeometricMixing(),
                           λ_mixing=Molly.MinimumMixing(),
                           scheduler=Molly.DefaultLambdaScheduler(), weight_special=1)

The Lennard-Jones 6-12 interaction between two atoms with a soft core potential, used for the appearing and disappearing of atoms.

See Gapsys et al. 2012. AbsoluteFESystem and RelativeFESystem convert each LennardJones to this type with lj_softcore=:gapsys. The potential energy is defined as

\[V(r_{ij}) = \left\{ \begin{array}{cl} \lambda \left( \frac{C^{(12)}}{r_{ij}^{12}} - \frac{C^{(6)}}{r_{ij}^{6}} \right), & \text{if} & r_{ij} \ge r_{LJ} \\ \lambda \left( (\frac{78C^{(12)}}{r_{LJ}^{14}}-\frac{21C^{(6)}}{r_{LJ}^{8}})r_{ij}^2 - (\frac{168C^{(12)}}{r_{LJ}^{13}}-\frac{48C^{(6)}}{r_{LJ}^{7}})r_{ij} + \frac{91C^{(12)}}{r_{LJ}^{12}}-\frac{28C^{(6)}}{r_{LJ}^{6}} \right), & \text{if} & r_{ij} \lt r_{LJ} \\ \end{array} \right.\]

and the force on each atom by

\[\vec{F}_i = \left\{ \begin{array}{cl} \lambda \left( \frac{12C^{(12)}}{r_{ij}^{13}} - \frac{6C^{(6)}}{r_{ij}^{7}} \right)\frac{\vec{r}_{ij}}{r_{ij}}, & \text{if} & r_{ij} \ge r_{LJ} \\ \lambda \left( (\frac{-156C^{(12)}}{r_{LJ}^{14}}+\frac{42C^{(6)}}{r_{LJ}^{8}})r_{ij} + (\frac{168C^{(12)}}{r_{LJ}^{13}}-\frac{48C^{(6)}}{r_{LJ}^{7}}) \right)\frac{\vec{r}_{ij}}{r_{ij}}, & \text{if} & r_{ij} \lt r_{LJ} \\ \end{array} \right.\]

where

\[r_{LJ} = \alpha \left( \frac{26C^{(12)}(1-\lambda)}{7C^{(6)}} \right)^{\frac{1}{6}}, \quad C^{(12)} = 4\varepsilon_{ij}\sigma_{ij}^{12}, \quad C^{(6)} = 4\varepsilon_{ij}\sigma_{ij}^{6}\]

and $\sigma_{ij}$ and $\varepsilon_{ij}$ are mixed with σ_mixing and ϵ_mixing. $\lambda$ is the steric coupling of the pair, which scheduler, a Molly.LambdaScheduler, gives from the λ of the atoms mixed with λ_mixing and their alchemical roles. If $\lambda$ is 1 this gives the standard LennardJones potential, and if it is 0 the interaction is turned off. $\alpha$ determines the strength of softening the function.

LinearBias(k, cv_target)

A linear bias on a collective variable (CV) towards a target value.

The potential energy is defined as

\[V(\boldsymbol{s}) = k |\boldsymbol{s} - \boldsymbol{s}_t|\]

where $s$ and $s_t$ are the system and target CV values respectively.

Arguments

• k: The energy constant for the bias. Must be compliant with the System energy units.
• cv_target: The target value of the collective variable.

LinearLambdaScheduler(; dual=true, LJindividual=false, LJspecial=false,
                      Cindividual=false, Cspecial=false, intraLJ=false,
                      Tscaled=true)

Lambda scheduler that scales the electrostatics and sterics of all alchemical atoms linearly and at the same time, from global_λ = 0 to 1.

A Molly.LambdaScheduler, see there for the keyword arguments.

LinearSchedule()

Schedule of a Molly.LinearLambdaScheduler, passed to Molly.LambdaScheduler.

ListHills(cvs, height, sigma, centers=typeof(float.(sigma))[],
          heights=fill(height, length(centers)))

Metadynamics memory for the CVs in the tuple cvs, storing every deposited hill explicitly, summed at each evaluation:

\[V(s) = \sum_{t' < t} h_{t'} \exp\left(-\frac{|s - s(t')|^2}{2\sigma^2}\right)\]

O(n_hills) per evaluation; see GridHills for an O(1) alternative.

sigma and the elements of centers are scalars for a single CV, or tuples with one entry per CV for several. Hills on periodic CVs such as torsions wrap around (see Molly.cv_period).

Arguments

• cvs: Tuple of CV descriptors, e.g. CalcDist; a single CV may be passed without the tuple.
• height: Default hill height, used unless Molly.add_hill! is given an explicit height.
• sigma: Gaussian width (standard deviation), in CV units.
• centers: CV values where hills have already been deposited. Defaults to an empty vector with the float element type of sigma; its element type sets the precision and units of every later deposit.
• heights: Deposited height per entry in centers; defaults to height.

LorentzMixing()

Mix parameters according to their arithmetic mean.

For example:

\[\sigma_{ij} = \frac{\sigma_i + \sigma_j}{2}\]

Give this when setting up a pairwise interaction such as LennardJones.

The reduced potentials and frame layout used as input to MBAR.

Returned by Molly.assemble_mbar_inputs, which describes the fields, and given to Molly.mbar_weights.

MBARResult(u, win_of, N_counts, states; kwargs...)
MBARResult(u, win_of, N_counts, kBT; names=nothing, kwargs...)
MBARResult(input::MBARInput, states; kwargs...)

Relative free energies of K thermodynamic states from MBAR, with their asymptotic covariance, see Shirts and Chodera 2008 Eq D6.

It is a leg source for ThermodynamicCycle: Molly.free_energy_difference and Molly.free_energy_difference_uncertainty return values in energy units.

Arguments

• u::AbstractMatrix - N×K reduced potentials, as for Molly.iterate_mbar. Per-frame shifts, as from Molly.assemble_mbar_inputs, do not change the result.
• win_of::AbstractVector - length-N vector giving the generating state index per frame.
• N_counts::AbstractVector - length-K sample counts per state.
• states::AbstractVector{<:ThermoState} - the K states, which give the thermal energy and the name of each state.
• kBT::AbstractVector - the thermal energy of each state, in energy units, when the states are not given.
• input::MBARInput - the output of Molly.assemble_mbar_inputs for states.

Keyword arguments

• names=nothing - state names, used to refer to states by name in a ThermodynamicCycleLeg.
• kwargs... - other keyword arguments, such as method, rtol, atol and max_iter, are passed to Molly.iterate_mbar.

Fields

• f::Vector{Float64} - relative free energies in units of each state's $k_B T$, with f[1] = 0.
• covariance::Matrix{Float64} - asymptotic covariance of f, with the gauge direction projected out; the variance of f[j] - f[i] is Θ[i,i] + Θ[j,j] - 2Θ[i,j].
• kBT::Vector - thermal energy of each state.
• names - state names, or nothing.

MTSIntegrator(; <keyword arguments>)

The rRESPA multiple time step integration algorithm.

See Tuckerman et al. 1992. This integrator allows different interactions to be evaluated at different frequencies, so more expensive interactions can be calculated less frequently.

The fraction arguments correspond to the number of times each interaction is applied per outer step. For example, if your system had two pairwise interactions and three specific interactions you could pass pi_fractions=(1, 1) and si_fractions=(2, 4). With a dt of 4 fs this means that the pairwise interactions are evaluated once per outer time step (every 4 fs), the first specific interaction twice (every 2 fs), and the second specific interaction four times (every 1 fs). All fractions should be multiples of previous fractions, and at least one fraction must be 1. When simulating, n_steps is the number of outer steps. Coupling, removing the center of mass motion and running loggers applies per outer time step.

With inner_step_neighbors=false and a neighbor finder with auto_update=true, the neighbors are checked before every force calculation that uses them at coordinates that have changed since they were last checked, and rebuilt if needed. With a neighbor finder with auto_update=false, the neighbors are found once per outer step on the schedule of the neighbor finder, after the coordinates are updated and before the interactions evaluated once per outer step. Interactions evaluated more often then use these neighbors throughout the next outer step, so if any pairwise interactions are evaluated more than once per outer step the neighbor finder should use a skin sufficient to cover atomic motion over an outer step.

Arguments

• dt::T: the outer time step of the simulation.
• pi_fractions: the number of times each pairwise interaction is evaluated per outer timestep, should be a tuple of integers with the same length as sys.pairwise_inters.
• si_fractions: the number of times each specific interaction is evaluated per outer timestep, should be a tuple of integers with the same length as sys.specific_inter_lists.
• gi_fractions: the number of times each general interaction is evaluated per outer timestep, should be a tuple of integers with the same length as sys.general_inters.
• coupling::C=nothing: the coupling which applies each outer step during the simulation.
• remove_CM_motion=1: remove the center of mass motion every this number of outer steps, set to false or 0 to not remove center of mass motion.
• inner_step_neighbors=false: whether to force recomputation of the neighbors at every inner step, useful when force calculation is slow compared to neighbor finding. If used, the auto_update and n_steps arguments of the neighbor finder are ignored.

MTSLangevinIntegrator(; <keyword arguments>)

The BAOAB-RESPA multiple time step algorithm for constant temperature dynamics.

See Tuckerman et al. 1992 and Lagardère et al. 2019. This integrator allows different interactions to be evaluated at different frequencies, so more expensive interactions can be calculated less frequently.

The fraction arguments correspond to the number of times each interaction is applied per outer step. For example, if your system had two pairwise interactions and three specific interactions you could pass pi_fractions=(1, 1) and si_fractions=(2, 4). With a dt of 4 fs this means that the pairwise interactions are evaluated once per outer time step (every 4 fs), the first specific interaction twice (every 2 fs), and the second specific interaction four times (every 1 fs). All fractions should be multiples of previous fractions, and at least one fraction must be 1. When simulating, n_steps is the number of outer steps. Coupling, removing the center of mass motion and running loggers applies per outer time step.

With inner_step_neighbors=false and a neighbor finder with auto_update=true, the neighbors are checked before every force calculation that uses them at coordinates that have changed since they were last checked, and rebuilt if needed. With a neighbor finder with auto_update=false, the neighbors are found once per outer step on the schedule of the neighbor finder, after the coordinates are updated and before the interactions evaluated once per outer step. Interactions evaluated more often then use these neighbors throughout the next outer step, so if any pairwise interactions are evaluated more than once per outer step the neighbor finder should use a skin sufficient to cover atomic motion over an outer step.

Arguments

• dt::T: the outer time step of the simulation.
• temperature::K: the equilibrium temperature of the simulation.
• friction::F: the friction coefficient of the simulation, e.g. 1.0u"ps^-1".
• pi_fractions: the number of times each pairwise interaction is evaluated per outer timestep, should be a tuple of integers with the same length as sys.pairwise_inters.
• si_fractions: the number of times each specific interaction is evaluated per outer timestep, should be a tuple of integers with the same length as sys.specific_inter_lists.
• gi_fractions: the number of times each general interaction is evaluated per outer timestep, should be a tuple of integers with the same length as sys.general_inters.
• coupling::C=nothing: the coupling which applies each outer step during the simulation.
• remove_CM_motion=1: remove the center of mass motion every this number of outer steps, set to false or 0 to not remove center of mass motion.
• inner_step_neighbors=false: whether to force recomputation of the neighbors at every inner step, useful when force calculation is slow compared to neighbor finding. If used, the auto_update and n_steps arguments of the neighbor finder are ignored.

Metadynamics(memory; n_steps=1, tempering=Molly.NoTempering(), logger=nothing)

A history-dependent bias potential for Metadynamics: a sum of Gaussians deposited at previously visited collective variable (CV) values.

\[V(\boldsymbol{s}) = \sum_{t' < t} h_{t'} \exp\left(-\frac{|\boldsymbol{s} - \boldsymbol{s}(t')|^2}{2\sigma^2}\right)\]

The memory argument (a Molly.AbstractMetadynamicsMemory) holds the CVs and controls how hills are stored: ListHills (exact, O(n_hills)) or GridHills (O(1), grid-based).

With the CVs in the memory, e.g. Metadynamics(ListHills((cv,), height, sigma)), the bias is an AtomsCalculators.jl calculator usable directly as a general_inters entry. Each CV is evaluated on the coordinates its own correction asks for, and hills on periodic CVs such as torsions wrap around (see Molly.cv_period). Hills are deposited into the memory passed in, which is also bias.memory, so biases built from one memory are walkers sharing their hills. Deposits into a memory aren't thread-safe, so walkers sharing a memory must run one after the other, not in parallel.

Hills are deposited when a simulator visits a state, on the first visit in each window of n_steps steps – no external logger is needed. Molecular dynamics integrators visit a state in their force calculation (visit_state=true in forces). MetropolisMonteCarlo, ReplicaExchangeMD, AWHSimulation and TSSSimulation aren't supported yet. The first hill is on the start configuration, at step 0, or at step 1 for integrators such as Langevin. A multiple time step integrator deposits once per outer step. Minimizers, loggers, couplers and other force or energy calls don't deposit. Chunked simulate! calls deposit as one continuous run, with or without init_step, as long as a bias uses one of the two throughout: mixing them, or visiting with a step number between runs, shifts the schedule. The number of hills deposited so far is bias.n_deposits and the step of the last deposit window is bias.last_deposit_step (outside a GridHills grid that window adds no hill).

Molly.add_hill! deposits a hill manually, on every n_steps-th call. Force or energy calls outside a simulator deposit only when given visit_state=true: with a step number, as forces and potential_energy pass (0 by default), they follow the step schedule; without one, as in AtomsCalculators.forces(sys, bias; visit_state=true), they count towards the same call counter as Molly.add_hill!.

tempering (a Molly.AbstractTempering) scales each deposited hill's height; the default Molly.NoTempering always deposits the full base height.

Arguments

• memory::AbstractMetadynamicsMemory: The CVs, and the storage and evaluation strategy for deposited hills.
• n_steps::Integer=1: Number of simulation steps between deposits into memory, starting with the start configuration. Molly.add_hill! calls and visiting force calls without a step number use a separate call counter, depositing on every n_steps-th such call without affecting the simulation schedule.
• tempering::AbstractTempering=NoTempering(): Scales the height of each deposited hill.
• logger=nothing: A MetadynamicsLogger to record each deposit for reweighting, available as bias.logger.

MetadynamicsLogger(; temperature=nothing, cv_points=nothing, filepath=nothing,
                   write_interval=100, overwrite=false)

Records each deposit of a Metadynamics bias for reweighting: the step, the CV value s, the bias V(s, t) it was sampled under and the offset c(t) (Tiwary & Parrinello, 2015). Pass it as Metadynamics(memory; logger=MetadynamicsLogger()); the records are in bias.logger, with values(bias.logger) the CV values.

temperature defaults to that of WellTempered and must be given otherwise. kT, bias.logger.kT, is fixed at the first deposit with the Boltzmann constant of its System; a later deposit with a different kT throws. c(t) integrates the bias over the grid of a GridHills, and over cv_points, a uniform grid of CV values, for other memories. See Molly.metadynamics_weights.

c(t) costs O(grid points or cv_points) per deposit, so 2D/3D grids want a larger n_steps.

With a filepath, the records are also appended to that text file, with the columns step, CV values, bias, offset, hill height and log weight, in batches of write_interval records. Call flush(bias.logger) after the run to write the last batch. An existing file is appended to as a new run with its own header, unless overwrite=true, which replaces it. Use one file per logger, also for walkers: a flush that finds its file changed by another writer throws. Molly.read_metadynamics_logger reads the file back.

MetropolisMonteCarlo(; <keyword arguments>)

A Monte Carlo simulator that uses the Metropolis algorithm to sample the configuration space.

A trial move with potential energy change $\Delta E$ is accepted with probability $\min(1, \exp(-\Delta E / k_B T))$. MonteCarloLogger can be used to record the acceptances. Not compatible with constraints, will print a warning and continue without applying constraints. Throws an error for history-dependent general interactions such as Metadynamics.

Arguments

• temperature::T: the temperature of the system.
• trial_moves::M: a function that performs the trial moves, called as trial_moves(sys; trial_args...), for example Molly.random_uniform_translation! or Molly.random_normal_translation!.
• trial_args::Dict=Dict(): a dictionary of keyword arguments to be passed to the trial move function.

Mie(; m, n, cutoff=NoCutoff(), use_neighbors=true, shortcut=Molly.LJZeroShortcut(),
    σ_mixing=Molly.LorentzMixing(), ϵ_mixing=Molly.GeometricMixing(), weight_special=1)

The Mie generalized interaction between two atoms.

When m equals 6 and n equals 12 this is equivalent to the Lennard-Jones interaction. The potential energy is defined as

\[V(r_{ij}) = C \varepsilon_{ij} \left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^n - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^m\right]\]

where

\[C = \frac{n}{n - m} \left( \frac{n}{m} \right) ^\frac{m}{n - m}\]

MinimumMixing()

Mix λ parameters by taking the minimum over the atoms in an interaction.

For an interaction between atoms $i$, $j$, $k$, ...:

\[\lambda_{ijk \dots} = \min(\lambda_i, \lambda_j, \lambda_k, \dots)\]

This works for pairwise and specific interactions with any number of atoms. Give this when setting up an alchemical interaction, where it is the default.

MissingResidueTemplateError(msg)

An error when a residue cannot be matched to a residue template.

MixingException(mixing, exceptions)

Provide mixing exceptions (NBFix) for specific pairs of atom types.

mixing is a mixing rule such as Molly.LorentzMixing, used for pairs that are not exceptions, and exceptions is a Molly.ExceptionList keyed by atom type. Give this when setting up a pairwise interaction such as LennardJones.

MolecularForceField(ff_files...; units=true, custom_residue_templates=nothing,
                    custom_renaming_scheme=nothing, float_type=Float64,
                    strictness=:warn)

A molecular force field.

Read one or more OpenMM force field XML files by passing them to the constructor. See the OpenMM documentation for how these files are formatted.

In order to assign force field parameters to the atoms in the simulation, the residues determined from a structure file are matched to templates provided by the force field file, as well as a template dictionary in XML format, which defines the standard topology (bonds) of the residues to be found in the simulation. At the moment, Molly provides a dictionary for all standard amino acids, nucleic acids and water for this purpose. If the system to be simulated contains other molecules, their template topologies must be defined either through CONECT records in the PDB file or by providing an extra custom template file to the custom_residue_templates keyword argument.

float_type should generally be Float64 since the float type of a System is determined later when creating the System. Behavior with unsupported files is determined by the strictness keyword argument. This can be :warn to emit warnings, :nowarn to suppress warnings or :error to error. Failures when reading in force field XML files throw a Molly.ForceFieldXMLError.

MolecularTopology(bond_is, bond_js, n_atoms)
MolecularTopology(atom_molecule_inds, molecule_atom_counts, bonded_atoms=[])

Topology information for a system.

Stores the index of the molecule each atom belongs to, the number of atoms in each molecule and the list of bonded atom pairs.

MollyCalculator(; <keyword arguments>)

A calculator for use with the AtomsCalculators.jl interface.

neighbors can optionally be given as a keyword argument when calling the calculation functions to save on computation when the neighbors are the same for multiple calls. In a similar way, n_threads can be given to determine the number of threads to use when running the calculation function. visit_state=true lets history-dependent general interactions such as Metadynamics update (see forces); pass the current step_n with it so they follow the step schedule. Note that this calculator is designed for using Molly in other contexts; if you want to use another calculator in Molly it can be given as general_inters when creating a System.

Not compatible with virial calculation. Not compatible with using atom properties such as σ and ϵ.

Arguments

• pairwise_inters::PI=(): the pairwise interactions in the system, i.e. interactions between all or most atom pairs such as electrostatics. Should be a Tuple or NamedTuple of PairwiseInteractions.
• specific_inter_lists::SI=(): the specific interactions in the system, i.e. interactions between specific atoms such as bonds or angles. Should be a Tuple or NamedTuple.
• general_inters::GI=(): the general interactions in the system, i.e. interactions involving all atoms such as implicit solvent. Each should implement the AtomsCalculators.jl interface. Should be a Tuple or NamedTuple.
• neighbor_finder::NF=nothing: the neighbor finder used to find close atoms and save on computation. nothing means that no neighbor finder is used and all pairs of atoms are considered, in which case use_neighbors should return false for the pairwise interactions.
• force_units::F=u"kJ * mol^-1 * nm^-1": the units of force of the system. Should be set to NoUnits if units are not being used.
• energy_units::E=u"kJ * mol^-1": the units of energy of the system. Should be set to NoUnits if units are not being used.
• k::K=Unitful.k or Unitful.k * Unitful.Na: the Boltzmann constant, which may be modified in some simulations. k is chosen based on the energy_units given.
• dims::Integer=3: the number of dimensions in the system.

MonteCarloBarostat(pressure, temperature, boundary; coupling_type=:isotropic,
                   n_steps=30, n_iterations=1,
                   scale_factor=0.01, scale_increment=1.1, max_volume_frac=0.3,
                   trial_find_neighbors=false) <: AbstractBarostat

The Monte Carlo barostat for controlling pressure.

See Chow and Ferguson 1995, Åqvist et al. 2004 and the OpenMM source code. At regular intervals a Monte Carlo step is attempted by scaling the coordinates and the bounding box by a randomly chosen amount. The step is accepted or rejected based on

\[\Delta G = \Delta E + \Delta W - N k_B T \ln \left( \frac{V + \Delta V}{V} \right)\]

where ΔE is the change in potential energy, ΔV is the change in volume, N is the number of molecules in the system, T is the equilibrium temperature and V is the system volume. ΔW is the work done by scaling the simulation box and its specific form changes depending on the type of scaling applied. In general and in the absence of shear stress:

\[\Delta W = (V + \Delta V) \cdot \sum w_i \cdot P_{i,i} \cdot \ln \left ( {\frac{V + \Delta V}{V}} \right )\]

where w_i is the proportional scaling along a specific box axis.

If ΔG ≤ 0 the step is always accepted, if ΔG > 0 the step is accepted with probability exp(-ΔG/kT).

The scale factor is modified over time to maintain an acceptance rate of around half. If the topology of the system is set then molecules are moved as a unit so properties such as bond lengths do not change.

The barostat assumes that the simulation is being run at a constant temperature but does not actively control the temperature. It should be used alongside a temperature coupling method such as the Langevin simulator or AndersenThermostat coupling. With a neighbor finder with auto_update=true, the default, the neighbor list is checked at the trial coordinates and rebuilt if the trial move could have invalidated it. Otherwise by default the neighbor list is not updated when calculating the energy of trial moves, which assumes that the volume changes are small. Set trial_find_neighbors to true to rebuild it for every trial move. Note that the barostat can change the bounding box of the system. Does not work with shear stresses, the anisotropic variant only applies independent linear scaling of the box vectors. If shear deformation is required the BerendsenBarostat or, preferably, the CRescaleBarostat should be used instead.

Not compatible with gradient calculation using Enzyme.

MonteCarloLogger(T=Float64)

A logger that records acceptances in a Monte Carlo simulation.

The logged quantities include the number of trial moves (n_trials), the number of successful acceptances (n_accept), an array named energy_rates which stores the value of $\frac{E}{k_B T}$ i.e. the argument of the Boltzmann factor for the states, and a BitVector named state_changed that stores whether a new state was accepted for the logged step.

MorseBond(; D, a, r0)

A Morse potential bond between two atoms.

The potential energy is defined as

\[V(r) = D(1 - e^{-a(r - r_0)})^2\]

MorseBondλ(; D, a, r0, λ_mixing=Molly.MinimumMixing(),
           scheduler=Molly.DefaultLambdaScheduler())

A Morse potential bond between two atoms in an alchemical system, the λ version of MorseBond.

AbsoluteFESystem and RelativeFESystem convert each MorseBond to this type. The potential energy is defined as

\[V(r) = w D(1 - e^{-a(r - r_0)})^2\]

where the weight $w$ depends on the alchemical roles of the atoms and on $\lambda$, the λ of the atoms mixed with λ_mixing. $w = 1$ if the bond involves an atom unique to one state or only environment atoms, otherwise $w = 1 - \lambda$ if it involves a core atom with the parameters of state A and $w = \lambda$ if it involves a core atom with the parameters of state B.

scheduler is a Molly.LambdaScheduler, which has to use dual topology. Single topology is not supported for this interaction, so RelativeFESystem throws an error if the scheduler has dual=false.

MullerBrown(; A=SVector(-200.0, -100.0, -170.0, 15.0)u"kJ * mol^-1",
            a=SVector(-1.0, -1.0, -6.5, 0.7)u"nm^-2",
            b=SVector(0.0, 0.0, 11.0, 0.6)u"nm^-2",
            c=SVector(-10.0, -10.0, -6.5, 0.7)u"nm^-2",
            x0=SVector(1.0, 0.0, -0.5, -1.0)u"nm", y0=SVector(0.0, 0.5, 1.5, 1.0)u"nm",
            force_units=u"kJ * mol^-1 * nm^-1", energy_units=u"kJ * mol^-1")

The Müller-Brown potential energy surface implemented as an AtomsCalculators.jl calculator.

See Müller and Brown 1979. The potential energy is defined as

\[V(x,y) = \sum_{k=1}^{4} A_k \exp[a_k(x-x_k^0)^2 + b_k(x-x_k^0)(y-y_k^0) + c_k(y-y_k^0)^2]\]

where A, a, b, c, x0, y0 are 4-element SVectors with standard defaults.

This potential is only compatible with 2D systems. It is often used for testing algorithms that find transition states or explore minimum energy pathways. There are 3 minima and 2 saddle points with the default parameters.

Not compatible with virial calculation.

NAMDLambdaScheduler(; dual=true, LJindividual=false, LJspecial=false,
                    Cindividual=false, Cspecial=false, intraLJ=false,
                    Tscaled=true)

Lambda scheduler with overlapping electrostatic and steric stages, similar to the separate electrostatic and van der Waals windows in NAMD.

Atoms that only exist in the first state lose their charges between global_λ = 0 and 0.5 and their Lennard-Jones interactions between 1/3 and 1. Atoms that only exist in the second state gain their Lennard-Jones interactions between 0 and 2/3 and their charges between 0.5 and 1. Core atoms are interpolated linearly.

A Molly.LambdaScheduler, see there for the keyword arguments.

NAMDSchedule()

Schedule of a Molly.NAMDLambdaScheduler, passed to Molly.LambdaScheduler.

NaNSimulationError(msg)

An error when a NaN is encountered during a simulation or setup.

NeighborList(n, list)
NeighborList()

Structure to contain neighbor lists, as returned by find_neighbors for most neighbor finders.

The first n entries of list are the (i, j, special) pairs, see neighbor_pairs.

NoCutoff()

Placeholder cutoff that does not alter the potential or force.

NoTempering()

Default tempering: every hill keeps its full base height (standard, constant-height Metadynamics).

NoseHoover(; <keyword arguments>)

The Nosé-Hoover chain integrator, an NVT simulator that couples the system to a chain of thermostats.

This follows the NoseHooverIntegrator in OpenMM, which uses the "middle" leapfrog scheme of Zhang et al. 2019 with the Nosé-Hoover chain of Martyna et al. 1992 propagated using the multiple time step and Yoshida-Suzuki factorization of Martyna et al. 1996. As with Langevin, the velocities are stored as in leapfrog, offset by half a time step behind the positions.

The equations of motion for a chain of $M$ thermostats are

\[\begin{aligned} \dot{\bf{r}}_i &= \bf{v}_i \\ \dot{\bf{v}}_i &= \frac{\bf{F}_i}{m_i} - v_{\xi_1} \bf{v}_i \\ \dot{\xi}_j &= v_{\xi_j} \\ \dot{v}_{\xi_j} &= G_j - v_{\xi_j} v_{\xi_{j+1}} \quad (j < M) \\ \dot{v}_{\xi_M} &= G_M \end{aligned}\]

where the forces on the thermostats are

\[G_1 = \frac{1}{Q_1} \left( \sum_i m_i \bf{v}_i^2 - N_\mathrm{df} k T \right), \qquad G_j = \frac{1}{Q_j} \left( Q_{j-1} v_{\xi_{j-1}}^2 - k T \right) \quad (j > 1)\]

and the thermostat masses are $Q_1 = N_\mathrm{df} k T \tau^2$ and $Q_j = k T \tau^2$ for $j > 1$. $T$ is the target temperature, $\tau$ is damping and $N_\mathrm{df}$ is the number of degrees of freedom of the system, which accounts for constraints and the removal of center of mass motion. $\tau$ is the inverse of the collision frequency used by OpenMM and $\tau_T / 2\pi$ for the oscillation period $\tau_T$ used by GROMACS.

Each time step $\Delta t$ consists of

\[\begin{aligned} \bf{v}_i &\leftarrow \bf{v}_i + \frac{\bf{F}_i}{m_i} \Delta t \\ \bf{r}_i &\leftarrow \bf{r}_i + \bf{v}_i \frac{\Delta t}{2} \\ \bf{v}_i &\leftarrow s \bf{v}_i \\ \bf{r}_i &\leftarrow \bf{r}_i + \bf{v}_i \frac{\Delta t}{2} \end{aligned}\]

with the forces calculated at the start of the step. If there are constraints, velocity constraints are applied after the first update and position constraints after the last update, with the velocities corrected by the constrained displacement divided by $\Delta t$. The velocity scale factor $s$ is found by propagating the chain over $\Delta t$ in n_mts multiple time steps, each split into n_yoshida_suzuki steps of length $\delta = w_k \Delta t / n_\mathrm{mts}$ with Yoshida-Suzuki weights $w_k$. Starting from $s = 1$, each of these steps is

\[\begin{aligned} v_{\xi_M} &\leftarrow v_{\xi_M} + \frac{\delta}{2} G_M \\ v_{\xi_j} &\leftarrow v_{\xi_j} e^{-\delta v_{\xi_{j+1}} / 2} + \frac{\delta}{2} G_j e^{-\delta v_{\xi_{j+1}} / 4} \quad (j = M - 1, \ldots, 1) \\ s &\leftarrow s e^{-\delta v_{\xi_1}} \\ \xi_j &\leftarrow \xi_j + \delta v_{\xi_j} \quad (j = 1, \ldots, M) \\ v_{\xi_j} &\leftarrow v_{\xi_j} e^{-\delta v_{\xi_{j+1}} / 2} + \frac{\delta}{2} G_j e^{-\delta v_{\xi_{j+1}} / 4} \quad (j = 1, \ldots, M - 1) \\ v_{\xi_M} &\leftarrow v_{\xi_M} + \frac{\delta}{2} G_M \end{aligned}\]

where the kinetic energy in $G_1$ is that of the velocities after the first position update scaled by $s^2$, $G_1$ is updated after $s$ and each $G_{j+1}$ is updated after $v_{\xi_j}$. The chain is not propagated if the kinetic energy is close to zero.

The quantity

\[H = \sum_i \frac{1}{2} m_i \bf{v}_i^2 + V(\bf{r}) + \sum_{j=1}^M \frac{1}{2} Q_j v_{\xi_j}^2 + N_\mathrm{df} k T \xi_1 + \sum_{j=2}^M k T \xi_j\]

is conserved, apart from the error arising from the finite time step, if no coupling is used. The thermostat terms are the heat bath energy. The following keyword arguments are given to the loggers each step and can be used in a GeneralObservableLogger, for example to log $H$:

• current_forces: the forces used to update the velocities during the step.
• heat_bath_energy: the heat bath energy, in the energy units of the system.
• thermostat_positions: a vector of the thermostat positions $\xi_j$.
• thermostat_velocities: a vector of the thermostat velocities $v_{\xi_j}$, as numbers in the inverse time units of dt.

The vectors are modified during the simulation so should be copied if stored.

The thermostat positions and velocities start at zero on each call to simulate!, so two calls to simulate! are not the same as one longer call.

Arguments

• dt::T: the time step of the simulation.
• temperature::K: the equilibrium temperature of the simulation.
• damping::D=100*dt: the time scale $\tau$ of the coupling to the thermostats.
• chain_length::Integer=3: the number of thermostats in the chain $M$.
• n_mts::Integer=3: the number of multiple time steps used to propagate the chain over each time step.
• n_yoshida_suzuki::Integer=7: the number of Yoshida-Suzuki steps used for each multiple time step, 1, 3, 5 or 7.
• coupling::C=nothing: the coupling which applies during the simulation.
• remove_CM_motion=1: remove the center of mass motion every this number of steps, set to false or 0 to not remove center of mass motion.

OpenFESchedule()

Schedule of a Molly.OpenFEScheduler, passed to Molly.LambdaScheduler.

OpenFEScheduler(; dual=false, LJindividual=false, LJspecial=false,
                Cindividual=true, Cspecial=false, intraLJ=false,
                Tscaled=false)

Lambda scheduler that reproduces the relative free energy setup of OpenFE.

The electrostatics and sterics follow the same stages as Molly.DefaultLambdaScheduler, but the defaults use single topology with individual scaling of the Coulomb parameters. As in OpenFE, the charges of the 1-4 interactions are scaled for each state, and the alchemical atoms do not contribute to the Lennard-Jones dispersion correction. Use lj_softcore=:gapsys and coulomb_softcore=:scaled to match OpenFE.

A Molly.LambdaScheduler, see there for the keyword arguments.

OverdampedLangevin(; <keyword arguments>)

Simulates the overdamped Langevin equation using the Euler-Maruyama method.

Not compatible with constraints, will print a warning and continue without applying constraints. Not compatible with virial calculation. Not compatible with gradient calculation using Enzyme.

Arguments

• dt::S: the time step of the simulation.
• temperature::K: the equilibrium temperature of the simulation.
• friction::F: the friction coefficient of the simulation, e.g. 1.0u"ps^-1".
• remove_CM_motion=1: remove the center of mass motion every this number of steps, set to false or 0 to not remove center of mass motion.

PME(dist_cutoff, atoms, boundary; error_tol=0.0005, order=5,
    ϵr=1.0, coulomb_const=Molly.coulomb_const, fixed_charges=true, mesh_dims=nothing,
    scheduler=Molly.DefaultLambdaScheduler(), float_type_high=Float64,
    grad_safe=false, n_threads=Threads.nthreads())

Particle mesh Ewald (PME) summation for long range electrostatics implemented as an AtomsCalculators.jl calculator.

Should be used alongside the CoulombEwald pairwise interaction, which provides the short range term, and the EwaldExclusion specific interaction, which provides the exclusions for bonded atoms. dist_cutoff, error_tol, ϵr and coulomb_const should match these interactions. SetupPME provides PME parameters when setting up a system from a file.

dist_cutoff is the cutoff distance for short range interactions. The Ewald splitting parameter is

\[\alpha = \frac{\sqrt{-\ln(2 \epsilon_{\mathrm{tol}})}}{r_c}\]

where $r_c$ is dist_cutoff and $\epsilon_{\mathrm{tol}}$ is error_tol. order is the order of the B-spline interpolation onto the grid, ϵr is the relative permittivity and coulomb_const is the Coulomb constant. fixed_charges should be set to false if the partial charges can change, for example when using a polarizable force field. mesh_dims gives the number of grid points in each dimension, overriding the value chosen from error_tol. The charge of each atom is scaled by its electrostatic coupling, which scheduler, a Molly.LambdaScheduler, gives from the λ and alchemical role of the atom; this has no effect for atoms that are not alchemical. float_type_high is the float type used to sum the reciprocal space energy and the charges for the constant self energy and neutralizing terms. grad_safe should be set to true if gradients are going to be calculated with Enzyme.jl. n_threads is used to pre-allocate memory on CPU and plan the FFTs.

This implementation is based on the implementation in OpenMM, which is based on the smooth PME algorithm from Essmann et al. 1995.

Only compatible with 3D systems. Not compatible with infinite boundaries.

PME_λ(dist_cutoff, atoms, boundary; error_tol=0.0005, order=5, ϵr=1.0,
      coulomb_const=Molly.coulomb_const, fixed_charges=true, mesh_dims=nothing,
      scheduler=Molly.GROMACSLambdaABFEScheduler(), λ=0.0, grad_safe=false,
      n_threads=Threads.nthreads())

Particle mesh Ewald (PME) summation for long range electrostatics in an alchemical system, the λ version of PME that mixes the two end states as GROMACS does.

AbsoluteFESystem and RelativeFESystem convert a PME to this type when the scheduler is Molly.GROMACSLambdaABFEScheduler or Molly.GROMACSLambdaRBFEScheduler. These are the only schedulers it is compatible with, and they have to use dual topology. The charges of the two end states are spread onto separate grids and the potential energy is defined as

\[E = (1 - \lambda_q) E_{\mathrm{PME}}(q^A) + \lambda_q E_{\mathrm{PME}}(q^B)\]

where $E_{\mathrm{PME}}(q)$ is the PME potential energy for atomic charges $q$. $q^A$ and $q^B$ are the charges of the atoms at a global λ of 0 and 1, where the atoms that are switched off in an end state have zero charge. $\lambda_q$ is the electrostatic coupling that the scheduler gives at the global λ λ, which is $\min(2\lambda, 1)$ for Molly.GROMACSLambdaABFEScheduler and $\lambda$ for Molly.GROMACSLambdaRBFEScheduler. The forces and the virial are mixed in the same way. The potential energy is linear in $\lambda_q$, whereas PME with a scheduler uses one grid on which the charge of each atom is scaled.

Molly.set_lambda! changes the global λ of a system without rebuilding the grids. fixed_charges is accepted for consistency with PME but has no effect, as the charges of the end states are always recalculated from the atoms. See PME for the other arguments, the interactions that this should be used alongside and the systems that it is compatible with.

A potential of mean force (PMF) along a collective variable.

Returned by Molly.pmf_with_uncertainty, which describes the fields, Molly.mbar_pmf and Molly.pmf.

PMFDeconvolution(state; grid, coupling=nothing, cv=nothing, target_temp=nothing,
                 target_pressure=nothing)

Prepare a sampled PMF deconvolution object for an AWHState or a TSSState, which is passed as pmf to AWHSimulation or TSSSimulation.

The backend is selected from the type of state. grid is (min, max, n_bins) for a 1D PMF, ((mins...), (maxs...), (n_bins...)) for an N-dimensional PMF, or a tuple of explicit bin edge vectors. coupling(xi, state_i) should return the dimensionless bias energy at PMF coordinate xi in the thermodynamic state with index state_i; when omitted, Molly attempts to infer matching BiasPotential interactions from the state Hamiltonians. cv can give a custom function for the PMF coordinate of a sample, in which case coupling must also be given. target_temp and target_pressure set the target ensemble of the PMF, and are only supported for AWH.

The estimator follows the sampled deconvolution form used by AWH, see Lindahl et al. 2014: each sampled PMF coordinate is accumulated in a self-normalized weighted histogram using the inverse time-dependent effective bias at the observed PMF bin. Direct deconvolution of final state free energies is not supported. Call Molly.pmf on the returned object to extract the current PMF estimate.

Base type for pairwise interactions.

An alias for NBodyInteraction{2}. Custom pairwise interactions should subtype this.

ParameterPlan(sys, params_dic)

The mapping from the entries of params_dic to the parameters of sys, built once and reused on every call to inject_gradients.

Atom keys are described in extract_parameters. It is an error to give an atom key whose atoms have different values, or an interaction key whose interactions have different values, which extract_parameters leaves out, or to use a single topology system.

Building the plan is where the dictionary keys are formed and looked up. A plan is valid for any dictionary with the same keys, so it can be reused across a gradient calculation and across finite differencing, but it must be rebuilt if the set of keys or the system changes.

sys may have units. The dictionary values are then plain numbers in the units that each parameter is stored with in sys, as returned by extract_parameters.

PeriodicFlatBottomBias(k, r_fb, cv_target)

A flat-bottomed square (harmonic) bias on a collective variable (CV) towards a target value.

The bias is zero when the value of the collective variable does not deviate from cv_target by more than r_fb, and is square (harmonic) outside this range.

This variant handles periodicity in the CV wrapping around the (-π, π) range.

The potential energy is defined as

\[V(\boldsymbol{s}) = \frac{1}{2} k (|\boldsymbol{s} - \boldsymbol{s}_t| - r_{fb})^2 H\]

where $s$ and $s_t$ are the system and target CV values respectively, $\boldsymbol{s} - \boldsymbol{s}_t$ is wrapped into the range $[-\pi, \pi)$, and

\[H = \left\{ \begin{array}{cl} 0 & \text{if} & |\boldsymbol{s} - \boldsymbol{s}_t| < r_{fb} \\ 1 & \text{if} & |\boldsymbol{s} - \boldsymbol{s}_t| \geq r_{fb} \\ \end{array} \right.\]

Arguments

• k: The energy constant for the bias. Must be compliant with the System energy units.
• r_fb: Width of flat-bottom potential well. Inside this region the bias potential is always 0.
• cv_target: The target value of the collective variable.

PeriodicTorsion(; periodicities, phases, ks, proper=true)

A periodic torsion angle between four atoms.

phases are in radians. The potential energy is defined as

\[V(\phi) = \sum_{n=1}^N k_n (1 + \cos(n \phi - \phi_{s,n}))\]

where ϕ is the angle between the planes defined by atoms (i, j, k) and (j, k, l).

Only compatible with 3D systems.

PeriodicTorsionλ(; periodicities, phases, ks, proper=true,
                 λ_mixing=Molly.MinimumMixing(),
                 scheduler=Molly.DefaultLambdaScheduler())

A periodic torsion angle between four atoms in an alchemical system, the λ version of PeriodicTorsion.

AbsoluteFESystem and RelativeFESystem convert each PeriodicTorsion to this type. phases and ϕ are as for PeriodicTorsion. The weights below depend on the alchemical roles of the atoms and on $\lambda$, the λ of the atoms mixed with λ_mixing. scheduler is a Molly.LambdaScheduler, which selects dual or single topology.

In dual topology the potential energy is defined as

\[V(\phi) = w \sum_{n=1}^N k_n (1 + \cos(n \phi - \phi_{s,n}))\]

where $w = 1 - \lambda$ if the torsion involves a core atom with the parameters of state A, $w = \lambda$ if it involves a core atom with the parameters of state B and $w = 1$ if it involves only unique and environment atoms.

In single topology the first half of the terms are those of state A and the second half are those of state B, and the potential energy is defined as

\[V(\phi) = w_A V_A(\phi) + w_B V_B(\phi)\]

where $V_A$ and $V_B$ are the PeriodicTorsion potential energies of the two halves. $w_A = 1 - \lambda$ and $w_B = \lambda$ if the torsion involves a core atom, and $w_A = w_B = 1$ otherwise.

Only compatible with 3D systems.

PolynomialCutoff(dist_activation, dist_cutoff)

Cutoff that multiplies the potential by a 5th-degree polynomial switching function between an activation distance and a cutoff distance, smoothly reducing it to zero at the cutoff distance. Equivalent to the switching function used in OpenMM.

\[\begin{aligned} V_c(r) &= \begin{cases} V(r), r \le r_a \\ S(r) V(r), r_a < r \le r_c \\ 0, r > r_c \end{cases} \\ F_c(r) &= \begin{cases} F(r), r \le r_a \\ S(r) F(r) - S'(r) V(r), r_a < r \le r_c \\ 0, r > r_c \end{cases} \\ S(r) &= 1 - 6t^5 + 15t^4 - 10t^3 \\ S'(r) &= \frac{-30t^4 + 60t^3 - 30t^2}{r_c - r_a} \\ t &= \frac{r - r_a}{r_c - r_a} \end{aligned}\]

ProductMixing()

Mix λ parameters according to their product.

\[\lambda_{ij} = \lambda_i \lambda_j\]

This only works for interactions between two atoms. Give this when setting up an alchemical interaction.

QuartersLambdaScheduler(; dual=true, LJindividual=false, LJspecial=false,
                        Cindividual=false, Cspecial=false, intraLJ=false,
                        Tscaled=true)

Lambda scheduler that removes the atoms of the first state before adding the atoms of the second state.

Atoms that only exist in the first state lose their charges between global_λ = 0 and 0.25 and their Lennard-Jones interactions between 0.25 and 0.5. Atoms that only exist in the second state gain their Lennard-Jones interactions between 0.5 and 0.75 and their charges between 0.75 and 1. Core atoms are interpolated linearly over the full range.

A Molly.LambdaScheduler, see there for the keyword arguments.

QuartersSchedule()

Schedule of a Molly.QuartersLambdaScheduler, passed to Molly.LambdaScheduler.

RBTorsion(; c0, c1, c2, c3, c4, c5)

A Ryckaert-Bellemans torsion angle between four atoms.

See Ryckaert and Bellemans 1975. The potential energy is defined as

\[V(\phi) = \sum_{n=0}^{5} C_n \cos^n(\psi)\]

where $\psi = \phi - 180^{\circ}$ and $\phi$ is the angle between the planes defined by atoms (i, j, k) and (j, k, l). Coefficients using the other sign convention, i.e. powers of $\cos(\phi)$ rather than $\cos(\psi)$, can be converted by multiplying $C_n$ by $(-1)^n$.

Only compatible with 3D systems.

RBTorsionλ(; c0, c1, c2, c3, c4, c5, λ_mixing=Molly.MinimumMixing(),
           scheduler=Molly.DefaultLambdaScheduler())

A Ryckaert-Bellemans torsion angle between four atoms in an alchemical system, the λ version of RBTorsion.

AbsoluteFESystem and RelativeFESystem convert each RBTorsion to this type. The potential energy is defined as

\[V(\phi) = w \sum_{n=0}^{5} C_n \cos^n(\psi)\]

where $\psi$ and the sign convention of the coefficients are as for RBTorsion, and the weight $w$ depends on the alchemical roles of the atoms and on $\lambda$, the λ of the atoms mixed with λ_mixing. $w = 1 - \lambda$ if the torsion involves a core atom with the parameters of state A, $w = \lambda$ if it involves a core atom with the parameters of state B and $w = 1$ if it involves only unique and environment atoms.

scheduler is a Molly.LambdaScheduler, which has to use dual topology. Single topology is not supported for this interaction, so RelativeFESystem throws an error if the scheduler has dual=false.

Only compatible with 3D systems.

ReaxFF(; force_field_file, version=ReaxFF2008(), float_type=Float64,
        array_type=Array, smoothness_correction=false,
        bond_dist_cutoff=0.5, hbond_dist_cutoff=0.75,
        valence_cutoff=0.001, valence_prod_cutoff=0.00001,
        hbond_bo_cutoff=0.01)

Reactive force field (ReaxFF) implemented as an AtomsCalculators.jl calculator.

See van Duin et al. 2001. force_field_file specifies the path to a ReaxFF parameter file. version selects the functional form; currently only ReaxFF2008 is supported. float_type sets the precision of the force-field parameters and internal ReaxFF calculations (Float64 strongly recommended). array_type sets CPU (Array) or GPU storage and the device to run on, e.g. CuArray after loading CUDA.jl.

bond_dist_cutoff and hbond_dist_cutoff set the distance cutoffs for bond construction and hydrogen-bond interactions, respectively. They accept lengths with units; unitless values are interpreted as nanometers.

valence_cutoff sets the bond-order threshold for angle and torsion terms. valence_prod_cutoff sets the bond-order product threshold for angle terms when smoothness_correction=false. hbond_bo_cutoff sets the donor–hydrogen bond-order threshold for hydrogen bonds. These three thresholds are dimensionless.

smoothness_correction=true enables the Furman-Wales smoothness corrections to bond orders and energy terms (Furman and Wales 2019). The default is false.

Use System(coord_file_path, reaxff; ...) to set up coordinates, neighbor finding, and the associated ElectronegativityEqualizationMethod charge equilibration automatically. Use a fresh ReaxFF instance for each system, since it stores a system-specific cached state.

This interaction is only compatible with 3D systems. The virial is not yet supported.

ReaxFF2008()

Selects the 2008 ReaxFF functional form used by ReaxFF (Chenoweth et al. 2008). This is the default and currently supported ReaxFF version.

RectangularBoundary(x, y)
RectangularBoundary(x)

Rectangular 2D bounding box defined by two side lengths.

If one length is given then both sides will have that length. Setting one or more values to Inf gives no boundary in that dimension.

ReplicaExchangeLogger(n_replicas)
ReplicaExchangeLogger(T, n_replicas)

A logger that records exchanges in a replica exchange simulation.

The logged quantities include the number of exchange attempts (n_attempts), number of successful exchanges (n_exchanges), exchanged replica indices (indices), exchange steps (steps) and the value of Δ i.e. the argument of Metropolis rate for the exchanges (deltas).

ReplicaExchangeMD(; dt, exchange_time)

A generalized simulator for Replica Exchange Molecular Dynamics (REMD).

See Sugita and Okamoto 1999. Handles Temperature REMD, Hamiltonian REMD, or multi-dimensional combinations automatically based on the definitions in the ThermoStates of the ReplicaSystem. Run with simulate! or simulate_remd!, with exchanges attempted by Molly.remd_exchange!.

Arguments

• dt::DT: the time step of the simulation.
• exchange_time::ET: the time interval between replica exchange attempts.

ReplicaSystem(thermo_states, replica_coords; <keyword arguments>)

A wrapper for replicas in a generalized ensemble or replica exchange simulation (REMD), for example with ReplicaExchangeMD.

Instead of instantiating completely disjoint System objects, ReplicaSystem automatically compiles an AlchemicalPartition from the provided ThermoState vector. This isolates shared, unperturbed topological and interactive components (e.g., bulk solvent) from state-specific perturbations. During an exchange attempt, cross-energies are evaluated efficiently by querying only the necessary subset of perturbed interactions, completely avoiding redundant evaluations of the shared system. The pairs are split between the shared and the state-specific parts by the neighbor finder, so the systems need a neighbor finder if the states share pairwise interactions, for example in temperature REMD, and an error is thrown otherwise.

Furthermore, upon a successful exchange, coordinates and velocities are not physically swapped in memory; the system simply updates internal pointers mapping the physical replica to its new thermodynamic state.

This is a sub-type of AbstractSystem from AtomsBase.jl and implements the interface described there, routing standard atomic property queries back to the unperturbed master system.

When using ReplicaSystem with CellListMapNeighborFinder, the number of threads used for both the simulation of replicas and the neighbor finder should be set to be the same. This can be done by passing number_of_batches=(min(n, 8), n) to CellListMapNeighborFinder during construction where n is the number of threads to be used per replica.

Arguments

• thermo_states::AbstractArray{<:ThermoState}: An array of thermodynamic states defining the replica ladder. Each state encapsulates its specific System interactions, integrator, inverse temperature (beta), and pressure (p). The number of states dictates n_replicas.
• replica_coords: The coordinates of the atoms in each replica. The number of elements must equal the length of thermo_states.
• replica_velocities=nothing: The velocities of the atoms in each replica. If not provided, they default to zero velocities using the system's units.
• replica_boundaries=nothing: The bounding box for each replica. If not provided, it defaults to duplicating the boundary of the reference system (the first ThermoState).
• replica_neighbor_finders=nothing: The neighbor finder of each thermodynamic state. If not provided, each state gets a copy of the neighbor finder of its system.
• replica_loggers=nothing: Logger collections for each thermodynamic state, recording the replica currently in that state. Stateful logger objects and TrajectoryWriter file paths cannot be shared between states.
• exchange_logger=nothing: The logger used to record replica exchange attempts. If nothing, a default ReplicaExchangeLogger is used.
• initial_step::Int=0: Absolute MD step for a new or resumed replica simulation.
• data::DA=nothing: Arbitrary data associated with the replica system.
• reuse_neighbors::Bool=true: Whether to reuse the active system's neighbor list when calculating energies for perturbed state differences. Generally improves performance.

SHAKE_RATTLE(; n_atoms, dist_tolerance=1e-8u"nm", vel_tolerance=1e-8u"nm^2 * ps^-1",
             dist_constraints=nothing, angle_constraints=nothing,
             gpu_block_size=128, max_iters=25, strictness=:warn)

Constrain distances during a simulation using the SHAKE and RATTLE algorithms.

Either or both of dist_constraints and angle_constraints must be given. The constraints should always be on the CPU, even for a GPU system; the System constructor moves the constraint data to the device of the system. SetupSHAKE_RATTLE provides SHAKE/RATTLE parameters when setting up a system from a file.

Velocity constraints will be imposed for simulators that integrate velocities such as VelocityVerlet. See Ryckaert et al. 1977 for SHAKE, Andersen 1983 for RATTLE, Elber et al. 2011 for a derivation of the linear system solved to satisfy the RATTLE algorithm, and Kräutler et al. 2001 for the M-SHAKE algorithm.

Arguments

• n_atoms: number of atoms in the system.
• dist_tolerance=1e-8u"nm": the tolerance used to end the iterative procedure when calculating position constraints, should have the same units as the coordinates.
• vel_tolerance=1e-8u"nm^2 * ps^-1": the tolerance used to end the iterative procedure when calculating velocity constraints, should have the same units as the velocities times the coordinates.
• dist_constraints=nothing: a vector of DistanceConstraint objects that define the distance constraints to be applied. If nothing, no distance constraints are applied.
• angle_constraints=nothing: a vector of AngleConstraint objects that define the angle constraints to be applied. If nothing, no angle constraints are applied.
• gpu_block_size=128: the number of threads per block to use for GPU calculations.
• max_iters=25: the maximum number of iterations to perform when doing SHAKE. If this number of iterations is reached, some constraints may not be satisfied.
• strictness=:warn: determines behavior when encountering possible problems, options are :warn to emit warnings, :nowarn to suppress warnings or :error to error. Note that this only applies to constraint setup.

SetupCoulomb(; cutoff, coulomb_const=Molly.coulomb_const)

Set up the Coulomb electrostatic interaction with the cutoff cutoff, for example DistanceCutoff, and the Coulomb constant coulomb_const.

Passed to the System constructor from files, where it creates a Coulomb pairwise interaction. Passing a cutoff as the non-bonded method is the same as passing SetupCoulomb with that cutoff and the default Coulomb constant.

SetupCoulombReactionField(; solvent_dielectric=78.3, solute_dielectric=1.0,
                          coulomb_const=Molly.coulomb_const, exclusion_correction=false)

Set up the Coulomb electrostatic interaction modified using the reaction field approximation between two atoms.

Passed to the System constructor from files, where it creates a CoulombReactionField pairwise interaction. If exclusion_correction is true, the reaction field is also applied to the excluded and special pairs and to each atom with itself, using the CoulombReactionFieldExclusion specific interaction, which requires a neighbor finder.

SetupEwald(; error_tol=0.0005, approximate_erfc=true,
           coulomb_const=138.93545764u"kJ * mol^-1 * nm")

Set up Ewald summation for long range electrostatics.

Passed to the System constructor from files, where it creates a Ewald general interaction, a CoulombEwald pairwise interaction and a EwaldExclusion specific interaction.

error_tol is the error tolerance for Ewald summation. approximate_erfc determines whether to use a fast approximation to the erfc function.

SetupImplicitSolventGBN2(; solvent_dielectric=Molly.gb_solvent_dielectric,
                         solute_dielectric=Molly.gb_solute_dielectric,
                         kappa=0.0u"nm^-1", offset=Molly.gbn2_offset,
                         dist_cutoff=0.0u"nm", probe_radius=Molly.gb_probe_radius,
                         sa_factor=Molly.gb_sa_factor, use_ACE=true,
                         neck_scale=Molly.gbn2_neck_scale, neck_cut=Molly.gbn2_neck_cut,
                         element_to_radius=Molly.mbondi2_element_to_radius,
                         element_to_screen=Molly.gbn2_element_to_screen,
                         element_to_screen_nucleic=Molly.gbn2_element_to_screen_nucleic,
                         atom_params=Molly.gbn2_atom_params,
                         atom_params_nucleic=Molly.gbn2_atom_params_nucleic,
                         data_d0=Molly.gbn2_data_d0, data_m0=Molly.gbn2_data_m0,
                         coulomb_const=nothing)

Set up the GBn2 solvation model.

Passed to the System constructor from files, where it creates an ImplicitSolventGBN2 general interaction. See ImplicitSolventGBN2 for argument descriptions. coulomb_const=nothing uses the Coulomb constant of the non-bonded method.

SetupImplicitSolventOBC(; use_OBC2=false,
                        solvent_dielectric=Molly.gb_solvent_dielectric,
                        solute_dielectric=Molly.gb_solute_dielectric, kappa=0.0u"nm^-1",
                        offset=Molly.obc_offset, dist_cutoff=0.0u"nm",
                        probe_radius=Molly.gb_probe_radius,
                        sa_factor=Molly.gb_sa_factor, use_ACE=true,
                        element_to_radius=Molly.mbondi2_element_to_radius,
                        element_to_screen=Molly.obc_element_to_screen,
                        coulomb_const=nothing)

Set up the Onufriev-Bashford-Case GBSA model.

Passed to the System constructor from files, where it creates an ImplicitSolventOBC general interaction. See ImplicitSolventOBC for argument descriptions. coulomb_const=nothing uses the Coulomb constant of the non-bonded method.

SetupLINCS(; dist_tolerance=1e-6u"nm", vel_tolerance=1e-6u"nm^2 * ps^-1",
           n_rec=4, n_iter=1, iter_vel_correction::Bool=false,
           gpu_block_size::Integer=128)

Set up constraints using the LINCS (LINear Constraint Solver) algorithm.

Passed to the System constructor from files, where it creates a set of LINCS constraints. See LINCS for argument descriptions.

SetupPME(; error_tol=0.0005, approximate_erfc=true, mesh_dims=nothing, order=5,
         coulomb_const=138.93545764u"kJ * mol^-1 * nm")

Set up the particle mesh Ewald (PME) summation for long range electrostatics.

Passed to the System constructor from files, where it creates a PME general interaction, a CoulombEwald pairwise interaction and a EwaldExclusion specific interaction.

error_tol is the error tolerance for Ewald summation, which sets the Ewald splitting parameter

\[\alpha = \frac{\sqrt{-\ln(2 \epsilon_{\mathrm{tol}})}}{r_c}\]

where $r_c$ is the cutoff distance and $\epsilon_{\mathrm{tol}}$ is error_tol, as in OpenMM. approximate_erfc determines whether to use a fast approximation to the erfc function. mesh_dims determines the number of PME grid points in each dimension and defaults to a value chosen from error_tol. order is the order of the B-spline interpolation onto the grid. coulomb_const is the Coulomb constant used by all three interactions.

SetupSHAKE_RATTLE(; dist_tolerance=1e-6u"nm", vel_tolerance=1e-6u"nm^2 * ps^-1",
                  gpu_block_size=128, max_iters=25)

Set up constraints using the SHAKE and RATTLE algorithms.

Passed to the System constructor from files, where it creates a set of SHAKE_RATTLE constraints. See SHAKE_RATTLE for argument descriptions.

ShiftedForceCutoff(dist_cutoff)

Cutoff that shifts the force to be continuous at a specified cutoff distance.

\[\begin{aligned} V_c(r) &= \begin{cases} V(r) - (r-r_c) V'(r_c) - V(r_c), r \le r_c \\ 0, r > r_c \end{cases} \\ F_c(r) &= \begin{cases} F(r) - F(r_c), r \le r_c \\ 0, r > r_c \end{cases} \end{aligned}\]

ShiftedPotentialCutoff(dist_cutoff)

Cutoff that shifts the potential to be continuous at a specified cutoff distance.

\[\begin{aligned} V_c(r) &= \begin{cases} V(r) - V(r_c), r \le r_c \\ 0, r > r_c \end{cases} \\ F_c(r) &= \begin{cases} F(r), r \le r_c \\ 0, r > r_c \end{cases} \end{aligned}\]

SoftSphere(; cutoff=NoCutoff(), use_neighbors=true, shortcut=Molly.LJZeroShortcut(),
           σ_mixing=Molly.LorentzMixing(), ϵ_mixing=Molly.GeometricMixing())

The soft-sphere potential.

The potential energy is defined as

\[V(r_{ij}) = 4\varepsilon_{ij} \left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}\]

SparsePairMatrix(n_atoms, pairs; listed=true, array_type=Array)

A symmetric n_atoms x n_atoms AbstractMatrix{Bool} that stores a list of atom pairs.

Listed pairs have the value listed, all other off-diagonal entries have the value !listed and the diagonal is false. This makes it a sparse alternative to the dense eligible and special matrices given to neighbor finders, which take memory proportional to n_atoms^2:

• SparsePairMatrix(n_atoms, excluded_pairs; listed=false) is an eligible matrix where every pair of different atoms is eligible apart from the excluded pairs.
• SparsePairMatrix(n_atoms, special_pairs; listed=true) is a special matrix where only the given pairs are special.

Neighbor finders build these matrices when they are given n_atoms, excluded_pairs and special_pairs, so usually they do not have to be constructed directly.

pairs is an iterable of (i, j) pairs, in either order and with duplicates allowed, or a sparse matrix whose stored true entries are the pairs. array_type is the array type used to store the lists, for example CuArray to use the matrix with a neighbor finder on the GPU.

For each atom the atoms it is paired with are stored in ascending order, so the memory used is proportional to the number of pairs and looking up an entry means scanning the partners of one atom. Indexing into a matrix stored on the GPU from the CPU is not allowed, as for other GPU arrays.

SpecificForce1Atoms(f1)

Force on one atom arising from an interaction such as a position restraint.

SpecificForce2Atoms(f1, f2)

Forces on two atoms arising from an interaction such as a bond potential.

SpecificForce3Atoms(f1, f2, f3)

Forces on three atoms arising from an interaction such as a bond angle potential.

SpecificForce4Atoms(f1, f2, f3, f4)

Forces on four atoms arising from an interaction such as a torsion potential.

SpecificForce5Atoms(f1, f2, f3, f4, f5)

Forces on five atoms arising from an interaction such as a CMAP torsion potential.

SquareBias(k, cv_target)

A harmonic bias on a collective variable (CV) towards a target value.

The potential energy is defined as

\[V(\boldsymbol{s}) = \frac{1}{2} k (\boldsymbol{s} - \boldsymbol{s}_t)^2\]

where $s$ and $s_t$ are the system and target CV values respectively.

Arguments

• k: The energy constant for the bias. Must be compliant with the System energy units.
• cv_target: The target value of the collective variable.

StatisticalInefficiency(series; maxlag=nothing)

Estimate the statistical inefficiency of a time series from its integrated autocorrelation time, truncated using the initial positive sequence.

maxlag sets the maximum lag considered. The fields are inefficiency, stride (the inefficiency rounded up, useful for subsampling the series), input_length, effective_size (the number of uncorrelated samples) and lag (the truncation lag used).

Uses initial positive sequence (IPS) on paired lags to choose L. Uses normalized ACF of the mean-removed series. Includes the (1 - τ/N) taper in the sum.

SteepestDescentMinimizer(; <keyword arguments>)

Steepest descent energy minimization.

Arguments

• step_size::D=0.01u"nm": the initial maximum displacement.
• max_steps::Int=1_000: the maximum number of steps.
• tol::F=1_000.0u"kJ * mol^-1 * nm^-1": the maximum force below which to finish minimization.
• constraint_bond_constant::K=500_000.0u"kJ * mol^-1 * nm^-2": the force constant for the harmonic bonds that are used instead of constraints during minimisation. Set to nothing to not use harmonic bonds and ignore constraints. Unused if the system does not have constraints.
• log_stream::L=devnull: stream to print minimization progress to.

StormerVerlet(; <keyword arguments>)

The Störmer-Verlet integrator.

Uses a velocity-based update for the first step since the previous step coordinates are not available. This means that two calls to simulate! will only be approximately the same as one longer call. Coupling methods are not supported.

Arguments

• dt::T: the time step of the simulation.

System(; <keyword arguments>)

A physical system to be simulated.

Properties unused in the simulation or in analysis can be left with their default values. The minimal required arguments are atoms, coords and boundary. atoms and coords should have the same length, along with velocities and atoms_data if these are provided. This is a sub-type of AbstractSystem from AtomsBase.jl and implements the interface described there.

The arguments have to be consistent with each other and errors are thrown when they are not. In particular coords and velocities should have the same number of dimensions as boundary, should use the same float type as float_type, and should be on the same device (CPU or GPU) as atoms. The same applies to the arrays in specific_inter_lists and in the neighbor finder. Cases that are allowed but are often mistakes, such as an interaction cutoff more than half the minimum box side, are reported according to strictness.

Arguments

• atoms::A: the atoms, or atom equivalents, in the system. Can be of any type but should be a bits type if the GPU is used.
• coords::C: the coordinates of the atoms in the system. Typically a vector of SVectors of 2 or 3 dimensions.
• boundary::B: the bounding box in which the simulation takes place.
• velocities::V=zero(coords) * u"ps^-1": the velocities of the atoms in the system.
• atoms_data::AD=[]: other data associated with the atoms, allowing the atoms to be bits types and hence work on the GPU.
• topology::TO=nothing: topological information about the system such as which atoms are in the same molecule.
• pairwise_inters::PI=(): the pairwise interactions in the system, i.e. interactions between all or most atom pairs such as electrostatics. Should be a Tuple or NamedTuple of PairwiseInteractions.
• specific_inter_lists::SI=(): the specific interactions in the system, i.e. interactions between specific atoms such as bonds or angles. Should be a Tuple or NamedTuple.
• general_inters::GI=(): the general interactions in the system, i.e. interactions involving all atoms such as implicit solvent. Each should implement the AtomsCalculators.jl interface. Should be a Tuple or NamedTuple.
• constraints::CN=(): the constraints for bonds and angles in the system. Should be a Tuple or NamedTuple.
• virtual_sites::VS=[]: the virtual sites present in the system; these are mass-less particles determined by the positions of other atoms.
• neighbor_finder::NF=nothing: the neighbor finder used to find close atoms and save on computation. nothing means that no neighbor finder is used and all pairs of atoms are considered, in which case use_neighbors should return false for the pairwise interactions, e.g. by giving use_neighbors=false to the built-in pairwise interactions.
• loggers::L=(): the loggers that record properties of interest during a simulation.
• force_units::F=u"kJ * mol^-1 * nm^-1": the units of force of the system. Should be set to NoUnits if units are not being used.
• energy_units::E=u"kJ * mol^-1": the units of energy of the system. Should be set to NoUnits if units are not being used.
• k::K=Unitful.k or Unitful.k * Unitful.Na: the Boltzmann constant, which may be modified in some simulations. k is chosen based on the energy_units given.
• float_type::T: the floating point type of the system, read from the boundary by default.
• float_type_high::TH: the floating point type used for accumulation where higher precision is useful, such as the potential energy and the virial. Float64 by default, or float_type on backends that do not support Float64 such as Metal.
• data::DA=nothing: arbitrary data associated with the system.
• grad_safe=false: should be set to true if the system is going to be used with Enzyme.jl.
• strictness=:warn: determines behavior when encountering possible problems, options are :warn to emit warnings, :nowarn to suppress warnings or :error to error.

System(coordinate_file, force_field; <keyword arguments>)

Read a coordinate file in a file format readable by Chemfiles and apply a force field, such as a MolecularForceField, to it.

Each residue is matched to a residue template of the force field, and a Molly.MissingResidueTemplateError is thrown if no template matches.

Arguments

• boundary=nothing: the bounding box used for simulation, read from the file by default.
• velocities=nothing: the velocities of the atoms in the system, set to zero by default.
• units::Bool=true: whether to use Unitful quantities. Should match the corresponding argument used for the force field argument.
• array_type=Array: the array type for the simulation, for example use CuArray or ROCArray for GPU support.
• float_type: the floating point type of the system, defaults to Float64 on CPU and Float32 on GPU.
• float_type_high: the floating point type used for accumulation where higher precision is useful, such as the potential energy and the virial. Float64 by default, or float_type on backends that do not support Float64 such as Metal.
• dist_cutoff=1.0u"nm": cutoff distance for the non-bonded interactions.
• dist_buffer=nothing: distance added to dist_cutoff when calculating neighbor lists. When neighbor_finder_auto_update=true this only affects performance and is chosen based on the neighbour finder type, when that is false this needs to be selected to avoid atoms crossing the buffer distance and defaults to 0.2u"nm".
• constraints=:none: which constraints to apply during the simulation, options are :none, :hbonds (bonds involving hydrogen), :allbonds (all bonds) and :hangles (all bonds plus H-X-H and H-O-X angles). Note that not all options may be supported depending on the bonding topology. Urey-Bradley bond terms are not treated as bonds when adding constraints.
• rigid_water=false: whether to constrain the bonds and angle in water molecules. Applied on top of constraints, so constraints=:hangles and rigid_water=false gives rigid water.
• constraint_algorithm=SetupLINCS(): Constraint algorithm to use for enforcing the constraints. Can be an instance of SetupLINCS or SetupSHAKE_RATTLE.
• nonbonded_method=SetupCoulombReactionField(): method for long range interaction summation, can be an instance of SetupCoulombReactionField (reaction field method), SetupPME (particle mesh Ewald summation), SetupEwald (Ewald summation, slow), SetupCoulomb or a cutoff like DistanceCutoff (short range only). For a cutoff, the cutoff distance must not exceed dist_cutoff. The Coulomb constant of the setup type is used for all the electrostatic interactions, including implicit solvent unless that sets its own.
• lj_cutoff=nothing: the cutoff method to use for the Lennard-Jones interaction, by default uses DistanceCutoff(dist_cutoff). The cutoff distance must not exceed dist_cutoff. If the long-range dispersion correction is being used, this must be a DistanceCutoff, ShiftedPotentialCutoff or PolynomialCutoff.
• dispersion_correction=nothing: whether to use the long-range dispersion correction, LJDispersionCorrection, or DExpDispersionCorrection for a double exponential force field. Defaults to the force field setting, which defaults to true for Lennard-Jones and false for double exponential.
• hydrogen_mass=false: if set to a number or unitful value (e.g. 2 or 2u"g/mol") then hydrogens are given the specified mass, which is transferred from the atom they are bonded to. This is applied before constraints are generated.
• center_coords::Bool=true: whether to center the coordinates in the simulation box.
• neighbor_finder_type=nothing: which type of neighbor finder to use. The default nothing chooses CellListMapNeighborFinder on CPU, GPUNeighborFinder on CUDA compatible GPUs and GPUCellListNeighborFinder on non-CUDA compatible GPUs, falling back to DistanceNeighborFinder there when the box is too small for a cell list. Nothing, the type of nothing, gives a system with no neighbor finder.
• neighbor_finder_auto_update=true: whether the neighbor finder rebuilds the neighbor list only when the atoms have moved far enough since the last rebuild that a pair within dist_cutoff could be missing from it, given the buffer distance dist_buffer. If true, neighbor_finder_n_steps is ignored.
• neighbor_finder_n_steps=10: with neighbor_finder_auto_update=false, the number of steps between neighbor finder updates. Can be tuned along with dist_buffer to ensure that particles do not cross the buffer distance during the update interval.
• launch_config=CUDALaunchConfig(): stored CUDA launch overrides for this system. Ignored on CPU and non-CUDA GPU backends.
• autotune_launch=true: whether to autotune CUDA launch parameters at the end of setup. This is a no-op on CPU and non-CUDA GPU backends.
• data=nothing: arbitrary data associated with the system.
• implicit_solvent=nothing: the implicit solvent model to use, can be nothing or an instance of SetupImplicitSolventOBC or SetupImplicitSolventGBN2.
• disulfide_bonds=true: whether or not to look for disulfide bonds between CYS residues in the structure file and add them to the topology. Uses geometric arguments to assign them.
• loggers=(): the loggers that record properties of interest during a simulation.
• n_threads=Threads.nthreads(): the number of threads the simulation will likely be run on, used for example when setting up PME and implicit solvent. Only relevant when running on CPU.
• grad_safe=false: should be set to true if the system is going to be used with Enzyme.jl.
• strictness=:warn: determines behavior when encountering possible problems, options are :warn to emit warnings, :nowarn to suppress warnings or :error to error.
• charge_types=:auto: how atoms are grouped into charge parameters for extract_parameters and inject_gradients when the force field takes charges from the residue templates, options are :auto, :atom_type, :template_type or :template_atom. See High-level gradient interface.

System(abstract_system; <keyword arguments>)

Convert an AtomsBase AbstractSystem to a Molly System.

The keyword arguments force_units and energy_units should be set as appropriate. Other keyword arguments are the same as for the main System constructor.

System(coord_file_path::AbstractString, reaxff::ReaxFF; eem_tolerance=1e-6,
       eem_max_iterations=200, eem_float_type=Float64, initial_charges=nothing,
       neighbor_skin=0.0, neighbor_auto_update=true, neighbor_refresh_n_steps=1,
       boundary=nothing, array_type=Array, velocities=nothing, units=true,
       loggers=(), strictness=:warn)

Construct a ReaxFF system from a coordinate file supported by Chemfiles.jl. The file must identify each atom's element, and every element must be present in the force field file. This constructor adds both ElectronegativityEqualizationMethod and ReaxFF to the system's general interactions and selects a neighbor finder automatically.

eem_float_type sets the charge-equilibration precision independently of the ReaxFF precision (Float64 recommended). eem_tolerance and eem_max_iterations control convergence; failure to converge raises an error. initial_charges contains the initial atomic partial charges in units of elementary charge and defaults to zeros. Its sum sets the total system charge that is conserved by charge equilibration.

neighbor_skin is a non-negative distance added to the largest ReaxFF distance cutoff. It accepts lengths with units; unitless values are interpreted as nanometers, regardless of units. Without a skin, the default, the neighbor list is rebuilt every step. With a positive skin and neighbor_auto_update=true, the list is rebuilt only when the atoms have moved far enough since the last rebuild that a pair within the ReaxFF cutoffs could be missing from it. With neighbor_auto_update=false, neighbor_refresh_n_steps sets the scheduled neighbor finder rebuild interval and must be at least one. An interval greater than one requires a positive skin distance. Select a skin large enough to cover atomic motion between rebuilds.

array_type must match that specified in the ReaxFF instance. The boundary is read from the file unless boundary is supplied, and velocities default to zero. Supplied boundaries and velocities must match the ReaxFF floating-point type and the units setting. With units=true, quantities carry units; with units=false, numerical values use nm, ps, g/mol, and kJ/mol units. loggers are passed to the system unchanged.

strictness determines behavior when encountering possible problems, such as a box smaller than twice the largest ReaxFF cutoff or a ReaxFF instance that has already been used to set up a system, options are :warn to emit warnings, :nowarn to suppress warnings or :error to error.

Use a fresh ReaxFF instance for each system, since it stores system-specific cached state.

System(crystal; <keyword arguments>)

Construct a System from a SimpleCrystals.jl Crystal struct.

Properties unused in the simulation or in analysis can be left with their default values. atoms, atoms_data, coords and boundary are automatically calculated from the Crystal struct. Extra atom parameters like σ have to be added manually after construction using the convenience constructor System(sys; <keyword arguments>).

System(sys; <keyword arguments>)

Convenience constructor for changing properties in a System.

The System is returned with the provided keyword arguments modified. Give deepcopy(sys) as the argument to make a new copy of the system.

Changing float_type, units, array_type or grad_safe converts the system, for example System(sys; float_type=Float64, units=false, array_type=Array) gives a unitless Float64 copy on the CPU. On a device without Float64 support, such as Metal, float_type defaults to Float32. Everything except atoms_data, topology and data is converted, including values given as keyword arguments. Units are stripped without conversion and can not be added. When the float type, units or device change, loggers are dropped with a warning unless given with the loggers keyword argument. On a new device, or with grad_safe=true, the neighbor finder is the one setup would choose. Custom types that can not be rebuilt from their converted fields should define a method of Molly.convert_component(conversion, x).

TSSHistoryForgetting(; alpha=0.19, n_epochs=16, phi=nothing)

Configure geometric history forgetting for TSS free-energy estimates, given as history_forgetting to TSSState.

At TSS iteration $t$ the history before approximately $\alpha t$ is dropped, so alpha is the forgotten fraction of the current history time. The retained history is grouped into logarithmically spaced epochs: n_epochs sets the target number of retained epochs, and phi optionally overrides the epoch growth factor inferred from alpha and n_epochs.

TSSSample

One TSS sample, passed to the sample loggers of a TSSSimulation by Molly.log_tss_sample!.

A sample is collected after each block of n_md_steps molecular dynamics steps, from the coordinates at the end of the block, so a cycle collects self_adjustment_steps samples per replica. As the system loggers log the initial configuration at the start step of a new simulation, the initial configuration of each replica is also passed as a sample, at that step, in the first simulate! call; its weights and log_den use the window estimator before any update.

Fields

• step::Int: the absolute MD step of the sample, on the same step count as the system loggers.
• replica_index::Int: the replica that produced the sample.
• system: the replica's active System, holding the sampled coordinates, boundary and velocities.
• active_state::Int: the global index of the thermodynamic state that propagated the coordinates.
• active_window::Int: the window of the cycle.
• state_indices: the global state indices of the window.
• reduced_potentials: the reduced potentials of the sampled coordinates in the window states, in the order of state_indices.
• log_den: the log normalization logsumexp(f .+ log_dens .- reduced_potentials), with the window estimator's free energies f and log sampling densities log_dens at the time of the sample.
• weights: the conditional probabilities of the window states given the sampled coordinates, from which the next active state is drawn.

system, state_indices, reduced_potentials and weights reference the simulation's own arrays, which are reused, so they are valid only during the call to Molly.log_tss_sample! and must not be modified. A logger that keeps data must copy it. Nothing is copied off the device for the logger: for a GPU system, system.coords is a device array.

TSSSampleLogger(observable::Function, T, n_steps)

A TSS sample logger that stores observable(sample)::T for each Molly.TSSSample whose step is a multiple of n_steps.

Samples are collected every n_md_steps MD steps, so with n_steps a multiple of n_md_steps it records every n_steps ÷ n_md_steps-th sample of each replica, at the same steps as a system logger with the same n_steps, including the initial configuration when the start step is a multiple of n_steps. The observable must copy any array it keeps from the sample. Access the stored values with values(logger).

TSSSimulation(state::TSSState; n_md_steps, n_cycles,
              self_adjustment_steps=1, log_freq=1000,
              n_replicas=nothing, first_states=nothing,
              first_windows=nothing, replica_active_states=nothing,
              loggers=(), replica_loggers=nothing,
              sample_loggers=(), replica_sample_loggers=nothing,
              pmf=nothing, frozen=false, initial_step=0)

A windowed TSS simulation driver.

See Predescu et al. 2024. Each cycle runs self_adjustment_steps blocks of n_md_steps molecular dynamics steps, collects TSS samples, updates window estimators unless frozen=true, and optionally records sampled PMF deconvolution diagnostics. initial_step sets the absolute MD step for a new or resumed simulation. loggers supplies the logger collection for a single-replica simulation. For multi-replica TSS, use replica_loggers to supply one logger collection per newly constructed replica active system.

sample_loggers supplies a collection of TSS sample loggers for a single-replica simulation, and replica_sample_loggers one collection per replica for multi-replica TSS, with distinct logger objects for each replica. Molly.log_tss_sample! is called with each logger for every collected Molly.TSSSample, in adaptive and frozen runs, and for the initial configuration at the steps where the system loggers log it. A replica's sample loggers are called by the task that runs that replica, in sample order, so with replica_parallel=:threads the replicas call their own loggers concurrently. Without sample loggers, no sample is constructed. Pass replica_rngs to simulate! to supply one RNG stream per replica; otherwise, each replica RNG is reseeded from the top-level rng.

TSSState(thermo_states; graph=nothing, first_state=1, first_window=nothing,
         gamma=nothing, initial_f=nothing, ETA=2.0, dens_reg=1e-6,
         reuse_neighbors=true, history_forgetting=nothing,
         adaptive_gamma=nothing, global_visit_control=true,
         visit_control_tolerance=1e-8,
         visit_control_max_iterations=1000,
         visit_control_damping=1.0, pi_regularization=1e-3)

Mutable state for a Times Square Sampling (TSS) simulation over a set of thermodynamic states.

See Predescu et al. 2024. thermo_states defines the global expanded ensemble. graph may be omitted for a single-window TSS run, or provided as a graph from Molly.tss_grid_graph for windowed TSS. history_forgetting can be a Molly.TSSHistoryForgetting. The state owns the active thermodynamic state, local window estimators, visit-control coupling state, and diagnostic histories. Loggers are attached by TSSSimulation, not by TSSState. History-dependent general interactions such as Metadynamics are not supported.

ThermoState(system::System, integrator; <keyword arguments>)

Thermodynamic state wrapper carrying the system, integrator, and derived thermodynamic properties (inverse temperature beta and pressure p). This serves as the definitive container for a single thermodynamic state across all generalized ensemble methods.

Arguments

• system::System: The simulation system used to evaluate energies.
• integrator: The integrator used to simulate the system.
• temperature=nothing: Explicit target temperature. If nothing, it is inferred from the integrator's thermostat or implicit temperature coupling.
• pressure=nothing: Explicit target pressure. If nothing, it is inferred from the integrator's barostat.
• name::AbstractString=nothing: A label for the state. If not provided, a default name based on the temperature and pressure is generated.

If system has loggers, construction emits a warning: generalized ensemble methods ignore ThermoState system loggers and use loggers attached to the appropriate generalized-ensemble object instead.

ThermodynamicCycle(legs::ThermodynamicCycleLeg...; correction=nothing)

A free energy composed from the free-energy differences of several legs.

Molly.free_energy returns $\sum_i c_i \Delta F_i$ plus correction, and Molly.free_energy_uncertainty its standard error.

correction is an analytic term in energy units, e.g. a standard-state correction (see Molly.standard_state_correction), the release of restraints or a long-range correction. nothing adds no term.

Each leg is converted to energy units with the temperature of its own states, so legs at different temperatures may be combined. A source raises an error when the two states of one leg have different temperatures, since their difference would need a single $k_B T$.

ThermodynamicCycleLeg(source, from => to; coefficient=1.0, name=Symbol(""))

One leg of a ThermodynamicCycle: the free-energy difference $F_\text{to} - F_\text{from}$ reported by source, multiplied by coefficient (usually ±1).

from and to are state indices or ThermoState names. source is any object that implements Molly.free_energy_difference and Molly.free_energy_difference_uncertainty, such as a TSSState, an AWHState, a Molly.MBARResult or a PMF.

TimeCorrelationLogger(observableA::Function, observableB::Function,
                        TA::DataType, TB::DataType,
                        observable_length::Integer, n_correlation::Integer)

A time correlation logger.

Estimates statistical correlations of normalized form

\[C(t)=\frac{\langle A_t\cdot B_0\rangle -\langle A\rangle\cdot \langle B\rangle}{\sqrt{\langle |A|^2\rangle\langle |B|^2\rangle}}\]

or unnormalized form

\[C(t)=\langle A_t\cdot B_0\rangle -\langle A \rangle\cdot \langle B\rangle\]

These can be used to estimate statistical error, or to compute transport coefficients from Green-Kubo type formulas. A and B are observables, functions of the form observable(sys::System, neighbors, step_n, buffers; n_threads::Integer). The return values of A and B can be of scalar or vector type (including Vector{SVector{...}}, like positions or velocities) and must implement dot.

n_correlation should typically be chosen so that dt * n_correlation > t_corr, where dt is the simulation timestep and t_corr is the decorrelation time for the considered system and observables. For the purpose of numerical stability, the logger internally records sums instead of running averages. The normalized and unnormalized form of the correlation function can be retrieved through values(logger::TimeCorrelationLogger; normalize::Bool).

Arguments

• observableA::Function: the function corresponding to observable A.
• observableB::Function: the function corresponding to observable B.
• TA::DataType: the type returned by observableA, supporting zero(TA).
• TB::DataType: the type returned by observableB, supporting zero(TB).
• observable_length::Integer: the length of the observables if they are vectors, or 1 if they are scalar-valued.
• n_correlation::Integer: the length of the computed correlation vector.

TrajectoryWriter(n_steps, filepath; format="", correction=:pbc, atom_inds=[],
                 excluded_res=String[], write_velocities=false, write_boundary=true,
                 overwrite=false, append_warn=true)

Write 3D structures to a file throughout a simulation.

Uses Chemfiles.jl to write to one of a variety of formats including DCD, XTC, PDB, CIF, MOL2, SDF, TRR and XYZ. The full list of file formats can be found in the Chemfiles docs. By default the format is guessed from the file extension but it can also be given as a string, e.g. format="DCD". BioStructures.jl is used to write to the PDB format.

The atom indices to be written can be given as a list or range to atom_inds, with all atoms being written by default. Residue names to be excluded can be given as excluded_res. Velocities can be written in addition to coordinates by setting write_velocities=true. Chemfiles does not support writing velocities to all file formats. The correction to be applied to the molecules is chosen with correction. :pbc, the default, keeps molecules whole, whereas :wrap wraps all atoms inside the simulation box regardless of connectivity.

The System should have atoms_data defined, and topology if bonding information is required. The file will be appended to, so should be deleted before simulation if it already exists. Setting overwrite=true deletes the file when the logger is constructed if it exists, rather than warning and appending to it. Setting append_warn=false appends to an existing file without the warning.

Not compatible with 2D systems. For the PDB format, the box size for the CRYST1 record is taken from the first snapshot; different box sizes at later snapshots will not be recorded. The CRYST1 record is not written for infinite boundaries.

Many formats store coordinates at reduced precision, for example 0.001 Å for PDB, about 0.001 nm for XTC and single precision for DCD. Energies calculated from frames read back, for example with EnsembleSystem, are affected by this, so a full precision format or CoordinatesLogger should be used if that matters.

TreeNeighborFinder(; n_atoms, dist_cutoff, excluded_pairs=(), special_pairs=(),
                   auto_update=true, n_steps=10, strictness=:warn)
TreeNeighborFinder(; eligible, dist_cutoff, special=nothing, auto_update=true,
                   n_steps=10, strictness=:warn)

Find close atoms by distance using a tree search.

dist_cutoff is the neighbor search distance, which should be the interaction cutoff distance plus a buffer distance since the list is reused until it is rebuilt. auto_update and n_steps determine when the list is rebuilt, as for DistanceNeighborFinder. The pairs of atoms that can interact are given by n_atoms, excluded_pairs and special_pairs, or by the eligible and special matrices, as for DistanceNeighborFinder. The pairs are stored on the CPU.

Cannot be used if one or more dimensions has infinite boundaries. Cannot be used with TriclinicBoundary.

TriclinicBoundary(v1, v2, v3; approx_images=true)
TriclinicBoundary(SVector(v1, v2, v3); approx_images=true)
TriclinicBoundary(SVector(l1, l2, l3), SVector(α, β, γ); approx_images=true)
TriclinicBoundary(arr; approx_images=true)

Triclinic 3D bounding box defined by 3 SVector{3} basis vectors or basis vector lengths and angles α/β/γ in radians.

The first basis vector must point along the x-axis and the second must lie in the xy plane. An approximation is used to find the closest periodic image when using the minimum image convention. The approximation is correct for distances shorter than half the shortest box height/width. Setting the keyword argument approx_images to false means the exact closest image is found, which is slower.

Not compatible with infinite boundaries.

UreyBradley(; kangle, θ0, kbond, r0)

An interaction between three atoms consisting of a harmonic bond angle and a harmonic bond between the outer atoms.

θ0 is in radians. The second atom is the middle atom. The potential energy is defined as

\[V(\theta, r) = \frac{1}{2} k_a (\theta - \theta_0)^2 + \frac{1}{2} k_b (r - r_0)^2\]

UreyBradleyλ(; kangle, θ0, kbond, r0, λ_mixing=Molly.MinimumMixing(),
             scheduler=Molly.DefaultLambdaScheduler())

An interaction between three atoms in an alchemical system consisting of a harmonic bond angle and a harmonic bond between the outer atoms, the λ version of UreyBradley.

AbsoluteFESystem and RelativeFESystem convert each UreyBradley to this type. The potential energy is defined as

\[V(\theta, r) = w \left( \frac{1}{2} k_a (\theta - \theta_0)^2 + \frac{1}{2} k_b (r - r_0)^2 \right)\]

where θ0 and the atom order are as for UreyBradley, and the weight $w$ and the parameters depend on the alchemical roles of the atoms and on $\lambda$, the λ of the atoms mixed with λ_mixing. scheduler is a Molly.LambdaScheduler, which selects dual or single topology.

In dual topology kangle, θ0, kbond and r0 are single values. $w = 1$ if the interaction involves an atom unique to one state or only environment atoms, otherwise $w = 1 - \lambda$ if it involves a core atom with the parameters of state A and $w = \lambda$ if it involves a core atom with the parameters of state B.

In single topology $w = 1$ and kangle, θ0, kbond and r0 are tuples of the values in states A and B. Each parameter $p$ is $(1 - \lambda) p^A + \lambda p^B$ if the interaction involves a core or unique atom and $p^A$ otherwise.

VelocityRescaleThermostat(temperature, coupling_const; n_steps=1) <: AbstractThermostat

The stochastic velocity rescaling thermostat.

See Bussi et al. 2007. In brief, acts like the BerendsenThermostat but adds a stochastic term, allowing correct sampling of isothermal ensembles.

Let $\Delta t$ be the simulation time step multiplied by n_steps, $\tau$ the coupling constant and $N_f$ the number of degrees of freedom used to calculate the instantaneous temperature of the system. Then $K = \frac{1}{2} \sum_i m_i v_i^2$ is the current kinetic energy and $\bar{K} = \frac{1}{2} N_f k_B T_0$ is the target kinetic energy for the reference temperature $T_0$.

Define $c = e^{-\Delta t/\tau}$. Draw $R \sim \mathcal{N}(0, 1)$ and $S \sim \chi^{2}_{N_f-1}$. Then

\[\lambda^2 = c + (1-c) \cdot \frac{\bar K}{N_f K} \cdot (R^2 + S)\;+\; 2 \cdot \sqrt{c(1-c) \frac{\bar K}{N_f K}} \cdot R, \qquad v' = \lambda\,v .\]

VelocityVerlet(; <keyword arguments>)

The velocity Verlet integrator.

Arguments

• dt::T: the time step of the simulation.
• coupling::C=nothing: the coupling which applies during the simulation.
• remove_CM_motion=1: remove the center of mass motion every this number of steps, set to false or 0 to not remove center of mass motion.

Verlet(; <keyword arguments>)

The leapfrog Verlet integrator.

This is a leapfrog integrator, so the velocities are offset by half a time step behind the positions.

Arguments

• dt::T: the time step of the simulation.
• coupling::C=nothing: the coupling which applies during the simulation.
• remove_CM_motion=1: remove the center of mass motion every this number of steps, set to false or 0 to not remove center of mass motion.

A virtual site, with coordinates determined by the coordinates of other atoms.

Created with OneParticleSite, TwoParticleAverageSite, ThreeParticleAverageSite, OutOfPlaneSite or LocalCoordinatesSite and given to a System as virtual_sites.

WaldmanHaglerMixing()

Mix σ and ϵ parameters according to the Waldman-Hagler rules.

See Waldman and Hagler 1993.

\[\begin{aligned} \sigma_{ij} &= \left( \frac{\sigma_i^6 + \sigma_j^6}{2} \right)^{\frac{1}{6}} \\ \varepsilon_{ij} &= 2 \sqrt{\varepsilon_i \varepsilon_j} \frac{\sigma_i^3 \sigma_j^3}{\sigma_i^6 + \sigma_j^6} \end{aligned}\]

Give this when setting up a pairwise interaction such as LennardJones.

WellTempered(bias_factor, temperature)

Standard well-tempered Metadynamics height decay (Barducci, Bussi & Parrinello, 2008):

\[h(\boldsymbol{s}) = h_0 \exp\left(-\frac{V_{bias}(\boldsymbol{s})}{k_B (\gamma - 1) T}\right)\]

where $h_0$ is the base height, $V_{bias}$ the bias accumulated so far at $\boldsymbol{s}$, $\gamma$ is bias_factor (> 1), $T$ is temperature and $k_B$ the Boltzmann constant of the System. Deposits become self-limiting, so the bias converges to an estimate of the free energy (scaled by bias_factor) instead of growing without bound.

Arguments

• bias_factor: Temperature boost $\gamma > 1$; larger tempers more slowly (approaching untempered Metadynamics as $\gamma \to \infty$).
• temperature: Temperature $T$ of the System.

Yukawa(; cutoff=NoCutoff(), use_neighbors=true, weight_special=1,
       coulomb_const=Molly.coulomb_const, kappa=1.0u"nm^-1")

The Yukawa electrostatic interaction between two atoms.

The potential energy is defined as

\[V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}} \exp(-\kappa r_{ij})\]

and the force on each atom by

\[\vec{F}_i = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}^2} \exp(-\kappa r_{ij})\left(\kappa r_{ij} + 1\right) \frac{\vec{r}_{ij}}{r_{ij}}\]

where $\kappa$ is kappa, the inverse of the screening length.

Alchemical role, in dual topology, of the copy of a core atom with the parameters of state A. Its interactions are scaled by 1 - global_λ.

Alchemical role, in dual topology, of the copy of a core atom with the parameters of state B. Its interactions are scaled by global_λ.

Alchemical role, in single topology, of a core atom present in both end states. Its parameters are interpolated between the A and B values.

Alchemical role of an atom only present in state A, including the decoupled atoms of an absolute transformation. It is decoupled from global_λ = 0 to 1, following the scheduler.

Alchemical role of an atom that is not alchemical, e.g. solvent, protein or ions. Its interactions are never scaled. The default of Atom.

Alchemical role of an atom only present in state B. It is coupled from global_λ = 0 to 1, following the scheduler.

The Coulomb constant, 1 / 4πϵ0, where ϵ0 is the vacuum permittivity.

values(logger)
values(logger::TimeCorrelationLogger; normalize::Bool=true)
values(logger::AverageObservableLogger; std::Bool=true)

Access the stored observations in a logger.

A_mixing(mixing, atom_i, atom_j, args...)

The rule for mixing A (Buckingham) parameters for a given mixing scheme.

Custom mixing schemes can implement this function.

AbsoluteFESystem(sys, global_λ, mapping; temp=298.0u"K",
                 scheduler=Molly.DefaultLambdaScheduler(dual=true), loggers=(),
                 lj_softcore=:gapsys, coulomb_softcore=:gapsys)

Set up an alchemical system for an absolute free energy calculation, in which the atoms in mapping are decoupled from the rest of sys.

The mapped atoms are fully coupled at global_λ = 0 and fully decoupled at global_λ = 1. The atom order of sys is kept. Interactions are replaced by their λ-dependent versions and random velocities are generated at temp.

Arguments

• sys: the System containing the atoms to decouple.
• global_λ: the global λ of the system, from 0 (coupled) to 1 (decoupled).
• mapping: a vector with the indices of the atoms to decouple.
• temp=298.0u"K": the temperature used to generate random velocities.
• scheduler=Molly.DefaultLambdaScheduler(dual=true): the lambda scheduler that turns global_λ into the couplings of the steric, electrostatic and bonded interactions, for example Molly.DefaultLambdaScheduler, Molly.GROMACSLambdaABFEScheduler or Molly.OpenFEScheduler. Absolute systems require dual topology (dual=true).
• loggers=(): the loggers that record properties of interest during a simulation.
• lj_softcore=:gapsys: the Lennard-Jones soft core, :beutler (Beutler et al. 1994), :gapsys (Gapsys et al. 2012) or :scaled, where the potential is scaled directly by λ. The 1-4 Lennard-Jones interactions of force fields with separate 1-4 parameters (e.g. CHARMM) use the same soft core. A DoubleExponential interaction always becomes a DoubleExponentialSoftCore, its own soft core.
• coulomb_softcore=:gapsys: the Coulomb soft core, :beutler, :gapsys or :scaled, where the potential is scaled directly by λ.