1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38 package ffx.potential.nonbonded.implicit;
39
40 import edu.rit.pj.IntegerForLoop;
41 import edu.rit.pj.ParallelRegion;
42 import ffx.crystal.Crystal;
43 import ffx.crystal.SymOp;
44 import ffx.numerics.atomic.AtomicDoubleArray3D;
45 import ffx.potential.bonded.Atom;
46
47 import java.util.List;
48 import java.util.logging.Level;
49 import java.util.logging.Logger;
50
51 import static ffx.numerics.multipole.GKSource.cn;
52 import static ffx.potential.parameters.MultipoleType.t000;
53 import static ffx.potential.parameters.MultipoleType.t001;
54 import static ffx.potential.parameters.MultipoleType.t002;
55 import static ffx.potential.parameters.MultipoleType.t010;
56 import static ffx.potential.parameters.MultipoleType.t011;
57 import static ffx.potential.parameters.MultipoleType.t020;
58 import static ffx.potential.parameters.MultipoleType.t100;
59 import static ffx.potential.parameters.MultipoleType.t101;
60 import static ffx.potential.parameters.MultipoleType.t110;
61 import static ffx.potential.parameters.MultipoleType.t200;
62 import static org.apache.commons.math3.util.FastMath.exp;
63 import static org.apache.commons.math3.util.FastMath.sqrt;
64
65
66
67
68
69
70
71 public class PermanentGKFieldRegion extends ParallelRegion {
72
73 private static final Logger logger = Logger.getLogger(PermanentGKFieldRegion.class.getName());
74
75 private static final double oneThird = 1.0 / 3.0;
76
77 private final double gkc;
78
79 private final double fc;
80
81 private final double fd;
82
83 private final double fq;
84
85 private final PermanentGKFieldLoop[] permanentGKFieldLoop;
86
87 protected Atom[] atoms;
88
89 private double[][][] globalMultipole;
90
91 private Crystal crystal;
92
93 private double[][][] sXYZ;
94
95 private int[][][] neighborLists;
96
97 private boolean[] use = null;
98
99 private double cut2;
100
101 private double[] born;
102
103 private AtomicDoubleArray3D sharedGKField;
104
105
106
107
108
109
110
111
112
113 public PermanentGKFieldRegion(int nt, double soluteDieletric, double solventDieletric, double gkc) {
114
115 fc = cn(0, soluteDieletric, solventDieletric);
116 fd = cn(1, soluteDieletric, solventDieletric);
117 fq = cn(2, soluteDieletric, solventDieletric);
118 this.gkc = gkc;
119
120 permanentGKFieldLoop = new PermanentGKFieldLoop[nt];
121 for (int i = 0; i < nt; i++) {
122 permanentGKFieldLoop[i] = new PermanentGKFieldLoop();
123 }
124 }
125
126 public void init(
127 Atom[] atoms,
128 double[][][] globalMultipole,
129 Crystal crystal,
130 double[][][] sXYZ,
131 int[][][] neighborLists,
132 boolean[] use,
133 double cut2,
134 double[] born,
135 AtomicDoubleArray3D sharedGKField) {
136 this.atoms = atoms;
137 this.globalMultipole = globalMultipole;
138 this.crystal = crystal;
139 this.sXYZ = sXYZ;
140 this.neighborLists = neighborLists;
141 this.use = use;
142 this.cut2 = cut2;
143 this.born = born;
144 this.sharedGKField = sharedGKField;
145 }
146
147 @Override
148 public void run() {
149 try {
150 int threadIndex = getThreadIndex();
151 int nAtoms = atoms.length;
152 execute(0, nAtoms - 1, permanentGKFieldLoop[threadIndex]);
153 } catch (Exception e) {
154 String message = "Fatal exception computing GK Energy in thread " + getThreadIndex() + "\n";
155 logger.log(Level.SEVERE, message, e);
156 }
157 }
158
159
160
161
162
163
164 private class PermanentGKFieldLoop extends IntegerForLoop {
165
166 private final double[][] a;
167 private final double[] gc;
168 private final double[] gux, guy, guz;
169 private final double[] gqxx, gqyy, gqzz;
170 private final double[] gqxy, gqxz, gqyz;
171 private final double[] dx_local;
172 private int threadID;
173 private double xi, yi, zi;
174 private double ci, uxi, uyi, uzi, qxxi, qxyi, qxzi, qyyi, qyzi, qzzi;
175 private double rbi;
176 private int iSymm;
177 private double[][] transOp;
178
179 PermanentGKFieldLoop() {
180 a = new double[4][3];
181 gc = new double[11];
182 gux = new double[11];
183 guy = new double[11];
184 guz = new double[11];
185 gqxx = new double[11];
186 gqyy = new double[11];
187 gqzz = new double[11];
188 gqxy = new double[11];
189 gqxz = new double[11];
190 gqyz = new double[11];
191 dx_local = new double[3];
192 transOp = new double[3][3];
193 }
194
195 @Override
196 public void run(int lb, int ub) {
197
198 threadID = getThreadIndex();
199
200 double[] x = sXYZ[0][0];
201 double[] y = sXYZ[0][1];
202 double[] z = sXYZ[0][2];
203
204 int nSymm = crystal.spaceGroup.symOps.size();
205 List<SymOp> symOps = crystal.spaceGroup.symOps;
206 for (iSymm = 0; iSymm < nSymm; iSymm++) {
207 SymOp symOp = symOps.get(iSymm);
208 crystal.getTransformationOperator(symOp, transOp);
209 for (int i = lb; i <= ub; i++) {
210 if (!use[i]) {
211 continue;
212 }
213 xi = x[i];
214 yi = y[i];
215 zi = z[i];
216 double[] multipolei = globalMultipole[0][i];
217 ci = multipolei[t000];
218 uxi = multipolei[t100];
219 uyi = multipolei[t010];
220 uzi = multipolei[t001];
221 qxxi = multipolei[t200] * oneThird;
222 qxyi = multipolei[t110] * oneThird;
223 qxzi = multipolei[t101] * oneThird;
224 qyyi = multipolei[t020] * oneThird;
225 qyzi = multipolei[t011] * oneThird;
226 qzzi = multipolei[t002] * oneThird;
227 rbi = born[i];
228 int[] list = neighborLists[iSymm][i];
229 for (int k : list) {
230 if (!use[k]) {
231 continue;
232 }
233 permanentGKField(i, k);
234 }
235
236
237 if (iSymm == 0) {
238 permanentGKField(i, i);
239 }
240 }
241 }
242 }
243
244 private void permanentGKField(int i, int k) {
245 dx_local[0] = sXYZ[iSymm][0][k] - xi;
246 dx_local[1] = sXYZ[iSymm][1][k] - yi;
247 dx_local[2] = sXYZ[iSymm][2][k] - zi;
248 final double r2 = crystal.image(dx_local);
249 if (r2 > cut2) {
250 return;
251 }
252 double xr = dx_local[0];
253 double yr = dx_local[1];
254 double zr = dx_local[2];
255 double xr2 = xr * xr;
256 double yr2 = yr * yr;
257 double zr2 = zr * zr;
258 final double rbk = born[k];
259 final double[] multipolek = globalMultipole[iSymm][k];
260 final double ck = multipolek[t000];
261 final double uxk = multipolek[t100];
262 final double uyk = multipolek[t010];
263 final double uzk = multipolek[t001];
264 final double qxxk = multipolek[t200] * oneThird;
265 final double qxyk = multipolek[t110] * oneThird;
266 final double qxzk = multipolek[t101] * oneThird;
267 final double qyyk = multipolek[t020] * oneThird;
268 final double qyzk = multipolek[t011] * oneThird;
269 final double qzzk = multipolek[t002] * oneThird;
270 final double rb2 = rbi * rbk;
271 final double expterm = exp(-r2 / (gkc * rb2));
272 final double expc = expterm / gkc;
273 final double expc1 = 1.0 - expc;
274 final double dexpc = -2.0 / (gkc * rb2);
275 final double expcdexpc = -expc * dexpc;
276 final double gf2 = 1.0 / (r2 + rb2 * expterm);
277 final double gf = sqrt(gf2);
278 final double gf3 = gf2 * gf;
279 final double gf5 = gf3 * gf2;
280 final double gf7 = gf5 * gf2;
281
282
283 a[0][0] = gf;
284 a[1][0] = -gf3;
285 a[2][0] = 3.0 * gf5;
286 a[3][0] = -15.0 * gf7;
287
288
289 a[0][1] = expc1 * a[1][0];
290 a[1][1] = expc1 * a[2][0];
291 a[2][1] = expc1 * a[3][0];
292
293 a[1][2] = expc1 * a[2][1] + expcdexpc * a[2][0];
294
295
296 a[0][1] = fc * a[0][1];
297 a[1][0] = fd * a[1][0];
298 a[1][1] = fd * a[1][1];
299 a[1][2] = fd * a[1][2];
300 a[2][0] = fq * a[2][0];
301 a[2][1] = fq * a[2][1];
302
303
304 gux[1] = xr * a[1][0];
305 guy[1] = yr * a[1][0];
306 guz[1] = zr * a[1][0];
307
308
309 gc[2] = xr * a[0][1];
310 gc[3] = yr * a[0][1];
311 gc[4] = zr * a[0][1];
312 gux[2] = a[1][0] + xr2 * a[1][1];
313 gux[3] = xr * yr * a[1][1];
314 gux[4] = xr * zr * a[1][1];
315 guy[2] = gux[3];
316 guy[3] = a[1][0] + yr2 * a[1][1];
317 guy[4] = yr * zr * a[1][1];
318 guz[2] = gux[4];
319 guz[3] = guy[4];
320 guz[4] = a[1][0] + zr2 * a[1][1];
321 gqxx[2] = xr * (2.0 * a[2][0] + xr2 * a[2][1]);
322 gqxx[3] = yr * xr2 * a[2][1];
323 gqxx[4] = zr * xr2 * a[2][1];
324 gqyy[2] = xr * yr2 * a[2][1];
325 gqyy[3] = yr * (2.0 * a[2][0] + yr2 * a[2][1]);
326 gqyy[4] = zr * yr2 * a[2][1];
327 gqzz[2] = xr * zr2 * a[2][1];
328 gqzz[3] = yr * zr2 * a[2][1];
329 gqzz[4] = zr * (2.0 * a[2][0] + zr2 * a[2][1]);
330 gqxy[2] = yr * (a[2][0] + xr2 * a[2][1]);
331 gqxy[3] = xr * (a[2][0] + yr2 * a[2][1]);
332 gqxy[4] = zr * xr * yr * a[2][1];
333 gqxz[2] = zr * (a[2][0] + xr2 * a[2][1]);
334 gqxz[3] = gqxy[4];
335 gqxz[4] = xr * (a[2][0] + zr2 * a[2][1]);
336 gqyz[2] = gqxy[4];
337 gqyz[3] = zr * (a[2][0] + yr2 * a[2][1]);
338 gqyz[4] = yr * (a[2][0] + zr2 * a[2][1]);
339
340
341 gux[5] = xr * (3.0 * a[1][1] + xr2 * a[1][2]);
342 gux[6] = yr * (a[1][1] + xr2 * a[1][2]);
343 gux[7] = zr * (a[1][1] + xr2 * a[1][2]);
344 gux[8] = xr * (a[1][1] + yr2 * a[1][2]);
345 gux[9] = zr * xr * yr * a[1][2];
346 gux[10] = xr * (a[1][1] + zr2 * a[1][2]);
347 guy[5] = yr * (a[1][1] + xr2 * a[1][2]);
348 guy[6] = xr * (a[1][1] + yr2 * a[1][2]);
349 guy[7] = gux[9];
350 guy[8] = yr * (3.0 * a[1][1] + yr2 * a[1][2]);
351 guy[9] = zr * (a[1][1] + yr2 * a[1][2]);
352 guy[10] = yr * (a[1][1] + zr2 * a[1][2]);
353 guz[5] = zr * (a[1][1] + xr2 * a[1][2]);
354 guz[6] = gux[9];
355 guz[7] = xr * (a[1][1] + zr2 * a[1][2]);
356 guz[8] = zr * (a[1][1] + yr2 * a[1][2]);
357 guz[9] = yr * (a[1][1] + zr2 * a[1][2]);
358 guz[10] = zr * (3.0 * a[1][1] + zr2 * a[1][2]);
359
360
361 double fix =
362 uxk * gux[2]
363 + uyk * gux[3]
364 + uzk * gux[4]
365 + 0.5
366 * (ck * gux[1]
367 + qxxk * gux[5]
368 + qyyk * gux[8]
369 + qzzk * gux[10]
370 + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9]))
371 + 0.5
372 * (ck * gc[2]
373 + qxxk * gqxx[2]
374 + qyyk * gqyy[2]
375 + qzzk * gqzz[2]
376 + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]));
377 double fiy =
378 uxk * guy[2]
379 + uyk * guy[3]
380 + uzk * guy[4]
381 + 0.5
382 * (ck * guy[1]
383 + qxxk * guy[5]
384 + qyyk * guy[8]
385 + qzzk * guy[10]
386 + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9]))
387 + 0.5
388 * (ck * gc[3]
389 + qxxk * gqxx[3]
390 + qyyk * gqyy[3]
391 + qzzk * gqzz[3]
392 + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]));
393 double fiz =
394 uxk * guz[2]
395 + uyk * guz[3]
396 + uzk * guz[4]
397 + 0.5
398 * (ck * guz[1]
399 + qxxk * guz[5]
400 + qyyk * guz[8]
401 + qzzk * guz[10]
402 + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9]))
403 + 0.5
404 * (ck * gc[4]
405 + qxxk * gqxx[4]
406 + qyyk * gqyy[4]
407 + qzzk * gqzz[4]
408 + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]));
409 double fkx =
410 uxi * gux[2]
411 + uyi * gux[3]
412 + uzi * gux[4]
413 - 0.5
414 * (ci * gux[1]
415 + qxxi * gux[5]
416 + qyyi * gux[8]
417 + qzzi * gux[10]
418 + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9]))
419 - 0.5
420 * (ci * gc[2]
421 + qxxi * gqxx[2]
422 + qyyi * gqyy[2]
423 + qzzi * gqzz[2]
424 + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]));
425 double fky =
426 uxi * guy[2]
427 + uyi * guy[3]
428 + uzi * guy[4]
429 - 0.5
430 * (ci * guy[1]
431 + qxxi * guy[5]
432 + qyyi * guy[8]
433 + qzzi * guy[10]
434 + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9]))
435 - 0.5
436 * (ci * gc[3]
437 + qxxi * gqxx[3]
438 + qyyi * gqyy[3]
439 + qzzi * gqzz[3]
440 + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]));
441 double fkz =
442 uxi * guz[2]
443 + uyi * guz[3]
444 + uzi * guz[4]
445 - 0.5
446 * (ci * guz[1]
447 + qxxi * guz[5]
448 + qyyi * guz[8]
449 + qzzi * guz[10]
450 + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9]))
451 - 0.5
452 * (ci * gc[4]
453 + qxxi * gqxx[4]
454 + qyyi * gqyy[4]
455 + qzzi * gqzz[4]
456 + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]));
457
458
459 if (i == k) {
460 fix *= 0.5;
461 fiy *= 0.5;
462 fiz *= 0.5;
463 fkx *= 0.5;
464 fky *= 0.5;
465 fkz *= 0.5;
466 }
467
468 double xc = fkx;
469 double yc = fky;
470 double zc = fkz;
471 fkx = xc * transOp[0][0] + yc * transOp[1][0] + zc * transOp[2][0];
472 fky = xc * transOp[0][1] + yc * transOp[1][1] + zc * transOp[2][1];
473 fkz = xc * transOp[0][2] + yc * transOp[1][2] + zc * transOp[2][2];
474
475 sharedGKField.add(threadID, i, fix, fiy, fiz);
476 sharedGKField.add(threadID, k, fkx, fky, fkz);
477 }
478 }
479 }