Skip to content

GitLab

Commit 0180a46c authored by Peter Eastman's avatar Peter Eastman
Browse files

Fixed compilation errors under Windows in reference Ewald summation code.

parent e702e8bb
Show whitespace changes
Inline Side-by-side
Showing with 95 additions and 21 deletions
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ 0180a46c
......@@ -32,6 +32,8 @@
#include "ReferenceLJCoulombIxn.h"
#include "ReferenceForce.h"
double erfc(double x);
/**---------------------------------------------------------------------------------------
ReferenceLJCoulombIxn constructor
......@@ -240,8 +242,8 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
int numRz = 60+1;
int kmax = std::max(numRx, std::max(numRy,numRz));
static const RealOpenMM alphaEwald = 3.123413;
RealOpenMM factorEwald = -1.0 / (4*alphaEwald*alphaEwald);
static const RealOpenMM alphaEwald = (RealOpenMM) 3.123413;
RealOpenMM factorEwald = -1 / (4*alphaEwald*alphaEwald);
RealOpenMM SQRT_PI = sqrt(PI);
......@@ -249,7 +251,7 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
static const RealOpenMM epsilon = 1.0;
static const RealOpenMM one = 1.0;
RealOpenMM recipCoeff = 4.0*M_PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon;
RealOpenMM recipCoeff = (RealOpenMM) 4.0*M_PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon;
RealOpenMM selfEwaldEnergy = 0.0;
RealOpenMM realSpaceEwaldEnergy = 0.0;
......@@ -275,9 +277,10 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
// setup K-vectors
d_complex eir[kmax][numberOfAtoms][3];
d_complex tab_xy[numberOfAtoms];
d_complex tab_qxyz[numberOfAtoms];
#define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
d_complex* eir = new d_complex[kmax*numberOfAtoms*3];
d_complex* tab_xy = new d_complex[numberOfAtoms];
d_complex* tab_qxyz = new d_complex[numberOfAtoms];
d_complex a,b,c;
int i,j,m;
......@@ -290,19 +293,19 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
for(i = 0; (i < numberOfAtoms); i++) {
for(m = 0; (m < 3); m++) {
eir[0][i][m].real() = 1;
eir[0][i][m].imag() = 0;
EIR(0, i, m).real() = 1;
EIR(0, i, m).imag() = 0;
}
for(m=0; (m<3); m++) {
eir[1][i][m].real() = cos(atomCoordinates[i][m]*recipBoxSize[m]);
eir[1][i][m].imag() = sin(atomCoordinates[i][m]*recipBoxSize[m]);
EIR(1, i, m).real() = cos(atomCoordinates[i][m]*recipBoxSize[m]);
EIR(1, i, m).imag() = sin(atomCoordinates[i][m]*recipBoxSize[m]);
}
for(j=2; (j<kmax); j++) {
for(m=0; (m<3); m++) {
eir[j][i][m] = eir[j-1][i][m] * eir[1][i][m];
EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
}
}
}
......@@ -322,13 +325,13 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
if(ry >= 0) {
for(int n = 0; n < numberOfAtoms; n++) {
tab_xy[n] = eir[rx][n][0] * eir[ry][n][1];
tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
}
}
else {
for(int n = 0; n < numberOfAtoms; n++) {
tab_xy[n]= eir[rx][n][0] * conj (eir[-ry][n][1]);
tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
}
}
......@@ -337,16 +340,16 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
if( rz >= 0) {
for( int n = 0; n < numberOfAtoms; n++) {
tab_qxyz[n].real() = atomParameters[n][QIndex] * (tab_xy[n] * eir[rz][n][2]).real();
tab_qxyz[n].imag() = atomParameters[n][QIndex] * (tab_xy[n] * eir[rz][n][2]).imag();
tab_qxyz[n].real() = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2)).real();
tab_qxyz[n].imag() = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2)).imag();
}
}
else {
for( int n = 0; n < numberOfAtoms; n++) {
tab_qxyz[n].real() = atomParameters[n][QIndex] * (tab_xy[n] * conj(eir[-rz][n][2])).real() ;
tab_qxyz[n].imag() = atomParameters[n][QIndex] * (tab_xy[n] * conj(eir[-rz][n][2])).imag() ;
tab_qxyz[n].real() = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2))).real() ;
tab_qxyz[n].imag() = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2))).imag() ;
}
}
......@@ -362,16 +365,16 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM kz = rz * recipBoxSize[2];
RealOpenMM k2 = kx * kx + ky * ky + kz * kz;
RealOpenMM ak = exp(k2*factorEwald) / k2;
RealOpenMM akv = 2.0 * ak * (1.0/k2 - factorEwald);
RealOpenMM akv = 2 * ak * (1/k2 - factorEwald);
recipEnergy += ak * ( cs * cs + ss * ss);
RealOpenMM vir = akv * ( cs * cs + ss * ss);
for(int n = 0; n < numberOfAtoms; n++) {
RealOpenMM force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
forces[n][0] += 2.0 * recipCoeff * force * kx ;
forces[n][1] += 2.0 * recipCoeff * force * ky ;
forces[n][2] += 2.0 * recipCoeff * force * kz ;
forces[n][0] += 2 * recipCoeff * force * kx ;
forces[n][1] += 2 * recipCoeff * force * ky ;
forces[n][2] += 2 * recipCoeff * force * kz ;
}
lowrz = 1 - numRz;
}
......@@ -442,6 +445,9 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
}
delete[] exclusionIndices;
delete[] eir;
delete[] tab_xy;
delete[] tab_qxyz;
// ***********************************************************************
......
platforms/reference/src/SimTKReference/erfc.cpp 0 → 100644
View file @ 0180a46c
/***************************
* erf.cpp
* author: Steve Strand
* written: 29-Jan-04
***************************/
#include <math.h>
static const double rel_error= 1E-12; //calculate 12 significant figures
//you can adjust rel_error to trade off between accuracy and speed
//but don't ask for > 15 figures (assuming usual 52 bit mantissa in a double)
double erf(double x)
//erf(x) = 2/sqrt(pi)*integral(exp(-t^2),t,0,x)
// = 2/sqrt(pi)*[x - x^3/3 + x^5/5*2! - x^7/7*3! + ...]
// = 1-erfc(x)
{
static const double two_sqrtpi= 1.128379167095512574; // 2/sqrt(pi)
if (fabs(x) > 2.2) {
return 1.0 - erfc(x); //use continued fraction when fabs(x) > 2.2
}
double sum= x, term= x, xsqr= x*x;
int j= 1;
do {
term*= xsqr/j;
sum-= term/(2*j+1);
++j;
term*= xsqr/j;
sum+= term/(2*j+1);
++j;
} while (fabs(term)/sum > rel_error);
return two_sqrtpi*sum;
}
double erfc(double x)
//erfc(x) = 2/sqrt(pi)*integral(exp(-t^2),t,x,inf)
// = exp(-x^2)/sqrt(pi) * [1/x+ (1/2)/x+ (2/2)/x+ (3/2)/x+ (4/2)/x+ ...]
// = 1-erf(x)
//expression inside [] is a continued fraction so '+' means add to denominator only
{
static const double one_sqrtpi= 0.564189583547756287; // 1/sqrt(pi)
if (fabs(x) < 2.2) {
return 1.0 - erf(x); //use series when fabs(x) < 2.2
}
if (signbit(x)) { //continued fraction only valid for x>0
return 2.0 - erfc(-x);
}
double a=1, b=x; //last two convergent numerators
double c=x, d=x*x+0.5; //last two convergent denominators
double q1, q2= b/d; //last two convergents (a/c and b/d)
double n= 1.0, t;
do {
t= a*n+b*x;
a= b;
b= t;
t= c*n+d*x;
c= d;
d= t;
n+= 0.5;
q1= q2;
q2= b/d;
} while (fabs(q1-q2)/q2 > rel_error);
return one_sqrtpi*exp(-x*x)*q2;
}
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