api/nblib/samples/methane-water-integration.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * This tests that sample code can run
  38  *
  39  * \author Victor Holanda <victor.holanda@cscs.ch>
  40  * \author Joe Jordan <ejjordan@kth.se>
  41  * \author Prashanth Kanduri <kanduri@cscs.ch>
  42  * \author Sebastian Keller <keller@cscs.ch>
  43  */
  44
  45 #include <cstdio>
  46
  47 #include "gromacs/utility/arrayref.h"
  48
  49 // The entire nblib public API can be included with a single header or individual components
  50 // can be included via their respective headers.
  51 #include "nblib/nblib.h"
  52
  53 using namespace nblib;
  54
  55 // Main function to write the MD program.
  56 int main(); // Keep the compiler happy
  57
  58 int main()
  59 {
  60     // Create the particles
  61     ParticleType Ow(ParticleTypeName("Ow"), Mass(15.999));
  62     ParticleType Hw(ParticleTypeName("Hw"), Mass(1.00784));
  63     ParticleType Cm(ParticleTypeName("Cm"), Mass(12.0107));
  64     ParticleType Hc(ParticleTypeName("Hc"), Mass(1.00784));
  65
  66     ParticleTypesInteractions interactions(CombinationRule::Geometric);
  67
  68     // Parameters from a GROMOS compatible force-field 2016H66
  69     // add non-bonded interactions for the particle types
  70     interactions.add(Ow.name(), C6(0.0026173456), C12(2.634129e-06));
  71     interactions.add(Hw.name(), C6(0.0), C12(0.0));
  72     interactions.add(Cm.name(), C6(0.01317904), C12(34.363044e-06));
  73     interactions.add(Hc.name(), C6(8.464e-05), C12(15.129e-09));
  74
  75     Molecule water(MoleculeName("Water"));
  76     Molecule methane(MoleculeName("Methane"));
  77
  78     water.addParticle(ParticleName("O"), Ow);
  79     water.addParticle(ParticleName("H1"), Hw);
  80     water.addParticle(ParticleName("H2"), Hw);
  81
  82     water.addExclusion(ParticleName("H1"), ParticleName("O"));
  83     water.addExclusion(ParticleName("H2"), ParticleName("O"));
  84
  85     methane.addParticle(ParticleName("C"), Cm);
  86     methane.addParticle(ParticleName("H1"), Hc);
  87     methane.addParticle(ParticleName("H2"), Hc);
  88     methane.addParticle(ParticleName("H3"), Hc);
  89     methane.addParticle(ParticleName("H4"), Hc);
  90
  91     methane.addExclusion(ParticleName("H1"), ParticleName("C"));
  92     methane.addExclusion(ParticleName("H2"), ParticleName("C"));
  93     methane.addExclusion(ParticleName("H3"), ParticleName("C"));
  94     methane.addExclusion(ParticleName("H4"), ParticleName("C"));
  95
  96     HarmonicBondType ohHarmonicBond(1, 1);
  97     HarmonicBondType hcHarmonicBond(2, 1);
  98
  99     HarmonicAngleType hohAngle(Degrees(120), 1);
 100     HarmonicAngleType hchAngle(Degrees(109.5), 1);
 101
 102     // add harmonic bonds for water
 103     water.addInteraction(ParticleName("O"), ParticleName("H1"), ohHarmonicBond);
 104     water.addInteraction(ParticleName("O"), ParticleName("H2"), ohHarmonicBond);
 105
 106     // add the angle for water
 107     water.addInteraction(ParticleName("H1"), ParticleName("O"), ParticleName("H2"), hohAngle);
 108
 109     // add harmonic bonds for methane
 110     methane.addInteraction(ParticleName("H1"), ParticleName("C"), hcHarmonicBond);
 111     methane.addInteraction(ParticleName("H2"), ParticleName("C"), hcHarmonicBond);
 112     methane.addInteraction(ParticleName("H3"), ParticleName("C"), hcHarmonicBond);
 113     methane.addInteraction(ParticleName("H4"), ParticleName("C"), hcHarmonicBond);
 114
 115     // add the angles for methane
 116     methane.addInteraction(ParticleName("H1"), ParticleName("C"), ParticleName("H2"), hchAngle);
 117     methane.addInteraction(ParticleName("H1"), ParticleName("C"), ParticleName("H3"), hchAngle);
 118     methane.addInteraction(ParticleName("H1"), ParticleName("C"), ParticleName("H4"), hchAngle);
 119     methane.addInteraction(ParticleName("H2"), ParticleName("C"), ParticleName("H3"), hchAngle);
 120     methane.addInteraction(ParticleName("H2"), ParticleName("C"), ParticleName("H4"), hchAngle);
 121     methane.addInteraction(ParticleName("H3"), ParticleName("C"), ParticleName("H4"), hchAngle);
 122
 123     // Define a box for the simulation
 124     Box box(6.05449);
 125
 126     // Define options for the non-bonded kernels
 127     NBKernelOptions options;
 128     // Use a simple cutoff rule for Coulomb
 129     options.coulombType = nblib::CoulombType::Cutoff;
 130     // Disable SIMD for this example
 131     options.nbnxmSimd = SimdKernels::SimdNo;
 132
 133     TopologyBuilder topologyBuilder;
 134
 135     // add molecules
 136     topologyBuilder.addMolecule(water, 2);
 137     topologyBuilder.addMolecule(methane, 1);
 138
 139     // add non-bonded interaction map
 140     topologyBuilder.addParticleTypesInteractions(interactions);
 141
 142     Topology topology = topologyBuilder.buildTopology();
 143
 144     // User defined coordinates.
 145     std::vector<Vec3> coordinates = {
 146         { 0.005, 0.600, 0.244 },    // Oxygen from water_1
 147         { -0.017, 0.690, 0.270 },   // Hydrogen_1 from water_1
 148         { 0.051, 0.610, 0.161 },    // Hydrogen_2 from water_1
 149         { 0.155, 0.341, 0.735 },    // Oxygen from water_2
 150         { 0.140, 0.284, 0.660 },    // Hydrogen_1 from water_2
 151         { 0.081, 0.402, 0.734 },    // Hydrogen_2 from water_2
 152         { -0.024, -0.222, -0.640 }, // Carbon from methane_1
 153         { -0.083, -0.303, -0.646 }, // Hydrogen_1 from methane_1
 154         { -0.080, -0.140, -0.642 }, // Hydrogen_2 from methane_1
 155         { 0.040, -0.221, -0.716 },  // Hydrogen_3 from methane_1
 156         { 0.027, -0.225, -0.553 }   // Hydrogen_4 from methane_1
 157     };
 158
 159     // User defined velocities.
 160     std::vector<Vec3> velocities = {
 161         { 0.1823, -0.4158, 0.487 },   // Oxygen from water_1
 162         { -1.7457, -0.5883, -0.460 }, // Hydrogen_1 from water_1
 163         { 2.5085, -0.1501, 1.762 },   // Hydrogen_2 from water_1
 164         { 0.6282, 0.4390, 0.001 },    // Oxygen from water_2
 165         { -0.3206, 0.0700, 0.4630 },  // Hydrogen_1 from water_2
 166         { -0.1556, -0.4529, 1.440 },  // Hydrogen_2 from water_2
 167         { 0.0, 0.0, 0.0 },            // Carbon from methane_1
 168         { 0.0, 0.0, 0.0 },            // Hydrogen_1 from methane_1
 169         { 0.0, 0.0, 0.0 },            // Hydrogen_2 from methane_1
 170         { 0.0, 0.0, 0.0 },            // Hydrogen_3 from methane_1
 171         { 0.0, 0.0, 0.0 },            // Hydrogen_4 from methane_1
 172     };
 173
 174     // Force buffer initialization for each particle.
 175     std::vector<Vec3> forces(topology.numParticles(), Vec3{ 0, 0, 0 });
 176
 177     SimulationState simulationState(coordinates, velocities, forces, box, topology);
 178
 179     // The non-bonded force calculator contains all the data needed to compute forces
 180     ForceCalculator forceCalculator(simulationState, options);
 181
 182     // The listed force calculator is also initialized with the required arguments
 183     ListedForceCalculator listedForceCalculator(topology.getInteractionData(),
 184                                                 topology.numParticles(), 4, box);
 185
 186     // Integrator is initialized with an array of inverse masses (constructed from topology) and
 187     // the bounding box
 188     LeapFrog integrator(topology, box);
 189
 190     // Print some diagnostic info
 191     printf("initial position of particle 0: x %4f y %4f z %4f\n", simulationState.coordinates()[0][0],
 192            simulationState.coordinates()[0][1], simulationState.coordinates()[0][2]);
 193
 194     // MD Loop
 195     int numSteps = 2;
 196
 197     for (auto i = 0; i < numSteps; i++)
 198     {
 199         zeroCartesianArray(simulationState.forces());
 200
 201         forceCalculator.compute(simulationState.coordinates(), simulationState.forces());
 202
 203         listedForceCalculator.compute(simulationState.coordinates(), simulationState.forces());
 204
 205         // Integrate with a time step of 1 fs, positions, velocities and forces
 206         integrator.integrate(1.0, simulationState.coordinates(), simulationState.velocities(),
 207                              simulationState.forces());
 208     }
 209
 210     printf("  final position of particle 9: x %4f y %4f z %4f\n", simulationState.coordinates()[9][0],
 211            simulationState.coordinates()[9][1], simulationState.coordinates()[9][2]);
 212
 213     return 0;
 214 } // main