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