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;
39
40 import edu.rit.pj.IntegerSchedule;
41 import edu.rit.pj.ParallelTeam;
42 import edu.rit.pj.reduction.SharedDouble;
43 import edu.rit.util.Range;
44 import ffx.crystal.Crystal;
45 import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
46 import ffx.numerics.atomic.AtomicDoubleArray3D;
47 import ffx.numerics.multipole.MultipoleTensor;
48 import ffx.potential.Platform;
49 import ffx.potential.bonded.Atom;
50 import ffx.potential.bonded.Bond;
51 import ffx.potential.bonded.LambdaInterface;
52 import ffx.potential.extended.ExtendedSystem;
53 import ffx.potential.nonbonded.pme.*;
54 import ffx.potential.parameters.AtomType;
55 import ffx.potential.parameters.ForceField;
56 import ffx.potential.parameters.ForceField.ELEC_FORM;
57 import ffx.potential.parameters.ForceField.ForceFieldType;
58 import ffx.potential.parameters.MultipoleType;
59 import ffx.potential.parameters.PolarizeType;
60 import ffx.potential.utils.EnergyException;
61 import ffx.utilities.Constants;
62 import ffx.utilities.FFXProperty;
63 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
64 import org.apache.commons.math3.linear.EigenDecomposition;
65
66 import java.util.ArrayList;
67 import java.util.Arrays;
68 import java.util.List;
69 import java.util.logging.Level;
70 import java.util.logging.Logger;
71
72 import static ffx.potential.nonbonded.pme.EwaldParameters.DEFAULT_EWALD_COEFFICIENT;
73 import static ffx.potential.parameters.ForceField.ELEC_FORM.PAM;
74 import static ffx.potential.parameters.ForceField.toEnumForm;
75 import static ffx.potential.parameters.MultipoleType.*;
76 import static ffx.utilities.Constants.ELEC_ANG_TO_DEBYE;
77 import static ffx.utilities.Constants.NS2SEC;
78 import static ffx.utilities.PropertyGroup.ElectrostaticsFunctionalForm;
79 import static ffx.utilities.PropertyGroup.LocalGeometryFunctionalForm;
80 import static java.lang.String.format;
81 import static java.lang.System.arraycopy;
82 import static java.util.Arrays.fill;
83 import static java.util.Collections.sort;
84 import static org.apache.commons.math3.util.FastMath.*;
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114 @SuppressWarnings("deprecation")
115 public class ParticleMeshEwald implements LambdaInterface {
116
117 private static final Logger logger = Logger.getLogger(ParticleMeshEwald.class.getName());
118
119
120
121 private static final int tensorCount = MultipoleTensor.tensorCount(3);
122
123
124
125
126 private final boolean lambdaTerm;
127
128
129
130 private final boolean nnTerm;
131 private final boolean reciprocalSpaceTerm;
132
133
134
135 private final ForceField forceField;
136 private static final double DEFAULT_POLAR_EPS = 1.0e-6;
137 @FFXProperty(name = "polar-eps", propertyGroup = ElectrostaticsFunctionalForm, defaultValue = "1.0e-6",
138 description = """
139 Sets the convergence criterion applied during computation of self-consistent induced dipoles.
140 The calculation is deemed to have converged when the rms change in Debyes in
141 the induced dipole at all polarizable sites is less than the value specified with this property.
142 The default value in the absence of the keyword is 1.0e-6 Debyes.
143 """)
144 private final double poleps;
145 private final boolean directFallback;
146
147
148
149 private final SCFPredictorParameters scfPredictorParameters;
150 private final EwaldParameters ewaldParameters;
151 private final ScaleParameters scaleParameters;
152 private final AlchemicalParameters alchemicalParameters;
153 private final RealSpaceNeighborParameters realSpaceNeighborParameters;
154 private final PMETimings pmeTimings;
155
156
157
158 private final int maxThreads;
159
160
161
162
163 private final ParallelTeam parallelTeam;
164
165
166
167
168
169
170
171
172 private final ParallelTeam sectionTeam;
173
174
175
176
177
178
179
180 private final ParallelTeam realSpaceTeam;
181
182
183
184
185
186
187
188 private final ParallelTeam fftTeam;
189 private final NeighborList neighborList;
190 private final InitializationRegion initializationRegion;
191 private final PermanentFieldRegion permanentFieldRegion;
192 private final InducedDipoleFieldRegion inducedDipoleFieldRegion;
193 private final InducedDipoleFieldReduceRegion inducedDipoleFieldReduceRegion;
194 private final ExpandInducedDipolesRegion expandInducedDipolesRegion;
195 private final DirectRegion directRegion;
196 private final SORRegion sorRegion;
197 private final OPTRegion optRegion;
198 private final PCGSolver pcgSolver;
199 private final ReciprocalSpace reciprocalSpace;
200 private final ReciprocalEnergyRegion reciprocalEnergyRegion;
201 private final RealSpaceEnergyRegion realSpaceEnergyRegion;
202 private final ReduceRegion reduceRegion;
203 private final GeneralizedKirkwood generalizedKirkwood;
204
205
206
207
208 public Polarization polarization;
209
210
211
212 public double[][][] coordinates;
213
214
215
216 public int[][][] neighborLists;
217
218
219
220 public double[][][] globalMultipole;
221
222
223
224 public double[][][] fractionalMultipole;
225
226
227
228 public double[][][] inducedDipole;
229 public double[][][] inducedDipoleCR;
230
231
232
233 public double[][] directDipole;
234 public double[][] directDipoleCR;
235 public double[][] directField;
236 public double[][] directFieldCR;
237
238
239
240 public double[][][] vacuumInducedDipole;
241 public double[][][] vacuumInducedDipoleCR;
242
243
244
245 public double[][] vacuumDirectDipole;
246 public double[][] vacuumDirectDipoleCR;
247
248
249
250 @FFXProperty(name = "electric", propertyGroup = LocalGeometryFunctionalForm, defaultValue = "332.063713",
251 description = """
252 Specifies a value for the so-called "electric constant" allowing conversion unit of electrostatic
253 potential energy values from electrons^2/Angstrom to kcal/mol. Internally, FFX stores a default value
254 for this constant as 332.063713 based on CODATA reference values. Since different force fields are
255 intended for use with slightly different values, this keyword allows overriding the default value.
256 """)
257 public double electric;
258
259
260
261 public double soluteDielectric;
262
263
264
265 protected Atom[] atoms;
266
267
268
269 protected int nAtoms;
270
271
272
273 protected int[][] ip11;
274 protected int[][] ip12;
275 protected int[][] ip13;
276
277
278
279
280 protected double totalMultipoleEnergy;
281
282
283
284
285 protected double permanentMultipoleEnergy;
286 protected double permanentRealSpaceEnergy;
287 protected double permanentSelfEnergy;
288 protected double permanentReciprocalEnergy;
289
290
291
292 protected double polarizationEnergy;
293 protected double inducedRealSpaceEnergy;
294 protected double inducedSelfEnergy;
295 protected double inducedReciprocalEnergy;
296 protected SCFAlgorithm scfAlgorithm;
297
298
299
300 private final SharedDouble shareddEdLambda;
301
302
303
304 private final SharedDouble sharedd2EdLambda2;
305
306
307
308 private final ELEC_FORM elecForm;
309
310
311
312 private Crystal crystal;
313
314
315
316 private int nSymm;
317
318
319
320 private boolean generalizedKirkwoodTerm;
321
322
323
324 private boolean gradient = false;
325
326
327
328 private int interactions;
329
330
331
332 private int gkInteractions;
333
334
335
336 private double solvationEnergy;
337
338
339
340 private LambdaMode lambdaMode = LambdaMode.OFF;
341
342
343
344 private double lambda = 1.0;
345
346
347
348 private double[][] localMultipole;
349 private MultipoleFrameDefinition[] frame;
350 private int[][] axisAtom;
351 private double[] ipdamp;
352 private double[] thole;
353 private double[] polarizability;
354
355
356
357 private int[][] mask12;
358 private int[][] mask13;
359 private int[][] mask14;
360 private int[][] mask15;
361
362
363
364 private boolean[] isSoft;
365
366
367
368 private int[] molecule;
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387 private boolean[] use;
388 private IntegerSchedule permanentSchedule;
389 private double[][] cartesianMultipolePhi;
390 private double[][] fracMultipolePhi;
391 private double[][] cartesianInducedDipolePhi;
392 private double[][] cartesianInducedDipolePhiCR;
393 private double[][] fractionalInducedDipolePhi;
394 private double[][] fractionalInducedDipolePhiCR;
395 private double[][] cartesianVacuumDipolePhi;
396 private double[][] cartesianVacuumDipolePhiCR;
397 private double[][] fractionalVacuumDipolePhi;
398 private double[][] fractionalVacuumDipolePhiCR;
399
400
401
402 private AtomicDoubleArrayImpl atomicDoubleArrayImpl;
403
404
405
406 private AtomicDoubleArray3D field;
407
408
409
410 private AtomicDoubleArray3D fieldCR;
411
412
413
414 private AtomicDoubleArray3D grad;
415
416
417
418 private AtomicDoubleArray3D torque;
419
420
421
422 private AtomicDoubleArray3D lambdaGrad;
423
424
425
426 private AtomicDoubleArray3D lambdaTorque;
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441 public ParticleMeshEwald(
442 Atom[] atoms,
443 int[] molecule,
444 ForceField forceField,
445 Crystal crystal,
446 NeighborList neighborList,
447 ELEC_FORM elecForm,
448 double ewaldCutoff,
449 double gkCutoff,
450 ParallelTeam parallelTeam) {
451 this.atoms = atoms;
452 this.molecule = molecule;
453 this.forceField = forceField;
454 this.crystal = crystal;
455 this.parallelTeam = parallelTeam;
456 this.neighborList = neighborList;
457 this.elecForm = elecForm;
458 neighborLists = neighborList.getNeighborList();
459 permanentSchedule = neighborList.getPairwiseSchedule();
460 nAtoms = atoms.length;
461 nSymm = crystal.spaceGroup.getNumberOfSymOps();
462 maxThreads = parallelTeam.getThreadCount();
463
464 electric = forceField.getDouble("ELECTRIC", Constants.DEFAULT_ELECTRIC);
465
466
467 soluteDielectric = forceField.getDouble("SOLUTE_DIELECTRIC", 1.0);
468 if (soluteDielectric > 1.0) {
469 electric = electric / soluteDielectric;
470 } else {
471 soluteDielectric = 1.0;
472 }
473
474 poleps = forceField.getDouble("POLAR_EPS", DEFAULT_POLAR_EPS);
475
476
477 lambdaTerm = forceField.getBoolean("ELEC_LAMBDATERM",
478 forceField.getBoolean("LAMBDATERM", false));
479
480 nnTerm = forceField.getBoolean("NNTERM", false);
481
482 if (nnTerm && lambdaTerm) {
483 logger.severe(" Use a neural network potential with alchemical simulations is not yet supported.");
484 }
485
486 double aewald = forceField.getDouble("EWALD_ALPHA", DEFAULT_EWALD_COEFFICIENT);
487 ewaldParameters = new EwaldParameters(ewaldCutoff, aewald);
488 scaleParameters = new ScaleParameters(elecForm, forceField);
489 reciprocalSpaceTerm = forceField.getBoolean("RECIPTERM", true);
490
491 SCFPredictor scfPredictor;
492 try {
493 String predictor = forceField.getString("SCF_PREDICTOR", "NONE");
494 predictor = predictor.replaceAll("-", "_").toUpperCase();
495 scfPredictor = SCFPredictor.valueOf(predictor);
496 } catch (Exception e) {
497 scfPredictor = SCFPredictor.NONE;
498 }
499
500 scfPredictorParameters = new SCFPredictorParameters(scfPredictor, nAtoms);
501 if (scfPredictor != SCFPredictor.NONE) {
502 scfPredictorParameters.init(forceField);
503 }
504
505 String algorithm = forceField.getString("SCF_ALGORITHM", "CG");
506 try {
507 algorithm = algorithm.replaceAll("-", "_").toUpperCase();
508 scfAlgorithm = SCFAlgorithm.valueOf(algorithm);
509 } catch (Exception e) {
510 scfAlgorithm = SCFAlgorithm.CG;
511 }
512
513
514 directFallback = forceField.getBoolean("DIRECT_SCF_FALLBACK", true);
515
516
517 String value = forceField.getString("ARRAY_REDUCTION", "MULTI");
518 try {
519 atomicDoubleArrayImpl = AtomicDoubleArrayImpl.valueOf(toEnumForm(value));
520 } catch (Exception e) {
521 atomicDoubleArrayImpl = AtomicDoubleArrayImpl.MULTI;
522 logger.info(
523 format(
524 " Unrecognized ARRAY-REDUCTION %s; defaulting to %s", value, atomicDoubleArrayImpl));
525 }
526 logger.fine(format(" PME using %s arrays.", atomicDoubleArrayImpl.toString()));
527
528 if (!scfAlgorithm.isSupported(Platform.FFX)) {
529
530
531 logger.fine(format(
532 " SCF algorithm %s is not supported by FFX reference implementation; falling back to CG!",
533 scfAlgorithm));
534 scfAlgorithm = SCFAlgorithm.CG;
535 }
536
537 pcgSolver = new PCGSolver(maxThreads, poleps, forceField, nAtoms);
538
539 alchemicalParameters = new AlchemicalParameters(forceField, lambdaTerm, nnTerm, polarization);
540 if (nnTerm) {
541 initNN();
542 }
543
544 String polar = forceField.getString("POLARIZATION", "MUTUAL");
545 if (elecForm == ELEC_FORM.FIXED_CHARGE) {
546 polar = "NONE";
547 }
548 boolean polarizationTerm = forceField.getBoolean("POLARIZETERM", true);
549 if (!polarizationTerm || polar.equalsIgnoreCase("NONE")) {
550 polarization = Polarization.NONE;
551 } else if (polar.equalsIgnoreCase("DIRECT")) {
552 polarization = Polarization.DIRECT;
553 } else {
554 polarization = Polarization.MUTUAL;
555 }
556
557 if (lambdaTerm) {
558 shareddEdLambda = new SharedDouble();
559 sharedd2EdLambda2 = new SharedDouble();
560 } else {
561 shareddEdLambda = null;
562 sharedd2EdLambda2 = null;
563 lambdaGrad = null;
564 lambdaTorque = null;
565 }
566
567 directRegion = new DirectRegion(maxThreads);
568 sorRegion = new SORRegion(maxThreads, forceField);
569 optRegion = new OPTRegion(maxThreads);
570
571 if (logger.isLoggable(Level.INFO)) {
572 StringBuilder sb = new StringBuilder();
573 sb.append("\n Electrostatics\n");
574 sb.append(format(" Polarization: %8s\n", polarization.toString()));
575 if (polarization == Polarization.MUTUAL) {
576 sb.append(format(" SCF Convergence Criteria: %8.3e\n", poleps));
577 sb.append(format(" SCF Predictor: %8s\n",
578 scfPredictorParameters.scfPredictor));
579 sb.append(format(" SCF Algorithm: %8s\n", scfAlgorithm));
580 if (scfAlgorithm == SCFAlgorithm.SOR) {
581 sb.append(format(" SOR Parameter: %8.3f\n", sorRegion.getSOR()));
582 } else {
583 sb.append(format(" CG Preconditioner Cut-Off: %8.3f\n",
584 pcgSolver.getPreconditionerCutoff()));
585 sb.append(format(" CG Preconditioner Ewald Coeff.: %8.3f\n",
586 pcgSolver.getPreconditionerEwald()));
587 sb.append(format(" CG Preconditioner Scale: %8.3f\n",
588 pcgSolver.getPreconditionerScale()));
589 sb.append(
590 format(" CG Preconditioner Mode: %15s\n", pcgSolver.getPreconditionerMode()));
591 }
592 }
593 if (ewaldParameters.aewald > 0.0) {
594 sb.append(" Particle-mesh Ewald\n");
595 sb.append(format(" Ewald Coefficient: %8.3f\n", ewaldParameters.aewald));
596 sb.append(format(" Particle Cut-Off: %8.3f (A)", ewaldParameters.off));
597 } else if (ewaldParameters.off != Double.POSITIVE_INFINITY) {
598 sb.append(format(" Electrostatics Cut-Off: %8.3f (A)\n",
599 ewaldParameters.off));
600 } else {
601 sb.append(" Electrostatics Cut-Off: NONE\n");
602 }
603
604 logger.info(sb.toString());
605 }
606
607
608 int sectionThreads;
609
610
611
612
613
614 int realSpaceThreads;
615
616
617
618
619
620 int reciprocalThreads;
621
622 boolean concurrent;
623 int realThreads = 1;
624 try {
625 realThreads = forceField.getInteger("PME_REAL_THREADS");
626 if (realThreads >= maxThreads || realThreads < 1) {
627 throw new Exception("pme-real-threads must be < ffx.nt and greater than 0");
628 }
629 concurrent = true;
630 } catch (Exception e) {
631 concurrent = false;
632 }
633 if (concurrent) {
634 sectionThreads = 2;
635 realSpaceThreads = realThreads;
636 reciprocalThreads = maxThreads - realThreads;
637 sectionTeam = new ParallelTeam(sectionThreads);
638 realSpaceTeam = new ParallelTeam(realSpaceThreads);
639 fftTeam = new ParallelTeam(reciprocalThreads);
640 } else {
641
642 sectionThreads = 1;
643 sectionTeam = new ParallelTeam(sectionThreads);
644 realSpaceTeam = parallelTeam;
645 fftTeam = parallelTeam;
646 }
647
648 realSpaceNeighborParameters = new RealSpaceNeighborParameters(maxThreads);
649 initializationRegion = new InitializationRegion(this, maxThreads, forceField);
650 expandInducedDipolesRegion = new ExpandInducedDipolesRegion(maxThreads);
651 initAtomArrays();
652
653
654
655
656
657
658 if (ewaldParameters.aewald > 0.0 && reciprocalSpaceTerm) {
659 reciprocalSpace = new ReciprocalSpace(this,
660 crystal.getUnitCell(), forceField, atoms,
661 ewaldParameters.aewald, fftTeam, parallelTeam);
662 reciprocalEnergyRegion = new ReciprocalEnergyRegion(maxThreads, ewaldParameters.aewald, electric);
663 } else {
664 reciprocalSpace = null;
665 reciprocalEnergyRegion = null;
666 }
667 permanentFieldRegion = new PermanentFieldRegion(realSpaceTeam, forceField, lambdaTerm);
668 inducedDipoleFieldRegion = new InducedDipoleFieldRegion(realSpaceTeam, forceField, lambdaTerm);
669 inducedDipoleFieldReduceRegion = new InducedDipoleFieldReduceRegion(maxThreads);
670 realSpaceEnergyRegion = new RealSpaceEnergyRegion(elecForm, maxThreads, forceField, lambdaTerm, electric);
671 reduceRegion = new ReduceRegion(maxThreads, forceField);
672
673 pmeTimings = new PMETimings(maxThreads);
674
675 if (lambdaTerm) {
676 logger.info(alchemicalParameters.toString());
677 }
678
679
680
681 generalizedKirkwoodTerm = forceField.getBoolean("GKTERM", false);
682 if (generalizedKirkwoodTerm || alchemicalParameters.doLigandGKElec) {
683 generalizedKirkwood = new GeneralizedKirkwood(forceField, atoms, this, crystal, parallelTeam, gkCutoff);
684 } else {
685 generalizedKirkwood = null;
686 }
687 }
688
689
690
691
692 public ELEC_FORM getElecForm() {
693 return elecForm;
694 }
695
696 public void computeInduceDipoleField() {
697 expandInducedDipoles();
698
699 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
700 reciprocalSpace.splineInducedDipoles(inducedDipole, inducedDipoleCR, use);
701 }
702 field.reset(parallelTeam);
703 fieldCR.reset(parallelTeam);
704 inducedDipoleFieldRegion.init(
705 atoms, crystal, use, molecule, ipdamp, thole, coordinates, realSpaceNeighborParameters,
706 inducedDipole, inducedDipoleCR, reciprocalSpaceTerm, reciprocalSpace, lambdaMode,
707 ewaldParameters, field, fieldCR, pmeTimings);
708 inducedDipoleFieldRegion.executeWith(sectionTeam);
709 pmeTimings.realSpaceSCFTotal = inducedDipoleFieldRegion.getRealSpaceSCFTotal();
710
711 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
712 reciprocalSpace.computeInducedPhi(
713 cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
714 fractionalInducedDipolePhi, fractionalInducedDipolePhiCR);
715 }
716
717 if (generalizedKirkwoodTerm) {
718 pmeTimings.gkEnergyTotal = -System.nanoTime();
719 generalizedKirkwood.computeInducedGKField();
720 pmeTimings.gkEnergyTotal += System.nanoTime();
721 logger.fine(
722 format(" Computed GK induced field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
723 }
724
725 inducedDipoleFieldReduceRegion.init(
726 atoms,
727 inducedDipole,
728 inducedDipoleCR,
729 generalizedKirkwoodTerm,
730 generalizedKirkwood,
731 ewaldParameters,
732 soluteDielectric,
733 cartesianInducedDipolePhi,
734 cartesianInducedDipolePhiCR,
735 field,
736 fieldCR);
737 inducedDipoleFieldReduceRegion.executeWith(parallelTeam);
738 }
739
740 public void destroy() {
741 if (fftTeam != null) {
742 try {
743 fftTeam.shutdown();
744 } catch (Exception ex) {
745 logger.warning(" Exception in shutting down fftTeam");
746 }
747 }
748 if (sectionTeam != null) {
749 try {
750 sectionTeam.shutdown();
751 } catch (Exception ex) {
752 logger.warning(" Exception in shutting down sectionTeam");
753 }
754 }
755 if (realSpaceTeam != null) {
756 try {
757 realSpaceTeam.shutdown();
758 } catch (Exception ex) {
759 logger.warning(" Exception in shutting down realSpaceTeam");
760 }
761 }
762 }
763
764
765
766
767
768
769
770
771 public double energy(boolean gradient, boolean print) {
772
773 this.gradient = gradient;
774
775
776 totalMultipoleEnergy = 0.0;
777 permanentMultipoleEnergy = 0.0;
778 permanentRealSpaceEnergy = 0.0;
779 permanentSelfEnergy = 0.0;
780 permanentReciprocalEnergy = 0.0;
781 polarizationEnergy = 0.0;
782 inducedRealSpaceEnergy = 0.0;
783 inducedSelfEnergy = 0.0;
784 inducedReciprocalEnergy = 0.0;
785 solvationEnergy = 0.0;
786
787
788 interactions = 0;
789 gkInteractions = 0;
790
791
792 pmeTimings.init();
793 if (reciprocalSpace != null) {
794 reciprocalSpace.initTimings();
795 }
796
797
798 if (lambdaTerm) {
799 shareddEdLambda.set(0.0);
800 sharedd2EdLambda2.set(0.0);
801 }
802
803 if (esvTerm) {
804 extendedSystem.initEsvPermElec();
805 extendedSystem.initEsvIndElec();
806 }
807
808 alchemicalParameters.doPermanentRealSpace = true;
809 alchemicalParameters.permanentScale = 1.0;
810 alchemicalParameters.doPolarization = true;
811 alchemicalParameters.polarizationScale = 1.0;
812
813
814 initializationRegion.init(lambdaTerm, extendedSystem,
815 atoms, coordinates, crystal, frame, axisAtom, globalMultipole,
816 dMultipoledTirationESV, dMultipoledTautomerESV,
817 polarizability, dPolardTitrationESV, dPolardTautomerESV,
818 thole, ipdamp, use, neighborLists,
819 realSpaceNeighborParameters.realSpaceLists,
820 alchemicalParameters.vaporLists,
821 grad, torque, lambdaGrad, lambdaTorque);
822 initializationRegion.executeWith(parallelTeam);
823
824
825 if (generalizedKirkwoodTerm || alchemicalParameters.doLigandGKElec) {
826 generalizedKirkwood.init();
827 }
828
829
830 double energy;
831 if (!lambdaTerm) {
832 lambdaMode = LambdaMode.OFF;
833 energy = computeEnergy(print);
834 if (nnTerm) {
835 lambdaMode = LambdaMode.VAPOR;
836 double temp = energy;
837 energy = nnElec();
838 logger.fine(format(" Vacuum energy: %20.8f", energy - temp));
839 }
840 } else {
841
842 lambdaMode = LambdaMode.CONDENSED;
843 energy = condensedEnergy();
844 if (logger.isLoggable(Level.FINE)) {
845 logger.fine(format(" Condensed energy: %20.8f", energy));
846 }
847
848
849 lambdaMode = LambdaMode.CONDENSED_NO_LIGAND;
850 double temp = energy;
851 energy = condensedNoLigandSCF();
852 if (logger.isLoggable(Level.FINE)) {
853 logger.fine(format(" Condensed no ligand energy: %20.8f", energy - temp));
854 }
855
856
857 if (alchemicalParameters.doLigandVaporElec) {
858 lambdaMode = LambdaMode.VAPOR;
859 temp = energy;
860 energy = ligandElec();
861 if (logger.isLoggable(Level.FINE)) {
862 logger.fine(format(" Vacuum energy: %20.8f", energy - temp));
863 }
864 }
865 }
866
867
868
869
870
871 if (gradient || lambdaTerm) {
872 reduceRegion.init(lambdaTerm, gradient,
873 atoms, coordinates, frame, axisAtom,
874 grad, torque, lambdaGrad, lambdaTorque);
875 reduceRegion.excuteWith(parallelTeam);
876 }
877
878
879 if (logger.isLoggable(Level.FINE)) {
880 pmeTimings.printRealSpaceTimings(maxThreads, realSpaceEnergyRegion);
881 if (ewaldParameters.aewald > 0.0 && reciprocalSpaceTerm) {
882 reciprocalSpace.printTimings();
883 }
884 }
885
886 return permanentMultipoleEnergy + polarizationEnergy;
887 }
888
889 public void expandInducedDipoles() {
890 if (nSymm > 1) {
891 expandInducedDipolesRegion.init(atoms, crystal, inducedDipole, inducedDipoleCR);
892 expandInducedDipolesRegion.executeWith(parallelTeam);
893 }
894 }
895
896 public int[][] getAxisAtoms() {
897 return axisAtom;
898 }
899
900 public double getCavitationEnergy() {
901 return generalizedKirkwood.getCavitationEnergy();
902 }
903
904 public double[][][] getCoordinates() {
905 return coordinates;
906 }
907
908 public double getDispersionEnergy() {
909 return generalizedKirkwood.getDispersionEnergy();
910 }
911
912 public double getEwaldCoefficient() {
913 return ewaldParameters.aewald;
914 }
915
916 public double getEwaldCutoff() {
917 return ewaldParameters.off;
918 }
919
920 public GeneralizedKirkwood getGK() {
921 return generalizedKirkwood;
922 }
923
924
925
926
927
928
929 public double getGKEnergy() {
930 return generalizedKirkwood.getGeneralizedKirkwoordEnergy();
931 }
932
933
934
935
936
937
938 public int getGKInteractions() {
939 return gkInteractions;
940 }
941
942
943
944
945
946
947 public int getInteractions() {
948 return interactions;
949 }
950
951
952
953
954
955
956 @Override
957 public double getLambda() {
958 return lambda;
959 }
960
961
962
963
964
965
966 @Override
967 public void setLambda(double lambda) {
968 assert (lambda >= 0.0 && lambda <= 1.0);
969 if (!lambdaTerm) {
970 return;
971 }
972 this.lambda = lambda;
973
974 initSoftCore();
975 alchemicalParameters.update(lambda);
976 if (generalizedKirkwoodTerm) {
977 generalizedKirkwood.setLambda(alchemicalParameters.polLambda);
978 generalizedKirkwood.setLambdaFunction(
979 alchemicalParameters.lPowPol,
980 alchemicalParameters.dlPowPol,
981 alchemicalParameters.d2lPowPol);
982 }
983 }
984
985 public double getPermanentEnergy() {
986 return permanentMultipoleEnergy;
987 }
988
989
990
991
992
993
994 public double getIndRealEnergy() {
995 return inducedRealSpaceEnergy;
996 }
997
998
999
1000
1001
1002
1003 public double getIndRecipEnergy() {
1004 return inducedReciprocalEnergy;
1005 }
1006
1007
1008
1009
1010
1011
1012 public double getIndSelfEnergy() {
1013 return inducedSelfEnergy;
1014 }
1015
1016
1017
1018
1019
1020
1021 public double getPermSelfEnergy() {
1022 return permanentSelfEnergy;
1023 }
1024
1025 public double getPermRealEnergy() {
1026 return permanentRealSpaceEnergy;
1027 }
1028
1029 public double getPermRecipEnergy() {
1030 return permanentReciprocalEnergy;
1031 }
1032
1033 public double getPolarEps() {
1034 return poleps;
1035 }
1036
1037 public int[][] getPolarization11() {
1038 return ip11;
1039 }
1040
1041 public int[][] getPolarization12() {
1042 return ip12;
1043 }
1044
1045 public int[][] getPolarization13() {
1046 return ip13;
1047 }
1048
1049
1050
1051
1052
1053
1054 public double getPolarizationEnergy() {
1055 return polarizationEnergy;
1056 }
1057
1058 public Polarization getPolarizationType() {
1059 return polarization;
1060 }
1061
1062 public ReciprocalSpace getReciprocalSpace() {
1063 return reciprocalSpace;
1064 }
1065
1066 public double getScale14() {
1067 return scaleParameters.m14scale;
1068 }
1069
1070
1071
1072
1073
1074
1075 public double getSolvationEnergy() {
1076 return solvationEnergy;
1077 }
1078
1079
1080
1081
1082
1083
1084
1085 public MultipoleType getMultipoleType(int i) {
1086 Atom atom = atoms[i];
1087 MultipoleType multipoleType = atom.getMultipoleType();
1088 if (esvTerm && extendedSystem.isTitrating(i)) {
1089 double[] multipole = multipoleType.getMultipole();
1090 double[] esvMultipole = new double[10];
1091 arraycopy(multipole, 0, esvMultipole, 0, multipole.length);
1092 double titrationLambda = extendedSystem.getTitrationLambda(i);
1093 double tautomerLambda = extendedSystem.getTautomerLambda(i);
1094 esvMultipole = extendedSystem.getTitrationUtils()
1095 .getMultipole(atom, titrationLambda, tautomerLambda, esvMultipole);
1096
1097 multipoleType = new MultipoleType(multipoleType, esvMultipole);
1098 }
1099 return multipoleType;
1100 }
1101
1102
1103
1104
1105
1106
1107
1108 public PolarizeType getPolarizeType(int i) {
1109 Atom atom = atoms[i];
1110 PolarizeType polarizeType = atom.getPolarizeType();
1111 if (polarizeType != null) {
1112 if (esvTerm && extendedSystem.isTitrating(i) && (extendedSystem.isTitratingHydrogen(i)
1113 || extendedSystem.isTitratingHeavy(i))) {
1114 double titrationLambda = extendedSystem.getTitrationLambda(i);
1115 double tautomerLambda = extendedSystem.getTautomerLambda(i);
1116 double esvPolarizability = extendedSystem.getTitrationUtils()
1117 .getPolarizability(atom, titrationLambda, tautomerLambda, polarizeType.polarizability);
1118 polarizeType = new PolarizeType(polarizeType, esvPolarizability);
1119 }
1120 }
1121 return polarizeType;
1122 }
1123
1124
1125
1126
1127 @Override
1128 public double getd2EdL2() {
1129 if (sharedd2EdLambda2 == null || !lambdaTerm) {
1130 return 0.0;
1131 }
1132 double d2EdL2 = sharedd2EdLambda2.get();
1133 if (generalizedKirkwoodTerm || alchemicalParameters.doLigandGKElec) {
1134 d2EdL2 += generalizedKirkwood.getd2EdL2();
1135 }
1136 return d2EdL2;
1137 }
1138
1139
1140
1141
1142 @Override
1143 public double getdEdL() {
1144 if (shareddEdLambda == null || !lambdaTerm) {
1145 return 0.0;
1146 }
1147 double dEdL = shareddEdLambda.get();
1148 if (generalizedKirkwoodTerm || alchemicalParameters.doLigandGKElec) {
1149 dEdL += generalizedKirkwood.getdEdL();
1150 }
1151 return dEdL;
1152 }
1153
1154
1155
1156
1157 @Override
1158 public void getdEdXdL(double[] gradient) {
1159 if (lambdaGrad == null || !lambdaTerm) {
1160 return;
1161 }
1162
1163 int index = 0;
1164 for (int i = 0; i < nAtoms; i++) {
1165 if (atoms[i].isActive()) {
1166 gradient[index++] += lambdaGrad.getX(i);
1167 gradient[index++] += lambdaGrad.getY(i);
1168 gradient[index++] += lambdaGrad.getZ(i);
1169 }
1170 }
1171 }
1172
1173 public void setAtoms(Atom[] atoms, int[] molecule) {
1174 if (lambdaTerm && atoms.length != nAtoms) {
1175 logger.severe(" Changing the number of atoms is not compatible with use of Lambda.");
1176 }
1177 this.atoms = atoms;
1178 this.molecule = molecule;
1179 nAtoms = atoms.length;
1180 initAtomArrays();
1181
1182 if (reciprocalSpace != null) {
1183 reciprocalSpace.setAtoms(atoms);
1184 }
1185
1186 if (generalizedKirkwood != null) {
1187 generalizedKirkwood.setAtoms(atoms);
1188 }
1189 }
1190
1191 public void setCrystal(Crystal crystal) {
1192
1193 int nSymmNew = crystal.spaceGroup.getNumberOfSymOps();
1194 if (nSymm < nSymmNew) {
1195 coordinates = new double[nSymmNew][3][nAtoms];
1196 globalMultipole = new double[nSymmNew][nAtoms][10];
1197 fractionalMultipole = new double[nSymmNew][nAtoms][10];
1198 inducedDipole = new double[nSymmNew][nAtoms][3];
1199 inducedDipoleCR = new double[nSymmNew][nAtoms][3];
1200 if (generalizedKirkwood != null) {
1201 generalizedKirkwood.setAtoms(atoms);
1202 }
1203 realSpaceNeighborParameters.allocate(nAtoms, nSymmNew);
1204 pcgSolver.allocateLists(nSymmNew, nAtoms);
1205 }
1206 nSymm = nSymmNew;
1207 neighborLists = neighborList.getNeighborList();
1208 this.crystal = crystal;
1209 if (reciprocalSpace != null) {
1210 reciprocalSpace.setCrystal(crystal.getUnitCell());
1211 }
1212 }
1213
1214 private void initAtomArrays() {
1215 if (localMultipole == null || localMultipole.length < nAtoms) {
1216 localMultipole = new double[nAtoms][10];
1217 frame = new MultipoleFrameDefinition[nAtoms];
1218 axisAtom = new int[nAtoms][];
1219 cartesianMultipolePhi = new double[nAtoms][tensorCount];
1220 fracMultipolePhi = new double[nAtoms][tensorCount];
1221 directDipole = new double[nAtoms][3];
1222 directDipoleCR = new double[nAtoms][3];
1223 directField = new double[nAtoms][3];
1224 directFieldCR = new double[nAtoms][3];
1225 vacuumDirectDipole = new double[nAtoms][3];
1226 vacuumDirectDipoleCR = new double[nAtoms][3];
1227 if (optRegion != null) {
1228 int optOrder = optRegion.optOrder;
1229 optRegion.optDipole = new double[optOrder + 1][nAtoms][3];
1230 optRegion.optDipoleCR = new double[optOrder + 1][nAtoms][3];
1231 }
1232
1233 cartesianInducedDipolePhi = new double[nAtoms][tensorCount];
1234 cartesianInducedDipolePhiCR = new double[nAtoms][tensorCount];
1235 fractionalInducedDipolePhi = new double[nAtoms][tensorCount];
1236 fractionalInducedDipolePhiCR = new double[nAtoms][tensorCount];
1237
1238 cartesianVacuumDipolePhi = new double[nAtoms][tensorCount];
1239 cartesianVacuumDipolePhiCR = new double[nAtoms][tensorCount];
1240 fractionalVacuumDipolePhi = new double[nAtoms][tensorCount];
1241 fractionalVacuumDipolePhiCR = new double[nAtoms][tensorCount];
1242
1243 mask12 = new int[nAtoms][];
1244 mask13 = new int[nAtoms][];
1245 mask14 = new int[nAtoms][];
1246 mask15 = new int[nAtoms][];
1247 ip11 = new int[nAtoms][];
1248 ip12 = new int[nAtoms][];
1249 ip13 = new int[nAtoms][];
1250 thole = new double[nAtoms];
1251 ipdamp = new double[nAtoms];
1252 polarizability = new double[nAtoms];
1253
1254 if (scfAlgorithm == SCFAlgorithm.CG) {
1255 pcgSolver.allocateVectors(nAtoms);
1256 }
1257 pcgSolver.allocateLists(nSymm, nAtoms);
1258
1259 if (scfPredictorParameters.scfPredictor != SCFPredictor.NONE) {
1260 int predictorOrder = scfPredictorParameters.predictorOrder;
1261 if (lambdaTerm) {
1262 scfPredictorParameters.predictorInducedDipole = new double[3][predictorOrder][nAtoms][3];
1263 scfPredictorParameters.predictorInducedDipoleCR =
1264 new double[3][predictorOrder][nAtoms][3];
1265 } else {
1266 scfPredictorParameters.predictorInducedDipole = new double[1][predictorOrder][nAtoms][3];
1267 scfPredictorParameters.predictorInducedDipoleCR =
1268 new double[1][predictorOrder][nAtoms][3];
1269 }
1270 }
1271
1272
1273
1274 grad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1275 torque = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1276 field = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1277 fieldCR = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1278 lambdaGrad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1279 lambdaTorque = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, maxThreads);
1280 isSoft = new boolean[nAtoms];
1281 use = new boolean[nAtoms];
1282
1283 coordinates = new double[nSymm][3][nAtoms];
1284 globalMultipole = new double[nSymm][nAtoms][10];
1285 fractionalMultipole = new double[nSymm][nAtoms][10];
1286 inducedDipole = new double[nSymm][nAtoms][3];
1287 inducedDipoleCR = new double[nSymm][nAtoms][3];
1288 vacuumInducedDipole = new double[nSymm][nAtoms][3];
1289 vacuumInducedDipoleCR = new double[nSymm][nAtoms][3];
1290
1291
1292 realSpaceNeighborParameters.allocate(nAtoms, nSymm);
1293
1294
1295 lambdaFactors = new LambdaFactors[maxThreads];
1296 for (int i = 0; i < maxThreads; i++) {
1297
1298
1299
1300
1301 if (lambdaTerm) {
1302
1303 lambdaFactors[i] = new LambdaFactorsOST();
1304 } else {
1305
1306 lambdaFactors[i] = LambdaDefaults;
1307 }
1308 }
1309 }
1310
1311
1312 fill(isSoft, false);
1313
1314
1315 fill(use, true);
1316
1317
1318 assignMultipoles();
1319
1320
1321 if (elecForm == PAM) {
1322 assignPolarizationGroups();
1323 }
1324
1325
1326 for (int i = 0; i < nAtoms; i++) {
1327 Atom ai = atoms[i];
1328
1329 if (elecForm == PAM) {
1330 PolarizeType polarizeType = ai.getPolarizeType();
1331 int index = ai.getIndex() - 1;
1332 thole[index] = polarizeType.thole;
1333 ipdamp[index] = polarizeType.pdamp;
1334 if (!(ipdamp[index] > 0.0)) {
1335 ipdamp[index] = Double.POSITIVE_INFINITY;
1336 } else {
1337 ipdamp[index] = 1.0 / ipdamp[index];
1338 }
1339 polarizability[index] = polarizeType.polarizability;
1340 }
1341
1342
1343 List<Atom> n12 = ai.get12List();
1344 mask12[i] = new int[n12.size()];
1345 int j = 0;
1346 for (Atom a12 : n12) {
1347 mask12[i][j++] = a12.getIndex() - 1;
1348 }
1349
1350
1351 List<Atom> n13 = ai.get13List();
1352 mask13[i] = new int[n13.size()];
1353 j = 0;
1354 for (Atom a13 : n13) {
1355 mask13[i][j++] = a13.getIndex() - 1;
1356 }
1357
1358
1359 List<Atom> n14 = ai.get14List();
1360 mask14[i] = new int[n14.size()];
1361 j = 0;
1362 for (Atom a14 : n14) {
1363 mask14[i][j++] = a14.getIndex() - 1;
1364 }
1365
1366
1367 List<Atom> n15 = ai.get15List();
1368 mask15[i] = new int[n15.size()];
1369 j = 0;
1370 for (Atom a15 : n15) {
1371 mask15[i][j++] = a15.getIndex() - 1;
1372 }
1373 }
1374 }
1375
1376
1377
1378
1379 private void initSoftCore() {
1380
1381 boolean rebuild = false;
1382
1383
1384 StringBuilder sb = new StringBuilder("\n Softcore Atoms:\n");
1385 int count = 0;
1386 for (int i = 0; i < nAtoms; i++) {
1387 Atom ai = atoms[i];
1388 boolean soft = ai.applyLambda();
1389 if (soft != isSoft[i]) {
1390 rebuild = true;
1391 }
1392 isSoft[i] = soft;
1393 if (soft) {
1394 count++;
1395 sb.append(ai).append("\n");
1396 }
1397 }
1398
1399
1400 if (alchemicalParameters.doLigandVaporElec
1401 && alchemicalParameters.vaporCrystal == null) {
1402 rebuild = true;
1403 }
1404
1405 if (!rebuild) {
1406 return;
1407 }
1408
1409 if (count > 0 && logger.isLoggable(Level.FINE)) {
1410 logger.fine(format(" Softcore atom count: %d", count));
1411 logger.fine(sb.toString());
1412 }
1413
1414
1415
1416 if (alchemicalParameters.doLigandVaporElec) {
1417 double vacuumOff = getVacuumOff();
1418 alchemicalParameters.vaporCrystal =
1419 new Crystal(3 * vacuumOff, 3 * vacuumOff, 3 * vacuumOff, 90.0, 90.0, 90.0, "P1");
1420 alchemicalParameters.vaporCrystal.setAperiodic(true);
1421 NeighborList vacuumNeighborList = new NeighborList(null,
1422 alchemicalParameters.vaporCrystal, atoms, vacuumOff, 2.0, parallelTeam);
1423 vacuumNeighborList.setIntermolecular(false);
1424
1425 alchemicalParameters.vaporLists = new int[1][nAtoms][];
1426 double[][] coords = new double[1][nAtoms * 3];
1427 for (int i = 0; i < nAtoms; i++) {
1428 coords[0][i * 3] = atoms[i].getX();
1429 coords[0][i * 3 + 1] = atoms[i].getY();
1430 coords[0][i * 3 + 2] = atoms[i].getZ();
1431 }
1432 boolean print = logger.isLoggable(Level.FINE);
1433 vacuumNeighborList.buildList(coords, alchemicalParameters.vaporLists, isSoft,
1434 true, print);
1435 alchemicalParameters.vaporPermanentSchedule = vacuumNeighborList.getPairwiseSchedule();
1436 alchemicalParameters.vaporEwaldSchedule = alchemicalParameters.vaporPermanentSchedule;
1437 alchemicalParameters.vacuumRanges = new Range[maxThreads];
1438 vacuumNeighborList.setDisableUpdates(
1439 forceField.getBoolean("DISABLE_NEIGHBOR_UPDATES", false));
1440 } else {
1441 alchemicalParameters.vaporCrystal = null;
1442 alchemicalParameters.vaporLists = null;
1443 alchemicalParameters.vaporPermanentSchedule = null;
1444 alchemicalParameters.vaporEwaldSchedule = null;
1445 alchemicalParameters.vacuumRanges = null;
1446 }
1447 }
1448
1449 private double getVacuumOff() {
1450 double maxr = 10.0;
1451 for (int i = 0; i < nAtoms; i++) {
1452 Atom ai = atoms[i];
1453 if (ai.applyLambda()) {
1454
1455
1456 for (int j = i + 1; j < nAtoms; j++) {
1457 Atom aj = atoms[j];
1458 if (aj.applyLambda()) {
1459 double dx = ai.getX() - aj.getX();
1460 double dy = ai.getY() - aj.getY();
1461 double dz = ai.getZ() - aj.getZ();
1462 double r = sqrt(dx * dx + dy * dy + dz * dz);
1463 maxr = max(r, maxr);
1464 }
1465 }
1466 }
1467 }
1468
1469 double vacuumOff = 2 * maxr;
1470 return vacuumOff;
1471 }
1472
1473
1474
1475
1476 private void initNN() {
1477
1478
1479 StringBuilder sb = new StringBuilder("\n NN Atoms:\n");
1480 boolean[] nnAtoms = new boolean[nAtoms];
1481
1482 int count = 0;
1483 for (int i = 0; i < nAtoms; i++) {
1484 Atom ai = atoms[i];
1485 nnAtoms[i] = ai.isNeuralNetwork();
1486 if (nnAtoms[i]) {
1487 count++;
1488 sb.append(ai).append("\n");
1489 }
1490 }
1491
1492 if (count > 0 && logger.isLoggable(Level.FINE)) {
1493 logger.fine(format(" Neural network atom count: %d", count));
1494 logger.fine(sb.toString());
1495 }
1496
1497
1498
1499 if (alchemicalParameters.doLigandVaporElec) {
1500 alchemicalParameters.vaporCrystal = new Crystal(10.0, 10.0, 10.0,
1501 90.0, 90.0, 90.0, "P1");
1502 alchemicalParameters.vaporCrystal.setAperiodic(true);
1503 NeighborList vacuumNeighborList = new NeighborList(null,
1504 alchemicalParameters.vaporCrystal, atoms, Double.POSITIVE_INFINITY, 2.0, parallelTeam);
1505 vacuumNeighborList.setIntermolecular(false);
1506
1507 alchemicalParameters.vaporLists = new int[1][nAtoms][];
1508 double[][] coords = new double[1][nAtoms * 3];
1509 for (int i = 0; i < nAtoms; i++) {
1510 coords[0][i * 3] = atoms[i].getX();
1511 coords[0][i * 3 + 1] = atoms[i].getY();
1512 coords[0][i * 3 + 2] = atoms[i].getZ();
1513 }
1514 boolean print = logger.isLoggable(Level.FINE);
1515 vacuumNeighborList.buildList(coords, alchemicalParameters.vaporLists, nnAtoms, true, print);
1516 alchemicalParameters.vaporPermanentSchedule = vacuumNeighborList.getPairwiseSchedule();
1517 alchemicalParameters.vaporEwaldSchedule = alchemicalParameters.vaporPermanentSchedule;
1518 alchemicalParameters.vacuumRanges = new Range[maxThreads];
1519 vacuumNeighborList.setDisableUpdates(forceField.getBoolean("DISABLE_NEIGHBOR_UPDATES", false));
1520 } else {
1521 alchemicalParameters.vaporCrystal = null;
1522 alchemicalParameters.vaporLists = null;
1523 alchemicalParameters.vaporPermanentSchedule = null;
1524 alchemicalParameters.vaporEwaldSchedule = null;
1525 alchemicalParameters.vacuumRanges = null;
1526 }
1527 }
1528
1529
1530
1531
1532
1533
1534
1535 private double condensedEnergy() {
1536 if (lambda < alchemicalParameters.polLambdaStart) {
1537
1538
1539
1540
1541
1542 alchemicalParameters.polarizationScale = 0.0;
1543 alchemicalParameters.doPolarization = false;
1544 } else if (lambda <= alchemicalParameters.polLambdaEnd) {
1545 alchemicalParameters.polarizationScale = alchemicalParameters.lPowPol;
1546 alchemicalParameters.doPolarization = true;
1547 } else {
1548 alchemicalParameters.polarizationScale = 1.0;
1549 alchemicalParameters.doPolarization = true;
1550 }
1551 alchemicalParameters.doPermanentRealSpace = true;
1552 alchemicalParameters.permanentScale = alchemicalParameters.lPowPerm;
1553 alchemicalParameters.dEdLSign = 1.0;
1554
1555 return computeEnergy(false);
1556 }
1557
1558
1559
1560
1561
1562
1563
1564 private double condensedNoLigandSCF() {
1565
1566 boolean skip = true;
1567 for (int i = 0; i < nAtoms; i++) {
1568 if (atoms[i].applyLambda()) {
1569 use[i] = false;
1570 } else {
1571 use[i] = true;
1572 skip = false;
1573 }
1574 }
1575
1576
1577 alchemicalParameters.doPermanentRealSpace = false;
1578 alchemicalParameters.permanentScale =
1579 1.0 - alchemicalParameters.lPowPerm;
1580 alchemicalParameters.dEdLSign = -1.0;
1581
1582
1583
1584 if (lambda <= alchemicalParameters.polLambdaEnd
1585 && alchemicalParameters.doNoLigandCondensedSCF) {
1586 alchemicalParameters.doPolarization = true;
1587 alchemicalParameters.polarizationScale =
1588 1.0 - alchemicalParameters.lPowPol;
1589 } else {
1590 alchemicalParameters.doPolarization = false;
1591 alchemicalParameters.polarizationScale = 0.0;
1592 }
1593
1594
1595 boolean gkBack = generalizedKirkwoodTerm;
1596 generalizedKirkwoodTerm = false;
1597
1598
1599
1600 double energy;
1601 if (skip) {
1602 energy = permanentMultipoleEnergy + polarizationEnergy + solvationEnergy;
1603 } else {
1604 energy = computeEnergy(false);
1605 Arrays.fill(use, true);
1606 }
1607
1608 generalizedKirkwoodTerm = gkBack;
1609
1610 return energy;
1611 }
1612
1613
1614
1615
1616
1617
1618 private double ligandElec() {
1619 for (int i = 0; i < nAtoms; i++) {
1620 use[i] = atoms[i].applyLambda();
1621 }
1622
1623
1624
1625 alchemicalParameters.doPermanentRealSpace = true;
1626 alchemicalParameters.permanentScale = 1.0 - alchemicalParameters.lPowPerm;
1627 alchemicalParameters.dEdLSign = -1.0;
1628 double lAlphaBack = alchemicalParameters.lAlpha;
1629 double dlAlphaBack = alchemicalParameters.dlAlpha;
1630 double d2lAlphaBack = alchemicalParameters.d2lAlpha;
1631 alchemicalParameters.lAlpha = 0.0;
1632 alchemicalParameters.dlAlpha = 0.0;
1633 alchemicalParameters.d2lAlpha = 0.0;
1634
1635
1636
1637 if (lambda <= alchemicalParameters.polLambdaEnd) {
1638 alchemicalParameters.doPolarization = true;
1639 alchemicalParameters.polarizationScale = 1.0 - alchemicalParameters.lPowPol;
1640 } else {
1641 alchemicalParameters.doPolarization = false;
1642 alchemicalParameters.polarizationScale = 0.0;
1643 }
1644
1645
1646 double offBack = ewaldParameters.off;
1647 double aewaldBack = ewaldParameters.aewald;
1648 ewaldParameters.setEwaldParameters(Double.MAX_VALUE, 0.0);
1649
1650 IntegerSchedule permanentScheduleBack = permanentSchedule;
1651 IntegerSchedule ewaldScheduleBack = realSpaceNeighborParameters.realSpaceSchedule;
1652 Range[] rangesBack = realSpaceNeighborParameters.realSpaceRanges;
1653 permanentSchedule = alchemicalParameters.vaporPermanentSchedule;
1654 realSpaceNeighborParameters.realSpaceSchedule = alchemicalParameters.vaporEwaldSchedule;
1655 realSpaceNeighborParameters.realSpaceRanges = alchemicalParameters.vacuumRanges;
1656
1657
1658 Crystal crystalBack = crystal;
1659 int nSymmBack = nSymm;
1660 int[][][] listsBack = neighborLists;
1661 neighborLists = alchemicalParameters.vaporLists;
1662 crystal = alchemicalParameters.vaporCrystal;
1663 nSymm = 1;
1664
1665
1666 boolean gkBack = generalizedKirkwoodTerm;
1667
1668
1669 SCFAlgorithm scfBack = scfAlgorithm;
1670 scfAlgorithm = SCFAlgorithm.SOR;
1671
1672 if (alchemicalParameters.doLigandGKElec) {
1673 generalizedKirkwoodTerm = true;
1674 generalizedKirkwood.setNeighborList(alchemicalParameters.vaporLists);
1675 generalizedKirkwood.setLambda(lambda);
1676 generalizedKirkwood.setCutoff(ewaldParameters.off);
1677 generalizedKirkwood.setCrystal(alchemicalParameters.vaporCrystal);
1678 generalizedKirkwood.setLambdaFunction(alchemicalParameters.polarizationScale,
1679 alchemicalParameters.dEdLSign * alchemicalParameters.dlPowPol,
1680 alchemicalParameters.dEdLSign * alchemicalParameters.d2lPowPol);
1681 } else {
1682 generalizedKirkwoodTerm = false;
1683 }
1684
1685 double energy = computeEnergy(false);
1686
1687
1688 ewaldParameters.setEwaldParameters(offBack, aewaldBack);
1689 neighborLists = listsBack;
1690 crystal = crystalBack;
1691 nSymm = nSymmBack;
1692 permanentSchedule = permanentScheduleBack;
1693 realSpaceNeighborParameters.realSpaceSchedule = ewaldScheduleBack;
1694 realSpaceNeighborParameters.realSpaceRanges = rangesBack;
1695 alchemicalParameters.lAlpha = lAlphaBack;
1696 alchemicalParameters.dlAlpha = dlAlphaBack;
1697 alchemicalParameters.d2lAlpha = d2lAlphaBack;
1698 generalizedKirkwoodTerm = gkBack;
1699 scfAlgorithm = scfBack;
1700
1701 fill(use, true);
1702
1703 return energy;
1704 }
1705
1706
1707
1708
1709
1710
1711 private double nnElec() {
1712 for (int i = 0; i < nAtoms; i++) {
1713 use[i] = atoms[i].isNeuralNetwork();
1714 }
1715
1716
1717
1718 alchemicalParameters.doPermanentRealSpace = true;
1719 alchemicalParameters.permanentScale = -1.0;
1720 alchemicalParameters.dEdLSign = 1.0;
1721 double lAlphaBack = alchemicalParameters.lAlpha;
1722 double dlAlphaBack = alchemicalParameters.dlAlpha;
1723 double d2lAlphaBack = alchemicalParameters.d2lAlpha;
1724 alchemicalParameters.lAlpha = 0.0;
1725 alchemicalParameters.dlAlpha = 0.0;
1726 alchemicalParameters.d2lAlpha = 0.0;
1727 alchemicalParameters.doPolarization = true;
1728 alchemicalParameters.polarizationScale = -1.0;
1729
1730
1731 double offBack = ewaldParameters.off;
1732 double aewaldBack = ewaldParameters.aewald;
1733 ewaldParameters.setEwaldParameters(Double.POSITIVE_INFINITY, 0.0);
1734
1735 IntegerSchedule permanentScheduleBack = permanentSchedule;
1736 IntegerSchedule ewaldScheduleBack = realSpaceNeighborParameters.realSpaceSchedule;
1737 Range[] rangesBack = realSpaceNeighborParameters.realSpaceRanges;
1738 permanentSchedule = alchemicalParameters.vaporPermanentSchedule;
1739 realSpaceNeighborParameters.realSpaceSchedule = alchemicalParameters.vaporEwaldSchedule;
1740 realSpaceNeighborParameters.realSpaceRanges = alchemicalParameters.vacuumRanges;
1741
1742
1743 Crystal crystalBack = crystal;
1744 int nSymmBack = nSymm;
1745 int[][][] listsBack = neighborLists;
1746 neighborLists = alchemicalParameters.vaporLists;
1747 crystal = alchemicalParameters.vaporCrystal;
1748 nSymm = 1;
1749
1750
1751 boolean gkBack = generalizedKirkwoodTerm;
1752
1753
1754 SCFAlgorithm scfBack = scfAlgorithm;
1755 scfAlgorithm = SCFAlgorithm.SOR;
1756 generalizedKirkwoodTerm = false;
1757
1758 double energy = computeEnergy(false);
1759
1760
1761 ewaldParameters.setEwaldParameters(offBack, aewaldBack);
1762 neighborLists = listsBack;
1763 crystal = crystalBack;
1764 nSymm = nSymmBack;
1765 permanentSchedule = permanentScheduleBack;
1766 realSpaceNeighborParameters.realSpaceSchedule = ewaldScheduleBack;
1767 realSpaceNeighborParameters.realSpaceRanges = rangesBack;
1768 alchemicalParameters.lAlpha = lAlphaBack;
1769 alchemicalParameters.dlAlpha = dlAlphaBack;
1770 alchemicalParameters.d2lAlpha = d2lAlphaBack;
1771 generalizedKirkwoodTerm = gkBack;
1772 scfAlgorithm = scfBack;
1773
1774 fill(use, true);
1775
1776 return energy;
1777 }
1778
1779
1780
1781
1782
1783
1784
1785 private double computeEnergy(boolean print) {
1786
1787 permanentMultipoleField();
1788
1789
1790 if (generalizedKirkwoodTerm) {
1791 pmeTimings.bornRadiiTotal -= System.nanoTime();
1792 generalizedKirkwood.setUse(use);
1793 generalizedKirkwood.computeBornRadii();
1794 pmeTimings.bornRadiiTotal += System.nanoTime();
1795 }
1796
1797
1798 if (polarization != Polarization.NONE && alchemicalParameters.doPolarization) {
1799
1800 if (generalizedKirkwoodTerm) {
1801 generalizedKirkwoodTerm = false;
1802
1803 selfConsistentField(logger.isLoggable(Level.FINE));
1804
1805 saveInducedDipolesToVacuumDipoles();
1806
1807 generalizedKirkwoodTerm = true;
1808
1809 permanentMultipoleField();
1810 }
1811
1812
1813 selfConsistentField(logger.isLoggable(Level.FINE));
1814
1815 if (esvTerm && polarization != Polarization.NONE) {
1816 for (int i = 0; i < nAtoms; i++) {
1817 if (extendedSystem.isTitrating(i)
1818 && (extendedSystem.isTitratingHydrogen(i) || extendedSystem.isTitratingHeavy(i))) {
1819 double dx = field.getX(i);
1820 double dy = field.getY(i);
1821 double dz = field.getZ(i);
1822 double dxCR = fieldCR.getX(i);
1823 double dyCR = fieldCR.getY(i);
1824 double dzCR = fieldCR.getZ(i);
1825
1826
1827 if (polarization == Polarization.MUTUAL) {
1828 dx += directField[i][0];
1829 dy += directField[i][1];
1830 dz += directField[i][2];
1831 dxCR += directFieldCR[i][0];
1832 dyCR += directFieldCR[i][1];
1833 dzCR += directFieldCR[i][2];
1834 }
1835 double fix = dx * dPolardTitrationESV[i] * dxCR;
1836 double fiy = dy * dPolardTitrationESV[i] * dyCR;
1837 double fiz = dz * dPolardTitrationESV[i] * dzCR;
1838 double titrdUdL = fix + fiy + fiz;
1839 double tautdUdL = 0.0;
1840 if (extendedSystem.isTautomerizing(i)) {
1841 fix = dx * dPolardTautomerESV[i] * dxCR;
1842 fiy = dy * dPolardTautomerESV[i] * dyCR;
1843 fiz = dz * dPolardTautomerESV[i] * dzCR;
1844 tautdUdL = fix + fiy + fiz;
1845 }
1846 extendedSystem.addIndElecDeriv(i, titrdUdL * electric * -0.5, tautdUdL * electric * -0.5);
1847 }
1848 }
1849 }
1850
1851 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
1852 if (gradient && polarization == Polarization.DIRECT) {
1853 reciprocalSpace.splineInducedDipoles(inducedDipole, inducedDipoleCR, use);
1854 field.reset(parallelTeam);
1855 fieldCR.reset(parallelTeam);
1856 inducedDipoleFieldRegion.init(
1857 atoms, crystal, use, molecule,
1858 ipdamp, thole, coordinates, realSpaceNeighborParameters,
1859 inducedDipole, inducedDipoleCR,
1860 reciprocalSpaceTerm, reciprocalSpace,
1861 lambdaMode, ewaldParameters,
1862 field, fieldCR, pmeTimings);
1863 inducedDipoleFieldRegion.executeWith(sectionTeam);
1864 reciprocalSpace.computeInducedPhi(
1865 cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
1866 fractionalInducedDipolePhi, fractionalInducedDipolePhiCR);
1867 }
1868 }
1869
1870 if (scfPredictorParameters.scfPredictor != SCFPredictor.NONE) {
1871 scfPredictorParameters.saveMutualInducedDipoles(lambdaMode,
1872 inducedDipole, inducedDipoleCR, directDipole, directDipoleCR);
1873 }
1874 }
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885 double[][][] inputDipole = inducedDipole;
1886 double[][][] inputDipoleCR = inducedDipoleCR;
1887 double[][] inputPhi = cartesianInducedDipolePhi;
1888 double[][] inputPhiCR = cartesianInducedDipolePhiCR;
1889 double[][] fracInputPhi = fractionalInducedDipolePhi;
1890 double[][] fracInputPhiCR = fractionalInducedDipolePhiCR;
1891 if (generalizedKirkwoodTerm) {
1892 inputDipole = vacuumInducedDipole;
1893 inputDipoleCR = vacuumInducedDipoleCR;
1894 inputPhi = cartesianVacuumDipolePhi;
1895 inputPhiCR = cartesianVacuumDipolePhiCR;
1896 fracInputPhi = fractionalVacuumDipolePhi;
1897 fracInputPhiCR = fractionalVacuumDipolePhiCR;
1898 }
1899
1900 double eself = 0.0;
1901 double erecip = 0.0;
1902 double eselfi = 0.0;
1903 double erecipi = 0.0;
1904 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
1905 reciprocalEnergyRegion.init(atoms, crystal, gradient, lambdaTerm, esvTerm, use,
1906 globalMultipole, fractionalMultipole, dMultipoledTirationESV, dMultipoledTautomerESV,
1907 cartesianMultipolePhi, fracMultipolePhi,
1908 polarization, inputDipole, inputDipoleCR, inputPhi, inputPhiCR, fracInputPhi,
1909 fracInputPhiCR, reciprocalSpace, alchemicalParameters, extendedSystem,
1910
1911 grad, torque, lambdaGrad, lambdaTorque, shareddEdLambda, sharedd2EdLambda2);
1912 reciprocalEnergyRegion.executeWith(parallelTeam);
1913 eself = reciprocalEnergyRegion.getPermanentSelfEnergy();
1914 erecip = reciprocalEnergyRegion.getPermanentReciprocalEnergy();
1915 eselfi = reciprocalEnergyRegion.getInducedDipoleSelfEnergy();
1916 erecipi = reciprocalEnergyRegion.getInducedDipoleReciprocalEnergy();
1917 interactions += nAtoms;
1918 }
1919
1920 pmeTimings.realSpaceEnergyTotal -= System.nanoTime();
1921 realSpaceEnergyRegion.init(atoms, crystal, extendedSystem, esvTerm, coordinates, frame, axisAtom,
1922 globalMultipole, dMultipoledTirationESV, dMultipoledTautomerESV,
1923 inputDipole, inputDipoleCR, use, molecule,
1924 ip11, mask12, mask13, mask14, mask15, isSoft, ipdamp, thole, realSpaceNeighborParameters,
1925 gradient, lambdaTerm, nnTerm, lambdaMode, polarization,
1926 ewaldParameters, scaleParameters, alchemicalParameters,
1927 pmeTimings.realSpaceEnergyTime,
1928
1929 grad, torque, lambdaGrad, lambdaTorque, shareddEdLambda, sharedd2EdLambda2);
1930 realSpaceEnergyRegion.executeWith(parallelTeam);
1931 double ereal = realSpaceEnergyRegion.getPermanentEnergy();
1932 double ereali = realSpaceEnergyRegion.getPolarizationEnergy();
1933 interactions += realSpaceEnergyRegion.getInteractions();
1934 pmeTimings.realSpaceEnergyTotal += System.nanoTime();
1935
1936 if (generalizedKirkwoodTerm) {
1937
1938
1939 double eGK = 0.0;
1940
1941
1942 AlchemicalParameters alchemicalParametersGK = new AlchemicalParameters(
1943 forceField, false, false, polarization);
1944
1945 alchemicalParametersGK.permanentScale = 0.0;
1946 alchemicalParametersGK.doPermanentRealSpace = false;
1947
1948 alchemicalParametersGK.polarizationScale = -1.0;
1949
1950
1951 AtomicDoubleArray3D gradGK = generalizedKirkwood.getGrad();
1952 AtomicDoubleArray3D torqueGK = generalizedKirkwood.getTorque();
1953
1954 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
1955 reciprocalEnergyRegion.init(
1956 atoms, crystal, gradient, false, esvTerm,
1957 use, globalMultipole, fractionalMultipole,
1958 dMultipoledTirationESV, dMultipoledTautomerESV,
1959 cartesianMultipolePhi, fracMultipolePhi,
1960 polarization,
1961 vacuumInducedDipole, vacuumInducedDipoleCR,
1962 cartesianVacuumDipolePhi, cartesianVacuumDipolePhiCR,
1963 fractionalVacuumDipolePhi, fractionalVacuumDipolePhiCR,
1964 reciprocalSpace, alchemicalParametersGK,
1965 extendedSystem,
1966 gradGK, torqueGK, null, null,
1967 shareddEdLambda, sharedd2EdLambda2);
1968 reciprocalEnergyRegion.executeWith(parallelTeam);
1969 eGK = reciprocalEnergyRegion.getInducedDipoleSelfEnergy()
1970 + reciprocalEnergyRegion.getInducedDipoleReciprocalEnergy();
1971 }
1972
1973 pmeTimings.gkEnergyTotal -= System.nanoTime();
1974 realSpaceEnergyRegion.init(atoms, crystal,
1975 extendedSystem, esvTerm,
1976 coordinates, frame, axisAtom,
1977 globalMultipole, dMultipoledTirationESV, dMultipoledTautomerESV,
1978 vacuumInducedDipole, vacuumInducedDipoleCR,
1979 use, molecule,
1980 ip11, mask12, mask13, mask14, mask15, isSoft, ipdamp, thole,
1981 realSpaceNeighborParameters,
1982 gradient, false, nnTerm, lambdaMode,
1983 polarization, ewaldParameters, scaleParameters, alchemicalParametersGK,
1984 pmeTimings.realSpaceEnergyTime,
1985
1986 gradGK,
1987 torqueGK,
1988 null,
1989 null,
1990 shareddEdLambda,
1991 sharedd2EdLambda2);
1992 realSpaceEnergyRegion.executeWith(parallelTeam);
1993 eGK += realSpaceEnergyRegion.getPolarizationEnergy();
1994
1995
1996 alchemicalParametersGK.polarizationScale = 1.0;
1997
1998 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
1999 reciprocalEnergyRegion.init(
2000 atoms, crystal, gradient, false, esvTerm, use,
2001 globalMultipole, fractionalMultipole,
2002 dMultipoledTirationESV, dMultipoledTautomerESV,
2003 cartesianMultipolePhi, fracMultipolePhi,
2004 polarization, inducedDipole, inducedDipoleCR,
2005 cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
2006 fractionalInducedDipolePhi, fractionalInducedDipolePhiCR,
2007 reciprocalSpace, alchemicalParametersGK, extendedSystem,
2008 gradGK, torqueGK, null, null,
2009 shareddEdLambda, sharedd2EdLambda2);
2010 reciprocalEnergyRegion.executeWith(parallelTeam);
2011 eGK += reciprocalEnergyRegion.getInducedDipoleSelfEnergy()
2012 + reciprocalEnergyRegion.getInducedDipoleReciprocalEnergy();
2013 }
2014
2015 pmeTimings.gkEnergyTotal -= System.nanoTime();
2016 realSpaceEnergyRegion.init(
2017 atoms,
2018 crystal,
2019 extendedSystem,
2020 esvTerm,
2021 coordinates,
2022 frame,
2023 axisAtom,
2024 globalMultipole,
2025 dMultipoledTirationESV,
2026 dMultipoledTautomerESV,
2027 inducedDipole,
2028 inducedDipoleCR,
2029 use,
2030 molecule,
2031 ip11,
2032 mask12,
2033 mask13,
2034 mask14,
2035 mask15,
2036 isSoft,
2037 ipdamp,
2038 thole,
2039 realSpaceNeighborParameters,
2040 gradient,
2041 false,
2042 nnTerm,
2043 lambdaMode,
2044 polarization,
2045 ewaldParameters,
2046 scaleParameters,
2047 alchemicalParametersGK,
2048 pmeTimings.realSpaceEnergyTime,
2049
2050 gradGK,
2051 torqueGK,
2052 null,
2053 null,
2054 shareddEdLambda,
2055 sharedd2EdLambda2);
2056 realSpaceEnergyRegion.executeWith(parallelTeam);
2057 eGK += realSpaceEnergyRegion.getPolarizationEnergy();
2058
2059
2060 solvationEnergy += generalizedKirkwood.solvationEnergy(eGK, gradient, print);
2061 if (gradient) {
2062
2063 generalizedKirkwood.reduce(grad, torque, lambdaGrad, lambdaTorque);
2064 }
2065 gkInteractions += generalizedKirkwood.getInteractions();
2066 pmeTimings.gkEnergyTotal += System.nanoTime();
2067 }
2068
2069
2070 permanentRealSpaceEnergy += ereal;
2071 permanentSelfEnergy += eself;
2072 permanentReciprocalEnergy += erecip;
2073 inducedRealSpaceEnergy += ereali;
2074 inducedSelfEnergy += eselfi;
2075 inducedReciprocalEnergy += erecipi;
2076 permanentMultipoleEnergy += eself + erecip + ereal;
2077 polarizationEnergy += eselfi + erecipi + ereali;
2078 totalMultipoleEnergy += ereal + eself + erecip + ereali + eselfi + erecipi;
2079
2080
2081 if (logger.isLoggable(Level.FINE)) {
2082 StringBuilder sb = new StringBuilder();
2083 sb.append(format("\n Global Cartesian PME, lambdaMode=%s\n", lambdaMode.toString()));
2084 sb.append(format(" Multipole Self-Energy: %16.8f\n", eself));
2085 sb.append(format(" Multipole Reciprocal: %16.8f\n", erecip));
2086 sb.append(format(" Multipole Real Space: %16.8f\n", ereal));
2087 sb.append(format(" Polarization Self-Energy:%16.8f\n", eselfi));
2088 sb.append(format(" Polarization Reciprocal: %16.8f\n", erecipi));
2089 sb.append(format(" Polarization Real Space: %16.8f\n", ereali));
2090 if (generalizedKirkwoodTerm) {
2091 sb.append(format(" Generalized Kirkwood: %16.8f\n", solvationEnergy));
2092 }
2093 logger.info(sb.toString());
2094 }
2095
2096 return permanentMultipoleEnergy + polarizationEnergy + solvationEnergy;
2097 }
2098
2099
2100
2101
2102 private void permanentMultipoleField() {
2103 try {
2104
2105 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2106 reciprocalSpace.computeBSplines();
2107 reciprocalSpace.splinePermanentMultipoles(globalMultipole, fractionalMultipole, use);
2108 }
2109
2110 field.reset(parallelTeam);
2111 fieldCR.reset(parallelTeam);
2112 permanentFieldRegion.init(atoms, crystal, coordinates,
2113 globalMultipole, inducedDipole, inducedDipoleCR, neighborLists,
2114 scaleParameters, use, molecule, ipdamp, thole, ip11,
2115 mask12, mask13, mask14, lambdaMode, reciprocalSpaceTerm, reciprocalSpace,
2116 ewaldParameters, pcgSolver, permanentSchedule, realSpaceNeighborParameters,
2117 field, fieldCR, pmeTimings);
2118
2119
2120 sectionTeam.execute(permanentFieldRegion);
2121
2122 pmeTimings.realSpacePermTotal = permanentFieldRegion.getRealSpacePermTotal();
2123
2124
2125 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2126 reciprocalSpace.computePermanentPhi(cartesianMultipolePhi, fracMultipolePhi);
2127 }
2128 } catch (RuntimeException e) {
2129 String message = "Fatal exception computing the permanent multipole field.\n";
2130 logger.log(Level.WARNING, message, e);
2131 throw e;
2132 } catch (Exception e) {
2133 String message = "Fatal exception computing the permanent multipole field.\n";
2134 logger.log(Level.SEVERE, message, e);
2135 }
2136 }
2137
2138 private void saveInducedDipolesToVacuumDipoles() {
2139 for (int i = 0; i < nAtoms; i++) {
2140 arraycopy(directDipole[i], 0, vacuumDirectDipole[i], 0, 3);
2141 arraycopy(directDipoleCR[i], 0, vacuumDirectDipoleCR[i], 0, 3);
2142 arraycopy(inducedDipole[0][i], 0, vacuumInducedDipole[0][i], 0, 3);
2143 arraycopy(inducedDipoleCR[0][i], 0, vacuumInducedDipoleCR[0][i], 0, 3);
2144 arraycopy(cartesianInducedDipolePhi[i], 0, cartesianVacuumDipolePhi[i], 0, tensorCount);
2145 arraycopy(cartesianInducedDipolePhiCR[i], 0, cartesianVacuumDipolePhiCR[i], 0,
2146 tensorCount);
2147 arraycopy(fractionalInducedDipolePhi[i], 0, fractionalVacuumDipolePhi[i], 0,
2148 tensorCount);
2149 arraycopy(fractionalInducedDipolePhiCR[i], 0, fractionalVacuumDipolePhiCR[i], 0,
2150 tensorCount);
2151 if (nSymm > 1) {
2152 for (int s = 1; s < nSymm; s++) {
2153 arraycopy(inducedDipole[s][i], 0, vacuumInducedDipole[s][i], 0, 3);
2154 arraycopy(inducedDipoleCR[s][i], 0, vacuumInducedDipoleCR[s][i], 0, 3);
2155 }
2156 }
2157 }
2158 }
2159
2160
2161
2162
2163 private int selfConsistentField(boolean print) {
2164 if (polarization == Polarization.NONE) {
2165 return -1;
2166 }
2167 long startTime = System.nanoTime();
2168
2169
2170 if (generalizedKirkwoodTerm) {
2171 pmeTimings.gkEnergyTotal = -System.nanoTime();
2172 generalizedKirkwood.computePermanentGKField();
2173 pmeTimings.gkEnergyTotal += System.nanoTime();
2174 logger.fine(
2175 format(" Computed GK permanent field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
2176 }
2177 directRegion.init(
2178 atoms,
2179 polarizability,
2180 globalMultipole,
2181 cartesianMultipolePhi,
2182 field,
2183 fieldCR,
2184 generalizedKirkwoodTerm,
2185 generalizedKirkwood,
2186 ewaldParameters,
2187 soluteDielectric,
2188 inducedDipole,
2189 inducedDipoleCR,
2190 directDipole,
2191 directDipoleCR,
2192 directField,
2193 directFieldCR
2194 );
2195 directRegion.executeWith(parallelTeam);
2196
2197
2198 if (polarization != Polarization.MUTUAL) {
2199 expandInducedDipoles();
2200 return 0;
2201 }
2202
2203
2204 if (scfPredictorParameters.scfPredictor != SCFPredictor.NONE) {
2205 switch (scfPredictorParameters.scfPredictor) {
2206 case ASPC -> scfPredictorParameters.aspcPredictor(lambdaMode, inducedDipole, inducedDipoleCR);
2207 case LS -> scfPredictorParameters.leastSquaresPredictor(lambdaMode, inducedDipole, inducedDipoleCR);
2208 case POLY -> scfPredictorParameters.polynomialPredictor(lambdaMode, inducedDipole, inducedDipoleCR);
2209 default -> {
2210 }
2211 }
2212 }
2213
2214
2215 expandInducedDipoles();
2216
2217
2218 try {
2219 return switch (scfAlgorithm) {
2220 case SOR -> scfBySOR(print, startTime);
2221 case EPT -> scfByEPT(print, startTime);
2222 default -> {
2223
2224 pcgSolver.init(atoms, coordinates, polarizability, ipdamp, thole,
2225 use, crystal, inducedDipole, inducedDipoleCR, directDipole, directDipoleCR,
2226 field, fieldCR, ewaldParameters, soluteDielectric, parallelTeam,
2227 realSpaceNeighborParameters.realSpaceSchedule,
2228 pmeTimings.realSpaceSCFTime);
2229 yield pcgSolver.scfByPCG(print, startTime, this);
2230 }
2231 };
2232 } catch (EnergyException ex) {
2233 if (directFallback) {
2234
2235 logger.warning(ex.toString());
2236
2237 if (generalizedKirkwoodTerm) {
2238 pmeTimings.gkEnergyTotal = -System.nanoTime();
2239 generalizedKirkwood.computePermanentGKField();
2240 pmeTimings.gkEnergyTotal += System.nanoTime();
2241 logger.fine(
2242 format(" Computed GK permanent field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
2243 }
2244 directRegion.init(
2245 atoms,
2246 polarizability,
2247 globalMultipole,
2248 cartesianMultipolePhi,
2249 field,
2250 fieldCR,
2251 generalizedKirkwoodTerm,
2252 generalizedKirkwood,
2253 ewaldParameters,
2254 soluteDielectric,
2255 inducedDipole,
2256 inducedDipoleCR,
2257 directDipole,
2258 directDipoleCR,
2259 directField,
2260 directFieldCR
2261 );
2262 directRegion.executeWith(parallelTeam);
2263 expandInducedDipoles();
2264 logger.info(" Direct induced dipoles computed due to SCF failure.");
2265 return 0;
2266 } else {
2267 throw ex;
2268 }
2269 }
2270 }
2271
2272
2273
2274
2275 private int scfBySOR(boolean print, long startTime) {
2276 long directTime = System.nanoTime() - startTime;
2277
2278
2279 StringBuilder sb = null;
2280 if (print) {
2281 sb = new StringBuilder("\n Self-Consistent Field\n Iter RMS Change (Debye) Time\n");
2282 }
2283 int completedSCFCycles = 0;
2284 int maxSCFCycles = 1000;
2285 double eps = 100.0;
2286 double previousEps;
2287 boolean done = false;
2288 while (!done) {
2289 long cycleTime = -System.nanoTime();
2290 try {
2291 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2292 reciprocalSpace.splineInducedDipoles(inducedDipole, inducedDipoleCR, use);
2293 }
2294 field.reset(parallelTeam);
2295 fieldCR.reset(parallelTeam);
2296 inducedDipoleFieldRegion.init(
2297 atoms, crystal, use, molecule, ipdamp, thole, coordinates,
2298 realSpaceNeighborParameters, inducedDipole, inducedDipoleCR,
2299 reciprocalSpaceTerm, reciprocalSpace, lambdaMode, ewaldParameters,
2300 field, fieldCR, pmeTimings);
2301 inducedDipoleFieldRegion.executeWith(sectionTeam);
2302 pmeTimings.realSpaceSCFTotal = inducedDipoleFieldRegion.getRealSpaceSCFTotal();
2303 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2304 reciprocalSpace.computeInducedPhi(
2305 cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
2306 fractionalInducedDipolePhi, fractionalInducedDipolePhiCR);
2307 }
2308
2309 if (generalizedKirkwoodTerm) {
2310
2311 pmeTimings.gkEnergyTotal = -System.nanoTime();
2312 generalizedKirkwood.computeInducedGKField();
2313 pmeTimings.gkEnergyTotal += System.nanoTime();
2314 logger.fine(
2315 format(" Computed GK induced field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
2316 }
2317
2318 sorRegion.init(
2319 atoms,
2320 polarizability,
2321 inducedDipole,
2322 inducedDipoleCR,
2323 directDipole,
2324 directDipoleCR,
2325 cartesianInducedDipolePhi,
2326 cartesianInducedDipolePhiCR,
2327 field,
2328 fieldCR,
2329 generalizedKirkwoodTerm,
2330 generalizedKirkwood,
2331 ewaldParameters);
2332 parallelTeam.execute(sorRegion);
2333
2334 expandInducedDipoles();
2335 } catch (Exception e) {
2336 String message = "Exception computing mutual induced dipoles.";
2337 logger.log(Level.SEVERE, message, e);
2338 }
2339 completedSCFCycles++;
2340 previousEps = eps;
2341 eps = sorRegion.getEps();
2342 eps = Constants.ELEC_ANG_TO_DEBYE * sqrt(eps / (double) nAtoms);
2343 cycleTime += System.nanoTime();
2344 if (print) {
2345 sb.append(format(" %4d %15.10f %7.4f\n", completedSCFCycles, eps, cycleTime * NS2SEC));
2346 }
2347
2348
2349 if (eps > previousEps) {
2350 if (sb != null) {
2351 logger.warning(sb.toString());
2352 }
2353 String message = format(" SCF convergence failure: (%10.5f > %10.5f)\n", eps, previousEps);
2354 throw new EnergyException(message);
2355 }
2356
2357
2358
2359 if (completedSCFCycles >= maxSCFCycles) {
2360 if (sb != null) {
2361 logger.warning(sb.toString());
2362 }
2363 String message = format(" Maximum SCF iterations reached: (%d)\n", completedSCFCycles);
2364 throw new EnergyException(message);
2365 }
2366
2367
2368 if (eps < poleps) {
2369 done = true;
2370 }
2371 }
2372 if (print) {
2373 sb.append(format(" Direct: %7.4f\n", NS2SEC * directTime));
2374 startTime = System.nanoTime() - startTime;
2375 sb.append(format(" Total: %7.4f", startTime * NS2SEC));
2376 logger.info(sb.toString());
2377 }
2378 return completedSCFCycles;
2379 }
2380
2381
2382
2383
2384 private int scfByEPT(boolean print, long startTime) {
2385
2386
2387 for (int i = 0; i < nAtoms; i++) {
2388 for (int j = 0; j < 3; j++) {
2389 optRegion.optDipole[0][i][j] = directDipole[i][j];
2390 optRegion.optDipoleCR[0][i][j] = directDipoleCR[i][j];
2391 }
2392 }
2393
2394 int optOrder = optRegion.optOrder;
2395
2396 for (int currentOptOrder = 1; currentOptOrder <= optOrder; currentOptOrder++) {
2397 try {
2398 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2399 reciprocalSpace.splineInducedDipoles(inducedDipole, inducedDipoleCR, use);
2400 }
2401 field.reset(parallelTeam);
2402 fieldCR.reset(parallelTeam);
2403 inducedDipoleFieldRegion.init(
2404 atoms, crystal, use, molecule, ipdamp, thole, coordinates, realSpaceNeighborParameters,
2405 inducedDipole, inducedDipoleCR, reciprocalSpaceTerm, reciprocalSpace, lambdaMode,
2406 ewaldParameters, field, fieldCR, pmeTimings);
2407 inducedDipoleFieldRegion.executeWith(sectionTeam);
2408 pmeTimings.realSpaceSCFTotal = inducedDipoleFieldRegion.getRealSpaceSCFTotal();
2409
2410 if (reciprocalSpaceTerm && ewaldParameters.aewald > 0.0) {
2411 reciprocalSpace.computeInducedPhi(
2412 cartesianInducedDipolePhi, cartesianInducedDipolePhiCR,
2413 fractionalInducedDipolePhi, fractionalInducedDipolePhiCR);
2414 }
2415
2416 if (generalizedKirkwoodTerm) {
2417
2418 pmeTimings.gkEnergyTotal = -System.nanoTime();
2419 generalizedKirkwood.computeInducedGKField();
2420 pmeTimings.gkEnergyTotal += System.nanoTime();
2421 logger.fine(
2422 format(" Computed GK induced field %8.3f (sec)", pmeTimings.gkEnergyTotal * 1.0e-9));
2423 }
2424
2425 optRegion.init(
2426 currentOptOrder,
2427 atoms,
2428 polarizability,
2429 inducedDipole,
2430 inducedDipoleCR,
2431 cartesianInducedDipolePhi,
2432 cartesianInducedDipolePhiCR,
2433 field,
2434 fieldCR,
2435 generalizedKirkwoodTerm,
2436 generalizedKirkwood,
2437 ewaldParameters,
2438 soluteDielectric);
2439 parallelTeam.execute(optRegion);
2440
2441 expandInducedDipoles();
2442
2443 } catch (Exception e) {
2444 String message = "Exception computing opt induced dipoles.";
2445 logger.log(Level.SEVERE, message, e);
2446 }
2447 }
2448
2449
2450 for (int i = 0; i < nAtoms; i++) {
2451 for (int j = 0; j < 3; j++) {
2452 inducedDipole[0][i][j] = 0.0;
2453 inducedDipoleCR[0][i][j] = 0.0;
2454 double sum = 0.0;
2455 double sump = 0.0;
2456 for (int k = 0; k <= optOrder; k++) {
2457 sum += optRegion.optDipole[k][i][j];
2458 sump += optRegion.optDipoleCR[k][i][j];
2459 inducedDipole[0][i][j] += optRegion.optCoefficients[k] * sum;
2460 inducedDipoleCR[0][i][j] += optRegion.optCoefficients[k] * sump;
2461 }
2462 }
2463 }
2464
2465
2466 expandInducedDipoles();
2467
2468 return optOrder;
2469 }
2470
2471
2472
2473
2474
2475
2476 public double getTotalMultipoleEnergy() {
2477 return permanentMultipoleEnergy + polarizationEnergy;
2478 }
2479
2480
2481
2482
2483
2484
2485 public void setPolarization(Polarization set) {
2486 this.polarization = set;
2487 }
2488
2489
2490
2491
2492 private void assignMultipoles() {
2493 if (forceField == null) {
2494 String message = "No force field is defined.\n";
2495 logger.log(Level.SEVERE, message);
2496 return;
2497 }
2498
2499 if (forceField.getForceFieldTypeCount(ForceFieldType.MULTIPOLE) < 1
2500 && forceField.getForceFieldTypeCount(ForceFieldType.CHARGE) < 1) {
2501 String message = "Force field has no permanent electrostatic types.\n";
2502 logger.log(Level.SEVERE, message);
2503 return;
2504 }
2505 if (nAtoms < 1) {
2506 String message = "No atoms are defined.\n";
2507 logger.log(Level.SEVERE, message);
2508 return;
2509 }
2510 for (int i = 0; i < nAtoms; i++) {
2511 Atom atom = atoms[i];
2512
2513 if (!assignMultipole(elecForm, atom, forceField, localMultipole[i], i, axisAtom, frame)) {
2514 logger.info(
2515 format("No MultipoleType could be assigned:\n %s --> %s", atom, atom.getAtomType()));
2516 StringBuilder sb = new StringBuilder();
2517 List<Bond> bonds = atom.getBonds();
2518 for (Bond bond12 : bonds) {
2519 Atom a2 = bond12.get1_2(atom);
2520 AtomType aType2 = a2.getAtomType();
2521 sb.append(format("\n 1-2 %s --> %s", a2, aType2));
2522 }
2523 for (Bond bond12 : bonds) {
2524 Atom atom2 = bond12.get1_2(atom);
2525 bonds = atom2.getBonds();
2526 for (Bond bond23 : bonds) {
2527 Atom a2 = bond23.get1_2(atom2);
2528 AtomType aType2 = a2.getAtomType();
2529 sb.append(format("\n 1-3 %s --> %s", a2, aType2));
2530 }
2531 }
2532
2533 List<MultipoleType> multipoleTypes = forceField.getMultipoleTypes(atom.getAtomType().getKey());
2534 if (multipoleTypes != null && !multipoleTypes.isEmpty()) {
2535 sb.append("\n Similar Multipole types:");
2536 for (MultipoleType multipoleType : multipoleTypes) {
2537 sb.append(format("\n %s", multipoleType));
2538 }
2539 }
2540
2541 logger.log(Level.SEVERE, sb.toString());
2542 }
2543 }
2544
2545
2546 StringBuilder sb = new StringBuilder();
2547 for (int i = 0; i < nAtoms; i++) {
2548 boolean flag = false;
2549 for (int j = 0; j < 10; j++) {
2550 if (Double.isNaN(localMultipole[i][j])) {
2551 flag = true;
2552 break;
2553 }
2554 }
2555 if (flag) {
2556 sb.append("\n").append(atoms[i].toString()).append("\n");
2557 sb.append(format("%d", i + 1));
2558 for (int j = 0; j < 10; j++) {
2559 sb.append(format(" %8.3f", localMultipole[i][j]));
2560 }
2561 sb.append("\n");
2562 }
2563 }
2564 if (!sb.isEmpty()) {
2565 String message = "Fatal exception: Error assigning multipoles. " + sb;
2566 logger.log(Level.SEVERE, message);
2567 }
2568 }
2569
2570
2571
2572
2573 protected void assignPolarizationGroups() {
2574
2575 List<Integer> group = new ArrayList<>();
2576 for (int i = 0; i < nAtoms; i++) {
2577 Atom a = atoms[i];
2578 if (a.getIndex() - 1 != i) {
2579 logger.info(format(" PME Index: %d Atom Index: %d\n Atom: %s",
2580 i, a.getIndex() - 1, a));
2581 logger.severe(" Atom indexing is not consistent in PME.");
2582 }
2583 }
2584 for (Atom ai : atoms) {
2585 group.clear();
2586 int index = ai.getIndex() - 1;
2587 group.add(index);
2588 PolarizeType polarizeType = ai.getPolarizeType();
2589 if (polarizeType != null) {
2590 if (polarizeType.polarizationGroup != null) {
2591 PolarizeType.growGroup(group, ai);
2592 sort(group);
2593 ip11[index] = new int[group.size()];
2594 int j = 0;
2595 for (int k : group) {
2596 ip11[index][j++] = k;
2597 }
2598 } else {
2599 ip11[index] = new int[group.size()];
2600 int j = 0;
2601 for (int k : group) {
2602 ip11[index][j++] = k;
2603 }
2604 }
2605 } else {
2606 String message = "The polarize keyword was not found for atom "
2607 + (index + 1) + " with type " + ai.getType();
2608 logger.severe(message);
2609 }
2610 }
2611
2612 int[] mask = new int[nAtoms];
2613 List<Integer> list = new ArrayList<>();
2614 List<Integer> keep = new ArrayList<>();
2615 for (int i = 0; i < nAtoms; i++) {
2616 mask[i] = -1;
2617 }
2618 for (int i = 0; i < nAtoms; i++) {
2619 list.clear();
2620 for (int j : ip11[i]) {
2621 list.add(j);
2622 mask[j] = i;
2623 }
2624 keep.clear();
2625 for (int j : list) {
2626 Atom aj = atoms[j];
2627 List<Bond> bonds = aj.getBonds();
2628 for (Bond b : bonds) {
2629 Atom ak = b.get1_2(aj);
2630 int k = ak.getIndex() - 1;
2631 if (mask[k] != i) {
2632 keep.add(k);
2633 }
2634 }
2635 }
2636 list.clear();
2637 for (int j : keep) {
2638 for (int k : ip11[j]) {
2639 list.add(k);
2640 }
2641 }
2642 sort(list);
2643 ip12[i] = new int[list.size()];
2644 int j = 0;
2645 for (int k : list) {
2646 ip12[i][j++] = k;
2647 }
2648 }
2649
2650
2651 for (int i = 0; i < nAtoms; i++) {
2652 mask[i] = -1;
2653 }
2654 for (int i = 0; i < nAtoms; i++) {
2655 for (int j : ip11[i]) {
2656 mask[j] = i;
2657 }
2658 for (int j : ip12[i]) {
2659 mask[j] = i;
2660 }
2661 list.clear();
2662 for (int j : ip12[i]) {
2663 for (int k : ip12[j]) {
2664 if (mask[k] != i) {
2665 if (!list.contains(k)) {
2666 list.add(k);
2667 }
2668 }
2669 }
2670 }
2671 ip13[i] = new int[list.size()];
2672 sort(list);
2673 int j = 0;
2674 for (int k : list) {
2675 ip13[i][j++] = k;
2676 }
2677 }
2678 }
2679
2680
2681
2682
2683
2684
2685
2686
2687 public void computeMoments(Atom[] activeAtoms, boolean forceEnergy) {
2688
2689 var netchg = 0.0;
2690 var netdpl = 0.0;
2691 var xdpl = 0.0;
2692 var ydpl = 0.0;
2693 var zdpl = 0.0;
2694 var xxqdp = 0.0;
2695 var xyqdp = 0.0;
2696 var xzqdp = 0.0;
2697 var yxqdp = 0.0;
2698 var yyqdp = 0.0;
2699 var yzqdp = 0.0;
2700 var zxqdp = 0.0;
2701 var zyqdp = 0.0;
2702 var zzqdp = 0.0;
2703
2704
2705 double xmid = 0.0;
2706 double ymid = 0.0;
2707 double zmid = 0.0;
2708 double totalMass = 0;
2709 for (Atom atom : activeAtoms) {
2710 var m = atom.getMass();
2711 totalMass += m;
2712 xmid = xmid + atom.getX() * m;
2713 ymid = ymid + atom.getY() * m;
2714 zmid = zmid + atom.getZ() * m;
2715 }
2716 if (totalMass > 0) {
2717 xmid /= totalMass;
2718 ymid /= totalMass;
2719 zmid /= totalMass;
2720 }
2721 int n = activeAtoms.length;
2722 double[] xcm = new double[n];
2723 double[] ycm = new double[n];
2724 double[] zcm = new double[n];
2725 int k = 0;
2726 for (Atom atom : activeAtoms) {
2727 xcm[k] = atom.getX() - xmid;
2728 ycm[k] = atom.getY() - ymid;
2729 zcm[k] = atom.getZ() - zmid;
2730 k++;
2731 }
2732
2733 if (forceEnergy) {
2734 energy(false, false);
2735 }
2736
2737
2738 k = 0;
2739 for (Atom atom : activeAtoms) {
2740 int i = atom.getIndex() - 1;
2741 double[] globalMultipolei = globalMultipole[0][i];
2742 double[] inducedDipolei = inducedDipole[0][i];
2743
2744 var ci = globalMultipolei[t000];
2745 var dix = globalMultipolei[t100];
2746 var diy = globalMultipolei[t010];
2747 var diz = globalMultipolei[t001];
2748 var uix = inducedDipolei[0];
2749 var uiy = inducedDipolei[1];
2750 var uiz = inducedDipolei[2];
2751
2752 netchg += ci;
2753 xdpl += xcm[k] * ci + dix + uix;
2754 ydpl += ycm[k] * ci + diy + uiy;
2755 zdpl += zcm[k] * ci + diz + uiz;
2756 xxqdp += xcm[k] * xcm[k] * ci + 2.0 * xcm[k] * (dix + uix);
2757 xyqdp += xcm[k] * ycm[k] * ci + xcm[k] * (diy + uiy) + ycm[k] * (dix + uix);
2758 xzqdp += xcm[k] * zcm[k] * ci + xcm[k] * (diz + uiz) + zcm[k] * (dix + uix);
2759 yxqdp += ycm[k] * xcm[k] * ci + ycm[k] * (dix + uix) + xcm[k] * (diy + uiy);
2760 yyqdp += ycm[k] * ycm[k] * ci + 2.0 * ycm[k] * (diy + uiy);
2761 yzqdp += ycm[k] * zcm[k] * ci + ycm[k] * (diz + uiz) + zcm[k] * (diy + uiy);
2762 zxqdp += zcm[k] * xcm[k] * ci + zcm[k] * (dix + uix) + xcm[k] * (diz + uiz);
2763 zyqdp += zcm[k] * ycm[k] * ci + zcm[k] * (diy + uiy) + ycm[k] * (diz + uiz);
2764 zzqdp += zcm[k] * zcm[k] * ci + 2.0 * zcm[k] * (diz + uiz);
2765 k++;
2766 }
2767
2768
2769 var qave = (xxqdp + yyqdp + zzqdp) / 3.0;
2770 xxqdp = 1.5 * (xxqdp - qave);
2771 xyqdp = 1.5 * xyqdp;
2772 xzqdp = 1.5 * xzqdp;
2773 yxqdp = 1.5 * yxqdp;
2774 yyqdp = 1.5 * (yyqdp - qave);
2775 yzqdp = 1.5 * yzqdp;
2776 zxqdp = 1.5 * zxqdp;
2777 zyqdp = 1.5 * zyqdp;
2778 zzqdp = 1.5 * (zzqdp - qave);
2779
2780
2781 for (Atom atom : activeAtoms) {
2782 int i = atom.getIndex() - 1;
2783 double[] globalMultipolei = globalMultipole[0][i];
2784 var qixx = globalMultipolei[t200];
2785 var qiyy = globalMultipolei[t020];
2786 var qizz = globalMultipolei[t002];
2787 var qixy = globalMultipolei[t110];
2788 var qixz = globalMultipolei[t101];
2789 var qiyz = globalMultipolei[t011];
2790 xxqdp += qixx;
2791 xyqdp += qixy;
2792 xzqdp += qixz;
2793 yxqdp += qixy;
2794 yyqdp += qiyy;
2795 yzqdp += qiyz;
2796 zxqdp += qixz;
2797 zyqdp += qiyz;
2798 zzqdp += qizz;
2799 }
2800
2801
2802 xdpl = xdpl * ELEC_ANG_TO_DEBYE;
2803 ydpl = ydpl * ELEC_ANG_TO_DEBYE;
2804 zdpl = zdpl * ELEC_ANG_TO_DEBYE;
2805 xxqdp = xxqdp * ELEC_ANG_TO_DEBYE;
2806 xyqdp = xyqdp * ELEC_ANG_TO_DEBYE;
2807 xzqdp = xzqdp * ELEC_ANG_TO_DEBYE;
2808 yxqdp = yxqdp * ELEC_ANG_TO_DEBYE;
2809 yyqdp = yyqdp * ELEC_ANG_TO_DEBYE;
2810 yzqdp = yzqdp * ELEC_ANG_TO_DEBYE;
2811 zxqdp = zxqdp * ELEC_ANG_TO_DEBYE;
2812 zyqdp = zyqdp * ELEC_ANG_TO_DEBYE;
2813 zzqdp = zzqdp * ELEC_ANG_TO_DEBYE;
2814
2815
2816 netdpl = sqrt(xdpl * xdpl + ydpl * ydpl + zdpl * zdpl);
2817 double[][] a = new double[3][3];
2818 a[0][0] = xxqdp;
2819 a[0][1] = xyqdp;
2820 a[0][2] = xzqdp;
2821 a[1][0] = yxqdp;
2822 a[1][1] = yyqdp;
2823 a[1][2] = yzqdp;
2824 a[2][0] = zxqdp;
2825 a[2][1] = zyqdp;
2826 a[2][2] = zzqdp;
2827 EigenDecomposition e = new EigenDecomposition(new Array2DRowRealMatrix(a));
2828
2829 var netqdp = e.getRealEigenvalues();
2830
2831 logger.info("\n Electric Moments\n");
2832 logger.info(format(" Total Electric Charge: %13.5f Electrons\n", netchg));
2833 logger.info(format(" Dipole Moment Magnitude: %13.5f Debye\n", netdpl));
2834 logger.info(format(" Dipole X,Y,Z-Components: %13.5f %13.5f %13.5f\n", xdpl, ydpl, zdpl));
2835 logger.info(format(" Quadrupole Moment Tensor: %13.5f %13.5f %13.5f", xxqdp, xyqdp, xzqdp));
2836 logger.info(format(" (Buckinghams) %13.5f %13.5f %13.5f", yxqdp, yyqdp, yzqdp));
2837 logger.info(format(" %13.5f %13.5f %13.5f\n", zxqdp, zyqdp, zzqdp));
2838 logger.info(
2839 format(
2840 " Principal Axes Quadrupole %13.5f %13.5f %13.5f\n", netqdp[2], netqdp[1], netqdp[0]));
2841 }
2842
2843
2844
2845
2846
2847
2848
2849
2850
2851
2852 private void log(int i, int k, double r, double eij) {
2853 logger.info(
2854 format(
2855 "%s %6d-%s %6d-%s %10.4f %10.4f",
2856 "ELEC",
2857 atoms[i].getIndex(),
2858 atoms[i].getAtomType().name,
2859 atoms[k].getIndex(),
2860 atoms[k].getAtomType().name,
2861 r,
2862 eij));
2863 }
2864
2865
2866
2867
2868
2869 private boolean esvTerm = false;
2870
2871
2872
2873
2874 private ExtendedSystem extendedSystem = null;
2875
2876
2877
2878
2879
2880
2881 public void attachExtendedSystem(ExtendedSystem system) {
2882
2883 esvTerm = true;
2884 extendedSystem = system;
2885
2886 setAtoms(extendedSystem.getExtendedAtoms(), extendedSystem.getExtendedMolecule());
2887
2888 if (dMultipoledTirationESV == null || dMultipoledTirationESV.length != nSymm
2889 || dMultipoledTirationESV[0].length != nAtoms) {
2890 dMultipoledTirationESV = new double[nSymm][nAtoms][10];
2891 dMultipoledTautomerESV = new double[nSymm][nAtoms][10];
2892 }
2893
2894 if (dPolardTitrationESV == null || dPolardTitrationESV.length != nAtoms) {
2895 dPolardTitrationESV = new double[nAtoms];
2896 dPolardTautomerESV = new double[nAtoms];
2897 }
2898
2899 logger.info(format(" Attached extended system (%d variables) to PME.\n", extendedSystem.getNumberOfVariables()));
2900 }
2901
2902
2903
2904
2905 double[] perAtomTitrationESV = null;
2906
2907
2908
2909
2910
2911 double[][][] dMultipoledTirationESV = null;
2912 double[][][] dMultipoledTautomerESV = null;
2913
2914
2915
2916
2917
2918 double[] dPolardTitrationESV = null;
2919 double[] dPolardTautomerESV = null;
2920
2921
2922
2923 LambdaFactors[] lambdaFactors = null;
2924
2925
2926
2927
2928
2929
2930
2931
2932
2933
2934
2935
2936 public class LambdaFactors {
2937
2938
2939
2940
2941 protected double lambdaProd = 1.0;
2942
2943
2944
2945 protected double lfAlpha = 0.0;
2946
2947
2948
2949 protected double dlfAlpha = 0.0;
2950
2951
2952
2953 protected double d2lfAlpha = 0.0;
2954
2955
2956
2957 protected double lfPowPerm = 1.0;
2958
2959
2960
2961 protected double dlfPowPerm = 0.0;
2962
2963
2964
2965 protected double d2lfPowPerm = 0.0;
2966
2967
2968
2969 protected double lfPowPol = 1.0;
2970
2971
2972
2973 protected double dlfPowPol = 0.0;
2974
2975
2976
2977 protected double d2lfPowPol = 0.0;
2978
2979
2980
2981 protected double dLpdL = 1.0;
2982
2983
2984
2985 protected double dLpdLi = 1.0;
2986
2987
2988
2989 protected double dLpdLk = 1.0;
2990
2991
2992
2993 protected int[] ik = new int[2];
2994
2995 public void print() {
2996 StringBuilder sb = new StringBuilder();
2997 if (this instanceof LambdaFactorsESV) {
2998 sb.append(format(" (QI-ESV) i,k:%d,%d", ik[0], ik[1]));
2999 } else {
3000 sb.append(" (QI-OST)");
3001 }
3002 sb.append(format(
3003 " lambda:%.2f lAlpha:%.2f,%.2f,%.2f lPowPerm:%.2f,%.2f,%.2f lPowPol:%.2f,%.2f,%.2f",
3004 lambdaProd,
3005 lfAlpha,
3006 dlfAlpha,
3007 d2lfAlpha,
3008 lfPowPerm,
3009 dlfPowPerm,
3010 d2lfPowPerm,
3011 lfPowPol,
3012 dlfPowPol,
3013 d2lfPowPol));
3014
3015 sb.append(format(
3016 "\n permExp:%.2f permAlpha:%.2f permWindow:%.2f,%.2f polExp:%.2f polWindow:%.2f,%.2f",
3017 alchemicalParameters.permLambdaExponent,
3018 alchemicalParameters.permLambdaAlpha,
3019 alchemicalParameters.permLambdaStart,
3020 alchemicalParameters.permLambdaEnd,
3021 alchemicalParameters.polLambdaExponent,
3022 alchemicalParameters.polLambdaStart,
3023 alchemicalParameters.polLambdaEnd));
3024 logger.info(sb.toString());
3025 }
3026
3027
3028
3029
3030
3031
3032
3033
3034 public void setFactors(int i, int k, LambdaMode mode) {
3035
3036 }
3037
3038
3039
3040
3041 public void setFactors() {
3042
3043 }
3044 }
3045
3046 public class LambdaFactorsOST extends LambdaFactors {
3047
3048 @Override
3049 public void setFactors() {
3050 lambdaProd = lambda;
3051 lfAlpha = alchemicalParameters.lAlpha;
3052 dlfAlpha = alchemicalParameters.dlAlpha;
3053 d2lfAlpha = alchemicalParameters.d2lAlpha;
3054 lfPowPerm = alchemicalParameters.permanentScale;
3055 dlfPowPerm =
3056 alchemicalParameters.dlPowPerm * alchemicalParameters.dEdLSign;
3057 d2lfPowPerm = alchemicalParameters.d2lPowPerm
3058 * alchemicalParameters.dEdLSign;
3059 lfPowPol = alchemicalParameters.polarizationScale;
3060 dlfPowPol =
3061 alchemicalParameters.dlPowPol * alchemicalParameters.dEdLSign;
3062 d2lfPowPol =
3063 alchemicalParameters.d2lPowPol * alchemicalParameters.dEdLSign;
3064 }
3065 }
3066
3067 public class LambdaFactorsESV extends LambdaFactors {
3068
3069 @Override
3070 public void setFactors(int i, int k, LambdaMode mode) {
3071 logger.info(format("Invoked Qi setFactors() method with i,k=%d,%d", i, k));
3072 ik[0] = i;
3073 ik[1] = k;
3074 final double L = lambda;
3075 lambdaProd = L * perAtomTitrationESV[i] * perAtomTitrationESV[k];
3076 final double esvli = perAtomTitrationESV[i];
3077 final double esvlk = perAtomTitrationESV[k];
3078 dLpdL = esvli * esvlk;
3079 dLpdLi = L * esvlk;
3080 dLpdLk = L * esvli;
3081
3082
3083 double permLambdaExponent = alchemicalParameters.permLambdaExponent;
3084 double permLambdaStart = alchemicalParameters.permLambdaStart;
3085 double permLambdaEnd = alchemicalParameters.permLambdaEnd;
3086
3087 double permWindow = 1.0 / (permLambdaEnd - permLambdaStart);
3088 double permLambda = (lambdaProd - permLambdaStart) * permWindow;
3089 lfPowPerm = pow(permLambda, alchemicalParameters.permLambdaExponent);
3090 dlfPowPerm = (permLambdaExponent < 1)
3091 ? 0.0 : permLambdaExponent * pow(permLambda, permLambdaExponent - 1) * permWindow;
3092 d2lfPowPerm = (permLambdaExponent < 2)
3093 ? 0.0 : permLambdaExponent
3094 * (permLambdaExponent - 1)
3095 * pow(permLambda, permLambdaExponent - 2)
3096 * permWindow
3097 * permWindow;
3098
3099
3100 double polLambdaExponent = alchemicalParameters.polLambdaExponent;
3101 double polLambdaStart = alchemicalParameters.polLambdaStart;
3102 double polLambdaEnd = alchemicalParameters.polLambdaEnd;
3103
3104 double polWindow = 1.0 / (polLambdaEnd - polLambdaStart);
3105 double polLambda = (lambdaProd - polLambdaStart) * polWindow;
3106 lfPowPol = pow(polLambda, polLambdaExponent);
3107 dlfPowPol = (polLambdaExponent < 1)
3108 ? 0.0 : polLambdaExponent * pow(polLambda, polLambdaExponent - 1) * polWindow;
3109 d2lfPowPol = (polLambdaExponent < 2)
3110 ? 0.0 : polLambdaExponent
3111 * (polLambdaExponent - 1)
3112 * pow(polLambda, polLambdaExponent - 2)
3113 * polWindow
3114 * polWindow;
3115
3116
3117 double permLambdaAlpha = alchemicalParameters.permLambdaAlpha;
3118 lfAlpha = permLambdaAlpha * (1.0 - permLambda) * (1.0 - permLambda);
3119 dlfAlpha = permLambdaAlpha * (1.0 - permLambda) * permWindow;
3120 d2lfAlpha = -permLambdaAlpha * permWindow * permWindow;
3121
3122
3123
3124
3125
3126
3127
3128
3129
3130
3131
3132 switch (mode) {
3133 case CONDENSED -> {
3134 lfPowPerm = 1.0 - lfPowPerm;
3135 dlfPowPerm = -dlfPowPerm;
3136 d2lfPowPerm = -d2lfPowPerm;
3137 if (polarization == Polarization.NONE || lambda < polLambdaStart) {
3138 lfPowPol = 0.0;
3139 dlfPowPol = 0.0;
3140 d2lfPowPol = 0.0;
3141 } else if (lambda <= polLambdaEnd) {
3142
3143 } else {
3144 lfPowPol = 1.0;
3145 dlfPowPol = 0.0;
3146 d2lfPowPol = 0.0;
3147 }
3148 }
3149 case CONDENSED_NO_LIGAND -> {
3150
3151 lfPowPerm = 1.0 - lfPowPerm;
3152 dlfPowPerm = -dlfPowPerm;
3153 d2lfPowPerm = -d2lfPowPerm;
3154 if (polarization == Polarization.NONE) {
3155 lfPowPol = 0.0;
3156 dlfPowPol = 0.0;
3157 d2lfPowPol = 0.0;
3158 } else if (lambda <= polLambdaEnd) {
3159 lfPowPol = 1.0 - lfPowPol;
3160 dlfPowPol = -dlfPowPol;
3161 d2lfPowPol = -d2lfPowPol;
3162 } else {
3163 lfPowPol = 0.0;
3164 dlfPowPol = 0.0;
3165 d2lfPowPol = 0.0;
3166 }
3167 }
3168 case VAPOR -> {
3169
3170 lfPowPerm = 1.0 - lfPowPerm;
3171 dlfPowPerm = -dlfPowPerm;
3172 d2lfPowPerm = -d2lfPowPerm;
3173 lfAlpha = 0.0;
3174 dlfAlpha = 0.0;
3175 d2lfAlpha = 0.0;
3176 if (polarization == Polarization.NONE || lambdaProd > polLambdaEnd) {
3177 lfPowPol = 0.0;
3178 dlfPowPol = 0.0;
3179 d2lfPowPol = 0.0;
3180 } else if (lambdaProd <= polLambdaEnd) {
3181 lfPowPol = 1.0 - lfPowPol;
3182 dlfPowPol = -dlfPowPol;
3183 d2lfPowPol = -d2lfPowPol;
3184 }
3185 }
3186 default -> {
3187 }
3188 }
3189 }
3190 }
3191
3192
3193
3194
3195
3196 public final LambdaFactors LambdaDefaults = new LambdaFactors();
3197 }