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.

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.

Generalized Born (GB) implicit solvent models augmented with the hydrophobic solvent accessible surface area (SA) term.

Custom GBSA methods should sub-type this abstract type.

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.

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.
• atom_type::T: the type of the atom.
• 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, res_name)

Data associated with an atom.

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

AtomType(type, class, element, charge, mass, σ, ϵ)

An atom type.

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; 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 = 1 - \frac{\kappa_T \delta t}{3 \tau} ( P_0 - P )\]

with the fractional change limited to max_scale_frac.

This barostat should be used with caution as it 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.

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

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

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(
           eligible=s.neighbor_finder.eligible, dist_cutoff=1.2u"nm",
           special=s.neighbor_finder.special, n_steps=10,
           x0=coords, unit_cell=boundary,
       )
CellListMapNeighborFinder{Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}, 3, Float64}
  Size of eligible matrix = (15954, 15954)
  n_steps = 10
  dist_cutoff = 1.2 nm

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

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

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

CoulombSoftCore(; cutoff, α, λ, p, use_neighbors, σ_mixing, weight_special, coulomb_const)

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(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 the true potential and zero between an activation point and a cutoff point, using a cubic Hermite spline.

DistanceConstraint(i, j, dist)

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

DistanceCutoff(dist_cutoff)

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

DistanceNeighborFinder(; eligible, dist_cutoff, special, n_steps)

Find close atoms by distance.

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

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

Gravity(; 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\]

ImplicitSolventGBN2(atoms, atoms_data, bonds)

GBn2 solvation model implemented as an AtomsCalculators.jl calculator.

Should be used along with a Coulomb or CoulombReactionField interaction.

ImplicitSolventOBC(atoms, atoms_data, bonds)

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

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

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

LennardJonesSoftCore(; cutoff, α, λ, p, use_neighbors, shortcut, σ_mixing, ϵ_mixing,
                     weight_special)

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.

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)

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)

Topology information for a system.

Stores the index of the molecule each atom belongs to and the number of atoms in each molecule.

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. 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. Each should implement the AtomsCalculators.jl interface. Typically a Tuple.
• 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.

MonteCarloAnisotropicBarostat(pressure, temperature, boundary; 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 anisotropic barostat for controlling pressure.

For 3D systems, pressure is a SVector of length 3 with components pressX, pressY, and pressZ representing the target pressure in each axis. For 2D systems, pressure is a SVector of length 2 with components pressX and pressY. To keep an axis fixed, set the corresponding pressure to nothing.

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 in a randomly selected axis. The step is accepted or rejected based on

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

where ΔE is the change in potential energy, P is the equilibrium pressure along the selected axis, Δ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. If ΔW ≤ 0 the step is always accepted, if ΔW > 0 the step is accepted with probability exp(-ΔW/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.

MonteCarloBarostat(pressure, temperature, boundary; 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 W = \Delta E + P \Delta V - N k_B T \ln \left( \frac{V + \Delta V}{V} \right)\]

where ΔE is the change in potential energy, P is the equilibrium pressure, Δ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. If ΔW ≤ 0 the step is always accepted, if ΔW > 0 the step is accepted with probability exp(-ΔW/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.

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.

MonteCarloMembraneBarostat(pressure, tension, temperature, boundary; n_steps=30,
                   n_iterations=1, scale_factor=0.01, scale_increment=1.1,
                   max_volume_frac=0.3, trial_find_neighbors=false,
                   xy_isotropy=false, z_axis_fixed=false, constant_volume=false)

The Monte Carlo membrane barostat for controlling pressure.

Set the xy_isotropy flag to true to scale the x and y axes isotropically. Set the z_axis_fixed flag to true to uncouple the z-axis and keep it fixed. Set the constant_volume flag to true to keep the system volume constant by scaling the z-axis accordingly. The z_axis_fixed and constant_volume flags cannot be true simultaneously.

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 in a randomly selected axis. The step is accepted or rejected based on

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

where ΔE is the change in potential energy, P is the equilibrium pressure along the selected axis, ΔV is the change in volume, γ is the surface tension, ΔA is the change in surface area, N is the number of molecules in the system, T is the equilibrium temperature and V is the system volume. If ΔW ≤ 0 the step is always accepted, if ΔW > 0 the step is accepted with probability exp(-ΔW/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.

This barostat is only available for 3D systems.

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.

NeighborList(n, list)
NeighborList()

Structure to contain neighbor lists.

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

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

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(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_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.
• replica_coords: the coordinates of the atoms in each replica.
• boundary::B: the bounding box in which the simulation takes place.
• n_replicas::Integer: the number of replicas of the system.
• 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). 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=(): 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). 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=(): 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. 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 the same for all replicas). Typically a Tuple. 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.

RescaleThermostat(temperature)

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

ResidueType(name, types, charges, indices)

A residue type.

SHAKE_RATTLE(constraints, n_atoms, dist_tolerance, vel_tolerance)

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

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

Not currently compatible with GPU simulation.

Arguments

• constraints: a vector of constraints to be imposed on the system.
• n_atoms::Integer: the number of atoms in the system.
• dist_tolerance: the tolerance used to end the iterative procedure when calculating position constraints, should have the same units as the coordinates.
• vel_tolerance: the tolerance used to end the iterative procedure when calculating velocity constraints, should have the same units as the velocities * the coordinates.

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

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

Write 3D output structures to a file in the PDB format throughout a simulation.

The System should have atoms_data defined.

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. 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. Each should implement the AtomsCalculators.jl interface. Typically a Tuple.
• constraints::CN=(): the constraints for bonds and angles in the system. Typically a Tuple.
• 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(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 implicit_solvent, kappa and rename_terminal_res keyword arguments are 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.
• gpu::Bool=false: whether to move the relevant parts of the system onto the 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.
• center_coords::Bool=true: whether to center the coordinates in the simulation box.
• use_cell_list::Bool=true: whether to use CellListMapNeighborFinder on CPU. If false, DistanceNeighborFinder is used.
• data=nothing: arbitrary data associated with the system.
• implicit_solvent=nothing: specify a string to add an implicit solvent model, options are "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".

System(sys; <keyword arguments>)

Convenience constructor for changing properties in a System.

A copy of the System is returned with the provided keyword arguments modified.

System(abstract_system; force_units=u"kJ * mol^-1 * nm^-1", energy_units=u"kJ * mol^-1")

Convert an AtomsBase AbstractSystem to a Molly System.

force_units and energy_units should be set as appropriate. To add properties not present in the AtomsBase interface (e.g. pair potentials) use the convenience constructor System(sys::System).

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

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.

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 able to simulate a cubic box, use CubicBoundary or small offsets instead. Not currently compatible with infinite boundaries.

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.

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.

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

This should only be used on systems containing just pairwise interactions, or where the specific interactions, general interactions and constraints do not contribute to the pressure.

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

This should only be used on systems containing just pairwise interactions, or where the specific interactions, general interactions and constraints do not contribute to the virial.

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(sys), 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, coupling, simulator, neighbors=nothing,
                step_n=0; n_threads=Threads.nthreads())

Apply a coupler to modify a simulation.

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

apply_loggers!(system, neighbors=nothing, step_n=0, 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; n_threads::Integer=Threads.nthreads())
apply_position_constraints!(sys, coord_storage, vel_storage, dt;
                            n_threads::Integer=Threads.nthreads())

Applies the system constraints to the coordinates.

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

apply_velocity_constraints!(sys; n_threads::Integer=Threads.nthreads())

Applies the system constraints to the velocities.

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

charge(atom)

The partial charge of an Atom.

charges(sys)

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

check_position_constraints(sys, constraints)

Checks if the position constraints are satisfied by the current coordinates of sys.

check_velocity_constraints(sys, constraints)

Checks if the velocity constraints are satisfied by the current velocities of sys.

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.

disable_constrained_interactions!(neighbor_finder, constraint_clusters)

Disables neighbor list interactions between atoms in a constraint.

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,
               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(sys), step_n=0; n_threads=Threads.nthreads())

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

hydrodynamic_radius(coords, boundary)

Calculate the hydrodynamic radius of a set of coordinates.

$R_{hyd}$ is defined by

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

inject_gradients(sys, params_dic)

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

Calculate the kinetic energy of a system.

log_property!(logger, system, 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.

molecule_centers(coords, boundary, topology)

Calculate the coordinates of the center of each molecule in a system.

Accounts for periodic boundary conditions by using the circular mean. If topology=nothing then the coordinates are returned.

Not currently compatible with TriclinicBoundary if the topology is set.

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.

potential_energy(system, neighbors=find_neighbors(sys), 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(sys, neighbors=find_neighbors(sys), step_n=0; n_threads=Threads.nthreads())

Calculate the pressure of a system.

The pressure is defined as

\[P = \frac{1}{V} \left( NkT - \frac{2}{D} \Xi \right)\]

where V is the system volume, N is the number of atoms, k is the Boltzmann constant, T is the system temperature, D is the number of dimensions and Ξ is the virial calculated using virial.

This should only be used on systems containing just pairwise interactions, or where the specific interactions, constraints and general interactions without virial defined do not contribute to the virial. 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))))

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

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)
random_velocities!(vels, sys, temp)

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

random_velocities(sys, temp)

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.

remd_exchange!(sys, sim, n, m; rng=Random.default_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.

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, scale_factor; ignore_molecules=false)

Scale the coordinates and bounding box of a system 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. Velocities are not scaled. If the topology of the system is set then atoms in the same molecule will be moved by the same amount according to the center of coordinates of the molecule. This can be disabled with ignore_molecules=true.

Not currently compatible with TriclinicBoundary if the topology is set.

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

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

run_loggers 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. Custom simulators should implement this function.

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

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

Not currently compatible with interactions that depend on step number.

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(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(sys, neighbors=find_neighbors(sys), step_n=0; n_threads=Threads.nthreads())
virial(inter, sys, neighbors, step_n; n_threads=Threads.nthreads())

Calculate the virial of a system or the virial resulting from a general interaction.

The virial is defined as

\[\Xi = -\frac{1}{2} \sum_{i,j>i} r_{ij} \cdot F_{ij}\]

Custom general interaction types can implement this function.

This should only be used on systems containing just pairwise interactions, or where the specific interactions, constraints and general interactions without virial defined do not contribute to the 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.