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 java.lang.String.format;
41 import static java.util.Arrays.fill;
42 import static java.util.Arrays.sort;
43 import static org.apache.commons.math3.util.FastMath.PI;
44 import static org.apache.commons.math3.util.FastMath.abs;
45 import static org.apache.commons.math3.util.FastMath.acos;
46 import static org.apache.commons.math3.util.FastMath.asin;
47 import static org.apache.commons.math3.util.FastMath.max;
48 import static org.apache.commons.math3.util.FastMath.min;
49 import static org.apache.commons.math3.util.FastMath.sqrt;
50
51 import edu.rit.pj.IntegerForLoop;
52 import edu.rit.pj.ParallelRegion;
53 import edu.rit.pj.reduction.SharedBooleanArray;
54 import edu.rit.pj.reduction.SharedDouble;
55 import ffx.numerics.atomic.AtomicDoubleArray3D;
56 import ffx.potential.bonded.Atom;
57 import ffx.potential.bonded.Residue;
58 import ffx.potential.parameters.VDWType;
59 import ffx.potential.utils.EnergyException;
60
61 import java.util.List;
62 import java.util.logging.Level;
63 import java.util.logging.Logger;
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82 public class SurfaceAreaRegion extends ParallelRegion {
83
84 private static final Logger logger = Logger.getLogger(SurfaceAreaRegion.class.getName());
85
86
87
88 private static final double delta = 1.0e-8;
89
90 private static final double delta2 = delta * delta;
91
92
93
94 private final Atom[] atoms;
95
96
97
98 private final int nAtoms;
99
100
101
102 private final int[][][] neighborLists;
103
104
105
106 private final double[] x, y, z;
107
108
109
110 private final boolean[] use;
111
112
113
114 private final InitLoop[] initLoop;
115
116
117
118 private final AtomOverlapLoop[] atomOverlapLoop;
119
120
121
122 private final SurfaceAreaLoop[] surfaceAreaLoop;
123
124
125
126 private final SharedDouble sharedSurfaceArea;
127
128
129
130 private final double probe;
131
132
133
134 private final int MAXARC = 1000;
135
136
137
138 private SharedBooleanArray skip;
139
140
141
142 private AtomicDoubleArray3D grad;
143
144
145
146 private int[][] overlaps;
147
148
149
150 private Integer[] overlapCounts;
151
152
153
154 private double[][] overlapDX;
155
156
157
158 private double[][] overlapDY;
159
160
161
162 private double[][] overlapDZ;
163
164
165
166 private double[][] overlapXY2;
167
168
169
170 private double[][] overlapR2;
171
172
173
174 private double[][] overlapR;
175
176
177
178 private IndexedDouble[][] gr;
179
180
181
182 private boolean[] buried;
183
184
185
186 private double[] area;
187
188
189
190 private double[] r;
191
192
193
194 private double surfaceTension;
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213 public SurfaceAreaRegion(
214 Atom[] atoms,
215 double[] x,
216 double[] y,
217 double[] z,
218 boolean[] use,
219 int[][][] neighborLists,
220 AtomicDoubleArray3D grad,
221 int nt,
222 double probe,
223 double surfaceTension) {
224 this.atoms = atoms;
225 this.nAtoms = atoms.length;
226 this.x = x;
227 this.y = y;
228 this.z = z;
229 this.use = use;
230 this.neighborLists = neighborLists;
231 this.probe = probe;
232 this.grad = grad;
233 this.surfaceTension = surfaceTension;
234
235 atomOverlapLoop = new AtomOverlapLoop[nt];
236 surfaceAreaLoop = new SurfaceAreaLoop[nt];
237 initLoop = new InitLoop[nt];
238 for (int i = 0; i < nt; i++) {
239 atomOverlapLoop[i] = new AtomOverlapLoop();
240 surfaceAreaLoop[i] = new SurfaceAreaLoop();
241 initLoop[i] = new InitLoop();
242 }
243 sharedSurfaceArea = new SharedDouble();
244 init();
245 }
246
247 @Override
248 public void finish() {
249 if (logger.isLoggable(Level.FINE)) {
250 int n = initLoop.length;
251 long initTime = 0;
252 long overlapTime = 0;
253 long cavTime = 0;
254 for (int i = 0; i < n; i++) {
255 initTime = max(initLoop[i].time, initTime);
256 overlapTime = max(atomOverlapLoop[i].time, overlapTime);
257 cavTime = max(surfaceAreaLoop[i].time, cavTime);
258 }
259 logger.fine(
260 format(
261 " Cavitation Init: %10.3f Overlap: %10.3f Cav: %10.3f",
262 initTime * 1e-9, overlapTime * 1e-9, cavTime * 1e-9));
263 }
264 }
265
266 public double getEnergy() {
267 return sharedSurfaceArea.get();
268 }
269
270 public double[] getArea() {
271 return area;
272 }
273
274 public double getResidueSurfaceArea(Residue residue) {
275 List<Atom> atoms = residue.getAtomList();
276 double residueArea = 0.0;
277 for (Atom atom : atoms) {
278 int i = atom.getXyzIndex() - 1;
279 residueArea += area[i];
280 }
281 return residueArea;
282 }
283
284 public final void init() {
285 if (overlapCounts == null || overlapCounts.length < nAtoms) {
286 overlapCounts = new Integer[nAtoms];
287 overlapDX = new double[nAtoms][MAXARC];
288 overlapDY = new double[nAtoms][MAXARC];
289 overlapDZ = new double[nAtoms][MAXARC];
290 overlapXY2 = new double[nAtoms][MAXARC];
291 overlapR2 = new double[nAtoms][MAXARC];
292 overlapR = new double[nAtoms][MAXARC];
293 gr = new IndexedDouble[nAtoms][MAXARC];
294 overlaps = new int[nAtoms][MAXARC];
295 buried = new boolean[nAtoms];
296 skip = new SharedBooleanArray(nAtoms);
297 area = new double[nAtoms];
298 r = new double[nAtoms];
299 }
300
301
302 for (int i = 0; i < nAtoms; i++) {
303 VDWType type = atoms[i].getVDWType();
304 double rmini = type.radius;
305 r[i] = rmini / 2.0;
306 if (r[i] != 0.0) {
307 r[i] = r[i] + probe;
308 }
309 skip.set(i, true);
310 }
311 }
312
313 @Override
314 public void run() {
315 try {
316 execute(0, nAtoms - 1, initLoop[getThreadIndex()]);
317 execute(0, nAtoms - 1, atomOverlapLoop[getThreadIndex()]);
318 execute(0, nAtoms - 1, surfaceAreaLoop[getThreadIndex()]);
319 } catch (Exception e) {
320 String message =
321 "Fatal exception computing Cavitation energy in thread " + getThreadIndex() + "\n";
322 logger.log(Level.SEVERE, message, e);
323 }
324 }
325
326 public void setSurfaceTension(double surfaceTension) {
327 this.surfaceTension = surfaceTension;
328 }
329
330 @Override
331 public void start() {
332 sharedSurfaceArea.set(0.0);
333 }
334
335 private static class IndexedDouble implements Comparable<IndexedDouble> {
336
337 public double value;
338 public int key;
339
340 IndexedDouble(double value, int key) {
341 this.value = value;
342 this.key = key;
343 }
344
345 @Override
346 public int compareTo(IndexedDouble d) {
347 return Double.compare(value, d.value);
348 }
349 }
350
351
352
353
354 private class InitLoop extends IntegerForLoop {
355
356 public long time;
357
358 @Override
359 public void finish() {
360 time += System.nanoTime();
361 }
362
363 @Override
364 public void run(int lb, int ub) throws Exception {
365
366
367 for (int i = lb; i <= ub; i++) {
368 buried[i] = false;
369 area[i] = 0.0;
370 overlapCounts[i] = 0;
371 double xr = x[i];
372 double yr = y[i];
373 double zr = z[i];
374 double rri = r[i];
375 final int[] list = neighborLists[0][i];
376 for (int k : list) {
377 double rrik = rri + r[k];
378 double dx = x[k] - xr;
379 double dy = y[k] - yr;
380 double dz = z[k] - zr;
381 double ccsq = dx * dx + dy * dy + dz * dz;
382 if (ccsq <= rrik * rrik) {
383 if (use[i] || use[k]) {
384 skip.set(k, false);
385 skip.set(i, false);
386 }
387 }
388 }
389 }
390 }
391
392 @Override
393 public void start() {
394 time = -System.nanoTime();
395 }
396 }
397
398
399
400
401
402
403 private class AtomOverlapLoop extends IntegerForLoop {
404
405 public long time;
406
407 @Override
408 public void finish() {
409 time += System.nanoTime();
410 }
411
412 @Override
413 public void run(int lb, int ub) {
414
415
416 for (int i = lb; i <= ub; i++) {
417 if (skip.get(i) || !use[i]) {
418 continue;
419 }
420 int[] list = neighborLists[0][i];
421 for (int k : list) {
422 if (k == i) {
423 continue;
424 }
425 pair(i, k);
426 if (!skip.get(k) || !use[k]) {
427 pair(k, i);
428 }
429 }
430 }
431 }
432
433 @Override
434 public void start() {
435 time = -System.nanoTime();
436 }
437
438
439
440
441
442
443
444 private void pair(int i, int k) {
445 double xi = x[i];
446 double yi = y[i];
447 double zi = z[i];
448 double rri = r[i];
449 double rplus = rri + r[k];
450 double dx = x[k] - xi;
451 double dy = y[k] - yi;
452 double dz = z[k] - zi;
453 if (abs(dx) >= rplus || abs(dy) >= rplus || abs(dz) >= rplus) {
454 return;
455 }
456
457
458
459 double xysq = dx * dx + dy * dy;
460 if (xysq < delta2) {
461 dx = delta;
462 dy = 0.0;
463 xysq = delta2;
464 }
465
466 double r2 = xysq + dz * dz;
467 double dr = sqrt(r2);
468 if (rplus - dr <= delta) {
469 return;
470 }
471 double rminus = rri - r[k];
472
473 synchronized (overlaps[i]) {
474
475 if (dr - abs(rminus) <= delta) {
476 if (rminus <= 0.0) {
477
478 buried[i] = true;
479 }
480 return;
481 }
482
483 int n = overlapCounts[i];
484 overlaps[i][n] = k;
485 overlapDX[i][n] = dx;
486 overlapDY[i][n] = dy;
487 overlapDZ[i][n] = dz;
488 overlapXY2[i][n] = xysq;
489 overlapR2[i][n] = r2;
490 overlapR[i][n] = dr;
491 double rri2 = 2.0 * rri;
492 gr[i][n] = new IndexedDouble((r2 + rplus * rminus) / (rri2 * overlapR[i][n]), n);
493 overlapCounts[i]++;
494 if (overlapCounts[i] >= MAXARC) {
495 throw new EnergyException(format(" Increase the value of MAXARC to (%d).", overlapCounts[i]));
496 }
497 }
498 }
499 }
500
501
502
503
504
505
506 private class SurfaceAreaLoop extends IntegerForLoop {
507
508
509 private static final double pix2 = 2.0 * PI;
510 private static final double pix4 = 4.0 * PI;
511 private static final double pid2 = PI / 2.0;
512 private static final double eps = 1.0e-8;
513 public long time;
514 private double localSurfaceEnergy;
515 private IndexedDouble[] arci;
516 private boolean[] omit;
517 private double[] xc;
518 private double[] yc;
519 private double[] zc;
520 private double[] dsq;
521 private double[] b;
522 private double[] bsq;
523 private double[] bg;
524 private double[] risq;
525 private double[] ri;
526 private double[] ther;
527 private double[] ider;
528 private double[] sign_yder;
529 private double[] arcf;
530 private double[] ex;
531 private double[] ux;
532 private double[] uy;
533 private double[] uz;
534 private int[] kent;
535 private int[] kout;
536 private int[] intag;
537 private int[] lt;
538 private int i;
539 private int j;
540 private int ib;
541 private int threadID;
542
543 SurfaceAreaLoop() {
544 allocateMemory(MAXARC);
545 }
546
547 @Override
548 public void finish() {
549 sharedSurfaceArea.addAndGet(localSurfaceEnergy);
550 time += System.nanoTime();
551 }
552
553 @Override
554 public void run(int lb, int ub) {
555
556 for (int ir = lb; ir <= ub; ir++) {
557 if (skip.get(ir) || !use[ir] || buried[i]) {
558 continue;
559 }
560 double rri = r[ir];
561 double rri2 = 2.0 * rri;
562 double rrisq = rri * rri;
563 surface(rri, rri2, rrisq, surfaceTension, false, ir);
564 if (area[ir] < 0.0) {
565 logger.log(Level.WARNING, format(" Negative surface area set to 0 for atom %d.", ir));
566 area[ir] = 0.0;
567 }
568 area[ir] *= rrisq * surfaceTension;
569 localSurfaceEnergy += area[ir];
570 }
571 }
572
573 @Override
574 public void start() {
575 time = -System.nanoTime();
576 threadID = getThreadIndex();
577 localSurfaceEnergy = 0.0;
578 fill(ider, 0);
579 fill(sign_yder, 0);
580 }
581
582
583
584
585
586
587
588
589
590
591
592 public void surface(double rri, double rri2, double rrisq, double wght, boolean moved, int ir) {
593
594 ib = 0;
595 int jb = 0;
596 double arcLength = 0.0;
597 double exang = 0.0;
598
599
600 if (overlapCounts[ir] == 0) {
601 area[ir] = pix4;
602 return;
603 }
604
605 if (overlapCounts[ir] == 1) {
606 double dx = overlapDX[ir][0];
607 double dy = overlapDY[ir][0];
608 double dz = overlapDZ[ir][0];
609 double r2 = overlapR2[ir][0];
610 double rr = overlapR[ir][0];
611 double arcsum = pix2;
612 ib = ib + 1;
613 arcLength += gr[ir][0].value * arcsum;
614 if (!moved) {
615 int k = overlaps[ir][0];
616 double t1 = arcsum * rrisq * (r2 - rrisq + r[k] * r[k]) / (rri2 * r2 * rr);
617 double gx = dx * t1 * wght;
618 double gy = dy * t1 * wght;
619 double gz = dz * t1 * wght;
620 grad.sub(threadID, ir, gx, gy, gz);
621 grad.add(threadID, k, gx, gy, gz);
622 }
623 area[ir] = ib * pix2 + exang + arcLength;
624 area[ir] = area[ir] % pix4;
625 return;
626 }
627
628
629
630
631
632 sort(gr[ir], 0, overlapCounts[ir]);
633 for (int j = 0; j < overlapCounts[ir]; j++) {
634 int k = gr[ir][j].key;
635 intag[j] = overlaps[ir][k];
636 xc[j] = overlapDX[ir][k];
637 yc[j] = overlapDY[ir][k];
638 zc[j] = overlapDZ[ir][k];
639 dsq[j] = overlapXY2[ir][k];
640 b[j] = overlapR[ir][k];
641 bsq[j] = overlapR2[ir][k];
642 omit[j] = false;
643 }
644
645
646 for (int i = 0; i < overlapCounts[ir]; i++) {
647 double gi = gr[ir][i].value * rri;
648 bg[i] = b[i] * gi;
649 risq[i] = rrisq - gi * gi;
650 ri[i] = sqrt(risq[i]);
651 ther[i] = pid2 - asin(min(1.0, max(-1.0, gr[ir][i].value)));
652 }
653
654
655 for (int k = 0; k < overlapCounts[ir] - 1; k++) {
656 if (omit[k]) {
657 continue;
658 }
659 double txk = xc[k];
660 double tyk = yc[k];
661 double tzk = zc[k];
662 double bk = b[k];
663 double therk = ther[k];
664 for (j = k + 1; j < overlapCounts[ir]; j++) {
665 if (omit[j]) {
666 continue;
667 }
668
669
670
671
672 double cc = (txk * xc[j] + tyk * yc[j] + tzk * zc[j]) / (bk * b[j]);
673
674
675 cc = acos(min(1.0, max(-1.0, cc)));
676 double td = therk + ther[j];
677
678 if (cc >= td) {
679 continue;
680 }
681
682 if (cc + ther[j] < therk) {
683 omit[j] = true;
684 continue;
685 }
686
687 if (cc > delta) {
688 if (pix2 - cc <= td) {
689 area[ir] = 0.0;
690 return;
691 }
692 }
693 }
694 }
695
696
697 for (int k = 0; k < overlapCounts[ir]; k++) {
698 if (omit[k]) {
699 continue;
700 }
701 boolean komit = omit[k];
702 omit[k] = true;
703 int narc = 0;
704 boolean top = false;
705 double txk = xc[k];
706 double tyk = yc[k];
707 double tzk = zc[k];
708 double dk = sqrt(dsq[k]);
709 double bsqk = bsq[k];
710 double bk = b[k];
711 double gk = gr[ir][k].value * rri;
712 double risqk = risq[k];
713 double rik = ri[k];
714 double therk = ther[k];
715
716
717 double t1 = tzk / (bk * dk);
718 double axx = txk * t1;
719 double axy = tyk * t1;
720 double axz = dk / bk;
721 double ayx = tyk / dk;
722 double ayy = txk / dk;
723 double azx = txk / bk;
724 double azy = tyk / bk;
725 double azz = tzk / bk;
726 for (int l = 0; l < overlapCounts[ir]; l++) {
727 if (omit[l]) {
728 continue;
729 }
730 double txl = xc[l];
731 double tyl = yc[l];
732 double tzl = zc[l];
733
734
735 double uxl = txl * axx + tyl * axy - tzl * axz;
736 double uyl = tyl * ayy - txl * ayx;
737 double uzl = txl * azx + tyl * azy + tzl * azz;
738 double cosine = min(1.0, max(-1.0, uzl / b[l]));
739 if (acos(cosine) < therk + ther[l]) {
740 double dsql = uxl * uxl + uyl * uyl;
741 double tb = uzl * gk - bg[l];
742 double txb = uxl * tb;
743 double tyb = uyl * tb;
744 double td = rik * dsql;
745 double tr2 = risqk * dsql - tb * tb;
746 tr2 = max(eps, tr2);
747 double tr = sqrt(tr2);
748 double txr = uxl * tr;
749 double tyr = uyl * tr;
750
751
752 tb = (txb + tyr) / td;
753 tb = min(1.0, max(-1.0, tb));
754 double tk1 = acos(tb);
755 if (tyb - txr < 0.0) {
756 tk1 = pix2 - tk1;
757 }
758 tb = (txb - tyr) / td;
759 tb = min(1.0, max(-1.0, tb));
760 double tk2 = acos(tb);
761 if (tyb + txr < 0.0) {
762 tk2 = pix2 - tk2;
763 }
764 double thec = (rrisq * uzl - gk * bg[l]) / (rik * ri[l] * b[l]);
765 double the = 0.0;
766 if (abs(thec) < 1.0) {
767 the = -acos(thec);
768 } else if (thec >= 1.0) {
769 the = 0.0;
770 } else if (thec <= -1.0) {
771 the = -PI;
772 }
773
774
775
776
777 cosine = min(1.0, max(-1.0, (uzl * gk - uxl * rik) / (b[l] * rri)));
778 double ti, tf;
779 if ((acos(cosine) - ther[l]) * (tk2 - tk1) <= 0.0) {
780 ti = tk2;
781 tf = tk1;
782 } else {
783 ti = tk1;
784 tf = tk2;
785 }
786 narc += 1;
787 if (narc > MAXARC) {
788 throw new EnergyException(format(" Increase value of MAXARC %d.", narc));
789 }
790 int narc1 = narc - 1;
791 if (tf <= ti) {
792 arcf[narc1] = tf;
793 arci[narc1] = new IndexedDouble(0.0, narc1);
794 tf = pix2;
795 lt[narc1] = l;
796 ex[narc1] = the;
797 top = true;
798 narc = narc + 1;
799 narc1 = narc - 1;
800 }
801 arcf[narc1] = tf;
802 arci[narc1] = new IndexedDouble(ti, narc1);
803 lt[narc1] = l;
804 ex[narc1] = the;
805 ux[l] = uxl;
806 uy[l] = uyl;
807 uz[l] = uzl;
808 }
809 }
810 omit[k] = komit;
811
812
813 if (narc <= 0) {
814 double arcsum = pix2;
815 ib += 1;
816 arcLength += gr[ir][k].value * arcsum;
817 if (!moved) {
818 int in = intag[k];
819 t1 = arcsum * rrisq * (bsqk - rrisq + r[in] * r[in]) / (rri2 * bsqk * bk);
820 grad.sub(threadID, ir, txk * t1 * wght, tyk * t1 * wght, tzk * t1 * wght);
821 grad.add(threadID, in, txk * t1 * wght, tyk * t1 * wght, tzk * t1 * wght);
822 }
823 continue;
824 }
825
826
827 sort(arci, 0, narc);
828 double arcsum = arci[0].value;
829 int mi = arci[0].key;
830 double t = arcf[mi];
831 int ni = mi;
832 for (j = 1; j < narc; j++) {
833 int m = arci[j].key;
834 if (t < arci[j].value) {
835 arcsum += (arci[j].value - t);
836 exang += ex[ni];
837 jb += 1;
838 if (jb >= MAXARC) {
839 throw new EnergyException(format("Increase the value of MAXARC (%d).", jb));
840 }
841 int l = lt[ni];
842 ider[l] += 1;
843 sign_yder[l] += 1;
844 kent[jb] = MAXARC * (l + 1) + (k + 1);
845 l = lt[m];
846 ider[l] += 1;
847 sign_yder[l] -= 1;
848 kout[jb] = MAXARC * (k + 1) + (l + 1);
849 }
850 double tt = arcf[m];
851 if (tt >= t) {
852 t = tt;
853 ni = m;
854 }
855 }
856 arcsum += (pix2 - t);
857 if (!top) {
858 exang += ex[ni];
859 jb = jb + 1;
860 int l = lt[ni];
861 ider[l] += 1;
862 sign_yder[l] += 1;
863 kent[jb] = MAXARC * (l + 1) + (k + 1);
864 l = lt[mi];
865 ider[l] += 1;
866 sign_yder[l] -= 1;
867 kout[jb] = MAXARC * (k + 1) + (l + 1);
868 }
869
870
871 for (int l = 0; l <= overlapCounts[ir]; l++) {
872 if (ider[l] == 0) {
873 continue;
874 }
875 double rcn = ider[l] * rrisq;
876 ider[l] = 0;
877 double uzl = uz[l];
878 double gl = gr[ir][l].value * rri;
879 double bgl = bg[l];
880 double bsql = bsq[l];
881 double risql = risq[l];
882 double wxlsq = bsql - uzl * uzl;
883 double wxl = sqrt(wxlsq);
884 double p = bgl - gk * uzl;
885 double v = risqk * wxlsq - p * p;
886 v = max(eps, v);
887 v = sqrt(v);
888 t1 = rri * (gk * (bgl - bsql) + uzl * (bgl - rrisq)) / (v * risql * bsql);
889 double deal = -wxl * t1;
890 double decl = -uzl * t1 - rri / v;
891 double dtkal = (wxlsq - p) / (wxl * v);
892 double dtkcl = (uzl - gk) / v;
893 double s = gk * b[l] - gl * uzl;
894 t1 = 2.0 * gk - uzl;
895 double t2 = rrisq - bgl;
896 double dtlal =
897 -(risql * wxlsq * b[l] * t1 - s * (wxlsq * t2 + risql * bsql))
898 / (risql * wxl * bsql * v);
899 double dtlcl = -(risql * b[l] * (uzl * t1 - bgl) - uzl * t2 * s) / (risql * bsql * v);
900 double gaca = rcn * (deal - (gk * dtkal - gl * dtlal) / rri) / wxl;
901 double gacb = (gk - uzl * gl / b[l]) * sign_yder[l] * rri / wxlsq;
902 sign_yder[l] = 0;
903 if (!moved) {
904 double faca = ux[l] * gaca - uy[l] * gacb;
905 double facb = uy[l] * gaca + ux[l] * gacb;
906 double facc = rcn * (decl - (gk * dtkcl - gl * dtlcl) / rri);
907 double dax = axx * faca - ayx * facb + azx * facc;
908 double day = axy * faca + ayy * facb + azy * facc;
909 double daz = azz * facc - axz * faca;
910 int in = intag[l];
911 grad.add(threadID, ir, dax * wght, day * wght, daz * wght);
912 grad.sub(threadID, in, dax * wght, day * wght, daz * wght);
913 }
914 }
915 arcLength += gr[ir][k].value * arcsum;
916 if (!moved) {
917 int in = intag[k];
918 t1 = arcsum * rrisq * (bsqk - rrisq + r[in] * r[in]) / (rri2 * bsqk * bk);
919 grad.sub(threadID, ir, txk * t1 * wght, tyk * t1 * wght, tzk * t1 * wght);
920 grad.add(threadID, in, txk * t1 * wght, tyk * t1 * wght, tzk * t1 * wght);
921 }
922 }
923 if (arcLength == 0.0) {
924 area[ir] = 0.0;
925 return;
926 }
927 if (jb == 0) {
928 area[ir] = ib * pix2 + exang + arcLength;
929 area[ir] = area[ir] % pix4;
930 return;
931 }
932
933
934 j = 0;
935 for (int k = 1; k <= jb; k++) {
936 if (kout[k] == -1) {
937 continue;
938 }
939 i = k;
940 boolean success = independentBoundaries(k, exang, jb, ir, arcLength);
941 if (success) {
942 return;
943 }
944 }
945 ib = ib + 1;
946 logger.log(Level.WARNING, format(" Connectivity error at atom %d", ir));
947 area[ir] = 0.0;
948 }
949
950 private void allocateMemory(int maxarc) {
951 arci = new IndexedDouble[maxarc];
952 arcf = new double[maxarc];
953 risq = new double[maxarc];
954 ri = new double[maxarc];
955 dsq = new double[maxarc];
956 bsq = new double[maxarc];
957 intag = new int[maxarc];
958 lt = new int[maxarc];
959 kent = new int[maxarc];
960 kout = new int[maxarc];
961 ider = new double[maxarc];
962 sign_yder = new double[maxarc];
963 xc = new double[maxarc];
964 yc = new double[maxarc];
965 zc = new double[maxarc];
966 b = new double[maxarc];
967 omit = new boolean[maxarc];
968 bg = new double[maxarc];
969 ther = new double[maxarc];
970 ex = new double[maxarc];
971 ux = new double[maxarc];
972 uy = new double[maxarc];
973 uz = new double[maxarc];
974 }
975
976
977
978
979
980
981
982
983
984
985 boolean independentBoundaries(int k, double exAngle, int jb, int ir, double arcLength) {
986 int m = kout[i];
987 kout[i] = -1;
988 j = j + 1;
989 for (int ii = 1; ii <= jb; ii++) {
990 if (m == kent[ii]) {
991 if (ii == k) {
992 ib++;
993 if (j == jb) {
994 area[ir] = ib * 2.0 * PI + exAngle + arcLength;
995 area[ir] = area[ir] % (4.0 * PI);
996 return true;
997 }
998 return false;
999 }
1000 i = ii;
1001 return independentBoundaries(k, exAngle, jb, ir, arcLength);
1002 }
1003 }
1004 return false;
1005 }
1006 }
1007 }