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.pme;
39
40 import static ffx.numerics.special.Erf.erfc;
41 import static org.apache.commons.math3.util.FastMath.exp;
42 import static org.apache.commons.math3.util.FastMath.min;
43 import static org.apache.commons.math3.util.FastMath.sqrt;
44
45 import edu.rit.pj.IntegerForLoop;
46 import edu.rit.pj.IntegerSchedule;
47 import edu.rit.pj.ParallelRegion;
48 import edu.rit.pj.ParallelSection;
49 import edu.rit.pj.ParallelTeam;
50 import ffx.crystal.Crystal;
51 import ffx.crystal.SymOp;
52 import ffx.numerics.atomic.AtomicDoubleArray3D;
53 import ffx.potential.bonded.Atom;
54 import ffx.potential.nonbonded.ReciprocalSpace;
55 import ffx.potential.parameters.ForceField;
56 import java.util.List;
57 import java.util.logging.Level;
58 import java.util.logging.Logger;
59
60
61
62
63
64
65
66
67
68
69
70
71 public class InducedDipoleFieldRegion extends ParallelRegion {
72
73 private static final Logger logger = Logger.getLogger(InducedDipoleFieldRegion.class.getName());
74
75 private final boolean intermolecularSoftcore;
76
77 private final boolean intramolecularSoftcore;
78
79 public double[][][] inducedDipole;
80
81 public double[][][] inducedDipoleCR;
82
83 private Atom[] atoms;
84
85 private Crystal crystal;
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104 private boolean[] use;
105
106 private int[] molecule;
107
108 private double[] ipdamp;
109 private double[] thole;
110
111 private double[][][] coordinates;
112
113
114
115 private int[][][] realSpaceLists;
116
117 private int[][] realSpaceCounts;
118
119 private IntegerSchedule realSpaceSchedule;
120
121 private ReciprocalSpace reciprocalSpace;
122
123 private boolean reciprocalSpaceTerm;
124
125 private LambdaMode lambdaMode = LambdaMode.OFF;
126
127 private double aewald;
128 private double an0, an1;
129 private long realSpaceSCFTotal;
130 private long[] realSpaceSCFTime;
131
132 private AtomicDoubleArray3D field;
133
134 private AtomicDoubleArray3D fieldCR;
135
136 private final InducedDipoleRealSpaceFieldSection inducedRealSpaceFieldSection;
137 private final InducedDipoleReciprocalFieldSection inducedReciprocalFieldSection;
138
139 public InducedDipoleFieldRegion(ParallelTeam parallelTeam,
140 ForceField forceField, boolean lambdaTerm) {
141 inducedRealSpaceFieldSection = new InducedDipoleRealSpaceFieldSection(parallelTeam);
142 inducedReciprocalFieldSection = new InducedDipoleReciprocalFieldSection();
143
144
145 if (lambdaTerm) {
146 intermolecularSoftcore = forceField.getBoolean("INTERMOLECULAR_SOFTCORE", false);
147 intramolecularSoftcore = forceField.getBoolean("INTRAMOLECULAR_SOFTCORE", false);
148 } else {
149 intermolecularSoftcore = false;
150 intramolecularSoftcore = false;
151 }
152 }
153
154
155
156
157
158
159 public void executeWith(ParallelTeam sectionTeam) {
160 try {
161 sectionTeam.execute(this);
162 } catch (RuntimeException e) {
163 String message = "Runtime exception computing the induced reciprocal space field.\n";
164 logger.log(Level.WARNING, message, e);
165 throw e;
166 } catch (Exception ex) {
167 String message = "Fatal exception computing the induced reciprocal space field.\n";
168 logger.log(Level.SEVERE, message, ex);
169 }
170 }
171
172 public long getRealSpaceSCFTotal() {
173 return realSpaceSCFTotal;
174 }
175
176 public void init(
177 Atom[] atoms,
178 Crystal crystal,
179 boolean[] use,
180 int[] molecule,
181 double[] ipdamp,
182 double[] thole,
183 double[][][] coordinates,
184 RealSpaceNeighborParameters realSpaceNeighborParameters,
185 double[][][] inducedDipole,
186 double[][][] inducedDipoleCR,
187 boolean reciprocalSpaceTerm,
188 ReciprocalSpace reciprocalSpace,
189 LambdaMode lambdaMode,
190 EwaldParameters ewaldParameters,
191 AtomicDoubleArray3D field,
192 AtomicDoubleArray3D fieldCR,
193 PMETimings pmeTimings) {
194
195 this.atoms = atoms;
196 this.crystal = crystal;
197 this.use = use;
198 this.molecule = molecule;
199 this.ipdamp = ipdamp;
200 this.thole = thole;
201 this.coordinates = coordinates;
202 this.realSpaceLists = realSpaceNeighborParameters.realSpaceLists;
203 this.realSpaceCounts = realSpaceNeighborParameters.realSpaceCounts;
204 this.realSpaceSchedule = realSpaceNeighborParameters.realSpaceSchedule;
205 this.inducedDipole = inducedDipole;
206 this.inducedDipoleCR = inducedDipoleCR;
207 this.reciprocalSpaceTerm = reciprocalSpaceTerm;
208 this.reciprocalSpace = reciprocalSpace;
209 this.lambdaMode = lambdaMode;
210 this.aewald = ewaldParameters.aewald;
211 this.an0 = ewaldParameters.an0;
212 this.an1 = ewaldParameters.an1;
213
214 this.field = field;
215 this.fieldCR = fieldCR;
216 this.realSpaceSCFTotal = pmeTimings.realSpaceSCFTotal;
217 this.realSpaceSCFTime = pmeTimings.realSpaceSCFTime;
218 }
219
220 @Override
221 public void run() {
222 try {
223 if (reciprocalSpaceTerm && aewald > 0.0) {
224 execute(inducedRealSpaceFieldSection, inducedReciprocalFieldSection);
225 } else {
226 execute(inducedRealSpaceFieldSection);
227 }
228 } catch (Exception e) {
229 String message = "Fatal exception computing the induced dipole field.\n";
230 logger.log(Level.SEVERE, message, e);
231 }
232 }
233
234 private class InducedDipoleRealSpaceFieldSection extends ParallelSection {
235
236 private final InducedDipoleRealSpaceFieldRegion polarizationRealSpaceFieldRegion;
237 private final ParallelTeam parallelTeam;
238
239 InducedDipoleRealSpaceFieldSection(ParallelTeam parallelTeam) {
240 this.parallelTeam = parallelTeam;
241 int nThreads = parallelTeam.getThreadCount();
242 polarizationRealSpaceFieldRegion = new InducedDipoleRealSpaceFieldRegion(nThreads);
243 }
244
245 @Override
246 public void run() {
247 try {
248 realSpaceSCFTotal -= System.nanoTime();
249 parallelTeam.execute(polarizationRealSpaceFieldRegion);
250 realSpaceSCFTotal += System.nanoTime();
251 } catch (Exception e) {
252 String message = "Fatal exception computing the real space field.\n";
253 logger.log(Level.SEVERE, message, e);
254 }
255 }
256 }
257
258 private class InducedDipoleReciprocalFieldSection extends ParallelSection {
259
260 @Override
261 public void run() {
262 reciprocalSpace.performConvolution();
263 }
264 }
265
266 private class InducedDipoleRealSpaceFieldRegion extends ParallelRegion {
267
268 private final InducedRealSpaceFieldLoop[] inducedRealSpaceFieldLoop;
269
270 InducedDipoleRealSpaceFieldRegion(int threadCount) {
271 inducedRealSpaceFieldLoop = new InducedRealSpaceFieldLoop[threadCount];
272 }
273
274 @Override
275 public void run() {
276 int threadIndex = getThreadIndex();
277 if (inducedRealSpaceFieldLoop[threadIndex] == null) {
278 inducedRealSpaceFieldLoop[threadIndex] = new InducedRealSpaceFieldLoop();
279 }
280 try {
281 int nAtoms = atoms.length;
282 execute(0, nAtoms - 1, inducedRealSpaceFieldLoop[threadIndex]);
283 } catch (Exception e) {
284 String message = "Fatal exception computing the induced real space field in thread " + getThreadIndex() + "\n";
285 logger.log(Level.SEVERE, message, e);
286 }
287 }
288
289 private class InducedRealSpaceFieldLoop extends IntegerForLoop {
290
291 private int threadID;
292 private double[][] ind, indCR;
293 private double[] x, y, z;
294
295 InducedRealSpaceFieldLoop() {}
296
297 @Override
298 public void finish() {
299 int threadIndex = getThreadIndex();
300 realSpaceSCFTime[threadIndex] += System.nanoTime();
301 }
302
303 @Override
304 public void run(int lb, int ub) {
305 final double[] dx = new double[3];
306 final double[][] transOp = new double[3][3];
307
308
309 int[][] lists = realSpaceLists[0];
310 int[] counts = realSpaceCounts[0];
311 for (int i = lb; i <= ub; i++) {
312 if (!use[i]) {
313 continue;
314 }
315 final int moleculei = molecule[i];
316 double fx = 0.0;
317 double fy = 0.0;
318 double fz = 0.0;
319 double px = 0.0;
320 double py = 0.0;
321 double pz = 0.0;
322 final double xi = x[i];
323 final double yi = y[i];
324 final double zi = z[i];
325 final double[] dipolei = ind[i];
326 final double uix = dipolei[0];
327 final double uiy = dipolei[1];
328 final double uiz = dipolei[2];
329 final double[] dipoleCRi = indCR[i];
330 final double pix = dipoleCRi[0];
331 final double piy = dipoleCRi[1];
332 final double piz = dipoleCRi[2];
333 final double pdi = ipdamp[i];
334 final double pti = thole[i];
335
336
337 final int[] list = lists[i];
338 final int npair = counts[i];
339 for (int j = 0; j < npair; j++) {
340 final int k = list[j];
341 if (!use[k]) {
342 continue;
343 }
344 boolean sameMolecule = (moleculei == molecule[k]);
345 if (lambdaMode == LambdaMode.VAPOR) {
346 if ((intermolecularSoftcore && !sameMolecule)
347 || (intramolecularSoftcore && sameMolecule)) {
348 continue;
349 }
350 }
351 final double pdk = ipdamp[k];
352 final double ptk = thole[k];
353 dx[0] = x[k] - xi;
354 dx[1] = y[k] - yi;
355 dx[2] = z[k] - zi;
356 final double r2 = crystal.image(dx);
357
358
359 final double r = sqrt(r2);
360 final double rr1 = 1.0 / r;
361 final double rr2 = rr1 * rr1;
362 final double ralpha = aewald * r;
363 final double exp2a = exp(-ralpha * ralpha);
364 final double bn0 = erfc(ralpha) * rr1;
365 final double bn1 = (bn0 + an0 * exp2a) * rr2;
366 final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
367 double scale3 = 1.0;
368 double scale5 = 1.0;
369 double damp = pdi * pdk;
370 final double pgamma = min(pti, ptk);
371 final double rdamp = r * damp;
372 damp = -pgamma * rdamp * rdamp * rdamp;
373 if (damp > -50.0) {
374 final double expdamp = exp(damp);
375 scale3 = 1.0 - expdamp;
376 scale5 = 1.0 - expdamp * (1.0 - damp);
377 }
378 double rr3 = rr1 * rr2;
379 double rr5 = 3.0 * rr3 * rr2;
380 rr3 *= (1.0 - scale3);
381 rr5 *= (1.0 - scale5);
382 final double xr = dx[0];
383 final double yr = dx[1];
384 final double zr = dx[2];
385 final double[] dipolek = ind[k];
386 final double ukx = dipolek[0];
387 final double uky = dipolek[1];
388 final double ukz = dipolek[2];
389 final double ukr = ukx * xr + uky * yr + ukz * zr;
390 final double bn2ukr = bn2 * ukr;
391 final double fimx = -bn1 * ukx + bn2ukr * xr;
392 final double fimy = -bn1 * uky + bn2ukr * yr;
393 final double fimz = -bn1 * ukz + bn2ukr * zr;
394 final double rr5ukr = rr5 * ukr;
395 final double fidx = -rr3 * ukx + rr5ukr * xr;
396 final double fidy = -rr3 * uky + rr5ukr * yr;
397 final double fidz = -rr3 * ukz + rr5ukr * zr;
398 fx += (fimx - fidx);
399 fy += (fimy - fidy);
400 fz += (fimz - fidz);
401 final double[] dipolepk = indCR[k];
402 final double pkx = dipolepk[0];
403 final double pky = dipolepk[1];
404 final double pkz = dipolepk[2];
405 final double pkr = pkx * xr + pky * yr + pkz * zr;
406 final double bn2pkr = bn2 * pkr;
407 final double pimx = -bn1 * pkx + bn2pkr * xr;
408 final double pimy = -bn1 * pky + bn2pkr * yr;
409 final double pimz = -bn1 * pkz + bn2pkr * zr;
410 final double rr5pkr = rr5 * pkr;
411 final double pidx = -rr3 * pkx + rr5pkr * xr;
412 final double pidy = -rr3 * pky + rr5pkr * yr;
413 final double pidz = -rr3 * pkz + rr5pkr * zr;
414 px += (pimx - pidx);
415 py += (pimy - pidy);
416 pz += (pimz - pidz);
417 final double uir = uix * xr + uiy * yr + uiz * zr;
418 final double bn2uir = bn2 * uir;
419 final double fkmx = -bn1 * uix + bn2uir * xr;
420 final double fkmy = -bn1 * uiy + bn2uir * yr;
421 final double fkmz = -bn1 * uiz + bn2uir * zr;
422 final double rr5uir = rr5 * uir;
423 final double fkdx = -rr3 * uix + rr5uir * xr;
424 final double fkdy = -rr3 * uiy + rr5uir * yr;
425 final double fkdz = -rr3 * uiz + rr5uir * zr;
426 field.add(threadID, k, fkmx - fkdx, fkmy - fkdy, fkmz - fkdz);
427 final double pir = pix * xr + piy * yr + piz * zr;
428 final double bn2pir = bn2 * pir;
429 final double pkmx = -bn1 * pix + bn2pir * xr;
430 final double pkmy = -bn1 * piy + bn2pir * yr;
431 final double pkmz = -bn1 * piz + bn2pir * zr;
432 final double rr5pir = rr5 * pir;
433 final double pkdx = -rr3 * pix + rr5pir * xr;
434 final double pkdy = -rr3 * piy + rr5pir * yr;
435 final double pkdz = -rr3 * piz + rr5pir * zr;
436 fieldCR.add(threadID, k, pkmx - pkdx, pkmy - pkdy, pkmz - pkdz);
437 }
438 field.add(threadID, i, fx, fy, fz);
439 fieldCR.add(threadID, i, px, py, pz);
440 }
441
442
443 List<SymOp> symOps = crystal.spaceGroup.symOps;
444 int nSymm = symOps.size();
445 for (int iSymm = 1; iSymm < nSymm; iSymm++) {
446 SymOp symOp = crystal.spaceGroup.getSymOp(iSymm);
447 crystal.getTransformationOperator(symOp, transOp);
448 lists = realSpaceLists[iSymm];
449 counts = realSpaceCounts[iSymm];
450 final double[] xs = coordinates[iSymm][0];
451 final double[] ys = coordinates[iSymm][1];
452 final double[] zs = coordinates[iSymm][2];
453 final double[][] inds = inducedDipole[iSymm];
454 final double[][] indCRs = inducedDipoleCR[iSymm];
455
456
457 for (int i = lb; i <= ub; i++) {
458 if (!use[i]) {
459 continue;
460 }
461 double fx = 0.0;
462 double fy = 0.0;
463 double fz = 0.0;
464 double px = 0.0;
465 double py = 0.0;
466 double pz = 0.0;
467 final double xi = x[i];
468 final double yi = y[i];
469 final double zi = z[i];
470 final double[] dipolei = ind[i];
471 final double uix = dipolei[0];
472 final double uiy = dipolei[1];
473 final double uiz = dipolei[2];
474 final double[] dipoleCRi = indCR[i];
475 final double pix = dipoleCRi[0];
476 final double piy = dipoleCRi[1];
477 final double piz = dipoleCRi[2];
478 final double pdi = ipdamp[i];
479 final double pti = thole[i];
480
481
482 final int[] list = lists[i];
483 final int npair = counts[i];
484 for (int j = 0; j < npair; j++) {
485 final int k = list[j];
486 if (!use[k]) {
487 continue;
488 }
489 double selfScale = 1.0;
490 if (i == k) {
491 selfScale = 0.5;
492 }
493 final double pdk = ipdamp[k];
494 final double ptk = thole[k];
495 dx[0] = xs[k] - xi;
496 dx[1] = ys[k] - yi;
497 dx[2] = zs[k] - zi;
498 final double r2 = crystal.image(dx);
499
500
501 final double r = sqrt(r2);
502 final double rr1 = 1.0 / r;
503 final double rr2 = rr1 * rr1;
504 final double ralpha = aewald * r;
505 final double exp2a = exp(-ralpha * ralpha);
506 final double bn0 = erfc(ralpha) * rr1;
507 final double bn1 = (bn0 + an0 * exp2a) * rr2;
508 final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
509 double scale3 = 1.0;
510 double scale5 = 1.0;
511 double damp = pdi * pdk;
512 final double pgamma = min(pti, ptk);
513 final double rdamp = r * damp;
514 damp = -pgamma * rdamp * rdamp * rdamp;
515 if (damp > -50.0) {
516 final double expdamp = exp(damp);
517 scale3 = 1.0 - expdamp;
518 scale5 = 1.0 - expdamp * (1.0 - damp);
519 }
520 double rr3 = rr1 * rr2;
521 double rr5 = 3.0 * rr3 * rr2;
522 rr3 *= (1.0 - scale3);
523 rr5 *= (1.0 - scale5);
524 final double xr = dx[0];
525 final double yr = dx[1];
526 final double zr = dx[2];
527 final double[] dipolek = inds[k];
528 final double ukx = dipolek[0];
529 final double uky = dipolek[1];
530 final double ukz = dipolek[2];
531 final double[] dipolepk = indCRs[k];
532 final double pkx = dipolepk[0];
533 final double pky = dipolepk[1];
534 final double pkz = dipolepk[2];
535 final double ukr = ukx * xr + uky * yr + ukz * zr;
536 final double bn2ukr = bn2 * ukr;
537 final double fimx = -bn1 * ukx + bn2ukr * xr;
538 final double fimy = -bn1 * uky + bn2ukr * yr;
539 final double fimz = -bn1 * ukz + bn2ukr * zr;
540 final double rr5ukr = rr5 * ukr;
541 final double fidx = -rr3 * ukx + rr5ukr * xr;
542 final double fidy = -rr3 * uky + rr5ukr * yr;
543 final double fidz = -rr3 * ukz + rr5ukr * zr;
544 fx += selfScale * (fimx - fidx);
545 fy += selfScale * (fimy - fidy);
546 fz += selfScale * (fimz - fidz);
547 final double pkr = pkx * xr + pky * yr + pkz * zr;
548 final double bn2pkr = bn2 * pkr;
549 final double pimx = -bn1 * pkx + bn2pkr * xr;
550 final double pimy = -bn1 * pky + bn2pkr * yr;
551 final double pimz = -bn1 * pkz + bn2pkr * zr;
552 final double rr5pkr = rr5 * pkr;
553 final double pidx = -rr3 * pkx + rr5pkr * xr;
554 final double pidy = -rr3 * pky + rr5pkr * yr;
555 final double pidz = -rr3 * pkz + rr5pkr * zr;
556 px += selfScale * (pimx - pidx);
557 py += selfScale * (pimy - pidy);
558 pz += selfScale * (pimz - pidz);
559 final double uir = uix * xr + uiy * yr + uiz * zr;
560 final double bn2uir = bn2 * uir;
561 final double fkmx = -bn1 * uix + bn2uir * xr;
562 final double fkmy = -bn1 * uiy + bn2uir * yr;
563 final double fkmz = -bn1 * uiz + bn2uir * zr;
564 final double rr5uir = rr5 * uir;
565 final double fkdx = -rr3 * uix + rr5uir * xr;
566 final double fkdy = -rr3 * uiy + rr5uir * yr;
567 final double fkdz = -rr3 * uiz + rr5uir * zr;
568 double xc = selfScale * (fkmx - fkdx);
569 double yc = selfScale * (fkmy - fkdy);
570 double zc = selfScale * (fkmz - fkdz);
571 double kx = (xc * transOp[0][0] + yc * transOp[1][0] + zc * transOp[2][0]);
572 double ky = (xc * transOp[0][1] + yc * transOp[1][1] + zc * transOp[2][1]);
573 double kz = (xc * transOp[0][2] + yc * transOp[1][2] + zc * transOp[2][2]);
574 field.add(threadID, k, kx, ky, kz);
575 final double pir = pix * xr + piy * yr + piz * zr;
576 final double bn2pir = bn2 * pir;
577 final double pkmx = -bn1 * pix + bn2pir * xr;
578 final double pkmy = -bn1 * piy + bn2pir * yr;
579 final double pkmz = -bn1 * piz + bn2pir * zr;
580 final double rr5pir = rr5 * pir;
581 final double pkdx = -rr3 * pix + rr5pir * xr;
582 final double pkdy = -rr3 * piy + rr5pir * yr;
583 final double pkdz = -rr3 * piz + rr5pir * zr;
584 xc = selfScale * (pkmx - pkdx);
585 yc = selfScale * (pkmy - pkdy);
586 zc = selfScale * (pkmz - pkdz);
587 kx = (xc * transOp[0][0] + yc * transOp[1][0] + zc * transOp[2][0]);
588 ky = (xc * transOp[0][1] + yc * transOp[1][1] + zc * transOp[2][1]);
589 kz = (xc * transOp[0][2] + yc * transOp[1][2] + zc * transOp[2][2]);
590 fieldCR.add(threadID, k, kx, ky, kz);
591 }
592 field.add(threadID, i, fx, fy, fz);
593 fieldCR.add(threadID, i, px, py, pz);
594 }
595 }
596 }
597
598 @Override
599 public IntegerSchedule schedule() {
600 return realSpaceSchedule;
601 }
602
603 @Override
604 public void start() {
605 threadID = getThreadIndex();
606 realSpaceSCFTime[threadID] -= System.nanoTime();
607 x = coordinates[0][0];
608 y = coordinates[0][1];
609 z = coordinates[0][2];
610 ind = inducedDipole[0];
611 indCR = inducedDipoleCR[0];
612 }
613 }
614 }
615 }