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