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-2024.
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;
39  
40  import static ffx.numerics.math.DoubleMath.sub;
41  import static ffx.potential.bonded.NamingUtils.renameAminoAcidToPDBStandard;
42  import static ffx.potential.bonded.NamingUtils.renameNucleicAcidToPDBStandard;
43  import static org.apache.commons.math3.util.FastMath.sqrt;
44  
45  import ffx.potential.bonded.AminoAcidUtils;
46  import ffx.potential.bonded.Atom;
47  import ffx.potential.bonded.Bond;
48  import ffx.potential.bonded.Molecule;
49  import ffx.potential.bonded.NucleicAcidUtils;
50  import ffx.potential.bonded.Polymer;
51  import ffx.potential.bonded.Residue;
52  import java.io.ByteArrayOutputStream;
53  import java.io.PrintStream;
54  import java.nio.charset.StandardCharsets;
55  import java.util.ArrayList;
56  import java.util.List;
57  import java.util.ListIterator;
58  import java.util.logging.Logger;
59  
60  /**
61   * The Utilities class provides methods to locate functional units of an organic system.
62   *
63   * @author Michael J. Schnieders
64   * @since 1.0
65   */
66  public final class Utilities {
67  
68    private static final Logger logger = Logger.getLogger(Utilities.class.getName());
69    /**
70     * The algorithms used to split arrays of atoms into a structural hierarchy use a bunch of List
71     * instances. Pooling them seems like it might be a performance win - although better algorithms
72     * probably exist. This is currently backed by ArrayLists.
73     */
74    private static final List<List<Atom>> atomListPool = new ArrayList<>();
75  
76    /**
77     * Finds the RMS deviation between the atoms of MolecularAssembly m1 and m2 provided they have the
78     * same number of atoms.
79     *
80     * @param m1 a {@link ffx.potential.MolecularAssembly} object.
81     * @param m2 a {@link ffx.potential.MolecularAssembly} object.
82     * @return a double.
83     */
84    public static double RMSCoordDev(MolecularAssembly m1, MolecularAssembly m2) {
85      if (m1 == null || m2 == null) {
86        return 0;
87      }
88      int n1 = m1.getAtomList().size();
89      int n2 = m2.getAtomList().size();
90      if (n1 != n2) {
91        return 0;
92      }
93      Atom a1, a2;
94      double[] d = new double[3];
95      double[] da = new double[3];
96      double[] db = new double[3];
97      double rms = 0;
98      ListIterator<Atom> li, lj;
99      for (li = m1.getAtomList().listIterator(), lj = m2.getAtomList().listIterator();
100         li.hasNext(); ) {
101       a1 = li.next();
102       a2 = lj.next();
103       a1.getXYZ(da);
104       a2.getXYZ(db);
105       sub(da, db, d);
106       rms += d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
107     }
108     return sqrt(rms / n1);
109   }
110 
111   /**
112    * This routine sub-divides a system into groups of ions, water, hetero molecules, and
113    * polynucleotides/polypeptides.
114    *
115    * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
116    * @param atoms a {@link java.util.List} object.
117    */
118   public static void biochemistry(MolecularAssembly molecularAssembly, List<Atom> atoms) {
119 
120     // logger.info(" Biochemistry called.");
121 
122     Atom atom, seed = null;
123     int num = 0;
124     int waterNum = 0;
125     int ionNum = 0;
126     int moleculeNum = 0;
127     List<String> segIDs = new ArrayList<>();
128     while (!atoms.isEmpty()) {
129       /*
130        Nitrogen is used to "seed" a backbone search because carbon can
131        be separated from the backbone by a sulfur (e.g., MET).
132       */
133       for (Atom a : atoms) {
134         seed = a;
135         if (seed.getAtomicNumber() == 7) {
136           break;
137         }
138       }
139       // If no nitrogen atoms remain, there are no nucleic or amino acids.
140       if (seed.getAtomicNumber() != 7) {
141         // logger.info(" Finished searching for proteins and nucleic acids.");
142         List<Atom> moleculeAtoms;
143         while (!atoms.isEmpty()) {
144           atom = atoms.get(0);
145           if (atom.getNumBonds() == 0) {
146             // A metal ion or noble gas
147             ionNum++;
148             Molecule ion = new Molecule(atom.getName() + "-" + ionNum);
149             ion.addMSNode(atom);
150             atoms.remove(0);
151             molecularAssembly.addMSNode(ion);
152             continue;
153           } else if (atom.getAtomicNumber() == 8 && isWaterOxygen(atom)) {
154             // Water
155             waterNum++;
156             Molecule water = new Molecule("Water-" + waterNum);
157             water.addMSNode(atom);
158             atoms.remove(0);
159             List<Bond> bonds = atom.getBonds();
160             for (Bond b : bonds) {
161               Atom hydrogen = b.get1_2(atom);
162               water.addMSNode(hydrogen);
163               atoms.remove(hydrogen);
164             }
165             molecularAssembly.addMSNode(water);
166             continue;
167           }
168           // All other molecules
169           moleculeNum++;
170           // logger.info(" Molecule: " + moleculeNum);
171           Molecule molecule = new Molecule("Molecule-" + moleculeNum);
172           moleculeAtoms = getAtomListFromPool();
173           collectAtoms(atoms.get(0), moleculeAtoms, true);
174           while (!moleculeAtoms.isEmpty()) {
175             atom = moleculeAtoms.get(0);
176             moleculeAtoms.remove(0);
177             molecule.addMSNode(atom);
178             atoms.remove(atom);
179           }
180           molecularAssembly.addMSNode(molecule);
181         }
182         break;
183       }
184 
185       List<Atom> backbone = findPolymer(seed, null);
186       if (!backbone.isEmpty()) {
187         seed = backbone.get(backbone.size() - 1);
188         backbone = findPolymer(seed, null);
189       }
190 
191       // logger.info(" Backbone length: " + backbone.size());
192       Character chainID = getChainID(num);
193       String segID = getSegID(chainID, segIDs);
194       Polymer c = new Polymer(chainID, segID, true);
195       if (backbone.size() > 2 && divideBackbone(backbone, c)) {
196         for (Atom a : c.getAtomList()) {
197           atoms.remove(a);
198         }
199         // logger.info(" Sequenced chain: " + c.getName());
200         molecularAssembly.addMSNode(c);
201         num++;
202       } else {
203         moleculeNum++;
204         // logger.info(" Molecule: " + moleculeNum);
205         Molecule hetero = new Molecule("Molecule-" + moleculeNum);
206         atom = backbone.get(0);
207         List<Atom> heteroAtomList = getAtomListFromPool();
208         collectAtoms(atom, heteroAtomList,  true);
209         for (Atom a : heteroAtomList) {
210           hetero.addMSNode(a);
211         }
212         for (Atom a : hetero.getAtomList()) {
213           atoms.remove(a);
214         }
215         molecularAssembly.addMSNode(hetero);
216       }
217     }
218   }
219 
220   /**
221    * Returns an atom bonded to the "end" atom, which is not equal to "other".
222    *
223    * @param end Atom
224    * @param other Atom
225    * @return Atom
226    */
227   public static Atom findSeed(Atom end, Atom other) {
228     for (Bond b : end.getBonds()) {
229       Atom seed = b.get1_2(end);
230       if (seed != other) {
231         return seed;
232       }
233     }
234     return null;
235   }
236 
237   /**
238    * Determine chainID for a given polymer number.
239    *
240    * @param i int
241    * @return Character
242    */
243   public static Character getChainID(int i) {
244     if (i > 35) {
245       i = i % 36;
246     }
247     char c;
248     if (i < 26) {
249       // 65 is 'A'. 90 is 'Z'.
250       c = (char) (i + 65);
251     } else {
252       i -= 26;
253       // 48 is '0'. 57 is '9'.
254       c = (char) (i + 48);
255     }
256     return c;
257   }
258 
259   /**
260    * Gets the stack trace for an exception and converts it to a String.
261    *
262    * @param ex An exception.
263    * @return A String of its stack trace.
264    */
265   public static String stackTraceToString(Throwable ex) {
266     ByteArrayOutputStream byteArrayOutputStream = new ByteArrayOutputStream();
267     try (PrintStream ps = new PrintStream(byteArrayOutputStream, true, StandardCharsets.UTF_8)) {
268       ex.printStackTrace(ps);
269       return byteArrayOutputStream.toString(StandardCharsets.UTF_8);
270     } catch (Exception e) {
271       logger.warning("Unable to convert stack trace to String");
272       return "";
273     }
274   }
275 
276   /**
277    * addAtomListToPool
278    *
279    * @param a a {@link java.util.List} object.
280    */
281   private static void addAtomListToPool(List<Atom> a) {
282     a.clear();
283     atomListPool.add(a);
284   }
285 
286   /**
287    * Collect all the atoms in a polymer opposite the end atom, and put them into the residue. This
288    * assumes that the "cap" is only linked to the rest of the polymer through the end atom.
289    *
290    * @param end Atom
291    * @param seed Atom
292    * @param residue Residue
293    */
294   private static void addCap(Atom end, Atom seed, Residue residue) {
295     List<Atom> cap = getAtomListFromPool();
296     cap.add(end);
297     collectAtoms(seed, cap);
298     // Assume the end & seed atoms are already part of the residue
299     cap.remove(0);
300     cap.remove(0);
301     for (Atom a : cap) {
302       residue.addMSNode(a);
303     }
304   }
305 
306   /**
307    * Add a phosphate and its bonded oxygens that are not bonded to a carbon to the specified
308    * residue.
309    *
310    * @param phosphate Atom
311    * @param residue Residue
312    */
313   private static void addPhosphate(Atom phosphate, Residue residue) {
314     if (phosphate == null) {
315       return;
316     }
317     residue.addMSNode(phosphate);
318     for (Bond b : phosphate.getBonds()) {
319       Atom oxygen = b.get1_2(phosphate);
320       // Add oxygens not bonded to a Carbon
321       if (numberOfBondsWith(oxygen, 6) == 0) {
322         residue.addMSNode(oxygen);
323         // Add hydrogens atoms for protonated oxygen groups
324         Atom hydrogen = findBondWith(oxygen, 1);
325         if (hydrogen != null) {
326           residue.addMSNode(hydrogen);
327         }
328       }
329     }
330   }
331 
332   private static Residue assignResidue(
333       List<Atom> backbone, int start, List<Atom> atoms, List<Atom> sidePolymer) {
334     int atomicnum;
335     char[] chars = {'S', 'P', 'O', 'N', 'C'};
336     // Bin number to atomic number => 0 = S, 1 = P, 2 = O, 3 = N, 4 = C
337     int[] bins = new int[5];
338 
339     for (Atom a : sidePolymer) {
340       atomicnum = a.getAtomicNumber();
341       switch (atomicnum) {
342         case 1:
343           // ignore hydrogens
344           break;
345         case 6:
346           // Carbon
347           bins[4]++;
348           break;
349         case 7:
350           // Nitrogen
351           bins[3]++;
352           break;
353         case 8:
354           // Oxygen
355           bins[2]++;
356           break;
357         case 15:
358           // Phosphorus
359           bins[1]++;
360           break;
361         case 16:
362           // Sulfur
363           bins[0]++;
364           break;
365         default:
366           return null;
367       }
368     }
369     StringBuilder key = new StringBuilder();
370     int atomCount = 0;
371     for (int i = 0; i < 5; i++) {
372       if (bins[i] != 0) {
373         atomCount += bins[i];
374         key.append(chars[i]);
375         key.append(bins[i]);
376       }
377     }
378     if (atomCount == 0) {
379       key.append("H"); // Glycine
380     }
381 
382     String resname = AminoAcidUtils.sidechainStoichiometry.get(key.toString());
383 
384     if (resname == null) {
385       resname = "Unknown";
386     } else {
387       resname = resname.intern();
388     }
389 
390     if (resname.equals("1") || resname.equals("2")) {
391       // Special case where atom string keys aren't unique
392       Atom alpha = backbone.get(start + 1);
393       Atom carbonyl = backbone.get(start + 2);
394       Atom beta = null;
395       List<Bond> alphabonds = alpha.getBonds();
396       for (Bond bond : alphabonds) {
397         beta = bond.get1_2(alpha);
398         // Don't want the peptide nitrogen or alpha hydrogen or carbonyl carbon
399         if (beta.getAtomicNumber() != 7 && beta.getAtomicNumber() != 1 && beta != carbonyl) {
400           break;
401         }
402         beta = null;
403       }
404       if (beta == null) {
405         return null;
406       }
407       List<Bond> betabonds = beta.getBonds();
408       Atom gamma;
409       int carboncount = 0;
410       for (Bond bond : betabonds) {
411         gamma = bond.get1_2(beta);
412         if (gamma.getAtomicNumber() == 6) {
413           carboncount++;
414         }
415       }
416       if (resname.equals("1")) {
417         if (carboncount == 2) {
418           resname = "PRO";
419         } else {
420           resname = "VAL";
421         }
422       } else {
423         if (carboncount == 2) {
424           resname = "LEU";
425         } else {
426           resname = "ILE";
427         }
428       }
429     } else if (resname.equals("3")) {
430       Atom c3 = backbone.get(start + 3);
431       int num = countCO(c3);
432       if (num == 2) {
433         resname = "A";
434       } else {
435         resname = "DG";
436       }
437     }
438     Residue residue = null;
439     try {
440       NucleicAcidUtils.NucleicAcid3.valueOf(resname.toUpperCase());
441       residue = new Residue(resname, Residue.ResidueType.NA);
442     } catch (Exception e) {
443       //
444     }
445     if (residue == null) {
446       try {
447         AminoAcidUtils.AminoAcid3.valueOf(resname.toUpperCase());
448         residue = new Residue(resname, Residue.ResidueType.AA);
449       } catch (Exception e) {
450         //
451       }
452     }
453     if (residue == null) {
454       residue = new Residue(resname, Residue.ResidueType.UNK);
455     }
456 
457     // Create the Residue group
458     for (Atom a : atoms) {
459       residue.addMSNode(a);
460       // logger.info(a.toString());
461     }
462     return residue;
463   }
464 
465   /**
466    * Given an array of bonded atoms, this function recursively collects all other connected atoms,
467    * without backtracking over atoms already in the list. Disulfide bonds are not crossed. (the
468    * intent is to search along a peptide backbone) Atoms preloaded into the List provide search
469    * termination.
470    *
471    * @param seed Atom
472    * @param atoms List
473    */
474   private static void collectAtoms(Atom seed, List<Atom> atoms) {
475     collectAtoms(seed, atoms, false);
476   }
477 
478   /**
479    * Given an array of bonded atoms, this function recursively collects all other connected atoms,
480    * without backtracking over atoms already in the list. Disulfide bonds can be crossed. (the
481    * intent is to not cross disulfide bonds in a protein, but cross disulfide bonds within molecules)
482    * Atoms preloaded into the List provide search termination.
483    *
484    * @param seed Atom
485    * @param atoms List
486    * @param searchDisulfide Allows crossing of disulfide bonds (e.g. molecules)
487    */
488   private static void collectAtoms(Atom seed, List<Atom> atoms, boolean searchDisulfide) {
489     if (seed == null) {
490       return;
491     }
492     atoms.add(seed);
493     for (Bond b : seed.getBonds()) {
494       Atom nextAtom = b.get1_2(seed);
495       if (nextAtom.getParent() != null) {
496         continue;
497       }
498       // If true, search across disulfide bonds.
499       if ((searchDisulfide || (nextAtom.getAtomicNumber() != 16 || seed.getAtomicNumber() != 16))
500           && !atoms.contains(nextAtom)) {
501         collectAtoms(nextAtom, atoms, searchDisulfide);
502       }
503     }
504   }
505 
506   /**
507    * countCO
508    *
509    * @param adjacent a {@link ffx.potential.bonded.Atom} object.
510    * @return a int.
511    */
512   private static int countCO(Atom adjacent) {
513     int total = 0;
514     for (Bond b : adjacent.getBonds()) {
515       Atom carbonyl = b.get1_2(adjacent);
516       if (carbonyl.getAtomicNumber() == 6) {
517         for (Bond b2 : carbonyl.getBonds()) {
518           Atom oxygen = b2.get1_2(carbonyl);
519           if (oxygen.getAtomicNumber() == 8) {
520             total++;
521           }
522         }
523       }
524     }
525     return total;
526   }
527 
528   /**
529    * divideBackbone
530    *
531    * @param backbone a {@link java.util.List} object.
532    * @param c a {@link ffx.potential.bonded.Polymer} object.
533    * @return a boolean.
534    */
535   private static boolean divideBackbone(List<Atom> backbone, Polymer c) {
536     int length = backbone.size();
537 
538     // Try to find a Phosphorus or Nitrogen in the backbone
539     int nitrogenCount = 0;
540     int phosphorusCount = 0;
541     for (Atom match : backbone) {
542       int atomicNumber = match.getAtomicNumber();
543       if (atomicNumber == 15) {
544         phosphorusCount++;
545       } else if (atomicNumber == 7) {
546         nitrogenCount++;
547       }
548     }
549 
550     // logger.info(format(" Divide Backbone: %d Nitrogen", nitrogenCount));
551     // logger.info(format(" Divide Backbone: %d Phosphorous", phosphorusCount));
552 
553     PolymerType type;
554     if (phosphorusCount >= nitrogenCount && phosphorusCount > 1) {
555       // logger.info(" Divide Backbone: Nucleic Acid");
556       type = PolymerType.NUCLEICACID;
557     } else if (nitrogenCount > phosphorusCount && nitrogenCount > 2) {
558       // logger.info(" Divide Backbone: Amino Acid");
559       type = PolymerType.AMINOACID;
560     } else {
561       // logger.info(" Divide Backbone: Unknown");
562       return false;
563     }
564     int start = -1;
565     for (int i = 0; i < length; i++) {
566       Residue res = patternMatch(i, backbone, type);
567       if (res != null) {
568         for (Atom a : res.getAtomList()) {
569           a.setParent(null);
570         }
571         if (!(res.getName().equals("Unknown"))) {
572           start = i;
573           // Want 5' to 3'
574           if (type == PolymerType.NUCLEICACID) {
575             Atom carbon5 = backbone.get(start + 1);
576             if (numberOfBondsWith(carbon5, 6) != 1) {
577               start = -1;
578             }
579           }
580           break;
581         }
582       }
583     }
584     if (start == -1) {
585       backbone = reverseAtomList(backbone);
586       for (int i = 0; i < length; i++) {
587         Residue res = patternMatch(i, backbone, type);
588         if (res != null) {
589           for (Atom a : res.getAtomList()) {
590             a.setParent(null);
591           }
592           if (!(res.getName().equals("Unknown"))) {
593             start = i;
594             break;
595           }
596         }
597       }
598     }
599     if (start == -1) {
600       return false;
601     }
602     // Potential Polypeptide
603     if (type == PolymerType.AMINOACID) {
604       Atom nitrogen, alpha, carbonyl = null;
605       Atom nitrogen2, carbonyl2;
606       Residue aa = null;
607       int lastRes = 0;
608       int firstRes = -1;
609       List<Residue> aaArray = new ArrayList<>();
610       while (start < length) {
611         aa = patternMatch(start, backbone, PolymerType.AMINOACID);
612         if (aa != null) {
613           if (firstRes == -1) {
614             firstRes = start;
615             carbonyl = findCarbonyl(backbone.get(start));
616           }
617           aaArray.add(aa);
618           lastRes = start;
619         }
620         start += 3;
621       }
622       // Make sure the first residue is found
623       aa = null;
624       if (carbonyl != null) {
625         alpha = findAlphaCarbon(carbonyl);
626         if (alpha != null) {
627           nitrogen = findBondWith(alpha, 7);
628           if (nitrogen != null) {
629             nitrogen2 = findBondWith(carbonyl, 7);
630             List<Atom> firstAtoms = getAtomListFromPool();
631             firstAtoms.add(nitrogen);
632             firstAtoms.add(alpha);
633             firstAtoms.add(carbonyl);
634             firstAtoms.add(nitrogen2);
635             aa = patternMatch(0, firstAtoms, PolymerType.AMINOACID);
636             addAtomListToPool(firstAtoms);
637             if (aa != null) {
638               addCap(alpha, nitrogen, aa);
639               aaArray.add(0, aa);
640             }
641           }
642         }
643       }
644       // Add the remaining atoms to the end of the Polymer
645       if (aa == null) {
646         nitrogen = backbone.get(firstRes);
647         alpha = backbone.get(firstRes + 1);
648         addCap(alpha, nitrogen, aaArray.get(0));
649       }
650       // Make sure the last residue is found
651       aa = null;
652       carbonyl = findCarbonyl(backbone.get(lastRes + 1));
653       if (carbonyl != null) {
654         nitrogen = findBondWith(carbonyl, 7);
655         if (nitrogen != null) {
656           alpha = findAlphaCarbon(nitrogen);
657           if (alpha != null) {
658             carbonyl2 = findCarbonyl(alpha);
659             if (carbonyl2 != null) {
660               List<Atom> lastAtoms = getAtomListFromPool();
661               lastAtoms.add(carbonyl);
662               lastAtoms.add(nitrogen);
663               lastAtoms.add(alpha);
664               lastAtoms.add(carbonyl2);
665               aa = patternMatch(1, lastAtoms, PolymerType.AMINOACID);
666               addAtomListToPool(lastAtoms);
667               if (aa != null) {
668                 addCap(alpha, carbonyl2, aa);
669                 aaArray.add(aa);
670               }
671             }
672           }
673         }
674       }
675       if (aa == null) {
676         carbonyl = backbone.get(lastRes + 2);
677         alpha = backbone.get(lastRes + 1);
678         addCap(alpha, carbonyl, aaArray.get(aaArray.size() - 1));
679       }
680       int index = 1;
681       for (Residue r : aaArray) {
682         r.setNumber(index++);
683         if(!renameAminoAcidToPDBStandard(r)){
684           return false;
685         }
686         c.addMSNode(r);
687       }
688       // Potential DNA/RNA
689     } else if (type == PolymerType.NUCLEICACID) {
690       Residue base;
691       int lastRes = 0;
692       boolean firstBase = true;
693       Atom phos = null;
694       Atom oxygen1 = null;
695       Atom phosphate1 = null;
696       List<Residue> na = new ArrayList<>();
697       while (start < length) {
698         base = patternMatch(start, backbone, PolymerType.NUCLEICACID);
699         if (base != null) {
700           phos = backbone.get(start - 1);
701           if (phos != null && phos.getAtomicNumber() == 15) {
702             addPhosphate(phos, base);
703           }
704           na.add(base);
705           if (firstBase) {
706             firstBase = false;
707             phosphate1 = backbone.get(start - 1);
708             oxygen1 = backbone.get(start);
709           }
710           lastRes = start;
711         }
712         start += 6;
713       }
714       // Make sure the first base is found
715       Atom o2, o3;
716       Atom c1, c2, c3;
717       if (phosphate1 != null && oxygen1 != null) {
718         o2 = findOtherOxygen(phosphate1, oxygen1);
719         if (o2 != null) {
720           c1 = findBondWith(o2, 6);
721           if (c1 != null) {
722             c2 = findCO(c1);
723             if (c2 != null) {
724               c3 = findC5(c2);
725               if (c3 != null) {
726                 o3 = findBondWith(c3, 8);
727                 if (o3 != null) {
728                   List<Atom> firstAtoms = getAtomListFromPool();
729                   firstAtoms.add(o3);
730                   firstAtoms.add(c3);
731                   firstAtoms.add(c2);
732                   firstAtoms.add(c1);
733                   firstAtoms.add(o2);
734                   firstAtoms.add(phosphate1);
735                   base = patternMatch(0, firstAtoms, type);
736                   if (base != null) {
737                     addCap(c3, o3, base);
738                     na.add(0, base);
739                   }
740                 }
741               }
742             }
743           }
744         }
745       }
746       // Make sure the last base is found
747       oxygen1 = backbone.get(lastRes + 4);
748       phosphate1 = backbone.get(lastRes + 5);
749       if (phosphate1 != null && oxygen1 != null) {
750         o2 = findOtherOxygen(phosphate1, oxygen1);
751         if (o2 != null) {
752           c1 = findBondWith(o2, 6);
753           if (c1 != null) {
754             c2 = findBondWith(c1, 6);
755             if (c2 != null) {
756               c3 = findCCO(c2);
757               if (c3 != null) {
758                 o3 = findBondWith(c3, 8);
759                 if (o3 != null) {
760                   List<Atom> lastAtoms = getAtomListFromPool();
761                   lastAtoms.add(phosphate1);
762                   lastAtoms.add(o2);
763                   lastAtoms.add(c1);
764                   lastAtoms.add(c2);
765                   lastAtoms.add(c3);
766                   lastAtoms.add(o3);
767                   base = patternMatch(1, lastAtoms, type);
768                   if (base != null) {
769                     addPhosphate(phosphate1, base);
770                     addCap(c3, o3, base);
771                     na.add(base);
772                   }
773                 }
774               }
775             }
776           }
777         }
778       }
779       int index = 1;
780       for (Residue r : na) {
781         r.setNumber(index++);
782         renameNucleicAcidToPDBStandard(r);
783         c.addMSNode(r);
784       }
785     } else {
786       return false;
787     }
788     return true;
789   }
790 
791   /**
792    * Returns a carbon that is bonded to the atom a, a carbonyl group, and a nitrogen. O=C-alpha-N
793    *
794    * @param a Atom
795    * @return Atom
796    */
797   private static Atom findAlphaCarbon(Atom a) {
798     for (Bond b : a.getBonds()) {
799       Atom alpha = b.get1_2(a);
800       if (alpha.getAtomicNumber() == 6 && findCO(alpha) != null && formsBondsWith(alpha, 7)) {
801         return alpha;
802       }
803     }
804     return null;
805   }
806 
807   /**
808    * Returns the first atom with the specified atomic number that bonds with atom a, or null
809    * otherwise.
810    *
811    * @param a Atom
812    * @param atomicNumber int
813    * @return Atom
814    */
815   private static Atom findBondWith(Atom a, int atomicNumber) {
816     for (Bond b : a.getBonds()) {
817       Atom other = b.get1_2(a);
818       if (other.getAtomicNumber() == atomicNumber) {
819         return other;
820       }
821     }
822     return null;
823   }
824 
825   /**
826    * Returns a carbon that is bonded to the adjacent atom, which bonds 1 carbon and 1 oxygen.
827    *
828    * @param adjacent Atom
829    * @return Atom
830    */
831   private static Atom findC5(Atom adjacent) {
832     for (Bond b : adjacent.getBonds()) {
833       Atom carbon = b.get1_2(adjacent);
834       if (carbon.getAtomicNumber() == 6
835           && numberOfBondsWith(carbon, 6) == 1
836           && numberOfBondsWith(carbon, 8) == 1) {
837         return carbon;
838       }
839     }
840     return null;
841   }
842 
843   /**
844    * Returns a carbon that is bonded to the adjacent atom, which bonds 2 carbons and 1 oxygen.
845    *
846    * @param adjacent Atom
847    * @return Atom
848    */
849   private static Atom findCCO(Atom adjacent) {
850     for (Bond b : adjacent.getBonds()) {
851       Atom carbon = b.get1_2(adjacent);
852       if (carbon.getAtomicNumber() == 6
853           && numberOfBondsWith(carbon, 6) == 2
854           && numberOfBondsWith(carbon, 8) >= 1) {
855         return carbon;
856       }
857     }
858     return null;
859   }
860 
861   /**
862    * Returns a carbon that is bonded to the adjacent atom, which bonds at least 1 oxygen.
863    *
864    * @param adjacent Atom
865    * @return Atom
866    */
867   private static Atom findCO(Atom adjacent) {
868     for (Bond b : adjacent.getBonds()) {
869       Atom carbon = b.get1_2(adjacent);
870       if (carbon.getAtomicNumber() == 6 && formsBondsWith(carbon, 8)) {
871         return carbon;
872       }
873     }
874     return null;
875   }
876 
877   /**
878    * Returns a carbon that is bonded to the adjacent atom and an oxygen.
879    *
880    * @param adjacent Atom
881    * @return Atom
882    */
883   private static Atom findCarbonyl(Atom adjacent) {
884     for (Bond b : adjacent.getBonds()) {
885       Atom carbonyl = b.get1_2(adjacent);
886       if (carbonyl.getAtomicNumber() == 6) {
887         for (Bond b2 : carbonyl.getBonds()) {
888           Atom oxygen = b2.get1_2(carbonyl);
889           if (oxygen.getAtomicNumber() == 8 && oxygen.getBonds().size() == 1) {
890             return carbonyl;
891           }
892         }
893       }
894     }
895     return null;
896   }
897 
898   /**
899    * Returns an oxygen atom that is bonded to atom p and a carbon, but is not atom o. O-P-X-C where
900    * X is the returned atom. This is useful for traversing a nucleic acid backbone.
901    *
902    * @param p Atom
903    * @param o Atom
904    * @return Atom
905    */
906   private static Atom findOtherOxygen(Atom p, Atom o) {
907     for (Bond b : p.getBonds()) {
908       Atom oxygen = b.get1_2(p);
909       if (oxygen.getAtomicNumber() == 8 && oxygen != o && formsBondsWith(oxygen, 6)) {
910         return oxygen;
911       }
912     }
913     return null;
914   }
915 
916   /**
917    * findPolymer
918    *
919    * @param currentAtom Atom
920    * @param path List
921    * @return List
922    */
923   private static List<Atom> findPolymer(Atom currentAtom, List<Atom> path) {
924 
925     // Atom has no bonds to follow
926     if (currentAtom.getBonds() == null) {
927       path = getAtomListFromPool();
928       path.add(currentAtom);
929       return path;
930     }
931 
932     // End of Recursion condition
933     if (currentAtom.getParent() != null) {
934       return null;
935     }
936 
937     // Only C,N,O,P in a DNA/RNA/protein backbone
938     int anum = currentAtom.getAtomicNumber();
939     if (anum != 6 && anum != 7 && anum != 8 && anum != 15) {
940       return null;
941     }
942 
943     // Allow the search to make it out of side chains, but not enter them.
944     if (path != null && path.size() > 6) {
945 
946       // Oxygen is only in the backbone for Nucleic Acids in a phosphate group
947       if (anum == 8) {
948         if (!formsBondsWith(currentAtom, 15)) {
949           // Found an oxygen not bonded to a phosphate.
950           return null;
951         }
952         // Nitrogen is only in the backbone in peptide bonds or at an N-terminus.
953       } else if (anum == 7) {
954         Atom alphaCarbon = findAlphaCarbon(currentAtom);
955         if (alphaCarbon == null) {
956           return null;
957         }
958         // Avoid more than 3 carbons in a row (phenyl groups, etc.)
959       } else if (anum == 6) {
960         Atom a;
961         int anum2, anum3;
962         int size = path.size();
963         a = path.get(size - 1);
964         anum2 = a.getAtomicNumber();
965         if (anum2 == 6) {
966           a = path.get(size - 2);
967           anum3 = a.getAtomicNumber();
968           if (anum3 == 6) {
969             return null;
970           }
971         }
972       }
973     }
974 
975     Atom previousAtom = null;
976     if (path != null) {
977       previousAtom = path.get(path.size() - 1);
978     }
979 
980     // Atoms with only one bond are at the end of a Polymer
981     List<Bond> bonds = currentAtom.getBonds();
982     if (bonds.size() == 1 && previousAtom != null) {
983       return null;
984     }
985 
986     // Initialization
987     if (path == null) {
988       path = getAtomListFromPool();
989       previousAtom = null;
990       // Or Continuation
991     } else {
992       List<Atom> pathclone = getAtomListFromPool();
993       pathclone.addAll(path);
994       path = pathclone;
995     }
996 
997     // Add the currentAtom to the growing path
998     path.add(currentAtom);
999 
1000     // Continue search in each bond direction, but no backtracking over previousAtom
1001     List<Atom> newPolymer, maxPolymer = getAtomListFromPool();
1002     for (Bond b : bonds) {
1003       Atom nextAtom = b.get1_2(currentAtom);
1004 
1005       // Check to avoid returning in the same direction and loops
1006       if (nextAtom != previousAtom && !path.contains(nextAtom)) {
1007         newPolymer = findPolymer(nextAtom, path);
1008         if (newPolymer != null) {
1009           // logger.info(" Next atom: " + nextAtom + " Size: " + newPolymer.size());
1010           // Check to see if the Polymers contain any of the same atoms,
1011           // and if so, use the shorter Polymer (avoids loops)
1012           // if (haveCommonAtom(newPolymer, maxPolymer)) {
1013           //    if (newPolymer.size() < maxPolymer.size()) {
1014           //        addAtomListToPool(maxPolymer);
1015           //        maxPolymer = newPolymer;
1016           //    }
1017           // } else if (newPolymer.size() > maxPolymer.size()) {
1018           if (newPolymer.size() > maxPolymer.size()) {
1019             addAtomListToPool(maxPolymer);
1020             maxPolymer = newPolymer;
1021           }
1022         }
1023       }
1024     }
1025 
1026     // Add the currentAtom to the longest discovered chain and return
1027     maxPolymer.add(0, currentAtom);
1028     return maxPolymer;
1029   }
1030 
1031   /**
1032    * True if Atom a forms a bond with another atom of the specified atomic Number.
1033    *
1034    * @param a Atom
1035    * @param atomicNumber int
1036    * @return boolean
1037    */
1038   private static boolean formsBondsWith(Atom a, int atomicNumber) {
1039     for (Bond b : a.getBonds()) {
1040       Atom other = b.get1_2(a);
1041       if (other.getAtomicNumber() == atomicNumber) {
1042         return true;
1043       }
1044     }
1045     return false;
1046   }
1047 
1048   /**
1049    * getAtomListFromPool
1050    *
1051    * @return a {@link java.util.List} object.
1052    */
1053   private static List<Atom> getAtomListFromPool() {
1054     if (atomListPool.isEmpty()) {
1055       return new ArrayList<>();
1056     }
1057     return atomListPool.remove(0);
1058   }
1059 
1060   /**
1061    * Convert possibly duplicate chainID into a unique segID.
1062    *
1063    * @param c chain ID just read.
1064    * @return a unique segID.
1065    */
1066   private static String getSegID(Character c, List<String> segIDs) {
1067     if (c == null || c.equals(' ')) {
1068       c = 'A';
1069     }
1070 
1071     // Loop through existing segIDs to find the first one that is unused.
1072     int m = 0;
1073     for (String segID : segIDs) {
1074       if (segID.endsWith(c.toString())) {
1075         m++;
1076       }
1077     }
1078 
1079     // If the count is greater than 0, then append it.
1080     String newSegID;
1081     if (m == 0) {
1082       newSegID = c.toString();
1083     } else {
1084       newSegID = c.toString() + m;
1085     }
1086 
1087     segIDs.add(newSegID);
1088     return newSegID;
1089   }
1090 
1091   /**
1092    * Returns true if the lists contain any atom in common.
1093    *
1094    * @param list1 List
1095    * @param list2 List
1096    * @return boolean
1097    */
1098   private static boolean haveCommonAtom(List<Atom> list1, List<Atom> list2) {
1099     if (list1 == null || list2 == null) {
1100       return false;
1101     }
1102     for (Atom a : list1) {
1103       if (list2.contains(a)) {
1104         return true;
1105       }
1106     }
1107     return false;
1108   }
1109 
1110   /**
1111    * Returns true if Atom a is water oxygen.
1112    *
1113    * @param a Atom
1114    * @return boolean
1115    */
1116   static boolean isWaterOxygen(Atom a) {
1117     if (a.getAtomicNumber() != 8) {
1118       return false;
1119     }
1120     for (Bond b : a.getBonds()) {
1121       Atom h = b.get1_2(a);
1122       if (h.getAtomicNumber() != 1) {
1123         return false;
1124       }
1125     }
1126     return true;
1127   }
1128 
1129   /**
1130    * This function returns the number of times atom "a" is bonded to an atom of the specified atomic
1131    * number.
1132    *
1133    * @param a Atom
1134    * @param atomicNumber int
1135    * @return int
1136    */
1137   private static int numberOfBondsWith(Atom a, int atomicNumber) {
1138     int total = 0;
1139     for (Bond b : a.getBonds()) {
1140       Atom other = b.get1_2(a);
1141       if (other.getAtomicNumber() == atomicNumber) {
1142         total++;
1143       }
1144     }
1145     return total;
1146   }
1147 
1148   /**
1149    * Check to see if a portion of the backbone matches that of a biological polymer, and if so
1150    * determine the respective residue
1151    *
1152    * @param start First atom.
1153    * @param backbone List of backbone atoms.
1154    * @param type Type of residue.
1155    * @return The residue.
1156    */
1157   private static Residue patternMatch(int start, List<Atom> backbone, PolymerType type) {
1158     int[] pattern;
1159 
1160     // Initialization
1161     if (type == PolymerType.AMINOACID) {
1162       pattern = AminoAcidUtils.AAPATTERN;
1163       if (backbone.size() < start + pattern.length) {
1164         return null;
1165       }
1166       // Check for correct Carbonyl placement
1167       Atom a = backbone.get(start + 1);
1168       if (formsBondsWith(a, 8)) {
1169         return null;
1170       }
1171       a = backbone.get(start + 2);
1172       if (!formsBondsWith(a, 8)) {
1173         return null;
1174       }
1175     } else if (type == PolymerType.NUCLEICACID) {
1176       pattern = NucleicAcidUtils.NAPATTERN;
1177       if (backbone.size() < start + pattern.length) {
1178         return null;
1179       }
1180     } else {
1181       return null;
1182     }
1183 
1184     int length = pattern.length;
1185     List<Atom> atoms = getAtomListFromPool();
1186     List<Atom> sidePolymer = getAtomListFromPool();
1187     for (int i = 0; i < length; i++) {
1188       Atom a = backbone.get(start + i);
1189       // Add backbone atoms to terminate sidePolymer search
1190       sidePolymer.add(a);
1191       if (a.getAtomicNumber() != pattern[i]) {
1192         return null;
1193       }
1194     }
1195     // Collect all the atoms in the Residue
1196     // Add the atom before and after the backbone pattern to
1197     // terminate the search, then remove them
1198     if (start > 0) {
1199       atoms.add(backbone.get(start - 1));
1200     }
1201     if (start + length < backbone.size()) {
1202       atoms.add(backbone.get(start + length));
1203     }
1204     collectAtoms(backbone.get(start), atoms);
1205     if (start > 0) {
1206       atoms.remove(0);
1207     }
1208     if (start + length < backbone.size()) {
1209       atoms.remove(0);
1210       // Collect Just Side chain atoms, then remove backbone termination
1211     }
1212 
1213     if (type == PolymerType.AMINOACID) {
1214       collectAtoms(sidePolymer.get(1), sidePolymer);
1215     } else {
1216       Atom seed = null;
1217       for (Atom a : atoms) {
1218         if (a.getAtomicNumber() == 7) {
1219           seed = a;
1220           break;
1221         }
1222       }
1223       if (seed != null && seed.getAtomicNumber() == 7) {
1224         sidePolymer.add(seed);
1225         collectAtoms(seed, sidePolymer);
1226       } else {
1227         return null;
1228       }
1229     }
1230 
1231     sidePolymer.subList(0, length + 1).clear();
1232 
1233     return assignResidue(backbone, start, atoms, sidePolymer);
1234   }
1235 
1236   /**
1237    * Returns an List with reversed ordering.
1238    *
1239    * @param atomList List
1240    * @return List
1241    */
1242   private static List<Atom> reverseAtomList(List<Atom> atomList) {
1243     List<Atom> reversedList = getAtomListFromPool();
1244     for (Atom a : atomList) {
1245       reversedList.add(0, a);
1246     }
1247     return reversedList;
1248   }
1249 
1250   /** An enumeration of recognized file types. */
1251   public enum FileType {
1252     XYZ,
1253     XPH,
1254     INT,
1255     ARC,
1256     PDB,
1257     CIF,
1258     ANY,
1259     SIM,
1260     UNK
1261   }
1262 
1263   /** An enumeration of recognized organic polymers. */
1264   public enum PolymerType {
1265     AMINOACID,
1266     NUCLEICACID,
1267     UNKNOWN
1268   }
1269 }