Molly API

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

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

ANIPotential(path; T=Float32, 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.

Note: the ANI-2x weights are Float32, so the energy/force paths run in Float32 internally regardless of the system's coordinate type.

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.
• num_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.

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.

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 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, use_neighbors, shortcut, σ_mixing, ϵ_mixing,
              λ_mixing, weight_special)

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)

Data associated with an atom.

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=1) <: 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, use_neighbors, shortcut, A_mixing, B_mixing,
           C_mixing, weight_special)

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=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=1) <: 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.

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=(), n_steps=10, x0=nothing,
                            number_of_batches=(0, 0), strictness=:warn)
CellListMapNeighborFinder(; eligible, dist_cutoff, boundary, special=nothing,
                            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 only updated every n_steps steps. 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=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, use_neighbors, weight_special, 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, use_neighbors=false, weight_special=1,
             coulomb_const=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 and error_tol should match these interactions.

The potential energy is defined as

\[V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0} \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, and $\epsilon_{\mathrm{tol}}$ is error_tol. 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=false,
                   λ_mixing=MinimumMixing(), scheduler=Molly.DefaultLambdaScheduler(),
                   weight_special=1, coulomb_const=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 Csoftcore=: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, use_neighbors, weight_special,
                        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} \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} - 1}{2\varepsilon_\mathrm{rf} + 1}, \quad c_\mathrm{rf} = \frac{1}{r_\mathrm{c}} \frac{3\varepsilon_\mathrm{rf}}{2\varepsilon_\mathrm{rf} + 1}\]

solvent_dielectric corresponds to $\varepsilon_\mathrm{rf}$. Setting solvent_dielectric=Inf gives conducting boundary conditions ($k_\mathrm{rf} = 1/(2r_\mathrm{c}^3)$, $c_\mathrm{rf} = 3/(2r_\mathrm{c})$).

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

CoulombReactionFieldScaled(; dist_cutoff, solvent_dielectric, use_neighbors,
                           λ_mixing, scheduler, weight_special, 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 Csoftcore=: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, use_neighbors, λ_mixing, scheduler, weight_special,
              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 Csoftcore=: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, α, use_neighbors, σ_mixing, ϵ_mixing, λ_mixing,
                       scheduler, weight_special, 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 Csoftcore=: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=false, σ_mixing=LorentzMixing(),
                            ϵ_mixing=GeometricMixing(), λ_mixing=MinimumMixing(),
                            scheduler=Molly.DefaultLambdaScheduler(), weight_special=1,
                            coulomb_const=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 Csoftcore=: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, α, use_neighbors,
                                     σ_mixing, ϵ_mixing, λ_mixing, scheduler,
                                     weight_special, 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 Csoftcore=: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.

CoulombSoftCoreGapsys(; cutoff, α, σQ, use_neighbors, λ_mixing, scheduler,
                      weight_special, 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 Csoftcore=: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=false, λ_mixing=MinimumMixing(),
                           scheduler=Molly.DefaultLambdaScheduler(), weight_special=1,
                           coulomb_const=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 Csoftcore=: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, α, σQ,
                                    use_neighbors, λ_mixing, scheduler,
                                    weight_special, 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 Csoftcore=: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.

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}\]

DPDInteraction(; a, γ, σ, r_c, dt, use_neighbors, key)

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.

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=(),
                       n_steps=10, array_type=nothing, strictness=:warn)
DistanceNeighborFinder(; eligible, dist_cutoff, special=nothing, 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 only updated every n_steps steps.

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, use_neighbors, α, β, shortcut, σ_mixing, ϵ_mixing,
                  weight_special)

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, use_neighbors, α, β, shortcut, σ_mixing,
                          ϵ_mixing, λ_mixing, scheduler, weight_special)

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.

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.

Ewald(dist_cutoff; error_tol=0.0005, 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 and error_tol 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. 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 and error_tol 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=(),
                          n_steps=10, array_type, max_neighbors=nothing,
                          output=:molly_pairs, ragged=true, strictness=:warn)
GPUCellListNeighborFinder(; eligible, dist_cutoff, special=nothing, n_steps=10,
                          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 dist_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 only updated every n_steps steps.

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=(),
                  n_steps=10, array_type, initialized=false, strictness=:warn)
GPUNeighborFinder(; eligible, dist_cutoff, special=nothing, 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.

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 only refreshed every n_steps steps. The buffer should be larger than the distance an atom can move in that time.
• 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.
• n_steps: the number of steps between Morton reorder and tile-list refresh passes.
• 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.

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.

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, G, use_neighbors)

The gravitational interaction between two atoms.

The potential energy is defined as

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

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=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=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.

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=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=gb_solvent_dielectric,
                    solute_dielectric=gb_solute_dielectric, kappa=0.0u"nm^-1",
                    offset=gbn2_offset, dist_cutoff=0.0u"nm", probe_radius=gb_probe_radius,
                    sa_factor=gb_sa_factor, use_ACE=true, neck_scale=gbn2_neck_scale,
                    neck_cut=gbn2_neck_cut, element_to_radius=mbondi2_element_to_radius,
                    element_to_screen=gbn2_element_to_screen,
                    element_to_screen_nucleic=gbn2_element_to_screen_nucleic,
                    atom_params=gbn2_atom_params,
                    atom_params_nucleic=gbn2_atom_params_nucleic, data_d0=gbn2_data_d0,
                    data_m0=gbn2_data_m0, 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. 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=gb_solvent_dielectric,
                   solute_dielectric=gb_solute_dielectric, kappa=0.0u"nm^-1",
                   offset=obc_offset, dist_cutoff=0.0u"nm", probe_radius=gb_probe_radius,
                   sa_factor=gb_sa_factor, use_ACE=true, use_OBC2=false,
                   element_to_radius=mbondi2_element_to_radius,
                   element_to_screen=obc_element_to_screen, 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. 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, dist_cutoff, σ_mixing=LorentzMixing(),
                       ϵ_mixing=GeometricMixing())

The long-range dispersion correction for the LennardJones interaction.

See Shirts et al. 2007. Approximately represents contributions from beyond the cutoff distance. Should be used alongside the LennardJones pairwise interaction when the long-range correction to the potential energy is required. 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, $V$ is the volume and $r_c$ is dist_cutoff. The forces are zero, with an isotropic virial contribution from the explicit volume dependence.

The number of atoms and atom σ and ϵ values are assumed not to change after setup (the box volume can change). Only compatible with 3D systems. Not compatible with cutoffs other than DistanceCutoff.

LJDispersionCorrectionλ(atoms, dist_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. 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)\]

where 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 and atom σ and ϵ values are assumed not to change after setup (the box volume can change). Only compatible with 3D systems. Not compatible with cutoffs other than DistanceCutoff.

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=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. This is a leapfrog-like integrator, so the velocities are 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.

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, use_neighbors, shortcut, σ_mixing, ϵ_mixing, weight_special)

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, scheduler)

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 LJsoftcore=: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, α, λ_mixing, scheduler)

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 LJsoftcore=: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, α, λ_mixing, scheduler)

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 LJsoftcore=: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, use_neighbors, shortcut, σ_mixing, ϵ_mixing, λ_mixing,
                   scheduler, weight_special)

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 LJsoftcore=: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, α, use_neighbors, shortcut, σ_mixing,
                            ϵ_mixing, λ_mixing, scheduler, weight_special)

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 LJsoftcore=: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, α, use_neighbors, shortcut, σ_mixing,
                           ϵ_mixing, λ_mixing, scheduler, weight_special)

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 LJsoftcore=: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.

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.

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.

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 n_steps argument of the neighbor finder is 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.

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 n_steps argument of the neighbor finder is ignored.

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.

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, use_neighbors, shortcut, σ_mixing, ϵ_mixing, weight_special)

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. 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=NoNeighborFinder(): the neighbor finder used to find close atoms and save on computation.
• 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. The neighbor list is not updated when making trial moves or after accepted moves. 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()
MonteCarloLogger(T)

A logger that records acceptances in a Monte Carlo simulation.

The logged quantities include the number of new selections (n_select), 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=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, a, b, c, x0, y0, force_units, energy_units)

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.

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.

NoNeighborFinder()

Placeholder neighbor finder that returns no neighbors.

When using this neighbor finder, ensure that use_neighbors for the interactions returns false.

NoseHoover(; <keyword arguments>)

The Nosé-Hoover integrator, a NVT simulator that extends velocity Verlet to control the temperature of the system.

See Evans and Holian 1985. The current implementation is limited to ergodic systems.

ζ, the thermostat friction variable, is reset on each call to simulate!. This means that two calls to simulate! will not be the same as one longer call. Not compatible with constraints, will print a warning and continue without applying constraints.

Arguments

• dt::T: the time step of the simulation.
• temperature::K: the equilibrium temperature of the simulation.
• damping::D=100*dt: the temperature damping time scale.
• 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.

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 LJsoftcore=:gapsys and Csoftcore=: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 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, 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 and error_tol 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. order is the order of the B-spline interpolation onto the grid and ϵr is the relative permittivity. 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 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,
      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.

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.

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)

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, λ_mixing=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.

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=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.

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.

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.

SetupCoulombReactionField(; solvent_dielectric, coulomb_const)

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.

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=gb_solvent_dielectric,
                         solute_dielectric=gb_solute_dielectric, kappa=0.0u"nm^-1",
                         offset=gbn2_offset, dist_cutoff=0.0u"nm",
                         probe_radius=gb_probe_radius, sa_factor=gb_sa_factor,
                         use_ACE=true, neck_scale=gbn2_neck_scale,
                         neck_cut=gbn2_neck_cut,
                         element_to_radius=mbondi2_element_to_radius,
                         element_to_screen=gbn2_element_to_screen,
                         element_to_screen_nucleic=gbn2_element_to_screen_nucleic,
                         atom_params=gbn2_atom_params,
                         atom_params_nucleic=gbn2_atom_params_nucleic,
                         data_d0=gbn2_data_d0, data_m0=gbn2_data_m0)

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.

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

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.

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,
         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. 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.

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, use_neighbors, shortcut, σ_mixing, ϵ_mixing)

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.

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=NoNeighborFinder(): the neighbor finder used to find close atoms and save on computation.
• 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=0.2u"nm": distance added to dist_cutoff when calculating classical neighbor lists every few steps.
• 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), or a cutoff like DistanceCutoff (short range only). For a cutoff, the cutoff distance must not exceed dist_cutoff.
• 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 Lennard-Jones dispersion correction is being used, this must be a DistanceCutoff.
• dispersion_correction=nothing: whether to use the long-range Lennard-Jones dispersion correction, LJDispersionCorrection. Defaults to the force field setting, which defaults to true.
• 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: which neighbor finder to use, default is 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. NoNeighborFinder can be used but in this case bonded atoms will not be excluded from the non-bonded interactions.
• neighbor_finder_n_steps=10: 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.

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(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.

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.

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.

TreeNeighborFinder(; n_atoms, dist_cutoff, excluded_pairs=(), special_pairs=(),
                   n_steps=10, strictness=:warn)
TreeNeighborFinder(; eligible, dist_cutoff, special=nothing, 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 only updated every n_steps steps. 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=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.

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.

Yukawa(; cutoff, use_neighbors, weight_special, coulomb_const, kappa)

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.

values(logger)
values(logger::TimeCorrelationLogger; normalize::Bool=true)
values(logger::AverageObservableLogger; std::Bool=true)

Access the stored observations in a logger.

AbsoluteFESystem(sys, global_λ, mapping; temp=298.0u"K",
                 scheduler=Molly.DefaultLambdaScheduler(dual=true), loggers=(),
                 LJsoftcore=:gapsys, Csoftcore=: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.
• LJsoftcore=: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.
• Csoftcore=:gapsys: the Coulomb soft core, :beutler, :gapsys or :scaled, where the potential is scaled directly by λ.

AutoCorrelationLogger(observable::Function, TA::DataType,
                        observable_length::Integer, n_correlation::Integer)

An autocorrelation logger, equivalent to a TimeCorrelationLogger in the case that observableA == observableB.

BoxLogger(n_steps)
BoxLogger(T, n_steps)

Log the boundary of a system throughout a simulation, represented as its box matrix. Useful for NPT simulations where the box dimensions fluctuate.

CoordinatesLogger(n_steps; dims=3)
CoordinatesLogger(T, n_steps; dims=3)

Log the coordinates throughout a simulation.

DensityLogger(n_steps)
DensityLogger(T, n_steps)

Log the density of a system throughout a simulation.

Not compatible with infinite boundaries.

ForcesLogger(n_steps; dims=3)
ForcesLogger(T, n_steps; dims=3)

Log the forces throughout a simulation.

The forces are those from the interactions and do not include forces applied by stochastic simulators such as Langevin.

KineticEnergyLogger(n_steps)
KineticEnergyLogger(T, n_steps)

Log the kinetic_energy of a system throughout a simulation.

LocalCoordinatesSite(atom_ind, atom_1, atom_2, atom_3, origin_weights, x_weights,
                     y_weights, local_position)

A virtual site defined by a position in a local coordinate system given by three atoms, matching the LocalCoordinatesSite of OpenMM.

Returns a VirtualSite defined by:

\[\begin{aligned} \mathbf{o} &= \sum_i w_i^o \mathbf{r}_i & \mathbf{x} &= \sum_i w_i^x \mathbf{r}_i & \mathbf{y} &= \sum_i w_i^y \mathbf{r}_i \\ \hat{\mathbf{x}} &= \frac{\mathbf{x}}{|\mathbf{x}|} & \hat{\mathbf{z}} &= \frac{\mathbf{x} \times \mathbf{y}}{|\mathbf{x} \times \mathbf{y}|} & \hat{\mathbf{y}} &= \hat{\mathbf{z}} \times \hat{\mathbf{x}} \end{aligned}\]

\[\mathbf{r} = \mathbf{o} + p_1 \hat{\mathbf{x}} + p_2 \hat{\mathbf{y}} + p_3 \hat{\mathbf{z}}\]

where $\mathbf{p}$ is local_position, $\sum_i w_i^o$ must equal 1 and $\sum_i w_i^x$ and $\sum_i w_i^y$ must equal 0.

Since the axes follow the atoms rather than the box, this places sites at a fixed distance and orientation, such as the lone pairs of a CHARMM force field.

Only compatible with 3D systems. Not compatible with virial calculation. Not compatible with gradient calculation using Enzyme.

OneParticleSite(atom_ind, atom_1)

A virtual site defined to have the same coordinates as another atom.

Returns a VirtualSite defined by:

\[\mathbf{r} = \mathbf{r}_1\]

This can be useful in alchemical simulations when multiple versions of an atom are required.

Not compatible with gradient calculation using Enzyme.

OutOfPlaneSite(atom_ind, atom_1, atom_2, atom_3, weight_12, weight_13, weight_cross)

A virtual site defined by the weighted average of the coordinates of three atoms and the cross product of their relative displacements.

Returns a VirtualSite defined by:

\[\mathbf{r} = \mathbf{r}_1 + w_{12} \mathbf{r}_{12} + w_{13} \mathbf{r}_{13} + w_{\mathrm{cross}} (\mathbf{r}_{12} \times \mathbf{r}_{13})\]

Only compatible with 3D systems. Not compatible with virial calculation. Not compatible with gradient calculation using Enzyme.

PotentialEnergyLogger(n_steps)
PotentialEnergyLogger(T, n_steps)

Log the potential_energy of a system throughout a simulation.

PressureLogger(n_steps)
PressureLogger(T, n_steps)

Log the pressure tensor of a system throughout a simulation.

For constrained systems on coordinate-scaling coupling steps, this records the pre-coupling pressure used for coupling, not a pressure recomputed with the post-coupling box.

This should only be used on 3-dimensional systems.

RelativeFESystem(sysA, sysB, global_λ, mapping, core_mapAB; temp=298.0u"K",
                 scheduler=Molly.DefaultLambdaScheduler(dual=true), loggers=(),
                 LJsoftcore=:gapsys, Csoftcore=:gapsys)

Set up a hybrid system for a relative free energy calculation that transforms system A (global_λ = 0) into system B (global_λ = 1).

Core atoms are interpolated between their parameters in A and B, unique A atoms are decoupled and unique B atoms are coupled. The environment, every atom of sysA that is neither core nor unique A, is taken from sysA. The atoms of the returned system are ordered as core, unique A, unique B and environment. In dual topology each core atom is added twice, with the parameters of A and as a massless copy with the parameters of B. Interactions are replaced by their λ-dependent versions and random velocities are generated at temp.

Arguments

• sysA: the System of state A, which also provides the environment.
• sysB: the System of state B, which provides the unique B atoms and the B parameters of the core atoms.
• global_λ: the global λ of the system, from 0 (state A) to 1 (state B).
• mapping: a dictionary of atom index vectors, where "core" and "unique_A" index into sysA and "unique_B" into sysB. An "env" entry is added to it. An atom unique to one end state cannot be bonded to an atom outside the mapping, as it is switched off at the other end state; such an atom belongs in "core", where only its parameters change, which is how a covalently bound ligand is set up.
• core_mapAB: a dictionary from each core atom index in sysA to its index in sysB.
• 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.GROMACSLambdaRBFEScheduler or Molly.OpenFEScheduler. dual=true uses dual topology (energy scaling) and dual=false single topology (parameter scaling), which does not support intraLJ=true.
• loggers=(): the loggers that record properties of interest during a simulation.
• LJsoftcore=: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.
• Csoftcore=:gapsys: the Coulomb soft core, :beutler, :gapsys or :scaled, where the potential is scaled directly by λ.

ScalarPressureLogger(n_steps)
ScalarPressureLogger(T, n_steps)

Log the scalar_pressure of a system throughout a simulation.

For constrained systems on coordinate-scaling coupling steps, this records the pre-coupling pressure used for coupling, not a pressure recomputed with the post-coupling box.

This should only be used on 3-dimensional systems.

ScalarVirialLogger(n_steps)
ScalarVirialLogger(T, n_steps)

Log the scalar_virial of a system throughout a simulation.

For constrained systems on coordinate-scaling coupling steps, this records the pre-coupling scalar virial for that step.

TemperatureLogger(n_steps)
TemperatureLogger(T, n_steps)

Log the temperature throughout a simulation.

ThreeParticleAverageSite(atom_ind, atom_1, atom_2, atom_3, weight_1, weight_2, weight_3)

A virtual site defined by the weighted average of the coordinates of three atoms.

Returns a VirtualSite defined by:

\[\mathbf{r} = w_1 \mathbf{r}_1 + w_2 \mathbf{r}_2 + w_3 \mathbf{r}_3\]

where $w_1 + w_2 + w_3$ must equal 1.

Not compatible with gradient calculation using Enzyme.

TotalEnergyLogger(n_steps)
TotalEnergyLogger(T, n_steps)

Log the total_energy of a system throughout a simulation.

TwoParticleAverageSite(atom_ind, atom_1, atom_2, weight_1, weight_2)

A virtual site defined by the weighted average of the coordinates of two atoms.

Returns a VirtualSite defined by:

\[\mathbf{r} = w_1 \mathbf{r}_1 + w_2 \mathbf{r}_2\]

where $w_1 + w_2$ must equal 1.

Not compatible with gradient calculation using Enzyme.

VelocitiesLogger(n_steps; dims=3)
VelocitiesLogger(T, n_steps; dims=3)

Log the velocities throughout a simulation.

VirialLogger(n_steps)
VirialLogger(T, n_steps)

Log the virial tensor of a system throughout a simulation.

For constrained systems on coordinate-scaling coupling steps, this records the pre-coupling virial for that step.

VolumeLogger(n_steps)
VolumeLogger(T, n_steps)

Log the volume of a system throughout a simulation.

Not compatible with infinite boundaries.

accelerations(system, neighbors=find_neighbors(system), step_n=0;
              n_threads=Threads.nthreads(), pairwise_inters=system.pairwise_inters,
              specific_inter_lists=system.specific_inter_lists,
              general_inters=system.general_inters, strictness=:warn)

Calculate the accelerations of all atoms in a system using the pairwise, specific and general interactions and Newton's second law of motion.

add_position_restraints(sys, k; atom_selector=Molly.is_any_atom,
                        restrain_coords=sys.coords)

Return a copy of a System with HarmonicPositionRestraints added to restrain the atoms.

The force constant k can be a single value or an array of equal length to the number of atoms in the system. The atom_selector function takes in each atom and atom data and determines whether to restrain that atom. For example, Molly.is_heavy_atom means non-hydrogen atoms are restrained.

ani2x_data_dir()

Path to the ANI-2x data directory: a lazily-downloaded artifact holding ani2x.h5 (the model weights) and 6mrr_ani2x.json (TorchANI reference energies). Requires Lux and HDF5. Load the potential with ANIPotential(joinpath(Molly.ani2x_data_dir(), "ani2x.h5")).

apply_coupling!(system, buffers, coupling, simulator, neighbors=nothing, step_n=0;
                n_threads=Threads.nthreads(), rng=Random.default_rng(),
                strictness=:warn)

Apply a coupler to modify a simulation.

Returns whether the coupling has invalidated the currently stored forces, for example by changing the coordinates. This information is useful for some simulators. If coupling is a tuple or named tuple then each coupler will be applied in turn. Custom couplers should implement this function.

apply_loggers!(system, neighbors=nothing, step_n=0, buffers=nothing, run_loggers=true;
               n_threads=Threads.nthreads(), strictness=:warn, kwargs...)

Run the loggers associated with a system.

run_loggers can be true, false or :skipstart, in which case the loggers are not run before the first step. Additional keyword arguments can be passed to the loggers if required.

apply_position_constraints!(sys, coord_storage; context=nothing,
                            n_threads=Threads.nthreads(), strictness=:warn)
apply_position_constraints!(sys, coord_storage, vel_storage, dt; context=nothing,
                            n_threads=Threads.nthreads(), strictness=:warn)

Apply the coordinate constraints to the system.

If vel_storage and dt are provided then velocity constraints are applied as well.

apply_velocity_constraints!(sys; context=nothing, n_threads=Threads.nthreads(),
                            strictness=:warn)

Apply the velocity constraints to the system.

array_type(sys)
array_type(arr)

The array type of a System, ReplicaSystem or array, for example Array for systems on CPU or CuArray for systems on a NVIDIA GPU.

assemble_mbar_inputs(coords_k, boundaries_k, states;
                     target_state=nothing, shift=false)

Assemble the reduced potentials matrix u (size N×K) for MBAR from per-window coordinates and boundaries.

The reduced potential of frame $n$ in state $k$ is

\[u_k(x_n) = \beta_k \left( U_k(x_n) + p_k V(x_n) \right)\]

where $\beta_k$ and $p_k$ are the inverse temperature and pressure of the ThermoState, and the pressure term is omitted for states without a pressure. Interactions shared between the states are evaluated once per frame where possible.

Arguments

• coords_k::Vector{<:Vector} - subsampled coordinates for each window k.
• boundaries_k::Vector{<:Vector} - subsampled boundaries for each window k (same lengths as coords_k[k]).
• states::AbstractVector{<:ThermoState} - thermodynamic states for each window.

Keyword arguments

• target_state::Union{Nothing, ThermoState} - if set, also compute u_target for that state.
• shift::Bool=false - if true, subtract per-frame row minima from u and return the shifts.

Returns

An MBARInput struct, which can be given to Molly.mbar_weights, with:

• u::Matrix{Float64} - N×K reduced potentials.
• u_target::Union{Vector{Float64}, Nothing} - reduced potential at target_state or nothing.
• N::Vector{Int} - sample counts per window.
• win_of::Vector{Int} - window index for each frame.
• shifts::Union{Vector{Float64},Nothing} - per-frame shifts when shift=true, else nothing.

bias_gradient(bias_type, cv_sim)

Calculate the gradient of a bias potential with respect to the value of a collective variable.

Custom bias types should implement this function along with potential_energy.

Arguments

• bias_type: a struct that defines the type of bias to be used, such as LinearBias.
• cv_sim: the value of a measured collective variable given the coordinates of a simulation.

bond_angle(coord_i, coord_j, coord_k, boundary)
bond_angle(vec_ji, vec_jk)

Calculate the bond or pseudo-bond angle in radians between three coordinates or two vectors.

The angle between j→i and j→k is returned in the range 0 to π.

box_center(boundary)

Calculate the center of a bounding box.

Dimensions with infinite length return zero.

calculate_cv(cv, coords, atoms, boundary, velocities; kwargs...)

Calculate the value of a collective variable (CV) with the current system state.

New CV types should implement this function. This function does not apply the molecule correction over the boundaries; BiasPotential makes the molecules whole beforehand when the correction field of the CV is :pbc. The gradient of this function with respect to coordinates, used to calculate forces, is by default calculated with automatic differentiation when Enzyme is imported. Alternatively, the Molly.cv_gradient function can be defined for a new CV type.

charge(atom)

The partial charge of an Atom.

Custom atom types should implement this function if charges are going to be used unless they have a charge field defined, which the function accesses by default.

charges(sys)

The partial charges of the atoms in a System or ReplicaSystem.

check_constraints(sys)
check_constraints(sys, constraints)

Check whether the coordinates and velocities of a system satisfy the coordinate and velocity constraints.

check_position_constraints(sys)
check_position_constraints(sys, constraints)

Check whether the coordinates of a system satisfy the position constraints.

check_velocity_constraints(sys)
check_velocity_constraints(sys, constraints)

Check whether the velocities of a system satisfy the velocity constraints.

compute_aevs(coords, species_indices, neighbors, boundary, aev_params, n_species)

Compute the Atomic Environment Vectors (AEVs) for all atoms, returning an (n_atoms, aev_length) matrix, as used by ANIPotential. neighbors is a NeighborList (or nothing for an all-pairs build). Requires Lux and HDF5.

cv_gradient(cv, coords, atoms, boundary, velocities; kwargs...)

Calculate the analytical gradient of a collective variable (CV) with respect to the system coordinates.

Returns a tuple containing the gradient (as an array of vectors) and the current value of the CV. When Enzyme is imported this defaults to using AD, but an explicit method can be provided for a given CV type and is defined for built-in CVs. The AD approach should work with and without units.

Explicit methods are defined for:

• CalcDist: distance between two atoms or groups (minimum, maximum, center of mass or single).
• CalcAngle: angle defined by three atoms.
• CalcRg: radius of gyration of a group of atoms.
• CalcRMSD: root-mean-square deviation from a reference structure using Kabsch alignment.
• CalcTorsion: torsion (dihedral) angle defined by four atoms.

density(sys)

The density of a System.

Returns zero density for infinite boundaries.

dipole_moment(sys)

The dipole moment of a system.

Calculated as the sum of the coordinates multiplied by the partial charges, $\boldsymbol{\mu} = \sum_i q_i \mathbf{r}_i$. Requires the charges on the atoms to be set.

displacements(coords, boundary)

Calculate the pairwise vector displacements of a set of coordinates, accounting for the periodic boundary conditions.

distances(coords, boundary)

Calculate the pairwise distances of a set of coordinates, accounting for the periodic boundary conditions.

extract_parameters(sys, ff=nothing)

A dictionary of the force field parameters of a system, keyed by strings, for use with inject_gradients.

An interaction contributes its parameters if it defines Molly.parameter_prefix and Molly.parameter_fields.

The values are plain numbers of the float type of the system. On a system with units they are in the units that each parameter is stored with in the system, for example kJ mol⁻¹ nm⁻² for a bond force constant in a system set up with the default units, and the units are re-attached by ParameterPlan and inject_gradients.

find_neighbors(system; n_threads=Threads.nthreads())
find_neighbors(system, neighbor_finder, current_neighbors=nothing, step_n=0,
               force_recompute=false; n_threads=Threads.nthreads())

Obtain a list of close atoms in a System.

Custom neighbor finders should implement this function. The neighbor list it returns should define Base.length, giving the number of pairs, and neighbor_pairs, giving the (i, j, special) pairs themselves, since that is what the pairwise force and energy loops use.

Returns a NeighborList for the classical neighbor finders and a GPUCellListNeighborList for GPUCellListNeighborFinder.

For GPUNeighborFinder, this returns nothing: the CUDA pairwise force and energy kernels build and cache their interacting tile list internally from the neighbor-finder metadata.

A neighbor finder may reuse the device buffers behind current_neighbors for the list it returns, as GPUCellListNeighborFinder does, so the list passed in should not be used afterwards.

float_type(sys)
float_type(boundary)

The float type a System, ReplicaSystem or bounding box uses.

Set with the float_type argument when constructing a system. When setting up a system from a file, the default is Float64 on CPU and Float32 on GPU. When setting up a system directly, the default is the float type of the boundary. See float_type_high for the higher precision accumulation type.

float_type_high(sys)

The float type a System or ReplicaSystem uses for accumulation where higher precision is useful, such as the potential_energy and the virial.

Set with the float_type_high argument when constructing a system, Float64 by default. See float_type for the type used for the coordinates, velocities and interaction parameters.

force(inter, vec_ij, atom_i, atom_j, force_units, special, coord_i, coord_j,
      boundary, velocity_i, velocity_j, step_n)
force(inter, coord_i, boundary, atom_i, force_units, velocity_i, step_n, data)
force(inter, coord_i, coord_j, boundary, atom_i, atom_j, force_units, velocity_i,
      velocity_j, step_n, data)
force(inter, coord_i, coord_j, coord_k, boundary, atom_i, atom_j, atom_k,
      force_units, velocity_i, velocity_j, velocity_k, step_n, data)
force(inter, coord_i, coord_j, coord_k, coord_l, boundary, atom_i, atom_j, atom_k,
      atom_l, force_units, velocity_i, velocity_j, velocity_k, velocity_l, step_n, data)
force(inter, coord_i, coord_j, coord_k, coord_l, coord_m, boundary, atom_i, atom_j,
      atom_k, atom_l, atom_m, force_units, velocity_i, velocity_j, velocity_k,
      velocity_l, velocity_m, step_n, data)

Calculate the force between atoms due to a given interaction type.

For pairwise interactions returns a single force vector and for specific interactions returns a type such as SpecificForce2Atoms. Custom pairwise and specific interaction types should implement this function. Pairwise interactions can use Molly.pairwise_force to apply a cutoff.

forces(system, neighbors=find_neighbors(system), step_n=0;
       n_threads=Threads.nthreads(), pairwise_inters=system.pairwise_inters,
       specific_inter_lists=system.specific_inter_lists,
       general_inters=system.general_inters, strictness=:warn)

Calculate the forces on all atoms in a system using the pairwise, specific and general interactions.

forces_virial(system, neighbors=find_neighbors(system), step_n=0;
              n_threads=Threads.nthreads(), pairwise_inters=system.pairwise_inters,
              specific_inter_lists=system.specific_inter_lists,
              general_inters=system.general_inters, strictness=:warn)

Calculate the forces on all atoms in a system and the virial using the pairwise, specific and general interactions. For constrained systems, constraint virial contributions are approximated using the same deterministic small-step constraint preview as virial.

Returns a tuple of the forces and the virial. This is faster than calling forces and virial separately.

has_ragged_neighbors(neighbors)

Whether a neighbor list also stores its neighbors per atom, default false.

See ragged_neighbors.

hydrodynamic_radius(coords, boundary)

Calculate the hydrodynamic radius of a set of coordinates.

$R_{hyd}$ is defined by

\[\frac{1}{R_{hyd}} = \frac{1}{2N^2}\sum_{i \neq j} \frac{1}{r_{ij}}\]

inject_gradients(sys, params_dic)
inject_gradients(sys, params_dic, plan)
inject_gradients(sys, params_dic, plan, coords)

A copy of sys with the parameters in params_dic injected, for use inside a function being differentiated with respect to params_dic.

plan is a ParameterPlan resolving the dictionary keys to the parameters of the system. Build it once outside the function being differentiated and pass it in; the two-argument form builds one on every call, which is convenient but much slower. coords replaces the coordinates of sys, so that a gradient can also be taken with respect to them.

sys may have units, in which case the values in params_dic are plain numbers in the units that each parameter is stored with in sys, as returned by extract_parameters. The unit of each field is re-attached when it is replaced. Gradients are only supported for systems without units.

is_any_atom(at, at_data)

Placeholder function that returns true, used to select any Atom.

is_heavy_atom(at, at_data)

Determine whether an Atom is a heavy atom, i.e. any element other than hydrogen.

is_on_gpu(sys)
is_on_gpu(arr)

Whether a System, ReplicaSystem or array type is on the GPU.

iterate_mbar(u, win_of, N_counts; rtol=1e-8, max_iter=10_000)

Solve the MBAR self-consistent equations, see Shirts and Chodera 2008 Eq C3.

The relative free energies satisfy

\[f_i = -\ln \sum_{n=1}^{N} \frac{\exp(-u_i(x_n))}{\sum_{k=1}^{K} N_k \exp(f_k - u_k(x_n))}\]

where the sum over $n$ runs over the frames from all states.

Arguments

• u::AbstractMatrix - N×K reduced potentials (rows = frames, cols = states).
• win_of::AbstractVector - length-N vector giving the generating state index per frame.
• N_counts::AbstractVector - length-K sample counts per state.

Keyword arguments

• rtol::Float64=1e-8 - relative convergence tolerance.
• max_iter::Int=10_000 - maximum iterations.

Returns

• f::Vector{Float64} - relative free energies per state (gauge-fixed to f[1]=0).
• logN::Vector{Float64} - log.(N_counts) for reuse downstream.

kinetic_energy(system; kin_tensor=nothing)

Calculate the kinetic energy of a system.

The scalar kinetic energy is defined as

\[K = \rm{Tr}\left[ \bf{K} \right]\]

where $\bf{K}$ is the kinetic energy tensor:

\[\bf{K} = \frac{1}{2} \sum_{i} m_i \bf{v_i} \otimes \bf{v_i}\]

kinetic_energy_tensor(system; kin_tensor=nothing)

Calculate the kinetic energy of a system in its tensorial form.

The kinetic energy tensor is defined as

\[\bf{K} = \frac{1}{2} \sum_{i} m_i \bf{v_i} \otimes \bf{v_i}\]

where $m_i$ is the mass and $\bf{v_i}$ is the velocity vector of atom $i$. See kinetic_energy for the scalar kinetic energy.

lambda(atom)

The λ of an Atom.

log_property!(logger, system, neighbors=nothing, step_n=0, buffers=nothing;
              n_threads=Threads.nthreads(), strictness=:warn, kwargs...)

Log a property of a system throughout a simulation.

Custom loggers should implement this function. Additional keyword arguments can be passed to the logger if required.

log_tss_sample!(logger, sample::TSSSample)

Record a TSS sample. Called by TSSSimulation once per collected sample, and for the initial configuration, for each logger in sample_loggers or replica_sample_loggers.

Custom sample loggers should implement this function. As for system loggers, the logger decides which samples to record, for example from sample.step; see TSSSampleLogger for a logger with a step interval.

logger_pressure_interval(logger)

Return the interval in simulation steps at which logger needs a valid pressure.

Custom pressure loggers should extend this function in addition to Molly.logger_virial_interval. The default Inf means that the logger does not require pressure-specific pre-coupling state.

logger_virial_interval(logger)

Return the interval in simulation steps at which logger needs a valid virial.

Custom loggers that need the total virial should extend this function. The default Inf means that the logger does not require virial calculation.

mass(atom)

The mass of an Atom.

Custom atom types should implement this function unless they have a mass field defined, which the function accesses by default. Virtual sites should have zero mass, and non-virtual sites should have non-zero mass.

masses(sys)

The masses of the atoms in a System or ReplicaSystem.

maxwell_boltzmann(atom_mass::Unitful.Mass, temp::Unitful.Temperature,
                  k::BoltzmannConstUnits=Unitful.k; rng=Random.default_rng())
maxwell_boltzmann(atom_mass::MolarMass, temp::Unitful.Temperature,
                  k_molar::MolarBoltzmannConstUnits=(Unitful.k * Unitful.Na);
                  rng=Random.default_rng())
maxwell_boltzmann(atom_mass::Real, temperature::Real,
                  k::Real=ustrip(u"u * nm^2 * ps^-2 * K^-1", Unitful.k);
                  rng=Random.default_rng())

Generate a random velocity along one dimension from the Maxwell-Boltzmann distribution, with optional custom Boltzmann constant.

mbar_pmf(u, u_target, f, N_counts, logN, observables, grid;
         shifts=nothing, zero=:min, kBT=nothing)

Compute a target-state PMF from solved MBAR quantities using the shared log-space PMF accumulator. Unlike Molly.pmf_with_uncertainty, this routine supports multi-dimensional observables and arbitrary grids, but it does not compute the large-sample covariance estimate.

mbar_weights(u, u_target, f, N_counts, logN; shifts=nothing, check=true)

Compute MBAR weights for sampled states and for a target state, see Shirts and Chodera 2008 Eq 13.

Arguments

• u::AbstractMatrix - N×K reduced potentials for sampled states.
• u_target::AbstractVector - length-N reduced potential for the target state.
• f::AbstractVector - length-K relative free energies.
• N_counts::AbstractVector - length-K sample counts per state.
• logN::AbstractVector - log.(N_counts).

Keyword arguments

• shifts::Union{Nothing, AbstractVector}=nothing - per-frame shifts previously subtracted from u, if any.
• check::Bool=true - perform basic normalization checks.

Returns

• W::Matrix{Float64} - N×K sampled-state weights.
• w_target::Vector{Float64} - length-N target-state weights.
• logD::Vector{Float64} - length-N log normalizers used in the weights.

mbar_weights(mbar_generator::MBARInput)

High-level MBAR wrapper that computes free energies and reweighting weights from a pre-assembled MBARInput struct.

Arguments

mbar_generator::MBARInput - result from Molly.assemble_mbar_inputs containing:

• u::AbstractMatrix - reduced potentials (N×K).
• u_target::Union{AbstractVector, Nothing} - reduced potentials at the target state.
• N::AbstractVector - sample counts per state.
• win_of::AbstractVector - window index for each frame.
• shifts::Union{AbstractVector,Nothing} - per-frame energy shifts (optional).

Returns

(W, w_target, logD) where:

• W::Matrix{Float64} - N×K sampled-state weights.
• w_target::Vector{Float64} - target-state weights.
• logD::Vector{Float64} - per-frame log normalizers.

This routine runs Molly.iterate_mbar to obtain relative free energies F_k and then calls the lower-level mbar_weights(u, u_target, f, N_counts, logN) using the contents of mbar_generator. All internal consistency checks are disabled for speed.

momentum(system)

The total momentum of a System.

This is the sum over atoms of the velocity multiplied by the mass, and is close to zero when the center of mass motion is removed with remove_CM_motion!.

neighbor_pairs(neighbors)

The (i, j, special) pairs of a neighbor list, as returned by find_neighbors.

This is what the pairwise force and energy loops iterate over, so a neighbor list type from a custom neighbor finder should define this along with Base.length. The result indexes into the list rather than copying it, so it should not be used after the list is passed back to find_neighbors, which may overwrite the pairs.

optimize_cuda_launch_config!(system)

Autotune and cache CUDA launch overrides for the given system.

On CUDA systems, this benchmarks a small launch-candidate set for Molly's tiled pairwise kernels, caches the result for the current device/system signature, and writes the winning values into system.launch_config. A fallback is provided for non-CUDA systems that does nothing. This is called when setting up a System from a file unless autotune_launch=false is given.

pairwise_force(inter, r, params)

Calculate the force magnitude between two atoms separated by distance r due to a pairwise interaction.

A positive force causes the atoms to move apart. This function is used in force to apply cutoff strategies by calculating the force at different values of r. Consequently, the parameters params should not include terms that depend on distance.

pairwise_pe(inter, r, params)

Calculate the potential energy between two atoms separated by distance r due to a pairwise interaction.

This function is used in potential_energy to apply cutoff strategies by calculating the potential energy at different values of r. Consequently, the parameters params should not include terms that depend on distance.

parameter_fields(T)

A tuple of (field, key_suffix) pairs naming the differentiable parameters of an interaction type, for example ((:k, "k"), (:r0, "r0")).

The key of a parameter is Molly.parameter_prefix followed by its suffix. Interactions whose parameters do not correspond one-to-one with fields, such as PeriodicTorsion, instead define Molly.n_parameters, Molly.parameter_keys, Molly.parameter_values and Molly.inject_parameters directly.

parameter_prefix(inter)
parameter_prefix(inter, inter_type)

The prefix that the parameters of an interaction have in a parameter dictionary, for example "inter_HB_CT/CT_".

Return nothing, the default, for an interaction with no differentiable parameters. Specific interactions are given the interaction type as a second argument, so their keys can depend on it. Define this along with Molly.parameter_fields to make a custom interaction work with extract_parameters and inject_gradients.

place_atoms(n_atoms, boundary; min_dist=nothing, max_attempts=100,
            rng=Random.default_rng())

Generate random coordinates.

Obtain n_atoms coordinates in bounding box boundary where no two points are closer than min_dist, accounting for periodic boundary conditions. The keyword argument max_attempts determines the number of failed tries after which to stop placing atoms. Cannot be used if one or more dimensions has infinite boundaries.

place_diatomics(n_molecules, boundary, bond_length; min_dist=nothing,
                max_attempts=100, aligned=false, rng=Random.default_rng())

Generate random diatomic molecule coordinates.

Obtain coordinates for n_molecules diatomics in bounding box boundary where no two points are closer than min_dist and the bond length is bond_length, accounting for periodic boundary conditions. The keyword argument max_attempts determines the number of failed tries after which to stop placing atoms. The keyword argument aligned determines whether the bonds all point the same direction (true) or random directions (false). Cannot be used if one or more dimensions has infinite boundaries.