Molly examples

The best examples for learning how the package works are in the Molly documentation section. Here we give further examples, showing what you can do with the package. Each is a self-contained block of code. Made something cool yourself? Make a PR to add it to this page.

Making and breaking bonds

There is an example of mutable atom properties in the main docs, but what if you want to make and break bonds during the simulation? In this case you can use a GeneralInteraction to make, break and apply the bonds. The partners of the atom can be stored in the atom type. We make a Logger to record when the bonds are present, allowing us to visualize them with the connection_frames keyword argument to visualize (this can take a while to plot).

using Molly
using Makie
using LinearAlgebra

struct BondableAtom
    mass::Float64
    σ::Float64
    ϵ::Float64
    partners::Set{Int}
end

struct BondableInteraction <: GeneralInteraction
    nl_only::Bool
    prob_formation::Float64
    prob_break::Float64
    dist_formation::Float64
    b0::Float64
    kb::Float64
end

function Molly.force!(forces, inter::BondableInteraction, s, i, j)
    i == j && return
    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r2 = sum(abs2, dr)
    if j in s.atoms[i].partners
        # Apply the force of a harmonic bond
        c = inter.kb * (norm(dr) - inter.b0)
        fdr = c * normalize(dr)
        forces[i] += fdr
        forces[j] -= fdr
        # Break bonds randomly
        if rand() < inter.prob_break
            delete!(s.atoms[i].partners, j)
            delete!(s.atoms[j].partners, i)
        end
    # Make bonds between close atoms randomly
    elseif r2 < inter.b0 * inter.dist_formation && rand() < inter.prob_formation
        push!(s.atoms[i].partners, j)
        push!(s.atoms[j].partners, i)
    end
    return nothing
end

struct BondLogger <: Logger
    n_steps::Int
    bonds::Vector{BitVector}
end

function Molly.log_property!(logger::BondLogger, s, step_n)
    if step_n % logger.n_steps == 0
        bonds = BitVector()
        for i in 1:length(s.coords)
            for j in 1:(i - 1)
                push!(bonds, j in s.atoms[i].partners)
            end
        end
        push!(logger.bonds, bonds)
    end
end

temperature = 0.01
timestep = 0.02
box_size = 10.0
n_steps = 2_000
n_atoms = 200

atoms = [BondableAtom(1.0, 0.1, 0.02, Set([])) for i in 1:n_atoms]
coords = [box_size .* rand(SVector{2}) for i in 1:n_atoms]
velocities = [velocity(1.0, temperature, dims=2) for i in 1:n_atoms]
general_inters = (SoftSphere(true), BondableInteraction(true, 0.1, 0.1, 1.1, 0.1, 2.0))

s = Simulation(
    simulator=VelocityVerlet(),
    atoms=atoms,
    general_inters=general_inters,
    coords=coords,
    velocities=velocities,
    temperature=temperature,
    box_size=box_size,
    neighbour_finder=DistanceNeighbourFinder(trues(n_atoms, n_atoms), 10, 2.0),
    thermostat=AndersenThermostat(5.0),
    loggers=Dict("coords" => CoordinateLogger(20, dims=2),
                    "bonds" => BondLogger(20, [])),
    timestep=timestep,
    n_steps=n_steps
)

simulate!(s)

connections = Tuple{Int, Int}[]
for i in 1:length(s.coords)
    for j in 1:(i - 1)
        push!(connections, (i, j))
    end
end

visualize(s.loggers["coords"],
            box_size,
            "sim_mutbond.gif",
            connections=connections,
            connection_frames=s.loggers["bonds"].bonds)