OptPmeRealSpace.C

Go to the documentation of this file.
00001 
00007 #include <string.h>
00008 #include "OptPmeRealSpace.h"
00009 #include <assert.h>
00010 
00011 OptPmeRealSpace::OptPmeRealSpace(PmeGrid grid,int natoms)
00012   : myGrid(grid), N(natoms) {
00013   int order = myGrid.order;
00014   M = new double[3*N*order];
00015   dM = new double[3*N*order];
00016 }
00017 
00018 OptPmeRealSpace::~OptPmeRealSpace() {
00019   delete [] M;
00020   delete [] dM;
00021 }
00022 
00023 
00024 // this should definitely help the compiler
00025 #define order 4
00026 
00027 void OptPmeRealSpace::fill_charges(double **q_arr, PmeParticle p[], int zstart, int zlen) {
00028   
00029   int i, j, k, l;
00030   int stride;
00031   int K1, K2, K3, dim2, dim3;
00032   int32  K3_1;
00033   double *Mi, *dMi;
00034   Mi = M; dMi = dM;
00035   K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3; dim2=myGrid.dim2; dim3=myGrid.dim3;
00036   K3_1 = K3 - 1;
00037   // order = myGrid.order;
00038   stride = 3*order;
00039 
00040   for (i=0; i<N; i++) {
00041     double q;
00042     int32 u1, u2, u2i, u3i;
00043     double fr1 = p[i].x;
00044     double fr2 = p[i].y;
00045     double fr3 = p[i].z;
00046     q = p[i].cg;
00047     u1  = (int)fr1;
00048     u2i = (int)fr2;
00049     u3i = (int)fr3;
00050     fr1 -= u1;
00051     fr2 -= u2i;
00052     fr3 -= u3i;
00053 
00054     // calculate b_spline for order = 4
00055     Mi[0] = ( ( (-1./6.) * fr1 + 0.5 ) * fr1 - 0.5 ) * fr1 + (1./6.);
00056     Mi[1] = ( ( 0.5 * fr1 - 1.0 ) * fr1 ) * fr1 + (2./3.);
00057     Mi[2] = ( ( -0.5 * fr1 + 0.5 ) * fr1 + 0.5 ) * fr1 + (1./6.);
00058     Mi[3] = (1./6.) * fr1 * fr1 * fr1;
00059     dMi[0] = ( -0.5 * fr1 + 1.0 )* fr1 - 0.5;
00060     dMi[1] = ( 1.5 * fr1 - 2.0 ) * fr1;
00061     dMi[2] = ( -1.5 * fr1 + 1.0 ) * fr1 + 0.5;
00062     dMi[3] = 0.5 * fr1 * fr1;
00063     Mi[4] = ( ( (-1./6.) * fr2 + 0.5 ) * fr2 - 0.5 ) * fr2 + (1./6.);
00064     Mi[5] = ( ( 0.5 * fr2 - 1.0 ) * fr2 ) * fr2 + (2./3.);
00065     Mi[6] = ( ( -0.5 * fr2 + 0.5 ) * fr2 + 0.5 ) * fr2 + (1./6.);
00066     Mi[7] = (1./6.) * fr2 * fr2 * fr2;
00067     dMi[4] = ( -0.5 * fr2 + 1.0 )* fr2 - 0.5;
00068     dMi[5] = ( 1.5 * fr2 - 2.0 ) * fr2;
00069     dMi[6] = ( -1.5 * fr2 + 1.0 ) * fr2 + 0.5;
00070     dMi[7] = 0.5 * fr2 * fr2;
00071     Mi[8] = ( ( (-1./6.) * fr3 + 0.5 ) * fr3 - 0.5 ) * fr3 + (1./6.);
00072     Mi[9] = ( ( 0.5 * fr3 - 1.0 ) * fr3 ) * fr3 + (2./3.);
00073     Mi[10] = ( ( -0.5 * fr3 + 0.5 ) * fr3 + 0.5 ) * fr3 + (1./6.);
00074     Mi[11] = (1./6.) * fr3 * fr3 * fr3;
00075     dMi[8] = ( -0.5 * fr3 + 1.0 )* fr3 - 0.5;
00076     dMi[9] = ( 1.5 * fr3 - 2.0 ) * fr3;
00077     dMi[10] = ( -1.5 * fr3 + 1.0 ) * fr3 + 0.5;
00078     dMi[11] = 0.5 * fr3 * fr3;
00079 
00080     u1  -= order;
00081     u2i -= order;
00082     u3i -= order;
00083     u3i += 1;
00084     for (j=0; j<order; j++) {
00085       double m1;
00086       int32 ind1;
00087       m1 = Mi[j]*q;
00088       u1++;
00089       //ind1 = (u1 + (u1 < 0 ? K1 : 0))*dim2;
00090       ind1 = (u1 + ((unsigned)u1 >>31)*K1)*dim2;
00091       u2 = u2i;
00092       for (k=0; k<order; k++) {
00093         double m1m2;
00094         int32 ind2;
00095         m1m2 = m1*Mi[order+k];
00096         u2++;
00097 
00098         //ind2 = ind1 + (u2 + (u2 < 0 ? K2 : 0));       
00099         ind2 = ind1 + (u2 + ((unsigned)u2 >>31)*K2);    
00100         //assert (ind2 < K1 * K2);      
00101         double *qline = q_arr[ind2];
00102         
00103         //if (qline == NULL)
00104         //printf ("qline [%d, %d] = NULL\n", (u1 + (u1 < 0 ? K1 : 0)), (u2 + (u2 < 0 ? K2 : 0)));
00105         
00106         //assert (qline != NULL);
00107         for (l=0; l<order; l++) {
00108           double m3;
00109           int32 ind;
00110           m3 = Mi[2*order + l];
00111           int32 u3 = u3i + l;
00112           ind = u3 - zstart;
00113           //if (ind >= K3) ind -= K3;
00114           ind = ind - ((unsigned)(K3_1 - ind) >>31)*K3;
00115 
00116           //if (ind < 0 || ind >= zlen) printf ("ind >= zlen  (%d,%d,%d)\n", ind, zlen, zstart); 
00117           
00118           //assert (ind < zlen && ind >= 0);
00119           qline[ind] += m1m2*m3; 
00120         }
00121       }
00122     }
00123     Mi += stride;
00124     dMi += stride;
00125   }
00126 }
00127 
00128 void OptPmeRealSpace::compute_forces(const double * const  * q_arr,
00129                                      const PmeParticle       p[], 
00130                                      Vector                  f[],
00131                                      int                     zstart, 
00132                                      int                     zlen,
00133                                      int                     istart,
00134                                      int                     iend) {
00135   
00136   int i, j, k, l, stride;
00137   double f1, f2, f3;
00138   double *Mi, *dMi;
00139   int K1, K2, K3, dim2;
00140   int32 K3_1;
00141 
00142   if (iend == 0)
00143     iend = N;
00144 
00145   K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3; dim2=myGrid.dim2;
00146   K3_1 = K3 - 1;
00147   // order = myGrid.order;
00148   stride=3*order;
00149   Mi = M + stride*istart;
00150   dMi = dM + stride*istart;
00151  
00152   for (i=istart; i<iend; i++) {
00153     double q;
00154     int32 u1, u2, u2i, u3i;
00155     q = p[i].cg;
00156     f1=f2=f3=0.0;
00157     u1 = (int)(p[i].x);
00158     u2i = (int)(p[i].y);
00159     u3i = (int)(p[i].z);
00160     u1 -= order;
00161     u2i -= order;
00162     u3i -= order;
00163     u3i += 1;
00164     for (j=0; j<order; j++) {
00165       double m1, d1;
00166       int32 ind1;
00167       m1=Mi[j]*q;
00168       d1=K1*dMi[j]*q;
00169       u1++;
00170       //ind1 = (u1 + (u1 < 0 ? K1 : 0))*dim2; 
00171       ind1 = (u1 + ((unsigned)u1 >>31)*K1)*dim2;
00172       u2 = u2i;
00173       for (k=0; k<order; k++) {
00174         double m2, d2, m1m2, m1d2, d1m2;
00175         int32 ind2;
00176         m2=Mi[order+k];
00177         d2=K2*dMi[order+k];
00178         m1m2=m1*m2;
00179         m1d2=m1*d2;
00180         d1m2=d1*m2;
00181         u2++;
00182         //ind2 = ind1 + (u2 + (u2 < 0 ? K2 : 0));
00183         ind2 = ind1 + (u2 + ((unsigned)u2 >>31)*K2);    
00184         const double *qline = q_arr[ind2];
00185         for (l=0; l<order; l++) {
00186           double term, m3, d3;
00187           int32 ind;
00188           m3=Mi[2*order+l];
00189           d3=K3*dMi[2*order+l];
00190           int32 u3 = u3i + l;
00191           ind = u3 - zstart;      
00192           //if (ind >= K3) ind -= K3;
00193           ind = ind - ((unsigned)(K3_1 -ind) >>31)*K3;
00194           
00195           //if (ind < 0 || ind >= zlen) printf ("ind >= zlen  (%d,%d,%d)\n", ind, zlen, zstart);          
00196           //assert (ind < zlen && ind >= 0);
00197 
00198           term = qline[ind];
00199           f1 -= d1m2 * m3 * term;
00200           f2 -= m1d2 * m3 * term;
00201           f3 -= m1m2 * d3 * term;
00202         }
00203       }
00204     }
00205     Mi += stride;
00206     dMi += stride;
00207     f[i].x = f1;
00208     f[i].y = f2;
00209     f[i].z = f3;
00210   }
00211 }

Generated on Tue Nov 21 01:17:14 2017 for NAMD by  doxygen 1.4.7