Skip to content

GitLab

Commit 6186c067 authored by Peter Eastman's avatar Peter Eastman
Browse files

Fixed compiler warnings

parent cd418100
Hide whitespace changes
Inline Side-by-side
Showing with 431 additions and 431 deletions
platforms/reference/src/SimTKReference/PME.cpp
View file @ 6186c067
......@@ -42,24 +42,24 @@ typedef int ivec[3];
struct pme
{
int natoms;
RealOpenMM ewaldcoeff;
int natoms;
RealOpenMM ewaldcoeff;
t_complex * grid; /* Memory for the grid we spread charges on.
t_complex * grid; /* Memory for the grid we spread charges on.
* Element (i,j,k) is accessed as:
* grid[i*ngrid[1]*ngrid[2] + j*ngrid[2] + k]
*/
int ngrid[3]; /* Total grid dimensions (all data is complex!) */
fftpack_t fftplan; /* Handle to fourier transform setup */
int ngrid[3]; /* Total grid dimensions (all data is complex!) */
fftpack_t fftplan; /* Handle to fourier transform setup */
int order; /* PME interpolation order. Almost always 4 */
int order; /* PME interpolation order. Almost always 4 */
/* Data for bspline interpolation, see the Essman PME paper */
RealOpenMM * bsplines_moduli[3]; /* 3 pointers, to x/y/z bspline moduli, each of length ngrid[x/y/z] */
RealOpenMM * bsplines_theta[3]; /* each of x/y/z has length order*natoms */
RealOpenMM * bsplines_dtheta[3]; /* each of x/y/z has length order*natoms */
RealOpenMM * bsplines_moduli[3]; /* 3 pointers, to x/y/z bspline moduli, each of length ngrid[x/y/z] */
RealOpenMM * bsplines_theta[3]; /* each of x/y/z has length order*natoms */
RealOpenMM * bsplines_dtheta[3]; /* each of x/y/z has length order*natoms */
ivec * particleindex; /* Array of length natoms. Each element is
ivec * particleindex; /* Array of length natoms. Each element is
* an ivec (3 ints) that specify the grid
* indices for that particular atom. Updated every step!
*/
......@@ -83,7 +83,7 @@ struct pme
* the central cell, i.e., in this case we would assume all coordinates fall in -10 nm < x,y,z < 20 nm.
*/
RealOpenMM epsilon_r; /* Dielectric coefficient to use, typically 1.0 */
RealOpenMM epsilon_r; /* Dielectric coefficient to use, typically 1.0 */
};
......@@ -95,128 +95,128 @@ struct pme
static void
pme_calculate_bsplines_moduli(pme_t pme)
{
int nmax;
int i,j,k,l,d;
int order;
int ndata;
RealOpenMM * data;
RealOpenMM * ddata;
RealOpenMM * bsplines_data;
RealOpenMM div;
RealOpenMM sc,ss,arg;
nmax = 0;
for(d=0;d<3;d++)
{
nmax = (pme->ngrid[d] > nmax) ? pme->ngrid[d] : nmax;
pme->bsplines_moduli[d] = (RealOpenMM *) malloc(sizeof(RealOpenMM)*pme->ngrid[d]);
}
order = pme->order;
/* temp storage in this routine */
data = (RealOpenMM *) malloc(sizeof(RealOpenMM)*order);
ddata = (RealOpenMM *) malloc(sizeof(RealOpenMM)*order);
bsplines_data = (RealOpenMM *) malloc(sizeof(RealOpenMM)*nmax);
data[order-1]=0;
data[1]=0;
data[0]=1;
for(k=3;k<order;k++)
{
div=1.0/(k-1.0);
data[k-1]=0;
for(l=1;l<(k-1);l++)
{
data[k-l-1]=div*(l*data[k-l-2]+(k-l)*data[k-l-1]);
}
data[0]=div*data[0];
}
/* differentiate */
ddata[0]=-data[0];
for(k=1;k<order;k++)
{
ddata[k]=data[k-1]-data[k];
}
div=1.0/(order-1);
data[order-1]=0;
for(l=1;l<(order-1);l++)
{
data[order-l-1]=div*(l*data[order-l-2]+(order-l)*data[order-l-1]);
}
data[0]=div*data[0];
for(i=0;i<nmax;i++)
{
bsplines_data[i]=0;
}
for(i=1;i<=order;i++)
{
bsplines_data[i]=data[i-1];
int nmax;
int i,j,k,l,d;
int order;
int ndata;
RealOpenMM * data;
RealOpenMM * ddata;
RealOpenMM * bsplines_data;
RealOpenMM div;
RealOpenMM sc,ss,arg;
nmax = 0;
for(d=0;d<3;d++)
{
nmax = (pme->ngrid[d] > nmax) ? pme->ngrid[d] : nmax;
pme->bsplines_moduli[d] = (RealOpenMM *) malloc(sizeof(RealOpenMM)*pme->ngrid[d]);
}
/* Evaluate the actual bspline moduli for X/Y/Z */
for(d=0;d<3;d++)
{
ndata = pme->ngrid[d];
for(i=0;i<ndata;i++)
{
sc=ss=0;
for(j=0;j<ndata;j++)
{
arg=(2.0*M_PI*i*j)/ndata;
sc+=bsplines_data[j]*cos(arg);
ss+=bsplines_data[j]*sin(arg);
}
pme->bsplines_moduli[d][i]=sc*sc+ss*ss;
}
for(i=0;i<ndata;i++)
{
if(pme->bsplines_moduli[d][i]<1.0e-7)
{
pme->bsplines_moduli[d][i]=(pme->bsplines_moduli[d][i-1]+pme->bsplines_moduli[d][i+1])*0.5;
}
}
}
/* Release temp storage */
free(data);
free(ddata);
free(bsplines_data);
order = pme->order;
/* temp storage in this routine */
data = (RealOpenMM *) malloc(sizeof(RealOpenMM)*order);
ddata = (RealOpenMM *) malloc(sizeof(RealOpenMM)*order);
bsplines_data = (RealOpenMM *) malloc(sizeof(RealOpenMM)*nmax);
data[order-1]=0;
data[1]=0;
data[0]=1;
for(k=3;k<order;k++)
{
div=(RealOpenMM) (1.0/(k-1.0));
data[k-1]=0;
for(l=1;l<(k-1);l++)
{
data[k-l-1]=div*(l*data[k-l-2]+(k-l)*data[k-l-1]);
}
data[0]=div*data[0];
}
/* differentiate */
ddata[0]=-data[0];
for(k=1;k<order;k++)
{
ddata[k]=data[k-1]-data[k];
}
div=(RealOpenMM) (1.0/(order-1));
data[order-1]=0;
for(l=1;l<(order-1);l++)
{
data[order-l-1]=div*(l*data[order-l-2]+(order-l)*data[order-l-1]);
}
data[0]=div*data[0];
for(i=0;i<nmax;i++)
{
bsplines_data[i]=0;
}
for(i=1;i<=order;i++)
{
bsplines_data[i]=data[i-1];
}
/* Evaluate the actual bspline moduli for X/Y/Z */
for(d=0;d<3;d++)
{
ndata = pme->ngrid[d];
for(i=0;i<ndata;i++)
{
sc=ss=0;
for(j=0;j<ndata;j++)
{
arg=(RealOpenMM) ((2.0*M_PI*i*j)/ndata);
sc+=bsplines_data[j]*cos(arg);
ss+=bsplines_data[j]*sin(arg);
}
pme->bsplines_moduli[d][i]=sc*sc+ss*ss;
}
for(i=0;i<ndata;i++)
{
if(pme->bsplines_moduli[d][i]<1.0e-7)
{
pme->bsplines_moduli[d][i]=(pme->bsplines_moduli[d][i-1]+pme->bsplines_moduli[d][i+1])/2;
}
}
}
/* Release temp storage */
free(data);
free(ddata);
free(bsplines_data);
}
static void
pme_update_grid_index_and_fraction(pme_t pme,
RealOpenMM ** atomCoordinates,
const RealOpenMM periodicBoxSize[3])
RealOpenMM ** atomCoordinates,
const RealOpenMM periodicBoxSize[3])
{
int i;
int d;
RealOpenMM t;
int ti;
for(i=0;i<pme->natoms;i++)
{
for(d=0;d<3;d++)
{
/* Index calculation (Look mom, no conditionals!):
int i;
int d;
RealOpenMM t;
int ti;
for(i=0;i<pme->natoms;i++)
{
for(d=0;d<3;d++)
{
/* Index calculation (Look mom, no conditionals!):
*
* Both for Cuda and modern CPUs it is nice to avoid conditionals, but we still need to apply periodic boundary conditions.
* Instead of having loops to add/subtract the box dimension, we do it this way:
*
* 1. First add the box size, to make sure this atom coordinate isnt -0.1 or something.
* 1. First add the box size, to make sure this atom coordinate isnt -0.1 or something.
* After this we assume all fractional box positions are *positive*.
* The reason for this is that we always want to round coordinates _down_ to get
* their grid index, and when taking the integer part of -3.4 we would get -3, not -4 as we want.
* Since we anyway need the grid indices to fall in the central box, it is more convenient
* to first manipulate the coordinates to be positive.
* 2. Convert to integer grid index
* 2. Convert to integer grid index
* Since we have added a whole box unit in step 1, this index might actually be larger than
* the grid dimension. Examples, assuming 10*10*10nm box and grid dimension 100*100*100 (spacing 0.1 nm):
*
......@@ -230,7 +230,7 @@ pme_update_grid_index_and_fraction(pme_t pme,
*
* The fraction is calculates as t-ti, which becomes { 0.43 , 0.35 , 0.7 }
*
* 3. Take the first integer index part (which can be larger than the grid) modulo the grid dimension
* 3. Take the first integer index part (which can be larger than the grid) modulo the grid dimension
*
* Now we get { 5 , 62 , 92 }
*
......@@ -240,14 +240,14 @@ pme_update_grid_index_and_fraction(pme_t pme,
* In practice, MD programs will apply PBC to reset particles inside the central box to avoid
* numerical problems, so this shouldnt cause any problems.
* (And, by adding 100.0 box lengths, we would lose a bit of numerical accuracy here!)
*/
t = (atomCoordinates[i][d] / periodicBoxSize[d] + 1.0)*pme->ngrid[d];
ti = (int) t;
pme->particlefraction[i][d] = t - ti;
pme->particleindex[i][d] = ti % pme->ngrid[d];
}
}
*/
t = (atomCoordinates[i][d] / periodicBoxSize[d] + 1)*pme->ngrid[d];
ti = (int) t;
pme->particlefraction[i][d] = t - ti;
pme->particleindex[i][d] = ti % pme->ngrid[d];
}
}
}
......@@ -259,96 +259,96 @@ pme_update_grid_index_and_fraction(pme_t pme,
static void
pme_update_bsplines(pme_t pme)
{
int i,j,k,l;
int order;
RealOpenMM dr,div;
RealOpenMM * data;
RealOpenMM * ddata;
int i,j,k,l;
int order;
RealOpenMM dr,div;
RealOpenMM * data;
RealOpenMM * ddata;
order = pme->order;
order = pme->order;
for(i=0; (i<pme->natoms); i++)
{
for(j=0; j<3; j++)
{
/* dr is relative offset from lower cell limit */
dr = pme->particlefraction[i][j];
data = &(pme->bsplines_theta[j][i*order]);
ddata = &(pme->bsplines_dtheta[j][i*order]);
data[order-1] = 0;
data[1] = dr;
data[0] = 1-dr;
for(k=3; k<order; k++)
{
div = 1.0/(k-1.0);
data[k-1] = div*dr*data[k-2];
for(l=1; l<(k-1); l++)
{
data[k-l-1] = div*((dr+l)*data[k-l-2]+(k-l-dr)*data[k-l-1]);
}
data[0] = div*(1-dr)*data[0];
}
/* differentiate */
ddata[0] = -data[0];
for(k=1; k<order; k++)
{
ddata[k] = data[k-1]-data[k];
}
div = 1.0/(order-1);
data[order-1] = div*dr*data[order-2];
for(l=1; l<(order-1); l++)
{
data[order-l-1] = div*((dr+l)*data[order-l-2]+(order-l-dr)*data[order-l-1]);
}
data[0] = div*(1-dr)*data[0];
}
{
for(j=0; j<3; j++)
{
/* dr is relative offset from lower cell limit */
dr = pme->particlefraction[i][j];
data = &(pme->bsplines_theta[j][i*order]);
ddata = &(pme->bsplines_dtheta[j][i*order]);
data[order-1] = 0;
data[1] = dr;
data[0] = 1-dr;
for(k=3; k<order; k++)
{
div = (RealOpenMM) (1.0/(k-1.0));
data[k-1] = div*dr*data[k-2];
for(l=1; l<(k-1); l++)
{
data[k-l-1] = div*((dr+l)*data[k-l-2]+(k-l-dr)*data[k-l-1]);
}
data[0] = div*(1-dr)*data[0];
}
/* differentiate */
ddata[0] = -data[0];
for(k=1; k<order; k++)
{
ddata[k] = data[k-1]-data[k];
}
div = (RealOpenMM) (1.0/(order-1));
data[order-1] = div*dr*data[order-2];
for(l=1; l<(order-1); l++)
{
data[order-l-1] = div*((dr+l)*data[order-l-2]+(order-l-dr)*data[order-l-1]);
}
data[0] = div*(1-dr)*data[0];
}
}
}
static void
pme_grid_spread_charge(pme_t pme,
RealOpenMM ** atomParameters)
RealOpenMM ** atomParameters)
{
static const int QIndex = 2; // atom charges are stored in atomParameters[atomID][2]
int order;
int i;
int ix,iy,iz;
int x0index,y0index,z0index;
int xindex,yindex,zindex;
int index;
RealOpenMM q;
RealOpenMM * thetax;
RealOpenMM * thetay;
RealOpenMM * thetaz;
order = pme->order;
int order;
int i;
int ix,iy,iz;
int x0index,y0index,z0index;
int xindex,yindex,zindex;
int index;
RealOpenMM q;
RealOpenMM * thetax;
RealOpenMM * thetay;
RealOpenMM * thetaz;
order = pme->order;
/* Reset the grid */
for(i=0;i<pme->ngrid[0]*pme->ngrid[1]*pme->ngrid[2];i++)
{
pme->grid[i].re = pme->grid[i].im = 0;
}
for(i=0;i<pme->ngrid[0]*pme->ngrid[1]*pme->ngrid[2];i++)
{
pme->grid[i].re = pme->grid[i].im = 0;
}
for(i=0;i<pme->natoms;i++)
{
q = atomParameters[i][QIndex];
for(i=0;i<pme->natoms;i++)
{
q = atomParameters[i][QIndex];
/* Grid index for the actual atom position */
x0index = pme->particleindex[i][0];
y0index = pme->particleindex[i][1];
z0index = pme->particleindex[i][2];
x0index = pme->particleindex[i][0];
y0index = pme->particleindex[i][1];
z0index = pme->particleindex[i][2];
/* Bspline factors for this atom in each dimension , calculated from fractional coordinates */
thetax = &(pme->bsplines_theta[0][i*order]);
thetay = &(pme->bsplines_theta[1][i*order]);
thetaz = &(pme->bsplines_theta[2][i*order]);
thetax = &(pme->bsplines_theta[0][i*order]);
thetay = &(pme->bsplines_theta[1][i*order]);
thetaz = &(pme->bsplines_theta[2][i*order]);
/* Loop over norder*norder*norder (typically 4*4*4) neighbor cells.
*
......@@ -368,90 +368,90 @@ pme_grid_spread_charge(pme_t pme,
* 3) When we parallelize things, we only need to communicate in one direction instead of two!
*/
for(ix=0;ix<order;ix++)
{
for(ix=0;ix<order;ix++)
{
/* Calculate index, apply PBC so we spread to index 0/1/2 when a particle is close to the upper limit of the grid */
xindex = (x0index + ix) % pme->ngrid[0];
xindex = (x0index + ix) % pme->ngrid[0];
for(iy=0;iy<order;iy++)
{
yindex = (y0index + iy) % pme->ngrid[1];
for(iy=0;iy<order;iy++)
{
yindex = (y0index + iy) % pme->ngrid[1];
for(iz=0;iz<order;iz++)
{
/* Can be optimized, but we keep it simple here */
zindex = (z0index + iz) % pme->ngrid[2];
for(iz=0;iz<order;iz++)
{
/* Can be optimized, but we keep it simple here */
zindex = (z0index + iz) % pme->ngrid[2];
/* Calculate index in the charge grid */
index = xindex*pme->ngrid[1]*pme->ngrid[2] + yindex*pme->ngrid[2] + zindex;
index = xindex*pme->ngrid[1]*pme->ngrid[2] + yindex*pme->ngrid[2] + zindex;
/* Add the charge times the bspline spread/interpolation factors to this grid position */
pme->grid[index].re += q*thetax[ix]*thetay[iy]*thetaz[iz];
}
}
}
}
pme->grid[index].re += q*thetax[ix]*thetay[iy]*thetaz[iz];
}
}
}
}
}
static void
pme_reciprocal_convolution(pme_t pme,
const RealOpenMM periodicBoxSize[3],
RealOpenMM * energy,
RealOpenMM pme_virial[3][3])
const RealOpenMM periodicBoxSize[3],
RealOpenMM * energy,
RealOpenMM pme_virial[3][3])
{
int kx,ky,kz;
int nx,ny,nz;
RealOpenMM mx,my,mz;
RealOpenMM mhx,mhy,mhz,m2;
RealOpenMM one_4pi_eps;
RealOpenMM virxx,virxy,virxz,viryy,viryz,virzz;
RealOpenMM bx,by,bz;
RealOpenMM d1,d2;
RealOpenMM eterm,vfactor,struct2,ets2;
RealOpenMM esum;
RealOpenMM factor;
RealOpenMM denom;
RealOpenMM boxfactor;
RealOpenMM maxkx,maxky,maxkz;
t_complex *ptr;
nx = pme->ngrid[0];
ny = pme->ngrid[1];
nz = pme->ngrid[2];
one_4pi_eps=ONE_4PI_EPS0/pme->epsilon_r;
factor=M_PI*M_PI/(pme->ewaldcoeff*pme->ewaldcoeff);
boxfactor = M_PI*periodicBoxSize[0]*periodicBoxSize[1]*periodicBoxSize[2];
esum = 0;
virxx = 0;
virxy = 0;
virxz = 0;
viryy = 0;
viryz = 0;
virzz = 0;
maxkx = (nx+1)/2;
maxky = (ny+1)/2;
maxkz = (nz+1)/2;
int kx,ky,kz;
int nx,ny,nz;
RealOpenMM mx,my,mz;
RealOpenMM mhx,mhy,mhz,m2;
RealOpenMM one_4pi_eps;
RealOpenMM virxx,virxy,virxz,viryy,viryz,virzz;
RealOpenMM bx,by,bz;
RealOpenMM d1,d2;
RealOpenMM eterm,vfactor,struct2,ets2;
RealOpenMM esum;
RealOpenMM factor;
RealOpenMM denom;
RealOpenMM boxfactor;
RealOpenMM maxkx,maxky,maxkz;
t_complex *ptr;
nx = pme->ngrid[0];
ny = pme->ngrid[1];
nz = pme->ngrid[2];
one_4pi_eps = (RealOpenMM) (ONE_4PI_EPS0/pme->epsilon_r);
factor = (RealOpenMM) (M_PI*M_PI/(pme->ewaldcoeff*pme->ewaldcoeff));
boxfactor = (RealOpenMM) (M_PI*periodicBoxSize[0]*periodicBoxSize[1]*periodicBoxSize[2]);
esum = 0;
virxx = 0;
virxy = 0;
virxz = 0;
viryy = 0;
viryz = 0;
virzz = 0;
maxkx = (RealOpenMM) ((nx+1)/2);
maxky = (RealOpenMM) ((ny+1)/2);
maxkz = (RealOpenMM) ((nz+1)/2);
for(kx=0;kx<nx;kx++)
{
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
mx = (kx<maxkx) ? kx : (kx-nx);
mx = (RealOpenMM) ((kx<maxkx) ? kx : (kx-nx));
mhx = mx/periodicBoxSize[0];
bx = boxfactor*pme->bsplines_moduli[0][kx];
for(ky=0;ky<ny;ky++)
{
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
my = (ky<maxky) ? ky : (ky-ny);
my = (RealOpenMM) ((ky<maxky) ? ky : (ky-ny));
mhy = my/periodicBoxSize[1];
by = pme->bsplines_moduli[1][ky];
for(kz=0;kz<nz;kz++)
{
for(kz=0;kz<nz;kz++)
{
/* If the net charge of the system is 0.0, there will not be any DC (direct current, zero frequency) component. However,
* we can still handle charged systems through a charge correction, in which case the DC
* component should be excluded from recprocal space. We will anyway run into problems below when dividing with the
......@@ -460,159 +460,159 @@ pme_reciprocal_convolution(pme_t pme,
* In cuda you could probably work around this by setting something to 0.0 instead, but the short story is that we
* should skip the zero frequency case!
*/
if (kx==0 && ky==0 && kz==0)
{
continue;
}
if (kx==0 && ky==0 && kz==0)
{
continue;
}
/* Calculate frequency. Grid indices in the upper half correspond to negative frequencies! */
mz = (kz<maxkz) ? kz : (kz-nz);
mhz = mz/periodicBoxSize[2];
mz = (RealOpenMM) ((kz<maxkz) ? kz : (kz-nz));
mhz = mz/periodicBoxSize[2];
/* Pointer to the grid cell in question */
ptr = pme->grid + kx*ny*nz + ky*nz + kz;
ptr = pme->grid + kx*ny*nz + ky*nz + kz;
/* Get grid data for this frequency */
d1 = ptr->re;
d2 = ptr->im;
d1 = ptr->re;
d2 = ptr->im;
/* Calculate the convolution - see the Essman/Darden paper for the equation! */
m2 = mhx*mhx+mhy*mhy+mhz*mhz;
bz = pme->bsplines_moduli[2][kz];
denom = m2*bx*by*bz;
m2 = mhx*mhx+mhy*mhy+mhz*mhz;
bz = pme->bsplines_moduli[2][kz];
denom = m2*bx*by*bz;
eterm = one_4pi_eps*exp(-factor*m2)/denom;
eterm = one_4pi_eps*exp(-factor*m2)/denom;
/* write back convolution data to grid */
ptr->re = d1*eterm;
ptr->im = d2*eterm;
ptr->re = d1*eterm;
ptr->im = d2*eterm;
struct2 = (d1*d1+d2*d2);
struct2 = (d1*d1+d2*d2);
/* Long-range PME contribution to the energy for this frequency */
ets2 = eterm*struct2;
esum += ets2;
/* Long-range PME contribution to the energy for this frequency */
ets2 = eterm*struct2;
esum += ets2;
/* PME long-range contribution to atomic virial. Since it is symmetric, we only calculate half the matrix inside this loop. */
vfactor = (factor*m2+1.0)*2.0/m2;
virxx += ets2*(vfactor*mhx*mhx-1.0);
/* PME long-range contribution to atomic virial. Since it is symmetric, we only calculate half the matrix inside this loop. */
vfactor = (factor*m2+1)*2/m2;
virxx += ets2*(vfactor*mhx*mhx-1);
virxy += ets2*vfactor*mhx*mhy;
virxz += ets2*vfactor*mhx*mhz;
viryy += ets2*(vfactor*mhy*mhy-1.0);
viryy += ets2*(vfactor*mhy*mhy-1);
viryz += ets2*vfactor*mhy*mhz;
virzz += ets2*(vfactor*mhz*mhz-1.0);
}
}
}
pme_virial[0][0] = 0.25*virxx;
pme_virial[1][1] = 0.25*viryy;
pme_virial[2][2] = 0.25*virzz;
pme_virial[0][1] = pme_virial[1][0] = 0.25*virxy;
pme_virial[0][2] = pme_virial[2][0] = 0.25*virxz;
pme_virial[1][2] = pme_virial[2][1] = 0.25*viryz;
virzz += ets2*(vfactor*mhz*mhz-1);
}
}
}
pme_virial[0][0] = (RealOpenMM) (0.25*virxx);
pme_virial[1][1] = (RealOpenMM) (0.25*viryy);
pme_virial[2][2] = (RealOpenMM) (0.25*virzz);
pme_virial[0][1] = pme_virial[1][0] = (RealOpenMM) (0.25*virxy);
pme_virial[0][2] = pme_virial[2][0] = (RealOpenMM) (0.25*virxz);
pme_virial[1][2] = pme_virial[2][1] = (RealOpenMM) (0.25*viryz);
/* The factor 0.5 is nothing special, but it is better to have it here than inside the loop :-) */
*energy = 0.5*esum;
*energy = 0.5*esum;
}
static void
pme_grid_interpolate_force(pme_t pme,
const RealOpenMM periodicBoxSize[3],
RealOpenMM ** atomParameters,
RealOpenMM ** forces)
const RealOpenMM periodicBoxSize[3],
RealOpenMM ** atomParameters,
RealOpenMM ** forces)
{
static const int QIndex = 2; // atom charges are stored in atomParameters[atomID][2]
int i;
int ix,iy,iz;
int x0index,y0index,z0index;
int xindex,yindex,zindex;
int index;
int order;
RealOpenMM q;
RealOpenMM * thetax;
RealOpenMM * thetay;
RealOpenMM * thetaz;
RealOpenMM * dthetax;
RealOpenMM * dthetay;
RealOpenMM * dthetaz;
RealOpenMM tx,ty,tz;
RealOpenMM dtx,dty,dtz;
RealOpenMM fx,fy,fz;
RealOpenMM gridvalue;
int nx,ny,nz;
nx = pme->ngrid[0];
ny = pme->ngrid[1];
nz = pme->ngrid[2];
order = pme->order;
int i;
int ix,iy,iz;
int x0index,y0index,z0index;
int xindex,yindex,zindex;
int index;
int order;
RealOpenMM q;
RealOpenMM * thetax;
RealOpenMM * thetay;
RealOpenMM * thetaz;
RealOpenMM * dthetax;
RealOpenMM * dthetay;
RealOpenMM * dthetaz;
RealOpenMM tx,ty,tz;
RealOpenMM dtx,dty,dtz;
RealOpenMM fx,fy,fz;
RealOpenMM gridvalue;
int nx,ny,nz;
nx = pme->ngrid[0];
ny = pme->ngrid[1];
nz = pme->ngrid[2];
order = pme->order;
/* This is almost identical to the charge spreading routine! */
for(i=0;i<pme->natoms;i++)
{
fx = fy = fz = 0;
for(i=0;i<pme->natoms;i++)
{
fx = fy = fz = 0;
q = atomParameters[i][QIndex];
q = atomParameters[i][QIndex];
/* Grid index for the actual atom position */
x0index = pme->particleindex[i][0];
y0index = pme->particleindex[i][1];
z0index = pme->particleindex[i][2];
x0index = pme->particleindex[i][0];
y0index = pme->particleindex[i][1];
z0index = pme->particleindex[i][2];
/* Bspline factors for this atom in each dimension , calculated from fractional coordinates */
thetax = &(pme->bsplines_theta[0][i*order]);
thetay = &(pme->bsplines_theta[1][i*order]);
thetaz = &(pme->bsplines_theta[2][i*order]);
dthetax = &(pme->bsplines_dtheta[0][i*order]);
dthetay = &(pme->bsplines_dtheta[1][i*order]);
dthetaz = &(pme->bsplines_dtheta[2][i*order]);
thetax = &(pme->bsplines_theta[0][i*order]);
thetay = &(pme->bsplines_theta[1][i*order]);
thetaz = &(pme->bsplines_theta[2][i*order]);
dthetax = &(pme->bsplines_dtheta[0][i*order]);
dthetay = &(pme->bsplines_dtheta[1][i*order]);
dthetaz = &(pme->bsplines_dtheta[2][i*order]);
/* See pme_grid_spread_charge() for comments about the order here, and only interpolation in one direction */
/* See pme_grid_spread_charge() for comments about the order here, and only interpolation in one direction */
/* Since we will add order^3 (typically 4*4*4=64) terms to the force on each particle, we use temporary fx/fy/fz
* variables, and only add it to memory forces[] at the end.
*/
for(ix=0;ix<order;ix++)
{
xindex = (x0index + ix) % pme->ngrid[0];
for(ix=0;ix<order;ix++)
{
xindex = (x0index + ix) % pme->ngrid[0];
/* Get both the bspline factor and its derivative with respect to the x coordinate! */
tx = thetax[ix];
dtx = dthetax[ix];
tx = thetax[ix];
dtx = dthetax[ix];
for(iy=0;iy<order;iy++)
{
yindex = (y0index + iy) % pme->ngrid[1];
for(iy=0;iy<order;iy++)
{
yindex = (y0index + iy) % pme->ngrid[1];
/* bspline + derivative wrt y */
ty = thetay[iy];
dty = dthetay[iy];
ty = thetay[iy];
dty = dthetay[iy];
for(iz=0;iz<order;iz++)
{
/* Can be optimized, but we keep it simple here */
zindex = (z0index + iz) % pme->ngrid[2];
for(iz=0;iz<order;iz++)
{
/* Can be optimized, but we keep it simple here */
zindex = (z0index + iz) % pme->ngrid[2];
/* bspline + derivative wrt z */
tz = thetaz[iz];
dtz = dthetaz[iz];
index = xindex*pme->ngrid[1]*pme->ngrid[2] + yindex*pme->ngrid[2] + zindex;
tz = thetaz[iz];
dtz = dthetaz[iz];
index = xindex*pme->ngrid[1]*pme->ngrid[2] + yindex*pme->ngrid[2] + zindex;
/* Get the fft+convoluted+ifft:d data from the grid, which must be real by definition */
/* Checking that the imaginary part is indeed zero might be a good check :-) */
gridvalue = pme->grid[index].re;
/* The d component of the force is calculated by taking the derived bspline in dimension d, normal bsplines in the other two */
fx += dtx*ty*tz*gridvalue;
fy += tx*dty*tz*gridvalue;
fz += tx*ty*dtz*gridvalue;
}
}
}
gridvalue = pme->grid[index].re;
/* The d component of the force is calculated by taking the derived bspline in dimension d, normal bsplines in the other two */
fx += dtx*ty*tz*gridvalue;
fy += tx*dty*tz*gridvalue;
fz += tx*ty*dtz*gridvalue;
}
}
}
/* Update memory force, note that we multiply by charge and some box stuff */
forces[i][0] -= q*fx*nx/periodicBoxSize[0];
forces[i][1] -= q*fy*ny/periodicBoxSize[1];
forces[i][2] -= q*fz*nz/periodicBoxSize[2];
}
forces[i][0] -= q*fx*nx/periodicBoxSize[0];
forces[i][1] -= q*fy*ny/periodicBoxSize[1];
forces[i][2] -= q*fz*nz/periodicBoxSize[2];
}
}
......@@ -621,10 +621,10 @@ pme_grid_interpolate_force(pme_t pme,
int
pme_init(pme_t * ppme,
RealOpenMM ewaldcoeff,
int natoms,
const int ngrid[3],
int pme_order,
RealOpenMM ewaldcoeff,
int natoms,
const int ngrid[3],
int pme_order,
RealOpenMM epsilon_r)
{
pme_t pme;
......@@ -635,27 +635,27 @@ pme_init(pme_t * ppme,
pme->order = pme_order;
pme->epsilon_r = epsilon_r;
pme->ewaldcoeff = ewaldcoeff;
pme->natoms = natoms;
pme->natoms = natoms;
for(d=0;d<3;d++)
{
pme->ngrid[d] = ngrid[d];
pme->bsplines_theta[d] = (RealOpenMM *)malloc(sizeof(RealOpenMM)*pme_order*natoms);
pme->bsplines_dtheta[d] = (RealOpenMM *)malloc(sizeof(RealOpenMM)*pme_order*natoms);
}
{
pme->ngrid[d] = ngrid[d];
pme->bsplines_theta[d] = (RealOpenMM *)malloc(sizeof(RealOpenMM)*pme_order*natoms);
pme->bsplines_dtheta[d] = (RealOpenMM *)malloc(sizeof(RealOpenMM)*pme_order*natoms);
}
pme->particlefraction = (rvec *)malloc(sizeof(rvec)*natoms);
pme->particleindex = (ivec *)malloc(sizeof(ivec)*natoms);
pme->particlefraction = (rvec *)malloc(sizeof(rvec)*natoms);
pme->particleindex = (ivec *)malloc(sizeof(ivec)*natoms);
/* Allocate charge grid storage */
pme->grid = (t_complex *)malloc(sizeof(t_complex)*ngrid[0]*ngrid[1]*ngrid[2]);
/* Allocate charge grid storage */
pme->grid = (t_complex *)malloc(sizeof(t_complex)*ngrid[0]*ngrid[1]*ngrid[2]);
fftpack_init_3d(&pme->fftplan,ngrid[0],ngrid[1],ngrid[2]);
fftpack_init_3d(&pme->fftplan,ngrid[0],ngrid[1],ngrid[2]);
/* Setup bspline moduli (see Essman paper) */
pme_calculate_bsplines_moduli(pme);
/* Setup bspline moduli (see Essman paper) */
pme_calculate_bsplines_moduli(pme);
*ppme = pme;
*ppme = pme;
return 0;
}
......@@ -665,12 +665,12 @@ pme_init(pme_t * ppme,
int pme_exec(pme_t pme,
RealOpenMM ** atomCoordinates,
RealOpenMM ** forces,
RealOpenMM ** atomParameters,
const RealOpenMM periodicBoxSize[3],
RealOpenMM * energy,
RealOpenMM pme_virial[3][3])
RealOpenMM ** atomCoordinates,
RealOpenMM ** forces,
RealOpenMM ** atomParameters,
const RealOpenMM periodicBoxSize[3],
RealOpenMM * energy,
RealOpenMM pme_virial[3][3])
{
/* Routine is called with coordinates in x, a box, and charges in q */
......@@ -679,28 +679,28 @@ int pme_exec(pme_t pme,
* and what its fractional offset in this grid cell is.
*/
/* Update charge grid indices and fractional offsets for each atom.
/* Update charge grid indices and fractional offsets for each atom.
* The indices/fractions are stored internally in the pme datatype
*/
pme_update_grid_index_and_fraction(pme,atomCoordinates,periodicBoxSize);
pme_update_grid_index_and_fraction(pme,atomCoordinates,periodicBoxSize);
/* Calculate bsplines (and their differentials) from current fractional coordinates, store in pme structure */
pme_update_bsplines(pme);
/* Calculate bsplines (and their differentials) from current fractional coordinates, store in pme structure */
pme_update_bsplines(pme);
/* Spread the charges on grid (using newly calculated bsplines in the pme structure) */
pme_grid_spread_charge(pme,atomParameters);
/* Spread the charges on grid (using newly calculated bsplines in the pme structure) */
pme_grid_spread_charge(pme,atomParameters);
/* do 3d-fft */
fftpack_exec_3d(pme->fftplan,FFTPACK_FORWARD,pme->grid,pme->grid);
/* do 3d-fft */
fftpack_exec_3d(pme->fftplan,FFTPACK_FORWARD,pme->grid,pme->grid);
/* solve in k-space */
pme_reciprocal_convolution(pme,periodicBoxSize,energy,pme_virial);
/* solve in k-space */
pme_reciprocal_convolution(pme,periodicBoxSize,energy,pme_virial);
/* do 3d-invfft */
fftpack_exec_3d(pme->fftplan,FFTPACK_BACKWARD,pme->grid,pme->grid);
/* do 3d-invfft */
fftpack_exec_3d(pme->fftplan,FFTPACK_BACKWARD,pme->grid,pme->grid);
/* Get the particle forces from the grid and bsplines in the pme structure */
pme_grid_interpolate_force(pme,periodicBoxSize,atomParameters,forces);
/* Get the particle forces from the grid and bsplines in the pme structure */
pme_grid_interpolate_force(pme,periodicBoxSize,atomParameters,forces);
return 0;
}
......@@ -710,24 +710,24 @@ int pme_exec(pme_t pme,
int
pme_destroy(pme_t pme)
{
int d;
int d;
free(pme->grid);
free(pme->grid);
for(d=0;d<3;d++)
{
free(pme->bsplines_moduli[d]);
free(pme->bsplines_theta[d]);
free(pme->bsplines_dtheta[d]);
}
for(d=0;d<3;d++)
{
free(pme->bsplines_moduli[d]);
free(pme->bsplines_theta[d]);
free(pme->bsplines_dtheta[d]);
}
free(pme->particlefraction);
free(pme->particleindex);
free(pme->particlefraction);
free(pme->particleindex);
fftpack_destroy(pme->fftplan);
fftpack_destroy(pme->fftplan);
/* destroy structure itself */
free(pme);
/* destroy structure itself */
free(pme);
return 0;
return 0;
}
platforms/reference/src/SimTKReference/ReferenceCustomBondIxn.cpp
View file @ 6186c067
......@@ -116,7 +116,7 @@ int ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices,
int atomBIndex = atomIndices[1];
ReferenceForce::getDeltaR( atomCoordinates[atomAIndex], atomCoordinates[atomBIndex], deltaR );
variables["r"] = deltaR[ReferenceForce::RIndex];
RealOpenMM dEdR = forceExpression.evaluate(variables);
RealOpenMM dEdR = (RealOpenMM) forceExpression.evaluate(variables);
dEdR = deltaR[ReferenceForce::RIndex] > zero ? (dEdR/deltaR[ReferenceForce::RIndex]) : zero;
forces[atomAIndex][0] += dEdR*deltaR[ReferenceForce::XIndex];
......@@ -127,7 +127,7 @@ int ReferenceCustomBondIxn::calculateBondIxn( int* atomIndices,
forces[atomBIndex][1] -= dEdR*deltaR[ReferenceForce::YIndex];
forces[atomBIndex][2] -= dEdR*deltaR[ReferenceForce::ZIndex];
RealOpenMM energy = energyExpression.evaluate(variables);
RealOpenMM energy = (RealOpenMM) energyExpression.evaluate(variables);
updateEnergy( energy, energiesByBond, twoI, atomIndices, energiesByAtom );
return ReferenceForce::DefaultReturn;
......
platforms/reference/src/SimTKReference/ReferenceCustomExternalIxn.cpp
View file @ 6186c067
......@@ -106,11 +106,11 @@ int ReferenceCustomExternalIxn::calculateForce( int atomIndex,
// ---------------------------------------------------------------------------------------
forces[atomIndex][0] -= forceExpressionX.evaluate(variables);
forces[atomIndex][1] -= forceExpressionY.evaluate(variables);
forces[atomIndex][2] -= forceExpressionZ.evaluate(variables);
forces[atomIndex][0] -= (RealOpenMM) forceExpressionX.evaluate(variables);
forces[atomIndex][1] -= (RealOpenMM) forceExpressionY.evaluate(variables);
forces[atomIndex][2] -= (RealOpenMM) forceExpressionZ.evaluate(variables);
if (energy != NULL)
*energy += energyExpression.evaluate(variables);
*energy += (RealOpenMM) energyExpression.evaluate(variables);
return ReferenceForce::DefaultReturn;
}
platforms/reference/src/SimTKReference/ReferenceLJCoulomb14.cpp
View file @ 6186c067
......@@ -157,9 +157,9 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
RealOpenMM dEdR = parameters[1]*( twelve*sig6 - six )*sig6;
if (cutoff)
dEdR += ONE_4PI_EPS0*parameters[2]*(inverseR-2.0f*krf*r2);
dEdR += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*(inverseR-2.0f*krf*r2));
else
dEdR += ONE_4PI_EPS0*parameters[2]*inverseR;
dEdR += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*inverseR);
dEdR *= inverseR*inverseR;
// accumulate forces
......@@ -172,9 +172,9 @@ int ReferenceLJCoulomb14::calculateBondIxn( int* atomIndices, RealOpenMM** atomC
RealOpenMM energy = parameters[1]*( sig6 - one )*sig6;
if (cutoff)
energy += ONE_4PI_EPS0*parameters[2]*(inverseR+krf*r2-crf);
energy += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*(inverseR+krf*r2-crf));
else
energy += ONE_4PI_EPS0*parameters[2]*inverseR;
energy += (RealOpenMM) (ONE_4PI_EPS0*parameters[2]*inverseR);
// accumulate energies
......
platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp
View file @ 6186c067
......@@ -206,7 +206,7 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
// **************************************************************************************
for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){
RealOpenMM selfEwaldEnergy = ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI;
RealOpenMM selfEwaldEnergy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI);
totalSelfEwaldEnergy -= selfEwaldEnergy;
if( energyByAtom ){
......@@ -373,8 +373,8 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM alphaR = alphaEwald * r;
RealOpenMM dEdR = ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
RealOpenMM dEdR = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR);
dEdR = (RealOpenMM) (dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
RealOpenMM sig = atomParameters[ii][SigIndex] + atomParameters[jj][SigIndex];
RealOpenMM sig2 = inverseR*sig;
......@@ -423,8 +423,8 @@ int ReferenceLJCoulombIxn::calculateEwaldIxn( int numberOfAtoms, RealOpenMM** at
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
RealOpenMM dEdR = ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erf(alphaR) - 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
RealOpenMM dEdR = (RealOpenMM) (ONE_4PI_EPS0 * atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR);
dEdR = (RealOpenMM) (dEdR * (erf(alphaR) - 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
// accumulate forces
......@@ -582,9 +582,9 @@ int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoo
RealOpenMM eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
RealOpenMM dEdR = eps*( twelve*sig6 - six )*sig6;
if (cutoff)
dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2);
dEdR += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR-2.0f*krf*r2));
else
dEdR += ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
dEdR += (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR);
dEdR *= inverseR*inverseR;
// accumulate forces
......@@ -601,9 +601,9 @@ int ReferenceLJCoulombIxn::calculateOneIxn( int ii, int jj, RealOpenMM** atomCoo
if( totalEnergy || energyByAtom ) {
if (cutoff)
energy = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
energy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf));
else
energy = ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR;
energy = (RealOpenMM) (ONE_4PI_EPS0*atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR);
energy += eps*(sig6-one)*sig6;
if( totalEnergy )
*totalEnergy += energy;
......
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