docs/nblib/guide-to-writing-MD-programs.rst

   1 Guide to Writing MD Programs
   2 ============================
   3
   4 The goal of NB-LIB’s is to enable researchers to programmatically define molecular simulations.
   5 Traditionally these have been performed using a collection of executables and a manual workflow followed by a “black-box” simulation engine.
   6 NB-LIB allows users to script a variety of novel simulation and analysis workflows at a more granular level.
   7
   8 Many possible use cases are facilitated by the flexibility that NB-LIB allows.
   9 These include customized update rules, defining custom forces, or orchestrating swarms of simulations.
  10 NB-LIB also allows for writing conventional MD simulations and analysis.
  11
  12 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.
  13
  14 Global Definitions
  15 ------------------
  16
  17 NB-LIB programs are written in C++ so its headers for I/O or advanced tasks must be included.
  18 In addition, one must include the headers for various capabilities and abstractions NB-LIB exposes as well.
  19 This can be directly copied from here.
  20 Finally, we use the namespace ``nblib`` for the data structures defined in the library.
  21 The last line in the block allows one to skip this specifier each time a function or a data structure is used.
  22
  23 .. code:: cpp
  24
  25    #include <cstdio>
  26
  27    #include "nblib/box.h"
  28    #include "nblib/forcecalculator.h"
  29    #include "nblib/integrator.h"
  30    #include "nblib/molecules.h"
  31    #include "nblib/nbkerneloptions.h"
  32    #include "nblib/particletype.h"
  33    #include "nblib/simulationstate.h"
  34    #include "nblib/topology.h"
  35
  36    using namespace nblib;
  37
  38 Define Particle Data
  39 --------------------
  40
  41 .. code:: cpp
  42
  43    // Parameters from a GROMOS compatible force-field 2016H66
  44
  45    struct OWaterAtom
  46    {
  47        ParticleName         name = "Ow";
  48        Mass                 mass = 15.999;
  49        C6                   c6   = 0.0026173456;
  50        C12                  c12  = 2.634129e-06;
  51    };
  52
  53    struct HwAtom
  54    {
  55        ParticleName         name = "Hw";
  56        Mass                 mass = 1.00784;
  57        C6                   c6   = 0.0;
  58        C12                  c12  = 0.0;
  59    };
  60
  61    struct CMethAtom
  62    {
  63        ParticleName         name = "Cm";
  64        Mass                 mass = 12.0107;
  65        C6                   c6   = 0.01317904;
  66        C12                  c12  = 34.363044e-06;
  67    };
  68
  69    struct HcAtom
  70    {
  71        ParticleName         name = "Hc";
  72        Mass                 mass = 1.00784;
  73        C6                   c6   = 8.464e-05;
  74        C12                  c12  = 15.129e-09;
  75    };
  76
  77 There can be as many structs of this kind as there are particle types in the system.
  78 Organizing the data like this is not strictly necessary, but is shown for the purpose of clarity.
  79 As shown here, there can be multiple particles that correspond to a single element as atomic mass can vary by molecular context.
  80 For example, the carbon atom in a carboxyl group would have different parameters from one in the methyl group.
  81 We can obtain the parameter set from any standard force-field, or generate new parameters to study new compounds or force fields.
  82 This example comes from the `2016H66 Parameter Set <https://pubs.acs.org/doi/10.1021/acs.jctc.6b00187>`__.
  83
  84 Defining Coordinates, Velocities and Force Buffers
  85 --------------------------------------------------
  86
  87 .. code:: cpp
  88
  89    std::vector<gmx::RVec> coordinates = {
  90        { 0.794, 1.439, 0.610 }, { 1.397, 0.673, 1.916 }, { 0.659, 1.080, 0.573 },
  91        { 1.105, 0.090, 3.431 }, { 1.741, 1.291, 3.432 }, { 1.936, 1.441, 5.873 },
  92        { 0.960, 2.246, 1.659 }, { 0.382, 3.023, 2.793 }, { 0.053, 4.857, 4.242 },
  93        { 2.655, 5.057, 2.211 }, { 4.114, 0.737, 0.614 }, { 5.977, 5.104, 5.217 },
  94    };
  95
  96    std::vector<gmx::RVec> velocities = {
  97        { 0.0055, -0.1400, 0.2127 }, { 0.0930, -0.0160, -0.0086 }, { 0.1678, 0.2476, -0.0660 },
  98        { 0.1591, -0.0934, -0.0835 }, { -0.0317, 0.0573, 0.1453 }, { 0.0597, 0.0013, -0.0462 },
  99        { 0.0484, -0.0357, 0.0168 }, { 0.0530, 0.0295, -0.2694 }, { -0.0550, -0.0896, 0.0494 },
 100        { -0.0799, -0.2534, -0.0079 }, { 0.0436, -0.1557, 0.1849 }, { -0.0214, 0.0446, 0.0758},
 101    };
 102
 103    std::vector<gmx::RVec> forces = {
 104        { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
 105        { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
 106        { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
 107        { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 }, { 0.0000, 0.0000, 0.0000 },
 108    };
 109
 110 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`__.
 111
 112   __ doxygen-ref-rvec_
 113
 114 Writing the MD Program
 115 ----------------------
 116
 117 As with as any basic C++ program, there needs to be a ``main()`` function.
 118
 119
 120 Define ParticleTypes
 121 ~~~~~~~~~~~~~~~~~~~~
 122
 123 .. code:: cpp
 124
 125    int main()
 126    {
 127        // Bring the parameter structs to scope
 128        OwAtom      owAtom;
 129        HwAtom      hwAtom;
 130        CMethAtom   cmethAtom;
 131        HcAtom      hcAtom;
 132
 133        // Create the particles
 134        ParticleType Ow(owAtom.name, owAtom.mass);
 135        ParticleType Hw(hwAtom.name, hwAtom.mass);
 136        ParticleType Cm(cmethAtom.name, cmethAtom.mass);
 137        ParticleType Hc(hcAtom.name, hcAtom.mass);
 138
 139 As before, the helper struct to define ``ParticleType`` data is not strictly needed, but is shown for clarity.
 140 The line ``ParticleType CMethAtom(ParticleName("Cm"), Mass(12.0107));`` would be sufficient.
 141
 142 Define Non-Bonded Interactions
 143 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 144
 145 .. code:: cpp
 146
 147    ParticleTypeInteractions interactions(CombinationRule::Geometric);
 148
 149    // add non-bonded interactions for the particle types
 150    interactions.add(owAtom.name, owAtom.c6, owAtom.c12);
 151    interactions.add(hwAtom.name, hwAtom.c6, hwAtom.c12);
 152    interactions.add(cmethAtom.name, cmethAtom.c6, cmethAtom.c12);
 153    interactions.add(hcAtom.name, hcAtom.c6, hcAtom.c12);
 154
 155 For the Lennard-Jones interactions, we define a ``ParticleTypeInteractions`` object.
 156 Each particle of the ``ParticleType`` interacts with each other based on the ``C6`` and ``C12`` parameters.
 157 These parameters of the two different particles are averaged using ``Geometric`` or ``LorentzBerthelot`` ``CombinationRule``.
 158 More details `here <http://manual.gromacs.org/documentation/2019/reference-manual/functions/nonbonded-interactions.html#the-lennard-jones-interaction>`__.
 159 By default ``CombinationRule::Geometric`` is selected.
 160
 161 We add the interaction parameters of each of the particle types into the ``ParticleTypeInteractions`` object.
 162 The result is a table that has interactions specified for all ``ParticleType`` pairs.
 163 The following matrix describes the pair-wise C6 parameter created using ``CombinationRule::Geometric``.
 164
 165 == ====== === ======= =======
 166 #  Ow     Hw  Cm      Hc
 167 == ====== === ======= =======
 168 Ow 0.0026 0.0 0.42    4.7e-4
 169 Hw 0.0    0.0 0.0     0.0
 170 Cm 0.42   0.0 0.013   1.05e-3
 171 Hc 4.7e-4 0.0 1.05e-3 8.5e-5
 172 == ====== === ======= =======
 173
 174 For a particular interaction pair, the user can also override the specified ``CombinationRule`` with custom parameters.
 175 The following overload would replace the parameters computed from a ``CombinationRule``  between ``Ow`` and ``Cm`` particle types.
 176
 177 .. code:: cpp
 178
 179    interactions.add("Ow", "Cm", 0.42, 42e-6);
 180
 181 To facilitate modular, reusable code, it is possible to combine multiple ``ParticleTypeInteractions`` objects.
 182 Assuming ``otherInteractions`` is defined, this can be done with ``interactions.merge(otherInteractions)``
 183
 184 Define Molecules
 185 ~~~~~~~~~~~~~~~~
 186
 187 .. code:: cpp
 188
 189    Molecule water("Water");
 190    Molecule methane("Methane");
 191
 192    water.addParticle(ParticleName("O"), Ow);
 193    water.addParticle(ParticleName("H1"), Hw);
 194    water.addParticle(ParticleName("H2"), Hw);
 195
 196    water.addExclusion("H1", "O");
 197    water.addExclusion("H2", "O");
 198
 199    methane.addParticle(ParticleName("C"), Cm);
 200    methane.addParticle(ParticleName("H1"), Hc);
 201    methane.addParticle(ParticleName("H2"), Hc);
 202    methane.addParticle(ParticleName("H3"), Hc);
 203    methane.addParticle(ParticleName("H4"), Hc);
 204
 205    methane.addExclusion("H1", "C");
 206    methane.addExclusion("H2", "C");
 207    methane.addExclusion("H3", "C");
 208    methane.addExclusion("H4", "C");
 209
 210 We begin declaring molecules with their constituent particles.
 211 A string identifier must uniquely identify a specific particle within the molecule.
 212 It is also possible to define partial charges on each particle for the computation of Coulomb interactions.
 213 ``water.addParticle(ParticleName("O"), Charge(-0.04), Ow);``
 214
 215 Adding exclusions ensures that non-bonded interactions are only computed when necessary.
 216 For example, if two  particles share a bond, the potential energy of the bond makes the non-bonded term negligible.
 217 Particle self-exclusions are enabled by default.
 218 We use the unique identifiers specified during ``addParticle()`` for this and the listed interactions later.
 219
 220 Define Listed Interactions
 221 ~~~~~~~~~~~~~~~~~~~~~~~~~~
 222
 223 Within a molecule, one can define interactions such as bonds, angles and dihedrals between the constituent particles.
 224 NB-LIB provides concrete implementations of several commonly used 2, 3 and 4 center interactions.
 225
 226 .. code:: cpp
 227
 228    HarmonicBondType ohHarmonicBond(1, 1);
 229    HarmonicBondType hcHarmonicBond(2, 1);
 230
 231    DefaultAngle hohAngle(Degrees(120), 1);
 232    DefaultAngle hchAngle(Degrees(109.5), 1);
 233
 234    //add harmonic bonds for water
 235    water.addInteraction("O", "H1", ohHarmonicBond);
 236    water.addInteraction("O", "H2", ohHarmonicBond);
 237
 238    // add the angle for water
 239    water.addInteraction("H1", "O", "H2", hohAngle);
 240
 241    // add harmonic bonds for methane
 242    methane.addInteraction("H1", "C", hcHarmonicBond);
 243    methane.addInteraction("H2", "C", hcHarmonicBond);
 244    methane.addInteraction("H3", "C", hcHarmonicBond);
 245    methane.addInteraction("H4", "C", hhcHarmonicBondc);
 246
 247    // add the angles for methane
 248    methane.addInteraction("H1", "C", "H2", hchAngle);
 249    methane.addInteraction("H1", "C", "H3", hchAngle);
 250    methane.addInteraction("H1", "C", "H4", hchAngle);
 251    methane.addInteraction("H2", "C", "H3", hchAngle);
 252    methane.addInteraction("H2", "C", "H4", hchAngle);
 253    methane.addInteraction("H3", "C", "H4", hchAngle);
 254
 255 Define Options for the Simulation and Non-Bonded Calculations
 256 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 257
 258 .. code:: cpp
 259
 260    // Define a box for the simulation
 261    Box box(6.05449);
 262
 263    // Define options for the non-bonded kernels
 264    NBKernelOptions options;
 265
 266 One can define the bounding box either with a single argument for a cube and 3 arguments to specify length, breadth and height separately.
 267
 268 ``NBKernelOptions`` contains a set of flags and configuration options for both hardware context and the relevant calculations for the simulation.
 269 The following table describes the possible options that can be set.
 270
 271 +----------------------+------+---------------------------------------+
 272 | Flag or Config       | Type | Implications                          |
 273 | Option               |      |                                       |
 274 +======================+======+=======================================+
 275 | ``useGpu``           | Bool | Use GPU for non-bonded computations   |
 276 |                      | ean  |                                       |
 277 +----------------------+------+---------------------------------------+
 278 | ``numThreads``       | Inte | Number of CPU threads to use          |
 279 |                      | ger  |                                       |
 280 +----------------------+------+---------------------------------------+
 281 | ``nbnxmSimd``        | Enum | Kernel SIMD type                      |
 282 |                      |      | (``SimdAuto``/``SimdNo``/``Simd4XM``/ |
 283 |                      |      | ``Simd2XMM``)                         |
 284 +----------------------+------+---------------------------------------+
 285 | ``ljCombination      | Enum | Lennard-Jones combination rule        |
 286 | Rule``               |      | (``Geometric``/``LorentzBerthelot``)  |
 287 +----------------------+------+---------------------------------------+
 288 | ``useHalfLJOptimizat | Bool | Enable i-cluster half-LJ optimization |
 289 | ion``                | ean  |                                       |
 290 +----------------------+------+---------------------------------------+
 291 | ``pairlistCutoff``   | Real | Specify pairlist and interaction      |
 292 |                      |      | cut-off                               |
 293 +----------------------+------+---------------------------------------+
 294 | ``computeVirialAndEn | Bool | Enable energy computations            |
 295 | ergy``               | ean  |                                       |
 296 +----------------------+------+---------------------------------------+
 297 | ``coulombType``      | Enum | Coulomb interaction function          |
 298 |                      |      | (``Pme``/``Cutoff``/``ReactionField`` |
 299 |                      |      | )                                     |
 300 +----------------------+------+---------------------------------------+
 301 | ``useTabulatedEwaldC | Bool | Use tabulated PME grid correction     |
 302 | orr``                | ean  | instead of analytical                 |
 303 +----------------------+------+---------------------------------------+
 304 | ``numIterations``    | Inte | Specify number of iterations for each |
 305 |                      | ger  | kernel                                |
 306 +----------------------+------+---------------------------------------+
 307 | ``cyclesPerPair``    | Bool | Enable printing cycles/pair instead   |
 308 |                      | ean  | of pairs/cycle                        |
 309 +----------------------+------+---------------------------------------+
 310 | ``timestep``         | Real | Specify the time step                 |
 311 +----------------------+------+---------------------------------------+
 312
 313 Define Topology and Simulation State
 314 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 315
 316 We build the system topology using the ``TopologyBuilder`` class.
 317 We add the ``Molecule`` objects that we defined previously along with the ``ParticleTypesInteractions`` using its public functions.
 318 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.
 319
 320 .. code:: cpp
 321
 322    TopologyBuilder topologyBuilder;
 323
 324    // add molecules
 325    topologyBuilder.addMolecule(water, 10);
 326    topologyBuilder.addMolecule(methane, 10);
 327
 328    // add non-bonded interaction map
 329    topologyBuilder.addParticleTypesInteractions(interactions);
 330
 331    Topology topology = topologyBuilder.buildTopology();
 332
 333 We now have all we need to fully describe our system using the ``SimulationState`` object.
 334 This is built using the topology, the box, and the particle coordinates and velocities.
 335 This object serves as a snapshot of the system that can be used for analysis or to start simulations from known states.
 336
 337 .. code:: cpp
 338
 339    SimulationState simulationState(coordinates, velocities, forces, box, topology);
 340
 341
 342
 343 Writing the MD Loop
 344 ~~~~~~~~~~~~~~~~~~~
 345
 346 Now that we have fully described our system and the problem, we need two entities to write an MD loop.
 347 The first is the ``ForceCalculator`` and the second is an Integrator.
 348 NB-LIB comes with a ``LeapFrog`` integrator but it is also possible for users to write custom integrators.
 349
 350 .. code:: cpp
 351
 352    // The force calculator contains all the data needed to compute forces
 353    ForceCalculator forceCalculator(simulationState, options);
 354
 355    // Integration requires masses, positions, and forces
 356    LeapFrog integrator(simulationState);
 357
 358    // Allocate a force buffer
 359    gmx::ArrayRef<gmx::RVec> userForces(topology.numParticles());
 360
 361    // MD Loop
 362    int numSteps = 100;
 363
 364    for (i = 0; i < numSteps; i++)
 365    {
 366      userForces = forceCalculator.compute();
 367
 368      // The forces are not automatically updated in case the user wants to add their own
 369      std::copy(userForces.begin(), userForces.end(), begin(simulationState.forces()));
 370
 371      // Integrate with a time step of 1 fs
 372      integrator.integrate(1.0);
 373    }
 374
 375    return 0;
 376    } // main
 377
 378 .. _doxygen-ref-rvec: ../doxygen/html-lib/namespacegmx.xhtml#a139c1919a9680de4ad1450f42e37d33b