Skip to content

GitLab

Commit e1e1e12d authored by peastman's avatar peastman
Browse files

Changed definition of polarizability to match CHARMM

parent 1aa2ccb7
Hide whitespace changes
Inline Side-by-side
Showing with 40 additions and 40 deletions
plugins/drude/platforms/cuda/src/CudaDrudeKernels.cpp
View file @ e1e1e12d
...@@ -126,9 +126,9 @@ void CudaCalcDrudeForceKernel::initialize(const System& system, const DrudeForce ...@@ -126,9 +126,9 @@ void CudaCalcDrudeForceKernel::initialize(const System& system, const DrudeForce
double a1 = (atoms[i][2] == -1 ? 1 : aniso12); double a1 = (atoms[i][2] == -1 ? 1 : aniso12);
double a2 = (atoms[i][3] == -1 || atoms[i][4] == -1 ? 1 : aniso34); double a2 = (atoms[i][3] == -1 || atoms[i][4] == -1 ? 1 : aniso34);
double a3 = 3-a1-a2; double a3 = 3-a1-a2;
double k3 = charge*charge/(polarizability*a3); double k3 = ONE_4PI_EPS0*charge*charge/(polarizability*a3);
double k1 = charge*charge/(polarizability*a1) - k3; double k1 = ONE_4PI_EPS0*charge*charge/(polarizability*a1) - k3;
double k2 = charge*charge/(polarizability*a2) - k3; double k2 = ONE_4PI_EPS0*charge*charge/(polarizability*a2) - k3;
if (atoms[i][2] == -1) { if (atoms[i][2] == -1) {
atoms[i][2] = 0; atoms[i][2] = 0;
k1 = 0; k1 = 0;
......
plugins/drude/platforms/cuda/tests/TestCudaDrudeForce.cpp
View file @ e1e1e12d
...@@ -76,9 +76,9 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg ...@@ -76,9 +76,9 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg
} }
void testSingleParticle() { void testSingleParticle() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
System system; System system;
system.addParticle(1.0); system.addParticle(1.0);
system.addParticle(1.0); system.addParticle(1.0);
...@@ -92,9 +92,9 @@ void testSingleParticle() { ...@@ -92,9 +92,9 @@ void testSingleParticle() {
} }
void testAnisotropicParticle() { void testAnisotropicParticle() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
const double a1 = 0.8; const double a1 = 0.8;
const double a2 = 1.1; const double a2 = 1.1;
const double k1 = k/a1; const double k1 = k/a1;
...@@ -124,9 +124,9 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) { ...@@ -124,9 +124,9 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) {
} }
void testThole() { void testThole() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
const double thole = 2.5; const double thole = 2.5;
System system; System system;
system.addParticle(1.0); system.addParticle(1.0);
...@@ -157,9 +157,9 @@ void testThole() { ...@@ -157,9 +157,9 @@ void testThole() {
} }
void testChangingParameters() { void testChangingParameters() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
Platform& platform = Platform::getPlatformByName("CUDA"); Platform& platform = Platform::getPlatformByName("CUDA");
// Create the system. // Create the system.
...@@ -184,9 +184,9 @@ void testChangingParameters() { ...@@ -184,9 +184,9 @@ void testChangingParameters() {
// Modify the parameters. // Modify the parameters.
const double k2 = 2.2; const double k2 = ONE_4PI_EPS0*2.2;
const double charge2 = 0.3; const double charge2 = 0.3;
const double alpha2 = charge2*charge2/k2; const double alpha2 = ONE_4PI_EPS0*charge2*charge2/k2;
drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1); drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
drude->updateParametersInContext(context); drude->updateParametersInContext(context);
state = context.getState(State::Energy); state = context.getState(State::Energy);
......
plugins/drude/platforms/opencl/src/OpenCLDrudeKernels.cpp
View file @ e1e1e12d
...@@ -131,9 +131,9 @@ void OpenCLCalcDrudeForceKernel::initialize(const System& system, const DrudeFor ...@@ -131,9 +131,9 @@ void OpenCLCalcDrudeForceKernel::initialize(const System& system, const DrudeFor
double a1 = (atoms[i][2] == -1 ? 1 : aniso12); double a1 = (atoms[i][2] == -1 ? 1 : aniso12);
double a2 = (atoms[i][3] == -1 || atoms[i][4] == -1 ? 1 : aniso34); double a2 = (atoms[i][3] == -1 || atoms[i][4] == -1 ? 1 : aniso34);
double a3 = 3-a1-a2; double a3 = 3-a1-a2;
double k3 = charge*charge/(polarizability*a3); double k3 = ONE_4PI_EPS0*charge*charge/(polarizability*a3);
double k1 = charge*charge/(polarizability*a1) - k3; double k1 = ONE_4PI_EPS0*charge*charge/(polarizability*a1) - k3;
double k2 = charge*charge/(polarizability*a2) - k3; double k2 = ONE_4PI_EPS0*charge*charge/(polarizability*a2) - k3;
if (atoms[i][2] == -1) { if (atoms[i][2] == -1) {
atoms[i][2] = 0; atoms[i][2] = 0;
k1 = 0; k1 = 0;
......
plugins/drude/platforms/opencl/tests/TestOpenCLDrudeForce.cpp
View file @ e1e1e12d
...@@ -76,9 +76,9 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg ...@@ -76,9 +76,9 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg
} }
void testSingleParticle() { void testSingleParticle() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
System system; System system;
system.addParticle(1.0); system.addParticle(1.0);
system.addParticle(1.0); system.addParticle(1.0);
...@@ -92,9 +92,9 @@ void testSingleParticle() { ...@@ -92,9 +92,9 @@ void testSingleParticle() {
} }
void testAnisotropicParticle() { void testAnisotropicParticle() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
const double a1 = 0.8; const double a1 = 0.8;
const double a2 = 1.1; const double a2 = 1.1;
const double k1 = k/a1; const double k1 = k/a1;
...@@ -124,9 +124,9 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) { ...@@ -124,9 +124,9 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) {
} }
void testThole() { void testThole() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
const double thole = 2.5; const double thole = 2.5;
System system; System system;
system.addParticle(1.0); system.addParticle(1.0);
...@@ -157,9 +157,9 @@ void testThole() { ...@@ -157,9 +157,9 @@ void testThole() {
} }
void testChangingParameters() { void testChangingParameters() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
Platform& platform = Platform::getPlatformByName("OpenCL"); Platform& platform = Platform::getPlatformByName("OpenCL");
// Create the system. // Create the system.
...@@ -184,9 +184,9 @@ void testChangingParameters() { ...@@ -184,9 +184,9 @@ void testChangingParameters() {
// Modify the parameters. // Modify the parameters.
const double k2 = 2.2; const double k2 = ONE_4PI_EPS0*2.2;
const double charge2 = 0.3; const double charge2 = 0.3;
const double alpha2 = charge2*charge2/k2; const double alpha2 = ONE_4PI_EPS0*charge2*charge2/k2;
drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1); drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
drude->updateParametersInContext(context); drude->updateParametersInContext(context);
state = context.getState(State::Energy); state = context.getState(State::Energy);
......
plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp
View file @ e1e1e12d
...@@ -136,9 +136,9 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include ...@@ -136,9 +136,9 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
RealOpenMM a1 = (p2 == -1 ? 1 : aniso12[i]); RealOpenMM a1 = (p2 == -1 ? 1 : aniso12[i]);
RealOpenMM a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34[i]); RealOpenMM a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34[i]);
RealOpenMM a3 = 3-a1-a2; RealOpenMM a3 = 3-a1-a2;
RealOpenMM k3 = charge[i]*charge[i]/(polarizability[i]*a3); RealOpenMM k3 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a3);
RealOpenMM k1 = charge[i]*charge[i]/(polarizability[i]*a1) - k3; RealOpenMM k1 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a1) - k3;
RealOpenMM k2 = charge[i]*charge[i]/(polarizability[i]*a2) - k3; RealOpenMM k2 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a2) - k3;
// Compute the isotropic force. // Compute the isotropic force.
......
plugins/drude/platforms/reference/tests/TestReferenceDrudeForce.cpp
View file @ e1e1e12d
...@@ -76,9 +76,9 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg ...@@ -76,9 +76,9 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg
} }
void testSingleParticle() { void testSingleParticle() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
System system; System system;
system.addParticle(1.0); system.addParticle(1.0);
system.addParticle(1.0); system.addParticle(1.0);
...@@ -92,9 +92,9 @@ void testSingleParticle() { ...@@ -92,9 +92,9 @@ void testSingleParticle() {
} }
void testAnisotropicParticle() { void testAnisotropicParticle() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
const double a1 = 0.8; const double a1 = 0.8;
const double a2 = 1.1; const double a2 = 1.1;
const double k1 = k/a1; const double k1 = k/a1;
...@@ -124,9 +124,9 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) { ...@@ -124,9 +124,9 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) {
} }
void testThole() { void testThole() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
const double thole = 2.5; const double thole = 2.5;
System system; System system;
system.addParticle(1.0); system.addParticle(1.0);
...@@ -157,9 +157,9 @@ void testThole() { ...@@ -157,9 +157,9 @@ void testThole() {
} }
void testChangingParameters() { void testChangingParameters() {
const double k = 1.5; const double k = ONE_4PI_EPS0*1.5;
const double charge = 0.1; const double charge = 0.1;
const double alpha = charge*charge/k; const double alpha = ONE_4PI_EPS0*charge*charge/k;
Platform& platform = Platform::getPlatformByName("Reference"); Platform& platform = Platform::getPlatformByName("Reference");
// Create the system. // Create the system.
...@@ -184,9 +184,9 @@ void testChangingParameters() { ...@@ -184,9 +184,9 @@ void testChangingParameters() {
// Modify the parameters. // Modify the parameters.
const double k2 = 2.2; const double k2 = ONE_4PI_EPS0*2.2;
const double charge2 = 0.3; const double charge2 = 0.3;
const double alpha2 = charge2*charge2/k2; const double alpha2 = ONE_4PI_EPS0*charge2*charge2/k2;
drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1); drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
drude->updateParametersInContext(context); drude->updateParametersInContext(context);
state = context.getState(State::Energy); state = context.getState(State::Energy);
......
wrappers/python/simtk/openmm/app/data/swm4ndp.xml
View file @ e1e1e12d
...@@ -30,7 +30,7 @@ ...@@ -30,7 +30,7 @@
<Atom type="swm4ndp-OD" charge="-1.71636" sigma="1" epsilon="0"/> <Atom type="swm4ndp-OD" charge="-1.71636" sigma="1" epsilon="0"/>
</NonbondedForce> </NonbondedForce>
<DrudeForce> <DrudeForce>
<Particle type1="swm4ndp-OD" type2="swm4ndp-O" charge="-1.71636" polarizability="7.040850e-6" thole="1.3"/> <Particle type1="swm4ndp-OD" type2="swm4ndp-O" charge="-1.71636" polarizability="0.000978253" thole="1.3"/>
</DrudeForce> </DrudeForce>
</ForceField> </ForceField>
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