Skip to content

GitLab

Unverified Commit 389f79e6 authored by peastman's avatar peastman Committed by GitHub
Browse files

Minor fixes to Nose Hoover code (#2670) 

* Fixed OpenCL compiler errors

* Fixed incorrect references to BAOAB

* Fixed a method that was incorrectly made public
parent 5f374e1d
Hide whitespace changes
Inline Side-by-side
Showing with 40 additions and 80 deletions
openmmapi/include/openmm/NoseHooverIntegrator.h
View file @ 389f79e6
......@@ -72,7 +72,7 @@ public:
explicit NoseHooverIntegrator(double temperature, double collisionFrequency, double stepSize,
int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 7);
virtual ~NoseHooverIntegrator();
~NoseHooverIntegrator();
/**
* Advance a simulation through time by taking a series of time steps.
*
......@@ -195,18 +195,6 @@ public:
* @param chainID the index of the Nose-Hoover chain thermostat (default=0).
*/
const NoseHooverChain& getThermostat(int chainID=0) const ;
/**
* This will be called by the Context when the user modifies aspects of the context state, such
* as positions, velocities, or parameters. This gives the Integrator a chance to discard cached
* information. This is <i>only</i> called when the user modifies information using methods of the Context
* object. It is <i>not</i> called when a ForceImpl object modifies state information in its updateContextState()
* method (unless the ForceImpl calls a Context method to perform the modification).
*
* @param changed this specifies what aspect of the Context was changed
*/
virtual void stateChanged(State::DataType changed) {
if (State::Positions == changed) forcesAreValid = false;
}
/**
* Return false, if this integrator was set up with the 'default constructor' that thermostats the whole system,
* true otherwise. Required for serialization.
......@@ -248,6 +236,18 @@ protected:
* cleanup. It will also get called again if the application calls reinitialize() on the Context.
*/
void cleanup();
/**
* This will be called by the Context when the user modifies aspects of the context state, such
* as positions, velocities, or parameters. This gives the Integrator a chance to discard cached
* information. This is <i>only</i> called when the user modifies information using methods of the Context
* object. It is <i>not</i> called when a ForceImpl object modifies state information in its updateContextState()
* method (unless the ForceImpl calls a Context method to perform the modification).
*
* @param changed this specifies what aspect of the Context was changed
*/
void stateChanged(State::DataType changed) {
if (State::Positions == changed) forcesAreValid = false;
}
/**
* Get the names of all Kernels used by this Integrator.
*/
......@@ -255,11 +255,11 @@ protected:
/**
* Compute the kinetic energy of the system at the current time.
*/
virtual double computeKineticEnergy();
double computeKineticEnergy();
/**
* Computing kinetic energy for this integrator does not require forces.
*/
bool kineticEnergyRequiresForce() const override;
bool kineticEnergyRequiresForce() const;
std::vector<NoseHooverChain> noseHooverChains;
std::vector<int> allAtoms;
......
platforms/common/src/CommonKernels.cpp
View file @ 389f79e6
......@@ -5791,15 +5791,14 @@ void CommonIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const No
}
/*
* Carry out the BAOAB integration, also known as the "LF-middle" scheme
* c.f. J. Phys. Chem. A 2019, 123, 6056−6079
* Carry out the LF-middle integration (c.f. J. Phys. Chem. A 2019, 123, 6056−6079)
*/
// B
// Velocity update
kernel1->execute(std::max(numAtoms, numPairs));
integration.applyVelocityConstraints(integrator.getConstraintTolerance());
// A
// Position update
kernel2->execute(numParticles);
// O
// Apply the thermostat
int numChains = integrator.getNumThermostats();
for(int chain = 0; chain < numChains; ++chain) {
const auto &thermostatChain = integrator.getThermostat(chain);
......@@ -5807,10 +5806,10 @@ void CommonIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const No
auto scaleFactors = propagateChain(context, thermostatChain, KEs, dt);
scaleVelocities(context, thermostatChain, scaleFactors);
}
// A
// Position update
kernel3->execute(numParticles);
integration.applyConstraints(integrator.getConstraintTolerance());
// B
// Apply constraint forces
kernel4->execute(numAtoms);
// Make sure any Drude-like particles have not wandered too far from home
if (numPairs > 0) kernelHardWall->execute(numPairs);
......
platforms/common/src/kernels/noseHooverChain.cc
View file @ 389f79e6
......@@ -18,12 +18,11 @@ KERNEL void propagateNoseHooverChain(GLOBAL mixed2* RESTRICT chainData, GLOBAL c
mixed wdt = ys * timeOverMTS;
chainData[chainLength-1].y += 0.5f * wdt * chainForces[chainLength-1];
for (int bead = chainLength - 2; bead >= 0; --bead) {
mixed aa = EXP(-0.25f * wdt * chainData[bead + 1].y);
mixed aa = exp(-0.25f * wdt * chainData[bead + 1].y);
chainData[bead].y = aa * (chainData[bead].y * aa + 0.5f * wdt * chainForces[bead]);
}
// update particle velocities
mixed aa = EXP(-wdt * chainData[0].y);
scale *= aa;
scale *= (mixed) exp(-wdt * chainData[0].y);;
// update the thermostat positions
for (int bead = 0; bead < chainLength; ++bead) {
chainData[bead].x += chainData[bead].y * wdt;
......@@ -32,7 +31,7 @@ KERNEL void propagateNoseHooverChain(GLOBAL mixed2* RESTRICT chainData, GLOBAL c
chainForces[0] = (scale * scale * KE2 - numDOFs * kT) / chainMasses[0];
// update thermostat velocities
for (int bead = 0; bead < chainLength - 1; ++bead) {
mixed aa = EXP(-0.25f * wdt * chainData[bead + 1].y);
mixed aa = exp(-0.25f * wdt * chainData[bead + 1].y);
chainData[bead].y = aa * (aa * chainData[bead].y + 0.5f * wdt * chainForces[bead]);
chainForces[bead + 1] = (chainMasses[bead] * chainData[bead].y * chainData[bead].y - kT) / chainMasses[bead + 1];
}
......
platforms/reference/include/ReferenceNoseHooverDynamics.h
View file @ 389f79e6
/* Portions copyright (c) 2006-2012 Stanford University and Simbios.
* Contributors: Pande Group
/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
* Contributors: Andy Simmonett, Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
......@@ -65,7 +65,7 @@ class ReferenceNoseHooverDynamics : public ReferenceDynamics {
/**---------------------------------------------------------------------------------------
BAstep = update the velocities and positions to start a BAOAB step using the leapfrog "middle" scheme
Perform the first half of a step using the leapfrog LF-Middle scheme
@param system the System to be integrated
@param atomCoordinates atom coordinates
......@@ -80,12 +80,12 @@ class ReferenceNoseHooverDynamics : public ReferenceDynamics {
--------------------------------------------------------------------------------------- */
void BAstep(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance, bool &forcesAreValid,
const std::vector<int> & allAtoms, const std::vector<std::tuple<int, int, double>> & allPairs, double maxPairDistance);
void step1(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance, bool &forcesAreValid,
const std::vector<int> & allAtoms, const std::vector<std::tuple<int, int, double>> & allPairs, double maxPairDistance);
/**---------------------------------------------------------------------------------------
ABstep = update the positions and velocities to end a BAOAB step using the leapfrog "middle" scheme
Perform the second half of a step using the leapfrog LF-Middle scheme
@param system the System to be integrated
@param atomCoordinates atom coordinates
......@@ -99,9 +99,9 @@ class ReferenceNoseHooverDynamics : public ReferenceDynamics {
@param maxPairDistance the maximum separation allowed for a Drude-like pair
--------------------------------------------------------------------------------------- */
void ABstep(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance, bool &forcesAreValid,
const std::vector<int> & allAtoms, const std::vector<std::tuple<int, int, double>> & allPairs, double maxPairDistance);
void step2(OpenMM::ContextImpl &context, const OpenMM::System& system, std::vector<OpenMM::Vec3>& atomCoordinates,
std::vector<OpenMM::Vec3>& velocities, std::vector<OpenMM::Vec3>& forces, std::vector<double>& masses, double tolerance, bool &forcesAreValid,
const std::vector<int> & allAtoms, const std::vector<std::tuple<int, int, double>> & allPairs, double maxPairDistance);
};
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 389f79e6
......@@ -2181,11 +2181,8 @@ void ReferenceIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const
dynamics->setReferenceConstraintAlgorithm(&extractConstraints(context));
prevStepSize = stepSize;
}
// BA
dynamics->BAstep(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
dynamics->step1(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
integrator.getAllThermostatedIndividualParticles(), integrator.getAllThermostatedPairs(), integrator.getMaximumPairDistance());
// O
int numChains = integrator.getNumThermostats();
for(int chain = 0; chain < numChains; ++chain) {
const auto &thermostatChain = integrator.getThermostat(chain);
......@@ -2193,10 +2190,8 @@ void ReferenceIntegrateNoseHooverStepKernel::execute(ContextImpl& context, const
std::pair<double, double> scaleFactors = propagateChain(context, thermostatChain, KEs, stepSize);
scaleVelocities(context, thermostatChain, scaleFactors);
}
// AB
dynamics->ABstep(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
dynamics->step2(context, context.getSystem(), posData, velData, forceData, masses, integrator.getConstraintTolerance(), forcesAreValid,
integrator.getAllThermostatedIndividualParticles(), integrator.getAllThermostatedPairs(), integrator.getMaximumPairDistance());
data.time += stepSize;
data.stepCount++;
}
......
platforms/reference/src/SimTKReference/ReferenceNoseHooverDynamics.cpp
View file @ 389f79e6
/* Portions copyright (c) 2006-2013 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group
/* Portions copyright (c) 2006-2020 Stanford University and Simbios.
* Contributors: Andy Simmonett, Peter Eastman, Pande Group
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
......@@ -37,17 +37,6 @@
using std::vector;
using namespace OpenMM;
/**---------------------------------------------------------------------------------------
ReferenceNoseHooverDynamics constructor
@param numberOfAtoms number of atoms
@param deltaT delta t for dynamics
@param friction friction coefficient
@param temperature temperature
--------------------------------------------------------------------------------------- */
ReferenceNoseHooverDynamics::ReferenceNoseHooverDynamics(int numberOfAtomsIn, double deltaT) :
ReferenceDynamics(numberOfAtomsIn, deltaT, 0.0) {
numberOfAtoms = numberOfAtomsIn;
......@@ -56,32 +45,10 @@ ReferenceNoseHooverDynamics::ReferenceNoseHooverDynamics(int numberOfAtomsIn, do
oldx.resize(numberOfAtoms);
}
/**---------------------------------------------------------------------------------------
ReferenceNoseHooverDynamics destructor
--------------------------------------------------------------------------------------- */
ReferenceNoseHooverDynamics::~ReferenceNoseHooverDynamics() {
}
/**---------------------------------------------------------------------------------------
Update -- driver routine for performing Velocity Verlet dynamics update of coordinates
and velocities
@param system the System to be integrated
@param atomCoordinates atom coordinates
@param velocities velocities
@param forces forces
@param masses atom masses
@param atomList list of all atoms not involved in a Drude-like pair
@param pairList list of all Drude-like pairs
@param maxPairDistance the maximum separation of any Drude-like pairs
--------------------------------------------------------------------------------------- */
void ReferenceNoseHooverDynamics::BAstep(OpenMM::ContextImpl &context, const OpenMM::System& system, vector<Vec3>& atomCoordinates,
void ReferenceNoseHooverDynamics::step1(OpenMM::ContextImpl &context, const OpenMM::System& system, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& masses, double tolerance, bool &forcesAreValid,
const std::vector<int> & atomList, const std::vector<std::tuple<int, int, double>> &pairList,
......@@ -145,7 +112,7 @@ void ReferenceNoseHooverDynamics::BAstep(OpenMM::ContextImpl &context, const Ope
}
void ReferenceNoseHooverDynamics::ABstep(OpenMM::ContextImpl &context, const OpenMM::System& system, vector<Vec3>& atomCoordinates,
void ReferenceNoseHooverDynamics::step2(OpenMM::ContextImpl &context, const OpenMM::System& system, vector<Vec3>& atomCoordinates,
vector<Vec3>& velocities,
vector<Vec3>& forces, vector<double>& masses, double tolerance, bool &forcesAreValid,
const std::vector<int> & atomList, const std::vector<std::tuple<int, int, double>> &pairList,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment