Molly API

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

A way to find near atoms to save on simulation time. Custom neighbor finders should sub-type this type.

AndersenThermostat(temperature, coupling_const)

Rescale random velocities according to the Andersen thermostat.

Atom(; <keyword arguments>)

An atom and its associated information. Properties unused in the simulation or in analysis can be left with their default values. The types used should be bits types if the GPU is going to be used.

Arguments

• index::Int: the index of the atom in the system.
• charge::C=0.0: the charge of the atom, used for electrostatic interactions.
• mass::M=0.0u"u": the mass of the atom.
• σ::S=0.0u"nm": the Lennard-Jones finite distance at which the inter-particle potential is zero.
• ϵ::E=0.0u"kJ * mol^-1": the Lennard-Jones depth of the potential well.
• solute::Bool=false: whether the atom is part of the solute.

AtomData(atom_type, atom_name, res_number, res_name)

Data associated with an atom. Storing this separately allows the atom types to be bits types and hence work on the GPU.

BerendsenThermostat(temperature, coupling_const)

The Berendsen thermostat. This thermostat should not be used in general as it can lead to simulation artifacts.

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

Find close atoms by distance, and store auxiliary arrays for in-place threading. x0 and unit_cell are optional initial coordinates and system unit cell that improve the first approximation of the cell list structure. The unit cell can be provided as a three-component vector of box sides on each direction, in which case the unit cell is considered OrthorhombicCell, or as a unit cell matrix, in which case the cell is considered a general TriclinicCell by the cell list algorithm.

Example

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

julia> box_size
3-element SVector{3, Quantity{Float64, 𝐋, Unitful.FreeUnits{(nm,), 𝐋, nothing}}} with indices SOneTo(3):
              5.676 nm
             5.6627 nm
             6.2963 nm

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

CoordinateLogger(n_steps; dims=3)

Log the coordinates throughout a simulation.

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

The Coulomb electrostatic interaction.

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

The Coulomb electrostatic interaction modified using the reaction field approximation.

DistanceCutoff(dist_cutoff)

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

DistanceNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff)

Find close atoms by distance.

DistanceVecNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff)

Find close atoms by distance in a GPU and differentiable safe manner.

Gravity(; G, nl_only)

The gravitational interaction.

HarmonicAngle(; th0, cth)

A harmonic bond angle between three atoms.

HarmonicBond(; b0, kb)

A harmonic bond between two atoms.

ImplicitSolventOBC(atoms, atoms_data, bonds)

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

An interaction between atoms that contributes to forces on the atoms.

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

A list of specific interactions between two atoms.

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

A list of specific interactions between three atoms.

InteractionList4Atoms(is, js, ks, ls, types, inters)
InteractionList4Atoms(inter_type)

A list of specific interactions between four atoms.

KineticEnergyLogger(n_steps)

Log the kinetic energy of the system throughout a simulation.

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

The Lennard-Jones 6-12 interaction. The potential is given by

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

and the force on each atom by

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

Mie(; m, n, cutoff, nl_only, lorentz_mixing, force_units, energy_units, skip_shortcut)

The Mie generalized interaction. When m equals 6 and n equals 12 this is equivalent to the Lennard-Jones interaction.

NeighborList()
NeighborList(n, list)

Structure to contain pre-allocated neighbor lists.

NeighborListVec(n, list)

Structure to contain neighbor lists for broadcasting.

NoCoupling()

Placeholder coupler that does nothing.

NoCutoff()

Placeholder cutoff that does not alter forces or potentials.

NoNeighborFinder()

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

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

An OpenMM atom type.

OpenMMForceField(ff_files...)
OpenMMForceField(T, ff_files...)
OpenMMForceField(atom_types, residue_types, bond_types, angle_types,
                    torsion_types, torsion_order, weight_14_coulomb,
                    weight_14_lj)

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

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

PeriodicTorsion(; periodicities, phases, ks, proper)

A periodic torsion angle between four atoms.

PeriodicTorsionType(periodicities, phases, ks, proper)

A periodic torsion type.

PotentialEnergyLogger(n_steps)

Log the potential energy of the system throughout a simulation.

RBTorsion(; f1, f2, f3, f4)

A Ryckaert-Bellemans torsion angle between four atoms.

RescaleThermostat(temperature)

The velocity rescaling thermostat. This thermostat should not be used in general as it can lead to simulation artifacts.

ShiftedForceCutoff(dist_cutoff)

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

ShiftedPotentialCutoff(dist_cutoff)

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

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

The soft-sphere potential.

SpecificForce2Atoms(f1, f2)

Forces on two atoms arising from an interaction.

SpecificForce3Atoms(f1, f2, f3)

Forces on three atoms arising from an interaction.

SpecificForce4Atoms(f1, f2, f3, f4)

Forces on four atoms arising from an interaction.

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

StormerVerlet(; <keyword arguments>)

The Störmer-Verlet integrator. In this case the velocities given to the simulator act as the previous step coordinates for the first step. Does not currently work with units or thermostats.

Arguments

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

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

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

System(; <keyword arguments>)

A physical system to be simulated. Properties unused in the simulation or in analysis can be left with their default values. atoms, atoms_data, coords and velocities should have the same length. This is a sub-type of AbstractSystem from AtomsBase.jl and implements the interface described there.

Arguments

• atoms::A: the atoms, or atom equivalents, in the system. Can be of any type but should be a bits type if the GPU is used.
• atoms_data::AD: other data associated with the atoms, allowing the atoms to be bits types and hence work on the GPU.
• pairwise_inters::PI=(): the pairwise interactions in the system, i.e. interactions between all or most atom pairs such as electrostatics. Typically a Tuple.
• specific_inter_lists::SI=(): the specific interactions in the system, i.e. interactions between specific atoms such as bonds or angles. Typically a Tuple.
• general_inters::GI=(): the general interactions in the system, i.e. interactions involving all atoms such as implicit solvent. Typically a Tuple.
• coords::C: the coordinates of the atoms in the system. Typically a vector of SVectors of 2 or 3 dimensions.
• velocities::V=zero(coords): the velocities of the atoms in the system.
• box_size::B: the size of the box in which the simulation takes place. Typically a SVector of 2 or 3 dimensions.
• neighbor_finder::NF=NoNeighborFinder(): the neighbor finder used to find close atoms and save on computation.
• loggers::L=Dict(): the loggers that record properties of interest during a simulation.
• force_units::F=u"kJ * mol^-1 * nm^-1": the units of force of the system.
• energy_units::E=u"kJ * mol^-1": the units of energy of the system.
• gpu_diff_safe::Bool: whether to use the code path suitable for the GPU and taking gradients. Defaults to isa(coords, CuArray).

System(coordinate_file, force_field; <keyword arguments>)

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

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

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

Arguments

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

TemperatureLogger(n_steps)
TemperatureLogger(T, n_steps)

Log the temperature throughout a simulation.

TotalEnergyLogger(n_steps)

Log the total energy of the system throughout a simulation.

TreeNeighborFinder(; nb_matrix, matrix_14, n_steps, dist_cutoff)

Find close atoms by distance using a tree search.

VelocityLogger(n_steps; dims=3)

Log the velocities throughout a simulation.

VelocityVerlet(; <keyword arguments>)

The velocity Verlet integrator.

Arguments

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

velocity(mass, temperature; dims=3)

Generate a random velocity from the Maxwell-Boltzmann distribution.

accelerations(system, neighbors=nothing; parallel=true)

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

apply_coupling!(system, simulator, coupling)

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

bond_angle(coord_i, coord_j, coord_k, box_size)
bond_angle(vec_ba, vec_bc)

Calculate the bond or pseudo-bond angle in radians between three coordinates or two vectors. The angle between B→A and B→C is returned in the range 0 to π.

born_radii_and_grad(inter, coords, box_size)

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

charge(atom)

The partial charge of an atom.

displacements(coords, box_size)

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

distances(coords, box_size)

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

extract_parameters(system, force_field)

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

find_neighbors(s::System,
                nf::CellListMapNeighborFinder,
                current_neighbors=nothing,
                step_n::Integer=0;
                parallel::Bool=true)

Find neighbors using CellListMap, without in-place updating. Should be called only the first time the cell lists are built. Modifies the mutable nf structure.

find_neighbors(system, neighbor_finder, current_neighbors=nothing,
                step_n=0; parallel=true)
find_neighbors(system; parallel=true)

Obtain a list of close atoms in a system. Custom neighbor finders should implement this function.

float_type(sys)

The float type a System uses.

force(inter::PairwiseInteraction, vec_ij, coord_i, coord_j,
      atom_i, atom_j, box_size)
force(inter::SpecificInteraction, coord_i, coord_j,
      box_size)
force(inter::SpecificInteraction, coord_i, coord_j,
      coord_k, box_size)
force(inter::SpecificInteraction, coord_i, coord_j,
      coord_k, coord_l, box_size)

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

forces(system, neighbors=nothing; parallel=true)
forces(system, coords, atoms, neighbors=nothing, neighbors_all=nothing)

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

forces(inter, system, neighbors=nothing)

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

inject_gradients(sys, params_dic)

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

is_gpu_diff_safe(sys)

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

kinetic_energy(s)

Compute the kinetic energy of the system.

log_property!(logger, system, neighbors=nothing, step_n=0)

Log a property of the system thoughout a simulation. Custom loggers should implement this function.

mass(atom)

The mass of an atom.

maxwell_boltzmann(mass, temperature)

Generate a random speed along one dimension from the Maxwell-Boltzmann distribution.

place_atoms(n_atoms, box_size, min_dist)

Obtain n_atoms 3D coordinates in a box with sides box_size where no two points are closer than min_dist, accounting for periodic boundary conditions.

place_diatomics(n_molecules, box_size, min_dist, bond_length)

Obtain 3D coordinates for n_molecules diatomics in a box with sides box_size where no two points are closer than min_dist and the bond length is bond_length, accounting for periodic boundary conditions.

potential_energy(s, neighbors=nothing)

Compute the potential energy of the system. If the interactions use neighbor lists, the neighbors should be computed first and passed to the function.

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

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

random_velocities!(sys, temp)

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

random_velocities(sys, temp)

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

rdf(coords, box_size; npoints=200)

Get the radial distribution function of a set of coordinates. This describes how density varies as a function of distance from each atom. Returns a list of distance bin centres and a list of the corresponding densities.

run_loggers!(system, neighbors=nothing, step_n=0)

Run the loggers associated with the system.

simulate!(system, simulator, n_steps; parallel=true)
simulate!(system, simulator; parallel=true)

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

temperature(system)

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

torsion_angle(coord_i, coord_j, coord_k, coord_l, box_size)
torsion_angle(vec_ab, vec_bc, vec_cd)

Calculate the torsion angle in radians defined by four coordinates or three vectors. The angle between the planes defined by atoms (i, j, k) and (j, k, l) is returned in the range -π to π.

total_energy(s, neighbors=nothing)

Compute the total energy of the system. If the interactions use neighbor lists, the neighbors should be computed first and passed to the function.

ustrip_vec(x)

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

vector(c1, c2, box_size)

Displacement between two coordinate values from c1 to c2, accounting for the bounding box. The minimum image convention is used, so the displacement is to the closest version of the coordinates accounting for the periodic boundaries.

vector1D(c1, c2, side_length)

Displacement between two 1D coordinate values from c1 to c2, accounting for the bounding box. The minimum image convention is used, so the displacement is to the closest version of the coordinate accounting for the periodic boundaries.

velocity_autocorr(vl, first_ind, last_ind)

Calculates the autocorrelation function of velocity from the velocity logger. This helps characterize the similarity between velocities observed at different time instances.

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

Visualize a simulation as an animation. This function is only available when GLMakie is imported. GLMakie v0.5 or later should be used. It can take a while to run, depending on the length and size of the simulation.

Arguments

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

wrap_coords(c, side_length)

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

wrap_coords_vec(c, box_size)

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