Commit 7db6258a authored by peastman's avatar peastman
Browse files

Created API for Gay-Berne force

parent 90ddfc31
...@@ -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-2015 Stanford University and the Authors. * * Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -47,6 +47,7 @@ ...@@ -47,6 +47,7 @@
#include "openmm/CustomNonbondedForce.h" #include "openmm/CustomNonbondedForce.h"
#include "openmm/CustomManyParticleForce.h" #include "openmm/CustomManyParticleForce.h"
#include "openmm/CustomTorsionForce.h" #include "openmm/CustomTorsionForce.h"
#include "openmm/GayBerneForce.h"
#include "openmm/GBSAOBCForce.h" #include "openmm/GBSAOBCForce.h"
#include "openmm/HarmonicAngleForce.h" #include "openmm/HarmonicAngleForce.h"
#include "openmm/HarmonicBondForce.h" #include "openmm/HarmonicBondForce.h"
...@@ -891,6 +892,41 @@ public: ...@@ -891,6 +892,41 @@ public:
virtual void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) = 0; virtual void copyParametersToContext(ContextImpl& context, const CustomManyParticleForce& force) = 0;
}; };
/**
* This kernel is invoked by GayBerneForce to calculate the forces acting on the system and the energy of the system.
*/
class CalcGayBerneForceKernel : public KernelImpl {
public:
static std::string Name() {
return "CalcGayBerneForce";
}
CalcGayBerneForceKernel(std::string name, const Platform& platform) : KernelImpl(name, platform) {
}
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the GayBerneForce this kernel will be used for
*/
virtual void initialize(const System& system, const GayBerneForce& force) = 0;
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
/**
* Copy changed parameters over to a context.
*
* @param context the context to copy parameters to
* @param force the GayBerneForce to copy the parameters from
*/
virtual void copyParametersToContext(ContextImpl& context, const GayBerneForce& force) = 0;
};
/** /**
* This kernel is invoked by VerletIntegrator to take one time step. * This kernel is invoked by VerletIntegrator to take one time step.
*/ */
......
...@@ -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) 2009-2015 Stanford University and the Authors. * * Portions copyright (c) 2009-2016 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -49,6 +49,7 @@ ...@@ -49,6 +49,7 @@
#include "openmm/CustomManyParticleForce.h" #include "openmm/CustomManyParticleForce.h"
#include "openmm/CustomNonbondedForce.h" #include "openmm/CustomNonbondedForce.h"
#include "openmm/Force.h" #include "openmm/Force.h"
#include "openmm/GayBerneForce.h"
#include "openmm/GBSAOBCForce.h" #include "openmm/GBSAOBCForce.h"
#include "openmm/HarmonicAngleForce.h" #include "openmm/HarmonicAngleForce.h"
#include "openmm/HarmonicBondForce.h" #include "openmm/HarmonicBondForce.h"
......
#ifndef OPENMM_GAYBERNEFORCE_H_
#define OPENMM_GAYBERNEFORCE_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "Context.h"
#include "Force.h"
#include <map>
#include <utility>
#include <vector>
#include "internal/windowsExport.h"
namespace OpenMM {
/**
* This class implements the Gay-Berne anisotropic potential. This is similar to a Lennard-Jones potential,
* but it represents the particles as ellipsoids rather than point particles. In addition to the standard
* sigma and epsilon parameters, each particle has three radii r_x, r_y, and r_z that give the size of the
* ellipsoid along each axis. It also has three scale factors e_x, e_y, and e_z that scale the strength
* of the interaction along each axis. You can think of this force as a Lennard-Jones interaction computed
* based on the distance between the nearest points on two ellipsoids. If two particles each have all their
* radii set to 0 and all their scale factors set 1, the interaction simplifies to a standard Lennard-Jones
* force between point particles.
*
* The orientation of a particle's ellipsoid is determined based on the positions of two other particles.
* The vector to the first particle sets the direction of the x axis. The vector to the second particle
* (after subtracting out any x component) sets the direction of the y axis. If the ellipsoid is axially
* symmetric (r_y=r_z and e_y=e_z), you can omit the second particle and define only an x axis direction.
* If the ellipsoid is a sphere (all three radii and all three scale factors are equal), both particles
* can be omitted.
*
* To determine the values of sigma and epsilon for an interaction, this class uses Lorentz-Berthelot
* combining rules: it takes the arithmetic mean of the sigmas and the geometric mean of the epsilons for
* the two interacting particles. You also can specify "exceptions", particular pairs of particles for
* which different values should be used.
*
* To use this class, create a GayBerneForce object, then call addParticle() once for each particle in the
* System to define its parameters. The number of particles for which you define parameters must be exactly
* equal to the number of particles in the System, or else an exception will be thrown when you try to
* create a Context. After a particle has been added, you can modify its force field parameters by calling
* setParticleParameters(). This will have no effect on Contexts that already exist unless you call
* updateParametersInContext().
*
* When using a cutoff, by default interactions are sharply truncated at the cutoff distance. Optionally
* you can instead use a switching function to make the interaction smoothly go to zero over a finite
* distance range. To enable this, call setUseSwitchingFunction(). You must also call setSwitchingDistance()
* to specify the distance at which the interaction should begin to decrease. The switching distance must be
* less than the cutoff distance.
*/
class OPENMM_EXPORT GayBerneForce : public Force {
public:
/**
* This is an enumeration of the different methods that may be used for handling long range nonbonded forces.
*/
enum NonbondedMethod {
/**
* No cutoff is applied to nonbonded interactions. The full set of N^2 interactions is computed exactly.
* This necessarily means that periodic boundary conditions cannot be used. This is the default.
*/
NoCutoff = 0,
/**
* Interactions beyond the cutoff distance are ignored.
*/
CutoffNonPeriodic = 1,
/**
* Periodic boundary conditions are used, so that each particle interacts only with the nearest periodic copy of
* each other particle. Interactions beyond the cutoff distance are ignored.
*/
CutoffPeriodic = 2
};
/**
* Create a GayBerneForce.
*/
GayBerneForce();
/**
* Get the number of particles for which force field parameters have been defined.
*/
int getNumParticles() const {
return particles.size();
}
/**
* Get the number of special interactions that should be calculated differently from other interactions.
*/
int getNumExceptions() const {
return exceptions.size();
}
/**
* Get the method used for handling long range interactions.
*/
NonbondedMethod getNonbondedMethod() const;
/**
* Set the method used for handling long range interactions.
*/
void setNonbondedMethod(NonbondedMethod method);
/**
* Get the cutoff distance (in nm) being used for interactions. If the NonbondedMethod in use
* is NoCutoff, this value will have no effect.
*
* @return the cutoff distance, measured in nm
*/
double getCutoffDistance() const;
/**
* Set the cutoff distance (in nm) being used for interactions. If the NonbondedMethod in use
* is NoCutoff, this value will have no effect.
*
* @param distance the cutoff distance, measured in nm
*/
void setCutoffDistance(double distance);
/**
* Get whether a switching function is applied to the interaction. If the nonbonded method is set
* to NoCutoff, this option is ignored.
*/
bool getUseSwitchingFunction() const;
/**
* Set whether a switching function is applied to the interaction. If the nonbonded method is set
* to NoCutoff, this option is ignored.
*/
void setUseSwitchingFunction(bool use);
/**
* Get the distance at which the switching function begins to reduce the interaction. This must be
* less than the cutoff distance.
*/
double getSwitchingDistance() const;
/**
* Set the distance at which the switching function begins to reduce the interaction. This must be
* less than the cutoff distance.
*/
void setSwitchingDistance(double distance);
/**
* Add the 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.
*
* @param sigma the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
* @param xparticle the index of the particle whose position defines the ellipsoid's x axis, or -1 if the ellipsoid is a sphere
* @param yparticle the index of the particle whose position defines the ellipsoid's y axis, or -1 if the ellipsoid is axially symmetric
* @param rx the radius of the ellipsoid along its x axis
* @param ry the radius of the ellipsoid along its y axis
* @param rz the radius of the ellipsoid along its z axis
* @param ex the factor by which epsilon is scaled along the ellipsoid's x axis
* @param ey the factor by which epsilon is scaled along the ellipsoid's y axis
* @param ez the factor by which epsilon is scaled along the ellipsoid's z axis
* @return the index of the particle that was added
*/
int addParticle(double sigma, double epsilon, int xparticle, int yparticle, double rx, double ry, double rz, double ex, double ey, double ez);
/**
* Get the parameters for a particle.
*
* @param index the index of the particle for which to get parameters
* @param[out] sigma the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
* @param[out] epsilon the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
* @param[out] xparticle the index of the particle whose position defines the ellipsoid's x axis, or -1 if the ellipsoid is a sphere
* @param[out] yparticle the index of the particle whose position defines the ellipsoid's y axis, or -1 if the ellipsoid is axially symmetric
* @param[out] rx the radius of the ellipsoid along its x axis
* @param[out] ry the radius of the ellipsoid along its y axis
* @param[out] rz the radius of the ellipsoid along its z axis
* @param[out] ex the factor by which epsilon is scaled along the ellipsoid's x axis
* @param[out] ey the factor by which epsilon is scaled along the ellipsoid's y axis
* @param[out] ez the factor by which epsilon is scaled along the ellipsoid's z axis
*/
void getParticleParameters(int index, double& sigma, double& epsilon, int& xparticle, int& yparticle, double& rx, double& ry, double& rz, double& ex, double& ey, double& ez) const;
/**
* Set the parameters for a particle.
*
* @param index the index of the particle for which to set parameters
* @param sigma the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
* @param xparticle the index of the particle whose position defines the ellipsoid's x axis, or -1 if the ellipsoid is a sphere
* @param yparticle the index of the particle whose position defines the ellipsoid's y axis, or -1 if the ellipsoid is axially symmetric
* @param rx the radius of the ellipsoid along its x axis
* @param ry the radius of the ellipsoid along its y axis
* @param rz the radius of the ellipsoid along its z axis
* @param ex the factor by which epsilon is scaled along the ellipsoid's x axis
* @param ey the factor by which epsilon is scaled along the ellipsoid's y axis
* @param ez the factor by which epsilon is scaled along the ellipsoid's z axis
*/
void setParticleParameters(int index, double sigma, double epsilon, int xparticle, int yparticle, double rx, double ry, double rz, double ex, double ey, double ez);
/**
* Add an interaction to the list of exceptions that should be calculated differently from other interactions. If
* epsilon is equal to 0, this will cause the interaction to be completely omitted from force and energy calculations.
*
* @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 sigma the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
* @param replace determines the behavior if there is already an exception for the same two particles. If true, the existing one is replaced. If false,
* an exception is thrown.
* @return the index of the exception that was added
*/
int addException(int particle1, int particle2, double sigma, double epsilon, bool replace = false);
/**
* 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[out] particle1 the index of the first particle involved in the interaction
* @param[out] particle2 the index of the second particle involved in the interaction
* @param[out] sigma the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
* @param[out] epsilon the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
*/
void getExceptionParameters(int index, int& particle1, int& particle2, double& sigma, double& epsilon) const;
/**
* Set the force field parameters for an interaction that should be calculated differently from others. If
* epsilon is 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 particle1 the index of the first particle involved in the interaction
* @param particle2 the index of the second particle involved in the interaction
* @param sigma the sigma parameter (corresponding to the van der Waals radius of the particle), measured in nm
* @param epsilon the epsilon parameter (corresponding to the well depth of the van der Waals interaction), measured in kJ/mol
*/
void setExceptionParameters(int index, int particle1, int particle2, double sigma, double epsilon);
/**
* Update the particle and exception parameters in a Context to match those stored in this Force object. This method
* provides an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
* Simply call setParticleParameters() and setExceptionParameters() to modify this object's parameters, then call
* updateParametersInContext() to copy them over to the Context.
*
* This method has several limitations. The only information it updates is the parameters of particles and exceptions.
* All other aspects of the Force (the nonbonded method, the cutoff distance, etc.) are unaffected and can only be
* changed by reinitializing the Context. Furthermore, only the sigma and epsilon values of an exception can be
* changed; the pair of particles involved in the exception cannot change. Finally, this method cannot be used to
* add new particles or exceptions, only to change the parameters of existing ones.
*/
void updateParametersInContext(Context& context);
/**
* Returns whether or not this force makes use of periodic boundary
* conditions.
*
* @returns true if force uses PBC and false otherwise
*/
bool usesPeriodicBoundaryConditions() const {
return nonbondedMethod == GayBerneForce::CutoffPeriodic;
}
protected:
ForceImpl* createImpl() const;
private:
class ParticleInfo;
class ExceptionInfo;
NonbondedMethod nonbondedMethod;
double cutoffDistance, switchingDistance;
bool useSwitchingFunction;
std::vector<ParticleInfo> particles;
std::vector<ExceptionInfo> exceptions;
std::map<std::pair<int, int>, int> exceptionMap;
};
/**
* This is an internal class used to record information about a particle.
* @private
*/
class GayBerneForce::ParticleInfo {
public:
int xparticle, yparticle;
double sigma, epsilon, rx, ry, rz, ex, ey, ez;
ParticleInfo() {
xparticle = yparticle = -1;
sigma = epsilon = rx = ry = rz = ex = ey = ez = 0.0;
}
ParticleInfo(double sigma, double epsilon, int xparticle, int yparticle, double rx, double ry, double rz, double ex, double ey, double ez) :
sigma(sigma), epsilon(epsilon), xparticle(xparticle), yparticle(yparticle), rx(rx), ry(ry), rz(rz), ex(ex), ey(ey), ez(ez) {
}
};
/**
* This is an internal class used to record information about an exception.
* @private
*/
class GayBerneForce::ExceptionInfo {
public:
int particle1, particle2;
double sigma, epsilon;
ExceptionInfo() {
particle1 = particle2 = -1;
sigma = epsilon = 0.0;
}
ExceptionInfo(int particle1, int particle2, double sigma, double epsilon) :
particle1(particle1), particle2(particle2), sigma(sigma), epsilon(epsilon) {
}
};
} // namespace OpenMM
#endif /*OPENMM_GAYBERNEFORCE_H_*/
#ifndef OPENMM_GAYBERNEFORCEIMPL_H_
#define OPENMM_GAYBERNEFORCEIMPL_H_
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "ForceImpl.h"
#include "openmm/GayBerneForce.h"
#include "openmm/Kernel.h"
#include <utility>
#include <set>
#include <string>
namespace OpenMM {
class System;
/**
* This is the internal implementation of GayBerneForce.
*/
class OPENMM_EXPORT GayBerneForceImpl : public ForceImpl {
public:
GayBerneForceImpl(const GayBerneForce& owner);
~GayBerneForceImpl();
void initialize(ContextImpl& context);
const GayBerneForce& getOwner() const {
return owner;
}
void updateContextState(ContextImpl& context) {
// This force field doesn't update the state directly.
}
double calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups);
std::map<std::string, double> getDefaultParameters() {
return std::map<std::string, double>(); // This force field doesn't define any parameters.
}
std::vector<std::string> getKernelNames();
void updateParametersInContext(ContextImpl& context);
private:
const GayBerneForce& owner;
Kernel kernel;
};
} // namespace OpenMM
#endif /*OPENMM_GAYBERNEFORCEIMPL_H_*/
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/Force.h"
#include "openmm/OpenMMException.h"
#include "openmm/GayBerneForce.h"
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/internal/GayBerneForceImpl.h"
#include <cmath>
#include <map>
#include <sstream>
#include <utility>
using namespace OpenMM;
using std::map;
using std::pair;
using std::set;
using std::string;
using std::stringstream;
using std::vector;
GayBerneForce::GayBerneForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), switchingDistance(-1.0), useSwitchingFunction(false) {
}
GayBerneForce::NonbondedMethod GayBerneForce::getNonbondedMethod() const {
return nonbondedMethod;
}
void GayBerneForce::setNonbondedMethod(NonbondedMethod method) {
nonbondedMethod = method;
}
double GayBerneForce::getCutoffDistance() const {
return cutoffDistance;
}
void GayBerneForce::setCutoffDistance(double distance) {
cutoffDistance = distance;
}
bool GayBerneForce::getUseSwitchingFunction() const {
return useSwitchingFunction;
}
void GayBerneForce::setUseSwitchingFunction(bool use) {
useSwitchingFunction = use;
}
double GayBerneForce::getSwitchingDistance() const {
return switchingDistance;
}
void GayBerneForce::setSwitchingDistance(double distance) {
switchingDistance = distance;
}
int GayBerneForce::addParticle(double sigma, double epsilon, int xparticle, int yparticle, double rx, double ry, double rz, double ex, double ey, double ez) {
particles.push_back(ParticleInfo(sigma, epsilon, xparticle, yparticle, rx, ry, rz, ex, ey, ez));
return particles.size()-1;
}
void GayBerneForce::getParticleParameters(int index, double& sigma, double& epsilon, int& xparticle, int& yparticle, double& rx, double& ry, double& rz, double& ex, double& ey, double& ez) const {
ASSERT_VALID_INDEX(index, particles);
sigma = particles[index].sigma;
epsilon = particles[index].epsilon;
xparticle = particles[index].xparticle;
yparticle = particles[index].yparticle;
rx = particles[index].rx;
ry = particles[index].ry;
rz = particles[index].rz;
ex = particles[index].ex;
ey = particles[index].ey;
ez = particles[index].ez;
}
void GayBerneForce::setParticleParameters(int index, double sigma, double epsilon, int xparticle, int yparticle, double rx, double ry, double rz, double ex, double ey, double ez) {
ASSERT_VALID_INDEX(index, particles);
particles[index].sigma = sigma;
particles[index].epsilon = epsilon;
particles[index].xparticle = xparticle;
particles[index].yparticle = yparticle;
particles[index].rx = rx;
particles[index].ry = ry;
particles[index].rz = rx;
particles[index].ex = ex;
particles[index].ey = ey;
particles[index].ez = ez;
}
int GayBerneForce::addException(int particle1, int particle2, double sigma, double epsilon, bool replace) {
map<pair<int, int>, int>::iterator iter = exceptionMap.find(pair<int, int>(particle1, particle2));
int newIndex;
if (iter == exceptionMap.end())
iter = exceptionMap.find(pair<int, int>(particle2, particle1));
if (iter != exceptionMap.end()) {
if (!replace) {
stringstream msg;
msg << "GayBerneForce: There is already an exception for particles ";
msg << particle1;
msg << " and ";
msg << particle2;
throw OpenMMException(msg.str());
}
exceptions[iter->second] = ExceptionInfo(particle1, particle2, sigma, epsilon);
newIndex = iter->second;
exceptionMap.erase(iter->first);
}
else {
exceptions.push_back(ExceptionInfo(particle1, particle2, sigma, epsilon));
newIndex = exceptions.size()-1;
}
exceptionMap[pair<int, int>(particle1, particle2)] = newIndex;
return newIndex;
}
void GayBerneForce::getExceptionParameters(int index, int& particle1, int& particle2, double& sigma, double& epsilon) const {
ASSERT_VALID_INDEX(index, exceptions);
particle1 = exceptions[index].particle1;
particle2 = exceptions[index].particle2;
sigma = exceptions[index].sigma;
epsilon = exceptions[index].epsilon;
}
void GayBerneForce::setExceptionParameters(int index, int particle1, int particle2, double sigma, double epsilon) {
ASSERT_VALID_INDEX(index, exceptions);
exceptions[index].particle1 = particle1;
exceptions[index].particle2 = particle2;
exceptions[index].sigma = sigma;
exceptions[index].epsilon = epsilon;
}
ForceImpl* GayBerneForce::createImpl() const {
return new GayBerneForceImpl(*this);
}
void GayBerneForce::updateParametersInContext(Context& context) {
dynamic_cast<GayBerneForceImpl&>(getImplInContext(context)).updateParametersInContext(getContextImpl(context));
}
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* Simbios, the NIH National Center for Physics-Based Simulation of *
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2016 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
#include "openmm/OpenMMException.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/internal/GayBerneForceImpl.h"
#include "openmm/kernels.h"
#include <set>
#include <sstream>
using namespace OpenMM;
using namespace std;
GayBerneForceImpl::GayBerneForceImpl(const GayBerneForce& owner) : owner(owner) {
}
GayBerneForceImpl::~GayBerneForceImpl() {
}
void GayBerneForceImpl::initialize(ContextImpl& context) {
kernel = context.getPlatform().createKernel(CalcGayBerneForceKernel::Name(), context);
// Check for errors in the specification of exceptions.
const System& system = context.getSystem();
if (owner.getNumParticles() != system.getNumParticles())
throw OpenMMException("GayBerneForce must have exactly as many particles as the System it belongs to.");
if (owner.getUseSwitchingFunction()) {
if (owner.getSwitchingDistance() < 0 || owner.getSwitchingDistance() >= owner.getCutoffDistance())
throw OpenMMException("GayBerneForce: Switching distance must satisfy 0 <= r_switch < r_cutoff");
}
for (int i = 0; i < owner.getNumParticles(); i++) {
int xparticle, yparticle;
double sigma, epsilon, rx, ry, rz, ex, ey, ez;
owner.getParticleParameters(i, sigma, epsilon, xparticle, yparticle, rx, ry, rz, ex, ey, ez);
if (xparticle < -1 || xparticle >= owner.getNumParticles()) {
stringstream msg;
msg << "GayBerneForce: Illegal particle index for xparticle: ";
msg << xparticle;
throw OpenMMException(msg.str());
}
if (yparticle < -1 || yparticle >= owner.getNumParticles()) {
stringstream msg;
msg << "GayBerneForce: Illegal particle index for a yparticle: ";
msg << yparticle;
throw OpenMMException(msg.str());
}
if (yparticle == -1 && (ry != rz || ey != ez))
throw OpenMMException("GayBerneForce: yparticle is -1 for a particle that is not axially symmetric");
if (xparticle == -1 && (rx != rz || ex != ez))
throw OpenMMException("GayBerneForce: xparticle is -1 for a particle that is not spherical");
if (xparticle == -1 && yparticle != -1)
throw OpenMMException("GayBerneForce: xparticle cannot be -1 if yparticle is not also -1");
}
vector<set<int> > exceptions(owner.getNumParticles());
for (int i = 0; i < owner.getNumExceptions(); i++) {
int particle1, particle2;
double sigma, epsilon;
owner.getExceptionParameters(i, particle1, particle2, sigma, epsilon);
if (particle1 < 0 || particle1 >= owner.getNumParticles()) {
stringstream msg;
msg << "GayBerneForce: Illegal particle index for an exception: ";
msg << particle1;
throw OpenMMException(msg.str());
}
if (particle2 < 0 || particle2 >= owner.getNumParticles()) {
stringstream msg;
msg << "GayBerneForce: 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 << "GayBerneForce: Multiple exceptions are specified for particles ";
msg << particle1;
msg << " and ";
msg << particle2;
throw OpenMMException(msg.str());
}
exceptions[particle1].insert(particle2);
exceptions[particle2].insert(particle1);
}
if (owner.getNonbondedMethod() == GayBerneForce::CutoffPeriodic) {
Vec3 boxVectors[3];
system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
double cutoff = owner.getCutoffDistance();
if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
throw OpenMMException("GayBerneForce: The cutoff distance cannot be greater than half the periodic box size.");
}
kernel.getAs<CalcGayBerneForceKernel>().initialize(context.getSystem(), owner);
}
double GayBerneForceImpl::calcForcesAndEnergy(ContextImpl& context, bool includeForces, bool includeEnergy, int groups) {
return kernel.getAs<CalcGayBerneForceKernel>().execute(context, includeForces, includeEnergy);
}
std::vector<std::string> GayBerneForceImpl::getKernelNames() {
std::vector<std::string> names;
names.push_back(CalcGayBerneForceKernel::Name());
return names;
}
void GayBerneForceImpl::updateParametersInContext(ContextImpl& context) {
kernel.getAs<CalcGayBerneForceKernel>().copyParametersToContext(context, owner);
}
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