Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

PmeBase.h

Go to the documentation of this file.
00001 
00007 #ifndef PME_BASE_H__
00008 #define PME_BASE_H__
00009 
00010 #include <math.h>
00011 #include <charm++.h>
00012 #include "MathArray.h"
00013 #include "Vector.h"
00014 #include "Lattice.h"
00015 
00016 #ifndef M_PI
00017 #define M_PI 3.14159265358979323846
00018 #endif
00019 
00020 struct PmeGrid {
00021   int K1, K2, K3;
00022   int dim2, dim3;
00023   int order;
00024   int block1, block2, block3;
00025   int xBlocks, yBlocks, zBlocks;
00026 };
00027 
00028 struct PmeParticle {
00029   double x, y, z;
00030   double cg;
00031 };
00032 
00033 void compute_b_spline(double frac[3], double *M, double *dM, int order);
00034 
00035 #define PME_MAX_EVALS 255 
00036 typedef MathArray<double,7> PmeReduction;
00037 
00038 #ifndef SQRT_PI
00039 #define SQRT_PI 1.7724538509055160273 /* mathematica 15 digits*/
00040 #endif
00041 
00042 
00043 static inline void scale_coordinates(PmeParticle p[], int N, Lattice lattice, PmeGrid grid) {
00044   Vector origin = lattice.origin();
00045   Vector recip1 = lattice.a_r();
00046   Vector recip2 = lattice.b_r();
00047   Vector recip3 = lattice.c_r();
00048   double ox = origin.x;
00049   double oy = origin.y;
00050   double oz = origin.z;
00051   double r1x = recip1.x;
00052   double r1y = recip1.y;
00053   double r1z = recip1.z;
00054   double r2x = recip2.x;
00055   double r2y = recip2.y;
00056   double r2z = recip2.z;
00057   double r3x = recip3.x;
00058   double r3y = recip3.y;
00059   double r3z = recip3.z;
00060   int K1 = grid.K1;
00061   int K2 = grid.K2;
00062   int K3 = grid.K3;
00063 
00064   for (int i=0; i<N; i++) {
00065     double px = p[i].x - ox;
00066     double py = p[i].y - oy;
00067     double pz = p[i].z - oz;
00068     double sx = px*r1x + py*r1y + pz*r1z;
00069     double sy = px*r2x + py*r2y + pz*r2z;
00070     double sz = px*r3x + py*r3y + pz*r3z;
00071     p[i].x = K1 * ( sx - floor(sx) );
00072     p[i].y = K2 * ( sy - floor(sy) );
00073     p[i].z = K3 * ( sz - floor(sz) );
00074     //  Check for rare rounding condition where K * ( 1 - epsilon ) == K
00075     //  which was observed with g++ on Intel x86 architecture.
00076     if ( p[i].x == K1 ) p[i].x = 0;
00077     if ( p[i].y == K2 ) p[i].y = 0;
00078     if ( p[i].z == K3 ) p[i].z = 0;
00079   }
00080 }
00081 
00082 
00083 static inline void scale_forces(Vector f[], int N, Lattice &lattice) {
00084   Vector recip1 = lattice.a_r();
00085   Vector recip2 = lattice.b_r();
00086   Vector recip3 = lattice.c_r();
00087   double r1x = recip1.x;
00088   double r1y = recip1.y;
00089   double r1z = recip1.z;
00090   double r2x = recip2.x;
00091   double r2y = recip2.y;
00092   double r2z = recip2.z;
00093   double r3x = recip3.x;
00094   double r3y = recip3.y;
00095   double r3z = recip3.z;
00096 
00097   for (int i=0; i<N; i++) {
00098     double f1 = f[i].x;
00099     double f2 = f[i].y;
00100     double f3 = f[i].z;
00101     f[i].x = f1*r1x + f2*r2x + f3*r3x;
00102     f[i].y = f1*r1y + f2*r2y + f3*r3y;
00103     f[i].z = f1*r1z + f2*r2z + f3*r3z;
00104   }
00105 }
00106 
00107 #endif

Generated on Sun Feb 12 04:07:56 2012 for NAMD by  doxygen 1.3.9.1