An fpga implementation of the Smooth Particle Mesh Ewald Reciprocal Sum Compute Engine (rsce)

The scaled and shifted fractional coordinates are computed in the get_scaled_fractionals() function in pmesh_kspace.c (called by the do_pmesh_kspace() in pmesh_kspace.c). The resulted x, y and z fractional coordinates for all particles are stored in the arrays: fr1[1..numatoms], fr2[1..numatoms], and fr3[1..numatoms] respectively. The array indexing for recip[] is adjusted by -4, so recip[4] is recip[0] and so on. As mentioned in section 3.2.4, since the simulation box is orthogonal, the x fractional coordinates fr1[i] only contains the corresponding x component of the Cartesian coordinates. After the fractional coordinates are derived, they are scaled and shifted. The shifting happens in the software implementation conforms the statement in [2]: “with m defined by m1’a1* + m2’a2* + m3’a3* , where mi’= mi for 0<mi<K/2 and mi’ = mi - Ki otherwise”. Examples of calculating the shifted and scaled coordinates are shown in Table 2 and Table 3. Furthermore, graphical representation of the scale and shift operation is shown in Figure 2.

w = ParticlePtr[n].x*recip[4]+ParticlePtr[n].y*recip[5]+ParticlePtr[n].z*recip[6];

fr1[n+1] = *nfft1 * (w - Nint(w) + .5);

w = ParticlePtr[n].x*recip[7]+ParticlePtr[n].y*recip[8]+ParticlePtr[n].z*recip[9];

fr2[n+1] = *nfft2 * (w - Nint(w) + .5);

w = ParticlePtr[n].x*recip[10]+ParticlePtr[n].y*recip[11]+ParticlePtr[n].z*recip[12];

fr3[n+1] = *nfft3 * (w - Nint(w) + .5);

}
In the pme.h definition, the Nint() function is defined as…

#define Nint(dumx) ( ((dumx) >=0.0) ? (int)((dumx) + 0.5) : (int)((dumx) - 0.5) )

Program output of the scaled and shifted fractional coordinates calculation (the variable “term” is calculated as (w-Nint(w)+0.5):

1: x=0.000000, w=0.000000, Nint(w)=0, term=0.500000, fr1=32.000000

2: x=60.810000, w=1.000000, Nint(w)=1, term=0.500000, fr1=32.000000

3: x=30.405000, w=0.500000, Nint(w)=1, term=0.000000, fr1=0.000000

4: x=3.921000, w=0.064480, Nint(w)=0, term=0.564480, fr1=36.126690

5: x=3.145000, w=0.051718, Nint(w)=0, term=0.551718, fr1=35.309982

6: x=4.572000, w=0.075185, Nint(w)=0, term=0.575185, fr1=36.811840

7: x=7.941000, w=0.130587, Nint(w)=0, term=0.630587, fr1=40.357573

8: x=7.480000, w=0.123006, Nint(w)=0, term=0.623006, fr1=39.872389

9: x=8.754000, w=0.143957, Nint(w)=0, term=0.643957, fr1=41.213222

10: x=4.168000, w=0.068541, Nint(w)=0, term=0.568541, fr1=36.386647

The computation of the B-Spline coefficients and the corresponding derivatives for all particles are initiated in the get_bspline_coeffs() function in bspline.c (called by the do_pmesh_kspace() function in pmesh_kspace.c). The actual calculations are done in fill_bspline() function in the source file bspline.c. The fractional part of the scaled and shifted fractional coordinates of a particle (w = fr1[n] - (int) fr1[n]) is input to the fill_bspline() function and the function returns with the B-Spline coefficients and derivates in two arrays of size “order” (the interpolation order). The B-Spline coefficients can be seen as influence/weights of a particle at the nearby grid points.

For n greater than 2:

The derivate is defined as:

{

int theta1_dim1, theta1_offset, theta2_dim1, theta2_offset,

dtheta2_dim1, dtheta2_offset, dtheta3_dim1, dtheta3_offset;

double *, double *);

static int n;

static double w;

--fr1; --fr2; --fr3;

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;

/* INPUT: */

/* numatoms: number of atoms */

/* fr1,fr2,fr3 the scaled and shifted fractional coords */

/* order: the order of spline interpolation */

/* OUTPUT */

/* theta1,theta2,theta3: the spline coeff arrays */

/* dtheta1,dtheta2,dtheta3: the 1st deriv of spline coeff arrays */

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

w = fr1[n] - ( int) fr1[n];

w = fr2[n] - ( int) fr2[n];

w = fr3[n] - ( int) fr3[n];

/* L100: */

}

return 0;

{

extern int diff(double *, double *, int *),

double *, double *, int *);

static int k;
/* Parameter adjustments */

--array;

/* compute standard b-spline recursion */

for (k = 3; k <= ( *order - 1); ++k) {

one_pass(&array[1], w, &k);

/* L10: */

}

/* perform standard b-spline differentiation */

} /* fill_bspline */

As an example, assume the scaled and shifted fractional coordinates (x, y, z) of a particle = (15.369, 8.850, 39.127), and the interpolation order=4. For the x dimension, (the distance to the grid point) w = 15.369 - 15 = 0.369.

{

--c;

c[2] = *x;

c[1] = 1. - *x;

return 0;

} /* init_ */
 init(&array[1], w, order);

c[4] = 0, c[2] = w = 0.369, c[1] = 1 - w = 0.631

{

static int j;

--c;

c[*k] = div * *x * c[*k - 1];

c[*k - j] = div * ((*x + j) * c[*k - j - 1] + (*k - j - *x) * c[*k - j]);

}

c[1] = div * (1 - *x) * c[1];

} /* one_pass */

one_pass(&array[1], w, &k);

}
k = 3 to order-1 (i.e. 3 to 3)

k = 3

c[3] = div * x * c[k-1] = 1/(3-1) * w * c[2] = 1/(3-1)*w*w = 0.0681

j = 1

c[3-1] = div * ((x+1) * c[3-1-1] + (3-1-x) * c [3-1])

c[2] = div * ((w+1) * c[1] + (2-w) * c[2])

c[2] = ( 1/(3-1) ) * ((w+1) * c[1] + (2-w) * c[2])

c[2] = ( 1/(3-1) ) * ((w+1) * (1-w) + (2-w) * w) = 0.733

c[1] = div * (1 – x) * c[1]

c[1] = div * (1 – w) * c[1]

c[1] = ( 1/(3-1) ) * (1 – w) * (1 – w) = 0.199

Differentiation

k = 4

c[4] = div * x * c[k-1] = 1/(4-1) * w * c[3]

j=1

c[4-1] = div * ((x+1) * c[4-1-1] + (4-1-x) * c [4-1])

c[3] = div * ((w+1) * c[2] + (3-w) * c[3])

j=2

c[4-2] = div * ((x+2) * c[4-2-1] + (4-2-x) * c [4-2])

c[2] = div * ((w+2) * c[1] + (2-w) * c[2])
c[1] = div * (1 – x) * c[1]

c[1] = div * (1 – w) * c[1]

c[1] = ( 1/(4-1) ) * (1 – w) * c[1]
c[3] = ( 1/(3-1) ) * w * c[2] = ( 1/(3-1) ) * w * w -- from previous pass

c[4] = ( 1/(4-1) ) * w * c[3] = ( 1/(4-1) ) * w * ( 1/(3-1) ) * w * w = 0.00837

c[2] = ( 1/2 ) * ((w+1) * (1-w) + (2-w) * w) -- from previous pass

c[3] = div * ((w+1) * c[2] + (3-w) * c[3])

c[3] = (1/3)*((w+1)*((1/2)*((w+1)*(1-w)+(2-w)*w)) + (3-w)*((1/2)*w*w))

c[3] = (1/3)*(1.369*(0.5*(1.369*0.631+1.631*0.369))+2.631*(0.5*0.369*0.369))

c[3] = (1/3)*(1.369*0.732839 + 2.631*0.068081 ) = 0.394

c[2] = div * ((w+2)*c[1] + (2-w)*c[2] )

c[2] = (1/3)*((w+2)*((1/2)*(1–w)*(1–w))+(2-w)*((1/2)*((w+1)*(1-w) + (2-w)*w)))

c[2] = (1/3)*(2.369*(0.5*0.631*0.631) + 1.631*(0.5*(1.369*0.631 + 1.631*0.369))

c[2] = (1/3)*(2.369*0.199081 + 1.631*0.732839) = 0.555
c[1] = div * (1 – w) * c[1]

c[1] = (1/3) * (1 – w) * ( (1/2) * (1 – w) * (1 – w) )

c[1] = (1/3) * (0.631) * ( 0.5 * 0.631 * 0.631 ) = 0.0419
Note: C[4] + c[3] + c[2] + c[1] = 1.0

Program output of the B-Spline calculation: (Please note that the sum of the derivates is zero and the sum of the coefficients is one).

After init: array[0]=0.000000

After init: array[1]=1.000000

After init: array[2]=0.000000

After init: array[3]=0.000000

After recursions: array[0]=0.000000

After recursions: array[1]=0.500000

After recursions: array[2]=0.500000

After recursions: array[3]=0.000000

Last recursion: array[0]=0.000000

Last recursion: array[1]=0.166667

Last recursion: array[2]=0.666667

Last recursion: array[3]=0.166667

fr1[1]=32.000000, w=0.000000, order=4

After init: array[0]=0.000000

After init: array[1]=0.691959

After init: array[2]=0.308041

After init: array[3]=0.000000

After recursions: array[0]=0.000000

After recursions: array[1]=0.239403

After recursions: array[2]=0.713152

After recursions: array[3]=0.047445

Last recursion: array[0]=0.000000

Last recursion: array[1]=0.055219

Last recursion: array[2]=0.586392

Last recursion: array[3]=0.353517

fr2[1]=34.308041, w=0.308041, order=4

theta2[5-8]=0.055219, 0.586392, 0.353517, 0.004872, sum=1.000000

dtheta2[5-8]=-0.239403, -0.473749, 0.665707, 0.047445, sum=0.000000

After init: array[0]=0.000000

After init: array[1]=0.573919

After init: array[2]=0.426081

After init: array[3]=0.000000

After recursions: array[0]=0.000000

After recursions: array[1]=0.164691

After recursions: array[2]=0.744536

After recursions: array[3]=0.090773

Last recursion: array[0]=0.000000

Last recursion: array[1]=0.031506

Last recursion: array[2]=0.523798

Last recursion: array[3]=0.431803

fr3[1]=33.426081, w=0.426081, order=4

theta3[5-8]=0.031506, 0.523798, 0.431803, 0.012892, sum=1.000000

dtheta3[5-8]=-0.164691, -0.579845, 0.653763, 0.090773, sum=-0.000000

:

After recursions: array[3]=0.345375

Last recursion: array[1]=0.000803

Last recursion: array[2]=0.262963

Last recursion: array[3]=0.640553

fr2[3]=34.831113, w=0.831113, order=4

theta2[13-16]=0.000803, 0.262963, 0.640553, 0.095682, sum=1.000000

dtheta2[13-16]=-0.014261, -0.626103, 0.294989, 0.345375, sum=0.000000

3.2.7 Composition of the Grid Charge Array

The grid charge array is composed in the fill_charge_grid() function in charge_grid.c (called by the do_pmesh_space() function in pmesh_space.c). During the simulation, the charge grid composing loop goes through all charged particle. For each charge, the counters k1, k2, and k3 are responsible to go through all the grid points that the charges have to interpolate to in the three dimensional space. In a 3D simulation system, each charge will interpolate to P*P*P grid points (P is the interpolation order). In the PME implementation, the index of charge array q starts at 1; that is, it is not 0-based. Also, there is a Nsign(i, j) function which is responsible for handling the wrap-around condition, for example, a point on the edge of the mesh (shifted and scaled reciprocal coordinate = 0.0) will be interpolated to grids located at 1, 62, 63, and 64, etc. An example is shown below:

{

int theta1_dim1, theta1_offset, theta2_dim1, theta2_offset,

static int ntot, i, j, k, n, i0, j0, k0;

extern /* Subroutine */ int clearq(double *, int *);

static 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;

--fr1; --fr2; --fr3;

q_dim2 = *nfftdim1;

q_dim3 = *nfftdim2;

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

q -= q_offset;

/* --------------------------------------------------------------------- */

/* INPUT: */

/* numatoms: number of atoms */

/* charge: the array of atomic charges */

/* theta1,theta2,theta3: the spline coeff arrays */

/* fr1,fr2,fr3 the scaled and shifted fractional coords */

/* nfft1,nfft2,nfft3: the charge grid dimensions */

/* nfftdim1,nfftdim2,nfftdim3: physical charge grid dims */

/* order: the order of spline interpolation */

/* OUTPUT: */

/* Q the charge grid */

/* --------------------------------------------------------------------- */

ntot = (*nfftdim1 * 2) * *nfftdim2 * *nfftdim3;
clearq(&q[q_offset], &ntot);
for (n = 1; n <= (*numatoms); ++n) {

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;

}

}

}

return 0;

/* below is to be used in charge_grid.c */

#define Nsign(i,j) ((j) >=0.0 ? Nabs((i)): -Nabs((i)) )
Based on the sign of variable j, this function returns a number that has the magnitude of variable i and the sign of variable j.
Program output shows the wrap around condition:

n=3, fr1=0.000000, fr2=34.831113, fr3=32.683046, q[62][32][30] = 0.000005

n=3, fr1=0.000000, fr2=34.831113, fr3=32.683046, q[63][32][30] = 0.000022

n=3, fr1=0.000000, fr2=34.831113, fr3=32.683046, q[64][32][30] = 0.000005

n=3, fr1=0.000000, fr2=34.831113, fr3=32.683046, q[1][32][30] = 0.000000
3.2.7.1 How fill_charge_grid() Function Works

Again, as an example, assume the scaled and shifted fractional coordinate (x, y, z) of a particle equals to (15.369, 8.850, 39.127) and the interpolation order=4. Therefore, fr1[n]=15.369, fr2[n]=8.850, and fr3[n]=39.127.

n=1

k0 = (int) fr3[] - order = 39 - 4 = 35

ith3 = 1

k = k0 + 1 + (nff3 – Nsign (nfft3, k0)/2 = 36+1+(64-64)/2 = 37

j0 = (int) fr2[] - order = 8 – 4 = 4

for (ith2=1, ith2<=4)

ith2=1

j0 = 5

prod = theta2[1+1*order] * theta3[1+1*order]*charge

i0 = 15 – 4 = 11

for (ith1=1; ith1<=4)

ith1=1

i0 = 12

q[13+(6+37*q_dim3)*q_dim2*2+1] = q[] + theta1[1+1*order]*prod

ith1=2

i0 = 13

q[13+(6+37*q_dim3)*q_dim2*2+1] = q[] + theta1[1+1*order]*prod

ith1=3

i0 = 14

q[15+(6+37*q_dim3)*q_dim2*2+1] = q[] + theta1[1+1*order]*prod

ith1=4

i0 = 15

q[16+(6+37*q_dim3)*q_dim2*2+1] = q[] + theta1[1+1*order]*prod

:
By observing the iteration, we know that ith3 goes from 1 to 4, ith2 goes from 1 to 4, ith1 goes from 1 to 4, k goes from 37 to 40, j goes from 6 to 9 and i goes from 13 to 16. The grid contribution result for the particle at (15.369, 8.850, 39.127) are shown in Table 4; there are 64 (4*4*4) grid points that the charges are interpolating to. The indexing of q[(i+(j+k*q_dim3)*q_dim2<<1)+1] array can be viewed as a 3-D mesh, as shown in Figure 3.