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.numerics.math; 39 40 import javax.annotation.Nullable; 41 42 import static java.lang.String.format; 43 import static org.apache.commons.math3.util.FastMath.abs; 44 import static org.apache.commons.math3.util.FastMath.acos; 45 import static org.apache.commons.math3.util.FastMath.max; 46 import static org.apache.commons.math3.util.FastMath.min; 47 import static org.apache.commons.math3.util.FastMath.sqrt; 48 49 import java.util.logging.Logger; 50 51 /** 52 * The FloatMath class is a simple math library that operates on 3-coordinate float arrays. 53 * 54 * <p>All methods are static and thread-safe. 55 * 56 * <p>Use instances of Float3 for convenience. 57 * 58 * @author Michael J. Schnieders 59 * @since 1.0 60 */ 61 public class FloatMath { 62 63 private static final Logger logger = Logger.getLogger(FloatMath.class.getName()); 64 65 /** 66 * Compute a * b + c and return the result in a new array. 67 * 68 * @param a First vector. 69 * @param b Scalar. 70 * @param c Second vector. 71 * @return Returns a * b + c. 72 */ 73 public static float[] fma(float[] a, float b, float[] c) { 74 return fma(a, b, c, new float[3]); 75 } 76 77 /** 78 * Compute a * b + c. 79 * 80 * @param a First vector. 81 * @param b Scalar. 82 * @param c Second vector. 83 * @param ret Result vector. 84 * @return Returns a * b + c in the vector ret. 85 */ 86 public static float[] fma(float[] a, float b, float[] c, float[] ret) { 87 ret[0] = Math.fma(a[0], b, c[0]); 88 ret[1] = Math.fma(a[1], b, c[1]); 89 ret[2] = Math.fma(a[2], b, c[2]); 90 return ret; 91 } 92 93 /** 94 * Finds the cross-product between two vectors. 95 * 96 * @param a First vector. 97 * @param b Second vector. 98 * @return Returns the cross-product. 99 */ 100 public static float[] X(float[] a, float[] b) { 101 return X(a, b, new float[3]); 102 } 103 104 /** 105 * Finds the cross-product between two vectors. 106 * 107 * @param a First vector. 108 * @param b Second vector. 109 * @param ret The cross-product of a x b. 110 * @return Returns the cross-product ret. 111 */ 112 public static float[] X(float[] a, float[] b, float[] ret) { 113 ret[0] = a[1] * b[2] - a[2] * b[1]; 114 ret[1] = a[2] * b[0] - a[0] * b[2]; 115 ret[2] = a[0] * b[1] - a[1] * b[0]; 116 return ret; 117 } 118 119 /** 120 * sum 121 * 122 * @param a an array of float. 123 * @param b an array of float. 124 * @return Returns the array ret. 125 */ 126 public static float[] add(float[] a, float[] b) { 127 return add(a, b, new float[3]); 128 } 129 130 /** 131 * sum 132 * 133 * @param a an array of float. 134 * @param b an array of float. 135 * @param ret an array of float. 136 * @return Returns the array ret. 137 */ 138 public static float[] add(float[] a, float[] b, float[] ret) { 139 ret[0] = a[0] + b[0]; 140 ret[1] = a[1] + b[1]; 141 ret[2] = a[2] + b[2]; 142 return ret; 143 } 144 145 /** 146 * angle 147 * 148 * @param i an array of float. 149 * @param j an array of float. 150 * @return Returns the angle. 151 */ 152 public static float angle(float[] i, float[] j) { 153 var x = dot(normalize(i), normalize(j)); 154 if (abs(x) > 1) { 155 logger.warning(format(" Angle: abs(dot) > 1 %10.6f", x)); 156 if (x > 0) { 157 x = 1; 158 } else { 159 x = -1; 160 } 161 } 162 return (float) acos(x); 163 } 164 165 /** 166 * Finds the angle formed by three atoms 167 * 168 * @param i Atom position vector. 169 * @param j Atom position vector (central atom). 170 * @param k Atom position vector. 171 * @return Returns the angle in the range [ -pi, pi ]. 172 */ 173 public static float bondAngle(float[] i, float[] j, float[] k) { 174 return angle(sub(i, j), sub(k, j)); 175 } 176 177 /** 178 * Finds the dihedral angle formed between 4 atoms 179 * 180 * @param a Atom position vector. 181 * @param b Atom position vector. 182 * @param c Atom position vector. 183 * @param d Atom position vector. 184 * @return The dihedral angle in the range [ -pi, pi ]. 185 */ 186 public static float dihedralAngle(float[] a, float[] b, float[] c, float[] d) { 187 var ba = sub(b, a); 188 var cb = sub(c, b); 189 var dc = sub(d, c); 190 var t = X(ba, cb); 191 var u = X(cb, dc); 192 var rt = dot(t, t); 193 var ru = dot(u, u); 194 var rtu = (float) sqrt(rt * ru); 195 if (rtu != 0.0) { 196 var rcb = length(cb); 197 var cosine = dot(t, u) / rtu; 198 var tu = X(t, u); 199 var sine = dot(cb, tu) / (rcb * rtu); 200 cosine = min(1.0f, max(-1.0f, cosine)); 201 var angle = (float) acos(cosine); 202 if (sine < 0.0) { 203 angle = -angle; 204 } 205 return angle; 206 } 207 return 0; 208 } 209 210 /** 211 * Finds the distance between two vectors. 212 * 213 * @param a First vector. 214 * @param b Second vector. 215 * @return Returns the distance between vectors a and b. 216 */ 217 public static float dist(float[] a, float[] b) { 218 return (float) sqrt(dist2(a, b)); 219 } 220 221 /** 222 * Finds the squared distance between two vectors. 223 * 224 * @param a First vector. 225 * @param b Second vector. 226 * @return Returns the squared distance between vectors a and b. 227 */ 228 public static float dist2(float[] a, float[] b) { 229 var dx = a[0] - b[0]; 230 var dy = a[1] - b[1]; 231 var dz = a[2] - b[2]; 232 return Math.fma(dx, dx, Math.fma(dy, dy, dz * dz)); 233 } 234 235 /** 236 * Finds the dot product between two vectors. 237 * 238 * @param a First vector. 239 * @param b Second vector. 240 * @return Returns the dot product of a and b. 241 */ 242 public static float dot(float[] a, float[] b) { 243 return Math.fma(a[0], b[0], Math.fma(a[1], b[1], a[2] * b[2])); 244 } 245 246 /** 247 * Finds the length of a vector. 248 * 249 * @param d A vector to find the length of. 250 * @return Returns the length of vector d. 251 */ 252 public static float length(float[] d) { 253 return (float) sqrt(length2(d)); 254 } 255 256 /** 257 * Finds the length of a vector squared. 258 * 259 * @param d A vector to find the length of squared. 260 * @return Returns the length of vector d squared. 261 */ 262 public static float length2(float[] d) { 263 return Math.fma(d[0], d[0], Math.fma(d[1], d[1], d[2] * d[2])); 264 } 265 266 /** 267 * logVector 268 * 269 * @param v an array of float. 270 */ 271 public static void log(float[] v) { 272 logger.info(toString(v)); 273 } 274 275 /** 276 * logVector. 277 * 278 * @param v an array of {@link float} objects. 279 * @param label a {@link String} object. 280 */ 281 public static void log(float[] v, String label) { 282 logger.info(toString(v, label)); 283 } 284 285 /** 286 * Normalizes a vector. 287 * 288 * @param n A vector to be normalized. 289 * @return Returns the normalized vector. 290 */ 291 public static float[] normalize(float[] n) { 292 return scale(n, (float) 1.0 / length(n), new float[3]); 293 } 294 295 /** 296 * Normalizes a vector. 297 * 298 * @param n A vector to be normalized. 299 * @param ret The normalized vector. 300 * @return Returns the normalized vector. 301 */ 302 public static float[] normalize(float[] n, float[] ret) { 303 return scale(n, (float) 1.0 / length(n), ret); 304 } 305 306 /** 307 * Scales a vector. 308 * 309 * @param n A vector to be scaled. 310 * @param a A scalar value. 311 * @return Returns the scaled vector. 312 */ 313 public static float[] scale(float[] n, float a) { 314 return scale(n, a, new float[3]); 315 } 316 317 /** 318 * Scales a vector. 319 * 320 * @param n A vector to be scaled. 321 * @param a A scalar value. 322 * @param ret The scaled Vector. 323 * @return Returns the array ret. 324 */ 325 public static float[] scale(float[] n, float a, float[] ret) { 326 ret[0] = n[0] * a; 327 ret[1] = n[1] * a; 328 ret[2] = n[2] * a; 329 return ret; 330 } 331 332 333 /** 334 * Squares values of a vector. 335 * 336 * @param n A vector to be squared. 337 * @return Returns the squared vector. 338 */ 339 public static float[] square(float[] n) { 340 return square(n, new float[3]); 341 } 342 343 /** 344 * Squares values of a vector. 345 * 346 * @param n A vector to be squared. 347 * @param ret The squared vector. 348 * @return Returns the array ret. 349 */ 350 public static float[] square(float[] n, float[] ret) { 351 ret[0] = n[0] * n[0]; 352 ret[1] = n[1] * n[1]; 353 ret[2] = n[2] * n[2]; 354 return ret; 355 } 356 357 /** 358 * Square root values of a vector. 359 * 360 * @param n A vector to determine square root. 361 * @return Returns the rooted vector. 362 */ 363 public static float[] squareRoot(float[] n) { 364 return squareRoot(n, new float[3]); 365 } 366 367 /** 368 * Square root values of a vector. 369 * 370 * @param n A vector to determine square root. 371 * @param ret The rooted vector. 372 * @return Returns the array ret. 373 */ 374 public static float[] squareRoot(float[] n, float[] ret) { 375 ret[0] = (float) sqrt(n[0]); 376 ret[1] = (float) sqrt(n[1]); 377 ret[2] = (float) sqrt(n[2]); 378 return ret; 379 } 380 381 /** 382 * Finds the difference between two vectors. 383 * 384 * @param a First vector 385 * @param b Second vector 386 * @return Returns the difference ret. 387 */ 388 public static float[] sub(float[] a, float[] b) { 389 return sub(a, b, new float[3]); 390 } 391 392 /** 393 * Finds the difference between two vectors 394 * 395 * @param a First vector 396 * @param b Second vector 397 * @param ret Return Values 398 * @return Returns the difference ret. 399 */ 400 public static float[] sub(float[] a, float[] b, float[] ret) { 401 ret[0] = a[0] - b[0]; 402 ret[1] = a[1] - b[1]; 403 ret[2] = a[2] - b[2]; 404 return ret; 405 } 406 407 /** 408 * logVector. 409 * 410 * @param v an array of float. 411 * @return Returns a String description of the vector. 412 */ 413 public static String toString(float[] v) { 414 StringBuilder sb = new StringBuilder("Vector ( "); 415 for (float d : v) { 416 sb.append(format("%g ", d)); 417 } 418 sb.append(")"); 419 return sb.toString(); 420 } 421 422 /** 423 * vectorToString. 424 * 425 * @param v an array of {@link float} objects. 426 * @param label a {@link String} object. 427 * @return Returns a String description of the vector. 428 */ 429 public static String toString(@Nullable float[] v, @Nullable String label) { 430 if (v == null) { 431 return null; 432 } 433 if (label == null) { 434 label = "v"; 435 } 436 StringBuilder sb = new StringBuilder(format(" %16s = [", label)); 437 for (float value : v) { 438 sb.append(format(" %16.8f", value)); 439 } 440 sb.append(" ]"); 441 return sb.toString(); 442 } 443 }