View Javadoc
1   //******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  //******************************************************************************
38  package ffx.potential.commands.test;
39  
40  import ffx.potential.MolecularAssembly;
41  import ffx.potential.Utilities;
42  import ffx.potential.bonded.Angle;
43  import ffx.potential.bonded.Atom;
44  import ffx.potential.bonded.Bond;
45  import ffx.potential.bonded.MSGroup;
46  import ffx.potential.bonded.MSNode;
47  import ffx.potential.bonded.Molecule;
48  import ffx.potential.bonded.Polymer;
49  import ffx.potential.bonded.Residue;
50  import ffx.potential.cli.PotentialCommand;
51  import ffx.potential.parameters.ForceField;
52  import ffx.potential.parsers.CIFFilter;
53  import ffx.potential.parsers.PDBFilter;
54  import ffx.potential.parsers.SystemFilter;
55  import ffx.potential.parsers.XYZFilter;
56  import ffx.utilities.FFXBinding;
57  import org.apache.commons.configuration2.CompositeConfiguration;
58  import org.openscience.cdk.AtomContainer;
59  import org.openscience.cdk.config.AtomTypeFactory;
60  import org.openscience.cdk.graph.rebond.RebondTool;
61  import org.openscience.cdk.interfaces.IAtom;
62  import org.openscience.cdk.interfaces.IBond;
63  import org.openscience.cdk.isomorphism.AtomMatcher;
64  import org.openscience.cdk.isomorphism.BondMatcher;
65  import org.openscience.cdk.isomorphism.Pattern;
66  import org.openscience.cdk.isomorphism.VentoFoggia;
67  import picocli.CommandLine.Command;
68  import picocli.CommandLine.Option;
69  import picocli.CommandLine.Parameters;
70  
71  import javax.vecmath.Point3d;
72  import java.io.BufferedWriter;
73  import java.io.File;
74  import java.io.FileWriter;
75  import java.util.ArrayList;
76  import java.util.Arrays;
77  import java.util.Comparator;
78  import java.util.List;
79  import java.util.logging.Level;
80  
81  import static ffx.numerics.math.DoubleMath.dihedralAngle;
82  import static ffx.potential.bonded.BondedUtils.intxyz;
83  import static ffx.potential.utils.Superpose.applyRotation;
84  import static ffx.potential.utils.Superpose.applyTranslation;
85  import static ffx.potential.utils.Superpose.calculateRotation;
86  import static ffx.potential.utils.Superpose.calculateTranslation;
87  import static ffx.potential.utils.Superpose.rmsd;
88  import static ffx.potential.utils.Superpose.superpose;
89  import static java.lang.Double.MAX_VALUE;
90  import static java.lang.String.format;
91  import static org.apache.commons.io.FilenameUtils.getFullPath;
92  import static org.apache.commons.io.FilenameUtils.getName;
93  import static org.apache.commons.io.FilenameUtils.removeExtension;
94  import static org.apache.commons.math3.util.FastMath.sqrt;
95  import static org.openscience.cdk.tools.periodictable.PeriodicTable.getSymbol;
96  
97  /**
98   * The ReorderAtoms script sorts the atoms within an XYZ file based on atomic weight.
99   * Usage: ffxc test.ReorderAtoms <filename>
100  */
101 @Command(description = " Reorder the atoms of an XYZ file based on atomic weight.",
102     name = "test.ReorderAtoms")
103 public class ReorderAtoms extends PotentialCommand {
104 
105   /**
106    * -a or --atomType Change atom types of first file to match second.
107    */
108   @Option(names = {"-a", "--atomType"}, paramLabel = "false", defaultValue = "false",
109       description = "Change atom types of first file to match second.")
110   private static boolean atomType = false;
111 
112   /**
113    * --sw or --swap Switch positions of equivalent atoms (untested for >2 atoms).
114    */
115   @Option(names = {"--sw", "--swap"}, paramLabel = "false", defaultValue = "false",
116       description = "Switch positions of equivalent atoms (untested for >2 atoms).")
117   private static boolean swap = false;
118 
119   /**
120    * -m or --match Match atom indices based on structure.
121    */
122   @Option(names = {"-m", "--match"}, paramLabel = "false", defaultValue = "false",
123       description = "Match atom indices based on structure.")
124   private static boolean match = false;
125 
126   /**
127    * -r or --reorder Reorder atoms based on environment.
128    */
129   @Option(names = {"-r", "--reorder"}, paramLabel = "false", defaultValue = "false",
130       description = "Reorder atoms based on environment.")
131   private static boolean reorder = false;
132 
133   /**
134    * The final argument(s) should be one or more filenames.
135    */
136   @Parameters(arity = "1..*", paramLabel = "files",
137       description = "Atomic coordinate files to reorder in XYZ/ARC format. If 2, change file 1 to match file 2.")
138   private List<String> filenames;
139 
140   /**
141    * Set bond tolerance (distance to determine if bond is made).
142    */
143   private double bondTolerance = 0.2;
144 
145   /**
146    * Constructor.
147    */
148   public ReorderAtoms() {
149     super();
150   }
151 
152   /**
153    * Constructor using Binding.
154    */
155   public ReorderAtoms(FFXBinding binding) {
156     super(binding);
157   }
158 
159   /**
160    * Constructor using args.
161    */
162   public ReorderAtoms(String[] args) {
163     super(args);
164   }
165 
166   @Override
167   public ReorderAtoms run() {
168     logger.warning(" The ReorderAtoms command is under development.\n " +
169         " Please verify the output.");
170 
171     System.setProperty("vdwterm", "false");
172     if (!init()) {
173       return null; // Preserve Groovy behavior on failed init
174     }
175 
176     if (reorder || swap) {
177       if (match || atomType) {
178         logger.warning(" Currently each flag (i.e., atomType, match, or reorder must be run individually. ");
179       }
180       for (String filename : filenames) {
181         potentialFunctions.openAll(filename);
182         SystemFilter systemFilter = potentialFunctions.getFilter();
183         int numModels = systemFilter.countNumModels();
184 
185         File saveLocation = new File(filename);
186         File saveFile = potentialFunctions.versionFile(saveLocation);
187         for (int i = 0; i < numModels; i++) {
188           MolecularAssembly assembly = systemFilter.getActiveMolecularSystem();
189           Atom[] atoms = assembly.getAtomArray();
190           List<Atom> atomList0 = assembly.getAtomList();
191           int nAtoms = atoms.length;
192           MolecularAssembly newAssembly = new MolecularAssembly(assembly.getName());
193           List<Bond> bondList = assembly.getBondList();
194 
195           // Reorder/swap based on environment comparisons
196           int atomListSize = atomList0.size();
197           for (int j = 0; j < atomListSize; j++) {
198             for (int k = 0; k < atomListSize; k++) {
199               if (j != k) {
200                 if (atomList0.get(j).getIndex() != atomList0.get(k).getIndex()) {
201                   AtomComparator ac = new AtomComparator();
202                   if (logger.isLoggable(Level.FINER)) {
203                     logger.finer(" Compare atom: " + atomList0.get(j) + " and " + atomList0.get(k));
204                   }
205                   int value = ac.compare(atomList0.get(j), atomList0.get(k));
206                   if (value == 0) {
207                     logger.info(" This " + atomList0.get(j) + " and " + atomList0.get(k) +
208                         " are \"equivalent\".");
209                     if (swap) {
210                       Atom temp = atomList0.get(j);
211                       atomList0.set(j, atomList0.get(k));
212                       atomList0.set(k, temp);
213                     }
214                   } else if (value < 0 && reorder) {
215                     if (logger.isLoggable(Level.FINER)) {
216                       logger.finer(" Swapped " + atomList0.get(j) + " and " + atomList0.get(k) + ".");
217                     }
218                     Atom temp = atomList0.get(j);
219                     atomList0.set(j, atomList0.get(k));
220                     atomList0.set(k, temp);
221                   }
222                 }
223               }
224             }
225           }
226 
227           int[] originalOrder = new int[nAtoms];
228           for (int j = 0; j < nAtoms; j++) {
229             originalOrder[j] = atoms[j].getIndex();
230             if (logger.isLoggable(Level.FINE)) {
231               logger.fine(format(" Original index for atom %3d: %3d New index: %3d", j, originalOrder[j],
232                   atomList0.get(j).getIndex()));
233             }
234           }
235 
236           // Build a new list of atoms in the new order
237           ArrayList<Atom> atomList = new ArrayList<>();
238           int atomIndex = 1;
239           for (int j = 0; j < nAtoms; j++) {
240             Atom a = atomList0.get(j);
241             double[] xyz = new double[]{a.getX(), a.getY(), a.getZ()};
242             Atom atom = new Atom(atomIndex++, a.getName(), a.getAtomType(), xyz);
243             atomList.add(atom);
244           }
245 
246           // Recreate bonds in the new ordering
247           for (Bond bond : bondList) {
248             Atom a1 = bond.getAtom(0);
249             Atom a2 = bond.getAtom(1);
250 
251             int index1 = -1;
252             int index2 = -1;
253             int counter = 0;
254             for (Atom atom : atomList0) {
255               int index = atom.getIndex();
256               if (index == a1.getIndex()) {
257                 index1 = counter;
258               } else if (index == a2.getIndex()) {
259                 index2 = counter;
260               }
261               counter++;
262             }
263             if (index1 == -1 || index2 == -1) {
264               logger.severe(format(" Index1 (%3d) and Index2 (%3d) suggests error for bond: %s", index1, index2, bond));
265             }
266             Atom newA1 = atomList.get(index1);
267             Atom newA2 = atomList.get(index2);
268             Bond b = new Bond(newA1, newA2);
269             b.setBondType(bond.getBondType());
270           }
271 
272           // Construct the force field for the reordered atoms
273           newAssembly.setForceField(assembly.getForceField());
274           newAssembly.setPotential(assembly.getPotentialEnergy());
275           newAssembly.setCrystal(assembly.getCrystal());
276           newAssembly.setFile(assembly.getFile());
277           Utilities.biochemistry(newAssembly, atomList);
278 
279           XYZFilter filter = new XYZFilter(saveFile, newAssembly, assembly.getForceField(),
280               assembly.getProperties());
281           if (numModels > 1) {
282             filter.writeFile(saveFile, true);
283           } else {
284             filter.writeFile(saveFile, false);
285           }
286           systemFilter.readNext(false, false);
287 
288           if (numModels > 1) {
289             logger.info(format(" Saved structure %d to " + saveFile.getName(), i + 1));
290           } else {
291             logger.info(format(" Saved to " + saveFile.getName()));
292           }
293         }
294       }
295     } else if (match && filenames != null && filenames.size() > 1) {
296       potentialFunctions.openAll(filenames.get(0));
297       SystemFilter systemFilter1 = potentialFunctions.getFilter();
298       potentialFunctions.openAll(filenames.get(1));
299       SystemFilter systemFilter2 = potentialFunctions.getFilter();
300       logger.info(" Matching atom order between files.");
301       MolecularAssembly assembly1 = systemFilter1.getActiveMolecularSystem();
302       MolecularAssembly assembly2 = systemFilter2.getActiveMolecularSystem();
303 
304       List<MSNode> molecules1 = assembly1.getMolecules();
305       List<MSNode> molecules2 = assembly2.getMolecules();
306       int numMol2 = molecules2.size();
307       boolean[] molDone = new boolean[numMol2];
308       int totalNumAtoms = 0;
309 
310       ArrayList<Atom> atomList = new ArrayList<>();
311       int atomIndex = 1;
312       File saveLocation = new File(filenames.get(0));
313       File saveFile = potentialFunctions.versionFile(saveLocation);
314       MolecularAssembly newAssembly = new MolecularAssembly(assembly1.getName());
315       newAssembly.setForceField(assembly1.getForceField());
316       newAssembly.setPotential(assembly1.getPotentialEnergy());
317       newAssembly.setCrystal(assembly1.getCrystal());
318       newAssembly.setFile(assembly1.getFile());
319 
320       for (MSNode mol1 : molecules1) {
321         List<Atom> atoms1List = mol1.getAtomList();
322         int numAtoms = atoms1List.size();
323         logger.info(format(" Molecule1 with %3d atoms: %s", numAtoms, mol1.toString()));
324 
325         // Determine atoms that are unique (i.e., different environments).
326         List<List<Integer>> equivalents = new ArrayList<>();
327         boolean[] equivalent = new boolean[numAtoms];
328         int numEquiv = 0;
329         for (int i = 0; i < numAtoms; i++) {
330           for (int j = 1; j < numAtoms; j++) {
331             if (i == j) {
332               continue;
333             }
334             AtomComparator ac = new AtomComparator();
335             int value = ac.compare(atoms1List.get(i), atoms1List.get(j));
336             if (value == 0) {
337               equivalent[i] = true;
338               equivalent[j] = true;
339               int iIndex = -1;
340               int jIndex = -1;
341               for (int k = 0; k < numEquiv; k++) {
342                 List<Integer> list = equivalents.get(k);
343                 if (list.contains(i)) {
344                   iIndex = k;
345                 }
346                 if (list.contains(j)) {
347                   jIndex = k;
348                 }
349                 if (iIndex >= 0 && jIndex >= 0) {
350                   break;
351                 }
352               }
353               if (iIndex == -1 && jIndex == -1) {
354                 List<Integer> ij = new ArrayList<>();
355                 ij.add(i);
356                 ij.add(j);
357                 equivalents.add(ij);
358                 numEquiv++;
359               } else if (iIndex >= 0 && jIndex == -1) {
360                 equivalents.get(iIndex).add(j);
361               } else if (iIndex == -1 && jIndex >= 0) {
362                 equivalents.get(jIndex).add(i);
363               }
364             }
365           }
366         }
367         int numUnique = 0;
368         for (boolean value : equivalent) {
369           if (!value) {
370             numUnique++;
371           }
372         }
373 
374         // Arrays for all and unique coordinates/masses
375         double[] xyz1 = new double[numAtoms * 3];
376         double[] mass = new double[numAtoms];
377         double[] compCoords1 = new double[numUnique * 3];
378         double[] compMass = new double[numUnique];
379         int index = 0;
380         int indexUnique = 0;
381         for (int i = 0; i < numAtoms; i++) {
382           Atom a = atoms1List.get(i);
383           if (!equivalent[i]) {
384             compCoords1[indexUnique * 3] = a.getX();
385             compCoords1[indexUnique * 3 + 1] = a.getY();
386             compCoords1[indexUnique * 3 + 2] = a.getZ();
387             compMass[indexUnique++] = a.getMass();
388           }
389           xyz1[index * 3] = a.getX();
390           xyz1[index * 3 + 1] = a.getY();
391           xyz1[index * 3 + 2] = a.getZ();
392           mass[index++] = a.getMass();
393         }
394 
395         for (int l = 0; l < numMol2; l++) {
396           MSNode mol2 = molecules2.get(l);
397           List<Atom> atoms2List = mol2.getAtomList();
398           int numAtoms2 = atoms2List.size();
399           logger.info(format(" Molecule2 with %3d atoms: %s", numAtoms2, mol2.toString()));
400           if (numAtoms != numAtoms2 || molDone[l]) {
401             continue;
402           }
403           double[] xyz2 = new double[numAtoms2 * 3];
404           double[] compCoords2 = new double[numUnique * 3];
405           index = 0;
406           indexUnique = 0;
407           for (int i = 0; i < numAtoms; i++) {
408             Atom a = atoms2List.get(i);
409             if (!equivalent[i]) {
410               compCoords2[indexUnique * 3] = a.getX();
411               compCoords2[indexUnique * 3 + 1] = a.getY();
412               compCoords2[indexUnique++ * 3 + 2] = a.getZ();
413             }
414             xyz2[index * 3] = a.getX();
415             xyz2[index * 3 + 1] = a.getY();
416             xyz2[index++ * 3 + 2] = a.getZ();
417           }
418 
419           // If no atoms are unique, fall back to all atoms
420           if (compCoords1.length == 0 || compCoords2.length == 0) {
421             logger.info(format(" Comparison Coordinates invalid (1: %3d 2: %3d) %3d of %3d unique atoms.", compCoords1.length, compCoords2.length, numUnique, numAtoms));
422             compCoords1 = xyz1;
423             compCoords2 = xyz2;
424             compMass = mass;
425           }
426 
427           double[] translate1 = calculateTranslation(compCoords1, compMass);
428           double[] translate2 = calculateTranslation(compCoords2, compMass);
429           applyTranslation(compCoords1, translate1);
430           applyTranslation(compCoords2, translate2);
431           // Translate all coordinates based on unique atoms
432           applyTranslation(xyz1, translate1);
433           applyTranslation(xyz2, translate2);
434           double[][] rotation = calculateRotation(compCoords1, compCoords2, compMass);
435           // Rotate all coordinates based on unique atoms
436           applyRotation(compCoords2, rotation);
437           applyRotation(xyz2, rotation);
438           double compRMSD = superpose(compCoords1, compCoords2, compMass);
439           logger.info(format(" Unique atom RMSD: %9.4f A", compRMSD));
440           if (compRMSD > 1.0) {
441             logger.warning(format(" Value of %9.3f A may lead to insufficient overlap. Double check produced ordering.", compRMSD));
442           }
443 
444           // Loop through atoms that are not unique and try to find optimal match.
445           for (int i = 0; i < numAtoms; i++) {
446             if (equivalent[i]) {
447               int ind = i * 3;
448               double[] coord1 = new double[]{xyz1[ind], xyz1[ind + 1], xyz1[ind + 2]};
449               double[] coord2 = new double[]{xyz2[ind], xyz2[ind + 1], xyz2[ind + 2]};
450               double value = rmsd(coord1, coord2, new double[]{mass[i]});
451               logger.info(format(" Atom %2d distance of %9.3f A", i, value));
452               for (List<Integer> list : equivalents) {
453                 if (list.contains(i)) {
454                   int minIndex = -1;
455                   double minValue = MAX_VALUE;
456                   for (Integer aIndex : list) {
457                     if (i != aIndex) {
458                       int ind2 = aIndex * 3;
459                       double[] tempCoord = new double[]{xyz1[ind2], xyz1[ind2 + 1], xyz1[ind2 + 2]};
460                       double value2 = rmsd(tempCoord, coord2, new double[]{mass[i]});
461                       if (value2 < minValue) {
462                         minIndex = aIndex;
463                         minValue = value2;
464                       }
465                     }
466                   }
467                   if (minValue < value) {
468                     // Update Atom
469                     Atom temp = atoms1List.get(i);
470                     atoms1List.set(i, atoms1List.get(minIndex));
471                     atoms1List.set(minIndex, temp);
472                     // Update Coords
473                     double[] tempCoord = new double[]{xyz1[ind], xyz1[ind + 1], xyz1[ind + 2]};
474                     int ind2 = minIndex * 3;
475                     xyz1[ind] = xyz1[ind2];
476                     xyz1[ind + 1] = xyz1[ind2 + 1];
477                     xyz1[ind + 2] = xyz1[ind2 + 2];
478                     xyz1[ind2] = tempCoord[0];
479                     xyz1[ind2 + 1] = tempCoord[1];
480                     xyz1[ind2 + 2] = tempCoord[2];
481                   }
482                   value = rmsd(new double[]{xyz1[ind], xyz1[ind + 1], xyz1[ind + 2]},
483                       new double[]{xyz2[ind], xyz2[ind + 1], xyz2[ind + 2]}, new double[]{mass[i]});
484                   logger.info(format(" Atom %2d distance updated to %9.3f A", i, value));
485                   break;
486                 }
487               }
488             }
489           }
490 
491           compRMSD = superpose(xyz1, xyz2, mass);
492           logger.info(format(" Final RMSD: %9.4f A", compRMSD));
493 
494           // Create a new set of Atoms for this molecule in the new order
495           for (int j = 0; j < numAtoms; j++) {
496             Atom a = atoms1List.get(j);
497             double[] xyz = new double[]{a.getX(), a.getY(), a.getZ()};
498             Atom atom = new Atom(atomIndex++, a.getName(), a.getAtomType(), xyz);
499             atomList.add(atom);
500           }
501 
502           // Create bonds for this molecule mapping indices to new positions with offset
503           List<Bond> bondList = mol1.getBondList();
504           for (Bond bond : bondList) {
505             Atom a1 = bond.getAtom(0);
506             Atom a2 = bond.getAtom(1);
507 
508             int index1 = -1;
509             int index2 = -1;
510             int counter = 0;
511             for (Atom atom : atoms1List) {
512               int aIndex = atom.getIndex();
513               if (aIndex == a1.getIndex()) {
514                 index1 = counter;
515               } else if (aIndex == a2.getIndex()) {
516                 index2 = counter;
517               }
518               counter++;
519             }
520             if (index1 == -1 || index2 == -1) {
521               logger.severe(format(" Index1 (%3d) and Index2 (%3d) suggests error for bond: %s", index1, index2, bond));
522             }
523             Atom newA1 = atomList.get(index1 + totalNumAtoms);
524             Atom newA2 = atomList.get(index2 + totalNumAtoms);
525             logger.info(format(" index1 %3d Index2 %3d", index1, index2));
526             if (!newA1.isBonded(newA2)) {
527               Bond b = new Bond(newA1, newA2);
528               b.setBondType(bond.getBondType());
529             }
530           }
531 
532           molDone[l] = true;
533           totalNumAtoms += numAtoms;
534         }
535         if (!Arrays.asList(molDone).contains(false)) {
536           break;
537         }
538       }
539 
540       Utilities.biochemistry(newAssembly, atomList);
541       XYZFilter filter = new XYZFilter(saveFile, newAssembly, newAssembly.getForceField(),
542           newAssembly.getProperties());
543       filter.writeFile(saveFile, true);
544       logger.info(format(" Saved to " + saveFile.getName()));
545     } else if (atomType && filenames != null && filenames.size() > 1) {
546       potentialFunctions.openAll(filenames.get(0));
547       SystemFilter systemFilter1 = potentialFunctions.getFilter();
548       potentialFunctions.openAll(filenames.get(1));
549       SystemFilter systemFilter2 = potentialFunctions.getFilter();
550       logger.info(" Changing atom types in file 1 to match file 2.");
551       // Atom types from desired file.
552       MolecularAssembly assembly1 = systemFilter1.getActiveMolecularSystem();
553       do {
554         MolecularAssembly assembly2 = systemFilter2.getActiveMolecularSystem();
555         Atom[] atoms1 = assembly1.getAtomArray();
556         Atom[] atoms2 = assembly2.getAtomArray();
557         ArrayList<ArrayList<Atom>> xyzatoms = new ArrayList<>();
558         MolecularAssembly outputAssembly = new MolecularAssembly(assembly1.getName());
559 
560         int numHydrogen = 0;
561         for (Atom atom : atoms2) {
562           if (atom.isHydrogen()) {
563             numHydrogen++;
564           }
565         }
566 
567         List<MSNode> entitiesInput = assembly2.getAllBondedEntities();
568         int[] molIndices = assembly2.getMoleculeNumbers();
569         int numEntitiesInput = entitiesInput.size();
570         int[] shiftIndices = new int[numEntitiesInput];
571         int count = 0;
572         for (int i = 0; i < numEntitiesInput; i++) {
573           for (int ind : molIndices) {
574             if (ind == i) {
575               count++;
576             }
577           }
578           shiftIndices[i] = count;
579         }
580         if (logger.isLoggable(Level.FINE)) {
581           logger.fine(format(" Number of entities in input: %d", numEntitiesInput));
582           for (MSNode entity : entitiesInput) {
583             logger.fine(format(" Entity: %s", entity.getName()));
584             int size = entity.getAtomList().size();
585             logger.fine(format("   Entity Size: %3d", size));
586             if (size > 0) {
587               logger.fine(format("   Entity First Atom: %s", entity.getAtomList().get(0).toString()));
588             } else {
589               logger.warning(" Entity did not contain atoms.");
590             }
591           }
592         }
593 
594         for (MSNode mol : entitiesInput) {
595           //noinspection unchecked
596           xyzatoms.add((ArrayList<Atom>) mol.getAtomList());
597         }
598         int atomIndex = 1;
599         // For each entity (protein, molecule, ion, etc) in the XYZ
600         logger.info(" Num Entities Input: " + numEntitiesInput);
601         for (int i = 0; i < numEntitiesInput; i++) {
602           MSNode mol = entitiesInput.get(i);
603           int numInputMolAtoms = xyzatoms.get(i).size();
604           int numMolHydrogen = 0;
605           for (Atom atom : xyzatoms.get(i)) {
606             if (atom.isHydrogen()) {
607               numMolHydrogen++;
608             }
609           }
610           if (logger.isLoggable(Level.FINE)) {
611             logger.fine(format(" Current entity number of atoms: %d (%d + %dH)", numInputMolAtoms,
612                 numInputMolAtoms - numMolHydrogen, numMolHydrogen));
613           }
614 
615           // Set up input (XYZ:assembly 2) file contents as CDK variable
616           AtomContainer xyzCDKAtoms = new AtomContainer();
617           for (Atom atom : xyzatoms.get(i)) {
618             String atomName = getSymbol(atom.getAtomType().atomicNumber);
619             xyzCDKAtoms.addAtom(new org.openscience.cdk.Atom(atomName, new Point3d(atom.getXYZ(null))));
620           }
621 
622           int zPrime;
623           int nAtoms = atoms1.length;
624           int nInputAtoms = atoms2.length;
625           if (nAtoms % nInputAtoms == 0) {
626             zPrime = (nAtoms / nInputAtoms);
627           } else if (nAtoms % (nInputAtoms - numHydrogen) == 0) {
628             zPrime = (nAtoms / (nInputAtoms - numHydrogen));
629           } else {
630             zPrime = 1;
631           }
632           logger.info(" Original ZPrime: " + zPrime);
633           // Add known input bonds; a limitation is that bonds all are given a Bond order of 1.
634           List<Bond> bonds = mol.getBondList();
635           IBond.Order order = IBond.Order.SINGLE;
636           int xyzBonds = bonds.size();
637           if (xyzBonds == 0) {
638             logger.warning(" No bonds detected in input structure. Please check input.\n " +
639                 "If correct, separate non-bonded entities into multiple CIFs.");
640           }
641           int shiftIndex;
642           for (Bond xyzBond : bonds) {
643             int atom0Index = xyzBond.getAtom(0).getXyzIndex();
644             int atom1Index = xyzBond.getAtom(1).getXyzIndex();
645             shiftIndex = (molIndices[atom0Index] == 0) ? 0 : shiftIndices[molIndices[atom0Index] - 1];
646             if (logger.isLoggable(Level.FINER)) {
647               logger.finer(format(" Bonded atom 1: %d, Bonded atom 2: %d", atom0Index, atom1Index));
648               logger.finer(format(" Atom0: %3d Atom1: %3d Shift: %3d final0,1: %3d, %3d MolIndex: %3d ShiftIndex: %3d",
649                   atom0Index, atom1Index, shiftIndex, atom0Index - shiftIndex - 1, atom1Index - shiftIndex - 1,
650                   molIndices[atom0Index], shiftIndices[molIndices[atom0Index]]));
651             }
652             xyzCDKAtoms.addBond(atom0Index - shiftIndex - 1, atom1Index - shiftIndex - 1, order);
653           }
654 
655           // Assign CDK atom types for the input molecule.
656           AtomTypeFactory factory = AtomTypeFactory.getInstance(
657               "org/openscience/cdk/config/data/jmol_atomtypes.txt", xyzCDKAtoms.getBuilder());
658 
659           // Assign atom types to CDK object.
660           for (IAtom atom : xyzCDKAtoms.atoms()) {
661             CIFFilter.setAtomTypes(factory, atom);
662             try {
663               factory.configure(atom);
664             } catch (Exception ex) {
665               logger.info(" Failed to configure atoms from CIF.\n" + ex + "\n" + Arrays.toString(
666                   ex.getStackTrace()));
667             }
668           }
669 
670           ArrayList<ArrayList<Integer>> zindices = new ArrayList<>();
671           int counter = 0;
672           // Bond atoms from CIF.
673           int cifBonds = CIFFilter.bondAtoms(atoms1, bondTolerance);
674           if (logger.isLoggable(Level.FINE)) {
675             logger.fine(
676                 format(" Created %d bonds between input atoms (%d in input).", cifBonds, xyzBonds));
677           }
678           List<Atom> atomPool = new ArrayList<>(Arrays.asList(atoms1));
679 
680           try {
681             while (!atomPool.isEmpty()) {
682               ArrayList<Atom> molecule = new ArrayList<>();
683               CIFFilter.collectAtoms(atomPool.get(0), molecule);
684               if (logger.isLoggable(Level.FINER)) {
685                 logger.finer(format(" Molecule (%d) Size: %d", counter, molecule.size()));
686               }
687               ArrayList<Integer> indices = new ArrayList<>();
688               while (!molecule.isEmpty()) {
689                 Atom atom = molecule.remove(0);
690                 indices.add(atom.getIndex());
691                 atomPool.remove(atom);
692               }
693 
694               if (logger.isLoggable(Level.FINER)) {
695                 logger.finer(format(
696                     " Molecule %d: %d atoms are ready and %d remain (%d atoms in input, %d atoms in CIF). ",
697                     counter + 1, indices.size(), atomPool.size(), nInputAtoms, nAtoms));
698               }
699               zindices.add(indices);
700               counter++;
701             }
702           } catch (Exception e) {
703             logger.severe(" Failed to separate copies within the asymmetric unit." + e + "\n" + Utilities.stackTraceToString(e));
704           }
705           zPrime = zindices.size();
706           logger.info(" Zprime: " + zPrime);
707           // Set up CIF contents as CDK variable
708           AtomContainer[] cifCDKAtomsArr = new AtomContainer[zPrime];
709           for (int j = 0; j < zPrime; j++) {
710             if (zPrime > 1) {
711               logger.info(format("\n Attempting entity %d of %d", j + 1, zPrime));
712             }
713             ArrayList<Integer> currentList = zindices.get(j);
714             int cifMolAtoms = currentList.size();
715             if (logger.isLoggable(Level.FINE)) {
716               logger.fine(format(" CIF atoms in current: %d", cifMolAtoms));
717             }
718             // Detect if CIF contains multiple copies (Z'>1)
719             if (cifMolAtoms % numInputMolAtoms == 0 || cifMolAtoms % (numInputMolAtoms - numMolHydrogen) == 0) {
720               cifCDKAtomsArr[j] = new AtomContainer();
721               for (Integer integer : currentList) {
722                 String atomName = getSymbol(atoms1[integer - 1].getAtomType().atomicNumber);
723                 cifCDKAtomsArr[j].addAtom(new org.openscience.cdk.Atom(atomName, new Point3d(atoms1[integer - 1].getXYZ(null))));
724               }
725               AtomContainer nullCDKAtoms = new AtomContainer();
726 
727               for (IAtom atom : cifCDKAtomsArr[j].atoms()) {
728                 if (atom.toString() == null) {
729                   nullCDKAtoms.addAtom(atom);
730                 }
731               }
732               cifCDKAtomsArr[j].remove(nullCDKAtoms);
733 
734               // Compute bonds for CIF molecule.
735               factory = AtomTypeFactory.getInstance(
736                   "org/openscience/cdk/config/data/jmol_atomtypes.txt", cifCDKAtomsArr[j].getBuilder());
737               for (IAtom atom : cifCDKAtomsArr[j].atoms()) {
738                 CIFFilter.setAtomTypes(factory, atom);
739                 try {
740                   factory.configure(atom);
741                 } catch (Exception ex) {
742                   logger.info(" Failed to configure CIF atoms.\n" + ex + "\n" + Utilities.stackTraceToString(ex));
743                 }
744               }
745 
746               RebondTool rebonder = new RebondTool(CIFFilter.MAX_COVALENT_RADIUS, CIFFilter.MIN_BOND_DISTANCE,
747                   bondTolerance);
748               try {
749                 rebonder.rebond(cifCDKAtomsArr[j]);
750               } catch (Exception ex) {
751                 logger.info("Failed to rebond CIF atoms.\n" + ex + "\n" + Utilities.stackTraceToString(ex));
752               }
753 
754               int cifMolBonds = cifCDKAtomsArr[j].getBondCount();
755               if (logger.isLoggable(Level.FINE)) {
756                 logger.fine(format(" Number of CIF bonds: %d (%d in input)", cifMolBonds, xyzBonds));
757               }
758               // Number of bonds matches.
759               // If cifMolBonds == 0 then ion or atom with implicit hydrogens (e.g. water, methane, etc.)
760               if (cifMolBonds != 0 && cifMolBonds % xyzBonds == 0) {
761                 Pattern pattern = VentoFoggia.findIdentical(
762                     xyzCDKAtoms, AtomMatcher.forElement(), BondMatcher.forAny());
763                 int[] p = pattern.match(cifCDKAtomsArr[j]);
764                 int pLength = p.length;
765                 if (p != null && pLength == numInputMolAtoms) {
766                   // Use matched atoms to update the positions of the input file atoms.
767                   for (int k = 0; k < pLength; k++) {
768                     if (logger.isLoggable(Level.FINEST)) {
769                       logger.finest(
770                           format(" %d input %s -> CIF %s", k, xyzCDKAtoms.getAtom(k).getSymbol(),
771                               cifCDKAtomsArr[j].getAtom(p[k]).getSymbol()));
772                     }
773                     Point3d point3d = cifCDKAtomsArr[j].getAtom(p[k]).getPoint3d();
774                     xyzatoms.get(i).get(k).setXYZ(new double[]{point3d.x, point3d.y, point3d.z});
775                   }
776                 } else {
777                   if (logger.isLoggable(Level.FINE)) {
778                     logger.fine(
779                         format(" Atoms from CIF (%d) and input (%d) structures don't match.", p.length,
780                             nAtoms));
781                   }
782                   continue;
783                 }
784               } else if ((xyzBonds - numMolHydrogen) == 0 || cifMolBonds % ((xyzBonds - numMolHydrogen)) == 0) {
785                 // Hydrogen most likely missing from file. If zero, then potentially water/methane (implicit hydrogen atoms).
786                 if (logger.isLoggable(Level.FINE)) {
787                   logger.info(" CIF may contain implicit hydrogen -- attempting to patch.");
788                 }
789                 // Match heavy atoms between CIF and input
790                 Pattern pattern = VentoFoggia.findSubstructure(
791                     cifCDKAtomsArr[j], AtomMatcher.forElement(), BondMatcher.forAny());
792                 int[] p = pattern.match(xyzCDKAtoms);
793                 int pLength = p.length;
794                 if (p != null && pLength == numInputMolAtoms - numMolHydrogen) {
795                   // Used matched atoms to update the positions of the input file atoms.
796                   for (int k = 0; k < pLength; k++) {
797                     Point3d point3d = cifCDKAtomsArr[j].getAtom(k).getPoint3d();
798                     xyzatoms.get(i).get(p[k]).setXYZ(new double[]{point3d.x, point3d.y, point3d.z});
799                   }
800                   // Add in hydrogen atoms
801                   Atom lastKnownAtom1 = null;
802                   int knownCount = 0;
803                   for (Atom hydrogen : xyzatoms.get(i)) {
804                     if (hydrogen.isHydrogen()) {
805                       Bond bond0 = hydrogen.getBonds().get(0);
806                       Atom atom1 = bond0.get1_2(hydrogen);
807                       double angle0_2 = 0;
808                       List<Angle> anglesList = hydrogen.getAngles(); // Same as bond
809                       Atom atom2 = null;
810                       switch (anglesList.size()) {
811                         case 0 ->
812                           // H-Cl No angles
813                           // Place hydrogen slightly off center of bonded atom (~1Å away).
814                             hydrogen.moveTo(new double[]{atom1.getX() - 0.6, atom1.getY() - 0.6,
815                                 atom1.getZ() - 0.6});
816                         case 1 -> {
817                           // H-O=C Need torsion
818                           // H-O-H
819                           for (Angle angle : anglesList) {
820                             atom2 = angle.get1_3(hydrogen);
821                             if (atom2 != null) {
822                               angle0_2 = angle.angleType.angle[0];
823                             }
824                             if (angle0_2 != 0) {
825                               //if angle0_2 is found then no need to look for another atom.
826                               break;
827                             }
828                           }
829                           assert atom2 != null;
830                           List<Bond> bonds2 = atom2.getBonds();
831                           Atom atom3 =
832                               (bonds2.size() > 1 && atom1 == bonds2.get(0).get1_2(atom2)) ? bonds2.get(
833                                   1).get1_2(atom2) : bonds2.get(0).get1_2(atom2);
834                           double diAng = dihedralAngle(hydrogen.getXYZ(null), atom1.getXYZ(null),
835                               atom2.getXYZ(null), atom3.getXYZ(null));
836                           if (atom1 != atom3) {
837                             intxyz(hydrogen, atom1, bond0.bondType.distance, atom2, angle0_2, atom3,
838                                 Math.toDegrees(diAng), 0);
839                           } else {
840                             // Likely water as Atom 3 is not unique. Since no hydrogen atoms
841                             //  are present, there isn't a third atom...
842                             double[] coord = new double[]{atom2.getX(), atom2.getY(), atom3.getZ()};
843                             double mag = sqrt(
844                                 coord[0] * coord[0] + coord[1] * coord[1] + coord[2] * coord[2]);
845                             coord[0] /= mag;
846                             coord[1] /= mag;
847                             coord[2] /= mag;
848                             hydrogen.moveTo(atom1.getX() - coord[0], atom1.getY() - coord[1],
849                                 atom1.getZ() - coord[2]);
850                           }
851                         }
852                         default -> {
853                           // H-C(-C)(-C)-C
854                           Atom atom2B = null;
855                           double angle0_2B = 0;
856                           Atom proposedAtom;
857                           double proposedAngle = 0;
858                           int chiral = 1;
859                           for (Angle angle : anglesList) {
860                             proposedAtom = angle.get1_3(hydrogen);
861                             if (proposedAtom != null && !proposedAtom.isHydrogen()) {
862                               proposedAngle = angle.angleType.angle[0];
863                             }
864                             if (proposedAngle != 0) {
865                               // If angle1_3 is found then no need to look for another atom.
866                               if (angle0_2 != 0) {
867                                 atom2B = proposedAtom;
868                                 angle0_2B = proposedAngle;
869                               } else {
870                                 atom2 = proposedAtom;
871                                 angle0_2 = proposedAngle;
872                               }
873                               proposedAngle = 0.0;
874                             }
875                             if (angle0_2 != 0 && angle0_2B != 0) {
876                               break;
877                             }
878                           }
879                           if (lastKnownAtom1 == null || lastKnownAtom1 != atom1) {
880                             lastKnownAtom1 = atom1;
881                             knownCount = 0;
882                           } else {
883                             knownCount++;
884                           }
885                           if (angle0_2B == 0.0) {
886                             // Hydrogen position depends on other hydrogen, use generic location
887                             chiral = 0;
888                             assert atom2 != null;
889                             List<Bond> bonds2 = atom2.getBonds();
890                             atom2B =
891                                 (atom1 == bonds2.get(0).get1_2(atom2)) ? bonds2.get(1).get1_2(atom2)
892                                     : bonds2.get(0).get1_2(atom2);
893                             // Evenly space out hydrogen
894                             angle0_2B = 180.0 - 120.0 * knownCount;
895                           } else if (anglesList.size() == 2) {
896                             // Trigonal hydrogen use equipartition between selected atoms
897                             chiral = 3;
898                           } else if (knownCount == 1) {
899                             // Already created one hydrogen (chiral = 1), use other.
900                             chiral = -1;
901                           }
902                           //TODO discern whether to use chiral = 1 or -1 when angles are to three heavy atoms.
903                           intxyz(hydrogen, atom1, bond0.bondType.distance, atom2, angle0_2, atom2B,
904                               angle0_2B, chiral);
905                         }
906                       }
907                     }
908                   }
909                 } else {
910                   logger.info(format(" Substructure match failed (%d vs expected %d)", pLength, numInputMolAtoms - numMolHydrogen));
911                 }
912               }
913 
914               // Build molecule for output assembly matching this entity
915               MSGroup molecule;
916               if (mol instanceof Polymer) {
917                 molecule = new Polymer(((Polymer) mol).getChainID(), mol.getName(), true);
918               } else {
919                 molecule = new Molecule(mol.getName());
920               }
921               List<Atom> newAtomList = new ArrayList<>();
922               for (Atom atom : xyzatoms.get(i)) {
923                 Atom molAtom;
924                 if (molecule instanceof Polymer) {
925                   molAtom = new Atom(atomIndex++, atom.getName(), atom.getAltLoc(), atom.getXYZ(null),
926                       atom.getResidueName(), atom.getResidueNumber(), atom.getChainID(), atom.getOccupancy(),
927                       atom.getTempFactor(), atom.getSegID());
928                   molAtom.setAtomType(atom.getAtomType());
929                 } else {
930                   molAtom = new Atom(atomIndex++, atom.getName(), atom.getAtomType(), atom.getXYZ(null));
931                   if (atom.getResidueName() != null) {
932                     molAtom.setResName(atom.getResidueName());
933                   }
934                 }
935                 newAtomList.add(molAtom);
936               }
937               List<Bond> bondList = mol.getBondList();
938               for (Bond bond : bondList) {
939                 Atom a1 = bond.getAtom(0);
940                 Atom a2 = bond.getAtom(1);
941                 if (logger.isLoggable(Level.FINE)) {
942                   logger.fine(format(" Bonded atom 1: %d, Bonded atom 2: %d", a1.getXyzIndex(), a2.getXyzIndex()));
943                 }
944                 int molIndex = a1.getXyzIndex();
945                 shiftIndex = (molIndices[molIndex] == 0) ? 0 : shiftIndices[molIndices[molIndex] - 1];
946                 Atom newA1 = newAtomList.get(a1.getIndex() - shiftIndex - 1);
947                 Atom newA2 = newAtomList.get(a2.getIndex() - shiftIndex - 1);
948                 Bond bond2 = new Bond(newA1, newA2);
949                 bond2.setBondType(bond.getBondType());
950               }
951               Residue res = null;
952               for (Atom atom : newAtomList) {
953                 if (molecule instanceof Polymer) {
954                   Polymer polyMol = (Polymer) mol;
955                   if (res == null) {
956                     res = new Residue(atom.getResidueName(), atom.getResidueNumber(), polyMol.getResidue(atom.getResidueNumber()).getResidueType());
957                   } else if (res.getResidueNumber() != atom.getResidueNumber()) {
958                     if (logger.isLoggable(Level.FINER)) {
959                       logger.finer(" Added Residue " + res.getResidueNumber() + " " + res.getName());
960                     }
961                     molecule.addMSNode(res);
962                     res = new Residue(atom.getResidueName(), atom.getResidueNumber(), polyMol.getResidue(atom.getResidueNumber()).getResidueType());
963                   }
964                   res.addMSNode(atom);
965                 } else {
966                   molecule.addMSNode(atom);
967                 }
968               }
969               if (molecule instanceof Polymer && res != null) {
970                 if (logger.isLoggable(Level.FINER)) {
971                   logger.finer(" Added Final Residue " + res.getResidueNumber() + " " + res.getName());
972                 }
973                 molecule.addMSNode(res);
974               }
975               outputAssembly.addMSNode(molecule);
976             } else {
977               if (logger.isLoggable(Level.INFO)) {
978                 logger.info(
979                     format(" Number of atoms in CIF (%d) molecule do not match input (%d + %dH = %d).",
980                         cifMolAtoms, nInputAtoms - numMolHydrogen, numMolHydrogen, nInputAtoms));
981               }
982             }
983           }
984         }
985         if (logger.isLoggable(Level.FINE)) {
986           logger.fine(format("\n Output Assembly Atoms: %d", outputAssembly.getAtomList().size()));
987         }
988         outputAssembly.setPotential(assembly2.getPotentialEnergy());
989         if (assembly1.getCrystal() != null) {
990           outputAssembly.setCrystal(assembly1.getCrystal());
991         }
992         outputAssembly.setForceField(assembly2.getForceField());
993         outputAssembly.setFile(assembly1.getFile());
994         outputAssembly.setName(assembly1.getName());
995 
996         if (outputAssembly.getAtomList().size() < 1) {
997           logger.info(" Atom types could not be matched. File could not be written.");
998         } else if (!writeOutputFile(outputAssembly, systemFilter1)) {
999           logger.info(" Output assembly file could not be written.");
1000         }
1001         // Finished with this block. Clean up/reset variables for next block.
1002         outputAssembly.destroy();
1003       } while (systemFilter1.readNext());
1004     } else {
1005       logger.info(helpString());
1006     }
1007     return this;
1008   }
1009 
1010   /**
1011    * Write molecular assembly to an output file.
1012    *
1013    * @param ma           Assembly to save
1014    * @param systemFilter SystemFilter (for path and name base)
1015    * @return True if write succeeded.
1016    */
1017   private static boolean writeOutputFile(MolecularAssembly ma, SystemFilter systemFilter) {
1018     File file = systemFilter.getFiles().get(0);
1019     String dir = getFullPath(file.getAbsolutePath()) + File.separator;
1020     String fileName = removeExtension(getName(file.getAbsolutePath()));
1021     String spacegroup;
1022     if (ma.getCrystal() != null) {
1023       spacegroup = ma.getCrystal().getUnitCell().spaceGroup.shortName;
1024     } else {
1025       spacegroup = null;
1026     }
1027     List<MSNode> entities = ma.getAllBondedEntities();
1028     File saveFile;
1029     if (systemFilter instanceof PDBFilter) {
1030       // Change name for different space groups. Input cannot handle multiple space groups in same file.
1031       if (entities.size() > 1) {
1032         fileName += "_z" + entities.size();
1033       }
1034       saveFile = new File(dir + fileName + ".pdb");
1035     } else {
1036       if (systemFilter.countNumModels() > 1) {
1037         saveFile = new File(dir + fileName + "_reorder.arc");
1038       } else {
1039         // If only structure then create XYZ.
1040         saveFile = XYZFilter.version(new File(dir + fileName + ".xyz"));
1041       }
1042     }
1043     ForceField ff = ma.getForceField();
1044     CompositeConfiguration properties = ma.getProperties();
1045     if (systemFilter instanceof PDBFilter) {
1046       PDBFilter pdbFilter = new PDBFilter(saveFile, ma, ff, properties);
1047       pdbFilter.writeFile(saveFile, true);
1048     } else {
1049       XYZFilter xyzFilter = new XYZFilter(saveFile, ma, ff, properties);
1050       xyzFilter.writeFile(saveFile, true);
1051     }
1052     logger.info("\n Saved output file:        " + saveFile.getAbsolutePath());
1053     File propertyFile = new File(saveFile.getParentFile(), removeExtension(saveFile.getName()) + ".properties");
1054     if (!propertyFile.exists()) {
1055       try {
1056         FileWriter fw = new FileWriter(propertyFile, false);
1057         BufferedWriter bw = new BufferedWriter(fw);
1058         ForceField forceField = ma.getForceField();
1059         String parameters = forceField.getString("parameters", "none");
1060         if (parameters != null && !parameters.equalsIgnoreCase("none")) {
1061           bw.write(format("parameters %s\n", parameters));
1062         }
1063         String forceFieldProperty = forceField.getString("forcefield", "none");
1064         if (forceFieldProperty != null && !forceFieldProperty.equalsIgnoreCase("none")) {
1065           bw.write(format("forcefield %s\n", forceFieldProperty));
1066         }
1067         String patch = forceField.getString("patch", "none");
1068         if (patch != null && !patch.equalsIgnoreCase("none")) {
1069           bw.write(format("patch %s\n", patch));
1070         }
1071         bw.write(format("spacegroup %s\n", spacegroup));
1072         if (entities.size() > 1) {
1073           bw.write("intermolecular-softcore true\n");
1074         }
1075         logger.info("\n Saved properties file: " + propertyFile.getAbsolutePath() + "\n");
1076         bw.close();
1077       } catch (Exception ex) {
1078         logger.info("Failed to write files.\n" + Utilities.stackTraceToString(ex));
1079       }
1080     } else {
1081       logger.info("\n Property file already exists:  " + propertyFile.getAbsolutePath() + "\n");
1082     }
1083     return true;
1084   }
1085 
1086   /**
1087    * Compare atoms based on connectivity and atom types to determine order.
1088    */
1089   private final class AtomComparator implements Comparator<Atom> {
1090     @Override
1091     public int compare(Atom a1, Atom a2) {
1092       int comp = Integer.compare(a1.getMoleculeNumber(), a2.getMoleculeNumber());
1093       if (logger.isLoggable(Level.FINER) && comp != 0) {
1094         logger.finer(format(" Different Molecule Number (%d vs %d)", a1.getMoleculeNumber(), a2.getMoleculeNumber()));
1095       }
1096       if (comp == 0) {
1097         comp = Integer.compare(a1.getResidueNumber(), a2.getResidueNumber());
1098         if (logger.isLoggable(Level.FINER) && comp != 0) {
1099           logger.finer(format(" Different Residue Numbers %d vs %d", a1.getResidueNumber(), a2.getResidueNumber()));
1100         }
1101         if (comp == 0) {
1102           // Want heavier atoms first so multiply by -1
1103           comp = -Double.compare(a1.getMass(), a2.getMass());
1104           if (logger.isLoggable(Level.FINER) && comp != 0) {
1105             logger.finer(format(" Different Masses %6.3f %6.3f", a1.getMass(), a2.getMass()));
1106           }
1107           if (comp == 0) {
1108             comp = Integer.compare(a1.getAtomType().type, a2.getAtomType().type);
1109             if (logger.isLoggable(Level.FINER) && comp != 0) {
1110               logger.finer(format(" Different Atom Types (%d vs %d)", a1.getAtomType().type, a2.getAtomType().type));
1111             }
1112             if (comp == 0) {
1113               comp = Integer.compare(a1.getBonds().size(), a2.getBonds().size());
1114               if (logger.isLoggable(Level.FINER) && comp != 0) {
1115                 logger.finer(format(" Different Bond Sizes (%d vs %d)", a1.getBonds().size(), a2.getBonds().size()));
1116               }
1117               if (comp == 0) {
1118                 comp = Integer.compare(a1.get12List().size(), a2.get12List().size());
1119                 if (logger.isLoggable(Level.FINER) && comp != 0) {
1120                   logger.finer(format(" Different 1-2 Size (%d vs %d).", a1.get12List().size(), a2.get12List().size()));
1121                 }
1122                 if (comp == 0) {
1123                   for (int i = 0; i < a1.get12List().size(); i++) {
1124                     Atom a12 = a1.get12List().get(i);
1125                     Atom a22 = a2.get12List().get(i);
1126                     comp = shallowCompare(a12, a22);
1127                     if (comp != 0) {
1128                       break;
1129                     }
1130                   }
1131                   if (logger.isLoggable(Level.FINER) && comp != 0) {
1132                     logger.finer(" Different 1-2 Shallow.");
1133                   }
1134                   if (comp == 0) {
1135                     comp = Integer.compare(a1.get13List().size(), a2.get13List().size());
1136                     if (logger.isLoggable(Level.FINER) && comp != 0) {
1137                       logger.finer(" Different 1-3 Size.");
1138                     }
1139                     if (comp == 0) {
1140                       for (int i = 0; i < a1.get13List().size(); i++) {
1141                         Atom a13 = a1.get13List().get(i);
1142                         Atom a23 = a2.get13List().get(i);
1143                         comp = shallowCompare(a13, a23);
1144                         if (comp != 0) {
1145                           break;
1146                         }
1147                       }
1148                       if (logger.isLoggable(Level.FINER) && comp != 0) {
1149                         logger.finer(" Different 1-3 Shallow.");
1150                       }
1151                       if (comp == 0) {
1152                         comp = Integer.compare(a1.get14List().size(), a2.get14List().size());
1153                         if (logger.isLoggable(Level.FINER) && comp != 0) {
1154                           logger.finer(" Different 1-4 Size.");
1155                         }
1156                         if (comp == 0) {
1157                           for (int i = 0; i < a1.get14List().size(); i++) {
1158                             Atom a14 = a1.get14List().get(i);
1159                             Atom a24 = a2.get14List().get(i);
1160                             comp = shallowCompare(a14, a24);
1161                             if (comp != 0) {
1162                               break;
1163                             }
1164                           }
1165                           if (logger.isLoggable(Level.FINER) && comp != 0) {
1166                             logger.finer(" Different 1-4 Shallow.");
1167                           }
1168                           if (comp == 0) {
1169                             comp = Integer.compare(a1.get15List().size(), a2.get15List().size());
1170                             if (logger.isLoggable(Level.FINER) && comp != 0) {
1171                               logger.finer(" Different 1-5 Size.");
1172                             }
1173                             if (comp == 0) {
1174                               for (int i = 0; i < a1.get15List().size(); i++) {
1175                                 Atom a15 = a1.get15List().get(i);
1176                                 Atom a25 = a2.get15List().get(i);
1177                                 comp = shallowCompare(a15, a25);
1178                                 if (comp != 0) {
1179                                   break;
1180                                 }
1181                               }
1182                               if (logger.isLoggable(Level.FINER) && comp != 0) {
1183                                 logger.finer(" Different 1-5 Shallow.");
1184                               }
1185                             }
1186                           }
1187                         }
1188                       }
1189                     }
1190                   }
1191                 }
1192               }
1193             }
1194           }
1195         }
1196       }
1197       return comp;
1198     }
1199 
1200     /**
1201      * Compare two atoms based on their bonded atoms' connections.
1202      *
1203      * @param a1 Atom 1
1204      * @param a2 Atom 2
1205      * @return 0 is uncertain, -1 if reverse, 1 if same.
1206      */
1207     private int shallowCompare(Atom a1, Atom a2) {
1208       int comp = Integer.compare(a1.getBonds().size(), a2.getBonds().size());
1209       if (logger.isLoggable(Level.FINER) && comp != 0) {
1210         logger.finer(" Different Shallow Bond Size.");
1211       }
1212       if (comp == 0) {
1213         comp = Integer.compare(a1.getAtomType().type, a2.getAtomType().type);
1214         if (logger.isLoggable(Level.FINER) && comp != 0) {
1215           logger.finer(" Different Shallow Atom Type.");
1216         }
1217         if (comp == 0) {
1218           comp = Integer.compare(a1.get12List().size(), a2.get12List().size());
1219           if (logger.isLoggable(Level.FINER) && comp != 0) {
1220             logger.finer(" Different Shallow 1-2 Size.");
1221           }
1222           if (comp == 0) {
1223             comp = Integer.compare(a1.get13List().size(), a2.get13List().size());
1224             if (logger.isLoggable(Level.FINER) && comp != 0) {
1225               logger.finer(" Different Shallow 1-3 Size.");
1226             }
1227             if (comp == 0) {
1228               comp = Integer.compare(a1.get14List().size(), a2.get14List().size());
1229               if (logger.isLoggable(Level.FINER) && comp != 0) {
1230                 logger.finer(" Different Shallow 1-4 Size.");
1231               }
1232               if (comp == 0) {
1233                 comp = Integer.compare(a1.get15List().size(), a2.get15List().size());
1234                 if (logger.isLoggable(Level.FINER) && comp != 0) {
1235                   logger.finer(" Different Shallow 1-5 Size.");
1236                 }
1237               }
1238             }
1239           }
1240         }
1241       }
1242       return comp;
1243     }
1244   }
1245 }