Guide to Writing MD Programs#

The goal of NB-LIB’s is to enable researchers to programmatically define molecular simulations. Traditionally these have been performed using a collection of executables and a manual workflow followed by a “black-box” simulation engine. NB-LIB allows users to script a variety of novel simulation and analysis workflows at a more granular level.

Many possible use cases are facilitated by the flexibility that NB-LIB allows. These include customized update rules, defining custom forces, or orchestrating swarms of simulations. NB-LIB also allows for writing conventional MD simulations and analysis.

This document goes over the steps to write MD programs using the API in NB-LIB that exposes features that are a part of the GROMACS package.

Global Definitions#

NB-LIB programs are written in C++ so its headers for I/O or advanced tasks must be included. In addition, one must include the headers for various capabilities and abstractions NB-LIB exposes as well. This can be directly copied from here. Finally, we use the namespace nblib for the data structures defined in the library. The last line in the block allows one to skip this specifier each time a function or a data structure is used.

#include <cstdio>

#include "nblib/box.h"
#include "nblib/forcecalculator.h"
#include "nblib/integrator.h"
#include "nblib/molecules.h"
#include "nblib/nbkerneloptions.h"
#include "nblib/particletype.h"
#include "nblib/simulationstate.h"
#include "nblib/topology.h"

using namespace nblib;

Define Particle Data#

// Parameters from a GROMOS compatible force-field 2016H66

struct OWaterAtom
{
    ParticleName         name = "Ow";
    Mass                 mass = 15.999;
    C6                   c6   = 0.0026173456;
    C12                  c12  = 2.634129e-06;
};

struct HwAtom
{
    ParticleName         name = "Hw";
    Mass                 mass = 1.00784;
    C6                   c6   = 0.0;
    C12                  c12  = 0.0;
};

struct CMethAtom
{
    ParticleName         name = "Cm";
    Mass                 mass = 12.0107;
    C6                   c6   = 0.01317904;
    C12                  c12  = 34.363044e-06;
};

struct HcAtom
{
    ParticleName         name = "Hc";
    Mass                 mass = 1.00784;
    C6                   c6   = 8.464e-05;
    C12                  c12  = 15.129e-09;
};

There can be as many structs of this kind as there are particle types in the system. Organizing the data like this is not strictly necessary, but is shown for the purpose of clarity. As shown here, there can be multiple particles that correspond to a single element as atomic mass can vary by molecular context. For example, the carbon atom in a carboxyl group would have different parameters from one in the methyl group. We can obtain the parameter set from any standard force-field, or generate new parameters to study new compounds or force fields. This example comes from the 2016H66 Parameter Set.

Defining Coordinates, Velocities and Force Buffers#

std::vector<gmx::RVec> coordinates = {
    { 0.794, 1.439, 0.610 }, { 1.397, 0.673, 1.916 }, { 0.659, 1.080, 0.573 },
    { 1.105, 0.090, 3.431 }, { 1.741, 1.291, 3.432 }, { 1.936, 1.441, 5.873 },
    { 0.960, 2.246, 1.659 }, { 0.382, 3.023, 2.793 }, { 0.053, 4.857, 4.242 },
    { 2.655, 5.057, 2.211 }, { 4.114, 0.737, 0.614 }, { 5.977, 5.104, 5.217 },
};

std::vector<gmx::RVec> velocities = {
    { 0.0055, -0.1400, 0.2127 }, { 0.0930, -0.0160, -0.0086 }, { 0.1678, 0.2476, -0.0660 },
    { 0.1591, -0.0934, -0.0835 }, { -0.0317, 0.0573, 0.1453 }, { 0.0597, 0.0013, -0.0462 },
    { 0.0484, -0.0357, 0.0168 }, { 0.0530, 0.0295, -0.2694 }, { -0.0550, -0.0896, 0.0494 },
    { -0.0799, -0.2534, -0.0079 }, { 0.0436, -0.1557, 0.1849 }, { -0.0214, 0.0446, 0.0758},
};

std::vector<gmx::RVec> forces = {
    { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
    { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
    { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
    { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
};

We can initialize coordinates for our particles using std::vector of gmx::RVec which is a specific data type for holding 3D vector quantities. Doxygen page on RVec here.

Writing the MD Program#

As with as any basic C++ program, there needs to be a main() function.

Define ParticleTypes#

int main()
{
    // Bring the parameter structs to scope
    OwAtom      owAtom;
    HwAtom      hwAtom;
    CMethAtom   cmethAtom;
    HcAtom      hcAtom;

    // Create the particles
    ParticleType Ow(owAtom.name, owAtom.mass);
    ParticleType Hw(hwAtom.name, hwAtom.mass);
    ParticleType Cm(cmethAtom.name, cmethAtom.mass);
    ParticleType Hc(hcAtom.name, hcAtom.mass);

As before, the helper struct to define ParticleType data is not strictly needed, but is shown for clarity. The line ParticleType CMethAtom(ParticleName("Cm"), Mass(12.0107)); would be sufficient.

Define Non-Bonded Interactions#

ParticleTypeInteractions interactions(CombinationRule::Geometric);

// add non-bonded interactions for the particle types
interactions.add(owAtom.name, owAtom.c6, owAtom.c12);
interactions.add(hwAtom.name, hwAtom.c6, hwAtom.c12);
interactions.add(cmethAtom.name, cmethAtom.c6, cmethAtom.c12);
interactions.add(hcAtom.name, hcAtom.c6, hcAtom.c12);

For the Lennard-Jones interactions, we define a ParticleTypeInteractions object. Each particle of the ParticleType interacts with each other based on the C6 and C12 parameters. These parameters of the two different particles are averaged using Geometric or LorentzBerthelot CombinationRule. More details here. By default CombinationRule::Geometric is selected.

We add the interaction parameters of each of the particle types into the ParticleTypeInteractions object. The result is a table that has interactions specified for all ParticleType pairs. The following matrix describes the pair-wise C6 parameter created using CombinationRule::Geometric.

#

Ow

Hw

Cm

Hc

Ow

0.0026

0.0

0.42

4.7e-4

Hw

0.0

0.0

0.0

0.0

Cm

0.42

0.0

0.013

1.05e-3

Hc

4.7e-4

0.0

1.05e-3

8.5e-5

For a particular interaction pair, the user can also override the specified CombinationRule with custom parameters. The following overload would replace the parameters computed from a CombinationRule between Ow and Cm particle types.

interactions.add("Ow", "Cm", 0.42, 42e-6);

To facilitate modular, reusable code, it is possible to combine multiple ParticleTypeInteractions objects. Assuming otherInteractions is defined, this can be done with interactions.merge(otherInteractions)

Define Molecules#

Molecule water("Water");
Molecule methane("Methane");

water.addParticle(ParticleName("O"), Ow);
water.addParticle(ParticleName("H1"), Hw);
water.addParticle(ParticleName("H2"), Hw);

water.addExclusion("H1", "O");
water.addExclusion("H2", "O");

methane.addParticle(ParticleName("C"), Cm);
methane.addParticle(ParticleName("H1"), Hc);
methane.addParticle(ParticleName("H2"), Hc);
methane.addParticle(ParticleName("H3"), Hc);
methane.addParticle(ParticleName("H4"), Hc);

methane.addExclusion("H1", "C");
methane.addExclusion("H2", "C");
methane.addExclusion("H3", "C");
methane.addExclusion("H4", "C");

We begin declaring molecules with their constituent particles. A string identifier must uniquely identify a specific particle within the molecule. It is also possible to define partial charges on each particle for the computation of Coulomb interactions. water.addParticle(ParticleName("O"), Charge(-0.04), Ow);

Adding exclusions ensures that non-bonded interactions are only computed when necessary. For example, if two particles share a bond, the potential energy of the bond makes the non-bonded term negligible. Particle self-exclusions are enabled by default. We use the unique identifiers specified during addParticle() for this and the listed interactions later.

Define Listed Interactions#

Within a molecule, one can define interactions such as bonds, angles and dihedrals between the constituent particles. NB-LIB provides concrete implementations of several commonly used 2, 3 and 4 center interactions.

HarmonicBondType ohHarmonicBond(1, 1);
HarmonicBondType hcHarmonicBond(2, 1);

DefaultAngle hohAngle(Degrees(120), 1);
DefaultAngle hchAngle(Degrees(109.5), 1);

//add harmonic bonds for water
water.addInteraction("O", "H1", ohHarmonicBond);
water.addInteraction("O", "H2", ohHarmonicBond);

// add the angle for water
water.addInteraction("H1", "O", "H2", hohAngle);

// add harmonic bonds for methane
methane.addInteraction("H1", "C", hcHarmonicBond);
methane.addInteraction("H2", "C", hcHarmonicBond);
methane.addInteraction("H3", "C", hcHarmonicBond);
methane.addInteraction("H4", "C", hhcHarmonicBondc);

// add the angles for methane
methane.addInteraction("H1", "C", "H2", hchAngle);
methane.addInteraction("H1", "C", "H3", hchAngle);
methane.addInteraction("H1", "C", "H4", hchAngle);
methane.addInteraction("H2", "C", "H3", hchAngle);
methane.addInteraction("H2", "C", "H4", hchAngle);
methane.addInteraction("H3", "C", "H4", hchAngle);

Define Options for the Simulation and Non-Bonded Calculations#

// Define a box for the simulation
Box box(6.05449);

// Define options for the non-bonded kernels
NBKernelOptions options;

One can define the bounding box either with a single argument for a cube and 3 arguments to specify length, breadth and height separately.

NBKernelOptions contains a set of flags and configuration options for both hardware context and the relevant calculations for the simulation. The following table describes the possible options that can be set.

Flag or Config Option

Type

Implications

useGpu

Bool ean

Use GPU for non-bonded computations

numThreads

Inte ger

Number of CPU threads to use

nbnxmSimd

Enum

Kernel SIMD type (SimdAuto/SimdNo/Simd4XM/ Simd2XMM)

ljCombination Rule

Enum

Lennard-Jones combination rule (Geometric/LorentzBerthelot)

useHalfLJOptimizat ion

Bool ean

Enable i-cluster half-LJ optimization

pairlistCutoff

Real

Specify pairlist and interaction cut-off

computeVirialAndEn ergy

Bool ean

Enable energy computations

coulombType

Enum

Coulomb interaction function (Pme/Cutoff/ReactionField )

useTabulatedEwaldC orr

Bool ean

Use tabulated PME grid correction instead of analytical

numIterations

Inte ger

Specify number of iterations for each kernel

cyclesPerPair

Bool ean

Enable printing cycles/pair instead of pairs/cycle

timestep

Real

Specify the time step

Define Topology and Simulation State#

We build the system topology using the TopologyBuilder class. We add the Molecule objects that we defined previously along with the ParticleTypesInteractions using its public functions. We get the actual Topology object complete with all exclusions, interaction maps and listed interaction data constructed based on the defined entities using the buildTopology()function.

TopologyBuilder topologyBuilder;

// add molecules
topologyBuilder.addMolecule(water, 10);
topologyBuilder.addMolecule(methane, 10);

// add non-bonded interaction map
topologyBuilder.addParticleTypesInteractions(interactions);

Topology topology = topologyBuilder.buildTopology();

We now have all we need to fully describe our system using the SimulationState object. This is built using the topology, the box, and the particle coordinates and velocities. This object serves as a snapshot of the system that can be used for analysis or to start simulations from known states.

SimulationState simulationState(coordinates, velocities, forces, box, topology);

Writing the MD Loop#

Now that we have fully described our system and the problem, we need two entities to write an MD loop. The first is the ForceCalculator and the second is an Integrator. NB-LIB comes with a LeapFrog integrator but it is also possible for users to write custom integrators.

// The force calculator contains all the data needed to compute forces
ForceCalculator forceCalculator(simulationState, options);

// Integration requires masses, positions, and forces
LeapFrog integrator(simulationState);

// Allocate a force buffer
gmx::ArrayRef<gmx::RVec> userForces(topology.numParticles());

// MD Loop
int numSteps = 100;

for (i = 0; i < numSteps; i++)
{
  userForces = forceCalculator.compute();

  // The forces are not automatically updated in case the user wants to add their own
  std::copy(userForces.begin(), userForces.end(), begin(simulationState.forces()));

  // Integrate with a time step of 1 fs
  integrator.integrate(1.0);
}

return 0;
} // main