Skip to content

GitLab

Commit 454eda76 authored by Saurabh Belsare's avatar Saurabh Belsare
Browse files

Added required headers. Still need to actually implement the functions

parent ad966d1c
Hide whitespace changes
Inline Side-by-side
Showing with 274 additions and 269 deletions
plugins/amoeba/openmmapi/include/openmm/AmoebaMultipoleForce.h
View file @ 454eda76
...@@ -316,15 +316,6 @@ public: ...@@ -316,15 +316,6 @@ public:
* This can be overridden by explicitly setting an alpha parameter and grid dimensions to use. * This can be overridden by explicitly setting an alpha parameter and grid dimensions to use.
*/ */
void setEwaldErrorTolerance(double tol); void setEwaldErrorTolerance(double tol);
/**
* Get the induced dipole moments of all particles.
*
* @param context the Context for which to get the induced dipoles
* @param[out] dipoles the induced dipole moment of particle i is stored into the i'th element
*/
void getInducedDipoles(Context& context, std::vector<Vec3>& dipoles);
/** /**
* Get the fixed dipole moments of all particles in the global reference frame. * Get the fixed dipole moments of all particles in the global reference frame.
* *
...@@ -332,6 +323,13 @@ public: ...@@ -332,6 +323,13 @@ public:
* @param[out] dipoles the fixed dipole moment of particle i is stored into the i'th element * @param[out] dipoles the fixed dipole moment of particle i is stored into the i'th element
*/ */
void getLabFramePermanentDipoles(Context& context, std::vector<Vec3>& dipoles); void getLabFramePermanentDipoles(Context& context, std::vector<Vec3>& dipoles);
/**
* Get the induced dipole moments of all particles.
*
* @param context the Context for which to get the induced dipoles
* @param[out] dipoles the induced dipole moment of particle i is stored into the i'th element
*/
void getInducedDipoles(Context& context, std::vector<Vec3>& dipoles);
/** /**
* Get the electrostatic potential. * Get the electrostatic potential.
......
plugins/amoeba/openmmapi/include/openmm/amoebaKernels.h
View file @ 454eda76
...@@ -59,7 +59,7 @@ public: ...@@ -59,7 +59,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the AmoebaBondForce this kernel will be used for * @param force the AmoebaBondForce this kernel will be used for
*/ */
...@@ -99,7 +99,7 @@ public: ...@@ -99,7 +99,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the AmoebaAngleForce this kernel will be used for * @param force the AmoebaAngleForce this kernel will be used for
*/ */
...@@ -139,7 +139,7 @@ public: ...@@ -139,7 +139,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the AmoebaInPlaneAngleForce this kernel will be used for * @param force the AmoebaInPlaneAngleForce this kernel will be used for
*/ */
...@@ -179,7 +179,7 @@ public: ...@@ -179,7 +179,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the PiTorsionForce this kernel will be used for * @param force the PiTorsionForce this kernel will be used for
*/ */
...@@ -219,7 +219,7 @@ public: ...@@ -219,7 +219,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the StretchBendForce this kernel will be used for * @param force the StretchBendForce this kernel will be used for
*/ */
...@@ -259,7 +259,7 @@ public: ...@@ -259,7 +259,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the OutOfPlaneBendForce this kernel will be used for * @param force the OutOfPlaneBendForce this kernel will be used for
*/ */
...@@ -299,7 +299,7 @@ public: ...@@ -299,7 +299,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the TorsionTorsionForce this kernel will be used for * @param force the TorsionTorsionForce this kernel will be used for
*/ */
...@@ -332,7 +332,7 @@ public: ...@@ -332,7 +332,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the MultipoleForce this kernel will be used for * @param force the MultipoleForce this kernel will be used for
*/ */
...@@ -348,6 +348,7 @@ public: ...@@ -348,6 +348,7 @@ public:
*/ */
virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0; virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
virtual void getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles) = 0;
virtual void getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles) = 0; virtual void getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles) = 0;
virtual void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid, virtual void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
...@@ -364,7 +365,7 @@ public: ...@@ -364,7 +365,7 @@ public:
/** /**
* Get the parameters being used for PME. * Get the parameters being used for PME.
* *
* @param alpha the separation parameter * @param alpha the separation parameter
* @param nx the number of grid points along the X axis * @param nx the number of grid points along the X axis
* @param ny the number of grid points along the Y axis * @param ny the number of grid points along the Y axis
...@@ -389,7 +390,7 @@ public: ...@@ -389,7 +390,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the GBSAOBCForce this kernel will be used for * @param force the GBSAOBCForce this kernel will be used for
*/ */
...@@ -429,7 +430,7 @@ public: ...@@ -429,7 +430,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the GBSAOBCForce this kernel will be used for * @param force the GBSAOBCForce this kernel will be used for
*/ */
...@@ -469,7 +470,7 @@ public: ...@@ -469,7 +470,7 @@ public:
/** /**
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the GBSAOBCForce this kernel will be used for * @param force the GBSAOBCForce this kernel will be used for
*/ */
......
plugins/amoeba/openmmapi/include/openmm/internal/AmoebaMultipoleForceImpl.h
View file @ 454eda76
...@@ -64,7 +64,7 @@ public: ...@@ -64,7 +64,7 @@ public:
/** /**
* Get the CovalentMap for an atom * Get the CovalentMap for an atom
* *
* @param force AmoebaMultipoleForce force reference * @param force AmoebaMultipoleForce force reference
* @param index the index of the atom for which to set parameters * @param index the index of the atom for which to set parameters
* @param minCovalentIndex minimum covalent index * @param minCovalentIndex minimum covalent index
...@@ -76,15 +76,13 @@ public: ...@@ -76,15 +76,13 @@ public:
/** /**
* Get the covalent degree for the CovalentEnd lists * Get the covalent degree for the CovalentEnd lists
* *
* @param force AmoebaMultipoleForce force reference * @param force AmoebaMultipoleForce force reference
* @param covalentDegree covalent degrees for the CovalentEnd lists * @param covalentDegree covalent degrees for the CovalentEnd lists
*/ */
static void getCovalentDegree(const AmoebaMultipoleForce& force, std::vector<int>& covalentDegree); static void getCovalentDegree(const AmoebaMultipoleForce& force, std::vector<int>& covalentDegree);
void getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
void getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles); void getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
void getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid, void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
std::vector< double >& outputElectrostaticPotential); std::vector< double >& outputElectrostaticPotential);
...@@ -92,7 +90,7 @@ public: ...@@ -92,7 +90,7 @@ public:
void getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments); void getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments);
void updateParametersInContext(ContextImpl& context); void updateParametersInContext(ContextImpl& context);
void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const; void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
private: private:
const AmoebaMultipoleForce& owner; const AmoebaMultipoleForce& owner;
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp
View file @ 454eda76
...@@ -69,19 +69,19 @@ void AmoebaMultipoleForce::setCutoffDistance(double distance) { ...@@ -69,19 +69,19 @@ void AmoebaMultipoleForce::setCutoffDistance(double distance) {
cutoffDistance = distance; cutoffDistance = distance;
} }
double AmoebaMultipoleForce::getAEwald() const { double AmoebaMultipoleForce::getAEwald() const {
return aewald; return aewald;
} }
void AmoebaMultipoleForce::setAEwald(double inputAewald) { void AmoebaMultipoleForce::setAEwald(double inputAewald) {
aewald = inputAewald; aewald = inputAewald;
} }
int AmoebaMultipoleForce::getPmeBSplineOrder() const { int AmoebaMultipoleForce::getPmeBSplineOrder() const {
return pmeBSplineOrder; return pmeBSplineOrder;
} }
void AmoebaMultipoleForce::getPmeGridDimensions(std::vector<int>& gridDimension) const { void AmoebaMultipoleForce::getPmeGridDimensions(std::vector<int>& gridDimension) const {
if (gridDimension.size() < 3) { if (gridDimension.size() < 3) {
gridDimension.resize(3); gridDimension.resize(3);
} }
...@@ -93,15 +93,15 @@ void AmoebaMultipoleForce::getPmeGridDimensions(std::vector<int>& gridDimension) ...@@ -93,15 +93,15 @@ void AmoebaMultipoleForce::getPmeGridDimensions(std::vector<int>& gridDimension)
gridDimension[0] = gridDimension[1] = gridDimension[2] = 0; gridDimension[0] = gridDimension[1] = gridDimension[2] = 0;
} }
return; return;
} }
void AmoebaMultipoleForce::setPmeGridDimensions(const std::vector<int>& gridDimension) { void AmoebaMultipoleForce::setPmeGridDimensions(const std::vector<int>& gridDimension) {
pmeGridDimension.resize(3); pmeGridDimension.resize(3);
pmeGridDimension[0] = gridDimension[0]; pmeGridDimension[0] = gridDimension[0];
pmeGridDimension[1] = gridDimension[1]; pmeGridDimension[1] = gridDimension[1];
pmeGridDimension[2] = gridDimension[2]; pmeGridDimension[2] = gridDimension[2];
return; return;
} }
void AmoebaMultipoleForce::getPMEParametersInContext(const Context& context, double& alpha, int& nx, int& ny, int& nz) const { void AmoebaMultipoleForce::getPMEParametersInContext(const Context& context, double& alpha, int& nx, int& ny, int& nz) const {
dynamic_cast<const AmoebaMultipoleForceImpl&>(getImplInContext(context)).getPMEParameters(alpha, nx, ny, nz); dynamic_cast<const AmoebaMultipoleForceImpl&>(getImplInContext(context)).getPMEParameters(alpha, nx, ny, nz);
...@@ -131,13 +131,13 @@ void AmoebaMultipoleForce::setEwaldErrorTolerance(double tol) { ...@@ -131,13 +131,13 @@ void AmoebaMultipoleForce::setEwaldErrorTolerance(double tol) {
ewaldErrorTol = tol; ewaldErrorTol = tol;
} }
int AmoebaMultipoleForce::addMultipole(double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole, int axisType, int AmoebaMultipoleForce::addMultipole(double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole, int axisType,
int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity) { int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity) {
multipoles.push_back(MultipoleInfo(charge, molecularDipole, molecularQuadrupole, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY, thole, dampingFactor, polarity)); multipoles.push_back(MultipoleInfo(charge, molecularDipole, molecularQuadrupole, axisType, multipoleAtomZ, multipoleAtomX, multipoleAtomY, thole, dampingFactor, polarity));
return multipoles.size()-1; return multipoles.size()-1;
} }
void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole, void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole,
int& axisType, int& multipoleAtomZ, int& multipoleAtomX, int& multipoleAtomY, double& thole, double& dampingFactor, double& polarity) const { int& axisType, int& multipoleAtomZ, int& multipoleAtomX, int& multipoleAtomY, double& thole, double& dampingFactor, double& polarity) const {
charge = multipoles[index].charge; charge = multipoles[index].charge;
...@@ -167,7 +167,7 @@ void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std ...@@ -167,7 +167,7 @@ void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std
polarity = multipoles[index].polarity; polarity = multipoles[index].polarity;
} }
void AmoebaMultipoleForce::setMultipoleParameters(int index, double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole, void AmoebaMultipoleForce::setMultipoleParameters(int index, double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole,
int axisType, int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity) { int axisType, int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity) {
multipoles[index].charge = charge; multipoles[index].charge = charge;
...@@ -230,8 +230,12 @@ void AmoebaMultipoleForce::getCovalentMaps(int index, std::vector< std::vector<i ...@@ -230,8 +230,12 @@ void AmoebaMultipoleForce::getCovalentMaps(int index, std::vector< std::vector<i
} }
} }
void AmoebaMultipoleForce::getInducedDipoles(Context& context, vector<Vec3>& dipoles) { void AmoebaMultipoleForce::getLabFramePermanentDipoles(Context& context, vector<Vec3>& dipoles) {
dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getInducedDipoles(getContextImpl(context), dipoles); dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getLabFramePermanentDipoles(getContextImpl(context), dipoles);
}
void AmoebaMultipoleForce::getLabFramePermanentDipoles(Context& context, vector<Vec3>& dipoles) {
dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getLabFramePermanentDipoles(getContextImpl(context), dipoles);
} }
void AmoebaMultipoleForce::getElectrostaticPotential(const std::vector< Vec3 >& inputGrid, Context& context, std::vector< double >& outputElectrostaticPotential) { void AmoebaMultipoleForce::getElectrostaticPotential(const std::vector< Vec3 >& inputGrid, Context& context, std::vector< double >& outputElectrostaticPotential) {
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForceImpl.cpp
View file @ 454eda76
...@@ -61,7 +61,7 @@ void AmoebaMultipoleForceImpl::initialize(ContextImpl& context) { ...@@ -61,7 +61,7 @@ void AmoebaMultipoleForceImpl::initialize(ContextImpl& context) {
double cutoff = owner.getCutoffDistance(); double cutoff = owner.getCutoffDistance();
if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2]) if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
throw OpenMMException("AmoebaMultipoleForce: The cutoff distance cannot be greater than half the periodic box size."); throw OpenMMException("AmoebaMultipoleForce: The cutoff distance cannot be greater than half the periodic box size.");
} }
double quadrupoleValidationTolerance = 1.0e-05; double quadrupoleValidationTolerance = 1.0e-05;
for (int ii = 0; ii < system.getNumParticles(); ii++) { for (int ii = 0; ii < system.getNumParticles(); ii++) {
...@@ -170,7 +170,7 @@ void AmoebaMultipoleForceImpl::getCovalentRange(const AmoebaMultipoleForce& forc ...@@ -170,7 +170,7 @@ void AmoebaMultipoleForceImpl::getCovalentRange(const AmoebaMultipoleForce& forc
*maxCovalentIndex = covalentList[ii]; *maxCovalentIndex = covalentList[ii];
} }
} }
} }
return; return;
} }
...@@ -179,10 +179,14 @@ void AmoebaMultipoleForceImpl::getCovalentDegree(const AmoebaMultipoleForce& for ...@@ -179,10 +179,14 @@ void AmoebaMultipoleForceImpl::getCovalentDegree(const AmoebaMultipoleForce& for
const int* CovalentDegrees = AmoebaMultipoleForceImpl::getCovalentDegrees(); const int* CovalentDegrees = AmoebaMultipoleForceImpl::getCovalentDegrees();
for (unsigned int kk = 0; kk < AmoebaMultipoleForce::CovalentEnd; kk++) { for (unsigned int kk = 0; kk < AmoebaMultipoleForce::CovalentEnd; kk++) {
covalentDegree[kk] = CovalentDegrees[kk]; covalentDegree[kk] = CovalentDegrees[kk];
} }
return; return;
} }
void AmoebaMultipoleForceImpl::getLabFramePermanentDipoles(ContextImpl& context, vector<Vec3>& dipoles) {
kernel.getAs<CalcAmoebaMultipoleForceKernel>().getLabFramePermanentDipoles(context, dipoles);
}
void AmoebaMultipoleForceImpl::getInducedDipoles(ContextImpl& context, vector<Vec3>& dipoles) { void AmoebaMultipoleForceImpl::getInducedDipoles(ContextImpl& context, vector<Vec3>& dipoles) {
kernel.getAs<CalcAmoebaMultipoleForceKernel>().getInducedDipoles(context, dipoles); kernel.getAs<CalcAmoebaMultipoleForceKernel>().getInducedDipoles(context, dipoles);
} }
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
View file @ 454eda76
...@@ -35,7 +35,7 @@ using namespace OpenMM; ...@@ -35,7 +35,7 @@ using namespace OpenMM;
AmoebaReferenceMultipoleForce::AmoebaReferenceMultipoleForce() : AmoebaReferenceMultipoleForce::AmoebaReferenceMultipoleForce() :
_nonbondedMethod(NoCutoff), _nonbondedMethod(NoCutoff),
_numParticles(0), _numParticles(0),
_electric(138.9354558456), _electric(138.9354558456),
_dielectric(1.0), _dielectric(1.0),
_mutualInducedDipoleConverged(0), _mutualInducedDipoleConverged(0),
...@@ -51,7 +51,7 @@ AmoebaReferenceMultipoleForce::AmoebaReferenceMultipoleForce() : ...@@ -51,7 +51,7 @@ AmoebaReferenceMultipoleForce::AmoebaReferenceMultipoleForce() :
AmoebaReferenceMultipoleForce::AmoebaReferenceMultipoleForce(NonbondedMethod nonbondedMethod) : AmoebaReferenceMultipoleForce::AmoebaReferenceMultipoleForce(NonbondedMethod nonbondedMethod) :
_nonbondedMethod(nonbondedMethod), _nonbondedMethod(nonbondedMethod),
_numParticles(0), _numParticles(0),
_electric(138.9354558456), _electric(138.9354558456),
_dielectric(1.0), _dielectric(1.0),
_mutualInducedDipoleConverged(0), _mutualInducedDipoleConverged(0),
...@@ -97,7 +97,7 @@ void AmoebaReferenceMultipoleForce::initialize() ...@@ -97,7 +97,7 @@ void AmoebaReferenceMultipoleForce::initialize()
_uScale[index++] = 1.0; _uScale[index++] = 1.0;
} }
AmoebaReferenceMultipoleForce::NonbondedMethod AmoebaReferenceMultipoleForce::getNonbondedMethod() const AmoebaReferenceMultipoleForce::NonbondedMethod AmoebaReferenceMultipoleForce::getNonbondedMethod() const
{ {
return _nonbondedMethod; return _nonbondedMethod;
} }
...@@ -107,7 +107,7 @@ void AmoebaReferenceMultipoleForce::setNonbondedMethod(AmoebaReferenceMultipoleF ...@@ -107,7 +107,7 @@ void AmoebaReferenceMultipoleForce::setNonbondedMethod(AmoebaReferenceMultipoleF
_nonbondedMethod = nonbondedMethod; _nonbondedMethod = nonbondedMethod;
} }
AmoebaReferenceMultipoleForce::PolarizationType AmoebaReferenceMultipoleForce::getPolarizationType() const AmoebaReferenceMultipoleForce::PolarizationType AmoebaReferenceMultipoleForce::getPolarizationType() const
{ {
return _polarizationType; return _polarizationType;
} }
...@@ -117,7 +117,7 @@ void AmoebaReferenceMultipoleForce::setPolarizationType(AmoebaReferenceMultipole ...@@ -117,7 +117,7 @@ void AmoebaReferenceMultipoleForce::setPolarizationType(AmoebaReferenceMultipole
_polarizationType = polarizationType; _polarizationType = polarizationType;
} }
int AmoebaReferenceMultipoleForce::getMutualInducedDipoleConverged() const int AmoebaReferenceMultipoleForce::getMutualInducedDipoleConverged() const
{ {
return _mutualInducedDipoleConverged; return _mutualInducedDipoleConverged;
} }
...@@ -127,7 +127,7 @@ void AmoebaReferenceMultipoleForce::setMutualInducedDipoleConverged(int mutualIn ...@@ -127,7 +127,7 @@ void AmoebaReferenceMultipoleForce::setMutualInducedDipoleConverged(int mutualIn
_mutualInducedDipoleConverged = mutualInducedDipoleConverged; _mutualInducedDipoleConverged = mutualInducedDipoleConverged;
} }
int AmoebaReferenceMultipoleForce::getMutualInducedDipoleIterations() const int AmoebaReferenceMultipoleForce::getMutualInducedDipoleIterations() const
{ {
return _mutualInducedDipoleIterations; return _mutualInducedDipoleIterations;
} }
...@@ -137,7 +137,7 @@ void AmoebaReferenceMultipoleForce::setMutualInducedDipoleIterations(int mutualI ...@@ -137,7 +137,7 @@ void AmoebaReferenceMultipoleForce::setMutualInducedDipoleIterations(int mutualI
_mutualInducedDipoleIterations = mutualInducedDipoleIterations; _mutualInducedDipoleIterations = mutualInducedDipoleIterations;
} }
RealOpenMM AmoebaReferenceMultipoleForce::getMutualInducedDipoleEpsilon() const RealOpenMM AmoebaReferenceMultipoleForce::getMutualInducedDipoleEpsilon() const
{ {
return _mutualInducedDipoleEpsilon; return _mutualInducedDipoleEpsilon;
} }
...@@ -147,7 +147,7 @@ void AmoebaReferenceMultipoleForce::setMutualInducedDipoleEpsilon(RealOpenMM mut ...@@ -147,7 +147,7 @@ void AmoebaReferenceMultipoleForce::setMutualInducedDipoleEpsilon(RealOpenMM mut
_mutualInducedDipoleEpsilon = mutualInducedDipoleEpsilon; _mutualInducedDipoleEpsilon = mutualInducedDipoleEpsilon;
} }
int AmoebaReferenceMultipoleForce::getMaximumMutualInducedDipoleIterations() const int AmoebaReferenceMultipoleForce::getMaximumMutualInducedDipoleIterations() const
{ {
return _maximumMutualInducedDipoleIterations; return _maximumMutualInducedDipoleIterations;
} }
...@@ -157,7 +157,7 @@ void AmoebaReferenceMultipoleForce::setMaximumMutualInducedDipoleIterations(int ...@@ -157,7 +157,7 @@ void AmoebaReferenceMultipoleForce::setMaximumMutualInducedDipoleIterations(int
_maximumMutualInducedDipoleIterations = maximumMutualInducedDipoleIterations; _maximumMutualInducedDipoleIterations = maximumMutualInducedDipoleIterations;
} }
RealOpenMM AmoebaReferenceMultipoleForce::getMutualInducedDipoleTargetEpsilon() const RealOpenMM AmoebaReferenceMultipoleForce::getMutualInducedDipoleTargetEpsilon() const
{ {
return _mutualInducedDipoleTargetEpsilon; return _mutualInducedDipoleTargetEpsilon;
} }
...@@ -172,7 +172,7 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector< ...@@ -172,7 +172,7 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector<
/* Setup for scaling maps: /* Setup for scaling maps:
* *
* _scaleMaps[particleIndex][ScaleType] = map, where map[covalentIndex] = scaleFactor * _scaleMaps[particleIndex][ScaleType] = map, where map[covalentIndex] = scaleFactor
* _maxScaleIndex[particleIndex] = max covalent index for particleIndex * _maxScaleIndex[particleIndex] = max covalent index for particleIndex
* *
* multipoleParticleCovalentInfo[ii][jj], jj =0,1,2,3 contains covalent indices (c12, c13, c14, c15) * multipoleParticleCovalentInfo[ii][jj], jj =0,1,2,3 contains covalent indices (c12, c13, c14, c15)
...@@ -208,8 +208,8 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector< ...@@ -208,8 +208,8 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector<
hit = 1; hit = 1;
} }
} }
} }
_scaleMaps[ii][P_SCALE][covalentIndex] = hit ? 0.5*_pScale[jj+1] : _pScale[jj+1]; _scaleMaps[ii][P_SCALE][covalentIndex] = hit ? 0.5*_pScale[jj+1] : _pScale[jj+1];
_scaleMaps[ii][M_SCALE][covalentIndex] = _mScale[jj+1]; _scaleMaps[ii][M_SCALE][covalentIndex] = _mScale[jj+1];
_maxScaleIndex[ii] = _maxScaleIndex[ii] < covalentIndex ? covalentIndex : _maxScaleIndex[ii]; _maxScaleIndex[ii] = _maxScaleIndex[ii] < covalentIndex ? covalentIndex : _maxScaleIndex[ii];
...@@ -231,7 +231,7 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector< ...@@ -231,7 +231,7 @@ void AmoebaReferenceMultipoleForce::setupScaleMaps(const vector< vector< vector<
} }
} }
RealOpenMM AmoebaReferenceMultipoleForce::getMultipoleScaleFactor(unsigned int particleI, unsigned int particleJ, ScaleType scaleType) const RealOpenMM AmoebaReferenceMultipoleForce::getMultipoleScaleFactor(unsigned int particleI, unsigned int particleJ, ScaleType scaleType) const
{ {
MapIntRealOpenMM scaleMap = _scaleMaps[particleI][scaleType]; MapIntRealOpenMM scaleMap = _scaleMaps[particleI][scaleType];
...@@ -243,13 +243,13 @@ RealOpenMM AmoebaReferenceMultipoleForce::getMultipoleScaleFactor(unsigned int p ...@@ -243,13 +243,13 @@ RealOpenMM AmoebaReferenceMultipoleForce::getMultipoleScaleFactor(unsigned int p
} }
} }
void AmoebaReferenceMultipoleForce::getDScaleAndPScale(unsigned int particleI, unsigned int particleJ, RealOpenMM& dScale, RealOpenMM& pScale) const void AmoebaReferenceMultipoleForce::getDScaleAndPScale(unsigned int particleI, unsigned int particleJ, RealOpenMM& dScale, RealOpenMM& pScale) const
{ {
dScale = getMultipoleScaleFactor(particleI, particleJ, D_SCALE); dScale = getMultipoleScaleFactor(particleI, particleJ, D_SCALE);
pScale = getMultipoleScaleFactor(particleI, particleJ, P_SCALE); pScale = getMultipoleScaleFactor(particleI, particleJ, P_SCALE);
} }
void AmoebaReferenceMultipoleForce::getMultipoleScaleFactors(unsigned int particleI, unsigned int particleJ, vector<RealOpenMM>& scaleFactors) const void AmoebaReferenceMultipoleForce::getMultipoleScaleFactors(unsigned int particleI, unsigned int particleJ, vector<RealOpenMM>& scaleFactors) const
{ {
scaleFactors[D_SCALE] = getMultipoleScaleFactor(particleI, particleJ, D_SCALE); scaleFactors[D_SCALE] = getMultipoleScaleFactor(particleI, particleJ, D_SCALE);
scaleFactors[P_SCALE] = getMultipoleScaleFactor(particleI, particleJ, P_SCALE); scaleFactors[P_SCALE] = getMultipoleScaleFactor(particleI, particleJ, P_SCALE);
...@@ -257,7 +257,7 @@ void AmoebaReferenceMultipoleForce::getMultipoleScaleFactors(unsigned int partic ...@@ -257,7 +257,7 @@ void AmoebaReferenceMultipoleForce::getMultipoleScaleFactors(unsigned int partic
scaleFactors[U_SCALE] = getMultipoleScaleFactor(particleI, particleJ, U_SCALE); scaleFactors[U_SCALE] = getMultipoleScaleFactor(particleI, particleJ, U_SCALE);
} }
RealOpenMM AmoebaReferenceMultipoleForce::normalizeRealVec(RealVec& vectorToNormalize) const RealOpenMM AmoebaReferenceMultipoleForce::normalizeRealVec(RealVec& vectorToNormalize) const
{ {
RealOpenMM norm = SQRT(vectorToNormalize.dot(vectorToNormalize)); RealOpenMM norm = SQRT(vectorToNormalize.dot(vectorToNormalize));
if (norm > 0.0) { if (norm > 0.0) {
...@@ -266,26 +266,26 @@ RealOpenMM AmoebaReferenceMultipoleForce::normalizeRealVec(RealVec& vectorToNorm ...@@ -266,26 +266,26 @@ RealOpenMM AmoebaReferenceMultipoleForce::normalizeRealVec(RealVec& vectorToNorm
return norm; return norm;
} }
void AmoebaReferenceMultipoleForce::initializeRealOpenMMVector(vector<RealOpenMM>& vectorToInitialize) const void AmoebaReferenceMultipoleForce::initializeRealOpenMMVector(vector<RealOpenMM>& vectorToInitialize) const
{ {
RealOpenMM zero = 0.0; RealOpenMM zero = 0.0;
vectorToInitialize.resize(_numParticles); vectorToInitialize.resize(_numParticles);
std::fill(vectorToInitialize.begin(), vectorToInitialize.end(), zero); std::fill(vectorToInitialize.begin(), vectorToInitialize.end(), zero);
} }
void AmoebaReferenceMultipoleForce::initializeRealVecVector(vector<RealVec>& vectorToInitialize) const void AmoebaReferenceMultipoleForce::initializeRealVecVector(vector<RealVec>& vectorToInitialize) const
{ {
vectorToInitialize.resize(_numParticles); vectorToInitialize.resize(_numParticles);
RealVec zeroVec(0.0, 0.0, 0.0); RealVec zeroVec(0.0, 0.0, 0.0);
std::fill(vectorToInitialize.begin(), vectorToInitialize.end(), zeroVec); std::fill(vectorToInitialize.begin(), vectorToInitialize.end(), zeroVec);
} }
void AmoebaReferenceMultipoleForce::copyRealVecVector(const vector<OpenMM::RealVec>& inputVector, vector<OpenMM::RealVec>& outputVector) const void AmoebaReferenceMultipoleForce::copyRealVecVector(const vector<OpenMM::RealVec>& inputVector, vector<OpenMM::RealVec>& outputVector) const
{ {
outputVector.resize(inputVector.size()); outputVector.resize(inputVector.size());
for (unsigned int ii = 0; ii < inputVector.size(); ii++) { for (unsigned int ii = 0; ii < inputVector.size(); ii++) {
outputVector[ii] = inputVector[ii]; outputVector[ii] = inputVector[ii];
} }
} }
void AmoebaReferenceMultipoleForce::loadParticleData(const vector<RealVec>& particlePositions, void AmoebaReferenceMultipoleForce::loadParticleData(const vector<RealVec>& particlePositions,
...@@ -295,9 +295,9 @@ void AmoebaReferenceMultipoleForce::loadParticleData(const vector<RealVec>& part ...@@ -295,9 +295,9 @@ void AmoebaReferenceMultipoleForce::loadParticleData(const vector<RealVec>& part
const vector<RealOpenMM>& tholes, const vector<RealOpenMM>& tholes,
const vector<RealOpenMM>& dampingFactors, const vector<RealOpenMM>& dampingFactors,
const vector<RealOpenMM>& polarity, const vector<RealOpenMM>& polarity,
vector<MultipoleParticleData>& particleData) const vector<MultipoleParticleData>& particleData) const
{ {
particleData.resize(_numParticles); particleData.resize(_numParticles);
for (unsigned int ii = 0; ii < _numParticles; ii++) { for (unsigned int ii = 0; ii < _numParticles; ii++) {
...@@ -343,8 +343,8 @@ void AmoebaReferenceMultipoleForce::zeroFixedMultipoleFields() ...@@ -343,8 +343,8 @@ void AmoebaReferenceMultipoleForce::zeroFixedMultipoleFields()
} }
void AmoebaReferenceMultipoleForce::checkChiralCenterAtParticle(MultipoleParticleData& particleI, int axisType, void AmoebaReferenceMultipoleForce::checkChiralCenterAtParticle(MultipoleParticleData& particleI, int axisType,
MultipoleParticleData& particleZ, MultipoleParticleData& particleX, MultipoleParticleData& particleZ, MultipoleParticleData& particleX,
MultipoleParticleData& particleY) const MultipoleParticleData& particleY) const
{ {
if (axisType == AmoebaMultipoleForce::ZThenX || particleY.particleIndex == -1) { if (axisType == AmoebaMultipoleForce::ZThenX || particleY.particleIndex == -1) {
...@@ -372,7 +372,7 @@ void AmoebaReferenceMultipoleForce::checkChiral(vector<MultipoleParticleData>& p ...@@ -372,7 +372,7 @@ void AmoebaReferenceMultipoleForce::checkChiral(vector<MultipoleParticleData>& p
const vector<int>& multipoleAtomXs, const vector<int>& multipoleAtomXs,
const vector<int>& multipoleAtomYs, const vector<int>& multipoleAtomYs,
const vector<int>& multipoleAtomZs, const vector<int>& multipoleAtomZs,
const vector<int>& axisTypes) const const vector<int>& axisTypes) const
{ {
for (unsigned int ii = 0; ii < _numParticles; ii++) { for (unsigned int ii = 0; ii < _numParticles; ii++) {
if (multipoleAtomYs[ii] > -1) { if (multipoleAtomYs[ii] > -1) {
...@@ -388,7 +388,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol ...@@ -388,7 +388,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol
const MultipoleParticleData& particleZ, const MultipoleParticleData& particleZ,
const MultipoleParticleData& particleX, const MultipoleParticleData& particleX,
MultipoleParticleData* particleY, MultipoleParticleData* particleY,
int axisType) const int axisType) const
{ {
// handle case where rotation matrix is identity (e.g. single ion) // handle case where rotation matrix is identity (e.g. single ion)
...@@ -403,25 +403,25 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol ...@@ -403,25 +403,25 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol
RealVec vectorX = particleX.position - particleI.position; RealVec vectorX = particleX.position - particleI.position;
normalizeRealVec(vectorZ); normalizeRealVec(vectorZ);
// branch based on axis type // branch based on axis type
if (axisType == AmoebaMultipoleForce::Bisector) { if (axisType == AmoebaMultipoleForce::Bisector) {
// bisector // bisector
// dx = dx1 + dx2 (in TINKER code) // dx = dx1 + dx2 (in TINKER code)
normalizeRealVec(vectorX); normalizeRealVec(vectorX);
vectorZ += vectorX; vectorZ += vectorX;
normalizeRealVec(vectorZ); normalizeRealVec(vectorZ);
} else if (axisType == AmoebaMultipoleForce::ZBisect) { } else if (axisType == AmoebaMultipoleForce::ZBisect) {
// z-bisect // z-bisect
// dx = dx1 + dx2 (in TINKER code) // dx = dx1 + dx2 (in TINKER code)
normalizeRealVec(vectorX); normalizeRealVec(vectorX);
vectorY = particleY->position - particleI.position; vectorY = particleY->position - particleI.position;
...@@ -429,13 +429,13 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol ...@@ -429,13 +429,13 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol
vectorX += vectorY; vectorX += vectorY;
normalizeRealVec(vectorX); normalizeRealVec(vectorX);
} else if (axisType == AmoebaMultipoleForce::ThreeFold) { } else if (axisType == AmoebaMultipoleForce::ThreeFold) {
// 3-fold // 3-fold
// dx = dx1 + dx2 + dx3 (in TINKER code) // dx = dx1 + dx2 + dx3 (in TINKER code)
normalizeRealVec(vectorX); normalizeRealVec(vectorX);
vectorY = particleY->position - particleI.position; vectorY = particleY->position - particleI.position;
...@@ -443,15 +443,15 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol ...@@ -443,15 +443,15 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol
vectorZ += vectorX + vectorY; vectorZ += vectorX + vectorY;
normalizeRealVec(vectorZ); normalizeRealVec(vectorZ);
} else if (axisType == AmoebaMultipoleForce::ZOnly) { } else if (axisType == AmoebaMultipoleForce::ZOnly) {
// z-only // z-only
vectorX = RealVec(0.1, 0.1, 0.1); vectorX = RealVec(0.1, 0.1, 0.1);
} }
RealOpenMM dot = vectorZ.dot(vectorX); RealOpenMM dot = vectorZ.dot(vectorX);
vectorX -= vectorZ*dot; vectorX -= vectorZ*dot;
...@@ -461,7 +461,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol ...@@ -461,7 +461,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol
RealVec rotationMatrix[3]; RealVec rotationMatrix[3];
rotationMatrix[0] = vectorX; rotationMatrix[0] = vectorX;
rotationMatrix[1] = vectorY; rotationMatrix[1] = vectorY;
rotationMatrix[2] = vectorZ; rotationMatrix[2] = vectorZ;
RealVec labDipole; RealVec labDipole;
for (int ii = 0; ii < 3; ii++) { for (int ii = 0; ii < 3; ii++) {
...@@ -471,7 +471,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol ...@@ -471,7 +471,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol
} }
} }
particleI.dipole = labDipole; particleI.dipole = labDipole;
RealOpenMM mPole[3][3]; RealOpenMM mPole[3][3];
RealOpenMM rPole[3][3] = { { 0.0, 0.0, 0.0 }, RealOpenMM rPole[3][3] = { { 0.0, 0.0, 0.0 },
{ 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 },
...@@ -488,7 +488,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol ...@@ -488,7 +488,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol
mPole[2][0] = particleI.quadrupole[QXZ]; mPole[2][0] = particleI.quadrupole[QXZ];
mPole[2][1] = particleI.quadrupole[QYZ]; mPole[2][1] = particleI.quadrupole[QYZ];
mPole[2][2] = particleI.quadrupole[QZZ]; mPole[2][2] = particleI.quadrupole[QZZ];
for (int ii = 0; ii < 3; ii++) { for (int ii = 0; ii < 3; ii++) {
for (int jj = ii; jj < 3; jj++) { for (int jj = ii; jj < 3; jj++) {
for (int kk = 0; kk < 3; kk++) { for (int kk = 0; kk < 3; kk++) {
...@@ -498,7 +498,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol ...@@ -498,7 +498,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle( Multipol
} }
} }
} }
particleI.quadrupole[QXX] = rPole[0][0]; particleI.quadrupole[QXX] = rPole[0][0];
particleI.quadrupole[QXY] = rPole[0][1]; particleI.quadrupole[QXY] = rPole[0][1];
particleI.quadrupole[QXZ] = rPole[0][2]; particleI.quadrupole[QXZ] = rPole[0][2];
...@@ -636,7 +636,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrix(vector<MultipoleParticle ...@@ -636,7 +636,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrix(vector<MultipoleParticle
const vector<int>& multipoleAtomXs, const vector<int>& multipoleAtomXs,
const vector<int>& multipoleAtomYs, const vector<int>& multipoleAtomYs,
const vector<int>& multipoleAtomZs, const vector<int>& multipoleAtomZs,
const vector<int>& axisTypes) const const vector<int>& axisTypes) const
{ {
for (unsigned int ii = 0; ii < _numParticles; ii++) { for (unsigned int ii = 0; ii < _numParticles; ii++) {
...@@ -649,19 +649,19 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrix(vector<MultipoleParticle ...@@ -649,19 +649,19 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrix(vector<MultipoleParticle
void AmoebaReferenceMultipoleForce::getAndScaleInverseRs(RealOpenMM dampI, RealOpenMM dampJ, void AmoebaReferenceMultipoleForce::getAndScaleInverseRs(RealOpenMM dampI, RealOpenMM dampJ,
RealOpenMM tholeI, RealOpenMM tholeJ, RealOpenMM tholeI, RealOpenMM tholeJ,
RealOpenMM r, vector<RealOpenMM>& rrI) const RealOpenMM r, vector<RealOpenMM>& rrI) const
{ {
RealOpenMM rI = 1.0/r; RealOpenMM rI = 1.0/r;
RealOpenMM r2I = rI*rI; RealOpenMM r2I = rI*rI;
rrI[0] = rI*r2I; rrI[0] = rI*r2I;
RealOpenMM constantFactor = 3.0; RealOpenMM constantFactor = 3.0;
for (unsigned int ii = 1; ii < rrI.size(); ii++) { for (unsigned int ii = 1; ii < rrI.size(); ii++) {
rrI[ii] = constantFactor*rrI[ii-1]*r2I; rrI[ii] = constantFactor*rrI[ii-1]*r2I;
constantFactor += 2.0; constantFactor += 2.0;
} }
RealOpenMM damp = dampI*dampJ; RealOpenMM damp = dampI*dampJ;
if (damp != 0.0) { if (damp != 0.0) {
RealOpenMM pgamma = tholeI < tholeJ ? tholeI : tholeJ; RealOpenMM pgamma = tholeI < tholeJ ? tholeI : tholeJ;
...@@ -669,7 +669,7 @@ void AmoebaReferenceMultipoleForce::getAndScaleInverseRs(RealOpenMM dampI, RealO ...@@ -669,7 +669,7 @@ void AmoebaReferenceMultipoleForce::getAndScaleInverseRs(RealOpenMM dampI, RealO
ratio = ratio*ratio*ratio; ratio = ratio*ratio*ratio;
damp = -pgamma*ratio; damp = -pgamma*ratio;
if (damp > -50.0) { if (damp > -50.0) {
RealOpenMM dampExp = EXP(damp); RealOpenMM dampExp = EXP(damp);
rrI[0] *= 1.0 - dampExp; rrI[0] *= 1.0 - dampExp;
...@@ -683,7 +683,7 @@ void AmoebaReferenceMultipoleForce::getAndScaleInverseRs(RealOpenMM dampI, RealO ...@@ -683,7 +683,7 @@ void AmoebaReferenceMultipoleForce::getAndScaleInverseRs(RealOpenMM dampI, RealO
void AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(const MultipoleParticleData& particleI, void AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(const MultipoleParticleData& particleI,
const MultipoleParticleData& particleJ, const MultipoleParticleData& particleJ,
RealOpenMM dScale, RealOpenMM pScale) RealOpenMM dScale, RealOpenMM pScale)
{ {
if (particleI.particleIndex == particleJ.particleIndex) if (particleI.particleIndex == particleJ.particleIndex)
...@@ -693,9 +693,9 @@ void AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(const Mu ...@@ -693,9 +693,9 @@ void AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(const Mu
RealOpenMM r = SQRT(deltaR.dot(deltaR)); RealOpenMM r = SQRT(deltaR.dot(deltaR));
vector<RealOpenMM> rrI(3); vector<RealOpenMM> rrI(3);
// get scaling factors, if needed // get scaling factors, if needed
getAndScaleInverseRs(particleI.dampingFactor, particleJ.dampingFactor, particleI.thole, particleJ.thole, r, rrI); getAndScaleInverseRs(particleI.dampingFactor, particleJ.dampingFactor, particleI.thole, particleJ.thole, r, rrI);
RealOpenMM rr3 = rrI[0]; RealOpenMM rr3 = rrI[0];
...@@ -710,8 +710,8 @@ void AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(const Mu ...@@ -710,8 +710,8 @@ void AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(const Mu
qDotDelta[1] = deltaR[0]*particleJ.quadrupole[QXY] + deltaR[1]*particleJ.quadrupole[QYY] + deltaR[2]*particleJ.quadrupole[QYZ]; qDotDelta[1] = deltaR[0]*particleJ.quadrupole[QXY] + deltaR[1]*particleJ.quadrupole[QYY] + deltaR[2]*particleJ.quadrupole[QYZ];
qDotDelta[2] = deltaR[0]*particleJ.quadrupole[QXZ] + deltaR[1]*particleJ.quadrupole[QYZ] + deltaR[2]*particleJ.quadrupole[QZZ]; qDotDelta[2] = deltaR[0]*particleJ.quadrupole[QXZ] + deltaR[1]*particleJ.quadrupole[QYZ] + deltaR[2]*particleJ.quadrupole[QZZ];
RealOpenMM dipoleDelta = particleJ.dipole.dot(deltaR); RealOpenMM dipoleDelta = particleJ.dipole.dot(deltaR);
RealOpenMM qdpoleDelta = qDotDelta.dot(deltaR); RealOpenMM qdpoleDelta = qDotDelta.dot(deltaR);
RealOpenMM factor = rr3*particleJ.charge - rr5*dipoleDelta + rr7*qdpoleDelta; RealOpenMM factor = rr3*particleJ.charge - rr5*dipoleDelta + rr7*qdpoleDelta;
RealVec field = deltaR*factor + particleJ.dipole*rr3 - qDotDelta*rr5_2; RealVec field = deltaR*factor + particleJ.dipole*rr3 - qDotDelta*rr5_2;
...@@ -719,17 +719,17 @@ void AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(const Mu ...@@ -719,17 +719,17 @@ void AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(const Mu
unsigned int particleIndex = particleI.particleIndex; unsigned int particleIndex = particleI.particleIndex;
_fixedMultipoleField[particleIndex] -= field*dScale; _fixedMultipoleField[particleIndex] -= field*dScale;
_fixedMultipoleFieldPolar[particleIndex] -= field*pScale; _fixedMultipoleFieldPolar[particleIndex] -= field*pScale;
// field at particle J due multipoles at particle I // field at particle J due multipoles at particle I
qDotDelta[0] = deltaR[0]*particleI.quadrupole[QXX] + deltaR[1]*particleI.quadrupole[QXY] + deltaR[2]*particleI.quadrupole[QXZ]; qDotDelta[0] = deltaR[0]*particleI.quadrupole[QXX] + deltaR[1]*particleI.quadrupole[QXY] + deltaR[2]*particleI.quadrupole[QXZ];
qDotDelta[1] = deltaR[0]*particleI.quadrupole[QXY] + deltaR[1]*particleI.quadrupole[QYY] + deltaR[2]*particleI.quadrupole[QYZ]; qDotDelta[1] = deltaR[0]*particleI.quadrupole[QXY] + deltaR[1]*particleI.quadrupole[QYY] + deltaR[2]*particleI.quadrupole[QYZ];
qDotDelta[2] = deltaR[0]*particleI.quadrupole[QXZ] + deltaR[1]*particleI.quadrupole[QYZ] + deltaR[2]*particleI.quadrupole[QZZ]; qDotDelta[2] = deltaR[0]*particleI.quadrupole[QXZ] + deltaR[1]*particleI.quadrupole[QYZ] + deltaR[2]*particleI.quadrupole[QZZ];
dipoleDelta = particleI.dipole.dot(deltaR); dipoleDelta = particleI.dipole.dot(deltaR);
qdpoleDelta = qDotDelta.dot(deltaR); qdpoleDelta = qDotDelta.dot(deltaR);
factor = rr3*particleI.charge + rr5*dipoleDelta + rr7*qdpoleDelta; factor = rr3*particleI.charge + rr5*dipoleDelta + rr7*qdpoleDelta;
field = deltaR*factor - particleI.dipole*rr3 - qDotDelta*rr5_2; field = deltaR*factor - particleI.dipole*rr3 - qDotDelta*rr5_2;
particleIndex = particleJ.particleIndex; particleIndex = particleJ.particleIndex;
_fixedMultipoleField[particleIndex] += field*dScale; _fixedMultipoleField[particleIndex] += field*dScale;
...@@ -748,7 +748,7 @@ void AmoebaReferenceMultipoleForce::calculateFixedMultipoleField(const vector<Mu ...@@ -748,7 +748,7 @@ void AmoebaReferenceMultipoleForce::calculateFixedMultipoleField(const vector<Mu
for (unsigned int jj = ii; jj < _numParticles; jj++) { for (unsigned int jj = ii; jj < _numParticles; jj++) {
// if site jj is less than max covalent scaling index then get/apply scaling constants // if site jj is less than max covalent scaling index then get/apply scaling constants
// otherwise add unmodified field and fieldPolar to particle fields // otherwise add unmodified field and fieldPolar to particle fields
RealOpenMM dScale, pScale; RealOpenMM dScale, pScale;
if (jj <= _maxScaleIndex[ii]) { if (jj <= _maxScaleIndex[ii]) {
...@@ -775,13 +775,13 @@ void AmoebaReferenceMultipoleForce::initializeInducedDipoles(vector<UpdateInduce ...@@ -775,13 +775,13 @@ void AmoebaReferenceMultipoleForce::initializeInducedDipoles(vector<UpdateInduce
} }
} }
void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxn(unsigned int particleI, void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxn(unsigned int particleI,
unsigned int particleJ, unsigned int particleJ,
RealOpenMM rr3, RealOpenMM rr3,
RealOpenMM rr5, RealOpenMM rr5,
const RealVec& deltaR, const RealVec& deltaR,
const vector<RealVec>& inducedDipole, const vector<RealVec>& inducedDipole,
vector<RealVec>& field) const vector<RealVec>& field) const
{ {
RealOpenMM dDotDelta = rr5*(inducedDipole[particleJ].dot(deltaR)); RealOpenMM dDotDelta = rr5*(inducedDipole[particleJ].dot(deltaR));
...@@ -790,7 +790,7 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxn(unsigned int p ...@@ -790,7 +790,7 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxn(unsigned int p
field[particleJ] += inducedDipole[particleI]*rr3 + deltaR*dDotDelta; field[particleJ] += inducedDipole[particleI]*rr3 + deltaR*dDotDelta;
} }
void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const MultipoleParticleData& particleI, void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const MultipoleParticleData& particleI,
const MultipoleParticleData& particleJ, const MultipoleParticleData& particleJ,
vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields) vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields)
{ {
...@@ -801,13 +801,13 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo ...@@ -801,13 +801,13 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
RealVec deltaR = particleJ.position - particleI.position; RealVec deltaR = particleJ.position - particleI.position;
RealOpenMM r = SQRT(deltaR.dot(deltaR)); RealOpenMM r = SQRT(deltaR.dot(deltaR));
vector<RealOpenMM> rrI(2); vector<RealOpenMM> rrI(2);
getAndScaleInverseRs(particleI.dampingFactor, particleJ.dampingFactor, getAndScaleInverseRs(particleI.dampingFactor, particleJ.dampingFactor,
particleI.thole, particleJ.thole, r, rrI); particleI.thole, particleJ.thole, r, rrI);
RealOpenMM rr3 = -rrI[0]; RealOpenMM rr3 = -rrI[0];
RealOpenMM rr5 = rrI[1]; RealOpenMM rr5 = rrI[1];
for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) { for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) {
calculateInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, rr3, rr5, deltaR, calculateInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, rr3, rr5, deltaR,
*updateInducedDipoleFields[ii].inducedDipoles, updateInducedDipoleFields[ii].inducedDipoleField); *updateInducedDipoleFields[ii].inducedDipoles, updateInducedDipoleFields[ii].inducedDipoleField);
...@@ -816,13 +816,13 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo ...@@ -816,13 +816,13 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
void AmoebaReferenceMultipoleForce::calculateInducedDipoleFields(const vector<MultipoleParticleData>& particleData, vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields) { void AmoebaReferenceMultipoleForce::calculateInducedDipoleFields(const vector<MultipoleParticleData>& particleData, vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields) {
// Initialize the fields to zero. // Initialize the fields to zero.
RealVec zeroVec(0.0, 0.0, 0.0); RealVec zeroVec(0.0, 0.0, 0.0);
for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++)
std::fill(updateInducedDipoleFields[ii].inducedDipoleField.begin(), updateInducedDipoleFields[ii].inducedDipoleField.end(), zeroVec); std::fill(updateInducedDipoleFields[ii].inducedDipoleField.begin(), updateInducedDipoleFields[ii].inducedDipoleField.end(), zeroVec);
// Add fields from all induced dipoles. // Add fields from all induced dipoles.
for (unsigned int ii = 0; ii < particleData.size(); ii++) for (unsigned int ii = 0; ii < particleData.size(); ii++)
for (unsigned int jj = ii; jj < particleData.size(); jj++) for (unsigned int jj = ii; jj < particleData.size(); jj++)
calculateInducedDipolePairIxns(particleData[ii], particleData[jj], updateInducedDipoleFields); calculateInducedDipolePairIxns(particleData[ii], particleData[jj], updateInducedDipoleFields);
...@@ -832,7 +832,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::updateInducedDipoleFields(const vector ...@@ -832,7 +832,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::updateInducedDipoleFields(const vector
vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields) vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields)
{ {
// Calculate the fields coming from induced dipoles. // Calculate the fields coming from induced dipoles.
calculateInducedDipoleFields(particleData, updateInducedDipoleFields); calculateInducedDipoleFields(particleData, updateInducedDipoleFields);
// Update the induced dipoles and calculate the convergence factor, maxEpsilon // Update the induced dipoles and calculate the convergence factor, maxEpsilon
...@@ -843,7 +843,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::updateInducedDipoleFields(const vector ...@@ -843,7 +843,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::updateInducedDipoleFields(const vector
*updateInducedDipoleFields[kk].fixedMultipoleField, *updateInducedDipoleFields[kk].fixedMultipoleField,
updateInducedDipoleFields[kk].inducedDipoleField, updateInducedDipoleFields[kk].inducedDipoleField,
*updateInducedDipoleFields[kk].inducedDipoles); *updateInducedDipoleFields[kk].inducedDipoles);
maxEpsilon = epsilon > maxEpsilon ? epsilon : maxEpsilon; maxEpsilon = epsilon > maxEpsilon ? epsilon : maxEpsilon;
} }
...@@ -878,11 +878,11 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesBySOR(const vector<Mult ...@@ -878,11 +878,11 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesBySOR(const vector<Mult
// loop until (1) induced dipoles are converged or // loop until (1) induced dipoles are converged or
// (2) iterations == max iterations or // (2) iterations == max iterations or
// (3) convergence factor (spsilon) increases // (3) convergence factor (spsilon) increases
while (!done) { while (!done) {
RealOpenMM epsilon = updateInducedDipoleFields(particleData, updateInducedDipoleField); RealOpenMM epsilon = updateInducedDipoleFields(particleData, updateInducedDipoleField);
epsilon = _polarSOR*_debye*SQRT(epsilon/(static_cast<RealOpenMM>(_numParticles))); epsilon = _polarSOR*_debye*SQRT(epsilon/(static_cast<RealOpenMM>(_numParticles)));
if (epsilon < getMutualInducedDipoleTargetEpsilon()) { if (epsilon < getMutualInducedDipoleTargetEpsilon()) {
...@@ -907,11 +907,11 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByDIIS(const vector<Mul ...@@ -907,11 +907,11 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByDIIS(const vector<Mul
int maxPrevious = 20; int maxPrevious = 20;
for (int iteration = 0; ; iteration++) { for (int iteration = 0; ; iteration++) {
// Compute the field from the induced dipoles. // Compute the field from the induced dipoles.
calculateInducedDipoleFields(particleData, updateInducedDipoleField); calculateInducedDipoleFields(particleData, updateInducedDipoleField);
// Record the current dipoles and the errors in them. // Record the current dipoles and the errors in them.
RealOpenMM maxEpsilon = 0; RealOpenMM maxEpsilon = 0;
prevErrors.push_back(vector<RealVec>()); prevErrors.push_back(vector<RealVec>());
prevErrors.back().resize(_numParticles); prevErrors.back().resize(_numParticles);
...@@ -933,9 +933,9 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByDIIS(const vector<Mul ...@@ -933,9 +933,9 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByDIIS(const vector<Mul
maxEpsilon = epsilon; maxEpsilon = epsilon;
} }
maxEpsilon = _debye*SQRT(maxEpsilon/_numParticles); maxEpsilon = _debye*SQRT(maxEpsilon/_numParticles);
// Decide whether to stop or continue iterating. // Decide whether to stop or continue iterating.
if (maxEpsilon < getMutualInducedDipoleTargetEpsilon()) if (maxEpsilon < getMutualInducedDipoleTargetEpsilon())
setMutualInducedDipoleConverged(true); setMutualInducedDipoleConverged(true);
if (maxEpsilon < getMutualInducedDipoleTargetEpsilon() || iteration == getMaximumMutualInducedDipoleIterations()) { if (maxEpsilon < getMutualInducedDipoleTargetEpsilon() || iteration == getMaximumMutualInducedDipoleIterations()) {
...@@ -943,9 +943,9 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByDIIS(const vector<Mul ...@@ -943,9 +943,9 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByDIIS(const vector<Mul
setMutualInducedDipoleIterations(iteration); setMutualInducedDipoleIterations(iteration);
return; return;
} }
// Select the new dipoles. // Select the new dipoles.
if (prevErrors.size() > maxPrevious) { if (prevErrors.size() > maxPrevious) {
prevErrors.erase(prevErrors.begin()); prevErrors.erase(prevErrors.begin());
for (int k = 0; k < numFields; k++) for (int k = 0; k < numFields; k++)
...@@ -972,9 +972,9 @@ void AmoebaReferenceMultipoleForce::computeDIISCoefficients(const vector<vector< ...@@ -972,9 +972,9 @@ void AmoebaReferenceMultipoleForce::computeDIISCoefficients(const vector<vector<
coefficients[0] = 1; coefficients[0] = 1;
return; return;
} }
// Create the DIIS matrix. // Create the DIIS matrix.
int rank = steps+1; int rank = steps+1;
Array2D<double> b(rank, rank); Array2D<double> b(rank, rank);
b[0][0] = 0; b[0][0] = 0;
...@@ -987,9 +987,9 @@ void AmoebaReferenceMultipoleForce::computeDIISCoefficients(const vector<vector< ...@@ -987,9 +987,9 @@ void AmoebaReferenceMultipoleForce::computeDIISCoefficients(const vector<vector<
sum += prevErrors[i][k].dot(prevErrors[j][k]); sum += prevErrors[i][k].dot(prevErrors[j][k]);
b[i+1][j+1] = b[j+1][i+1] = sum; b[i+1][j+1] = b[j+1][i+1] = sum;
} }
// Solve using SVD. Since the right hand side is (-1, 0, 0, 0, ...), this is simpler than the general case. // Solve using SVD. Since the right hand side is (-1, 0, 0, 0, ...), this is simpler than the general case.
JAMA::SVD<double> svd(b); JAMA::SVD<double> svd(b);
Array2D<double> u, v; Array2D<double> u, v;
svd.getU(u); svd.getU(u);
...@@ -1002,7 +1002,7 @@ void AmoebaReferenceMultipoleForce::computeDIISCoefficients(const vector<vector< ...@@ -1002,7 +1002,7 @@ void AmoebaReferenceMultipoleForce::computeDIISCoefficients(const vector<vector<
for (int j = 0; j < effectiveRank; j++) for (int j = 0; j < effectiveRank; j++)
d -= u[0][j]*v[i][j]/s[j]; d -= u[0][j]*v[i][j]/s[j];
coefficients[i-1] = d; coefficients[i-1] = d;
} }
} }
void AmoebaReferenceMultipoleForce::calculateInducedDipoles(const vector<MultipoleParticleData>& particleData) void AmoebaReferenceMultipoleForce::calculateInducedDipoles(const vector<MultipoleParticleData>& particleData)
...@@ -1014,12 +1014,12 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipoles(const vector<Multipo ...@@ -1014,12 +1014,12 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipoles(const vector<Multipo
calculateFixedMultipoleField(particleData); calculateFixedMultipoleField(particleData);
// initialize inducedDipoles // initialize inducedDipoles
// if polarization type is 'Direct', then return after initializing; otherwise // if polarization type is 'Direct', then return after initializing; otherwise
// converge induced dipoles. // converge induced dipoles.
for (unsigned int ii = 0; ii < _numParticles; ii++) { for (unsigned int ii = 0; ii < _numParticles; ii++) {
_fixedMultipoleField[ii] *= particleData[ii].polarity; _fixedMultipoleField[ii] *= particleData[ii].polarity;
_fixedMultipoleFieldPolar[ii] *= particleData[ii].polarity; _fixedMultipoleFieldPolar[ii] *= particleData[ii].polarity;
} }
_inducedDipole.resize(_numParticles); _inducedDipole.resize(_numParticles);
...@@ -1045,7 +1045,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const Mu ...@@ -1045,7 +1045,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const Mu
const MultipoleParticleData& particleK, const MultipoleParticleData& particleK,
const vector<RealOpenMM>& scalingFactors, const vector<RealOpenMM>& scalingFactors,
vector<RealVec>& forces, vector<RealVec>& forces,
vector<RealVec>& torque) const vector<RealVec>& torque) const
{ {
unsigned int iIndex = particleI.particleIndex; unsigned int iIndex = particleI.particleIndex;
unsigned int kIndex = particleK.particleIndex; unsigned int kIndex = particleK.particleIndex;
...@@ -1438,9 +1438,9 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP ...@@ -1438,9 +1438,9 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
const MultipoleParticleData& particleV, const MultipoleParticleData& particleV,
MultipoleParticleData* particleW, MultipoleParticleData* particleW,
int axisType, const Vec3& torque, int axisType, const Vec3& torque,
vector<RealVec>& forces) const vector<RealVec>& forces) const
{ {
static const int U = 0; static const int U = 0;
static const int V = 1; static const int V = 1;
static const int W = 2; static const int W = 2;
...@@ -1454,17 +1454,17 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP ...@@ -1454,17 +1454,17 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
static const int VS = 10; static const int VS = 10;
static const int WS = 11; static const int WS = 11;
static const int LastVectorIndex = 12; static const int LastVectorIndex = 12;
static const int X = 0; static const int X = 0;
static const int Y = 1; static const int Y = 1;
static const int Z = 2; static const int Z = 2;
static const int I = 3; static const int I = 3;
RealOpenMM norms[LastVectorIndex]; RealOpenMM norms[LastVectorIndex];
RealOpenMM angles[LastVectorIndex][2]; RealOpenMM angles[LastVectorIndex][2];
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
// get coordinates of this atom and the z & x axis atoms // get coordinates of this atom and the z & x axis atoms
// compute the vector between the atoms and 1/sqrt(d2), d2 is distance between // compute the vector between the atoms and 1/sqrt(d2), d2 is distance between
// this atom and the axis atom // this atom and the axis atom
...@@ -1486,12 +1486,12 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP ...@@ -1486,12 +1486,12 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
vectorW = vectorU.cross(vectorV); vectorW = vectorU.cross(vectorV);
} }
norms[W] = normalizeRealVec(vectorW); norms[W] = normalizeRealVec(vectorW);
RealVec vectorUV, vectorUW, vectorVW; RealVec vectorUV, vectorUW, vectorVW;
vectorUV = vectorV.cross(vectorU); vectorUV = vectorV.cross(vectorU);
vectorUW = vectorW.cross(vectorU); vectorUW = vectorW.cross(vectorU);
vectorVW = vectorW.cross(vectorV); vectorVW = vectorW.cross(vectorV);
norms[UV] = normalizeRealVec(vectorUV); norms[UV] = normalizeRealVec(vectorUV);
norms[UW] = normalizeRealVec(vectorUW); norms[UW] = normalizeRealVec(vectorUW);
norms[VW] = normalizeRealVec(vectorVW); norms[VW] = normalizeRealVec(vectorVW);
...@@ -1501,7 +1501,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP ...@@ -1501,7 +1501,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
angles[UV][0] = vectorU.dot(vectorV); angles[UV][0] = vectorU.dot(vectorV);
angles[UV][1] = SQRT(1.0 - angles[UV][0]*angles[UV][0]); angles[UV][1] = SQRT(1.0 - angles[UV][0]*angles[UV][0]);
angles[UW][0] = vectorU.dot(vectorW); angles[UW][0] = vectorU.dot(vectorW);
angles[UW][1] = SQRT(1.0 - angles[UW][0]*angles[UW][0]); angles[UW][1] = SQRT(1.0 - angles[UW][0]*angles[UW][0]);
...@@ -1515,26 +1515,26 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP ...@@ -1515,26 +1515,26 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
dphi *= -1.0; dphi *= -1.0;
// branch based on axis type // branch based on axis type
if (axisType == AmoebaMultipoleForce::ZThenX || axisType == AmoebaMultipoleForce::Bisector) { if (axisType == AmoebaMultipoleForce::ZThenX || axisType == AmoebaMultipoleForce::Bisector) {
RealOpenMM factor1; RealOpenMM factor1;
RealOpenMM factor2; RealOpenMM factor2;
RealOpenMM factor3; RealOpenMM factor3;
RealOpenMM factor4; RealOpenMM factor4;
RealOpenMM half = 0.5; RealOpenMM half = 0.5;
factor1 = dphi[V]/(norms[U]*angles[UV][1]); factor1 = dphi[V]/(norms[U]*angles[UV][1]);
factor2 = dphi[W]/(norms[U]); factor2 = dphi[W]/(norms[U]);
factor3 = -dphi[U]/(norms[V]*angles[UV][1]); factor3 = -dphi[U]/(norms[V]*angles[UV][1]);
if (axisType == AmoebaMultipoleForce::Bisector) { if (axisType == AmoebaMultipoleForce::Bisector) {
factor2 *= half; factor2 *= half;
factor4 = half*dphi[W]/(norms[V]); factor4 = half*dphi[W]/(norms[V]);
} else { } else {
factor4 = 0.0; factor4 = 0.0;
} }
for (int ii = 0; ii < 3; ii++) { for (int ii = 0; ii < 3; ii++) {
double forceU = vectorUV[ii]*factor1 + factor2*vectorUW[ii]; double forceU = vectorUV[ii]*factor1 + factor2*vectorUW[ii];
forces[particleU.particleIndex][ii] -= forceU; forces[particleU.particleIndex][ii] -= forceU;
...@@ -1547,7 +1547,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP ...@@ -1547,7 +1547,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
} else if (axisType == AmoebaMultipoleForce::ZBisect) { } else if (axisType == AmoebaMultipoleForce::ZBisect) {
RealVec vectorR = vectorV + vectorW; RealVec vectorR = vectorV + vectorW;
RealVec vectorS = vectorU.cross(vectorR); RealVec vectorS = vectorU.cross(vectorR);
norms[R] = normalizeRealVec(vectorR); norms[R] = normalizeRealVec(vectorR);
...@@ -1574,7 +1574,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP ...@@ -1574,7 +1574,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForceForParticle(const MultipoleP
angles[WS][0] = vectorW.dot(vectorS); angles[WS][0] = vectorW.dot(vectorS);
angles[WS][1] = SQRT(1.0 - angles[WS][0]*angles[WS][0]); angles[WS][1] = SQRT(1.0 - angles[WS][0]*angles[WS][0]);
RealVec t1 = vectorV - vectorS*angles[VS][0]; RealVec t1 = vectorV - vectorS*angles[VS][0];
RealVec t2 = vectorW - vectorS*angles[WS][0]; RealVec t2 = vectorW - vectorS*angles[WS][0];
...@@ -1656,7 +1656,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForce(vector<MultipoleParticleDat ...@@ -1656,7 +1656,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForce(vector<MultipoleParticleDat
const vector<int>& multipoleAtomZs, const vector<int>& multipoleAtomZs,
const vector<int>& axisTypes, const vector<int>& axisTypes,
vector<RealVec>& torques, vector<RealVec>& torques,
vector<RealVec>& forces) const vector<RealVec>& forces) const
{ {
// map torques to forces // map torques to forces
...@@ -1666,7 +1666,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForce(vector<MultipoleParticleDat ...@@ -1666,7 +1666,7 @@ void AmoebaReferenceMultipoleForce::mapTorqueToForce(vector<MultipoleParticleDat
mapTorqueToForceForParticle(particleData[ii], mapTorqueToForceForParticle(particleData[ii],
particleData[multipoleAtomZs[ii]], particleData[multipoleAtomXs[ii]], particleData[multipoleAtomZs[ii]], particleData[multipoleAtomXs[ii]],
multipoleAtomYs[ii] > -1 ? &particleData[multipoleAtomYs[ii]] : NULL, multipoleAtomYs[ii] > -1 ? &particleData[multipoleAtomYs[ii]] : NULL,
axisTypes[ii], torques[ii], forces); axisTypes[ii], torques[ii], forces);
} }
} }
} }
...@@ -1680,7 +1680,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostatic(const vector<Mu ...@@ -1680,7 +1680,7 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostatic(const vector<Mu
vector<RealOpenMM> scaleFactors(LAST_SCALE_TYPE_INDEX); vector<RealOpenMM> scaleFactors(LAST_SCALE_TYPE_INDEX);
for (unsigned int kk = 0; kk < scaleFactors.size(); kk++) { for (unsigned int kk = 0; kk < scaleFactors.size(); kk++) {
scaleFactors[kk] = 1.0; scaleFactors[kk] = 1.0;
} }
// main loop over particle pairs // main loop over particle pairs
...@@ -1878,7 +1878,7 @@ void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments(const ...@@ -1878,7 +1878,7 @@ void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments(const
} }
// convert the quadrupole from traced to traceless form // convert the quadrupole from traced to traceless form
outputMultipoleMoments.resize(13); outputMultipoleMoments.resize(13);
RealOpenMM qave = (xxqdp + yyqdp + zzqdp)/3.0; RealOpenMM qave = (xxqdp + yyqdp + zzqdp)/3.0;
outputMultipoleMoments[4] = 0.5*(xxqdp-qave); outputMultipoleMoments[4] = 0.5*(xxqdp-qave);
...@@ -1901,7 +1901,7 @@ void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments(const ...@@ -1901,7 +1901,7 @@ void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments(const
outputMultipoleMoments[7] = outputMultipoleMoments[5]; outputMultipoleMoments[7] = outputMultipoleMoments[5];
outputMultipoleMoments[10] = outputMultipoleMoments[6]; outputMultipoleMoments[10] = outputMultipoleMoments[6];
outputMultipoleMoments[11] = outputMultipoleMoments[9]; outputMultipoleMoments[11] = outputMultipoleMoments[9];
RealOpenMM debye = 4.80321; RealOpenMM debye = 4.80321;
outputMultipoleMoments[0] = netchg; outputMultipoleMoments[0] = netchg;
...@@ -1910,16 +1910,16 @@ void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments(const ...@@ -1910,16 +1910,16 @@ void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments(const
outputMultipoleMoments[1] = dpl[0]; outputMultipoleMoments[1] = dpl[0];
outputMultipoleMoments[2] = dpl[1]; outputMultipoleMoments[2] = dpl[1];
outputMultipoleMoments[3] = dpl[2]; outputMultipoleMoments[3] = dpl[2];
debye *= 3.0; debye *= 3.0;
for (unsigned int ii = 4; ii < 13; ii++) { for (unsigned int ii = 4; ii < 13; ii++) {
outputMultipoleMoments[ii] *= 100.0*debye; outputMultipoleMoments[ii] *= 100.0*debye;
} }
} }
RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPotentialForParticleGridPoint(const MultipoleParticleData& particleI, const RealVec& gridPoint) const RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPotentialForParticleGridPoint(const MultipoleParticleData& particleI, const RealVec& gridPoint) const
{ {
RealVec deltaR = particleI.position - gridPoint; RealVec deltaR = particleI.position - gridPoint;
getPeriodicDelta(deltaR); getPeriodicDelta(deltaR);
...@@ -1961,7 +1961,7 @@ void AmoebaReferenceMultipoleForce::calculateElectrostaticPotential(const vector ...@@ -1961,7 +1961,7 @@ void AmoebaReferenceMultipoleForce::calculateElectrostaticPotential(const vector
{ {
// setup, including calculating induced dipoles // setup, including calculating induced dipoles
// initialize potential // initialize potential
// calculate contribution of each particle to potential at grid point // calculate contribution of each particle to potential at grid point
// apply prefactor // apply prefactor
...@@ -1988,14 +1988,14 @@ void AmoebaReferenceMultipoleForce::calculateElectrostaticPotential(const vector ...@@ -1988,14 +1988,14 @@ void AmoebaReferenceMultipoleForce::calculateElectrostaticPotential(const vector
} }
AmoebaReferenceMultipoleForce::UpdateInducedDipoleFieldStruct::UpdateInducedDipoleFieldStruct(vector<OpenMM::RealVec>& inputFixed_E_Field, vector<OpenMM::RealVec>& inputInducedDipoles) : AmoebaReferenceMultipoleForce::UpdateInducedDipoleFieldStruct::UpdateInducedDipoleFieldStruct(vector<OpenMM::RealVec>& inputFixed_E_Field, vector<OpenMM::RealVec>& inputInducedDipoles) :
fixedMultipoleField(&inputFixed_E_Field), inducedDipoles(&inputInducedDipoles) { fixedMultipoleField(&inputFixed_E_Field), inducedDipoles(&inputInducedDipoles) {
inducedDipoleField.resize(fixedMultipoleField->size()); inducedDipoleField.resize(fixedMultipoleField->size());
} }
AmoebaReferenceGeneralizedKirkwoodMultipoleForce::AmoebaReferenceGeneralizedKirkwoodMultipoleForce(AmoebaReferenceGeneralizedKirkwoodForce* amoebaReferenceGeneralizedKirkwoodForce) : AmoebaReferenceGeneralizedKirkwoodMultipoleForce::AmoebaReferenceGeneralizedKirkwoodMultipoleForce(AmoebaReferenceGeneralizedKirkwoodForce* amoebaReferenceGeneralizedKirkwoodForce) :
AmoebaReferenceMultipoleForce(NoCutoff), AmoebaReferenceMultipoleForce(NoCutoff),
_amoebaReferenceGeneralizedKirkwoodForce(amoebaReferenceGeneralizedKirkwoodForce), _amoebaReferenceGeneralizedKirkwoodForce(amoebaReferenceGeneralizedKirkwoodForce),
_gkc(2.455) _gkc(2.455)
{ {
RealOpenMM solventDielectric = _amoebaReferenceGeneralizedKirkwoodForce->getSolventDielectric(); RealOpenMM solventDielectric = _amoebaReferenceGeneralizedKirkwoodForce->getSolventDielectric();
...@@ -2016,33 +2016,33 @@ AmoebaReferenceGeneralizedKirkwoodMultipoleForce::AmoebaReferenceGeneralizedKirk ...@@ -2016,33 +2016,33 @@ AmoebaReferenceGeneralizedKirkwoodMultipoleForce::AmoebaReferenceGeneralizedKirk
for (unsigned int ii = 0; ii < _scaledRadii.size(); ii++) { for (unsigned int ii = 0; ii < _scaledRadii.size(); ii++) {
_scaledRadii[ii] *= _atomicRadii[ii]; _scaledRadii[ii] *= _atomicRadii[ii];
} }
} }
AmoebaReferenceGeneralizedKirkwoodMultipoleForce::~AmoebaReferenceGeneralizedKirkwoodMultipoleForce() AmoebaReferenceGeneralizedKirkwoodMultipoleForce::~AmoebaReferenceGeneralizedKirkwoodMultipoleForce()
{ {
delete _amoebaReferenceGeneralizedKirkwoodForce; delete _amoebaReferenceGeneralizedKirkwoodForce;
}; };
int AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getIncludeCavityTerm() const int AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getIncludeCavityTerm() const
{ {
return _includeCavityTerm; return _includeCavityTerm;
}; };
RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getProbeRadius() const RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getProbeRadius() const
{ {
return _probeRadius; return _probeRadius;
}; };
RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getSurfaceAreaFactor() const RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getSurfaceAreaFactor() const
{ {
return _surfaceAreaFactor; return _surfaceAreaFactor;
}; };
RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getDielectricOffset() const RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::getDielectricOffset() const
{ {
return _dielectricOffset; return _dielectricOffset;
}; };
void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::zeroFixedMultipoleFields() void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::zeroFixedMultipoleFields()
{ {
this->AmoebaReferenceMultipoleForce::zeroFixedMultipoleFields(); this->AmoebaReferenceMultipoleForce::zeroFixedMultipoleFields();
...@@ -2057,7 +2057,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateFixedMultipoleFi ...@@ -2057,7 +2057,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateFixedMultipoleFi
this->AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(particleI, particleJ, dScale, pScale); this->AmoebaReferenceMultipoleForce::calculateFixedMultipoleFieldPairIxn(particleI, particleJ, dScale, pScale);
// get deltaR, R2, and R between 2 atoms // get deltaR, R2, and R between 2 atoms
RealVec deltaR = particleJ.position - particleI.position; RealVec deltaR = particleJ.position - particleI.position;
RealOpenMM r = SQRT(deltaR.dot(deltaR)); RealOpenMM r = SQRT(deltaR.dot(deltaR));
RealOpenMM rb2 = _bornRadii[particleI.particleIndex]*_bornRadii[particleJ.particleIndex]; RealOpenMM rb2 = _bornRadii[particleI.particleIndex]*_bornRadii[particleJ.particleIndex];
...@@ -2067,36 +2067,36 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateFixedMultipoleFi ...@@ -2067,36 +2067,36 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateFixedMultipoleFi
RealOpenMM uxi = particleI.dipole[0]; RealOpenMM uxi = particleI.dipole[0];
RealOpenMM uyi = particleI.dipole[1]; RealOpenMM uyi = particleI.dipole[1];
RealOpenMM uzi = particleI.dipole[2]; RealOpenMM uzi = particleI.dipole[2];
RealOpenMM qxxi = particleI.quadrupole[QXX]; RealOpenMM qxxi = particleI.quadrupole[QXX];
RealOpenMM qxyi = particleI.quadrupole[QXY]; RealOpenMM qxyi = particleI.quadrupole[QXY];
RealOpenMM qxzi = particleI.quadrupole[QXZ]; RealOpenMM qxzi = particleI.quadrupole[QXZ];
RealOpenMM qyyi = particleI.quadrupole[QYY]; RealOpenMM qyyi = particleI.quadrupole[QYY];
RealOpenMM qyzi = particleI.quadrupole[QYZ]; RealOpenMM qyzi = particleI.quadrupole[QYZ];
RealOpenMM qzzi = particleI.quadrupole[QZZ]; RealOpenMM qzzi = particleI.quadrupole[QZZ];
RealOpenMM xr = deltaR[0]; RealOpenMM xr = deltaR[0];
RealOpenMM yr = deltaR[1]; RealOpenMM yr = deltaR[1];
RealOpenMM zr = deltaR[2]; RealOpenMM zr = deltaR[2];
RealOpenMM ck = particleJ.charge; RealOpenMM ck = particleJ.charge;
RealOpenMM xr2 = xr*xr; RealOpenMM xr2 = xr*xr;
RealOpenMM yr2 = yr*yr; RealOpenMM yr2 = yr*yr;
RealOpenMM zr2 = zr*zr; RealOpenMM zr2 = zr*zr;
RealOpenMM r2 = xr2 + yr2 + zr2; RealOpenMM r2 = xr2 + yr2 + zr2;
RealOpenMM uxk = particleJ.dipole[0]; RealOpenMM uxk = particleJ.dipole[0];
RealOpenMM uyk = particleJ.dipole[1]; RealOpenMM uyk = particleJ.dipole[1];
RealOpenMM uzk = particleJ.dipole[2]; RealOpenMM uzk = particleJ.dipole[2];
RealOpenMM qxxk = particleJ.quadrupole[QXX]; RealOpenMM qxxk = particleJ.quadrupole[QXX];
RealOpenMM qxyk = particleJ.quadrupole[QXY]; RealOpenMM qxyk = particleJ.quadrupole[QXY];
RealOpenMM qxzk = particleJ.quadrupole[QXZ]; RealOpenMM qxzk = particleJ.quadrupole[QXZ];
RealOpenMM qyyk = particleJ.quadrupole[QYY]; RealOpenMM qyyk = particleJ.quadrupole[QYY];
RealOpenMM qyzk = particleJ.quadrupole[QYZ]; RealOpenMM qyzk = particleJ.quadrupole[QYZ];
RealOpenMM qzzk = particleJ.quadrupole[QZZ]; RealOpenMM qzzk = particleJ.quadrupole[QZZ];
RealOpenMM expterm = EXP(-r2/(_gkc*rb2)); RealOpenMM expterm = EXP(-r2/(_gkc*rb2));
RealOpenMM expc = expterm/_gkc; RealOpenMM expc = expterm/_gkc;
RealOpenMM dexpc = -2.0/(_gkc*rb2); RealOpenMM dexpc = -2.0/(_gkc*rb2);
...@@ -2284,11 +2284,11 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateFixedMultipoleFi ...@@ -2284,11 +2284,11 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateFixedMultipoleFi
void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipolePairGkIxn(const MultipoleParticleData& particleI, void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipolePairGkIxn(const MultipoleParticleData& particleI,
const MultipoleParticleData& particleJ, const MultipoleParticleData& particleJ,
const vector<OpenMM::RealVec>& inputFields, const vector<OpenMM::RealVec>& inputFields,
vector<OpenMM::RealVec>& outputFields) const vector<OpenMM::RealVec>& outputFields) const
{ {
RealOpenMM a[3][3]; RealOpenMM a[3][3];
RealVec deltaR = particleJ.position - particleI.position; RealVec deltaR = particleJ.position - particleI.position;
RealOpenMM r = SQRT(deltaR.dot(deltaR)); RealOpenMM r = SQRT(deltaR.dot(deltaR));
...@@ -2306,7 +2306,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipolePai ...@@ -2306,7 +2306,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipolePai
RealOpenMM r2 = xr2 + yr2 + zr2; RealOpenMM r2 = xr2 + yr2 + zr2;
RealOpenMM expterm = EXP(-r2/(_gkc*rb2)); RealOpenMM expterm = EXP(-r2/(_gkc*rb2));
RealOpenMM expc = expterm /_gkc; RealOpenMM expc = expterm /_gkc;
RealOpenMM gf2 = 1.0/(r2+rb2*expterm); RealOpenMM gf2 = 1.0/(r2+rb2*expterm);
RealOpenMM gf = SQRT(gf2); RealOpenMM gf = SQRT(gf2);
...@@ -2317,7 +2317,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipolePai ...@@ -2317,7 +2317,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipolePai
RealVec duks = inputFields[jIndex]; RealVec duks = inputFields[jIndex];
// reaction potential auxiliary terms // reaction potential auxiliary terms
a[1][0] = -gf3; a[1][0] = -gf3;
a[2][0] = 3.0*gf5; a[2][0] = 3.0*gf5;
...@@ -2325,22 +2325,22 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipolePai ...@@ -2325,22 +2325,22 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipolePai
RealOpenMM expc1 = 1.0 - expc; RealOpenMM expc1 = 1.0 - expc;
a[1][1] = expc1*a[2][0]; a[1][1] = expc1*a[2][0];
// unweighted dipole reaction potential gradient tensor // unweighted dipole reaction potential gradient tensor
RealVec gux, guy, guz; RealVec gux, guy, guz;
gux[0] = (a[1][0] + xr2*a[1][1]); gux[0] = (a[1][0] + xr2*a[1][1]);
gux[1] = xr*yr*a[1][1]; gux[1] = xr*yr*a[1][1];
gux[2] = xr*zr*a[1][1]; gux[2] = xr*zr*a[1][1];
guy[0] = gux[1]; guy[0] = gux[1];
guy[1] = (a[1][0] + yr2*a[1][1]); guy[1] = (a[1][0] + yr2*a[1][1]);
guy[2] = yr*zr*a[1][1]; guy[2] = yr*zr*a[1][1];
guz[0] = gux[2]; guz[0] = gux[2];
guz[1] = guy[2]; guz[1] = guy[2];
guz[2] = (a[1][0] + zr2*a[1][1]); guz[2] = (a[1][0] + zr2*a[1][1]);
outputFields[iIndex][0] += _fd*duks.dot(gux); outputFields[iIndex][0] += _fd*duks.dot(gux);
outputFields[iIndex][1] += _fd*duks.dot(guy); outputFields[iIndex][1] += _fd*duks.dot(guy);
outputFields[iIndex][2] += _fd*duks.dot(guz); outputFields[iIndex][2] += _fd*duks.dot(guz);
...@@ -2387,8 +2387,8 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipoles(c ...@@ -2387,8 +2387,8 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateInducedDipoles(c
for (unsigned int ii = 0; ii < _numParticles; ii++) { for (unsigned int ii = 0; ii < _numParticles; ii++) {
_fixedMultipoleField[ii] *= particleData[ii].polarity; _fixedMultipoleField[ii] *= particleData[ii].polarity;
_fixedMultipoleFieldPolar[ii] *= particleData[ii].polarity; _fixedMultipoleFieldPolar[ii] *= particleData[ii].polarity;
_gkField[ii] *= particleData[ii].polarity; _gkField[ii] *= particleData[ii].polarity;
_inducedDipole[ii] = _fixedMultipoleField[ii]; _inducedDipole[ii] = _fixedMultipoleField[ii];
_inducedDipolePolar[ii] = _fixedMultipoleFieldPolar[ii]; _inducedDipolePolar[ii] = _fixedMultipoleFieldPolar[ii];
...@@ -2418,7 +2418,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodPa ...@@ -2418,7 +2418,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodPa
const MultipoleParticleData& particleJ, const MultipoleParticleData& particleJ,
vector<RealVec>& forces, vector<RealVec>& forces,
vector<RealVec>& torques, vector<RealVec>& torques,
vector<RealOpenMM>& dBorn) const vector<RealOpenMM>& dBorn) const
{ {
RealOpenMM e,ei; RealOpenMM e,ei;
...@@ -2469,9 +2469,9 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodPa ...@@ -2469,9 +2469,9 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodPa
szi = _inducedDipoleS[iIndex][2] + _inducedDipolePolarS[iIndex][2]; szi = _inducedDipoleS[iIndex][2] + _inducedDipolePolarS[iIndex][2];
// decide whether to compute the current interaction // decide whether to compute the current interaction
RealVec deltaR = particleJ.position - particleI.position; RealVec deltaR = particleJ.position - particleI.position;
RealOpenMM r = SQRT(deltaR.dot(deltaR)); RealOpenMM r = SQRT(deltaR.dot(deltaR));
xr = deltaR[0]; xr = deltaR[0];
yr = deltaR[1]; yr = deltaR[1];
...@@ -3863,7 +3863,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrosta ...@@ -3863,7 +3863,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrosta
vector<RealOpenMM> scaleFactors(LAST_SCALE_TYPE_INDEX); vector<RealOpenMM> scaleFactors(LAST_SCALE_TYPE_INDEX);
for (unsigned int kk = 0; kk < scaleFactors.size(); kk++) { for (unsigned int kk = 0; kk < scaleFactors.size(); kk++) {
scaleFactors[kk] = 1.0; scaleFactors[kk] = 1.0;
} }
RealOpenMM eDiffEnergy = 0.0; RealOpenMM eDiffEnergy = 0.0;
for (unsigned int ii = 0; ii < particleData.size(); ii++) { for (unsigned int ii = 0; ii < particleData.size(); ii++) {
...@@ -3889,7 +3889,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrosta ...@@ -3889,7 +3889,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateElectrosta
} }
void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateGrycukChainRulePairIxn(const MultipoleParticleData& particleI, const MultipoleParticleData& particleJ, void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateGrycukChainRulePairIxn(const MultipoleParticleData& particleI, const MultipoleParticleData& particleJ,
const vector<RealOpenMM>& dBorn, vector<RealVec>& forces) const const vector<RealOpenMM>& dBorn, vector<RealVec>& forces) const
{ {
unsigned int iIndex = particleI.particleIndex; unsigned int iIndex = particleI.particleIndex;
unsigned int jIndex = particleJ.particleIndex; unsigned int jIndex = particleJ.particleIndex;
...@@ -3937,7 +3937,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateGrycukChainRuleP ...@@ -3937,7 +3937,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateGrycukChainRuleP
uik = r + sk; uik = r + sk;
uik4 = uik*uik; uik4 = uik*uik;
uik4 = uik4*uik4; uik4 = uik4*uik4;
de -= 0.25*M_PI*(sk2+4.0*sk*r+r2)/ (r2*uik4); de -= 0.25*M_PI*(sk2+4.0*sk*r+r2)/ (r2*uik4);
RealOpenMM dbr = term*de/r; RealOpenMM dbr = term*de/r;
de = dbr*dBorn[iIndex]; de = dbr*dBorn[iIndex];
...@@ -3946,7 +3946,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateGrycukChainRuleP ...@@ -3946,7 +3946,7 @@ void AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateGrycukChainRuleP
forces[jIndex] += deltaR*de; forces[jIndex] += deltaR*de;
} }
RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateCavityTermEnergyAndForces(vector<RealOpenMM>& dBorn) const RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateCavityTermEnergyAndForces(vector<RealOpenMM>& dBorn) const
{ {
RealOpenMM energy = 0.0; RealOpenMM energy = 0.0;
...@@ -3973,7 +3973,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -3973,7 +3973,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
const MultipoleParticleData& particleJ, const MultipoleParticleData& particleJ,
RealOpenMM pscale, RealOpenMM dscale, RealOpenMM pscale, RealOpenMM dscale,
vector<RealVec>& forces, vector<RealVec>& forces,
vector<RealVec>& torques) const vector<RealVec>& torques) const
{ {
static const RealOpenMM uscale = 1.0; static const RealOpenMM uscale = 1.0;
...@@ -4022,7 +4022,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4022,7 +4022,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
RealOpenMM ddsc7_3 = 0.0; RealOpenMM ddsc7_3 = 0.0;
// apply Thole polarization damping to scale factors // apply Thole polarization damping to scale factors
RealOpenMM damp = particleI.dampingFactor*particleJ.dampingFactor; RealOpenMM damp = particleI.dampingFactor*particleJ.dampingFactor;
if (damp != 0.0) { if (damp != 0.0) {
pgamma = particleJ.thole > particleI.thole ? particleI.thole : particleJ.thole; pgamma = particleJ.thole > particleI.thole ? particleI.thole : particleJ.thole;
...@@ -4059,9 +4059,9 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4059,9 +4059,9 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
psc3 = scale3*pscale; psc3 = scale3*pscale;
psc5 = scale5*pscale; psc5 = scale5*pscale;
psc7 = scale7*pscale; psc7 = scale7*pscale;
// construct auxiliary vectors for permanent terms // construct auxiliary vectors for permanent terms
RealOpenMM dixr1 = particleI.dipole[1]*zr - particleI.dipole[2]*yr; RealOpenMM dixr1 = particleI.dipole[1]*zr - particleI.dipole[2]*yr;
RealOpenMM dixr2 = particleI.dipole[2]*xr - particleI.dipole[0]*zr; RealOpenMM dixr2 = particleI.dipole[2]*xr - particleI.dipole[0]*zr;
RealOpenMM dixr3 = particleI.dipole[0]*yr - particleI.dipole[1]*xr; RealOpenMM dixr3 = particleI.dipole[0]*yr - particleI.dipole[1]*xr;
...@@ -4087,14 +4087,14 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4087,14 +4087,14 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
RealOpenMM rxqkr3 = xr*qkr2 - yr*qkr1; RealOpenMM rxqkr3 = xr*qkr2 - yr*qkr1;
// get intermediate variables for permanent energy terms // get intermediate variables for permanent energy terms
RealOpenMM sc3 = particleI.dipole[0]*xr + particleI.dipole[1]*yr + particleI.dipole[2]*zr; RealOpenMM sc3 = particleI.dipole[0]*xr + particleI.dipole[1]*yr + particleI.dipole[2]*zr;
RealOpenMM sc4 = particleJ.dipole[0]*xr + particleJ.dipole[1]*yr + particleJ.dipole[2]*zr; RealOpenMM sc4 = particleJ.dipole[0]*xr + particleJ.dipole[1]*yr + particleJ.dipole[2]*zr;
RealOpenMM sc5 = qir1*xr + qir2*yr + qir3*zr; RealOpenMM sc5 = qir1*xr + qir2*yr + qir3*zr;
RealOpenMM sc6 = qkr1*xr + qkr2*yr + qkr3*zr; RealOpenMM sc6 = qkr1*xr + qkr2*yr + qkr3*zr;
// construct auxiliary vectors for induced terms // construct auxiliary vectors for induced terms
RealOpenMM dixuk1 = particleI.dipole[1]*_inducedDipoleS[jIndex][2] - particleI.dipole[2]*_inducedDipoleS[jIndex][1]; RealOpenMM dixuk1 = particleI.dipole[1]*_inducedDipoleS[jIndex][2] - particleI.dipole[2]*_inducedDipoleS[jIndex][1];
RealOpenMM dixuk2 = particleI.dipole[2]*_inducedDipoleS[jIndex][0] - particleI.dipole[0]*_inducedDipoleS[jIndex][2]; RealOpenMM dixuk2 = particleI.dipole[2]*_inducedDipoleS[jIndex][0] - particleI.dipole[0]*_inducedDipoleS[jIndex][2];
RealOpenMM dixuk3 = particleI.dipole[0]*_inducedDipoleS[jIndex][1] - particleI.dipole[1]*_inducedDipoleS[jIndex][0]; RealOpenMM dixuk3 = particleI.dipole[0]*_inducedDipoleS[jIndex][1] - particleI.dipole[1]*_inducedDipoleS[jIndex][0];
...@@ -4158,7 +4158,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4158,7 +4158,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
RealOpenMM rxqkuip1 = yr*qkuip3 - zr*qkuip2; RealOpenMM rxqkuip1 = yr*qkuip3 - zr*qkuip2;
RealOpenMM rxqkuip2 = zr*qkuip1 - xr*qkuip3; RealOpenMM rxqkuip2 = zr*qkuip1 - xr*qkuip3;
RealOpenMM rxqkuip3 = xr*qkuip2 - yr*qkuip1; RealOpenMM rxqkuip3 = xr*qkuip2 - yr*qkuip1;
// get intermediate variables for induction energy terms // get intermediate variables for induction energy terms
RealOpenMM sci1 = _inducedDipoleS[iIndex][0]*particleJ.dipole[0] + _inducedDipoleS[iIndex][1]*particleJ.dipole[1] + RealOpenMM sci1 = _inducedDipoleS[iIndex][0]*particleJ.dipole[0] + _inducedDipoleS[iIndex][1]*particleJ.dipole[1] +
...@@ -4216,7 +4216,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4216,7 +4216,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
RealOpenMM gfi6 = -rr7*(sci3*psc7+scip3*dsc7); RealOpenMM gfi6 = -rr7*(sci3*psc7+scip3*dsc7);
// get the induced force // get the induced force
RealOpenMM ftm2i1 = gfi1*xr + 0.5* RealOpenMM ftm2i1 = gfi1*xr + 0.5*
(- rr3*particleJ.charge*(_inducedDipoleS[iIndex][0]*psc3+_inducedDipolePolarS[iIndex][0]*dsc3) (- rr3*particleJ.charge*(_inducedDipoleS[iIndex][0]*psc3+_inducedDipolePolarS[iIndex][0]*dsc3)
+ rr5*sc4*(_inducedDipoleS[iIndex][0]*psc5+_inducedDipolePolarS[iIndex][0]*dsc5) + rr5*sc4*(_inducedDipoleS[iIndex][0]*psc5+_inducedDipolePolarS[iIndex][0]*dsc5)
...@@ -4231,7 +4231,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4231,7 +4231,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
+ 0.5*gfi4*((qkui1-qiuk1)*psc5 + 0.5*gfi4*((qkui1-qiuk1)*psc5
+ (qkuip1-qiukp1)*dsc5) + (qkuip1-qiukp1)*dsc5)
+ gfi5*qir1 + gfi6*qkr1; + gfi5*qir1 + gfi6*qkr1;
RealOpenMM ftm2i2 = gfi1*yr + 0.5* RealOpenMM ftm2i2 = gfi1*yr + 0.5*
(- rr3*particleJ.charge*(_inducedDipoleS[iIndex][1]*psc3+_inducedDipolePolarS[iIndex][1]*dsc3) (- rr3*particleJ.charge*(_inducedDipoleS[iIndex][1]*psc3+_inducedDipolePolarS[iIndex][1]*dsc3)
+ rr5*sc4*(_inducedDipoleS[iIndex][1]*psc5+_inducedDipolePolarS[iIndex][1]*dsc5) + rr5*sc4*(_inducedDipoleS[iIndex][1]*psc5+_inducedDipolePolarS[iIndex][1]*dsc5)
...@@ -4261,7 +4261,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4261,7 +4261,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
+ 0.5*gfi4*((qkui3-qiuk3)*psc5 + 0.5*gfi4*((qkui3-qiuk3)*psc5
+ (qkuip3-qiukp3)*dsc5) + (qkuip3-qiukp3)*dsc5)
+ gfi5*qir3 + gfi6*qkr3; + gfi5*qir3 + gfi6*qkr3;
// intermediate values needed for partially excluded interactions // intermediate values needed for partially excluded interactions
RealOpenMM fridmp1 = 0.5*(rr3*((gli1+gli6)*pscale RealOpenMM fridmp1 = 0.5*(rr3*((gli1+gli6)*pscale
...@@ -4283,7 +4283,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4283,7 +4283,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
+ rr7*(gli3*pscale+glip3*dscale)*ddsc7_3); + rr7*(gli3*pscale+glip3*dscale)*ddsc7_3);
// get the induced-induced derivative terms // get the induced-induced derivative terms
RealOpenMM findmp1 = 0.5*uscale*(scip2*rr3*ddsc3_1 RealOpenMM findmp1 = 0.5*uscale*(scip2*rr3*ddsc3_1
- rr5*ddsc5_1*(sci3*scip4+scip3*sci4)); - rr5*ddsc5_1*(sci3*scip4+scip3*sci4));
...@@ -4309,12 +4309,12 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4309,12 +4309,12 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
ftm2i1 -= fdir1 - findmp1; ftm2i1 -= fdir1 - findmp1;
ftm2i2 -= fdir2 - findmp2; ftm2i2 -= fdir2 - findmp2;
ftm2i3 -= fdir3 - findmp3; ftm2i3 -= fdir3 - findmp3;
} }
// now perform the torque calculation // now perform the torque calculation
// intermediate terms for torque between multipoles i and j // intermediate terms for torque between multipoles i and j
RealOpenMM gti2 = 0.5*(sci4*psc5+scip4*dsc5)*rr5; RealOpenMM gti2 = 0.5*(sci4*psc5+scip4*dsc5)*rr5;
RealOpenMM gti3 = 0.5*(sci3*psc5+scip3*dsc5)*rr5; RealOpenMM gti3 = 0.5*(sci3*psc5+scip3*dsc5)*rr5;
RealOpenMM gti4 = gfi4; RealOpenMM gti4 = gfi4;
...@@ -4346,7 +4346,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4346,7 +4346,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
RealOpenMM ttm3i3 = -rr3*(dkxui3*psc3+dkxuip3*dsc3)*0.5 RealOpenMM ttm3i3 = -rr3*(dkxui3*psc3+dkxuip3*dsc3)*0.5
+ gti3*dkxr3 - gti4*((uixqkr3+rxqkui3)*psc5 + gti3*dkxr3 - gti4*((uixqkr3+rxqkui3)*psc5
+(uixqkrp3+rxqkuip3)*dsc5)*0.5 - gti6*rxqkr3; +(uixqkrp3+rxqkuip3)*dsc5)*0.5 - gti6*rxqkr3;
// update force and torque // update force and torque
RealVec force, torqueI, torqueJ; RealVec force, torqueI, torqueJ;
...@@ -4610,7 +4610,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED ...@@ -4610,7 +4610,7 @@ RealOpenMM AmoebaReferenceGeneralizedKirkwoodMultipoleForce::calculateKirkwoodED
+ gti3*dkxr3 - gti4*((uixqkr3+rxqkui3)*psc5 + gti3*dkxr3 - gti4*((uixqkr3+rxqkui3)*psc5
+(uixqkrp3+rxqkuip3)*dsc5)*0.5 - gti6*rxqkr3; +(uixqkrp3+rxqkuip3)*dsc5)*0.5 - gti6*rxqkr3;
// update force and torque // update force and torque
force[0] += ftm2i1; force[0] += ftm2i1;
force[1] += ftm2i2; force[1] += ftm2i2;
...@@ -4635,13 +4635,13 @@ const RealOpenMM AmoebaReferencePmeMultipoleForce::SQRT_PI = 1.77245385091; ...@@ -4635,13 +4635,13 @@ const RealOpenMM AmoebaReferencePmeMultipoleForce::SQRT_PI = 1.77245385091;
AmoebaReferencePmeMultipoleForce::AmoebaReferencePmeMultipoleForce() : AmoebaReferencePmeMultipoleForce::AmoebaReferencePmeMultipoleForce() :
AmoebaReferenceMultipoleForce(PME), AmoebaReferenceMultipoleForce(PME),
_cutoffDistance(1.0), _cutoffDistanceSquared(1.0), _cutoffDistance(1.0), _cutoffDistanceSquared(1.0),
_pmeGridSize(0), _totalGridSize(0), _alphaEwald(0.0) _pmeGridSize(0), _totalGridSize(0), _alphaEwald(0.0)
{ {
_fftplan = NULL; _fftplan = NULL;
_pmeGrid = NULL; _pmeGrid = NULL;
_pmeGridDimensions = IntVec(-1, -1, -1); _pmeGridDimensions = IntVec(-1, -1, -1);
} }
AmoebaReferencePmeMultipoleForce::~AmoebaReferencePmeMultipoleForce() AmoebaReferencePmeMultipoleForce::~AmoebaReferencePmeMultipoleForce()
{ {
...@@ -4652,19 +4652,19 @@ AmoebaReferencePmeMultipoleForce::~AmoebaReferencePmeMultipoleForce() ...@@ -4652,19 +4652,19 @@ AmoebaReferencePmeMultipoleForce::~AmoebaReferencePmeMultipoleForce()
delete _pmeGrid; delete _pmeGrid;
} }
}; };
RealOpenMM AmoebaReferencePmeMultipoleForce::getCutoffDistance() const RealOpenMM AmoebaReferencePmeMultipoleForce::getCutoffDistance() const
{ {
return _cutoffDistance; return _cutoffDistance;
}; };
void AmoebaReferencePmeMultipoleForce::setCutoffDistance(RealOpenMM cutoffDistance) void AmoebaReferencePmeMultipoleForce::setCutoffDistance(RealOpenMM cutoffDistance)
{ {
_cutoffDistance = cutoffDistance; _cutoffDistance = cutoffDistance;
_cutoffDistanceSquared = cutoffDistance*cutoffDistance; _cutoffDistanceSquared = cutoffDistance*cutoffDistance;
}; };
RealOpenMM AmoebaReferencePmeMultipoleForce::getAlphaEwald() const RealOpenMM AmoebaReferencePmeMultipoleForce::getAlphaEwald() const
{ {
return _alphaEwald; return _alphaEwald;
}; };
...@@ -4674,7 +4674,7 @@ void AmoebaReferencePmeMultipoleForce::setAlphaEwald(RealOpenMM alphaEwald) ...@@ -4674,7 +4674,7 @@ void AmoebaReferencePmeMultipoleForce::setAlphaEwald(RealOpenMM alphaEwald)
_alphaEwald = alphaEwald; _alphaEwald = alphaEwald;
}; };
void AmoebaReferencePmeMultipoleForce::getPmeGridDimensions(vector<int>& pmeGridDimensions) const void AmoebaReferencePmeMultipoleForce::getPmeGridDimensions(vector<int>& pmeGridDimensions) const
{ {
pmeGridDimensions.resize(3); pmeGridDimensions.resize(3);
...@@ -4687,7 +4687,7 @@ void AmoebaReferencePmeMultipoleForce::getPmeGridDimensions(vector<int>& pmeGrid ...@@ -4687,7 +4687,7 @@ void AmoebaReferencePmeMultipoleForce::getPmeGridDimensions(vector<int>& pmeGrid
void AmoebaReferencePmeMultipoleForce::setPmeGridDimensions(vector<int>& pmeGridDimensions) void AmoebaReferencePmeMultipoleForce::setPmeGridDimensions(vector<int>& pmeGridDimensions)
{ {
if ((pmeGridDimensions[0] == _pmeGridDimensions[0]) && if ((pmeGridDimensions[0] == _pmeGridDimensions[0]) &&
(pmeGridDimensions[1] == _pmeGridDimensions[1]) && (pmeGridDimensions[1] == _pmeGridDimensions[1]) &&
(pmeGridDimensions[2] == _pmeGridDimensions[2])) (pmeGridDimensions[2] == _pmeGridDimensions[2]))
return; return;
...@@ -4761,9 +4761,9 @@ void AmoebaReferencePmeMultipoleForce::initializePmeGrid() ...@@ -4761,9 +4761,9 @@ void AmoebaReferencePmeMultipoleForce::initializePmeGrid()
for (int jj = 0; jj < _totalGridSize; jj++) { for (int jj = 0; jj < _totalGridSize; jj++) {
_pmeGrid[jj].re = _pmeGrid[jj].im = 0.0; _pmeGrid[jj].re = _pmeGrid[jj].im = 0.0;
} }
} }
void AmoebaReferencePmeMultipoleForce::getPeriodicDelta(RealVec& deltaR) const void AmoebaReferencePmeMultipoleForce::getPeriodicDelta(RealVec& deltaR) const
{ {
deltaR -= _periodicBoxVectors[2]*FLOOR(deltaR[2]*_recipBoxVectors[2][2]+0.5); deltaR -= _periodicBoxVectors[2]*FLOOR(deltaR[2]*_recipBoxVectors[2][2]+0.5);
deltaR -= _periodicBoxVectors[1]*FLOOR(deltaR[1]*_recipBoxVectors[1][1]+0.5); deltaR -= _periodicBoxVectors[1]*FLOOR(deltaR[1]*_recipBoxVectors[1][1]+0.5);
...@@ -4773,8 +4773,8 @@ void AmoebaReferencePmeMultipoleForce::getPeriodicDelta(RealVec& deltaR) const ...@@ -4773,8 +4773,8 @@ void AmoebaReferencePmeMultipoleForce::getPeriodicDelta(RealVec& deltaR) const
void AmoebaReferencePmeMultipoleForce::getDampedInverseDistances(const MultipoleParticleData& particleI, void AmoebaReferencePmeMultipoleForce::getDampedInverseDistances(const MultipoleParticleData& particleI,
const MultipoleParticleData& particleJ, const MultipoleParticleData& particleJ,
RealOpenMM dscale, RealOpenMM pscale, RealOpenMM r, RealOpenMM dscale, RealOpenMM pscale, RealOpenMM r,
RealVec& dampedDInverseDistances, RealVec& dampedDInverseDistances,
RealVec& dampedPInverseDistances) const RealVec& dampedPInverseDistances) const
{ {
RealVec scaleFactor = RealVec(1.0, 1.0, 1.0); RealVec scaleFactor = RealVec(1.0, 1.0, 1.0);
...@@ -4792,8 +4792,8 @@ void AmoebaReferencePmeMultipoleForce::getDampedInverseDistances(const Multipole ...@@ -4792,8 +4792,8 @@ void AmoebaReferencePmeMultipoleForce::getDampedInverseDistances(const Multipole
scaleFactor[0] = 1.0 - expdamp; scaleFactor[0] = 1.0 - expdamp;
scaleFactor[1] = 1.0 - expdamp*(1.0-damp); scaleFactor[1] = 1.0 - expdamp*(1.0-damp);
scaleFactor[2] = 1.0 - expdamp*(1.0-damp+(0.6f*damp*damp)); scaleFactor[2] = 1.0 - expdamp*(1.0-damp+(0.6f*damp*damp));
} }
} }
RealVec dampedDScale = scaleFactor*dscale; RealVec dampedDScale = scaleFactor*dscale;
RealOpenMM r2 = r*r; RealOpenMM r2 = r*r;
...@@ -4808,7 +4808,7 @@ void AmoebaReferencePmeMultipoleForce::getDampedInverseDistances(const Multipole ...@@ -4808,7 +4808,7 @@ void AmoebaReferencePmeMultipoleForce::getDampedInverseDistances(const Multipole
dampedPInverseDistances = dampedDInverseDistances; dampedPInverseDistances = dampedDInverseDistances;
} else { } else {
RealVec dampedPScale = scaleFactor*pscale; RealVec dampedPScale = scaleFactor*pscale;
dampedPInverseDistances[0] = (1.0-dampedPScale[0])/r3; dampedPInverseDistances[0] = (1.0-dampedPScale[0])/r3;
dampedPInverseDistances[1] = 3.0*(1.0-dampedPScale[1])/r5; dampedPInverseDistances[1] = 3.0*(1.0-dampedPScale[1])/r5;
dampedPInverseDistances[2] = 15.0*(1.0-dampedPScale[2])/r7; dampedPInverseDistances[2] = 15.0*(1.0-dampedPScale[2])/r7;
} }
...@@ -4976,7 +4976,7 @@ void AmoebaReferencePmeMultipoleForce::calculateFixedMultipoleFieldPairIxn(const ...@@ -4976,7 +4976,7 @@ void AmoebaReferencePmeMultipoleForce::calculateFixedMultipoleFieldPairIxn(const
RealVec qj = RealVec(qxJ.dot(deltaR), qyJ.dot(deltaR), qzJ.dot(deltaR)); RealVec qj = RealVec(qxJ.dot(deltaR), qyJ.dot(deltaR), qzJ.dot(deltaR));
RealOpenMM qjr = qj.dot(deltaR); RealOpenMM qjr = qj.dot(deltaR);
RealVec fim = qj*(2.0*bn2) - particleJ.dipole*bn1 - deltaR*(bn1*particleJ.charge - bn2*djr+bn3*qjr); RealVec fim = qj*(2.0*bn2) - particleJ.dipole*bn1 - deltaR*(bn1*particleJ.charge - bn2*djr+bn3*qjr);
RealVec fjm = qi*(-2.0*bn2) - particleI.dipole*bn1 + deltaR*(bn1*particleI.charge + bn2*dir+bn3*qir); RealVec fjm = qi*(-2.0*bn2) - particleI.dipole*bn1 + deltaR*(bn1*particleI.charge + bn2*dir+bn3*qir);
...@@ -5105,7 +5105,7 @@ void AmoebaReferencePmeMultipoleForce::computeBSplinePoint(vector<RealOpenMM4>& ...@@ -5105,7 +5105,7 @@ void AmoebaReferencePmeMultipoleForce::computeBSplinePoint(vector<RealOpenMM4>&
/** /**
* Compute b-spline coefficients. * Compute b-spline coefficients.
*/ */
void AmoebaReferencePmeMultipoleForce::computeAmoebaBsplines(const vector<MultipoleParticleData>& particleData) void AmoebaReferencePmeMultipoleForce::computeAmoebaBsplines(const vector<MultipoleParticleData>& particleData)
{ {
// get the B-spline coefficients for each multipole site // get the B-spline coefficients for each multipole site
...@@ -5127,7 +5127,7 @@ void AmoebaReferencePmeMultipoleForce::computeAmoebaBsplines(const vector<Multip ...@@ -5127,7 +5127,7 @@ void AmoebaReferencePmeMultipoleForce::computeAmoebaBsplines(const vector<Multip
_thetai[jj][ii*AMOEBA_PME_ORDER+kk] = thetaiTemp[kk]; _thetai[jj][ii*AMOEBA_PME_ORDER+kk] = thetaiTemp[kk];
} }
} }
// Record the grid point. // Record the grid point.
_iGrid[ii] = igrid; _iGrid[ii] = igrid;
...@@ -5136,7 +5136,7 @@ void AmoebaReferencePmeMultipoleForce::computeAmoebaBsplines(const vector<Multip ...@@ -5136,7 +5136,7 @@ void AmoebaReferencePmeMultipoleForce::computeAmoebaBsplines(const vector<Multip
void AmoebaReferencePmeMultipoleForce::transformMultipolesToFractionalCoordinates(const vector<MultipoleParticleData>& particleData) { void AmoebaReferencePmeMultipoleForce::transformMultipolesToFractionalCoordinates(const vector<MultipoleParticleData>& particleData) {
// Build matrices for transforming the dipoles and quadrupoles. // Build matrices for transforming the dipoles and quadrupoles.
RealVec a[3]; RealVec a[3];
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
...@@ -5172,7 +5172,7 @@ void AmoebaReferencePmeMultipoleForce::transformMultipolesToFractionalCoordinate ...@@ -5172,7 +5172,7 @@ void AmoebaReferencePmeMultipoleForce::transformMultipolesToFractionalCoordinate
void AmoebaReferencePmeMultipoleForce::transformPotentialToCartesianCoordinates(const vector<RealOpenMM>& fphi, vector<RealOpenMM>& cphi) const { void AmoebaReferencePmeMultipoleForce::transformPotentialToCartesianCoordinates(const vector<RealOpenMM>& fphi, vector<RealOpenMM>& cphi) const {
// Build a matrix for transforming the potential. // Build a matrix for transforming the potential.
RealVec a[3]; RealVec a[3];
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
...@@ -5194,9 +5194,9 @@ void AmoebaReferencePmeMultipoleForce::transformPotentialToCartesianCoordinates( ...@@ -5194,9 +5194,9 @@ void AmoebaReferencePmeMultipoleForce::transformPotentialToCartesianCoordinates(
b[i][j] += a[index1[i]][index2[j]]*a[index2[i]][index1[j]]; b[i][j] += a[index1[i]][index2[j]]*a[index2[i]][index1[j]];
} }
} }
// Transform the potential. // Transform the potential.
for (int i = 0; i < _numParticles; i++) { for (int i = 0; i < _numParticles; i++) {
cphi[10*i] = fphi[20*i]; cphi[10*i] = fphi[20*i];
cphi[10*i+1] = a[0][0]*fphi[20*i+1] + a[0][1]*fphi[20*i+2] + a[0][2]*fphi[20*i+3]; cphi[10*i+1] = a[0][0]*fphi[20*i+1] + a[0][1]*fphi[20*i+2] + a[0][2]*fphi[20*i+3];
...@@ -5210,18 +5210,18 @@ void AmoebaReferencePmeMultipoleForce::transformPotentialToCartesianCoordinates( ...@@ -5210,18 +5210,18 @@ void AmoebaReferencePmeMultipoleForce::transformPotentialToCartesianCoordinates(
} }
} }
void AmoebaReferencePmeMultipoleForce::spreadFixedMultipolesOntoGrid(const vector<MultipoleParticleData>& particleData) void AmoebaReferencePmeMultipoleForce::spreadFixedMultipolesOntoGrid(const vector<MultipoleParticleData>& particleData)
{ {
transformMultipolesToFractionalCoordinates(particleData); transformMultipolesToFractionalCoordinates(particleData);
// Clear the grid. // Clear the grid.
for (int gridIndex = 0; gridIndex < _totalGridSize; gridIndex++) for (int gridIndex = 0; gridIndex < _totalGridSize; gridIndex++)
_pmeGrid[gridIndex] = t_complex(0, 0); _pmeGrid[gridIndex] = t_complex(0, 0);
// Loop over atoms and spread them on the grid. // Loop over atoms and spread them on the grid.
for (int atomIndex = 0; atomIndex < _numParticles; atomIndex++) { for (int atomIndex = 0; atomIndex < _numParticles; atomIndex++) {
RealOpenMM atomCharge = _transformed[atomIndex].charge; RealOpenMM atomCharge = _transformed[atomIndex].charge;
RealVec atomDipole = RealVec(_transformed[atomIndex].dipole[0], RealVec atomDipole = RealVec(_transformed[atomIndex].dipole[0],
...@@ -5405,19 +5405,19 @@ void AmoebaReferencePmeMultipoleForce::computeFixedPotentialFromGrid() ...@@ -5405,19 +5405,19 @@ void AmoebaReferencePmeMultipoleForce::computeFixedPotentialFromGrid()
void AmoebaReferencePmeMultipoleForce::spreadInducedDipolesOnGrid(const vector<RealVec>& inputInducedDipole, void AmoebaReferencePmeMultipoleForce::spreadInducedDipolesOnGrid(const vector<RealVec>& inputInducedDipole,
const vector<RealVec>& inputInducedDipolePolar) { const vector<RealVec>& inputInducedDipolePolar) {
// Create the matrix to convert from Cartesian to fractional coordinates. // Create the matrix to convert from Cartesian to fractional coordinates.
RealVec cartToFrac[3]; RealVec cartToFrac[3];
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) for (int j = 0; j < 3; j++)
cartToFrac[j][i] = _pmeGridDimensions[j]*_recipBoxVectors[i][j]; cartToFrac[j][i] = _pmeGridDimensions[j]*_recipBoxVectors[i][j];
// Clear the grid. // Clear the grid.
for (int gridIndex = 0; gridIndex < _totalGridSize; gridIndex++) for (int gridIndex = 0; gridIndex < _totalGridSize; gridIndex++)
_pmeGrid[gridIndex] = t_complex(0, 0); _pmeGrid[gridIndex] = t_complex(0, 0);
// Loop over atoms and spread them on the grid. // Loop over atoms and spread them on the grid.
for (int atomIndex = 0; atomIndex < _numParticles; atomIndex++) { for (int atomIndex = 0; atomIndex < _numParticles; atomIndex++) {
RealVec inducedDipole = RealVec(inputInducedDipole[atomIndex][0]*cartToFrac[0][0] + inputInducedDipole[atomIndex][1]*cartToFrac[0][1] + inputInducedDipole[atomIndex][2]*cartToFrac[0][2], RealVec inducedDipole = RealVec(inputInducedDipole[atomIndex][0]*cartToFrac[0][0] + inputInducedDipole[atomIndex][1]*cartToFrac[0][1] + inputInducedDipole[atomIndex][2]*cartToFrac[0][2],
inputInducedDipole[atomIndex][0]*cartToFrac[1][0] + inputInducedDipole[atomIndex][1]*cartToFrac[1][1] + inputInducedDipole[atomIndex][2]*cartToFrac[1][2], inputInducedDipole[atomIndex][0]*cartToFrac[1][0] + inputInducedDipole[atomIndex][1]*cartToFrac[1][1] + inputInducedDipole[atomIndex][2]*cartToFrac[1][2],
...@@ -5432,7 +5432,7 @@ void AmoebaReferencePmeMultipoleForce::spreadInducedDipolesOnGrid(const vector<R ...@@ -5432,7 +5432,7 @@ void AmoebaReferencePmeMultipoleForce::spreadInducedDipolesOnGrid(const vector<R
int y = (gridPoint[1]+iy) % _pmeGridDimensions[1]; int y = (gridPoint[1]+iy) % _pmeGridDimensions[1];
for (int iz = 0; iz < AMOEBA_PME_ORDER; iz++) { for (int iz = 0; iz < AMOEBA_PME_ORDER; iz++) {
int z = (gridPoint[2]+iz) % _pmeGridDimensions[2]; int z = (gridPoint[2]+iz) % _pmeGridDimensions[2];
RealOpenMM4 t = _thetai[0][atomIndex*AMOEBA_PME_ORDER+ix]; RealOpenMM4 t = _thetai[0][atomIndex*AMOEBA_PME_ORDER+ix];
RealOpenMM4 u = _thetai[1][atomIndex*AMOEBA_PME_ORDER+iy]; RealOpenMM4 u = _thetai[1][atomIndex*AMOEBA_PME_ORDER+iy];
RealOpenMM4 v = _thetai[2][atomIndex*AMOEBA_PME_ORDER+iz]; RealOpenMM4 v = _thetai[2][atomIndex*AMOEBA_PME_ORDER+iz];
...@@ -5654,7 +5654,7 @@ void AmoebaReferencePmeMultipoleForce::computeInducedPotentialFromGrid() ...@@ -5654,7 +5654,7 @@ void AmoebaReferencePmeMultipoleForce::computeInducedPotentialFromGrid()
} }
RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceFixedMultipoleForceAndEnergy(const vector<MultipoleParticleData>& particleData, RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceFixedMultipoleForceAndEnergy(const vector<MultipoleParticleData>& particleData,
vector<RealVec>& forces, vector<RealVec>& torques) const vector<RealVec>& forces, vector<RealVec>& torques) const
{ {
RealOpenMM multipole[10]; RealOpenMM multipole[10];
const int deriv1[] = {1, 4, 7, 8, 10, 15, 17, 13, 14, 19}; const int deriv1[] = {1, 4, 7, 8, 10, 15, 17, 13, 14, 19};
...@@ -5734,7 +5734,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceFixedMultipol ...@@ -5734,7 +5734,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceFixedMultipol
*/ */
RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipoleForceAndEnergy(AmoebaReferenceMultipoleForce::PolarizationType polarizationType, RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipoleForceAndEnergy(AmoebaReferenceMultipoleForce::PolarizationType polarizationType,
const vector<MultipoleParticleData>& particleData, const vector<MultipoleParticleData>& particleData,
vector<RealVec>& forces, vector<RealVec>& torques) const vector<RealVec>& forces, vector<RealVec>& torques) const
{ {
RealOpenMM multipole[10]; RealOpenMM multipole[10];
...@@ -5820,7 +5820,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipole ...@@ -5820,7 +5820,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipole
f[0] += (inducedDipole[k]+inducedDipolePolar[k])*_phi[20*i+j1]; f[0] += (inducedDipole[k]+inducedDipolePolar[k])*_phi[20*i+j1];
f[1] += (inducedDipole[k]+inducedDipolePolar[k])*_phi[20*i+j2]; f[1] += (inducedDipole[k]+inducedDipolePolar[k])*_phi[20*i+j2];
f[2] += (inducedDipole[k]+inducedDipolePolar[k])*_phi[20*i+j3]; f[2] += (inducedDipole[k]+inducedDipolePolar[k])*_phi[20*i+j3];
if (polarizationType == AmoebaReferenceMultipoleForce::Mutual) if (polarizationType == AmoebaReferenceMultipoleForce::Mutual)
{ {
f[0] += (inducedDipole[k]*_phip[10*i+j1] + inducedDipolePolar[k]*_phid[10*i+j1]); f[0] += (inducedDipole[k]*_phip[10*i+j1] + inducedDipolePolar[k]*_phid[10*i+j1]);
...@@ -5899,13 +5899,13 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector ...@@ -5899,13 +5899,13 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields) vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields)
{ {
// Initialize the fields to zero. // Initialize the fields to zero.
RealVec zeroVec(0.0, 0.0, 0.0); RealVec zeroVec(0.0, 0.0, 0.0);
for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++)
std::fill(updateInducedDipoleFields[ii].inducedDipoleField.begin(), updateInducedDipoleFields[ii].inducedDipoleField.end(), zeroVec); std::fill(updateInducedDipoleFields[ii].inducedDipoleField.begin(), updateInducedDipoleFields[ii].inducedDipoleField.end(), zeroVec);
// Add fields from direct space interactions. // Add fields from direct space interactions.
for (unsigned int ii = 0; ii < particleData.size(); ii++) { for (unsigned int ii = 0; ii < particleData.size(); ii++) {
for (unsigned int jj = ii + 1; jj < particleData.size(); jj++) { for (unsigned int jj = ii + 1; jj < particleData.size(); jj++) {
calculateDirectInducedDipolePairIxns(particleData[ii], particleData[jj], updateInducedDipoleFields); calculateDirectInducedDipolePairIxns(particleData[ii], particleData[jj], updateInducedDipoleFields);
...@@ -5916,7 +5916,7 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector ...@@ -5916,7 +5916,7 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
calculateReciprocalSpaceInducedDipoleField(updateInducedDipoleFields); calculateReciprocalSpaceInducedDipoleField(updateInducedDipoleFields);
// self ixn // self ixn
RealOpenMM term = (4.0/3.0)*(_alphaEwald*_alphaEwald*_alphaEwald)/SQRT_PI; RealOpenMM term = (4.0/3.0)*(_alphaEwald*_alphaEwald*_alphaEwald)/SQRT_PI;
for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) { for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) {
...@@ -5924,7 +5924,7 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector ...@@ -5924,7 +5924,7 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
vector<RealVec>& field = updateInducedDipoleFields[ii].inducedDipoleField; vector<RealVec>& field = updateInducedDipoleFields[ii].inducedDipoleField;
for (unsigned int jj = 0; jj < particleData.size(); jj++) { for (unsigned int jj = 0; jj < particleData.size(); jj++) {
field[jj] += inducedDipoles[jj]*term; field[jj] += inducedDipoles[jj]*term;
} }
} }
} }
...@@ -5932,7 +5932,7 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxn(unsig ...@@ -5932,7 +5932,7 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxn(unsig
RealOpenMM preFactor1, RealOpenMM preFactor2, RealOpenMM preFactor1, RealOpenMM preFactor2,
const RealVec& delta, const RealVec& delta,
const vector<RealVec>& inducedDipole, const vector<RealVec>& inducedDipole,
vector<RealVec>& field) const vector<RealVec>& field) const
{ {
// field at i due induced dipole at j // field at i due induced dipole at j
...@@ -5952,7 +5952,7 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons ...@@ -5952,7 +5952,7 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
{ {
// compute the real space portion of the Ewald summation // compute the real space portion of the Ewald summation
RealOpenMM uscale = 1.0; RealOpenMM uscale = 1.0;
RealVec deltaR = particleJ.position - particleI.position; RealVec deltaR = particleJ.position - particleI.position;
...@@ -6013,10 +6013,10 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons ...@@ -6013,10 +6013,10 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
calculateDirectInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, preFactor1, preFactor2, deltaR, calculateDirectInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, preFactor1, preFactor2, deltaR,
*updateInducedDipoleFields[ii].inducedDipoles, *updateInducedDipoleFields[ii].inducedDipoles,
updateInducedDipoleFields[ii].inducedDipoleField); updateInducedDipoleFields[ii].inducedDipoleField);
} }
} }
RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeSelfEnergy(const vector<MultipoleParticleData>& particleData) const RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeSelfEnergy(const vector<MultipoleParticleData>& particleData) const
{ {
RealOpenMM cii = 0.0; RealOpenMM cii = 0.0;
RealOpenMM dii = 0.0; RealOpenMM dii = 0.0;
...@@ -6044,7 +6044,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeSelfEnergy(const vector ...@@ -6044,7 +6044,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeSelfEnergy(const vector
} }
void AmoebaReferencePmeMultipoleForce::calculatePmeSelfTorque(const vector<MultipoleParticleData>& particleData, void AmoebaReferencePmeMultipoleForce::calculatePmeSelfTorque(const vector<MultipoleParticleData>& particleData,
vector<RealVec>& torques) const vector<RealVec>& torques) const
{ {
RealOpenMM term = (2.0/3.0)*(_electric/_dielectric)*(_alphaEwald*_alphaEwald*_alphaEwald)/SQRT_PI; RealOpenMM term = (2.0/3.0)*(_electric/_dielectric)*(_alphaEwald*_alphaEwald*_alphaEwald)/SQRT_PI;
for (unsigned int ii = 0; ii < _numParticles; ii++) { for (unsigned int ii = 0; ii < _numParticles; ii++) {
...@@ -6058,11 +6058,11 @@ void AmoebaReferencePmeMultipoleForce::calculatePmeSelfTorque(const vector<Multi ...@@ -6058,11 +6058,11 @@ void AmoebaReferencePmeMultipoleForce::calculatePmeSelfTorque(const vector<Multi
} }
} }
RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeDirectElectrostaticPairIxn(const MultipoleParticleData& particleI, RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeDirectElectrostaticPairIxn(const MultipoleParticleData& particleI,
const MultipoleParticleData& particleJ, const MultipoleParticleData& particleJ,
const vector<RealOpenMM>& scalingFactors, const vector<RealOpenMM>& scalingFactors,
vector<RealVec>& forces, vector<RealVec>& forces,
vector<RealVec>& torques) const vector<RealVec>& torques) const
{ {
unsigned int iIndex = particleI.particleIndex; unsigned int iIndex = particleI.particleIndex;
...@@ -6483,7 +6483,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector ...@@ -6483,7 +6483,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector
vector<RealOpenMM> scaleFactors(LAST_SCALE_TYPE_INDEX); vector<RealOpenMM> scaleFactors(LAST_SCALE_TYPE_INDEX);
for (unsigned int kk = 0; kk < scaleFactors.size(); kk++) { for (unsigned int kk = 0; kk < scaleFactors.size(); kk++) {
scaleFactors[kk] = 1.0; scaleFactors[kk] = 1.0;
} }
// loop over particle pairs for direct space interactions // loop over particle pairs for direct space interactions
...@@ -6508,6 +6508,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector ...@@ -6508,6 +6508,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector
energy += computeReciprocalSpaceInducedDipoleForceAndEnergy(getPolarizationType(), particleData, forces, torques); energy += computeReciprocalSpaceInducedDipoleForceAndEnergy(getPolarizationType(), particleData, forces, torques);
energy += computeReciprocalSpaceFixedMultipoleForceAndEnergy(particleData, forces, torques); energy += computeReciprocalSpaceFixedMultipoleForceAndEnergy(particleData, forces, torques);
energy += calculatePmeSelfEnergy(particleData); energy += calculatePmeSelfEnergy(particleData);
return energy; return energy;
} }
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