# PmeBase.inl

Go to the documentation of this file.
```00001
00007 #ifndef PME_BASE_INL__
00008 #define PME_BASE_INL__
00009
00010 #include "PmeBase.h"
00011 #include <math.h>
00012 #include <charm++.h>
00013 #include "Vector.h"
00014 #include "Lattice.h"
00015
00016
00017 static inline void scale_coordinates(PmeParticle p[], int N, Lattice lattice, PmeGrid grid) {
00018   Vector origin = lattice.origin();
00019   Vector recip1 = lattice.a_r();
00020   Vector recip2 = lattice.b_r();
00021   Vector recip3 = lattice.c_r();
00022   double ox = origin.x;
00023   double oy = origin.y;
00024   double oz = origin.z;
00025   double r1x = recip1.x;
00026   double r1y = recip1.y;
00027   double r1z = recip1.z;
00028   double r2x = recip2.x;
00029   double r2y = recip2.y;
00030   double r2z = recip2.z;
00031   double r3x = recip3.x;
00032   double r3y = recip3.y;
00033   double r3z = recip3.z;
00034   int K1 = grid.K1;
00035   int K2 = grid.K2;
00036   int K3 = grid.K3;
00037   double shift1 = ((K1 + grid.order - 1)/2)/(double)K1;
00038   double shift2 = ((K2 + grid.order - 1)/2)/(double)K2;
00039   double shift3 = ((K3 + grid.order - 1)/2)/(double)K3;
00040
00041   for (int i=0; i<N; i++) {
00042     double px = p[i].x - ox;
00043     double py = p[i].y - oy;
00044     double pz = p[i].z - oz;
00045     double sx = shift1 + px*r1x + py*r1y + pz*r1z;
00046     double sy = shift2 + px*r2x + py*r2y + pz*r2z;
00047     double sz = shift3 + px*r3x + py*r3y + pz*r3z;
00048     p[i].x = K1 * ( sx - floor(sx) );
00049     p[i].y = K2 * ( sy - floor(sy) );
00050     p[i].z = K3 * ( sz - floor(sz) );
00051     //  Check for rare rounding condition where K * ( 1 - epsilon ) == K
00052     //  which was observed with g++ on Intel x86 architecture.
00053     if ( p[i].x == K1 ) p[i].x = 0;
00054     if ( p[i].y == K2 ) p[i].y = 0;
00055     if ( p[i].z == K3 ) p[i].z = 0;
00056   }
00057 }
00058
00059
00060 static inline void scale_forces(Vector f[], int N, Lattice &lattice) {
00061   Vector recip1 = lattice.a_r();
00062   Vector recip2 = lattice.b_r();
00063   Vector recip3 = lattice.c_r();
00064   double r1x = recip1.x;
00065   double r1y = recip1.y;
00066   double r1z = recip1.z;
00067   double r2x = recip2.x;
00068   double r2y = recip2.y;
00069   double r2z = recip2.z;
00070   double r3x = recip3.x;
00071   double r3y = recip3.y;
00072   double r3z = recip3.z;
00073
00074   for (int i=0; i<N; i++) {
00075     double f1 = f[i].x;
00076     double f2 = f[i].y;
00077     double f3 = f[i].z;
00078     f[i].x = f1*r1x + f2*r2x + f3*r3x;
00079     f[i].y = f1*r1y + f2*r2y + f3*r3y;
00080     f[i].z = f1*r1z + f2*r2z + f3*r3z;
00081   }
00082 }
00083
00084
00085 template <class REAL>
00086 static inline void compute_b_spline(REAL * __restrict frac, REAL *M, REAL *dM, int order) {
00087   int j, n;
00088   REAL x,y,z,x1,y1,z1, div;
00089   REAL * __restrict Mx, * __restrict My, * __restrict Mz, * __restrict dMx, * __restrict dMy, * __restrict dMz;
00090   Mx=M-1; My=M+order-1; Mz=M+2*order-1;
00091   dMx=dM-1; dMy =dM+order-1; dMz=dM+2*order-1;
00092   x=frac[0];
00093   y=frac[1];
00094   z=frac[2];
00095   x1=1.0-x; y1=1.0-y; z1=1.0-z;
00096   /* Do n=3 case first */
00097   Mx[1]=.5*x1*x1;
00098   Mx[2]=x1*x + .5;
00099   Mx[3]=0.5*x*x;
00100   Mx[order]=0.0;
00101   My[1]=.5*y1*y1;
00102   My[2]=y1*y + .5;
00103   My[3]=0.5*y*y;
00104   My[order]=0.0;
00105   Mz[1]=.5*z1*z1;
00106   Mz[2]=z1*z + .5;
00107   Mz[3]=0.5*z*z;
00108   Mz[order]=0.0;
00109
00110   /* Recursively fill in the rest.  */
00111   for (n=4; n<=order-1; n++) {
00112     REAL div=1.0/(n-1);
00113     int j;
00114     Mx[n] = x*div*Mx[n-1];
00115     My[n] = y*div*My[n-1];
00116     Mz[n] = z*div*Mz[n-1];
00117     for (j=1; j<=n-2; j++) {
00118       Mx[n-j] = ((x+j)*Mx[n-j-1] + (n-x-j)*Mx[n-j])*div;
00119       My[n-j] = ((y+j)*My[n-j-1] + (n-y-j)*My[n-j])*div;
00120       Mz[n-j] = ((z+j)*Mz[n-j-1] + (n-z-j)*Mz[n-j])*div;
00121     }
00122     Mx[1] *= (1.0-x)*div;
00123     My[1] *= (1.0-y)*div;
00124     Mz[1] *= (1.0-z)*div;
00125   }
00126   /* Now do the derivatives  */
00127   dMx[1]=-Mx[1]; dMy[1]=-My[1]; dMz[1]=-Mz[1];
00128   for (j=2; j <= order; j++) {
00129     dMx[j] = Mx[j-1] - Mx[j];
00130     dMy[j] = My[j-1] - My[j];
00131     dMz[j] = Mz[j-1] - Mz[j];
00132   }
00133   /* Now finish the job!    */
00134   div=1.0/(order-1);
00135   Mx[order] = x*div*Mx[order-1];
00136   My[order] = y*div*My[order-1];
00137   Mz[order] = z*div*Mz[order-1];
00138   for (j=1; j<=order-2; j++) {
00139     Mx[order-j] = ((x+j)*Mx[order-j-1] + (order-x-j)*Mx[order-j])*div;
00140     My[order-j] = ((y+j)*My[order-j-1] + (order-y-j)*My[order-j])*div;
00141     Mz[order-j] = ((z+j)*Mz[order-j-1] + (order-z-j)*Mz[order-j])*div;
00142   }
00143   Mx[1] *= (1.0-x)*div;
00144   My[1] *= (1.0-y)*div;
00145   Mz[1] *= (1.0-z)*div;
00146 }
00147
00148
00149 #endif
00150
```

Generated on Fri May 25 01:17:14 2018 for NAMD by  1.4.7