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.algorithms.commands;
39
40 import ffx.algorithms.cli.AlgorithmsCommand;
41 import ffx.numerics.Potential;
42 import ffx.numerics.math.RunningStatistics;
43 import ffx.potential.bonded.Atom;
44 import ffx.potential.bonded.Bond;
45 import ffx.potential.parsers.SystemFilter;
46 import ffx.potential.utils.ProgressiveAlignmentOfCrystals;
47 import ffx.utilities.FFXBinding;
48 import picocli.CommandLine.Command;
49 import picocli.CommandLine.Option;
50 import picocli.CommandLine.Parameters;
51
52 import java.util.ArrayList;
53 import java.util.Collections;
54 import java.util.List;
55 import java.util.logging.Level;
56
57 import static ffx.utilities.StringUtils.parseAtomRanges;
58 import static java.lang.String.format;
59 import static org.apache.commons.io.FilenameUtils.concat;
60 import static org.apache.commons.io.FilenameUtils.getBaseName;
61 import static org.apache.commons.io.FilenameUtils.getFullPath;
62
63 @Command(description = " Determine the RMSD for crystal polymorphs using the Progressive Alignment of Crystals (PAC) algorithm.",
64 name = "SuperposeCrystals")
65 public class SuperposeCrystals extends AlgorithmsCommand {
66
67 @Option(names = {"--ac", "--alchemicalAtoms"}, paramLabel = "", defaultValue = "",
68 description = "Atom indices to be excluded from both crystals (e.g. 1-24,32-65). Use if molecular identity and atom ordering are the same for both crystals (otherwise use \"--ac1\" and \"--ac2\".")
69 private String excludeAtoms = "";
70
71 @Option(names = {"--ac1", "--alchemicalAtoms1"}, paramLabel = "", defaultValue = "",
72 description = "Atom indices to be excluded in the first crystal (e.g. 1-24,32-65).")
73 private String excludeAtomsA = "";
74
75 @Option(names = {"--ac2", "--alchemicalAtoms2"}, paramLabel = "", defaultValue = "",
76 description = "Atom indices to be excluded in the second crystal (e.g. 1-24,32-65).")
77 private String excludeAtomsB = "";
78
79 @Option(names = {"--ca", "--carbonAlphas"}, paramLabel = "false", defaultValue = "false",
80 description = "Consider only alpha carbons for proteins.")
81 private static boolean alphaCarbons;
82
83 @Option(names = {"--cd", "--createDirectories"}, paramLabel = "false", defaultValue = "false",
84 description = "Create subdirectories for free energy simulations.")
85 private static boolean createFE;
86
87 @Option(names = {"--as", "--autoSymmetry"}, paramLabel = "false", defaultValue = "false",
88 description = "Automatically generate symmetry operators.")
89 private static boolean autoSym;
90
91 @Option(names = {"--st", "--symTolerance"}, paramLabel = "0.5", defaultValue = "0.5",
92 description = "Add atoms to torsion based symmetries that are beneath this cutoff.")
93 private static double symTolerance;
94
95 @Option(names = {"--ws", "--writeSym"}, paramLabel = "false", defaultValue = "false",
96 description = "Write the created symmetry operators in the corresponding Properties/Key file.")
97 private static boolean writeSym;
98
99 @Option(names = {"--gc", "--gyrationComponents"}, paramLabel = "false", defaultValue = "false",
100 description = "Display components for radius of gyration for final clusters.")
101 private static boolean gyrationComponents;
102
103 @Option(names = {"--ht", "--hitTolerance"}, paramLabel = "-1.0", defaultValue = "-1.0",
104 description = "Sum comparisons that attain a value lower than this tolerance.")
105 private double hitTol;
106
107 @Option(names = {"--if", "--inflationFactor"}, paramLabel = "5.0", defaultValue = "5.0",
108 description = "Inflation factor used to determine replicates expansion (IF * nAU in replicates).")
109 private double inflationFactor;
110
111 @Option(names = {"--ih", "--includeHydrogen"}, paramLabel = "false", defaultValue = "false",
112 description = "Include hydrogen atoms.")
113 private boolean includeHydrogen;
114
115 @Option(names = {"--in", "--inertia"}, paramLabel = "false", defaultValue = "false",
116 description = "Display moments of inertia for final clusters.")
117 private static boolean inertia;
118
119 @Option(names = {"-l", "--linkage"}, paramLabel = "1", defaultValue = "1",
120 description = "Single (0), Average (1), or Complete (2) coordinate linkage for molecule prioritization.")
121 private int linkage;
122
123 @Option(names = {"--lm", "--lowMemory"}, paramLabel = "false", defaultValue = "false",
124 description = "Reduce memory usage at the cost of efficiency.")
125 private static boolean lowMemory;
126
127 @Option(names = {"--mt", "--moleculeTolerance"}, paramLabel = "0.0015", defaultValue = "0.0015",
128 description = "Tolerance to determine if two AUs are different.")
129 private double matchTol;
130
131 @Option(names = {"--mw", "--massWeighted"}, paramLabel = "false", defaultValue = "false",
132 description = "Use mass-weighted atomic coordinates for alignment.")
133 private static boolean massWeighted;
134
135 @Option(names = {"--na", "--numAU"}, paramLabel = "20", defaultValue = "20",
136 description = "Number of asymmetric units included to calculate the RMSD.")
137 private int numAU;
138
139 @Option(names = {"--pc", "--prioritizeCrystals"}, paramLabel = "0", defaultValue = "0",
140 description = "Prioritize crystals based on high density (0), low density (1), or file order (2).")
141 private int crystalPriority;
142
143 @Option(names = {"--ps", "--printSymOp"}, paramLabel = "-1.0", defaultValue = "-1.0",
144 description = "Print optimal SymOp to align input crystals (print out atom deviations above value).")
145 private static double printSym;
146
147 @Option(names = {"-r", "--restart"}, paramLabel = "false", defaultValue = "false",
148 description = "Restart from a previously written RMSD matrix (if one exists).")
149 private static boolean restart;
150
151 @Option(names = {"-s", "--save"}, paramLabel = "-1.0", defaultValue = "-1.0",
152 description = "Save structures less than or equal to this cutoff.")
153 private static double save;
154
155 @Option(names = {"--sc", "--saveClusters"}, paramLabel = "0", defaultValue = "0",
156 description = "Save files for the superposed crystal clusters (1=PDB, 2=XYZ).")
157 private static int saveClusters;
158
159 @Option(names = {"--sm", "--saveMachineLearning"}, paramLabel = "false", defaultValue = "false",
160 description = "Final structures for each comparison will be written out with RMSD in a CSV.")
161 private static boolean machineLearning;
162
163 @Option(names = {"--th", "--thorough"}, paramLabel = "false", defaultValue = "false",
164 description = "More intensive, less efficient version of PAC.")
165 private static boolean thorough;
166
167 @Option(names = {"-w", "--write"}, paramLabel = "false", defaultValue = "false",
168 description = "Write out the PAC RMSD matrix.")
169 private static boolean write;
170
171 @Option(names = {"--zp", "--zPrime"}, paramLabel = "-1", defaultValue = "-1",
172 description = "Z'' for both crystals (assumes same value).")
173 private int zPrime;
174
175 @Option(names = {"--zp1", "--zPrime1"}, paramLabel = "-1", defaultValue = "-1",
176 description = "Z'' for crystal 1 (default will try to autodetect).")
177 private int zPrime1;
178
179 @Option(names = {"--zp2", "--zPrime2"}, paramLabel = "-1", defaultValue = "-1",
180 description = "Z'' for crystal 2 (default will try to autodetect).")
181 private int zPrime2;
182
183 @Parameters(arity = "1..2", paramLabel = "files",
184 description = "Atomic coordinate file(s) to compare in XYZ format.")
185 List<String> filenames = null;
186
187 private final List<Atom> previouslyIncludedBase = new ArrayList<>();
188 private final List<Integer> queueIndices = new ArrayList<>();
189 public RunningStatistics runningStatistics;
190 public final StringBuilder symOpsA = new StringBuilder();
191 public final StringBuilder symOpsB = new StringBuilder();
192
193 public SuperposeCrystals() {
194 super();
195 }
196
197 public SuperposeCrystals(FFXBinding binding) {
198 super(binding);
199 }
200
201 public SuperposeCrystals(String[] args) {
202 super(args);
203 }
204
205 @Override
206 public SuperposeCrystals run() {
207 System.setProperty("vdwterm", "false");
208
209 if (!init()) {
210 return this;
211 }
212
213 if (filenames == null) {
214 logger.info(helpString());
215 return this;
216 }
217
218 algorithmFunctions.openAll(filenames.get(0));
219 SystemFilter baseFilter = algorithmFunctions.getFilter();
220 SystemFilter targetFilter;
221
222 boolean isSymmetric = false;
223 int numFiles = filenames.size();
224 if (numFiles == 1) {
225 logger.info(
226 "\n PAC will be applied between all pairs of conformations within the supplied file.\n");
227 isSymmetric = true;
228 algorithmFunctions.openAll(filenames.get(0));
229 targetFilter = algorithmFunctions.getFilter();
230 } else {
231 logger.info(
232 "\n PAC will compare all conformations in the first file to all those in the second file.\n");
233 algorithmFunctions.openAll(filenames.get(1));
234 targetFilter = algorithmFunctions.getFilter();
235 String filenameA = getBaseName(filenames.get(0));
236 String filenameB = getBaseName(filenames.get(1));
237 symOpsA.append(format("\n Add the following symop to properties/key file of %s (Apply on %s to generate %s):\nsymop", filenameA, filenameB, filenameA));
238 symOpsB.append(format("\n\n Add the following symop to properties/key file of %s (Apply on %s to generate %s):\nsymop", filenameB, filenameA, filenameB));
239 }
240
241 String filename = filenames.get(0);
242 String pacFilename = concat(getFullPath(filename), getBaseName(filename) + ".dst");
243
244 if (zPrime > 0 && zPrime % 1 == 0) {
245 zPrime1 = zPrime;
246 zPrime2 = zPrime;
247 }
248
249 if (printSym < 0 && writeSym || autoSym) {
250 logger.info("\n Printing distance for atoms greater than large distance of 10.0 Å");
251 printSym = 10.0;
252 }
253
254 if (createFE) {
255 lowMemory = true;
256 }
257
258 if (autoSym) {
259 boolean loggerLevelChanged = false;
260 if (logger.isLoggable(Level.INFO) && !logger.isLoggable(Level.FINE)) {
261 loggerLevelChanged = true;
262 logger.setLevel(Level.WARNING);
263 }
264
265 final Atom[] atomsBase = baseFilter.getActiveMolecularSystem().getAtomArray();
266 final Atom[] atomsTarget = targetFilter.getActiveMolecularSystem().getAtomArray();
267 final int nAtomsBase = atomsBase.length;
268 final int nAtomsTarget = atomsTarget.length;
269
270 List<List<Atom>> smallAtomGroups = new ArrayList<>();
271 List<Double> smallGroupsRMSD = new ArrayList<>();
272
273 List<Atom> mapBaseAtoms = new ArrayList<>();
274 List<Integer> mapBaseIndices = new ArrayList<>();
275 List<Atom> mapTargetAtoms = new ArrayList<>();
276 List<Integer> mapTargetIndices = new ArrayList<>();
277
278 if (!excludeAtoms.isEmpty()) {
279 List<Integer> alchAtoms = parseAtomRanges("AutoSym", excludeAtoms, nAtomsBase);
280 for (int i = 0; i < nAtomsBase - 1; i++) {
281 if (!alchAtoms.contains(i + 1)) {
282 logger.warning(format(" Include Atoms for both:\n %s\n %s.", atomsBase[i].toString(), atomsTarget[i].toString()));
283 mapBaseIndices.add(i + 1);
284 mapBaseAtoms.add(atomsBase[i]);
285 mapTargetIndices.add(i + 1);
286 mapTargetAtoms.add(atomsTarget[i]);
287 } else {
288 previouslyIncludedBase.add(atomsBase[i]);
289 }
290 }
291 } else if (!excludeAtomsA.isEmpty() || !excludeAtomsB.isEmpty()) {
292 logger.warning(format(" Independent alchemical atoms has not been implemented."));
293 }
294
295 assert (mapBaseIndices.size() == mapBaseAtoms.size());
296 assert (mapTargetIndices.size() == mapTargetAtoms.size());
297 if (mapBaseAtoms.size() != mapTargetAtoms.size()) {
298 logger.warning(format(" Error in atom selection. Base size of %3d does not match Target size of %3d.", mapBaseAtoms.size(), mapTargetAtoms.size()));
299 if (logger.isLoggable(Level.FINE)) {
300 StringBuilder sb = new StringBuilder();
301 for (Integer ind : mapBaseIndices) {
302 sb.append(ind).append(",");
303 }
304 sb.append("\n");
305 logger.fine(" Base indices:" + sb.toString());
306 sb.setLength(0);
307 for (Integer ind : mapTargetIndices) {
308 sb.append(ind).append(",");
309 }
310 sb.append("\n");
311 logger.fine(" Target indices:" + sb.toString());
312 }
313 return this;
314 }
315
316 for (Atom a : atomsBase) {
317 if (previouslyIncludedBase.contains(a)) {
318 continue;
319 }
320 queueIndices.clear();
321 addAtomIndex(a);
322
323 ArrayList<Atom> queueAtoms = new ArrayList<>();
324 for (Integer index : queueIndices) {
325 queueAtoms.add(atomsBase[index - 1]);
326 }
327
328 if (queueAtoms.size() < 3) {
329 boolean groupFound = false;
330 for (Atom atom : queueAtoms) {
331 Bond[] bonds = atom.getBonds().toArray(new Bond[0]);
332 for (Bond bond : bonds) {
333 Atom atomGrouped = bond.get1_2(atom);
334 for (List<Atom> atomList : smallAtomGroups) {
335 if (atomList.contains(atomGrouped)) {
336 groupFound = true;
337 for (Atom atomNeedingGroup : queueAtoms) {
338 atomList.add(atomNeedingGroup);
339 }
340 }
341 if (groupFound) {
342 break;
343 }
344 }
345 if (groupFound) {
346 break;
347 }
348 }
349 if (groupFound) {
350 break;
351 }
352 }
353 } else {
354 smallAtomGroups.add(queueAtoms);
355 }
356 }
357
358 for (int i = 0; i < nAtomsBase; i++) {
359 Atom atom = atomsBase[i];
360 boolean found = false;
361 for (List<Atom> list : smallAtomGroups) {
362 if (list.contains(atom)) {
363 found = true;
364 break;
365 }
366 }
367 if (!found) {
368 boolean groupFound = false;
369 Bond[] bonds = atom.getBonds().toArray(new Bond[0]);
370 for (Bond bond : bonds) {
371 Atom atomGrouped = bond.get1_2(atom);
372 for (List<Atom> atomList : smallAtomGroups) {
373 if (atomList.contains(atomGrouped)) {
374 groupFound = true;
375 atomList.add(atom);
376 }
377 if (groupFound) {
378 break;
379 }
380 }
381 if (groupFound) {
382 break;
383 }
384 }
385 }
386 }
387
388 StringBuilder sb = new StringBuilder(" Small Groups:\n");
389 for (List<Atom> group : smallAtomGroups) {
390 ArrayList<Integer> includeIndices = new ArrayList<>();
391 sb.append(" New Group: ");
392 for (Atom atom : group) {
393 sb.append(atom.getIndex() + ", ");
394 includeIndices.add(atom.getIndex());
395 }
396 sb.append("\n");
397
398 excludeAtoms = "";
399
400 for (int i = 0; i < nAtomsBase; i++) {
401 int ind = i + 1;
402 if (!includeIndices.contains(ind) && excludeAtoms.length() == 0) {
403 excludeAtoms = String.valueOf(ind);
404 } else if (!includeIndices.contains(ind)) {
405 excludeAtoms += "," + ind;
406 }
407 }
408
409 if (excludeAtoms == null || excludeAtoms.isEmpty()) {
410 continue;
411 } else {
412 excludeAtomsA = excludeAtoms;
413 excludeAtomsB = excludeAtoms;
414 }
415 if (logger.isLoggable(Level.FINE)) {
416 logger.fine("A: " + excludeAtomsA);
417 logger.fine("B: " + excludeAtomsB);
418 }
419
420 ProgressiveAlignmentOfCrystals pac = new ProgressiveAlignmentOfCrystals(baseFilter, targetFilter,
421 isSymmetric);
422 runningStatistics =
423 pac.comparisons(numAU, inflationFactor, matchTol, hitTol, zPrime1, zPrime2, excludeAtomsA, excludeAtomsB,
424 alphaCarbons, includeHydrogen, massWeighted, crystalPriority, thorough, saveClusters, save,
425 restart, write, machineLearning, inertia, gyrationComponents, linkage, printSym,
426 lowMemory, createFE, false, pacFilename, new StringBuilder(""), new StringBuilder(""));
427 double min = runningStatistics.getMin();
428 smallGroupsRMSD.add(min);
429 }
430
431 logger.info(sb.toString());
432 int numGroups = smallAtomGroups.size();
433 if (logger.isLoggable(Level.FINE)) {
434 logger.fine(format(" Number of minimal groups: %3d", numGroups));
435 }
436
437 int[][] smallIndexGroups = new int[numGroups][];
438 for (int i = 0; i < numGroups; i++) {
439 List<Atom> currentAtoms = smallAtomGroups.get(i);
440 int currentSize = currentAtoms.size();
441 smallIndexGroups[i] = new int[currentSize];
442 for (int j = 0; j < currentSize; j++) {
443 smallIndexGroups[i][j] = currentAtoms.get(j).getIndex();
444 }
445 }
446
447 ArrayList<ArrayList<Integer>> finalGroups = new ArrayList<>();
448 ArrayList<Integer> acceptedAtoms = new ArrayList<>();
449 ArrayList<Integer> queuedIndices = new ArrayList<>();
450 ArrayList<Integer> usedGroups = new ArrayList<>();
451
452 for (int i = 0; i < numGroups; i++) {
453 List<Atom> currentAtomGroup = smallAtomGroups.get(i);
454 int[] currentIndexGroup = smallIndexGroups[i];
455 int numAtomsInGroup = currentAtomGroup.size();
456 acceptedAtoms.clear();
457 queuedIndices.clear();
458 if (logger.isLoggable(Level.FINE)) {
459 logger.fine(format(" Minimal group number: %3d Size: %3d RMSD: %9.3f", i, numAtomsInGroup, smallGroupsRMSD.get(i)));
460 }
461 if (!usedGroups.contains(i)) {
462 for (int index : currentIndexGroup) {
463 acceptedAtoms.add(index);
464 }
465 usedGroups.add(i);
466
467 if (smallGroupsRMSD.get(i) > symTolerance) {
468 if (logger.isLoggable(Level.FINE)) {
469 logger.fine(format(" Run (group) %3d added outright to final.", i));
470 }
471 } else {
472 for (int j = 0; j < acceptedAtoms.size(); j++) {
473 Atom a1 = atomsBase[acceptedAtoms.get(j) - 1];
474 if (logger.isLoggable(Level.FINE)) {
475 logger.fine(format(" Current group: %3d Atom in Group: %3d of %3d (was %3d==%3d) Atom Index: %3d", i, j, acceptedAtoms.size(), numAtomsInGroup, currentIndexGroup.length, a1.getIndex()));
476 }
477 Bond[] b1 = a1.getBonds().toArray(new Bond[0]);
478 for (Bond b : b1) {
479 Atom a2 = b.get1_2(a1);
480 int a2Index = a2.getIndex();
481 boolean currentIndexGroupContains = false;
482 for (int idx : currentIndexGroup) {
483 if (idx == a2Index) {
484 currentIndexGroupContains = true;
485 break;
486 }
487 }
488 if (!currentIndexGroupContains) {
489 for (int k = 0; k < numGroups; k++) {
490 int[] currentIndexGroup2 = smallIndexGroups[k];
491 boolean currentIndexGroup2Contains = false;
492 for (int idx : currentIndexGroup2) {
493 if (idx == a2Index) {
494 currentIndexGroup2Contains = true;
495 break;
496 }
497 }
498 if (!usedGroups.contains(k) && smallGroupsRMSD.get(k) < symTolerance && currentIndexGroup2Contains) {
499 for (int ind : currentIndexGroup2) {
500 queuedIndices.add(ind);
501 }
502
503 excludeAtoms = "";
504 for (int l = 0; l < nAtomsBase; l++) {
505 int ind = l + 1;
506 boolean exclude = !acceptedAtoms.contains(ind) && !queuedIndices.contains(ind);
507 if (exclude) {
508 if (excludeAtoms.length() == 0) {
509 excludeAtoms = String.valueOf(ind);
510 } else {
511 excludeAtoms += "," + ind;
512 }
513 }
514 }
515
516 excludeAtomsA = excludeAtoms;
517 excludeAtomsB = excludeAtoms;
518
519 if (logger.isLoggable(Level.FINE)) {
520 logger.fine(" Excluded atoms A: " + excludeAtomsA);
521 logger.fine(" Excluded atoms B: " + excludeAtomsB);
522 }
523
524 ProgressiveAlignmentOfCrystals pac = new ProgressiveAlignmentOfCrystals(baseFilter, targetFilter,
525 isSymmetric);
526 runningStatistics =
527 pac.comparisons(numAU, inflationFactor, matchTol, hitTol, zPrime1, zPrime2, excludeAtomsA, excludeAtomsB,
528 alphaCarbons, includeHydrogen, massWeighted, crystalPriority, thorough, saveClusters, save,
529 restart, write, machineLearning, inertia, gyrationComponents, linkage, printSym,
530 lowMemory, createFE, false, pacFilename, new StringBuilder(""), new StringBuilder(""));
531 double min = runningStatistics.getMin();
532 if (min > symTolerance) {
533 if (logger.isLoggable(Level.FINE)) {
534 logger.fine(format(" Run %3d Group %3d failed with RMSD %9.3f > tolerance (%9.3f).", i, k, min, symTolerance));
535 }
536 queuedIndices.clear();
537 } else {
538 for (Integer value : queuedIndices) {
539 if (!acceptedAtoms.contains(value)) {
540 acceptedAtoms.add(value);
541 }
542 }
543 if (logger.isLoggable(Level.FINE)) {
544 logger.fine(format(" Run %3d Group %3d added to accepted with RMSD: %9.3f.", i, k, min));
545 }
546 usedGroups.add(k);
547 }
548 }
549 }
550 }
551 }
552 }
553 }
554 ArrayList<Integer> clone = new ArrayList<>(acceptedAtoms.size());
555 clone.addAll(acceptedAtoms);
556 finalGroups.add(clone);
557 if (logger.isLoggable(Level.FINE)) {
558 logger.fine(format(" Run %3d accepted atoms (size: %3d) added to finalGroups (size: %3d)", i, acceptedAtoms.size(), finalGroups.size()));
559 }
560 }
561 }
562
563 if (logger.isLoggable(Level.FINE)) {
564 logger.fine(format(" Final Num Groups: %3d", finalGroups.size()));
565 }
566
567 if (loggerLevelChanged) {
568 logger.setLevel(Level.INFO);
569 }
570
571 for (ArrayList<Integer> group : finalGroups) {
572 int groupSize = group.size();
573 if (logger.isLoggable(Level.FINE)) {
574 logger.fine(format(" Final Size: %3d", groupSize));
575 }
576 if (groupSize <= 0) {
577 logger.warning(" Final group of size zero was encountered. Check selection logic.");
578 continue;
579 }
580 queuedIndices.clear();
581 for (Integer ind : group) {
582 queuedIndices.add(ind);
583 }
584 excludeAtoms = "";
585 for (int l = 0; l < nAtomsBase; l++) {
586 int ind = l + 1;
587 if (!queuedIndices.contains(ind) && excludeAtoms.length() == 0) {
588 excludeAtoms = String.valueOf(ind);
589 } else if (!queuedIndices.contains(ind)) {
590 excludeAtoms += "," + ind;
591 }
592 }
593 excludeAtomsA = excludeAtoms;
594 excludeAtomsB = excludeAtoms;
595 if (logger.isLoggable(Level.FINE)) {
596 logger.fine("A: " + excludeAtomsA);
597 logger.fine("B: " + excludeAtomsB);
598 }
599
600 ProgressiveAlignmentOfCrystals pac = new ProgressiveAlignmentOfCrystals(baseFilter, targetFilter,
601 isSymmetric);
602 runningStatistics =
603 pac.comparisons(numAU, inflationFactor, matchTol, hitTol, zPrime1, zPrime2, excludeAtomsA, excludeAtomsB,
604 alphaCarbons, includeHydrogen, massWeighted, crystalPriority, thorough, saveClusters, save,
605 restart, write, machineLearning, inertia, gyrationComponents, linkage, printSym,
606 lowMemory, createFE, writeSym, pacFilename, symOpsA, symOpsB);
607 }
608 } else {
609 if (excludeAtoms != null && !excludeAtoms.isEmpty()) {
610 excludeAtomsA = excludeAtoms;
611 excludeAtomsB = excludeAtoms;
612 }
613 ProgressiveAlignmentOfCrystals pac = new ProgressiveAlignmentOfCrystals(baseFilter, targetFilter,
614 isSymmetric);
615 runningStatistics =
616 pac.comparisons(numAU, inflationFactor, matchTol, hitTol, zPrime1, zPrime2, excludeAtomsA, excludeAtomsB,
617 alphaCarbons, includeHydrogen, massWeighted, crystalPriority, thorough, saveClusters, save,
618 restart, write, machineLearning, inertia, gyrationComponents, linkage, printSym,
619 lowMemory, createFE, writeSym, pacFilename, symOpsA, symOpsB);
620 }
621
622 if (printSym >= 0.0) {
623 int finalCharA = symOpsA.length() - 2;
624 if (symOpsA.substring(finalCharA).equals("\\\n")) {
625 symOpsA.setLength(finalCharA);
626 }
627 int finalCharB = symOpsB.length() - 2;
628 if (symOpsB.substring(finalCharB).equals("\\\n")) {
629 symOpsB.setLength(finalCharB);
630 }
631 if (symOpsA.indexOf("\\") >= 0 && symOpsB.indexOf("\\") >= 0) {
632 logger.info(symOpsA.toString() + symOpsB.toString());
633 } else if (logger.isLoggable(Level.FINE)) {
634 logger.fine(" Following generated SymOp strings did not contain symmetry operators." + symOpsA.toString() + symOpsB.toString());
635 }
636 }
637
638 baseFilter.closeReader();
639 targetFilter.closeReader();
640 return this;
641 }
642
643 private void addAtomIndex(Atom aAdded) {
644 if (!previouslyIncludedBase.contains(aAdded)) {
645 queueIndices.add(aAdded.getIndex());
646 previouslyIncludedBase.add(aAdded);
647 int numBondsAdded = aAdded.getNumBonds();
648 int numHAdded = aAdded.getNumberOfBondedHydrogen();
649 int nonHbonds = numBondsAdded - numHAdded;
650
651 Bond[] bonds4Added = aAdded.getBonds().toArray(new Bond[0]);
652 for (Bond b : bonds4Added) {
653 Atom aBound2Added = b.get1_2(aAdded);
654 if (nonHbonds == 1) {
655 addAtomIndex(aBound2Added);
656 for (Bond bAdded : bonds4Added) {
657 Atom a3 = bAdded.get1_2(aAdded);
658 if (a3.isHydrogen()) {
659 addAtomIndex(a3);
660 }
661 }
662 }
663
664 int numBondsB2A = aBound2Added.getNumBonds();
665 int numHB2A = aBound2Added.getNumberOfBondedHydrogen();
666 int nonHB2A = numBondsB2A - numHB2A;
667 Bond[] bondsB2A = aBound2Added.getBonds().toArray(new Bond[0]);
668 if (nonHB2A == 1) {
669 addAtomIndex(aBound2Added);
670 for (Bond b2a : bondsB2A) {
671 Atom aBond2B2A = b2a.get1_2(aBound2Added);
672 if (aBond2B2A.isHydrogen()) {
673 addAtomIndex(aBond2B2A);
674 }
675 }
676 } else if (numBondsAdded == 1 || numBondsB2A == 1) {
677 addAtomIndex(aBound2Added);
678 } else if (nonHB2A == 2) {
679 addAtomIndex(aBound2Added);
680 for (Bond b2a : bondsB2A) {
681 Atom aBond2B2A = b2a.get1_2(aBound2Added);
682 if (aBond2B2A.isHydrogen()) {
683 addAtomIndex(aBond2B2A);
684 }
685 }
686 } else {
687 for (Bond b2a : bondsB2A) {
688 Atom aBond2B2A = b2a.get1_2(aBound2Added);
689 if (!aBond2B2A.isHydrogen() && nonHbonds == 1) {
690 int numBonds2B2A = aBond2B2A.getNumBonds();
691 if ((numBonds2B2A == 1) && (aBond2B2A.isBonded(aBound2Added) || aBond2B2A.isBonded(aAdded))) {
692 addAtomIndex(aBond2B2A);
693 }
694 }
695 }
696 }
697 }
698 }
699 }
700
701 @Override
702 public List<Potential> getPotentials() {
703 return Collections.emptyList();
704 }
705 }