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.utils; 39 40 import static ffx.numerics.math.DoubleMath.dist2; 41 import static java.lang.String.format; 42 import static java.lang.System.arraycopy; 43 import static org.apache.commons.math3.util.FastMath.max; 44 import static org.apache.commons.math3.util.FastMath.sqrt; 45 46 import com.github.quickhull3d.QuickHull3D; 47 import ffx.potential.bonded.Atom; 48 import ffx.utilities.Constants; 49 import java.util.Arrays; 50 import java.util.logging.Logger; 51 import java.util.stream.IntStream; 52 53 /** 54 * This ConvexHullOps class uses the QuickHull3D package by John E. Lloyd to construct and operate on 55 * 3D convex hulls: the minimal convex polyhedron that contains all points in a set of points. This 56 * is especially useful for max-dist operations, as the most distant points in a set are guaranteed 57 * to be part of the convex polyhedron. 58 * 59 * <p>The QuickHull3D package website is at quickhull3d.github.io/quickhull3d/ The algorithm it uses 60 * is described in Barber, Dobkin, and Huhdanpaa, "The Quickhull Algorithm for Convex Hulls" (ACM 61 * Transactions on Mathematical Software, Vol. 22, No. 4, December 1996), and the code is based on 62 * the C package qhull. 63 * 64 * @author Jacob M. Litman 65 * @author Michael J. Schnieders 66 * @since 1.0.0 67 */ 68 public class ConvexHullOps { 69 70 private static final Logger logger = Logger.getLogger(ConvexHullOps.class.getName()); 71 72 /** 73 * Constructs a convex hull from a set of atoms. 74 * 75 * @param atoms Atoms to build a convex hull for. 76 * @return A QuickHull3D implementation of convex hulls. 77 */ 78 public static QuickHull3D constructHull(Atom[] atoms) { 79 int nAts = atoms.length; 80 if (nAts < 4) { 81 throw new IllegalArgumentException( 82 format(" 3D convex hull ill-defined for less than 4 points, found %d", nAts)); 83 } 84 double[] xyz = new double[nAts * 3]; 85 for (int i = 0; i < nAts; i++) { 86 Atom at = atoms[i]; 87 int i3 = 3 * i; 88 xyz[i3] = at.getX(); 89 xyz[i3 + 1] = at.getY(); 90 xyz[i3 + 2] = at.getZ(); 91 } 92 return new QuickHull3D(xyz); 93 } 94 95 /** 96 * UNTESTED: Identifies atoms forming the convex hull. 97 * 98 * @param quickHull3D A QuickHull3D. 99 * @param allAtoms Atoms used in building the QuickHull3D. 100 * @return Atoms forming the convex hull. 101 */ 102 public static Atom[] identifyHullAtoms(QuickHull3D quickHull3D, Atom[] allAtoms) { 103 int[] indices = quickHull3D.getVertexPointIndices(); 104 return Arrays.stream(indices).mapToObj((int i) -> allAtoms[i]).toArray(Atom[]::new); 105 } 106 107 /** 108 * Find the maximum pairwise distance between vertex points on a convex hull. 109 * 110 * @param quickHull3D A QuickHull3D object. 111 * @return Maximum vertex-vertex distance. 112 */ 113 public static double maxDist(QuickHull3D quickHull3D) { 114 long time = -System.nanoTime(); 115 int nVerts = quickHull3D.getNumVertices(); 116 if (nVerts < 2) { 117 return 0; 118 } 119 double[] vertPoints = new double[3 * nVerts]; 120 quickHull3D.getVertices(vertPoints); 121 double maxDist = IntStream.range(0, nVerts).parallel().mapToDouble((int i) -> { 122 double[] xyz = new double[3]; 123 arraycopy(vertPoints, 3 * i, xyz, 0, 3); 124 double mij = 0; 125 for (int j = i + 1; j < nVerts; j++) { 126 double[] xyzJ = new double[3]; 127 arraycopy(vertPoints, 3 * j, xyzJ, 0, 3); 128 double distIJ = dist2(xyz, xyzJ); 129 mij = max(mij, distIJ); 130 } 131 return mij; 132 }).max().getAsDouble(); 133 maxDist = sqrt(maxDist); 134 time += System.nanoTime(); 135 if (time > 1E9) { 136 logger.warning(format(" Required %12.6g sec to find max distance on a convex hull." 137 + " It may be time to further optimize this!", Constants.NS2SEC * time)); 138 } 139 return maxDist; 140 } 141 142 /** 143 * Maximum pairwise distance between atoms in an array. Uses either the convex hull method (more 144 * than 10 atoms), or a brute-force loop (10 atoms or fewer). 145 * 146 * @param atoms Atoms to check max pairwise distance for. 147 * @return Max pairwise distance in Angstroms, or 0 (0 or 1 atoms given). 148 */ 149 public static double maxDist(Atom[] atoms) { 150 int nAts = atoms.length; 151 if (nAts < 2) { 152 return 0; 153 } else if (nAts > 10) { 154 return maxDist(constructHull(atoms)); 155 } else { 156 double maxDist = 0; 157 for (int i = 0; i < nAts - 1; i++) { 158 Atom atI = atoms[i]; 159 double[] xyzI = new double[3]; 160 xyzI = atI.getXYZ(xyzI); 161 for (int j = i + 1; j < nAts; j++) { 162 Atom atJ = atoms[j]; 163 double[] xyzJ = new double[3]; 164 xyzJ = atJ.getXYZ(xyzJ); 165 double dist = dist2(xyzI, xyzJ); 166 maxDist = max(dist, maxDist); 167 } 168 } 169 return sqrt(maxDist); 170 } 171 } 172 }