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.Tr_0_0_0;
41 import static ffx.crystal.SymOp.ZERO_ROTATION;
42 import static ffx.crystal.SymOp.applyCartesianSymOp;
43 import static ffx.crystal.SymOp.asMatrixString;
44 import static ffx.crystal.SymOp.invertSymOp;
45 import static ffx.numerics.math.ScalarMath.mod;
46 import static ffx.numerics.math.DoubleMath.dist2;
47 import static ffx.numerics.math.MatrixMath.mat4Mat4;
48 import static ffx.potential.parsers.DistanceMatrixFilter.writeDistanceMatrixRow;
49 import static ffx.potential.utils.StructureMetrics.momentsOfInertia;
50 import static ffx.potential.utils.StructureMetrics.radiusOfGyration;
51 import static ffx.potential.utils.StructureMetrics.radiusOfGyrationComponents;
52 import static ffx.potential.utils.Superpose.applyRotation;
53 import static ffx.potential.utils.Superpose.applyTranslation;
54 import static ffx.potential.utils.Superpose.calculateRotation;
55 import static ffx.potential.utils.Superpose.calculateTranslation;
56 import static ffx.potential.utils.Superpose.rmsd;
57 import static ffx.potential.utils.Superpose.rotate;
58 import static ffx.potential.utils.Superpose.translate;
59 import static ffx.utilities.StringUtils.parseAtomRanges;
60 import static ffx.utilities.Resources.logResources;
61 import static ffx.utilities.StringUtils.writeAtomRanges;
62 import static java.lang.String.format;
63 import static java.lang.System.arraycopy;
64 import static java.util.Arrays.copyOf;
65 import static java.util.Arrays.fill;
66 import static java.util.Arrays.sort;
67 import static org.apache.commons.io.FilenameUtils.getName;
68 import static org.apache.commons.math3.util.FastMath.abs;
69 import static org.apache.commons.math3.util.FastMath.sqrt;
70 import static org.apache.commons.math3.util.FastMath.max;
71 import static org.apache.commons.math3.util.FastMath.ceil;
72 import static org.apache.commons.math3.util.FastMath.min;
73 import static org.apache.commons.math3.util.FastMath.cbrt;
74 import static org.apache.commons.math3.util.FastMath.PI;
75
76 import edu.rit.mp.DoubleBuf;
77 import edu.rit.pj.Comm;
78 import ffx.crystal.Crystal;
79 import ffx.crystal.ReplicatesCrystal;
80 import ffx.crystal.SymOp;
81 import ffx.numerics.math.RunningStatistics;
82 import ffx.potential.MolecularAssembly;
83 import ffx.potential.Utilities;
84 import ffx.potential.bonded.Atom;
85 import ffx.potential.bonded.Bond;
86 import ffx.potential.bonded.MSNode;
87 import ffx.potential.parameters.ForceField;
88 import ffx.potential.parsers.DistanceMatrixFilter;
89 import ffx.potential.parsers.PDBFilter;
90 import ffx.potential.parsers.SystemFilter;
91 import ffx.potential.parsers.XYZFilter;
92 import ffx.utilities.DoubleIndexPair;
93 import ffx.utilities.IndexIndexPair;
94
95 import java.io.BufferedWriter;
96 import java.io.File;
97 import java.io.FileWriter;
98 import java.util.ArrayList;
99 import java.util.Arrays;
100 import java.util.List;
101 import java.util.logging.Level;
102 import java.util.logging.Logger;
103 import java.util.stream.IntStream;
104
105 import org.apache.commons.configuration2.CompositeConfiguration;
106 import org.apache.commons.io.FileUtils;
107 import org.apache.commons.io.FilenameUtils;
108 import org.apache.commons.lang3.ArrayUtils;
109
110
111
112
113
114
115
116
117 public class ProgressiveAlignmentOfCrystals {
118
119
120
121
122 private static final Logger logger = Logger.getLogger(
123 ProgressiveAlignmentOfCrystals.class.getName());
124
125
126
127 private final SystemFilter baseFilter;
128
129
130
131 private final int baseSize;
132
133
134
135 private final String baseLabel;
136
137
138
139 private final SystemFilter targetFilter;
140
141
142
143 private final int targetSize;
144
145
146
147 private final String targetLabel;
148
149
150
151 private final int numWorkItems;
152
153
154
155
156 private final boolean isSymmetric;
157
158
159
160 private String rmsdLabel;
161
162
163
164 public final double[] distRow;
165
166
167
168
169 private int restartRow = 0;
170
171
172
173
174
175 private int restartColumn = 0;
176
177
178
179 private final Comm world;
180
181
182
183 private final boolean useMPI;
184
185
186
187 private final int numProc;
188
189
190
191 private final int rank;
192
193
194
195
196 private final double[][] distances;
197
198
199
200 private final DoubleBuf[] buffers;
201
202
203
204 private final double[] myDistances;
205
206
207
208 private final DoubleBuf myBuffer;
209
210
211
212 private double[] massN;
213
214
215
216 private double massSum = 0;
217
218
219
220 private int[] comparisonAtoms;
221
222
223
224 private DoubleIndexPair[] baseAUDist;
225
226
227
228 private DoubleIndexPair[] baseAUDist_2;
229
230
231
232 private DoubleIndexPair[] targetAUDist;
233
234
235
236 private DoubleIndexPair[] targetAUDist_2;
237
238
239
240 private final List<Integer> tempListIndices = new ArrayList<>();
241
242
243
244 private Integer[] uniqueBase;
245
246
247
248 private Integer[] uniqueTarget;
249
250
251
252
253 private double[] baseXYZoriginal;
254
255
256
257
258
259 private double[] baseXYZ;
260
261
262
263
264 private double[] targetXYZoriginal;
265
266
267
268
269
270 private double[] targetXYZ;
271
272
273
274
275 private double[][][] baseCoM;
276
277
278
279
280 private double[][] baseAUoriginal;
281
282
283
284 private double[] baseAU;
285
286
287
288 private double[] base3AUs;
289
290
291
292
293 private double[][] baseNAUs;
294
295
296
297 private double[][][] bestBaseNAUs;
298
299
300
301 private double[][] targetCoM;
302
303
304
305
306 private double[][] targetAUoriginal;
307
308
309
310 private double[] targetAU;
311
312
313
314 private double[] target3AUs;
315
316
317
318 private double[] targetNAUs;
319
320
321
322
323 private double[][][] bestTargetNAUs;
324
325
326
327 private DoubleIndexPair[] pairedAUs;
328
329
330
331 private int numberOfHits = 0;
332
333
334
335 private File[] fileCache;
336
337
338
339 private String[] nameCache;
340
341
342
343 private ArrayList<ArrayList<Bond>> bondCache;
344
345
346
347 private Atom[][] atomCache;
348
349
350
351
352 private ForceField[] forceFieldCache;
353
354
355
356 private Crystal[] crystalCache;
357
358
359
360
361 private final double[] gyrations = new double[2];
362
363
364
365 private double[][] bestBaseMandV = new double[3][4];
366
367
368
369 private double[][] bestTargetMandV = new double[3][4];
370
371
372
373 private double[] bestBaseRg = new double[3];
374
375
376
377 private double[] bestTargetRg = new double[3];
378
379
380
381 private SymOp baseSymOp;
382
383
384
385 private SymOp[][] bestBaseSymOp;
386
387
388
389 private SymOp targetSymOp;
390
391
392
393 private SymOp[][] bestTargetSymOp;
394
395
396
397 private ArrayList<SymOp> baseSymOps = new ArrayList<>();
398
399
400
401 private ArrayList<Integer> baseDistMap = new ArrayList<>();
402
403
404
405 private ArrayList<SymOp> targetSymOps = new ArrayList<>();
406
407
408
409 private ArrayList<Integer> targetDistMap = new ArrayList<>();
410
411
412
413 private double printSym;
414
415
416
417 private final StringBuilder stringBuilder = new StringBuilder();
418
419
420
421
422 private SymOp baseTransformSymOp;
423
424
425
426
427 private SymOp targetTransformSymOp;
428
429
430
431 private SymOp[][] bestBaseTransformSymOp;
432
433
434
435 private SymOp[][] bestTargetTransformSymOp;
436
437
438
439 private static final double MATCH_TOLERANCE = 1.0E-12;
440
441
442
443
444 private static final double SYM_TOLERANCE = 2.0E-8;
445
446
447
448
449
450
451
452
453 public ProgressiveAlignmentOfCrystals(SystemFilter baseFilter, SystemFilter targetFilter,
454 boolean isSymmetric) {
455
456 this.baseFilter = baseFilter;
457 this.targetFilter = targetFilter;
458
459 this.isSymmetric = isSymmetric;
460
461
462 baseSize = baseFilter.countNumModels();
463 baseLabel = getName(baseFilter.getFile().getAbsolutePath());
464 targetSize = targetFilter.countNumModels();
465 targetLabel = getName(targetFilter.getFile().getAbsolutePath());
466
467 assert !isSymmetric || (baseSize == targetSize);
468
469 if (logger.isLoggable(Level.FINER)) {
470 logger.finer(format("\n Conformations for %s: %d", baseLabel, baseSize));
471 logger.finer(format(" Conformations for %s: %d", targetLabel, targetSize));
472 }
473
474
475 distRow = new double[targetSize];
476
477 CompositeConfiguration properties = baseFilter.getActiveMolecularSystem().getProperties();
478 useMPI = properties.getBoolean("pj.use.mpi", true);
479
480 if (useMPI) {
481 world = Comm.world();
482
483 numProc = world.size();
484
485 rank = world.rank();
486 } else {
487 world = null;
488 numProc = 1;
489 rank = 0;
490 }
491
492
493
494
495 int extra = targetSize % numProc;
496 int paddedTargetSize = targetSize;
497 if (extra != 0) {
498 paddedTargetSize = targetSize - extra + numProc;
499 }
500 numWorkItems = paddedTargetSize / numProc;
501
502 if (numProc > 1) {
503 logger.info(format(" Number of MPI Processes: %d", numProc));
504 logger.info(format(" Rank of this MPI Process: %d", rank));
505 logger.info(format(" Work per process per row: %d", numWorkItems));
506 }
507
508
509 fill(distRow, -1.0);
510
511
512 distances = new double[numProc][numWorkItems];
513
514
515 for (int i = 0; i < numProc; i++) {
516 fill(distances[i], -2.0);
517 }
518
519
520 buffers = new DoubleBuf[numProc];
521 for (int i = 0; i < numProc; i++) {
522 buffers[i] = DoubleBuf.buffer(distances[i]);
523 }
524
525
526 myDistances = distances[rank];
527 myBuffer = buffers[rank];
528 }
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561 private double compare(final File file1, final String name1, final List<Bond> bondList1,
562 final Atom[] atoms1, final ForceField forceField1, final File file2, final String name2,
563 final List<Bond> bondList2, final Atom[] atoms2, final ForceField forceField2, final int z1,
564 final int z2, final int compareAtomsSize, final int nAU, final int baseSearchValue,
565 final int targetSearchValue, final double matchTol, final int compNum, final boolean strict,
566 final int saveClusters, final boolean machineLearning, final int linkage, final boolean inertia,
567 final boolean gyrationComponents) {
568
569
570 boolean useSave = saveClusters > 0;
571 boolean useSym = printSym >= 0.0;
572 int nCoords = compareAtomsSize * 3;
573
574 int nBaseMols = baseXYZ.length / nCoords;
575 int nTargetMols = targetXYZ.length / nCoords;
576
577 if (logger.isLoggable(Level.FINER)) {
578 stringBuilder.append(format("\n Base: %s Target: %s", name1, name2));
579 stringBuilder.append(format("""
580
581 Comparing %3d of %5d in base sphere.
582 Comparing %3d of %5d in target sphere.
583 """, nAU, nBaseMols, nAU, nTargetMols));
584 }
585
586 if (useSym) {
587
588 if (bestBaseSymOp == null || bestBaseSymOp.length != z1 || bestBaseSymOp[0].length != z2) {
589 bestBaseSymOp = new SymOp[z1][z2];
590 }
591 if (bestTargetSymOp == null || bestTargetSymOp.length != z1
592 || bestTargetSymOp[0].length != z2) {
593 bestTargetSymOp = new SymOp[z1][z2];
594 }
595 if (bestBaseTransformSymOp == null || bestBaseTransformSymOp.length != z1
596 || bestBaseTransformSymOp[0].length != z2) {
597 bestBaseTransformSymOp = new SymOp[z1][z2];
598 }
599 if (bestTargetTransformSymOp == null || bestTargetTransformSymOp.length != z1
600 || bestTargetTransformSymOp[0].length != z2) {
601 bestTargetTransformSymOp = new SymOp[z1][z2];
602 }
603 }
604
605
606 if (logger.isLoggable(Level.FINER)) {
607 stringBuilder.append(" Prioritize target system.\n");
608
609 stringBuilder.append(" Search conformations of each crystal:\n");
610 }
611
612 arraycopy(baseXYZ, baseAUDist[0].index() * nCoords, baseAU, 0, nCoords);
613
614 tempListIndices.clear();
615
616 numberUniqueAUs(baseXYZ, baseAUDist, baseAU.clone(), nCoords, baseSearchValue, strict, nBaseMols,
617 massN, matchTol);
618 int baseLength = tempListIndices.size();
619 if (uniqueBase == null || uniqueBase.length != baseLength) {
620 uniqueBase = new Integer[baseLength];
621 }
622 tempListIndices.toArray(uniqueBase);
623
624
625 arraycopy(targetXYZ, targetAUDist[0].index() * nCoords, targetAU, 0, nCoords);
626 tempListIndices.clear();
627
628 numberUniqueAUs(targetXYZ, targetAUDist, targetAU.clone(), nCoords, targetSearchValue, strict,
629 nTargetMols, massN, matchTol);
630 int targetLength = tempListIndices.size();
631 if (uniqueTarget == null || uniqueTarget.length != targetLength) {
632 uniqueTarget = new Integer[targetLength];
633 }
634 tempListIndices.toArray(uniqueTarget);
635 final boolean reverse = baseLength > targetLength;
636 if (logger.isLoggable(Level.FINE)) {
637 stringBuilder.append(format("""
638
639 %d conformations detected out of %d in base crystal.
640 %d conformations detected out of %d in target crystal.
641 """, baseLength, baseSearchValue, targetLength, targetSearchValue));
642 }
643
644 if (useSym) {
645
646 if (z1 > 1 || z2 > 1) {
647 stringBuilder.append(
648 "\n Utilizing Sym Ops with Z''>1. Attempting to map all unique structures.\n");
649 }
650 }
651
652
653
654 double min1 = Double.MAX_VALUE;
655 double max1 = Double.MIN_VALUE;
656 int minBIndex = -1;
657 int minTIndex = -1;
658 int maxBIndex = -1;
659 int maxTIndex = -1;
660
661
662 IndexIndexPair[][] minZinds = new IndexIndexPair[z1][z2];
663 for (IndexIndexPair[] row : minZinds) {
664 fill(row, new IndexIndexPair(-1, -1));
665 }
666
667
668 double[][] minZdiffs = new double[z1][z2];
669 for (double[] row : minZdiffs) {
670 fill(row, Double.MAX_VALUE);
671 }
672
673 double[][] minDiffs = new double[baseLength][targetLength];
674 for (double[] row : minDiffs) {
675 fill(row, Double.MAX_VALUE);
676 }
677
678 for (int i = 0; i < baseLength; i++) {
679
680 int bIndex = uniqueBase[i];
681
682 int baseInd = baseAUDist[bIndex].index();
683
684 int currZ1 = baseDistMap.get(baseInd) % z1;
685
686 double[] baseXYZ = new double[nCoords];
687
688 for (int j = 0; j < targetLength; j++) {
689
690 arraycopy(this.baseXYZ, baseInd * nCoords, baseXYZ, 0, nCoords);
691
692 int tIndex = uniqueTarget[j];
693
694 int targetInd = targetAUDist[tIndex].index();
695
696 int currZ2 = targetDistMap.get(targetInd) % z2;
697
698 double[] targetXYZ = new double[nCoords];
699
700 arraycopy(this.targetXYZ, targetInd * nCoords, targetXYZ, 0, nCoords);
701
702 double[] tempTranB = calculateTranslation(baseXYZ, massN);
703 applyTranslation(baseXYZ, tempTranB);
704
705 double[] tempTranT = calculateTranslation(targetXYZ, massN);
706 applyTranslation(targetXYZ, tempTranT);
707
708 double[][] tempRotT = calculateRotation(baseXYZ, targetXYZ, massN);
709 applyRotation(targetXYZ, tempRotT);
710
711 double value = rmsd(baseXYZ, targetXYZ, massN);
712
713 if (logger.isLoggable(Level.FINEST)) {
714 stringBuilder.append(format(
715 "\n Comp %3d: Base %2d IND %3d (Z''=%2d) \t Target %2d IND %3d (Z''=%2d) Diff: %8.4f\n",
716 j + i * targetLength, bIndex, baseInd, currZ1, tIndex, targetInd, currZ2, value));
717 if (useSym) {
718 stringBuilder.append(" Base Temp Translation: ").append(Arrays.toString(tempTranB))
719 .append("\n");
720 stringBuilder.append(" Target Temp Translation: ").append(Arrays.toString(tempTranT));
721 stringBuilder.append(matrixToString(tempRotT, j + i * targetLength, "Temp Rot"))
722 .append("\n");
723 }
724 }
725
726 minDiffs[i][j] = value;
727
728 if (nAU == 1 && !strict) {
729
730 arraycopy(baseXYZ, 0, baseNAUs[0], 0, nAU * nCoords);
731 arraycopy(targetXYZ, 0, targetNAUs, 0, nAU * nCoords);
732 arraycopy(baseXYZ, 0, bestBaseNAUs[currZ1][currZ2], 0, nAU * nCoords);
733 arraycopy(targetXYZ, 0, bestTargetNAUs[currZ1][currZ2], 0, nAU * nCoords);
734 }
735 if (useSym && value < minZdiffs[currZ1][currZ2]) {
736 minZdiffs[currZ1][currZ2] = value;
737 minZinds[currZ1][currZ2] = new IndexIndexPair(baseInd, targetInd);
738 arraycopy(baseXYZ, 0, baseAU, 0, nCoords);
739 arraycopy(targetXYZ, 0, targetAU, 0, nCoords);
740 if (logger.isLoggable(Level.FINER)) {
741 stringBuilder.append(
742 format("\n Base AU: %2d of %2d" + "\t\tvs\t\t Target AU: %2d of %2d \t (%9.4f)",
743 currZ1 + 1, z1, currZ2 + 1, z2, minZdiffs[currZ1][currZ2]));
744 if (logger.isLoggable(Level.FINER)) {
745 stringBuilder.append("\n Base Translation: ").append(Arrays.toString(tempTranB))
746 .append("\n");
747 stringBuilder.append(" Target Translation: ").append(Arrays.toString(tempTranT));
748 stringBuilder.append(matrixToString(tempRotT, currZ1, "Target Rotation"));
749 }
750 }
751 baseTransformSymOp = new SymOp(ZERO_ROTATION, tempTranB);
752 bestBaseTransformSymOp[currZ1][currZ2] = baseTransformSymOp;
753
754 targetTransformSymOp = new SymOp(ZERO_ROTATION, tempTranT);
755 targetTransformSymOp = targetTransformSymOp.append(new SymOp(tempRotT, Tr_0_0_0));
756 bestTargetTransformSymOp[currZ1][currZ2] = targetTransformSymOp;
757 SymOp base = baseSymOps.get(baseDistMap.get(minZinds[currZ1][currZ2].sortedIndex()));
758 baseSymOp = new SymOp(base.rot, base.tr);
759 SymOp target = targetSymOps.get(
760 targetDistMap.get(minZinds[currZ1][currZ2].referenceIndex()));
761 targetSymOp = new SymOp(target.rot, target.tr);
762 bestBaseSymOp[currZ1][currZ2] = baseSymOp;
763 bestTargetSymOp[currZ1][currZ2] = targetSymOp;
764
765 if (logger.isLoggable(Level.FINEST)) {
766 stringBuilder.append("\n Base Sym Op Translation: ")
767 .append(Arrays.toString(bestBaseSymOp[currZ1][currZ2].tr));
768 stringBuilder.append(
769 matrixToString(bestBaseSymOp[currZ1][currZ2].rot, currZ2 + currZ1 * z2,
770 "Base Sym Op Rot")).append("\n");
771 stringBuilder.append(" Base Transform Translation: ")
772 .append(Arrays.toString(bestBaseTransformSymOp[currZ1][currZ2].tr));
773 stringBuilder.append(
774 matrixToString(bestBaseTransformSymOp[currZ1][currZ2].rot, currZ2 + currZ1 * z2,
775 "Base Transform Rot")).append("\n");
776 stringBuilder.append(" Target Sym Op Translation: ")
777 .append(Arrays.toString(bestTargetSymOp[currZ1][currZ2].tr));
778 stringBuilder.append(
779 matrixToString(bestTargetSymOp[currZ1][currZ2].rot, currZ2 + currZ1 * z2,
780 "Target Sym Op Rot")).append("\n");
781 stringBuilder.append(" Target Transform Translation: ")
782 .append(Arrays.toString(bestTargetTransformSymOp[currZ1][currZ2].tr));
783 stringBuilder.append(
784 matrixToString(bestTargetTransformSymOp[currZ1][currZ2].rot, currZ2 + currZ1 * z2,
785 "Target Transform Rot")).append("\n");
786 }
787 if (logger.isLoggable(Level.FINER)) {
788 stringBuilder.append(format(
789 " \n Saved %3d: Base %2d IND %3d (Z''=%2d) \t Target %2d IND %3d (Z''=%2d) Diff: %8.4f\n",
790 j + i * targetLength, bIndex, baseInd, currZ1, tIndex, targetInd, currZ2, value));
791 }
792 }
793
794 if (minDiffs[i][j] < min1) {
795 min1 = minDiffs[i][j];
796 minBIndex = bIndex;
797 minTIndex = tIndex;
798 }
799
800 if (minDiffs[i][j] > max1) {
801 max1 = minDiffs[i][j];
802 maxBIndex = bIndex;
803 maxTIndex = tIndex;
804 }
805 }
806 }
807
808 int baseBestZ = baseDistMap.get(baseAUDist[minBIndex].index()) % z1;
809 int targetBestZ = targetDistMap.get(targetAUDist[minTIndex].index()) % z2;
810 if (logger.isLoggable(Level.FINER)) {
811 stringBuilder.append(
812 format("\n Base min index: %d (Z''=%d) Target min Index %d (Z''=%d)", minBIndex, baseBestZ,
813 minTIndex, targetBestZ));
814 }
815
816 int maxLength = max(baseLength, targetLength);
817 int[] baseTargetMap = new int[maxLength];
818
819 if (reverse) {
820 for (int i = 0; i < baseLength; i++) {
821 double minValue = Double.MAX_VALUE;
822 for (int j = 0; j < targetLength; j++) {
823 if (minValue > minDiffs[i][j]) {
824 minValue = minDiffs[i][j];
825 baseTargetMap[i] = uniqueTarget[j];
826 }
827 }
828 }
829 } else {
830 for (int i = 0; i < targetLength; i++) {
831 double minValue = Double.MAX_VALUE;
832 for (int j = 0; j < baseLength; j++) {
833 if (minValue > minDiffs[j][i]) {
834 minValue = minDiffs[j][i];
835 baseTargetMap[i] = uniqueBase[j];
836 }
837 }
838 }
839 }
840
841
842 if (nAU == 1 && !strict) {
843
844
845 double baseGyration = radiusOfGyration(bestBaseNAUs[baseBestZ][targetBestZ].clone(), massN);
846 gyrations[0] = baseGyration;
847 double targetGyration = radiusOfGyration(bestTargetNAUs[baseBestZ][targetBestZ].clone(),
848 massN);
849 gyrations[1] = targetGyration;
850 if (useSym) {
851
852 final int numAtomsPerMol1 = atoms1.length/z1;
853 final int numAtomsPerMol2 = atoms2.length/z2;
854
855 boolean multisym = z1 > 1 || z2 > 1;
856 StringBuilder alchemicalAtoms = null;
857 StringBuilder alchemicalAtoms2 = null;
858 StringBuilder separation = new StringBuilder();
859 for (int i = 0; i < z1; i++) {
860 for (int j = 0; j < z2; j++) {
861 separation.append("\n Atom Separation (A) Description");
862 if (multisym) {
863 separation.append(format(" (Z1'' =%2d Z2'' =%2d)\n Base: Target: Distance:\n", i, j));
864 } else {
865 separation.append(" \n");
866 }
867 for (int k = 0; k < compareAtomsSize; k++) {
868 final int index = k * 3;
869 final double value = rmsd(
870 new double[]{bestBaseNAUs[i][j][index], bestBaseNAUs[i][j][index + 1],
871 bestBaseNAUs[i][j][index + 2]},
872 new double[]{bestTargetNAUs[i][j][index], bestTargetNAUs[i][j][index + 1],
873 bestTargetNAUs[i][j][index + 2]}, massN);
874 if (logger.isLoggable(Level.INFO)) {
875 if (printSym < value) {
876
877 separation.append(
878 format(" %4d %4d %14.6f\n", i * numAtomsPerMol1 + comparisonAtoms[k] + 1,
879 j * numAtomsPerMol2 + comparisonAtoms[k] + 1, value));
880
881 if (alchemicalAtoms == null) {
882 alchemicalAtoms = new StringBuilder();
883 if (multisym) {
884 alchemicalAtoms.append(format(" %3d %3d ", i, j));
885 }
886 alchemicalAtoms.append(i * numAtomsPerMol1 + comparisonAtoms[k] + 1);
887 alchemicalAtoms2 = new StringBuilder();
888 if (multisym) {
889 alchemicalAtoms2.append(format(" %3d %3d ", i, j));
890 }
891 alchemicalAtoms2.append(j * numAtomsPerMol2 + comparisonAtoms[k] + 1);
892 } else if (alchemicalAtoms.charAt(alchemicalAtoms.length() - 1) == '\n') {
893 alchemicalAtoms.append(format(" %3d %3d ", i, j))
894 .append(i * numAtomsPerMol1 + comparisonAtoms[k] + 1);
895 alchemicalAtoms2.append(format(" %3d %3d ", i, j))
896 .append(j * numAtomsPerMol2 + comparisonAtoms[k] + 1);
897 } else {
898 alchemicalAtoms.append(",")
899 .append(i * numAtomsPerMol1 + comparisonAtoms[k] + 1);
900 alchemicalAtoms2.append(",")
901 .append(j * numAtomsPerMol2 + comparisonAtoms[k] + 1);
902 }
903 }
904 }
905 }
906 if (useSave) {
907 if (machineLearning) {
908 saveAssembly(file1, name1, bondList1, atoms1, forceField1, bestBaseNAUs[i][j],
909 comparisonAtoms, nAU, "_c1", 0.000, compNum, saveClusters);
910 saveAssembly(file2, name2, bondList2, atoms2, forceField2, bestTargetNAUs[i][j],
911 comparisonAtoms, nAU, "_c2", min1, compNum, saveClusters);
912 } else {
913 saveAssembly(file1, name1, bondList1, atoms1, forceField1, bestBaseNAUs[i][j],
914 comparisonAtoms, nAU, "_c1", compNum, saveClusters);
915 saveAssembly(file2, name2, bondList2, atoms2, forceField2, bestTargetNAUs[i][j],
916 comparisonAtoms, nAU, "_c2", compNum, saveClusters);
917 }
918 }
919 if (alchemicalAtoms != null) {
920 alchemicalAtoms.append("\n");
921 alchemicalAtoms2.append("\n");
922 }
923 }
924 }
925 if (alchemicalAtoms != null) {
926 stringBuilder.append(separation);
927 stringBuilder.append("\n Suggested Alchemical Atoms:\n");
928
929 if (multisym) {
930 stringBuilder.append(" Base: (").append(baseLabel).append(")\n");
931 }
932 stringBuilder.append(alchemicalAtoms).append("\n");
933 if (multisym) {
934 stringBuilder.append(" Target: (").append(targetLabel).append(")\n")
935 .append(alchemicalAtoms2).append("\n");
936 }
937 }
938 stringBuilder.append(format(" Max RMSD: %9.4f B Index: %d T Index: %d", max1, maxBIndex, maxTIndex));
939 }
940 if (useSave && !useSym) {
941 if (machineLearning) {
942 saveAssembly(file1, name1, bondList1, atoms1, forceField1,
943 bestBaseNAUs[baseBestZ][targetBestZ], comparisonAtoms, nAU, "_c1", 0.000, compNum,
944 saveClusters);
945 saveAssembly(file2, name2, bondList2, atoms2, forceField2,
946 bestTargetNAUs[baseBestZ][targetBestZ], comparisonAtoms, nAU, "_c2", min1, compNum,
947 saveClusters);
948 } else {
949 saveAssembly(file1, name1, bondList1, atoms1, forceField1,
950 bestBaseNAUs[baseBestZ][targetBestZ], comparisonAtoms, nAU, "_c1", compNum, saveClusters);
951 saveAssembly(file2, name2, bondList2, atoms2, forceField2,
952 bestTargetNAUs[baseBestZ][targetBestZ], comparisonAtoms, nAU, "_c2", compNum, saveClusters);
953 }
954 }
955 if (logger.isLoggable(Level.FINER) && useSym) {
956 for (int i = 0; i < z1; i++) {
957 for (int j = 0; j < z2; j++) {
958 printSym(compareAtomsSize, file2, name2, bondList2, atoms2, forceField2, saveClusters, j);
959 }
960 }
961 }
962 return min1;
963 }
964 if (logger.isLoggable(Level.FINER)) {
965 stringBuilder.append(
966 " Minimum RMSD_1 Between Unique Base and Target AUs:\n i j (bInd tInd) RMSD_1\n");
967 for (int i = 0; i < baseLength; i++) {
968 for (int j = 0; j < targetLength; j++) {
969 stringBuilder.append(
970 format(" %d %d (%4d %4d) %4.4f\n", i, j, baseAUDist[uniqueBase[i]].index(),
971 targetAUDist[uniqueTarget[j]].index(), minDiffs[i][j]));
972 }
973 }
974 }
975
976
977 double bestRMSD = Double.MAX_VALUE;
978
979 for (double[] row : minZdiffs) {
980 fill(row, Double.MAX_VALUE);
981 }
982 if (logger.isLoggable(Level.FINE)) {
983 stringBuilder.append(
984 format("\n Trial RMSD_1 (%8s) RMSD_3 (%8s) %8s G(r1) G(r2)\n", rmsdLabel,
985 rmsdLabel, rmsdLabel));
986 }
987 baseAUDist_2 = new DoubleIndexPair[nBaseMols];
988 targetAUDist_2 = new DoubleIndexPair[nTargetMols];
989
990
991 int currentComparison = 1;
992 for (int l = 0; l < baseLength; l++) {
993 final int bIndex = uniqueBase[l];
994 final int baseInd = baseAUDist[bIndex].index();
995 final int currZ1 = baseDistMap.get(baseInd) % z1;
996
997 if (l > 0) {
998
999 arraycopy(baseXYZoriginal, 0, baseXYZ, 0, baseXYZoriginal.length);
1000 }
1001
1002 prioritizeReplicates(baseXYZ, baseCoM[0], compareAtomsSize, baseAUDist_2,
1003 baseInd, linkage);
1004
1005 int baseAUIndex = baseAUDist_2[0].index() * nCoords;
1006 arraycopy(baseXYZ, baseAUIndex, baseAU, 0, nCoords);
1007 for (int i = 0; i < 4; i++) {
1008 for (int j = 0; j < nAU; j++) {
1009 final int molIndex = baseAUDist_2[j].index() * nCoords;
1010 arraycopy(baseXYZ, molIndex, baseNAUs[i], j * nCoords, nCoords);
1011 }
1012 }
1013
1014 double[] translation = calculateTranslation(baseAU, massN);
1015 applyTranslation(baseAU, translation);
1016 applyTranslation(baseNAUs[1], translation);
1017 applyTranslation(baseNAUs[2], translation);
1018 applyTranslation(baseNAUs[3], translation);
1019 applyTranslation(baseXYZ, translation);
1020
1021 centerOfMass(baseCoM[1], baseXYZ, massN, massSum, compareAtomsSize);
1022 if (useSym) {
1023
1024 baseTransformSymOp = new SymOp(ZERO_ROTATION, Tr_0_0_0);
1025
1026 final SymOp base = baseSymOps.get(baseInd);
1027 baseSymOp = new SymOp(base.rot, base.tr);
1028 baseTransformSymOp = baseTransformSymOp.append(new SymOp(ZERO_ROTATION, translation));
1029 if (logger.isLoggable(Level.FINEST)) {
1030 stringBuilder.append(format("\n Base %d (%2d) to Origin Translation: ", l, bIndex))
1031 .append(Arrays.toString(translation)).append("\n");
1032 }
1033 }
1034
1035
1036 if (logger.isLoggable(Level.FINEST)) {
1037 stringBuilder.append("\n Base 3 Conformations:\n");
1038 }
1039 for (int i = 0; i < 3; i++) {
1040 int baseIndex = baseAUDist_2[i].index() * nCoords;
1041 arraycopy(baseXYZ, baseIndex, base3AUs, i * nCoords, nCoords);
1042 }
1043 translation = calculateTranslation(base3AUs, massN);
1044 applyTranslation(base3AUs, translation);
1045 applyTranslation(baseNAUs[2], translation);
1046 applyTranslation(baseNAUs[3], translation);
1047 applyTranslation(baseXYZ, translation);
1048
1049 centerOfMass(baseCoM[2], baseXYZ, massN, massSum, compareAtomsSize);
1050 if (useSym) {
1051 baseTransformSymOp = baseTransformSymOp.append(new SymOp(ZERO_ROTATION, translation));
1052 if (logger.isLoggable(Level.FINEST)) {
1053 stringBuilder.append(format("\n Base %d (%2d) 2nd Translation: ", l, bIndex))
1054 .append(Arrays.toString(translation)).append("\n");
1055 }
1056 }
1057
1058
1059 final double maxDist = baseAUDist[baseAUDist.length - 1].doubleValue();
1060 if (logger.isLoggable(Level.FINEST)) {
1061 stringBuilder.append(format("\n System 1 farthest distance: %4.8f", sqrt(maxDist)));
1062 }
1063 for (int i = 0; i < nAU; i++) {
1064 final int molIndex = baseAUDist_2[i].index() * nCoords;
1065 arraycopy(baseXYZ, molIndex, baseNAUs[3], i * nCoords, nCoords);
1066 }
1067 translation = calculateTranslation(baseNAUs[3], massN);
1068
1069 applyTranslation(baseNAUs[3], translation);
1070 applyTranslation(baseXYZ, translation);
1071 if (useSym) {
1072 baseTransformSymOp = baseTransformSymOp.append(new SymOp(ZERO_ROTATION, translation));
1073 if (logger.isLoggable(Level.FINEST)) {
1074 stringBuilder.append(format("\n Base %d (%2d) Final Translation: ", l, bIndex))
1075 .append(Arrays.toString(translation)).append("\n");
1076 }
1077 }
1078
1079 for (int m = 0; m < targetLength; m++) {
1080 final int tIndex = uniqueTarget[m];
1081 final int targetInd = targetAUDist[tIndex].index();
1082 final int currZ2 = targetDistMap.get(targetInd) % z2;
1083 if (strict || !reverse && baseTargetMap[m] == bIndex
1084 || reverse && baseTargetMap[l] == tIndex) {
1085
1086 arraycopy(targetXYZoriginal, 0, targetXYZ, 0, targetXYZoriginal.length);
1087
1088
1089 prioritizeReplicates(targetXYZ, targetCoM, compareAtomsSize,
1090 targetAUDist_2, targetInd, linkage);
1091
1092
1093 arraycopy(targetXYZ, targetAUDist_2[0].index() * nCoords, targetAU, 0, nCoords);
1094
1095
1096 translation = calculateTranslation(targetAU, massN);
1097
1098 applyTranslation(targetAU, translation);
1099
1100 applyTranslation(targetXYZ, translation);
1101
1102
1103 double[][] rotation = calculateRotation(baseAU, targetAU, massN);
1104 applyRotation(targetAU, rotation);
1105 applyRotation(targetXYZ, rotation);
1106
1107 centerOfMass(targetCoM, targetXYZ, massN, massSum, compareAtomsSize);
1108
1109 if (useSym) {
1110
1111 targetTransformSymOp = new SymOp(ZERO_ROTATION, Tr_0_0_0);
1112
1113 final SymOp target = targetSymOps.get(targetInd);
1114 targetSymOp = new SymOp(target.rot, target.tr);
1115
1116
1117 targetTransformSymOp = targetTransformSymOp.append(new SymOp(ZERO_ROTATION, translation))
1118 .append(new SymOp(rotation, Tr_0_0_0));
1119
1120 if (logger.isLoggable(Level.FINEST)) {
1121 stringBuilder.append(matrixToString(baseSymOp.rot, bIndex, "Base Sym Rot"));
1122 stringBuilder.append(format("\n Base %d Sym Op Translation: ", bIndex))
1123 .append(Arrays.toString(baseSymOp.tr)).append("\n");
1124 stringBuilder.append(matrixToString(targetSymOp.rot, tIndex, "Target Sym Rot"));
1125 stringBuilder.append(format("\n Target %d Sym Op Translation: ", tIndex))
1126 .append(Arrays.toString(targetSymOp.tr)).append("\n");
1127 stringBuilder.append(format("\n Trans target %d (%2d) sys to origin: \n\t", m, tIndex))
1128 .append(Arrays.toString(translation)).append("\n");
1129 stringBuilder.append(matrixToString(rotation, tIndex, "1st Rot"));
1130 }
1131 }
1132
1133
1134
1135
1136 pairEntities(max(3, nAU), 1);
1137
1138 double checkRMSD1 = -3.0;
1139 double n1RMSD = -4.0;
1140 if (logger.isLoggable(Level.FINE)) {
1141 checkRMSD1 = rmsd(baseAU, targetAU, massN);
1142 for (int i = 0; i < nAU; i++) {
1143 final int offset = i * nCoords;
1144 final int molIndex = pairedAUs[i].index() * nCoords;
1145 arraycopy(targetXYZ, molIndex, targetNAUs, offset, nCoords);
1146 }
1147 n1RMSD = rmsd(baseNAUs[1], targetNAUs, massN);
1148 if (logger.isLoggable(Level.FINEST) && useSave) {
1149 saveAssembly(file1, name1, bondList1, atoms1, forceField1, baseAU, comparisonAtoms, 1,
1150 "_c1_1", compNum, saveClusters);
1151 saveAssembly(file1, name1, bondList1, atoms1, forceField1, baseNAUs[1],
1152 comparisonAtoms, nAU, "_c1_1N", compNum, saveClusters);
1153 saveAssembly(file2, name2, bondList2, atoms2, forceField2, targetAU, comparisonAtoms,
1154 1, "_c2_1", compNum, saveClusters);
1155 saveAssembly(file2, name2, bondList2, atoms2, forceField2, targetNAUs, comparisonAtoms,
1156 nAU, "_c2_1N", compNum, saveClusters);
1157 }
1158 }
1159
1160
1161 for (int i = 0; i < 3; i++) {
1162 final int targetIndex = pairedAUs[i].index() * nCoords;
1163 final int coordIndex = i * nCoords;
1164 arraycopy(targetXYZ, targetIndex, target3AUs, coordIndex, nCoords);
1165 }
1166
1167 translation = calculateTranslation(target3AUs, massN);
1168 applyTranslation(target3AUs, translation);
1169 applyTranslation(targetNAUs, translation);
1170 applyTranslation(targetXYZ, translation);
1171
1172
1173
1174 rotation = calculateRotation(base3AUs, target3AUs, massN);
1175 applyRotation(target3AUs, rotation);
1176 applyRotation(targetNAUs, rotation);
1177 applyRotation(targetXYZ, rotation);
1178
1179 centerOfMass(targetCoM, targetXYZ, massN, massSum, compareAtomsSize);
1180 if (useSym) {
1181 applyTranslation(targetAU, translation);
1182 applyRotation(targetAU, rotation);
1183 targetTransformSymOp = targetTransformSymOp.append(new SymOp(ZERO_ROTATION, translation))
1184 .append(new SymOp(rotation, Tr_0_0_0));
1185 if (logger.isLoggable(Level.FINEST)) {
1186 stringBuilder.append(
1187 format("\n Target %d (%2d) 2nd Translation/Rotation: ", m, tIndex))
1188 .append(Arrays.toString(translation)).append("\n");
1189 printSym(compareAtomsSize, file2, name2, bondList2, atoms2, forceField2, saveClusters, currZ2);
1190 }
1191 }
1192
1193 double checkRMSD2 = -5.0;
1194 double n3RMSD = -6.0;
1195 if (logger.isLoggable(Level.FINE)) {
1196 checkRMSD2 = rmsd(base3AUs, target3AUs, massN);
1197 n3RMSD = rmsd(baseNAUs[2], targetNAUs, massN);
1198 if (logger.isLoggable(Level.FINEST) && useSave) {
1199 saveAssembly(file1, name1, bondList1, atoms1, forceField1, base3AUs, comparisonAtoms,
1200 3, "_c1_3", compNum, saveClusters);
1201 saveAssembly(file1, name1, bondList1, atoms1, forceField1, baseNAUs[2],
1202 comparisonAtoms, nAU, "_c1_3N", compNum, saveClusters);
1203 saveAssembly(file2, name2, bondList2, atoms2, forceField2, target3AUs, comparisonAtoms,
1204 3, "_c2_3", compNum, saveClusters);
1205 saveAssembly(file2, name2, bondList2, atoms2, forceField2, targetNAUs, comparisonAtoms,
1206 nAU, "_c2_3N", compNum, saveClusters);
1207 }
1208 }
1209
1210
1211
1212 pairEntities(nAU, 2);
1213
1214 for (int i = 0; i < nAU; i++) {
1215 final int offset = i * nCoords;
1216 final int molIndex = pairedAUs[i].index() * nCoords;
1217 arraycopy(targetXYZ, molIndex, targetNAUs, offset, nCoords);
1218 }
1219
1220 translation = calculateTranslation(targetNAUs, massN);
1221 applyTranslation(targetNAUs, translation);
1222 applyTranslation(targetXYZ, translation);
1223 rotation = calculateRotation(baseNAUs[3], targetNAUs, massN);
1224 if (useSym) {
1225 applyTranslation(targetAU, translation);
1226 applyRotation(targetAU, rotation);
1227 targetTransformSymOp = targetTransformSymOp.append(new SymOp(ZERO_ROTATION, translation))
1228 .append(new SymOp(rotation, Tr_0_0_0));
1229
1230 if (logger.isLoggable(Level.FINEST)) {
1231 stringBuilder.append(matrixToString(rotation, bIndex, "Target System Final Rotation"));
1232 stringBuilder.append(format("\n Target %d Final Translation: ", bIndex))
1233 .append(Arrays.toString(translation)).append("\n");
1234 printSym(compareAtomsSize, file2, name2, bondList2, atoms2, forceField2, saveClusters, currZ2);
1235 }
1236 }
1237 if (logger.isLoggable(Level.FINEST) && useSave) {
1238 saveAssembly(file1, name1, bondList1, atoms1, forceField1, baseNAUs[3], comparisonAtoms,
1239 nAU, "_c1_N", compNum, saveClusters);
1240 saveAssembly(file2, name2, bondList2, atoms2, forceField2, targetNAUs, comparisonAtoms,
1241 nAU, "_c2_N", compNum, saveClusters);
1242 }
1243 applyRotation(targetNAUs, rotation);
1244 applyRotation(targetXYZ, rotation);
1245 final double rmsdSymOp = rmsd(baseNAUs[3], targetNAUs, massN);
1246
1247 if (logger.isLoggable(Level.FINEST) && useSave) {
1248 saveAssembly(file1, name1, bondList1, atoms1, forceField1, baseNAUs[3], comparisonAtoms,
1249 nAU, "_c1_N2", compNum, saveClusters);
1250 saveAssembly(file2, name2, bondList2, atoms2, forceField2, targetNAUs, comparisonAtoms,
1251 nAU, "_c2_N2", compNum, saveClusters);
1252 }
1253
1254 final double baseGyration = radiusOfGyration(baseNAUs[3], massN);
1255 final double targetGyration = radiusOfGyration(targetNAUs, massN);
1256
1257 if (logger.isLoggable(Level.FINE)) {
1258 final int totalComparisons = (strict) ? baseLength * targetLength : maxLength;
1259 String output = format(" %2d of %2d: %7.4f (%8.4f) %7.4f (%8.4f) %8.4f %8.4f %8.4f",
1260 currentComparison, totalComparisons, checkRMSD1, n1RMSD, checkRMSD2, n3RMSD,
1261 rmsdSymOp, baseGyration, targetGyration);
1262
1263 if (logger.isLoggable(Level.FINER)) {
1264 if (reverse) {
1265 output += format(" b: %2d t: %2d bt: %2d", baseAUDist[bIndex].index(),
1266 targetAUDist[uniqueTarget[m]].index(), targetAUDist[baseTargetMap[m]].index());
1267 } else {
1268 output += format(" b: %2d t: %2d tb: %2d", baseAUDist[bIndex].index(),
1269 targetAUDist[uniqueTarget[m]].index(), baseAUDist[baseTargetMap[m]].index());
1270 }
1271 }
1272 stringBuilder.append(output).append("\n");
1273 }
1274
1275 if (useSym && rmsdSymOp < minZdiffs[currZ1][currZ2]) {
1276 minZdiffs[currZ1][currZ2] = rmsdSymOp;
1277 bestTargetTransformSymOp[currZ1][currZ2] = new SymOp(targetTransformSymOp.rot,
1278 targetTransformSymOp.tr);
1279 bestTargetSymOp[currZ1][currZ2] = new SymOp(targetSymOp.rot, targetSymOp.tr);
1280 bestBaseSymOp[currZ1][currZ2] = new SymOp(baseSymOp.rot, baseSymOp.tr);
1281 bestBaseTransformSymOp[currZ1][currZ2] = new SymOp(baseTransformSymOp.rot,
1282 baseTransformSymOp.tr);
1283 }
1284 if (rmsdSymOp < bestRMSD) {
1285
1286 gyrations[0] = baseGyration;
1287 gyrations[1] = targetGyration;
1288 bestRMSD = rmsdSymOp;
1289 baseBestZ = currZ1;
1290 targetBestZ = currZ2;
1291 arraycopy(baseNAUs[3], 0, bestBaseNAUs[currZ1][currZ2], 0, nAU * nCoords);
1292 arraycopy(targetNAUs, 0, bestTargetNAUs[currZ1][currZ2], 0, nAU * nCoords);
1293 }
1294 currentComparison++;
1295 }
1296 }
1297 }
1298
1299 double finalRMSD;
1300 if (bestRMSD < Double.MAX_VALUE) {
1301 finalRMSD = bestRMSD;
1302 } else {
1303 stringBuilder.append(" This RMSD was filtered out! Try the --st flag or increasing --if.\n");
1304
1305 finalRMSD = -7.0;
1306 }
1307 if (nAU == 1 && abs(min1 - finalRMSD) > MATCH_TOLERANCE) {
1308 stringBuilder.append(
1309 format(" WARNING: Single molecule overlay (%8.4f) does not match PAC RMSD_1 (%8.4f)", min1,
1310 finalRMSD));
1311 }
1312 if (useSave) {
1313 if (machineLearning) {
1314 saveAssembly(file1, name1, bondList1, atoms1, forceField1,
1315 bestBaseNAUs[baseBestZ][targetBestZ], comparisonAtoms, nAU, "_c1", 0.000, compNum, saveClusters);
1316 saveAssembly(file2, name2, bondList2, atoms2, forceField2,
1317 bestTargetNAUs[baseBestZ][targetBestZ], comparisonAtoms, nAU, "_c2", finalRMSD, compNum,
1318 saveClusters);
1319 } else {
1320 saveAssembly(file1, name1, bondList1, atoms1, forceField1,
1321 bestBaseNAUs[baseBestZ][targetBestZ], comparisonAtoms, nAU, "_c1", compNum, saveClusters);
1322 saveAssembly(file2, name2, bondList2, atoms2, forceField2,
1323 bestTargetNAUs[baseBestZ][targetBestZ], comparisonAtoms, nAU, "_c2", compNum, saveClusters);
1324 }
1325 }
1326
1327 if (inertia) {
1328 bestBaseMandV = momentsOfInertia(bestBaseNAUs[baseBestZ][targetBestZ], massN, false, false,
1329 true);
1330 bestTargetMandV = momentsOfInertia(bestTargetNAUs[baseBestZ][targetBestZ], massN, false, false,
1331 true);
1332 }
1333
1334 if (gyrationComponents) {
1335 bestBaseRg = radiusOfGyrationComponents(bestBaseNAUs[baseBestZ][targetBestZ], massN, true);
1336 bestTargetRg = radiusOfGyrationComponents(bestTargetNAUs[baseBestZ][targetBestZ], massN, true);
1337 }
1338
1339 return finalRMSD;
1340 }
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375 public RunningStatistics comparisons(final int nAU, final double inflationFactor, double matchTol, final double hitTol, final int zPrime,
1376 final int zPrime2, final String excludeAtomsA, final String excludeAtomsB, final boolean alphaCarbons,
1377 final boolean includeHydrogen, final boolean massWeighted, final int crystalPriority, final boolean strict,
1378 int saveClusters, final double save, final boolean restart, final boolean write, final boolean machineLearning, boolean inertia,
1379 final boolean gyrationComponents, int linkage, final double printSym, boolean lowMemory,
1380 final String pacFileName) {
1381 this.printSym = printSym;
1382
1383
1384
1385
1386 RunningStatistics runningStatistics;
1387 if (restart) {
1388 runningStatistics = readMatrix(pacFileName);
1389 if (runningStatistics == null) {
1390 runningStatistics = new RunningStatistics();
1391 }
1392 } else {
1393 runningStatistics = new RunningStatistics();
1394 File file = new File(pacFileName);
1395 if (file.exists() && file.delete()) {
1396 logger.info(format(" PAC RMSD file (%s) was deleted.", pacFileName));
1397 logger.info(" To restart from a previous run, use the '-r' flag.");
1398 }
1399 }
1400
1401
1402 if(save >= 0.0){
1403 lowMemory = true;
1404 }
1405
1406
1407
1408
1409
1410 for (int row = 0; row < restartRow; row++) {
1411 baseFilter.readNext(false, false, row + 1 >= restartRow);
1412 }
1413
1414 MolecularAssembly baseAssembly = baseFilter.getActiveMolecularSystem();
1415
1416
1417 Atom[] atoms1 = baseAssembly.getAtomArray();
1418 final int nAtoms = atoms1.length;
1419
1420
1421 MolecularAssembly targetAssembly = targetFilter.getActiveMolecularSystem();
1422 Atom[] atoms2 = targetAssembly.getAtomArray();
1423 final int nAtoms2 = atoms2.length;
1424
1425
1426 ArrayList<Integer> unique = new ArrayList<>(parseAtomRanges("Base Assembly", excludeAtomsA, nAtoms));
1427 if (invalidAtomSelection(unique, atoms1, alphaCarbons, includeHydrogen)) {
1428 logger.warning("\n No atoms were selected for the PAC RMSD in first crystal.");
1429 return null;
1430 }
1431 int[] comparisonAtoms = unique.stream().mapToInt(i -> i).toArray();
1432 List<MSNode> bondedEntities = baseAssembly.getAllBondedEntities();
1433
1434
1435 int z1;
1436 if(zPrime > 0){
1437 z1 = zPrime;
1438 }else{
1439 z1 = guessZPrime(unique, baseAssembly.getMoleculeNumbers(), bondedEntities.size());
1440 }
1441
1442
1443 unique = new ArrayList<>(parseAtomRanges("target", excludeAtomsB, nAtoms2));
1444 if (invalidAtomSelection(unique, atoms2, alphaCarbons, includeHydrogen)) {
1445 logger.warning("\n No atoms were selected for the PAC RMSD in second crystal.");
1446 return null;
1447 }
1448 int[] comparisonAtoms2 = unique.stream().mapToInt(i -> i).toArray();
1449
1450 List<MSNode> bondedEntities2 = targetAssembly.getAllBondedEntities();
1451 int z2;
1452 if(zPrime2 > 0){
1453 z2 = zPrime2;
1454 }else{
1455 logger.info(format(" Z2' info: %2d Num Entities %2d water.", bondedEntities.size(), baseAssembly.getWater().size()));
1456 z2 = guessZPrime(unique, targetAssembly.getMoleculeNumbers(), bondedEntities2.size());
1457 }
1458
1459
1460 int compareAtomsSize = comparisonAtoms.length;
1461 int compareAtomsSize2 = comparisonAtoms2.length;
1462
1463
1464
1465
1466 if (z1 > 1) {
1467 compareAtomsSize /= z1;
1468 }
1469 if (z2 > 1) {
1470 compareAtomsSize2 /= z2;
1471 }
1472 if(z1 != z2){
1473 logger.warning(format(" Z'1 (%2d) does not equal Z'2 (%2d).", z1, z2));
1474 }
1475
1476
1477 final boolean useSym = printSym >= 0.0;
1478 if (useSym && (z1 > 1 || z2 > 1)) {
1479 matchTol = SYM_TOLERANCE;
1480 }
1481 if (useSym && z1 != z2) {
1482 logger.warning(
1483 " Comparisons to determine symmetry should have the same number of molecules in the asymmetric unit.");
1484 }
1485
1486 Crystal baseCrystal = baseAssembly.getCrystal().getUnitCell();
1487 Crystal targetCrystal = targetFilter.getActiveMolecularSystem().getCrystal().getUnitCell();
1488
1489 int baseSearchValue = (baseCrystal.spaceGroup.respectsChirality()) ? z1 : 2 * z1;
1490 int targetSearchValue = (targetCrystal.spaceGroup.respectsChirality()) ? z2 : 2 * z2;
1491
1492
1493 final int nCoords = compareAtomsSize * 3;
1494
1495 if (this.comparisonAtoms == null) {
1496 this.comparisonAtoms = new int[compareAtomsSize];
1497 arraycopy(comparisonAtoms, 0, this.comparisonAtoms, 0, compareAtomsSize);
1498 }
1499
1500 int massIndex = 0;
1501 final double[] mass = new double[compareAtomsSize];
1502 for (Integer value : this.comparisonAtoms) {
1503 Atom atom = atoms1[value];
1504 double m = atom.getMass();
1505 mass[massIndex++] = (massWeighted) ? m : 1.0;
1506 }
1507
1508 massIndex = 0;
1509 final double[] mass2 = new double[compareAtomsSize2];
1510 for (Integer value : comparisonAtoms2) {
1511 if (massIndex == compareAtomsSize2) {
1512
1513 break;
1514 }
1515 Atom atom = atoms2[value];
1516 final double m = atom.getMass();
1517 mass2[massIndex++] = (massWeighted) ? m : 1.0;
1518 }
1519
1520 if (!Arrays.equals(mass, mass2)) {
1521 logger.warning(" Mass arrays do not match. Check atom ordering or size.");
1522 if (logger.isLoggable(Level.FINER)) {
1523 for (int i = 0; i < compareAtomsSize; i++) {
1524 if (mass[i] != mass2[i]) {
1525 logger.finer(format(" Index: %d Mass 1: %4.4f Mass 2: %4.4f", i, mass[i], mass2[i]));
1526 }
1527 }
1528 }
1529 }
1530
1531
1532 if (massN == null) {
1533 massN = (nAU > 3) ? new double[compareAtomsSize * nAU] : new double[compareAtomsSize * 3];
1534 for (int i = 0; i < nAU; i++) {
1535 arraycopy(mass, 0, massN, i * compareAtomsSize, compareAtomsSize);
1536 }
1537 }
1538
1539 if (massSum == 0) {
1540 for (int i = 0; i < compareAtomsSize; i++) {
1541 massSum += massN[i];
1542 }
1543 }
1544
1545 if (machineLearning) {
1546 saveClusters = 1;
1547 }
1548 if (saveClusters > 2) {
1549 saveClusters = 0;
1550 logger.info(" Save flag specified incorrectly (1:PDB; 2:XYZ). Not saving files.");
1551 }
1552
1553 if (logger.isLoggable(Level.FINER)) {
1554 if (linkage == 0) {
1555 logger.finer(" Single linkage will be used.");
1556 } else if (linkage == 2) {
1557 logger.finer(" Complete linkage will be used.");
1558 } else if (linkage == 1) {
1559 logger.finer(" Average linkage will be used.");
1560 }
1561 }
1562 if (linkage != 1 && linkage != 0 && linkage != 2) {
1563 logger.warning(
1564 "Prioritization method specified incorrectly (--pm {0, 1, 2}). Using default of average linkage.");
1565 linkage = 1;
1566 }
1567
1568
1569 if (compareAtomsSize == 0 || compareAtomsSize2 == 0) {
1570 logger.severe(" No atoms were selected for comparison.");
1571 return runningStatistics;
1572 }
1573 logger.info(format("\n %d atoms will be used for the PAC RMSD out of %d in first crystal.",
1574 compareAtomsSize * z1, nAtoms));
1575 logger.info(format(" %d atoms will be used for the PAC RMSD out of %d in second crystal.\n",
1576 compareAtomsSize2 * z2, nAtoms2));
1577
1578
1579 rmsdLabel = format("RMSD_%d", nAU);
1580
1581
1582 double minTime = Double.MAX_VALUE;
1583 double maxTime = -Double.MIN_VALUE;
1584
1585 if (!lowMemory) {
1586
1587 int size = (int) ceil((double) targetSize / (double) numProc);
1588 atomCache = new Atom[size][nAtoms2];
1589 fileCache = new File[size];
1590 nameCache = new String[size];
1591 bondCache = new ArrayList<>();
1592 for (int i = 0; i < size; i++) {
1593 bondCache.add(new ArrayList<>());
1594 }
1595 forceFieldCache = new ForceField[size];
1596 crystalCache = new Crystal[size];
1597
1598
1599 for (int column = restartColumn; column < targetSize; column++) {
1600 final int targetRank = column % numProc;
1601 if (targetRank == rank) {
1602 final int assemblyNum = column / numProc;
1603 MolecularAssembly currentAssembly = targetFilter.getActiveMolecularSystem();
1604 Atom[] arrayAtom = currentAssembly.getAtomArray();
1605 final int atomSize = arrayAtom.length;
1606 for (int i = 0; i < atomSize; i++) {
1607 final double[] xyz = new double[3];
1608 arrayAtom[i].getXYZ(xyz);
1609 atomCache[assemblyNum][i] = new Atom(i, arrayAtom[i].getName(),
1610 arrayAtom[i].getAtomType(), xyz);
1611 }
1612 List<Bond> currentBonds = currentAssembly.getBondList();
1613 List<Bond> currentBondCache = bondCache.get(assemblyNum);
1614 for (Bond b : currentBonds) {
1615 currentBondCache.add(new Bond(b.getAtom(0), b.getAtom(1)));
1616 }
1617 ForceField currentForcefield = currentAssembly.getForceField();
1618 fileCache[assemblyNum] = new File(currentAssembly.getFile().getName());
1619 forceFieldCache[assemblyNum] = new ForceField(currentForcefield.getProperties());
1620 nameCache[assemblyNum] = currentAssembly.getName();
1621 Crystal cXtal = currentAssembly.getCrystal().getUnitCell();
1622 crystalCache[assemblyNum] = new Crystal(cXtal.a, cXtal.b, cXtal.c, cXtal.alpha, cXtal.beta,
1623 cXtal.gamma, cXtal.spaceGroup.pdbName);
1624 }
1625 targetFilter.readNext(false, false, (column + 1) % numProc == rank);
1626 }
1627 }
1628
1629 if (logger.isLoggable(Level.FINER)) {
1630 logResources();
1631 }
1632
1633
1634 File saveLocation1 = new File(FilenameUtils.removeExtension(baseAssembly.getFile().getAbsoluteFile().getAbsolutePath()) + "_lt1_" + save + ".arc");
1635 File saveLocation2 = new File(FilenameUtils.removeExtension(targetAssembly.getFile().getAbsoluteFile().getAbsolutePath()) + "_lt2_" + save + ".arc");
1636 if(save >= 0.0) {
1637 logger.info(format(" Saving trajectories less than %7.4f to :\n %s\n %s", save, saveLocation1.getName(), saveLocation2.getName()));
1638 }
1639
1640 for (int row = restartRow; row < baseSize; row++) {
1641
1642 fill(myDistances, -8.0);
1643 int myIndex = 0;
1644
1645 baseAssembly = baseFilter.getActiveMolecularSystem();
1646 baseCrystal = baseAssembly.getCrystal().getUnitCell();
1647 atoms1 = baseAssembly.getAtomArray();
1648
1649 final double[] reducedBaseCoords = reduceSystem(atoms1, comparisonAtoms);
1650 final double baseDensity = baseCrystal.getUnitCell().getDensity(baseAssembly.getMass());
1651 if (baseCrystal == null || baseCrystal.aperiodic()) {
1652 logger.warning(" " + baseAssembly.getName() + " does not have a crystal. Consider using Superpose command.\n");
1653 continue;
1654 }
1655
1656
1657
1658
1659
1660 if (logger.isLoggable(Level.FINER)) {
1661 logger.finer(format(" Unit Cell Volume: (Base) %4.2f (Target) %4.2f", baseCrystal.volume,
1662 targetCrystal.volume));
1663 logger.finer(format(" Unit Cell Symm Ops: (Base) %d (Target) %d", baseCrystal.getNumSymOps(),
1664 targetCrystal.getNumSymOps()));
1665 logger.finer(format(" Z'': (Base) %d (Target) %d", z1, z2));
1666 }
1667
1668
1669
1670
1671
1672
1673
1674 final int[] baseLMN = new int[3];
1675 double inflationFactorOrig = inflationFactor;
1676 baseXYZoriginal = generateInflatedSphere(baseCrystal.getUnitCell(), reducedBaseCoords, z1, nAU,
1677 linkage, mass, baseLMN, strict, baseSymOps, baseDistMap, inflationFactor);
1678 if (inflationFactorOrig != inflationFactor) {
1679 logger.warning(format(
1680 " Default replicates crystal was too small for comparison. Increasing inflation factor from %9.3f to %9.3f",
1681 inflationFactor, inflationFactorOrig));
1682 }
1683 int nBaseCoords = baseXYZoriginal.length;
1684 baseXYZ = new double[nBaseCoords];
1685 arraycopy(baseXYZoriginal, 0, baseXYZ, 0, nBaseCoords);
1686
1687
1688 int nBaseMols = nBaseCoords / (compareAtomsSize * 3);
1689 if (baseAUDist == null || baseAUDist.length != nBaseMols) {
1690 baseAUDist = new DoubleIndexPair[nBaseMols];
1691 }
1692 if (baseCoM == null || baseCoM.length != nBaseMols) {
1693 baseCoM = new double[3][nBaseMols][3];
1694 }
1695 centerOfMass(baseCoM[0], baseXYZ, massN, massSum, compareAtomsSize);
1696 prioritizeReplicates(baseXYZ, baseCoM[0], compareAtomsSize, baseAUDist, 0,
1697 linkage);
1698
1699 for (int column = restartColumn; column < targetSize; column++) {
1700 final int targetRank = column % numProc;
1701 if (targetRank == rank) {
1702 stringBuilder.setLength(0);
1703 long time = -System.nanoTime();
1704 final int assemblyNum = column / numProc;
1705 if (!lowMemory) {
1706 targetCrystal = crystalCache[assemblyNum];
1707 atoms2 = atomCache[assemblyNum];
1708 } else {
1709 targetAssembly = targetFilter.getActiveMolecularSystem();
1710 targetCrystal = targetAssembly.getCrystal().getUnitCell();
1711 atoms2 = targetAssembly.getAtomArray();
1712 }
1713 if (targetCrystal == null || targetCrystal.aperiodic()) {
1714 if (!lowMemory) {
1715 logger.warning(" " + nameCache[assemblyNum] + " does not have a crystal. Consider using Superpose command.\n");
1716 } else {
1717 logger.warning(" " + targetAssembly.getName() + " does not have a crystal. Consider using Superpose command.\n");
1718 }
1719 continue;
1720 }
1721 double rmsd;
1722
1723 boolean densityCheck;
1724 if (isSymmetric && row == column) {
1725 stringBuilder.append(
1726 format("\n Comparing Model %4d (%s) of %s\n with Model %4d (%s) of %s\n",
1727 row + 1, baseCrystal.toShortString(), baseLabel, column + 1,
1728 targetCrystal.toShortString(), targetLabel));
1729
1730 rmsd = 0.0;
1731
1732 stringBuilder.append(format("\n PAC %s: %12s %7.4f A\n", rmsdLabel, "", rmsd));
1733 } else if (isSymmetric && row > column) {
1734
1735 rmsd = -10.0;
1736 } else {
1737 stringBuilder.append(
1738 format("\n Comparing Model %4d (%s) of %s\n with Model %4d (%s) of %s\n",
1739 row + 1, baseCrystal.toShortString(), baseLabel, column + 1,
1740 targetCrystal.toShortString(), targetLabel));
1741
1742
1743 final double[] reducedTargetCoords = reduceSystem(atoms2, comparisonAtoms2);
1744
1745
1746
1747 final int[] targetLMN = new int[3];
1748 inflationFactorOrig = inflationFactor;
1749 targetXYZoriginal = generateInflatedSphere(targetCrystal.getUnitCell(),
1750 reducedTargetCoords, z2, nAU, linkage, mass2, targetLMN, strict, targetSymOps,
1751 targetDistMap, inflationFactor);
1752 if (inflationFactorOrig != inflationFactor) {
1753 logger.warning(format(
1754 " Default replicates crystal was too small for comparison. Increasing inflation factor from %9.3f to %9.3f",
1755 inflationFactor, inflationFactorOrig));
1756 }
1757 int nTargetCoords = targetXYZoriginal.length;
1758 targetXYZ = new double[nTargetCoords];
1759 arraycopy(targetXYZoriginal, 0, targetXYZ, 0, nTargetCoords);
1760
1761
1762 double densityMass = 0;
1763 if (!lowMemory) {
1764 for (Atom atom : atomCache[assemblyNum]) {
1765 densityMass += atom.getMass();
1766 }
1767 } else {
1768 densityMass = targetFilter.getActiveMolecularSystem().getMass();
1769 }
1770 final double targetDensity = targetCrystal.getUnitCell().getDensity(densityMass);
1771 if (logger.isLoggable(Level.FINER)) {
1772 stringBuilder.append(
1773 format("\n Base Density: %4.4f Target Density: %4.4f\n", baseDensity,
1774 targetDensity));
1775 }
1776
1777 densityCheck =
1778 (crystalPriority == 1) ? baseDensity < targetDensity : baseDensity > targetDensity;
1779
1780
1781 int nTargetMols = targetXYZ.length / (compareAtomsSize2 * 3);
1782 if (targetAUDist == null || targetAUDist.length != nTargetMols) {
1783 targetAUDist = new DoubleIndexPair[nTargetMols];
1784 }
1785 if (targetCoM == null || targetCoM.length != nTargetMols) {
1786 targetCoM = new double[nTargetMols][3];
1787 }
1788
1789 centerOfMass(targetCoM, targetXYZ, massN, massSum, compareAtomsSize);
1790 prioritizeReplicates(targetXYZ, targetCoM, compareAtomsSize,
1791 targetAUDist, 0, linkage);
1792
1793
1794 File file1;
1795 File file2;
1796 String name1;
1797 String name2;
1798 List<Bond> bondList1;
1799 List<Bond> bondList2;
1800 ForceField forceField1;
1801 ForceField forceField2;
1802 arraycopy(baseXYZoriginal, 0, baseXYZ, 0, nBaseCoords);
1803 if (densityCheck || crystalPriority == 2) {
1804
1805 atoms1 = baseAssembly.getAtomArray();
1806 file1 = baseAssembly.getFile();
1807 name1 = baseAssembly.getName();
1808 bondList1 = baseAssembly.getBondList();
1809 forceField1 = baseAssembly.getForceField();
1810 if (!lowMemory) {
1811 atoms2 = atomCache[assemblyNum];
1812 file2 = fileCache[assemblyNum];
1813 name2 = nameCache[assemblyNum];
1814 bondList2 = bondCache.get(assemblyNum);
1815 forceField2 = forceFieldCache[assemblyNum];
1816 } else {
1817 MolecularAssembly mobileAssembly = targetFilter.getActiveMolecularSystem();
1818 atoms2 = mobileAssembly.getAtomArray();
1819 file2 = mobileAssembly.getFile();
1820 name2 = mobileAssembly.getName();
1821 bondList2 = mobileAssembly.getBondList();
1822 forceField2 = mobileAssembly.getForceField();
1823 }
1824 } else {
1825
1826 double[] tempXYZ = new double[nBaseCoords];
1827 arraycopy(baseXYZ, 0, tempXYZ, 0, nBaseCoords);
1828 baseXYZ = new double[nTargetCoords];
1829 arraycopy(targetXYZ, 0, baseXYZ, 0, nTargetCoords);
1830 targetXYZ = new double[nBaseCoords];
1831 arraycopy(tempXYZ, 0, targetXYZ, 0, nBaseCoords);
1832 arraycopy(baseXYZoriginal, 0, tempXYZ, 0, nBaseCoords);
1833 baseXYZoriginal = new double[nTargetCoords];
1834 arraycopy(targetXYZoriginal, 0, baseXYZoriginal, 0, nTargetCoords);
1835 targetXYZoriginal = new double[nBaseCoords];
1836 arraycopy(tempXYZ, 0, targetXYZoriginal, 0, nBaseCoords);
1837 tempXYZ = new double[nCoords];
1838 arraycopy(reducedBaseCoords, 0, tempXYZ, 0, nCoords);
1839 arraycopy(reducedTargetCoords, 0, reducedBaseCoords, 0, nCoords);
1840 arraycopy(tempXYZ, 0, reducedTargetCoords, 0, nCoords);
1841 ArrayList<SymOp> also = new ArrayList<>(baseSymOps);
1842 baseSymOps = new ArrayList<>(targetSymOps);
1843 targetSymOps = new ArrayList<>(also);
1844 ArrayList<Integer> ali = new ArrayList<>(baseDistMap);
1845 baseDistMap = new ArrayList<>(targetDistMap);
1846 targetDistMap = new ArrayList<>(ali);
1847 DoubleIndexPair[] tempDIP = new DoubleIndexPair[nBaseMols];
1848 arraycopy(baseAUDist, 0, tempDIP, 0, nBaseMols);
1849 baseAUDist = new DoubleIndexPair[nTargetMols];
1850 arraycopy(targetAUDist, 0, baseAUDist, 0, nTargetMols);
1851 targetAUDist = new DoubleIndexPair[nBaseMols];
1852 arraycopy(tempDIP, 0, targetAUDist, 0, nBaseMols);
1853 double[][] tempDD = new double[nBaseMols][3];
1854 arraycopy(baseCoM[0], 0, tempDD, 0, nBaseMols);
1855 baseCoM = new double[3][nTargetMols][3];
1856 arraycopy(targetCoM, 0, baseCoM[0], 0, nTargetMols);
1857 targetCoM = new double[nBaseMols][3];
1858 arraycopy(tempDD, 0, targetCoM, 0, nBaseMols);
1859 atoms2 = baseAssembly.getAtomArray();
1860 file2 = baseAssembly.getFile();
1861 name2 = baseAssembly.getName();
1862 bondList2 = baseAssembly.getBondList();
1863 forceField2 = baseAssembly.getForceField();
1864 int[] tempAtoms = comparisonAtoms.clone();
1865 comparisonAtoms = comparisonAtoms2.clone();
1866 comparisonAtoms2 = tempAtoms.clone();
1867 int temp = z1;
1868 z1 = z2;
1869 z2 = temp;
1870 temp = baseSearchValue;
1871 baseSearchValue = targetSearchValue;
1872 targetSearchValue = temp;
1873 temp = nBaseCoords;
1874 nBaseCoords = nTargetCoords;
1875 nTargetCoords = temp;
1876 temp = nBaseMols;
1877 nBaseMols = nTargetMols;
1878 nTargetMols = temp;
1879 if (!lowMemory) {
1880 atoms1 = atomCache[assemblyNum];
1881 file1 = fileCache[assemblyNum];
1882 name1 = nameCache[assemblyNum];
1883 bondList1 = bondCache.get(assemblyNum);
1884 forceField1 = forceFieldCache[assemblyNum];
1885 } else {
1886 MolecularAssembly staticAssembly = targetFilter.getActiveMolecularSystem();
1887 atoms1 = staticAssembly.getAtomArray();
1888 file1 = staticAssembly.getFile();
1889 name1 = staticAssembly.getName();
1890 bondList1 = staticAssembly.getBondList();
1891 forceField1 = staticAssembly.getForceField();
1892 }
1893 }
1894
1895 if (baseAU == null) {
1896 baseAU = new double[nCoords];
1897 }
1898 if (targetAU == null) {
1899 targetAU = new double[nCoords];
1900 }
1901 if (pairedAUs == null || pairedAUs.length != max(3, nAU)) {
1902 pairedAUs = new DoubleIndexPair[max(3, nAU)];
1903 }
1904 if (base3AUs == null) {
1905 base3AUs = new double[nCoords * 3];
1906 }
1907 if (target3AUs == null) {
1908 target3AUs = new double[nCoords * 3];
1909 }
1910 if (baseNAUs == null) {
1911 baseNAUs = new double[4][nAU * nCoords];
1912 }
1913 if (targetNAUs == null) {
1914 targetNAUs = new double[nAU * nCoords];
1915 }
1916 if (bestBaseNAUs == null) {
1917 bestBaseNAUs = new double[z1][z2][nAU * nCoords];
1918 }
1919 if (bestTargetNAUs == null) {
1920 bestTargetNAUs = new double[z1][z2][nAU * nCoords];
1921 }
1922 if (useSym) {
1923 if (baseSymOps == null) {
1924 baseSymOps = new ArrayList<>();
1925 }
1926 if (targetSymOps == null) {
1927 targetSymOps = new ArrayList<>();
1928 }
1929 if (baseAUoriginal == null || baseAUoriginal.length != z1) {
1930 baseAUoriginal = new double[z1][nCoords];
1931 }
1932 if (targetAUoriginal == null || targetAUoriginal.length != z2) {
1933 targetAUoriginal = new double[z2][nCoords];
1934 }
1935 }
1936
1937 if (useSym) {
1938 for (int i = 0; i < z1; i++) {
1939 arraycopy(reducedBaseCoords, i * compareAtomsSize * 3, baseAUoriginal[i], 0,
1940 compareAtomsSize * 3);
1941 }
1942 for (int i = 0; i < z2; i++) {
1943 arraycopy(reducedTargetCoords, i * compareAtomsSize2 * 3, targetAUoriginal[i], 0,
1944 compareAtomsSize2 * 3);
1945 }
1946 }
1947
1948 if (logger.isLoggable(Level.FINER)) {
1949 logger.finer(
1950 format(" %3d Molecules in Base Crystal \t %3d Molecules in Target Crystal",
1951 (nBaseCoords / 3) / compareAtomsSize, (nTargetCoords / 3) / compareAtomsSize));
1952 }
1953
1954 if (logger.isLoggable(Level.FINEST) && saveClusters > 0) {
1955 saveAssembly(file1, name1, bondList1, atoms1, forceField1, baseXYZoriginal,
1956 comparisonAtoms, nBaseCoords / 3, "_c1_rep", -1, saveClusters);
1957 saveAssembly(file2, name2, bondList2, atoms2, forceField2, targetXYZoriginal,
1958 comparisonAtoms2, nTargetCoords / 3, "_c2_rep", -1, saveClusters);
1959 }
1960
1961 rmsd = compare(file1, name1, bondList1, atoms1, forceField1, file2, name2, bondList2,
1962 atoms2, forceField2, z1, z2, compareAtomsSize, nAU, baseSearchValue,
1963 targetSearchValue, matchTol, row * targetSize + column, strict, saveClusters,
1964 machineLearning, linkage, inertia, gyrationComponents);
1965 time += System.nanoTime();
1966 final double timeSec = time * 1.0e-9;
1967
1968 if (timeSec < minTime) {
1969 minTime = timeSec;
1970 }
1971 if (timeSec > maxTime) {
1972 maxTime = timeSec;
1973 }
1974
1975 final double avgGyration = (gyrations[0] + gyrations[1]) / 2;
1976 final double diff = max(gyrations[0], gyrations[1]) - avgGyration;
1977 stringBuilder.append(
1978 format("\n PAC %7s: %12s %7.4f A (%5.3f sec) G(r) %7.4f A +/- %7.4f\n", rmsdLabel,
1979 "", rmsd, timeSec, avgGyration, diff));
1980 if (inertia) {
1981 stringBuilder.append("""
1982
1983 Moments of Inertia and Principle Axes for first crystal:
1984 Moments (amu Ang^2): \t X-, Y-, and Z-Components of Axes:
1985 """);
1986 for (int i = 1; i < 4; i++) {
1987 stringBuilder.append(
1988 format(" %16.3f %12.6f %12.6f %12.6f\n", bestBaseMandV[i - 1][0],
1989 bestBaseMandV[0][i], bestBaseMandV[1][i], bestBaseMandV[2][i]));
1990 }
1991 stringBuilder.append("""
1992 Moments of Inertia and Principle Axes for second crystal:
1993 Moments (amu Ang^2): \t X-, Y-, and Z-Components of Axes:
1994 """);
1995 for (int i = 1; i < 4; i++) {
1996 stringBuilder.append(
1997 format(" %16.3f %12.6f %12.6f %12.6f\n", bestTargetMandV[i - 1][0],
1998 bestTargetMandV[0][i], bestTargetMandV[1][i], bestTargetMandV[2][i]));
1999 }
2000 }
2001 if (gyrationComponents) {
2002
2003 double Rmax = max(max(bestBaseRg[0], bestBaseRg[1]), bestBaseRg[2]);
2004 Rmax *= Rmax;
2005 double Rmed = max(min(bestBaseRg[0], bestBaseRg[1]), bestBaseRg[2]);
2006 Rmed *= Rmed;
2007 double Rmin = min(min(bestBaseRg[0], bestBaseRg[1]), bestBaseRg[2]);
2008 Rmin *= Rmin;
2009 double Rg = Rmax + Rmed + Rmin;
2010
2011 double asphericity = Rmax - 0.5 * (Rmin + Rmed);
2012 double acylindricity = Rmed - Rmin;
2013 double anisotropy = (asphericity * asphericity + 0.75 * acylindricity * acylindricity);
2014 double totalMass = Arrays.stream(mass).sum();
2015
2016 double volume = totalMass / max(baseCrystal.getDensity(massSum),
2017 targetCrystal.getDensity(massSum));
2018 double targetRadius = cbrt(0.75 / PI * volume);
2019 stringBuilder.append(format(" Base Radius Metric: %7.4f", avgGyration / targetRadius));
2020 if (logger.isLoggable(Level.FINER)) {
2021 stringBuilder.append(
2022 format("\n Target Base Radius: %9.3f Ã…\n", targetRadius));
2023 stringBuilder.append(
2024 format(" Asphericity: (%6.3f Ã…^2) %9.3f\n", asphericity, asphericity / Rmax));
2025 stringBuilder.append(format(" Acylindricity: (%6.3f Ã…^2) %9.3f\n", acylindricity,
2026 acylindricity / Rmax));
2027 stringBuilder.append(
2028 format(" Anisotropy: (%6.3f Ã…) %9.3f\n", sqrt(sqrt(anisotropy)),
2029 anisotropy / (Rg * Rg)));
2030 stringBuilder.append(
2031 format(" Ryz %9.3fÃ… Rxz %9.3fÃ… Rxy %9.3fÃ…\n", bestBaseRg[0], bestBaseRg[1],
2032 bestBaseRg[2]));
2033 }
2034
2035 Rmax = max(max(bestTargetRg[0], bestTargetRg[1]), bestTargetRg[2]);
2036 Rmax *= Rmax;
2037 Rmed = max(min(bestTargetRg[0], bestTargetRg[1]), bestTargetRg[2]);
2038 Rmed *= Rmed;
2039 Rmin = min(min(bestTargetRg[0], bestTargetRg[1]), bestTargetRg[2]);
2040 Rmin *= Rmin;
2041 Rg = Rmax + Rmed + Rmin;
2042
2043 asphericity = Rmax - 0.5 * (Rmin + Rmed);
2044 acylindricity = Rmed - Rmin;
2045 anisotropy = (asphericity * asphericity + 0.75 * acylindricity * acylindricity);
2046 totalMass = Arrays.stream(mass).sum();
2047 volume = totalMass / targetDensity;
2048 targetRadius = cbrt(0.75 / PI * volume);
2049 if (logger.isLoggable(Level.FINER)) {
2050 stringBuilder.append(
2051 format("\n Target Target Radius: %9.3f Ã…\n", targetRadius));
2052 stringBuilder.append(
2053 format(" Asphericity: (%6.3f Ã…) %9.3f\n", asphericity, asphericity / Rmax));
2054 stringBuilder.append(format(" Acylindricity: (%6.3f Ã…) %9.3f\n", acylindricity,
2055 acylindricity / Rmax));
2056 stringBuilder.append(
2057 format(" Anisotropy: (%6.3f Ã…) %9.3f\n", sqrt(sqrt(anisotropy)),
2058 anisotropy / (Rg * Rg)));
2059 stringBuilder.append(
2060 format(" Ryz %9.3fÃ… Rxz %9.3fÃ… Rxy %9.3fÃ…\n", bestTargetRg[0],
2061 bestTargetRg[1], bestTargetRg[2]));
2062 }
2063 }
2064 if (logger.isLoggable(Level.FINER)) {
2065 stringBuilder.append(
2066 format("\n Gyration Crystal 1 (%s): %7.4f Crystal 2 (%s): %7.4f.\n", name1,
2067 gyrations[0], name2, gyrations[1]));
2068 }
2069 if (useSym) {
2070 StringBuilder sbSO = new StringBuilder(
2071 format("\n Sym Op to move %s onto %s:\nsymop ", name2, name1));
2072 StringBuilder sbInv = new StringBuilder(
2073 format("\n Inverted Sym Op to move %s onto %s:\nsymop ", name1, name2));
2074 final int[] mol1 = baseAssembly.getMoleculeNumbers();
2075 final int[] mol2 = targetAssembly.getMoleculeNumbers();
2076 for (int i = 0; i < z1; i++) {
2077 for (int j = 0; j < z2; j++) {
2078 if (bestTargetTransformSymOp[i][j] != null) {
2079 bestTargetTransformSymOp[i][j] = bestTargetTransformSymOp[i][j].prepend(
2080 bestTargetSymOp[i][j]);
2081
2082
2083
2084 bestBaseTransformSymOp[i][j] = bestBaseTransformSymOp[i][j].prepend(
2085 bestBaseSymOp[i][j]);
2086 bestTargetTransformSymOp[i][j] = bestTargetTransformSymOp[i][j].append(
2087 invertSymOp(bestBaseTransformSymOp[i][j]));
2088
2089 if (strict || logger.isLoggable(Level.FINEST) || z1 != z2 || i == j) {
2090 SymOp inverted = invertSymOp(bestTargetTransformSymOp[i][j]);
2091 if (logger.isLoggable(Level.FINEST)) {
2092 stringBuilder.append(
2093 format("\n Sym Op to move %s (%2d) onto %s (%2d):\n", name2, j, name1,
2094 i)).append(bestTargetTransformSymOp[i][j]).append("\n");
2095 stringBuilder.append(
2096 format("\n\n Inverted Sym Op to move %s (%2d) onto %s (%2d):\n", name1,
2097 i, name2, j)).append(inverted).append("\n");
2098 }
2099 ArrayList<Integer> mol1list = new ArrayList<>();
2100 ArrayList<Integer> mol2list = new ArrayList<>();
2101 for (int k = 0; k < nAtoms; k++) {
2102 if ((z1 == 1 || i == mol1[k]) && ArrayUtils.contains(comparisonAtoms, k)) {
2103 mol1list.add(k);
2104 }
2105 if ((z2 == 1 || j == mol2[k]) && ArrayUtils.contains(comparisonAtoms2, k)) {
2106 mol2list.add(k);
2107 }
2108 }
2109 final int[] mol1arr = mol1list.stream().mapToInt(Integer::intValue).toArray();
2110 final int[] mol2arr = mol2list.stream().mapToInt(Integer::intValue).toArray();
2111 sbSO.append(format(" %s %s", writeAtomRanges(mol1arr), writeAtomRanges(mol2arr)))
2112 .append(asMatrixString(bestTargetTransformSymOp[i][j]));
2113 sbInv.append(format(" %s %s", writeAtomRanges(mol2arr), writeAtomRanges(mol1arr))).append(asMatrixString(inverted));
2114 if (i + 1 < z1 || j + 1 < z2) {
2115 sbSO.append("\\\n");
2116 sbInv.append("\\\n");
2117 }
2118
2119 if (logger.isLoggable(Level.FINER)) {
2120 stringBuilder.append(format(
2121 "\n Sym Op Application from Scratch: Target: %s (%2d) onto Base: %s (%2d)",
2122 name2, j, name1, i));
2123 final double[] symMol = new double[nCoords];
2124 for (int k = 0; k < compareAtomsSize; k++) {
2125 final int l = k * 3;
2126 final double[] xyz = new double[]{targetAUoriginal[j][l],
2127 targetAUoriginal[j][l + 1], targetAUoriginal[j][l + 2]};
2128 applyCartesianSymOp(xyz, xyz, bestTargetTransformSymOp[i][j]);
2129 symMol[l] = xyz[0];
2130 symMol[l + 1] = xyz[1];
2131 symMol[l + 2] = xyz[2];
2132 stringBuilder.append(format(
2133 "\n %8.3f %8.3f %8.3f moved to %8.3f %8.3f %8.3f compared to %8.3f %8.3f %8.3f",
2134 targetAUoriginal[j][l], targetAUoriginal[j][l + 1],
2135 targetAUoriginal[j][l + 2], symMol[l], symMol[l + 1], symMol[l + 2],
2136 baseAUoriginal[i][l], baseAUoriginal[i][l + 1],
2137 baseAUoriginal[i][l + 2]));
2138 }
2139 if (saveClusters > 0) {
2140 saveAssembly(file2, name2, bondList2, atoms2, forceField2, symMol,
2141 comparisonAtoms, 1, "_moved", 0, saveClusters);
2142 }
2143
2144 final double value = rmsd(baseAUoriginal[i], symMol, massN);
2145 final double[] symMol2 = new double[nCoords];
2146 stringBuilder.append("\n\n Sym Op Inverse Application:");
2147 for (int k = 0; k < compareAtomsSize; k++) {
2148 final int l = k * 3;
2149 final double[] xyz = new double[]{symMol[l], symMol[l + 1], symMol[l + 2]};
2150 applyCartesianSymOp(xyz, xyz, inverted);
2151 symMol2[l] = xyz[0];
2152 symMol2[l + 1] = xyz[1];
2153 symMol2[l + 2] = xyz[2];
2154 stringBuilder.append(format(
2155 "\n %8.3f %8.3f %8.3f moved to %8.3f %8.3f %8.3f compared to %8.3f %8.3f %8.3f",
2156 symMol[l], symMol[l + 1], symMol[l + 2], symMol2[l], symMol2[l + 1],
2157 symMol2[l + 2], targetAUoriginal[j][l], targetAUoriginal[j][l + 1],
2158 targetAUoriginal[j][l + 2]));
2159 }
2160 stringBuilder.append(format("\n\n Sym Op RMSD: %8.8f \n\n", value));
2161 }
2162 }
2163 } else {
2164 if (logger.isLoggable(Level.INFO)) {
2165 logger.info(format(" Symmetry Operator %2d %2d was filtered out. "
2166 + "Try using a single asymmetric unit (--na 1).", i, j));
2167 }
2168 }
2169 }
2170 }
2171 stringBuilder.append(sbSO.toString()).append("\n").append(sbInv.toString())
2172 .append("\n");
2173 }
2174
2175 if (!densityCheck && crystalPriority != 2) {
2176 baseXYZ = new double[nTargetCoords];
2177 arraycopy(targetXYZ, 0, baseXYZ, 0, nTargetCoords);
2178 baseXYZoriginal = new double[nTargetCoords];
2179 arraycopy(targetXYZoriginal, 0, baseXYZoriginal, 0, nTargetCoords);
2180 arraycopy(reducedTargetCoords, 0, reducedBaseCoords, 0, nCoords);
2181 baseSymOps = new ArrayList<>(targetSymOps);
2182 baseDistMap = new ArrayList<>(targetDistMap);
2183 int[] tempAtoms = comparisonAtoms.clone();
2184 comparisonAtoms = comparisonAtoms2.clone();
2185 comparisonAtoms2 = tempAtoms.clone();
2186 DoubleIndexPair[] tempDIP = new DoubleIndexPair[nBaseMols];
2187 arraycopy(baseAUDist, 0, tempDIP, 0, nBaseMols);
2188 baseAUDist = new DoubleIndexPair[nTargetMols];
2189 arraycopy(targetAUDist, 0, baseAUDist, 0, nTargetMols);
2190 targetAUDist = new DoubleIndexPair[nBaseMols];
2191 arraycopy(tempDIP, 0, targetAUDist, 0, nBaseMols);
2192 final double[][] tempDD = new double[nBaseMols][3];
2193 arraycopy(baseCoM[0], 0, tempDD, 0, nBaseMols);
2194 baseCoM = new double[3][nTargetMols][3];
2195 arraycopy(targetCoM, 0, baseCoM[0], 0, nTargetMols);
2196 targetCoM = new double[nBaseMols][3];
2197 arraycopy(tempDD, 0, targetCoM, 0, nBaseMols);
2198 int temp = z1;
2199 z1 = z2;
2200 z2 = temp;
2201 temp = baseSearchValue;
2202 baseSearchValue = targetSearchValue;
2203 targetSearchValue = temp;
2204 nBaseCoords = nTargetCoords;
2205 nBaseMols = nTargetMols;
2206 }
2207 if (hitTol >= 0.0 && hitTol > rmsd) {
2208 numberOfHits++;
2209 }
2210 }
2211 myDistances[myIndex] = rmsd;
2212 myIndex++;
2213 if (!stringBuilder.isEmpty()) {
2214 logger.info(stringBuilder.toString());
2215 }
2216 if (save >= 0.0 && rmsd < save) {
2217 XYZFilter xyzFilter = new XYZFilter(saveLocation1.getAbsoluteFile(), baseAssembly, baseAssembly.getForceField(),
2218 baseAssembly.getProperties());
2219 xyzFilter.writeFile(saveLocation1.getAbsoluteFile(), true);
2220 XYZFilter xyzFilter2 = new XYZFilter(saveLocation2.getAbsoluteFile(), targetAssembly, targetAssembly.getForceField(),
2221 targetAssembly.getProperties());
2222 xyzFilter2.writeFile(saveLocation2.getAbsoluteFile(), true);
2223 }
2224 }
2225 if (lowMemory) {
2226 targetFilter.readNext(false, false, (column + 1) % numProc == rank);
2227 }
2228 }
2229 restartColumn = 0;
2230 if (lowMemory) {
2231 targetFilter.readNext(true, false, 0 == rank);
2232 }
2233 baseFilter.readNext(false, false);
2234
2235
2236 gatherRMSDs(row, runningStatistics);
2237
2238
2239 if (rank == 0 && write) {
2240 final int firstColumn = (isSymmetric) ? row : 0;
2241 writeDistanceMatrixRow(pacFileName, distRow, firstColumn);
2242 }
2243 }
2244
2245 if (minTime < Double.MAX_VALUE) {
2246 logger.info(format("\n Minimum PAC Time: %7.4f", minTime));
2247 }
2248 if (logger.isLoggable(Level.FINE) && maxTime > Double.MIN_VALUE && maxTime != minTime) {
2249 logger.info(format(" Maximum PAC Time: %7.4f", maxTime));
2250 }
2251
2252 baseFilter.closeReader();
2253 targetFilter.closeReader();
2254
2255 if (rank == 0 || logger.isLoggable(Level.FINER)) {
2256 final double min = runningStatistics.getMin();
2257 if (Double.MAX_VALUE - min < matchTol) {
2258 logger.warning(" SuperposeCrystals was not able to compare the provided structures.");
2259 if (logger.isLoggable(Level.FINE)) {
2260 logger.info(format(" RMSD Minimum: %8.6f", min));
2261 }
2262 } else {
2263 if(hitTol >= 0.0){
2264 logger.info(format(" Number of \"Hits\": %8d", numberOfHits));
2265 }
2266 logger.info(format(" RMSD Minimum: %8.6f", min));
2267 logger.info(format(" RMSD Maximum: %8.6f", runningStatistics.getMax()));
2268 logger.info(format(" RMSD Mean: %8.6f", runningStatistics.getMean()));
2269 final double variance = runningStatistics.getVariance();
2270 if (!Double.isNaN(variance)) {
2271 logger.info(format(" RMSD Variance: %8.6f", variance));
2272 }
2273 }
2274 }
2275
2276
2277 return runningStatistics;
2278 }
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288 private static boolean addLooseUnequal(List<Double> values, final double value, final double tol) {
2289 boolean found = false;
2290 for (Double dbl : values) {
2291 if (abs(dbl - value) < tol) {
2292 found = true;
2293 break;
2294 }
2295 }
2296 if (!found) {
2297 values.add(value);
2298 }
2299
2300 return !found;
2301 }
2302
2303
2304
2305
2306
2307
2308
2309
2310 public void addRotation(final double[][] rot, final double[][] totalTransform, final boolean prepend) {
2311 double[][] transform = new double[][]{{rot[0][0], rot[0][1], rot[0][2], 0.0},
2312 {rot[1][0], rot[1][1], rot[1][2], 0.0}, {rot[2][0], rot[2][1], rot[2][2], 0.0},
2313 {0.0, 0.0, 0.0, 1.0}};
2314 transform =
2315 (prepend) ? mat4Mat4(totalTransform, transform) : mat4Mat4(transform, totalTransform);
2316 final int length = totalTransform.length;
2317 for (int i = 0; i < length; i++) {
2318 arraycopy(transform[i], 0, totalTransform[i], 0, totalTransform[i].length);
2319 }
2320 }
2321
2322
2323
2324
2325
2326
2327
2328
2329 public static void addTranslation(final double[] translation, final double[][] totalTransform,
2330 final boolean prepend) {
2331 double[][] transform = new double[][]{{1.0, 0.0, 0.0, translation[0]},
2332 {0.0, 1.0, 0.0, translation[1]}, {0.0, 0.0, 1.0, translation[2]}, {0.0, 0.0, 0.0, 1.0}};
2333 transform =
2334 (prepend) ? mat4Mat4(totalTransform, transform) : mat4Mat4(transform, totalTransform);
2335 final int length = totalTransform.length;
2336 for (int i = 0; i < length; i++) {
2337 arraycopy(transform[i], 0, totalTransform[i], 0, totalTransform[i].length);
2338 }
2339 }
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
2351 private static void centerOfMass(final double[][] centersOfMass, final double[] coords, final double[] mass,
2352 final double massSum, final int nAtoms) {
2353 final int nAU = coords.length / (nAtoms * 3);
2354 for (int i = 0; i < nAU; i++) {
2355 final int auIndex = i * nAtoms * 3;
2356 final double[] comI = new double[3];
2357 for (int j = 0; j < nAtoms; j++) {
2358 final int atomIndex = auIndex + j * 3;
2359 final double m = mass[j];
2360 for (int k = 0; k < 3; k++) {
2361 comI[k] += coords[atomIndex + k] * m;
2362 }
2363 }
2364 for (int j = 0; j < 3; j++) {
2365 comI[j] /= massSum;
2366 }
2367 centersOfMass[i] = comI;
2368 }
2369 }
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382
2383
2384 private static boolean checkInflatedSphere(final double[] xyz, final double[] massN, final double massSum,
2385 final int compareAtomsSize, final int nAU, final int linkage, final int[] startLMN, final ArrayList<SymOp> symOps,
2386 final ArrayList<Integer> indexOrder) {
2387
2388 final double[][] centerOfMass = new double[xyz.length / 3][3];
2389 centerOfMass(centerOfMass, xyz, massN, massSum, compareAtomsSize);
2390 DoubleIndexPair[] auDist = new DoubleIndexPair[xyz.length / 3 / compareAtomsSize];
2391 prioritizeReplicates(xyz, centerOfMass, compareAtomsSize, auDist, 0, linkage);
2392
2393 boolean redo = false;
2394 for (int i = 0; i < nAU; i++) {
2395 if (logger.isLoggable(Level.FINE)) {
2396 logger.fine(format(
2397 " i: %3d\t auDist Index: %3d \t indexOrder: %3d\t indget(au): %3d\t ReplicatesVector: %2d (%2d) %2d (%2d) %2d (%2d)",
2398 i, auDist[i].index(), indexOrder.get(i), indexOrder.get(auDist[i].index()),
2399 symOps.get(indexOrder.get(auDist[i].index())).replicatesVector[0] + 1, startLMN[0],
2400 symOps.get(indexOrder.get(auDist[i].index())).replicatesVector[1] + 1, startLMN[1],
2401 symOps.get(indexOrder.get(auDist[i].index())).replicatesVector[2] + 1, startLMN[2]));
2402 }
2403 final int[] lmn = symOps.get(indexOrder.get(auDist[i].index())).replicatesVector;
2404 if (lmn[0] == 0 || lmn[0] == startLMN[0] - 1) {
2405 redo = true;
2406 }
2407 if (lmn[1] == 0 || lmn[1] == startLMN[1] - 1) {
2408 redo = true;
2409 }
2410 if (lmn[2] == 0 || lmn[2] == startLMN[2] - 1) {
2411 redo = true;
2412 }
2413 if (redo) {
2414 break;
2415 }
2416 }
2417 if (redo && logger.isLoggable(Level.FINE)) {
2418 logger.fine(" REDO EXPANSION.");
2419 }
2420 return redo;
2421 }
2422
2423
2424
2425
2426
2427
2428
2429 private static void copyArrayValues(double[][] newArray, double[][] oldArray) {
2430 final int arrayLength = newArray.length;
2431 for (int i = 0; i < arrayLength; i++) {
2432 arraycopy(oldArray[i], 0, newArray[i], 0, newArray[i].length);
2433 }
2434 }
2435
2436
2437
2438
2439
2440
2441
2442
2443
2444
2445
2446 private static void determineComparableAtoms(Atom[] atoms, ArrayList<Integer> indices,
2447 ArrayList<Integer> unique, boolean alphaCarbons, boolean includeHydrogen) {
2448 final int nAtoms = atoms.length;
2449 for (int i = 0; i < nAtoms; i++) {
2450 if (!unique.contains(i)) {
2451 final Atom atom = atoms[i];
2452 if (alphaCarbons) {
2453 final String atomName = atom.getName();
2454 final int atomAtNum = atom.getAtomicNumber();
2455 final boolean proteinCheck = atomName.equals("CA") && atomAtNum == 6;
2456 final boolean aminoCheck = (atomName.equals("N1") || atomName.equals("N9")) && atomAtNum == 7;
2457 if (proteinCheck || aminoCheck) {
2458 indices.add(i);
2459 }
2460 } else if (includeHydrogen || !atom.isHydrogen()) {
2461 indices.add(i);
2462 }
2463
2464 atom.setActive(true);
2465 }
2466 }
2467 }
2468
2469
2470
2471
2472
2473
2474
2475 private void gatherRMSDs(int row, RunningStatistics runningStatistics) {
2476 if (useMPI) {
2477 try {
2478 if (logger.isLoggable(Level.FINER)) {
2479 logger.finer(" Receiving MPI results.");
2480 }
2481
2482 world.gather(0, myBuffer, buffers);
2483 if (rank == 0) {
2484 for (int workItem = 0; workItem < numWorkItems; workItem++) {
2485 for (int proc = 0; proc < numProc; proc++) {
2486 final int column = numProc * workItem + proc;
2487
2488 if (column < targetSize) {
2489 distRow[column] = distances[proc][workItem];
2490 if (!isSymmetric) {
2491 runningStatistics.addValue(distRow[column]);
2492 } else if (column >= row) {
2493
2494 runningStatistics.addValue(distRow[column]);
2495 }
2496 if (logger.isLoggable(Level.FINER)) {
2497 logger.finer(format(" %d %d %14.8f", row, column, distances[proc][workItem]));
2498 }
2499 }
2500 }
2501 }
2502 } else {
2503 for (int workItem = 0; workItem < numWorkItems; workItem++) {
2504 final int column = numProc * workItem + rank;
2505
2506 if (column < targetSize) {
2507 distRow[column] = distances[rank][workItem];
2508 if (!isSymmetric) {
2509 runningStatistics.addValue(distRow[column]);
2510 } else if (column >= row) {
2511
2512 runningStatistics.addValue(distRow[column]);
2513 }
2514 if (logger.isLoggable(Level.FINER)) {
2515 logger.finer(format(" %d %d %14.8f", row, column, distances[rank][workItem]));
2516 }
2517 }
2518 }
2519 }
2520 } catch (Exception e) {
2521 logger.severe(" Exception collecting distance values." + e);
2522 }
2523 } else {
2524 for (int i = 0; i < targetSize; i++) {
2525 distRow[i] = myDistances[i];
2526 if (!isSymmetric) {
2527 runningStatistics.addValue(distRow[i]);
2528 } else if (i >= row) {
2529
2530 runningStatistics.addValue(distRow[i]);
2531 }
2532 }
2533 }
2534 }
2535
2536
2537
2538
2539
2540
2541
2542
2543
2544
2545
2546
2547
2548
2549
2550
2551
2552
2553
2554
2555 private static double[] generateInflatedSphere(final Crystal unitCell,
2556 final double[] reducedCoords, final int zPrime, final int nAU, final int linkage,
2557 final double[] mass, int[] lmn, final boolean strict, ArrayList<SymOp> symOps,
2558 ArrayList<Integer> indexOrder, double inflationFactor) {
2559 symOps.clear();
2560 indexOrder.clear();
2561
2562 final int nCoords = reducedCoords.length;
2563
2564 final int nAtoms = nCoords / 3;
2565
2566 final int zAtoms = nAtoms / zPrime;
2567
2568 final double[] x = new double[nAtoms];
2569 final double[] y = new double[nAtoms];
2570 final double[] z = new double[nAtoms];
2571 for (int i = 0; i < nAtoms; i++) {
2572 final int atomIndex = i * 3;
2573 x[i] = reducedCoords[atomIndex];
2574 y[i] = reducedCoords[atomIndex + 1];
2575 z[i] = reducedCoords[atomIndex + 2];
2576 }
2577
2578 final double volume = unitCell.volume / unitCell.getNumSymOps() / zPrime;
2579 final double radius = cbrt(0.75 / PI * volume * max(3, nAU) * inflationFactor);
2580 Crystal replicatesCrystal = ReplicatesCrystal.replicatesCrystalFactory(unitCell, radius * 2.0,
2581 lmn);
2582
2583
2584
2585
2586 final int nSymm = replicatesCrystal.getNumSymOps();
2587 final double[][] xS = new double[nSymm][nAtoms];
2588 final double[][] yS = new double[nSymm][nAtoms];
2589 final double[][] zS = new double[nSymm][nAtoms];
2590
2591 final int numEntities = nSymm * zPrime;
2592
2593 final double[][] cartCenterOfMass = new double[numEntities][3];
2594
2595
2596 List<SymOp> inflatedSymOps = replicatesCrystal.spaceGroup.symOps;
2597 for (int iSym = 0; iSym < nSymm; iSym++) {
2598 final SymOp symOp = inflatedSymOps.get(iSym);
2599
2600 final double[][] rot = new double[3][3];
2601 replicatesCrystal.getTransformationOperator(symOp, rot);
2602
2603 replicatesCrystal.applySymOp(nAtoms, x, y, z, xS[iSym], yS[iSym], zS[iSym], symOp);
2604 for (int zp = 0; zp < zPrime; zp++) {
2605 final int symIndex = zp * zAtoms;
2606
2607 final double[] centerOfMass = new double[3];
2608 double zMass = 0.0;
2609 for (int i = 0; i < zAtoms; i++) {
2610 final double m = mass[i];
2611 final int coordIndex = symIndex + i;
2612 centerOfMass[0] += xS[iSym][coordIndex] * m;
2613 centerOfMass[1] += yS[iSym][coordIndex] * m;
2614 centerOfMass[2] += zS[iSym][coordIndex] * m;
2615 zMass += m;
2616 }
2617 centerOfMass[0] /= zMass;
2618 centerOfMass[1] /= zMass;
2619 centerOfMass[2] /= zMass;
2620
2621 final double[] translate = moveIntoCrystal(replicatesCrystal, centerOfMass);
2622
2623 final double[] trans = symOp.tr.clone();
2624 replicatesCrystal.toCartesianCoordinates(trans, trans);
2625 final int translateLength = translate.length;
2626 for (int i = 0; i < translateLength; i++) {
2627 trans[i] += translate[i];
2628 }
2629 symOps.add(new SymOp(rot, trans, symOp.replicatesVector));
2630 for (int i = 0; i < zAtoms; i++) {
2631 final int coordIndex = symIndex + i;
2632 xS[iSym][coordIndex] += translate[0];
2633 centerOfMass[0] += translate[0];
2634 yS[iSym][coordIndex] += translate[1];
2635 centerOfMass[1] += translate[1];
2636 zS[iSym][coordIndex] += translate[2];
2637 centerOfMass[2] += translate[2];
2638 }
2639
2640
2641 cartCenterOfMass[iSym * zPrime + zp] = centerOfMass;
2642 }
2643 }
2644
2645
2646
2647
2648 final DoubleIndexPair[] auDist = new DoubleIndexPair[numEntities];
2649 final double[] cartCenter = new double[3];
2650
2651
2652
2653 final double[] fracCenter = {0.5, 0.5, 0.5};
2654 replicatesCrystal.toCartesianCoordinates(fracCenter, cartCenter);
2655
2656 if (logger.isLoggable(Level.FINER)) {
2657 if (logger.isLoggable(Level.FINEST)) {
2658 logger.finer(" Replicates crystal " + replicatesCrystal);
2659 }
2660 logger.finer(format(" Replicates Volume: %8.4f", replicatesCrystal.volume));
2661 logger.finer(
2662 format(" Expanded Crystal Center: %14.8f %14.8f %14.8f", cartCenter[0], cartCenter[1],
2663 cartCenter[2]));
2664 }
2665 logger.fine(
2666 format(" Expanded Crystal Center: %14.8f %14.8f %14.8f", cartCenter[0], cartCenter[1],
2667 cartCenter[2]));
2668
2669 for (int i = 0; i < numEntities; i++) {
2670
2671 auDist[i] = new DoubleIndexPair(i, dist2(cartCenter, cartCenterOfMass[i]));
2672 }
2673
2674
2675
2676 sort(auDist);
2677 for (DoubleIndexPair molsDist : auDist) {
2678 indexOrder.add(molsDist.index());
2679 }
2680
2681 if (logger.isLoggable(Level.FINEST)) {
2682 logger.finest("\n Copy SymOp Distance");
2683 }
2684 double[] systemCoords = new double[numEntities * zAtoms * 3];
2685 for (int n = 0; n < nSymm; n++) {
2686 for (int zp = 0; zp < zPrime; zp++) {
2687 final int index = n * zPrime + zp;
2688
2689 final int iSym = auDist[index].index();
2690 final double distance = auDist[index].doubleValue();
2691 if (logger.isLoggable(Level.FINEST) && n < 30) {
2692 logger.finest(format(" %4d %5d %14.8f", index, iSym, sqrt(distance)));
2693 }
2694
2695
2696 for (int i = 0; i < zAtoms; i++) {
2697 final int symIndex = index * zAtoms * 3;
2698 final int atomIndex = symIndex + i * 3;
2699 final int conversion = iSym / zPrime;
2700 final int shift = i + (iSym % zPrime) * zAtoms;
2701 systemCoords[atomIndex] = xS[conversion][shift];
2702 systemCoords[atomIndex + 1] = yS[conversion][shift];
2703 systemCoords[atomIndex + 2] = zS[conversion][shift];
2704 }
2705 }
2706 }
2707 final double massSum = Arrays.stream(mass).sum();
2708 if (logger.isLoggable(Level.FINE)) {
2709 logger.fine(" Checking replicates crystal.");
2710 logger.fine(format(
2711 " SysCoords: %3d (%3d)(%3d), reducedCoords Size: %3d (%3d)\n mass Size %3d, massSum: %6.3f, nAU: %2d, nCoords: %3d\n "
2712 + "auDist Size: %3d, l: %2d, m: %2d, n: %2d, numEntities: %3d, nAtoms: %3d, zPrime: %3d, zAtoms: %3d",
2713 systemCoords.length, systemCoords.length / 3, systemCoords.length / 3 / zAtoms,
2714 reducedCoords.length, reducedCoords.length / 3, mass.length, massSum, nAU, nCoords,
2715 auDist.length, lmn[0], lmn[1], lmn[2], numEntities, nAtoms, zPrime, zAtoms));
2716 }
2717 if (strict && checkInflatedSphere(systemCoords, mass, massSum, zAtoms, nAU, linkage, lmn, symOps,
2718 indexOrder)) {
2719
2720
2721 inflationFactor += 5;
2722 systemCoords = generateInflatedSphere(unitCell, reducedCoords, zPrime, nAU, linkage, mass, lmn,
2723 true, symOps, indexOrder, inflationFactor);
2724 }
2725 return systemCoords;
2726 }
2727
2728
2729
2730
2731
2732
2733
2734
2735 private static int guessZPrime(final ArrayList<Integer> unique, final int[] molNum, final int numEntities) {
2736 final int[] molSize = new int[numEntities];
2737 final int nAtoms = molNum.length;
2738 for(int i = 0; i < nAtoms; i++){
2739 if(unique.contains(i)){
2740 molSize[molNum[i]]++;
2741 }
2742 }
2743 int z;
2744 switch (numEntities) {
2745 case 2 -> {
2746 if (molSize[0] == molSize[1]) {
2747 z = 2;
2748 } else {
2749 z = 1;
2750 }
2751 }
2752 case 3 -> {
2753 if (molSize[0] == molSize[1] && molSize[1] == molSize[2]) {
2754 z = 3;
2755 } else {
2756 z = 1;
2757 }
2758 }
2759 case 4 -> {
2760 if (molSize[0] == molSize[1] && molSize[1] == molSize[2] && molSize[2] == molSize[3]) {
2761 z = 4;
2762 } else if (molSize[0] == molSize[2] && molSize[1] == molSize[3]) {
2763 z = 2;
2764 } else {
2765 z = 1;
2766 }
2767 }
2768 case 5 -> {
2769 if (molSize[0] == molSize[1] && molSize[1] == molSize[2] && molSize[2] == molSize[3] && molSize[3] == molSize[4]) {
2770 z = 5;
2771 } else {
2772 z = 1;
2773 }
2774 }
2775 case 6 -> {
2776 if (molSize[0] == molSize[1] && molSize[1] == molSize[2] && molSize[2] == molSize[3] && molSize[3] == molSize[4] && molSize[4] == molSize[5]) {
2777 z = 6;
2778 } else if (molSize[0] == molSize[2] && molSize[2] == molSize[4] && molSize[1] == molSize[3] && molSize[3] == molSize[5]) {
2779 z = 3;
2780 } else if (molSize[0] == molSize[3] && molSize[1] == molSize[4] && molSize[2] == molSize[5]) {
2781 z = 2;
2782 } else {
2783 z = 1;
2784 }
2785 }
2786 default -> z = 1;
2787 }
2788 if (logger.isLoggable(Level.FINER)) {
2789 logger.finer(format(" Number of species in asymmetric unit (Z Prime): %d", z));
2790 }
2791 return z;
2792 }
2793
2794
2795
2796
2797
2798
2799
2800
2801
2802
2803 private static boolean invalidAtomSelection(final ArrayList<Integer> indices, final Atom[] atoms,
2804 final boolean alphaCarbons, final boolean includeHydrogen) {
2805 final ArrayList<Integer> unique = new ArrayList<>();
2806 for (Integer value : indices) {
2807 if (!unique.contains(value)) {
2808 unique.add(value);
2809 }
2810 }
2811 indices.clear();
2812 determineComparableAtoms(atoms, indices, unique, alphaCarbons, includeHydrogen);
2813 return indices.isEmpty();
2814 }
2815
2816
2817
2818
2819
2820
2821
2822
2823
2824 public static String matrixToString(final double[][] matrix, final int index, String description) {
2825 StringBuilder value = new StringBuilder(format("\n %s %d: ", description, index));
2826 for (double[] doubles : matrix) {
2827 value.append(format("\n\t%s", Arrays.toString(doubles)));
2828 }
2829 value.append("\n");
2830 return value.toString();
2831 }
2832
2833
2834
2835
2836
2837
2838
2839
2840
2841 private static double[] moveIntoCrystal(final Crystal crystal, final double[] com) {
2842
2843 final double[] translate = new double[3];
2844 final double[] currentCoM = new double[3];
2845 currentCoM[0] = com[0];
2846 currentCoM[1] = com[1];
2847 currentCoM[2] = com[2];
2848
2849
2850 crystal.toFractionalCoordinates(com, translate);
2851 translate[0] = mod(translate[0], 1.0);
2852 translate[1] = mod(translate[1], 1.0);
2853 translate[2] = mod(translate[2], 1.0);
2854 crystal.toCartesianCoordinates(translate, translate);
2855
2856
2857 com[0] = translate[0];
2858 com[1] = translate[1];
2859 com[2] = translate[2];
2860
2861 translate[0] -= currentCoM[0];
2862 translate[1] -= currentCoM[1];
2863 translate[2] -= currentCoM[2];
2864
2865 return translate;
2866 }
2867
2868
2869
2870
2871
2872
2873
2874
2875
2876
2877
2878
2879
2880
2881
2882 private void numberUniqueAUs(double[] allCoords, DoubleIndexPair[] auDist, double[] auCoords,
2883 int nCoords, int upperLimit, boolean strict, int nAUinReplicates, double[] massN,
2884 double matchTol) {
2885
2886
2887 final List<Double> uniqueDiffs = new ArrayList<>();
2888 final ArrayList<double[]> tempListXYZ = new ArrayList<>();
2889 tempListIndices.add(0);
2890 uniqueDiffs.add(rmsd(auCoords, auCoords, massN));
2891 tempListXYZ.add(auCoords);
2892
2893 int index = 1;
2894
2895 final boolean useSym = printSym >= 0.0;
2896 final int numConfCheck = (strict || upperLimit <= 0 || useSym) ? nAUinReplicates : upperLimit;
2897 if (logger.isLoggable(Level.FINEST)) {
2898 logger.finest("RMSD Differences in Replicates Crystal:");
2899 }
2900 final double[] target = new double[nCoords];
2901 arraycopy(tempListXYZ.get(0), 0, target, 0, nCoords);
2902 translate(target, massN);
2903 while ((uniqueDiffs.size() < numConfCheck)) {
2904 final double[] baseCheckMol = new double[nCoords];
2905 final int auIndex = auDist[index].index();
2906 arraycopy(allCoords, auIndex * nCoords, baseCheckMol, 0, nCoords);
2907 translate(baseCheckMol, massN);
2908 rotate(target, baseCheckMol, massN);
2909 final double value = rmsd(target, baseCheckMol, massN);
2910 if (addLooseUnequal(uniqueDiffs, value, matchTol) || useSym) {
2911 if (logger.isLoggable(Level.FINEST)) {
2912 logger.finest(format(" Sorted Index: %4d Dist: %18.16f", auIndex, value));
2913 }
2914 tempListIndices.add(index);
2915 tempListXYZ.add(baseCheckMol);
2916 }
2917 index++;
2918 if (index >= auDist.length) {
2919 break;
2920 }
2921 }
2922 if (logger.isLoggable(Level.FINER)) {
2923 logger.finer(" RMSD_1 from 1st AU:\n i RMSD AU Index");
2924 final int numUnique = uniqueDiffs.size();
2925 for (int i = 0; i < numUnique; i++) {
2926 logger.finer(format(" %d %4.4f %4d", i, uniqueDiffs.get(i),
2927 auDist[tempListIndices.get(i)].index()));
2928 }
2929 }
2930 }
2931
2932
2933
2934
2935
2936
2937
2938 private void pairEntities(int desiredAUs, int comparisonNum) {
2939 tempListIndices.clear();
2940
2941 for (DoubleIndexPair doubleIndexPair : targetAUDist_2) {
2942
2943
2944 tempListIndices.add(doubleIndexPair.index());
2945 }
2946
2947
2948 for (int i = 0; i < desiredAUs; i++) {
2949 double minDist = Double.MAX_VALUE;
2950 Integer minIndex = -1;
2951 final double[] baseCoMCurr = baseCoM[comparisonNum][baseAUDist_2[i].index()];
2952 for (Integer target : tempListIndices) {
2953 final double dist = dist2(baseCoMCurr, targetCoM[target]);
2954 if (dist < minDist) {
2955 minIndex = target;
2956 minDist = dist;
2957 }
2958 if (abs(minDist) < MATCH_TOLERANCE) {
2959
2960 if (logger.isLoggable(Level.FINEST)) {
2961 stringBuilder.append("\n \tExit out of match loop.\n");
2962 }
2963 break;
2964 }
2965 }
2966 pairedAUs[i] = new DoubleIndexPair(minIndex, minDist);
2967 if (!tempListIndices.remove(minIndex)) {
2968 logger.warning(format(" Index value of %d was not found (%4.4f).", minIndex, minDist));
2969 }
2970 if (logger.isLoggable(Level.FINEST)) {
2971 stringBuilder.append(format("\n Base position: %d: %8.4f %8.4f %8.4f\n", i,
2972 baseCoM[comparisonNum][baseAUDist_2[i].index()][0],
2973 baseCoM[comparisonNum][baseAUDist_2[i].index()][1],
2974 baseCoM[comparisonNum][baseAUDist_2[i].index()][2]));
2975 stringBuilder.append(
2976 format(" Match Distance: %d: %8.4f\n", i, sqrt(pairedAUs[i].doubleValue())));
2977 int pairedIndex = pairedAUs[i].index();
2978 stringBuilder.append(
2979 format(" Target position: %d: %8.4f %8.4f %8.4f\n", i, targetCoM[pairedIndex][0],
2980 targetCoM[pairedIndex][1], targetCoM[pairedIndex][2]));
2981 }
2982 }
2983 if (logger.isLoggable(Level.FINER)) {
2984 stringBuilder.append(" Distance between pairs:\n"
2985 + " Index Base Index Target Index Match Index Distance\n");
2986 for (int i = 0; i < desiredAUs; i++) {
2987 stringBuilder.append(format(" %5d %10d %14d %14d %10.4f\n", i, baseAUDist_2[i].index(),
2988 targetAUDist_2[i].index(), pairedAUs[i].index(), sqrt(pairedAUs[i].doubleValue())));
2989 }
2990 }
2991 }
2992
2993
2994
2995
2996
2997
2998
2999
3000
3001
3002
3003
3004
3005 private void printSym(int compareAtomsSize, File file, String name, List<Bond> bondList,
3006 Atom[] atoms, ForceField forceField, int saveClusters, int currZ2) {
3007
3008
3009
3010 final double[][] tempTransform = new double[4][4];
3011 copyArrayValues(tempTransform, targetTransformSymOp.asMatrix());
3012 addTranslation(targetSymOp.tr, tempTransform, true);
3013 addRotation(targetSymOp.rot, tempTransform, true);
3014 final double[] bestTranslation = new double[]{tempTransform[0][3] / tempTransform[3][3],
3015 tempTransform[1][3] / tempTransform[3][3], tempTransform[2][3] / tempTransform[3][3]};
3016 final SymOp symOp = new SymOp(copyOf(tempTransform, 3), bestTranslation);
3017 final double[] newMol = new double[compareAtomsSize * 3];
3018 StringBuilder printString = new StringBuilder();
3019 String title = format("\n Print Current Symmetry Operator (Z'=%2d):\n"
3020 + " \t Original Coords \t\t Current Coords \t==\t Applied Sym Op Coords", currZ2);
3021 final double[] originalAU = targetAUoriginal[currZ2];
3022 for (int i = 0; i < compareAtomsSize; i++) {
3023 final int k = i * 3;
3024 final double[] xyz = new double[]{originalAU[k], originalAU[k + 1], originalAU[k + 2]};
3025 applyCartesianSymOp(xyz, xyz, symOp);
3026 newMol[k] = xyz[0];
3027 newMol[k + 1] = xyz[1];
3028 newMol[k + 2] = xyz[2];
3029 printString.append(
3030 format("\n %9.3f %9.3f %9.3f to %9.3f %9.3f %9.3f to %9.3f %9.3f %9.3f", originalAU[k],
3031 originalAU[k + 1], originalAU[k + 2], targetAU[k], targetAU[k + 1], targetAU[k + 2],
3032 newMol[k], newMol[k + 1], newMol[k + 2]));
3033 }
3034 title += format(" RMSD: %9.3f (should be 0.000)", rmsd(newMol, targetAU, massN));
3035 stringBuilder.append(title).append(printString);
3036 saveAssembly(file, name, bondList, atoms, forceField, newMol, comparisonAtoms, 1, "_moveMethod",
3037 0, saveClusters);
3038 stringBuilder.append("\n").append(symOp).append("\n");
3039 }
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049
3050
3051
3052 private static void prioritizeReplicates(final double[] coordsXYZ, final double[][] centerOfMasses, final int nAtoms,
3053 final DoubleIndexPair[] auDist, final int index, final int linkage) {
3054
3055
3056 final int nCoords = nAtoms * 3;
3057 final int length = coordsXYZ.length;
3058 final int nMols = length / nCoords;
3059 if (auDist.length != nMols) {
3060 logger.warning(" Number of molecules does not match distance sort array length.");
3061 }
3062 if (linkage == 0) {
3063
3064 final int centerIndex = index * nCoords;
3065 for (int i = 0; i < nMols; i++) {
3066 double tempDist = Double.MAX_VALUE;
3067 final int molIndex = i * nCoords;
3068 for (int j = 0; j < nAtoms; j++) {
3069 final int centerAtomIndex = centerIndex + j * 3;
3070 final double[] centerXYZ = {coordsXYZ[centerAtomIndex], coordsXYZ[centerAtomIndex + 1],
3071 coordsXYZ[centerAtomIndex + 2]};
3072 for (int k = 0; k < nAtoms; k++) {
3073 final int atomIndex = molIndex + k * 3;
3074 final double[] xyz = {coordsXYZ[atomIndex], coordsXYZ[atomIndex + 1],
3075 coordsXYZ[atomIndex + 2]};
3076 final double currDist = dist2(centerXYZ, xyz);
3077 if (currDist < tempDist) {
3078 tempDist = currDist;
3079 }
3080 if (abs(tempDist) < MATCH_TOLERANCE) {
3081 break;
3082 }
3083 }
3084 if (abs(tempDist) < MATCH_TOLERANCE) {
3085 break;
3086 }
3087 }
3088 auDist[i] = new DoubleIndexPair(i, tempDist);
3089 }
3090
3091 sort(auDist);
3092
3093 } else if (linkage == 1) {
3094
3095 final double[] coordCenter = centerOfMasses[index];
3096 for (int i = 0; i < nMols; i++) {
3097 final double[] moleculeCenter = centerOfMasses[i];
3098 auDist[i] = new DoubleIndexPair(i, dist2(coordCenter, moleculeCenter));
3099 }
3100
3101 sort(auDist);
3102
3103
3104
3105
3106 final DoubleIndexPair[] molDists = new DoubleIndexPair[nMols];
3107 int auIndex0 = auDist[0].index();
3108 final int auIndex1 = auDist[1].index();
3109 molDists[0] = new DoubleIndexPair(auIndex0, auDist[0].doubleValue());
3110
3111 double[] avgCenter = new double[3];
3112 avgCenter[0] = (centerOfMasses[auIndex0][0] + centerOfMasses[auIndex1][0]) / 2;
3113 avgCenter[1] = (centerOfMasses[auIndex0][1] + centerOfMasses[auIndex1][1]) / 2;
3114 avgCenter[2] = (centerOfMasses[auIndex0][2] + centerOfMasses[auIndex1][2]) / 2;
3115
3116 for (int i = 1; i < nMols; i++) {
3117 auIndex0 = auDist[i].index();
3118 final double[] moleculeCenter = centerOfMasses[auIndex0];
3119 molDists[i] = new DoubleIndexPair(auIndex0, dist2(avgCenter, moleculeCenter));
3120 }
3121 sort(molDists);
3122
3123
3124
3125 final DoubleIndexPair[] molDists3 = new DoubleIndexPair[nMols];
3126 molDists3[0] = new DoubleIndexPair(molDists[0].index(), molDists[0].doubleValue());
3127 molDists3[1] = new DoubleIndexPair(molDists[1].index(), molDists[1].doubleValue());
3128 avgCenter = new double[3];
3129 final int bound = min(3, molDists.length);
3130 for (int i = 0; i < bound; i++) {
3131 final int molIndex = molDists[i].index();
3132 avgCenter[0] += centerOfMasses[molIndex][0];
3133 avgCenter[1] += centerOfMasses[molIndex][1];
3134 avgCenter[2] += centerOfMasses[molIndex][2];
3135 }
3136 for (int i = 0; i < 3; i++) {
3137 avgCenter[i] /= 3;
3138 }
3139
3140 for (int i = 2; i < nMols; i++) {
3141 auIndex0 = molDists[i].index();
3142 final double[] moleculeCenter = centerOfMasses[auIndex0];
3143 molDists3[i] = new DoubleIndexPair(auIndex0, dist2(avgCenter, moleculeCenter));
3144 }
3145
3146 sort(molDists3);
3147 arraycopy(molDists3, 0, auDist, 0, nMols);
3148 } else if (linkage == 2) {
3149
3150 final int centerIndex = index * nCoords;
3151 for (int i = 0; i < nMols; i++) {
3152 double tempDist = 0.0;
3153 final int molIndex = i * nCoords;
3154 for (int j = 0; j < nAtoms; j++) {
3155 final int centerAtomIndex = centerIndex + j * 3;
3156 final double[] centerXYZ = {coordsXYZ[centerAtomIndex], coordsXYZ[centerAtomIndex + 1],
3157 coordsXYZ[centerAtomIndex + 2]};
3158 for (int k = 0; k < nAtoms; k++) {
3159 final int atomIndex = molIndex + k * 3;
3160 final double[] xyz = {coordsXYZ[atomIndex], coordsXYZ[atomIndex + 1],
3161 coordsXYZ[atomIndex + 2]};
3162 final double currDist = dist2(centerXYZ, xyz);
3163 if (currDist > tempDist) {
3164 tempDist = currDist;
3165 }
3166 if (abs(tempDist) < MATCH_TOLERANCE) {
3167 break;
3168 }
3169 }
3170 if (abs(tempDist) < MATCH_TOLERANCE) {
3171 break;
3172 }
3173 }
3174 auDist[i] = new DoubleIndexPair(i, tempDist);
3175 }
3176
3177 sort(auDist);
3178 } else {
3179 logger.severe(format(" Linkage value of %d was unrecognized. Try 0 (, 1, or 2.", linkage));
3180 }
3181 }
3182
3183
3184
3185
3186
3187
3188
3189 private RunningStatistics readMatrix(final String filename) {
3190 restartRow = 0;
3191 restartColumn = 0;
3192
3193 final DistanceMatrixFilter distanceMatrixFilter = new DistanceMatrixFilter();
3194 final RunningStatistics runningStatistics = distanceMatrixFilter.readDistanceMatrix(filename, baseSize,
3195 targetSize);
3196
3197 if (runningStatistics != null && runningStatistics.getCount() > 0) {
3198 restartRow = distanceMatrixFilter.getRestartRow();
3199 restartColumn = distanceMatrixFilter.getRestartColumn();
3200
3201 if (isSymmetric) {
3202
3203 if (restartRow == baseSize && restartColumn == 1) {
3204 logger.info(format(" Complete symmetric distance matrix found (%d x %d).", restartRow,
3205 restartRow));
3206 } else {
3207 restartColumn = 0;
3208 logger.info(format(
3209 " Incomplete symmetric distance matrix found.\n Restarting at row %d, column %d.",
3210 restartRow + 1, restartColumn + 1));
3211 }
3212 } else if (restartRow == baseSize && restartColumn == targetSize) {
3213 logger.info(format(" Complete distance matrix found (%d x %d).", restartRow, restartColumn));
3214 } else {
3215 restartColumn = 0;
3216 logger.info(format(" Incomplete distance matrix found.\n Restarting at row %d, column %d.",
3217 restartRow + 1, restartColumn + 1));
3218 }
3219 }
3220
3221 return runningStatistics;
3222 }
3223
3224
3225
3226
3227
3228
3229
3230
3231 private static double[] reduceSystem(final Atom[] atoms, final int[] comparisonAtoms) {
3232
3233 final double[] reducedCoords = new double[comparisonAtoms.length * 3];
3234 int coordIndex = 0;
3235 for (Integer value : comparisonAtoms) {
3236 final Atom atom = atoms[value];
3237 reducedCoords[coordIndex++] = atom.getX();
3238 reducedCoords[coordIndex++] = atom.getY();
3239 reducedCoords[coordIndex++] = atom.getZ();
3240 }
3241 return reducedCoords;
3242 }
3243
3244
3245
3246
3247
3248
3249
3250
3251
3252
3253
3254
3255
3256
3257
3258 private void saveAssembly(final File file, String name, final List<Bond> bondList, final Atom[] atoms,
3259 final ForceField forceField0, final double[] coords, final int[] comparisonAtoms, final int nAU,
3260 final String description, final int compNum, final int saveClusters) {
3261
3262 final String fileName = FilenameUtils.removeExtension(file.getName());
3263 File saveLocation;
3264 if (saveClusters == 2) {
3265 saveLocation = new File(fileName + description + ".xyz");
3266 } else if (saveClusters == 1) {
3267 saveLocation = new File(fileName + description + ".pdb");
3268 } else {
3269 return;
3270 }
3271
3272 final MolecularAssembly currentAssembly = new MolecularAssembly(name);
3273 final ArrayList<Atom> newAtomList = new ArrayList<>();
3274 int atomIndex = 1;
3275 final int nCoords = coords.length / (nAU);
3276
3277 try {
3278 for (int n = 0; n < nAU; n++) {
3279
3280
3281
3282
3283
3284 final ArrayList<Atom> atomList = new ArrayList<>();
3285
3286 int atomValue = 0;
3287
3288
3289
3290 final int[] atomSelection = new int[nCoords / 3];
3291 arraycopy(comparisonAtoms, 0, atomSelection, 0, nCoords / 3);
3292 for (final Integer i : atomSelection) {
3293 final Atom a = atoms[i];
3294
3295 final double[] xyz = new double[3];
3296 final int coordIndex = n * nCoords + atomValue;
3297 xyz[0] = coords[coordIndex];
3298 xyz[1] = coords[coordIndex + 1];
3299 xyz[2] = coords[coordIndex + 2];
3300 final Atom atom = new Atom(atomIndex++, a.getName(), a.getAtomType(), xyz);
3301 atomList.add(atom);
3302 atomValue += 3;
3303 }
3304
3305 if (saveClusters == 2) {
3306
3307 for (final Bond bond : bondList) {
3308 final Atom a1 = bond.getAtom(0);
3309 final Atom a2 = bond.getAtom(1);
3310
3311 final int a1Ind = a1.getIndex() - 1;
3312 final int a2Ind = a2.getIndex() - 1;
3313 if (IntStream.of(comparisonAtoms).anyMatch(x -> x == a1Ind) && IntStream.of(
3314 comparisonAtoms).anyMatch(x -> x == a2Ind)) {
3315 final Atom newA1 = atomList.get(a1Ind);
3316 final Atom newA2 = atomList.get(a2Ind);
3317 final Bond b = new Bond(newA1, newA2);
3318 b.setBondType(bond.getBondType());
3319 }
3320 }
3321 }
3322 newAtomList.addAll(atomList);
3323 }
3324 } catch (Exception exception) {
3325 stringBuilder.append("\n Error saving moved coordinates to PDB.\n").append(exception).append("\n");
3326 }
3327
3328 if (logger.isLoggable(Level.FINEST)) {
3329 final int newSize = newAtomList.size();
3330 stringBuilder.append(format("\n Save Num AU: %3d", nAU));
3331 stringBuilder.append(format("\n Save Num Active: %3d", nCoords));
3332 stringBuilder.append(format("\n Save Num Coords: %3d", coords.length));
3333 stringBuilder.append(format("\n Save Assembly Length: %3d", newSize));
3334 if (newSize > 0) {
3335 final Atom zero = newAtomList.get(0);
3336 stringBuilder.append(
3337 format("\n Save Assembly First Atom: %3s %9.3f %9.3f %9.3f\n", zero.getName(),
3338 zero.getX(), zero.getY(), zero.getZ()));
3339 } else {
3340 stringBuilder.append("\n WARNING: New list containing new atoms was empty.\n");
3341 }
3342 }
3343
3344
3345 final ForceField forceField = new ForceField(forceField0.getProperties());
3346
3347
3348 forceField.clearProperty("a-axis");
3349 forceField.clearProperty("b-axis");
3350 forceField.clearProperty("c-axis");
3351 forceField.clearProperty("alpha");
3352 forceField.clearProperty("beta");
3353 forceField.clearProperty("gamma");
3354 forceField.clearProperty("spacegroup");
3355 currentAssembly.setForceField(forceField);
3356
3357
3358
3359 Utilities.biochemistry(currentAssembly, newAtomList);
3360
3361 currentAssembly.setFile(file);
3362 if (saveClusters == 2) {
3363 final File key = new File(fileName + ".key");
3364 final File properties = new File(fileName + ".properties");
3365 if (key.exists()) {
3366 final File keyComparison = new File(fileName + description + ".key");
3367 try {
3368 if (keyComparison.createNewFile()) {
3369 FileUtils.copyFile(key, keyComparison);
3370 } else {
3371 stringBuilder.append("\n Could not create properties file.\n");
3372 }
3373 } catch (Exception ex) {
3374
3375 if (logger.isLoggable(Level.FINER)) {
3376 stringBuilder.append(ex);
3377 }
3378 }
3379 } else if (properties.exists()) {
3380 final File propertiesComparison = new File(fileName + description + ".properties");
3381 try {
3382 if (propertiesComparison.createNewFile()) {
3383 FileUtils.copyFile(properties, propertiesComparison);
3384 } else {
3385 stringBuilder.append("\n Could not create properties file.\n");
3386 }
3387 } catch (Exception ex) {
3388
3389 stringBuilder.append(
3390 "\n Neither key nor properties file detected therefore only creating XYZ.\n");
3391 if (logger.isLoggable(Level.FINER)) {
3392 stringBuilder.append(ex);
3393 }
3394 }
3395 }
3396 final XYZFilter xyzFilter = new XYZFilter(saveLocation, currentAssembly, forceField,
3397 currentAssembly.getProperties());
3398 xyzFilter.writeFile(saveLocation, true);
3399 } else {
3400 final PDBFilter pdbfilter = new PDBFilter(saveLocation, currentAssembly, forceField,
3401 currentAssembly.getProperties());
3402 pdbfilter.setModelNumbering(compNum);
3403 pdbfilter.writeFile(saveLocation, true, false, false);
3404 }
3405 currentAssembly.destroy();
3406 }
3407
3408
3409
3410
3411
3412
3413
3414
3415
3416
3417
3418 private void saveAssembly(final File file, final String name, final List<Bond> bondList, final Atom[] atoms,
3419 final ForceField forceField, double[] coords, final int[] comparisonAtoms, final int nAU, final String description,
3420 final double finalRMSD, final int compNum, final int saveClusters) {
3421 saveAssembly(file, name, bondList, atoms, forceField, coords, comparisonAtoms, nAU, description,
3422 compNum, saveClusters);
3423 final String fileName = FilenameUtils.removeExtension(file.getName());
3424 try {
3425
3426 File csv = PDBFilter.version(new File(fileName + description + ".csv"));
3427 if (csv.createNewFile()) {
3428 BufferedWriter bw = new BufferedWriter(new FileWriter(csv, false));
3429 bw.append("rms\n");
3430 bw.append(format("%4.4f\n", finalRMSD));
3431 bw.close();
3432 } else {
3433 logger.warning(" Could not create CSV file \"" + csv.getName() + "\"");
3434 }
3435 } catch (Exception ex) {
3436 logger.info(ex.toString());
3437 }
3438 }
3439 }
3440
3441
3442
3443
3444
3445
3446
3447
3448
3449
3450
3451
3452
3453
3454
3455
3456
3457
3458
3459
3460
3461
3462
3463
3464
3465
3466
3467
3468
3469
3470
3471
3472
3473
3474
3475
3476
3477
3478
3479
3480
3481
3482
3483
3484
3485
3486
3487
3488
3489
3490
3491
3492
3493
3494
3495
3496
3497
3498
3499
3500
3501
3502
3503
3504
3505
3506
3507
3508
3509
3510
3511
3512
3513
3514
3515
3516
3517
3518
3519
3520
3521
3522
3523
3524
3525
3526
3527
3528
3529
3530
3531
3532
3533
3534
3535
3536
3537
3538
3539
3540
3541
3542
3543
3544
3545
3546
3547
3548
3549
3550
3551
3552
3553
3554
3555
3556
3557
3558
3559
3560
3561
3562
3563
3564
3565
3566
3567
3568
3569
3570
3571
3572
3573
3574
3575
3576
3577
3578
3579
3580
3581
3582
3583
3584
3585
3586
3587
3588
3589
3590
3591
3592
3593
3594
3595
3596
3597
3598
3599
3600
3601
3602
3603
3604
3605
3606
3607
3608
3609
3610
3611
3612
3613
3614
3615
3616
3617
3618
3619
3620
3621
3622
3623
3624
3625
3626
3627
3628
3629
3630
3631
3632
3633
3634
3635
3636
3637
3638
3639
3640
3641
3642
3643
3644
3645
3646
3647
3648
3649
3650
3651
3652
3653
3654
3655
3656
3657
3658
3659
3660
3661
3662
3663
3664
3665
3666
3667
3668
3669
3670
3671
3672
3673
3674
3675
3676
3677
3678
3679
3680
3681
3682
3683
3684
3685
3686
3687
3688
3689
3690
3691
3692
3693
3694
3695
3696
3697
3698
3699
3700
3701
3702
3703
3704
3705
3706
3707
3708
3709
3710
3711
3712
3713
3714
3715
3716
3717
3718
3719
3720
3721
3722
3723
3724
3725
3726
3727
3728
3729
3730
3731
3732
3733
3734
3735
3736
3737
3738
3739
3740
3741
3742
3743
3744
3745
3746
3747
3748
3749
3750
3751
3752
3753
3754
3755
3756
3757
3758
3759
3760
3761
3762
3763
3764
3765
3766
3767
3768
3769
3770
3771
3772
3773
3774
3775
3776
3777
3778
3779
3780
3781
3782
3783
3784
3785
3786
3787
3788
3789
3790
3791
3792
3793
3794
3795
3796
3797
3798
3799
3800
3801
3802
3803
3804
3805
3806
3807
3808
3809
3810
3811
3812
3813
3814
3815
3816
3817
3818
3819
3820
3821
3822
3823
3824
3825
3826
3827
3828
3829
3830
3831
3832
3833
3834
3835
3836
3837
3838
3839
3840
3841
3842
3843
3844
3845
3846
3847
3848
3849
3850
3851
3852
3853
3854
3855
3856
3857
3858
3859
3860
3861
3862
3863
3864