1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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
57
58
59
60
61
62
63
64
65
66
67
68
69
70 public class StructureMetrics {
71
72 private static final Logger logger = Logger.getLogger(StructureMetrics.class.getName());
73
74
75
76
77
78
79
80 public static double radiusOfGyration(Atom[] atoms) {
81 Double3 radius = new Double3(radiusOfGyrationComponents(atoms));
82 return radius.length();
83 }
84
85
86
87
88
89
90
91
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
102
103
104
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
114
115
116
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
139
140
141
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
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
163
164
165
166
167
168
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
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
191 Double3 radius = new Double3();
192 if (pmp) {
193
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
225 Array2DRowRealMatrix cMatrix = new Array2DRowRealMatrix(tensor, false);
226 EigenDecomposition eigenDecomposition = new EigenDecomposition(cMatrix);
227
228 radius.set(eigenDecomposition.getRealEigenvalues()).scaleI(1.0 / nAtoms).sqrtI();
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
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
266
267
268
269
270
271
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
295
296
297
298
299
300
301
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
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
324
325
326
327
328
329
330
331
332
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
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
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
395 Array2DRowRealMatrix cMatrix = new Array2DRowRealMatrix(tensor, false);
396 EigenDecomposition eigenDecomposition = new EigenDecomposition(cMatrix);
397
398 moment = eigenDecomposition.getRealEigenvalues();
399 vec = eigenDecomposition.getV().getData();
400
401
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
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
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
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
458 if (print) {
459 logger.info(format("\n Center of Mass Coordinates: %8.4f %8.4f %8.4f", xcm, ycm, zcm));
460
461 double[] angles = new Rotation(vec, 1.0E-7).getAngles(
462 RotationOrder.XYZ, RotationConvention.VECTOR_OPERATOR);
463 double radian = 180 / PI;
464
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 }