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