NAMD
PmeBase.inl
Go to the documentation of this file.
1
7 #ifndef PME_BASE_INL__
8 #define PME_BASE_INL__
9
10 #include "PmeBase.h"
11 #include <math.h>
12 #include <charm++.h>
13 #include "Vector.h"
14 #include "Lattice.h"
15
16
17 static inline void scale_coordinates(PmeParticle p[], int N, Lattice lattice, PmeGrid grid) {
18  Vector origin = lattice.origin();
19  Vector recip1 = lattice.a_r();
20  Vector recip2 = lattice.b_r();
21  Vector recip3 = lattice.c_r();
22  double ox = origin.x;
23  double oy = origin.y;
24  double oz = origin.z;
25  double r1x = recip1.x;
26  double r1y = recip1.y;
27  double r1z = recip1.z;
28  double r2x = recip2.x;
29  double r2y = recip2.y;
30  double r2z = recip2.z;
31  double r3x = recip3.x;
32  double r3y = recip3.y;
33  double r3z = recip3.z;
34  int K1 = grid.K1;
35  int K2 = grid.K2;
36  int K3 = grid.K3;
37  double shift1 = ((K1 + grid.order - 1)/2)/(double)K1;
38  double shift2 = ((K2 + grid.order - 1)/2)/(double)K2;
39  double shift3 = ((K3 + grid.order - 1)/2)/(double)K3;
40
41  for (int i=0; i<N; i++) {
42  double px = p[i].x - ox;
43  double py = p[i].y - oy;
44  double pz = p[i].z - oz;
45  double sx = shift1 + px*r1x + py*r1y + pz*r1z;
46  double sy = shift2 + px*r2x + py*r2y + pz*r2z;
47  double sz = shift3 + px*r3x + py*r3y + pz*r3z;
48  p[i].x = K1 * ( sx - floor(sx) );
49  p[i].y = K2 * ( sy - floor(sy) );
50  p[i].z = K3 * ( sz - floor(sz) );
51  // Check for rare rounding condition where K * ( 1 - epsilon ) == K
52  // which was observed with g++ on Intel x86 architecture.
53  if ( p[i].x == K1 ) p[i].x = 0;
54  if ( p[i].y == K2 ) p[i].y = 0;
55  if ( p[i].z == K3 ) p[i].z = 0;
56  }
57 }
58
59
60 static inline void scale_forces(Vector f[], int N, Lattice &lattice) {
61  Vector recip1 = lattice.a_r();
62  Vector recip2 = lattice.b_r();
63  Vector recip3 = lattice.c_r();
64  double r1x = recip1.x;
65  double r1y = recip1.y;
66  double r1z = recip1.z;
67  double r2x = recip2.x;
68  double r2y = recip2.y;
69  double r2z = recip2.z;
70  double r3x = recip3.x;
71  double r3y = recip3.y;
72  double r3z = recip3.z;
73
74  for (int i=0; i<N; i++) {
75  double f1 = f[i].x;
76  double f2 = f[i].y;
77  double f3 = f[i].z;
78  f[i].x = f1*r1x + f2*r2x + f3*r3x;
79  f[i].y = f1*r1y + f2*r2y + f3*r3y;
80  f[i].z = f1*r1z + f2*r2z + f3*r3z;
81  }
82 }
83
84
85 template <class REAL>
86 static inline void compute_b_spline(REAL * __restrict frac, REAL *M, REAL *dM, int order) {
87  int j, n;
88  REAL x,y,z,x1,y1,z1, div;
89  REAL * __restrict Mx, * __restrict My, * __restrict Mz, * __restrict dMx, * __restrict dMy, * __restrict dMz;
90  Mx=M-1; My=M+order-1; Mz=M+2*order-1;
91  dMx=dM-1; dMy =dM+order-1; dMz=dM+2*order-1;
92  x=frac[0];
93  y=frac[1];
94  z=frac[2];
95  x1=1.0-x; y1=1.0-y; z1=1.0-z;
96  /* Do n=3 case first */
97  Mx[1]=.5*x1*x1;
98  Mx[2]=x1*x + .5;
99  Mx[3]=0.5*x*x;
100  Mx[order]=0.0;
101  My[1]=.5*y1*y1;
102  My[2]=y1*y + .5;
103  My[3]=0.5*y*y;
104  My[order]=0.0;
105  Mz[1]=.5*z1*z1;
106  Mz[2]=z1*z + .5;
107  Mz[3]=0.5*z*z;
108  Mz[order]=0.0;
109
110  /* Recursively fill in the rest. */
111  for (n=4; n<=order-1; n++) {
112  REAL div=1.0/(n-1);
113  int j;
114  Mx[n] = x*div*Mx[n-1];
115  My[n] = y*div*My[n-1];
116  Mz[n] = z*div*Mz[n-1];
117  for (j=1; j<=n-2; j++) {
118  Mx[n-j] = ((x+j)*Mx[n-j-1] + (n-x-j)*Mx[n-j])*div;
119  My[n-j] = ((y+j)*My[n-j-1] + (n-y-j)*My[n-j])*div;
120  Mz[n-j] = ((z+j)*Mz[n-j-1] + (n-z-j)*Mz[n-j])*div;
121  }
122  Mx[1] *= (1.0-x)*div;
123  My[1] *= (1.0-y)*div;
124  Mz[1] *= (1.0-z)*div;
125  }
126  /* Now do the derivatives */
127  dMx[1]=-Mx[1]; dMy[1]=-My[1]; dMz[1]=-Mz[1];
128  for (j=2; j <= order; j++) {
129  dMx[j] = Mx[j-1] - Mx[j];
130  dMy[j] = My[j-1] - My[j];
131  dMz[j] = Mz[j-1] - Mz[j];
132  }
133  /* Now finish the job! */
134  div=1.0/(order-1);
135  Mx[order] = x*div*Mx[order-1];
136  My[order] = y*div*My[order-1];
137  Mz[order] = z*div*Mz[order-1];
138  for (j=1; j<=order-2; j++) {
139  Mx[order-j] = ((x+j)*Mx[order-j-1] + (order-x-j)*Mx[order-j])*div;
140  My[order-j] = ((y+j)*My[order-j-1] + (order-y-j)*My[order-j])*div;
141  Mz[order-j] = ((z+j)*Mz[order-j-1] + (order-z-j)*Mz[order-j])*div;
142  }
143  Mx[1] *= (1.0-x)*div;
144  My[1] *= (1.0-y)*div;
145  Mz[1] *= (1.0-z)*div;
146 }
147
148
149 #endif
150
static void scale_coordinates(PmeParticle p[], int N, Lattice lattice, PmeGrid grid)
Definition: PmeBase.inl:17
Vector a_r() const
Definition: Lattice.h:268
double x
Definition: PmeBase.h:26
Definition: Vector.h:64
int K2
Definition: PmeBase.h:18
int K1
Definition: PmeBase.h:18
Vector c_r() const
Definition: Lattice.h:270
BigReal z
Definition: Vector.h:66
double z
Definition: PmeBase.h:26
double y
Definition: PmeBase.h:26
Vector origin() const
Definition: Lattice.h:262
Vector b_r() const
Definition: Lattice.h:269
#define order
Definition: PmeRealSpace.C:235
int order
Definition: PmeBase.h:20
gridSize z
static void scale_forces(Vector f[], int N, Lattice &lattice)
Definition: PmeBase.inl:60
BigReal x
Definition: Vector.h:66
int K3
Definition: PmeBase.h:18
BigReal y
Definition: Vector.h:66
gridSize y
gridSize x
static void compute_b_spline(REAL *__restrict frac, REAL *M, REAL *dM, int order)
Definition: PmeBase.inl:86