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 / 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 DISP:
781 tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
782 dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
783 surfaceAreaRegion = null;
784 chandlerCavitation = null;
785 break;
786 case SEV_DISP:
787 tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
788 double[] radii = new double[nAtoms];
789 int index = 0;
790 for (Atom atom : atoms) {
791 radii[index] = atom.getVDWType().radius / 2.0;
792 index++;
793 }
794 ConnollyRegion connollyRegion = new ConnollyRegion(atoms, radii, threadCount);
795
796
797 double wiggle = forceField.getDouble("WIGGLE", ConnollyRegion.DEFAULT_WIGGLE);
798 connollyRegion.setWiggle(wiggle);
799 chandlerCavitation = new ChandlerCavitation(atoms, connollyRegion, forceField);
800 dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
801 surfaceAreaRegion = null;
802 break;
803 case GAUSS_DISP:
804 dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
805 surfaceAreaRegion = null;
806 tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
807 GaussVol gaussVol = new GaussVol(atoms, forceField, parallelTeam);
808 chandlerCavitation = new ChandlerCavitation(atoms, gaussVol, forceField);
809 break;
810 case BORN_CAV_DISP:
811 tensionDefault = DEFAULT_CAVDISP_SURFACE_TENSION;
812 surfaceAreaRegion = null;
813 chandlerCavitation = null;
814 dispersionRegion = new DispersionRegion(threadCount, atoms, forceField);
815 break;
816 case HYDROPHOBIC_PMF:
817 case BORN_SOLV:
818 case NONE:
819 default:
820 tensionDefault = DEFAULT_CAVONLY_SURFACE_TENSION;
821 surfaceAreaRegion = null;
822 dispersionRegion = null;
823 chandlerCavitation = null;
824 break;
825 }
826
827 gkc = forceField.getDouble("GKC", DEFAULT_GKC);
828
829 surfaceTension = forceField.getDouble("SURFACE_TENSION", tensionDefault);
830 solventPressue = forceField.getDouble("SOLVENT_PRESSURE", DEFAULT_SOLVENT_PRESSURE);
831
832 double crossOver = forceField.getDouble("CROSS_OVER", DEFAULT_CROSSOVER);
833 if (chandlerCavitation != null) {
834
835 chandlerCavitation.setCrossOver(crossOver);
836
837 chandlerCavitation.setSurfaceTension(surfaceTension);
838 chandlerCavitation.setSolventPressure(solventPressue);
839 crossOver = chandlerCavitation.getCrossOver();
840 }
841 if (surfaceAreaRegion != null) {
842 surfaceAreaRegion.setSurfaceTension(surfaceTension);
843 }
844 double dispersionOffset = forceField.getDouble("DISPERSION_OFFSET", DEFAULT_DISPERSION_OFFSET);
845 if (dispersionRegion != null) {
846 dispersionRegion.setDispersionOffest(dispersionOffset);
847 }
848
849 initializationRegion = new InitializationRegion(threadCount);
850 bornRadiiRegion = new BornRadiiRegion(threadCount, nAtoms, forceField, neckCorrection,
851 tanhCorrection, perfectHCTScale);
852 permanentGKFieldRegion = new PermanentGKFieldRegion(threadCount, soluteDielectric, solventDielectric, gkc);
853 inducedGKFieldRegion = new InducedGKFieldRegion(threadCount, soluteDielectric, solventDielectric, gkc);
854 if (!perfectRadii) {
855 bornGradRegion = new BornGradRegion(threadCount, neckCorrection, tanhCorrection, perfectHCTScale);
856 } else {
857
858 bornGradRegion = null;
859 }
860 boolean gkQI = forceField.getBoolean("GK_QI", false);
861 gkEnergyRegion = new GKEnergyRegion(threadCount, polarization, nonPolarModel, surfaceTension,
862 probe, electric, soluteDielectric, solventDielectric, gkc, gkQI);
863
864 if (gkQI) {
865 logger.info(" Continuum Solvation (QI)");
866 } else {
867 logger.info(" Continuum Solvation");
868 }
869 logger.info(format(" Radii: %8s", soluteRadiiType));
870 if (cutoff != Double.POSITIVE_INFINITY) {
871 logger.info(format(" Generalized Kirkwood Cut-Off: %8.4f (A)", cutoff));
872 } else {
873 logger.info(" Generalized Kirkwood Cut-Off: NONE");
874 }
875 logger.info(format(" GKC: %8.4f", gkc));
876 logger.info(format(" Solvent Dielectric: %8.4f", solventDielectric));
877 logger.info(format(" Solute Dielectric: %8.4f", soluteDielectric));
878
879 if (perfectRadii) {
880 logger.info(format(" Use Perfect Born Radii: %8B", perfectRadii));
881 } else {
882 logger.info(format(" Descreen with vdW Radii: %8B", descreenVDW));
883 logger.info(format(" Descreen with Hydrogen Atoms: %8B", descreenHydrogen));
884 logger.info(format(" Descreen Offset: %8.4f", descreenOffset));
885 if (neckCorrection) {
886 logger.info(format(" Use Neck Correction: %8B", neckCorrection));
887 logger.info(format(" Sneck Scale Factor: %8.4f", sneck));
888 logger.info(format(" Chemically Aware Sneck: %8B", chemicallyAwareSneck));
889 }
890 if (tanhCorrection) {
891 logger.info(format(" Use Tanh Correction %8B", tanhCorrection));
892 logger.info(format(" Beta0: %8.4f", beta0));
893 logger.info(format(" Beta1: %8.4f", beta1));
894 logger.info(format(" Beta2: %8.4f", beta2));
895 }
896 if (perfectHCTScale) {
897 logger.info(format(" GaussVol HCT Scale Factors: %8B", perfectHCTScale));
898 }
899 logger.info(format(" General HCT Scale Factor: %8.4f",
900 forceField.getDouble("HCT-SCALE", DEFAULT_HCT_SCALE)));
901 if (elementHCTScale) {
902 logger.info(format(" Element-Specific HCT Scale Factors: %8B", elementHCTScale));
903 Integer[] elementHCTkeyset = elementHCTScaleFactors.keySet().toArray(new Integer[0]);
904 Arrays.sort(elementHCTkeyset);
905 for (Integer key : elementHCTkeyset) {
906 logger.info(format(" HCT-Element # %d: %8.4f", key, elementHCTScaleFactors.get(key)));
907 }
908 }
909 }
910
911 logger.info(
912 format(" Non-Polar Model: %10s",
913 nonPolarModel.toString().replace('_', '-')));
914
915 if (nonPolarModel.equals(NonPolarModel.GAUSS_DISP)) {
916 logger.info(format(" GaussVol Radii Offset: %2.4f",
917 forceField.getDouble("GAUSSVOL_RADII_OFFSET", 0.0)));
918 logger.info(format(" GaussVol Radii Scale: %2.4f",
919 forceField.getDouble("GAUSSVOL_RADII_SCALE", 1.15)));
920 }
921
922 if (dispersionRegion != null) {
923 logger.info(format(" Dispersion Integral Offset: %8.4f (A)",
924 dispersionRegion.getDispersionOffset()));
925 }
926
927 if (surfaceAreaRegion != null) {
928 logger.info(format(" Cavitation Probe Radius: %8.4f (A)", probe));
929 logger.info(
930 format(" Cavitation Surface Tension: %8.4f (Kcal/mol/A^2)", surfaceTension));
931 } else if (chandlerCavitation != null) {
932 logger.info(
933 format(" Cavitation Solvent Pressure: %8.4f (Kcal/mol/A^3)", solventPressue));
934 logger.info(
935 format(" Cavitation Surface Tension: %8.4f (Kcal/mol/A^2)", surfaceTension));
936 logger.info(format(" Cavitation Cross-Over Radius: %8.4f (A)", crossOver));
937 }
938
939
940 if (logger.isLoggable(Level.FINE)) {
941 logger.fine(" Base Radii Descreen Radius Overlap Scale Neck Scale");
942 for (int i = 0; i < nAtoms; i++) {
943 logger.info(
944 format(" %s %8.6f %8.6f %5.3f %5.3f",
945 atoms[i].toString(), baseRadius[i], descreenRadius[i], overlapScale[i],
946 neckScale[i]));
947 }
948 }
949 }
950
951
952
953
954
955
956 public boolean getUsePerfectRadii() {
957 return perfectRadii;
958 }
959
960
961
962
963
964
965
966 public double[] getPerfectRadii() {
967 bornRadiiRegion.init(atoms, crystal, sXYZ, neighborLists, baseRadius, descreenRadius,
968 overlapScale, neckScale, descreenOffset, use, cut2, nativeEnvironmentApproximation, born);
969 return bornRadiiRegion.getPerfectRadii();
970 }
971
972
973
974
975
976
977
978 public static NonPolarModel getNonPolarModel(String nonpolarModel) {
979 try {
980 return NonPolarModel.valueOf(toEnumForm(nonpolarModel));
981 } catch (IllegalArgumentException ex) {
982 logger.warning(" Unrecognized nonpolar model requested; defaulting to NONE.");
983 return NonPolarModel.NONE;
984 }
985 }
986
987
988
989
990 public void computeBornRadii() {
991
992 applySoluteRadii();
993
994
995 if (tanhCorrection) {
996 BornTanhRescaling.setBeta0(beta0);
997 BornTanhRescaling.setBeta1(beta1);
998 BornTanhRescaling.setBeta2(beta2);
999 }
1000
1001 try {
1002 bornRadiiRegion.init(
1003 atoms,
1004 crystal,
1005 sXYZ,
1006 neighborLists,
1007 baseRadius,
1008 descreenRadius,
1009 overlapScale,
1010 neckScale,
1011 descreenOffset,
1012 use,
1013 cut2,
1014 nativeEnvironmentApproximation,
1015 born);
1016 parallelTeam.execute(bornRadiiRegion);
1017 } catch (Exception e) {
1018 String message = "Fatal exception computing Born radii.";
1019 logger.log(Level.SEVERE, message, e);
1020 }
1021 }
1022
1023
1024
1025
1026 public void computeInducedGKField() {
1027 try {
1028 fieldGK.reset(parallelTeam);
1029 fieldGKCR.reset(parallelTeam);
1030 inducedGKFieldRegion.init(atoms, inducedDipole, inducedDipoleCR, crystal, sXYZ,
1031 neighborLists, use, cut2, born, fieldGK, fieldGKCR);
1032 parallelTeam.execute(inducedGKFieldRegion);
1033 fieldGK.reduce(parallelTeam);
1034 fieldGKCR.reduce(parallelTeam);
1035 } catch (Exception e) {
1036 String message = "Fatal exception computing induced GK field.";
1037 logger.log(Level.SEVERE, message, e);
1038 }
1039 }
1040
1041
1042
1043
1044 public void computePermanentGKField() {
1045 try {
1046 fieldGK.reset(parallelTeam);
1047 permanentGKFieldRegion.init(
1048 atoms, globalMultipole, crystal, sXYZ, neighborLists, use, cut2, born, fieldGK);
1049 parallelTeam.execute(permanentGKFieldRegion);
1050 fieldGK.reduce(parallelTeam);
1051 } catch (Exception e) {
1052 String message = "Fatal exception computing permanent GK field.";
1053 logger.log(Level.SEVERE, message, e);
1054 }
1055 }
1056
1057
1058
1059
1060
1061
1062 public double[] getBaseRadii() {
1063 return baseRadius;
1064 }
1065
1066
1067
1068
1069
1070
1071 public double[] getDescreenRadii() {
1072 return descreenRadius;
1073 }
1074
1075
1076
1077
1078
1079
1080 public double getCavitationEnergy() {
1081 return cavitationEnergy;
1082 }
1083
1084
1085
1086
1087
1088
1089 public ChandlerCavitation getChandlerCavitation() {
1090 return chandlerCavitation;
1091 }
1092
1093
1094
1095
1096
1097
1098 public double getCutoff() {
1099 return cutoff;
1100 }
1101
1102
1103
1104
1105
1106
1107 public void setCutoff(double cutoff) {
1108 this.cutoff = cutoff;
1109 this.cut2 = cutoff * cutoff;
1110 }
1111
1112
1113
1114
1115
1116
1117 public double getDielecOffset() {
1118 return dOffset;
1119 }
1120
1121
1122
1123
1124
1125
1126 public double getDescreenOffset() {
1127 return descreenOffset;
1128 }
1129
1130
1131
1132
1133
1134
1135 public double getDispersionEnergy() {
1136 return dispersionEnergy;
1137 }
1138
1139 public DispersionRegion getDispersionRegion() {
1140 return dispersionRegion;
1141 }
1142
1143 public AtomicDoubleArray3D getFieldGK() {
1144 return fieldGK;
1145 }
1146
1147 public AtomicDoubleArray3D getFieldGKCR() {
1148 return fieldGKCR;
1149 }
1150
1151 public SurfaceAreaRegion getSurfaceAreaRegion() {
1152 return surfaceAreaRegion;
1153 }
1154
1155
1156
1157
1158
1159
1160 public double getGeneralizedKirkwoordEnergy() {
1161 return gkEnergy;
1162 }
1163
1164
1165
1166
1167
1168
1169 public double getGeneralizedKirkwoordPermanentEnergy() {
1170 return gkPermanentEnergy;
1171 }
1172
1173
1174
1175
1176
1177
1178 public double getGeneralizedKirkwoordPolariztionEnergy() {
1179 return gkPolarizationEnergy;
1180 }
1181
1182 public AtomicDoubleArray3D getGrad() {
1183 return grad;
1184 }
1185
1186
1187
1188
1189
1190
1191 public int getInteractions() {
1192 return gkEnergyRegion.getInteractions();
1193 }
1194
1195
1196
1197
1198 @Override
1199 public double getLambda() {
1200 return lambda;
1201 }
1202
1203
1204
1205
1206
1207
1208 @Override
1209 public void setLambda(double lambda) {
1210 if (lambdaTerm) {
1211 this.lambda = lambda;
1212 } else {
1213
1214 this.lambda = 1.0;
1215 lPow = 1.0;
1216 dlPow = 0.0;
1217 dl2Pow = 0.0;
1218 }
1219 }
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232 public boolean getNativeEnvironmentApproximation() {
1233 return nativeEnvironmentApproximation;
1234 }
1235
1236
1237
1238
1239
1240
1241 public NonPolarModel getNonPolarModel() {
1242 return nonPolarModel;
1243 }
1244
1245
1246
1247
1248
1249
1250 public double[] getOverlapScale() {
1251 return overlapScale;
1252 }
1253
1254
1255
1256
1257
1258
1259 public double[] getNeckScale() {
1260 return neckScale;
1261 }
1262
1263 public boolean getTanhCorrection() {
1264 return tanhCorrection;
1265 }
1266
1267
1268
1269
1270
1271
1272 public double getProbeRadius() {
1273 return probe;
1274 }
1275
1276
1277
1278
1279
1280
1281 public double getSolventPermittivity() {
1282 return solventDielectric;
1283 }
1284
1285
1286
1287
1288
1289
1290 public double getSolutePermittivity() {
1291 return soluteDielectric;
1292 }
1293
1294
1295
1296
1297
1298
1299 public double getSurfaceTension() {
1300 return surfaceTension;
1301 }
1302
1303 public AtomicDoubleArray3D getTorque() {
1304 return torque;
1305 }
1306
1307
1308
1309
1310
1311
1312 @Override
1313 public double getd2EdL2() {
1314 if (lambdaTerm) {
1315 return dl2Pow * solvationEnergy;
1316 } else {
1317 return 0.0;
1318 }
1319 }
1320
1321
1322
1323
1324 @Override
1325 public double getdEdL() {
1326 if (lambdaTerm) {
1327 return dlPow * solvationEnergy;
1328 }
1329 return 0.0;
1330 }
1331
1332
1333
1334
1335
1336
1337 @Override
1338 public void getdEdXdL(double[] gradient) {
1339
1340 }
1341
1342 public void init() {
1343 initializationRegion.init(this, atoms, lambdaTerm, grad, torque, bornRadiiChainRule);
1344 initializationRegion.executeWith(parallelTeam);
1345 }
1346
1347 public void reduce(
1348 AtomicDoubleArray3D g,
1349 AtomicDoubleArray3D t,
1350 AtomicDoubleArray3D lg,
1351 AtomicDoubleArray3D lt) {
1352 grad.reduce(parallelTeam);
1353 torque.reduce(parallelTeam);
1354 for (int i = 0; i < nAtoms; i++) {
1355 g.add(0, i, lPow * grad.getX(i), lPow * grad.getY(i), lPow * grad.getZ(i));
1356 t.add(0, i, lPow * torque.getX(i), lPow * torque.getY(i), lPow * torque.getZ(i));
1357 if (lambdaTerm) {
1358 lg.add(0, i, dlPow * grad.getX(i), dlPow * grad.getY(i), dlPow * grad.getZ(i));
1359 lt.add(0, i, dlPow * torque.getX(i), dlPow * torque.getY(i), dlPow * torque.getZ(i));
1360 }
1361 }
1362 }
1363
1364 public AtomicDoubleArray getSelfEnergy() {
1365
1366 return gkEnergyRegion.getSelfEnergy();
1367 }
1368
1369 public double[] getBorn() {
1370 return bornRadiiRegion.getBorn();
1371 }
1372
1373
1374
1375
1376
1377
1378 public void setElementHCTScaleFactors(HashMap<Integer, Double> elementHCT) {
1379 for (Integer element : elementHCT.keySet()) {
1380 if (elementHCTScaleFactors.containsKey(element)) {
1381 elementHCTScaleFactors.replace(element, elementHCT.get(element));
1382 } else {
1383 logger.info("No HCT scale factor set for element: " + element);
1384 }
1385 }
1386 initAtomArrays();
1387 }
1388
1389
1390
1391
1392
1393
1394 public void setAtoms(Atom[] atoms) {
1395 this.atoms = atoms;
1396 nAtoms = atoms.length;
1397 maxNumAtoms = max(nAtoms, maxNumAtoms);
1398 initAtomArrays();
1399 }
1400
1401
1402
1403
1404
1405
1406 public void setCrystal(Crystal crystal) {
1407 this.crystal = crystal;
1408 }
1409
1410
1411
1412
1413
1414
1415 public void setNeighborList(int[][][] neighbors) {
1416 this.neighborLists = neighbors;
1417 }
1418
1419
1420
1421
1422
1423
1424 public void setUse(boolean[] use) {
1425 this.use = use;
1426 }
1427
1428
1429
1430
1431
1432
1433
1434
1435 public double solvationEnergy(boolean gradient, boolean print) {
1436 return solvationEnergy(0.0, gradient, print);
1437 }
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447 public double solvationEnergy(double gkInducedCorrectionEnergy, boolean gradient, boolean print) {
1448 cavitationEnergy = 0.0;
1449 dispersionEnergy = 0.0;
1450 gkEnergy = 0.0;
1451
1452 try {
1453
1454 gkTime = -System.nanoTime();
1455 gkEnergyRegion.init(
1456 atoms,
1457 globalMultipole,
1458 inducedDipole,
1459 inducedDipoleCR,
1460 crystal,
1461 sXYZ,
1462 neighborLists,
1463 use,
1464 cut2,
1465 baseRadius,
1466 born,
1467 gradient,
1468 parallelTeam,
1469 grad,
1470 torque,
1471 bornRadiiChainRule);
1472 parallelTeam.execute(gkEnergyRegion);
1473 gkTime += System.nanoTime();
1474
1475
1476 switch (nonPolarModel) {
1477 case CAV:
1478 cavitationTime = -System.nanoTime();
1479 parallelTeam.execute(surfaceAreaRegion);
1480 cavitationEnergy = surfaceAreaRegion.getEnergy();
1481 cavitationTime += System.nanoTime();
1482 break;
1483 case CAV_DISP:
1484 dispersionTime = -System.nanoTime();
1485 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1486 parallelTeam.execute(dispersionRegion);
1487 dispersionEnergy = dispersionRegion.getEnergy();
1488 dispersionTime += System.nanoTime();
1489 cavitationTime = -System.nanoTime();
1490 parallelTeam.execute(surfaceAreaRegion);
1491 cavitationEnergy = surfaceAreaRegion.getEnergy();
1492 cavitationTime += System.nanoTime();
1493 break;
1494 case DISP:
1495 dispersionTime = -System.nanoTime();
1496 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1497 parallelTeam.execute(dispersionRegion);
1498 dispersionEnergy = dispersionRegion.getEnergy();
1499 dispersionTime += System.nanoTime();
1500 break;
1501 case SEV_DISP:
1502 dispersionTime = -System.nanoTime();
1503 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1504 parallelTeam.execute(dispersionRegion);
1505 dispersionEnergy = dispersionRegion.getEnergy();
1506 dispersionTime += System.nanoTime();
1507 cavitationTime = -System.nanoTime();
1508 cavitationTime = -System.nanoTime();
1509 double[][] positions = new double[nAtoms][3];
1510 for (int i = 0; i < nAtoms; i++) {
1511 positions[i][0] = atoms[i].getX();
1512 positions[i][1] = atoms[i].getY();
1513 positions[i][2] = atoms[i].getZ();
1514 }
1515 chandlerCavitation.energyAndGradient(positions, grad);
1516 cavitationEnergy = chandlerCavitation.getEnergy();
1517 cavitationTime += System.nanoTime();
1518 break;
1519 case GAUSS_DISP:
1520 dispersionTime = -System.nanoTime();
1521 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1522 parallelTeam.execute(dispersionRegion);
1523 dispersionEnergy = dispersionRegion.getEnergy();
1524 dispersionTime += System.nanoTime();
1525 cavitationTime = -System.nanoTime();
1526 positions = new double[nAtoms][3];
1527 for (int i = 0; i < nAtoms; i++) {
1528 positions[i][0] = atoms[i].getX();
1529 positions[i][1] = atoms[i].getY();
1530 positions[i][2] = atoms[i].getZ();
1531 }
1532 chandlerCavitation.energyAndGradient(positions, grad);
1533 cavitationEnergy = chandlerCavitation.getEnergy();
1534 cavitationTime += System.nanoTime();
1535 break;
1536 case BORN_CAV_DISP:
1537 dispersionTime = -System.nanoTime();
1538 dispersionRegion.init(atoms, crystal, use, neighborLists, x, y, z, cut2, gradient, grad);
1539 parallelTeam.execute(dispersionRegion);
1540 dispersionEnergy = dispersionRegion.getEnergy();
1541 dispersionTime += System.nanoTime();
1542 break;
1543 case HYDROPHOBIC_PMF:
1544 case BORN_SOLV:
1545 case NONE:
1546 default:
1547 break;
1548 }
1549 } catch (Exception e) {
1550 String message = "Fatal exception computing the continuum solvation energy.";
1551 logger.log(Level.SEVERE, message, e);
1552 }
1553
1554
1555 if (gradient && !perfectRadii) {
1556 try {
1557 gkTime -= System.nanoTime();
1558 bornGradRegion.init(
1559 atoms, crystal, sXYZ, neighborLists,
1560 baseRadius, descreenRadius, overlapScale, neckScale, descreenOffset,
1561 bornRadiiRegion.getUnscaledBornIntegral(), use, cut2,
1562 nativeEnvironmentApproximation, born, grad, bornRadiiChainRule);
1563 bornGradRegion.executeWith(parallelTeam);
1564 gkTime += System.nanoTime();
1565 } catch (Exception e) {
1566 String message = "Fatal exception computing Born radii chain rule term.";
1567 logger.log(Level.SEVERE, message, e);
1568 }
1569 }
1570
1571 gkEnergy = gkEnergyRegion.getEnergy() + gkInducedCorrectionEnergy;
1572 gkPermanentEnergy = gkEnergyRegion.getPermanentEnergy();
1573 gkPolarizationEnergy = gkEnergy - gkPermanentEnergy;
1574
1575
1576
1577
1578
1579 solvationEnergy = cavitationEnergy + dispersionEnergy + gkEnergy;
1580
1581 if (print) {
1582 logger.info(format(" Generalized Kirkwood%16.8f %10.3f", gkEnergy, gkTime * 1e-9));
1583 switch (nonPolarModel) {
1584 case CAV:
1585 logger.info(format(" Cavitation %16.8f %10.3f", cavitationEnergy, cavitationTime * 1e-9));
1586 break;
1587 case DISP:
1588 logger.info(format(" Dispersion %16.8f %10.3f", dispersionEnergy, dispersionTime * 1e-9));
1589 break;
1590 case CAV_DISP:
1591 case SEV_DISP:
1592 case GAUSS_DISP:
1593 logger.info(format(" Cavitation %16.8f %10.3f", cavitationEnergy, cavitationTime * 1e-9));
1594 logger.info(format(" Dispersion %16.8f %10.3f", dispersionEnergy, dispersionTime * 1e-9));
1595 break;
1596 case BORN_CAV_DISP:
1597 logger.info(format(" Dispersion %16.8f %10.3f", dispersionEnergy, dispersionTime * 1e-9));
1598 break;
1599 case HYDROPHOBIC_PMF:
1600 case BORN_SOLV:
1601 case NONE:
1602 default:
1603 break;
1604 }
1605 }
1606
1607 if (lambdaTerm) {
1608 return lPow * solvationEnergy;
1609 } else {
1610 return solvationEnergy;
1611 }
1612 }
1613
1614 private void initAtomArrays() {
1615 sXYZ = particleMeshEwald.coordinates;
1616
1617 x = sXYZ[0][0];
1618 y = sXYZ[0][1];
1619 z = sXYZ[0][2];
1620
1621 globalMultipole = particleMeshEwald.globalMultipole;
1622 inducedDipole = particleMeshEwald.inducedDipole;
1623 inducedDipoleCR = particleMeshEwald.inducedDipoleCR;
1624 neighborLists = particleMeshEwald.neighborLists;
1625
1626 if (grad == null) {
1627 int threadCount = parallelTeam.getThreadCount();
1628 grad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
1629 torque = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, threadCount);
1630 bornRadiiChainRule = atomicDoubleArrayFactory(atomicDoubleArrayImpl, threadCount, nAtoms);
1631 } else {
1632 grad.alloc(nAtoms);
1633 torque.alloc(nAtoms);
1634 bornRadiiChainRule.alloc(nAtoms);
1635 }
1636
1637 fieldGK.alloc(nAtoms);
1638 fieldGKCR.alloc(nAtoms);
1639
1640 if (baseRadius == null || baseRadius.length < nAtoms) {
1641 baseRadius = new double[nAtoms];
1642 descreenRadius = new double[nAtoms];
1643 overlapScale = new double[nAtoms];
1644 neckScale = new double[nAtoms];
1645 born = new double[nAtoms];
1646 use = new boolean[nAtoms];
1647 }
1648
1649 fill(use, true);
1650
1651 setSoluteRadii(forceField, atoms, soluteRadiiType);
1652 applySoluteRadii();
1653
1654 if (dispersionRegion != null) {
1655 dispersionRegion.allocate(atoms);
1656 }
1657
1658 if (surfaceAreaRegion != null) {
1659 surfaceAreaRegion.init();
1660 }
1661
1662 if (chandlerCavitation != null) {
1663 GaussVol gaussVol = chandlerCavitation.getGaussVol();
1664 try {
1665 gaussVol.allocate(atoms);
1666 } catch (Exception e) {
1667 logger.severe(" " + e);
1668 }
1669 }
1670 }
1671
1672
1673
1674
1675
1676
1677
1678 public void udpateSoluteParameters(int i) {
1679 Atom atom = atoms[i];
1680 SoluteType soluteType = atom.getSoluteType();
1681 if (soluteType == null) {
1682 logger.severe(format(" No SoluteType for atom %s.", atom));
1683 return;
1684 }
1685
1686
1687 baseRadius[i] = soluteType.gkDiameter * 0.5 * bondiScale;
1688
1689 overlapScale[i] = hctScale;
1690
1691 if (elementHCTScale) {
1692 int atomicNumber = atom.getAtomicNumber();
1693 if (elementHCTScaleFactors.containsKey(atomicNumber)) {
1694 overlapScale[i] = elementHCTScaleFactors.get(atomicNumber);
1695 } else {
1696 logger.fine(format(" No element-specific HCT scale factor for %s", atom));
1697 overlapScale[i] = hctScale;
1698 }
1699 }
1700
1701 descreenRadius[i] = baseRadius[i];
1702
1703
1704 AtomType atomType = atom.getAtomType();
1705 if (radiiOverride.containsKey(atomType.type)) {
1706 double override = radiiOverride.get(atomType.type);
1707
1708 baseRadius[i] = baseRadius[i] * override / bondiScale;
1709 logger.fine(format(
1710 " Scaling %s (atom type %d) to %7.4f (Bondi factor %7.4f)",
1711 atom, atomType.type, baseRadius[i], override));
1712 descreenRadius[i] = baseRadius[i];
1713 }
1714
1715
1716 if (descreenVDW) {
1717 descreenRadius[i] = atom.getVDWType().radius / 2.0;
1718 }
1719
1720
1721 if (!descreenHydrogen && atom.getAtomicNumber() == 1) {
1722 overlapScale[i] = 0.0;
1723 }
1724
1725
1726
1727 if (!atom.isHydrogen()) {
1728
1729 if (!neckCorrection || overlapScale[i] == 0.0) {
1730 neckScale[i] = 0.0;
1731 } else {
1732 if (chemicallyAwareSneck) {
1733
1734 int numBoundHeavyAtoms = 0;
1735 for (Atom boundAtom : atom.get12List()) {
1736 if (!boundAtom.isHydrogen()) {
1737 numBoundHeavyAtoms++;
1738 }
1739 }
1740
1741 if (numBoundHeavyAtoms == 0) {
1742
1743 neckScale[i] = 1.0;
1744 } else {
1745 neckScale[i] = atom.getSoluteType().sneck * (5.0 - numBoundHeavyAtoms) / 4.0;
1746 }
1747 } else {
1748
1749 neckScale[i] = atom.getSoluteType().sneck;
1750 }
1751 }
1752
1753
1754
1755
1756 for (Atom boundAtom : atom.get12List()) {
1757 if (boundAtom.isHydrogen()) {
1758 int hydrogenIndex = boundAtom.getIndex();
1759 neckScale[hydrogenIndex - 1] = neckScale[i];
1760 }
1761 }
1762 }
1763 }
1764
1765
1766
1767
1768 public void applySoluteRadii() {
1769
1770 for (int i = 0; i < nAtoms; i++) {
1771 udpateSoluteParameters(i);
1772 }
1773
1774
1775 if (perfectHCTScale) {
1776 double[][] coords = new double[nAtoms][3];
1777 int index = 0;
1778 for (Atom atom : atoms) {
1779 coords[index][0] = atom.getX();
1780 coords[index][1] = atom.getY();
1781 coords[index][2] = atom.getZ();
1782 index++;
1783 }
1784 GaussVol gaussVol = new GaussVol(atoms, forceField, parallelTeam);
1785 gaussVol.computeVolumeAndSA(coords);
1786 double[] selfVolumesFractions = gaussVol.getSelfVolumeFractions();
1787 for (int i = 0; i < nAtoms; i++) {
1788
1789 overlapScale[i] = selfVolumesFractions[i] * hctScale;
1790
1791
1792 if (!descreenHydrogen && atoms[i].getAtomicNumber() == 1) {
1793 overlapScale[i] = 0.0;
1794 }
1795 }
1796 }
1797 }
1798
1799 public void setSneck(double sneck_input) {
1800 this.sneck = sneck_input;
1801 initAtomArrays();
1802 }
1803
1804 public double[] getTanhBetas() {
1805 double[] betas = {beta0, beta1, beta2};
1806 return betas;
1807 }
1808
1809 public void setTanhBetas(double[] betas) {
1810 this.beta0 = betas[0];
1811 this.beta1 = betas[1];
1812 this.beta2 = betas[2];
1813 }
1814
1815 void setLambdaFunction(double lPow, double dlPow, double dl2Pow) {
1816 if (lambdaTerm) {
1817 this.lPow = lPow;
1818 this.dlPow = dlPow;
1819 this.dl2Pow = dl2Pow;
1820 } else {
1821
1822 this.lambda = 1.0;
1823 this.lPow = 1.0;
1824 this.dlPow = 0.0;
1825 this.dl2Pow = 0.0;
1826 }
1827 }
1828
1829 public enum NonPolarModel {
1830 CAV,
1831 CAV_DISP,
1832 DISP,
1833 SEV_DISP,
1834 GAUSS_DISP,
1835 HYDROPHOBIC_PMF,
1836 BORN_CAV_DISP,
1837 BORN_SOLV,
1838 NONE
1839 }
1840 }