Skip to content

GitLab

Unverified Commit 50efef13 authored by Andy Simmonett's avatar Andy Simmonett
Browse files

Start implementing changes from review - CUDA broken

parent 3251d4e4
Hide whitespace changes
Inline Side-by-side
Showing with 643 additions and 248 deletions
olla/include/openmm/kernels.h
View file @ 50efef13
......@@ -64,7 +64,7 @@
#include "openmm/VariableLangevinIntegrator.h"
#include "openmm/VariableVerletIntegrator.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/VelocityVerletIntegrator.h"
#include "openmm/NoseHooverIntegrator.h"
#include <iosfwd>
#include <set>
#include <string>
......@@ -1075,7 +1075,7 @@ public:
* @param system the System this kernel will be applied to
* @param integrator the VelocityVerletIntegrator this kernel will be used for
*/
virtual void initialize(const System& system, const VelocityVerletIntegrator& integrator) = 0;
virtual void initialize(const System& system, const NoseHooverIntegrator& integrator) = 0;
/**
* Execute the kernel.
*
......@@ -1084,14 +1084,14 @@ public:
* @param forcesAreValid a reference to the parent integrator's boolean for keeping
* track of the validity of the current forces.
*/
virtual void execute(ContextImpl& context, const VelocityVerletIntegrator& integrator, bool &forcesAreValid) = 0;
virtual void execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) = 0;
/**
* Compute the kinetic energy.
*
* @param context the context in which to execute this kernel
* @param integrator the VelocityVerletIntegrator this kernel is being used for
*/
virtual double computeKineticEnergy(ContextImpl& context, const VelocityVerletIntegrator& integrator) = 0;
virtual double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) = 0;
};
/**
......@@ -1346,11 +1346,11 @@ public:
*
* @param context the context in which to execute this kernel
* @param noseHooverChain the object describing the chain to be propagated.
* @param kineticEnergy the kineticEnergy of the particles being thermostated by this chain.
* @param kineticEnergy the {center of mass, relative} kineticEnergies of the particles being thermostated by this chain.
* @param timeStep the time step used by the integrator.
* @return the velocity scale factor to apply to the particles associated with this heat bath.
*/
virtual double propagateChain(ContextImpl& context, const NoseHooverChain &noseHooverChain, double kineticEnergy, double timeStep) = 0;
virtual std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> kineticEnergy, double timeStep) = 0;
/**
* Execute the kernal that computes the total (kinetic + potential) heat bath energy.
*
......@@ -1367,15 +1367,15 @@ public:
* @param noseHooverChain the chain whose energy is to be determined.
* @param downloadValue whether the computed value should be downloaded and returned.
*/
virtual double computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) = 0;
virtual std::pair<double, double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) = 0;
/**
* Execute the kernel that scales the velocities of particles associated with a nose hoover chain
*
* @param context the context in which to execute this kernel
* @param noseHooverChain the chain whose energy is to be determined.
* @param scaleFactor the multiplicative factor by which velocities are scaled.
* @param scaleFactor the multiplicative factor by which {absolute, relative} velocities are scaled.
*/
virtual void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, double scaleFactor) = 0;
virtual void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactor) = 0;
};
/**
......
openmmapi/include/OpenMM.h
View file @ 50efef13
......@@ -75,7 +75,7 @@
#include "openmm/VariableVerletIntegrator.h"
#include "openmm/Vec3.h"
#include "openmm/VerletIntegrator.h"
#include "openmm/VelocityVerletIntegrator.h"
#include "openmm/NoseHooverIntegrator.h"
#include "openmm/VirtualSite.h"
#include "openmm/Platform.h"
#include "openmm/serialization/XmlSerializer.h"
......
openmmapi/include/openmm/NoseHooverChain.h
View file @ 50efef13
......@@ -49,6 +49,12 @@ namespace OpenMM {
*
* where the total number of timesteps used to propagate the chain in each step is
* the number of MTS steps multiplied by the number of terms in the Yoshida-Suzuki decomposition.
*
* Two types of NHC may be created. The first is a simple thermostat that couples with a given subset
* of the atoms within a system, controling their absolute motion. The second is more elaborate and
* can thermostat tethered pairs of atoms and in this case two thermostats are created: one that controls
* the absolute center of mass velocity of each pair and another that controls their motion relative to
* one another.
*/
class OPENMM_EXPORT NoseHooverChain {
......@@ -56,23 +62,29 @@ public:
/**
* Create a NoseHooverChain.
*
* @param defaultTemperature the default temperature of the heat bath (in Kelvin)
* @param defaultCollisionFrequency the default collision frequency (in 1/ps)
* @param defaultNumDOFs the default number of degrees of freedom in the particles that
* interact with this chain
* @param defaultChainLength the default length of (number of particles in) this heat bath
* @param defaultNumMTS the default number of multi time steps used to propagate this chain
* @param defaultNumYoshidaSuzuki the default number of Yoshida Suzuki steps used to propagate this chain (1, 3, or 5).
* @param defaultChainID the default chain id used to distinguish this Nose-Hoover chain from others that may
* be used to control a different set of particles, e.g. for Drude oscillators
* @param thermostatedAtoms the list of atoms to be handled by this thermostat
* @param parentAtoms the list of parent atoms, if this thermostat handles Drude particles. Empty vector otherwise.
* @param defaultTemperature the default temperature of the heat bath for absolute motion (in Kelvin)
* @param defaultRelativeTemperature the default temperature of the heat bath for relative motion(in Kelvin).
* This is only used if the list of thermostated pairs is not empty.
* @param defaultCollisionFrequency the default collision frequency for absolute motion (in 1/ps)
* @param defaultRelativeCollisionFrequency the default collision frequency for relative motion(in 1/ps).
* This is only used if the list of thermostated pairs is not empty.
* @param defaultNumDOFs the default number of degrees of freedom in the particles that
* interact with this chain
* @param defaultChainLength the default length of (number of particles in) this heat bath
* @param defaultNumMTS the default number of multi time steps used to propagate this chain
* @param defaultNumYoshidaSuzuki the default number of Yoshida Suzuki steps used to propagate this chain (1, 3, or 5).
* @param defaultChainID the default chain id used to distinguish this Nose-Hoover chain from others that may
* be used to control a different set of particles, e.g. for Drude oscillators
* @param thermostatedAtoms the list of atoms to be handled by this thermostat
* @param thermostatedPairs the list of connected pairs to be thermostated; their absolute center of mass motion will
* be thermostated independently from their motion relative to one another.
*/
NoseHooverChain(double defaultTemperature, double defaultCollisionFrequency, int defaultNumDOFs, int defaultChainLength,
NoseHooverChain(double defaultTemperature, double defaultRelativeTemperature, double defaultCollisionFrequency,
double defaultRelativeCollisionFrequency, int defaultNumDOFs, int defaultChainLength,
int defaultNumMTS, int defaultNumYoshidaSuzuki, int defaultChainID,
const std::vector<int>& thermostatedAtoms, const std::vector<int>& parentAtoms);
const std::vector<int>& thermostatedAtoms, const std::vector< std::pair< int, int > > &thermostatedPairs);
/**
* Get the default temperature of the heat bath (in Kelvin).
* Get the default temperature of the heat bath for treating absolute particle motion (in Kelvin).
*
* @return the default temperature of the heat bath, measured in Kelvin.
*/
......@@ -80,8 +92,8 @@ public:
return defaultTemp;
}
/**
* Set the default temperature of the heat bath. This will affect any new Contexts you create,
* but not ones that already exist.
* Set the default temperature of the heat bath for treating absolute particle motion.
* This will affect any new Contexts you create, but not ones that already exist.
*
* @param temperature the default temperature of the heat bath (in Kelvin)
*/
......@@ -89,7 +101,25 @@ public:
defaultTemp = temperature;
}
/**
* Get the default collision frequency (in 1/ps).
* Get the default temperature of the heat bath for treating relative particle motion (in Kelvin).
*
* @return the default temperature of the heat bath, measured in Kelvin.
*/
double getDefaultRelativeTemperature() const {
return defaultRelativeTemp;
}
/**
* Set the default temperature of the heat bath for treating relative motion if this thermostat has
* been set up to treat connected pairs of atoms. This will affect any new Contexts you create,
* but not ones that already exist.
*
* @param temperature the default temperature of the heat bath for relative motion (in Kelvin)
*/
void setDefaultRelativeTemperature(double temperature) {
defaultRelativeTemp = temperature;
}
/**
* Get the default collision frequency for treating absolute particle motion (in 1/ps).
*
* @return the default collision frequency, measured in 1/ps.
*/
......@@ -97,14 +127,32 @@ public:
return defaultFreq;
}
/**
* Set the default collision frequency. This will affect any new Contexts you create,
* but not those that already exist.
* Set the default collision frequency for treating absolute particle motion.
* This will affect any new Contexts you create, but not those that already exist.
*
* @param frequency the default collision frequency (in 1/ps)
*/
void setDefaultCollisionFrequency(double frequency) {
defaultFreq = frequency;
}
/**
* Get the default collision frequency for treating relative particle motion (in 1/ps).
*
* @return the default collision frequency, measured in 1/ps.
*/
double getDefaultRelativeCollisionFrequency() const {
return defaultRelativeFreq;
}
/**
* Set the default collision frequency for treating relative particle motion if this thermostat has
* been set up to handle connected pairs of atoms. This will affect any new Contexts you create,
* but not those that already exist.
*
* @param frequency the default collision frequency (in 1/ps)
*/
void setDefaultRelativeCollisionFrequency(double frequency) {
defaultRelativeFreq = frequency;
}
/**
* Get the default number of degrees of freedom in the particles controled by this heat bath.
*
......@@ -209,26 +257,22 @@ public:
thermostatedAtoms = atomIDs;
}
/**
* In case this thermostat handles the kinetic energy of Drude particles relative to the
* atoms that they are attached to, get the atom ids of all parent atoms.
* Get the list of any connected pairs to be handled by this thermostat.
* If this is a regular thermostat, returns an empty vector.
*
* @returns ids of all parent atoms
* @returns list of connected pairs.
*/
const std::vector<int>& getParentAtoms() const {
return parentAtoms;
const std::vector< std::pair< int, int > >& getThermostatedPairs() const {
return thermostatedPairs;
}
/**
* In case this thermostat handles the kinetic energy of Drude particles
* set the atom IDs of all parent atoms.
*
* @param parentIDs
* @param pairIDs the list of connected pairs to thermostat.
*/
void setParentAtoms(const std::vector<int>& parentIDs){
if (not parentIDs.size() == thermostatedAtoms.size()){
throw OpenMMException("The number of parent atoms has to be the same as the number of thermostated atoms, or zero.");
}
parentAtoms = parentIDs;
void setThermostatedPairs(const std::vector< std::pair< int, int > >& pairIDs){
thermostatedPairs = pairIDs;
}
/**
* Get the default weights used in the Yoshida Suzuki multi time step decomposition (dimensionless)
......@@ -246,11 +290,12 @@ public:
return false;
}
private:
double defaultTemp, defaultFreq; //, defaultTimeStep;
double defaultTemp, defaultFreq, defaultRelativeTemp, defaultRelativeFreq;
int defaultNumDOFs, defaultChainLength, defaultNumMTS, defaultNumYS;
// The suffix used to distinguish NH chains, e.g. for Drude particles vs. regular particles.
int defaultChainID;
std::vector<int> thermostatedAtoms, parentAtoms;
std::vector<int> thermostatedAtoms;
std::vector<std::pair<int, int>> thermostatedPairs;
};
} // namespace OpenMM
......
openmmapi/include/openmm/VelocityVerletIntegrator.h → openmmapi/include/openmm/NoseHooverIntegrator.h
View file @ 50efef13
#ifndef OPENMM_VELOCITYVERLETINTEGRATOR_H_
#define OPENMM_VELOCITYVERLETINTEGRATOR_H_
#ifndef OPENMM_NOSEHOOVERINTEGRATOR_H_
#define OPENMM_NOSEHOOVERINTEGRATOR_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
......@@ -47,16 +47,32 @@ class System;
* thermostats, using the velocity Verlet propagation algorithm.
*/
class OPENMM_EXPORT VelocityVerletIntegrator : public Integrator {
class OPENMM_EXPORT NoseHooverIntegrator : public Integrator {
public:
/**
* Create a VelocityVerletIntegrator.
* Create a NoseHooverIntegrator. This version creates a bare velocity Verlet integrator
* with no thermostats; any thermostats should be added by calling addThermostat.
*
* @param stepSize the step size with which to integrate the system (in picoseconds)
*/
explicit VelocityVerletIntegrator(double stepSize);
explicit NoseHooverIntegrator(double stepSize);
/**
* Create a NoseHooverIntegrator.
*
* @param stepSize the step size with which to integrate the system (in picoseconds)
* @param system the system to be thermostated. Note: this must be setup, i.e. all
* particles should have been added, before calling this function.
* @param temperature the target temperature for the system.
* @param collisionFrequency the frequency of the interaction with the heat bath (in 1/ps).
* @param chainLength the number of beads in the Nose-Hoover chain.
* @param numMTS the number of step in the multiple time step chain propagation algorithm.
* @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
* used in the chain propagation algorithm (must be 1, 3, or 5).
*/
explicit NoseHooverIntegrator(double stepSize, System &system, double temperature, int collisionFrequnency,
int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 3);
virtual ~VelocityVerletIntegrator();
virtual ~NoseHooverIntegrator();
/**
* Advance a simulation through time by taking a series of time steps.
*
......@@ -64,7 +80,7 @@ public:
*/
void step(int steps);
/**
* Add a Nose-Hoover Chain thermostat to control the temperature of the system
* Add a simple Nose-Hoover Chain thermostat to control the temperature of the full system
*
* @param system the system to be thermostated. Note: this must be setup, i.e. all
* particles should have been added, before calling this function.
......@@ -75,28 +91,38 @@ public:
* @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
* used in the chain propagation algorithm (must be 1, 3, or 5).
*/
int addNoseHooverChainThermostat(System& system, double temperature, double collisionFrequency,
int chainLength, int numMTS, int numYoshidaSuzuki);
int addThermostat(System& system, double temperature, double collisionFrequency,
int chainLength, int numMTS, int numYoshidaSuzuki);
/**
* Add a Nose-Hoover Chain thermostat to control the temperature of the system
* Add a Nose-Hoover Chain thermostat to control the temperature of a collection of atoms and/or pairs of
* connected atoms within the full system. A list of atoms defining the atoms to be thermostated is
* provided and the thermostat will only control members of that list. Additionally a list of pairs of
* connected atoms may be provided; in this case both the center of mass absolute motion of each pair is
* controlled as well as their motion relative to each other, which is independently thermostated.
*
* @param system the system to be thermostated. Note: this must be setup, i.e. all
* particles should have been added, before calling this function.
* @param mask list of particle ids to be thermostated.
* @param parents either an empty list of a list describing the parent atoms that each thermostated
* atom is connected to.
* @param temperature the target temperature for the system.
* @param collisionFrequency the frequency of the interaction with the heat bath (in 1/ps).
* @param thermostatedParticles list of particle ids to be thermostated.
* @param thermostatedPairs a list of pairs of connected atoms whose absolute center of mass motion
* and motion relative to one another will be independently thermostated.
* @param temperature the target temperature for each pair's absolute of center of mass motion.
* @param relativeTemperature the target temperature for each pair's relative motion.
* @param collisionFrequency the frequency of the interaction with the heat bath for the
* pairs' center of mass motion (in 1/ps).
* @param relativeCollisionFrequency the frequency of the interaction with the heat bath for the
* pairs' relative motion (in 1/ps).
* @param chainLength the number of beads in the Nose-Hoover chain.
* @param numMTS the number of step in the multiple time step chain propagation algorithm.
* @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
* used in the chain propagation algorithm (must be 1, 3, or 5).
*/
int addMaskedNoseHooverChainThermostat(System& system, const std::vector<int>& mask, const std::vector<int>& parents,
double temperature, double collisionFrequency,
int chainLength, int numMTS, int numYoshidaSuzuki);
int addSubsystemThermostat(System& system, const std::vector<int>& thermostatedParticles,
const std::vector< std::pair< int, int> >& thermostatedPairs,
double temperature, double relativeTemperature,
double collisionFrequency, double relativeCollisionFrequency,
int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 3);
/**
* Get the temperature of the i-th chain (in Kelvin).
* Get the temperature of the i-th chain for controling absolute particle motion (in Kelvin).
*
* @param chainID the index of the Nose-Hoover chain (default=0).
*
......@@ -104,14 +130,29 @@ public:
*/
double getTemperature(int chainID=0) const;
/**
* set the temperature of the i-th chain.
* set the (absolute motion) temperature of the i-th chain.
*
* @param temperature the temperature for the Nose-Hoover chain thermostat (in Kelvin).
* @param chainID The id of the Nose-Hoover chain for which the temperature is set (default=0).
*/
void setTemperature(double temperature, int chainID=0);
/**
* Get the collision frequency of the i-th chain (in 1/picosecond).
* Get the temperature of the i-th chain for controling pairs' relative particle motion (in Kelvin).
*
* @param chainID the index of the Nose-Hoover chain (default=0).
*
* @return the temperature.
*/
double getRelativeTemperature(int chainID=0) const;
/**
* set the (relative pair motion) temperature of the i-th chain.
*
* @param temperature the temperature for the Nose-Hoover chain thermostat (in Kelvin).
* @param chainID The id of the Nose-Hoover chain for which the temperature is set (default=0).
*/
void setRelativeTemperature(double temperature, int chainID=0);
/**
* Get the collision frequency for absolute motion of the i-th chain (in 1/picosecond).
*
* @param chainID the index of the Nose-Hoover chain (default=0).
*
......@@ -119,12 +160,27 @@ public:
*/
double getCollisionFrequency(int chainID=0) const;
/**
* Set the collision frequency of the i-th chain.
* Set the collision frequency for absolute motion of the i-th chain.
*
* @param frequency the collision frequency in picosecond.
* @param chainID the index of the Nose-Hoover chain (default=0).
*/
void setCollisionFrequency(double frequency, int chainID=0);
/**
* Get the collision frequency for pairs' relative motion of the i-th chain (in 1/picosecond).
*
* @param chainID the index of the Nose-Hoover chain (default=0).
*
* @return the collision frequency.
*/
double getRelativeCollisionFrequency(int chainID=0) const;
/**
* Set the collision frequency for pairs' relative motion of the i-th chain.
*
* @param frequency the collision frequency in picosecond.
* @param chainID the index of the Nose-Hoover chain (default=0).
*/
void setRelativeCollisionFrequency(double frequency, int chainID=0);
/**
* Compute the total (potential + kinetic) heat bath energy for all heat baths
* associated with this integrator, at the current time.
......@@ -154,11 +210,11 @@ protected:
* Advance any Nose-Hoover chains associated with this integrator and determine
* scale factor for the velocities.
*
* @param kineticEnergy the kinetic energy of the system that the chain is thermostating
* @param kineticEnergy the {absolute, relative} kinetic energies of the system that the chain is thermostating
* @param chainID id of the Nose-Hoover-Chain
* @return the scale factor to be applied to the velocities of the particles thermostated by the chain.
*/
double propagateChain(double kineticEnergy, int chainID=0);
std::pair<double, double> propagateChain(std::pair<double, double> kineticEnergy, int chainID=0);
/**
* This will be called by the Context when it is created. It informs the Integrator
* of what context it will be integrating, and gives it a chance to do any necessary initialization.
......@@ -186,4 +242,4 @@ protected:
} // namespace OpenMM
#endif /*OPENMM_VELOCITYVERLETINTEGRATOR_H_*/
#endif /*OPENMM_NOSEHOOVERINTEGRATOR_H_*/
openmmapi/src/NoseHooverChain.cpp
View file @ 50efef13
......@@ -34,13 +34,15 @@
using namespace OpenMM;
NoseHooverChain::NoseHooverChain(double defaultTemperature, double defaultCollisionFrequency,
NoseHooverChain::NoseHooverChain(double defaultTemperature, double defaultRelativeTemperature, double defaultCollisionFrequency,
double defaultRelativeCollisionFrequency,
int defaultNumDOFs, int defaultChainLength, int defaultNumMTS,
int defaultNumYoshidaSuzuki, int defaultChainID,
const std::vector<int>& thermostatedAtoms, const std::vector<int>& parentAtoms):
defaultTemp(defaultTemperature), defaultFreq(defaultCollisionFrequency), defaultNumDOFs(defaultNumDOFs),
const std::vector<int>& thermostatedAtoms, const std::vector<std::pair<int, int> > &thermostatedPairs):
defaultTemp(defaultTemperature), defaultRelativeTemp(defaultRelativeTemperature), defaultFreq(defaultCollisionFrequency),
defaultRelativeFreq(defaultRelativeCollisionFrequency), defaultNumDOFs(defaultNumDOFs),
defaultChainLength(defaultChainLength), defaultNumMTS(defaultNumMTS), defaultNumYS(defaultNumYoshidaSuzuki),
defaultChainID(defaultChainID), thermostatedAtoms(thermostatedAtoms), parentAtoms(parentAtoms)
defaultChainID(defaultChainID), thermostatedAtoms(thermostatedAtoms), thermostatedPairs(thermostatedPairs)
{}
......
openmmapi/src/VelocityVerletIntegrator.cpp → openmmapi/src/NoseHooverIntegrator.cpp
View file @ 50efef13
......@@ -29,7 +29,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/VelocityVerletIntegrator.h"
#include "openmm/NoseHooverIntegrator.h"
#include "openmm/Context.h"
#include "openmm/Force.h"
#include "openmm/System.h"
......@@ -46,24 +46,28 @@ using namespace OpenMM;
using std::string;
using std::vector;
VelocityVerletIntegrator::VelocityVerletIntegrator(double stepSize):
NoseHooverIntegrator::NoseHooverIntegrator(double stepSize):
forcesAreValid(false)
{
setStepSize(stepSize);
setConstraintTolerance(1e-5);
}
NoseHooverIntegrator::NoseHooverIntegrator(double stepSize, System &system, double temperature, int collisionFrequnency,
int chainLength, int numMTS, int numYoshidaSuzuki) : forcesAreValid(false) {
setStepSize(stepSize);
setConstraintTolerance(1e-5);
addThermostat(system, temperature, collisionFrequnency, chainLength, numMTS, numYoshidaSuzuki);
}
VelocityVerletIntegrator::~VelocityVerletIntegrator() {}
NoseHooverIntegrator::~NoseHooverIntegrator() {}
double VelocityVerletIntegrator::propagateChain(double kineticEnergy, int chainID) {
std::pair<double, double> NoseHooverIntegrator::propagateChain(std::pair<double, double> kineticEnergy, int chainID) {
return nhcKernel.getAs<NoseHooverChainKernel>().propagateChain(*context, noseHooverChains.at(chainID), kineticEnergy, getStepSize());
}
int VelocityVerletIntegrator::addNoseHooverChainThermostat(System& system, double temperature, double collisionFrequency,
int NoseHooverIntegrator::addThermostat(System& system, double temperature, double collisionFrequency,
int chainLength, int numMTS, int numYoshidaSuzuki) {
int numForces = system.getNumForces();
std::vector<int> thermostatedParticles;
std::vector<int> parentParticles;
for(int particle = 0; particle < system.getNumParticles(); ++particle) {
double mass = system.getParticleMass(particle);
if ( (mass > 0) && (mass < 0.8) ){
......@@ -75,51 +79,62 @@ int VelocityVerletIntegrator::addNoseHooverChainThermostat(System& system, doubl
}
}
return addMaskedNoseHooverChainThermostat(system, thermostatedParticles, parentParticles, temperature, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
return addSubsystemThermostat(system, thermostatedParticles, std::vector<std::pair<int, int>>(), temperature, temperature,
collisionFrequency, collisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
}
int VelocityVerletIntegrator::addMaskedNoseHooverChainThermostat(System& system, const std::vector<int>& thermostatedParticles, const std::vector<int>& parentParticles,
double temperature, double collisionFrequency,
int chainLength, int numMTS, int numYoshidaSuzuki){
const auto & mask = thermostatedParticles; // just an alias
int NoseHooverIntegrator::addSubsystemThermostat(System& system, const std::vector<int>& thermostatedParticles,
const std::vector< std::pair< int, int> > &thermostatedPairs,
double temperature, double relativeTemperature,
double collisionFrequency, double relativeCollisionFrequency,
int chainLength, int numMTS, int numYoshidaSuzuki) {
if (context) {
throw OpenMMException("Nose-Hoover chains cannot be added after binding this integrator to a context.");
}
if ((parentParticles.size() != mask.size()) && (parentParticles.size() != 0)) {
throw OpenMMException("The number of parent particles has to be either equal to the number of thermostated particles (to thermostat relative motion)"
" or zero (to thermostat absolute motion).");
}
// figure out the number of DOFs
int nDOF = 3*mask.size();
int nDOF = 3*(thermostatedParticles.size() + thermostatedPairs.size());
for (int constraintNum = 0; constraintNum < system.getNumConstraints(); constraintNum++) {
int particle1, particle2;
double distance;
system.getConstraintParameters(constraintNum, particle1, particle2, distance);
bool particle1_in_mask = (std::find(mask.begin(), mask.end(), particle1) != mask.end());
bool particle2_in_mask = (std::find(mask.begin(), mask.end(), particle2) != mask.end());
bool particle1_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(), thermostatedParticles.end(), particle1)
!= thermostatedParticles.end())) ||
(std::find_if(thermostatedPairs.begin(), thermostatedPairs.end(),
[&particle1](const std::pair<int, int>& pair){ return pair.first == particle1 || pair.second == particle1;})
!= thermostatedPairs.end());
bool particle2_in_thermostatedParticles = ((std::find(thermostatedParticles.begin(), thermostatedParticles.end(), particle2)
!= thermostatedParticles.end())) ||
(std::find_if(thermostatedPairs.begin(), thermostatedPairs.end(),
[&particle2](const std::pair<int, int>& pair){ return pair.first == particle2 || pair.second == particle2;})
!= thermostatedPairs.end());
if ((system.getParticleMass(particle1) > 0) && (system.getParticleMass(particle2) > 0)){
if ((particle1_in_mask && !particle2_in_mask) || (!particle1_in_mask && particle2_in_mask)){
if ((particle1_in_thermostatedParticles && !particle2_in_thermostatedParticles) ||
(!particle1_in_thermostatedParticles && particle2_in_thermostatedParticles)){
throw OpenMMException("Cannot add only one of particles " + std::to_string(particle1) + " and " + std::to_string(particle2)
+ " to NoseHooverChain, because they are connected by a constraint.");
}
if (particle1_in_mask && particle2_in_mask){
if (particle1_in_thermostatedParticles && particle2_in_thermostatedParticles){
nDOF -= 1;
}
}
}
int numForces = system.getNumForces();
if (parentParticles.size() == 0){ // remove 3 degrees of freedom from thermostats that act on absolute motions
if (thermostatedPairs.size() == 0){ // remove 3 degrees of freedom from thermostats that act on absolute motions
for (int forceNum = 0; forceNum < numForces; ++forceNum) {
if (dynamic_cast<CMMotionRemover*>(&system.getForce(forceNum))) nDOF -= 3;
// if (dynamic_cast<CMMotionRemover*>(&system.getForce(forceNum))) nDOF -= 3;
}
}
// make sure that thermostats do not overlap
// make sure that thermostats do not overlap
for (auto &other_nhc: noseHooverChains) {
for (auto &particle: mask){
bool isParticleInOtherChain = (std::find(other_nhc.getThermostatedAtoms().begin(),
other_nhc.getThermostatedAtoms().end(),
particle) == other_nhc.getThermostatedAtoms().begin());
for (auto &particle: thermostatedParticles){
bool isParticleInOtherChain = (std::find(other_nhc.getThermostatedAtoms().begin(),
other_nhc.getThermostatedAtoms().end(),
particle) != other_nhc.getThermostatedAtoms().end()) ||
(std::find_if(other_nhc.getThermostatedPairs().begin(),
other_nhc.getThermostatedPairs().end(),
[&particle](const std::pair<int, int>& pair){ return pair.first == particle || pair.second == particle;})
!= other_nhc.getThermostatedPairs().end());
if (isParticleInOtherChain){
throw OpenMMException("Found particle " + std::to_string(particle) + "in a different NoseHooverChain, "
"but particles can only be thermostated by one thermostat.");
......@@ -127,16 +142,19 @@ int VelocityVerletIntegrator::addMaskedNoseHooverChainThermostat(System& system,
}
}
// create and add new chain
int chainID = noseHooverChains.size();
auto nhc = NoseHooverChain(temperature, collisionFrequency, nDOF, chainLength,
numMTS, numYoshidaSuzuki, chainID,
thermostatedParticles, parentParticles);
// create and add new chain
int chainID = 0;
for(const auto &c : noseHooverChains) {
// We need to account for the fact that a NHC with pairs thermostated has both a relative and an absolute chain
chainID += c.getThermostatedPairs().size() ? 2 : 1;
}
auto nhc = NoseHooverChain(temperature, relativeTemperature, collisionFrequency, relativeCollisionFrequency,
nDOF, chainLength, numMTS, numYoshidaSuzuki, chainID, thermostatedParticles, thermostatedPairs);
noseHooverChains.push_back(nhc);
return chainID;
}
double VelocityVerletIntegrator::getTemperature(int chainID) const {
double NoseHooverIntegrator::getTemperature(int chainID) const {
if (chainID >= noseHooverChains.size()) {
throw OpenMMException("Cannot get temperature for chainID " + std::to_string(chainID)
+ ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
......@@ -145,7 +163,7 @@ double VelocityVerletIntegrator::getTemperature(int chainID) const {
return noseHooverChains.at(chainID).getDefaultTemperature();
}
void VelocityVerletIntegrator::setTemperature(double temperature, int chainID){
void NoseHooverIntegrator::setTemperature(double temperature, int chainID){
if (chainID >= noseHooverChains.size()) {
throw OpenMMException("Cannot set temperature for chainID " + std::to_string(chainID)
+ ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
......@@ -155,15 +173,35 @@ void VelocityVerletIntegrator::setTemperature(double temperature, int chainID){
}
double VelocityVerletIntegrator::getCollisionFrequency(int chainID) const {
double NoseHooverIntegrator::getRelativeTemperature(int chainID) const {
if (chainID >= noseHooverChains.size()) {
throw OpenMMException("Cannot get relative temperature for chainID " + std::to_string(chainID)
+ ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
);
}
return noseHooverChains.at(chainID).getDefaultRelativeTemperature();
}
void NoseHooverIntegrator::setRelativeTemperature(double temperature, int chainID){
if (chainID >= noseHooverChains.size()) {
throw OpenMMException("Cannot set relative temperature for chainID " + std::to_string(chainID)
+ ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
);
}
noseHooverChains.at(chainID).setDefaultRelativeTemperature(temperature);
}
double NoseHooverIntegrator::getCollisionFrequency(int chainID) const {
if (chainID >= noseHooverChains.size()) {
throw OpenMMException("Cannot get collision frequency for chainID " + std::to_string(chainID)
+ ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
);
}
return noseHooverChains.at(chainID).getDefaultCollisionFrequency();
}
void VelocityVerletIntegrator::setCollisionFrequency(double frequency, int chainID){
}
void NoseHooverIntegrator::setCollisionFrequency(double frequency, int chainID){
if (chainID >= noseHooverChains.size()) {
throw OpenMMException("Cannot set collision frequency for chainID " + std::to_string(chainID)
+ ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
......@@ -172,12 +210,31 @@ void VelocityVerletIntegrator::setCollisionFrequency(double frequency, int chain
noseHooverChains.at(chainID).setDefaultCollisionFrequency(frequency);
}
double VelocityVerletIntegrator::computeKineticEnergy() {
double NoseHooverIntegrator::getRelativeCollisionFrequency(int chainID) const {
if (chainID >= noseHooverChains.size()) {
throw OpenMMException("Cannot get relative collision frequency for chainID " + std::to_string(chainID)
+ ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
);
}
return noseHooverChains.at(chainID).getDefaultRelativeCollisionFrequency();
}
void NoseHooverIntegrator::setRelativeCollisionFrequency(double frequency, int chainID){
if (chainID >= noseHooverChains.size()) {
throw OpenMMException("Cannot set relative collision frequency for chainID " + std::to_string(chainID)
+ ". Only " + std::to_string(noseHooverChains.size()) + " have been registered with this integrator."
);
}
noseHooverChains.at(chainID).setDefaultRelativeCollisionFrequency(frequency);
}
double NoseHooverIntegrator::computeKineticEnergy() {
double kE = 0.0;
if(noseHooverChains.size()) {
for (const auto &nhc: noseHooverChains){
if (nhc.getParentAtoms().size() == 0) {
kE += nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, true);
if (nhc.getThermostatedPairs().size() == 0) {
auto KEs = nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, true);
kE += std::get<0>(KEs);
}
}
} else {
......@@ -186,7 +243,7 @@ double VelocityVerletIntegrator::computeKineticEnergy() {
return kE;
}
double VelocityVerletIntegrator::computeHeatBathEnergy() {
double NoseHooverIntegrator::computeHeatBathEnergy() {
double energy = 0;
for(auto &nhc : noseHooverChains) {
energy += nhcKernel.getAs<NoseHooverChainKernel>().computeHeatBathEnergy(*context, nhc);
......@@ -194,7 +251,7 @@ double VelocityVerletIntegrator::computeHeatBathEnergy() {
return energy;
}
void VelocityVerletIntegrator::initialize(ContextImpl& contextRef) {
void NoseHooverIntegrator::initialize(ContextImpl& contextRef) {
if (owner != NULL && &contextRef.getOwner() != owner)
throw OpenMMException("This Integrator is already bound to a context");
context = &contextRef;
......@@ -206,22 +263,22 @@ void VelocityVerletIntegrator::initialize(ContextImpl& contextRef) {
forcesAreValid = false;
}
void VelocityVerletIntegrator::cleanup() {
void NoseHooverIntegrator::cleanup() {
vvKernel = Kernel();
nhcKernel = Kernel();
}
vector<string> VelocityVerletIntegrator::getKernelNames() {
vector<string> NoseHooverIntegrator::getKernelNames() {
std::vector<std::string> names;
names.push_back(NoseHooverChainKernel::Name());
names.push_back(IntegrateVelocityVerletStepKernel::Name());
return names;
}
void VelocityVerletIntegrator::step(int steps) {
void NoseHooverIntegrator::step(int steps) {
if (context == NULL)
throw OpenMMException("This Integrator is not bound to a context!");
double scale, kineticEnergy;
std::pair<double, double> scale, kineticEnergy;
for (int i = 0; i < steps; ++i) {
context->updateContextState();
for(auto &nhc : noseHooverChains) {
......
platforms/cuda/include/CudaKernels.h
View file @ 50efef13
......@@ -1380,7 +1380,7 @@ public:
* @param system the System this kernel will be applied to
* @param integrator the VelocityVerletIntegrator this kernel will be used for
*/
void initialize(const System& system, const VelocityVerletIntegrator& integrator);
void initialize(const System& system, const NoseHooverIntegrator& integrator);
/**
* Execute the kernel.
*
......@@ -1389,14 +1389,14 @@ public:
* @param forcesAreValid a reference to the parent integrator's boolean for keeping
* track of the validity of the current forces.
*/
void execute(ContextImpl& context, const VelocityVerletIntegrator& integrator, bool &forcesAreValid);
void execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid);
/**
* Compute the kinetic energy.
*
* @param context the context in which to execute this kernel
* @param integrator the VelocityVerletIntegrator this kernel is being used for
*/
double computeKineticEnergy(ContextImpl& context, const VelocityVerletIntegrator& integrator);
double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator);
private:
CudaContext& cu;
CUfunction kernel1, kernel2, kernel3;
......
platforms/cuda/src/CudaKernels.cpp
View file @ 50efef13
......@@ -7076,7 +7076,7 @@ double CudaIntegrateVerletStepKernel::computeKineticEnergy(ContextImpl& context,
return cu.getIntegrationUtilities().computeKineticEnergy(0.5*integrator.getStepSize());
}
void CudaIntegrateVelocityVerletStepKernel::initialize(const System& system, const VelocityVerletIntegrator& integrator) {
void CudaIntegrateVelocityVerletStepKernel::initialize(const System& system, const NoseHooverIntegrator& integrator) {
cu.getPlatformData().initializeContexts(system);
cu.setAsCurrent();
map<string, string> defines;
......@@ -7086,7 +7086,7 @@ void CudaIntegrateVelocityVerletStepKernel::initialize(const System& system, con
kernel3 = cu.getKernel(module, "integrateVelocityVerletPart3");
}
void CudaIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const VelocityVerletIntegrator& integrator, bool &forcesAreValid) {
void CudaIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
cu.setAsCurrent();
CudaIntegrationUtilities& integration = cu.getIntegrationUtilities();
int numAtoms = cu.getNumAtoms();
......@@ -7136,7 +7136,7 @@ void CudaIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const
cu.reorderAtoms();
}
double CudaIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VelocityVerletIntegrator& integrator) {
double CudaIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) {
return cu.getIntegrationUtilities().computeKineticEnergy(0);
}
......@@ -8512,7 +8512,7 @@ double CudaNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& contex
bool useDouble = cu.getUseDoublePrecision() || cu.getUseMixedPrecision();
int chainID = nhc.getDefaultChainID();
const auto & parents = nhc.getParentAtoms();
const auto & parents = nhc.getThermostatedPairs();
const auto & thermostatedAtoms = nhc.getThermostatedAtoms();
auto nAtoms = cu.getPaddedNumAtoms();
if (!masks.count(chainID)) {
......
platforms/cuda/tests/TestCudaVelocityVerletIntegrator.cpp
View file @ 50efef13
......@@ -30,7 +30,7 @@
* -------------------------------------------------------------------------- */
#include "CudaTests.h"
#include "TestVelocityVerletIntegrator.h"
#include "TestNoseHooverIntegrator.h"
void runPlatformTests() {
}
platforms/reference/include/ReferenceKernels.h
View file @ 50efef13
......@@ -1150,7 +1150,7 @@ public:
* @param system the System this kernel will be applied to
* @param integrator the VelocityVerletIntegrator this kernel will be used for
*/
void initialize(const System& system, const VelocityVerletIntegrator& integrator);
void initialize(const System& system, const NoseHooverIntegrator& integrator);
/**
* Execute the kernel.
*
......@@ -1159,14 +1159,14 @@ public:
* @param forcesAreValid a reference to the parent integrator's boolean for keeping
* track of the validity of the current forces.
*/
void execute(ContextImpl& context, const VelocityVerletIntegrator& integrator, bool &forcesAreValid);
void execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid);
/**
* Compute the kinetic energy.
*
* @param context the context in which to execute this kernel
* @param integrator the VelocityVerletIntegrator this kernel is being used for
*/
double computeKineticEnergy(ContextImpl& context, const VelocityVerletIntegrator& integrator);
double computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator);
private:
ReferencePlatform::PlatformData& data;
ReferenceVelocityVerletDynamics* dynamics;
......@@ -1446,11 +1446,11 @@ public:
*
* @param context the context in which to execute this kernel
* @param noseHooverChain the object describing the chain to be propagated.
* @param kineticEnergy the kineticEnergy of the particles being thermostated by this chain.
* @param kineticEnergy the {absolute, relative} kinetic energies of the particles being thermostated by this chain.
* @param timeStep the time step used by the integrator.
* @return the velocity scale factor to apply to the particles associated with this heat bath.
* @return the velocity scale factors to apply to the {absolute, relative} motion of particles associated with this heat bath.
*/
virtual double propagateChain(ContextImpl& context, const NoseHooverChain &nhc, double kineticEnergy, double timeStep);
virtual std::pair<double, double> propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergy, double timeStep);
/**
* Execute the kernal that computes the total (kinetic + potential) heat bath energy.
*
......@@ -1468,16 +1468,16 @@ public:
* @param downloadValue whether the computed value should be downloaded and returned.
*
*/
virtual double computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
virtual std::pair<double, double> computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue);
/**
* Execute the kernel that scales the velocities of particles associated with a nose hoover chain
*
* @param context the context in which to execute this kernel
* @param noseHooverChain the chain whose energy is to be determined.
* @param scaleFactor the multiplicative factor by which velocities are scaled.
* @param scaleFactor the multiplicative factor by which {absolute, relative} velocities are scaled.
*/
virtual void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, double scaleFactor);
virtual void scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactor);
private:
......
platforms/reference/src/ReferenceKernels.cpp
View file @ 50efef13
......@@ -2150,14 +2150,14 @@ ReferenceIntegrateVelocityVerletStepKernel::~ReferenceIntegrateVelocityVerletSte
delete dynamics;
}
void ReferenceIntegrateVelocityVerletStepKernel::initialize(const System& system, const VelocityVerletIntegrator& integrator) {
void ReferenceIntegrateVelocityVerletStepKernel::initialize(const System& system, const NoseHooverIntegrator& integrator) {
int numParticles = system.getNumParticles();
masses.resize(numParticles);
for (int i = 0; i < numParticles; ++i)
masses[i] = system.getParticleMass(i);
}
void ReferenceIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const VelocityVerletIntegrator& integrator, bool &forcesAreValid) {
void ReferenceIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, const NoseHooverIntegrator& integrator, bool &forcesAreValid) {
double stepSize = integrator.getStepSize();
vector<Vec3>& posData = extractPositions(context);
vector<Vec3>& velData = extractVelocities(context);
......@@ -2175,7 +2175,7 @@ void ReferenceIntegrateVelocityVerletStepKernel::execute(ContextImpl& context, c
data.stepCount++;
}
double ReferenceIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const VelocityVerletIntegrator& integrator) {
double ReferenceIntegrateVelocityVerletStepKernel::computeKineticEnergy(ContextImpl& context, const NoseHooverIntegrator& integrator) {
return computeShiftedKineticEnergy(context, masses, 0);
}
......@@ -2470,8 +2470,10 @@ void ReferenceNoseHooverChainKernel::initialize() {
//SimTKOpenMMUtilities::setRandomNumberSeed((unsigned int) thermostat.getRandomNumberSeed());
}
double ReferenceNoseHooverChainKernel::propagateChain(ContextImpl& context, const NoseHooverChain &nhc, double kineticEnergy, double timeStep) {
if (kineticEnergy < 1e-8) return 1.0; // (catches the problem of zero velocities in the first dynamics step, where we have nothing to scale)
std::pair<double, double> ReferenceNoseHooverChainKernel::propagateChain(ContextImpl& context, const NoseHooverChain &nhc, std::pair<double, double> kineticEnergy, double timeStep) {
double absKE = kineticEnergy.first;
double relKE = kineticEnergy.second;
if (absKE < 1e-8) return {1.0, 1.0}; // (catches the problem of zero velocities in the first dynamics step, where we have nothing to scale)
// Get the variables describing the NHC
int chainLength = nhc.getDefaultChainLength();
int chainID = nhc.getDefaultChainID();
......@@ -2497,10 +2499,34 @@ double ReferenceNoseHooverChainKernel::propagateChain(ContextImpl& context, cons
if (chainVelocities.size() < chainLength){
chainVelocities.resize(chainLength, 0);
}
double scale = chainPropagator->propagate(kineticEnergy, chainVelocities, chainPositions, numDOFs,
temperature, collisionFrequency, timeStep,
double absScale = chainPropagator->propagate(absKE, chainVelocities, chainPositions, numDOFs,
temperature, collisionFrequency, timeStep,
numMTS, nhc.getDefaultYoshidaSuzukiWeights());
double relScale = 0;
int numPairs = nhc.getThermostatedPairs().size();
if(numPairs){
double relativeTemperature = nhc.getDefaultRelativeTemperature();
double relativeCollisionFrequency = nhc.getDefaultRelativeCollisionFrequency();
if (allChainPositions.size() < chainID+2){
allChainPositions.resize(chainID+2);
}
if (allChainVelocities.size() < chainID+2){
allChainVelocities.resize(chainID+2);
}
auto& chainPositions = allChainPositions.at(chainID+1);
auto& chainVelocities = allChainVelocities.at(chainID+1);
if (chainPositions.size() < chainLength){
chainPositions.resize(chainLength, 0);
}
if (chainVelocities.size() < chainLength){
chainVelocities.resize(chainLength, 0);
}
relScale = chainPropagator->propagate(relKE, chainVelocities, chainPositions, 3*numPairs,
relativeTemperature, relativeCollisionFrequency, timeStep,
numMTS, nhc.getDefaultYoshidaSuzukiWeights());
return scale;
}
return {absScale, relScale};
}
double ReferenceNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& context, const NoseHooverChain &nhc) {
......@@ -2527,9 +2553,9 @@ double ReferenceNoseHooverChainKernel::computeHeatBathEnergy(ContextImpl& contex
return kineticEnergy + potentialEnergy;
}
double ReferenceNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) {
const std::vector<int>& mask = noseHooverChain.getThermostatedAtoms();
const std::vector<int>& parents = noseHooverChain.getParentAtoms();
std::pair<double, double> ReferenceNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& context, const NoseHooverChain &noseHooverChain, bool downloadValue) {
const std::vector<int>& atomsList = noseHooverChain.getThermostatedAtoms();
const std::vector<std::pair<int,int>>& pairsList = noseHooverChain.getThermostatedPairs();
std::vector<Vec3>& velocities = extractVelocities(context);
const System& system = context.getSystem();
int numParticles = system.getNumParticles();
......@@ -2537,44 +2563,64 @@ double ReferenceNoseHooverChainKernel::computeMaskedKineticEnergy(ContextImpl& c
for (int i = 0; i < numParticles; ++i)
masses[i] = system.getParticleMass(i);
double ke = 0;
if (parents.size() == 0){
// scale absolute velocities
for (auto m: mask){
ke += 0.5 * masses[m] * velocities[m].dot(velocities[m]);
}
} else {
// scale velocities relative to parent
assert(parents.size() == mask.size());
Vec3 velocity;
double mass;
for (int i=0; i < mask.size(); i++){
velocity = velocities[mask[i]] - velocities[parents[i]];
mass = masses[mask[i]] * masses[parents[i]] / (masses[mask[i]] + masses[parents[i]]);
ke += 0.5 * mass * velocity.dot(velocity);
}
double comKE = 0;
double relKE = 0;
// kinetic energy of individual atoms
for (const auto &m: atomsList){
comKE += 0.5 * masses[m] * velocities[m].dot(velocities[m]);
}
// center of mass kinetic energy of pairs
for (const auto &p: pairsList){
double m1 = masses[p.first];
double m2 = masses[p.second];
Vec3 v1 = velocities[p.first];
Vec3 v2 = velocities[p.second];
double invMass = 1.0 / (m1 + m2);
double redMass = m1 * m2 * invMass;
double fracM1 = m1 * invMass;
double fracM2 = m2 * invMass;
Vec3 comVelocity = fracM1 * v1 + fracM2 * v2;
Vec3 relVelocity = v2 - v1;
comKE += 0.5 * (m1 + m2) * comVelocity.dot(comVelocity);
relKE += 0.5 * redMass * relVelocity.dot(relVelocity);
}
// We ignore the downloadValue argument here and always return the correct value
return ke;
return {comKE, relKE};
}
void ReferenceNoseHooverChainKernel::scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, double scaleFactor) {
const std::vector<int>& mask = noseHooverChain.getThermostatedAtoms();
const std::vector<int>& parents = noseHooverChain.getParentAtoms();
void ReferenceNoseHooverChainKernel::scaleVelocities(ContextImpl& context, const NoseHooverChain &noseHooverChain, std::pair<double, double> scaleFactors) {
const auto& atoms = noseHooverChain.getThermostatedAtoms();
const auto& pairs = noseHooverChain.getThermostatedPairs();
std::vector<Vec3>& velocities = extractVelocities(context);
double absScale = scaleFactors.first;
double relScale = scaleFactors.second;
if (parents.size() == 0){
// scale absolute velocities
for (auto m: mask){
velocities[m] *= scaleFactor;
}
} else {
// scale velocities relative to parent
assert(parents.size() == mask.size());
for (int i=0; i < mask.size(); i++){
velocities[mask[i]] = velocities[parents[i]] + scaleFactor * (velocities[mask[i]] - velocities[parents[i]]);
}
const System& system = context.getSystem();
int numParticles = system.getNumParticles();
std::vector<double> masses(numParticles);
for (int i = 0; i < numParticles; ++i)
masses[i] = system.getParticleMass(i);
// scale absolute velocities
for (const auto &a: atoms){
velocities[a] *= absScale;
}
// scale relative velocities and absolute center of mass velocities for each pair
for (const auto &p: pairs){
int p1 = p.first;
int p2 = p.second;
double m1 = masses[p.first];
double m2 = masses[p.second];
Vec3 v1 = velocities[p.first];
Vec3 v2 = velocities[p.second];
double invMass = 1.0 / (m1 + m2);
double fracM1 = m1 * invMass;
double fracM2 = m2 * invMass;
Vec3 comVelocity = fracM1 * v1 + fracM2 * v2;
Vec3 relVelocity = v2 - v1;
velocities[p1] = absScale * comVelocity - relScale * relVelocity * fracM2;
velocities[p2] = absScale * comVelocity + relScale * relVelocity * fracM1;
}
}
......
platforms/reference/tests/TestReferenceVelocityVerletIntegrator.cpp
View file @ 50efef13
......@@ -30,7 +30,7 @@
* -------------------------------------------------------------------------- */
#include "ReferenceTests.h"
#include "TestVelocityVerletIntegrator.h"
#include "TestNoseHooverIntegrator.h"
void runPlatformTests() {
}
plugins/drude/openmmapi/include/OpenMMDrude.h
View file @ 50efef13
......@@ -35,6 +35,6 @@
#include "openmm/DrudeForce.h"
#include "openmm/DrudeLangevinIntegrator.h"
#include "openmm/DrudeSCFIntegrator.h"
#include "openmm/DrudeVelocityVerletIntegrator.h"
#include "openmm/DrudeNIntegrator.h"
#endif /*OPENMM_DRUDE_H_*/
plugins/drude/openmmapi/include/openmm/DrudeVelocityVerletIntegrator.h → plugins/drude/openmmapi/include/openmm/DrudeNoseHooverIntegrator.h
View file @ 50efef13
#ifndef OPENMM_DRUDEVELOCITYVERLETINTEGRATOR_H_
#define OPENMM_DRUDEVELOCITYVERLETINTEGRATOR_H_
#ifndef OPENMM_DRUDENOSEHOOVERINTEGRATOR_H_
#define OPENMM_DRUDENOSEHOOVERINTEGRATOR_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
......@@ -32,7 +32,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/VelocityVerletIntegrator.h"
#include "openmm/NoseHooverIntegrator.h"
#include "openmm/Kernel.h"
#include "openmm/internal/windowsExportDrude.h"
......@@ -49,33 +49,28 @@ namespace OpenMM {
* particles.
*/
class OPENMM_EXPORT_DRUDE DrudeVelocityVerletIntegrator : public VelocityVerletIntegrator {
class OPENMM_EXPORT_DRUDE DrudeNoseHooverIntegrator : public NoseHooverIntegrator {
public:
/**
* Create a DrudeVelocityVerletIntegrator.
* Create a DrudeNoseHooverIntegrator.
*
* @param stepSize the step size with which to integrator the system (in picoseconds)
*/
DrudeVelocityVerletIntegrator(double stepSize);
virtual ~DrudeVelocityVerletIntegrator();
/**
* Add a Nose-Hoover Chain thermostat to control the temperature of the system
*
* @param system the system to be thermostated. Note: this must be setup, i.e. all
* particles should have been added, before calling this function.
* @param temperature the target temperature for the system.
* @param collisionFrequency the frequency of the interaction with the heat bath (in 1/ps).
* @param temperature the target temperature for the drude particles..
* @param collisionFrequency the frequency of the interaction of the drude particles with the heat bath (in 1/ps).
* @param drudeTemperature the target temperature for the Drude particles, relative to their parent atom.
* @param collisionFrequency the frequency of the system's interaction with the heat bath (in 1/ps).
* @param drudeCollisionFrequency the frequency of the drude particles' interaction with the heat bath (in 1/ps).
* @param chainLength the number of beads in the Nose-Hoover chain.
* @param numMTS the number of step in the multiple time step chain propagation algorithm.
* @param numYoshidaSuzuki the number of terms in the Yoshida-Suzuki multi time step decomposition
* used in the chain propagation algorithm (must be 1, 3, or 5).
*/
int addDrudeNoseHooverChainThermostat(System& system, double temperature, double collisionFrequency,
double drudeTemperature, double drudeCollisionFrequency,
int chainLength, int numMTS, int numYoshidaSuzuki);
DrudeNoseHooverIntegrator(double stepSize, System &system, double temperature, double drudeTemperature,
double collisionFrequency, double drudeCollisionFrequency,
int chainLength = 3, int numMTS = 3, int numYoshidaSuzuki = 3);
virtual ~DrudeNoseHooverIntegrator();
/**
* This will be called by the Context when it is created. It informs the Integrator
* of what context it will be integrating, and gives it a chance to do any necessary initialization.
......@@ -94,4 +89,4 @@ public:
} // namespace OpenMM
#endif /*OPENMM_DRUDEVELOCITYVERLETINTEGRATOR_H_*/
#endif /*OPENMM_DRUDENOSEHOOVERINTEGRATOR_H_*/
plugins/drude/openmmapi/src/DrudeVelocityVerletIntegrator.cpp → plugins/drude/openmmapi/src/DrudeNoseHooverIntegrator.cpp
View file @ 50efef13
......@@ -29,7 +29,7 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/DrudeVelocityVerletIntegrator.h"
#include "openmm/DrudeNoseHooverIntegrator.h"
#include "openmm/Context.h"
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
......@@ -46,13 +46,10 @@ using namespace OpenMM;
using std::string;
using std::vector;
DrudeVelocityVerletIntegrator::DrudeVelocityVerletIntegrator(double stepSize) : VelocityVerletIntegrator(stepSize) { }
DrudeVelocityVerletIntegrator::~DrudeVelocityVerletIntegrator() { }
int DrudeVelocityVerletIntegrator::addDrudeNoseHooverChainThermostat(System& system, double temperature, double collisionFrequency,
double drudeTemperature, double drudeCollisionFrequency,
int chainLength, int numMTS, int numYoshidaSuzuki) {
DrudeNoseHooverIntegrator::DrudeNoseHooverIntegrator(double stepSize, System &system, double temperature, double drudeTemperature,
double collisionFrequency, double drudeCollisionFrequency,
int chainLength, int numMTS, int numYoshidaSuzuki) :
NoseHooverIntegrator(stepSize) {
const DrudeForce* drudeForce = NULL;
for (int i = 0; i < system.getNumForces(); i++)
if (dynamic_cast<const DrudeForce*>(&system.getForce(i)) != NULL) {
......@@ -63,30 +60,31 @@ int DrudeVelocityVerletIntegrator::addDrudeNoseHooverChainThermostat(System& sys
}
if (drudeForce == NULL)
throw OpenMMException("The System does not contain a DrudeForce");
std::set<int> realParticlesSet;
std::set<int> realParticlesSet;
vector<int> realParticles, drudeParticles, drudeParents;
vector<std::pair<int, int>> drudePairs;
for (int i = 0; i < system.getNumParticles(); i++) {
if (system.getParticleMass(i) > 0.0) realParticlesSet.insert(i);
if (system.getParticleMass(i) > 0.0) realParticlesSet.insert(i);
}
for (int i = 0; i < drudeForce->getNumParticles(); i++) {
int p, p1, p2, p3, p4;
double charge, polarizability, aniso12, aniso34;
drudeForce->getParticleParameters(i, p, p1, p2, p3, p4, charge, polarizability, aniso12, aniso34);
realParticlesSet.erase(p);
drudeParticles.push_back(p);
drudeParents.push_back(p1);
realParticlesSet.erase(p1);
drudePairs.push_back({p,p1});
}
for(const auto &p : realParticlesSet) realParticles.push_back(p);
addMaskedNoseHooverChainThermostat(system, realParticles, vector<int>(), temperature, collisionFrequency,
chainLength, numMTS, numYoshidaSuzuki);
addMaskedNoseHooverChainThermostat(system, drudeParticles, drudeParents, drudeTemperature, drudeCollisionFrequency,
chainLength, numMTS, numYoshidaSuzuki);
return noseHooverChains.size() - 1;
addSubsystemThermostat(system, realParticles, drudePairs, temperature, drudeTemperature, collisionFrequency,
drudeCollisionFrequency, chainLength, numMTS, numYoshidaSuzuki);
}
void DrudeVelocityVerletIntegrator::initialize(ContextImpl& contextRef) {
VelocityVerletIntegrator::initialize(contextRef);
DrudeNoseHooverIntegrator::~DrudeNoseHooverIntegrator() { }
void DrudeNoseHooverIntegrator::initialize(ContextImpl& contextRef) {
NoseHooverIntegrator::initialize(contextRef);
if (owner != NULL && &contextRef.getOwner() != owner)
throw OpenMMException("This Integrator is already bound to a context");
const DrudeForce* drudeForce = NULL;
......@@ -113,17 +111,17 @@ void DrudeVelocityVerletIntegrator::initialize(ContextImpl& contextRef) {
owner = &contextRef.getOwner();
}
double DrudeVelocityVerletIntegrator::computeDrudeKineticEnergy() {
double DrudeNoseHooverIntegrator::computeDrudeKineticEnergy() {
double kE = 0.0;
for (const auto &nhc: noseHooverChains){
if (nhc.getParentAtoms().size() != 0) {
kE += nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, true);
if (nhc.getThermostatedPairs().size() != 0) {
kE += nhcKernel.getAs<NoseHooverChainKernel>().computeMaskedKineticEnergy(*context, nhc, true).second;
}
}
return kE;
}
double DrudeVelocityVerletIntegrator::computeTotalKineticEnergy() {
double DrudeNoseHooverIntegrator::computeTotalKineticEnergy() {
return vvKernel.getAs<IntegrateVelocityVerletStepKernel>().computeKineticEnergy(*context, *this);
}
plugins/drude/tests/TestDrudeNoseHoover.h 0 → 100644
View file @ 50efef13
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2019 Stanford University and the Authors. *
* Authors: Andreas Krämer and Andrew C. Simmonett *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/NoseHooverChain.h"
#include "openmm/CMMotionRemover.h"
#include "openmm/DrudeNoseHooverIntegrator.h"
#include "openmm/Context.h"
#include "openmm/State.h"
#include "openmm/HarmonicBondForce.h"
#include "openmm/VirtualSite.h"
#include "openmm/NonbondedForce.h"
#include "openmm/CustomExternalForce.h"
#include "openmm/System.h"
#include "openmm/DrudeForce.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <numeric>
using namespace OpenMM;
using namespace std;
//extern "C" OPENMM_EXPORT void registerDrudeReferenceKernelFactories();
extern "C" OPENMM_EXPORT void registerKernelFactories();
Platform& initializePlatform(int argc, char* argv[]);
void testWaterBox(Platform& platform){
// Create a box of SWM4-NDP water molecules. This involves constraints, virtual sites,
// and Drude particles.
const int gridSize = 3;
const int numMolecules = gridSize*gridSize*gridSize;
const double spacing = 0.8;
const double boxSize = spacing*(gridSize+1);
const double temperature = 300.0;
const double temperatureDrude = 10.0;
System system;
NonbondedForce* nonbonded = new NonbondedForce();
DrudeForce* drude = new DrudeForce();
CMMotionRemover* cmm = new CMMotionRemover(1);
system.addForce(cmm);
system.addForce(nonbonded);
system.addForce(drude);
system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
nonbonded->setCutoffDistance(1.2);
nonbonded->setSwitchingDistance(0.8);
for (int i = 0; i < numMolecules; i++) {
int startIndex = system.getNumParticles();
system.addParticle(15.6); // O
system.addParticle(0.4); // D
system.addParticle(1.0); // H1
system.addParticle(1.0); // H2
system.addParticle(0.0); // M
nonbonded->addParticle(1.71636, 0.318395, 0.21094*4.184);
nonbonded->addParticle(-1.71636, 1, 0);
nonbonded->addParticle(0.55733, 1, 0);
nonbonded->addParticle(0.55733, 1, 0);
nonbonded->addParticle(-1.11466, 1, 0);
for (int j = 0; j < 5; j++)
for (int k = 0; k < j; k++)
nonbonded->addException(startIndex+j, startIndex+k, 0, 1, 0);
system.addConstraint(startIndex, startIndex+2, 0.09572);
system.addConstraint(startIndex, startIndex+3, 0.09572);
system.addConstraint(startIndex+2, startIndex+3, 0.15139);
system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, ONE_4PI_EPS0*1.71636*1.71636/(100000*4.184), 1, 1);
}
vector<Vec3> positions;
for (int i = 0; i < gridSize; i++)
for (int j = 0; j < gridSize; j++)
for (int k = 0; k < gridSize; k++) {
Vec3 pos(i*spacing, j*spacing, k*spacing);
positions.push_back(pos);
positions.push_back(pos);
positions.push_back(pos+Vec3(0.09572, 0, 0));
positions.push_back(pos+Vec3(-0.023999, 0.092663, 0));
positions.push_back(pos);
}
int numStandardDof = 3*3*numMolecules - system.getNumConstraints();
int numDrudeDof = 3*numMolecules;
int numDof = numStandardDof+numDrudeDof;
// Simulate it and check the temperature.
int chainLength = 4;
int numMTS = 3;
int numYS = 3;
double frequency = 50.0;
double frequencyDrude = 25.0;
int randomSeed = 100;
DrudeNoseHooverIntegrator integ(0.0005, system, temperature, temperatureDrude,
frequency, frequencyDrude, chainLength, numMTS, numYS);;
Context context(system, integ, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temperature, randomSeed);
std::vector<Vec3> velocities = context.getState(State::Velocities).getVelocities();
for (int i = 0; i < numMolecules; i++){
Vec3 noize;
for (int j = 0; j < 3; j++){
noize[j] = float(((i+18311)*(j+18253) * 313419097822414) % 18313) / float(18313);
noize[j] *= sqrt(3 * BOLTZ * temperatureDrude / 0.4);
}
velocities[5*i+1] = velocities[5*i] + noize;
}
context.setVelocities(velocities);
context.applyConstraints(1e-6);
// Equilibrate.
integ.step(5000);
// Compute the internal and center of mass temperatures.
double totalKE = 0;
int numSteps = 100000;
double meanTemp = 0.0;
double meanDrudeTemp = 0.0;
double meanConserved = 0.0;
for (int i = 0; i < numSteps; i++) {
integ.step(1);
State state = context.getState(State::Energy);
double KE = state.getKineticEnergy();
double PE = state.getPotentialEnergy();
double fullKE = integ.computeTotalKineticEnergy();
double drudeKE = integ.computeDrudeKineticEnergy();
KE = fullKE - drudeKE;
double temp = KE/(0.5*numStandardDof*BOLTZ);
double drudeTemp = drudeKE/(0.5*numDrudeDof*BOLTZ);
meanTemp = (i*meanTemp + temp)/(i+1);
meanDrudeTemp = (i*meanDrudeTemp + drudeTemp)/(i+1);
double heatBathEnergy = integ.computeHeatBathEnergy();
double conserved = PE + fullKE + heatBathEnergy;
meanConserved = (i*meanConserved + conserved)/(i+1);
#define DEBUG 1
#if DEBUG
if(i%10 == 0)
std::cout << std::setw(6) << i
<< std::setprecision(8) << std::setw(16) << KE
<< std::setprecision(8) << std::setw(16) << drudeKE
<< std::setprecision(8) << std::setw(16) << meanTemp
<< std::setprecision(8) << std::setw(16) << meanDrudeTemp
<< std::setprecision(8) << std::setw(16) << conserved
<< std::endl;
#endif
totalKE += KE;
// ASSERT(fabs(meanConserved - conserved) < 0.5);
}
totalKE /= numSteps;
// ASSERT_USUALLY_EQUAL_TOL(temperature, meanTemp, 0.02);
// ASSERT_USUALLY_EQUAL_TOL(temperatureDrude, meanDrudeTemp, 0.10);
}
int main(int argc, char* argv[]) {
try {
Platform& platform = initializePlatform(argc, argv);
//registerDrudeReferenceKernelFactories();
//registerKernelFactories();
testWaterBox(platform);
}
catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
tests/TestVelocityVerletIntegrator.h → tests/TestNoseHooverIntegrator.h
View file @ 50efef13
......@@ -35,7 +35,7 @@
#include "openmm/HarmonicBondForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/System.h"
#include "openmm/VelocityVerletIntegrator.h"
#include "openmm/NoseHooverIntegrator.h"
#include "openmm/VirtualSite.h"
#include "SimTKOpenMMRealType.h"
#include "sfmt/SFMT.h"
......@@ -51,7 +51,7 @@ void testVVSingleBond() {
System system;
system.addParticle(2.0);
system.addParticle(2.0);
VelocityVerletIntegrator integrator(0.01);
NoseHooverIntegrator integrator(0.01);
HarmonicBondForce* forceField = new HarmonicBondForce();
forceField->addBond(0, 1, 1.5, 1);
system.addForce(forceField);
......@@ -86,7 +86,7 @@ void testVVConstraints() {
const int numParticles = 8;
const int numConstraints = 5;
System system;
VelocityVerletIntegrator integrator(0.001);
NoseHooverIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
......@@ -138,7 +138,7 @@ void testVVConstraints() {
void testVVConstrainedClusters() {
const int numParticles = 7;
System system;
VelocityVerletIntegrator integrator(0.001);
NoseHooverIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) {
......@@ -205,7 +205,7 @@ void testVVConstrainedMasslessParticles() {
vector<Vec3> positions(2);
positions[0] = Vec3(-1, 0, 0);
positions[1] = Vec3(1, 0, 0);
VelocityVerletIntegrator integrator(0.01);
NoseHooverIntegrator integrator(0.01);
bool failed = false;
try {
// This should throw an exception.
......@@ -250,7 +250,7 @@ void testThreeParticleVirtualSite() {
forceField->addParticle(2, params);
params[0] = 0.4;
forceField->addParticle(3, params);
VelocityVerletIntegrator integrator(0.002);
NoseHooverIntegrator integrator(0.002);
Context context(system, integrator, platform);
vector<Vec3> positions(4);
positions[0] = Vec3(0, 0, 0);
......
tests/TestNoseHooverThermostat.h
View file @ 50efef13
......@@ -31,7 +31,7 @@
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/NoseHooverChain.h"
#include "openmm/VelocityVerletIntegrator.h"
#include "openmm/NoseHooverIntegrator.h"
#include "openmm/Context.h"
#include "openmm/State.h"
#include "openmm/HarmonicBondForce.h"
......@@ -65,8 +65,8 @@ void testHarmonicOscillator() {
auto harmonic_restraint = new CustomExternalForce("0.5*(x^2+y^2+z^2)");
harmonic_restraint->addParticle(0);
system.addForce(harmonic_restraint);
VelocityVerletIntegrator integrator(0.001);
integrator.addNoseHooverChainThermostat(system, temperature, frequency, chain_length, mts, ys);
NoseHooverIntegrator integrator(0.001);
integrator.addThermostat(system, temperature, frequency, chain_length, mts, ys);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocities(velocities);
......@@ -156,15 +156,16 @@ void testDimerBox(bool constrain=true) {
}
}
}
VelocityVerletIntegrator integrator(0.001);
bool simpleConstruct = true;
double temperature = 300; // kelvin
double collisionFrequency = 25; // 1/ps
int numMTS = 3;
int numYS = 3;
int chainLength = 5;
integrator.addNoseHooverChainThermostat(system, temperature, collisionFrequency,
chainLength, numMTS, numYS);
auto integrator = simpleConstruct ? NoseHooverIntegrator(0.001, system, temperature, collisionFrequency, chainLength, numMTS, numYS)
: NoseHooverIntegrator(0.001);
if (!simpleConstruct)
integrator.addThermostat(system, temperature, collisionFrequency, chainLength, numMTS, numYS);
Context context(system, integrator, platform);
context.setPositions(positions);
context.setVelocitiesToTemperature(temperature);
......@@ -212,7 +213,7 @@ void testDimerBox(bool constrain=true) {
void testCheckpoints() {
double timeStep = 0.001;
VelocityVerletIntegrator integrator(timeStep), newIntegrator(timeStep);
NoseHooverIntegrator integrator(timeStep), newIntegrator(timeStep);
System system;
double mass = 1;
system.addParticle(mass);
......@@ -223,17 +224,11 @@ void testCheckpoints() {
system.addForce(force);
double kineticEnergy = 1e6;
double temperature=300, collisionFrequency=1, chainLength=3, numMTS=3, numYS=3;
integrator.addMaskedNoseHooverChainThermostat(system, std::vector<int>(1,0), std::vector<int>(), temperature, collisionFrequency,
chainLength, numMTS, numYS);
newIntegrator.addMaskedNoseHooverChainThermostat(system, std::vector<int>(1,0), std::vector<int>(), temperature, collisionFrequency,
chainLength, numMTS, numYS);
chainLength = 10;
integrator.addMaskedNoseHooverChainThermostat(system, std::vector<int>(1,1), std::vector<int>(1,0),
temperature, collisionFrequency,
chainLength, numMTS, numYS);
newIntegrator.addMaskedNoseHooverChainThermostat(system, std::vector<int>(1,1), std::vector<int>(1,0),
temperature, collisionFrequency,
chainLength, numMTS, numYS);
integrator.addSubsystemThermostat(system, std::vector<int>(), std::vector<std::pair<int,int>>{{0,1}}, temperature, temperature, collisionFrequency, collisionFrequency,
chainLength, numMTS, numYS);
newIntegrator.addSubsystemThermostat(system, std::vector<int>(), std::vector<std::pair<int,int>>{{0,1}}, temperature, temperature, collisionFrequency, collisionFrequency,
chainLength, numMTS, numYS);
Context context(system, integrator, platform);
Context newContext(system, newIntegrator, platform);
std::vector<Vec3> positions(2);
......
wrappers/generateWrappers.py
View file @ 50efef13
......@@ -78,7 +78,7 @@ class WrapperGenerator:
'Vec3 OpenMM::LocalCoordinatesSite::getXWeights',
'Vec3 OpenMM::LocalCoordinatesSite::getYWeights',
'std::vector<double> OpenMM::NoseHooverChain::getDefaultYoshidaSuzukiWeights',
'virtual void OpenMM::VelocityVerletIntegrator::stateChanged',
'virtual void OpenMM::NoseHooverIntegrator::stateChanged',
]
self.hideClasses = ['Kernel', 'KernelImpl', 'KernelFactory', 'ContextImpl', 'SerializationNode', 'SerializationProxy']
self.nodeByID={}
......@@ -168,6 +168,7 @@ class CHeaderGenerator(WrapperGenerator):
'std::vector< std::string >': 'OpenMM_StringArray',
'std::vector< Vec3 >': 'OpenMM_Vec3Array',
'std::vector< std::pair< int, int > >': 'OpenMM_BondArray',
'const std::vector< std::pair< int, int > >': 'OpenMM_BondArray',
'std::map< std::string, double >': 'OpenMM_ParameterArray',
'std::map< std::string, std::string >': 'OpenMM_PropertyArray',
'std::vector< double >': 'OpenMM_DoubleArray',
......
wrappers/python/src/swig_doxygen/swigInputConfig.py
View file @ 50efef13
......@@ -168,6 +168,7 @@ UNITS = {
("*", "getCutoffDistance") : ("unit.nanometers", ()),
("*", "getSwitchingDistance") : ("unit.nanometers", ()),
("*", "getDefaultCollisionFrequency") : ("1/unit.picosecond", ()),
("*", "getDefaultRelativeCollisionFrequency") : ("1/unit.picosecond", ()),
("*", "getDefaultPeriodicBoxVectors")
: (None, ('unit.nanometer', 'unit.nanometer', 'unit.nanometer')),
("*", "getDefaultPressure") : ("unit.bar", ()),
......@@ -182,6 +183,7 @@ UNITS = {
("*", "setDefaultSurfaceTension") : (None, ("unit.bar*unit.nanometer",)),
("*", "getDefaultTemperature") : ("unit.kelvin", ()),
("*", "setDefaultTemperature") : (None, ("unit.kelvin",)),
("*", "getDefaultRelativeTemperature") : ("unit.kelvin", ()),
("*", "getErrorTolerance") : (None, ()),
("*", "getEwaldErrorTolerance") : (None, ()),
("*", "getFriction") : ("1/unit.picosecond", ()),
......@@ -468,7 +470,7 @@ UNITS = {
("MonteCarloMembraneBarostat", "getXYMode") : (None, ()),
("MonteCarloMembraneBarostat", "getZMode") : (None, ()),
("DrudeLangevinIntegrator", "getDrudeFriction") : ("1/unit.picosecond", ()),
("NoseHooverChain", "getParentAtoms") : (None, ()),
("NoseHooverChain", "getThermostatedPairs") : (None, ()),
("NoseHooverChain", "getThermostatedAtoms") : (None, ()),
("NoseHooverChain", "getDefaultYoshidaSuzukiWeights") : (None, ()),
("DrudeSCFIntegrator", "getMinimizationErrorTolerance") : ("unit.kilojoules_per_mole/unit.nanometer", ()),
......
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