View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
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   * SurfaceAreaRegion performs an analytical computation of the weighted solvent accessible surface
67   * area of each atom and the first derivatives of the area with respect to Cartesian coordinates
68   *
69   * <p>Literature references:
70   *
71   * <p>T. J. Richmond, "Solvent Accessible Surface Area and Excluded Volume in Proteins", Journal of
72   * Molecular Biology, 178, 63-89 (1984)
73   *
74   * <p>L. Wesson and D. Eisenberg, "Atomic Solvation Parameters Applied to Molecular Dynamics of
75   * Proteins in Solution", Protein Science, 1, 227-235 (1992)
76   *
77   * <p>This was ported from the TINKER "surface.f" code.
78   *
79   * @author Michael J. Schnieders
80   * @since 1.0
81   */
82  public class SurfaceAreaRegion extends ParallelRegion {
83  
84    private static final Logger logger = Logger.getLogger(SurfaceAreaRegion.class.getName());
85    /**
86     * Tolerance used in the tests for sphere overlaps and for colinearity.
87     */
88    private static final double delta = 1.0e-8;
89  
90    private static final double delta2 = delta * delta;
91    /**
92     * Array of atoms.
93     */
94    private final Atom[] atoms;
95    /**
96     * Number of atoms.
97     */
98    private final int nAtoms;
99    /**
100    * Atomic neighbor lists.
101    */
102   private final int[][][] neighborLists;
103   /**
104    * Atomic coordinates.
105    */
106   private final double[] x, y, z;
107   /**
108    * Per-atom flag to indicate the atom is used.
109    */
110   private final boolean[] use;
111   /**
112    * Initialization loops.
113    */
114   private final InitLoop[] initLoop;
115   /**
116    * Atom overlap loops.
117    */
118   private final AtomOverlapLoop[] atomOverlapLoop;
119   /**
120    * Surface area loops.
121    */
122   private final SurfaceAreaLoop[] surfaceAreaLoop;
123   /**
124    * Total surface area.
125    */
126   private final SharedDouble sharedSurfaceArea;
127   /**
128    * Radius of the probe sphere.
129    */
130   private final double probe;
131   /**
132    * Maximum number of arcs.
133    */
134   private final int MAXARC = 1000;
135   /**
136    * Atoms to skip in the area calculation.
137    */
138   private SharedBooleanArray skip;
139   /**
140    * Atomic gradient array.
141    */
142   private AtomicDoubleArray3D grad;
143   /**
144    * Atom i overlaps. intag1[atom index][overlap index].
145    */
146   private int[][] overlaps;
147   /**
148    * Number of overlaps for each atom.
149    */
150   private Integer[] overlapCounts;
151   /**
152    * X-separation for atom i with overlap k.
153    */
154   private double[][] overlapDX;
155   /**
156    * Y-separation for atom i with overlap k.
157    */
158   private double[][] overlapDY;
159   /**
160    * Z-separation for atom i with overlap k.
161    */
162   private double[][] overlapDZ;
163   /**
164    * DX*DX + DY*DY for atom i with overlap k.
165    */
166   private double[][] overlapXY2;
167   /**
168    * R^2 for atom i with overlap k.
169    */
170   private double[][] overlapR2;
171   /**
172    * R for atom i with overlap k.
173    */
174   private double[][] overlapR;
175   /**
176    * Degree of overlap for atom i with overlap k.
177    */
178   private IndexedDouble[][] gr;
179   /**
180    * Per atom flag to indicate if the atom is buried and has no surface area.
181    */
182   private boolean[] buried;
183   /**
184    * Accessible surface area of each atom.
185    */
186   private double[] area;
187   /**
188    * Atomic radii.
189    */
190   private double[] r;
191   /**
192    * Weight assigned to each atom's area; if set to 1.0, return is actual area in square Angstroms
193    */
194   private double surfaceTension;
195 
196   /**
197    * This class is a port of the Cavitation code in TINKER.
198    *
199    * <p>GK implicit solvents are moving to use Gaussian based definitions of surface area for
200    * efficiency.
201    *
202    * @param atoms          Atom array.
203    * @param x              X-coordinate array.
204    * @param y              Y-coordinate array.
205    * @param z              Z-coordinate array.
206    * @param use            Specifies if the atom is used in the potential.
207    * @param neighborLists  Neighbor-list array.
208    * @param grad           Gradient array.
209    * @param nt             Number of threads.
210    * @param probe          Solvent probe radius.
211    * @param surfaceTension Surface tension.
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     // Set the sphere radii.
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    * Initialize arrays for Cavitation calculation.
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       // Set the "skip" array to exclude all inactive atoms
366       // that do not overlap any of the current active atoms
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    * Compute Cavitation energy for a range of atoms.
400    *
401    * @since 1.0
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       // Find overlaps with the current sphere.
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      * Find overlaps of atom i by atom k.
440      *
441      * @param i Index of atom i.
442      * @param k Index of atom k.
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       // Check for overlap of spheres by testing center to center
458       // distance against sum and difference of radii.
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         // Check for a completely buried "ir" sphere.
475         if (dr - abs(rminus) <= delta) {
476           if (rminus <= 0.0) {
477             // SA for this atom is zero.
478             buried[i] = true;
479           }
480           return;
481         }
482         // Calculate overlap parameters between "i" and "ir" sphere.
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    * Compute Cavitation energy for a range of atoms.
503    *
504    * @since 1.0
505    */
506   private class SurfaceAreaLoop extends IntegerForLoop {
507 
508     // Set pi multiples, overlap criterion and tolerances.
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       // Compute the area and derivatives of current "ir" sphere
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      * Calculate surface area.
584      *
585      * @param rri   Radius.
586      * @param rri2  Diameter.
587      * @param rrisq Radius squared.
588      * @param wght  Surface tension.
589      * @param moved Atom has been moved.
590      * @param ir    Atom index.
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       // Case where no other spheres overlap the current sphere.
600       if (overlapCounts[ir] == 0) {
601         area[ir] = pix4;
602         return;
603       }
604       // Case where only one sphere overlaps the current sphere.
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        General case where more than one sphere intersects the
629        current sphere; sort intersecting spheres by their degree of
630        overlap with the current main sphere
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       // Radius of the each circle on the surface of the "ir" sphere.
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       // Find boundary of inaccessible area on "ir" sphere.
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            Check to see if J circle is intersecting K circle;
670            get distance between circle centers and sum of radii.
671           */
672           double cc = (txk * xc[j] + tyk * yc[j] + tzk * zc[j]) / (bk * b[j]);
673 
674           // Check acos FORTRAN vs. Java.
675           cc = acos(min(1.0, max(-1.0, cc)));
676           double td = therk + ther[j];
677           // Check to see if circles enclose separate regions
678           if (cc >= td) {
679             continue;
680           }
681           // Check for circle J completely inside circle K
682           if (cc + ther[j] < therk) {
683             omit[j] = true;
684             continue;
685           }
686           // Check for circles essentially parallel.
687           if (cc > delta) {
688             if (pix2 - cc <= td) {
689               area[ir] = 0.0;
690               return;
691             }
692           }
693         }
694       }
695 
696       // Find T value of circle intersections.
697       for (int k = 0; k < overlapCounts[ir]; k++) {
698         if (omit[k]) {
699           continue; // goto 110
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         // Rotation matrix elements.
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           // Rotate spheres so K vector collinear with z-axis.
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             // Get T values of intersection for K circle.
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              See if "tk1" is entry or exit point; check t=0
775              point; "ti" is exit point, "tf" is entry point.
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         // Special case; K circle without intersections.
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         // General case; sum up arclength and set connectivity code.
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         // Calculate the surface area derivatives.
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       // Find number of independent boundaries and check connectivity.
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      * Find number of independent boundaries and check connectivity.
978      *
979      * @param k         Atom index.
980      * @param exAngle   Ex angle.
981      * @param jb        Upper limit.
982      * @param ir        Atom index.
983      * @param arcLength Arc length.
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 }