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.ParallelTeam;
41 import ffx.crystal.Crystal;
42 import ffx.numerics.atomic.AtomicDoubleArray;
43 import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
44 import ffx.numerics.atomic.AtomicDoubleArray3D;
45 import ffx.potential.bonded.Atom;
46 import ffx.potential.bonded.LambdaInterface;
47 import ffx.potential.nonbonded.implicit.BornGradRegion;
48 import ffx.potential.nonbonded.implicit.BornRadiiRegion;
49 import ffx.potential.nonbonded.implicit.BornTanhRescaling;
50 import ffx.potential.nonbonded.implicit.ChandlerCavitation;
51 import ffx.potential.nonbonded.implicit.ConnollyRegion;
52 import ffx.potential.nonbonded.implicit.DispersionRegion;
53 import ffx.potential.nonbonded.implicit.GKEnergyRegion;
54 import ffx.potential.nonbonded.implicit.GaussVol;
55 import ffx.potential.nonbonded.implicit.InducedGKFieldRegion;
56 import ffx.potential.nonbonded.implicit.InitializationRegion;
57 import ffx.potential.nonbonded.implicit.PermanentGKFieldRegion;
58 import ffx.potential.nonbonded.implicit.SurfaceAreaRegion;
59 import ffx.potential.nonbonded.pme.Polarization;
60 import ffx.potential.parameters.AtomType;
61 import ffx.potential.parameters.ForceField;
62 import ffx.potential.parameters.SoluteType;
63 import ffx.potential.parameters.SoluteType.SOLUTE_RADII_TYPE;
64 import ffx.utilities.Constants;
65 import ffx.utilities.FFXProperty;
66
67 import java.util.Arrays;
68 import java.util.HashMap;
69 import java.util.logging.Level;
70 import java.util.logging.Logger;
71
72 import static ffx.potential.nonbonded.implicit.DispersionRegion.DEFAULT_DISPERSION_OFFSET;
73 import static ffx.potential.parameters.ForceField.toEnumForm;
74 import static ffx.potential.parameters.SoluteType.setSoluteRadii;
75 import static ffx.utilities.Constants.dWater;
76 import static ffx.utilities.PropertyGroup.ImplicitSolvent;
77 import static java.lang.String.format;
78 import static java.util.Arrays.fill;
79 import static org.apache.commons.math3.util.FastMath.max;
80
81
82
83
84
85
86
87
88
89
90
91
92 public class GeneralizedKirkwood implements LambdaInterface {
93
94
95
96
97 public static final double DEFAULT_DIELECTRIC_OFFSET = 0.09;
98
99 private static final Logger logger = Logger.getLogger(GeneralizedKirkwood.class.getName());
100
101
102
103 private static final double DEFAULT_SOLUTE_SCALE = 1.0;
104
105
106
107
108
109
110
111 private final double dOffset = DEFAULT_DIELECTRIC_OFFSET;
112
113
114
115 private final ForceField forceField;
116
117
118
119 private final Polarization polarization;
120
121
122
123
124 private final ParticleMeshEwald particleMeshEwald;
125
126
127
128 private final ParallelTeam parallelTeam;
129
130
131
132 private final InitializationRegion initializationRegion;
133
134
135
136 private final BornRadiiRegion bornRadiiRegion;
137
138
139
140 private final PermanentGKFieldRegion permanentGKFieldRegion;
141
142
143
144 private final InducedGKFieldRegion inducedGKFieldRegion;
145
146
147
148 private final GKEnergyRegion gkEnergyRegion;
149
150
151
152 private final BornGradRegion bornGradRegion;
153
154
155
156 private final DispersionRegion dispersionRegion;
157
158
159
160 private final SurfaceAreaRegion surfaceAreaRegion;
161
162
163
164 private final ChandlerCavitation chandlerCavitation;
165
166
167
168
169 private final HashMap<Integer, Double> radiiOverride = new HashMap<>();
170
171
172
173 public double electric;
174
175
176
177 private final double bondiScale;
178
179
180
181
182 private Atom[] atoms;
183
184
185
186 private int nAtoms;
187
188
189
190 private double[] x, y, z;
191
192
193
194 private double[] baseRadius;
195
196
197
198 private double[] descreenRadius;
199
200
201
202
203
204
205
206
207 private double[] overlapScale;
208
209
210
211
212 private double[] neckScale;
213
214
215
216 private final boolean perfectHCTScale;
217
218
219
220
221 @FFXProperty(name = "solvent-dielectric", propertyGroup = ImplicitSolvent, defaultValue = "78.3", description = """
222 The dielectric constant used for the solvent in generalized Kirkwood calculations.
223 The default of 78.3 corresponds to water.
224 """)
225 private final double solventDielectric;
226
227
228
229 @FFXProperty(name = "solute-dielectric", propertyGroup = ImplicitSolvent, defaultValue = "1.0", description = """
230 The dielectric constant used for the solute(s) in generalized Kirkwood calculations.
231 The default of 1.0 is consistent with all solute dielectric response arising from either
232 polarization via induced dipoles and/or permanent dipole realignment.
233 """)
234 private final double soluteDielectric;
235
236
237
238
239
240 private static final double DEFAULT_HCT_SCALE = 0.72;
241
242
243
244
245 @FFXProperty(name = "hct-scale", propertyGroup = ImplicitSolvent, defaultValue = "0.72", description =
246 "The default overlap scale factor for Hawkins-Cramer-Truhlar pairwise descreening.")
247 private final double hctScale;
248
249
250
251
252 @FFXProperty(name = "element-hct-scale", clazz = Boolean.class, propertyGroup = ImplicitSolvent, defaultValue = "true",
253 description =
254 "Flag to turn on element specific overlap scale factors for Hawkins-Cramer-Truhlar pairwise descreening.")
255 private final boolean elementHCTScale;
256
257
258
259
260 private final HashMap<Integer, Double> elementHCTScaleFactors;
261
262
263
264
265 @FFXProperty(name = "descreen-vdw", propertyGroup = ImplicitSolvent, defaultValue = "true", description =
266 "If true, the descreening size of each atom is based on its force field van der Waals radius.")
267 private final boolean descreenVDW;
268
269
270
271
272 @FFXProperty(name = "descreen-hydrogen", propertyGroup = ImplicitSolvent, defaultValue = "false", description =
273 "If true, hydrogen atoms are contribute to the pairwise descreening integrals.")
274 private final boolean descreenHydrogen;
275
276 private static final double DEFAULT_DESCREEN_OFFSET = 0.3;
277
278
279
280 @FFXProperty(name = "descreen-offset", propertyGroup = ImplicitSolvent, defaultValue = "0.3", description =
281 "Offset applied to the pairwise descreening integral to improve stability at small separation.")
282 private final double descreenOffset;
283
284
285
286
287 @FFXProperty(name = "neck-correction", propertyGroup = ImplicitSolvent, defaultValue = "true", description =
288 "Apply a neck correction during descreening.")
289 private final boolean neckCorrection;
290
291
292
293
294 private static final double DEFAULT_NECK_SCALE = 0.1350;
295
296
297
298 @FFXProperty(name = "neck-scale", propertyGroup = ImplicitSolvent, defaultValue = "0.1350", description =
299 "The overlap scale factor to use during the descreening neck correction.")
300 private double sneck;
301
302
303
304
305
306 @FFXProperty(name = "chemically-aware-neck-scale", propertyGroup = ImplicitSolvent, defaultValue = "true",
307 description = """
308 If the neck descreening correction is being used, apply a smaller overlap scale
309 factors as the number of bonded heavy atoms increases.
310 """)
311 private final boolean chemicallyAwareSneck;
312
313
314
315
316
317 @FFXProperty(name = "tanh-correction", propertyGroup = ImplicitSolvent, defaultValue = "true", description = """
318 If the neck descreening correction is being used, apply a smaller overlap scale
319 factors as the number of bonded heavy atoms increases.
320 """)
321 private final boolean tanhCorrection;
322
323
324
325
326 public static final double DEFAULT_TANH_BETA0 = 0.9563;
327
328
329
330 public static final double DEFAULT_TANH_BETA1 = 0.2578;
331
332
333
334 public static final double DEFAULT_TANH_BETA2 = 0.0810;
335
336
337
338
339 @FFXProperty(name = "tanh-beta0", propertyGroup = ImplicitSolvent, defaultValue = "0.9563", description =
340 "The coefficient beta0 for tanh rescaling of descreening integrals.")
341 private double beta0;
342
343
344
345
346 @FFXProperty(name = "tanh-beta1", propertyGroup = ImplicitSolvent, defaultValue = "0.2578", description =
347 "The coefficient beta1 for tanh rescaling of descreening integrals.")
348 private double beta1;
349
350
351
352
353 @FFXProperty(name = "tanh-beta2", propertyGroup = ImplicitSolvent, defaultValue = "0.0810", description =
354 "The coefficient beta2 for tanh rescaling of descreening integrals.")
355 private double beta2;
356
357
358
359
360 public static final double DEFAULT_GKC = 2.455;
361
362
363
364 @FFXProperty(name = "gkc", propertyGroup = ImplicitSolvent, defaultValue = "2.455", description =
365 "The Generalized Kirkwood cross-term parameter.")
366 public final double gkc;
367
368 private static final NonPolarModel DEFAULT_NONPOLAR_MODEL = NonPolarModel.GAUSS_DISP;
369
370
371
372
373 @FFXProperty(name = "nonpolar-model", clazz = String.class, propertyGroup = ImplicitSolvent,
374 defaultValue = "gauss-disp", description = """
375 [CAV / CAV-DISP / DISP / GAUSS-DISP / SEV-DISP / NONE ]
376 The non-polar contribution to the implicit solvent.
377 """)
378 private final NonPolarModel nonPolarModel;
379
380
381
382
383
384 private static final double DEFAULT_CAVONLY_SURFACE_TENSION = 0.0049;
385
386
387
388
389 public final double probe;
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405 public static final double DEFAULT_SOLVENT_PRESSURE = 0.0334;
406
407
408
409
410
411
412
413 public static final double DEFAULT_CAVDISP_SURFACE_TENSION = 0.080;
414
415
416
417
418
419
420
421
422
423
424 public static final double DEFAULT_CROSSOVER =
425 3.0 * DEFAULT_CAVDISP_SURFACE_TENSION / DEFAULT_SOLVENT_PRESSURE;
426
427
428
429
430 @FFXProperty(name = "surface-tension", propertyGroup = ImplicitSolvent, defaultValue = "0.080", description =
431 "The cavitation surface tension coefficient (kcal/mol/A^2).")
432 private final double surfaceTension;
433
434
435
436
437 @FFXProperty(name = "solvent-pressure", propertyGroup = ImplicitSolvent, defaultValue = "0.0334", description =
438 "The solvent pressure for nonpolar models with an explicit volume term (kcal/mol/A^3).")
439 private final double solventPressue;
440
441
442
443
444 @FFXProperty(name = "gk-radius", clazz = String.class, propertyGroup = ImplicitSolvent, defaultValue = "solute",
445 description = """
446 [SOLUTE / VDW / CONSENSUS]
447 The base atomic radii to use for generalized Kirkwood calculations. The default is to use solute radii,
448 which were fit to experimental solvation free energy differences. Alternatively, force field
449 specific van der Waals radii (vdw) or consensus Bondi radii (consensus) can be chosen.
450 """)
451 private SOLUTE_RADII_TYPE soluteRadiiType;
452
453
454
455
456 private double[] born;
457
458
459
460 private boolean[] use = null;
461
462
463
464 private Crystal crystal;
465
466
467
468 private double[][][] sXYZ;
469
470
471
472 private double[][][] globalMultipole;
473
474
475
476 private double[][][] inducedDipole;
477
478
479
480 private double[][][] inducedDipoleCR;
481
482
483
484 private AtomicDoubleArrayImpl atomicDoubleArrayImpl;
485
486
487
488 private AtomicDoubleArray3D grad;
489
490
491
492 private AtomicDoubleArray3D torque;
493
494
495
496 private AtomicDoubleArray bornRadiiChainRule;
497
498
499
500 private final AtomicDoubleArray3D fieldGK;
501
502
503
504 private final AtomicDoubleArray3D fieldGKCR;
505
506
507
508 private int[][][] neighborLists;
509
510
511
512
513
514 private int maxNumAtoms;
515
516
517
518 private double cutoff;
519
520
521
522 private double cut2;
523
524
525
526 private boolean lambdaTerm;
527
528
529
530 private double lambda = 1.0;
531
532
533
534
535 private double lPow = 1.0;
536
537
538
539 private double dlPow = 0.0;
540
541
542
543 private double dl2Pow = 0.0;
544
545
546
547 private double solvationEnergy = 0.0;
548
549
550
551 private double gkEnergy = 0.0;
552
553
554
555 private double gkPermanentEnergy = 0.0;
556
557
558
559 private double gkPolarizationEnergy = 0.0;
560
561
562
563 private double dispersionEnergy = 0.0;
564
565
566
567 private double cavitationEnergy = 0.0;
568
569
570
571 private long gkTime = 0;
572
573
574
575 private long dispersionTime = 0;
576
577
578
579 private long cavitationTime = 0;
580
581
582
583 private final boolean nativeEnvironmentApproximation;
584
585
586
587 private final boolean perfectRadii;
588
589
590
591
592
593
594
595
596
597
598 public GeneralizedKirkwood(ForceField forceField, Atom[] atoms, ParticleMeshEwald particleMeshEwald,
599 Crystal crystal, ParallelTeam parallelTeam, double gkCutoff) {
600 this.forceField = forceField;
601 this.atoms = atoms;
602 this.particleMeshEwald = particleMeshEwald;
603 this.crystal = crystal;
604 this.parallelTeam = parallelTeam;
605 this.electric = electric;
606 this.cutoff = gkCutoff;
607 nAtoms = atoms.length;
608 maxNumAtoms = nAtoms;
609 polarization = particleMeshEwald.polarization;
610
611
612 electric = forceField.getDouble("ELECTRIC", Constants.DEFAULT_ELECTRIC);
613
614 solventDielectric = forceField.getDouble("SOLVENT_DIELECTRIC", dWater);
615
616 soluteDielectric = forceField.getDouble("SOLUTE_DIELECTRIC", 1.0);
617
618
619 String value = forceField.getString("ARRAY_REDUCTION", "MULTI");
620 try {
621 atomicDoubleArrayImpl = AtomicDoubleArrayImpl.valueOf(toEnumForm(value));
622 } catch (Exception e) {
623 atomicDoubleArrayImpl = AtomicDoubleArrayImpl.MULTI;
624 logger.info(
625 format(" Unrecognized ARRAY-REDUCTION %s; defaulting to %s", value,
626 atomicDoubleArrayImpl));
627 }
628
629 String gkRadius = forceField.getString("GK_RADIUS", "SOLUTE");
630 try {
631 soluteRadiiType = SOLUTE_RADII_TYPE.valueOf(gkRadius.trim().toUpperCase());
632 } catch (Exception e) {
633 soluteRadiiType = SOLUTE_RADII_TYPE.SOLUTE;
634 }
635
636
637 if (soluteRadiiType != SOLUTE_RADII_TYPE.SOLUTE) {
638 bondiScale = forceField.getDouble("SOLUTE_SCALE", DEFAULT_SOLUTE_SCALE);
639 } else {
640
641 bondiScale = forceField.getDouble("SOLUTE_SCALE", 1.0);
642 }
643
644 hctScale = forceField.getDouble("HCT_SCALE", DEFAULT_HCT_SCALE);
645 elementHCTScale = forceField.getBoolean("ELEMENT_HCT_SCALE", true);
646 descreenVDW = forceField.getBoolean("DESCREEN_VDW", true);
647 descreenHydrogen = forceField.getBoolean("DESCREEN_HYDROGEN", false);
648 if (descreenVDW && !descreenHydrogen) {
649 perfectHCTScale = forceField.getBoolean("PERFECT_HCT_SCALE", false);
650 } else {
651 perfectHCTScale = false;
652 }
653 descreenOffset = forceField.getDouble("DESCREEN_OFFSET", DEFAULT_DESCREEN_OFFSET);
654
655 neckCorrection = forceField.getBoolean("NECK_CORRECTION", true);
656 double sn = forceField.getDouble("SNECK", DEFAULT_NECK_SCALE);
657 sneck = forceField.getDouble("NECK_SCALE", sn);
658 chemicallyAwareSneck = forceField.getBoolean("CHEMICALLY_AWARE_SNECK", true);
659 tanhCorrection = forceField.getBoolean("TANH_CORRECTION", true);
660 double b0 = forceField.getDouble("BETA0", DEFAULT_TANH_BETA0);
661 double b1 = forceField.getDouble("BETA1", DEFAULT_TANH_BETA1);
662 double b2 = forceField.getDouble("BETA2", DEFAULT_TANH_BETA2);
663 beta0 = forceField.getDouble("TANH_BETA0", b0);
664 beta1 = forceField.getDouble("TANH_BETA1", b1);
665 beta2 = forceField.getDouble("TANH_BETA2", b2);
666
667 BornTanhRescaling.setBeta0(beta0);
668 BornTanhRescaling.setBeta1(beta1);
669 BornTanhRescaling.setBeta2(beta2);
670
671
672 HashMap<Integer, Double> DEFAULT_HCT_ELEMENTS = new HashMap<>();
673
674 DEFAULT_HCT_ELEMENTS.put(1, 0.7200);
675 DEFAULT_HCT_ELEMENTS.put(6, 0.6950);
676 DEFAULT_HCT_ELEMENTS.put(7, 0.7673);
677 DEFAULT_HCT_ELEMENTS.put(8, 0.7965);
678 DEFAULT_HCT_ELEMENTS.put(15, 0.6117);
679 DEFAULT_HCT_ELEMENTS.put(16, 0.7204);
680
681
682 elementHCTScaleFactors = new HashMap<>();
683 elementHCTScaleFactors.put(1, DEFAULT_HCT_ELEMENTS.get(1));
684 elementHCTScaleFactors.put(6, DEFAULT_HCT_ELEMENTS.get(6));
685 elementHCTScaleFactors.put(7, DEFAULT_HCT_ELEMENTS.get(7));
686 elementHCTScaleFactors.put(8, DEFAULT_HCT_ELEMENTS.get(8));
687 elementHCTScaleFactors.put(15, DEFAULT_HCT_ELEMENTS.get(15));
688 elementHCTScaleFactors.put(16, DEFAULT_HCT_ELEMENTS.get(16));
689
690
691
692
693
694 String[] tmpElemScaleFactors = forceField.getProperties().getStringArray("hct-element");
695 for (String elemScale : tmpElemScaleFactors) {
696 String[] singleScale = elemScale.trim().split(" +");
697 if (elementHCTScaleFactors.containsKey(Integer.parseInt(singleScale[0]))) {
698 elementHCTScaleFactors.replace(Integer.parseInt(singleScale[0]), Double.parseDouble(singleScale[1]));
699 } else {
700 elementHCTScaleFactors.put(Integer.parseInt(singleScale[0]), Double.parseDouble(singleScale[1]));
701 }
702 }
703
704 perfectRadii = forceField.getBoolean("PERFECT_RADII", false);
705
706
707 String radiiProp = forceField.getString("GK_RADIIOVERRIDE", null);
708 if (radiiProp != null) {
709 String[] tokens = radiiProp.split("A");
710 for (String token : tokens) {
711 if (!token.contains("r")) {
712 logger.severe("Invalid radius override.");
713 }
714 int separator = token.indexOf("r");
715 int type = Integer.parseInt(token.substring(0, separator));
716 double factor = Double.parseDouble(token.substring(separator + 1));
717 logger.info(format(" (GK) Scaling AtomType %d with bondi factor %.2f", type, factor));
718 radiiOverride.put(type, factor);
719 }
720 }
721
722 NonPolarModel npModel;
723 try {
724 String cavModel = forceField.getString("CAVMODEL", DEFAULT_NONPOLAR_MODEL.name());
725 cavModel = forceField.getString("NONPOLAR_MODEL", cavModel);
726 npModel = getNonPolarModel(cavModel.toUpperCase());
727 } catch (Exception ex) {
728 npModel = NonPolarModel.NONE;
729 logger.warning(format(" Error parsing non-polar model (set to NONE) %s", ex));
730 }
731 nonPolarModel = npModel;
732
733 int threadCount = parallelTeam.getThreadCount();
734 fieldGK = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
735 fieldGKCR = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
736
737 nativeEnvironmentApproximation =
738 forceField.getBoolean("NATIVE_ENVIRONMENT_APPROXIMATION", false);
739 probe = forceField.getDouble("PROBE_RADIUS", 1.4);
740 cut2 = cutoff * cutoff;
741 lambdaTerm = forceField.getBoolean("GK_LAMBDATERM", forceField.getBoolean("LAMBDATERM", false));
742
743
744
745 double polLambdaExp = forceField.getDouble("POLARIZATION_LAMBDA_EXPONENT", 3.0);
746 if (polLambdaExp == 0.0) {
747 lambdaTerm = false;
748 logger.info(" GK lambda term set to false.");
749 }
750
751
752 if (!lambdaTerm
753 && particleMeshEwald.getPolarizationType() != Polarization.NONE) {
754 if (forceField.getBoolean("ELEC_LAMBDATERM",
755 forceField.getBoolean("LAMBDATERM", false))) {
756 logger.info(" If PME includes polarization and is a function of lambda, GK must also.");
757 lambdaTerm = true;
758 }
759 }
760
761 initAtomArrays();
762
763 double tensionDefault;
764 switch (nonPolarModel) {
765 case CAV:
766 tensionDefault = DEFAULT_CAVONLY_SURFACE_TENSION;
767 surfaceAreaRegion = new SurfaceAreaRegion(atoms, x, y, z, use,
768 neighborLists, grad, threadCount, probe, tensionDefault);
769 dispersionRegion = null;
770 chandlerCavitation = null;
771 break;
772 case CAV_DISP:
773 tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
774 surfaceAreaRegion = new SurfaceAreaRegion(atoms, x, y, z, use,
775 neighborLists, grad, threadCount, probe, tensionDefault);
776 dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
777 chandlerCavitation = null;
778 break;
779 case DISP:
780 tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
781 dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
782 surfaceAreaRegion = null;
783 chandlerCavitation = null;
784 break;
785 case SEV_DISP:
786 tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
787 double[] radii = new double[nAtoms];
788 int index = 0;
789 for (Atom atom : atoms) {
790 radii[index] = atom.getVDWType().radius / 2.0;
791 index++;
792 }
793 ConnollyRegion connollyRegion = new ConnollyRegion(atoms, radii, threadCount);
794
795
796 double wiggle = forceField.getDouble("WIGGLE", ConnollyRegion.DEFAULT_WIGGLE);
797 connollyRegion.setWiggle(wiggle);
798 chandlerCavitation = new ChandlerCavitation(atoms, connollyRegion, forceField);
799 dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
800 surfaceAreaRegion = null;
801 break;
802 case GAUSS_DISP:
803 dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
804 surfaceAreaRegion = null;
805 tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
806 GaussVol gaussVol = new GaussVol(atoms, forceField, parallelTeam);
807 chandlerCavitation = new ChandlerCavitation(atoms, gaussVol, forceField);
808 break;
809 case BORN_CAV_DISP:
810 tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
811 surfaceAreaRegion = null;
812 chandlerCavitation = null;
813 dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
814 break;
815 case HYDROPHOBIC_PMF:
816 case BORN_SOLV:
817 case NONE:
818 default:
819 tensionDefault = DEFAULT_CAVONLY_SURFACE_TENSION;
820 surfaceAreaRegion = null;
821 dispersionRegion = null;
822 chandlerCavitation = null;
823 break;
824 }
825
826 gkc = forceField.getDouble("GKC", DEFAULT_GKC);
827
828 surfaceTension = forceField.getDouble("SURFACE_TENSION", tensionDefault);
829 solventPressue = forceField.getDouble("SOLVENT_PRESSURE", DEFAULT_SOLVENT_PRESSURE);
830
831 double crossOver = forceField.getDouble("CROSS_OVER", DEFAULT_CROSSOVER);
832 if (chandlerCavitation != null) {
833
834 chandlerCavitation.setCrossOver(crossOver);
835
836 chandlerCavitation.setSurfaceTension(surfaceTension);
837 chandlerCavitation.setSolventPressure(solventPressue);
838 crossOver = chandlerCavitation.getCrossOver();
839 }
840 if (surfaceAreaRegion != null) {
841 surfaceAreaRegion.setSurfaceTension(surfaceTension);
842 }
843 double dispersionOffset = forceField.getDouble("DISPERSION_OFFSET", DEFAULT_DISPERSION_OFFSET);
844 if (dispersionRegion != null) {
845 dispersionRegion.setDispersionOffest(dispersionOffset);
846 }
847
848 initializationRegion = new InitializationRegion(threadCount);
849 bornRadiiRegion = new BornRadiiRegion(threadCount, nAtoms, forceField, neckCorrection,
850 tanhCorrection, perfectHCTScale);
851 permanentGKFieldRegion = new PermanentGKFieldRegion(threadCount, soluteDielectric, solventDielectric, gkc);
852 inducedGKFieldRegion = new InducedGKFieldRegion(threadCount, soluteDielectric, solventDielectric, gkc);
853 if (!perfectRadii) {
854 bornGradRegion = new BornGradRegion(threadCount, neckCorrection, tanhCorrection, perfectHCTScale);
855 } else {
856
857 bornGradRegion = null;
858 }
859 boolean gkQI = forceField.getBoolean("GK_QI", false);
860 gkEnergyRegion = new GKEnergyRegion(threadCount, polarization, nonPolarModel, surfaceTension,
861 probe, electric, soluteDielectric, solventDielectric, gkc, gkQI);
862
863 if (gkQI) {
864 logger.info(" Continuum Solvation (QI)");
865 } else {
866 logger.info(" Continuum Solvation");
867 }
868 logger.info(format(" Radii: %8s", soluteRadiiType));
869 if (cutoff != Double.POSITIVE_INFINITY) {
870 logger.info(format(" Generalized Kirkwood Cut-Off: %8.4f (A)", cutoff));
871 } else {
872 logger.info(" Generalized Kirkwood Cut-Off: NONE");
873 }
874 logger.info(format(" GKC: %8.4f", gkc));
875 logger.info(format(" Solvent Dielectric: %8.4f", solventDielectric));
876 logger.info(format(" Solute Dielectric: %8.4f", soluteDielectric));
877
878 if (perfectRadii) {
879 logger.info(format(" Use Perfect Born Radii: %8B", perfectRadii));
880 } else {
881 logger.info(format(" Descreen with vdW Radii: %8B", descreenVDW));
882 logger.info(format(" Descreen with Hydrogen Atoms: %8B", descreenHydrogen));
883 logger.info(format(" Descreen Offset: %8.4f", descreenOffset));
884 if (neckCorrection) {
885 logger.info(format(" Use Neck Correction: %8B", neckCorrection));
886 logger.info(format(" Sneck Scale Factor: %8.4f", sneck));
887 logger.info(format(" Chemically Aware Sneck: %8B", chemicallyAwareSneck));
888 }
889 if (tanhCorrection) {
890 logger.info(format(" Use Tanh Correction %8B", tanhCorrection));
891 logger.info(format(" Beta0: %8.4f", beta0));
892 logger.info(format(" Beta1: %8.4f", beta1));
893 logger.info(format(" Beta2: %8.4f", beta2));
894 }
895 if (perfectHCTScale) {
896 logger.info(format(" GaussVol HCT Scale Factors: %8B", perfectHCTScale));
897 }
898 logger.info(format(" General HCT Scale Factor: %8.4f",
899 forceField.getDouble("HCT-SCALE", DEFAULT_HCT_SCALE)));
900 if (elementHCTScale) {
901 logger.info(format(" Element-Specific HCT Scale Factors: %8B", elementHCTScale));
902 Integer[] elementHCTkeyset = elementHCTScaleFactors.keySet().toArray(new Integer[0]);
903 Arrays.sort(elementHCTkeyset);
904 for (Integer key : elementHCTkeyset) {
905 logger.info(format(" HCT-Element # %d: %8.4f", key, elementHCTScaleFactors.get(key)));
906 }
907 }
908 }
909
910 logger.info(
911 format(" Non-Polar Model: %10s",
912 nonPolarModel.toString().replace('_', '-')));
913
914 if (nonPolarModel.equals(NonPolarModel.GAUSS_DISP)) {
915 logger.info(format(" GaussVol Radii Offset: %2.4f",
916 forceField.getDouble("GAUSSVOL_RADII_OFFSET", 0.0)));
917 logger.info(format(" GaussVol Radii Scale: %2.4f",
918 forceField.getDouble("GAUSSVOL_RADII_SCALE", 1.15)));
919 }
920
921 if (dispersionRegion != null) {
922 logger.info(format(" Dispersion Integral Offset: %8.4f (A)",
923 dispersionRegion.getDispersionOffset()));
924 }
925
926 if (surfaceAreaRegion != null) {
927 logger.info(format(" Cavitation Probe Radius: %8.4f (A)", probe));
928 logger.info(
929 format(" Cavitation Surface Tension: %8.4f (Kcal/mol/A^2)", surfaceTension));
930 } else if (chandlerCavitation != null) {
931 logger.info(
932 format(" Cavitation Solvent Pressure: %8.4f (Kcal/mol/A^3)", solventPressue));
933 logger.info(
934 format(" Cavitation Surface Tension: %8.4f (Kcal/mol/A^2)", surfaceTension));
935 logger.info(format(" Cavitation Cross-Over Radius: %8.4f (A)", crossOver));
936 }
937
938
939 if (logger.isLoggable(Level.FINE)) {
940 logger.fine(" Base Radii Descreen Radius Overlap Scale Neck Scale");
941 for (int i = 0; i < nAtoms; i++) {
942 logger.info(
943 format(" %s %8.6f %8.6f %5.3f %5.3f",
944 atoms[i].toString(), baseRadius[i], descreenRadius[i], overlapScale[i],
945 neckScale[i]));
946 }
947 }
948 }
949
950
951
952
953
954
955 public boolean getUsePerfectRadii() {
956 return perfectRadii;
957 }
958
959
960
961
962
963
964
965 public double[] getPerfectRadii() {
966 bornRadiiRegion.init(atoms, crystal, sXYZ, neighborLists, baseRadius, descreenRadius,
967 overlapScale, neckScale, descreenOffset, use, cut2, nativeEnvironmentApproximation, born);
968 return bornRadiiRegion.getPerfectRadii();
969 }
970
971
972
973
974
975
976
977 public static NonPolarModel getNonPolarModel(String nonpolarModel) {
978 try {
979 return NonPolarModel.valueOf(toEnumForm(nonpolarModel));
980 } catch (IllegalArgumentException ex) {
981 logger.warning(" Unrecognized nonpolar model requested; defaulting to NONE.");
982 return NonPolarModel.NONE;
983 }
984 }
985
986
987
988
989 public void computeBornRadii() {
990
991 applySoluteRadii();
992
993
994 if (tanhCorrection) {
995 BornTanhRescaling.setBeta0(beta0);
996 BornTanhRescaling.setBeta1(beta1);
997 BornTanhRescaling.setBeta2(beta2);
998 }
999
1000 try {
1001 bornRadiiRegion.init(
1002 atoms,
1003 crystal,
1004 sXYZ,
1005 neighborLists,
1006 baseRadius,
1007 descreenRadius,
1008 overlapScale,
1009 neckScale,
1010 descreenOffset,
1011 use,
1012 cut2,
1013 nativeEnvironmentApproximation,
1014 born);
1015 parallelTeam.execute(bornRadiiRegion);
1016 } catch (Exception e) {
1017 String message = "Fatal exception computing Born radii.";
1018 logger.log(Level.SEVERE, message, e);
1019 }
1020 }
1021
1022
1023
1024
1025 public void computeInducedGKField() {
1026 try {
1027 fieldGK.reset(parallelTeam);
1028 fieldGKCR.reset(parallelTeam);
1029 inducedGKFieldRegion.init(atoms, inducedDipole, inducedDipoleCR, crystal, sXYZ,
1030 neighborLists, use, cut2, born, fieldGK, fieldGKCR);
1031 parallelTeam.execute(inducedGKFieldRegion);
1032 fieldGK.reduce(parallelTeam);
1033 fieldGKCR.reduce(parallelTeam);
1034 } catch (Exception e) {
1035 String message = "Fatal exception computing induced GK field.";
1036 logger.log(Level.SEVERE, message, e);
1037 }
1038 }
1039
1040
1041
1042
1043 public void computePermanentGKField() {
1044 try {
1045 fieldGK.reset(parallelTeam);
1046 permanentGKFieldRegion.init(
1047 atoms, globalMultipole, crystal, sXYZ, neighborLists, use, cut2, born, fieldGK);
1048 parallelTeam.execute(permanentGKFieldRegion);
1049 fieldGK.reduce(parallelTeam);
1050 } catch (Exception e) {
1051 String message = "Fatal exception computing permanent GK field.";
1052 logger.log(Level.SEVERE, message, e);
1053 }
1054 }
1055
1056
1057
1058
1059
1060
1061
1062
1063 public double[] getBaseRadii() {
1064 return baseRadius;
1065 }
1066
1067
1068
1069
1070
1071
1072
1073
1074 public double[] getDescreenRadii() {
1075 return descreenRadius;
1076 }
1077
1078
1079
1080
1081
1082
1083 public double getCavitationEnergy() {
1084 return cavitationEnergy;
1085 }
1086
1087
1088
1089
1090
1091
1092 public ChandlerCavitation getChandlerCavitation() {
1093 return chandlerCavitation;
1094 }
1095
1096
1097
1098
1099
1100
1101 public double getCutoff() {
1102 return cutoff;
1103 }
1104
1105
1106
1107
1108
1109
1110 public void setCutoff(double cutoff) {
1111 this.cutoff = cutoff;
1112 this.cut2 = cutoff * cutoff;
1113 }
1114
1115
1116
1117
1118
1119
1120 public double getDielecOffset() {
1121 return dOffset;
1122 }
1123
1124
1125
1126
1127
1128
1129 public double getDescreenOffset() {
1130 return descreenOffset;
1131 }
1132
1133
1134
1135
1136
1137
1138 public double getDispersionEnergy() {
1139 return dispersionEnergy;
1140 }
1141
1142 public DispersionRegion getDispersionRegion() {
1143 return dispersionRegion;
1144 }
1145
1146 public AtomicDoubleArray3D getFieldGK() {
1147 return fieldGK;
1148 }
1149
1150 public AtomicDoubleArray3D getFieldGKCR() {
1151 return fieldGKCR;
1152 }
1153
1154 public SurfaceAreaRegion getSurfaceAreaRegion() {
1155 return surfaceAreaRegion;
1156 }
1157
1158
1159
1160
1161
1162
1163 public double getGeneralizedKirkwoordEnergy() {
1164 return gkEnergy;
1165 }
1166
1167
1168
1169
1170
1171
1172 public double getGeneralizedKirkwoordPermanentEnergy() {
1173 return gkPermanentEnergy;
1174 }
1175
1176
1177
1178
1179
1180
1181 public double getGeneralizedKirkwoordPolariztionEnergy() {
1182 return gkPolarizationEnergy;
1183 }
1184
1185 public AtomicDoubleArray3D getGrad() {
1186 return grad;
1187 }
1188
1189
1190
1191
1192
1193
1194 public int getInteractions() {
1195 return gkEnergyRegion.getInteractions();
1196 }
1197
1198
1199
1200
1201 @Override
1202 public double getLambda() {
1203 return lambda;
1204 }
1205
1206
1207
1208
1209
1210
1211 @Override
1212 public void setLambda(double lambda) {
1213 if (lambdaTerm) {
1214 this.lambda = lambda;
1215 } else {
1216
1217 this.lambda = 1.0;
1218 lPow = 1.0;
1219 dlPow = 0.0;
1220 dl2Pow = 0.0;
1221 }
1222 }
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235 public boolean getNativeEnvironmentApproximation() {
1236 return nativeEnvironmentApproximation;
1237 }
1238
1239
1240
1241
1242
1243
1244 public NonPolarModel getNonPolarModel() {
1245 return nonPolarModel;
1246 }
1247
1248
1249
1250
1251
1252
1253 public double[] getOverlapScale() {
1254 return overlapScale;
1255 }
1256
1257
1258
1259
1260
1261
1262 public double[] getNeckScale() {
1263 return neckScale;
1264 }
1265
1266
1267
1268
1269
1270
1271
1272
1273 public boolean getTanhCorrection() {
1274 return tanhCorrection;
1275 }
1276
1277
1278
1279
1280
1281
1282 public double getProbeRadius() {
1283 return probe;
1284 }
1285
1286
1287
1288
1289
1290
1291 public double getSolventPermittivity() {
1292 return solventDielectric;
1293 }
1294
1295
1296
1297
1298
1299
1300 public double getSolutePermittivity() {
1301 return soluteDielectric;
1302 }
1303
1304
1305
1306
1307
1308
1309 public double getSurfaceTension() {
1310 return surfaceTension;
1311 }
1312
1313 public AtomicDoubleArray3D getTorque() {
1314 return torque;
1315 }
1316
1317
1318
1319
1320
1321
1322 @Override
1323 public double getd2EdL2() {
1324 if (lambdaTerm) {
1325 return dl2Pow * solvationEnergy;
1326 } else {
1327 return 0.0;
1328 }
1329 }
1330
1331
1332
1333
1334 @Override
1335 public double getdEdL() {
1336 if (lambdaTerm) {
1337 return dlPow * solvationEnergy;
1338 }
1339 return 0.0;
1340 }
1341
1342
1343
1344
1345
1346
1347 @Override
1348 public void getdEdXdL(double[] gradient) {
1349
1350 }
1351
1352 public void init() {
1353 initializationRegion.init(this, atoms, lambdaTerm, grad, torque, bornRadiiChainRule);
1354 initializationRegion.executeWith(parallelTeam);
1355 }
1356
1357 public void reduce(
1358 AtomicDoubleArray3D g,
1359 AtomicDoubleArray3D t,
1360 AtomicDoubleArray3D lg,
1361 AtomicDoubleArray3D lt) {
1362 grad.reduce(parallelTeam);
1363 torque.reduce(parallelTeam);
1364 for (int i = 0; i < nAtoms; i++) {
1365 g.add(0, i, lPow * grad.getX(i), lPow * grad.getY(i), lPow * grad.getZ(i));
1366 t.add(0, i, lPow * torque.getX(i), lPow * torque.getY(i), lPow * torque.getZ(i));
1367 if (lambdaTerm) {
1368 lg.add(0, i, dlPow * grad.getX(i), dlPow * grad.getY(i), dlPow * grad.getZ(i));
1369 lt.add(0, i, dlPow * torque.getX(i), dlPow * torque.getY(i), dlPow * torque.getZ(i));
1370 }
1371 }
1372 }
1373
1374 public AtomicDoubleArray getSelfEnergy() {
1375
1376 return gkEnergyRegion.getSelfEnergy();
1377 }
1378
1379 public double[] getBorn() {
1380 return bornRadiiRegion.getBorn();
1381 }
1382
1383
1384
1385
1386
1387
1388 public void setElementHCTScaleFactors(HashMap<Integer, Double> elementHCT) {
1389 for (Integer element : elementHCT.keySet()) {
1390 if (elementHCTScaleFactors.containsKey(element)) {
1391 elementHCTScaleFactors.replace(element, elementHCT.get(element));
1392 } else {
1393 logger.info("No HCT scale factor set for element: " + element);
1394 }
1395 }
1396 initAtomArrays();
1397 }
1398
1399
1400
1401
1402
1403
1404 public void setAtoms(Atom[] atoms) {
1405 this.atoms = atoms;
1406 nAtoms = atoms.length;
1407 maxNumAtoms = max(nAtoms, maxNumAtoms);
1408 initAtomArrays();
1409 }
1410
1411
1412
1413
1414
1415
1416 public void setCrystal(Crystal crystal) {
1417 this.crystal = crystal;
1418 }
1419
1420
1421
1422
1423
1424
1425 public void setNeighborList(int[][][] neighbors) {
1426 this.neighborLists = neighbors;
1427 }
1428
1429
1430
1431
1432
1433
1434 public void setUse(boolean[] use) {
1435 this.use = use;
1436 }
1437
1438
1439
1440
1441
1442
1443
1444
1445 public double solvationEnergy(boolean gradient, boolean print) {
1446 return solvationEnergy(0.0, gradient, print);
1447 }
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457 public double solvationEnergy(double gkInducedCorrectionEnergy, boolean gradient, boolean print) {
1458 cavitationEnergy = 0.0;
1459 dispersionEnergy = 0.0;
1460 gkEnergy = 0.0;
1461
1462 try {
1463
1464 gkTime = -System.nanoTime();
1465 gkEnergyRegion.init(
1466 atoms,
1467 globalMultipole,
1468 inducedDipole,
1469 inducedDipoleCR,
1470 crystal,
1471 sXYZ,
1472 neighborLists,
1473 use,
1474 cut2,
1475 baseRadius,
1476 born,
1477 gradient,
1478 parallelTeam,
1479 grad,
1480 torque,
1481 bornRadiiChainRule);
1482 parallelTeam.execute(gkEnergyRegion);
1483 gkTime += System.nanoTime();
1484
1485
1486 switch (nonPolarModel) {
1487 case CAV:
1488 cavitationTime = -System.nanoTime();
1489 parallelTeam.execute(surfaceAreaRegion);
1490 cavitationEnergy = surfaceAreaRegion.getEnergy();
1491 cavitationTime += System.nanoTime();
1492 break;
1493 case CAV_DISP:
1494 dispersionTime = -System.nanoTime();
1495 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1496 parallelTeam.execute(dispersionRegion);
1497 dispersionEnergy = dispersionRegion.getEnergy();
1498 dispersionTime += System.nanoTime();
1499 cavitationTime = -System.nanoTime();
1500 parallelTeam.execute(surfaceAreaRegion);
1501 cavitationEnergy = surfaceAreaRegion.getEnergy();
1502 cavitationTime += System.nanoTime();
1503 break;
1504 case DISP:
1505 dispersionTime = -System.nanoTime();
1506 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1507 parallelTeam.execute(dispersionRegion);
1508 dispersionEnergy = dispersionRegion.getEnergy();
1509 dispersionTime += System.nanoTime();
1510 break;
1511 case SEV_DISP:
1512 dispersionTime = -System.nanoTime();
1513 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1514 parallelTeam.execute(dispersionRegion);
1515 dispersionEnergy = dispersionRegion.getEnergy();
1516 dispersionTime += System.nanoTime();
1517 cavitationTime = -System.nanoTime();
1518 cavitationTime = -System.nanoTime();
1519 double[][] positions = new double[nAtoms][3];
1520 for (int i = 0; i < nAtoms; i++) {
1521 positions[i][0] = atoms[i].getX();
1522 positions[i][1] = atoms[i].getY();
1523 positions[i][2] = atoms[i].getZ();
1524 }
1525 chandlerCavitation.energyAndGradient(positions, grad);
1526 cavitationEnergy = chandlerCavitation.getEnergy();
1527 cavitationTime += System.nanoTime();
1528 break;
1529 case GAUSS_DISP:
1530 dispersionTime = -System.nanoTime();
1531 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1532 parallelTeam.execute(dispersionRegion);
1533 dispersionEnergy = dispersionRegion.getEnergy();
1534 dispersionTime += System.nanoTime();
1535 cavitationTime = -System.nanoTime();
1536 positions = new double[nAtoms][3];
1537 for (int i = 0; i < nAtoms; i++) {
1538 positions[i][0] = atoms[i].getX();
1539 positions[i][1] = atoms[i].getY();
1540 positions[i][2] = atoms[i].getZ();
1541 }
1542 chandlerCavitation.energyAndGradient(positions, grad);
1543 cavitationEnergy = chandlerCavitation.getEnergy();
1544 cavitationTime += System.nanoTime();
1545 break;
1546 case BORN_CAV_DISP:
1547 dispersionTime = -System.nanoTime();
1548 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1549 parallelTeam.execute(dispersionRegion);
1550 dispersionEnergy = dispersionRegion.getEnergy();
1551 dispersionTime += System.nanoTime();
1552 break;
1553 case HYDROPHOBIC_PMF:
1554 case BORN_SOLV:
1555 case NONE:
1556 default:
1557 break;
1558 }
1559 } catch (Exception e) {
1560 String message = "Fatal exception computing the continuum solvation energy.";
1561 logger.log(Level.SEVERE, message, e);
1562 }
1563
1564
1565 if (gradient && !perfectRadii) {
1566 try {
1567 gkTime -= System.nanoTime();
1568 bornGradRegion.init(
1569 atoms, crystal, sXYZ, neighborLists,
1570 baseRadius, descreenRadius, overlapScale, neckScale, descreenOffset,
1571 bornRadiiRegion.getUnscaledBornIntegral(), use, cut2,
1572 nativeEnvironmentApproximation, born, grad, bornRadiiChainRule);
1573 bornGradRegion.executeWith(parallelTeam);
1574 gkTime += System.nanoTime();
1575 } catch (Exception e) {
1576 String message = "Fatal exception computing Born radii chain rule term.";
1577 logger.log(Level.SEVERE, message, e);
1578 }
1579 }
1580
1581 gkEnergy = gkEnergyRegion.getEnergy() + gkInducedCorrectionEnergy;
1582 gkPermanentEnergy = gkEnergyRegion.getPermanentEnergy();
1583 gkPolarizationEnergy = gkEnergy - gkPermanentEnergy;
1584
1585
1586
1587
1588
1589 solvationEnergy = cavitationEnergy + dispersionEnergy + gkEnergy;
1590
1591 if (print) {
1592 logger.info(format(" Generalized Kirkwood%16.8f %10.3f", gkEnergy, gkTime * 1e-9));
1593 switch (nonPolarModel) {
1594 case CAV:
1595 logger.info(format(" Cavitation %16.8f %10.3f", cavitationEnergy, cavitationTime * 1e-9));
1596 break;
1597 case DISP:
1598 logger.info(format(" Dispersion %16.8f %10.3f", dispersionEnergy, dispersionTime * 1e-9));
1599 break;
1600 case CAV_DISP:
1601 case SEV_DISP:
1602 case GAUSS_DISP:
1603 logger.info(format(" Cavitation %16.8f %10.3f", cavitationEnergy, cavitationTime * 1e-9));
1604 logger.info(format(" Dispersion %16.8f %10.3f", dispersionEnergy, dispersionTime * 1e-9));
1605 break;
1606 case BORN_CAV_DISP:
1607 logger.info(format(" Dispersion %16.8f %10.3f", dispersionEnergy, dispersionTime * 1e-9));
1608 break;
1609 case HYDROPHOBIC_PMF:
1610 case BORN_SOLV:
1611 case NONE:
1612 default:
1613 break;
1614 }
1615 }
1616
1617 if (lambdaTerm) {
1618 return lPow * solvationEnergy;
1619 } else {
1620 return solvationEnergy;
1621 }
1622 }
1623
1624 private void initAtomArrays() {
1625 sXYZ = particleMeshEwald.coordinates;
1626
1627 x = sXYZ[0][0];
1628 y = sXYZ[0][1];
1629 z = sXYZ[0][2];
1630
1631 globalMultipole = particleMeshEwald.globalMultipole;
1632 inducedDipole = particleMeshEwald.inducedDipole;
1633 inducedDipoleCR = particleMeshEwald.inducedDipoleCR;
1634 neighborLists = particleMeshEwald.neighborLists;
1635
1636 if (grad == null) {
1637 int threadCount = parallelTeam.getThreadCount();
1638 grad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
1639 torque = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
1640 bornRadiiChainRule = atomicDoubleArrayImpl.createInstance(threadCount, nAtoms);
1641 } else {
1642 grad.alloc(nAtoms);
1643 torque.alloc(nAtoms);
1644 bornRadiiChainRule.alloc(nAtoms);
1645 }
1646
1647 fieldGK.alloc(nAtoms);
1648 fieldGKCR.alloc(nAtoms);
1649
1650 if (baseRadius == null || baseRadius.length < nAtoms) {
1651 baseRadius = new double[nAtoms];
1652 descreenRadius = new double[nAtoms];
1653 overlapScale = new double[nAtoms];
1654 neckScale = new double[nAtoms];
1655 born = new double[nAtoms];
1656 use = new boolean[nAtoms];
1657 }
1658
1659 fill(use, true);
1660
1661 setSoluteRadii(forceField, atoms, soluteRadiiType);
1662 applySoluteRadii();
1663
1664 if (dispersionRegion != null) {
1665 dispersionRegion.allocate(atoms);
1666 }
1667
1668 if (surfaceAreaRegion != null) {
1669 surfaceAreaRegion.init();
1670 }
1671
1672 if (chandlerCavitation != null) {
1673 GaussVol gaussVol = chandlerCavitation.getGaussVol();
1674 try {
1675 gaussVol.allocate(atoms);
1676 } catch (Exception e) {
1677 logger.severe(" " + e);
1678 }
1679 }
1680 }
1681
1682
1683
1684
1685
1686
1687
1688 public void udpateSoluteParameters(int i) {
1689 Atom atom = atoms[i];
1690 SoluteType soluteType = atom.getSoluteType();
1691 if (soluteType == null) {
1692 logger.severe(format(" No SoluteType for atom %s.", atom));
1693 return;
1694 }
1695
1696
1697 baseRadius[i] = soluteType.gkDiameter * 0.5 * bondiScale;
1698
1699 overlapScale[i] = hctScale;
1700
1701 if (elementHCTScale) {
1702 int atomicNumber = atom.getAtomicNumber();
1703 if (elementHCTScaleFactors.containsKey(atomicNumber)) {
1704 overlapScale[i] = elementHCTScaleFactors.get(atomicNumber);
1705 } else {
1706 logger.fine(format(" No element-specific HCT scale factor for %s", atom));
1707 overlapScale[i] = hctScale;
1708 }
1709 }
1710
1711 descreenRadius[i] = baseRadius[i];
1712
1713
1714 AtomType atomType = atom.getAtomType();
1715 if (radiiOverride.containsKey(atomType.type)) {
1716 double override = radiiOverride.get(atomType.type);
1717
1718 baseRadius[i] = baseRadius[i] * override / bondiScale;
1719 logger.fine(format(
1720 " Scaling %s (atom type %d) to %7.4f (Bondi factor %7.4f)",
1721 atom, atomType.type, baseRadius[i], override));
1722 descreenRadius[i] = baseRadius[i];
1723 }
1724
1725
1726 if (descreenVDW) {
1727 descreenRadius[i] = atom.getVDWType().radius / 2.0;
1728 }
1729
1730
1731 if (!descreenHydrogen && atom.getAtomicNumber() == 1) {
1732 overlapScale[i] = 0.0;
1733 }
1734
1735
1736
1737 if (!atom.isHydrogen()) {
1738
1739 if (!neckCorrection || overlapScale[i] == 0.0) {
1740 neckScale[i] = 0.0;
1741 } else {
1742 if (chemicallyAwareSneck) {
1743
1744 int numBoundHeavyAtoms = 0;
1745 for (Atom boundAtom : atom.get12List()) {
1746 if (!boundAtom.isHydrogen()) {
1747 numBoundHeavyAtoms++;
1748 }
1749 }
1750
1751 if (numBoundHeavyAtoms == 0) {
1752
1753 neckScale[i] = 1.0;
1754 } else {
1755 neckScale[i] = atom.getSoluteType().sneck * (5.0 - numBoundHeavyAtoms) / 4.0;
1756 }
1757 } else {
1758
1759 neckScale[i] = atom.getSoluteType().sneck;
1760 }
1761 }
1762
1763
1764
1765
1766 for (Atom boundAtom : atom.get12List()) {
1767 if (boundAtom.isHydrogen()) {
1768 int hydrogenIndex = boundAtom.getIndex();
1769 neckScale[hydrogenIndex - 1] = neckScale[i];
1770 }
1771 }
1772 }
1773 }
1774
1775
1776
1777
1778 public void applySoluteRadii() {
1779
1780 for (int i = 0; i < nAtoms; i++) {
1781 udpateSoluteParameters(i);
1782 }
1783
1784
1785 if (perfectHCTScale) {
1786 double[][] coords = new double[nAtoms][3];
1787 int index = 0;
1788 for (Atom atom : atoms) {
1789 coords[index][0] = atom.getX();
1790 coords[index][1] = atom.getY();
1791 coords[index][2] = atom.getZ();
1792 index++;
1793 }
1794 GaussVol gaussVol = new GaussVol(atoms, forceField, parallelTeam);
1795 gaussVol.computeVolumeAndSA(coords);
1796 double[] selfVolumesFractions = gaussVol.getSelfVolumeFractions();
1797 for (int i = 0; i < nAtoms; i++) {
1798
1799 overlapScale[i] = selfVolumesFractions[i] * hctScale;
1800
1801
1802 if (!descreenHydrogen && atoms[i].getAtomicNumber() == 1) {
1803 overlapScale[i] = 0.0;
1804 }
1805 }
1806 }
1807 }
1808
1809 public void setSneck(double sneck_input) {
1810 this.sneck = sneck_input;
1811 initAtomArrays();
1812 }
1813
1814 public double[] getTanhBetas() {
1815 double[] betas = {beta0, beta1, beta2};
1816 return betas;
1817 }
1818
1819 public void setTanhBetas(double[] betas) {
1820 this.beta0 = betas[0];
1821 this.beta1 = betas[1];
1822 this.beta2 = betas[2];
1823 }
1824
1825 void setLambdaFunction(double lPow, double dlPow, double dl2Pow) {
1826 if (lambdaTerm) {
1827 this.lPow = lPow;
1828 this.dlPow = dlPow;
1829 this.dl2Pow = dl2Pow;
1830 } else {
1831
1832 this.lambda = 1.0;
1833 this.lPow = 1.0;
1834 this.dlPow = 0.0;
1835 this.dl2Pow = 0.0;
1836 }
1837 }
1838
1839 public enum NonPolarModel {
1840 CAV,
1841 CAV_DISP,
1842 DISP,
1843 SEV_DISP,
1844 GAUSS_DISP,
1845 HYDROPHOBIC_PMF,
1846 BORN_CAV_DISP,
1847 BORN_SOLV,
1848 NONE
1849 }
1850 }