Skip to content

GitLab

Commit 308534a7 authored by Peter Eastman's avatar Peter Eastman
Browse files

Fixed bug with ZBisect axis type in AmoebaMultipoleForce

parent b465c982
Hide whitespace changes
Inline Side-by-side
Showing with 194 additions and 3 deletions
plugins/amoeba/platforms/cuda/src/kernels/multipoles.cu
View file @ 308534a7
......@@ -141,7 +141,7 @@ extern "C" __global__ void computeLabFrameMoments(real4* __restrict__ posq, int4
// Check the chirality and see whether it needs to be reversed
bool reverse = false;
if (axisType != 0 && particles.x >= 0 && particles.y >=0 && particles.z >= 0) {
if (axisType == 0 && particles.x >= 0 && particles.y >=0 && particles.z >= 0) {
real4 posY = posq[particles.y];
real delta[4][3];
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp
View file @ 308534a7
......@@ -3131,6 +3131,98 @@ void testTriclinic() {
ASSERT_EQUAL_TOL(expectedPotential[i], potential[i], 1e-4);
}
void testZBisect() {
System system;
for (int i = 0; i < 7; i++)
system.addParticle(1.0);
system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
AmoebaMultipoleForce* force = new AmoebaMultipoleForce();
system.addForce(force);
force->setNonbondedMethod(AmoebaMultipoleForce::PME);
force->setCutoffDistance(1.2);
double charge[] = {-1.01875, 0, 0, 0, -0.51966, 0.25983, 0.25983};
double dipole[7][3] = {
{0.06620218576365969, 0.056934176095985306, 0.06298584667720743},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0.007556121391156931},
{-0.05495981592297553, 0, -0.0030787530116780605},
{-0.05495981592297553, 0, -0.0030787530116780605}};
double quadrupole[7][9] = {
{-0.0004042865090302655, 0.0010450291005955025, -0.0010871640586155112, 0.0010450291005955025, 0.0002512789255424535, -0.0009504541350087216, -0.0010871640586155112, -0.0009504541350087216, 0.00015300758348781198},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.00035403072392177476, 0.009334284009749387, 0, 0.009334284009749387, -0.00039025708016361216, 0, 0, 0, 3.622635624183737e-05},
{-3.42848251678095e-05, 0.07374084367702016, -1.894859653979126e-06, 0.07374084367702016, -0.00010024087598069868, 0, -1.894859653979126e-06, 0, 0.00013452570114850818},
{-3.42848251678095e-05, 0.07374084367702016, -1.894859653979126e-06, 0.07374084367702016, -0.00010024087598069868, 0, -1.894859653979126e-06, 0, 0.00013452570114850818}
};
int axis[7][4] = {
{2, 2, 1, 3},
{5, -1, -1, -1},
{5, -1, -1, -1},
{5, -1, -1, -1},
{1, 5, 6, -1},
{0, 4, 6, -1},
{0, 4, 5, -1}
};
double thole = 0.39;
double damping[] = {0.33178695365189015, 0.33178695365189015, 0.33178695365189015, 0.33178695365189015, 0.306987653777382, 0.2813500172269554, 0.2813500172269554};
double polarity[] = {0.001334, 0.001334, 0.001334, 0.001334, 0.000837, 0.000496, 0.000496};
for (int i = 0; i < 7; i++) {
vector<double> d, q;
for (int j = 0; j < 3; j++)
d.push_back(dipole[i][j]);
for (int j = 0; j < 9; j++)
q.push_back(quadrupole[i][j]);
force->addMultipole(charge[i], d, q, axis[i][0], axis[i][1], axis[i][2], axis[i][3], thole, damping[i], polarity[i]);
}
for (int i = 0; i < 4; i++) {
vector<int> map;
if (i != 0) map.push_back(0);
force->setCovalentMap(i, AmoebaMultipoleForce::Covalent12, map);
map.clear();
if (i != 1) map.push_back(1);
if (i != 2) map.push_back(2);
if (i != 3) map.push_back(3);
force->setCovalentMap(i, AmoebaMultipoleForce::Covalent13, map);
map.clear();
map.push_back(0);
map.push_back(1);
map.push_back(2);
map.push_back(3);
force->setCovalentMap(i, AmoebaMultipoleForce::PolarizationCovalent11, map);
}
for (int i = 4; i < 7; i++) {
vector<int> map;
if (i != 4) map.push_back(4);
force->setCovalentMap(i, AmoebaMultipoleForce::Covalent12, map);
map.clear();
if (i != 5) map.push_back(5);
if (i != 6) map.push_back(6);
force->setCovalentMap(i, AmoebaMultipoleForce::Covalent13, map);
map.clear();
map.push_back(4);
map.push_back(5);
map.push_back(6);
force->setCovalentMap(i, AmoebaMultipoleForce::PolarizationCovalent11, map);
}
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
vector<Vec3> positions;
positions.push_back(Vec3(-0.06317711175870899, -0.04905009196658128, 0.0767217));
positions.push_back(Vec3(-0.049166918626451395, -0.20747614470348363, 0.03979849999999996));
positions.push_back(Vec3(-0.19317150000000005, -0.05811762921948427, 0.1632788999999999));
positions.push_back(Vec3(0.04465103038516016, -0.018345116763806235, 0.18531239999999993));
positions.push_back(Vec3(0.005630299999999998, 0.40965770000000035, 0.5731495));
positions.push_back(Vec3(0.036148100000000016, 0.3627041999999996, 0.49299430000000033));
positions.push_back(Vec3(0.07781149999999992, 0.4178183000000004, 0.6355703000000004));
context.setPositions(positions);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(-84.1532, state.getPotentialEnergy(), 0.01);
}
int main(int argc, char* argv[]) {
try {
std::cout << "TestCudaAmoebaMultipoleForce running test..." << std::endl;
......@@ -3183,6 +3275,10 @@ int main(int argc, char* argv[]) {
// triclinic box of water
testTriclinic();
// test the ZBisect axis type.
testZBisect();
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
View file @ 308534a7
......@@ -359,7 +359,7 @@ void AmoebaReferenceMultipoleForce::checkChiralCenterAtParticle(MultipoleParticl
MultipoleParticleData& particleY) const
{
if (axisType == AmoebaMultipoleForce::ZThenX || particleY.particleIndex == -1) {
if (axisType != AmoebaMultipoleForce::ZThenX || particleY.particleIndex == -1) {
return;
}
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaMultipoleForce.cpp
View file @ 308534a7
......@@ -2938,6 +2938,98 @@ void testTriclinic() {
ASSERT_EQUAL_TOL(expectedEnergy, state.getPotentialEnergy(), 1e-2);
}
void testZBisect() {
System system;
for (int i = 0; i < 7; i++)
system.addParticle(1.0);
system.setDefaultPeriodicBoxVectors(Vec3(4, 0, 0), Vec3(0, 4, 0), Vec3(0, 0, 4));
AmoebaMultipoleForce* force = new AmoebaMultipoleForce();
system.addForce(force);
force->setNonbondedMethod(AmoebaMultipoleForce::PME);
force->setCutoffDistance(1.2);
double charge[] = {-1.01875, 0, 0, 0, -0.51966, 0.25983, 0.25983};
double dipole[7][3] = {
{0.06620218576365969, 0.056934176095985306, 0.06298584667720743},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0},
{0, 0, 0.007556121391156931},
{-0.05495981592297553, 0, -0.0030787530116780605},
{-0.05495981592297553, 0, -0.0030787530116780605}};
double quadrupole[7][9] = {
{-0.0004042865090302655, 0.0010450291005955025, -0.0010871640586155112, 0.0010450291005955025, 0.0002512789255424535, -0.0009504541350087216, -0.0010871640586155112, -0.0009504541350087216, 0.00015300758348781198},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0},
{0.00035403072392177476, 0.009334284009749387, 0, 0.009334284009749387, -0.00039025708016361216, 0, 0, 0, 3.622635624183737e-05},
{-3.42848251678095e-05, 0.07374084367702016, -1.894859653979126e-06, 0.07374084367702016, -0.00010024087598069868, 0, -1.894859653979126e-06, 0, 0.00013452570114850818},
{-3.42848251678095e-05, 0.07374084367702016, -1.894859653979126e-06, 0.07374084367702016, -0.00010024087598069868, 0, -1.894859653979126e-06, 0, 0.00013452570114850818}
};
int axis[7][4] = {
{2, 2, 1, 3},
{5, -1, -1, -1},
{5, -1, -1, -1},
{5, -1, -1, -1},
{1, 5, 6, -1},
{0, 4, 6, -1},
{0, 4, 5, -1}
};
double thole = 0.39;
double damping[] = {0.33178695365189015, 0.33178695365189015, 0.33178695365189015, 0.33178695365189015, 0.306987653777382, 0.2813500172269554, 0.2813500172269554};
double polarity[] = {0.001334, 0.001334, 0.001334, 0.001334, 0.000837, 0.000496, 0.000496};
for (int i = 0; i < 7; i++) {
vector<double> d, q;
for (int j = 0; j < 3; j++)
d.push_back(dipole[i][j]);
for (int j = 0; j < 9; j++)
q.push_back(quadrupole[i][j]);
force->addMultipole(charge[i], d, q, axis[i][0], axis[i][1], axis[i][2], axis[i][3], thole, damping[i], polarity[i]);
}
for (int i = 0; i < 4; i++) {
vector<int> map;
if (i != 0) map.push_back(0);
force->setCovalentMap(i, AmoebaMultipoleForce::Covalent12, map);
map.clear();
if (i != 1) map.push_back(1);
if (i != 2) map.push_back(2);
if (i != 3) map.push_back(3);
force->setCovalentMap(i, AmoebaMultipoleForce::Covalent13, map);
map.clear();
map.push_back(0);
map.push_back(1);
map.push_back(2);
map.push_back(3);
force->setCovalentMap(i, AmoebaMultipoleForce::PolarizationCovalent11, map);
}
for (int i = 4; i < 7; i++) {
vector<int> map;
if (i != 4) map.push_back(4);
force->setCovalentMap(i, AmoebaMultipoleForce::Covalent12, map);
map.clear();
if (i != 5) map.push_back(5);
if (i != 6) map.push_back(6);
force->setCovalentMap(i, AmoebaMultipoleForce::Covalent13, map);
map.clear();
map.push_back(4);
map.push_back(5);
map.push_back(6);
force->setCovalentMap(i, AmoebaMultipoleForce::PolarizationCovalent11, map);
}
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
vector<Vec3> positions;
positions.push_back(Vec3(-0.06317711175870899, -0.04905009196658128, 0.0767217));
positions.push_back(Vec3(-0.049166918626451395, -0.20747614470348363, 0.03979849999999996));
positions.push_back(Vec3(-0.19317150000000005, -0.05811762921948427, 0.1632788999999999));
positions.push_back(Vec3(0.04465103038516016, -0.018345116763806235, 0.18531239999999993));
positions.push_back(Vec3(0.005630299999999998, 0.40965770000000035, 0.5731495));
positions.push_back(Vec3(0.036148100000000016, 0.3627041999999996, 0.49299430000000033));
positions.push_back(Vec3(0.07781149999999992, 0.4178183000000004, 0.6355703000000004));
context.setPositions(positions);
State state = context.getState(State::Energy);
ASSERT_EQUAL_TOL(-84.1532, state.getPotentialEnergy(), 0.01);
}
int main(int numberOfArguments, char* argv[]) {
try {
......@@ -2987,7 +3079,10 @@ int main(int numberOfArguments, char* argv[]) {
// triclinic box of water
testTriclinic();
// test the ZBisect axis type.
testZBisect();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
......
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