Commit dd1a959b authored by Peter Eastman's avatar Peter Eastman
Browse files

Long range dispersion correction is enabled by default

parent 5228a07e
...@@ -46,7 +46,7 @@ using std::string; ...@@ -46,7 +46,7 @@ using std::string;
using std::stringstream; using std::stringstream;
using std::vector; using std::vector;
NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), rfDielectric(78.3), ewaldErrorTol(5e-4), useDispersionCorrection(false) { NonbondedForce::NonbondedForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), rfDielectric(78.3), ewaldErrorTol(5e-4), useDispersionCorrection(true) {
} }
NonbondedForce::NonbondedMethod NonbondedForce::getNonbondedMethod() const { NonbondedForce::NonbondedMethod NonbondedForce::getNonbondedMethod() const {
......
...@@ -95,10 +95,6 @@ void NonbondedForceImpl::initialize(ContextImpl& context) { ...@@ -95,10 +95,6 @@ void NonbondedForceImpl::initialize(ContextImpl& context) {
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("NonbondedForce: The cutoff distance cannot be greater than half the periodic box size."); throw OpenMMException("NonbondedForce: The cutoff distance cannot be greater than half the periodic box size.");
} }
else {
if (owner.getUseDispersionCorrection())
throw OpenMMException("NonbondedForce: Dispersion correction may only be used with periodic boundary conditions.");
}
dynamic_cast<CalcNonbondedForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner); dynamic_cast<CalcNonbondedForceKernel&>(kernel.getImpl()).initialize(context.getSystem(), owner);
} }
...@@ -176,6 +172,9 @@ int NonbondedForceImpl::findZero(const NonbondedForceImpl::ErrorFunction& f, int ...@@ -176,6 +172,9 @@ int NonbondedForceImpl::findZero(const NonbondedForceImpl::ErrorFunction& f, int
} }
double NonbondedForceImpl::calcDispersionCorrection(const System& system, const NonbondedForce& force) { double NonbondedForceImpl::calcDispersionCorrection(const System& system, const NonbondedForce& force) {
if (force.getNonbondedMethod() == NonbondedForce::NoCutoff || force.getNonbondedMethod() == NonbondedForce::CutoffNonPeriodic)
return 0.0;
// Identify all particle classes (defined by sigma and epsilon), and count the number of // Identify all particle classes (defined by sigma and epsilon), and count the number of
// particles in each class. // particles in each class.
......
...@@ -631,14 +631,14 @@ void testDispersionCorrection() { ...@@ -631,14 +631,14 @@ void testDispersionCorrection() {
Context context(system, integrator, platform); Context context(system, integrator, platform);
context.setPositions(positions); context.setPositions(positions);
double energy1 = context.getState(State::Energy).getPotentialEnergy(); double energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(true); nonbonded->setUseDispersionCorrection(false);
context.reinitialize(); context.reinitialize();
context.setPositions(positions); context.setPositions(positions);
double energy2 = context.getState(State::Energy).getPotentialEnergy(); double energy2 = context.getState(State::Energy).getPotentialEnergy();
double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9; double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3; double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize); double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4); ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
// Now modify half the particles to be different, and see if it is still correct. // Now modify half the particles to be different, and see if it is still correct.
...@@ -648,11 +648,11 @@ void testDispersionCorrection() { ...@@ -648,11 +648,11 @@ void testDispersionCorrection() {
numType2++; numType2++;
} }
int numType1 = numParticles-numType2; int numType1 = numParticles-numType2;
nonbonded->setUseDispersionCorrection(false); nonbonded->setUseDispersionCorrection(true);
context.reinitialize(); context.reinitialize();
context.setPositions(positions); context.setPositions(positions);
energy1 = context.getState(State::Energy).getPotentialEnergy(); energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(true); nonbonded->setUseDispersionCorrection(false);
context.reinitialize(); context.reinitialize();
context.setPositions(positions); context.setPositions(positions);
energy2 = context.getState(State::Energy).getPotentialEnergy(); energy2 = context.getState(State::Energy).getPotentialEnergy();
...@@ -667,7 +667,7 @@ void testDispersionCorrection() { ...@@ -667,7 +667,7 @@ void testDispersionCorrection() {
term1 /= (numParticles*(numParticles+1))/2; term1 /= (numParticles*(numParticles+1))/2;
term2 /= (numParticles*(numParticles+1))/2; term2 /= (numParticles*(numParticles+1))/2;
expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize); expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4); ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
} }
int main() { int main() {
......
...@@ -636,14 +636,14 @@ void testDispersionCorrection() { ...@@ -636,14 +636,14 @@ void testDispersionCorrection() {
Context context(system, integrator, platform); Context context(system, integrator, platform);
context.setPositions(positions); context.setPositions(positions);
double energy1 = context.getState(State::Energy).getPotentialEnergy(); double energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(true); nonbonded->setUseDispersionCorrection(false);
context.reinitialize(); context.reinitialize();
context.setPositions(positions); context.setPositions(positions);
double energy2 = context.getState(State::Energy).getPotentialEnergy(); double energy2 = context.getState(State::Energy).getPotentialEnergy();
double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9; double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3; double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize); double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4); ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
// Now modify half the particles to be different, and see if it is still correct. // Now modify half the particles to be different, and see if it is still correct.
...@@ -653,11 +653,11 @@ void testDispersionCorrection() { ...@@ -653,11 +653,11 @@ void testDispersionCorrection() {
numType2++; numType2++;
} }
int numType1 = numParticles-numType2; int numType1 = numParticles-numType2;
nonbonded->setUseDispersionCorrection(false); nonbonded->setUseDispersionCorrection(true);
context.reinitialize(); context.reinitialize();
context.setPositions(positions); context.setPositions(positions);
energy1 = context.getState(State::Energy).getPotentialEnergy(); energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(true); nonbonded->setUseDispersionCorrection(false);
context.reinitialize(); context.reinitialize();
context.setPositions(positions); context.setPositions(positions);
energy2 = context.getState(State::Energy).getPotentialEnergy(); energy2 = context.getState(State::Energy).getPotentialEnergy();
...@@ -672,7 +672,7 @@ void testDispersionCorrection() { ...@@ -672,7 +672,7 @@ void testDispersionCorrection() {
term1 /= (numParticles*(numParticles+1))/2; term1 /= (numParticles*(numParticles+1))/2;
term2 /= (numParticles*(numParticles+1))/2; term2 /= (numParticles*(numParticles+1))/2;
expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize); expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4); ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
} }
int main() { int main() {
......
...@@ -384,14 +384,14 @@ void testDispersionCorrection() { ...@@ -384,14 +384,14 @@ void testDispersionCorrection() {
Context context(system, integrator, platform); Context context(system, integrator, platform);
context.setPositions(positions); context.setPositions(positions);
double energy1 = context.getState(State::Energy).getPotentialEnergy(); double energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(true); nonbonded->setUseDispersionCorrection(false);
context.reinitialize(); context.reinitialize();
context.setPositions(positions); context.setPositions(positions);
double energy2 = context.getState(State::Energy).getPotentialEnergy(); double energy2 = context.getState(State::Energy).getPotentialEnergy();
double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9; double term1 = (0.5*pow(1.1, 12)/pow(cutoff, 9))/9;
double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3; double term2 = (0.5*pow(1.1, 6)/pow(cutoff, 3))/3;
double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize); double expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4); ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
// Now modify half the particles to be different, and see if it is still correct. // Now modify half the particles to be different, and see if it is still correct.
...@@ -401,11 +401,11 @@ void testDispersionCorrection() { ...@@ -401,11 +401,11 @@ void testDispersionCorrection() {
numType2++; numType2++;
} }
int numType1 = numParticles-numType2; int numType1 = numParticles-numType2;
nonbonded->setUseDispersionCorrection(false); nonbonded->setUseDispersionCorrection(true);
context.reinitialize(); context.reinitialize();
context.setPositions(positions); context.setPositions(positions);
energy1 = context.getState(State::Energy).getPotentialEnergy(); energy1 = context.getState(State::Energy).getPotentialEnergy();
nonbonded->setUseDispersionCorrection(true); nonbonded->setUseDispersionCorrection(false);
context.reinitialize(); context.reinitialize();
context.setPositions(positions); context.setPositions(positions);
energy2 = context.getState(State::Energy).getPotentialEnergy(); energy2 = context.getState(State::Energy).getPotentialEnergy();
...@@ -420,7 +420,7 @@ void testDispersionCorrection() { ...@@ -420,7 +420,7 @@ void testDispersionCorrection() {
term1 /= (numParticles*(numParticles+1))/2; term1 /= (numParticles*(numParticles+1))/2;
term2 /= (numParticles*(numParticles+1))/2; term2 /= (numParticles*(numParticles+1))/2;
expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize); expected = 8*M_PI*numParticles*numParticles*(term1-term2)/(boxSize*boxSize*boxSize);
ASSERT_EQUAL_TOL(expected, energy2-energy1, 1e-4); ASSERT_EQUAL_TOL(expected, energy1-energy2, 1e-4);
} }
int main() { int main() {
......
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