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