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.implicit;
39
40 import static ffx.numerics.atomic.AtomicDoubleArray.atomicDoubleArrayFactory;
41 import static ffx.numerics.math.DoubleMath.add;
42 import static ffx.numerics.math.DoubleMath.length2;
43 import static ffx.numerics.math.DoubleMath.scale;
44 import static ffx.numerics.math.DoubleMath.sub;
45 import static java.lang.Double.compare;
46 import static java.lang.String.format;
47 import static org.apache.commons.math3.util.FastMath.PI;
48 import static org.apache.commons.math3.util.FastMath.cbrt;
49 import static org.apache.commons.math3.util.FastMath.exp;
50 import static org.apache.commons.math3.util.FastMath.pow;
51
52 import edu.rit.pj.IntegerForLoop;
53 import edu.rit.pj.ParallelRegion;
54 import edu.rit.pj.ParallelTeam;
55 import edu.rit.pj.reduction.SharedDouble;
56 import ffx.numerics.atomic.AtomicDoubleArray;
57 import ffx.numerics.atomic.AtomicDoubleArray.AtomicDoubleArrayImpl;
58 import ffx.numerics.atomic.AtomicDoubleArray3D;
59 import ffx.potential.bonded.Atom;
60 import ffx.potential.parameters.ForceField;
61
62 import java.util.ArrayList;
63 import java.util.Collections;
64 import java.util.List;
65 import java.util.logging.Level;
66 import java.util.logging.Logger;
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82 public class GaussVol {
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110 private static final Logger logger = Logger.getLogger(GaussVol.class.getName());
111
112
113
114 private static final double offset = 0.005;
115
116
117
118 private static final double KFC = 2.2269859253;
119
120 private static final double PFC = 2.5;
121
122
123
124 private static final double sphereConversion = KFC;
125
126
127
128 private static final int MAX_ORDER = 16;
129
130
131
132 private static final double ANG3 = 1.0;
133
134
135
136 private static final double VOLMINA = 0.01 * ANG3;
137
138
139
140 private static final double VOLMINB = 0.1 * ANG3;
141
142
143
144 private static final double MIN_GVOL = Double.MIN_VALUE;
145
146
147
148
149
150 public static final double DEFAULT_GAUSSVOL_RADII_OFFSET = 0.0;
151
152
153
154
155
156
157
158
159 public static final double DEFAULT_GAUSSVOL_RADII_SCALE = 1.15;
160
161 private static final double RMIN_TO_SIGMA = 1.0 / pow(2.0, 1.0 / 6.0);
162
163
164 private final GaussVolRegion gaussVolRegion;
165
166 private enum GAUSSVOL_MODE {COMPUTE_TREE, RESCAN_TREE}
167
168 private final AtomicDoubleArray3D grad;
169 private final SharedDouble totalVolume;
170 private final SharedDouble energy;
171 private final AtomicDoubleArray gradV;
172 private final AtomicDoubleArray freeVolume;
173 private final AtomicDoubleArray selfVolume;
174 private final ParallelTeam parallelTeam;
175
176 private final double vdwRadiiOffset;
177 private final double vdwRadiiScale;
178 private final boolean includeHydrogen;
179 private final boolean useSigma;
180 private static final double FOUR_THIRDS_PI = 4.0 / 3.0 * PI;
181
182
183
184
185 private final Atom[] atoms;
186
187
188
189 private final int nAtoms;
190
191
192
193 private double[] radii;
194
195
196
197 private double[] volumes;
198
199
200
201
202 private final double[] radiiOffset;
203
204
205
206 private final double[] volumeOffset;
207
208
209
210 private final double[] gammas;
211
212
213
214
215 private final boolean[] ishydrogen;
216
217
218
219 private final double[] selfVolumeFraction;
220
221
222
223 private final GaussianOverlapTree tree;
224
225
226
227 private int maximumDepth = 0;
228
229
230
231 private int totalNumberOfOverlaps = 0;
232
233
234
235 private double surfaceArea;
236
237
238
239 private double volume;
240
241
242
243 private double[] volumeGradient;
244
245
246
247 private double[] surfaceAreaGradient;
248
249
250
251
252
253
254
255
256 public GaussVol(Atom[] atoms, ForceField forceField, ParallelTeam parallelTeam) {
257 this.atoms = atoms;
258 nAtoms = atoms.length;
259 tree = new GaussianOverlapTree(nAtoms);
260 radii = new double[nAtoms];
261 volumes = new double[nAtoms];
262 gammas = new double[nAtoms];
263 ishydrogen = new boolean[nAtoms];
264 radiiOffset = new double[nAtoms];
265 volumeOffset = new double[nAtoms];
266 selfVolumeFraction = new double[nAtoms];
267
268 vdwRadiiOffset = forceField.getDouble("GAUSSVOL_RADII_OFFSET", DEFAULT_GAUSSVOL_RADII_OFFSET);
269 vdwRadiiScale = forceField.getDouble("GAUSSVOL_RADII_SCALE", DEFAULT_GAUSSVOL_RADII_SCALE);
270 includeHydrogen = forceField.getBoolean("GAUSSVOL_HYDROGEN", false);
271 useSigma = forceField.getBoolean("GAUSSVOL_USE_SIGMA", false);
272
273 for (int i = 0; i < nAtoms; i++) {
274 updateAtom(i);
275 }
276
277 this.parallelTeam = parallelTeam;
278 int nThreads = 1;
279 if (parallelTeam != null) {
280 nThreads = parallelTeam.getThreadCount();
281 gaussVolRegion = new GaussVolRegion(nThreads);
282 } else {
283 gaussVolRegion = null;
284 }
285
286 totalVolume = new SharedDouble();
287 energy = new SharedDouble();
288 AtomicDoubleArrayImpl atomicDoubleArrayImpl = AtomicDoubleArrayImpl.MULTI;
289 grad = new AtomicDoubleArray3D(atomicDoubleArrayImpl, nAtoms, nThreads);
290 gradV = atomicDoubleArrayFactory(atomicDoubleArrayImpl, nThreads, nAtoms);
291 freeVolume = atomicDoubleArrayFactory(atomicDoubleArrayImpl, nThreads, nAtoms);
292 selfVolume = atomicDoubleArrayFactory(atomicDoubleArrayImpl, nThreads, nAtoms);
293 }
294
295
296
297
298
299
300 public double[] getRadii() {
301 return radii;
302 }
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324 private static double overlapGaussianAlpha(
325 GaussianVca g1, GaussianVca g2, GaussianVca g12, double[] dVdr, double[] dVdV, double[] sfp) {
326 double[] c1 = g1.c;
327 double[] c2 = g2.c;
328 double[] dist = new double[3];
329
330 sub(c2, c1, dist);
331 double d2 = length2(dist);
332 double a12 = g1.a + g2.a;
333 double deltai = 1.0 / a12;
334
335
336 double df = (g1.a) * (g2.a) * deltai;
337 double ef = exp(-df * d2);
338 double gvol = ((g1.v * g2.v) / pow(PI / df, 1.5)) * ef;
339
340
341 double dgvol = -2.0 * df * gvol;
342 dVdr[0] = dgvol;
343
344
345 double dgvolv = g1.v > 0.0 ? gvol / g1.v : 0.0;
346 dVdV[0] = dgvolv;
347
348
349
350
351
352
353 double[] c1a = new double[3];
354 double[] c2a = new double[3];
355 scale(c1, g1.a * deltai, c1a);
356 scale(c2, g2.a * deltai, c2a);
357 add(c1a, c2a, g12.c);
358 g12.a = a12;
359 g12.v = gvol;
360
361
362 double[] sp = new double[1];
363 double s = switchingFunction(gvol, VOLMINA, VOLMINB, sp);
364 sfp[0] = sp[0] * gvol + s;
365 return s * gvol;
366 }
367
368
369
370
371
372
373
374
375
376
377 private static double switchingFunction(
378 double gvol, double volmina, double volminb, double[] sp) {
379 if (gvol > volminb) {
380 sp[0] = 0.0;
381 return 1.0;
382 } else if (gvol < volmina) {
383 sp[0] = 0.0;
384 return 0.0;
385 }
386 double swd = 1.0 / (volminb - volmina);
387 double swu = (gvol - volmina) * swd;
388 double swu2 = swu * swu;
389 double swu3 = swu * swu2;
390 sp[0] = swd * 30.0 * swu2 * (1.0 - 2.0 * swu + swu2);
391 return swu3 * (10.0 - 15.0 * swu + 6.0 * swu2);
392 }
393
394
395
396
397
398
399 public double[] getSelfVolumeFractions() {
400 return selfVolumeFraction;
401 }
402
403
404
405
406
407
408
409 public double computeVolumeAndSA(double[][] positions) {
410
411 if (parallelTeam == null || parallelTeam.getThreadCount() == 1) {
412
413 for (int i = 0; i < nAtoms - 1; i++) {
414 updateAtom(i);
415 }
416
417
418 computeTree(positions);
419
420
421 computeVolume(totalVolume, energy, grad, gradV, freeVolume, selfVolume);
422 } else {
423
424 try {
425 gaussVolRegion.init(GAUSSVOL_MODE.COMPUTE_TREE, positions);
426 parallelTeam.execute(gaussVolRegion);
427 } catch (Exception e) {
428 logger.severe(" Exception evaluating GaussVol " + e);
429 }
430 }
431
432 if (volumeGradient == null || volumeGradient.length != nAtoms * 3) {
433 volumeGradient = new double[nAtoms * 3];
434 surfaceAreaGradient = new double[nAtoms * 3];
435 }
436
437 double selfVolumeSum = 0.0;
438 double noOverlapVolume = 0.0;
439 for (int i = 0; i < nAtoms; i++) {
440 double si = selfVolume.get(i);
441 selfVolumeSum += si;
442 double x = grad.getX(i);
443 double y = grad.getY(i);
444 double z = grad.getZ(i);
445 if (ishydrogen[i]) {
446 x = 0.0;
447 y = 0.0;
448 z = 0.0;
449 selfVolumeFraction[i] = 0.0;
450 } else {
451 double vi = volumes[i];
452 noOverlapVolume += vi;
453 selfVolumeFraction[i] = si / vi;
454 }
455
456 int index = i * 3;
457 volumeGradient[index++] = x;
458 volumeGradient[index++] = y;
459 volumeGradient[index] = z;
460
461
462 index = i * 3;
463 surfaceAreaGradient[index++] = x;
464 surfaceAreaGradient[index++] = y;
465 surfaceAreaGradient[index] = z;
466 }
467
468
469 if (logger.isLoggable(Level.FINE)) {
470 double meanHCT = selfVolumeSum / noOverlapVolume;
471 logger.fine(format("\n Mean Self-Volume Fraction: %16.8f", meanHCT));
472 logger.fine(format(" Cube Root of the Fraction: %16.8f", cbrt(meanHCT)));
473 if (logger.isLoggable(Level.FINEST)) {
474 logger.finest(" Atomic Self-volume Fractions");
475 for (int i = 0; i < nAtoms; i++) {
476 logger.finest(format(" %6d %16.8f", i + 1, selfVolumeFraction[i]));
477 }
478 }
479 }
480
481
482 volume = totalVolume.get();
483
484
485 double[] radiiBak = this.radii;
486 double[] volumeBak = this.volumes;
487
488
489 this.radii = radiiOffset;
490 this.volumes = volumeOffset;
491
492
493
494
495 if (parallelTeam == null || parallelTeam.getThreadCount() == 1) {
496
497
498 rescanTreeVolumes(positions);
499
500 computeVolume(totalVolume, energy, grad, gradV, freeVolume, selfVolume);
501 } else {
502
503 try {
504 gaussVolRegion.init(GAUSSVOL_MODE.RESCAN_TREE, positions);
505 parallelTeam.execute(gaussVolRegion);
506 } catch (Exception e) {
507 logger.severe(" Exception evaluating GaussVol " + e);
508 }
509 }
510
511
512 this.radii = radiiBak;
513 this.volumes = volumeBak;
514
515 double selfVolumeOffsetSum = 0;
516 for (int i = 0; i < nAtoms; i++) {
517 selfVolumeOffsetSum += selfVolume.get(i);
518 double x = grad.getX(i);
519 double y = grad.getY(i);
520 double z = grad.getZ(i);
521 if (ishydrogen[i]) {
522 x = 0.0;
523 y = 0.0;
524 z = 0.0;
525 }
526
527 int index = i * 3;
528 surfaceAreaGradient[index] = (x - surfaceAreaGradient[index++]) / offset;
529 surfaceAreaGradient[index] = (y - surfaceAreaGradient[index++]) / offset;
530 surfaceAreaGradient[index] = (z - surfaceAreaGradient[index]) / offset;
531 }
532
533
534 surfaceArea = (selfVolumeOffsetSum - selfVolumeSum) / offset;
535
536 return volume;
537 }
538
539
540
541
542
543
544 public int getMaximumDepth() {
545 return maximumDepth;
546 }
547
548
549
550
551
552
553 public double getSurfaceArea() {
554 return surfaceArea;
555 }
556
557 public double[] getSurfaceAreaGradient() {
558 return surfaceAreaGradient;
559 }
560
561
562
563
564
565
566 public int getTotalNumberOfOverlaps() {
567 return totalNumberOfOverlaps;
568 }
569
570
571
572
573
574
575 public double getVolume() {
576 return volume;
577 }
578
579 public double[] getVolumeGradient() {
580 return volumeGradient;
581 }
582
583 public void allocate(Atom[] atoms) throws Exception {
584 if (atoms.length != this.atoms.length) {
585 throw new Exception(" allocate: number of atoms does not match");
586 }
587 }
588
589
590
591
592
593
594 private void computeTree(double[][] positions) {
595 tree.computeOverlapTreeR(positions, radii, volumes, gammas, ishydrogen);
596 }
597
598
599
600
601
602
603
604
605
606
607
608
609 private void computeVolume(
610 SharedDouble totalVolume,
611 SharedDouble totalEnergy,
612 AtomicDoubleArray3D grad,
613 AtomicDoubleArray gradV,
614 AtomicDoubleArray freeVolume,
615 AtomicDoubleArray selfVolume) {
616
617 totalVolume.set(0.0);
618 totalEnergy.set(0.0);
619 grad.reset(0, 0, nAtoms - 1);
620 gradV.reset(0, 0, nAtoms - 1);
621 freeVolume.reset(0, 0, nAtoms - 1);
622 selfVolume.reset(0, 0, nAtoms - 1);
623
624
625 tree.computeVolume2R(totalVolume, totalEnergy, grad, gradV, freeVolume, selfVolume);
626
627
628 grad.reduce(0, nAtoms - 1);
629 gradV.reduce(0, nAtoms - 1);
630 freeVolume.reduce(0, nAtoms - 1);
631 selfVolume.reduce(0, nAtoms - 1);
632
633 for (int i = 0; i < nAtoms; ++i) {
634 if (volumes[i] > 0) {
635 gradV.set(0, i, gradV.get(i) / volumes[i]);
636 }
637 }
638 }
639
640
641
642
643
644
645 private void rescanTreeVolumes(double[][] positions) {
646 tree.rescanTreeV(positions, radii, volumes, gammas, ishydrogen);
647 }
648
649
650
651
652 private void rescanTreeGammas() {
653 tree.rescanTreeG(gammas);
654 }
655
656
657
658
659
660
661 private void getStats(int[] nov) {
662 if (nov.length != nAtoms) {
663 return;
664 }
665
666 for (int i = 0; i < nAtoms; i++) {
667 nov[i] = 0;
668 }
669 for (int atom = 0; atom < nAtoms; atom++) {
670 int slot = atom + 1;
671 nov[atom] = tree.nChildrenUnderSlotR(slot);
672 }
673 }
674
675
676
677
678 private void printTree() {
679 tree.printTree();
680 }
681
682
683
684
685 private static class GaussianVca {
686
687
688 public double v;
689
690 public double a;
691
692 public double[] c = new double[3];
693 }
694
695
696
697
698
699
700
701
702
703
704
705 private static class GaussianOverlap implements Comparable<GaussianOverlap> {
706
707
708
709
710 public int level;
711
712
713
714 public GaussianVca g;
715
716
717
718 public double volume;
719
720
721
722 public int atom;
723
724
725
726 double dvv1;
727
728
729
730 double[] dv1;
731
732
733
734 double gamma1i;
735
736
737
738 double selfVolume;
739
740
741
742 double sfp;
743
744
745
746 int parentIndex;
747
748
749
750 int childrenStartIndex;
751
752
753
754 int childrenCount;
755
756 public GaussianOverlap() {
757 g = new GaussianVca();
758 dv1 = new double[3];
759 }
760
761 public GaussianOverlap(GaussianVca g, double volume, double selfVolume, int atom) {
762 this.g = g;
763 this.volume = volume;
764 this.selfVolume = selfVolume;
765 this.atom = atom;
766
767 dv1 = new double[3];
768 }
769
770
771
772
773
774
775
776 @Override
777 public int compareTo(GaussianOverlap o) {
778 return compare(volume, o.volume);
779 }
780
781
782
783
784 public void printOverlap() {
785 logger.info(format(" Gaussian Overlap %d: Atom: %d, Parent: %d, ChildrenStartIndex: %d, ChildrenCount: %d,"
786 + "Volume: %6.3f, Gamma: %6.3f, Gauss.a: %6.3f, Gauss.v: %6.3f, Gauss.center (%6.3f,%6.3f,%6.3f),"
787 + "dedx: %6.3f, dedy: %6.3f, dedz: %6.3f, sfp: %6.3f",
788 level, atom, parentIndex, childrenStartIndex, childrenCount, volume, gamma1i,
789 g.a, g.v, g.c[0], g.c[1], g.c[2], dv1[0], dv1[1], dv1[2], sfp));
790 }
791 }
792
793 public void updateAtom(int i) {
794 Atom atom = atoms[i];
795 ishydrogen[i] = atom.isHydrogen();
796 if (includeHydrogen) {
797 ishydrogen[i] = false;
798 }
799 radii[i] = atom.getVDWType().radius / 2.0;
800 if (useSigma) {
801 radii[i] *= RMIN_TO_SIGMA;
802 }
803 radii[i] += vdwRadiiOffset;
804 radii[i] *= vdwRadiiScale;
805 volumes[i] = FOUR_THIRDS_PI * pow(radii[i], 3);
806 gammas[i] = 1.0;
807 radiiOffset[i] = radii[i] + offset;
808 volumeOffset[i] = FOUR_THIRDS_PI * pow(radiiOffset[i], 3);
809 }
810
811
812
813
814 private class GaussianOverlapTree {
815
816
817
818
819 int nAtoms;
820
821
822
823 List<GaussianOverlap> overlaps;
824
825
826
827
828
829
830 GaussianOverlapTree(int nAtoms) {
831 this.nAtoms = nAtoms;
832 overlaps = Collections.synchronizedList(new ArrayList<>(nAtoms + 1));
833 }
834
835
836
837
838
839
840
841
842
843
844 void initOverlapTree(
845 double[][] pos, double[] radii, double[] volumes, double[] gammas, boolean[] ishydrogen) {
846
847
848 overlaps = Collections.synchronizedList(new ArrayList<>(nAtoms + 1));
849
850
851 GaussianOverlap overlap = new GaussianOverlap();
852 overlap.level = 0;
853 overlap.volume = 0;
854 overlap.dvv1 = 0.;
855 overlap.selfVolume = 0;
856 overlap.sfp = 1.;
857 overlap.gamma1i = 0.;
858 overlap.parentIndex = -1;
859 overlap.atom = -1;
860 overlap.childrenStartIndex = 1;
861 overlap.childrenCount = nAtoms;
862 overlaps.add(0, overlap);
863
864
865 for (int iat = 0; iat < nAtoms; iat++) {
866 overlap = new GaussianOverlap();
867 double a = sphereConversion / (radii[iat] * radii[iat]);
868 double vol = ishydrogen[iat] ? 0 : volumes[iat];
869 overlap.level = 1;
870 overlap.g.v = vol;
871 overlap.g.a = a;
872 overlap.g.c = pos[iat];
873 overlap.volume = vol;
874 overlap.dvv1 = 1;
875 overlap.selfVolume = 0;
876 overlap.sfp = 1;
877 overlap.gamma1i = gammas[iat];
878 overlap.parentIndex = 0;
879 overlap.atom = iat;
880 overlap.childrenStartIndex = -1;
881 overlap.childrenCount = -1;
882 overlaps.add(iat + 1, overlap);
883 }
884 }
885
886
887
888
889
890
891
892
893 int addChildren(int parentIndex, List<GaussianOverlap> childrenOverlaps) {
894
895
896 int startIndex = overlaps.size();
897
898 int noverlaps = childrenOverlaps.size();
899
900
901 GaussianOverlap root = overlaps.get(parentIndex);
902
903
904 root.childrenStartIndex = startIndex;
905 root.childrenCount = noverlaps;
906
907
908 Collections.sort(childrenOverlaps);
909
910 int rootLevel = root.level;
911 int nextLevel = rootLevel + 1;
912
913
914 for (GaussianOverlap child : childrenOverlaps) {
915 child.level = nextLevel;
916
917
918 child.parentIndex = parentIndex;
919
920
921 child.childrenStartIndex = -1;
922 child.childrenCount = -1;
923
924
925 overlaps.add(child);
926 }
927
928 return startIndex;
929 }
930
931
932
933
934
935
936
937
938 void computeChildren(int rootIndex, List<GaussianOverlap> childrenOverlaps) {
939 int parentIndex;
940 int siblingStart, siblingCount;
941
942
943 childrenOverlaps.clear();
944
945
946 GaussianOverlap root = overlaps.get(rootIndex);
947
948
949 parentIndex = root.parentIndex;
950
951
952 if (parentIndex < 0) {
953 throw new IllegalArgumentException(" Cannot compute children of master node!");
954 }
955
956 if (root.level >= MAX_ORDER) {
957 return;
958 }
959
960 GaussianOverlap parent = overlaps.get(parentIndex);
961
962
963 siblingStart = parent.childrenStartIndex;
964 siblingCount = parent.childrenCount;
965
966
967 if (siblingStart < 0 || siblingCount < 0) {
968 throw new IllegalArgumentException(
969 format(" Parent %s of overlap %s has no sibilings.", parent, root));
970 }
971
972
973 if (rootIndex < siblingStart && rootIndex > siblingStart + siblingCount - 1) {
974 throw new IllegalArgumentException(
975 format(" Node %s is somehow not the child of its parent %s", root, parent));
976 }
977
978
979
980 for (int slotj = rootIndex + 1; slotj < siblingStart + siblingCount; slotj++) {
981
982 GaussianVca g12 = new GaussianVca();
983
984 GaussianOverlap sibling = overlaps.get(slotj);
985 double gvol;
986 double[] dVdr = new double[1];
987 double[] dVdV = new double[1];
988 double[] sfp = new double[1];
989
990
991 int atom2 = sibling.atom;
992 GaussianVca g1 = root.g;
993
994
995 GaussianVca g2 = overlaps.get(atom2 + 1).g;
996 gvol = overlapGaussianAlpha(g1, g2, g12, dVdr, dVdV, sfp);
997
998
999
1000
1001
1002 if (gvol > MIN_GVOL) {
1003 GaussianOverlap ov = new GaussianOverlap(g12, gvol, 0.0, atom2);
1004
1005
1006 sub(g2.c, g1.c, ov.dv1);
1007 scale(ov.dv1, -dVdr[0], ov.dv1);
1008
1009
1010 ov.dvv1 = dVdV[0];
1011 ov.sfp = sfp[0];
1012 ov.gamma1i = root.gamma1i + overlaps.get(atom2 + 1).gamma1i;
1013 childrenOverlaps.add(ov);
1014 }
1015 }
1016 }
1017
1018
1019
1020
1021
1022
1023 private void computeAndAddChildrenR(int root) {
1024 List<GaussianOverlap> childrenOverlaps = Collections.synchronizedList(new ArrayList<>());
1025 computeChildren(root, childrenOverlaps);
1026 int nOverlaps = childrenOverlaps.size();
1027 if (nOverlaps > 0) {
1028 int startSlot = addChildren(root, childrenOverlaps);
1029 for (int ichild = startSlot; ichild < startSlot + nOverlaps; ichild++) {
1030 computeAndAddChildrenR(ichild);
1031 totalNumberOfOverlaps++;
1032 }
1033 }
1034 }
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045 void computeOverlapTreeR(
1046 double[][] pos, double[] radii, double[] volumes, double[] gammas, boolean[] ishydrogen) {
1047 initOverlapTree(pos, radii, volumes, gammas, ishydrogen);
1048 for (int slot = 1; slot <= nAtoms; slot++) {
1049 computeAndAddChildrenR(slot);
1050 totalNumberOfOverlaps++;
1051 }
1052 }
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074 void computeVolumeUnderSlot2R(
1075 int slot,
1076 double[] psi1i,
1077 double[] f1i,
1078 double[] p1i,
1079 double[] psip1i,
1080 double[] fp1i,
1081 double[] pp1i,
1082 double[] energy1i,
1083 double[] fenergy1i,
1084 double[] penergy1i,
1085 int threadID,
1086 AtomicDoubleArray3D dr,
1087 AtomicDoubleArray dv,
1088 AtomicDoubleArray freeVolume,
1089 AtomicDoubleArray selfVolume) {
1090
1091 GaussianOverlap ov = overlaps.get(slot);
1092
1093
1094 if (ov.level >= maximumDepth) {
1095 maximumDepth = ov.level;
1096 }
1097
1098
1099
1100 double cf = ov.level % 2 == 0 ? -1.0 : 1.0;
1101
1102 double volcoeff = ov.level > 0 ? cf : 0;
1103
1104 double volcoeffp = ov.level > 0 ? volcoeff / (double) ov.level : 0;
1105
1106 int atom = ov.atom;
1107 double ai = overlaps.get(atom + 1).g.a;
1108 double a1i = ov.g.a;
1109 double a1 = a1i - ai;
1110
1111
1112 psi1i[0] = volcoeff * ov.volume;
1113 f1i[0] = volcoeff * ov.sfp;
1114
1115
1116 psip1i[0] = volcoeffp * ov.volume;
1117 fp1i[0] = volcoeffp * ov.sfp;
1118
1119
1120 energy1i[0] = volcoeffp * ov.gamma1i * ov.volume;
1121 fenergy1i[0] = volcoeffp * ov.sfp * ov.gamma1i;
1122
1123
1124 if (ov.childrenStartIndex >= 0) {
1125 for (int sloti = ov.childrenStartIndex;
1126 sloti < ov.childrenStartIndex + ov.childrenCount;
1127 sloti++) {
1128 double[] psi1it = new double[1];
1129 double[] f1it = new double[1];
1130 double[] p1it = new double[3];
1131 double[] psip1it = new double[1];
1132 double[] fp1it = new double[1];
1133 double[] pp1it = new double[3];
1134 double[] energy1it = new double[1];
1135 double[] fenergy1it = new double[1];
1136 double[] penergy1it = new double[3];
1137 computeVolumeUnderSlot2R(
1138 sloti,
1139 psi1it,
1140 f1it,
1141 p1it,
1142 psip1it,
1143 fp1it,
1144 pp1it,
1145 energy1it,
1146 fenergy1it,
1147 penergy1it,
1148 threadID,
1149 dr,
1150 dv,
1151 freeVolume,
1152 selfVolume);
1153 psi1i[0] += psi1it[0];
1154 f1i[0] += f1it[0];
1155 add(p1i, p1it, p1i);
1156
1157 psip1i[0] += psip1it[0];
1158 fp1i[0] += fp1it[0];
1159 add(pp1i, pp1it, pp1i);
1160
1161 energy1i[0] += energy1it[0];
1162 fenergy1i[0] += fenergy1it[0];
1163 add(penergy1i, penergy1it, penergy1i);
1164 }
1165 }
1166
1167
1168 if (ov.level > 0) {
1169
1170 freeVolume.add(threadID, atom, psi1i[0]);
1171 selfVolume.add(threadID, atom, psip1i[0]);
1172
1173
1174 double c2 = ai / a1i;
1175
1176
1177 double[] work1 = new double[3];
1178 double[] work2 = new double[3];
1179 double[] work3 = new double[3];
1180 scale(penergy1i, c2, work1);
1181 scale(ov.dv1, -fenergy1i[0], work2);
1182 add(work1, work2, work3);
1183 dr.add(threadID, atom, work3[0], work3[1], work3[2]);
1184
1185
1186 dv.add(threadID, atom, ov.g.v * fenergy1i[0]);
1187
1188
1189 c2 = a1 / a1i;
1190
1191
1192 scale(ov.dv1, f1i[0], work1);
1193 scale(pp1i, c2, work2);
1194 add(work1, work2, p1i);
1195
1196
1197 scale(ov.dv1, fp1i[0], work1);
1198 scale(pp1i, c2, work2);
1199 add(work1, work2, pp1i);
1200
1201
1202 scale(ov.dv1, fenergy1i[0], work1);
1203 scale(penergy1i, c2, work2);
1204 add(work1, work2, penergy1i);
1205
1206
1207 f1i[0] = ov.dvv1 * f1i[0];
1208 fp1i[0] = ov.dvv1 * fp1i[0];
1209 fenergy1i[0] = ov.dvv1 * fenergy1i[0];
1210 }
1211 }
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223 void computeVolume2R(
1224 SharedDouble volume,
1225 SharedDouble energy,
1226 AtomicDoubleArray3D dr,
1227 AtomicDoubleArray dv,
1228 AtomicDoubleArray freeVolume,
1229 AtomicDoubleArray selfvolume) {
1230
1231 double[] psi1i = new double[1];
1232 double[] f1i = new double[1];
1233 double[] p1i = new double[3];
1234 double[] psip1i = new double[1];
1235 double[] fp1i = new double[1];
1236 double[] pp1i = new double[3];
1237 double[] energy1i = new double[1];
1238 double[] fenergy1i = new double[1];
1239 double[] penergy1i = new double[3];
1240
1241 int threadID = 0;
1242
1243 computeVolumeUnderSlot2R(
1244 0,
1245 psi1i,
1246 f1i,
1247 p1i,
1248 psip1i,
1249 fp1i,
1250 pp1i,
1251 energy1i,
1252 fenergy1i,
1253 penergy1i,
1254 threadID,
1255 dr,
1256 dv,
1257 freeVolume,
1258 selfvolume);
1259
1260 volume.addAndGet(psi1i[0]);
1261 energy.addAndGet(energy1i[0]);
1262 }
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273 void rescanTreeV(
1274 double[][] pos, double[] radii, double[] volumes, double[] gammas, boolean[] ishydrogen) {
1275 initRescanTreeV(pos, radii, volumes, gammas, ishydrogen);
1276 rescanR(0);
1277 }
1278
1279
1280
1281
1282
1283
1284 void rescanR(int slot) {
1285 int parentIndex;
1286
1287
1288 GaussianOverlap ov = overlaps.get(slot);
1289
1290
1291 parentIndex = ov.parentIndex;
1292 if (parentIndex > 0) {
1293 GaussianVca g12 = new GaussianVca();
1294 double[] dVdr = new double[1];
1295 double[] dVdV = new double[1];
1296 double[] sfp = new double[1];
1297
1298 int atom = ov.atom;
1299 GaussianOverlap parent = overlaps.get(parentIndex);
1300 GaussianVca g1 = parent.g;
1301
1302
1303 GaussianVca g2 = overlaps.get(atom + 1).g;
1304 double gvol = overlapGaussianAlpha(g1, g2, g12, dVdr, dVdV, sfp);
1305 ov.g = g12;
1306 ov.volume = gvol;
1307
1308
1309
1310 sub(g2.c, g1.c, ov.dv1);
1311 scale(ov.dv1, -dVdr[0], ov.dv1);
1312
1313
1314 ov.dvv1 = dVdV[0];
1315 ov.sfp = sfp[0];
1316 ov.gamma1i = parent.gamma1i + overlaps.get(atom + 1).gamma1i;
1317 }
1318
1319
1320 for (int slotChild = ov.childrenStartIndex;
1321 slotChild < ov.childrenStartIndex + ov.childrenCount;
1322 slotChild++) {
1323 rescanR(slotChild);
1324 }
1325 }
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336 void initRescanTreeV(
1337 double[][] pos, double[] radii, double[] volumes, double[] gammas, boolean[] ishydrogen) {
1338 int slot = 0;
1339 GaussianOverlap ov = overlaps.get(slot);
1340 ov.level = 0;
1341 ov.volume = 0.0;
1342 ov.dv1 = new double[3];
1343 ov.dvv1 = 0.0;
1344 ov.selfVolume = 0.0;
1345 ov.sfp = 1.0;
1346 ov.gamma1i = 0.0;
1347
1348 slot = 1;
1349 for (int iat = 0; iat < nAtoms; iat++, slot++) {
1350 double a = sphereConversion / (radii[iat] * radii[iat]);
1351 double vol = ishydrogen[iat] ? 0.0 : volumes[iat];
1352 ov = overlaps.get(slot);
1353 ov.level = 1;
1354 ov.g.v = vol;
1355 ov.g.a = a;
1356 ov.g.c = pos[iat];
1357 ov.volume = vol;
1358 ov.dv1 = new double[3];
1359 ov.dvv1 = 1.0;
1360 ov.selfVolume = 0.0;
1361 ov.sfp = 1.0;
1362 ov.gamma1i = gammas[iat];
1363 }
1364 }
1365
1366
1367
1368
1369
1370
1371 void rescanGammaR(int slot) {
1372
1373
1374 GaussianOverlap go = overlaps.get(slot);
1375
1376
1377 int parentIndex = go.parentIndex;
1378 if (parentIndex > 0) {
1379 int atom = go.atom;
1380 GaussianOverlap parent = overlaps.get(parentIndex);
1381 go.gamma1i = parent.gamma1i + overlaps.get(atom + 1).gamma1i;
1382 }
1383
1384
1385 for (int slotChild = go.childrenStartIndex;
1386 slotChild < go.childrenStartIndex + go.childrenCount;
1387 slotChild++) {
1388 rescanGammaR(slotChild);
1389 }
1390 }
1391
1392
1393
1394
1395
1396
1397 void rescanTreeG(double[] gammas) {
1398
1399 int slot = 0;
1400 GaussianOverlap ov = overlaps.get(slot);
1401 ov.gamma1i = 0.;
1402
1403 slot = 1;
1404 for (int iat = 0; iat < nAtoms; iat++, slot++) {
1405 ov = overlaps.get(slot);
1406 ov.gamma1i = gammas[iat];
1407 }
1408
1409 rescanGammaR(0);
1410 }
1411
1412
1413
1414
1415 void printTree() {
1416
1417
1418 for (int i = 1; i <= nAtoms; i++) {
1419 printTreeR(i);
1420 }
1421 }
1422
1423
1424
1425
1426
1427
1428 void printTreeR(int slot) {
1429 GaussianOverlap ov = overlaps.get(slot);
1430 logger.info(format("tg: %d ", slot));
1431 ov.printOverlap();
1432 for (int i = ov.childrenStartIndex; i < ov.childrenStartIndex + ov.childrenCount; i++) {
1433 printTreeR(i);
1434 }
1435 }
1436
1437
1438
1439
1440
1441
1442
1443 int nChildrenUnderSlotR(int slot) {
1444 int n = 0;
1445 if (overlaps.get(slot).childrenCount > 0) {
1446 n += overlaps.get(slot).childrenCount;
1447
1448 for (int i = 0; i < overlaps.get(slot).childrenCount; i++) {
1449 n += nChildrenUnderSlotR(overlaps.get(slot).childrenStartIndex + i);
1450 }
1451 }
1452 return n;
1453 }
1454 }
1455
1456
1457
1458
1459 private class GaussVolRegion extends ParallelRegion {
1460
1461 private final InitGaussVol[] initGaussVol;
1462 private final GaussianOverlapTree[] localTree;
1463 private final ComputeTreeLoop[] computeTreeLoops;
1464 private final ComputeVolumeLoop[] computeVolumeLoops;
1465 private final RescanTreeLoop[] rescanTreeLoops;
1466 private final ReductionLoop[] reductionLoops;
1467 private GAUSSVOL_MODE mode = GAUSSVOL_MODE.COMPUTE_TREE;
1468 private double[][] coordinates = null;
1469
1470
1471
1472
1473
1474
1475 public GaussVolRegion(int nThreads) {
1476 initGaussVol = new InitGaussVol[nThreads];
1477 localTree = new GaussianOverlapTree[nThreads];
1478 computeTreeLoops = new ComputeTreeLoop[nThreads];
1479 computeVolumeLoops = new ComputeVolumeLoop[nThreads];
1480 rescanTreeLoops = new RescanTreeLoop[nThreads];
1481 reductionLoops = new ReductionLoop[nThreads];
1482 }
1483
1484 public void init(GAUSSVOL_MODE mode, double[][] coordinates) {
1485 this.mode = mode;
1486 this.coordinates = coordinates;
1487
1488
1489 totalVolume.set(0.0);
1490 energy.set(0.0);
1491 grad.reset(parallelTeam);
1492 gradV.reset(parallelTeam, 0, nAtoms - 1);
1493 freeVolume.reset(parallelTeam, 0, nAtoms - 1);
1494 selfVolume.reset(parallelTeam, 0, nAtoms - 1);
1495 }
1496
1497 @Override
1498 public void run() throws Exception {
1499 int threadIndex = getThreadIndex();
1500 if (computeTreeLoops[threadIndex] == null) {
1501 initGaussVol[threadIndex] = new InitGaussVol();
1502 localTree[threadIndex] = new GaussianOverlapTree(nAtoms);
1503 computeTreeLoops[threadIndex] = new ComputeTreeLoop();
1504 computeVolumeLoops[threadIndex] = new ComputeVolumeLoop();
1505 rescanTreeLoops[threadIndex] = new RescanTreeLoop();
1506 reductionLoops[threadIndex] = new ReductionLoop();
1507 }
1508 try {
1509 if (mode == GAUSSVOL_MODE.COMPUTE_TREE) {
1510 execute(0, nAtoms - 1, initGaussVol[threadIndex]);
1511 execute(1, nAtoms, computeTreeLoops[threadIndex]);
1512 execute(1, nAtoms, computeVolumeLoops[threadIndex]);
1513 execute(0, nAtoms - 1, reductionLoops[threadIndex]);
1514 } else {
1515 execute(1, nAtoms, rescanTreeLoops[threadIndex]);
1516 execute(1, nAtoms, computeVolumeLoops[threadIndex]);
1517 execute(0, nAtoms - 1, reductionLoops[threadIndex]);
1518 }
1519 } catch (RuntimeException ex) {
1520 logger.warning("Runtime exception computing tree in thread: " + threadIndex);
1521 throw ex;
1522 } catch (Exception e) {
1523 String message = "Fatal exception computing tree in thread: " + threadIndex + "\n";
1524 logger.log(Level.SEVERE, message, e);
1525 }
1526 }
1527
1528 private class InitGaussVol extends IntegerForLoop {
1529
1530 @Override
1531 public void run(int first, int last) throws Exception {
1532 for (int i = first; i <= last; i++) {
1533
1534 updateAtom(i);
1535 }
1536 }
1537 }
1538
1539
1540
1541
1542 private class ComputeTreeLoop extends IntegerForLoop {
1543
1544 @Override
1545 public void run(int first, int last) throws Exception {
1546
1547 int threadIndex = getThreadIndex();
1548 for (int slot = first; slot <= last; slot++) {
1549 localTree[threadIndex].computeAndAddChildrenR(slot);
1550 }
1551 }
1552
1553 @Override
1554 public void start() {
1555 localTree[getThreadIndex()].initOverlapTree(
1556 coordinates, radii, volumes, gammas, ishydrogen);
1557 }
1558 }
1559
1560
1561
1562
1563 private class ComputeVolumeLoop extends IntegerForLoop {
1564
1565 @Override
1566 public void run(int first, int last) throws Exception {
1567 int threadIndex = getThreadIndex();
1568 for (int slot = first; slot <= last; slot++) {
1569 double[] psi1i = new double[1];
1570 double[] f1i = new double[1];
1571 double[] p1i = new double[3];
1572 double[] psip1i = new double[1];
1573 double[] fp1i = new double[1];
1574 double[] pp1i = new double[3];
1575 double[] energy1i = new double[1];
1576 double[] fenergy1i = new double[1];
1577 double[] penergy1i = new double[3];
1578 localTree[threadIndex].computeVolumeUnderSlot2R(
1579 slot,
1580 psi1i,
1581 f1i,
1582 p1i,
1583 psip1i,
1584 fp1i,
1585 pp1i,
1586 energy1i,
1587 fenergy1i,
1588 penergy1i,
1589 threadIndex,
1590 grad,
1591 gradV,
1592 freeVolume,
1593 selfVolume);
1594 totalVolume.addAndGet(psi1i[0]);
1595 energy.addAndGet(energy1i[0]);
1596 }
1597 }
1598 }
1599
1600
1601
1602
1603 private class RescanTreeLoop extends IntegerForLoop {
1604
1605 @Override
1606 public void run(int first, int last) throws Exception {
1607
1608 int threadIndex = getThreadIndex();
1609 for (int slot = first; slot <= last; slot++) {
1610 localTree[threadIndex].rescanR(slot);
1611 }
1612 }
1613
1614 @Override
1615 public void start() {
1616 localTree[getThreadIndex()].initRescanTreeV(
1617 coordinates, radii, volumes, gammas, ishydrogen);
1618 }
1619 }
1620
1621
1622
1623
1624 private class ReductionLoop extends IntegerForLoop {
1625
1626 int threadID;
1627
1628 @Override
1629 public void run(int lb, int ub) {
1630 grad.reduce(lb, ub);
1631 gradV.reduce(lb, ub);
1632 freeVolume.reduce(lb, ub);
1633 selfVolume.reduce(lb, ub);
1634 for (int i = lb; i <= ub; i++) {
1635 if (volumes[i] > 0) {
1636 gradV.set(0, i, gradV.get(i) / volumes[i]);
1637 }
1638 }
1639 }
1640
1641 @Override
1642 public void start() {
1643 threadID = getThreadIndex();
1644 }
1645 }
1646 }
1647 }