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.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 }