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 ffx.crystal.SymOp.asMatrixString;
41 import static ffx.potential.parsers.DistanceMatrixFilter.writeDistanceMatrixRow;
42 import static ffx.utilities.StringUtils.writeAtomRanges;
43 import static java.lang.String.format;
44 import static java.lang.System.arraycopy;
45 import static java.util.Arrays.fill;
46 import static org.apache.commons.io.FilenameUtils.concat;
47 import static org.apache.commons.io.FilenameUtils.getBaseName;
48 import static org.apache.commons.io.FilenameUtils.getFullPath;
49 import static org.apache.commons.math3.util.FastMath.sqrt;
50
51 import edu.rit.mp.DoubleBuf;
52 import edu.rit.pj.Comm;
53 import ffx.crystal.SymOp;
54 import ffx.numerics.math.Double3;
55 import ffx.numerics.math.RunningStatistics;
56 import ffx.potential.AssemblyState;
57 import ffx.potential.ForceFieldEnergy;
58 import ffx.potential.MolecularAssembly;
59 import ffx.potential.bonded.Atom;
60 import ffx.potential.parsers.DistanceMatrixFilter;
61 import ffx.potential.parsers.SystemFilter;
62 import ffx.potential.parsers.XYZFilter;
63 import java.io.File;
64 import java.util.ArrayList;
65 import java.util.logging.Logger;
66 import org.apache.commons.math3.linear.Array2DRowRealMatrix;
67 import org.apache.commons.math3.linear.EigenDecomposition;
68
69 public class Superpose {
70
71
72
73
74 private static final Logger logger = Logger.getLogger(Superpose.class.getName());
75
76 private final SystemFilter baseFilter;
77 private final SystemFilter targetFilter;
78 private final boolean isSymmetric;
79
80 private final int baseSize;
81 private final int targetSize;
82
83 private int restartRow;
84 private int restartColumn;
85 private final double[] distRow;
86
87
88
89
90 private final Comm world;
91
92
93
94 private final int numProc;
95
96
97
98 private final int rank;
99
100
101
102
103 private final double[][] distances;
104
105
106
107 private final DoubleBuf[] buffers;
108
109
110
111 private final double[] myDistance;
112
113
114
115 private final DoubleBuf myBuffer;
116
117
118 public Superpose(SystemFilter baseFilter, SystemFilter targetFilter, boolean isSymmetric) {
119 this.baseFilter = baseFilter;
120 this.targetFilter = targetFilter;
121 this.isSymmetric = isSymmetric;
122
123
124 baseSize = baseFilter.countNumModels();
125 targetSize = targetFilter.countNumModels();
126
127
128 distRow = new double[targetSize];
129 fill(distRow, -1.0);
130
131 world = Comm.world();
132
133 numProc = world.size();
134
135 rank = world.rank();
136 if (numProc > 1) {
137 logger.info(format(" Number of MPI Processes: %d", numProc));
138 logger.info(format(" Rank of this MPI Process: %d", rank));
139 }
140
141 distances = new double[numProc][1];
142
143
144 for (int i = 0; i < numProc; i++) {
145 fill(distances[i], -1.0);
146 }
147
148
149 buffers = new DoubleBuf[numProc];
150 for (int i = 0; i < numProc; i++) {
151 buffers[i] = DoubleBuf.buffer(distances[i]);
152 }
153
154
155 myDistance = distances[rank];
156 myBuffer = buffers[rank];
157 }
158
159
160
161
162
163
164
165
166
167
168
169 public void calculateRMSDs(int[] usedIndices, boolean dRMSD, boolean verbose, boolean restart,
170 boolean write, boolean saveSnapshots, boolean printSym) {
171
172 String filename = baseFilter.getFile().getAbsolutePath();
173 String targetFilename = targetFilter.getFile().getAbsolutePath();
174
175
176 File baseOutputFile = null;
177 File targetOutputFile = null;
178 SystemFilter baseOutputFilter = null;
179 SystemFilter targetOutputFilter = null;
180 if (saveSnapshots) {
181 String baseOutputName = concat(getFullPath(filename),
182 getBaseName(filename) + "_superposed.arc");
183 String targetOutputName = concat(getFullPath(filename),
184 getBaseName(targetFilename) + "_superposed.arc");
185 baseOutputFile = SystemFilter.version(new File(baseOutputName));
186 targetOutputFile = SystemFilter.version(new File(targetOutputName));
187 MolecularAssembly baseAssembly = baseFilter.getActiveMolecularSystem();
188 MolecularAssembly targetAssembly = targetFilter.getActiveMolecularSystem();
189 baseOutputFilter = new XYZFilter(baseOutputFile, baseAssembly,
190 baseAssembly.getForceField(),
191 baseAssembly.getProperties());
192 targetOutputFilter = new XYZFilter(targetOutputFile, targetAssembly,
193 targetAssembly.getForceField(),
194 targetAssembly.getProperties());
195 }
196
197 String matrixFilename = concat(getFullPath(filename), getBaseName(filename) + ".dst");
198 RunningStatistics runningStatistics;
199 if (restart) {
200
201 runningStatistics = readMatrix(matrixFilename, isSymmetric, baseSize, targetSize);
202 if (runningStatistics == null) {
203 runningStatistics = new RunningStatistics();
204 }
205 } else {
206 runningStatistics = new RunningStatistics();
207 File file = new File(matrixFilename);
208 if (file.exists() && file.delete()) {
209 logger.info(format(" RMSD file (%s) was deleted.", matrixFilename));
210 logger.info(" To restart from a previous run, use the '-r' flag.");
211 }
212 }
213
214
215 int nUsed = usedIndices.length;
216 int nUsedVars = nUsed * 3;
217
218
219 MolecularAssembly baseMolecularAssembly = baseFilter.getActiveMolecularSystem();
220 ForceFieldEnergy baseForceFieldEnergy = baseMolecularAssembly.getPotentialEnergy();
221 int nVars = baseForceFieldEnergy.getNumberOfVariables();
222 double[] baseCoords = new double[nVars];
223 baseForceFieldEnergy.getCoordinates(baseCoords);
224 double[] baseUsedCoords = new double[nUsedVars];
225
226
227 Atom[] atoms = baseMolecularAssembly.getAtomArray();
228 double[] mass = new double[nUsed];
229 for (int i = 0; i < nUsed; i++) {
230 mass[i] = atoms[usedIndices[i]].getMass();
231 }
232
233
234 MolecularAssembly targetMolecularAssembly = targetFilter.getActiveMolecularSystem();
235 ForceFieldEnergy targetForceFieldEnergy = targetMolecularAssembly.getPotentialEnergy();
236 double[] targetCoords = new double[nVars];
237 double[] targetUsedCoords = new double[nUsedVars];
238
239
240 for (int row = 0; row < restartRow; row++) {
241 baseFilter.readNext(false, false);
242 }
243
244
245
246
247 int paddedTargetSize = targetSize;
248 int extra = targetSize % numProc;
249 if (extra != 0) {
250 paddedTargetSize = targetSize - extra + numProc;
251 logger.fine(format(" Target size %d vs. Padded size %d", targetSize, paddedTargetSize));
252 }
253
254
255 for (int row = restartRow; row < baseSize; row++) {
256
257
258 myDistance[0] = -1.0;
259
260 if (row == restartRow) {
261 if (dRMSD) {
262 logger.info(
263 "\n Coordinate RMSD\n Snapshots Original After Translation After Rotation dRMSD");
264 } else if (verbose) {
265 logger.info(
266 "\n Coordinate RMSD\n Snapshots Original After Translation After Rotation");
267 }
268 }
269
270
271 for (int column = restartColumn; column < paddedTargetSize; column++) {
272
273
274 if (column < targetSize) {
275 int targetRank = column % numProc;
276 if (targetRank == rank) {
277 if (isSymmetric && row == column) {
278
279 myDistance[0] = 0.0;
280 if (verbose) {
281 logger.info(format(" %6d %6d %s %8.5f",
282 row + 1, column + 1, "Diagonal", myDistance[0]));
283 }
284 } else if (isSymmetric && row > column) {
285
286 myDistance[0] = -1.0;
287 } else {
288
289
290
291 baseForceFieldEnergy.getCoordinates(baseCoords);
292
293
294 AssemblyState origStateB = new AssemblyState(targetMolecularAssembly);
295
296
297 targetForceFieldEnergy.getCoordinates(targetCoords);
298
299 extractCoordinates(usedIndices, baseCoords, baseUsedCoords);
300 extractCoordinates(usedIndices, targetCoords, targetUsedCoords);
301
302 double origRMSD = rmsd(baseUsedCoords, targetUsedCoords, mass);
303
304
305 double[] baseTranslation = calculateTranslation(baseUsedCoords, mass);
306 applyTranslation(baseCoords, baseTranslation);
307 double[] targetTranslation = calculateTranslation(targetUsedCoords, mass);
308 applyTranslation(targetCoords, targetTranslation);
309
310 extractCoordinates(usedIndices, baseCoords, baseUsedCoords);
311 extractCoordinates(usedIndices, targetCoords, targetUsedCoords);
312 double translatedRMSD = rmsd(baseUsedCoords, targetUsedCoords, mass);
313
314
315 double[][] rotation = calculateRotation(baseUsedCoords, targetUsedCoords, mass);
316 applyRotation(targetCoords, rotation);
317
318 extractCoordinates(usedIndices, targetCoords, targetUsedCoords);
319 double rotatedRMSD = rmsd(baseUsedCoords, targetUsedCoords, mass);
320
321 if (dRMSD) {
322 double disRMSD = calcDRMSD(baseUsedCoords, targetUsedCoords, nUsed * 3);
323 logger.info(format(" %6d %6d %8.5f %8.5f %8.5f %8.5f",
324 row + 1, column + 1, origRMSD, translatedRMSD, rotatedRMSD, disRMSD));
325 } else if (verbose) {
326 logger.info(format(" %6d %6d %8.5f %8.5f %8.5f",
327 row + 1, column + 1, origRMSD, translatedRMSD, rotatedRMSD));
328 }
329
330
331 myDistance[0] = rotatedRMSD;
332
333 if(printSym){
334 StringBuilder sbSO = new StringBuilder(
335 format("\n Sym Op to move %s onto %s:\nsymop ", targetFilename, filename));
336 StringBuilder sbInv = new StringBuilder(
337 format("\n Inverted Sym Op to move %s onto %s:\nsymop ", filename, targetFilename));
338 SymOp bestBaseSymOp = new SymOp(SymOp.ZERO_ROTATION, SymOp.Tr_0_0_0).append(new SymOp(SymOp.ZERO_ROTATION, baseTranslation).append(SymOp.invertSymOp(new SymOp(rotation, targetTranslation))));
339 double[] inverseBaseTranslation = new double[]{-baseTranslation[0], -baseTranslation[1], -baseTranslation[2]};
340 SymOp bestTargetSymOp = new SymOp(SymOp.ZERO_ROTATION, SymOp.Tr_0_0_0).append(new SymOp(SymOp.ZERO_ROTATION, targetTranslation).append(new SymOp(rotation, inverseBaseTranslation)));
341 ArrayList<Integer> mol1List = new ArrayList<>();
342 ArrayList<Integer> mol2List = new ArrayList<>();
343 Atom[] atomArr1 = baseMolecularAssembly.getAtomArray();
344 Atom[] atomArr2 = targetMolecularAssembly.getAtomArray();
345 int nAtoms = atomArr1.length;
346 for(int i = 0; i < nAtoms; i++){
347 if(atomArr1[i].isActive()){
348 mol1List.add(i);
349 }
350 }
351 nAtoms = atomArr2.length;
352 for(int i = 0; i < nAtoms; i++){
353 if(atomArr2[i].isActive()){
354 mol2List.add(i);
355 }
356 }
357 int[] mol1arr = mol1List.stream().mapToInt(Integer::intValue).toArray();
358 int[] mol2arr = mol2List.stream().mapToInt(Integer::intValue).toArray();
359 sbSO.append(format(" %s %s", writeAtomRanges(mol2arr), writeAtomRanges(mol1arr)))
360 .append(asMatrixString(bestBaseSymOp));
361 sbInv.append(format(" %s %s", writeAtomRanges(mol2arr), writeAtomRanges(mol1arr))).append(asMatrixString(bestTargetSymOp));
362 logger.info(format(" %s\n %s", sbSO, sbInv));
363 }
364
365
366 if (saveSnapshots && numProc == 1) {
367 MolecularAssembly molecularAssembly = targetOutputFilter.getActiveMolecularSystem();
368 molecularAssembly.getPotentialEnergy().setCoordinates(targetCoords);
369 targetOutputFilter.writeFile(targetOutputFile, true);
370 molecularAssembly = baseOutputFilter.getActiveMolecularSystem();
371 molecularAssembly.getPotentialEnergy().setCoordinates(baseCoords);
372 baseOutputFilter.writeFile(baseOutputFile, true);
373 origStateB.revertState();
374 }
375 }
376 }
377
378
379 targetFilter.readNext(false, false);
380 }
381
382
383 if ((column + 1) % numProc == 0) {
384 gatherRMSDs(row, column, runningStatistics);
385 }
386 }
387
388
389 restartColumn = 0;
390
391
392 targetFilter.readNext(true, false);
393
394
395 baseFilter.readNext(false, false);
396
397
398 if (rank == 0 && write) {
399 int firstColumn = 0;
400 if (isSymmetric) {
401 firstColumn = row;
402 }
403 writeDistanceMatrixRow(matrixFilename, distRow, firstColumn);
404 }
405 }
406
407 baseFilter.closeReader();
408 targetFilter.closeReader();
409
410 logger.info(format(" RMSD Minimum: %8.6f", runningStatistics.getMin()));
411 logger.info(format(" RMSD Maximum: %8.6f", runningStatistics.getMax()));
412 logger.info(format(" RMSD Mean: %8.6f", runningStatistics.getMean()));
413 double variance = runningStatistics.getVariance();
414 if (!Double.isNaN(variance)) {
415 logger.info(format(" RMSD Variance: %8.6f", variance));
416 }
417 }
418
419
420
421
422
423
424
425
426 private void gatherRMSDs(int row, int column, RunningStatistics runningStatistics) {
427 try {
428 logger.finer(" Receiving results.");
429 world.allGather(myBuffer, buffers);
430 for (int i = 0; i < numProc; i++) {
431 int c = (column + 1) - numProc + i;
432 if (c < targetSize) {
433 distRow[c] = distances[i][0];
434 if (!isSymmetric) {
435 runningStatistics.addValue(distRow[c]);
436 } else if (c > row) {
437
438 runningStatistics.addValue(distRow[c]);
439 }
440 logger.finer(format(" %d %d %16.8f", row, c, distances[i][0]));
441 }
442 }
443 } catch (Exception e) {
444 logger.severe(" Exception collecting distance values." + e);
445 }
446 }
447
448
449
450
451
452
453
454
455
456
457 private RunningStatistics readMatrix(String filename, boolean isSymmetric, int expectedRows,
458 int expectedColumns) {
459 restartRow = 0;
460 restartColumn = 0;
461
462 DistanceMatrixFilter distanceMatrixFilter = new DistanceMatrixFilter();
463 RunningStatistics runningStatistics = distanceMatrixFilter.readDistanceMatrix(
464 filename, expectedRows, expectedColumns);
465
466 if (runningStatistics != null && runningStatistics.getCount() > 0) {
467 restartRow = distanceMatrixFilter.getRestartRow();
468 restartColumn = distanceMatrixFilter.getRestartColumn();
469
470 if (isSymmetric) {
471
472 if (restartRow == expectedRows && restartColumn == 1) {
473 logger.info(format(" Complete symmetric distance matrix found (%d x %d).", restartRow,
474 restartRow));
475 } else {
476 restartColumn = 0;
477 logger.info(format(
478 " Incomplete symmetric distance matrix found.\n Restarting at row %d, column %d.",
479 restartRow + 1, restartColumn + 1));
480 }
481 } else if (restartRow == expectedRows && restartColumn == expectedColumns) {
482 logger.info(format(" Complete distance matrix found (%d x %d).", restartRow, restartColumn));
483 } else {
484 restartColumn = 0;
485 logger.info(format(" Incomplete distance matrix found.\n Restarting at row %d, column %d.",
486 restartRow + 1, restartColumn + 1));
487 }
488 }
489
490 return runningStatistics;
491 }
492
493
494
495
496
497
498
499
500 public static void extractCoordinates(int[] usedIndices, double[] x, double[] xUsed) {
501 int nUsed = usedIndices.length;
502 for (int u = 0; u < nUsed; u++) {
503 int u3 = 3 * u;
504 int i3 = 3 * usedIndices[u];
505 arraycopy(x, i3, xUsed, u3, 3);
506 }
507 }
508
509
510
511
512
513
514
515
516
517 public static double calcDRMSD(double[] xUsed, double[] x2Used, int nUsed) {
518 double disRMSD = 0.0;
519 int counter = 0;
520 for (int i = 0; i < nUsed; i = i + 3) {
521 Double3 xi = new Double3(xUsed[i], xUsed[i + 1], xUsed[i + 2]);
522 Double3 x2i = new Double3(x2Used[i], x2Used[i + 1], x2Used[i + 2]);
523 for (int j = i + 3; j < nUsed; j = j + 3) {
524 Double3 xj = new Double3(xUsed[j], xUsed[j + 1], xUsed[j + 2]);
525 Double3 x2j = new Double3(x2Used[j], x2Used[j + 1], x2Used[j + 2]);
526 double dis1 = xi.sub(xj).length();
527 double dis2 = x2i.sub(x2j).length();
528 double diff = dis1 - dis2;
529 disRMSD += diff * diff;
530 counter++;
531 }
532 }
533 disRMSD = disRMSD / counter;
534 return sqrt(disRMSD);
535 }
536
537
538
539
540
541
542
543 public static void applyRotation(double[] x2, double[][] rot) {
544 int n = x2.length / 3;
545
546
547 for (int i = 0; i < n; i++) {
548 int k = i * 3;
549 double xrot = x2[k] * rot[0][0] + x2[k + 1] * rot[0][1] + x2[k + 2] * rot[0][2];
550 double yrot = x2[k] * rot[1][0] + x2[k + 1] * rot[1][1] + x2[k + 2] * rot[1][2];
551 double zrot = x2[k] * rot[2][0] + x2[k + 1] * rot[2][1] + x2[k + 2] * rot[2][2];
552 x2[k] = xrot;
553 x2[k + 1] = yrot;
554 x2[k + 2] = zrot;
555 }
556 }
557
558
559
560
561
562
563
564 public static void applyTranslation(double[] x, final double[] translation) {
565 int n = x.length / 3;
566 for (int i = 0; i < n; i++) {
567 int k = i * 3;
568 for (int j = 0; j < 3; j++) {
569 x[k + j] += translation[j];
570 }
571 }
572 }
573
574
575
576
577
578
579
580
581
582
583 public static double[][] calculateRotation(double[] x1, double[] x2, double[] mass) {
584
585
586 double xxyx = 0.0;
587 double xxyy = 0.0;
588 double xxyz = 0.0;
589 double xyyx = 0.0;
590 double xyyy = 0.0;
591 double xyyz = 0.0;
592 double xzyx = 0.0;
593 double xzyy = 0.0;
594 double xzyz = 0.0;
595
596 int n = x1.length / 3;
597 for (int i = 0; i < n; i++) {
598 int k = i * 3;
599 double weigh = mass[i];
600 xxyx = xxyx + weigh * x1[k] * x2[k];
601 xxyy = xxyy + weigh * x1[k + 1] * x2[k];
602 xxyz = xxyz + weigh * x1[k + 2] * x2[k];
603 xyyx = xyyx + weigh * x1[k] * x2[k + 1];
604 xyyy = xyyy + weigh * x1[k + 1] * x2[k + 1];
605 xyyz = xyyz + weigh * x1[k + 2] * x2[k + 1];
606 xzyx = xzyx + weigh * x1[k] * x2[k + 2];
607 xzyy = xzyy + weigh * x1[k + 1] * x2[k + 2];
608 xzyz = xzyz + weigh * x1[k + 2] * x2[k + 2];
609 }
610
611 double[][] c = new double[4][4];
612 c[0][0] = xxyx + xyyy + xzyz;
613 c[0][1] = xzyy - xyyz;
614 c[1][0] = c[0][1];
615 c[1][1] = xxyx - xyyy - xzyz;
616 c[0][2] = xxyz - xzyx;
617 c[2][0] = c[0][2];
618 c[1][2] = xxyy + xyyx;
619 c[2][1] = c[1][2];
620 c[2][2] = xyyy - xzyz - xxyx;
621 c[0][3] = xyyx - xxyy;
622 c[3][0] = c[0][3];
623 c[1][3] = xzyx + xxyz;
624 c[3][1] = c[1][3];
625 c[2][3] = xyyz + xzyy;
626 c[3][2] = c[2][3];
627 c[3][3] = xzyz - xxyx - xyyy;
628
629
630 Array2DRowRealMatrix cMatrix = new Array2DRowRealMatrix(c, false);
631 EigenDecomposition eigenDecomposition = new EigenDecomposition(cMatrix);
632
633
634 double[] q = eigenDecomposition.getEigenvector(0).toArray();
635
636
637 double q02 = q[0] * q[0];
638 double q12 = q[1] * q[1];
639 double q22 = q[2] * q[2];
640 double q32 = q[3] * q[3];
641 double[][] rot = new double[3][3];
642 rot[0][0] = q02 + q12 - q22 - q32;
643 rot[1][0] = 2.0 * (q[1] * q[2] - q[0] * q[3]);
644 rot[2][0] = 2.0 * (q[1] * q[3] + q[0] * q[2]);
645 rot[0][1] = 2.0 * (q[1] * q[2] + q[0] * q[3]);
646 rot[1][1] = q02 - q12 + q22 - q32;
647 rot[2][1] = 2.0 * (q[2] * q[3] - q[0] * q[1]);
648 rot[0][2] = 2.0 * (q[1] * q[3] - q[0] * q[2]);
649 rot[1][2] = 2.0 * (q[2] * q[3] + q[0] * q[1]);
650 rot[2][2] = q02 - q12 - q22 + q32;
651
652 return rot;
653 }
654
655
656
657
658
659
660
661
662 public static double[] calculateTranslation(double[] x, final double[] mass) {
663 double xmid = 0.0;
664 double ymid = 0.0;
665 double zmid = 0.0;
666 double norm = 0.0;
667 int n = x.length / 3;
668
669 for (int i = 0; i < n; i++) {
670 int k = 3 * i;
671 double weigh = mass[i];
672 xmid += x[k] * weigh;
673 ymid += x[k + 1] * weigh;
674 zmid += x[k + 2] * weigh;
675 norm += weigh;
676 }
677 xmid /= norm;
678 ymid /= norm;
679 zmid /= norm;
680
681 return new double[] {-xmid, -ymid, -zmid};
682 }
683
684
685
686
687
688
689
690
691
692 public static double rmsd(double[] x1, double[] x2, double[] mass) {
693 double rmsfit = 0.0;
694 double norm = 0.0;
695 int n = x1.length / 3;
696 for (int i = 0; i < n; i++) {
697 int k = 3 * i;
698 double weigh = mass[i];
699 double xr = x1[k] - x2[k];
700 double yr = x1[k + 1] - x2[k + 1];
701 double zr = x1[k + 2] - x2[k + 2];
702 double dist2 = xr * xr + yr * yr + zr * zr;
703 norm = norm + weigh;
704 double rmsterm = dist2 * weigh;
705 rmsfit = rmsfit + rmsterm;
706 }
707 return sqrt(rmsfit / norm);
708 }
709
710
711
712
713
714
715
716
717 public static void rotate(double[] x1, double[] x2, double[] mass) {
718 double[][] rotation = calculateRotation(x1, x2, mass);
719 applyRotation(x2, rotation);
720 }
721
722
723
724
725
726
727
728
729
730 public static void translate(double[] x1, double[] mass1, double[] x2, double[] mass2) {
731
732 translate(x1, mass1);
733
734 translate(x2, mass2);
735 }
736
737
738
739
740
741
742
743 public static void translate(double[] x, final double[] mass) {
744 double[] translation = calculateTranslation(x, mass);
745 applyTranslation(x, translation);
746 }
747
748
749
750
751
752
753
754
755
756
757
758 public static double superpose(double[] x1, double[] x2, double[] mass) {
759 translate(x1, mass);
760 translate(x2, mass);
761 rotate(x1, x2, mass);
762 return rmsd(x1, x2, mass);
763 }
764 }