using namespace nblib;
-// Main function to write the MD program.
-int main(); // Keep the compiler happy
-
int main()
{
// Create the particles
HarmonicBondType ohHarmonicBond(1, 1);
HarmonicBondType hcHarmonicBond(2, 1);
- HarmonicAngleType hohAngle(Degrees(120), 1);
- HarmonicAngleType hchAngle(Degrees(109.5), 1);
+ HarmonicAngle hohAngle(1, Degrees(120));
+ HarmonicAngle hchAngle(1, Degrees(109.5));
// add harmonic bonds for water
water.addInteraction(ParticleName("O"), ParticleName("H1"), ohHarmonicBond);
SimulationState simulationState(coordinates, velocities, forces, box, topology);
// The non-bonded force calculator contains all the data needed to compute forces
- ForceCalculator forceCalculator(simulationState, options);
+ auto forceCalculator = setupGmxForceCalculatorCpu(simulationState.topology(), options);
+
+ // build the pair list
+ forceCalculator->updatePairlist(simulationState.coordinates(), simulationState.box());
// The listed force calculator is also initialized with the required arguments
- ListedForceCalculator listedForceCalculator(topology.getInteractionData(),
- topology.numParticles(), 4, box);
+ ListedForceCalculator listedForceCalculator(
+ topology.getInteractionData(), topology.numParticles(), 4, box);
// Integrator is initialized with an array of inverse masses (constructed from topology) and
// the bounding box
- LeapFrog integrator(topology, box);
+ LeapFrog integrator(simulationState.topology(), simulationState.box());
// Print some diagnostic info
- printf("initial position of particle 0: x %4f y %4f z %4f\n", simulationState.coordinates()[0][0],
- simulationState.coordinates()[0][1], simulationState.coordinates()[0][2]);
+ printf("initial position of particle 0: x %4f y %4f z %4f\n",
+ simulationState.coordinates()[0][0],
+ simulationState.coordinates()[0][1],
+ simulationState.coordinates()[0][2]);
// MD Loop
int numSteps = 2;
{
zeroCartesianArray(simulationState.forces());
- forceCalculator.compute(simulationState.coordinates(), simulationState.forces());
+ forceCalculator->compute(
+ simulationState.coordinates(), simulationState.box(), simulationState.forces());
- listedForceCalculator.compute(simulationState.coordinates(), simulationState.forces());
+ listedForceCalculator.compute(
+ simulationState.coordinates(), simulationState.forces(), gmx::ArrayRef<real>{});
// Integrate with a time step of 1 fs, positions, velocities and forces
- integrator.integrate(1.0, simulationState.coordinates(), simulationState.velocities(),
- simulationState.forces());
+ integrator.integrate(
+ 1.0, simulationState.coordinates(), simulationState.velocities(), simulationState.forces());
}
- printf(" final position of particle 9: x %4f y %4f z %4f\n", simulationState.coordinates()[9][0],
- simulationState.coordinates()[9][1], simulationState.coordinates()[9][2]);
+ printf(" final position of particle 9: x %4f y %4f z %4f\n",
+ simulationState.coordinates()[9][0],
+ simulationState.coordinates()[9][1],
+ simulationState.coordinates()[9][2]);
return 0;
} // main