Molly documentation

This documentation will first introduce the main features of the package with some examples, then will give details on each component of a simulation. There are further examples in the Molly examples section. For more information on specific types or functions, see the Molly API section or call ?function_name in Julia. The Differentiable simulation with Molly section describes taking gradients through simulations and the Free energy calculation section covers an approach to estimate free energies.

The package takes a modular approach to molecular simulation. To run a simulation you create a System object and call simulate! on it. The different components of the system and simulation can be used as defined by the package, or you can define your own versions. An important principle of the package is that your custom components, particularly force and potential energy functions, should be easy to define and just as performant as the built-in versions.

Simulation basics

Let's look at the simulation of a fluid acting under the Lennard-Jones potential to start with. First, we'll need some atoms with the relevant parameters defined.

using Molly

n_atoms = 100
atom_mass = 10.0u"g/mol"
atoms = [Atom(mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]

See the Unitful.jl documentation for more information on the unit annotations. Molly re-exports Unitful.jl and StaticArrays.jl since they are often required to run simulations. You can use your own atom types in Molly, provided that the mass function is defined and any fields required by the interactions are present. Next, we'll need some starting coordinates and velocities.

boundary = CubicBoundary(2.0u"nm") # Periodic boundary conditions with a 2 nm cube
coords = place_atoms(n_atoms, boundary; min_dist=0.3u"nm") # Random placement without clashing

temp = 100.0u"K"
velocities = [random_velocity(atom_mass, temp) for i in 1:n_atoms]

We store the coordinates and velocities as static arrays for performance. They can be of 2 or 3 dimensions and of any number type, e.g. Float64 or Float32. Setting individual dimensions in a CubicBoundary to Inf * u"nm" makes the simulation have no boundary in that dimension. You can also use a TriclinicBoundary. Simulations in 2 dimensions should use a RectangularBoundary.

The vector function calculates the vector from one coordinate to another accounting for periodic boundary conditions by using the minimum image convention:

vector(coords[1], coords[2], boundary)

3-element SVector{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}} with indices SOneTo(3):
 -0.8619698558762985 nm
 -0.9475020064484192 nm
  0.8359421827141784 nm

Now we can define our pairwise interactions, i.e. those between most or all atom pairs, as a tuple. Because we have defined the relevant parameters for the atoms, we can use the built-in LennardJones type. Pairwise interactions use a neighbor list by default, which needs a neighbor finder as shown in later examples. This small system does not need one, so we give use_neighbors=false to consider all pairs of atoms.

pairwise_inters = (LennardJones(use_neighbors=false),) # Don't forget the trailing comma!

Finally, we can define the system and run the simulation. We use an AndersenThermostat to keep a constant temperature, and we log the temperature and coordinates every 10 steps. Periodic boundary conditions are automatically used with the cubic box we defined earlier.

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    velocities=velocities,
    pairwise_inters=pairwise_inters,
    loggers=(
        temp=TemperatureLogger(10),
        coords=CoordinatesLogger(10),
    ),
)

simulator = VelocityVerlet(
    dt=0.002u"ps",
    coupling=AndersenThermostat(temp, 1.0u"ps"),
)

simulate!(sys, simulator, 1_000)

atoms, coords and boundary are the minimum required properties to define a System, though you would generally want to add interactions to a System to do something useful with it. The data keyword argument can give arbitrary data to the System that can be accessed with sys.data, for example metadata or properties to be used with a custom simulator. System implements the AbstractSystem interface from AtomsBase.jl in cases where units are provided. Showing sys in the REPL gives a summary of the system:

julia> sys
Molly System
  n atoms: 100
  n molecules: -
  n dimensions: 3
  n pairwise interactions: 1
  n specific interaction lists: 0
  n general interactions: 0
  n constrained atoms: 0
  n virtual sites: 0
  n loggers: 2
  neighbor finder: Nothing
  total mass: 1000.0 g mol^-1
  Degrees of freedom: 297
  Array type: Array
  Float type: general Float64, high Float64
  Units: energy kJ mol^-1, force kJ nm^-1 mol^-1, length nm, time ps, mass g mol^-1
  Boundary: CubicBoundary{3, Float64, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}}(Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}[2.0 nm, 2.0 nm, 2.0 nm])

Various functions can be used on a System:

potential_energy(sys) # -40.91077282719628 kJ mol^-1
kinetic_energy(sys)   # 119.47758705080862 kJ mol^-1
total_energy(sys)     #  78.56681422361234 kJ mol^-1

forces(sys)
accelerations(sys)
virial(sys)
pressure(sys)
dipole_moment(sys)

masses(sys)
density(sys) # 207.56738339673083 kg m^-3
temperature(sys) # 96.76667184796673 K
random_velocities(sys, 300.0u"K")

float_type(sys) # Float64
is_on_gpu(sys)  # false

# Access properties
sys.atoms
sys.coords
sys.boundary
sys.velocities
sys.topology
sys.pairwise_inters
sys.constraints
sys.virtual_sites
sys.virtual_site_flags
sys.neighbor_finder
sys.loggers
sys.total_mass

# Write out the system to a file, which requires atoms_data such as the element
sys_data = System(sys; atoms_data=[AtomData(element="Ar") for _ in 1:length(sys)])
write_structure("file.pdb", sys_data)
write_structure("file.mol2", sys_data)

# AtomsBase.jl interface
import AtomsBase
AtomsBase.mass(sys, 5) # 10.0 u
AtomsBase.position(sys, 10)

# Define a new system with certain properties changed
System(sys; coords=(sys.coords .* 0.5))

By default the simulation is run in parallel on the number of threads available to Julia, but this behaviour can be changed by giving the keyword argument n_threads to simulate!. For example, n_threads=1 uses no parallelization. n_threads=4 will limit the simulation to 4 threads. Sometimes setting n_threads higher than the number of available threads can improve performance.

The values stored by the loggers can be accessed using values, e.g. values(sys.loggers.coords). An animation of the stored coordinates can be saved by using visualize, which is available when GLMakie.jl is imported.

using GLMakie
visualize(sys.loggers.coords, boundary, "sim_lj.mp4")

A Molly System can also be constructed from an AtomsBase.jl system. For example:

System(ab_sys; force_units=u"kJ/Å", energy_units=u"kJ")

GPU acceleration

To run simulations on the GPU you will need to have a GPU available and then load the appropriate package:

Modern GPUs can run simulations of over 100,000 atoms, as seen in the examples. Metal/Apple Silicon devices can only run with 32 bit precision, so be sure to use Float32 in this case. Non-CUDA backends are less well-tested with Molly than CUDA. The tiled neighbor finder GPUNeighborFinder is CUDA-specific, so on other GPU backends the O(N) cell list GPUCellListNeighborFinder should be used instead.

Simulation setup is similar to above, but with the coordinates, velocities and atoms moved to the GPU. This example also shows setting up a simulation to run with Float32, which gives much better performance on GPUs. Of course, you will need to determine whether this level of numerical accuracy is appropriate in your case. Here is an example script for an NVIDIA GPU using CUDA:

using Molly
using CUDA

n_atoms = 100
atom_mass = 10.0f0u"g/mol"
boundary = CubicBoundary(2.0f0u"nm")
temp = 100.0f0u"K"
atoms = CuArray([Atom(mass=atom_mass, σ=0.3f0u"nm", ϵ=0.2f0u"kJ * mol^-1") for i in 1:n_atoms])
coords = CuArray(place_atoms(n_atoms, boundary; min_dist=0.3u"nm"))
velocities = CuArray([random_velocity(atom_mass, temp) for i in 1:n_atoms])
simulator = VelocityVerlet(dt=0.002f0u"ps")

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    velocities=velocities,
    pairwise_inters=(LennardJones(use_neighbors=false),),
    loggers=(
        temp=TemperatureLogger(typeof(1.0f0u"K"), 10),
        coords=CoordinatesLogger(typeof(1.0f0u"nm"), 10),
    ),
)

simulate!(deepcopy(sys), simulator, 20) # Compile function
simulate!(sys, simulator, 1_000)

To use another GPU package, just swap out CUDA for your desired package and CuArray for your desired array type. The device to run on can be changed with device!, e.g. device!(1). There are two GPU code paths currently: a fast path specific to CUDA and a slower path using KernelAbstractions.jl that is suitable for all backends.

The number of GPU threads used for the GPU kernels can be tuned with the environmental variables MOLLY_GPUNTHREADS_PAIRWISE, MOLLY_GPUNTHREADS_SPECIFIC, MOLLY_GPUNTHREADS_DISTANCENF, MOLLY_GPUNTHREADS_CELLLIST, MOLLY_GPUNTHREADS_COPY, MOLLY_GPUNTHREADS_IMPLICIT and MOLLY_GPUTILES_IMPLICIT. In general these should only be changed if GPU memory errors occur on smaller GPUs.

For the CUDA fast path, used with GPUNeighborFinder, users can explicitly call optimize_cuda_launch_config! on a system prior to a simulation. This will benchmark various launch configurations for the pairwise force and energy kernels and store the optimal parameters in the system, caching them for systems of the same type, size and box on the same device. Users can also manually override the kernel parameters by setting the environment variables MOLLY_CUDA_FORCE_BLOCK_Y, MOLLY_CUDA_ENERGY_BLOCK_Y, MOLLY_CUDA_TILE_THREADS_X, MOLLY_CUDA_TILE_THREADS_Y, and MOLLY_CUDA_FORCE_MAXREGS.

Simulating diatomic molecules

If we want to define specific interactions between atoms, for example bonds, we can do this as well. Using the same definitions as the first example, let's set up the coordinates so that paired atoms are 1 Å apart.

atoms = [Atom(mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]

coords = place_atoms(n_atoms ÷ 2, boundary; min_dist=0.3u"nm")
for i in 1:length(coords)
    push!(coords, coords[i] .+ [0.1, 0.0, 0.0]u"nm")
end

velocities = [random_velocity(atom_mass, temp) for i in 1:n_atoms]

We could have used place_diatomics instead here. Now we can use the built-in interaction list and bond types to place harmonic bonds between paired atoms.

bonds = InteractionList2Atoms(
    collect(1:(n_atoms ÷ 2)),           # First atom indices
    collect((1 + n_atoms ÷ 2):n_atoms), # Second atom indices
    [HarmonicBond(k=300_000.0u"kJ * mol^-1 * nm^-2", r0=0.1u"nm") for i in 1:(n_atoms ÷ 2)],
)

specific_inter_lists = (bonds,) # Don't forget the trailing comma!

This time we are also going to use a neighbor list to speed up the Lennard-Jones calculation since we don't care about interactions beyond a certain distance. We can use the built-in DistanceNeighborFinder. The arguments are the number of atoms, the pairs of atoms excluded from the non-bonded interactions and the distance cutoff to be classed as a neighbor. The neighbor list is reused until the atoms have moved far enough that a pair within the interaction cutoff could be missing from it, so we use a distance cutoff for the neighbor list that is larger than the cutoff for the interaction.

# All pairs apart from bonded pairs are eligible for non-bonded interactions
neighbor_finder = DistanceNeighborFinder(
    n_atoms=n_atoms,
    excluded_pairs=[(i, i + (n_atoms ÷ 2)) for i in 1:(n_atoms ÷ 2)],
    dist_cutoff=1.0u"nm",
)

cutoff = DistanceCutoff(0.8u"nm") # Less than half the box side
pairwise_inters = (LennardJones(cutoff=cutoff),)

Since there is a neighbor finder this time, we leave out use_neighbors=false so that the interaction uses the neighbor list. Now we can simulate as before.

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    velocities=velocities,
    pairwise_inters=pairwise_inters,
    specific_inter_lists=specific_inter_lists,
    neighbor_finder=neighbor_finder,
    loggers=(
        temp=TemperatureLogger(10),
        coords=CoordinatesLogger(10),
    ),
)

simulator = VelocityVerlet(
    dt=0.002u"ps",
    coupling=AndersenThermostat(temp, 1.0u"ps"),
)
simulate!(sys, simulator, 1_000)

This time when we view the trajectory we can add lines to show the bonds.

visualize(
    sys.loggers.coords,
    boundary,
    "sim_diatomic.mp4";
    connections=[(i, i + (n_atoms ÷ 2)) for i in 1:(n_atoms ÷ 2)],
)

The neighbors can be found using find_neighbors, e.g. find_neighbors(sys), which returns a NeighborList for the classical neighbor finders, a GPUCellListNeighborList for GPUCellListNeighborFinder and nothing for GPUNeighborFinder, whose CUDA kernels manage their tile list internally.

Simulating gravity

Molly is geared primarily to molecular simulation, but can also be used to simulate other physical systems. Let's set up a gravitational simulation. This example also shows a 2D simulation and no specified units.

atoms = [Atom(mass=1.0f0), Atom(mass=1.0f0)]
coords = [SVector(0.3f0, 0.5f0), SVector(0.7f0, 0.5f0)]
velocities = [SVector(0.0f0, 1.0f0), SVector(0.0f0, -1.0f0)]
pairwise_inters = (Gravity(use_neighbors=false, G=1.5f0),)
simulator = VelocityVerlet(dt=0.002f0)
boundary = RectangularBoundary(1.0f0)

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    velocities=velocities,
    pairwise_inters=pairwise_inters,
    loggers=(coords=CoordinatesLogger(Float32, 10; dims=2),),
    force_units=NoUnits,
    energy_units=NoUnits,
)

simulate!(sys, simulator, 2_000)

When we view the simulation we can use some extra options:

visualize(
    sys.loggers.coords,
    boundary,
    "sim_gravity.mp4";
    trails=4,
    framerate=15,
    color=[:orange, :lightgreen],
)

Simulating a protein

The recommended way to run a macromolecular simulation is to read in a force field in OpenMM XML format to a MolecularForceField and then read in a coordinate file in a format supported by Chemfiles.jl. Files for common force fields can be found at OpenMM and OpenMM force fields. This sets up a system in the same data structures as above and that is simulated in the same way. Here we carry out an energy minimization, simulate with a Langevin integrator in the NPT ensemble and use a TrajectoryWriter to write the trajectory as a DCD file (or another file format, by changing the file extension).

data_dir = joinpath(dirname(pathof(Molly)), "..", "data")
T = Float32 # Float32 is much faster on GPU
ff = MolecularForceField(
    joinpath(data_dir, "force_fields", "ff99SBildn.xml"),
    joinpath(data_dir, "force_fields", "tip3p_standard.xml"),
)

sys = System(
    joinpath(data_dir, "6mrr_equil.pdb"),
    ff;
    nonbonded_method=SetupPME(),
    loggers=(
        energy=TotalEnergyLogger(10),
        writer=TrajectoryWriter(10, "traj_6mrr_5ps.dcd"),
    ),
    array_type=Array, # CuArray for CUDA GPU
    float_type=T, # By default Float32 on GPU and Float64 on CPU
)

minimizer = SteepestDescentMinimizer()
simulate!(sys, minimizer)

temp = T(298.0)u"K"
random_velocities!(sys, temp)
simulator = Langevin(
    dt=T(0.001)u"ps",
    temperature=temp,
    friction=T(1.0)u"ps^-1",
    coupling=MonteCarloBarostat(T(1.0)u"bar", temp, sys.boundary),
)

simulate!(sys, simulator, 5_000)

The above 5 ps simulation looks something like this when you view it in PyMOL:

The system setup procedure is tested against OpenMM, following their template matching procedure to assign force field parameters to the structures read from the structure file. Some margin in residue and atom naming is allowed, as the naming present in the structure files is queried against a renaming dictionary that contains common alternative names present in PDB files, which you can consult in pdbNames.xml. You can extend this dictionary yourself, if you really want to use your own naming; or you can build an additional renaming dictionary following the same structure as the one mentioned above, and pass it as a keyword argument custom_renaming_scheme when you build your MolecularForceField. The bonding topology of the system is automatically inferred for standard residues (protein and nucleic acids, plus water). If your simulation contains other types of molecules, you must provide the topology for them. You can do this either by using a structure file format with explicit bond definitions, such as Mol2 or mmCIF, by defining the appropriate CONECT records in a PDB file, or by providing a custom topology template in a format equivalent to the one found in residues.xml.

If a residue in the structure file cannot be matched to a residue template then the error message names the residue, gives its atoms and diagnoses the most likely cause by comparing it to the closest template, for example:

could not match residue TRP (residue number 2 of chain "A") to any of the residue
templates in the force field. The set of atoms is similar to TRP, but the residue is
missing 1 C atom. The residue has 23 atoms: N, H, CA, ...

Common causes are missing heavy atoms or hydrogens, missing bond information for non-standard residues, non-standard atom or residue naming, and using a force field that does not cover the molecules in the file. A template with the same name as the residue is used if it matches, otherwise the templates are tried in alphabetical order; if more than one of them matches and they give different parameters then this is reported according to strictness and the first is used.

To run on the GPU, set array_type=GPUArrayType, where GPUArrayType is the array type for your GPU backend (for example CuArray for NVIDIA or ROCArray for AMD). The floating point type can be set with float_type. A boundary given as a keyword argument is not converted to this type, so on GPU, where Float32 is the default, give it in Float32, e.g. boundary=CubicBoundary(10.0f0u"nm"). Certain quantities such as the potential_energy and the virial are accumulated using a higher precision type, as are the large constant terms of Ewald and PME summation. This can be set with float_type_high but should generally be left as the default Float64. The nonbonded method can be selected using the nonbonded_method keyword argument to System. This can can be an instance of SetupCoulombReactionField (reaction field method), SetupPME (particle mesh Ewald summation), SetupEwald (Ewald summation, slow), SetupCoulomb or a cutoff like DistanceCutoff (short range only). Setting exclusion_correction=true also applies the reaction field to the excluded and special pairs and to each atom with itself, using the CoulombReactionFieldExclusion specific interaction. For a system set up without files, Molly.reaction_field_exclusions gives this list of interactions. To run with constraints, use the constraints (:none, :hbonds, :allbonds or :hangles) and rigid_water keyword arguments. SHAKE_RATTLE (constraint_algorithm=SetupSHAKE_RATTLE()) only supports clusters of constraints around a central atom, such as :hbonds and rigid water, so use the default SetupLINCS with :allbonds. Neither algorithm supports angle constraints that share atoms with other constraints, so :hangles currently fails for proteins and most small molecules. Note that rigid_water defaults to false, whereas OpenMM makes water rigid by default, so set rigid_water=true to reproduce OpenMM behavior. Hydrogen mass repartitioning can be used by setting for example hydrogen_mass=2, and is applied before constraints are generated. Unlike OpenMM, it is also applied to the hydrogens of rigid water.

You can use an implicit solvent method by giving the implicit_solvent keyword argument. The options are instances of SetupImplicitSolventOBC and SetupImplicitSolventGBN2, corresponding to the Onufriev-Bashford-Case GBSA model and the GB-Neck2 model. Other options detailed in the docstring for System include overriding the boundary dimensions in the file (boundary) and modifying the non-bonded interaction and neighbor list cutoff distances (dist_cutoff and dist_buffer). The strictness keyword argument determines behavior when encountering possible problems and can be set to :error or :nowarn rather than the default :warn. It can be set globally with the MOLLY_STRICTNESS environmental variable.

Harmonic position restraints can be added to a System for equilibration using Molly.add_position_restraints:

sys_res = Molly.add_position_restraints(
    sys,
    100_000.0u"kJ * mol^-1 * nm^-2";
    atom_selector=Molly.is_heavy_atom,
)

Supported OpenMM force field XML tags

See the OpenMM documentation for the available tags. The following tags are supported:

• <AtomTypes>: both atom types and atom classes are supported
• <Residues>: <VirtualSite> tags are supported except for type="localCoords", and the override attribute on <Residue> tags is supported
• <Patches>: patches that apply to one or two residue templates are supported, patches that apply to more than two residue templates and multiple single residue patches acting on one residue template are not supported
• <HarmonicBondForce>
• <HarmonicAngleForce>
• <PeriodicTorsionForce>: both <Proper> and <Improper> tags are supported
• <CMAPTorsionForce>
• <NonbondedForce>: <UseAttributeFromResidue> tags other than <UseAttributeFromResidue name="charge"/> are not supported, useDispersionCorrection is supported and is true by default
• <CustomTorsionForce>: the special case where energy="k*(theta-theta0)^2" is supported as it is used to define improper torsions in some force fields
• <CustomNonbondedForce>: the special case of the double exponential potential is supported with bondCutoff="3", global parameters alpha and beta, and per-particle parameters sigma and epsilon; its long-range dispersion correction, DExpDispersionCorrection, is off by default but can be turned on with dispersion_correction=true in System or the useDispersionCorrection attribute of <NonbondedForce>
• <LennardJonesForce>: <NBFixPair> tags and sigma14/epsilon14 attributes in <Atom> tags are supported, useDispersionCorrection is supported and is true by default
• <AmoebaUreyBradleyForce>: this defines the bond-like part of a general Urey-Bradley interaction and is not specific to the AMOEBA force field, note that k is equivalent to the k/2 harmonic bond factor
• <Include>

The following tags are not yet supported and in general will be ignored rather than throwing an error when reading in a MolecularForceField:

• <RBTorsionForce>
• <GBSAOBCForce>
• <CustomBondForce>
• <CustomAngleForce>
• <CustomTorsionForce>: apart from the special case described above
• <CustomNonbondedForce>: apart from the special case described above
• <CustomGBForce>
• <CustomHbondForce>
• <CustomManyParticleForce>
• <Script>: since this is not supported, XML files that do further setup with Python code in <Script> may not behave as expected
• Any of the other polarisable force field tags, i.e. <DrudeForce>, <HippoNonbondedForce> and those starting with Amoeba other than <AmoebaUreyBradleyForce>.

In general, custom forces should be implemented as described in Forces and energies.

Structure file formats

Coordinate files can be in any format where Chemfiles.jl can read the positions, atoms and bonds. Where Chemfiles.jl cannot read the residues, it is assumed that the whole file represents one residue. For example, water can be loaded from various formats:

using Molly, Test

data_dir = joinpath(dirname(pathof(Molly)), "..", "data")
ff = MolecularForceField(joinpath(data_dir, "force_fields", "tip3p_standard.xml"))

water_pdb  = System(joinpath(data_dir, "water_formats", "water.pdb" ), ff)
water_cif  = System(joinpath(data_dir, "water_formats", "water.cif" ), ff)
water_mol2 = System(joinpath(data_dir, "water_formats", "water.mol2"), ff)
water_sdf  = System(joinpath(data_dir, "water_formats", "water.sdf" ), ff) # Residue inferred for SDF

@test potential_energy(water_pdb ) ≈ potential_energy(water_cif) ≈
      potential_energy(water_mol2) ≈ potential_energy(water_sdf)

ReaxFF

ReaxFF is a reactive empirical force-field that allows the description of chemical bonding and reactive events without expensive quantum chemical calculations. (van Duin et al. 2001).

The reactive description is enabled through a continuous bond order formalism, where bond orders are calculated from interatomic distances and then corrected based on each atom's local coordination environment. A fixed bonding topology is therefore not required. Partial charges for all atoms are also updated each MD step using a dynamic charge equilibration scheme, allowing the charge distribution to respond to changes in geometry as bonds form and break.

Force-field parameters are assigned by chemical element and their combinations. Parameters sets are generally developed for particular chemical applications, so careful selection of a force-field file is important for the environments and reactions you intend to simulate. A selection of literature ReaxFF force-fields is provided in the ReaxFF force-field directory.

For implementation notes, see ReaxFF implementation details.

Simulating methane combustion

This section introduces the setup of ReaxFF MD simulations in Molly by adapting the procedure for high-temperature oxidation of methane in the canonical ensemble (NVT) as described in Chenoweth et al. (2008). This involves simulating one methane molecule in a box with 100 oxygen molecules using periodic boundary conditions, starting with low-temperature equilibration at 100 K and high-temperature equilibration at 2500 K, followed by 1 ns of production NVT simulation at 2500 K. This example walks through selecting a ReaxFF force-field file, constructing a Molly system compatible with ReaxFF simulation, minimizing the initial structure, controlling the temperature, and recording a reactive trajectory.

Computational setup

We first load Molly and CUDA and locate the bundled ReaxFF data, containing literature force-field files and example systems. Setting array_type = CuArray selects GPU execution; use Array for CPU execution. If running on CPU, the CUDA import can be omitted. We use Float64 to set internal ReaxFF calculations to double-precision, which is highly recommended as single-precision arithmetic can lead to significant energy drifts in long NVE simulations (see discussion in Senftle et al. (2016)).

using Molly
using CUDA

data_dir = joinpath(dirname(pathof(Molly)), "..", "data", "ReaxFF")
float_type = Float64
array_type = CuArray

Load the force field

We construct a ReaxFF interaction using the Chenoweth C/H/O parameter set for hydrocarbon oxidation. ReaxFF force fields in general are trained for highly specific applications, and should not be assumed to be transferable to other systems outside of their parameterization domain, even if they contain the same elements. The force-field must contain parameters for every element that will be present in the simulated system.

The force_field_file argument selects the parameter file, while float_type and array_type carry over our precision and CPU/GPU choices.

reaxff = ReaxFF(
    force_field_file=joinpath(data_dir, "force_fields", "Chenoweth_2008_CHO.ff"),
    float_type=float_type,
    array_type=array_type
)

Construct the system

The methane + oxygen system is loaded from the supplied PDB file. In principle any structure file that is supported by Chemfiles can be loaded, but they must contain element information for every atom since ReaxFF requires atomic coordinates and elements. In our case, the PDB file must contain valid element names in the column Element symbol (columns 77-78).

A periodic cubic box with side length 2.5 nm is chosen to give a reasonable system density. Molly sets up the charge equilibration and ReaxFF interactions using the selected force field by constructing a general interactions tuple with ElectronegativityEqualizationMethod and ReaxFF interactions attached.

eem_float_type selects the numerical precision of the EEM solver (Float64 highly recommended), while eem_tolerance controls its convergence tolerance (typically 1e-6 is used).

The loggers record total energy and write trajectory frames every 100 steps. These loggers run during minimization, equilibration, and production.

sys = System(
    joinpath(data_dir, "data", "methane_oxidation.pdb"),
    reaxff;
    boundary=CubicBoundary(float_type(2.5)u"nm"),
    array_type=array_type,
    eem_float_type=Float64,
    eem_tolerance=1e-6,
    loggers=(
        total_energy=TotalEnergyLogger(100),
        writer=TrajectoryWriter(100, "methane_oxidation_ReaxFF_NVT.dcd"; correction=:pbc),
    )
)

Minimize the initial structure

Before starting dynamics, we relax the initial structure using energy minimization via steepest-descent. Setting run_loggers=true includes the minimization energies and trajectory frames in the recorded output.

minimizer = SteepestDescentMinimizer()
simulate!(sys, minimizer; run_loggers=true)

Equilibrate at low temperature

Atomic velocities are initialized from the Maxwell–Boltzmann distribution at a temperature of 100 K, followed by 10 ps of molecular dynamics using velocity Verlet and a Berendsen thermostat. We note that the Berendsen thermostat does not generate the correct canonical distribution, but we use it following the procedure in the paper.

The timestep, dt=0.1u"fs", used in ReaxFF is typically small, with simulations reported in literature using the range 0.1 - 0.25 fs. We select the coupling time, 0.1u"ps", of the thermostat to match the settings in the Chenoweth et al. paper.

# low-temperature equilibration
temp = 100.0u"K"
random_velocities!(sys, temp)

simulator = VelocityVerlet(
    dt=0.1u"fs",
    coupling=BerendsenThermostat(temp, 0.1u"ps")
)
simulate!(sys, simulator, 100_000)      # 10 ps

Heat to the reaction temperature

We raise the thermostat target to 2500 K and equilibrate for 100 ps, retaining a timestep of 0.1 fs and thermostat coupling of 0.1 ps. Although Chenoweth et al. turned off C–O and H–O bond parameters during this equilibration section to prevent reactions, this example keeps the full force field active, so reactions can occur during equilibration.

# high-temperature equilibration
temp = 2500.0u"K"
simulator = VelocityVerlet(
    dt=0.1u"fs",
    coupling=BerendsenThermostat(temp, 0.1u"ps")
)
simulate!(sys, simulator, 1_000_000)    # 100 ps

Run production dynamics

We run production MD at 2500 K for 1 ns with the same 0.1 fs timestep. The thermostat coupling is increased to 0.5 ps to match the paper’s production settings, weakening the temperature coupling.

# start production
simulator = VelocityVerlet(
    dt=0.1u"fs",
    coupling=BerendsenThermostat(temp, 0.5u"ps")
)
simulate!(sys, simulator, 10_000_000)   # 1 ns

Inspect the output

The trajectory is written to methane_oxidation_ReaxFF_NVT.dcd in the working directory. The 100-step logging interval corresponds to saving a frame every 10 fs. The trajectory and energy histories also include minimization, equilibration, and production.

The recorded total energies can be accessed through the system’s logger:

energies = values(sys.loggers.total_energy)

Simulating a solvated protein

This example simulates solvated ubiquitin (PDB ID: 1UBQ) in the isothermal-isobaric ensemble (NPT) using the Monti et al. ProtReaxFF force-field (Monti et al. 2013).

It runs on a NVIDIA GPU at 300 K and 1 bar, using a Nosé–Hoover thermostat and an isotropic Monte Carlo barostat. Since ReaxFF does not currently implement the virial, the Monte Carlo barostat is required for pressure coupling. In addition, when using the Monte Carlo barostat without a neighbor finder skin, which is the default for ReaxFF, the trial_find_neighbors field must be set to true in order to refresh the neighbors when a trial attempt changes the atomic coordinates. With a skin given by neighbor_skin, the neighbor finder rebuilds the list only when needed, including for trial attempts.

using Molly
using CUDA

data_dir = joinpath(dirname(pathof(Molly)), "..", "data", "ReaxFF")
float_type = Float64
array_type = CuArray

reaxff = ReaxFF(
    force_field_file=joinpath(data_dir, "force_fields", "Monti_2013_CHONSMgPNaCuCl.ff"),
    float_type=float_type,
    array_type=array_type
)

sys = System(
    joinpath(data_dir, "data", "1UBQ_solvated.pdb"),
    reaxff;
    boundary=CubicBoundary(float_type(6.503)u"nm"),
    array_type=array_type,
    eem_float_type=Float64,
    eem_tolerance=1e-6,
    loggers=(
        total_energy=TotalEnergyLogger(1_000),
        writer=TrajectoryWriter(1_000, "1UBQ_ReaxFF_NPT.dcd"; correction=:pbc),
    )
)

minimizer = SteepestDescentMinimizer()
simulate!(sys, minimizer; run_loggers=true)

temp = 300.0u"K"
random_velocities!(sys, temp)

barostat = MonteCarloBarostat(
    1.0u"bar",
    temp,
    sys.boundary;
    n_steps=30,
    trial_find_neighbors=true,
)

simulator = NoseHoover(
    dt=0.25u"fs",
    temperature=temp,
    damping=0.1u"ps",
    coupling=barostat
)

simulate!(sys, simulator, 4_000_000)    # 1 ns

Furman and Wales smoothness correction

ReaxFF uses several cutoffs in its functional forms for bond-orders and distances to limit the number of interactions that require evaluating, therefore saving computation. However, abrupt truncation at these various cutoffs introduces discontinuities in the potential energy and forces, which can result in poor energy conservation and other issues. This implementation provides an optional flag to turn on the smoothness correction based on Furman and Wales (2019), which applies seventh-order polynomial tapering functions to bond orders and selected bonded energy contributions near their cutoffs.

Enable the correction by setting smoothness_correction=true when constructing the ReaxFF interaction, before creating the System. The default is false.

reaxff = ReaxFF(
    force_field_file=joinpath(data_dir, "force_fields", "Chenoweth_2008_CHO.ff"),
    float_type=float_type,
    array_type=array_type,
    smoothness_correction=true,
)

In this implementation, this smoothness correction smooths the uncorrected bond-order components and applies tapering to the valence-angle, penalty, conjugation, torsion, hydrogen-bond, and C–O triple-bond correction ReaxFF energy terms. The same force-field parameter file is used, but enabling the correction changes the evaluated energies and forces near the affected cutoffs. More details about the implementation of this correction can be found in the Furman-Wales smoothness correction section of the developer documentation.

Enhanced sampling

Molly has the ReplicaSystem struct and simulators such as ReplicaExchangeMD to carry out replica exchange molecular dynamics (REMD). On CPU these are run in parallel by dividing up the number of available threads. On GPU each device is driven by its own worker process from Distributed.jl, so add one worker per GPU with addprocs and load Molly and the GPU package on the workers with @everywhere before calling simulate!. All visible devices are used by default, and the gpu_devices keyword argument selects others, e.g. gpu_devices=collect(CUDA.devices())[1:1] with one worker. For example, to run temperature REMD on a protein with 4 replicas and attempt exchanges every 2.5 ps:

using Molly
using Statistics

data_dir = joinpath(dirname(pathof(Molly)), "..", "data")
ff = MolecularForceField(
    joinpath(data_dir, "force_fields", "ff99SBildn.xml"),
    joinpath(data_dir, "force_fields", "tip3p_standard.xml"),
)

sys = System(joinpath(data_dir, "6mrr_equil.pdb"), ff)

minimizer = SteepestDescentMinimizer()
simulate!(sys, minimizer)

dt            = 0.0005u"ps"
n_replicas    = 4
temps         = [240.0u"K", 280.0u"K", 320.0u"K", 360.0u"K"]
thermo_states = ThermoState[]

for temp in temps
        intg = Langevin(dt=dt, temperature=temp, friction=1.0u"ps^-1")
        push!(thermo_states, ThermoState(sys, intg; temperature=temp))
end

replica_loggers = [(temp=TemperatureLogger(10), coords=CoordinatesLogger(10)) for i in 1:n_replicas]

rep_sys = ReplicaSystem(
    thermo_states,
    [copy(sys.coords) for _ in 1:n_replicas];
    replica_loggers=replica_loggers,
)

sim = ReplicaExchangeMD(dt=dt, exchange_time=2.5u"ps")

simulate!(rep_sys, sim, 40_000; assign_velocities=true)

println(rep_sys.exchange_logger.n_attempts)
# 12

for id in 1:n_replicas
        mean_temp = mean(values(rep_sys.replica_loggers[id].temp)[(end - 10):end])
        # Given physical coordinates swap thermal states, they should average out across the ladder bounds
        println(mean_temp)
end
# 242.68166672295013 K
# 285.37217289488996 K
# 314.0608014652237 K
# 362.86427110336984 K

ReplicaSystem records the absolute MD step, and subsequent simulate! calls resume from rep_sys.current_step. Set initial_step in the constructor when restoring replica coordinates and velocities from a checkpoint.

The accelerated weight histogram method (AWHState and AWHSimulation) has also been implemented in Molly.jl, allowing enhanced sampling and on-the-fly estimation of free energies along arbitrary collective variables and alchemical transformations. A more detailed overview of this can be found in the Free energy calculation section.

As described in the next section, bias potentials can also be added to standard simulators for enhanced sampling.

Biased simulation

Molly allows users to bias simulations along one or several collective variables (CVs). A bias potential needs to be defined for every CV using the BiasPotential struct, which specifies the CV function and the functional form of the bias potential:

struct BiasPotential{C, B}
    cv_type::C
    bias_type::B
end

BiasPotential is a general interaction in Molly, and the potential energy and forces that result from the bias are calculated with AtomsCalculators.potential_energy and AtomsCalculators.forces! methods that take the BiasPotential struct as input. In the force calculation, the gradient of the bias potential with respect to the CV and the gradient of the CV function with respect to the system coordinates are calculated in two separate steps. Either calculation can be performed with an explicitly defined gradient function or with automatic differentiation.

The available CV types are:

• CalcDist, the distance between two sets of atoms. The way the distance is calculated is set by giving CalcMinDist, CalcMaxDist, CalcCMDist or CalcSingleDist as the dist_type argument.
• CalcRg, the radius of gyration of a set of atoms.
• CalcRMSD, the RMSD of a set of atoms to a target structure.
• CalcTorsion, the torsion angle defined by four atoms.

Currently, CV calculation is always done on the CPU. Other CV functions can be added by the user. Every CV type needs to have its own struct, an associated method of the calculate_cv function and potentially a method for the Molly.cv_gradient function.

To define your own CV function, first define the struct:

struct MyCV
    # Properties of the CV, e.g. indices of atoms over which the CV should be calculated
    correction::Symbol # Set to :pbc to make molecules whole, or :wrap to not
    has_virial::Bool   # Read when the virial is needed, e.g. with a barostat
end

If has_virial is true, a method of Molly.calculate_virial!(virial, cv::MyCV, coords, forces, atoms, boundary) that adds the contribution of the bias forces to the virial matrix also has to be defined. If it is false, the bias does not contribute to the virial.

Next, you need to define a method for the calculate_cv function. The method should look as follows, taking properties of the system rather than the system itself as input:

function Molly.calculate_cv(cv::MyCV, coords, atoms, boundary, velocities; kwargs...)
    # Function to calculate the CV value given a system configuration and properties of the CV
end

The gradient of calculate_cv is by default calculated with automatic differentiation if Enzyme is imported. However, it is also possible to manually add a method to the Molly.cv_gradient function to calculate the gradient without using automatic differentiation:

function Molly.cv_gradient(cv_type::MyCV, coords, atoms, boundary, velocities; kwargs...)
    # Function to calculate the gradient of the CV function with respect to the coordinates of the system
end

The system can be biased along the chosen CV using different bias potentials. The bias potential is a function of the system's current CV value and the target CV value and maps e.g. the difference between these two values to a potential energy. The available bias potential types are:

• LinearBias
• SquareBias
• FlatBottomSquareBias
• PeriodicFlatBottomBias, for periodic CVs such as torsion angles

You can define your own type of bias potential by first defining a new struct:

struct MyBias
    # Properties of the bias, including a force constant and the CV target value
end

The functional form of the bias potential is then defined by adding a method to the potential_energy function. The method should take cv_sim (the output of calculate_cv) as input:

function Molly.potential_energy(bias_fn::MyBias, cv_sim; kwargs...)
    pe = ... # Determined by the target value of the CV and cv_sim, the system's current CV value
    return pe
end

Finally, you need to add a method for MyBias to the Molly.bias_gradient function. The method must return the derivative of the potential_energy method for MyBias with respect to cv_sim:

function Molly.bias_gradient(bias_fn::MyBias, cv_sim; kwargs...)
    bias_grad = ... # Derivative of bias with respect to cv_sim, potentially calculated with autodiff
    return bias_grad
end

Both methods are given the keyword arguments of the force or energy call, such as step_n and visit_state (see General interactions).

Once the CV and the bias function have been defined, a BiasPotential can be constructed and added as a general interaction in the system setup, e.g. general_inters = (BiasPotential(cv_type, bias_type),). A system can be biased along multiple CVs simultaneously by adding several BiasPotentials to general_inters.

Below, we show an example of how to set up a simulation in which the distance between two atoms is biased in a simulation of a Lennard-Jones fluid:

using Molly
using Enzyme

n_atoms = 100
boundary = CubicBoundary(10.0u"nm")
temp = 298.0u"K"
atom_mass = 10.0u"g/mol"

atoms = [Atom(mass=atom_mass, σ=0.3u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]
coords = place_atoms(n_atoms, boundary; min_dist=0.3u"nm")
velocities = [random_velocity(atom_mass, temp) for i in 1:n_atoms]

pairwise_inters = (LennardJones(use_neighbors=false),)

# Bias distance between atoms 1 and 2
define_cv = CalcDist([1], [2])

# Harmonic bias potential with a target distance of 1.5 nm
define_bias = SquareBias(400u"kJ * mol^-1 * nm^-2", 1.5u"nm")

general_inters = (BiasPotential(define_cv, define_bias),)

simulator = VelocityVerlet(
    dt=0.002u"ps",
    coupling=AndersenThermostat(temp, 1.0u"ps"),
)

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    velocities=velocities,
    pairwise_inters=pairwise_inters,
    general_inters=general_inters,
)

simulate!(sys, simulator, 100_000)

See also this example.

Monte Carlo sampling

Molly has the MetropolisMonteCarlo simulator to carry out Monte Carlo sampling with Metropolis selection rates. For example, to perform simulated annealing on charged particles to form a crystal lattice:

n_atoms = 100
atoms = [Atom(mass=10.0u"g/mol", charge=1.0) for i in 1:n_atoms]
boundary = RectangularBoundary(4.0u"nm")

coords = place_atoms(n_atoms, boundary; min_dist=0.2u"nm")
pairwise_inters = (Coulomb(use_neighbors=false),)

temperatures = [1198.0, 798.0, 398.0, 198.0, 98.0, 8.0]u"K"
sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    pairwise_inters=pairwise_inters,
    loggers=(
        coords=CoordinatesLogger(n_atoms, dims=2),
        montecarlo=MonteCarloLogger(),
    ),
)

trial_args = Dict(:shift_size => 0.1u"nm")
for t in temperatures
    sim = MetropolisMonteCarlo(;
        temperature=t,
        trial_moves=Molly.random_uniform_translation!,
        trial_args=trial_args,
    )

    simulate!(sys, sim, 10_000)
end

println(sys.loggers.montecarlo.n_accept)
# 15234

visualize(sys.loggers.coords, boundary, "sim_montecarlo.gif")

trial_moves should be a function that takes a System as its argument and optional keyword arguments trial_args. It should modify the coordinates as appropriate, accounting for any boundary conditions. Molly.random_uniform_translation! and Molly.random_normal_translation! are provided as common trial move functions. MonteCarloLogger records various properties throughout the simulation.

Units

Molly is fairly opinionated about using Unitful.jl units as shown above: you don't have to use them, but it is better if you do. Any consistent unit scheme can be used, or no units at all. Molly is most strict about the mixture of molar and non-molar types. For example, if your energy and force units are molar then your atom masses should be g/mol or similar. If you are not using units then no quantities can have Unitful annotations and you are responsible for ensuring a consistent unit system. Whilst you occasionally may run into friction with dimension mismatches, using units has the major advantages of catching whole classes of errors and letting you physically interpret the numbers in your system. The performance overhead of using units is minimal. Units are not fully compatible with Enzyme, though simple cases may work.

All your interaction types need to return the same units of force and energy or the simulation will not run. By default these are kJ * mol^-1 * nm^-1 for force and kJ * mol^-1 for energy, but this can be changed using the force_units and energy_units arguments to System and some interactions. These arguments should be NoUnits if you are not using units. If you need to strip units for downstream analysis, use the ustrip or ustrip_vec functions. It should be noted that charges are stored as dimensionless, i.e. 1.0 represents an atomic charge of +1. It is possible that you may run into issues when using different but valid units of the same dimension together, e.g. 1.0u"nm" and 10.0u"Å". In this case, try using the same units throughout.

The float type of a System is read from the boundary by default and can be set with the float_type argument. The coordinates, velocities and boundary should all use this float type; mixing Float32 and Float64 data silently loses the precision you asked for on CPU and can fail to compile on GPU, so it is reported according to strictness. The parameters of atoms and interactions should also use this float type but are not checked. Note that some defaults are Float64, such as coulomb_const in Coulomb and CoulombReactionField, so give these explicitly when using Float32.

A System can be converted to another float type, to unitless values or to another device by giving float_type, units or array_type to the copy constructor. For example, a system simulated with units in Float32 on the GPU can be evaluated in Float64 without units on the CPU, which is suitable for Enzyme:

sys_cpu = System(sys; float_type=Float64, units=false, array_type=Array, grad_safe=true)

The atoms, coordinates, velocities, boundary, interactions, constraints, virtual sites, neighbor finder and Boltzmann constant are converted, and loggers are dropped with a warning unless given with the loggers keyword argument. Units are stripped without conversion, so the numbers are in the units each value is stored with. Converting to a higher precision does not recover the precision lost when the original values were stored. Interaction types are rebuilt from their converted fields by calling the constructor of the same name with the fields in order. Custom types where this does not work, for example because they contain buffers, should define a method of Molly.convert_component(conversion, x); see the existing methods in the Molly source for examples.

Atom types

Molly has a built-in Atom type with a few properties commonly used in molecular simulation defined. The mass and charge functions can be used on an Atom.

Custom atom types can be used just as effectively, and are the way to provide per-atom parameters that the built-in Atom does not have. A custom atom type needs to have:

• The mass, either as a mass field or by defining a method for mass.
• The charge, either as a charge field or by defining a method for charge. Only required by interactions that use charges.
• A field for every other atom parameter read by the interactions in the system. Molly.required_atom_fields(inter) gives these for a built-in interaction, for example (:σ, :ϵ, :λ) for LennardJones and (:charge,) for Coulomb.

For example, an atom type that works with LennardJones and Coulomb and carries an extra property:

struct MyAtom{T, M, S, E}
    mass::M
    charge::T
    σ::S
    ϵ::E
    λ::T
    hydrophobicity::T # Read by a custom interaction
end

Custom atom types should generally be bits types, i.e. isbitstype(MyAtom) should be true, to work on the GPU. The list of atoms passed to the System constructor should be concretely typed. To differentiate through a simulation, a custom atom type also needs a Base.zero method, which is used to allocate the gradient; differentiating with respect to the atom parameters themselves currently requires the built-in Atom type.

Additional non-bits type data for the atoms that is not directly used when calculating the interactions can be passed to the System constructor with the atoms_data keyword argument. For example the built-in AtomData type contains fields such as the atom name that are useful when writing trajectories. Setting up a System from a structure file and a force field always uses the built-in Atom type.

Forces and energies

Interactions define how different parts of the system interact. The force on each particle in the system is derived from the potential corresponding to the interaction.

\[\vec{F}_i = -\sum_j \frac{dV_{ij}(r_{ij})}{dr_{ij}}\frac{\vec{r}_{ij}}{r_{ij}}\]

In Molly there are three types of interactions:

• Pairwise interactions are present between all or most atom pairs, and account for example for non-bonded terms in molecular mechanics force fields. All pairwise interactions are a subtype of PairwiseInteraction.
• Specific interactions are present between specific atoms, and account for example for bonded terms in molecular mechanics force fields.
• General interactions are a free-form interaction type that can access the whole system and outputs forces for all atoms. This is useful for neural network potentials, implicit solvent models and other cases that require maximum flexibility. General interactions should be compatible with the AtomsCalculators.jl interface.

The available pairwise interactions are:

• LennardJones
• LennardJonesSoftCoreBeutler
• LennardJonesSoftCoreGapsys
• AshbaughHatch
• SoftSphere
• Mie
• Buckingham
• Coulomb
• CoulombSoftCoreBeutler
• CoulombSoftCoreGapsys
• CoulombReactionField
• CoulombSoftCoreBeutlerReactionField
• CoulombSoftCoreGapsysReactionField
• CoulombEwald
• CoulombSoftCoreBeutlerEwald
• CoulombSoftCoreGapsysEwald
• Yukawa
• DoubleExponential
• Gravity
• DPDInteraction

The available specific interactions (1-5 atoms) are:

• HarmonicPositionRestraint - 1 atom
• HarmonicBond - 2 atoms
• MorseBond - 2 atoms
• FENEBond - 2 atoms
• LennardJones14 - 2 atoms
• LennardJones14SoftCoreBeutler - 2 atoms
• LennardJones14SoftCoreGapsys - 2 atoms
• EwaldExclusion - 2 atoms
• CoulombReactionFieldExclusion - 2 atoms
• HarmonicAngle - 3 atoms
• CosineAngle - 3 atoms
• UreyBradley - 3 atoms
• PeriodicTorsion - 4 atoms
• RBTorsion - 4 atoms
• HarmonicTorsion - 4 atoms
• CMAPTorsion - 5 atoms

The available general interactions are:

• Ewald
• PME
• LJDispersionCorrection
• DExpDispersionCorrection
• ImplicitSolventOBC
• ImplicitSolventGBN2
• ANIPotential
• MullerBrown
• ASECalculator

In addition there are lambda-scaled variants for many of the above interactions (λ or Scaled suffixes). Some interactions combine instances of the above. For example, particle mesh Ewald summation uses the CoulombEwald pairwise interaction, the EwaldExclusion specific interaction and the PME general interaction, allowing the short range terms to use the neighbors. Similarly, the reaction field with exclusion_correction uses the CoulombReactionField pairwise interaction and the CoulombReactionFieldExclusion specific interaction for the excluded pairs and the self terms.

By default, functions like potential_energy and forces use the interactions in the system. However, different interactions can be passed as keyword arguments. This can be useful for isolating the contribution of individual interactions:

# Evaluate the potential energy of the system using the interactions in the system
potential_energy(sys)

# Evaluate the potential energy of the system using only the Lennard-Jones potential
potential_energy(
    sys;
    pairwise_inters=(LennardJones(),),
    specific_inter_lists=(),
    general_inters=(),
)

Pairwise interactions

Some pairwise interactions define mixing functions which determine how the parameters from each atom are combined. For example, the default σ_mixing for LennardJones is Molly.LorentzMixing, which calculates (atom_i.σ + atom_j.σ) / 2. Other mixing functions are available, such as Molly.GeometricMixing and Molly.WaldmanHaglerMixing. Custom mixing functions can be used instead by defining the struct and appropriate functions such as Molly.σ_mixing and Molly.ϵ_mixing. For example:

struct MyMixing end

Molly.σ_mixing(m::MyMixing, atom_i, atom_j, args...) = (atom_i.σ + atom_j.σ) / 2
Molly.ϵ_mixing(m::MyMixing, atom_i, atom_j, args...) = (atom_i.ϵ + atom_j.ϵ) / 2

Use args... as above since some interactions pass further arguments. The atom_type field of the atoms is available, allowing features like changing the weight of solute-solvent interactions. With GPUNeighborFinder only the atom fields used by an interaction are passed between GPU threads, so a custom mixing function that reads other fields, such as atom_type, has to list them, e.g. Molly.mixing_atom_fields(::MyMixing) = (:atom_type,). Shortcut functions are also available to skip the interaction entirely, with similar logic. Specific exceptions can be given for pairs of atom types by using Molly.MixingException and Molly.ExceptionList:

d = Dict((1, 2) => 0.3u"nm") # Pairs of atom type 1 and 2 have a σ value of 0.3 nm
exceptions = Molly.ExceptionList(d)
σ_mix = Molly.MixingException(MyMixing(), exceptions) # Give as e.g. σ_mixing argument to LennardJones

To define your own pairwise interaction, first define the struct, which must be a subtype of PairwiseInteraction:

struct MyPairwiseInter{C} <: PairwiseInteraction
    # Any properties, e.g. cutoffs and interaction parameters
    cutoff::C
    param::Float64
end

You can also define a use_neighbors method, which determines whether the neighbor list is used to omit distant atoms (true) or whether all atom pairs are always considered (false):

Molly.use_neighbors(inter::MyPairwiseInter) = false

This is true by default, in which case you must specify a neighbor finder when setting up the System. For built-in interactions this function accesses the use_neighbors field of the struct, which is set with the use_neighbors keyword argument and is true by default. To work on the GPU the struct should be a bits type, i.e. isbitstype(MyPairwiseInter) should be true.

Next, you need to define a method for the force and Molly.pairwise_force functions acting between a pair of atoms.

This has a set series of arguments:

function Molly.force(inter::MyPairwiseInter,
                     vec_ij,
                     atom_i,
                     atom_j,
                     force_units,
                     special,
                     coord_i,
                     coord_j,
                     boundary,
                     velocity_i,
                     velocity_j,
                     step_n,
                     args...)
    r = norm(vec_ij)
    params = (inter.param,)
    f = Molly.force_cutoff(inter.cutoff, inter, r, params)
    # Multiply the force magnitude by the normalised vector between the atoms
    fdr = f * vec_ij / r
    return fdr
end

function Molly.pairwise_force(::MyPairwiseInter, r, params)
    # Replace this with your force calculation
    # A positive force causes the atoms to move apart
    f = 0.0
    return f
end

The logic here is that any parameter mixing, shortcutting or setup is done in force then force_cutoff is used to apply the cutoff by calling Molly.pairwise_force, which calculates the force magnitude solely from interatomic distance r. params should not include any distance-based information, since Molly.pairwise_force may be called for other distances to apply the cutoff.

Most of the arguments to force will generally not be used but are passed to allow maximum flexibility. You can use args... to indicate unused further arguments, e.g. Molly.force(inter::MyPairwiseInter, vec_ij, args...). args... should be added to the end, since more arguments may be allowed in future. vec_ij is the vector between the closest images of atoms i and j accounting for the periodic boundary conditions. Atom properties can be accessed, e.g. atom_i.σ. force_units can be useful for returning a zero force under certain conditions. If you want to use velocity_i and velocity_j you should set Molly.pairwise_uses_velocity(::MyPairwiseInter) = true, otherwise the velocities are passed as nothing for speed. step_n is the step number in the simulator, allowing time-dependent interactions. Beware that this step counter starts from 1 every time simulate! is called unless you give the init_step argument, and can also be 0 to calculate forces before the first step. It also doesn't work with simulate_remd!.

Typically the force function is where most computation time is spent during the simulation, so consider optimising this function if you want high performance. One nice feature of Molly is that this function will work on both the CPU and the GPU. If you need a different version of the function on GPU, you can define Molly.force_gpu and Molly.potential_energy_gpu. The returned values should have the float type of the system. Virial computation is done automatically when required using the force function. The force function should not have side effects like updating atom properties to avoid race conditions; use a general interaction if you need this.

The argument special is a Bool determining whether the atom pair interaction should be treated as special. This is specified during neighbor finder construction. When simulating molecules, for example, non-bonded interactions for atoms in a 1-4 bonding arrangement (i-x-x-j) are often weighted by a factor such as 0.5. For interactions where this is relevant, special can be used to apply this weighting after calling force_cutoff in the force function. It can have a variety of uses depending on the context, for example if you have multiple interactions and want to exclude certain atom pairs from one of the interactions only.

To use your custom interaction in a simulation, add it to the list of pairwise interactions:

pairwise_inters = (MyPairwiseInter(NoCutoff(), 1.0),) # Don't forget the trailing comma!

Then create a System and simulate as above. Note that you can also use a named tuple instead of a tuple if you want to access interactions by name:

pairwise_inters = (MyPairwiseInter=MyPairwiseInter(NoCutoff(), 1.0),)

If you wish to calculate potential energies or log the energy throughout a simulation, you will need to define methods for the potential_energy and Molly.pairwise_pe functions. These have the same arguments and logic as force, except the fifth argument to potential_energy is the energy units not the force units, and should return a single value corresponding to the potential energy:

function Molly.potential_energy(inter::MyPairwiseInter,
                                vec_ij,
                                atom_i,
                                atom_j,
                                energy_units,
                                special,
                                coord_i,
                                coord_j,
                                boundary,
                                velocity_i,
                                velocity_j,
                                step_n,
                                args...)
    r = norm(vec_ij)
    params = (inter.param,)
    return Molly.pe_cutoff(inter.cutoff, inter, r, params)
end

function Molly.pairwise_pe(::MyPairwiseInter, r, params)
    # Replace this with your potential energy calculation
    pe = 0.0
    return pe
end

Specific interactions

To define your own specific interaction, first define the struct:

struct MySpecificInter
    # Properties, e.g. a bond distance corresponding to the energy minimum
end

Next, you need to define a method for the force function. The form of this will depend on whether the interaction involves 1, 2, 3, 4 or 5 atoms. For example in the 2 atom case:

function Molly.force(inter::MySpecificInter,
                     coord_i,
                     coord_j,
                     boundary,
                     atom_i,
                     atom_j,
                     force_units,
                     velocity_i,
                     velocity_j,
                     step_n,
                     data)
    dr = vector(coord_i, coord_j, boundary)

    # Replace this with your force calculation
    # A positive force causes the atoms to move apart
    f = 0.0

    fdr = f * normalize(dr)
    return SpecificForce2Atoms(-fdr, fdr)
end

Again, most of these arguments are rarely used and can be replaced with args.... The 3 atom case would define Molly.force(inter::MySpecificInter, coord_i, coord_j, coord_k, boundary, atom_i, atom_j, atom_k, force_units, velocity_i, velocity_j, velocity_k, step_n, data) and return SpecificForce3Atoms(f1, f2, f3). Virial computation is done automatically when required using the force function.

To use your custom interaction, add it to the specific interaction lists along with the atom indices:

specific_inter_lists = (
    InteractionList2Atoms(
        [1, 3],
        [2, 4],
        [MySpecificInter(), MySpecificInter()],
        ["", ""], # Interaction type information, optional
        nothing,  # Interaction data, optional
    ), # Don't forget the trailing comma!
)

Giving interaction data, including arrays, as the last argument to the interaction list means that it can be accessed in the force and potential_energy functions via the data argument. Interaction lists can be indexed and iterated over, with each entry being a NamedTuple of the atom indices, the interaction and the interaction type. For 3 atom interactions use InteractionList3Atoms and pass 3 sets of indices. If using the GPU, the inner list of indices and interactions should be moved to the GPU with CuArray. The number in the interaction list and the return type from force must match, e.g. InteractionList3Atoms must always return SpecificForce3Atoms from the corresponding force function. If some atoms are required in the interaction for force calculation but have no force applied to them by the interaction, give a zero force vector for those atoms. Again a method for potential_energy with the same arguments, except the seventh argument is the energy units not the force units, can be defined:

function Molly.potential_energy(inter::MySpecificInter,
                                coord_i,
                                coord_j,
                                boundary,
                                atom_i,
                                atom_j,
                                energy_units,
                                velocity_i,
                                velocity_j,
                                step_n,
                                data)
    # Example harmonic bond interaction
    dr = vector(coord_i, coord_j, boundary)
    r = norm(dr)
    return (inter.k / 2) * (r - inter.r0) ^ 2
end

General interactions

To define your own general interaction, first define the struct:

struct MyGeneralInter
    # Properties, e.g. a neural network model
end

Next, you need to define a method for the AtomsCalculators.forces! function (note this is different to the force function above, and to the forces function).

import AtomsCalculators

function AtomsCalculators.forces!(fs,
                                  sys,
                                  inter::MyGeneralInter;
                                  neighbors=nothing,
                                  step_n=0,
                                  n_threads=Threads.nthreads(),
                                  kwargs...)
    # kwargs... is required, neighbors/step_n/n_threads can be omitted if not used

    # Calculate the forces on all atoms using the interaction and the system
    # Add the forces to the existing forces, which have the same shape as the coordinates
    # For example, a neural network might do something like this
    fs .+= inter.model(sys.coords, sys.atoms)
    return fs
end

The neighbors calculated from the neighbor list are available in this function, but may or may not be used depending on the context. You could carry out your own neighbor finding in this function if required. Since the forces are available for mutation in the function, general interactions can be used to modify the forces in arbitrary ways, though it is up to you to write an appropriate potential energy function if energies are required. In cases where the forces are not simply added to, the order of general interactions can matter. General interactions are applied after pairwise and specific interactions by default. Define Molly.evaluate_before_pairwise(::MyGeneralInter) = true to evaluate an interaction before them. On the GPU, a general interaction that only adds to the forces, as above, can define Molly.concurrent_with_pairwise(::MyGeneralInter) = true so that it is calculated at the same time as the pairwise and specific interactions, which is faster. Note that this function calculates forces not accelerations; if you have a neural network that calculates accelerations you should multiply these by the masses of the system, masses(sys), to get the forces according to F=ma.

If you want to use the virial in your simulation, for example for to control pressure, you will also need to calculate the contribution to the virial in the AtomsCalculators.forces! function. To do this, the function should take the additional keyword arguments needs_vir and buffers that will be passed by the simulator. needs_vir is a Bool that determines whether the virial is needed that call and buffers.virial contains the virial tensor. For example, the following could go in the function above:

    if needs_vir
        # The virial is a 3x3 SMatrix with units of energy
        buffers.virial .+= rand(SMatrix{3, 3, Float64}) * sys.energy_units
    end

A method for the AtomsCalculators.potential_energy function that takes the same arguments and returns a single value can also be defined. You can add the AtomsCalculators.@generate_interface annotation to these functions as per the AtomsCalculators.jl docs.

A general interaction that changes with the simulation history, such as Metadynamics, can use the visit_state keyword argument. This is given to forces and potential_energy and defaults to false. It is true for a molecular dynamics simulator force call, and false for minimizers and loggers, so that only visited states update the interaction. MetropolisMonteCarlo never passes visit_state=true, and ReplicaExchangeMD, AWHSimulation and TSSSimulation evaluate the energies of other thermodynamic states with frozen copies of the interactions, so none of these are compatible with such interactions.

To use your custom interaction in a simulation, add it to the list of general interactions:

general_inters = (MyGeneralInter(),) # Don't forget the trailing comma!

general_inters=general_inters can be given as a keyword argument when setting up the System.

Cutoffs

As most atomic interactions decay rapidly with the interatomic distance, a cutoff radius beyond which the interaction is ignored is often used. There are a number of approaches to avoid discontinuities in the force/potential energy at the cutoff distance. Since the cutoff strategy is independent of the interaction for which is used, each interaction is defined without including a specific cutoff. See Pairwise interactions for how to make a pairwise interaction compatible with arbitrary cutoffs. The corresponding interaction constructor has a cutoff field (default NoCutoff for the built-in interactions) which is then used via dispatch to apply the chosen cutoff, e.g. SoftSphere(cutoff=ShiftedPotentialCutoff(1.2u"nm")).

The available cutoffs are:

• NoCutoff
• DistanceCutoff
• ShiftedPotentialCutoff
• ShiftedForceCutoff
• CubicSplineCutoff
• PolynomialCutoff

The following built-in interactions (and their lambda-scaled variants) can use a cutoff:

• LennardJones
• LennardJonesSoftCoreBeutler
• LennardJonesSoftCoreGapsys
• AshbaughHatch
• SoftSphere
• Mie
• Buckingham
• Coulomb
• CoulombSoftCoreBeutler
• CoulombSoftCoreGapsys
• Yukawa
• DoubleExponential
• Gravity

In addition, CoulombReactionField, CoulombEwald and the implicit solvent and Ewald general interactions have a dist_cutoff argument for a cutoff distance. Cutoffs combine well with neighbor finders, provided that the cutoff distance is less than the neighbor distance (see Neighbor finders for more).

To define a custom cutoff, you should first decide whether it has one cutoff point (the force/potential energy is modified up to the cutoff distance) or two (the force/potential energy is modified between an activation distance and the cutoff distance). Then, create a struct that subtypes Molly.AbstractCutoff and cutoff functions:

struct MyCutoff{P, D} <: Molly.AbstractCutoff{1} # One cutoff point
    dist_cutoff::D # This field is required
end

MyCutoff(dist_cutoff::D) where D = MyCutoff{1, D}(dist_cutoff)

function Molly.pe_apply_cutoff(::MyCutoff, inter, r, params)
    # Modify the potential energy, pairwise_pe with other values of r can be used
    pe_r = Molly.pairwise_pe(inter, r, params)
    return pe_r
end

function Molly.force_apply_cutoff(::MyCutoff, inter, r, params)
    # Modify the force magnitude, pairwise_force with other values of r can be used
    f_r = Molly.pairwise_force(inter, r, params)
    return f_r
end

The two cutoff point case is similar except it should sub-type Molly.AbstractCutoff{2} and the field dist_activation is also required in the struct.

Boundaries

Molly allows the use of various periodic and infinite boundary conditions. The available 3D boundaries are:

• CubicBoundary
• TriclinicBoundary

The available 2D boundaries are:

• RectangularBoundary

Some examples of using boundaries:

CubicBoundary(2.0u"nm"                             ) # Periodic cube with 2 nm sides
CubicBoundary(2.0u"nm"   , 2.0u"nm"   , 2.0u"nm"   ) # Periodic cube with 2 nm sides
CubicBoundary(4.0u"nm"   , 5.0u"nm"   , 6.0u"nm"   ) # Periodic cuboid
CubicBoundary(2.0u"nm"   , 2.0u"nm"   , Inf * u"nm") # Infinite boundary in z direction
CubicBoundary(Inf * u"nm", Inf * u"nm", Inf * u"nm") # Infinite boundary, no periodicity
CubicBoundary(Inf * u"nm"                          ) # Infinite boundary, no periodicity

RectangularBoundary(2.0u"nm"             ) # Periodic square
RectangularBoundary(4.0u"nm", 5.0u"nm"   ) # Periodic rectangle
RectangularBoundary(2.0u"nm", Inf * u"nm") # Infinite boundary in y direction

# Periodic triclinic from basis vectors
TriclinicBoundary(SVector(
    SVector(2.2      , 0.0      , 0.0      )u"nm",
    SVector(1.0      , 1.7320508, 0.0      )u"nm",
    SVector(1.37888  , 0.5399122, 1.0233204)u"nm",
))

# Periodic triclinic from basis vector lengths and angles α/β/γ
b = TriclinicBoundary(
    SVector(2.2, 2.0, 1.8)u"nm",
    deg2rad.(SVector(50.0, 40.0, 60.0)),
)

# Volume of bounding box
volume(b) # 3.8993746318188633 nm^3

# Random coordinate uniformly distributed within boundary
random_coord(b) # SVector(2.651062310435411, 2.1702306804433973, 0.9518105403051831)u"nm"

# Wrap a coordinate back into the boundary if it is outside
wrap_coords(SVector(1.0, 1.0, 1.0)u"nm", b) # SVector(3.2, 1.0, 1.0)u"nm"

The Molly.box_center, AtomsBase.n_dimensions, float_type, place_atoms and place_diatomics functions are also available for boundaries.

The appropriate boundary to use will depend on your simulation. For example, having different lengths in each dimension would usually only make sense in a situation where forces or restraints depended on the dimension.

Simulators

Simulators define what type of simulation is run. This could be anything from a simple energy minimization to complicated replica exchange MD. The available simulators are:

• SteepestDescentMinimizer
• FIREMinimizer
• LBFGSMinimizer
• VelocityVerlet
• DPDVelocityVerlet
• Verlet
• StormerVerlet
• Langevin
• LangevinSplitting
• OverdampedLangevin
• NoseHoover
• MTSIntegrator
• MTSLangevinIntegrator
• ReplicaExchangeMD
• MetropolisMonteCarlo

Many of these require a time step dt as an argument. Many also remove the center of mass motion every time step, which can be tuned with the remove_CM_motion argument (false or a number of steps). The third argument to simulate! can be an Integer number of steps or a simulation time, e.g. simulate!(sys, simulator, 10.0u"ns"), in which case the number of steps is calculated from dt.

The statistical ensemble that is sampled depends on the simulator and on any Coupling that is applied:

The available thermostats and barostats are listed under Coupling. Not every coupler samples the corresponding ensemble correctly, for example the BerendsenThermostat and BerendsenBarostat do not, so read the docstrings before using them.

Common options when calling simulate! with these simulators include:

• show_progress to decide whether to show a progress bar for the simulation. This is true by default in the REPL/IJulia/Pluto, otherwise false, and can be set globally with the environmental variable MOLLY_SHOW_PROGRESS.
• shortcut to provide a struct that has a Molly.shortcut_sim method defined, checking at the end of each step when to stop the simulation early.
• init_step to specify the step number before the first step is taken, useful for time-dependent potentials. Default 0.

The LangevinSplitting simulator can be used to define a variety of integrators such as velocity Verlet (splitting "BAB"), the Langevin implementation in Langevin ("BAOA"), and symplectic Euler integrators ("AB" and "BA").

The NoseHoover simulator gives the energy of its thermostats to the loggers as the heat_bath_energy keyword argument, so the quantity conserved by the integrator can be logged to check the simulation:

function conserved_energy(sys, neighbors, step_n, buffers; heat_bath_energy, kwargs...)
    return kinetic_energy(sys) + potential_energy(sys, neighbors, step_n) + heat_bath_energy
end

loggers = (conserved=GeneralObservableLogger(conserved_energy, typeof(1.0u"kJ * mol^-1"), 100),)

To define your own simulator, first define a struct:

struct MySimulator
    # Any properties, e.g. the time step or coupling methods
end

Then, define a method for simulate! that carries out the simulation. This example shows some helper functions that you can use:

function Molly.simulate!(sys,
                         sim::MySimulator,
                         n_steps::Integer;
                         n_threads::Integer=Threads.nthreads(),
                         run_loggers=true)
    # Find neighbors like this
    neighbors = find_neighbors(sys, sys.neighbor_finder; n_threads=n_threads)

    # Storage used by the force calculation, couplers and loggers
    buffers = Molly.init_buffers!(sys, n_threads)

    # Calculate accelerations like this
    # Includes appropriate unit conversion for when the force units are per mol
    # See src/simulators.jl for more efficient code that reuses memory
    # visit_state=true lets history-dependent interactions such as Metadynamics update
    accels_t = accelerations(sys, neighbors, 0; n_threads=n_threads, visit_state=true)

    for step_n in 1:n_steps
        # Example velocity Verlet update
        sys.coords .+= sys.velocities .* sim.dt .+ accels_t .* sim.dt^2 ./ 2

        # Ensure coordinates stay within the simulation box like this
        sys.coords .= wrap_coords.(sys.coords, (sys.boundary,))

        # Find new neighbors like this, after the coordinates change and before
        #   anything uses them
        # The neighbor finder decides whether to rebuild the list based on how far the
        #   atoms have moved, or on the step number with auto_update=false
        neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n;
                                   n_threads=n_threads)

        accels_t_dt = accelerations(sys, neighbors, step_n; n_threads=n_threads,
                                    visit_state=true)
        sys.velocities .+= (accels_t .+ accels_t_dt) .* sim.dt ./ 2
        accels_t = accels_t_dt

        # Remove center of mass motion like this
        remove_CM_motion!(sys)

        # Apply coupling like this
        recompute_forces = Molly.apply_coupling!(sys, buffers, sim.coupling, sim, neighbors,
                                                 step_n; n_threads=n_threads)

        # A coupler such as a barostat can change the coordinates, in which case the
        #   neighbors are found again regardless of the step number
        # A neighbor finder with auto_update=true would also notice the change without
        #   force_recompute=true, and only rebuild the list if needed
        if recompute_forces
            neighbors = find_neighbors(sys, sys.neighbor_finder, neighbors, step_n, true;
                                       n_threads=n_threads)
            accels_t = accelerations(sys, neighbors, step_n; n_threads=n_threads)
        end

        # Apply the loggers like this
        # Computed quantities can also be given as keyword arguments to apply_loggers!
        Molly.apply_loggers!(sys, neighbors, step_n, buffers, run_loggers; n_threads=n_threads)
    end

    return sys
end

See more in the source code, for example how to apply constraints to coordinates and velocities. To use your custom simulator, give it as the second argument when calling simulate!.

To define your own replica exchange simulator, first define a struct:

struct MyREMDSimulator
    dt # Time step
    exchange_time # Time between exchanges
    simulators # A list of simulators to use for each replica e.g. Langevin
    # Other properties of the simulation
end

Then define the function that carries out the exchange, Molly.remd_exchange!:

function Molly.remd_exchange!(sys::ReplicaSystem,
                              sim::MyREMDSimulator,
                              n::Integer,
                              m::Integer;
                              n_threads::Integer=Threads.nthreads(),
                              rng=Random.default_rng())
    # Attempt to exchange the replicas with index n and m
    # First define Δ for the REMD scheme
    make_exchange = Δ <= 0 || rand(rng) < exp(-Δ) # Example of Metropolis acceptance rate
    if make_exchange
        # Exchange coordinates and velocities of replicas
        # Also scale the velocities to match the temperature if required
    end

    return Δ, make_exchange
end

The above function returns Δ, the argument of the acceptance rate that is logged by ReplicaExchangeLogger, and a boolean indicating whether the exchange was successful.

Then, define a method for the simulate! function to perform the parallel simulation. This does any initial setup such as assigning velocities then uses simulate_remd! to run the simulation:

function Molly.simulate!(sys::ReplicaSystem,
                         sim::MyREMDSimulator,
                         n_steps::Integer;
                         rng=Random.default_rng(),
                         n_threads::Integer=Threads.nthreads())
    # Do any initial setup if necessary
    simulate_remd!(sys, sim, n_steps; rng=rng, n_threads=n_threads)
end

Under the hood there are multiple implementations for the forces function, used by accelerations, and for potential_energy: a version geared towards CPUs and parallelism, a CUDA version, and a version for other GPU backends. You can define different versions of a simulator for CPU, CUDA and generic GPU systems by dispatching on System{D, Array} or System{D, CuArray} and System{D, AT} where AT <: AbstractGPUArray respectively. This also applies to coupling methods, neighbor finders and analysis functions. You do not have to define different versions though: you may only intend to use the simulator one way, or one version may be performant in all cases.

Coupling

Some simulators can be modified by adding coupling methods to allow properties like temperature and pressure to be controlled during a simulation. The available couplers are:

• ImmediateThermostat
• VelocityRescaleThermostat
• AndersenThermostat
• BerendsenThermostat
• BerendsenBarostat
• CRescaleBarostat
• MonteCarloBarostat

Currently the VelocityVerlet, Verlet, Langevin and NoseHoover simulators support coupling methods, with the default being nothing. StormerVerlet does not support coupling methods because its recurrence stores coordinate history rather than advancing velocities as primary state. Couplers are given to the coupling keyword argument during simulator construction:

temp = 300.0u"K"
press = 1.0u"bar"
thermostat = AndersenThermostat(temp, 1.0u"ps")
barostat = MonteCarloBarostat(press, temp, sys.boundary)

# Velocity Verlet with Andersen thermostat
VelocityVerlet(dt=0.001u"ps", coupling=thermostat)

# Velocity Verlet with Andersen thermostat and Monte Carlo barostat
VelocityVerlet(dt=0.001u"ps", coupling=(thermostat, barostat))

The appropriate coupling to use will depend on the situation. For example, the Langevin simulator controls temperature so does not require a thermostat. The MonteCarloBarostat for controlling pressure assumes 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 the AndersenThermostat coupling. In general, something should be implemented as a coupling method if it can be applied "on top of" a standard simulator, but should be its own simulator if it requires a change to the integration scheme or what is integrated (e.g. NoseHoover).

To define your own coupling method, first define the struct:

struct MyCoupler
    # Any properties, e.g. a target temperature or coupling constant
end

Then, define the function that implements the coupling every time step:

function Molly.apply_coupling!(sys, buffers, coupling::MyCoupler, sim, neighbors, step_n;
                               n_threads=Threads.nthreads(), rng=Random.default_rng(), kwargs...)
    # Do something to the simulation, e.g. scale the velocities
    # Return whether the coupling has invalidated the currently stored forces,
    #   for example by changing the coordinates
    recompute_forces = false
    return recompute_forces
end

The functions random_velocity, maxwell_boltzmann and temperature may be useful here. To use your custom coupler, give it as the coupling argument to the simulator as above.

Note that as some coupling methods need the virial to work, a custom coupler can be accompanied by a function that tells Molly the interval of simulation steps at which the virial is required. This helps save computation. For example:

# The virial is needed every c.n_steps
Molly.needs_virial(c::MyCoupler) = c.n_steps

# The virial is not needed, this is the default
Molly.needs_virial(c::MyCoupler) = Inf

The use of the virial tensor allows for non-isotropic pressure control. Molly follows the definition in LAMMPS, $\bf{W} = \sum_i \bf{r_i} \otimes \bf{f_i}$, taking into account pairwise and specific interactions (of two or more atoms), the contribution of the Ewald and PME methods, dispersion corrections and bias potentials. The virial tensor of GROMACS is $-\bf{W} / 2$, but the pressure tensor has the same definition in all three packages; see virial and pressure for the differences. Direct calls to virial, scalar_virial, pressure and scalar_pressure approximate constraint contributions with a deterministic small-step constraint preview. Contributions from implicit solvent methods are ignored. Use a non-virial barostat such as the MonteCarloBarostat to control pressure in cases where some virial contributions are missing.

During supported constrained simulations, Molly can add constraint contributions to the total virial for steps where a barostat or virial/pressure logger requests it. For the initial simulation step, the same preview convention is used so that interactions and constraints both contribute to the logged virial/pressure. If a coordinate-scaling coupling method changes the box on a constrained step, virial and pressure loggers record the pre-coupling virial/pressure for that step, matching the state used by the coupling method. Other state loggers, such as BoxLogger, continue to record the current post-coupling state. The virial is compatible with virtual sites apart from OutOfPlaneSite and LocalCoordinatesSite. As described previously, custom general interactions should implement virial calculation if required.

Loggers

Loggers record properties of the simulation to allow monitoring and analysis. The available loggers are:

• GeneralObservableLogger
• TemperatureLogger
• CoordinatesLogger
• VelocitiesLogger
• TotalEnergyLogger
• KineticEnergyLogger
• PotentialEnergyLogger
• ForcesLogger
• VolumeLogger
• DensityLogger
• VirialLogger
• ScalarVirialLogger
• PressureLogger
• ScalarPressureLogger
• DisplacementsLogger
• TrajectoryWriter
• TimeCorrelationLogger
• AutoCorrelationLogger
• AverageObservableLogger
• ReplicaExchangeLogger
• MonteCarloLogger

Many of the loggers can be initialised with just the number of steps between recorded values, e.g. CoordinatesLogger(10). An optional first argument is the type of the recorded value; the above is equivalent to CoordinatesLogger(typeof(1.0u"nm"), 10) but if the simulation did not use units then CoordinatesLogger(Float64, 10) would be required. If the simulation is in 2D, giving dims=2 as a keyword argument is required for some loggers. A logger's history can be accessed with values(my_logger).

To define your own logger, first define the struct and a method for values to access the stored values:

struct MyLogger
    n_steps::Int
    history::Vector{Float64}
    # Any other properties
end

Base.values(logger::MyLogger) = logger.history

Then, define the logging function that is called every step by the simulator:

function Molly.log_property!(logger::MyLogger,
                             sys,
                             neighbors,
                             step_n,
                             buffers;
                             n_threads=Threads.nthreads(),
                             kwargs...)
    if step_n % logger.n_steps == 0
        # Record some property or carry out some action
    end
end

The use of n_steps is optional and is an example of how to record a property periodically throughout the simulation. If a custom logger needs the total virial for constrained simulations, define Molly.logger_virial_interval so the simulator can compute the virial on the same interval:

Molly.logger_virial_interval(logger::MyLogger) = logger.n_steps

If it specifically logs pressure, also define Molly.logger_pressure_interval:

Molly.logger_pressure_interval(logger::MyLogger) = logger.n_steps

Both functions return Inf by default, meaning no virial or pressure state is needed. To use your custom logger, add it to the named tuple of loggers given when creating the System:

loggers = (mylogger=MyLogger(10, []),) # Don't forget the trailing comma!

In addition to being run at the end of each step, loggers are run before the first step, i.e. at step 0. This means that a logger that records a value every step for a simulation with 100 steps will end up with 101 values. Running loggers before the first step can be disabled by giving run_loggers=:skipstart as a keyword argument to simulate!, which can be useful when splitting up simulations into multiple simulate! calls. For example, this runs the loggers 301 times:

simulate!(sys, simulator, 100) # Default run_loggers=true
simulate!(sys, simulator, 100; run_loggers=:skipstart)
simulate!(sys, simulator, 100; run_loggers=:skipstart)

Running loggers can be disabled entirely with run_loggers=false, which is the default for the energy minimizers such as SteepestDescentMinimizer.

Many times, a logger will just record an observation to an Array containing a record of past observations. For this purpose, you can use the GeneralObservableLogger without defining a custom logging function. Define your observation function as

function my_observable(sys::System, neighbors, step_n, buffers; n_threads::Integer, kwargs...)
    # Probe the system for some desired property
    return observation
end

Keyword arguments current_forces and current_potential_energy can also be used here to avoid recomputing values that are passed from the simulator:

function my_pe_observable(sys::System, neighbors, step_n, buffers; n_threads::Integer,
                          current_potential_energy=nothing, kwargs...)
    if isnothing(current_potential_energy)
        # Compute potential energy
        return potential_energy(sys, neighbors; n_threads=n_threads)
    else
        # Potential energy was passed from simulator, reuse
        return current_potential_energy
    end
end

These keyword arguments are also available in Molly.log_property!. Which values are passed depends on the simulator being used, for example SteepestDescentMinimizer passes current_potential_energy because it uses it for minimization. Note that loggers are called after Molly.apply_coupling!, so the coordinates may have changed since the potential energy or forces were computed.

A logger which records the property every n_steps can be constructed through

my_logger = GeneralObservableLogger(my_observable, T, n_steps)

where T = typeof(observation) is the type of the return value for my_observable. AverageObservableLogger is similar but records a running average rather than storing observations.

The TimeCorrelationLogger logger can be used to compute correlation functions of the form

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

where A and B are scalar or vector centered observables and the brackets are ensemble averages. This includes the computations of autocorrelation functions, which can be used to gather insight into the dynamical properties of the system, for instance using Green-Kubo formulas, or the statistical properties of a sampling method.

Let's look at a simple example, computing the velocity autocorrelation function for a simple system consisting of diatomic molecules defined by HarmonicBond potentials between pairs of atoms, and an additional SoftSphere potential between all pairs of atoms. Let's start by defining the system.

n_atoms = 400
atom_mass = 10.0u"g/mol"
atoms = [Atom(mass=atom_mass, σ=0.2u"nm", ϵ=0.2u"kJ * mol^-1") for i in 1:n_atoms]

# Initialization
boundary = CubicBoundary(6.0u"nm")
coords = place_diatomics(n_atoms ÷ 2, boundary, 0.2u"nm"; min_dist=0.2u"nm")

temp = 50.0u"K"
velocities = [random_velocity(atom_mass, temp) .* 0.01 for i in 1:n_atoms]

# Interaction potentials
pairwise_inters = (SoftSphere(cutoff=DistanceCutoff(0.6u"nm")),)

bonds = [HarmonicBond(k=10000u"kJ * mol^-1 * nm^-2", r0=0.2u"nm") for i in 1:(n_atoms ÷ 2)]
specific_inter_lists = (InteractionList2Atoms(
    collect(1:2:n_atoms),
    collect(2:2:n_atoms),
    bonds,
),)

# Define system
nf = DistanceNeighborFinder(n_atoms=n_atoms, dist_cutoff=0.6u"nm")

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    velocities=velocities,
    pairwise_inters=pairwise_inters,
    specific_inter_lists=specific_inter_lists,
    neighbor_finder=nf,
)

We leave the loggers empty until we thermalize the system using Langevin dynamics.

dt = 0.002u"ps"
simulator = LangevinSplitting(
    dt=dt,
    temperature=temp,
    friction=10.0u"g * mol^-1 * ps^-1",
    splitting="BAOAB",
)
simulate!(sys, simulator, 10_000)
@show temperature(sys)

temperature(sys) = 48.76795299825687 K

Next we define our correlation logger, add it to the system's loggers and run a longer simulation. Note that we need to redeclare the system when adding a logger.

# Velocity observable
# args and kwargs because more complex observables may require neighbors and parallelism
V(s::System, args...; kwargs...) = s.velocities
V_Type = eltype(sys.velocities)
logger = TimeCorrelationLogger(V, V, V_Type, V_Type, n_atoms, 1_000)

sys = System(
    atoms=atoms,
    coords=sys.coords,
    boundary=boundary,
    velocities=sys.velocities,
    pairwise_inters=pairwise_inters,
    specific_inter_lists=specific_inter_lists,
    neighbor_finder=nf,
    loggers=(velocity_autocorrelation=logger,)
)
simulate!(sys, simulator, 100_000)

Check the output:

show(sys.loggers)

(velocity_autocorrelation = AutoCorrelationLogger with n_correlation 1000 and 100001 samples collected for observable V,)

Note we also could have used the convenience function AutoCorrelationLogger to define our logger since the two observables we are correlating are the same.

using GLMakie
f = Figure()
ax = Axis(f[1, 1], xlabel="Time / ps", ylabel="Correlation")
lines!((1:1000) .* ustrip(dt), values(sys.loggers.velocity_autocorrelation))

As expected, the velocities are highly correlated at small time offsets and the correlation decays rapidly. The oscillatory behavior is due to the harmonic bond interactions.

Constraints

Molly implements SHAKE/RATTLE (see SHAKE_RATTLE) and LINCS (see LINCS) to perform constrained molecular dynamics. These methods are useful for maintaining bond lengths and angles during a simulation, often allowing the use of longer time steps and therefore more efficient use of computing resources.

SHAKE/RATTLE

The constraints satisfied by SHAKE are solely on the atomic coordinates:

\[\begin{aligned} \vec{r}_{ij} \cdot \vec{r}_{ij} = d^2_{ij} \end{aligned}\]

whereas RATTLE also constrains the velocities:

\[\begin{aligned} \vec{r}_{ij} \cdot \vec{v}_{ij} = 0 \end{aligned}\]

Here $\vec{r}_{ij}$ is the vector between atoms i and j in a constraint, $d_{ij}$ is the bond length to be maintained and $\vec{v}_{ij}$ is the difference in the velocity vectors for atoms i and j. SHAKE was originally derived for the Verlet integration scheme (Ryckaert et al. 1977) with RATTLE extending SHAKE to work for velocity Verlet where the velocities are also integrated (Andersen 1983).

LINCS

LINCS (LINear Constraint Solver) is a non-iterative constraint algorithm that uses matrix expansion to approximate the inverse of the constraint coupling matrix (Hess et al. 1997). It is typically faster than SHAKE/RATTLE for large systems. The implementation in Molly includes explicit velocity constraints. The key parameters controlling accuracy are n_rec, the order of the matrix expansion for coupling matrix inversion (default 4), and n_iter, the number of outer correction iterations for rotational lengthening (default 1). Higher values of either improve accuracy at the cost of performance.

Currently, constraints are supported by the following simulators:

• VelocityVerlet
• DPDVelocityVerlet
• Verlet
• StormerVerlet
• Langevin
• NoseHoover
• MTSIntegrator
• MTSLangevinIntegrator

Due to the nature of the velocity treatment in each integrator, the velocities stored in the system are expected to satisfy the constraints only for VelocityVerlet, DPDVelocityVerlet, MTSIntegrator and MTSLangevinIntegrator (in other cases the simulator behaviour is still correct).

The following simulators automatically use harmonic bonds in place of constraints, where the force constant can be adjusted by changing constraint_bond_constant:

• SteepestDescentMinimizer
• FIREMinimizer
• LBFGSMinimizer

Simulators incompatible with constraints will print a warning and continue without applying constraints when used with systems containing constraints.

Molly supports DistanceConstraint and AngleConstraint on CPU and GPU. Distance constraints fix the distance between two atoms. Angle constraints are defined for sets of three atoms (e.g. water) and restrict the angle and the two bond lengths.

Constraints can be added when setting up a system from a file as described above. To add constraints to a system manually, use something like the following:

dist_constraints = [
    DistanceConstraint(1, 2, 1.0u"nm"),
    DistanceConstraint(2, 3, 0.11u"nm"),
]
angle_constraints = [AngleConstraint(4, 5, 6, deg2rad(104.5), 0.1u"nm", 0.1u"nm")]

shake = SHAKE_RATTLE(
    n_atoms=6,
    dist_constraints=dist_constraints,
    angle_constraints=angle_constraints,
)
# SHAKE_RATTLE with 0 2-atom clusters, 1 3-atom clusters, 0 4-atom clusters and 1 angle clusters

constraints=(shake,) can then be given when setting up a System.

Alternatively, using LINCS:

masses = [1.0u"g/mol", 1.0u"g/mol", 1.0u"g/mol", 1.0u"g/mol", 1.0u"g/mol", 1.0u"g/mol"]

lincs = LINCS(
    masses=masses,
    dist_constraints=dist_constraints,
    angle_constraints=angle_constraints,
)
# LINCS with 2 distance and 1 angle constraints (n_rec=4, n_iter=1) (implicitly 5 distance constraints)

constraints=(lincs,) can then be given when setting up a System.

The DistanceConstraints, AngleConstraints and masses given to SHAKE_RATTLE and LINCS should always be on the CPU, even for a GPU system. The System constructor moves the constraint data to the device of the system, so no manual transfer is required.

See this example for more.

This diagram demonstrates the four allowed constraint types:

• Single bond between two atoms (e.g. a hydrogen molecule).
• 3 atoms, 2 DistanceConstraints, angle is free (e.g. the hydrogens on either carbon of ethelene).
• 3 atoms, 1 AngleConstraint, implemented as 3 distance constraints (e.g. a water molecule).
• 3 atoms around 1 central atom, 3 DistanceConstraints (e.g. an ammonia molecule).

These constraints provide enough flexibility to constrain all hydrogen atoms in organic molecules as well as water molecules.

For SHAKE/RATTLE, all velocity constraints and diatomic distance constraints are solved analytically while larger constraints are linearized and solved iteratively via matrix inverse. The direct matrix inverse does not scale well beyond clusters with 3 constraints and is not implemented. Other methods can be used to solve larger constraint clusters, these are not yet supported by Molly. LINCS uses a matrix expansion approach that scales better with system size and does not have the same cluster size limitations.

Virtual sites

Virtual sites are massless particles whose coordinates are defined by the coordinates of other atoms. One use case is to carry partial charge at a location separate from the atom centers, as seen in four-point water models like TIP4P. Molly allows virtual sites to be defined in the following ways:

• OneParticleSite: defined to have the same coordinates as another atom, can be useful in alchemical simulations when multiple versions of an atom are required.
• TwoParticleAverageSite: defined by the weighted average of the coordinates of two atoms.
• ThreeParticleAverageSite: defined by the weighted average of the coordinates of three atoms.
• OutOfPlaneSite: defined by the weighted average of the coordinates of three atoms and the cross product of their relative displacements.
• LocalCoordinatesSite: defined by a position in a local coordinate system given by three atoms, matching the site of the same name in OpenMM. This places a site at a fixed distance and orientation whatever the bond lengths are, as used for the lone pairs of a CHARMM force field.

Virtual sites should have an entry in the atom, coordinate and velocity arrays. They can be involved in any interaction type, with the forces being distributed back to the parent atoms automatically after all forces have been calculated. forces, accelerations and sys.velocities are zero for virtual site atoms since they are not integrated. They share all the non-bonded exclusions of, and are excluded from, the first of their parent atoms. The parent atoms must not be virtual sites themselves. They cannot participate in constraints. Virtual sites apart from OutOfPlaneSite and LocalCoordinatesSite are compatible with virial calculation, since the coordinates of the other sites are linear in the coordinates of their parent atoms.

A virtual site can be set up manually, for example for a molecule of TIP4P water:

using Molly

coords = [
    SVector(0.5558, 1.9020, 1.2139), # O
    SVector(0.4860, 1.9198, 1.2769), # H1
    SVector(0.5800, 1.9881, 1.1799), # H2
    SVector(1.0   , 1.0   , 1.0   ), # Virtual site, initial coordinates not used
] * u"nm"

# Virtual sites should have zero mass
atoms = [
    Atom(mass=15.99943u"g/mol", charge=0.0      , σ=0.316555u"nm", ϵ=0.749279u"kJ * mol^-1"),
    Atom(mass=1.007947u"g/mol", charge=0.525868 , σ=0.0u"nm"     , ϵ=0.0u"kJ * mol^-1"     ),
    Atom(mass=1.007947u"g/mol", charge=0.525868 , σ=0.0u"nm"     , ϵ=0.0u"kJ * mol^-1"     ),
    Atom(mass=0.0u"g/mol"     , charge=-1.051736, σ=0.0u"nm"     , ϵ=0.0u"kJ * mol^-1"     ),
]

boundary = CubicBoundary(3.0u"nm")

virtual_sites = [
    # Atom 4 is a virtual site defined by atoms 1/2/3
    ThreeParticleAverageSite(4, 1, 2, 3, 0.820314, 0.089843, 0.089843),
]

sys = System(
    atoms=atoms,
    coords=coords,
    boundary=boundary,
    virtual_sites=virtual_sites,
)

sys.virtual_site_flags # [false, false, false, true]

Virtual sites are set up appropriately when reading in a structure file. The virtual site atoms need to be present in the structure file to assign the correct residue template, which should have the virtual sites defined by <VirtualSite> entries.

Neighbor finders

Neighbor finders find close atoms periodically throughout the simulation, saving on computation time by allowing the force calculation between distant atoms to be omitted. The available neighbor finders are:

• CellListMapNeighborFinder
• GPUNeighborFinder
• GPUCellListNeighborFinder
• DistanceNeighborFinder
• TreeNeighborFinder

The recommended neighbor finder is CellListMapNeighborFinder on CPU, GPUNeighborFinder on NVIDIA GPUs and GPUCellListNeighborFinder on other GPUs, falling back to DistanceNeighborFinder there when the box is too small for a cell list.

The find_neighbors function returns an object that can be indexed into, and all pairs can be extracted with neighbor_pairs:

nbs = find_neighbors(sys)
for i in eachindex(nbs)
    println(nbs[i])
end

neighbor_pairs(nbs)

Pairwise interactions use the neighbor list by default, i.e. use_neighbors returns true, so a System with pairwise interactions needs a neighbor finder unless the interactions are given use_neighbors=false. A pairwise interaction with use_neighbors=false considers all pairs of atoms, without the excluded and special pairs of the neighbor finder.

The pairs of atoms that interact are given to a neighbor finder as the number of atoms along with lists of the excluded pairs, which do not interact through the pairwise interactions, for example bonded atoms, and the special pairs, which have scaled interactions, for example 1-4 atoms:

dist_cutoff = 1.0u"nm" # Interaction cutoff distance
dist_buffer = 0.2u"nm" # Buffer distance
neighbor_finder = DistanceNeighborFinder(
    n_atoms=n_atoms,
    excluded_pairs=[(1, 2), (2, 3), (1, 3)],
    special_pairs=[(1, 4)],
    dist_cutoff=(dist_cutoff + dist_buffer),
)

Every neighbor finder takes these arguments. The pairs can be given as any iterable of (i, j) pairs, or as a sparse matrix whose true entries are the pairs. They are stored as SparsePairMatrixs, so the memory used is proportional to the number of pairs and grows linearly with the number of atoms. For a system on the GPU using DistanceNeighborFinder, GPUNeighborFinder or GPUCellListNeighborFinder, give the array type of the system as array_type, for example array_type=CuArray, so that the pairs are stored on the GPU. Systems set up from a file store the pairs in this way.

The dist_cutoff of a neighbor finder is the distance used to search for neighbors, and is not the same as the interaction cutoff distance (see Cutoffs). Since the neighbor list is reused until it is rebuilt, dist_cutoff should be the interaction cutoff distance plus a buffer distance, as above. By default the list is rebuilt when the atoms have moved far enough that a pair within the interaction cutoff could be missing from it, see Updating the neighbor list, so the buffer distance only affects how often the list is rebuilt. When setting up a System from a file the buffer is added automatically and the two distances are given separately as the dist_cutoff and dist_buffer keyword arguments.

Alternatively, the pairs can be described by an eligible matrix, which is false for pairs excluded from the pairwise interactions, and a special matrix, which is true for special pairs:

# All pairs are eligible apart from each atom with itself
eligible = trues(n_atoms, n_atoms)
for i in 1:n_atoms
    eligible[i, i] = false
end
neighbor_finder = DistanceNeighborFinder(
    eligible=eligible,
    special=falses(n_atoms, n_atoms),
    dist_cutoff=(dist_cutoff + dist_buffer),
)

Dense matrices take memory proportional to the square of the number of atoms, which limits the size of system that can be simulated, especially on the GPU. They can be useful for small systems or patterns of pairs that are not sparse. A SparsePairMatrix can also be given as the eligible and special matrices.

GPUNeighborFinder follows a different CUDA-specific path based on the tiled GPU strategy of Eastman and Pande 2010. Instead of materializing a conventional neighbor list, it stores sparse excluded and special pairs and lets the CUDA pairwise kernels reorder atoms and cache a compact list of interacting 32x32 tiles internally, building the exclusion masks of the few tiles that need them from the sparse pairs. The memory it uses grows linearly with the number of atoms, and for large systems the interacting tiles are found by searching a tree of bounding boxes, so the time taken also grows close to linearly. Accordingly, find_neighbors returns nothing for GPUNeighborFinder, and the force and energy kernels decide when to refresh the tile list. With auto_update=false it should still be called at the same points as for the other neighbor finders, since it then marks when the tile list is rebuilt. When using it, set dist_cutoff to the interaction cutoff distance plus a buffer distance as above.

GPUCellListNeighborFinder is an O(N) cell list that does materialize a neighbor list, returning a GPUCellListNeighborList, which also gives the neighbors of each atom as a padded matrix via ragged_neighbors. It runs on any GPU backend and is the best option on GPUs other than NVIDIA ones, where the tiled kernels of GPUNeighborFinder are not available. Three-dimensional CubicBoundary and TriclinicBoundary systems are supported, as long as opposite box faces are at least three times dist_cutoff apart so that the grid has at least three cells along every box axis. Like GPUNeighborFinder it stores the exclusions and special pairs as SparsePairMatrixs on the device, converting dense eligible and special matrices at construction rather than keeping them, since a dense mask is n_atoms^2 bytes on the device. Most of the memory it uses for a large system is the per-atom neighbor matrix and the pair list. The matrix can be left out with ragged=false when only the pairs are needed, as for the pairwise interactions of a System, which roughly doubles the number of atoms that fit on a GPU; systems set up from a file do this.

Updating the neighbor list

The built-in neighbor finders have three ways of deciding when to rebuild the neighbor list:

• auto_update=true, the default: the list is rebuilt only when the atoms have moved far enough since it was built that a pair within the interaction cutoff could be missing from it, so no interacting pairs are missed whatever the buffer distance. n_steps is ignored.
• auto_update=false with n_steps=1: the list is rebuilt every step without tracking how far the atoms have moved, which is suitable when there is no buffer, as for ReaxFF by default.
• auto_update=false with n_steps greater than one: the list is rebuilt every n_steps steps. Interacting pairs are missed if an atom moves further than about half the buffer distance in that time, which nothing checks, so the buffer distance must be chosen carefully.

With auto_update=true, the neighbor finder stores the coordinates that the list was built from. Each time it is called, the two largest displacements d₁ and d₂ of the atoms since then are found, and the list is rebuilt if d₁ + d₂ is larger than the buffer distance dist_cutoff - r_c, where r_c is the largest cutoff of the interactions that use the list. A pair of atoms that was further than dist_cutoff apart when the list was built can not then be within r_c of each other. This is less strict than the common rule of rebuilding when any atom has moved more than half the buffer distance. When a barostat changes the box, the displacements are measured relative to the scaling of the box and the buffer distance is reduced by the amount the box shrinks, so scaling the coordinates does not by itself cause a rebuild. The interaction cutoff r_c is read from the cutoff or dist_cutoff property of the pairwise interactions for which use_neighbors is true, and from Molly.neighbor_cutoff for general interactions that read the neighbor list, such as ANIPotential and ReaxFF. If one of these interactions has no cutoff or its cutoff is not known, or if the buffer distance is zero, the list is rebuilt every time without tracking the displacements. Custom general interactions that read the pairs of the neighbor list should define Molly.neighbor_cutoff.

Since no interacting pairs are missed with auto_update=true whatever the buffer distance, it can be chosen for speed alone: a larger buffer means fewer rebuilds but more pairs in the list. Molly.neighbor_update_counts gives the number of times a neighbor finder has built a list and checked the displacements, which helps to choose it. GPUNeighborFinder makes the check in the force and energy kernels each time they run.

Analysis

Molly contains some tools for analysing the results of simulations. Functions that may be useful for analysis include:

• visualize
• rdf
• distances
• displacements
• rmsd
• radius_gyration
• hydrodynamic_radius
• bond_angle
• torsion_angle
• momentum
• dipole_moment

Julia is a language well-suited to implementing all kinds of analysis for molecular simulations.

Randomness

Many functions can take in a random number generator to allow reproducible randomness. For example:

julia> using Molly, Random

julia> random_velocity(10.0u"g/mol", 300.0u"K")
3-element SVector{3, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(nm, ps^-1), 𝐋 𝐓^-1, nothing}}} with indices SOneTo(3):
    0.5749142510427172 nm ps^-1
   -0.5373536210638031 nm ps^-1
 -0.044898786241729376 nm ps^-1

julia> random_velocity(10.0u"g/mol", 300.0u"K"; rng=Xoshiro(10))
3-element SVector{3, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(nm, ps^-1), 𝐋 𝐓^-1, nothing}}} with indices SOneTo(3):
    0.1703351361585439 nm ps^-1
 -0.027009534351612234 nm ps^-1
   0.23543457751795477 nm ps^-1

julia> random_velocity(10.0u"g/mol", 300.0u"K"; rng=Xoshiro(10))
3-element SVector{3, Quantity{Float64, 𝐋 𝐓^-1, Unitful.FreeUnits{(nm, ps^-1), 𝐋 𝐓^-1, nothing}}} with indices SOneTo(3):
    0.1703351361585439 nm ps^-1
 -0.027009534351612234 nm ps^-1
   0.23543457751795477 nm ps^-1

This may not apply across Julia versions, though you can use StableRNGs.jl. It also does not apply across different backends such as CPU and GPU.

Functions that generate randomness for many atoms at once may use the Philox4x32-10 algorithm from PhiloxRNG.jl internally, with the given rng only being used to seed Philox. Philox generates the same raw random numbers on all backends for the same rng state. Still, whole simulations are not expected to match bit for bit across backends due to floating point rounding differences.

Other functions, such as random_velocity and the Monte Carlo methods, draw from the given rng directly. Some functions require Random.default_rng() for thread safety, and will error if given a different rng.

Performance tips

Here is a checklist to ensure that you are getting the optimal performance from your simulations:

• On CPU, you should tune the n_threads argument to simulate!. If running on a single thread, it should be 1. Otherwise you should try various values, including larger than the number of threads available to Julia (which balances the load appropriately). Make sure to start Julia with as many threads as possible using -t. Generally, Float32 is not much faster than Float64 on CPU so Float64 is the default float type.
• On GPU, using Float32 will give vastly better performance and is the default float type. You can try changing the number of threads for each kernel as described in the GPU acceleration section, but the defaults are generally suitable for modern hardware. Multiple simulations can be run on different GPUs using device!. It is not currently possible to split one simulation onto multiple devices.
• If you run a simulation using CUDA GPUs with GPUNeighborFinder, the optimize_cuda_launch_config! function can be used to automatically test several launch parameters for the CUDA kernels and select the most performant ones. This is done automatically when setting a system up from a file.
• Run a short simulate! call once to ensure JIT compilation. You can run it on deepcopy(sys) if you don't want to affect sys, though beware of side effects like writing out trajectory files and consider using run_loggers=false.
• Make sure all arrays, such as coordinates and velocities, are concretely typed.
• In general, using units doesn't slow things down as described in the Units section, but you could try running without units.

Troubleshooting

There are many reasons that setting up and running a simulation can go wrong:

• Errors in parameterising a system with a force field are often due to discrepancies between the residues found in the structure and those found in the force field file. See the OpenMM FAQ for more.
• NaNs can arise in simulations, usually due to the system "blowing up" after encountering a massive force. Reasons for a massive force include errors in assigning force field parameters, a time step that is too large, incorrectly set boundaries and a lack of energy minimisation. You can try writing out the trajectory every step or reproducing exactly the same simulation (see Randomness) to find the issue. See the OpenMM FAQ for more. Note that the actual error thrown could appear after the point where the problem originated. For example, errors about constraint convergence could indicate earlier problems with the coordinates.
• GPU memory can be exceeded. In this case, make sure you are running with Float32 and the GPU does not have other processes running. If you still run out of memory, the system may be too large. Molly does not currently support splitting a simulation over multiple GPUs.
• Adding an interaction causes a huge slowdown. You may have introduced a type instability or be allocating memory in the inner force/potential energy calls. See the Julia performance tips for more. Make sure to avoid hidden type promotions, for example 1 / 2 is a Float64 which will promote Float32s it interacts with. You can use T(1 / 2) to avoid this, where T is a function type parameter.
• All array types must be on the same device. For example, if the coordinates are a CuArray then the atoms also need to be a CuArray. This is checked during setup.
• Minimum box side (...) is less than 2 * dist_cutoff indicates that the box is too small for the specified non-bonded cutoff distance. Each box side should be twice the cutoff distance, otherwise a particle would see multiple periodic versions of the same atom (though it wouldn't in Molly since we use the minimum image convention).
• Unitful errors like DimensionError: 1.0 and 1.0 nm are not dimensionally compatible usually indicate that a quantity with physical dimensions has not been given an appropriate unit. Whilst frustrating, fixing these errors often prevents silent unit conversion errors later on.
• Enzyme errors when taking gradients are common, see the Differentiable simulation with Molly section.

Do open an issue if you run into problems.