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 int main()
  56 {
  57     // Create the particles
  58     ParticleType Ow(ParticleTypeName("Ow"), Mass(15.999));
  59     ParticleType Hw(ParticleTypeName("Hw"), Mass(1.00784));
  60     ParticleType Cm(ParticleTypeName("Cm"), Mass(12.0107));
  61     ParticleType Hc(ParticleTypeName("Hc"), Mass(1.00784));
  62
  63     ParticleTypesInteractions interactions(CombinationRule::Geometric);
  64
  65     // Parameters from a GROMOS compatible force-field 2016H66
  66     // add non-bonded interactions for the particle types
  67     interactions.add(Ow.name(), C6(0.0026173456), C12(2.634129e-06));
  68     interactions.add(Hw.name(), C6(0.0), C12(0.0));
  69     interactions.add(Cm.name(), C6(0.01317904), C12(34.363044e-06));
  70     interactions.add(Hc.name(), C6(8.464e-05), C12(15.129e-09));
  71
  72     Molecule water(MoleculeName("Water"));
  73     Molecule methane(MoleculeName("Methane"));
  74
  75     water.addParticle(ParticleName("O"), Ow);
  76     water.addParticle(ParticleName("H1"), Hw);
  77     water.addParticle(ParticleName("H2"), Hw);
  78
  79     water.addExclusion(ParticleName("H1"), ParticleName("O"));
  80     water.addExclusion(ParticleName("H2"), ParticleName("O"));
  81
  82     methane.addParticle(ParticleName("C"), Cm);
  83     methane.addParticle(ParticleName("H1"), Hc);
  84     methane.addParticle(ParticleName("H2"), Hc);
  85     methane.addParticle(ParticleName("H3"), Hc);
  86     methane.addParticle(ParticleName("H4"), Hc);
  87
  88     methane.addExclusion(ParticleName("H1"), ParticleName("C"));
  89     methane.addExclusion(ParticleName("H2"), ParticleName("C"));
  90     methane.addExclusion(ParticleName("H3"), ParticleName("C"));
  91     methane.addExclusion(ParticleName("H4"), ParticleName("C"));
  92
  93     HarmonicBondType ohHarmonicBond(1, 1);
  94     HarmonicBondType hcHarmonicBond(2, 1);
  95
  96     HarmonicAngle hohAngle(1, Degrees(120));
  97     HarmonicAngle hchAngle(1, Degrees(109.5));
  98
  99     // add harmonic bonds for water
 100     water.addInteraction(ParticleName("O"), ParticleName("H1"), ohHarmonicBond);
 101     water.addInteraction(ParticleName("O"), ParticleName("H2"), ohHarmonicBond);
 102
 103     // add the angle for water
 104     water.addInteraction(ParticleName("H1"), ParticleName("O"), ParticleName("H2"), hohAngle);
 105
 106     // add harmonic bonds for methane
 107     methane.addInteraction(ParticleName("H1"), ParticleName("C"), hcHarmonicBond);
 108     methane.addInteraction(ParticleName("H2"), ParticleName("C"), hcHarmonicBond);
 109     methane.addInteraction(ParticleName("H3"), ParticleName("C"), hcHarmonicBond);
 110     methane.addInteraction(ParticleName("H4"), ParticleName("C"), hcHarmonicBond);
 111
 112     // add the angles for methane
 113     methane.addInteraction(ParticleName("H1"), ParticleName("C"), ParticleName("H2"), hchAngle);
 114     methane.addInteraction(ParticleName("H1"), ParticleName("C"), ParticleName("H3"), hchAngle);
 115     methane.addInteraction(ParticleName("H1"), ParticleName("C"), ParticleName("H4"), hchAngle);
 116     methane.addInteraction(ParticleName("H2"), ParticleName("C"), ParticleName("H3"), hchAngle);
 117     methane.addInteraction(ParticleName("H2"), ParticleName("C"), ParticleName("H4"), hchAngle);
 118     methane.addInteraction(ParticleName("H3"), ParticleName("C"), ParticleName("H4"), hchAngle);
 119
 120     // Define a box for the simulation
 121     Box box(6.05449);
 122
 123     // Define options for the non-bonded kernels
 124     NBKernelOptions options;
 125     // Use a simple cutoff rule for Coulomb
 126     options.coulombType = nblib::CoulombType::Cutoff;
 127     // Disable SIMD for this example
 128     options.nbnxmSimd = SimdKernels::SimdNo;
 129
 130     TopologyBuilder topologyBuilder;
 131
 132     // add molecules
 133     topologyBuilder.addMolecule(water, 2);
 134     topologyBuilder.addMolecule(methane, 1);
 135
 136     // add non-bonded interaction map
 137     topologyBuilder.addParticleTypesInteractions(interactions);
 138
 139     Topology topology = topologyBuilder.buildTopology();
 140
 141     // User defined coordinates.
 142     std::vector<Vec3> coordinates = {
 143         { 0.005, 0.600, 0.244 },    // Oxygen from water_1
 144         { -0.017, 0.690, 0.270 },   // Hydrogen_1 from water_1
 145         { 0.051, 0.610, 0.161 },    // Hydrogen_2 from water_1
 146         { 0.155, 0.341, 0.735 },    // Oxygen from water_2
 147         { 0.140, 0.284, 0.660 },    // Hydrogen_1 from water_2
 148         { 0.081, 0.402, 0.734 },    // Hydrogen_2 from water_2
 149         { -0.024, -0.222, -0.640 }, // Carbon from methane_1
 150         { -0.083, -0.303, -0.646 }, // Hydrogen_1 from methane_1
 151         { -0.080, -0.140, -0.642 }, // Hydrogen_2 from methane_1
 152         { 0.040, -0.221, -0.716 },  // Hydrogen_3 from methane_1
 153         { 0.027, -0.225, -0.553 }   // Hydrogen_4 from methane_1
 154     };
 155
 156     // User defined velocities.
 157     std::vector<Vec3> velocities = {
 158         { 0.1823, -0.4158, 0.487 },   // Oxygen from water_1
 159         { -1.7457, -0.5883, -0.460 }, // Hydrogen_1 from water_1
 160         { 2.5085, -0.1501, 1.762 },   // Hydrogen_2 from water_1
 161         { 0.6282, 0.4390, 0.001 },    // Oxygen from water_2
 162         { -0.3206, 0.0700, 0.4630 },  // Hydrogen_1 from water_2
 163         { -0.1556, -0.4529, 1.440 },  // Hydrogen_2 from water_2
 164         { 0.0, 0.0, 0.0 },            // Carbon from methane_1
 165         { 0.0, 0.0, 0.0 },            // Hydrogen_1 from methane_1
 166         { 0.0, 0.0, 0.0 },            // Hydrogen_2 from methane_1
 167         { 0.0, 0.0, 0.0 },            // Hydrogen_3 from methane_1
 168         { 0.0, 0.0, 0.0 },            // Hydrogen_4 from methane_1
 169     };
 170
 171     // Force buffer initialization for each particle.
 172     std::vector<Vec3> forces(topology.numParticles(), Vec3{ 0, 0, 0 });
 173
 174     SimulationState simulationState(coordinates, velocities, forces, box, topology);
 175
 176     // The non-bonded force calculator contains all the data needed to compute forces
 177     auto forceCalculator = setupGmxForceCalculatorCpu(simulationState.topology(), options);
 178
 179     // build the pair list
 180     forceCalculator->updatePairlist(simulationState.coordinates(), simulationState.box());
 181
 182     // The listed force calculator is also initialized with the required arguments
 183     ListedForceCalculator listedForceCalculator(
 184             topology.getInteractionData(), 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(simulationState.topology(), simulationState.box());
 189
 190     // Print some diagnostic info
 191     printf("initial position of particle 0: x %4f y %4f z %4f\n",
 192            simulationState.coordinates()[0][0],
 193            simulationState.coordinates()[0][1],
 194            simulationState.coordinates()[0][2]);
 195
 196     // MD Loop
 197     int numSteps = 2;
 198
 199     for (auto i = 0; i < numSteps; i++)
 200     {
 201         zeroCartesianArray(simulationState.forces());
 202
 203         forceCalculator->compute(
 204                 simulationState.coordinates(), simulationState.box(), simulationState.forces());
 205
 206         listedForceCalculator.compute(simulationState.coordinates(), simulationState.forces());
 207
 208         // Integrate with a time step of 1 fs, positions, velocities and forces
 209         integrator.integrate(
 210                 1.0, simulationState.coordinates(), simulationState.velocities(), simulationState.forces());
 211     }
 212
 213     printf("  final position of particle 9: x %4f y %4f z %4f\n",
 214            simulationState.coordinates()[9][0],
 215            simulationState.coordinates()[9][1],
 216            simulationState.coordinates()[9][2]);
 217
 218     return 0;
 219 } // main