JSoC '22 : My contributions to Molly.jl

During the summer of 2022 I had the opportunity to contribute to Molly.jl as part of the Julia Season of Contribution 2022.

As a college student interested in molecular simulation, I was exposed to Julia around the second year of my undergraduate studies. Going through various Julia forums and actively following the discourse chemistry and materials science section allowed me to find a thread related to a "Chemistry Birds of a Feather" session that would be hosted in JuliaCon 2021. I signed up for it immediately. After listening to the ideas proposed in the BoF, I felt as though I could help contribute to the growing chemistry community as well. Joining the JuliaMolSim Slack channel (now a Zulip channel), allowed me to understand what the community demanded and how I could contribute. Around November of 2021, I reached out to Joe Greener, the developer of Molly.jl, a julia-based molecular dynamics (MD) package, since I had some experience working with MD on LAMMPS. From November to March, I was able to make a few contributions and hence decided to sign up for JSoC '22 so that I could spend the summer contributing a bit more.

What is Molecular Dynamics?

Molecular Dynamics (MD) is a simulation technique that allows us to observe the dynamic evolution of a molecular system. It works by numerically integrating equations of motion (such as Newton's law) and plots out a trajectory of motion for atoms and molecules. It is able to model a variety of physical systems from gases to proteins, macromolecules and solids. It's a particularly handy simulation tool for researchers working with complex fluids, soft matter and biological macromolecules.

What is Molly.jl?

Molly.jl is a pure Julia package for molecular dynamics and other atomistic and molecular simulations. It has a focus on being differentiable and fast. It has several features but misses out on some of the capabilities present in popular MD software such as LAMMPS and GROMACS. One of these features is bond and angle constraint algorithms. This is an area I felt I could contribute to.

My contributions

I started by contributing to some other features to Molly.jl initially such as implementing an analysis tool (velocity autocorrelation) and adding a bonded interaction to the list of potentials.

The initial challenge was adding a "constraint" variable to the "System" struct of Molly.jl without it affecting the rest of the ecosystem. This was a bit challenging for me since I had a very limited prior experience contributing to open source software.

After adding a constraint variable, I started my work on adding bond constraints with the standard "SHAKE" scheme, following implementations from other open-source software as well as the original research texts relating to the algorithm.

Apart from constraints, I was also interested in working on improving the visualization aspects of julia chemistry packages, but wasn't able to get very far working on them.

An example

This is a simple example where the lengths of bonds are constrained by a SHAKE constraint algorithm.

Start by defining the atom properties:

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

    coords = place_atoms(n_atoms ÷ 2, boundary, min_dist=0.3u"nm")

    for i in 1:length(coords)
        push!(coords, coords[i].+ [0.15, 0.0, 0.0]u"nm")
    end

    temp = 100.0u"K"

    velocities = [velocity(atom_mass, temp) for i in 1:n_atoms]
    #=
    bonds = InteractionList2Atoms(
        collect(1:(n_atoms÷2)),
        collect((1+n_atoms÷2):n_atoms),
        repeat([""], n_atoms÷2),
        [],
       )
    =#

    nb_matrix = trues(n_atoms, n_atoms)

    for i in (1:n_atoms÷2)
        nb_matrix[i, i + (n_atoms ÷ 2)] = false
        nb_matrix[i + (n_atoms ÷ 2), i] = false
    end

    neighbor_finder = DistanceNeighborFinder(nb_matrix=nb_matrix, n_steps=10, dist_cutoff=1.5u"nm")

Now for the constraints, we can define bond lengths and

    bond_lengths = [0.1u"nm" for i in 1:(n_atoms÷2)]

    sh = SHAKE(bond_lengths, collect(1:(n_atoms÷2)), collect((1+n_atoms÷2):n_atoms))

    constraint_list = (sh,)

Now we can setup the system and the simulator:

    sys = System(atoms=atoms,
                 pairwise_inters=(LennardJones(nl_only=true),),

                            constraints=constraint_list,
                            coords=coords,
                            velocities=velocities,
                            boundary=boundary,
                            neighbor_finder=neighbor_finder,
                            loggers=Dict(
                                         "temp" => TemperatureLogger(10),
                                         "coords" => CoordinateLogger(10),
                                        )
                           )


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

    simulate!(sys, simulator, 1_000)

For visualizing the output, we can use GLMakie:

    using GLMakie

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