Skip to content

GitLab

"openmmapi/src/NoseHooverIntegrator.cpp" did not exist on "925b00ec7ee5ba38e3a72f1da8a88df5a3d51bd5"
Commit 0c96184b authored by Lee-Ping Wang's avatar Lee-Ping Wang
Browse files

Revert to revision 3571.

parent f028c07d
Hide whitespace changes
Inline Side-by-side
Showing with 61 additions and 81 deletions
plugins/amoeba/openmmapi/include/openmm/AmoebaMultipoleForce.h
View file @ 0c96184b
......@@ -319,21 +319,16 @@ public:
* <i>lowest nonvanishing moment</i> has a well defined value. This means that if the system has a net
* nonzero charge, the dipole and quadrupole moments are not well defined and should be ignored. If the
* net charge is zero, the dipole moment is well defined (and really represents a dipole density), but
* the quadrupole moment is still undefined and should be ignored. Specify "compute = True" if you want
* to do an energy and force calculation prior to getting the moments; if compute = False, it will be faster
* but we will lose accuracy. Specify "order" for the maximum order computed (order = 2 means we go up
* through quadrupoles.)
* the quadrupole moment is still undefined and should be ignored.
*
* @param context context
* @param compute Whether to compute energies and forces first
* @param order Maximum order to go through, above 2 for quadrupole has no effect
* @param outputMultipoleMonents (charge,
dipole_x, dipole_y, dipole_z,
quadrupole_xx, quadrupole_xy, quadrupole_xz,
quadrupole_yx, quadrupole_yy, quadrupole_yz,
quadrupole_zx, quadrupole_zy, quadrupole_zz)
*/
void getSystemMultipoleMoments(Context& context, bool compute, int order, std::vector< double >& outputMultipoleMoments);
void getSystemMultipoleMoments(Context& context, std::vector< double >& outputMultipoleMoments);
/**
* Update the multipole parameters in a Context to match those stored in this Force object. This method
* provides an efficient method to update certain parameters in an existing Context without needing to reinitialize it.
......
plugins/amoeba/openmmapi/include/openmm/amoebaKernels.h
View file @ 0c96184b
......@@ -351,7 +351,7 @@ public:
virtual void getElectrostaticPotential( ContextImpl& context, const std::vector< Vec3 >& inputGrid,
std::vector< double >& outputElectrostaticPotential ) = 0;
virtual void getSystemMultipoleMoments( ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMonents ) = 0;
virtual void getSystemMultipoleMoments( ContextImpl& context, std::vector< double >& outputMultipoleMonents ) = 0;
/**
* Copy changed parameters over to a context.
*
......
plugins/amoeba/openmmapi/include/openmm/internal/AmoebaMultipoleForceImpl.h
View file @ 0c96184b
......@@ -85,7 +85,7 @@ public:
void getElectrostaticPotential( ContextImpl& context, const std::vector< Vec3 >& inputGrid,
std::vector< double >& outputElectrostaticPotential );
void getSystemMultipoleMoments( ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMonents );
void getSystemMultipoleMoments( ContextImpl& context, std::vector< double >& outputMultipoleMonents );
void updateParametersInContext(ContextImpl& context);
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp
View file @ 0c96184b
......@@ -230,8 +230,8 @@ void AmoebaMultipoleForce::getElectrostaticPotential( const std::vector< Vec3 >&
dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getElectrostaticPotential(getContextImpl(context), inputGrid, outputElectrostaticPotential);
}
void AmoebaMultipoleForce::getSystemMultipoleMoments(Context& context, bool compute, int order, std::vector< double >& outputMultipoleMonents ){
dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getSystemMultipoleMoments(getContextImpl(context), compute, order, outputMultipoleMonents);
void AmoebaMultipoleForce::getSystemMultipoleMoments(Context& context, std::vector< double >& outputMultipoleMonents ){
dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getSystemMultipoleMoments(getContextImpl(context), outputMultipoleMonents);
}
ForceImpl* AmoebaMultipoleForce::createImpl() {
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForceImpl.cpp
View file @ 0c96184b
......@@ -188,8 +188,8 @@ void AmoebaMultipoleForceImpl::getElectrostaticPotential( ContextImpl& context,
kernel.getAs<CalcAmoebaMultipoleForceKernel>().getElectrostaticPotential(context, inputGrid, outputElectrostaticPotential);
}
void AmoebaMultipoleForceImpl::getSystemMultipoleMoments( ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMonents ){
kernel.getAs<CalcAmoebaMultipoleForceKernel>().getSystemMultipoleMoments(context, compute, order, outputMultipoleMonents);
void AmoebaMultipoleForceImpl::getSystemMultipoleMoments( ContextImpl& context, std::vector< double >& outputMultipoleMonents ){
kernel.getAs<CalcAmoebaMultipoleForceKernel>().getSystemMultipoleMoments(context, outputMultipoleMonents);
}
void AmoebaMultipoleForceImpl::updateParametersInContext(ContextImpl& context) {
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 0c96184b
......@@ -1596,7 +1596,7 @@ void CudaCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl&
}
template <class T, class T4, class M4>
void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextImpl& context, int order, vector<double>& outputMultipoleMoments) {
void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextImpl& context, vector<double>& outputMultipoleMoments) {
// Compute the local coordinates relative to the center of mass.
int numAtoms = cu.getNumAtoms();
vector<T4> posq;
......@@ -1641,66 +1641,56 @@ void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextIm
double zyqdp = 0.0;
double zzqdp = 0.0;
vector<T> labDipoleVec, inducedDipoleVec, quadrupoleVec;
if (order >= 1) {
labFrameDipoles->download(labDipoleVec);
inducedDipole->download(inducedDipoleVec);
}
if (order >= 2) {
labFrameQuadrupoles->download(quadrupoleVec);
}
labFrameDipoles->download(labDipoleVec);
inducedDipole->download(inducedDipoleVec);
labFrameQuadrupoles->download(quadrupoleVec);
for (int i = 0; i < numAtoms; i++) {
totalCharge += posqLocal[i].w;
double netDipoleX = (labDipoleVec[3*i] + inducedDipoleVec[3*i]);
double netDipoleY = (labDipoleVec[3*i+1] + inducedDipoleVec[3*i+1]);
double netDipoleZ = (labDipoleVec[3*i+2] + inducedDipoleVec[3*i+2]);
if (order >= 1) {
xdpl += posqLocal[i].x*posqLocal[i].w + netDipoleX;
ydpl += posqLocal[i].y*posqLocal[i].w + netDipoleY;
zdpl += posqLocal[i].z*posqLocal[i].w + netDipoleZ;
}
if (order >= 2) {
xxqdp += posqLocal[i].x*posqLocal[i].x*posqLocal[i].w + 2*posqLocal[i].x*netDipoleX;
xyqdp += posqLocal[i].x*posqLocal[i].y*posqLocal[i].w + posqLocal[i].x*netDipoleY + posqLocal[i].y*netDipoleX;
xzqdp += posqLocal[i].x*posqLocal[i].z*posqLocal[i].w + posqLocal[i].x*netDipoleZ + posqLocal[i].z*netDipoleX;
yxqdp += posqLocal[i].y*posqLocal[i].x*posqLocal[i].w + posqLocal[i].y*netDipoleX + posqLocal[i].x*netDipoleY;
yyqdp += posqLocal[i].y*posqLocal[i].y*posqLocal[i].w + 2*posqLocal[i].y*netDipoleY;
yzqdp += posqLocal[i].y*posqLocal[i].z*posqLocal[i].w + posqLocal[i].y*netDipoleZ + posqLocal[i].z*netDipoleY;
zxqdp += posqLocal[i].z*posqLocal[i].x*posqLocal[i].w + posqLocal[i].z*netDipoleX + posqLocal[i].x*netDipoleZ;
zyqdp += posqLocal[i].z*posqLocal[i].y*posqLocal[i].w + posqLocal[i].z*netDipoleY + posqLocal[i].y*netDipoleZ;
zzqdp += posqLocal[i].z*posqLocal[i].z*posqLocal[i].w + 2*posqLocal[i].z*netDipoleZ;
}
double netDipoleX = (labDipoleVec[3*i] + inducedDipoleVec[3*i]);
double netDipoleY = (labDipoleVec[3*i+1] + inducedDipoleVec[3*i+1]);
double netDipoleZ = (labDipoleVec[3*i+2] + inducedDipoleVec[3*i+2]);
xdpl += posqLocal[i].x*posqLocal[i].w + netDipoleX;
ydpl += posqLocal[i].y*posqLocal[i].w + netDipoleY;
zdpl += posqLocal[i].z*posqLocal[i].w + netDipoleZ;
xxqdp += posqLocal[i].x*posqLocal[i].x*posqLocal[i].w + 2*posqLocal[i].x*netDipoleX;
xyqdp += posqLocal[i].x*posqLocal[i].y*posqLocal[i].w + posqLocal[i].x*netDipoleY + posqLocal[i].y*netDipoleX;
xzqdp += posqLocal[i].x*posqLocal[i].z*posqLocal[i].w + posqLocal[i].x*netDipoleZ + posqLocal[i].z*netDipoleX;
yxqdp += posqLocal[i].y*posqLocal[i].x*posqLocal[i].w + posqLocal[i].y*netDipoleX + posqLocal[i].x*netDipoleY;
yyqdp += posqLocal[i].y*posqLocal[i].y*posqLocal[i].w + 2*posqLocal[i].y*netDipoleY;
yzqdp += posqLocal[i].y*posqLocal[i].z*posqLocal[i].w + posqLocal[i].y*netDipoleZ + posqLocal[i].z*netDipoleY;
zxqdp += posqLocal[i].z*posqLocal[i].x*posqLocal[i].w + posqLocal[i].z*netDipoleX + posqLocal[i].x*netDipoleZ;
zyqdp += posqLocal[i].z*posqLocal[i].y*posqLocal[i].w + posqLocal[i].z*netDipoleY + posqLocal[i].y*netDipoleZ;
zzqdp += posqLocal[i].z*posqLocal[i].z*posqLocal[i].w + 2*posqLocal[i].z*netDipoleZ;
}
// Convert the quadrupole from traced to traceless form.
if (order >= 2) {
double qave = (xxqdp + yyqdp + zzqdp)/3;
xxqdp = 1.5*(xxqdp-qave);
xyqdp = 1.5*xyqdp;
xzqdp = 1.5*xzqdp;
yxqdp = 1.5*yxqdp;
yyqdp = 1.5*(yyqdp-qave);
yzqdp = 1.5*yzqdp;
zxqdp = 1.5*zxqdp;
zyqdp = 1.5*zyqdp;
zzqdp = 1.5*(zzqdp-qave);
double qave = (xxqdp + yyqdp + zzqdp)/3;
xxqdp = 1.5*(xxqdp-qave);
xyqdp = 1.5*xyqdp;
xzqdp = 1.5*xzqdp;
yxqdp = 1.5*yxqdp;
yyqdp = 1.5*(yyqdp-qave);
yzqdp = 1.5*yzqdp;
zxqdp = 1.5*zxqdp;
zyqdp = 1.5*zyqdp;
zzqdp = 1.5*(zzqdp-qave);
// Add the traceless atomic quadrupoles to the total quadrupole moment.
for (int i = 0; i < numAtoms; i++) {
xxqdp = xxqdp + 3*quadrupoleVec[5*i];
xyqdp = xyqdp + 3*quadrupoleVec[5*i+1];
xzqdp = xzqdp + 3*quadrupoleVec[5*i+2];
yxqdp = yxqdp + 3*quadrupoleVec[5*i+1];
yyqdp = yyqdp + 3*quadrupoleVec[5*i+3];
yzqdp = yzqdp + 3*quadrupoleVec[5*i+4];
zxqdp = zxqdp + 3*quadrupoleVec[5*i+2];
zyqdp = zyqdp + 3*quadrupoleVec[5*i+4];
zzqdp = zzqdp + -3*(quadrupoleVec[5*i]+quadrupoleVec[5*i+3]);
}
}
double debye = 48.0321; // If we use 4.80321 then the result is ten times too small.
for (int i = 0; i < numAtoms; i++) {
xxqdp = xxqdp + 3*quadrupoleVec[5*i];
xyqdp = xyqdp + 3*quadrupoleVec[5*i+1];
xzqdp = xzqdp + 3*quadrupoleVec[5*i+2];
yxqdp = yxqdp + 3*quadrupoleVec[5*i+1];
yyqdp = yyqdp + 3*quadrupoleVec[5*i+3];
yzqdp = yzqdp + 3*quadrupoleVec[5*i+4];
zxqdp = zxqdp + 3*quadrupoleVec[5*i+2];
zyqdp = zyqdp + 3*quadrupoleVec[5*i+4];
zzqdp = zzqdp + -3*(quadrupoleVec[5*i]+quadrupoleVec[5*i+3]);
}
double debye = 4.80321;
outputMultipoleMoments.resize(13);
outputMultipoleMoments[0] = totalCharge;
outputMultipoleMoments[1] = xdpl*debye;
......@@ -1717,15 +1707,14 @@ void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextIm
outputMultipoleMoments[12] = zzqdp*debye;
}
void CudaCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, bool compute, int order, vector<double>& outputMultipoleMoments) {
if (compute)
context.calcForcesAndEnergy(false, false, -1);
void CudaCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, vector<double>& outputMultipoleMoments) {
context.calcForcesAndEnergy(false, false, -1);
if (cu.getUseDoublePrecision())
computeSystemMultipoleMoments<double, double4, double4>(context, order, outputMultipoleMoments);
computeSystemMultipoleMoments<double, double4, double4>(context, outputMultipoleMoments);
else if (cu.getUseMixedPrecision())
computeSystemMultipoleMoments<float, float4, double4>(context, order, outputMultipoleMoments);
computeSystemMultipoleMoments<float, float4, double4>(context, outputMultipoleMoments);
else
computeSystemMultipoleMoments<float, float4, float4>(context, order, outputMultipoleMoments);
computeSystemMultipoleMoments<float, float4, float4>(context, outputMultipoleMoments);
}
void CudaCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaMultipoleForce& force) {
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h
View file @ 0c96184b
......@@ -341,15 +341,13 @@ public:
* Get the system multipole moments
*
* @param context context
* @param compute Whether to compute energies and forces first
* @param order Maximum order to go through, above 2 for quadrupole has no effect
* @param outputMultipoleMoments (charge,
* dipole_x, dipole_y, dipole_z,
* quadrupole_xx, quadrupole_xy, quadrupole_xz,
* quadrupole_yx, quadrupole_yy, quadrupole_yz,
* quadrupole_zx, quadrupole_zy, quadrupole_zz )
*/
void getSystemMultipoleMoments(ContextImpl& context, bool compute, int order, std::vector<double>& outputMultipoleMoments);
void getSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
/**
* Copy changed parameters over to a context.
*
......@@ -370,7 +368,7 @@ private:
const char* getSortKey() const {return "value.y";}
};
void initializeScaleFactors();
template <class T, class T4, class M4> void computeSystemMultipoleMoments(ContextImpl& context, int order, std::vector<double>& outputMultipoleMoments);
template <class T, class T4, class M4> void computeSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
int numMultipoles, maxInducedIterations;
int fixedFieldThreads, inducedFieldThreads, electrostaticsThreads;
double inducedEpsilon;
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp
View file @ 0c96184b
......@@ -2009,7 +2009,7 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::No
context.setPositions(positions);
if( testName == "testSystemMultipoleMoments" ){
amoebaMultipoleForce->getSystemMultipoleMoments(context, 1, 2, outputMultipoleMoments );
amoebaMultipoleForce->getSystemMultipoleMoments(context, outputMultipoleMoments );
} else if( testName == "testMultipoleGridPotential" ){
amoebaMultipoleForce->getElectrostaticPotential( inputGrid, context, outputGridPotential );
} else {
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 0c96184b
......@@ -708,7 +708,7 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextI
return;
}
void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMoments){
void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments){
// retrieve masses
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 0c96184b
......@@ -380,8 +380,6 @@ public:
* Get the system multipole moments.
*
* @param context context
* @param compute Whether to compute energies and forces first
* @param order Maximum order to go through, above 2 for quadrupole has no effect
* @param outputMultipoleMonents vector of multipole moments:
(charge,
dipole_x, dipole_y, dipole_z,
......@@ -389,7 +387,7 @@ public:
quadrupole_yx, quadrupole_yy, quadrupole_yz,
quadrupole_zx, quadrupole_zy, quadrupole_zz )
*/
void getSystemMultipoleMoments(ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMoments);
void getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments);
/**
* Copy changed parameters over to a context.
*
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaMultipoleForce.cpp
View file @ 0c96184b
......@@ -1915,7 +1915,7 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::No
context.setPositions(positions);
if( testName == "testSystemMultipoleMoments" ){
amoebaMultipoleForce->getSystemMultipoleMoments( context, 1, 2, outputMultipoleMoments );
amoebaMultipoleForce->getSystemMultipoleMoments( context, outputMultipoleMoments );
} else if( testName == "testMultipoleGridPotential" ){
amoebaMultipoleForce->getElectrostaticPotential( inputGrid, context, outputGridPotential );
} else {
......
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