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-2025.
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 java.lang.String.format;
41  import static org.apache.commons.math3.util.FastMath.PI;
42  
43  import ffx.numerics.math.Double3;
44  import ffx.potential.bonded.Atom;
45  
46  import java.util.Arrays;
47  import java.util.logging.Logger;
48  
49  import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
50  import org.apache.commons.math3.geometry.euclidean.threed.RotationConvention;
51  import org.apache.commons.math3.geometry.euclidean.threed.RotationOrder;
52  import org.apache.commons.math3.linear.Array2DRowRealMatrix;
53  import org.apache.commons.math3.linear.EigenDecomposition;
54  
55  /**
56   * Structure Metrics contains functionality to calculate characteristics of coordinate systems.
57   * <p>
58   * This includes: Gyrate computes the radius of gyration of a molecular system from its atomic
59   * coordinates.
60   * <p>
61   * Inertia computes the principal moments of inertia for the system, and optionally translates the
62   * center of mass to the origin and rotates the principal axes onto the global axes. Reference:
63   * Herbert Goldstein, "Classical Mechanics, 2nd Edition", Addison-Wesley, Reading, MA, 1980; see the
64   * Euler angle xyz convention in Appendix B
65   *
66   * @author Michael J. Schnieders
67   * @author Aaron J. Nessler
68   * @since 1.0
69   */
70  public class StructureMetrics {
71  
72    private static final Logger logger = Logger.getLogger(StructureMetrics.class.getName());
73  
74    /**
75     * Compute the radius of gyration for all atoms in the supplied array.
76     *
77     * @param atoms Atom array.
78     * @return The radius of gyration.
79     */
80    public static double radiusOfGyration(Atom[] atoms) {
81      Double3 radius = new Double3(radiusOfGyrationComponents(atoms));
82      return radius.length();
83    }
84  
85    /**
86     * Compute the average radius of gyration.
87     *
88     * @param x Array of atomic X-coordinates.
89     * @param y Array of atomic X-coordinates.
90     * @param z Array of atomic X-coordinates.
91     * @return The radius of gyration.
92     */
93    public static double radiusOfGyration(double[] x, double[] y, double[] z, double[] mass) {
94      assert (x.length == y.length);
95      assert (y.length == z.length);
96      Double3 radius = new Double3(radiusOfGyrationComponents(x, y, z, mass, true));
97      return radius.length();
98    }
99  
100   /**
101    * Compute the average radius of gyration.
102    *
103    * @param xyz Array of atomic coordinates (xyz = [X0, Y0, Z0, X1, Y1, Z1, ...].
104    * @return The radius of gyration.
105    */
106   public static double radiusOfGyration(double[] xyz, double[] mass) {
107     assert (xyz.length % 3 == 0);
108     Double3 radius = new Double3(radiusOfGyrationComponents(xyz, mass, true));
109     return radius.length();
110   }
111 
112   /**
113    * Compute the radius of gyration for all atoms in the supplied array.
114    *
115    * @param atoms Atom array.
116    * @return The radius of gyration.
117    */
118   public static double[] radiusOfGyrationComponents(Atom[] atoms) {
119     int nAtoms = atoms.length;
120     double[] x = new double[nAtoms];
121     double[] y = new double[nAtoms];
122     double[] z = new double[nAtoms];
123     double[] mass = new double[nAtoms];
124 
125     int index = 0;
126     for (int i = 0; i < nAtoms; i++) {
127       mass[i] = atoms[i].getMass();
128       x[index] = atoms[i].getX();
129       y[index] = atoms[i].getY();
130       z[index] = atoms[i].getZ();
131       index++;
132     }
133 
134     return radiusOfGyrationComponents(x, y, z, mass, false);
135   }
136 
137   /**
138    * Compute the components that make up the radius of gyration about yz-, xz-, xy-planes.
139    *
140    * @param xyz Coordinates for calculation
141    * @return radius of gyration along axes
142    */
143   public static double[] radiusOfGyrationComponents(double[] xyz, double[] mass, boolean pmp) {
144     assert (xyz.length % 3 == 0);
145     int nAtoms = xyz.length / 3;
146     // Find the centroid of the atomic coordinates.
147     double[] x = new double[nAtoms];
148     double[] y = new double[nAtoms];
149     double[] z = new double[nAtoms];
150 
151     for (int i = 0; i < nAtoms; i++) {
152       int index = i * 3;
153       x[i] = xyz[index++];
154       y[i] = xyz[index++];
155       z[i] = xyz[index];
156     }
157 
158     return radiusOfGyrationComponents(x, y, z, mass, pmp);
159   }
160 
161   /**
162    * Compute the components that make up the radius of gyration tensor about yz-, xz-, xy-planes.
163    *
164    * @param x   Coordinates for calculation
165    * @param y   Coordinates for calculation
166    * @param z   Coordinates for calculation
167    * @param pmp Principal moment plane
168    * @return radius of gyration about planes (yz, xz, xy).
169    */
170   public static double[] radiusOfGyrationComponents(double[] x, double[] y, double[] z,
171                                                     double[] mass, boolean pmp) {
172     assert (x.length == y.length);
173     assert (y.length == z.length);
174     assert (x.length <= mass.length);
175     double massSum = Arrays.stream(mass).sum();
176 
177     // Find the centroid of the atomic coordinates.
178     int nAtoms = x.length;
179     double xc = 0.0;
180     double yc = 0.0;
181     double zc = 0.0;
182     for (int i = 0; i < nAtoms; i++) {
183       xc += x[i] * mass[i];
184       yc += y[i] * mass[i];
185       zc += z[i] * mass[i];
186     }
187     Double3 centroid = new Double3(xc, yc, zc);
188     centroid.scaleI(1.0 / massSum);
189 
190     // Compute the radius of gyration.
191     Double3 radius = new Double3();
192     if (pmp) {
193       // gyration tensor about principal planes
194       double xterm;
195       double yterm;
196       double zterm;
197       double xx = 0.0;
198       double xy = 0.0;
199       double xz = 0.0;
200       double yy = 0.0;
201       double yz = 0.0;
202       double zz = 0.0;
203       for (int i = 0; i < nAtoms; i++) {
204         xterm = x[i] - xc;
205         yterm = y[i] - yc;
206         zterm = z[i] - zc;
207         xx += xterm * xterm;
208         xy += xterm * yterm;
209         xz += xterm * zterm;
210         yy += yterm * yterm;
211         yz += yterm * zterm;
212         zz += zterm * zterm;
213       }
214       double[][] tensor = new double[3][3];
215       tensor[0][0] = xx;
216       tensor[0][1] = xy;
217       tensor[0][2] = xz;
218       tensor[1][0] = xy;
219       tensor[1][1] = yy;
220       tensor[1][2] = yz;
221       tensor[2][0] = xz;
222       tensor[2][1] = yz;
223       tensor[2][2] = zz;
224       // Diagonalize the matrix.
225       Array2DRowRealMatrix cMatrix = new Array2DRowRealMatrix(tensor, false);
226       EigenDecomposition eigenDecomposition = new EigenDecomposition(cMatrix);
227       // Extract the quaternions.
228       radius.set(eigenDecomposition.getRealEigenvalues()).scaleI(1.0 / nAtoms).sqrtI();
229 //      vec = eigenDecomposition.getV().getData();
230 //
231 //      // Select the direction for each principal moment axis
232 //      double dot = 0.0;
233 //      for (int i = 0; i < 2; i++) {
234 //        for (int j = 0; j < nAtoms; j++) {
235 //          xterm = vec[i][0] * (x[j] - xc);
236 //          yterm = vec[i][1] * (y[j] - yc);
237 //          zterm = vec[i][2] * (z[j] - zc);
238 //          dot = xterm + yterm + zterm;
239 //          if (dot < 0.0) {
240 //            for (int k = 0; k < 3; k++) {
241 //              vec[i][k] = -vec[i][k];
242 //            }
243 //          }
244 //          if (dot != 0.0) {
245 //            break;
246 //          }
247 //        }
248 //        if (dot != 0.0) {
249 //          break;
250 //        }
251 //      }
252     } else {
253       for (int i = 0; i < nAtoms; i++) {
254         Double3 xyz = new Double3(x[i], y[i], z[i]);
255         xyz.subI(centroid);
256         radius.addI(xyz.squareI());
257       }
258       radius.scaleI(1.0 / nAtoms).sqrtI();
259     }
260 
261     return radius.get();
262   }
263 
264   /**
265    * Compute the moments of inertia for all atoms in the supplied array.
266    *
267    * @param atoms Atom array.
268    * @param moved Move coordinates
269    * @param print Display values to user
270    * @param pma   Use principal moment axes.
271    * @return The moments of inertia.
272    */
273   public static double[][] momentsOfInertia(Atom[] atoms, boolean moved, boolean print,
274                                             boolean pma) {
275     double[] mass = new double[atoms.length];
276     int nAtoms = atoms.length;
277     double[] x = new double[nAtoms];
278     double[] y = new double[nAtoms];
279     double[] z = new double[nAtoms];
280 
281     int index = 0;
282     for (Atom atom : atoms) {
283       mass[index] = atom.getMass();
284       x[index] = atom.getX();
285       y[index] = atom.getY();
286       z[index] = atom.getZ();
287       index++;
288     }
289 
290     return momentsOfInertia(x, y, z, mass, moved, print, pma);
291   }
292 
293   /**
294    * Compute the moments of inertia.
295    *
296    * @param xyz   Array of atomic coordinates (xyz = [X0, Y0, Z0, X1, Y1, Z1, ...].
297    * @param mass  Mass of atoms
298    * @param moved Move from original coordinates to selection
299    * @param print Display values to user
300    * @param pma   Use principal moment axes.
301    * @return The radius of gyration.
302    */
303   public static double[][] momentsOfInertia(double[] xyz, double[] mass, boolean moved,
304                                             boolean print, boolean pma) {
305     assert (xyz.length % 3 == 0);
306     int nAtoms = xyz.length / 3;
307     // Find the centroid of the atomic coordinates.
308     double[] x = new double[nAtoms];
309     double[] y = new double[nAtoms];
310     double[] z = new double[nAtoms];
311 
312     for (int i = 0; i < nAtoms; i++) {
313       int index = i * 3;
314       x[i] = xyz[index++];
315       y[i] = xyz[index++];
316       z[i] = xyz[index];
317     }
318 
319     return momentsOfInertia(x, y, z, mass, moved, print, pma);
320   }
321 
322   /**
323    * Compute the moments of inertia
324    *
325    * @param x     Array of atomic X-coordinates.
326    * @param y     Array of atomic X-coordinates.
327    * @param z     Array of atomic X-coordinates.
328    * @param mass  mass of atoms
329    * @param moved Move coordinates to principal axes
330    * @param print Print out values to screen
331    * @param pma   Report moments of inertia to principal axes.
332    * @return The moment of inertia.
333    */
334   public static double[][] momentsOfInertia(double[] x, double[] y, double[] z, double[] mass,
335                                             boolean moved, boolean print, boolean pma) {
336     assert (x.length == y.length);
337     assert (y.length == z.length);
338 
339     // Find the centroid of the atomic coordinates.
340     double total = 0.0;
341     double xcm = 0.0;
342     double ycm = 0.0;
343     double zcm = 0.0;
344     int nAtoms = x.length;
345     for (int i = 0; i < nAtoms; i++) {
346       double massValue = mass[i];
347       total += massValue;
348       xcm += x[i] * massValue;
349       ycm += y[i] * massValue;
350       zcm += z[i] * massValue;
351     }
352     xcm /= total;
353     ycm /= total;
354     zcm /= total;
355 
356     // Compute and then diagonalize the inertia tensor
357     double xx = 0.0;
358     double xy = 0.0;
359     double xz = 0.0;
360     double yy = 0.0;
361     double yz = 0.0;
362     double zz = 0.0;
363 
364     double xterm;
365     double yterm;
366     double zterm;
367 
368     for (int i = 0; i < nAtoms; i++) {
369       double massValue = mass[i];
370       xterm = x[i] - xcm;
371       yterm = y[i] - ycm;
372       zterm = z[i] - zcm;
373       xx += xterm * xterm * massValue;
374       xy += xterm * yterm * massValue;
375       xz += xterm * zterm * massValue;
376       yy += yterm * yterm * massValue;
377       yz += yterm * zterm * massValue;
378       zz += zterm * zterm * massValue;
379     }
380     double[][] tensor = new double[3][3];
381     tensor[0][0] = yy + zz;
382     tensor[0][1] = -xy;
383     tensor[0][2] = -xz;
384     tensor[1][0] = -xy;
385     tensor[1][1] = xx + zz;
386     tensor[1][2] = -yz;
387     tensor[2][0] = -xz;
388     tensor[2][1] = -yz;
389     tensor[2][2] = xx + yy;
390 
391     double[] moment;
392     double[][] vec;
393     if (pma) {
394       // Diagonalize the matrix.
395       Array2DRowRealMatrix cMatrix = new Array2DRowRealMatrix(tensor, false);
396       EigenDecomposition eigenDecomposition = new EigenDecomposition(cMatrix);
397       // Extract the quaternions.
398       moment = eigenDecomposition.getRealEigenvalues();
399       vec = eigenDecomposition.getV().getData();
400 
401       // Select the direction for each principal moment axis
402       double dot = 0.0;
403       for (int i = 0; i < 2; i++) {
404         for (int j = 0; j < nAtoms; j++) {
405           xterm = vec[i][0] * (x[j] - xcm);
406           yterm = vec[i][1] * (y[j] - ycm);
407           zterm = vec[i][2] * (z[j] - zcm);
408           dot = xterm + yterm + zterm;
409           if (dot < 0.0) {
410             for (int k = 0; k < 3; k++) {
411               vec[i][k] = -vec[i][k];
412             }
413           }
414           if (dot != 0.0) {
415             break;
416           }
417         }
418         if (dot != 0.0) {
419           break;
420         }
421       }
422 
423       // Moment axes must give a right-handed coordinate system.
424       xterm = vec[0][0] * (vec[1][1] * vec[2][2] - vec[2][1] * vec[1][2]);
425       yterm = vec[0][1] * (vec[2][0] * vec[1][2] - vec[1][0] * vec[2][2]);
426       zterm = vec[0][2] * (vec[1][0] * vec[2][1] - vec[2][0] * vec[1][1]);
427       dot = xterm + yterm + zterm;
428       if (dot < 0.0) {
429         for (int i = 0; i < 3; i++) {
430           vec[2][i] = -vec[2][i];
431         }
432       }
433 
434       // Principal moment axes form rows of Euler rotation matrix.
435       if (moved) {
436         double[][] a = new double[3][3];
437         for (int i = 0; i < 3; i++) {
438           for (int j = 0; j < 3; j++) {
439             a[j][i] = vec[j][i];
440           }
441         }
442         // Translate to origin, then apply Euler rotation matrix.
443         for (int i = 0; i < nAtoms; i++) {
444           xterm = x[i] - xcm;
445           yterm = y[i] - ycm;
446           zterm = z[i] - zcm;
447           x[i] = a[0][0] * xterm + a[1][0] * yterm + a[2][0] * zterm;
448           y[i] = a[0][1] * xterm + a[1][1] * yterm + a[2][1] * zterm;
449           z[i] = a[0][2] * xterm + a[1][2] * yterm + a[2][2] * zterm;
450         }
451       }
452     } else {
453       vec = new double[][]{{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
454       moment = new double[]{tensor[0][0], tensor[1][1], tensor[2][2]};
455     }
456 
457     // print the center of mass and Euler angle values
458     if (print) {
459       logger.info(format("\n Center of Mass Coordinates: %8.4f %8.4f %8.4f", xcm, ycm, zcm));
460       // invert vec
461       double[] angles = new Rotation(vec, 1.0E-7).getAngles(
462           RotationOrder.XYZ, RotationConvention.VECTOR_OPERATOR);
463       double radian = 180 / PI;
464       // Convert to degrees
465       for (int i = 0; i < 3; i++) {
466         angles[i] *= radian;
467       }
468       logger.info(format(" Euler Angles (Phi/Theta/Psi): %8.3f %8.3f %8.3f", angles[0], angles[1], angles[2]));
469       logger.info(
470           " Moments of Inertia and Principle Axes:\n  Moments (amu Ang^2): \t X-, Y-, and Z-Components of Axes:");
471       for (int i = 0; i < 3; i++) {
472         logger.info(format("  %16.3f %12.6f %12.6f %12.6f", moment[i], vec[i][0], vec[i][1], vec[i][2]));
473       }
474     }
475 
476     double[][] momentsAndVectors = new double[3][4];
477     for (int i = 0; i < 3; i++) {
478       for (int j = 0; j < 4; j++) {
479         if (j == 0) {
480           momentsAndVectors[i][j] = moment[i];
481         } else {
482           momentsAndVectors[i][j] = vec[i][j - 1];
483         }
484       }
485     }
486 
487     return momentsAndVectors;
488   }
489 }