Skip to content

GitLab

Commit 21d67074 authored by Rossen Apostolov's avatar Rossen Apostolov
Browse files

Fixed a typo in Madelung energy formula.

parent fd64ac97
Hide whitespace changes
Inline Side-by-side
Showing with 91 additions and 93 deletions
platforms/reference/tests/TestReferenceEwald.cpp
View file @ 21d67074
...@@ -85,109 +85,20 @@ void testEwaldExact() { ...@@ -85,109 +85,20 @@ void testEwaldExact() {
const vector<Vec3>& forces = state.getForces(); const vector<Vec3>& forces = state.getForces();
// The potential energy of an ion in a crystal is // The potential energy of an ion in a crystal is
// E = - (M*e/ 4*pi*epsilon0*a0), // E = - (M*e^2/ 4*pi*epsilon0*a0),
// where // where
// M : Madelung constant (dimensionless, for FCC cells such as NaCl it is 1.7476) // M : Madelung constant (dimensionless, for FCC cells such as NaCl it is 1.7476)
// e : 1.6022 × 10−19 C // e : 1.6022 × 10−19 C
// 4*pi*epsilon0 : 1.112 × 10−10 C²/(J m) // 4*pi*epsilon0 : 1.112 × 10−10 C²/(J m)
// a0 : 0.282 x 10-9 m (perfect cell) // a0 : 0.282 x 10-9 m (perfect cell)
// Therefore // Therefore
double exactEnergy = -8.9290; // eV double exactEnergy = -14.3061e-19; // J (per ion pair)
exactEnergy = -4.4645; // eV per 1 ion exactEnergy = -7.1531e-19; // J per ion
exactEnergy = -7.15292069e-19; // J per ion
exactEnergy = -430.820; // kJ/mol per ion exactEnergy = -430.820; // kJ/mol per ion
ASSERT_EQUAL_TOL(exactEnergy * numParticles , state.getPotentialEnergy(), 100*TOL); ASSERT_EQUAL_TOL(exactEnergy * numParticles , state.getPotentialEnergy(), 100*TOL);
// Gromacs result // Gromacs result
// ASSERT_EQUAL_TOL(-430494.0, state.getPotentialEnergy(), 10*TOL); // ASSERT_EQUAL_TOL(-430494.0, state.getPotentialEnergy(), 10*TOL);
}
void testEwald2Ions() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
system.setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
nonbonded->setEwaldErrorTolerance(TOL);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(3.048000,2.764000,3.156000);
positions[1] = Vec3(2.809000,2.888000,2.571000);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 10*TOL);
}
void testWaterSystem() {
ReferencePlatform platform;
System system;
static int numParticles;
numParticles = 648;
for (int i = 0 ; i < numParticles ; i++)
{
system.addParticle(1.0);
}
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0 ; i < numParticles/3 ; i++)
{
nonbonded->addParticle(-0.82, 1, 0);
nonbonded->addParticle(0.41, 1, 0);
nonbonded->addParticle(0.41, 1, 0);
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 0.8;
nonbonded->setCutoffDistance(cutoff);
system.setPeriodicBoxVectors(Vec3(1.86206, 0, 0), Vec3(0, 1.86206, 0), Vec3(0, 0, 1.86206));
nonbonded->setEwaldErrorTolerance(TOL);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
#include "water.dat"
context.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state1.getForces();
// for (int i = 0 ; i < 42 ; i++)
// cout << "f [" << i << " : ]" << forces[i] << endl;
// cout << "PotentialEnergy: " << state1.getPotentialEnergy() << endl;
// Take a small step in the direction of the energy gradient.
double norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state1.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-3;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = state1.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
// See whether the potential energy changed by the expected amount.
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
State state2 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state1.getPotentialEnergy())/delta, 0.01)
} }
void testEwaldPME() { void testEwaldPME() {
...@@ -294,14 +205,101 @@ void testEwaldPME() { ...@@ -294,14 +205,101 @@ void testEwaldPME() {
} }
void testEwald2Ions() {
ReferencePlatform platform;
System system;
system.addParticle(1.0);
system.addParticle(1.0);
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
nonbonded->addParticle(1.0, 1, 0);
nonbonded->addParticle(-1.0, 1, 0);
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
const double cutoff = 2.0;
nonbonded->setCutoffDistance(cutoff);
system.setPeriodicBoxVectors(Vec3(6, 0, 0), Vec3(0, 6, 0), Vec3(0, 0, 6));
nonbonded->setEwaldErrorTolerance(TOL);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(2);
positions[0] = Vec3(3.048000,2.764000,3.156000);
positions[1] = Vec3(2.809000,2.888000,2.571000);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
ASSERT_EQUAL_VEC(Vec3(-123.711, 64.1877, -302.716), forces[0], 10*TOL);
ASSERT_EQUAL_VEC(Vec3(123.711, -64.1877, 302.716), forces[1], 10*TOL);
ASSERT_EQUAL_TOL(-217.276, state.getPotentialEnergy(), 10*TOL);
}
void testWaterSystem() {
ReferencePlatform platform;
System system;
static int numParticles;
numParticles = 648;
for (int i = 0 ; i < numParticles ; i++)
{
system.addParticle(1.0);
}
VerletIntegrator integrator(0.01);
NonbondedForce* nonbonded = new NonbondedForce();
for (int i = 0 ; i < numParticles/3 ; i++)
{
nonbonded->addParticle(-0.82, 1, 0);
nonbonded->addParticle(0.41, 1, 0);
nonbonded->addParticle(0.41, 1, 0);
}
nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
const double cutoff = 0.8;
nonbonded->setCutoffDistance(cutoff);
system.setPeriodicBoxVectors(Vec3(1.86206, 0, 0), Vec3(0, 1.86206, 0), Vec3(0, 0, 1.86206));
nonbonded->setEwaldErrorTolerance(TOL);
system.addForce(nonbonded);
Context context(system, integrator, platform);
vector<Vec3> positions(numParticles);
#include "water.dat"
context.setPositions(positions);
State state1 = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state1.getForces();
// for (int i = 0 ; i < 42 ; i++)
// cout << "f [" << i << " : ]" << forces[i] << endl;
// cout << "PotentialEnergy: " << state1.getPotentialEnergy() << endl;
// Take a small step in the direction of the energy gradient.
double norm = 0.0;
for (int i = 0; i < numParticles; ++i) {
Vec3 f = state1.getForces()[i];
norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
}
norm = std::sqrt(norm);
const double delta = 1e-3;
double step = delta/norm;
for (int i = 0; i < numParticles; ++i) {
Vec3 p = positions[i];
Vec3 f = state1.getForces()[i];
positions[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
}
context.setPositions(positions);
// See whether the potential energy changed by the expected amount.
nonbonded->setNonbondedMethod(NonbondedForce::Ewald);
State state2 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state1.getPotentialEnergy())/delta, 0.01)
}
int main() { int main() {
try { try {
testEwaldExact(); testEwaldExact();
testEwaldPME();
// testEwald2Ions(); // testEwald2Ions();
// testWaterSystem(); // testWaterSystem();
testEwaldPME();
} }
catch(const exception& e) { catch(const exception& e) {
cout << "exception: " << e.what() << endl; cout << "exception: " << e.what() << 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