Molly API

Generalized Born (GB) implicit solvent models augmented with the hydrophobic solvent accessible surface area (SA) term. Custom GBSA methods should sub-type this type.

AndersenThermostat(temperature, coupling_const)

Rescale random velocities according to the Andersen thermostat.

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: the index of the atom in the system.
• charge::C=0.0: the charge of the atom, used for electrostatic interactions.
• mass::M=1.0u"u": the mass of the atom.
• σ::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.
• solute::Bool=false: whether the atom is part of the solute.

AtomData(atom_type, atom_name, res_number, res_name)

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.

BerendsenThermostat(temperature, coupling_const)

The Berendsen thermostat. This thermostat should be used with caution as it can lead to simulation artifacts. The scaling factor for the velocities each step is

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

Buckingham(; cutoff, nl_only, weight_14, force_units, energy_units)

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.

CellListMapNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff, x0, unit_cell)

Find close atoms by distance and store auxiliary arrays for in-place threading. x0 and unit_cell are optional initial coordinates and system unit cell that improve the first approximation of the cell list structure. Can not be used if one or more dimensions has infinite boundaries.

Example

julia> coords
15954-element Vector{SVector{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}}:
 [2.5193063341012127 nm, 3.907448346081021 nm, 4.694954671434135 nm]
 [2.4173958848835233 nm, 3.916034913604175 nm, 4.699661024574953 nm]
 ⋮
 [1.818842280373283 nm, 5.592152965227421 nm, 4.992100424805031 nm]
 [1.7261366568663976 nm, 5.610326185704369 nm, 5.084523386833478 nm]

julia> boundary
CubicBoundary{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}(Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}[5.676 nm, 5.6627 nm, 6.2963 nm])

julia> neighbor_finder = CellListMapNeighborFinder(
           nb_matrix=s.neighbor_finder.nb_matrix, matrix_14=s.neighbor_finder.matrix_14, 
           n_steps=10, dist_cutoff=1.2u"nm",
           x0=coords, unit_cell=boundary,
       )
CellListMapNeighborFinder{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, 3, Float64}
  Size of nb_matrix = (15954, 15954)
  n_steps = 10
  dist_cutoff = 1.2 nm

CosineAngle(; k, θ0)

A cosine bond angle between three atoms. The potential energy is defined as

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

Coulomb(; cutoff, nl_only, weight_14, coulomb_const, force_units, energy_units)

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

CoulombReactionField(; dist_cutoff, solvent_dielectric, nl_only, weight_14,
                        coulomb_const, force_units, energy_units)

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

CoulombSoftCore(; cutoff, α, λ, p, nl_only, lorentz_mixing, weight_14,
                coulomb_const, force_units, energy_units)

The Coulomb electrostatic interaction between two atoms with a soft core. The potential energy is defined as

\[V(r_{ij}) = \frac{q_i q_j}{4 \pi \varepsilon_0 (r_{ij}^6 + \alpha \sigma_{ij}^6 \lambda^p)^{\frac{1}{6}}}\]

If $\alpha$ or $\lambda$ are zero this gives the standard Coulomb potential.

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

Cubic 3D bounding box defined by 3 side lengths. Setting one or more values to Inf gives no boundary in that dimension.

CubicSplineCutoff(dist_activation, dist_cutoff)

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

DistanceCutoff(dist_cutoff)

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

DistanceNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff)

Find close atoms by distance.

DistanceVecNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff)

Find close atoms by distance in a GPU and Zygote compatible manner.

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

GeneralObservableLogger(observable::Function, T, n_steps)

A logger which holds a record of regularly sampled observations of the system. observable should return an object of type T and support the method observable(s::System, neighbors; n_threads::Integer)::T.

Gravity(; G, nl_only)

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. Not currently compatible with automatic differentiation using Zygote.

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

ImplicitSolventGBN2(atoms, atoms_data, bonds)

GBn2 solvation model. Should be used along with a Coulomb or CoulombReactionField interaction.

ImplicitSolventOBC(atoms, atoms_data, bonds)

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

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

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

InteractionList2Atoms(is, js, types, inters)
InteractionList2Atoms(inter_type)

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

InteractionList3Atoms(is, js, ks, types, inters)
InteractionList3Atoms(inter_type)

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

InteractionList4Atoms(is, js, ks, ls, types, inters)
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. 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.
• remove_CM_motion::Bool=true: whether to remove the center of mass motion every time step.

LangevinSplitting(; <keyword arguments>)

The Langevin simulator using a general splitting scheme, consisting 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 automatic differentiation using Zygote.

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::Bool=true: whether to remove the center of mass motion every time step.

LennardJones(; cutoff, nl_only, lorentz_mixing, weight_14, weight_solute_solvent,
             force_units, energy_units, skip_shortcut)

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

LennardJonesSoftCore(; cutoff, α, λ, p, nl_only, lorentz_mixing, weight_14,
                     weight_solute_solvent, force_units, energy_units, skip_shortcut)

The Lennard-Jones 6-12 interaction between two atoms with a soft core. The potential energy is defined as

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

and the force on each atom by

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

where

\[r_{ij}^{\text{sc}} = \left(r_{ij}^6 + \alpha \sigma_{ij}^6 \lambda^p \right)^{1/6}\]

If $\alpha$ or $\lambda$ are zero this gives the standard LennardJones potential.

MetropolisMonteCarlo(; <keyword arguments>)

A Monte Carlo simulator that uses the Metropolis algorithm to sample the configuration space. simulate! for this simulator accepts an optional keyword argument log_states::Bool=true which determines whether to run the loggers or not (for example, during equilibration).

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, nl_only, lorentz_mixing, force_units, energy_units, skip_shortcut)

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

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

NeighborList()
NeighborList(n, list)

Structure to contain pre-allocated neighbor lists.

NeighborListVec(close, all)

Structure to contain neighbor lists for broadcasting. Each component may be present or absent depending on the interactions in the system.

NoCoupling()

Placeholder coupler that does nothing.

NoCutoff()

Placeholder cutoff that does not alter forces or potentials.

NoNeighborFinder()

Placeholder neighbor finder that returns no neighbors. When using this neighbor finder, ensure that nl_only for the interactions is set to false.

OpenMMAtomType(type, element, mass, σ, ϵ)

An OpenMM atom type.

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

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

OpenMMResidueType(name, types, charges, indices)

An OpenMM residue type.

A pairwise interaction that will apply to all or most atom pairs. Custom pairwise interactions should sub-type this type.

PeriodicTorsion(; periodicities, phases, ks, proper)

A periodic torsion angle between four atoms. The potential energy is defined as

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

PeriodicTorsionType(periodicities, phases, ks, proper)

A periodic torsion type.

RBTorsion(; f1, f2, f3, f4)

A Ryckaert-Bellemans torsion angle between four atoms.

RectangularBoundary(x, y)
RectangularBoundary(arr)

Rectangular 2D bounding box defined by 2 side lengths. 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. atoms, atoms_data and the elements in replica_coords and replica_velocities should have the same length. The number of elements in replica_coords, 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. 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.
• atoms_data::AD: other data associated with the atoms, allowing the atoms to be bits types and hence work on the GPU.
• pairwise_inters::PI=(): the pairwise interactions in the system, i.e. interactions between all or most atom pairs such as electrostatics (to be used if same for all replicas). Typically a Tuple. 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::SI=(): the specific interactions in the system, i.e. interactions between specific atoms such as bonds or angles (to be used if same for all replicas). Typically a Tuple. 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::GI=(): the general interactions in the system, i.e. interactions involving all atoms such as implicit solvent (to be used if same for all replicas). Typically a Tuple. 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 same for all replicas). Typically a Tuple.
• replica_constraints=[() for _ in 1:n_replicas]: the constraints for bonds and angles in each replica. This is only used if no value is passed to the argument replica_constraints.
• n_replicas::Integer: the number of replicas of the system.
• replica_coords: the coordinates of the atoms in each replica.
• replica_velocities=[zero(replica_coords[1]) * u"ps^-1" for _ in 1:n_replicas]: the velocities of the atoms in each replica.
• boundary::B: the bounding box in which the simulation takes place.
• neighbor_finder::NF=NoNeighborFinder(): the neighbor finder used to find close atoms and save on computation.
• exchange_logger::EL=ReplicaExchangeLogger(n_replicas): the logger used to record the exchange of replicas.
• replica_loggers=[() for _ in 1:n_replicas]: the loggers for each replica 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: the Boltzmann constant, which may be modified in some simulations.
• gpu_diff_safe::Bool: whether to use the code path suitable for the GPU and taking gradients. Defaults to isa(replica_coords[1], CuArray).

RescaleThermostat(temperature)

The velocity rescaling thermostat that immediately rescales the velocities to match a target temperature. This thermostat should be used with caution as it can lead to simulation artifacts. The scaling factor for the velocities each step is

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

SHAKE(dists, is, js)

Constrains a set of bonds to defined distances.

ShiftedForceCutoff(dist_cutoff)

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

ShiftedPotentialCutoff(dist_cutoff)

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

SoftSphere(; cutoff, nl_only, lorentz_mixing, force_units, energy_units, skip_shortcut)

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.

A specific interaction between sets of specific atoms, e.g. a bond angle. Custom specific interactions should sub-type this type.

SteepestDescentMinimizer(; <keyword arguments>)

Steepest descent energy minimization. Not currently compatible with automatic differentiation using Zygote.

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.
• run_loggers::Bool=false: whether to run the loggers during 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 every time step.

Arguments

• dt::T: the time step of the simulation.
• coupling::C=NoCoupling(): the coupling which applies during the simulation.

StructureWriter(n_steps, filepath, excluded_res=String[])

Write 3D output structures to the PDB file format throughout a 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. atoms, atoms_data, coords and velocities should have the same length. 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.
• atoms_data::AD: other data associated with the atoms, allowing the atoms to be bits types and hence work on the GPU.
• pairwise_inters::PI=(): the pairwise interactions in the system, i.e. interactions between all or most atom pairs such as electrostatics. Typically a Tuple.
• specific_inter_lists::SI=(): the specific interactions in the system, i.e. interactions between specific atoms such as bonds or angles. Typically a Tuple.
• general_inters::GI=(): the general interactions in the system, i.e. interactions involving all atoms such as implicit solvent. Typically a Tuple.
• constraints::CN=(): the constraints for bonds and angles in the system. Typically a Tuple.
• coords::C: the coordinates of the atoms in the system. Typically a vector of SVectors of 2 or 3 dimensions.
• velocities::V=zero(coords) * u"ps^-1": the velocities of the atoms in the system.
• boundary::B: the bounding box in which the simulation takes place.
• 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: the Boltzmann constant, which may be modified in some simulations.
• gpu_diff_safe::Bool: whether to use the code path suitable for the GPU and taking gradients. Defaults to isa(coords, CuArray).

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.

Arguments

• velocities=nothing: the velocities of the atoms in the system, set to zero by default.
• boundary=nothing: the bounding box used for simulation, read from the file by default.
• loggers=(): the loggers that record properties of interest during a simulation.
• units::Bool=true: whether to use Unitful quantities.
• gpu::Bool=false: whether to move the relevant parts of the system onto the GPU.
• gpu_diff_safe::Bool: whether to use the code path suitable for the GPU and taking gradients, defaults to the value of gpu.
• dist_cutoff=1.0u"nm": cutoff distance for long-range interactions.
• dist_neighbors=1.2u"nm": cutoff distance for the neighbor list, should be greater than dist_cutoff.
• implicit_solvent=nothing: specify a string to add an implicit solvent model, options are "obc1", "obc2" and "gbn2".
• center_coords::Bool=true: whether to center the coordinates in the simulation box.

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. Not currently compatible with automatic differentiation using Zygote.

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.

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

A time correlation logger, which 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.

TreeNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff)

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(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. Not currently compatible with automatic differentiation using Zygote.

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::Bool=true: whether to remove the center of mass motion every time step.

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::Bool=true: whether to remove the center of mass motion every time step.

n_dimensions(boundary)

Number of dimensions of a System, ReplicaSystem or bounding box.

velocity(mass, temperature; dims=3)
velocity(mass, temperature, k; dims=3)

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

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

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

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

Log the coordinates throughout a simulation.

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

Log the forces throughout a simulation.

KineticEnergyLogger(n_steps)
KineticEnergyLogger(T, n_steps)

Log the kinetic_energy of the system throughout a simulation.

PotentialEnergyLogger(n_steps)
PotentialEnergyLogger(T, n_steps)

Log the potential_energy of the 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 the system throughout a simulation.

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

Log the velocities throughout a simulation.

accelerations(system, neighbors=nothing; n_threads=Threads.nthreads())

Calculate the accelerations of all atoms using the pairwise, specific and general interactions and Newton's second law of motion. If the interactions use neighbor lists, the neighbors should be computed first and passed to the function.

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_constraints!(system, old_coords, dt)

Applies all the bond and angle constraints associated with the System.

apply_constraints!(sys, constraint, old_coords, dt)

Updates the coordinates and/or velocities of a System based on the constraints.

apply_coupling!(system, coupling, simulator)

Apply a coupler to modify a simulation. Custom couplers should implement this function.

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

born_radii_and_grad(inter, coords, boundary)

Calculate Born radii and gradients of Born radii and surface area overlap with respect to atomic distance. Custom GBSA methods should implement this function.

box_center(boundary)

Calculate the center of a bounding box. Dimensions with infinite length return zero.

box_volume(boundary)

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

charge(atom)

The partial charge of an Atom.

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(system, force_field)

Form a Dict of all parameters in a System, allowing gradients to be tracked.

find_neighbors(system; n_threads=Threads.nthreads())
find_neighbors(system, neighbor_finder, current_neighbors=nothing,
                step_n=0; 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::PairwiseInteraction, vec_ij, coord_i, coord_j,
      atom_i, atom_j, boundary)
force(inter::PairwiseInteraction, vec_ij, coord_i, coord_j,
      atom_i, atom_j, boundary, weight_14)
force(inter::SpecificInteraction, coord_i, coord_j,
      boundary)
force(inter::SpecificInteraction, coord_i, coord_j,
      coord_k, boundary)
force(inter::SpecificInteraction, coord_i, coord_j,
      coord_k, coord_l, boundary)

Calculate the force between atoms due to a given interation type. For PairwiseInteractions returns a single force vector and for SpecificInteractions returns a type such as SpecificForce2Atoms. Custom pairwise and specific interaction types should implement this function.

forces(system, neighbors=nothing; n_threads=Threads.nthreads())

Calculate the forces on all atoms in the system using the pairwise, specific and general interactions. If the interactions use neighbor lists, the neighbors should be computed first and passed to the function.

forces(inter, system, neighbors=nothing)

Calculate the forces on all atoms in the system arising from a general interaction. If the interaction uses neighbor lists, the neighbors should be computed first and passed to the function. Custom general interaction types should implement this function.

inject_gradients(sys, params_dic)

Add parameters from a dictionary to a System. Allows gradients for individual parameters to be tracked. Returns atoms, pairwise interactions, specific interaction lists and general interactions.

is_any_atom(at, at_data)

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

is_gpu_diff_safe(sys)

Whether a System or ReplicaSystem uses the code path suitable for the GPU and for taking gradients.

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)

Whether a System or ReplicaSystem is on the GPU.

kinetic_energy(s)

Calculate the kinetic energy of the system.

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

Log a property of the system thoughout 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(mass, temperature; rng=Random.GLOBAL_RNG)
maxwell_boltzmann(mass, temperature, k; rng=Random.GLOBAL_RNG)

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

place_atoms(n_atoms, boundary; min_dist=nothing, max_attempts=100)

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)

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.

potential_energy(s, neighbors=nothing)

Calculate the potential energy of the system using the pairwise, specific and general interactions. If the interactions use neighbor lists, the neighbors should be computed first and passed to the function. Not currently compatible with automatic differentiation using Zygote.

potential_energy(inter::PairwiseInteraction, vec_ij, coord_i, coord_j,
                 atom_i, atom_j, boundary)
potential_energy(inter::SpecificInteraction, coords_i, coords_j,
                 boundary)
potential_energy(inter::SpecificInteraction, coords_i, coords_j,
                 coords_k, boundary)
potential_energy(inter::SpecificInteraction, coords_i, coords_j,
                 coords_k, coords_l, boundary)
potential_energy(inter, system, neighbors=nothing)

Calculate the potential energy due to a given interation type. Custom interaction types should implement this function.

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.

rand_coord(boundary)

Generate a random coordinate uniformly distributed within a bounding box.

random_normal_translation!(sys::System; shift_size=oneunit(eltype(eltype(sys.coords))))

Performs a random translation of the coordinates of a randomly selected atom in a System. The translation is generated using a uniformly choosen 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))))

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)

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

random_velocities(sys, temp)

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

rdf(coords, boundary; npoints=200)

Calculate the radial distribution function of a set of coordinates. 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.

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

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.

run_loggers!(system, neighbors=nothing, step_n=0; n_threads=Threads.nthreads(), kwargs...)

Run the loggers associated with the system. Ignored for gradient calculation during automatic differentiation. Additional keyword arguments can be passed to the loggers if required.

simulate!(system, simulator, n_steps; n_threads=Threads.nthreads())
simulate!(system, simulator; n_threads=Threads.nthreads())

Run a simulation on a system according to the rules of the given simulator. Custom simulators should implement this function.

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

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

temperature(system)

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(s, neighbors=nothing)

Calculate the total energy of the system as the sum of the kinetic_energy and the potential_energy. If the interactions use neighbor lists, the neighbors should be computed first and passed to the function. Not currently compatible with automatic differentiation using Zygote.

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

velocity_autocorr(vl, first_ind, last_ind)

Calculate the autocorrelation function of velocity from the velocity logger. This characterizes the similarity between velocities observed at different time instances.

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

Visualize a simulation as an animation. This function is only available when GLMakie is imported. GLMakie v0.6 or later should be used. 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.

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.