int fft_back(double *array, double *fftable, double *ffwork, int *nfft1

{

int isign;

int *, double *, int *, int *, double *,

int *, double *, int *);
--array; --fftable; --ffwork;
isign = -1;

pubz3d(&isign, nfft1, nfft2, nfft3, &array[1], nfftdim1, nfftdim2, &

fftable[1], nfftable, &ffwork[1], nffwork);

return 0;

} /* fft_back */

3.2.9 Computation of the Reciprocal Energy, EER

The reciprocal energy is calculated in the scalar_sum() function in charge_grid.c (called by do_pmesh_kspace() in pmesh_kspace.c). The Ewald reciprocal energy is calculated according to the equation 4.7 [2]:

{

int q_dim2, q_dim3, q_offset;

int k;

int k1;

int k2, k3, m1, m2, m3;

double struc2, pi, energy;

int indtop, nf1, nf2, nf3;

double fac;

int nff, ind, jnd;

double msq;
q_dim2 = *nfftdim1;

q_dim3 = *nfftdim2;

q_offset = (q_dim2 * (q_dim3 + 1) + 1 << 1) + 1;

q -= q_offset;

recip -= 4;

--bsp_mod1;

--bsp_mod2;

--bsp_mod3;

--vir;
indtop = *nfft1 * *nfft2 * *nfft3;

pi = 3.14159265358979323846;

fac = pi*pi / (*ewaldcof * *ewaldcof);

nff = *nfft1 * *nfft2;

nf1 = *nfft1 / 2;

if (nf1 << 1 < *nfft1) {

++nf1;

}

nf2 = *nfft2 / 2;

++nf2;

nf3 = *nfft3 / 2;

if (nf3 << 1 < *nfft3) {

++nf3;

for (k = 1; k <= 6; ++k)

vir[k] = 0.;

#endif
for (ind = 1; ind <= (indtop - 1); ++ind) {

k3 = ind / nff + 1;

jnd = ind - (k3 - 1) * nff;

k2 = jnd / *nfft1 + 1;

k1 = jnd - (k2 - 1) * *nfft1 + 1;

m1 = k1 - 1;

if (k1 > nf1) {

m1 = k1 - 1 - *nfft1;

}

m2 = k2 - 1;

m2 = k2 - 1 - *nfft2;

}

m3 = k3 - 1;

m3 = k3 - 1 - *nfft3;

}

mhat1 = recip[4] * m1 + recip[7] * m2 + recip[10] * m3;

mhat3 = recip[6] * m1 + recip[9] * m2 + recip[12] * m3;

msq = mhat1 * mhat1 + mhat2 * mhat2 + mhat3 * mhat3;

denom = pi * *volume * bsp_mod1[k1] * bsp_mod2[k2] * bsp_mod3[k3] * msq;

vterm = (fac * msq + 1.) * 2. / msq;

energy += eterm * struc2;

vir[1] += eterm * struc2 * (vterm * mhat1 * mhat1 - 1.);

vir[2] += eterm * struc2 * (vterm * mhat1 * mhat2);

vir[3] += eterm * struc2 * (vterm * mhat1 * mhat3);

vir[4] += eterm * struc2 * (vterm * mhat2 * mhat2 - 1.);

vir[5] += eterm * struc2 * (vterm * mhat2 * mhat3);

vir[6] += eterm * struc2 * (vterm * mhat3 * mhat3 - 1.);

#endif

eterm * q[(k1 + (k2 + k3 * q_dim3) * q_dim2 << 1) + 1];

eterm * q[(k1 + (k2 + k3 * q_dim3) * q_dim2 << 1) + 2];

}

*eer = energy * .5;

for (k = 1; k <= 6; ++k)

vir[k] *= .5;

#endif
return 0;

} /* scalar_sum */

3.2.10 Computation of F(Q) using 3D-FFT and Q Array Update

Similar to the IFFT operation described in the section 3.2.8, the 3D FFT is done in fft_forward() function in the fftcall.c. It invokes the dynamic library libpubfft.a to perform the 3D-FFT operation. The transformed elements are stored in the original charge array; hence, this is termed in-place FFT operation. The 3D FFT is performed on the updated Q charge array (in step 3.2.9) to obtain the convolution (θ_rec * Q) which is necessary to calculate the reciprocal forces (as shown in equation 4.9).

The following calculation shows the derivation of the convolution result. As described in the section 3.2.9, the updated Q array contains the value of B∙C∙F^-1(Q).

(θ_rec * Q)

= F[F^-1(θ_rec)∙F^-1(Q)]

= F[B∙C∙F^-1(Q)]
In the fft_forward() function in fftcalls.c…
int fft_forward(double *array, double *fftable, double *ffwork, int *nfft1, int *nfft2, int *nfft3, int *nfftdim1, int *nfftdim2, int *nfftdim3, int *nfftable, int *nffwork)

{

int isign;

int *, double *, int *, int *, double *,

int *, double *, int *);
--array;

--fftable;

--ffwork;
isign = 1;

pubz3d(&isign, nfft1, nfft2, nfft3, &array[1], nfftdim1, nfftdim2, &

fftable[1], nfftable, &ffwork[1], nffwork);

return 0;

The reciprocal forces for all particles are calculated in the grad_sum() function in charge_grid.c (called by do_pmesh_kspace() in pmesh_kspace.c).

This function goes through all particles. For each particle, it identifies all the grid points that this particle has been interpolated to. Then, it calculates the contribution of the reciprocal force exerted on this charged particle by each of these grid points. For a 3-D simulation space, there are three components of the forces, namely, the x, y and z directional forces. The convolution (θ_rec * Q) is already resided in the Q charge array which is calculated in the operation described in the section 3.2.10. The derivates of the B-Spline coefficients (dtheta array) which are needed in the calculation of the ∂Q/∂r_αi are already computed in the program step described in the section 3.2.6.

For x direction: ∂Q/∂r_αi = q*dthetax*thetay*thetaz.

For y direction, ∂Q/∂r_αi = q*thetax*dthetay*thetaz.

For z direction, ∂Q/∂r_αi = q*thetax*thetay*dthetaz.
In the grad_sum( ) function in charge_grid.c…
int grad_sum( int *numatoms, PmeParticlePtr ParticlePtr,

double *recip, double *theta1, double *theta2, double

*theta3, double *dtheta1, double *dtheta2, double *

dtheta3, PmeVectorPtr rfparticle, double *

fr1, double *fr2, double *fr3, int *order, int *nfft1,

int *nfft2, int *nfft3, int *nfftdim1, int *nfftdim2,

int *nfftdim3, double *q)

{

int theta1_dim1, theta1_offset, theta2_dim1, theta2_offset;

int dtheta2_dim1, dtheta2_offset, dtheta3_dim1, dtheta3_offset;

int q_dim2, q_dim3, q_offset;
double term; int i, j, k, n; double f1, f2;

int i0, j0, k0; double f3; int ith1, ith2, ith3;

theta1_dim1 = *order;

theta1_offset = theta1_dim1 + 1;

theta1 -= theta1_offset;

theta2_dim1 = *order;

theta2_offset = theta2_dim1 + 1;

theta2 -= theta2_offset;

theta3_dim1 = *order;

theta3_offset = theta3_dim1 + 1;

theta3 -= theta3_offset;

dtheta1_dim1 = *order;

dtheta1_offset = dtheta1_dim1 + 1;

dtheta1 -= dtheta1_offset;

dtheta2_dim1 = *order;

dtheta2_offset = dtheta2_dim1 + 1;

dtheta2 -= dtheta2_offset;

dtheta3_dim1 = *order;

dtheta3_offset = dtheta3_dim1 + 1;

dtheta3 -= dtheta3_offset;

--fr1;

--fr3;
q_dim2 = *nfftdim1;

q_dim3 = *nfftdim2;

q_offset = (q_dim2 * (q_dim3 + 1) + 1 << 1) + 1;

q -= q_offset;

for (n = 1; n <= (*numatoms); ++n) {

f1 = 0.;

f2 = 0.;

k0 = ( int) fr3[n] - *order;

for (ith3 = 1; ith3 <= ( *order); ++ith3) {

++k0;

j0 = ( int) fr2[n] - *order;

for (ith2 = 1; ith2 <=(*order); ++ith2) {

++j0;

i0 = ( int) fr1[n] - *order;

for (ith1 = 1; ith1 <= (*order); ++ith1) {

++i0;

term = ParticlePtr[n-1].cg * q[(i + (j + k * q_dim3) * q_dim2 << 1) + 1];

f1 -= *nfft1 * term * dtheta1[ith1 + n * dtheta1_dim1] *

theta2[ith2 + n * theta2_dim1] * theta3[ith3 + n * theta3_dim1];
f2 -= *nfft2 * term * theta1[ith1 + n * theta1_dim1] *

dtheta2[ith2 + n * dtheta2_dim1] * theta3[ith3 + n * theta3_dim1];

theta2[ith2 + n * theta2_dim1] * dtheta3[ith3 + n * dtheta3_dim1];

}

}

}

rfparticle[n-1].y += recip[5] * f1 + recip[8] * f2 + recip[11] * f3;

rfparticle[n-1].z += recip[6] * f1 + recip[9] * f2 + recip[12] * f3;

printf("n-1=%0d, f.x=%0f, f.y=%0f, f.z=%0f\n",

n-1, rfparticle[n-1].x, rfparticle[n-1].y, rfparticle[n-1].z);

}

return 0;

The adjustment for bonded interaction is done is the adjust_recip() function in the utility_ser2.c. After this step, the reciprocal energy and force calculations are complete.

/* interactions between bonded particles, eg in water excluded H1-O H2-O H1-H2 */

/* interactions */
int adjust_recip( int *numwats, PmeParticlePtr particlelist, double *ewaldcof,

double *ene, PmeVectorPtr afparticle, double *vir)

{

int i2, i3;

PmeParticle particle1, particle2;

*ene = 0.;

vir[i] = 0.;

}
for (n = 1; n <= (*numwats); ++n) {

ilo = (n - 1) * 3 ;

particle1= particlelist[ilo];

i2 = ilo + 1;

particle2= particlelist[i2];

get_adj_pair(ilo, i2, particle1, particle2, ewaldcof, ene, afparticle, &vir[1]);

i2 = ilo + 2;

particle2= particlelist[i2];

get_adj_pair(ilo, i2, particle1, particle2,ewaldcof, ene, afparticle, &vir[1]);

i2 = ilo + 1;

particle1= particlelist[i2];

i3 = ilo + 2;

particle2= particlelist[i3];

get_adj_pair(i2, i3, particle1, particle2, ewaldcof, ene, afparticle, &vir[1]);

return 0;

1. 
   T. Darden, D York, L Pedersen. A Smooth Particle Mesh method. Journal of Chemistry Physics 1995.
   
2. 
   Smooth Particle Mesh Ewald Implementation source code written by A. Toukmaji of Duke University. Homepage: http://www.ee.duke.edu/~ayt/