An fpga implementation of the Smooth Particle Mesh Ewald Reciprocal Sum Compute Engine (rsce)
3.2 Program Flow In this section, the program flow to calculate the reciprocal energy and forces are described. The PME program starts with the main() program in pme_test.c and it does the following:
In the following subsections, each of the above steps will be explained in more detail. 3.2.1 Memory Allocation The memory allocation is done in the calc_recip_sum() function in recip_sum2.c: The following data arrays are allocated for reciprocal sum calculation:
The code that performs the memory allocation: In the calc_recip_sum() function of recip_sum2.c… dtheta1=dvector(0,numatoms*order); /*used in spline interpolation */ dtheta2=dvector(0,numatoms*order); dtheta3=dvector(0,numatoms*order); theta1=dvector(0,numatoms*order); theta2=dvector(0,numatoms*order); theta3=dvector(0,numatoms*order); bsp_mod1=dvector(0,nfft); bsp_mod2=dvector(0,nfft); bsp_mod3=dvector(0,nfft); : fr1=dvector(0,numatoms); /* fractional coordinates */ fr2=dvector(0,numatoms); fr3=dvector(0,numatoms); : q = dvector(0,siz_q); 3.2.2 Computation of Modulus of the IDFT of B-Spline Coef. The modulus of the IDFT of the B-Spline Coefficients represent the inverse of the B(m1, m2, m3) array in equation (4.8) of the SPME paper. As shown in equation 4.7 in [2], the B(m1, m2, m3) array is necessary in the calculation of the reciprocal energy: Since the computation of the modulus of the IDFT of B-Spline coefficients is not related to the locations of the charges, this computation can be pre-computed at the beginning of the simulation. This step is divided into two sub-steps:
Please refer to the appendix C for the SPME paper [2] for more information on the Cardinal Spline interpolation and the derivation of bi(mi). All the related source codes are attached here for reference:
{ int i_1; int nmax; extern int fill_bspline(double *, int *, double *, double *); int i;
double w, array[MAXORD]; extern int dftmod(double *, double *, int *) ; double darray[MAXORD], bsp_arr[MAXN]; /* Parameter adjustments */ --bsp_mod1; --bsp_mod2; --bsp_mod3; /* this routine loads the moduli of the inverse DFT of the B splines */ /* bsp_mod1-3 hold these values, nfft1-3 are the grid dimensions, */ /* Order is the order of the B spline approx. */ if (*order > MAXORD) { printf( "Error:order too large! check on MAXORD(pmesh_kspace.c) \n"); exit(2);
} /* Computing MAX */ i_1 = max(*nfft2,*nfft1); nmax = max(*nfft3,i_1); if (nmax > MAXN) { printf("Error: nfft1-3 too large! check on MAXN(pmesh_kspace.c)\n"); exit(3);
} w = 0.;
fill_bspline(&w, order, array, darray); // Mn(k) for (i = 1; i <=nmax; ++i) { bsp_arr[i - 1] = 0.; } i_1 = *order + 1; for (i = 2; i <= i_1; ++i) { bsp_arr[i - 1] = array[i - 2]; //only the first “order” of bsp_arrs //contains non-zero values }
dftmod(&bsp_mod2[1], bsp_arr, nfft2); dftmod(&bsp_mod3[1], bsp_arr, nfft3); return 0; } /* load_bsp_moduli */ In pmesh_kspace.c… int fill_bspline(double *w, int *order, double *array, double *darray) { extern int diff(double *, double *, int *), init(double *, double *, int *), one_pass(double *, double *, int *); static int k;
/* do linear case */ init(&array[1], w, order); /* compute standard b-spline recursion */ for (k = 3; k <= ( *order - 1); ++k) { one_pass(&array[1], w, &k); }
diff(&array[1], &darray[1], order); /* one more recursion */ one_pass(&array[1], w, order); return 0; } /* fill_bspline */
{ --c; c[*order] = 0.; c[2] = *x; c[1] = 1. - *x; return 0; } /* init_ */
{ static int j; static double div; --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]; return 0; } /* one_pass */ int diff(double *c, double *d, int *order) { static int j; --c; --d;
for (j = 2; j <= (*order); ++j) { d[j] = c[j - 1] - c[j]; } return 0; } /* diff_ */ In dftmod() function in pmesh_kspace.c… int dftmod(double *bsp_mod, double *bsp_arr, int *nfft) { double cos(double), sin(double); double tiny; int j, k; double twopi, arg, sum1, sum2;
--bsp_mod; --bsp_arr;
/* storing it into bsp_mod */ twopi = 6.2831853071795862; tiny = 1.e-7; for (k = 1; k <=(*nfft); ++k) { sum1 = 0.; sum2 = 0.;
arg = twopi * (k - 1) * (j - 1) / *nfft; sum1 += bsp_arr[j] * cos(arg); sum2 += bsp_arr[j] * sin(arg); /* L250: */ } bsp_mod[k] = sum1*sum1 + sum2*sum2; } for (k = 1; k <=(*nfft); ++k) { if (bsp_mod[k] < tiny) { bsp_mod[k] = (bsp_mod[k - 1] + bsp_mod[k + 1]) * .5; } } return 0; } /* dftmod_ */ Program output for calculating the B-Spline coefficients of the grid point (Interpolation order is 4 and grid size is 64): Inside load_bsp_moduli() function in pmesh_kspace.c Before fill_bspline: array[0]=0.000000 Before fill_bspline: array[1]=0.000000 Before fill_bspline: array[2]=0.000000 Before fill_bspline: array[3]=0.000000 Before fill_bspline: array[4]=0.000000 Before fill_bspline: array[5]=0.000000 order = 4 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 After fill_bspline: array[0]=0.166667 After fill_bspline: array[1]=0.666667 After fill_bspline: array[2]=0.166667 After fill_bspline: array[3]=0.000000 After fill_bspline: array[4]=0.000000 After fill_bspline: array[5]=0.000000 bsp_ar[0]=0.000000 bsp_ar[1]=0.166667 bsp_ar[2]=0.666667 bsp_ar[3]=0.166667 bsp_ar[4]=0.000000 : :
bsp_ar[62]=0.000000 bsp_ar[63]=0.000000 Program output for calculating the modulus of IDFT of the B-Spline coefficients: bsp_ar[0]=0.000000 bsp_mod1[0]=0.000000 bsp_mod2[0]=0.000000 bsp_mod3[0]=0.000000 bsp_ar[1]=0.166667 bsp_mod1[1]=1.000000 bsp_mod2[1]=1.000000 bsp_mod3[1]=1.000000 bsp_ar[2]=0.666667 bsp_mod1[2]=0.996792 bsp_mod2[2]=0.996792 bsp_mod3[2]=0.996792 bsp_ar[3]=0.166667 bsp_mod1[3]=0.987231 : bsp_mod1[61]=0.949897 bsp_mod2[61]=0.949897 bsp_mod3[61]=0.949897 : bsp_mod1[64]=0.996792 bsp_mod2[64]=0.996792 bsp_mod3[64]=0.996792 3.2.3 Obtaining the Initial Coordinates The loading of the initial x, y, and z coordinates is done in the main() function in pme_test.c. The (x, y, z) co-ordinates of all particles are obtained from the input pdb file. The value cgo is the charge of Oxygen and cgh is the charge of Hydrogen; their charge magnitudes are hard coded to constant numbers. The value cgh is ½ of cgo, therefore, the simulation environment is neutral. Also observed the way it derives the magnitude of the charge. In main() function in pme_test.c… /* used for water system, hydrogen and oxygen charges */ factor = sqrt(332.17752); cgh = factor * .417; cgo = cgh * -2.; In main() function in pme_test.c… numwats = numatoms / 3; n = 0;
for (i = 1; i <= numwats;i++) { fscanf(infile,"%lf %lf %lf",&ParticlePtr[n].x,&ParticlePtr[n].y,&ParticlePtr[n].z); fscanf(infile,"%lf %lf %lf",&ParticlePtr[n+1].x,&ParticlePtr[n+1].y,&ParticlePtr[n+1].z); fscanf(infile,"%lf %lf %lf",&ParticlePtr[n+2].x,&ParticlePtr[n+2].y,&ParticlePtr[n+2].z); ParticlePtr[n].cg = cgo; ParticlePtr[n + 1].cg = cgh; ParticlePtr[n + 2].cg = cgh; n += 3;
} 3.2.4. Construction of the Reciprocal Lattice Vectors The construction of reciprocal lattice vectors is done in the main() function in pme_test.c. By observing the code, only array elements recip[0], recip[4], and recip[8] have a value of 1/box; others have a value of zero. This implies that only orthogonal simulation box is supported. In main() function in pme_test.c… volume = box * box * box; reclng[0] = box; reclng[1] = box; reclng[2] = box; for (i = 1; i <= 3; i++) { for (j = 1; j <= 3; ++j) { recip[i + j * 3 - 4] = 0.; } recip[i + i * 3 - 4] = 1. / box; } Sample output (array index is adjusted to be zero-based): recip[0] = 0.016445 recip[1] = 0.000000 recip[2] = 0.000000 recip[3] = 0.000000 recip[4] = 0.016445 recip[5] = 0.000000 recip[6] = 0.000000 recip[7] = 0.000000 recip[8] = 0.016445 Directory: ~pc ~pc -> The Tablet War: Apple v s The Rest ~pc -> From: object-oriented analysis and design, Grady Booch, Addison-Wesley, 1998 ~pc -> Analysis of an Industry Price War: The Tablet price war ~pc -> Biography of Pok Chi Lau Home address: 2600 Download 1.53 Mb. Share with your friends: |