2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * This tests that sample code can run
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>
47 #include "gromacs/utility/arrayref.h"
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"
53 using namespace nblib;
55 // Main function to write the MD program.
56 int main(); // Keep the compiler happy
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));
66 ParticleTypesInteractions interactions(CombinationRule::Geometric);
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));
75 Molecule water(MoleculeName("Water"));
76 Molecule methane(MoleculeName("Methane"));
78 water.addParticle(ParticleName("O"), Ow);
79 water.addParticle(ParticleName("H1"), Hw);
80 water.addParticle(ParticleName("H2"), Hw);
82 water.addExclusion(ParticleName("H1"), ParticleName("O"));
83 water.addExclusion(ParticleName("H2"), ParticleName("O"));
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);
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"));
96 HarmonicBondType ohHarmonicBond(1, 1);
97 HarmonicBondType hcHarmonicBond(2, 1);
99 HarmonicAngle hohAngle(Degrees(120), 1);
100 HarmonicAngle hchAngle(Degrees(109.5), 1);
102 // add harmonic bonds for water
103 water.addInteraction(ParticleName("O"), ParticleName("H1"), ohHarmonicBond);
104 water.addInteraction(ParticleName("O"), ParticleName("H2"), ohHarmonicBond);
106 // add the angle for water
107 water.addInteraction(ParticleName("H1"), ParticleName("O"), ParticleName("H2"), hohAngle);
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);
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);
123 // Define a box for the simulation
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;
133 TopologyBuilder topologyBuilder;
136 topologyBuilder.addMolecule(water, 2);
137 topologyBuilder.addMolecule(methane, 1);
139 // add non-bonded interaction map
140 topologyBuilder.addParticleTypesInteractions(interactions);
142 Topology topology = topologyBuilder.buildTopology();
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
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
174 // Force buffer initialization for each particle.
175 std::vector<Vec3> forces(topology.numParticles(), Vec3{ 0, 0, 0 });
177 SimulationState simulationState(coordinates, velocities, forces, box, topology);
179 // The non-bonded force calculator contains all the data needed to compute forces
180 ForceCalculator forceCalculator(simulationState, options);
182 // The listed force calculator is also initialized with the required arguments
183 ListedForceCalculator listedForceCalculator(
184 topology.getInteractionData(), topology.numParticles(), 4, box);
186 // Integrator is initialized with an array of inverse masses (constructed from topology) and
188 LeapFrog integrator(topology, box);
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]);
199 for (auto i = 0; i < numSteps; i++)
201 zeroCartesianArray(simulationState.forces());
203 forceCalculator.compute(simulationState.coordinates(), simulationState.forces());
205 listedForceCalculator.compute(simulationState.coordinates(), simulationState.forces());
207 // Integrate with a time step of 1 fs, positions, velocities and forces
208 integrator.integrate(
209 1.0, simulationState.coordinates(), simulationState.velocities(), simulationState.forces());
212 printf(" final position of particle 9: x %4f y %4f z %4f\n",
213 simulationState.coordinates()[9][0],
214 simulationState.coordinates()[9][1],
215 simulationState.coordinates()[9][2]);