NAMD
GridForceGrid.inl
Go to the documentation of this file.
1 
7 #ifndef GRIDFORCEGRID_INL
8 #define GRIDFORCEGRID_INL
9 
10 #include "GridForceGrid.h"
11 
12 inline int GridforceFullBaseGrid::compute_VdV(Position pos, float &V, Vector &dV) const
13 {
14  //SimParameters *simParams = Node::Object()->simParameters;
15  int inds[3];
16  Vector g, dg;
17  Vector gapscale = Vector(1, 1, 1);
18 
19  int err = get_inds(pos, inds, dg, gapscale);
20  if (err) {
21  return -1;
22  }
23 
24  DebugM(1, "gapscale = " << gapscale << "\n");
25  DebugM(1, "dg = " << dg << "\n");
26  DebugM(1, "ind + dg = " << inds[0]+dg[0] << " " << inds[1]+dg[1] << " " << inds[2]+dg[2] << "\n");
27  DebugM(3, "compute_VdV: generation = " << generation << "\n" << endi);
28 
29  // Pass to subgrid if one exists here
30  for (int i = 0; i < numSubgrids; i++) {
31  if (((inds[0] >= subgrids[i]->pmin[0] && inds[0] <= subgrids[i]->pmax[0]) || subgrids[i]->cont[0]) &&
32  ((inds[1] >= subgrids[i]->pmin[1] && inds[1] <= subgrids[i]->pmax[1]) || subgrids[i]->cont[1]) &&
33  ((inds[2] >= subgrids[i]->pmin[2] && inds[2] <= subgrids[i]->pmax[2]) || subgrids[i]->cont[2]))
34  {
35  return subgrids[i]->compute_VdV(pos, V, dV);
36  }
37  }
38 
39  // Compute b
40  float b[64]; // Matrix of values at 8 box corners
41  compute_b(b, inds, gapscale);
42  for (int j = 0; j < 64; j++) DebugM(1, "b[" << j << "] = " << b[j] << "\n" << endi);
43 
44  // Compute a
45  float a[64];
46  compute_a(a, b);
47  for (int j = 0; j < 64; j++) DebugM(1, "a[" << j << "] = " << a[j] << "\n" << endi);
48 
49  // Calculate powers of x, y, z for later use
50  // e.g. x[2] = x^2
51  float x[4], y[4], z[4];
52  x[0] = 1; y[0] = 1; z[0] = 1;
53  for (int j = 1; j < 4; j++) {
54  x[j] = x[j-1] * dg.x;
55  y[j] = y[j-1] * dg.y;
56  z[j] = z[j-1] * dg.z;
57  }
58 
59  V = compute_V(a, x, y, z);
60  dV = Tensor::diagonal(gapscale) * (compute_dV(a, x, y, z) * inv);
61 
62  return 0;
63 }
64 
65 
66 inline int GridforceLiteGrid::compute_VdV(Position pos, float &V, Vector &dV) const
67 {
68  int inds[3];
69  Vector g, dg;
70 
71  int err = get_inds(pos, inds, dg);
72  if (err) {
73  return -1;
74  }
75 
76  float wts[4][8];
77  float results[4];
78 
79  // compute_wts(wts, dg);
80  // wts[0][0] = (1-dg.x) * (1-dg.y) * (1-dg.z);
81  // wts[0][1] = (1-dg.x) * (1-dg.y) * dg.z;
82  // wts[0][2] = (1-dg.x) * dg.y * (1-dg.z);
83  // wts[0][3] = (1-dg.x) * dg.y * dg.z;
84  // wts[0][4] = dg.x * (1-dg.y) * (1-dg.z);
85  // wts[0][5] = dg.x * (1-dg.y) * dg.z;
86  // wts[0][6] = dg.x * dg.y * (1-dg.z);
87  // wts[0][7] = dg.x * dg.y * dg.z;
88 
89  int i = 1;
90  wts[i][0] = -(1-dg.y) * (1-dg.z);
91  wts[i][1] = -(1-dg.y) * dg.z;
92  wts[i][2] = - dg.y * (1-dg.z);
93  wts[i][3] = - dg.y * dg.z;
94  for (int j=0; j<4; j++) wts[i][j+4] = -wts[i][j];
95 
96  i = 2;
97  wts[i][0] = -(1-dg.x) * (1-dg.z);
98  wts[i][1] = -(1-dg.x) * dg.z;
99  wts[i][2] = -wts[i][0];
100  wts[i][3] = -wts[i][1];
101  wts[i][4] = - dg.x * (1-dg.z);
102  wts[i][5] = - dg.x * dg.z;
103  wts[i][6] = -wts[i][4];
104  wts[i][7] = -wts[i][5];
105 
106  i = 3;
107  wts[i][0] = - (1-dg.x) * (1-dg.y);
108  wts[i][1] = -wts[i][0];
109  wts[i][2] = - (1-dg.x) * dg.y ;
110  wts[i][3] = -wts[i][2];
111  wts[i][4] = - dg.x * (1-dg.y);
112  wts[i][5] = -wts[i][4];
113  wts[i][6] = - dg.x * dg.y ;
114  wts[i][7] = -wts[i][6];
115 
116  i = 0;
117  for (int j=0; j<4; j++) wts[i][j] = (1-dg.x) * wts[i+1][j+4];
118  for (int j=0; j<4; j++) wts[i][j+4] = dg.x * wts[i+1][j+4];
119 
120  for (i = 0; i < 4; i++) {
121  results[i] = linear_interpolate(inds[0], inds[1], inds[2], 0, wts[i]);
122  }
123 
124  V = results[0];
125  dV = Vector(results[1], results[2], results[3]) * inv;
126 
127  return 0;
128 }
129 
130 
131 inline int GridforceFullBaseGrid::get_inds(Position pos, int *inds, Vector &dg, Vector &gapscale) const
132 {
133  Vector p = pos - origin;
134  Vector g;
135 
136  g = inv * p;
137 
138  for (int i = 0; i < 3; i++) {
139  inds[i] = (int)floor(g[i]);
140  dg[i] = g[i] - inds[i];
141  }
142 
143  for (int i = 0; i < 3; i++) {
144  if (inds[i] < 0 || inds[i] >= k[i]-1) {
145  if (cont[i]) inds[i] = k[i]-1;
146  else return -1; // Outside potential and grid is not continuous
147  }
148  if (cont[i] && inds[i] == k[i]-1) {
149  // Correct for non-unit spacing between continuous grid images
150  gapscale[i] *= gapinv[i];
151  if (g[i] < 0.0) dg[i] = 1.0 + g[i]*gapinv[i]; // = (gap[i] + g[i]) * gapinv[i]
152  else dg[i] = (g[i] - inds[i]) * gapinv[i];
153  }
154  }
155 
156  return 0;
157 }
158 
159 
160 inline float GridforceFullBaseGrid::compute_V(float *a, float *x, float *y, float *z) const
161 {
162  float V = 0.0;
163  long int ind = 0;
164  for (int l = 0; l < 4; l++) {
165  for (int k = 0; k < 4; k++) {
166  for (int j = 0; j < 4; j++) {
167  V += a[ind] * x[j] * y[k] * z[l];
168  ind++;
169  }
170  }
171  }
172  return V;
173 }
174 
175 
176 inline Vector GridforceFullBaseGrid::compute_dV(float *a, float *x, float *y, float *z) const
177 {
178  Vector dV = 0;
179  long int ind = 0;
180  for (int l = 0; l < 4; l++) {
181  for (int k = 0; k < 4; k++) {
182  for (int j = 0; j < 4; j++) {
183  if (j > 0) dV.x += a[ind] * j * x[j-1] * y[k] * z[l]; // dV/dx
184  if (k > 0) dV.y += a[ind] * k * x[j] * y[k-1] * z[l]; // dV/dy
185  if (l > 0) dV.z += a[ind] * l * x[j] * y[k] * z[l-1]; // dV/dz
186  ind++;
187  }
188  }
189  }
190  return dV;
191 }
192 
193 
194 inline Vector GridforceFullBaseGrid::compute_d2V(float *a, float *x, float *y, float *z) const
195 {
196  Vector d2V = 0;
197  int ind = 0;
198  for (int l = 0; l < 4; l++) {
199  for (int k = 0; k < 4; k++) {
200  for (int j = 0; j < 4; j++) {
201  if (j > 0 && k > 0) d2V.x += a[ind] * j * k * x[j-1] * y[k-1] * z[l]; // d2V/dxdy
202  if (j > 0 && l > 0) d2V.y += a[ind] * j * l * x[j-1] * y[k] * z[l-1]; // d2V/dxdz
203  if (k > 0 && l > 0) d2V.z += a[ind] * k * l * x[j] * y[k-1] * z[l-1]; // d2V/dydz
204  ind++;
205  }
206  }
207  }
208  return d2V;
209 }
210 
211 
212 inline float GridforceFullBaseGrid::compute_d3V(float *a, float *x, float *y, float *z) const
213 {
214  float d3V = 0.0;
215  long int ind = 0;
216  for (int l = 0; l < 4; l++) {
217  for (int k = 0; k < 4; k++) {
218  for (int j = 0; j < 4; j++) {
219  if (j > 0 && k > 0 && l > 0) d3V += a[ind] * j * k * l * x[j-1] * y[k-1] * z[l-1]; // d3V/dxdydz
220  ind++;
221  }
222  }
223  }
224  return d3V;
225 }
226 
227 
228 inline void GridforceFullBaseGrid::compute_a(float *a, float *b) const
229 {
230  // Static sparse 64x64 matrix times vector ... nicer looking way than this?
231  a[0] = b[0];
232  a[1] = b[8];
233  a[2] = -3*b[0] + 3*b[1] - 2*b[8] - b[9];
234  a[3] = 2*b[0] - 2*b[1] + b[8] + b[9];
235  a[4] = b[16];
236  a[5] = b[32];
237  a[6] = -3*b[16] + 3*b[17] - 2*b[32] - b[33];
238  a[7] = 2*b[16] - 2*b[17] + b[32] + b[33];
239  a[8] = -3*b[0] + 3*b[2] - 2*b[16] - b[18];
240  a[9] = -3*b[8] + 3*b[10] - 2*b[32] - b[34];
241  a[10] = 9*b[0] - 9*b[1] - 9*b[2] + 9*b[3] + 6*b[8] + 3*b[9] - 6*b[10] - 3*b[11]
242  + 6*b[16] - 6*b[17] + 3*b[18] - 3*b[19] + 4*b[32] + 2*b[33] + 2*b[34] + b[35];
243  a[11] = -6*b[0] + 6*b[1] + 6*b[2] - 6*b[3] - 3*b[8] - 3*b[9] + 3*b[10] + 3*b[11]
244  - 4*b[16] + 4*b[17] - 2*b[18] + 2*b[19] - 2*b[32] - 2*b[33] - b[34] - b[35];
245  a[12] = 2*b[0] - 2*b[2] + b[16] + b[18];
246  a[13] = 2*b[8] - 2*b[10] + b[32] + b[34];
247  a[14] = -6*b[0] + 6*b[1] + 6*b[2] - 6*b[3] - 4*b[8] - 2*b[9] + 4*b[10] + 2*b[11]
248  - 3*b[16] + 3*b[17] - 3*b[18] + 3*b[19] - 2*b[32] - b[33] - 2*b[34] - b[35];
249  a[15] = 4*b[0] - 4*b[1] - 4*b[2] + 4*b[3] + 2*b[8] + 2*b[9] - 2*b[10] - 2*b[11]
250  + 2*b[16] - 2*b[17] + 2*b[18] - 2*b[19] + b[32] + b[33] + b[34] + b[35];
251  a[16] = b[24];
252  a[17] = b[40];
253  a[18] = -3*b[24] + 3*b[25] - 2*b[40] - b[41];
254  a[19] = 2*b[24] - 2*b[25] + b[40] + b[41];
255  a[20] = b[48];
256  a[21] = b[56];
257  a[22] = -3*b[48] + 3*b[49] - 2*b[56] - b[57];
258  a[23] = 2*b[48] - 2*b[49] + b[56] + b[57];
259  a[24] = -3*b[24] + 3*b[26] - 2*b[48] - b[50];
260  a[25] = -3*b[40] + 3*b[42] - 2*b[56] - b[58];
261  a[26] = 9*b[24] - 9*b[25] - 9*b[26] + 9*b[27] + 6*b[40] + 3*b[41] - 6*b[42] - 3*b[43]
262  + 6*b[48] - 6*b[49] + 3*b[50] - 3*b[51] + 4*b[56] + 2*b[57] + 2*b[58] + b[59];
263  a[27] = -6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 3*b[40] - 3*b[41] + 3*b[42] + 3*b[43]
264  - 4*b[48] + 4*b[49] - 2*b[50] + 2*b[51] - 2*b[56] - 2*b[57] - b[58] - b[59];
265  a[28] = 2*b[24] - 2*b[26] + b[48] + b[50];
266  a[29] = 2*b[40] - 2*b[42] + b[56] + b[58];
267  a[30] = -6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 4*b[40] - 2*b[41] + 4*b[42] + 2*b[43]
268  - 3*b[48] + 3*b[49] - 3*b[50] + 3*b[51] - 2*b[56] - b[57] - 2*b[58] - b[59];
269  a[31] = 4*b[24] - 4*b[25] - 4*b[26] + 4*b[27] + 2*b[40] + 2*b[41] - 2*b[42] - 2*b[43]
270  + 2*b[48] - 2*b[49] + 2*b[50] - 2*b[51] + b[56] + b[57] + b[58] + b[59];
271  a[32] = -3*b[0] + 3*b[4] - 2*b[24] - b[28];
272  a[33] = -3*b[8] + 3*b[12] - 2*b[40] - b[44];
273  a[34] = 9*b[0] - 9*b[1] - 9*b[4] + 9*b[5] + 6*b[8] + 3*b[9] - 6*b[12] - 3*b[13]
274  + 6*b[24] - 6*b[25] + 3*b[28] - 3*b[29] + 4*b[40] + 2*b[41] + 2*b[44] + b[45];
275  a[35] = -6*b[0] + 6*b[1] + 6*b[4] - 6*b[5] - 3*b[8] - 3*b[9] + 3*b[12] + 3*b[13]
276  - 4*b[24] + 4*b[25] - 2*b[28] + 2*b[29] - 2*b[40] - 2*b[41] - b[44] - b[45];
277  a[36] = -3*b[16] + 3*b[20] - 2*b[48] - b[52];
278  a[37] = -3*b[32] + 3*b[36] - 2*b[56] - b[60];
279  a[38] = 9*b[16] - 9*b[17] - 9*b[20] + 9*b[21] + 6*b[32] + 3*b[33] - 6*b[36] - 3*b[37]
280  + 6*b[48] - 6*b[49] + 3*b[52] - 3*b[53] + 4*b[56] + 2*b[57] + 2*b[60] + b[61];
281  a[39] = -6*b[16] + 6*b[17] + 6*b[20] - 6*b[21] - 3*b[32] - 3*b[33] + 3*b[36] + 3*b[37]
282  - 4*b[48] + 4*b[49] - 2*b[52] + 2*b[53] - 2*b[56] - 2*b[57] - b[60] - b[61];
283  a[40] = 9*b[0] - 9*b[2] - 9*b[4] + 9*b[6] + 6*b[16] + 3*b[18] - 6*b[20] - 3*b[22]
284  + 6*b[24] - 6*b[26] + 3*b[28] - 3*b[30] + 4*b[48] + 2*b[50] + 2*b[52] + b[54];
285  a[41] = 9*b[8] - 9*b[10] - 9*b[12] + 9*b[14] + 6*b[32] + 3*b[34] - 6*b[36] - 3*b[38]
286  + 6*b[40] - 6*b[42] + 3*b[44] - 3*b[46] + 4*b[56] + 2*b[58] + 2*b[60] + b[62];
287  a[42] = -27*b[0] + 27*b[1] + 27*b[2] - 27*b[3] + 27*b[4] - 27*b[5] - 27*b[6] + 27*b[7]
288  - 18*b[8] - 9*b[9] + 18*b[10] + 9*b[11] + 18*b[12] + 9*b[13] - 18*b[14] - 9*b[15]
289  - 18*b[16] + 18*b[17] - 9*b[18] + 9*b[19] + 18*b[20] - 18*b[21] + 9*b[22] - 9*b[23]
290  - 18*b[24] + 18*b[25] + 18*b[26] - 18*b[27] - 9*b[28] + 9*b[29] + 9*b[30] - 9*b[31]
291  - 12*b[32] - 6*b[33] - 6*b[34] - 3*b[35] + 12*b[36] + 6*b[37] + 6*b[38] + 3*b[39]
292  - 12*b[40] - 6*b[41] + 12*b[42] + 6*b[43] - 6*b[44] - 3*b[45] + 6*b[46] + 3*b[47]
293  - 12*b[48] + 12*b[49] - 6*b[50] + 6*b[51] - 6*b[52] + 6*b[53] - 3*b[54] + 3*b[55]
294  - 8*b[56] - 4*b[57] - 4*b[58] - 2*b[59] - 4*b[60] - 2*b[61] - 2*b[62] - b[63];
295  a[43] = 18*b[0] - 18*b[1] - 18*b[2] + 18*b[3] - 18*b[4] + 18*b[5] + 18*b[6] - 18*b[7]
296  + 9*b[8] + 9*b[9] - 9*b[10] - 9*b[11] - 9*b[12] - 9*b[13] + 9*b[14] + 9*b[15]
297  + 12*b[16] - 12*b[17] + 6*b[18] - 6*b[19] - 12*b[20] + 12*b[21] - 6*b[22] + 6*b[23]
298  + 12*b[24] - 12*b[25] - 12*b[26] + 12*b[27] + 6*b[28] - 6*b[29] - 6*b[30] + 6*b[31]
299  + 6*b[32] + 6*b[33] + 3*b[34] + 3*b[35] - 6*b[36] - 6*b[37] - 3*b[38] - 3*b[39]
300  + 6*b[40] + 6*b[41] - 6*b[42] - 6*b[43] + 3*b[44] + 3*b[45] - 3*b[46] - 3*b[47]
301  + 8*b[48] - 8*b[49] + 4*b[50] - 4*b[51] + 4*b[52] - 4*b[53] + 2*b[54] - 2*b[55]
302  + 4*b[56] + 4*b[57] + 2*b[58] + 2*b[59] + 2*b[60] + 2*b[61] + b[62] + b[63];
303  a[44] = -6*b[0] + 6*b[2] + 6*b[4] - 6*b[6] - 3*b[16] - 3*b[18] + 3*b[20] + 3*b[22]
304  - 4*b[24] + 4*b[26] - 2*b[28] + 2*b[30] - 2*b[48] - 2*b[50] - b[52] - b[54];
305  a[45] = -6*b[8] + 6*b[10] + 6*b[12] - 6*b[14] - 3*b[32] - 3*b[34] + 3*b[36] + 3*b[38]
306  - 4*b[40] + 4*b[42] - 2*b[44] + 2*b[46] - 2*b[56] - 2*b[58] - b[60] - b[62];
307  a[46] = 18*b[0] - 18*b[1] - 18*b[2] + 18*b[3] - 18*b[4] + 18*b[5] + 18*b[6] - 18*b[7]
308  + 12*b[8] + 6*b[9] - 12*b[10] - 6*b[11] - 12*b[12] - 6*b[13] + 12*b[14] + 6*b[15]
309  + 9*b[16] - 9*b[17] + 9*b[18] - 9*b[19] - 9*b[20] + 9*b[21] - 9*b[22] + 9*b[23]
310  + 12*b[24] - 12*b[25] - 12*b[26] + 12*b[27] + 6*b[28] - 6*b[29] - 6*b[30] + 6*b[31]
311  + 6*b[32] + 3*b[33] + 6*b[34] + 3*b[35] - 6*b[36] - 3*b[37] - 6*b[38] - 3*b[39]
312  + 8*b[40] + 4*b[41] - 8*b[42] - 4*b[43] + 4*b[44] + 2*b[45] - 4*b[46] - 2*b[47]
313  + 6*b[48] - 6*b[49] + 6*b[50] - 6*b[51] + 3*b[52] - 3*b[53] + 3*b[54] - 3*b[55]
314  + 4*b[56] + 2*b[57] + 4*b[58] + 2*b[59] + 2*b[60] + b[61] + 2*b[62] + b[63];
315  a[47] = -12*b[0] + 12*b[1] + 12*b[2] - 12*b[3] + 12*b[4] - 12*b[5] - 12*b[6] + 12*b[7]
316  - 6*b[8] - 6*b[9] + 6*b[10] + 6*b[11] + 6*b[12] + 6*b[13] - 6*b[14] - 6*b[15]
317  - 6*b[16] + 6*b[17] - 6*b[18] + 6*b[19] + 6*b[20] - 6*b[21] + 6*b[22] - 6*b[23]
318  - 8*b[24] + 8*b[25] + 8*b[26] - 8*b[27] - 4*b[28] + 4*b[29] + 4*b[30] - 4*b[31]
319  - 3*b[32] - 3*b[33] - 3*b[34] - 3*b[35] + 3*b[36] + 3*b[37] + 3*b[38] + 3*b[39]
320  - 4*b[40] - 4*b[41] + 4*b[42] + 4*b[43] - 2*b[44] - 2*b[45] + 2*b[46] + 2*b[47]
321  - 4*b[48] + 4*b[49] - 4*b[50] + 4*b[51] - 2*b[52] + 2*b[53] - 2*b[54] + 2*b[55]
322  - 2*b[56] - 2*b[57] - 2*b[58] - 2*b[59] - b[60] - b[61] - b[62] - b[63];
323  a[48] = 2*b[0] - 2*b[4] + b[24] + b[28];
324  a[49] = 2*b[8] - 2*b[12] + b[40] + b[44];
325  a[50] = -6*b[0] + 6*b[1] + 6*b[4] - 6*b[5] - 4*b[8] - 2*b[9] + 4*b[12] + 2*b[13]
326  - 3*b[24] + 3*b[25] - 3*b[28] + 3*b[29] - 2*b[40] - b[41] - 2*b[44] - b[45];
327  a[51] = 4*b[0] - 4*b[1] - 4*b[4] + 4*b[5] + 2*b[8] + 2*b[9] - 2*b[12] - 2*b[13]
328  + 2*b[24] - 2*b[25] + 2*b[28] - 2*b[29] + b[40] + b[41] + b[44] + b[45];
329  a[52] = 2*b[16] - 2*b[20] + b[48] + b[52];
330  a[53] = 2*b[32] - 2*b[36] + b[56] + b[60];
331  a[54] = -6*b[16] + 6*b[17] + 6*b[20] - 6*b[21] - 4*b[32] - 2*b[33] + 4*b[36] + 2*b[37]
332  - 3*b[48] + 3*b[49] - 3*b[52] + 3*b[53] - 2*b[56] - b[57] - 2*b[60] - b[61];
333  a[55] = 4*b[16] - 4*b[17] - 4*b[20] + 4*b[21] + 2*b[32] + 2*b[33] - 2*b[36] - 2*b[37]
334  + 2*b[48] - 2*b[49] + 2*b[52] - 2*b[53] + b[56] + b[57] + b[60] + b[61];
335  a[56] = -6*b[0] + 6*b[2] + 6*b[4] - 6*b[6] - 4*b[16] - 2*b[18] + 4*b[20] + 2*b[22]
336  - 3*b[24] + 3*b[26] - 3*b[28] + 3*b[30] - 2*b[48] - b[50] - 2*b[52] - b[54];
337  a[57] = -6*b[8] + 6*b[10] + 6*b[12] - 6*b[14] - 4*b[32] - 2*b[34] + 4*b[36] + 2*b[38]
338  - 3*b[40] + 3*b[42] - 3*b[44] + 3*b[46] - 2*b[56] - b[58] - 2*b[60] - b[62];
339  a[58] = 18*b[0] - 18*b[1] - 18*b[2] + 18*b[3] - 18*b[4] + 18*b[5] + 18*b[6] - 18*b[7]
340  + 12*b[8] + 6*b[9] - 12*b[10] - 6*b[11] - 12*b[12] - 6*b[13] + 12*b[14] + 6*b[15]
341  + 12*b[16] - 12*b[17] + 6*b[18] - 6*b[19] - 12*b[20] + 12*b[21] - 6*b[22] + 6*b[23]
342  + 9*b[24] - 9*b[25] - 9*b[26] + 9*b[27] + 9*b[28] - 9*b[29] - 9*b[30] + 9*b[31]
343  + 8*b[32] + 4*b[33] + 4*b[34] + 2*b[35] - 8*b[36] - 4*b[37] - 4*b[38] - 2*b[39]
344  + 6*b[40] + 3*b[41] - 6*b[42] - 3*b[43] + 6*b[44] + 3*b[45] - 6*b[46] - 3*b[47]
345  + 6*b[48] - 6*b[49] + 3*b[50] - 3*b[51] + 6*b[52] - 6*b[53] + 3*b[54] - 3*b[55]
346  + 4*b[56] + 2*b[57] + 2*b[58] + b[59] + 4*b[60] + 2*b[61] + 2*b[62] + b[63];
347  a[59] = -12*b[0] + 12*b[1] + 12*b[2] - 12*b[3] + 12*b[4] - 12*b[5] - 12*b[6] + 12*b[7]
348  - 6*b[8] - 6*b[9] + 6*b[10] + 6*b[11] + 6*b[12] + 6*b[13] - 6*b[14] - 6*b[15]
349  - 8*b[16] + 8*b[17] - 4*b[18] + 4*b[19] + 8*b[20] - 8*b[21] + 4*b[22] - 4*b[23]
350  - 6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 6*b[28] + 6*b[29] + 6*b[30] - 6*b[31]
351  - 4*b[32] - 4*b[33] - 2*b[34] - 2*b[35] + 4*b[36] + 4*b[37] + 2*b[38] + 2*b[39]
352  - 3*b[40] - 3*b[41] + 3*b[42] + 3*b[43] - 3*b[44] - 3*b[45] + 3*b[46] + 3*b[47]
353  - 4*b[48] + 4*b[49] - 2*b[50] + 2*b[51] - 4*b[52] + 4*b[53] - 2*b[54] + 2*b[55]
354  - 2*b[56] - 2*b[57] - b[58] - b[59] - 2*b[60] - 2*b[61] - b[62] - b[63];
355  a[60] = 4*b[0] - 4*b[2] - 4*b[4] + 4*b[6] + 2*b[16] + 2*b[18] - 2*b[20] - 2*b[22]
356  + 2*b[24] - 2*b[26] + 2*b[28] - 2*b[30] + b[48] + b[50] + b[52] + b[54];
357  a[61] = 4*b[8] - 4*b[10] - 4*b[12] + 4*b[14] + 2*b[32] + 2*b[34] - 2*b[36] - 2*b[38]
358  + 2*b[40] - 2*b[42] + 2*b[44] - 2*b[46] + b[56] + b[58] + b[60] + b[62];
359  a[62] = -12*b[0] + 12*b[1] + 12*b[2] - 12*b[3] + 12*b[4] - 12*b[5] - 12*b[6] + 12*b[7]
360  - 8*b[8] - 4*b[9] + 8*b[10] + 4*b[11] + 8*b[12] + 4*b[13] - 8*b[14] - 4*b[15]
361  - 6*b[16] + 6*b[17] - 6*b[18] + 6*b[19] + 6*b[20] - 6*b[21] + 6*b[22] - 6*b[23]
362  - 6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 6*b[28] + 6*b[29] + 6*b[30] - 6*b[31]
363  - 4*b[32] - 2*b[33] - 4*b[34] - 2*b[35] + 4*b[36] + 2*b[37] + 4*b[38] + 2*b[39]
364  - 4*b[40] - 2*b[41] + 4*b[42] + 2*b[43] - 4*b[44] - 2*b[45] + 4*b[46] + 2*b[47]
365  - 3*b[48] + 3*b[49] - 3*b[50] + 3*b[51] - 3*b[52] + 3*b[53] - 3*b[54] + 3*b[55]
366  - 2*b[56] - b[57] - 2*b[58] - b[59] - 2*b[60] - b[61] - 2*b[62] - b[63];
367  a[63] = 8*b[0] - 8*b[1] - 8*b[2] + 8*b[3] - 8*b[4] + 8*b[5] + 8*b[6] - 8*b[7]
368  + 4*b[8] + 4*b[9] - 4*b[10] - 4*b[11] - 4*b[12] - 4*b[13] + 4*b[14] + 4*b[15]
369  + 4*b[16] - 4*b[17] + 4*b[18] - 4*b[19] - 4*b[20] + 4*b[21] - 4*b[22] + 4*b[23]
370  + 4*b[24] - 4*b[25] - 4*b[26] + 4*b[27] + 4*b[28] - 4*b[29] - 4*b[30] + 4*b[31]
371  + 2*b[32] + 2*b[33] + 2*b[34] + 2*b[35] - 2*b[36] - 2*b[37] - 2*b[38] - 2*b[39]
372  + 2*b[40] + 2*b[41] - 2*b[42] - 2*b[43] + 2*b[44] + 2*b[45] - 2*b[46] - 2*b[47]
373  + 2*b[48] - 2*b[49] + 2*b[50] - 2*b[51] + 2*b[52] - 2*b[53] + 2*b[54] - 2*b[55]
374  + b[56] + b[57] + b[58] + b[59] + b[60] + b[61] + b[62] + b[63];
375 }
376 
377 
378 inline int GridforceLiteGrid::get_inds(Position pos, int *inds, Vector &dg) const
379 {
380  Vector p = pos - origin;
381  Vector g;
382 
383  g = inv * p;
384 
385  for (int i = 0; i < 3; i++) {
386  inds[i] = (int)floor(g[i]);
387  dg[i] = g[i] - inds[i];
388  }
389 
390  for (int i = 0; i < 3; i++) {
391  if (inds[i] < 0 || inds[i] >= k[i]-1) {
392  return -1; // Outside potential and grid is not continuous
393  }
394  }
395 
396  return 0;
397 }
398 
399 
400 inline void GridforceLiteGrid::compute_wts(float *wts, const Vector &dg) const
401 {
402  wts[0] = (1-dg.x) * (1-dg.y) * (1-dg.z);
403  wts[1] = (1-dg.x) * (1-dg.y) * dg.z;
404  wts[2] = (1-dg.x) * dg.y * (1-dg.z);
405  wts[3] = (1-dg.x) * dg.y * dg.z;
406  wts[4] = dg.x * (1-dg.y) * (1-dg.z);
407  wts[5] = dg.x * (1-dg.y) * dg.z;
408  wts[6] = dg.x * dg.y * (1-dg.z);
409  wts[7] = dg.x * dg.y * dg.z;
410  DebugM(2, "dg = " << dg << "\n" << endi);
411 }
412 
413 
414 inline float GridforceLiteGrid::linear_interpolate(int i0, int i1, int i2, int i3, const float *wts) const
415 {
416 #ifdef DEBUGM
417  float vals[8];
418  vals[0] = get_grid(i0, i1, i2, i3);
419  vals[1] = get_grid(i0, i1, i2+1, i3);
420  vals[2] = get_grid(i0, i1+1, i2, i3);
421  vals[3] = get_grid(i0, i1+1, i2+1, i3);
422  vals[4] = get_grid(i0+1, i1, i2, i3);
423  vals[5] = get_grid(i0+1, i1, i2+1, i3);
424  vals[6] = get_grid(i0+1, i1+1, i2, i3);
425  vals[7] = get_grid(i0+1, i1+1, i2+1, i3);
426 
427  switch (i3) {
428  case 0:
429  DebugM(2, "V\n" << endi);
430  break;
431  case 1:
432  DebugM(2, "dV/dx\n" << endi);
433  break;
434  case 2:
435  DebugM(2, "dV/dy\n" << endi);
436  break;
437  case 3:
438  DebugM(2, "dV/dz\n" << endi);
439  break;
440  }
441 
442  for (int i = 0; i < 8; i++) {
443  DebugM(2, "vals[" << i << "] = " << vals[i] << " wts[" << i << "] = " << wts[i] << "\n" << endi);
444  }
445 #endif
446 
447  float result =
448  wts[0] * get_grid(i0, i1, i2, i3) +
449  wts[1] * get_grid(i0, i1, i2+1, i3) +
450  wts[2] * get_grid(i0, i1+1, i2, i3) +
451  wts[3] * get_grid(i0, i1+1, i2+1, i3) +
452  wts[4] * get_grid(i0+1, i1, i2, i3) +
453  wts[5] * get_grid(i0+1, i1, i2+1, i3) +
454  wts[6] * get_grid(i0+1, i1+1, i2, i3) +
455  wts[7] * get_grid(i0+1, i1+1, i2+1, i3);
456 
457  DebugM(2, "result = " << result << "\n" << endi);
458 
459  return result;
460 }
461 
462 
463 inline Position GridforceGrid::wrap_position(const Position &pos, const Lattice &lattice)
464 {
465  // Wrap 'pos' about grid center, using periodic cell information in 'lattice'
466  // Position pos_wrapped = pos;
467  // Position center = get_center();
468  // pos_wrapped += lattice.wrap_delta(pos);
469  // pos_wrapped += lattice.delta(pos_wrapped, center) - (pos_wrapped - center);
470 
471  Position pos_wrapped = pos + lattice.wrap_delta(pos - get_center() + lattice.origin());
472 
473  return pos_wrapped;
474 }
475 
476 #endif
int compute_VdV(Position pos, float &V, Vector &dV) const
float compute_V(float *a, float *x, float *y, float *z) const
int get_inds(Position pos, int *inds, Vector &dg) const
static Tensor diagonal(const Vector &v1)
Definition: Tensor.h:37
float linear_interpolate(int i0, int i1, int i2, int i3, const float *wts) const
Definition: Vector.h:64
virtual void compute_b(float *b, int *inds, Vector gapscale) const =0
Vector wrap_delta(const Position &pos1) const
Definition: Lattice.h:206
Vector compute_dV(float *a, float *x, float *y, float *z) const
int compute_VdV(Position pos, float &V, Vector &dV) const
#define DebugM(x, y)
Definition: Debug.h:59
BigReal z
Definition: Vector.h:66
void compute_a(float *a, float *b) const
Vector origin() const
Definition: Lattice.h:262
void compute_wts(float *wts, const Vector &dg) const
float get_grid(int i0, int i1, int i2, int i3) const
float compute_d3V(float *a, float *x, float *y, float *z) const
GridforceFullSubGrid ** subgrids
virtual Position get_center(void) const =0
gridSize z
BigReal x
Definition: Vector.h:66
int get_inds(Position pos, int *inds, Vector &dg, Vector &gapscale) const
BigReal y
Definition: Vector.h:66
Position wrap_position(const Position &pos, const Lattice &lattice)
gridSize y
infostream & endi(infostream &s)
Definition: InfoStream.C:38
gridSize x
Vector compute_d2V(float *a, float *x, float *y, float *z) const
for(int i=0;i< n1;++i)