Molly API

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.

Not currently compatible with TriclinicBoundary. Not currently compatible with virial calculation.

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.

AndersenThermostat(temperature, coupling_const)

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.

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

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

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 can not 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.

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-λ_{ij}) &, r_{ij}\leq 2^{1/6}σ \\ λ_{ij}V_{\text{LJ}}(r_{ij}) &, 2^{1/6}σ \leq r_{ij} \end{cases}\]

and the force on each atom by

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

where

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

If $\lambda$ is one 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::Int=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.

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; 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)

The Berendsen barostat for controlling pressure.

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 are :isotropic, :semiisotropic and :anisotropic.

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

BerendsenThermostat(temperature, coupling_const)

The Berendsen thermostat for controlling temperature.

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.

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

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.

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

The stochastic cell rescale barostat.

See Bernetti and Bussi 2020 and Del Tatto et al. 2022. In brief, this is an extension of the Berendsen barostat 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 dt}{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 scaling factor $\mu$ is a matrix, allowing non-isotropic pressure control. Available options are :isotropic, :semiisotropic and :anisotropic.

CellListMapNeighborFinder(; eligible, dist_cutoff, special, n_steps, x0,
                          unit_cell, dims)

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

This is the recommended neighbor finder on CPU. x0 and unit_cell are optional initial coordinates and system unit cell that improve the first approximation of the cell list structure. The number of dimensions dims is inferred from unit_cell or x0, or assumed to be 3 otherwise.

Can not be used if one or more dimensions has infinite boundaries.

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

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

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. dist_cutoff and error_tol should match the general interaction.

dist_cutoff is the cutoff distance for short range interactions. approximate_erfc determines whether to use a fast approximation to the erfc function.

CoulombReactionField(; dist_cutoff, solvent_dielectric, use_neighbors, weight_special,
                        coulomb_const)

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

CoulombSoftCoreBeutler(; cutoff, α, λ, use_neighbors, σ_mixing, ϵ_mixing,
                       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. The potential energy is defined as

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

and the force on each atom by

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

where

\[r_{Q} = (\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}})+r_{ij}^6)\]

and

\[C^{(12)} = 4\epsilon\sigma^{12} C^{(6)} = 4\epsilon\sigma^{6}\]

If $\lambda$ is 1.0, this gives the standard Coulomb potential and means the atom is fully turned on. If $\lambda$ is zero the interaction is turned off. $\alpha$ determines the strength of softening the function.

CoulombSoftCoreGapsys(; cutoff, α, λ, σQ, use_neighbors, 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. The potential energy is defined as

\[V(r_{ij}) = \left\{ \begin{array}{cl} \frac{1}{4\pi\epsilon_0} \frac{q_iq_j}{r_{ij}}, & \text{if} & r \ge r_{LJ} \\ \frac{1}{4\pi\epsilon_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 \lt r_{LJ} \\ \end{array} \right.\]

and the force on each atom by

\[\vec{F}_i = \left\{ \begin{array}{cl} \frac{1}{4\pi\epsilon_0} \frac{q_iq_j}{r_{ij}^2}\frac{\vec{r_{ij}}}{r_{ij}}, & \text{if} & r \ge r_{LJ} \\ \frac{1}{4\pi\epsilon_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 \lt r_{LJ} \\ \end{array} \right.\]

where

\[r_{Q} = \alpha(1-\lambda)^{1/6}(1+σ_Q|qi*qj|)\]

If $\lambda$ is 1.0, this gives the standard Coulomb potential and means the atom is fully turned on. If $\lambda$ is zero the interaction is turned off. $\alpha$ determines the strength of softening the function.

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

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.

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.

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(; eligible, dist_cutoff, special, n_steps)

Find close atoms by distance.

This is the recommended neighbor finder on non-NVIDIA GPUs.

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, eligible=nothing, special=nothing)

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

Should be used alongside the CoulombEwald pairwise interaction, which provide the short range term. dist_cutoff and error_tol should match CoulombEwald.

dist_cutoff is the cutoff distance for short range interactions. eligible indicates pairs eligible for short range interaction, and can be a matrix like the neighbor list or nothing to indicate that all pairs are eligible. special should also be given where relevant, as these interactions are excluded from long range calculation.

This algorithm is O(N^2) and in general PME should be used instead. Only compatible with 3D systems and CubicBoundary. Runs on the CPU, even for GPU systems.

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

GPUNeighborFinder(; eligible, dist_cutoff, special, n_steps_reorder, initialized)

Use the non-bonded forces/potential energy algorithm from Eastman and Pande 2010 to avoid calculating a neighbor list.

This is the recommended neighbor finder on NVIDIA GPUs.

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; n_threads::Integer)::T.

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

HamiltonianREMD(; <keyword arguments>)

A simulator for a parallel Hamiltonian replica exchange MD (H-REMD) simulation on a ReplicaSystem.

The replicas are expected to have different Hamiltonians, i.e. different interactions. When calling simulate!, the assign_velocities keyword argument determines whether to assign random velocities at the appropriate temperature for each replica.

Arguments

• dt::DT: the time step of the simulation.
• temperature::T: the temperatures of the simulation.
• simulators::ST: individual simulators for simulating each replica.
• exchange_time::ET: the time interval between replica exchange attempts.

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

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

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

ImmediateThermostat(temperature)

The immediate velocity rescaling thermostat for controlling temperature.

Velocities are immediately rescaled to match a target temperature. 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)

GBn2 solvation model implemented as an AtomsCalculators.jl calculator.

Should be used along with a Coulomb interaction.

Not currently compatible with virial calculation.

ImplicitSolventOBC(atoms, atoms_data, bonds)

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

Should be used along with a Coulomb interaction. The keyword argument use_OBC2 determines whether to use parameter set I (false, the default) or II (true).

Not currently compatible with virial calculation.

InteractionList1Atoms(is, inters)
InteractionList1Atoms(is, inters, types)
InteractionList1Atoms(inter_type)

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

InteractionList2Atoms(is, js, inters)
InteractionList2Atoms(is, js, inters, types)
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)
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)
InteractionList4Atoms(inter_type)

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

Langevin(; <keyword arguments>)

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

See Zhang et al. 2019. This is a leapfrog 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.
• coupling::C=NoCoupling(): 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 currently 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.
• 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}^{6}}{r_{ij}^{6}}\right)^2 -\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{6}\right] \vec{r}_{ij} \end{aligned}\]

The potential energy does not include the long range dispersion correction present in some other implementations that approximately represents contributions from beyond the cutoff distance.

LennardJonesSoftCoreBeutler(; cutoff, α, λ, use_neighbors, shortcut, σ_mixing,
                            ϵ_mixing, 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. The potential energy is defined as

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

and the force on each atom by

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

where

\[r_{LJ} = (\frac{\alpha(1-\lambda)C^{(12)}}{C^{(6)}}+r^6)^{1/6}\]

and

\[C^{(12)} = 4\epsilon\sigma^{12} C^{(6)} = 4\epsilon\sigma^{6}\]

If $\lambda$ is 1.0, this gives the standard LennardJones potential and means the atom is fully turned on. If $\lambda$ is zero the interaction is turned off. $\alpha$ determines the strength of softening the function.

LennardJonesSoftCoreGapsys(; cutoff, α, λ, use_neighbors, shortcut, σ_mixing,
                           ϵ_mixing, 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. The potential energy is defined as

\[V(r_{ij}) = \left\{ \begin{array}{cl} \frac{C^{(12)}}{r_{ij}^{12}} - \frac{C^{(6)}}{r_{ij}^{6}}, & \text{if} & r \ge r_{LJ} \\ (\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}}, & \text{if} & r \lt r_{LJ} \\ \end{array} \right.\]

and the force on each atom by

\[\vec{F}_i = \left\{ \begin{array}{cl} (\frac{12C^{(12)}}{r_{ij}^{13}} - \frac{6C^{(6)}}{r_{ij}^{7}})\frac{\vec{r_{ij}}}{r_{ij}}, & \text{if} & r \ge r_{LJ} \\ ((\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}}))\frac{\vec{r_{ij}}}{r_{ij}}, & \text{if} & r \lt r_{LJ} \\ \end{array} \right.\]

where

\[r_{LJ} = \alpha(\frac{26C^{(12)}(1-\lambda)}{7C^{(6)}})^{\frac{1}{6}}\]

and

\[C^{(12)} = 4\epsilon\sigma^{12} C^{(6)} = 4\epsilon\sigma^{6}\]

If $\lambda$ is 1.0, this gives the standard LennardJones potential and means the atom is fully turned on. If $\lambda$ is zero the interaction is turned off. $\alpha$ determines the strength of softening the function.

MetropolisMonteCarlo(; <keyword arguments>)

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

Arguments

• temperature::T: the temperature of the system.
• trial_moves::M: a function that performs the trial moves.
• trial_args::Dict: a dictionary of 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}\]

MolecularForceField(ff_files...; units=true)
MolecularForceField(T, ff_files...; units=true)
MolecularForceField(atom_types, residue_types, bond_types, angle_types,
                    torsion_types, torsion_order, weight_14_coulomb,
                    weight_14_lj, attributes_from_residue)

A molecular force field.

Read one or more OpenMM force field XML files by passing them to the constructor.

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 currently compatible with virial calculation. Not currently 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)

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 currently work with shear stresses, the anisotropic variant only applies independent linear scaling of the box vectors. If shear deformation is required the BerendsenBarostat or, preferrably, the CRescaleBarostat should be used instead.

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

MullerBrown(; A, a, b, c, x0, y0, force_units, energy_units)

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

The potential energy is defined as

\[V(x,y) = \sum_{n=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 currently compatible with virial calculation.

NeighborList(n, list)
NeighborList()

Structure to contain neighbor lists.

NoCoupling()

Placeholder coupler that does nothing.

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.

Not currently 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=NoCoupling(): 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.

OverdampedLangevin(; <keyword arguments>)

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

Not currently 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 of 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.

PME(dist_cutoff, atoms, boundary; error_tol=0.0005, order=5,
    ϵr=1.0, fixed_charges=true, eligible=nothing, special=nothing,
    grad_safe=false, n_threads=Threads.nthreads())

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

Should be used alongside the CoulombEwald pairwise interaction, which provide the short range term. dist_cutoff and error_tol should match CoulombEwald.

dist_cutoff is the cutoff distance for short range interactions. eligible indicates pairs eligible for short range interaction, and can be a matrix like the neighbor list or nothing to indicate that all pairs are eligible. special should also be given where relevant, as these interactions are excluded from long range calculation. fixed_charges should be set to false if the partial charges can change, for example when using a polarizable force field. 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.

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.

Base type for pairwise interactions.

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

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

Only compatible with 3D systems.

RBTorsion(; f1, f2, f3, f4)

A Ryckaert-Bellemans torsion angle between four atoms.

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

ReplicaSystem(; <keyword arguments>)

A wrapper for replicas in a replica exchange simulation.

Each individual replica is a System. Properties unused in the simulation or in analysis can be left with their default values. The minimal required arguments are atoms, replica_coords, boundary and n_replicas. atoms and the elements in replica_coords should have the same length, along with atoms_data and the elements in replica_velocities if these are provided. The number of elements in replica_coords, replica_boundaries, replica_velocities, replica_loggers and the interaction arguments replica_pairwise_inters, replica_specific_inter_lists, replica_general_inters and replica_constraints should be equal to n_replicas if used. This is a sub-type of AbstractSystem from AtomsBase.jl and implements the interface described there.

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 nbatches=(min(n, 8), n) to CellListMapNeighborFinder during construction where n is the number of threads to be used per replica.

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.
• replica_coords: the coordinates of the atoms in each replica.
• n_replicas::Integer: the number of replicas of the system.
• boundary=nothing: the bounding box in which the simulation takes place. This is only used if no value is passed to the argument replica_pairwise_inters.
• replica_boundaries=nothing: the bounding box for each replica.
• replica_velocities=[zero(replica_coords[1]) * u"ps^-1" for _ in 1:n_replicas]: the velocities of the atoms in each replica.
• 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 (to be used if the same for all replicas). This is only used if no value is passed to the argument replica_topology.
• replica_topology=[nothing for _ in 1:n_replicas]: the topological information for each replica.
• pairwise_inters=(): the pairwise interactions in the system, i.e. interactions between all or most atom pairs such as electrostatics (to be used if the same for all replicas). Should be a Tuple or NamedTuple of PairwiseInteractions. This is only used if no value is passed to the argument replica_pairwise_inters.
• replica_pairwise_inters=[() for _ in 1:n_replicas]: the pairwise interactions for each replica.
• specific_inter_lists=(): the specific interactions in the system, i.e. interactions between specific atoms such as bonds or angles (to be used if the same for all replicas). Should be a Tuple or NamedTuple. This is only used if no value is passed to the argument replica_specific_inter_lists.
• replica_specific_inter_lists=[() for _ in 1:n_replicas]: the specific interactions in each replica.
• general_inters=(): the general interactions in the system, i.e. interactions involving all atoms such as implicit solvent (to be used if the same for all replicas). Each should implement the AtomsCalculators.jl interface. Should be a Tuple or NamedTuple. This is only used if no value is passed to the argument replica_general_inters.
• replica_general_inters=[() for _ in 1:n_replicas]: the general interactions for each replica.
• constraints::CN=(): the constraints for bonds and angles in the system (to be used if the same for all replicas). Should be a Tuple or NamedTuple. This is only used if no value is passed to the argument replica_constraints.
• replica_constraints=[() for _ in 1:n_replicas]: the constraints for bonds and angles in each replica.
• neighbor_finder::NF=NoNeighborFinder(): the neighbor finder used to find close atoms and save on computation. It is duplicated for each replica.
• replica_loggers=[() for _ in 1:n_replicas]: the loggers for each replica that record properties of interest during a simulation.
• exchange_logger::EL=ReplicaExchangeLogger(n_replicas): the logger used to record the exchange of replicas.
• 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.
• data::DA=nothing: arbitrary data associated with the replica system.

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)

Constrain distances during a simulation using the SHAKE and RATTLE algorithms. Either or both of dist_constraints and angle_constraints must be given.

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 Krautler et al. 2000 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 if iterations is reached, some constraints may not be satisfied.

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

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.

SteepestDescentMinimizer(; <keyword arguments>)

Steepest descent energy minimization.

Arguments

• step_size::D=0.01u"nm": the initial maximum displacement.
• max_steps::Int=1000: the maximum number of steps.
• tol::F=1000.0u"kJ * mol^-1 * nm^-1": the maximum force below which to finish minimization.
• log_stream::L=devnull: stream to print minimization progress to.

StormerVerlet(; <keyword arguments>)

The Störmer-Verlet integrator.

The velocity calculation is accurate to O(dt).

Does not currently work with coupling methods that alter the velocity. Does not currently remove the center of mass motion.

Arguments

• dt::T: the time step of the simulation.
• coupling::C=NoCoupling(): the coupling which applies during 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.

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.
• 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.
• data::DA=nothing: arbitrary data associated with the system.

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.

System(coordinate_file, force_field; <keyword arguments>)

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

Atom names should exactly match residue templates - no searching of residue templates is carried out.

System(coordinate_file, topology_file; <keyword arguments>)
System(T, coordinate_file, topology_file; <keyword arguments>)

Read a Gromacs coordinate file and a Gromacs topology file with all includes collapsed into one file.

Gromacs file reading should be considered experimental. The rename_terminal_res keyword argument is not available when reading Gromacs files.

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.
• loggers=(): the loggers that record properties of interest during a simulation.
• units::Bool=true: whether to use Unitful quantities.
• array_type=Array: the array type for the simulation, for example use CuArray or ROCArray for GPU support.
• dist_cutoff=1.0u"nm": cutoff distance for long-range interactions.
• dist_buffer=0.2u"nm": distance added to dist_cutoff when calculating neighbors every few steps. Not relevant if GPUNeighborFinder is used since the neighbors are calculated each step.
• constraints=:none: which constraints to apply during the simulation, options are :none, :hbonds (bonds involving hydrogen), :allbonds 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.
• 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.
• nonbonded_method=:none: method for long range interaction summation, options are :none (short range only), :cutoff (reaction field method), :pme (particle mesh Ewald summation) and :ewald (Ewald summation, slow).
• ewald_error_tol=0.0005: the error tolerance for Ewald summation, used when nonbonded_method is :pme or :ewald.
• approximate_pme=true: whether to use a fast approximation to the erfc function, used when nonbonded_method is :pme.
• 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 DistanceNeighborFinder on non-CUDA compatible GPUs.
• data=nothing: arbitrary data associated with the system.
• implicit_solvent=:none: the implicit solvent model to use, options are :none, :obc1, :obc2 and :gbn2.
• kappa=0.0u"nm^-1": the kappa value for the implicit solvent model if one is used.
• rename_terminal_res=true: whether to rename the first and last residues to match the appropriate atom templates, for example the first (N-terminal) residue could be changed from "MET" to "NMET".
• grad_safe=false: should be set to true if the system is going to be used with Enzyme.jl and nonbonded_method is :pme.

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 paramaters like σ have to be added manually after construction using the convenience constructor System(sys; <keyword arguments>).

TemperatureREMD(; <keyword arguments>)

A simulator for a parallel temperature replica exchange MD (T-REMD) simulation on a ReplicaSystem.

See Sugita and Okamoto 1999. The corresponding ReplicaSystem should have the same number of replicas as the number of temperatures in the simulator. When calling simulate!, the assign_velocities keyword argument determines whether to assign random velocities at the appropriate temperature for each replica.

Arguments

• dt::DT: the time step of the simulation.
• temperatures::TP: the temperatures corresponding to the replicas.
• simulators::ST: individual simulators for simulating each replica.
• exchange_time::ET: the time interval between replica exchange attempts.

ThermoState(name::AbstractString, β, p, system)
ThermoState(system::System, β, p; name::Union{Nothing, AbstractString}=nothing)

Thermodynamic state wrapper carrying inverse temperature β = 1/kBT, pressure p, and the System used to evaluate energies.

Fields:

• name::String - label for the state.
• β - inverse temperature with units compatible with 1/system.energy_units.
• p - pressure Quantity or nothing.
• system::System - simulation system used to compute potential energy.

The second constructor checks unit consistency for β and p and sets a default name when not provided.

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; 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)

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.

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(; eligible, dist_cutoff, special, n_steps)

Find close atoms by distance using a tree search.

Can not be used if one or more dimensions has infinite boundaries. Can not 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 currently 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\]

VelocityRescaleThermostat(temperature, coupling_const; n_steps=1)

The stochastic velocity rescaling thermostat.

See Bussi et al. 2007. In brief, acts like the BerendsenThermostat but adds an stochastic term, allowing correct sampling of isothermal ensembles.

Let Δt be the simulation timestep, Nf the kinetic DOFs used to calculate the instantaneous temperature of the system. Then, $K = \frac{1}{2} \cdot \sum m \cdot v^2$ is the current kinetic energy, and $K̄ = \frac{1}{2} Nf k_B T_0$ is the target kinetic energy for a reference temperature T_0.

Define $c = e^{-Δt/τ}$. Draw $R \sim 𝒩(0,1)$ and $S \sim \chi^{2}_{Nf-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=NoCoupling(): 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=NoCoupling(): 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.

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

\[F(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 r_{ij}^2} \exp(-\kappa r_{ij})\left(\kappa r_{ij} + 1\right) \vec{r}_{ij}\]

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

Access the stored observations in a logger.

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

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

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.

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.

This should only be used on 3-dimensional systems.

ScalarPressureLogger(n_steps)
ScalarPressureLogger(T, n_steps)

Log the scalar_pressure of a system throughout a simulation.

This should only be used on 3-dimensional systems.

ScalarVirialLogger(n_steps)
ScalarVirialLogger(T, n_steps)

Log the scalar_virial tensor of a system throughout a simulation.

TemperatureLogger(n_steps)
TemperatureLogger(T, n_steps)

Log the temperature throughout a simulation.

TotalEnergyLogger(n_steps)
TotalEnergyLogger(T, n_steps)

Log the total_energy of a system throughout a simulation.

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.

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

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=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, is_heavy_atom means non-hydrogen atoms are restrained.

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

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, buffers, neighbors=nothing, step_n=0, run_loggers=true;
               n_threads=Threads.nthreads(), kwargs...)

Run the loggers associated with a system.

run_loggers can be true, false or :skipzero, in which case the loggers are not run before the first step. Additional keyword arguments can be passed to the loggers if required. Ignored for gradient calculation during automatic differentiation.

apply_position_constraints!(sys, coord_storage)
apply_position_constraints!(sys, coord_storage, vel_storage, dt)

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)

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.

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::Vector{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

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

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.

charge(atom)

The partial charge of an Atom.

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.

density(sys)

The density of a System.

Returns zero density for infinite boundaries.

dipole_moment(sys)

The dipole moment μ of a system.

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.

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.

float_type(sys)
float_type(boundary)

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

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)
force(inter, coord_i, coord_j, boundary, atom_i, atom_j, force_units, velocity_i,
      velocity_j, step_n)
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)
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)

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.

forces(system, neighbors=find_neighbors(system), step_n=0;
       n_threads=Threads.nthreads())

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

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

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

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

is_any_atom(at, at_data)

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

is_heavy_atom(at, at_data)

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

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

log_property!(logger, system, buffers, neighbors=nothing, step_n=0;
              n_threads=Threads.nthreads(), 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.

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.

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

pairwise_force(inter, r, params)

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

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.

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. Can not 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). Can not be used if one or more dimensions has infinite boundaries.

pmf_with_uncertainty(u, u_target, f, N_counts, logN, R_k;
                     shifts=nothing, nbins=nothing, edges=nothing, kBT=nothing,
                     zero=:min, rmin=nothing, rmax=nothing)

Estimate a 1D PMF along a scalar CV and its large-sample uncertainty using MBAR, see Shirts and Chodera 2008 Eq D8.

Arguments

• u::AbstractMatrix - N×K reduced potentials.
• u_target::AbstractVector - length-N reduced potentials at the target state.
• f::AbstractVector - length-K relative free energies.
• N_counts::AbstractVector - length-K sample counts.
• logN::AbstractVector - log.(N_counts).
• R_k::Vector{<:AbstractVector} - CV values per window, concatenating to length N.

Keyword arguments

• shifts=nothing - per-frame shifts used on u, if any.
• nbins=nothing, edges=nothing - bin count or explicit bin edges.
• rmin=nothing, rmax=nothing - bounds when edges is omitted.
• zero=:min - PMF gauge: :min or :last.
• kBT=nothing - if set, also return dimensional F_energy = F*kBT.

Returns

PMF struct with:

• centers - the center of the histogram bins used to sample the CV.
• widths - the width of the histogram bins used to sample the CV.
• edges - the edges of the histogram bins used to sample the CV.
• F - PMF in kBT units.
• F_energy - PMF in energy units, if provided.
• sigma_F - uncertainty of the PMF in kBT units.
• sigma_F_energy - uncertainty of the PMF in energy units, if provided.
• p - probability density along the CV.
• var_p - variance of the probability density along the CV.

pmf_with_uncertainty(coords_k, boundaries_k, states, target_state, CV;
                     shift=false)

High-level PMF wrapper that builds MBAR inputs from trajectories, solves MBAR, and computes the PMF along CV.

Arguments

• coords_k::AbstractVector - coordinates per window.
• boundaries_k::AbstractVector - boundaries per window.
• states::Vector{ThermoState} - thermodynamic states per window.
• target_state::ThermoState - state to reweight to.
• CV::AbstractVector - CV values per window.

Keyword arguments

• shift::Bool=false - subtract per-frame minima from u for stability.

Returns

Same PMF struct as the lower-level pmf_with_uncertainty.

potential_energy(system, neighbors=find_neighbors(system), step_n=0;
                 n_threads=Threads.nthreads())

Calculate the potential energy of a system using the pairwise, specific and general interactions.

potential_energy(inter, vec_ij, atom_i, atom_j, energy_units, special, coord_i, coord_j,
                 boundary, velocity_i, velocity_j, step_n)
potential_energy(inter, coord_i, boundary, atom_i, energy_units, velocity_i, step_n)
potential_energy(inter, coord_i, coord_j, boundary, atom_i, atom_j, energy_units,
                 velocity_i, velocity_j, step_n)
potential_energy(inter, coord_i, coord_j, coord_k, boundary, atom_i, atom_j, atom_k,
                 energy_units, velocity_i, velocity_j, velocity_k, step_n)
potential_energy(inter, coord_i, coord_j, coord_k, coord_l, boundary, atom_i, atom_j,
                 atom_k, atom_l, energy_units, velocity_i, velocity_j, velocity_k,
                 velocity_l, step_n)

Calculate the potential energy due to a given interaction type.

Custom interaction types should implement this function.

pressure(system, neighbors=find_neighbors(system), step_n=0, buffers=nothing;
         recompute=true, n_threads=Threads.nthreads())

Calculate the pressure tensor of the system.

The pressure is defined as

\[\bf{P} = \frac{ 2 \cdot \bf{K} + \bf{W} }{V}\]

where $V$ is the system volume, $\bf{K}$ is the kinetic energy tensor and $\bf{W}$ is the virial tensor.

To calculate the scalar pressure, see scalar_pressure.

Not compatible with infinite boundaries.

radius_gyration(coords, atoms)

Calculate the radius of gyration of a set of coordinates.

Assumes the coordinates do not cross the bounding box, i.e. all coordinates correspond to the same periodic image.

random_coord(boundary; rng=Random.default_rng())

Generate a random coordinate uniformly distributed within a bounding box.

random_normal_translation!(sys::System; shift_size=oneunit(eltype(eltype(sys.coords))),
                           rng=Random.default_rng())

Performs a random translation of the coordinates of a randomly selected atom in a System.

The translation is generated using a uniformly chosen direction and length selected from the standard normal distribution i.e. with mean 0 and standard deviation 1, scaled by shift_size which should have appropriate length units.

random_uniform_translation!(sys::System; shift_size=oneunit(eltype(eltype(sys.coords))),
                            rng=Random.default_rng())

Performs a random translation of the coordinates of a randomly selected atom in a System.

The translation is generated using a uniformly selected direction and uniformly selected length in range [0, 1) scaled by shift_size which should have appropriate length units.

random_velocities!(sys, temp; rng=Random.default_rng())
random_velocities!(vels, sys, temp; rng=Random.default_rng())

Set the velocities of a System, or a vector, to random velocities generated from the Maxwell-Boltzmann distribution.

random_velocities(sys, temp; rng=Random.default_rng())

Generate random velocities from the Maxwell-Boltzmann distribution for a System.

random_velocity(atom_mass::Union{Unitful.Mass, MolarMass}, temp::Unitful.Temperature;
                dims=3, rng=Random.default_rng())
random_velocity(atom_mass::Union{Unitful.Mass, MolarMass}, temp::Unitful.Temperature,
                k::Union{BoltzmannConstUnits, MolarBoltzmannConstUnits};
                dims=3, rng=Random.default_rng())
random_velocity(atom_mass::Real, temp::Real, k::Real=ustrip(u"u * nm^2 * ps^-2 * K^-1", Unitful.k);
                dims=3, rng=Random.default_rng())

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

rdf(coords, boundary; npoints=200)

Calculate the radial distribution function of a set of coordinates.

This function is only available when KernelDensity is imported. This describes how density varies as a function of distance from each atom. Returns a list of distance bin centers and a list of the corresponding densities.

read_frame!(ens_sys::EnsembleSystem, frame_idx::Integer)

Read a frame from an EnsembleSystem and return a System representing the frame.

remd_exchange!(sys, sim, n, m; n_threads=Threads.nthreads(), rng=Random.default_rng())

Attempt an exchange of replicas n and m in a ReplicaSystem during a REMD simulation.

Successful exchanges should exchange coordinates and velocities as appropriate. Returns acceptance quantity Δ and a Bool indicating whether the exchange was successful.

remove_CM_motion!(system)

Remove the center of mass motion from a System.

rmsd(coords_1, coords_2)

Calculate the root-mean-square deviation (RMSD) of two sets of 3D coordinates after superimposition by the Kabsch algorithm.

Assumes the coordinates do not cross the bounding box, i.e. all coordinates in each set correspond to the same periodic image.

scalar_pressure(system, neighbors=find_neighbors(system), step_n=0, buffers=nothing;
                recompute=true, n_threads=Threads.nthreads())

Calculate the pressure of the system as a scalar.

This is the trace of the pressure tensor.

scalar_virial(system, neighbors=find_neighbors(system), step_n=0;
              n_threads=Threads.nthreads())

Calculate the virial of the system as a scalar.

This is the trace of the virial tensor.

scale_boundary(boundary, scale_factor)

Scale the sides of a bounding box by a scaling factor.

The scaling factor can be a single number or a SVector of the appropriate number of dimensions corresponding to the scaling factor for each axis. For a 3D bounding box the volume scales as the cube of the scaling factor.

scale_coords!(sys::System{<:Any, AT}, μ::SMatrix{D,D};
              rotate::Bool=true,
              ignore_molecules::Bool=false,
              scale_velocities::Bool=false)

Rigid-molecular barostat update with optional rotation.

• Box: B′ = μ * B
• Positions: r′ = μ * r (implemented via COM affine + optional rotation of internal offsets)
• Velocities: v′ = μ⁻¹ * v (applied when scale_velocities=true)

simulate!(system, simulator, n_steps; <keyword arguments>)
simulate!(system, simulator; <keyword arguments>)

Run a simulation on a system according to the rules of the given simulator.

Custom simulators should implement this function. Constraints are applied during minimization, which can lead to issues.

Arguments

• n_threads=Threads.nthreads(): the number of threads to run the simulation on, only relevant when running on CPU.
• run_loggers: whether to run the loggers during the simulation. Can be true, false or :skipzero, in which case the loggers are not run before the first step. run_loggers is true by default except for SteepestDescentMinimizer, where it is false.
• rng=Random.default_rng(): the random number generator used for the simulation. Setting this allows reproducible stochastic simulations.

simulate_remd!(sys, remd_sim, n_steps; n_threads=Threads.nthreads(),
               run_loggers=true, rng=Random.default_rng())

Run a REMD simulation on a ReplicaSystem using a REMD simulator.

Not currently compatible with interactions that depend on step number.

statistical_inefficiency(series::AbstractVector; maxlag::Union{Nothing, Int}=nothing)

Integrated autocorrelation time estimator with IPS truncation and finite-sample taper.

Returns a StatisticalInefficiency struct with:

• inefficiency: statistical inefficiency
• stride: ceil(Int, g)
• input_length: input length
• effective_size: floor(N / stride)
• lag: truncation lag used

Notes:

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

temperature(system; kin_tensor=nothing, recompute=true)

Calculate the temperature of a system from the kinetic energy of the atoms.

torsion_angle(coord_i, coord_j, coord_k, coord_l, boundary)
torsion_angle(vec_ij, vec_jk, vec_kl)

Calculate the torsion angle in radians defined by four coordinates or three vectors.

The angle between the planes defined by atoms (i, j, k) and (j, k, l) is returned in the range -π to π.

total_energy(system, neighbors=find_neighbors(sys); n_threads=Threads.nthreads())

Calculate the total energy of a system as the sum of the kinetic_energy and the potential_energy.

use_neighbors(inter)

Whether a pairwise interaction uses the neighbor list, default false.

Custom pairwise interactions can define a method for this function. For built-in interactions such as LennardJones this function accesses the use_neighbors field of the struct.

ustrip_vec(x)
ustrip_vec(u, x)

Broadcasted form of ustrip from Unitful.jl, allowing e.g. ustrip_vec.(coords).

vector(c1, c2, boundary)

Displacement between two coordinate values from c1 to c2, accounting for periodic boundary conditions.

The minimum image convention is used, so the displacement is to the closest version of the coordinates accounting for the periodic boundaries. For the TriclinicBoundary an approximation is used to find the closest version by default.

vector_1D(c1, c2, side_length)

Displacement between two 1D coordinate values from c1 to c2, accounting for periodic boundary conditions in a CubicBoundary or RectangularBoundary.

The minimum image convention is used, so the displacement is to the closest version of the coordinate accounting for the periodic boundaries.

virial(system, neighbors=find_neighbors(system), step_n=0;
       n_threads=Threads.nthreads())

Calculate the virial tensor of the system.

The virial, in its most general form, is defined as:

\[\bf{W} = \sum_i \bf{r_i} \otimes \bf{f_i}\]

where $\bf{r_i}$ and $\bf{f_i}$ are the position and force vectors, respectively, acting on the i$^{th}$ atom. In Molly.jl, we implement the virial definition used in LAMMPS, and take into account pairwise and specific interactions, as well as the K-space contribution of the Ewald and PME methods, computed as indicated in the Essmann et al. 1995. Contributions from constraints are ignored.

To calculate the scalar virial, see scalar_virial.

visualize(coord_logger, boundary, out_filepath; <keyword arguments>)

Visualize a simulation as an animation.

This function is only available when GLMakie is imported. It can take a while to run, depending on the length of the simulation and the number of atoms.

Arguments

• connections=Tuple{Int, Int}[]: pairs of atoms indices to link with bonds.
• connection_frames: the frames in which bonds are shown. Should be a list of the same length as the number of frames, where each item is a list of Bools of the same length as connections. Defaults to always true.
• trails::Integer=0: the number of preceding frames to show as transparent trails.
• framerate::Integer=30: the frame rate of the animation.
• color=:purple: the color of the atoms. Can be a single color or a list of colors of the same length as the number of atoms.
• connection_color=:orange: the color of the bonds. Can be a single color or a list of colors of the same length as connections.
• markersize=0.05: the size of the atom markers, in the units of the data.
• linewidth=2.0: the width of the bond lines.
• transparency=true: whether transparency is active on the plot.
• show_boundary::Bool=true: whether to show the bounding box as lines.
• boundary_linewidth=2.0: the width of the boundary lines.
• boundary_color=:black: the color of the boundary lines.
• kwargs...: other keyword arguments are passed to the point plotting function.

volume(sys)
volume(boundary)

Calculate the volume (3D) or area (2D) of a System or bounding box.

Returns infinite volume for infinite boundaries.

wrap_coord_1D(c, side_length)

Ensure a 1D coordinate is within the bounding box and return the coordinate.

wrap_coords(c, boundary)

Ensure a coordinate is within the bounding box and return the coordinate.

write_structure(filepath, sys; format="", correction=:pbc, atom_inds=[],
                excluded_res=String[], write_velocities=false, write_boundary=true)

Write the 3D structure of a system to a file.

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 overwritten if it already exists.

Not compatible with 2D systems. The CRYST1 record is not written for infinite boundaries.