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.potential.parameters.MultipoleType.getRotationMatrix;
41 import static ffx.potential.parameters.MultipoleType.rotateMultipole;
42 import static ffx.potential.parameters.MultipoleType.t000;
43 import static ffx.potential.parameters.MultipoleType.t001;
44 import static ffx.potential.parameters.MultipoleType.t002;
45 import static ffx.potential.parameters.MultipoleType.t010;
46 import static ffx.potential.parameters.MultipoleType.t011;
47 import static ffx.potential.parameters.MultipoleType.t020;
48 import static ffx.potential.parameters.MultipoleType.t100;
49 import static ffx.potential.parameters.MultipoleType.t101;
50 import static ffx.potential.parameters.MultipoleType.t110;
51 import static ffx.potential.parameters.MultipoleType.t200;
52 import static org.apache.commons.math3.util.FastMath.max;
53
54 import edu.rit.pj.IntegerForLoop;
55 import edu.rit.pj.IntegerSchedule;
56 import edu.rit.pj.ParallelRegion;
57 import edu.rit.pj.ParallelTeam;
58 import ffx.crystal.Crystal;
59 import ffx.crystal.SymOp;
60 import ffx.numerics.atomic.AtomicDoubleArray3D;
61 import ffx.potential.bonded.Atom;
62 import ffx.potential.extended.ExtendedSystem;
63 import ffx.potential.nonbonded.ParticleMeshEwald;
64 import ffx.potential.parameters.ForceField;
65 import ffx.potential.parameters.MultipoleType;
66 import ffx.potential.parameters.MultipoleType.MultipoleFrameDefinition;
67 import ffx.potential.parameters.PolarizeType;
68 import java.util.List;
69 import java.util.logging.Level;
70 import java.util.logging.Logger;
71
72
73
74
75
76
77
78
79 public class InitializationRegion extends ParallelRegion {
80
81 private static final Logger logger = Logger.getLogger(InitializationRegion.class.getName());
82
83 private final ParticleMeshEwald particleMeshEwald;
84
85
86
87
88
89 private final boolean rotateMultipoles;
90
91 private final boolean useCharges;
92
93 private final boolean useDipoles;
94
95 private final boolean useQuadrupoles;
96
97
98
99 private final InitializationLoop[] initializationLoop;
100
101
102
103 private final RotateMultipolesLoop[] rotateMultipolesLoop;
104
105
106
107
108 private boolean lambdaTerm;
109
110 private ExtendedSystem esvSystem;
111
112
113
114 private boolean esvTerm;
115
116
117 private Atom[] atoms;
118
119 private double[][][] coordinates;
120
121 private Crystal crystal;
122
123 private MultipoleFrameDefinition[] frame;
124
125 private int[][] axisAtom;
126
127 private double[][][] globalMultipole;
128
129 private double[][][] titrationMultipole;
130 private double[][][] tautomerMultipole;
131
132 private double[] polarizability;
133 private double[] titrationPolarizability;
134 private double[] tautomerPolarizability;
135 private double[] thole;
136 private double[] ipdamp;
137
138
139
140
141 private boolean[] use;
142
143
144
145 private int[][][] neighborLists;
146
147
148
149 private int[][][] realSpaceLists;
150
151 private int[][][] vaporLists;
152
153 private AtomicDoubleArray3D grad;
154
155 private AtomicDoubleArray3D torque;
156
157 private AtomicDoubleArray3D lambdaGrad;
158
159 private AtomicDoubleArray3D lambdaTorque;
160
161 public InitializationRegion(ParticleMeshEwald particleMeshEwald, int maxThreads,
162 ForceField forceField) {
163 initializationLoop = new InitializationLoop[maxThreads];
164 rotateMultipolesLoop = new RotateMultipolesLoop[maxThreads];
165 useCharges = forceField.getBoolean("USE_CHARGES", true);
166 useDipoles = forceField.getBoolean("USE_DIPOLES", true);
167 useQuadrupoles = forceField.getBoolean("USE_QUADRUPOLES", true);
168 rotateMultipoles = forceField.getBoolean("ROTATE_MULTIPOLES", true);
169 this.particleMeshEwald = particleMeshEwald;
170 }
171
172
173
174
175
176
177 public void executeWith(ParallelTeam parallelTeam) {
178 try {
179 parallelTeam.execute(this);
180 } catch (RuntimeException e) {
181 String message = "RuntimeException expanding coordinates and rotating multipoles.\n";
182 logger.log(Level.WARNING, message, e);
183 throw e;
184 } catch (Exception e) {
185 String message = "Fatal exception expanding coordinates and rotating multipoles.\n";
186 logger.log(Level.SEVERE, message, e);
187 }
188 }
189
190 public void init(
191 boolean lambdaTerm,
192 ExtendedSystem esvSystem,
193 Atom[] atoms,
194 double[][][] coordinates,
195 Crystal crystal,
196 MultipoleFrameDefinition[] frame,
197 int[][] axisAtom,
198 double[][][] globalMultipole,
199 double[][][] titrationMultipole,
200 double[][][] tautomerMultipole,
201 double[] polarizability,
202 double[] titrationPolarizability,
203 double[] tautomerPolarizability,
204 double[] thole,
205 double[] ipdamp,
206 boolean[] use,
207 int[][][] neighborLists,
208 int[][][] realSpaceLists,
209 int[][][] vaporLists,
210 AtomicDoubleArray3D grad,
211 AtomicDoubleArray3D torque,
212 AtomicDoubleArray3D lambdaGrad,
213 AtomicDoubleArray3D lambdaTorque) {
214 this.lambdaTerm = lambdaTerm;
215 this.esvSystem = esvSystem;
216 if (esvSystem != null) {
217 this.esvTerm = true;
218 }
219 this.atoms = atoms;
220 this.coordinates = coordinates;
221 this.crystal = crystal;
222 this.frame = frame;
223 this.axisAtom = axisAtom;
224 this.globalMultipole = globalMultipole;
225 this.titrationMultipole = titrationMultipole;
226 this.tautomerMultipole = tautomerMultipole;
227 this.polarizability = polarizability;
228 this.titrationPolarizability = titrationPolarizability;
229 this.tautomerPolarizability = tautomerPolarizability;
230 this.thole = thole;
231 this.ipdamp = ipdamp;
232 this.use = use;
233 this.neighborLists = neighborLists;
234 this.realSpaceLists = realSpaceLists;
235 this.vaporLists = vaporLists;
236 this.grad = grad;
237 this.torque = torque;
238 this.lambdaGrad = lambdaGrad;
239 this.lambdaTorque = lambdaTorque;
240 }
241
242 @Override
243 public void run() {
244
245 int nAtoms = atoms.length;
246 int threadIndex = getThreadIndex();
247 if (initializationLoop[threadIndex] == null) {
248 initializationLoop[threadIndex] = new InitializationLoop();
249 rotateMultipolesLoop[threadIndex] = new RotateMultipolesLoop();
250 }
251 try {
252 execute(0, nAtoms - 1, initializationLoop[threadIndex]);
253 execute(0, nAtoms - 1, rotateMultipolesLoop[threadIndex]);
254 } catch (Exception e) {
255 String message = "Fatal exception initializing coordinates in thread: " + threadIndex + "\n";
256 logger.log(Level.SEVERE, message, e);
257 }
258 }
259
260 private class InitializationLoop extends IntegerForLoop {
261
262 private final double[] in = new double[3];
263 private final double[] out = new double[3];
264 private double[] x;
265 private double[] y;
266 private double[] z;
267 private int threadID;
268
269 @Override
270 public void run(int lb, int ub) {
271 grad.reset(threadID, lb, ub);
272 torque.reset(threadID, lb, ub);
273 if (lambdaTerm) {
274 lambdaGrad.reset(threadID, lb, ub);
275 lambdaTorque.reset(threadID, lb, ub);
276 }
277
278
279 for (int i = lb; i <= ub; i++) {
280 Atom atom = atoms[i];
281 x[i] = atom.getX();
282 y[i] = atom.getY();
283 z[i] = atom.getZ();
284 use[i] = atom.getUse();
285
286
287
288
289
290
291 int size = neighborLists[0][i].length;
292 if (vaporLists != null) {
293 size = max(size, vaporLists[0][i].length);
294 }
295 if (realSpaceLists[0][i] == null || realSpaceLists[0][i].length < size) {
296 realSpaceLists[0][i] = new int[size];
297 }
298 }
299
300
301 List<SymOp> symOps = crystal.spaceGroup.symOps;
302 int nSymm = symOps.size();
303 for (int iSymm = 1; iSymm < nSymm; iSymm++) {
304 SymOp symOp = symOps.get(iSymm);
305 double[] xs = coordinates[iSymm][0];
306 double[] ys = coordinates[iSymm][1];
307 double[] zs = coordinates[iSymm][2];
308 for (int i = lb; i <= ub; i++) {
309 in[0] = x[i];
310 in[1] = y[i];
311 in[2] = z[i];
312 crystal.applySymOp(in, out, symOp);
313 xs[i] = out[0];
314 ys[i] = out[1];
315 zs[i] = out[2];
316 int size = neighborLists[iSymm][i].length;
317 if (realSpaceLists[iSymm][i] == null || realSpaceLists[iSymm][i].length < size) {
318 realSpaceLists[iSymm][i] = new int[size];
319 }
320 }
321 }
322 }
323
324 @Override
325 public IntegerSchedule schedule() {
326 return IntegerSchedule.fixed();
327 }
328
329 @Override
330 public void start() {
331 x = coordinates[0][0];
332 y = coordinates[0][1];
333 z = coordinates[0][2];
334 threadID = getThreadIndex();
335 }
336 }
337
338 private class RotateMultipolesLoop extends IntegerForLoop {
339
340
341 private final double[] localOrigin = new double[3];
342 private final double[][] frameCoords = new double[4][3];
343 private final double[][] rotmat = new double[3][3];
344 private final double[] tempDipole = new double[3];
345 private final double[][] tempQuadrupole = new double[3][3];
346 private final double[] dipole = new double[3];
347 private final double[][] quadrupole = new double[3][3];
348
349 @Override
350 public void run(int lb, int ub) {
351 List<SymOp> symOps = crystal.spaceGroup.symOps;
352 int nSymm = symOps.size();
353 for (int iSymm = 0; iSymm < nSymm; iSymm++) {
354 final double[] x = coordinates[iSymm][0];
355 final double[] y = coordinates[iSymm][1];
356 final double[] z = coordinates[iSymm][2];
357 for (int ii = lb; ii <= ub; ii++) {
358 Atom atom = atoms[ii];
359 double chargeScale = 1.0;
360 double dipoleScale = 1.0;
361 double quadrupoleScale = 1.0;
362 double polarizabilityScale = 1.0;
363 if (!useCharges) {
364 chargeScale = 0.0;
365 }
366 if (!useDipoles) {
367 dipoleScale = 0.0;
368 }
369 if (!useQuadrupoles) {
370 quadrupoleScale = 0.0;
371 }
372
373 double elecScale = 1.0;
374 if (!atom.getElectrostatics()) {
375 elecScale = 0.0;
376 }
377
378
379 MultipoleType multipoleType = particleMeshEwald.getMultipoleType(ii);
380 double[] in = multipoleType.getMultipole();
381
382 frame[ii] = multipoleType.frameDefinition;
383
384 axisAtom[ii] = atom.getAxisAtomIndices();
385
386 if (rotateMultipoles) {
387
388 localOrigin[0] = x[ii];
389 localOrigin[1] = y[ii];
390 localOrigin[2] = z[ii];
391
392
393 int[] referenceSites = axisAtom[ii];
394 int nSites = 0;
395 if (referenceSites != null) {
396 nSites = referenceSites.length;
397 }
398 for (int i = 0; i < nSites; i++) {
399 int index = referenceSites[i];
400 frameCoords[i][0] = x[index];
401 frameCoords[i][1] = y[index];
402 frameCoords[i][2] = z[index];
403 }
404
405
406 tempDipole[0] = in[t100];
407 tempDipole[1] = in[t010];
408 tempDipole[2] = in[t001];
409
410 tempQuadrupole[0][0] = in[t200];
411 tempQuadrupole[1][1] = in[t020];
412 tempQuadrupole[2][2] = in[t002];
413 tempQuadrupole[0][1] = in[t110];
414 tempQuadrupole[0][2] = in[t101];
415 tempQuadrupole[1][2] = in[t011];
416 tempQuadrupole[1][0] = in[t110];
417 tempQuadrupole[2][0] = in[t101];
418 tempQuadrupole[2][1] = in[t011];
419
420
421 boolean needsChiralInversion = false;
422
423
424
425
426
427
428
429
430
431
432
433
434
435 getRotationMatrix(frame[ii], localOrigin, frameCoords, rotmat);
436 rotateMultipole(rotmat, tempDipole, tempQuadrupole, dipole, quadrupole);
437
438 double[] out = globalMultipole[iSymm][ii];
439
440 out[t000] = in[0] * chargeScale * elecScale;
441
442 out[t100] = dipole[0] * dipoleScale * elecScale;
443 out[t010] = dipole[1] * dipoleScale * elecScale;
444 out[t001] = dipole[2] * dipoleScale * elecScale;
445
446 out[t200] = quadrupole[0][0] * quadrupoleScale * elecScale;
447 out[t020] = quadrupole[1][1] * quadrupoleScale * elecScale;
448 out[t002] = quadrupole[2][2] * quadrupoleScale * elecScale;
449 out[t110] = quadrupole[0][1] * quadrupoleScale * elecScale;
450 out[t101] = quadrupole[0][2] * quadrupoleScale * elecScale;
451 out[t011] = quadrupole[1][2] * quadrupoleScale * elecScale;
452
453 if (esvTerm && esvSystem.isTitrating(ii)) {
454 double[] esvMultipoleTitrDot = new double[10];
455 double[] esvMultipoleTautDot = new double[10];
456 double titrationLambda = esvSystem.getTitrationLambda(ii);
457 double tautomerLambda = esvSystem.getTautomerLambda(ii);
458 esvMultipoleTitrDot = esvSystem.getTitrationUtils()
459 .getMultipoleTitrationDeriv(atom, titrationLambda,
460 tautomerLambda, esvMultipoleTitrDot);
461 esvMultipoleTautDot = esvSystem.getTitrationUtils()
462 .getMultipoleTautomerDeriv(atom, titrationLambda,
463 tautomerLambda, esvMultipoleTautDot);
464 double tempCharge = esvMultipoleTitrDot[0];
465
466 tempDipole[0] = esvMultipoleTitrDot[t100];
467 tempDipole[1] = esvMultipoleTitrDot[t010];
468 tempDipole[2] = esvMultipoleTitrDot[t001];
469
470 tempQuadrupole[0][0] = esvMultipoleTitrDot[t200];
471 tempQuadrupole[1][1] = esvMultipoleTitrDot[t020];
472 tempQuadrupole[2][2] = esvMultipoleTitrDot[t002];
473 tempQuadrupole[0][1] = esvMultipoleTitrDot[t110];
474 tempQuadrupole[0][2] = esvMultipoleTitrDot[t101];
475 tempQuadrupole[1][2] = esvMultipoleTitrDot[t011];
476 tempQuadrupole[1][0] = esvMultipoleTitrDot[t110];
477 tempQuadrupole[2][0] = esvMultipoleTitrDot[t101];
478 tempQuadrupole[2][1] = esvMultipoleTitrDot[t011];
479 if (needsChiralInversion) {
480 tempDipole[1] = -tempDipole[1];
481 tempQuadrupole[0][1] = -tempQuadrupole[0][1];
482 tempQuadrupole[1][0] = -tempQuadrupole[1][0];
483 tempQuadrupole[1][2] = -tempQuadrupole[1][2];
484 tempQuadrupole[2][1] = -tempQuadrupole[2][1];
485 }
486 rotateMultipole(rotmat, tempDipole, tempQuadrupole, dipole, quadrupole);
487 out = titrationMultipole[iSymm][ii];
488 out[t000] = tempCharge * chargeScale * elecScale;
489
490 out[t100] = dipole[0] * dipoleScale * elecScale;
491 out[t010] = dipole[1] * dipoleScale * elecScale;
492 out[t001] = dipole[2] * dipoleScale * elecScale;
493
494 out[t200] = quadrupole[0][0] * quadrupoleScale * elecScale;
495 out[t020] = quadrupole[1][1] * quadrupoleScale * elecScale;
496 out[t002] = quadrupole[2][2] * quadrupoleScale * elecScale;
497 out[t110] = quadrupole[0][1] * quadrupoleScale * elecScale;
498 out[t101] = quadrupole[0][2] * quadrupoleScale * elecScale;
499 out[t011] = quadrupole[1][2] * quadrupoleScale * elecScale;
500
501 tempCharge = esvMultipoleTautDot[0];
502
503 tempDipole[0] = esvMultipoleTautDot[t100];
504 tempDipole[1] = esvMultipoleTautDot[t010];
505 tempDipole[2] = esvMultipoleTautDot[t001];
506
507 tempQuadrupole[0][0] = esvMultipoleTautDot[t200];
508 tempQuadrupole[1][1] = esvMultipoleTautDot[t020];
509 tempQuadrupole[2][2] = esvMultipoleTautDot[t002];
510 tempQuadrupole[0][1] = esvMultipoleTautDot[t110];
511 tempQuadrupole[0][2] = esvMultipoleTautDot[t101];
512 tempQuadrupole[1][2] = esvMultipoleTautDot[t011];
513 tempQuadrupole[1][0] = esvMultipoleTautDot[t110];
514 tempQuadrupole[2][0] = esvMultipoleTautDot[t101];
515 tempQuadrupole[2][1] = esvMultipoleTautDot[t011];
516 if (needsChiralInversion) {
517 tempDipole[1] = -tempDipole[1];
518 tempQuadrupole[0][1] = -tempQuadrupole[0][1];
519 tempQuadrupole[1][0] = -tempQuadrupole[1][0];
520 tempQuadrupole[1][2] = -tempQuadrupole[1][2];
521 tempQuadrupole[2][1] = -tempQuadrupole[2][1];
522 }
523 rotateMultipole(rotmat, tempDipole, tempQuadrupole, dipole, quadrupole);
524 out = tautomerMultipole[iSymm][ii];
525 out[t000] = tempCharge * chargeScale * elecScale;
526
527 out[t100] = dipole[0] * dipoleScale * elecScale;
528 out[t010] = dipole[1] * dipoleScale * elecScale;
529 out[t001] = dipole[2] * dipoleScale * elecScale;
530
531 out[t200] = quadrupole[0][0] * quadrupoleScale * elecScale;
532 out[t020] = quadrupole[1][1] * quadrupoleScale * elecScale;
533 out[t002] = quadrupole[2][2] * quadrupoleScale * elecScale;
534 out[t110] = quadrupole[0][1] * quadrupoleScale * elecScale;
535 out[t101] = quadrupole[0][2] * quadrupoleScale * elecScale;
536 out[t011] = quadrupole[1][2] * quadrupoleScale * elecScale;
537 }
538 } else {
539
540
541 double[] out = globalMultipole[iSymm][ii];
542 out[t000] = in[t000] * chargeScale * elecScale;
543 out[t100] = in[t100] * dipoleScale * elecScale;
544 out[t010] = in[t010] * dipoleScale * elecScale;
545 out[t001] = in[t001] * dipoleScale * elecScale;
546 out[t200] = in[t200] * quadrupoleScale * elecScale;
547 out[t020] = in[t020] * quadrupoleScale * elecScale;
548 out[t002] = in[t002] * quadrupoleScale * elecScale;
549 out[t110] = in[t110] * quadrupoleScale * elecScale;
550 out[t101] = in[t101] * quadrupoleScale * elecScale;
551 out[t011] = in[t011] * quadrupoleScale * elecScale;
552
553 if (esvTerm && esvSystem.isTitrating(ii)) {
554 double[] esvMultipoleTitrDot = new double[10];
555 double titrationLambda = esvSystem.getTitrationLambda(ii);
556 double tautomerLambda = esvSystem.getTautomerLambda(ii);
557 esvMultipoleTitrDot = esvSystem.getTitrationUtils()
558 .getMultipoleTitrationDeriv(atom, titrationLambda,
559 tautomerLambda, esvMultipoleTitrDot);
560 in = esvMultipoleTitrDot;
561 out = titrationMultipole[iSymm][ii];
562 out[t000] = in[t000] * chargeScale * elecScale;
563 out[t100] = in[t100] * dipoleScale * elecScale;
564 out[t010] = in[t010] * dipoleScale * elecScale;
565 out[t001] = in[t001] * dipoleScale * elecScale;
566 out[t200] = in[t200] * quadrupoleScale * elecScale;
567 out[t020] = in[t020] * quadrupoleScale * elecScale;
568 out[t002] = in[t002] * quadrupoleScale * elecScale;
569 out[t110] = in[t110] * quadrupoleScale * elecScale;
570 out[t101] = in[t101] * quadrupoleScale * elecScale;
571 out[t011] = in[t011] * quadrupoleScale * elecScale;
572
573 double[] esvMultipoleTautDot = new double[10];
574 esvMultipoleTautDot = esvSystem.getTitrationUtils()
575 .getMultipoleTautomerDeriv(atom, titrationLambda,
576 tautomerLambda, esvMultipoleTautDot);
577 in = esvMultipoleTautDot;
578 out = tautomerMultipole[iSymm][ii];
579 out[t000] = in[t000] * chargeScale * elecScale;
580 out[t100] = in[t100] * dipoleScale * elecScale;
581 out[t010] = in[t010] * dipoleScale * elecScale;
582 out[t001] = in[t001] * dipoleScale * elecScale;
583 out[t200] = in[t200] * quadrupoleScale * elecScale;
584 out[t020] = in[t020] * quadrupoleScale * elecScale;
585 out[t002] = in[t002] * quadrupoleScale * elecScale;
586 out[t110] = in[t110] * quadrupoleScale * elecScale;
587 out[t101] = in[t101] * quadrupoleScale * elecScale;
588 out[t011] = in[t011] * quadrupoleScale * elecScale;
589 }
590 }
591
592
593 PolarizeType polarizeType = particleMeshEwald.getPolarizeType(ii);
594 if (polarizeType != null) {
595 polarizability[ii] = polarizeType.polarizability * polarizabilityScale * elecScale;
596 if (esvTerm && esvSystem.isTitrating(ii) && (esvSystem.isTitratingHydrogen(ii)
597 || esvSystem.isTitratingHeavy(ii))) {
598 titrationPolarizability[ii] = 0.0;
599 tautomerPolarizability[ii] = 0.0;
600 double titrationLambda = esvSystem.getTitrationLambda(ii);
601 double tautomerLambda = esvSystem.getTautomerLambda(ii);
602 titrationPolarizability[ii] = esvSystem.getTitrationUtils()
603 .getPolarizabilityTitrationDeriv(atom, titrationLambda, tautomerLambda);
604 tautomerPolarizability[ii] = esvSystem.getTitrationUtils()
605 .getPolarizabilityTautomerDeriv(atom, titrationLambda, tautomerLambda);
606 }
607 thole[ii] = polarizeType.thole;
608 ipdamp[ii] = polarizeType.pdamp;
609 if (!(ipdamp[ii] > 0.0)) {
610 ipdamp[ii] = Double.POSITIVE_INFINITY;
611 } else {
612 ipdamp[ii] = 1.0 / ipdamp[ii];
613 }
614 } else {
615 polarizability[ii] = 0.0;
616 thole[ii] = 0.0;
617 ipdamp[ii] = Double.POSITIVE_INFINITY;
618 }
619 }
620 }
621 }
622
623 @Override
624 public IntegerSchedule schedule() {
625 return IntegerSchedule.fixed();
626 }
627 }
628 }