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