Skip to content

GitLab

Commit a138663d authored by Peter Eastman's avatar Peter Eastman
Browse files

Changed how kernels get initialized. StandardMMForceField now has you set the... 

Changed how kernels get initialized.  StandardMMForceField now has you set the 1-4 interactions explicitly.  CMMotionRemover allows the removal frequency to be set.
parent bf93e42f
Show whitespace changes
Inline Side-by-side
Showing with 343 additions and 443 deletions
olla/include/kernels.h
View file @ a138663d
......@@ -32,8 +32,16 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "AndersenThermostat.h"
#include "BrownianIntegrator.h"
#include "CMMotionRemover.h"
#include "GBSAOBCForceField.h"
#include "KernelImpl.h"
#include "LangevinIntegrator.h"
#include "StandardMMForceField.h"
#include "Stream.h"
#include "System.h"
#include "VerletIntegrator.h"
#include <set>
#include <string>
#include <vector>
......@@ -56,50 +64,28 @@ public:
CalcStandardMMForceFieldKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel, setting up the values of all the force field parameters.
* Initialize the kernel.
*
* @param bondIndices the two atoms connected by each bond term
* @param bondParameters the force parameters (length, k) for each bond term
* @param angleIndices the three atoms connected by each angle term
* @param angleParameters the force parameters (angle, k) for each angle term
* @param periodicTorsionIndices the four atoms connected by each periodic torsion term
* @param periodicTorsionParameters the force parameters (k, phase, periodicity) for each periodic torsion term
* @param rbTorsionIndices the four atoms connected by each Ryckaert-Bellemans torsion term
* @param rbTorsionParameters the coefficients (in order of increasing powers) for each Ryckaert-Bellemans torsion term
* @param bonded14Indices each element contains the indices of two atoms whose nonbonded interactions should be reduced since
* they form a bonded 1-4 pair
* @param lj14Scale the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
* @param coulomb14Scale the factor by which Coulomb interactions should be reduced for bonded 1-4 pairs
* @param system the System this kernel will be applied to
* @param force the StandardMMForceField this kernel will be used for
* @param exclusions the i'th element lists the indices of all atoms with which the i'th atom should not interact through
* nonbonded forces. Bonded 1-4 pairs are also included in this list, since they should be omitted from
* the standard nonbonded calculation.
* @param nonbondedParameters the nonbonded force parameters (charge, sigma, epsilon) for each atom
* @param nonbondedMethod the method to use for handling long range nonbonded interactions
* @param nonbondedCutoff the cutoff distance for nonbonded interactions (if nonbondedMethod involves a cutoff)
* @param periodicBoxSize the size of the periodic box (if nonbondedMethod involves a periodic boundary conditions)
*/
virtual void initialize(const std::vector<std::vector<int> >& bondIndices, const std::vector<std::vector<double> >& bondParameters,
const std::vector<std::vector<int> >& angleIndices, const std::vector<std::vector<double> >& angleParameters,
const std::vector<std::vector<int> >& periodicTorsionIndices, const std::vector<std::vector<double> >& periodicTorsionParameters,
const std::vector<std::vector<int> >& rbTorsionIndices, const std::vector<std::vector<double> >& rbTorsionParameters,
const std::vector<std::vector<int> >& bonded14Indices, double lj14Scale, double coulomb14Scale,
const std::vector<std::set<int> >& exclusions, const std::vector<std::vector<double> >& nonbondedParameters,
NonbondedMethod nonbondedMethod, double nonbondedCutoff, double periodicBoxSize[3]) = 0;
virtual void initialize(const System& system, const StandardMMForceField& force, const std::vector<std::set<int> >& exclusions) = 0;
/**
* Execute the kernel to calculate the forces.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param forces a Stream of type Double3 containing the force (x, y, z) on each atom. On entry, this contains the forces that
* have been calculated so far. The kernel should add its own forces to the values already in the stream.
* @param context the context in which to execute this kernel
*/
virtual void executeForces(const Stream& positions, Stream& forces) = 0;
virtual void executeForces(OpenMMContextImpl& context) = 0;
/**
* Execute the kernel to calculate the energy.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param context the context in which to execute this kernel
* @return the potential energy due to the StandardMMForceField
*/
virtual double executeEnergy(const Stream& positions) = 0;
virtual double executeEnergy(OpenMMContextImpl& context) = 0;
};
/**
......@@ -113,28 +99,25 @@ public:
CalcGBSAOBCForceFieldKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel, setting up the values of all the force field parameters.
* Initialize the kernel.
*
* @param atomParameters the force parameters (charge, atomic radius, scaling factor) for each atom
* @param solventDielectric the dielectric constant of the solvent
* @param soluteDielectric the dielectric constant of the solute
* @param system the System this kernel will be applied to
* @param force the GBSAOBCForceField this kernel will be used for
*/
virtual void initialize(const std::vector<std::vector<double> >& atomParameters, double solventDielectric, double soluteDielectric) = 0;
virtual void initialize(const System& system, const GBSAOBCForceField& force) = 0;
/**
* Execute the kernel to calculate the forces.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param forces a Stream of type Double3 containing the force (x, y, z) on each atom. On entry, this contains the forces that
* have been calculated so far. The kernel should add its own forces to the values already in the stream.
* @param context the context in which to execute this kernel
*/
virtual void executeForces(const Stream& positions, Stream& forces) = 0;
virtual void executeForces(OpenMMContextImpl& context) = 0;
/**
* Execute the kernel to calculate the energy.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param context the context in which to execute this kernel
* @return the potential energy due to the GBSAOBCForceField
*/
virtual double executeEnergy(const Stream& positions) = 0;
virtual double executeEnergy(OpenMMContextImpl& context) = 0;
};
/**
......@@ -148,23 +131,19 @@ public:
IntegrateVerletStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel, setting up all parameters related to integrator.
* Initialize the kernel.
*
* @param masses the mass of each atom
* @param constraintIndices each element contains the indices of two atoms whose distance should be constrained
* @param constraintLengths the required distance between each pair of constrained atoms
* @param system the System this kernel will be applied to
* @param integrator the VerletIntegrator this kernel will be used for
*/
virtual void initialize(const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
const std::vector<double>& constraintLengths) = 0;
virtual void initialize(const System& system, const VerletIntegrator& integrator) = 0;
/**
* Execute the kernel.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @param forces a Stream of type Double3 containing the force (x, y, z) on each atom
* @param stepSize the integration step size
* @param context the context in which to execute this kernel
* @param integrator the VerletIntegrator this kernel is being used for
*/
virtual void execute(Stream& positions, Stream& velocities, const Stream& forces, double stepSize) = 0;
virtual void execute(OpenMMContextImpl& context, const VerletIntegrator& integrator) = 0;
};
/**
......@@ -178,25 +157,19 @@ public:
IntegrateLangevinStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel, setting up all parameters related to integrator.
* Initialize the kernel.
*
* @param masses the mass of each atom
* @param constraintIndices each element contains the indices of two atoms whose distance should be constrained
* @param constraintLengths the required distance between each pair of constrained atoms
* @param system the System this kernel will be applied to
* @param integrator the LangevinIntegrator this kernel will be used for
*/
virtual void initialize(const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
const std::vector<double>& constraintLengths) = 0;
virtual void initialize(const System& system, const LangevinIntegrator& integrator) = 0;
/**
* Execute the kernel.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @param forces a Stream of type Double3 containing the force (x, y, z) on each atom
* @param temperature the temperature of the heat bath
* @param friction the friction coefficient coupling the system to the heat bath
* @param stepSize the integration step size
* @param context the context in which to execute this kernel
* @param integrator the LangevinIntegrator this kernel is being used for
*/
virtual void execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize) = 0;
virtual void execute(OpenMMContextImpl& context, const LangevinIntegrator& integrator) = 0;
};
/**
......@@ -210,25 +183,19 @@ public:
IntegrateBrownianStepKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel, setting up all parameters related to integrator.
* Initialize the kernel.
*
* @param masses the mass of each atom
* @param constraintIndices each element contains the indices of two atoms whose distance should be constrained
* @param constraintLengths the required distance between each pair of constrained atoms
* @param system the System this kernel will be applied to
* @param integrator the BrownianIntegrator this kernel will be used for
*/
virtual void initialize(const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
const std::vector<double>& constraintLengths) = 0;
virtual void initialize(const System& system, const BrownianIntegrator& integrator) = 0;
/**
* Execute the kernel.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @param forces a Stream of type Double3 containing the force (x, y, z) on each atom
* @param temperature the temperature of the heat bath
* @param friction the friction coefficient coupling the system to the heat bath
* @param stepSize the integration step size
* @param context the context in which to execute this kernel
* @param integrator the BrownianIntegrator this kernel is being used for
*/
virtual void execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize) = 0;
virtual void execute(OpenMMContextImpl& context, const BrownianIntegrator& integrator) = 0;
};
/**
......@@ -242,20 +209,18 @@ public:
ApplyAndersenThermostatKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel, setting up the values of unchanging parameters.
* Initialize the kernel.
*
* @param masses the mass of each atom
* @param system the System this kernel will be applied to
* @param thermostat the AndersenThermostat this kernel will be used for
*/
virtual void initialize(const std::vector<double>& masses) = 0;
virtual void initialize(const System& system, const AndersenThermostat& thermostat) = 0;
/**
* Execute the kernel.
*
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @param temperature the temperature of the heat bath
* @param collisionFrequency the frequency at which atom collide with particles in the heat bath
* @param stepSize the integration step size
* @param context the context in which to execute this kernel
*/
virtual void execute(Stream& velocities, double temperature, double collisionFrequency, double stepSize) = 0;
virtual void execute(OpenMMContextImpl& context) = 0;
};
/**
......@@ -269,18 +234,17 @@ public:
CalcKineticEnergyKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel, setting up the atomic masses.
* Initialize the kernel.
*
* @param masses the mass of each atom
* @param system the System this kernel will be applied to
*/
virtual void initialize(const std::vector<double>& masses) = 0;
virtual void initialize(const System& system) = 0;
/**
* Execute the kernel.
*
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @return the kinetic energy of the system
* @param context the context in which to execute this kernel
*/
virtual double execute(const Stream& velocities) = 0;
virtual double execute(OpenMMContextImpl& context) = 0;
};
/**
......@@ -294,17 +258,18 @@ public:
RemoveCMMotionKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel, setting up the atomic masses.
* Initialize the kernel.
*
* @param masses the mass of each atom
* @param system the System this kernel will be applied to
* @param force the CMMotionRemover this kernel will be used for
*/
virtual void initialize(const std::vector<double>& masses) = 0;
virtual void initialize(const System& system, const CMMotionRemover& force) = 0;
/**
* Execute the kernel.
*
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @param context the context in which to execute this kernel
*/
virtual void execute(Stream& velocities) = 0;
virtual void execute(OpenMMContextImpl& context) = 0;
};
} // namespace OpenMM
......
openmmapi/include/CMMotionRemover.h
View file @ a138663d
......@@ -48,9 +48,23 @@ public:
/**
* Create a CMMotionRemover.
*/
CMMotionRemover();
CMMotionRemover(int frequency = 1);
/**
* Get the frequency (in time steps) at which center of mass motion should be removed
*/
int getFrequency() const {
return frequency;
}
/**
* Set the frequency (in time steps) at which center of mass motion should be removed
*/
void setFrequency(int freq) {
frequency = freq;
}
protected:
ForceImpl* createImpl();
private:
int frequency;
};
} // namespace OpenMM
......
openmmapi/include/StandardMMForceField.h
View file @ a138663d
......@@ -80,8 +80,9 @@ public:
* @param numAngles the number of harmonic bond angle terms in the potential function
* @param numPeriodicTorsions the number of periodic torsion terms in the potential function
* @param numRBTorsions the number of Ryckaert-Bellemans torsion terms in the potential function
* @param numNonbonded14 the number of nonbonded 1-4 terms in the potential function
*/
StandardMMForceField(int numAtoms, int numBonds, int numAngles, int numPeriodicTorsions, int numRBTorsions);
StandardMMForceField(int numAtoms, int numBonds, int numAngles, int numPeriodicTorsions, int numRBTorsions, int numNonbonded14);
/**
* Get the number of atoms in the system.
*/
......@@ -112,10 +113,16 @@ public:
int getNumRBTorsions() const {
return rbTorsions.size();
}
/**
* Get the number of nonbonded 1-4 terms in the potential function
*/
int getNumNonbonded14() const {
return nb14s.size();
}
/**
* Get the method used for handling long range nonbonded interactions.
*/
NonbondedMethod getNonbondedMethod();
NonbondedMethod getNonbondedMethod() const;
/**
* Set the method used for handling long range nonbonded interactions.
*/
......@@ -124,7 +131,7 @@ public:
* Get the cutoff distance (in nm) being used for nonbonded interactions. If the NonbondedMethod in use
* does not use cutoffs, this value will have no effect.
*/
double getCutoffDistance();
double getCutoffDistance() const;
/**
* Set the cutoff distance (in nm) being used for nonbonded interactions. If the NonbondedMethod in use
* does not use cutoffs, this value will have no effect.
......@@ -138,7 +145,7 @@ public:
* @param y on exit, this contains the width of the periodic box along the y axis
* @param z on exit, this contains the width of the periodic box along the z axis
*/
void getPeriodicBoxSize(double& x, double& y, double& z);
void getPeriodicBoxSize(double& x, double& y, double& z) const;
/**
* Set the dimensions of the periodic box (in nm). If the NonbondedMethod in use does not use periodic
* boundary conditions, these values will have no effect.
......@@ -266,6 +273,28 @@ public:
* @param c5 the coefficient of the 5th order term
*/
void setRBTorsionParameters(int index, int atom1, int atom2, int atom3, int atom4, double c0, double c1, double c2, double c3, double c4, double c5);
/**
* Get the force field parameters for a nonbonded 1-4 term.
*
* @param index the index of the interaction for which to get parameters
* @param atom1 the index of the first atom involved in the interaction
* @param atom2 the index of the second atom involved in the interaction
* @param charge the scaled product of the atomic charges (i.e. the strength of the Coulomb interaction), measured in units of the proton charge
* @param radius the van der Waals radius of the atom (sigma in the Lennard Jones potential), measured in nm
* @param depth the well depth of the van der Waals interaction (epsilon in the Lennard Jones potential), measured in kJ/mol
*/
void getNonbonded14Parameters(int index, int& atom1, int& atom2, double& charge, double& radius, double& depth) const;
/**
* Set the force field parameters for a nonbonded 1-4 term.
*
* @param index the index of the interaction for which to get parameters
* @param atom1 the index of the first atom involved in the interaction
* @param atom2 the index of the second atom involved in the interaction
* @param charge the scaled product of the atomic charges (i.e. the strength of the Coulomb interaction), measured in units of the proton charge
* @param radius the van der Waals radius of the atom (sigma in the Lennard Jones potential), measured in nm
* @param depth the well depth of the van der Waals interaction (epsilon in the Lennard Jones potential), measured in kJ/mol
*/
void setNonbonded14Parameters(int index, int atom1, int atom2, double charge, double radius, double depth);
protected:
ForceImpl* createImpl();
private:
......@@ -274,6 +303,7 @@ private:
class AngleInfo;
class PeriodicTorsionInfo;
class RBTorsionInfo;
class NB14Info;
NonbondedMethod nonbondedMethod;
double cutoffDistance;
double periodicBoxSize[3];
......@@ -291,6 +321,7 @@ private:
std::vector<AngleInfo> angles;
std::vector<PeriodicTorsionInfo> periodicTorsions;
std::vector<RBTorsionInfo> rbTorsions;
std::vector<NB14Info> nb14s;
#if defined(_MSC_VER)
#pragma warning(pop)
......@@ -347,6 +378,16 @@ public:
}
};
class StandardMMForceField::NB14Info {
public:
int atom1, atom2;
double charge, radius, depth;
NB14Info() {
atom1 = atom2 = -1;
charge = radius = depth = 0.0;
}
};
} // namespace OpenMM
#endif /*OPENMM_STANDARDMMFORCEFIELD_H_*/
openmmapi/src/AndersenThermostatImpl.cpp
View file @ a138663d
......@@ -44,16 +44,11 @@ AndersenThermostatImpl::AndersenThermostatImpl(AndersenThermostat& owner) : owne
void AndersenThermostatImpl::initialize(OpenMMContextImpl& context) {
kernel = context.getPlatform().createKernel(ApplyAndersenThermostatKernel::Name(), context);
const System& system = context.getSystem();
vector<double> masses(system.getNumAtoms());
for (int i = 0; i < system.getNumAtoms(); ++i)
masses[i] = system.getAtomMass(i);
dynamic_cast<ApplyAndersenThermostatKernel&>(kernel.getImpl()).initialize(masses);
dynamic_cast<ApplyAndersenThermostatKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
}
void AndersenThermostatImpl::updateContextState(OpenMMContextImpl& context) {
dynamic_cast<ApplyAndersenThermostatKernel&>(kernel.getImpl()).execute(context.getVelocities(), context.getParameter(AndersenThermostat::Temperature),
context.getParameter(AndersenThermostat::CollisionFrequency), context.getIntegrator().getStepSize());
dynamic_cast<ApplyAndersenThermostatKernel&>(kernel.getImpl()).execute(context);
}
std::map<std::string, double> AndersenThermostatImpl::getDefaultParameters() {
......
openmmapi/src/BrownianIntegrator.cpp
View file @ a138663d
......@@ -62,7 +62,7 @@ void BrownianIntegrator::initialize(OpenMMContextImpl& contextRef) {
constraintIndices[i].push_back(atom2);
constraintLengths[i] = distance;
}
dynamic_cast<IntegrateBrownianStepKernel&>(kernel.getImpl()).initialize(masses, constraintIndices, constraintLengths);
dynamic_cast<IntegrateBrownianStepKernel&>(kernel.getImpl()).initialize(system, *this);
}
vector<string> BrownianIntegrator::getKernelNames() {
......@@ -75,8 +75,7 @@ void BrownianIntegrator::step(int steps) {
for (int i = 0; i < steps; ++i) {
context->updateContextState();
context->calcForces();
dynamic_cast<IntegrateBrownianStepKernel&>(kernel.getImpl()).execute(context->getPositions(), context->getVelocities(), context->getForces(), temperature,
friction, getStepSize());
dynamic_cast<IntegrateBrownianStepKernel&>(kernel.getImpl()).execute(*context, *this);
context->setTime(context->getTime()+getStepSize());
}
}
openmmapi/src/CMMotionRemover.cpp
View file @ a138663d
......@@ -34,7 +34,7 @@
using namespace OpenMM;
CMMotionRemover::CMMotionRemover() {
CMMotionRemover::CMMotionRemover(int frequency) : frequency(frequency) {
}
ForceImpl* CMMotionRemover::createImpl() {
......
openmmapi/src/CMMotionRemoverImpl.cpp
View file @ a138663d
......@@ -45,14 +45,11 @@ CMMotionRemoverImpl::CMMotionRemoverImpl(CMMotionRemover& owner) : owner(owner)
void CMMotionRemoverImpl::initialize(OpenMMContextImpl& context) {
kernel = context.getPlatform().createKernel(RemoveCMMotionKernel::Name(), context);
const System& system = context.getSystem();
vector<double> masses(system.getNumAtoms());
for (int i = 0; i < system.getNumAtoms(); ++i)
masses[i] = system.getAtomMass(i);
dynamic_cast<RemoveCMMotionKernel&>(kernel.getImpl()).initialize(masses);
dynamic_cast<RemoveCMMotionKernel&>(kernel.getImpl()).initialize(system, owner);
}
void CMMotionRemoverImpl::updateContextState(OpenMMContextImpl& context) {
dynamic_cast<RemoveCMMotionKernel&>(kernel.getImpl()).execute(context.getVelocities());
dynamic_cast<RemoveCMMotionKernel&>(kernel.getImpl()).execute(context);
}
std::vector<std::string> CMMotionRemoverImpl::getKernelNames() {
......
openmmapi/src/GBSAOBCForceFieldImpl.cpp
View file @ a138663d
......@@ -42,23 +42,15 @@ GBSAOBCForceFieldImpl::GBSAOBCForceFieldImpl(GBSAOBCForceField& owner) : owner(o
void GBSAOBCForceFieldImpl::initialize(OpenMMContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcGBSAOBCForceFieldKernel::Name(), context);
vector<vector<double> > atomParameters(owner.getNumAtoms());
for (int i = 0; i < owner.getNumAtoms(); ++i) {
double charge, radius, scalingFactor;
owner.getAtomParameters(i, charge, radius, scalingFactor);
atomParameters[i].push_back(charge);
atomParameters[i].push_back(radius);
atomParameters[i].push_back(scalingFactor);
}
dynamic_cast<CalcGBSAOBCForceFieldKernel&>(kernel.getImpl()).initialize(atomParameters, owner.getSolventDielectric(), owner.getSoluteDielectric());
dynamic_cast<CalcGBSAOBCForceFieldKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
}
void GBSAOBCForceFieldImpl::calcForces(OpenMMContextImpl& context, Stream& forces) {
dynamic_cast<CalcGBSAOBCForceFieldKernel&>(kernel.getImpl()).executeForces(context.getPositions(), forces);
dynamic_cast<CalcGBSAOBCForceFieldKernel&>(kernel.getImpl()).executeForces(context);
}
double GBSAOBCForceFieldImpl::calcEnergy(OpenMMContextImpl& context) {
return dynamic_cast<CalcGBSAOBCForceFieldKernel&>(kernel.getImpl()).executeEnergy(context.getPositions());
return dynamic_cast<CalcGBSAOBCForceFieldKernel&>(kernel.getImpl()).executeEnergy(context);
}
std::vector<std::string> GBSAOBCForceFieldImpl::getKernelNames() {
......
openmmapi/src/LangevinIntegrator.cpp
View file @ a138663d
......@@ -48,21 +48,7 @@ LangevinIntegrator::LangevinIntegrator(double temperature, double frictionCoeff,
void LangevinIntegrator::initialize(OpenMMContextImpl& contextRef) {
context = &contextRef;
kernel = context->getPlatform().createKernel(IntegrateLangevinStepKernel::Name(), contextRef);
const System& system = context->getSystem();
vector<double> masses(system.getNumAtoms());
vector<std::vector<int> > constraintIndices(system.getNumConstraints());
vector<double> constraintLengths(system.getNumConstraints());
for (int i = 0; i < system.getNumAtoms(); ++i)
masses[i] = system.getAtomMass(i);
for (int i = 0; i < system.getNumConstraints(); ++i) {
int atom1, atom2;
double distance;
system.getConstraintParameters(i, atom1, atom2, distance);
constraintIndices[i].push_back(atom1);
constraintIndices[i].push_back(atom2);
constraintLengths[i] = distance;
}
dynamic_cast<IntegrateLangevinStepKernel&>(kernel.getImpl()).initialize(masses, constraintIndices, constraintLengths);
dynamic_cast<IntegrateLangevinStepKernel&>(kernel.getImpl()).initialize(contextRef.getSystem(), *this);
}
vector<string> LangevinIntegrator::getKernelNames() {
......@@ -75,8 +61,7 @@ void LangevinIntegrator::step(int steps) {
for (int i = 0; i < steps; ++i) {
context->updateContextState();
context->calcForces();
dynamic_cast<IntegrateLangevinStepKernel&>(kernel.getImpl()).execute(context->getPositions(), context->getVelocities(), context->getForces(), temperature,
friction, getStepSize());
dynamic_cast<IntegrateLangevinStepKernel&>(kernel.getImpl()).execute(*context, *this);
context->setTime(context->getTime()+getStepSize());
}
}
openmmapi/src/OpenMMContextImpl.cpp
View file @ a138663d
......@@ -70,7 +70,7 @@ OpenMMContextImpl::OpenMMContextImpl(OpenMMContext& owner, System& system, Integ
vector<double> masses(system.getNumAtoms());
for (size_t i = 0; i < masses.size(); ++i)
masses[i] = system.getAtomMass(i);
dynamic_cast<CalcKineticEnergyKernel&>(kineticEnergyKernel.getImpl()).initialize(masses);
dynamic_cast<CalcKineticEnergyKernel&>(kineticEnergyKernel.getImpl()).initialize(system);
for (size_t i = 0; i < forceImpls.size(); ++i)
forceImpls[i]->initialize(*this);
integrator.initialize(*this);
......@@ -107,7 +107,7 @@ void OpenMMContextImpl::calcForces() {
}
double OpenMMContextImpl::calcKineticEnergy() {
return dynamic_cast<CalcKineticEnergyKernel&>(kineticEnergyKernel.getImpl()).execute(velocities);
return dynamic_cast<CalcKineticEnergyKernel&>(kineticEnergyKernel.getImpl()).execute(*this);
}
double OpenMMContextImpl::calcPotentialEnergy() {
......
openmmapi/src/StandardMMForceField.cpp
View file @ a138663d
......@@ -36,13 +36,13 @@
using namespace OpenMM;
StandardMMForceField::StandardMMForceField(int numAtoms, int numBonds, int numAngles, int numPeriodicTorsions, int numRBTorsions) :
atoms(numAtoms), bonds(numBonds), angles(numAngles), periodicTorsions(numPeriodicTorsions), rbTorsions(numRBTorsions),
StandardMMForceField::StandardMMForceField(int numAtoms, int numBonds, int numAngles, int numPeriodicTorsions, int numRBTorsions, int numNonbonded14) :
atoms(numAtoms), bonds(numBonds), angles(numAngles), periodicTorsions(numPeriodicTorsions), rbTorsions(numRBTorsions), nb14s(numNonbonded14),
nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
periodicBoxSize[0] = periodicBoxSize[1] = periodicBoxSize[2] = 2.0;
}
StandardMMForceField::NonbondedMethod StandardMMForceField::getNonbondedMethod() {
StandardMMForceField::NonbondedMethod StandardMMForceField::getNonbondedMethod() const {
return nonbondedMethod;
}
......@@ -50,7 +50,7 @@ void StandardMMForceField::setNonbondedMethod(NonbondedMethod method) {
nonbondedMethod = method;
}
double StandardMMForceField::getCutoffDistance() {
double StandardMMForceField::getCutoffDistance() const {
return cutoffDistance;
}
......@@ -58,7 +58,7 @@ void StandardMMForceField::setCutoffDistance(double distance) {
cutoffDistance = distance;
}
void StandardMMForceField::getPeriodicBoxSize(double& x, double& y, double& z) {
void StandardMMForceField::getPeriodicBoxSize(double& x, double& y, double& z) const {
x = periodicBoxSize[0];
y = periodicBoxSize[1];
z = periodicBoxSize[2];
......@@ -158,6 +158,22 @@ void StandardMMForceField::setRBTorsionParameters(int index, int atom1, int atom
rbTorsions[index].c[5] = c5;
}
void StandardMMForceField::getNonbonded14Parameters(int index, int& atom1, int& atom2, double& charge, double& radius, double& depth) const {
atom1 = nb14s[index].atom1;
atom2 = nb14s[index].atom2;
charge = nb14s[index].charge;
radius = nb14s[index].radius;
depth = nb14s[index].depth;
}
void StandardMMForceField::setNonbonded14Parameters(int index, int atom1, int atom2, double charge, double radius, double depth) {
nb14s[index].atom1 = atom1;
nb14s[index].atom2 = atom2;
nb14s[index].charge = charge;
nb14s[index].radius = radius;
nb14s[index].depth = depth;
}
ForceImpl* StandardMMForceField::createImpl() {
return new StandardMMForceFieldImpl(*this);
}
openmmapi/src/StandardMMForceFieldImpl.cpp
View file @ a138663d
......@@ -46,91 +46,26 @@ StandardMMForceFieldImpl::~StandardMMForceFieldImpl() {
void StandardMMForceFieldImpl::initialize(OpenMMContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcStandardMMForceFieldKernel::Name(), context);
vector<vector<int> > bondIndices(owner.getNumBonds());
vector<vector<double> > bondParameters(owner.getNumBonds());
vector<vector<int> > angleIndices(owner.getNumAngles());
vector<vector<double> > angleParameters(owner.getNumAngles());
vector<vector<int> > periodicTorsionIndices(owner.getNumPeriodicTorsions());
vector<vector<double> > periodicTorsionParameters(owner.getNumPeriodicTorsions());
vector<vector<int> > rbTorsionIndices(owner.getNumRBTorsions());
vector<vector<double> > rbTorsionParameters(owner.getNumRBTorsions());
vector<vector<int> > bonded14Indices;
vector<set<int> > exclusions(owner.getNumAtoms());
vector<vector<double> > nonbondedParameters(owner.getNumAtoms());
vector<vector<int> > bondIndices(owner.getNumBonds());
set<pair<int, int> > bonded14set;
for (int i = 0; i < owner.getNumBonds(); ++i) {
int atom1, atom2;
double length, k;
owner.getBondParameters(i, atom1, atom2, length, k);
bondIndices[i].push_back(atom1);
bondIndices[i].push_back(atom2);
bondParameters[i].push_back(length);
bondParameters[i].push_back(k);
}
for (int i = 0; i < owner.getNumAngles(); ++i) {
int atom1, atom2, atom3;
double angle, k;
owner.getAngleParameters(i, atom1, atom2, atom3, angle, k);
angleIndices[i].push_back(atom1);
angleIndices[i].push_back(atom2);
angleIndices[i].push_back(atom3);
angleParameters[i].push_back(angle);
angleParameters[i].push_back(k);
}
for (int i = 0; i < owner.getNumPeriodicTorsions(); ++i) {
int atom1, atom2, atom3, atom4, periodicity;
double phase, k;
owner.getPeriodicTorsionParameters(i, atom1, atom2, atom3, atom4, periodicity, phase, k);
periodicTorsionIndices[i].push_back(atom1);
periodicTorsionIndices[i].push_back(atom2);
periodicTorsionIndices[i].push_back(atom3);
periodicTorsionIndices[i].push_back(atom4);
periodicTorsionParameters[i].push_back(k);
periodicTorsionParameters[i].push_back(phase);
periodicTorsionParameters[i].push_back(periodicity);
}
for (int i = 0; i < owner.getNumRBTorsions(); ++i) {
int atom1, atom2, atom3, atom4;
double c0, c1, c2, c3, c4, c5;
owner.getRBTorsionParameters(i, atom1, atom2, atom3, atom4, c0, c1, c2, c3, c4, c5);
rbTorsionIndices[i].push_back(atom1);
rbTorsionIndices[i].push_back(atom2);
rbTorsionIndices[i].push_back(atom3);
rbTorsionIndices[i].push_back(atom4);
rbTorsionParameters[i].push_back(c0);
rbTorsionParameters[i].push_back(c1);
rbTorsionParameters[i].push_back(c2);
rbTorsionParameters[i].push_back(c3);
rbTorsionParameters[i].push_back(c4);
rbTorsionParameters[i].push_back(c5);
}
for (int i = 0; i < owner.getNumAtoms(); ++i) {
double charge, radius, depth;
owner.getAtomParameters(i, charge, radius, depth);
nonbondedParameters[i].push_back(charge);
nonbondedParameters[i].push_back(radius);
nonbondedParameters[i].push_back(depth);
}
set<pair<int, int> > bonded14set;
findExclusions(bondIndices, exclusions, bonded14set);
bonded14Indices.resize(bonded14set.size());
int index = 0;
for (set<pair<int, int> >::const_iterator iter = bonded14set.begin(); iter != bonded14set.end(); ++iter) {
bonded14Indices[index].push_back(iter->first);
bonded14Indices[index++].push_back(iter->second);
}
double boxSize[3];
owner.getPeriodicBoxSize(boxSize[0], boxSize[1], boxSize[2]);
dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).initialize(bondIndices, bondParameters, angleIndices, angleParameters,
periodicTorsionIndices, periodicTorsionParameters, rbTorsionIndices, rbTorsionParameters, bonded14Indices, 0.5, 1.0/1.2, exclusions,
nonbondedParameters, CalcStandardMMForceFieldKernel::NonbondedMethod(owner.getNonbondedMethod()), owner.getCutoffDistance(), boxSize);
dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner, exclusions);
}
void StandardMMForceFieldImpl::calcForces(OpenMMContextImpl& context, Stream& forces) {
dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).executeForces(context.getPositions(), forces);
dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).executeForces(context);
}
double StandardMMForceFieldImpl::calcEnergy(OpenMMContextImpl& context) {
return dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).executeEnergy(context.getPositions());
return dynamic_cast<CalcStandardMMForceFieldKernel&>(kernel.getImpl()).executeEnergy(context);
}
std::vector<std::string> StandardMMForceFieldImpl::getKernelNames() {
......
openmmapi/src/VerletIntegrator.cpp
View file @ a138663d
......@@ -46,21 +46,7 @@ VerletIntegrator::VerletIntegrator(double stepSize) {
void VerletIntegrator::initialize(OpenMMContextImpl& contextRef) {
context = &contextRef;
kernel = context->getPlatform().createKernel(IntegrateVerletStepKernel::Name(), contextRef);
const System& system = context->getSystem();
vector<double> masses(system.getNumAtoms());
vector<std::vector<int> > constraintIndices(system.getNumConstraints());
vector<double> constraintLengths(system.getNumConstraints());
for (int i = 0; i < system.getNumAtoms(); ++i)
masses[i] = system.getAtomMass(i);
for (int i = 0; i < system.getNumConstraints(); ++i) {
int atom1, atom2;
double distance;
system.getConstraintParameters(i, atom1, atom2, distance);
constraintIndices[i].push_back(atom1);
constraintIndices[i].push_back(atom2);
constraintLengths[i] = distance;
}
dynamic_cast<IntegrateVerletStepKernel&>(kernel.getImpl()).initialize(masses, constraintIndices, constraintLengths);
dynamic_cast<IntegrateVerletStepKernel&>(kernel.getImpl()).initialize(contextRef.getSystem(), *this);
}
vector<string> VerletIntegrator::getKernelNames() {
......@@ -73,7 +59,7 @@ void VerletIntegrator::step(int steps) {
for (int i = 0; i < steps; ++i) {
context->updateContextState();
context->calcForces();
dynamic_cast<IntegrateVerletStepKernel&>(kernel.getImpl()).execute(context->getPositions(), context->getVelocities(), context->getForces(), getStepSize());
dynamic_cast<IntegrateVerletStepKernel&>(kernel.getImpl()).execute(*context, *this);
context->setTime(context->getTime()+getStepSize());
}
}
platforms/cuda/include/CudaPlatform.h
View file @ a138663d
......@@ -68,6 +68,7 @@ public:
_gpuContext* gpu;
bool removeCM;
bool useOBC;
int cmMotionFrequency;
};
} // namespace OpenMM
......
platforms/cuda/src/CudaKernels.cpp
View file @ a138663d
......@@ -35,6 +35,7 @@
#include "CudaStreamImpl.h"
#include "LangevinIntegrator.h"
#include "ReferencePlatform.h"
#include "internal/OpenMMContextImpl.h"
#include "kernels/gputypes.h"
#include "kernels/cudaKernels.h"
#include <cmath>
......@@ -47,19 +48,13 @@ using namespace std;
CudaCalcStandardMMForceFieldKernel::~CudaCalcStandardMMForceFieldKernel() {
}
void CudaCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >& bondIndices, const vector<vector<double> >& bondParameters,
const vector<vector<int> >& angleIndices, const vector<vector<double> >& angleParameters,
const vector<vector<int> >& periodicTorsionIndices, const vector<vector<double> >& periodicTorsionParameters,
const vector<vector<int> >& rbTorsionIndices, const vector<vector<double> >& rbTorsionParameters,
const vector<vector<int> >& bonded14Indices, double lj14Scale, double coulomb14Scale,
const vector<set<int> >& exclusions, const vector<vector<double> >& nonbondedParameters,
NonbondedMethod nonbondedMethod, double nonbondedCutoff, double periodicBoxSize[3]) {
numAtoms = nonbondedParameters.size();
numBonds = bondIndices.size();
numAngles = angleIndices.size();
numPeriodicTorsions = periodicTorsionIndices.size();
numRBTorsions = rbTorsionIndices.size();
num14 = bonded14Indices.size();
void CudaCalcStandardMMForceFieldKernel::initialize(const System& system, const StandardMMForceField& force, const std::vector<std::set<int> >& exclusions) {
numAtoms = force.getNumAtoms();
numBonds = force.getNumBonds();
numAngles = force.getNumAngles();
numPeriodicTorsions = force.getNumPeriodicTorsions();
numRBTorsions = force.getNumRBTorsions();
num14 = force.getNumNonbonded14();
const float RadiansToDegrees = 180.0/3.14159265;
_gpuContext* gpu = data.gpu;
......@@ -71,10 +66,10 @@ void CudaCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >&
vector<float> length(numBonds);
vector<float> k(numBonds);
for (int i = 0; i < numBonds; i++) {
atom1[i] = bondIndices[i][0];
atom2[i] = bondIndices[i][1];
length[i] = (float) bondParameters[i][0];
k[i] = (float) bondParameters[i][1];
double lengthValue, kValue;
force.getBondParameters(i, atom1[i], atom2[i], lengthValue, kValue);
length[i] = (float) lengthValue;
k[i] = (float) kValue;
}
gpuSetBondParameters(gpu, atom1, atom2, length, k);
}
......@@ -88,11 +83,10 @@ void CudaCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >&
vector<float> angle(numAngles);
vector<float> k(numAngles);
for (int i = 0; i < numAngles; i++) {
atom1[i] = angleIndices[i][0];
atom2[i] = angleIndices[i][1];
atom3[i] = angleIndices[i][2];
angle[i] = (float) (angleParameters[i][0]*RadiansToDegrees);
k[i] = (float) angleParameters[i][1];
double angleValue, kValue;
force.getAngleParameters(i, atom1[i], atom2[i], atom3[i], angleValue, kValue);
angle[i] = (float) (angleValue*RadiansToDegrees);
k[i] = (float) kValue;
}
gpuSetBondAngleParameters(gpu, atom1, atom2, atom3, angle, k);
}
......@@ -108,13 +102,10 @@ void CudaCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >&
vector<float> phase(numPeriodicTorsions);
vector<int> periodicity(numPeriodicTorsions);
for (int i = 0; i < numPeriodicTorsions; i++) {
atom1[i] = periodicTorsionIndices[i][0];
atom2[i] = periodicTorsionIndices[i][1];
atom3[i] = periodicTorsionIndices[i][2];
atom4[i] = periodicTorsionIndices[i][3];
k[i] = (float) periodicTorsionParameters[i][0];
phase[i] = (float) (periodicTorsionParameters[i][1]*RadiansToDegrees);
periodicity[i] = (int) periodicTorsionParameters[i][2];
double kValue, phaseValue;
force.getPeriodicTorsionParameters(i, atom1[i], atom2[i], atom3[i], atom4[i], periodicity[i], phaseValue, kValue);
k[i] = (float) kValue;
phase[i] = (float) (phaseValue*RadiansToDegrees);
}
gpuSetDihedralParameters(gpu, atom1, atom2, atom3, atom4, k, phase, periodicity);
}
......@@ -133,16 +124,14 @@ void CudaCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >&
vector<float> c4(numRBTorsions);
vector<float> c5(numRBTorsions);
for (int i = 0; i < numRBTorsions; i++) {
atom1[i] = rbTorsionIndices[i][0];
atom2[i] = rbTorsionIndices[i][1];
atom3[i] = rbTorsionIndices[i][2];
atom4[i] = rbTorsionIndices[i][3];
c0[i] = (float) rbTorsionParameters[i][0];
c1[i] = (float) rbTorsionParameters[i][1];
c2[i] = (float) rbTorsionParameters[i][2];
c3[i] = (float) rbTorsionParameters[i][3];
c4[i] = (float) rbTorsionParameters[i][4];
c5[i] = (float) rbTorsionParameters[i][5];
double c[6];
force.getRBTorsionParameters(i, atom1[i], atom2[i], atom3[i], atom4[i], c[0], c[1], c[2], c[3], c[4], c[5]);
c0[i] = (float) c[0];
c1[i] = (float) c[1];
c2[i] = (float) c[2];
c3[i] = (float) c[3];
c4[i] = (float) c[4];
c5[i] = (float) c[5];
}
gpuSetRbDihedralParameters(gpu, atom1, atom2, atom3, atom4, c0, c1, c2, c3, c4, c5);
}
......@@ -157,10 +146,12 @@ void CudaCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >&
vector<char> symbol;
vector<vector<int> > exclusionList(numAtoms);
for (int i = 0; i < numAtoms; i++) {
double charge, radius, depth;
force.getAtomParameters(i, charge, radius, depth);
atom[i] = i;
q[i] = (float) nonbondedParameters[i][0];
c6[i] = (float) (4*nonbondedParameters[i][2]*pow(nonbondedParameters[i][1], 6.0));
c12[i] = (float) (4*nonbondedParameters[i][2]*pow(nonbondedParameters[i][1], 12.0));
q[i] = (float) charge;
c6[i] = (float) (4*depth*pow(radius, 6.0));
c12[i] = (float) (4*depth*pow(radius, 12.0));
exclusionList[i] = vector<int>(exclusions[i].begin(), exclusions[i].end());
exclusionList[i].push_back(i);
}
......@@ -177,20 +168,19 @@ void CudaCalcStandardMMForceFieldKernel::initialize(const vector<vector<int> >&
vector<float> q1(num14);
vector<float> q2(num14);
for (int i = 0; i < num14; i++) {
atom1[i] = bonded14Indices[i][0];
atom2[i] = bonded14Indices[i][1];
double sig = 0.5*(nonbondedParameters[atom1[i]][1]+nonbondedParameters[atom2[i]][1]);
double eps = sqrt(nonbondedParameters[atom1[i]][2]*nonbondedParameters[atom2[i]][2]);
c6[i] = (float) (4*eps*pow(sig, 6.0)*lj14Scale);
c12[i] = (float) (4*eps*pow(sig, 12.0)*lj14Scale);
q1[i] = (float) nonbondedParameters[atom1[i]][0];
q2[i] = (float) nonbondedParameters[atom2[i]][0];
double charge, sig, eps;
force.getNonbonded14Parameters(i, atom1[i], atom2[i], charge, sig, eps);
c6[i] = (float) (4*eps*pow(sig, 6.0));
c12[i] = (float) (4*eps*pow(sig, 12.0));
float q = (float) std::sqrt(charge);
q1[i] = q;
q2[i] = q;
}
gpuSetLJ14Parameters(gpu, 138.935485f, (float) coulomb14Scale, atom1, atom2, c6, c12, q1, q2);
gpuSetLJ14Parameters(gpu, 138.935485f, 1.0f, atom1, atom2, c6, c12, q1, q2);
}
}
void CudaCalcStandardMMForceFieldKernel::executeForces(const Stream& positions, Stream& forces) {
void CudaCalcStandardMMForceFieldKernel::executeForces(OpenMMContextImpl& context) {
_gpuContext* gpu = data.gpu;
if (data.useOBC) {
kCalculateCDLJObcGbsaForces1(gpu);
......@@ -205,84 +195,90 @@ void CudaCalcStandardMMForceFieldKernel::executeForces(const Stream& positions,
kReduceBornSumAndForces(gpu);
}
double CudaCalcStandardMMForceFieldKernel::executeEnergy(const Stream& positions) {
double CudaCalcStandardMMForceFieldKernel::executeEnergy(OpenMMContextImpl& context) {
// We don't currently have GPU kernels to calculate energy, so instead we have the reference
// platform do it. This is VERY slow.
LangevinIntegrator integrator(0.0, 1.0, 0.0);
ReferencePlatform platform;
OpenMMContext context(system, integrator, platform);
OpenMMContext refContext(system, integrator, platform);
const Stream& positions = context.getPositions();
double* posData = new double[positions.getSize()*3];
positions.saveToArray(posData);
vector<Vec3> pos(positions.getSize());
for (int i = 0; i < pos.size(); i++)
pos[i] = Vec3(posData[3*i], posData[3*i+1], posData[3*i+2]);
delete[] posData;
context.setPositions(pos);
return context.getState(State::Energy).getPotentialEnergy();
refContext.setPositions(pos);
return refContext.getState(State::Energy).getPotentialEnergy();
}
CudaCalcGBSAOBCForceFieldKernel::~CudaCalcGBSAOBCForceFieldKernel() {
}
void CudaCalcGBSAOBCForceFieldKernel::initialize(const vector<vector<double> >& atomParameters, double solventDielectric, double soluteDielectric) {
int numAtoms = atomParameters.size();
void CudaCalcGBSAOBCForceFieldKernel::initialize(const System& system, const GBSAOBCForceField& force) {
int numAtoms = system.getNumAtoms();
_gpuContext* gpu = data.gpu;
vector<int> atom(numAtoms);
vector<float> radius(numAtoms);
vector<float> scale(numAtoms);
for (int i = 0; i < numAtoms; i++) {
double charge, atomRadius, scalingFactor;
force.getAtomParameters(i, charge, atomRadius, scalingFactor);
atom[i] = i;
radius[i] = (float) atomParameters[i][1];
scale[i] = (float) atomParameters[i][2];
radius[i] = (float) atomRadius;
scale[i] = (float) scalingFactor;
}
gpuSetObcParameters(gpu, soluteDielectric, solventDielectric, atom, radius, scale);
gpuSetObcParameters(gpu, force.getSoluteDielectric(), force.getSolventDielectric(), atom, radius, scale);
data.useOBC = true;
}
void CudaCalcGBSAOBCForceFieldKernel::executeForces(const Stream& positions, Stream& forces) {
void CudaCalcGBSAOBCForceFieldKernel::executeForces(OpenMMContextImpl& context) {
}
double CudaCalcGBSAOBCForceFieldKernel::executeEnergy(const Stream& positions) {
double CudaCalcGBSAOBCForceFieldKernel::executeEnergy(OpenMMContextImpl& context) {
}
//CudaIntegrateVerletStepKernel::~CudaIntegrateVerletStepKernel() {
//}
//
//void CudaIntegrateVerletStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
//void CudaIntegrateVerletStepKernel::initialize(const System& system, const VerletIntegrator& integrator) {
//}
//
//void CudaIntegrateVerletStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double stepSize) {
//void CudaIntegrateVerletStepKernel::execute(OpenMMContextImpl& context, const VerletIntegrator& integrator) {
//}
CudaIntegrateLangevinStepKernel::~CudaIntegrateLangevinStepKernel() {
}
void CudaIntegrateLangevinStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
const vector<double>& constraintLengths) {
void CudaIntegrateLangevinStepKernel::initialize(const System& system, const LangevinIntegrator& integrator) {
// Set masses.
_gpuContext* gpu = data.gpu;
vector<float> mass(masses.size());
for (int i = 0; i < (int) mass.size(); i++)
mass[i] = (float) masses[i];
int numAtoms = system.getNumAtoms();
vector<float> mass(numAtoms);
for (int i = 0; i < numAtoms; i++)
mass[i] = (float) system.getAtomMass(i);
gpuSetMass(gpu, mass);
// Set constraints.
int numConstraints = constraintLengths.size();
int numConstraints = system.getNumConstraints();
vector<int> atom1(numConstraints);
vector<int> atom2(numConstraints);
vector<float> distance(numConstraints);
vector<float> invMass1(numConstraints);
vector<float> invMass2(numConstraints);
for (int i = 0; i < numConstraints; i++) {
atom1[i] = constraintIndices[i][0];
atom2[i] = constraintIndices[i][1];
distance[i] = (float) constraintLengths[i];
invMass1[i] = 1.0f/mass[atom1[i]];
invMass2[i] = 1.0f/mass[atom2[i]];
int atom1Index, atom2Index;
double constraintDistance;
system.getConstraintParameters(i, atom1Index, atom2Index, constraintDistance);
atom1[i] = atom1Index;
atom2[i] = atom2Index;
distance[i] = (float) constraintDistance;
invMass1[i] = 1.0f/mass[atom1Index];
invMass2[i] = 1.0f/mass[atom2Index];
}
gpuSetShakeParameters(gpu, atom1, atom2, distance, invMass1, invMass2);
gpuBuildThreadBlockWorkList(gpu);
......@@ -297,8 +293,11 @@ void CudaIntegrateLangevinStepKernel::initialize(const vector<double>& masses, c
prevStepSize = -1.0;
}
void CudaIntegrateLangevinStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize) {
void CudaIntegrateLangevinStepKernel::execute(OpenMMContextImpl& context, const LangevinIntegrator& integrator) {
_gpuContext* gpu = data.gpu;
double temperature = integrator.getTemperature();
double friction = integrator.getFriction();
double stepSize = integrator.getStepSize();
if (temperature != prevTemp || friction != prevFriction || stepSize != prevStepSize) {
// Initialize the GPU parameters.
......@@ -312,8 +311,11 @@ void CudaIntegrateLangevinStepKernel::execute(Stream& positions, Stream& velocit
}
kUpdatePart1(gpu);
kApplyFirstShake(gpu);
if (data.removeCM)
if (data.removeCM) {
int step = context.getTime()/stepSize;
if (step%data.cmMotionFrequency == 0)
gpu->bCalculateCM = true;
}
kUpdatePart2(gpu);
kApplySecondShake(gpu);
}
......@@ -321,30 +323,33 @@ void CudaIntegrateLangevinStepKernel::execute(Stream& positions, Stream& velocit
//CudaIntegrateBrownianStepKernel::~CudaIntegrateBrownianStepKernel() {
//}
//
//void CudaIntegrateBrownianStepKernel::initialize(const vector<double>& masses, const vector<vector<int> >& constraintIndices,
// const vector<double>& constraintLengths) {
//void CudaIntegrateBrownianStepKernel::initialize(const System& system, const BrownianIntegrator& integrator) {
//}
//
//void CudaIntegrateBrownianStepKernel::execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize) {
//void CudaIntegrateBrownianStepKernel::execute(OpenMMContextImpl& context, const BrownianIntegrator& integrator) {
//}
//
//CudaApplyAndersenThermostatKernel::~CudaApplyAndersenThermostatKernel() {
//}
//
//void CudaApplyAndersenThermostatKernel::initialize(const vector<double>& masses) {
//void CudaApplyAndersenThermostatKernel::initialize(const System& system, const AndersenThermostat& thermostat) {
//}
//
//void CudaApplyAndersenThermostatKernel::execute(Stream& velocities, double temperature, double collisionFrequency, double stepSize) {
//void CudaApplyAndersenThermostatKernel::execute(OpenMMContextImpl& context) {
//}
void CudaCalcKineticEnergyKernel::initialize(const vector<double>& masses) {
this->masses = masses;
void CudaCalcKineticEnergyKernel::initialize(const System& system) {
int numAtoms = system.getNumAtoms();
masses.resize(numAtoms);
for (size_t i = 0; i < numAtoms; ++i)
masses[i] = system.getAtomMass(i);
}
double CudaCalcKineticEnergyKernel::execute(const Stream& velocities) {
double CudaCalcKineticEnergyKernel::execute(OpenMMContextImpl& context) {
// We don't currently have a GPU kernel to do this, so we retrieve the velocities and calculate the energy
// on the CPU.
const Stream& velocities = context.getVelocities();
double* v = new double[velocities.getSize()*3];
velocities.saveToArray(v);
double energy = 0.0;
......@@ -354,9 +359,10 @@ double CudaCalcKineticEnergyKernel::execute(const Stream& velocities) {
return 0.5*energy;
}
void CudaRemoveCMMotionKernel::initialize(const vector<double>& masses) {
void CudaRemoveCMMotionKernel::initialize(const System& system, const CMMotionRemover& force) {
data.removeCM = true;
data.cmMotionFrequency = force.getFrequency();
}
void CudaRemoveCMMotionKernel::execute(Stream& velocities) {
void CudaRemoveCMMotionKernel::execute(OpenMMContextImpl& context) {
}
platforms/cuda/src/CudaKernels.h
View file @ a138663d
......@@ -54,47 +54,28 @@ public:
}
~CudaCalcStandardMMForceFieldKernel();
/**
* Initialize the kernel, setting up the values of all the force field parameters.
* Initialize the kernel.
*
* @param bondIndices the two atoms connected by each bond term
* @param bondParameters the force parameters (length, k) for each bond term
* @param angleIndices the three atoms connected by each angle term
* @param angleParameters the force parameters (angle, k) for each angle term
* @param periodicTorsionIndices the four atoms connected by each periodic torsion term
* @param periodicTorsionParameters the force parameters (k, phase, periodicity) for each periodic torsion term
* @param rbTorsionIndices the four atoms connected by each Ryckaert-Bellemans torsion term
* @param rbTorsionParameters the coefficients (in order of increasing powers) for each Ryckaert-Bellemans torsion term
* @param bonded14Indices each element contains the indices of two atoms whose nonbonded interactions should be reduced since
* they form a bonded 1-4 pair
* @param lj14Scale the factor by which van der Waals interactions should be reduced for bonded 1-4 pairs
* @param coulomb14Scale the factor by which Coulomb interactions should be reduced for bonded 1-4 pairs
* @param system the System this kernel will be applied to
* @param force the StandardMMForceField this kernel will be used for
* @param exclusions the i'th element lists the indices of all atoms with which the i'th atom should not interact through
* nonbonded forces. Bonded 1-4 pairs are also included in this list, since they should be omitted from
* the standard nonbonded calculation.
* @param nonbondedParameters the nonbonded force parameters (charge, sigma, epsilon) for each atom
*/
void initialize(const std::vector<std::vector<int> >& bondIndices, const std::vector<std::vector<double> >& bondParameters,
const std::vector<std::vector<int> >& angleIndices, const std::vector<std::vector<double> >& angleParameters,
const std::vector<std::vector<int> >& periodicTorsionIndices, const std::vector<std::vector<double> >& periodicTorsionParameters,
const std::vector<std::vector<int> >& rbTorsionIndices, const std::vector<std::vector<double> >& rbTorsionParameters,
const std::vector<std::vector<int> >& bonded14Indices, double lj14Scale, double coulomb14Scale,
const std::vector<std::set<int> >& exclusions, const std::vector<std::vector<double> >& nonbondedParameters,
NonbondedMethod nonbondedMethod, double nonbondedCutoff, double periodicBoxSize[3]);
void initialize(const System& system, const StandardMMForceField& force, const std::vector<std::set<int> >& exclusions);
/**
* Execute the kernel to calculate the forces.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param forces a Stream of type Double3 containing the force (x, y, z) on each atom. On entry, this contains the forces that
* have been calculated so far. The kernel should add its own forces to the values already in the stream.
* @param context the context in which to execute this kernel
*/
void executeForces(const Stream& positions, Stream& forces);
void executeForces(OpenMMContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param context the context in which to execute this kernel
* @return the potential energy due to the StandardMMForceField
*/
double executeEnergy(const Stream& positions);
double executeEnergy(OpenMMContextImpl& context);
private:
CudaPlatform::PlatformData& data;
int numAtoms, numBonds, numAngles, numPeriodicTorsions, numRBTorsions, num14;
......@@ -110,28 +91,25 @@ public:
}
~CudaCalcGBSAOBCForceFieldKernel();
/**
* Initialize the kernel, setting up the values of all the force field parameters.
* Initialize the kernel.
*
* @param atomParameters the force parameters (charge, atomic radius, scaling factor) for each atom
* @param solventDielectric the dielectric constant of the solvent
* @param soluteDielectric the dielectric constant of the solute
* @param system the System this kernel will be applied to
* @param force the GBSAOBCForceField this kernel will be used for
*/
void initialize(const std::vector<std::vector<double> >& atomParameters, double solventDielectric, double soluteDielectric);
void initialize(const System& system, const GBSAOBCForceField& force);
/**
* Execute the kernel to calculate the forces.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param forces a Stream of type Double3 containing the force (x, y, z) on each atom. On entry, this contains the forces that
* have been calculated so far. The kernel should add its own forces to the values already in the stream.
* @param context the context in which to execute this kernel
*/
void executeForces(const Stream& positions, Stream& forces);
void executeForces(OpenMMContextImpl& context);
/**
* Execute the kernel to calculate the energy.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param context the context in which to execute this kernel
* @return the potential energy due to the GBSAOBCForceField
*/
double executeEnergy(const Stream& positions);
double executeEnergy(OpenMMContextImpl& context);
private:
CudaPlatform::PlatformData& data;
};
......@@ -146,23 +124,19 @@ private:
// }
// ~CudaIntegrateVerletStepKernel();
// /**
// * Initialize the kernel, setting up all parameters related to integrator.
// * Initialize the kernel.
// *
// * @param masses the mass of each atom
// * @param constraintIndices each element contains the indices of two atoms whose distance should be constrained
// * @param constraintLengths the required distance between each pair of constrained atoms
// * @param system the System this kernel will be applied to
// * @param integrator the VerletIntegrator this kernel will be used for
// */
// void initialize(const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
// const std::vector<double>& constraintLengths);
// void initialize(const System& system, const VerletIntegrator& integrator);
// /**
// * Execute the kernel.
// *
// * @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
// * @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
// * @param forces a Stream of type Double3 containing the force (x, y, z) on each atom
// * @param stepSize the integration step size
// * @param context the context in which to execute this kernel
// * @param integrator the VerletIntegrator this kernel is being used for
// */
// void execute(Stream& positions, Stream& velocities, const Stream& forces, double stepSize);
// void execute(OpenMMContextImpl& context, const VerletIntegrator& integrator);
//private:
// CudaVerletDynamics* dynamics;
// CudaShakeAlgorithm* shake;
......@@ -182,25 +156,19 @@ public:
}
~CudaIntegrateLangevinStepKernel();
/**
* Initialize the kernel, setting up all parameters related to integrator.
* Initialize the kernel, setting up the atomic masses.
*
* @param masses the mass of each atom
* @param constraintIndices each element contains the indices of two atoms whose distance should be constrained
* @param constraintLengths the required distance between each pair of constrained atoms
* @param system the System this kernel will be applied to
* @param integrator the LangevinIntegrator this kernel will be used for
*/
void initialize(const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
const std::vector<double>& constraintLengths);
void initialize(const System& system, const LangevinIntegrator& integrator);
/**
* Execute the kernel.
*
* @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @param forces a Stream of type Double3 containing the force (x, y, z) on each atom
* @param temperature the temperature of the heat bath
* @param friction the friction coefficient coupling the system to the heat bath
* @param stepSize the integration step size
* @param context the context in which to execute this kernel
* @param integrator the LangevinIntegrator this kernel is being used for
*/
void execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize);
void execute(OpenMMContextImpl& context, const LangevinIntegrator& integrator);
private:
CudaPlatform::PlatformData& data;
double prevTemp, prevFriction, prevStepSize;
......@@ -216,25 +184,19 @@ private:
// }
// ~CudaIntegrateBrownianStepKernel();
// /**
// * Initialize the kernel, setting up all parameters related to integrator.
// * Initialize the kernel.
// *
// * @param masses the mass of each atom
// * @param constraintIndices each element contains the indices of two atoms whose distance should be constrained
// * @param constraintLengths the required distance between each pair of constrained atoms
// * @param system the System this kernel will be applied to
// * @param integrator the BrownianIntegrator this kernel will be used for
// */
// void initialize(const std::vector<double>& masses, const std::vector<std::vector<int> >& constraintIndices,
// const std::vector<double>& constraintLengths);
// void initialize(const System& system, const BrownianIntegrator& integrator);
// /**
// * Execute the kernel.
// *
// * @param positions a Stream of type Double3 containing the position (x, y, z) of each atom
// * @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
// * @param forces a Stream of type Double3 containing the force (x, y, z) on each atom
// * @param temperature the temperature of the heat bath
// * @param friction the friction coefficient coupling the system to the heat bath
// * @param stepSize the integration step size
// * @param context the context in which to execute this kernel
// * @param integrator the BrownianIntegrator this kernel is being used for
// */
// void execute(Stream& positions, Stream& velocities, const Stream& forces, double temperature, double friction, double stepSize);
// void execute(OpenMMContextImpl& context, const BrownianIntegrator& integrator);
//private:
// CudaBrownianDynamics* dynamics;
// CudaShakeAlgorithm* shake;
......@@ -254,20 +216,18 @@ private:
// }
// ~CudaApplyAndersenThermostatKernel();
// /**
// * Initialize the kernel, setting up the values of unchanging parameters.
// * Initialize the kernel.
// *
// * @param masses the mass of each atom
// * @param system the System this kernel will be applied to
// * @param thermostat the AndersenThermostat this kernel will be used for
// */
// void initialize(const std::vector<double>& masses);
// void initialize(const System& system, const AndersenThermostat& thermostat);
// /**
// * Execute the kernel.
// *
// * @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
// * @param temperature the temperature of the heat bath
// * @param collisionFrequency the frequency at which atom collide with particles in the heat bath
// * @param stepSize the integration step size
// * @param context the context in which to execute this kernel
// */
// void execute(Stream& velocities, double temperature, double collisionFrequency, double stepSize);
// void execute(OpenMMContextImpl& context);
//private:
// CudaAndersenThermostat* thermostat;
// RealOpenMM* masses;
......@@ -281,18 +241,17 @@ public:
CudaCalcKineticEnergyKernel(std::string name, const Platform& platform) : CalcKineticEnergyKernel(name, platform) {
}
/**
* Initialize the kernel, setting up the atomic masses.
* Initialize the kernel.
*
* @param masses the mass of each atom
* @param system the System this kernel will be applied to
*/
void initialize(const std::vector<double>& masses);
void initialize(const System& system);
/**
* Execute the kernel.
*
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @return the kinetic energy of the system
* @param context the context in which to execute this kernel
*/
double execute(const Stream& velocities);
double execute(OpenMMContextImpl& context);
private:
std::vector<double> masses;
};
......@@ -307,15 +266,16 @@ public:
/**
* Initialize the kernel, setting up the atomic masses.
*
* @param masses the mass of each atom
* @param system the System this kernel will be applied to
* @param force the CMMotionRemover this kernel will be used for
*/
void initialize(const std::vector<double>& masses);
void initialize(const System& system, const CMMotionRemover& force);
/**
* Execute the kernel.
*
* @param velocities a Stream of type Double3 containing the velocity (x, y, z) of each atom
* @param context the context in which to execute this kernel
*/
void execute(Stream& velocities);
void execute(OpenMMContextImpl& context);
private:
CudaPlatform::PlatformData& data;
};
......
platforms/cuda/tests/TestCudaCMMotionRemover.cpp
View file @ a138663d
......@@ -63,7 +63,7 @@ void testMotionRemoval() {
CudaPlatform platform;
System system(numAtoms, 0);
LangevinIntegrator integrator(0.0, 1e-5, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 1, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 1, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, i+1);
forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
......
platforms/cuda/tests/TestCudaGBSAOBCForceField.cpp
View file @ a138663d
......@@ -58,7 +58,7 @@ void testSingleAtom() {
GBSAOBCForceField* forceField = new GBSAOBCForceField(1);
forceField->setAtomParameters(0, 0.5, 0.15, 1);
system.addForce(forceField);
system.addForce(new StandardMMForceField(1, 0, 0, 0, 0));
system.addForce(new StandardMMForceField(1, 0, 0, 0, 0, 0));
OpenMMContext context(system, integrator, platform);
vector<Vec3> positions(1);
positions[0] = Vec3(0, 0, 0);
......@@ -81,7 +81,7 @@ void testForce() {
for (int i = 0; i < numAtoms; ++i)
forceField->setAtomParameters(i, i%2 == 0 ? -1 : 1, 0.15, 1);
system.addForce(forceField);
system.addForce(new StandardMMForceField(numAtoms, 0, 0, 0, 0));
system.addForce(new StandardMMForceField(numAtoms, 0, 0, 0, 0, 0));
OpenMMContext context(system, integrator, platform);
// Set random positions for all the atoms.
......
platforms/cuda/tests/TestCudaLangevinIntegrator.cpp
View file @ a138663d
......@@ -55,7 +55,7 @@ void testSingleBond() {
system.setAtomMass(0, 2.0);
system.setAtomMass(1, 2.0);
LangevinIntegrator integrator(0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 1, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(2, 1, 0, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1.5, 1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
......@@ -99,7 +99,7 @@ void testTemperature() {
CudaPlatform platform;
System system(numAtoms, 0);
LangevinIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 2.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
......@@ -135,7 +135,7 @@ void testConstraints() {
CudaPlatform platform;
System system(numAtoms, numConstraints);
LangevinIntegrator integrator(temp, 2.0, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(numAtoms, 0, 0, 0, 0, 0);
for (int i = 0; i < numAtoms; ++i) {
system.setAtomMass(i, 10.0);
forceField->setAtomParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
......
platforms/cuda/tests/TestCudaStandardMMForceField.cpp
View file @ a138663d
......@@ -52,7 +52,7 @@ void testBonds() {
CudaPlatform platform;
System system(3, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 2, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(3, 2, 0, 0, 0, 0);
forceField->setBondParameters(0, 0, 1, 1.5, 0.8);
forceField->setBondParameters(1, 1, 2, 1.2, 0.7);
system.addForce(forceField);
......@@ -74,7 +74,7 @@ void testAngles() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 2, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 2, 0, 0, 0);
forceField->setAngleParameters(0, 0, 1, 2, PI_M/3, 1.1);
forceField->setAngleParameters(1, 1, 2, 3, PI_M/2, 1.2);
system.addForce(forceField);
......@@ -99,7 +99,7 @@ void testPeriodicTorsions() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 1, 0);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 1, 0, 0);
forceField->setPeriodicTorsionParameters(0, 0, 1, 2, 3, 2, PI_M/3, 1.1);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
......@@ -122,7 +122,7 @@ void testRBTorsions() {
CudaPlatform platform;
System system(4, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 0, 1);
StandardMMForceField* forceField = new StandardMMForceField(4, 0, 0, 0, 1, 0);
forceField->setRBTorsionParameters(0, 0, 1, 2, 3, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6);
system.addForce(forceField);
OpenMMContext context(system, integrator, platform);
......@@ -155,7 +155,7 @@ void testCoulomb() {
CudaPlatform platform;
System system(2, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0, 0);
forceField->setAtomParameters(0, 0.5, 1, 0);
forceField->setAtomParameters(1, -1.5, 1, 0);
system.addForce(forceField);
......@@ -176,7 +176,7 @@ void testLJ() {
CudaPlatform platform;
System system(2, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(2, 0, 0, 0, 0, 0);
forceField->setAtomParameters(0, 0, 1.2, 1);
forceField->setAtomParameters(1, 0, 1.4, 2);
system.addForce(forceField);
......@@ -199,7 +199,7 @@ void testExclusionsAnd14() {
CudaPlatform platform;
System system(5, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0, 2);
forceField->setBondParameters(0, 0, 1, 1, 0);
forceField->setBondParameters(1, 1, 2, 1, 0);
forceField->setBondParameters(2, 2, 3, 1, 0);
......@@ -217,6 +217,8 @@ void testExclusionsAnd14() {
}
forceField->setAtomParameters(0, 0, 1.5, 1);
forceField->setAtomParameters(i, 0, 1.5, 1);
forceField->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
positions[i] = Vec3(r, 0, 0);
OpenMMContext context(system, integrator, platform);
context.setPositions(positions);
......@@ -242,6 +244,8 @@ void testExclusionsAnd14() {
forceField->setAtomParameters(0, 2, 1.5, 0);
forceField->setAtomParameters(i, 2, 1.5, 0);
forceField->setNonbonded14Parameters(0, 0, 3, i == 3 ? 4/1.2 : 0, 1.5, 0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
OpenMMContext context2(system, integrator, platform);
context2.setPositions(positions);
state = context2.getState(State::Forces | State::Energy);
......@@ -266,7 +270,7 @@ void testCutoff() {
CudaPlatform platform;
System system(3, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 0, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(3, 0, 0, 0, 0, 0);
forceField->setAtomParameters(0, 1.0, 1, 0);
forceField->setAtomParameters(1, 1.0, 1, 0);
forceField->setAtomParameters(2, 1.0, 1, 0);
......@@ -299,7 +303,7 @@ void testCutoff14() {
CudaPlatform platform;
System system(5, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(5, 4, 0, 0, 0, 2);
forceField->setBondParameters(0, 0, 1, 1, 0);
forceField->setBondParameters(1, 1, 2, 1, 0);
forceField->setBondParameters(2, 2, 3, 1, 0);
......@@ -323,6 +327,8 @@ void testCutoff14() {
for (int j = 1; j < 5; ++j)
forceField->setAtomParameters(j, 0, 1.5, 0);
forceField->setAtomParameters(i, 0, 1.5, 1);
forceField->setNonbonded14Parameters(0, 0, 3, 0, 1.5, i == 3 ? 0.5 : 0.0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0);
context.reinitialize();
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
......@@ -349,6 +355,8 @@ void testCutoff14() {
const double q = 0.7;
forceField->setAtomParameters(0, q, 1.5, 0);
forceField->setAtomParameters(i, q, 1.5, 0);
forceField->setNonbonded14Parameters(0, 0, 3, i == 3 ? q*q/1.2 : 0, 1.5, 0);
forceField->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0);
context.reinitialize();
context.setPositions(positions);
state = context.getState(State::Forces | State::Energy);
......@@ -376,7 +384,7 @@ void testPeriodic() {
CudaPlatform platform;
System system(3, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
StandardMMForceField* forceField = new StandardMMForceField(3, 1, 0, 0, 0);
StandardMMForceField* forceField = new StandardMMForceField(3, 1, 0, 0, 0, 0);
forceField->setAtomParameters(0, 1.0, 1, 0);
forceField->setAtomParameters(1, 1.0, 1, 0);
forceField->setAtomParameters(2, 1.0, 1, 0);
......
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