Skip to content

GitLab

Commit 29e3fa57 authored by Peter Eastman's avatar Peter Eastman
Browse files

Redesigned the API for specifying exclusions and 1-4 interactions.

parent afd06645
Show whitespace changes
Inline Side-by-side
Showing with 316 additions and 195 deletions
CMakeLists.txt
View file @ 29e3fa57
...@@ -47,7 +47,7 @@ ENDIF(${CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT}) ...@@ -47,7 +47,7 @@ ENDIF(${CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT})
INCLUDE(Dart) INCLUDE(Dart)
# SUBDIRS (platforms/reference/tests platforms/cuda) SUBDIRS (tests)
# The source is organized into subdirectories, but we handle them all from # The source is organized into subdirectories, but we handle them all from
# this CMakeLists file rather than letting CMake visit them as SUBDIRS. # this CMakeLists file rather than letting CMake visit them as SUBDIRS.
......
olla/include/kernels.h
View file @ 29e3fa57
...@@ -226,11 +226,8 @@ public: ...@@ -226,11 +226,8 @@ public:
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the NonbondedForce this kernel will be used for * @param force the NonbondedForce this kernel will be used for
* @param exclusions the i'th element lists the indices of all particles with which the i'th particle 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.
*/ */
virtual void initialize(const System& system, const NonbondedForce& force, const std::vector<std::set<int> >& exclusions) = 0; virtual void initialize(const System& system, const NonbondedForce& force) = 0;
/** /**
* Execute the kernel to calculate the forces. * Execute the kernel to calculate the forces.
* *
......
openmmapi/include/NonbondedForce.h
View file @ 29e3fa57
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008 Stanford University and the Authors. * * Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -35,6 +35,7 @@ ...@@ -35,6 +35,7 @@
#include "Force.h" #include "Force.h"
#include "Vec3.h" #include "Vec3.h"
#include <map> #include <map>
#include <set>
#include <vector> #include <vector>
#include "internal/windowsExport.h" #include "internal/windowsExport.h"
...@@ -47,15 +48,21 @@ namespace OpenMM { ...@@ -47,15 +48,21 @@ namespace OpenMM {
* calculated with the Lorentz-Bertelot combining rule: it uses the arithmetic mean of the sigmas and the * calculated with the Lorentz-Bertelot combining rule: it uses the arithmetic mean of the sigmas and the
* geometric mean of the epsilons for the two interacting particles. * geometric mean of the epsilons for the two interacting particles.
* *
* If the System also contains a HarmonicBondForce, nonbonded interactions are automatically excluded between * To use this class, create a NonbondedForce object, then call addParticle() once for each particle in the
* particles which are separated by three or fewer bonds. Most molecular force fields omit Coulomb and * System to define its parameters. The number of particles for which you define nonbonded parameters must
* Lennard-Jones interactions between particles separated by one or two bonds, while using modified parameters * be exactly equal to the number of particles in the System, or else an exception will be thrown when you
* for those separated by three bonds (known as "1-4 interactions"). This class lets you provide a list of * try to create an OpenMMContext. After a particle has been added, you can modify its force field parameters
* 1-4 interactions to include in the potential, along with the parameters to use for each one. * by calling setParticleParameters().
* *
* When creating a NonbondedForce, you specify the number of particles and the number of 1-4 interactions as * NonbondedForce also lets you specify "exceptions", particular pairs of particles whose interactions should be
* arguments to the constructor. You then loop over them and call setParticleParameters() for each particle and * computed based on different parameters than those defined for the individual particles. This can be used to
* setNonbond14Parameters() for each 1-4 interaction. * completely exclude certain interactions from the force calculation, or to alter how they interact with each other.
*
* Many molecular force fields omit Coulomb and Lennard-Jones interactions between particles separated by one
* or two bonds, while using modified parameters for those separated by three bonds (known as "1-4 interactions").
* This class provides a convenience method for this case called createExceptionsFromBonds(). You pass to it
* a list of bonds and the scale factors to use for 1-4 interactions. It identifies all pairs of particles which
* are separated by 1, 2, or 3 bonds, then automatically creates exceptions for them.
*/ */
class OPENMM_EXPORT NonbondedForce : public Force { class OPENMM_EXPORT NonbondedForce : public Force {
...@@ -81,28 +88,26 @@ public: ...@@ -81,28 +88,26 @@ public:
*/ */
CutoffPeriodic = 2, CutoffPeriodic = 2,
/** /**
* Using Ewald summation method for long-range electrostatics. * Periodic boundary conditions are used, and Ewald summation is used to compute the interaction of each particle
* with all periodic copies of every other particle.
*/ */
Ewald = 3 Ewald = 3
}; };
/** /**
* Create a NonbondedForce. * Create a NonbondedForce.
*
* @param numParticles the number of particles in the system
* @param numNonbonded14 the number of nonbonded 1-4 terms in the potential function
*/ */
NonbondedForce(int numParticles, int numNonbonded14); NonbondedForce();
/** /**
* Get the number of particles in the system. * Get the number of particles for which force field parameters have been defined.
*/ */
int getNumParticles() const { int getNumParticles() const {
return particles.size(); return particles.size();
} }
/** /**
* Get the number of nonbonded 1-4 terms in the potential function * Get the number of special interactions that should be calculated differently from other interactions.
*/ */
int getNumNonbonded14() const { int getNumExceptions() const {
return nb14s.size(); return exceptions.size();
} }
/** /**
* Get the method used for handling long range nonbonded interactions. * Get the method used for handling long range nonbonded interactions.
...@@ -146,13 +151,25 @@ public: ...@@ -146,13 +151,25 @@ public:
* @param c the vector defining the third edge of the periodic box * @param c the vector defining the third edge of the periodic box
*/ */
void setPeriodicBoxVectors(Vec3 a, Vec3 b, Vec3 c); void setPeriodicBoxVectors(Vec3 a, Vec3 b, Vec3 c);
/**
* Add the nonbonded force parameters for a particle. This should be called once for each particle
* in the System. When it is called for the i'th time, it specifies the parameters for the i'th particle.
* For calculating the Lennard-Jones interaction between two particles, the arithmetic mean of the sigmas
* and the geometric mean of the epsilons for the two interacting particles is used (the Lorentz-Bertelot
* combining rule).
*
* @param charge the charge of the particle, measured in units of the proton charge
* @param sigma the sigma parameter of the Lennard-Jones potential (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the epsilon parameter of the Lennard-Jones potential (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
*/
void addParticle(double charge, double sigma, double epsilon);
/** /**
* Get the nonbonded force parameters for a particle. * Get the nonbonded force parameters for a particle.
* *
* @param index the index of the particle for which to get parameters * @param index the index of the particle for which to get parameters
* @param charge the charge of the particle, measured in units of the proton charge * @param charge the charge of the particle, measured in units of the proton charge
* @param sigma the van der Waals radius of the particle (sigma in the Lennard Jones potential), measured in nm * @param sigma the sigma parameter of the Lennard-Jones potential (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the well depth of the van der Waals interaction (epsilon in the Lennard Jones potential), measured in kJ/mol * @param epsilon the epsilon parameter of the Lennard-Jones potential (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
*/ */
void getParticleParameters(int index, double& charge, double& sigma, double& epsilon) const; void getParticleParameters(int index, double& charge, double& sigma, double& epsilon) const;
/** /**
...@@ -162,40 +179,71 @@ public: ...@@ -162,40 +179,71 @@ public:
* *
* @param index the index of the particle for which to set parameters * @param index the index of the particle for which to set parameters
* @param charge the charge of the particle, measured in units of the proton charge * @param charge the charge of the particle, measured in units of the proton charge
* @param sigma the van der Waals radius of the particle (sigma in the Lennard Jones potential), measured in nm * @param sigma the sigma parameter of the Lennard-Jones potential (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the well depth of the van der Waals interaction (epsilon in the Lennard Jones potential), measured in kJ/mol * @param epsilon the epsilon parameter of the Lennard-Jones potential (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
*/ */
void setParticleParameters(int index, double charge, double sigma, double epsilon); void setParticleParameters(int index, double charge, double sigma, double epsilon);
/** /**
* Get the force field parameters for a nonbonded 1-4 term. * Add an interaction to the list of exceptions that should be calculated differently from other interactions.
* If chargeProd and epsilon are both equal to 0, this will cause the interaction to be completely omitted from
* force and energy calculations.
*
* In many cases, you can use createExceptionsFromBonds() rather than adding each exception explicitly.
*
* @param particle1 the index of the first particle involved in the interaction
* @param particle2 the index of the second particle involved in the interaction
* @param chargeProd the scaled product of the atomic charges (i.e. the strength of the Coulomb interaction), measured in units of the proton charge squared
* @param sigma the sigma parameter of the Lennard-Jones potential (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the epsilon parameter of the Lennard-Jones potential (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
*/
void addException(int particle1, int particle2, double chargeProd, double sigma, double epsilon);
/**
* Get the force field parameters for an interaction that should be calculated differently from others.
* *
* @param index the index of the interaction for which to get parameters * @param index the index of the interaction for which to get parameters
* @param particle1 the index of the first particle involved in the interaction * @param particle1 the index of the first particle involved in the interaction
* @param particle2 the index of the second particle involved in the interaction * @param particle2 the index of the second particle involved in the interaction
* @param chargeProd the scaled product of the atomic charges (i.e. the strength of the Coulomb interaction), measured in units of the proton charge squared * @param chargeProd the scaled product of the atomic charges (i.e. the strength of the Coulomb interaction), measured in units of the proton charge squared
* @param sigma the van der Waals radius of the particle (sigma in the Lennard Jones potential), measured in nm * @param sigma the sigma parameter of the Lennard-Jones potential (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the well depth of the van der Waals interaction (epsilon in the Lennard Jones potential), measured in kJ/mol * @param epsilon the epsilon parameter of the Lennard-Jones potential (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
*/ */
void getNonbonded14Parameters(int index, int& particle1, int& particle2, double& chargeProd, double& sigma, double& epsilon) const; void getExceptionParameters(int index, int& particle1, int& particle2, double& chargeProd, double& sigma, double& epsilon) const;
/** /**
* Set the force field parameters for a nonbonded 1-4 term. * Set the force field parameters for an interaction that should be calculated differently from others.
* If chargeProd and epsilon are both equal to 0, this will cause the interaction to be completely omitted from
* force and energy calculations.
* *
* @param index the index of the interaction for which to get parameters * @param index the index of the interaction for which to get parameters
* @param particle1 the index of the first particle involved in the interaction * @param particle1 the index of the first particle involved in the interaction
* @param particle2 the index of the second particle involved in the interaction * @param particle2 the index of the second particle involved in the interaction
* @param chargeProd the scaled product of the atomic charges (i.e. the strength of the Coulomb interaction), measured in units of the proton charge squared * @param chargeProd the scaled product of the atomic charges (i.e. the strength of the Coulomb interaction), measured in units of the proton charge squared
* @param sigma the van der Waals radius of the particle (sigma in the Lennard Jones potential), measured in nm * @param sigma the sigma parameter of the Lennard-Jones potential (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the well depth of the van der Waals interaction (epsilon in the Lennard Jones potential), measured in kJ/mol * @param epsilon the epsilon parameter of the Lennard-Jones potential (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
*/ */
void setNonbonded14Parameters(int index, int particle1, int particle2, double chargeProd, double sigma, double epsilon); void setExceptionParameters(int index, int particle1, int particle2, double chargeProd, double sigma, double epsilon);
/**
* Identify exceptions based on the molecular topology. Particles which are separated by one or two bonds are set
* to not interact at all, while pairs of particles separated by three bonds (known as "1-4 interactions") have
* their Coulomb and Lennard-Jones interactions reduced by a fixed factor.
*
* @param bonds the set of bonds based on which to construct exceptions. Each element specifies the indices of
* two particles that are bonded to each other.
* @param coulomb14Scale pairs of particles separated by three bonds will have the strength of their Coulomb interaction
* multiplied by this factor
* @param lj14Scale pairs of particles separated by three bonds will have the strength of their Lennard-Jones interaction
* multiplied by this factor
*/
void createExceptionsFromBonds(const std::vector<std::pair<int, int> >& bonds, double coulomb14Scale, double lj14Scale);
protected: protected:
ForceImpl* createImpl(); ForceImpl* createImpl();
private: private:
class ParticleInfo; class ParticleInfo;
class NB14Info; class ExceptionInfo;
NonbondedMethod nonbondedMethod; NonbondedMethod nonbondedMethod;
double cutoffDistance; double cutoffDistance;
Vec3 periodicBoxVectors[3]; Vec3 periodicBoxVectors[3];
void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
// Retarded visual studio compiler complains about being unable to // Retarded visual studio compiler complains about being unable to
// export private stl class members. // export private stl class members.
...@@ -205,7 +253,7 @@ private: ...@@ -205,7 +253,7 @@ private:
#pragma warning(disable:4251) #pragma warning(disable:4251)
#endif #endif
std::vector<ParticleInfo> particles; std::vector<ParticleInfo> particles;
std::vector<NB14Info> nb14s; std::vector<ExceptionInfo> exceptions;
#if defined(_MSC_VER) #if defined(_MSC_VER)
#pragma warning(pop) #pragma warning(pop)
#endif #endif
...@@ -218,16 +266,22 @@ public: ...@@ -218,16 +266,22 @@ public:
ParticleInfo() { ParticleInfo() {
charge = sigma = epsilon = 0.0; charge = sigma = epsilon = 0.0;
} }
ParticleInfo(double charge, double sigma, double epsilon) :
charge(charge), sigma(sigma), epsilon(epsilon) {
}
}; };
class NonbondedForce::NB14Info { class NonbondedForce::ExceptionInfo {
public: public:
int particle1, particle2; int particle1, particle2;
double chargeProd, sigma, epsilon; double chargeProd, sigma, epsilon;
NB14Info() { ExceptionInfo() {
particle1 = particle2 = -1; particle1 = particle2 = -1;
chargeProd = sigma = epsilon = 0.0; chargeProd = sigma = epsilon = 0.0;
} }
ExceptionInfo(int particle1, int particle2, double chargeProd, double sigma, double epsilon) :
particle1(particle1), particle2(particle2), chargeProd(chargeProd), sigma(sigma), epsilon(epsilon) {
}
}; };
} // namespace OpenMM } // namespace OpenMM
......
openmmapi/include/internal/NonbondedForceImpl.h
View file @ 29e3fa57
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008 Stanford University and the Authors. * * Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -63,8 +63,6 @@ public: ...@@ -63,8 +63,6 @@ public:
} }
std::vector<std::string> getKernelNames(); std::vector<std::string> getKernelNames();
private: private:
void findExclusions(const std::vector<std::vector<int> >& bondIndices, std::vector<std::set<int> >& exclusions, std::set<std::pair<int, int> >& bonded14Indices) const;
void addExclusionsToSet(const std::vector<std::set<int> >& bonded12, std::set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const;
NonbondedForce& owner; NonbondedForce& owner;
Kernel kernel; Kernel kernel;
}; };
......
openmmapi/src/NonbondedForce.cpp
View file @ 29e3fa57
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008 Stanford University and the Authors. * * Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -33,11 +33,15 @@ ...@@ -33,11 +33,15 @@
#include "OpenMMException.h" #include "OpenMMException.h"
#include "NonbondedForce.h" #include "NonbondedForce.h"
#include "internal/NonbondedForceImpl.h" #include "internal/NonbondedForceImpl.h"
#include <cmath>
#include <utility>
using namespace OpenMM; using namespace OpenMM;
using std::pair;
using std::set;
using std::vector;
NonbondedForce::NonbondedForce(int numParticles, int numNonbonded14) : particles(numParticles), nb14s(numNonbonded14), NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
nonbondedMethod(NoCutoff), cutoffDistance(1.0) {
periodicBoxVectors[0] = Vec3(2, 0, 0); periodicBoxVectors[0] = Vec3(2, 0, 0);
periodicBoxVectors[1] = Vec3(0, 2, 0); periodicBoxVectors[1] = Vec3(0, 2, 0);
periodicBoxVectors[2] = Vec3(0, 0, 2); periodicBoxVectors[2] = Vec3(0, 0, 2);
...@@ -77,6 +81,10 @@ void NonbondedForce::setPeriodicBoxVectors(Vec3 a, Vec3 b, Vec3 c) { ...@@ -77,6 +81,10 @@ void NonbondedForce::setPeriodicBoxVectors(Vec3 a, Vec3 b, Vec3 c) {
periodicBoxVectors[2] = c; periodicBoxVectors[2] = c;
} }
void NonbondedForce::addParticle(double charge, double sigma, double epsilon) {
particles.push_back(ParticleInfo(charge, sigma, epsilon));
}
void NonbondedForce::getParticleParameters(int index, double& charge, double& sigma, double& epsilon) const { void NonbondedForce::getParticleParameters(int index, double& charge, double& sigma, double& epsilon) const {
charge = particles[index].charge; charge = particles[index].charge;
sigma = particles[index].sigma; sigma = particles[index].sigma;
...@@ -89,22 +97,73 @@ void NonbondedForce::setParticleParameters(int index, double charge, double sigm ...@@ -89,22 +97,73 @@ void NonbondedForce::setParticleParameters(int index, double charge, double sigm
particles[index].epsilon = epsilon; particles[index].epsilon = epsilon;
} }
void NonbondedForce::getNonbonded14Parameters(int index, int& particle1, int& particle2, double& chargeProd, double& sigma, double& epsilon) const { void NonbondedForce::addException(int particle1, int particle2, double chargeProd, double sigma, double epsilon) {
particle1 = nb14s[index].particle1; exceptions.push_back(ExceptionInfo(particle1, particle2, chargeProd, sigma, epsilon));
particle2 = nb14s[index].particle2; }
chargeProd = nb14s[index].chargeProd; void NonbondedForce::getExceptionParameters(int index, int& particle1, int& particle2, double& chargeProd, double& sigma, double& epsilon) const {
sigma = nb14s[index].sigma; particle1 = exceptions[index].particle1;
epsilon = nb14s[index].epsilon; particle2 = exceptions[index].particle2;
chargeProd = exceptions[index].chargeProd;
sigma = exceptions[index].sigma;
epsilon = exceptions[index].epsilon;
} }
void NonbondedForce::setNonbonded14Parameters(int index, int particle1, int particle2, double chargeProd, double sigma, double epsilon) { void NonbondedForce::setExceptionParameters(int index, int particle1, int particle2, double chargeProd, double sigma, double epsilon) {
nb14s[index].particle1 = particle1; exceptions[index].particle1 = particle1;
nb14s[index].particle2 = particle2; exceptions[index].particle2 = particle2;
nb14s[index].chargeProd = chargeProd; exceptions[index].chargeProd = chargeProd;
nb14s[index].sigma = sigma; exceptions[index].sigma = sigma;
nb14s[index].epsilon = epsilon; exceptions[index].epsilon = epsilon;
} }
ForceImpl* NonbondedForce::createImpl() { ForceImpl* NonbondedForce::createImpl() {
return new NonbondedForceImpl(*this); return new NonbondedForceImpl(*this);
} }
void NonbondedForce::createExceptionsFromBonds(const vector<pair<int, int> >& bonds, double coulomb14Scale, double lj14Scale) {
// Find particles separated by 1, 2, or 3 bonds.
vector<set<int> > exclusions(particles.size());
vector<set<int> > bonded12(exclusions.size());
for (int i = 0; i < (int) bonds.size(); ++i) {
bonded12[bonds[i].first].insert(bonds[i].second);
bonded12[bonds[i].second].insert(bonds[i].first);
}
for (int i = 0; i < (int) exclusions.size(); ++i)
addExclusionsToSet(bonded12, exclusions[i], i, i, 2);
// Find particles separated by 1 or 2 bonds and create the exceptions.
for (int i = 0; i < (int) exclusions.size(); ++i) {
set<int> bonded13;
addExclusionsToSet(bonded12, bonded13, i, i, 1);
for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
if (*iter < i) {
if (bonded13.find(*iter) == bonded13.end()) {
// This is a 1-4 interaction.
const ParticleInfo& particle1 = particles[*iter];
const ParticleInfo& particle2 = particles[i];
const double chargeProd = coulomb14Scale*particle1.charge*particle2.charge;
const double sigma = 0.5*(particle1.sigma+particle2.sigma);
const double epsilon = lj14Scale*std::sqrt(particle1.epsilon*particle2.epsilon);
addException(*iter, i, chargeProd, sigma, epsilon);
}
else {
// This interaction should be completely excluded.
addException(*iter, i, 0.0, 1.0, 0.0);
}
}
}
}
void NonbondedForce::addExclusionsToSet(const vector<set<int> >& bonded12, set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const {
for (set<int>::const_iterator iter = bonded12[fromParticle].begin(); iter != bonded12[fromParticle].end(); ++iter) {
if (*iter != baseParticle)
exclusions.insert(*iter);
if (currentLevel > 0)
addExclusionsToSet(bonded12, exclusions, baseParticle, *iter, currentLevel-1);
}
}
openmmapi/src/NonbondedForceImpl.cpp
View file @ 29e3fa57
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008 Stanford University and the Authors. * * Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -29,14 +29,17 @@ ...@@ -29,14 +29,17 @@
* USE OR OTHER DEALINGS IN THE SOFTWARE. * * USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */ * -------------------------------------------------------------------------- */
#include "OpenMMException.h"
#include "internal/OpenMMContextImpl.h" #include "internal/OpenMMContextImpl.h"
#include "internal/NonbondedForceImpl.h" #include "internal/NonbondedForceImpl.h"
#include "kernels.h" #include "kernels.h"
#include <sstream>
using namespace OpenMM; using namespace OpenMM;
using std::pair; using std::pair;
using std::vector; using std::vector;
using std::set; using std::set;
using std::stringstream;
NonbondedForceImpl::NonbondedForceImpl(NonbondedForce& owner) : owner(owner) { NonbondedForceImpl::NonbondedForceImpl(NonbondedForce& owner) : owner(owner) {
} }
...@@ -47,40 +50,40 @@ NonbondedForceImpl::~NonbondedForceImpl() { ...@@ -47,40 +50,40 @@ NonbondedForceImpl::~NonbondedForceImpl() {
void NonbondedForceImpl::initialize(OpenMMContextImpl& context) { void NonbondedForceImpl::initialize(OpenMMContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcNonbondedForceKernel::Name(), context); kernel = context.getPlatform().createKernel(CalcNonbondedForceKernel::Name(), context);
// See if the system contains a HarmonicBondForce. If so, use it to identify exclusions. // Check for errors in the specification of exceptions.
System& system = context.getSystem(); System& system = context.getSystem();
vector<set<int> > exclusions(owner.getNumParticles()); if (owner.getNumParticles() != system.getNumParticles())
vector<vector<int> > bondIndices; throw OpenMMException("NonbondedForce must have exactly as many particles as the System it belongs to.");
set<pair<int, int> > bonded14set; vector<set<int> > exceptions(owner.getNumParticles());
for (int i = 0; i < system.getNumForces(); i++) { for (int i = 0; i < owner.getNumExceptions(); i++) {
if (dynamic_cast<HarmonicBondForce*>(&system.getForce(i)) != NULL) {
const HarmonicBondForce& force = dynamic_cast<const HarmonicBondForce&>(system.getForce(i));
bondIndices.resize(force.getNumBonds());
for (int j = 0; j < force.getNumBonds(); ++j) {
int particle1, particle2; int particle1, particle2;
double length, k; double chargeProd, sigma, epsilon;
force.getBondParameters(j, particle1, particle2, length, k); owner.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
bondIndices[j].push_back(particle1); if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
bondIndices[j].push_back(particle2); stringstream msg;
msg << "NonbondedForce: Illegal particle index for an exception: ";
msg << particle1;
throw OpenMMException(msg.str());
} }
break; if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
stringstream msg;
msg << "NonbondedForce: Illegal particle index for an exception: ";
msg << particle2;
throw OpenMMException(msg.str());
} }
if (exceptions[particle1].count(particle2) > 0 || exceptions[particle2].count(particle1) > 0) {
stringstream msg;
msg << "NonbondedForce: Multiple exceptions are specified for particles ";
msg << particle1;
msg << " and ";
msg << particle2;
throw OpenMMException(msg.str());
} }
exceptions[particle1].insert(particle2);
// Also treat constrained distances as bonds. exceptions[particle2].insert(particle1);
int numBonds = bondIndices.size();
bondIndices.resize(numBonds+system.getNumConstraints());
for (int j = 0; j < system.getNumConstraints(); j++) {
int particle1, particle2;
double distance;
system.getConstraintParameters(j, particle1, particle2, distance);
bondIndices[j+numBonds].push_back(particle1);
bondIndices[j+numBonds].push_back(particle2);
} }
findExclusions(bondIndices, exclusions, bonded14set); dynamic_cast<CalcNonbondedForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
dynamic_cast<CalcNonbondedForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner, exclusions);
} }
void NonbondedForceImpl::calcForces(OpenMMContextImpl& context, Stream& forces) { void NonbondedForceImpl::calcForces(OpenMMContextImpl& context, Stream& forces) {
...@@ -96,30 +99,3 @@ std::vector<std::string> NonbondedForceImpl::getKernelNames() { ...@@ -96,30 +99,3 @@ std::vector<std::string> NonbondedForceImpl::getKernelNames() {
names.push_back(CalcNonbondedForceKernel::Name()); names.push_back(CalcNonbondedForceKernel::Name());
return names; return names;
} }
void NonbondedForceImpl::findExclusions(const vector<vector<int> >& bondIndices, vector<set<int> >& exclusions, set<pair<int, int> >& bonded14Indices) const {
vector<set<int> > bonded12(exclusions.size());
for (int i = 0; i < (int) bondIndices.size(); ++i) {
bonded12[bondIndices[i][0]].insert(bondIndices[i][1]);
bonded12[bondIndices[i][1]].insert(bondIndices[i][0]);
}
for (int i = 0; i < (int) exclusions.size(); ++i)
addExclusionsToSet(bonded12, exclusions[i], i, i, 2);
for (int i = 0; i < (int) exclusions.size(); ++i) {
set<int> bonded13;
addExclusionsToSet(bonded12, bonded13, i, i, 1);
for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter)
if (*iter < i && bonded13.find(*iter) == bonded13.end())
bonded14Indices.insert(pair<int, int> (*iter, i));
}
}
void NonbondedForceImpl::addExclusionsToSet(const vector<set<int> >& bonded12, set<int>& exclusions, int baseParticle, int fromParticle, int currentLevel) const {
for (set<int>::const_iterator iter = bonded12[fromParticle].begin(); iter != bonded12[fromParticle].end(); ++iter) {
if (*iter != baseParticle)
exclusions.insert(*iter);
if (currentLevel > 0)
addExclusionsToSet(bonded12, exclusions, baseParticle, *iter, currentLevel-1);
}
}
platforms/brook/src/BrookCalcNonbondedForceKernel.cpp
View file @ 29e3fa57
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008 Stanford University and the Authors. * * Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Peter Eastman, Mark Friedrichs * * Authors: Peter Eastman, Mark Friedrichs *
* Contributors: * * Contributors: *
* * * *
...@@ -129,7 +129,7 @@ int BrookCalcNonbondedForceKernel::setLog( FILE* log ){ ...@@ -129,7 +129,7 @@ int BrookCalcNonbondedForceKernel::setLog( FILE* log ){
* *
*/ */
void BrookCalcNonbondedForceKernel::initialize( const System& system, const NonbondedForce& force, const std::vector<std::set<int> >& exclusions ){ void BrookCalcNonbondedForceKernel::initialize( const System& system, const NonbondedForce& force ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -174,6 +174,21 @@ void BrookCalcNonbondedForceKernel::initialize( const System& system, const Nonb ...@@ -174,6 +174,21 @@ void BrookCalcNonbondedForceKernel::initialize( const System& system, const Nonb
nonbondedParameters[ii].push_back( depth ); nonbondedParameters[ii].push_back( depth );
} }
// Go through the exclusions.
std::vector<std::set<int> >& exclusions(_numberOfParticles);
std::vector<int> nb14s;
for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
exclusions[particle1].insert(particle2);
exclusions[particle2].insert(particle1);
exclusions.push_back(pair<int, int>(particle1, particle2));
if (chargeProd != 0.0 || epsilon != 0.0)
nb14s.push_back(i);
}
brookNonBonded.setup( _numberOfParticles, nonbondedParameters, exclusions, getPlatform() ); brookNonBonded.setup( _numberOfParticles, nonbondedParameters, exclusions, getPlatform() );
_openMMBrookInterface.setTriggerForceKernel( this ); _openMMBrookInterface.setTriggerForceKernel( this );
_openMMBrookInterface.setTriggerEnergyKernel( this ); _openMMBrookInterface.setTriggerEnergyKernel( this );
...@@ -188,7 +203,7 @@ void BrookCalcNonbondedForceKernel::initialize( const System& system, const Nonb ...@@ -188,7 +203,7 @@ void BrookCalcNonbondedForceKernel::initialize( const System& system, const Nonb
// nonbonded 14 ixns // nonbonded 14 ixns
initialize14Interactions( system, force ); initialize14Interactions( system, force, nb14s );
} }
...@@ -200,7 +215,7 @@ void BrookCalcNonbondedForceKernel::initialize( const System& system, const Nonb ...@@ -200,7 +215,7 @@ void BrookCalcNonbondedForceKernel::initialize( const System& system, const Nonb
* *
*/ */
void BrookCalcNonbondedForceKernel::initialize14Interactions( const System& system, const NonbondedForce& force ){ void BrookCalcNonbondedForceKernel::initialize14Interactions( const System& system, const NonbondedForce& force, const std::vector<int>& nb14s ){
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -216,7 +231,7 @@ void BrookCalcNonbondedForceKernel::initialize14Interactions( const System& syst ...@@ -216,7 +231,7 @@ void BrookCalcNonbondedForceKernel::initialize14Interactions( const System& syst
// create _brookBondParameters object containing particle indices/parameters // create _brookBondParameters object containing particle indices/parameters
int numberOf14Forces = force.getNumNonbonded14(); int numberOf14Forces = nb14s.size();
if( numberOf14Forces > 0 ){ if( numberOf14Forces > 0 ){
_brookBondParameters = new BrookBondParameters( BondName, NumberOfParticlesInBond, NumberOfParametersInBond, numberOf14Forces, getLog() ); _brookBondParameters = new BrookBondParameters( BondName, NumberOfParticlesInBond, NumberOfParametersInBond, numberOf14Forces, getLog() );
...@@ -229,7 +244,7 @@ void BrookCalcNonbondedForceKernel::initialize14Interactions( const System& syst ...@@ -229,7 +244,7 @@ void BrookCalcNonbondedForceKernel::initialize14Interactions( const System& syst
int particles[NumberOfParticlesInBond]; int particles[NumberOfParticlesInBond];
double parameters[NumberOfParametersInBond]; double parameters[NumberOfParametersInBond];
force.getNonbonded14Parameters( ii, particle1, particle2, charge, radius, depth ); force.getExceptionParameters( nb14s[ii], particle1, particle2, charge, radius, depth );
//(void) fprintf( log, "%s idx=%d [%d %d] [%f %f %f]\n", methodName.c_str(), ii, particle1, particle2, charge, radius, depth ); //(void) fprintf( log, "%s idx=%d [%d %d] [%f %f %f]\n", methodName.c_str(), ii, particle1, particle2, charge, radius, depth );
particles[0] = particle1; particles[0] = particle1;
......
platforms/brook/src/BrookCalcNonbondedForceKernel.h
View file @ 29e3fa57
...@@ -54,12 +54,9 @@ class BrookCalcNonbondedForceKernel : public CalcNonbondedForceKernel { ...@@ -54,12 +54,9 @@ class BrookCalcNonbondedForceKernel : public CalcNonbondedForceKernel {
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the NonbondedForce this kernel will be used for * @param force the NonbondedForce this kernel will be used for
* @param exclusions the i'th element lists the indices of all particles with which the i'th particle 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.
*/ */
void initialize( const System& system, const NonbondedForce& force, const std::vector<std::set<int> >& exclusions ); void initialize( const System& system, const NonbondedForce& force );
/** /**
* Initialize the 14 ixns * Initialize the 14 ixns
......
platforms/brook/tests/TestBrookCMMotionRemover.cpp
View file @ 29e3fa57
...@@ -88,10 +88,10 @@ void testMotionRemoval( FILE* log ) { ...@@ -88,10 +88,10 @@ void testMotionRemoval( FILE* log ) {
bonds->setBondParameters(0, 2, 3, 2.0, 0.5); bonds->setBondParameters(0, 2, 3, 2.0, 0.5);
system.addForce(bonds); system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(numberOfParticles, 0); NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < numberOfParticles; ++i) { for (int i = 0; i < numberOfParticles; ++i) {
system.setParticleMass(i, (double) (i+1) ); system.setParticleMass(i, (double) (i+1) );
nonbonded->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0); nonbonded->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
} }
system.addForce(nonbonded); system.addForce(nonbonded);
......
platforms/brook/tests/TestBrookLangevinIntegrator.cpp
View file @ 29e3fa57
...@@ -217,10 +217,10 @@ void testLangevinTemperature( FILE* log ){ ...@@ -217,10 +217,10 @@ void testLangevinTemperature( FILE* log ){
System system(numberOfParticles, 0); System system(numberOfParticles, 0);
LangevinIntegrator integrator(temp, 0.2, 0.002); LangevinIntegrator integrator(temp, 0.2, 0.002);
NonbondedForce* forceField = new NonbondedForce(numberOfParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numberOfParticles; ++i){ for (int i = 0; i < numberOfParticles; ++i){
system.setParticleMass(i, mass ); system.setParticleMass(i, mass );
forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0); forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
} }
system.addForce(forceField); system.addForce(forceField);
...@@ -315,10 +315,10 @@ void testLangevinConstraints( FILE* log ){ ...@@ -315,10 +315,10 @@ void testLangevinConstraints( FILE* log ){
System system( numParticles, numConstraints ); System system( numParticles, numConstraints );
LangevinIntegrator integrator( temp, 2.0, 0.001 ); LangevinIntegrator integrator( temp, 2.0, 0.001 );
integrator.setConstraintTolerance(1e-5); integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, mass); system.setParticleMass(i, mass);
forceField->setParticleParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0); forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
} }
for (int i = 0; i < numConstraints; ++i){ for (int i = 0; i < numConstraints; ++i){
system.setConstraintParameters(i, 2*i, 2*i+1, 1.0); system.setConstraintParameters(i, 2*i, 2*i+1, 1.0);
......
platforms/brook/tests/TestBrookNonbondedForce.cpp
View file @ 29e3fa57
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2008 Stanford University and the Authors. * * Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Authors: Mark Friedrichs * * Authors: Mark Friedrichs *
* Contributors: * * Contributors: *
* * * *
...@@ -73,9 +73,9 @@ void testBrookCoulomb( FILE* log ){ ...@@ -73,9 +73,9 @@ void testBrookCoulomb( FILE* log ){
// int index, double charge, double radius, double depth // int index, double charge, double radius, double depth
NonbondedForce* forceField = new NonbondedForce( numberOfParticles, 0 ); NonbondedForce* forceField = new NonbondedForce();
forceField->setParticleParameters(0, 0.5, 1, 0); forceField->addParticle(0.5, 1, 0);
forceField->setParticleParameters(1, -1.5, 1, 0); forceField->addParticle(-1.5, 1, 0);
system.addForce(forceField); system.addForce(forceField);
//(void) fprintf( log, "%s: Calling context\n", methodName.c_str() ); //(void) fprintf( log, "%s: Calling context\n", methodName.c_str() );
...@@ -141,9 +141,9 @@ void testBrookLJ( FILE* log ){ ...@@ -141,9 +141,9 @@ void testBrookLJ( FILE* log ){
// int index, double charge, double radius, double depth // int index, double charge, double radius, double depth
NonbondedForce* forceField = new NonbondedForce( numberOfParticles, 0 ); NonbondedForce* forceField = new NonbondedForce();
forceField->setParticleParameters(0, 0, 1.2, 1); forceField->addParticle(0, 1.2, 1);
forceField->setParticleParameters(1, 0, 1.4, 2); forceField->addParticle(0, 1.4, 2);
system.addForce(forceField); system.addForce(forceField);
//(void) fprintf( log, "%s: Calling context\n", methodName.c_str() ); //(void) fprintf( log, "%s: Calling context\n", methodName.c_str() );
...@@ -208,13 +208,25 @@ void testBrookExclusionsAnd14( FILE* log ){ ...@@ -208,13 +208,25 @@ void testBrookExclusionsAnd14( FILE* log ){
// int index, double charge, double radius, double depth // int index, double charge, double radius, double depth
HarmonicBondForce* bonds = new HarmonicBondForce(4); NonbondedForce* nonbonded = new NonbondedForce();
bonds->setBondParameters(0, 0, 1, 1, 0); for (int i = 0; i < numberOfParticles; i++)
bonds->setBondParameters(1, 1, 2, 1, 0); nonbonded->addParticle(0, 1.5, 0)
bonds->setBondParameters(2, 2, 3, 1, 0); vector<pair<int, int> > bonds;
bonds->setBondParameters(3, 3, 4, 1, 0); bonds.push_back(pair<int, int>(0, 1));
system.addForce(bonds); bonds.push_back(pair<int, int>(1, 2));
NonbondedForce* nonbonded = new NonbondedForce( numberOfParticles, 2); bonds.push_back(pair<int, int>(2, 3));
bonds.push_back(pair<int, int>(3, 4));
nonbonded->createExceptionsFromBonds(bonds, 0.0, 0.0);
int first14, second14;
for (int i = 0; i < nonbonded->getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
nonbonded->getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
if ((particle1 == 0 && particle2 == 3) || (particle1 == 3 && particle2 == 0))
first14 = i;
if ((particle1 == 1 && particle2 == 4) || (particle1 == 4 && particle2 == 1))
second14 = i;
}
system.addForce(nonbonded); system.addForce(nonbonded);
//(void) fprintf( log, "%s: Calling context\n", methodName.c_str() ); //(void) fprintf( log, "%s: Calling context\n", methodName.c_str() );
...@@ -240,8 +252,8 @@ void testBrookExclusionsAnd14( FILE* log ){ ...@@ -240,8 +252,8 @@ void testBrookExclusionsAnd14( FILE* log ){
} }
nonbonded->setParticleParameters(0, 0, 1.5, 1); nonbonded->setParticleParameters(0, 0, 1.5, 1);
nonbonded->setParticleParameters(ii, 0, 1.5, 1); nonbonded->setParticleParameters(ii, 0, 1.5, 1);
nonbonded->setNonbonded14Parameters(0, 0, 3, 0, 1.5, ii == 3 ? 0.5 : 0.0); nonbonded->setExceptionParameters(first14, 0, 3, 0, 1.5, ii == 3 ? 0.5 : 0.0);
nonbonded->setNonbonded14Parameters(1, 1, 4, 0, 1.5, 0.0); nonbonded->setExceptionParameters(second14, 1, 4, 0, 1.5, 0.0);
positions[ii] = Vec3(r, 0, 0); positions[ii] = Vec3(r, 0, 0);
context.reinitialize(); context.reinitialize();
...@@ -284,8 +296,8 @@ void testBrookExclusionsAnd14( FILE* log ){ ...@@ -284,8 +296,8 @@ void testBrookExclusionsAnd14( FILE* log ){
nonbonded->setParticleParameters( 0, 2, 1.5, 0 ); nonbonded->setParticleParameters( 0, 2, 1.5, 0 );
nonbonded->setParticleParameters( ii, 2, 1.5, 0 ); nonbonded->setParticleParameters( ii, 2, 1.5, 0 );
nonbonded->setNonbonded14Parameters( 0, 0, 3, ii == 3 ? 4/1.2 : 0, 1.5, 0 ); nonbonded->setExceptionParameters( first14, 0, 3, ii == 3 ? 4/1.2 : 0, 1.5, 0 );
nonbonded->setNonbonded14Parameters( 1, 1, 4, 0, 1.5, 0 ); nonbonded->setExceptionParameters( second14, 1, 4, 0, 1.5, 0 );
context.reinitialize(); context.reinitialize();
......
platforms/brook/tests/TestBrookVerletIntegrator.cpp
View file @ 29e3fa57
...@@ -183,10 +183,10 @@ void testVerletConstraints( FILE* log ){ ...@@ -183,10 +183,10 @@ void testVerletConstraints( FILE* log ){
System system(numParticles, numConstraints); System system(numParticles, numConstraints);
VerletIntegrator integrator(0.001); VerletIntegrator integrator(0.001);
integrator.setConstraintTolerance(1e-5); integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, mass ); system.setParticleMass(i, mass );
forceField->setParticleParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0); forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
} }
for (int i = 0; i < numConstraints; ++i){ for (int i = 0; i < numConstraints; ++i){
system.setConstraintParameters(i, 2*i, 2*i+1, 1.0); system.setConstraintParameters(i, 2*i, 2*i+1, 1.0);
......
platforms/cuda/src/CudaKernels.cpp
View file @ 29e3fa57
...@@ -238,14 +238,27 @@ double CudaCalcRBTorsionForceKernel::executeEnergy(OpenMMContextImpl& context) { ...@@ -238,14 +238,27 @@ double CudaCalcRBTorsionForceKernel::executeEnergy(OpenMMContextImpl& context) {
CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() { CudaCalcNonbondedForceKernel::~CudaCalcNonbondedForceKernel() {
} }
void CudaCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force, const std::vector<std::set<int> >& exclusions) { void CudaCalcNonbondedForceKernel::initialize(const System& system, const NonbondedForce& force) {
if (data.primaryKernel == NULL) if (data.primaryKernel == NULL)
data.primaryKernel = this; data.primaryKernel = this;
data.hasNonbonded = true; data.hasNonbonded = true;
numParticles = force.getNumParticles(); numParticles = force.getNumParticles();
num14 = force.getNumNonbonded14();
_gpuContext* gpu = data.gpu; _gpuContext* gpu = data.gpu;
// Identify which exceptions are 1-4 interactions.
vector<pair<int, int> > exclusions;
vector<int> nb14s;
for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2;
double chargeProd, sigma, epsilon;
force.getExceptionParameters(i, particle1, particle2, chargeProd, sigma, epsilon);
exclusions.push_back(pair<int, int>(particle1, particle2));
if (chargeProd != 0.0 || epsilon != 0.0)
nb14s.push_back(i);
}
num14 = nb14s.size();
// Initialize nonbonded interactions. // Initialize nonbonded interactions.
{ {
...@@ -262,9 +275,12 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -262,9 +275,12 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
q[i] = (float) charge; q[i] = (float) charge;
c6[i] = (float) (4*depth*pow(radius, 6.0)); c6[i] = (float) (4*depth*pow(radius, 6.0));
c12[i] = (float) (4*depth*pow(radius, 12.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); exclusionList[i].push_back(i);
} }
for (int i = 0; i < exclusions.size(); i++) {
exclusionList[exclusions[i].first].push_back(exclusions[i].second);
exclusionList[exclusions[i].second].push_back(exclusions[i].first);
}
CudaNonbondedMethod method = NO_CUTOFF; CudaNonbondedMethod method = NO_CUTOFF;
if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) { if (force.getNonbondedMethod() != NonbondedForce::NoCutoff) {
gpuSetNonbondedCutoff(gpu, force.getCutoffDistance(), 78.3f); gpuSetNonbondedCutoff(gpu, force.getCutoffDistance(), 78.3f);
...@@ -299,7 +315,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon ...@@ -299,7 +315,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
vector<float> q2(num14); vector<float> q2(num14);
for (int i = 0; i < num14; i++) { for (int i = 0; i < num14; i++) {
double charge, sig, eps; double charge, sig, eps;
force.getNonbonded14Parameters(i, particle1[i], particle2[i], charge, sig, eps); force.getExceptionParameters(nb14s[i], particle1[i], particle2[i], charge, sig, eps);
c6[i] = (float) (4*eps*pow(sig, 6.0)); c6[i] = (float) (4*eps*pow(sig, 6.0));
c12[i] = (float) (4*eps*pow(sig, 12.0)); c12[i] = (float) (4*eps*pow(sig, 12.0));
q1[i] = (float) charge; q1[i] = (float) charge;
......
platforms/cuda/src/CudaKernels.h
View file @ 29e3fa57
...@@ -215,11 +215,8 @@ public: ...@@ -215,11 +215,8 @@ public:
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the NonbondedForce this kernel will be used for * @param force the NonbondedForce this kernel will be used for
* @param exclusions the i'th element lists the indices of all particles with which the i'th particle 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.
*/ */
void initialize(const System& system, const NonbondedForce& force, const std::vector<std::set<int> >& exclusions); void initialize(const System& system, const NonbondedForce& force);
/** /**
* Execute the kernel to calculate the forces. * Execute the kernel to calculate the forces.
* *
......
platforms/cuda/tests/TestCudaAndersenThermostat.cpp
View file @ 29e3fa57
...@@ -55,10 +55,10 @@ void testTemperature() { ...@@ -55,10 +55,10 @@ void testTemperature() {
CudaPlatform platform; CudaPlatform platform;
System system(numParticles, 0); System system(numParticles, 0);
VerletIntegrator integrator(0.01); VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 2.0); system.setParticleMass(i, 2.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0); forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
} }
system.addForce(forceField); system.addForce(forceField);
AndersenThermostat* thermstat = new AndersenThermostat(temp, collisionFreq); AndersenThermostat* thermstat = new AndersenThermostat(temp, collisionFreq);
...@@ -93,10 +93,10 @@ void testRandomSeed() { ...@@ -93,10 +93,10 @@ void testRandomSeed() {
CudaPlatform platform; CudaPlatform platform;
System system(numParticles, 0); System system(numParticles, 0);
VerletIntegrator integrator(0.01); VerletIntegrator integrator(0.01);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 2.0); system.setParticleMass(i, 2.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0); forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
} }
system.addForce(forceField); system.addForce(forceField);
AndersenThermostat* thermostat = new AndersenThermostat(temp, collisionFreq); AndersenThermostat* thermostat = new AndersenThermostat(temp, collisionFreq);
......
platforms/cuda/tests/TestCudaBrownianIntegrator.cpp
View file @ 29e3fa57
...@@ -132,10 +132,10 @@ void testConstraints() { ...@@ -132,10 +132,10 @@ void testConstraints() {
System system(numParticles, numConstraints); System system(numParticles, numConstraints);
BrownianIntegrator integrator(temp, 2.0, 0.001); BrownianIntegrator integrator(temp, 2.0, 0.001);
integrator.setConstraintTolerance(1e-5); integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 10.0); system.setParticleMass(i, 10.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0); forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
} }
system.setConstraintParameters(0, 0, 1, 1.0); system.setConstraintParameters(0, 0, 1, 1.0);
system.setConstraintParameters(1, 1, 2, 1.0); system.setConstraintParameters(1, 1, 2, 1.0);
...@@ -178,10 +178,10 @@ void testRandomSeed() { ...@@ -178,10 +178,10 @@ void testRandomSeed() {
CudaPlatform platform; CudaPlatform platform;
System system(numParticles, 0); System system(numParticles, 0);
BrownianIntegrator integrator(temp, 2.0, 0.001); BrownianIntegrator integrator(temp, 2.0, 0.001);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 2.0); system.setParticleMass(i, 2.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0); forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
} }
system.addForce(forceField); system.addForce(forceField);
vector<Vec3> positions(numParticles); vector<Vec3> positions(numParticles);
......
platforms/cuda/tests/TestCudaCMMotionRemover.cpp
View file @ 29e3fa57
...@@ -67,10 +67,10 @@ void testMotionRemoval(Integrator& integrator) { ...@@ -67,10 +67,10 @@ void testMotionRemoval(Integrator& integrator) {
HarmonicBondForce* bonds = new HarmonicBondForce(1); HarmonicBondForce* bonds = new HarmonicBondForce(1);
bonds->setBondParameters(0, 2, 3, 2.0, 0.5); bonds->setBondParameters(0, 2, 3, 2.0, 0.5);
system.addForce(bonds); system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(numParticles, 0); NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, i+1); system.setParticleMass(i, i+1);
nonbonded->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0); nonbonded->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
} }
system.addForce(nonbonded); system.addForce(nonbonded);
CMMotionRemover* remover = new CMMotionRemover(); CMMotionRemover* remover = new CMMotionRemover();
......
platforms/cuda/tests/TestCudaEwald.cpp
View file @ 29e3fa57
...@@ -58,9 +58,9 @@ void testEwald() { ...@@ -58,9 +58,9 @@ void testEwald() {
CudaPlatform platform; CudaPlatform platform;
System system(2, 0); System system(2, 0);
VerletIntegrator integrator(0.01); VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce(2, 0); NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setParticleParameters(0, 1.0, 1, 0); nonbonded->addParticle(1.0, 1, 0);
nonbonded->setParticleParameters(1, -1.0, 1, 0); nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald); nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0; const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff); nonbonded->setCutoffDistance(cutoff);
...@@ -90,10 +90,10 @@ void testPeriodic() { ...@@ -90,10 +90,10 @@ void testPeriodic() {
HarmonicBondForce* bonds = new HarmonicBondForce(1); HarmonicBondForce* bonds = new HarmonicBondForce(1);
bonds->setBondParameters(0, 0, 1, 1, 0); bonds->setBondParameters(0, 0, 1, 1, 0);
system.addForce(bonds); system.addForce(bonds);
NonbondedForce* nonbonded = new NonbondedForce(3, 0); NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->setParticleParameters(0, 1.0, 1, 0); nonbonded->addParticle(1.0, 1, 0);
nonbonded->setParticleParameters(1, 1.0, 1, 0); nonbonded->addParticle(1.0, 1, 0);
nonbonded->setParticleParameters(2, 1.0, 1, 0); nonbonded->addParticle(1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic); nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 2.0; const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff); nonbonded->setCutoffDistance(cutoff);
......
platforms/cuda/tests/TestCudaGBSAOBCForce.cpp
View file @ 29e3fa57
...@@ -57,9 +57,9 @@ void testSingleParticle() { ...@@ -57,9 +57,9 @@ void testSingleParticle() {
system.setParticleMass(0, 2.0); system.setParticleMass(0, 2.0);
LangevinIntegrator integrator(0, 0.1, 0.01); LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* gbsa = new GBSAOBCForce(1); GBSAOBCForce* gbsa = new GBSAOBCForce(1);
NonbondedForce* nonbonded = new NonbondedForce(1, 0); NonbondedForce* nonbonded = new NonbondedForce();
gbsa->setParticleParameters(0, 0.5, 0.15, 1); gbsa->setParticleParameters(0, 0.5, 0.15, 1);
nonbonded->setParticleParameters(0, 0.5, 1, 0); nonbonded->addParticle(0.5, 1, 0);
system.addForce(gbsa); system.addForce(gbsa);
system.addForce(nonbonded); system.addForce(nonbonded);
OpenMMContext context(system, integrator, platform); OpenMMContext context(system, integrator, platform);
...@@ -80,11 +80,11 @@ void testCutoffAndPeriodic() { ...@@ -80,11 +80,11 @@ void testCutoffAndPeriodic() {
System system(2, 0); System system(2, 0);
LangevinIntegrator integrator(0, 0.1, 0.01); LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* gbsa = new GBSAOBCForce(2); GBSAOBCForce* gbsa = new GBSAOBCForce(2);
NonbondedForce* nonbonded = new NonbondedForce(2, 0); NonbondedForce* nonbonded = new NonbondedForce();
gbsa->setParticleParameters(0, -1, 0.15, 1); gbsa->setParticleParameters(0, -1, 0.15, 1);
nonbonded->setParticleParameters(0, -1, 1, 0); nonbonded->addParticle(-1, 1, 0);
gbsa->setParticleParameters(1, 1, 0.15, 1); gbsa->setParticleParameters(1, 1, 0.15, 1);
nonbonded->setParticleParameters(1, 1, 1, 0); nonbonded->addParticle(1, 1, 0);
const double cutoffDistance = 3.0; const double cutoffDistance = 3.0;
const double boxSize = 10.0; const double boxSize = 10.0;
nonbonded->setCutoffDistance(cutoffDistance); nonbonded->setCutoffDistance(cutoffDistance);
...@@ -131,11 +131,11 @@ void testForce(int numParticles, NonbondedForce::NonbondedMethod method) { ...@@ -131,11 +131,11 @@ void testForce(int numParticles, NonbondedForce::NonbondedMethod method) {
System system(numParticles, 0); System system(numParticles, 0);
LangevinIntegrator integrator(0, 0.1, 0.01); LangevinIntegrator integrator(0, 0.1, 0.01);
GBSAOBCForce* gbsa = new GBSAOBCForce(numParticles); GBSAOBCForce* gbsa = new GBSAOBCForce(numParticles);
NonbondedForce* nonbonded = new NonbondedForce(numParticles, 0); NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
double charge = i%2 == 0 ? -1 : 1; double charge = i%2 == 0 ? -1 : 1;
gbsa->setParticleParameters(i, charge, 0.15, 1); gbsa->setParticleParameters(i, charge, 0.15, 1);
nonbonded->setParticleParameters(i, charge, 1, 0); nonbonded->addParticle(charge, 1, 0);
} }
nonbonded->setNonbondedMethod(method); nonbonded->setNonbondedMethod(method);
nonbonded->setCutoffDistance(3.0); nonbonded->setCutoffDistance(3.0);
......
platforms/cuda/tests/TestCudaLangevinIntegrator.cpp
View file @ 29e3fa57
...@@ -100,10 +100,10 @@ void testTemperature() { ...@@ -100,10 +100,10 @@ void testTemperature() {
CudaPlatform platform; CudaPlatform platform;
System system(numParticles, 0); System system(numParticles, 0);
LangevinIntegrator integrator(temp, 2.0, 0.01); LangevinIntegrator integrator(temp, 2.0, 0.01);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 2.0); system.setParticleMass(i, 2.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0); forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
} }
system.addForce(forceField); system.addForce(forceField);
OpenMMContext context(system, integrator, platform); OpenMMContext context(system, integrator, platform);
...@@ -137,10 +137,10 @@ void testConstraints() { ...@@ -137,10 +137,10 @@ void testConstraints() {
System system(numParticles, numConstraints); System system(numParticles, numConstraints);
LangevinIntegrator integrator(temp, 2.0, 0.01); LangevinIntegrator integrator(temp, 2.0, 0.01);
integrator.setConstraintTolerance(1e-5); integrator.setConstraintTolerance(1e-5);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 10.0); system.setParticleMass(i, 10.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0); forceField->addParticle((i%2 == 0 ? 0.2 : -0.2), 0.5, 5.0);
} }
system.setConstraintParameters(0, 0, 1, 1.0); system.setConstraintParameters(0, 0, 1, 1.0);
system.setConstraintParameters(1, 1, 2, 1.0); system.setConstraintParameters(1, 1, 2, 1.0);
...@@ -183,10 +183,10 @@ void testRandomSeed() { ...@@ -183,10 +183,10 @@ void testRandomSeed() {
CudaPlatform platform; CudaPlatform platform;
System system(numParticles, 0); System system(numParticles, 0);
LangevinIntegrator integrator(temp, 2.0, 0.01); LangevinIntegrator integrator(temp, 2.0, 0.01);
NonbondedForce* forceField = new NonbondedForce(numParticles, 0); NonbondedForce* forceField = new NonbondedForce();
for (int i = 0; i < numParticles; ++i) { for (int i = 0; i < numParticles; ++i) {
system.setParticleMass(i, 2.0); system.setParticleMass(i, 2.0);
forceField->setParticleParameters(i, (i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0); forceField->addParticle((i%2 == 0 ? 1.0 : -1.0), 1.0, 5.0);
} }
system.addForce(forceField); system.addForce(forceField);
vector<Vec3> positions(numParticles); vector<Vec3> positions(numParticles);
......
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