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.bonded;
39  
40  import ffx.numerics.math.DoubleMath;
41  import ffx.potential.bonded.AminoAcidUtils.ResiduePosition;
42  import ffx.potential.bonded.BondedUtils.MissingAtomTypeException;
43  import ffx.potential.bonded.BondedUtils.MissingHeavyAtomException;
44  import ffx.potential.parameters.AtomType;
45  import ffx.potential.parameters.ForceField;
46  
47  import java.util.Arrays;
48  import java.util.HashMap;
49  import java.util.List;
50  import java.util.Map;
51  import java.util.logging.Level;
52  import java.util.logging.Logger;
53  
54  import static ffx.numerics.math.DoubleMath.dist;
55  import static ffx.potential.bonded.AminoAcidUtils.ResiduePosition.FIRST_RESIDUE;
56  import static ffx.potential.bonded.AminoAcidUtils.ResiduePosition.LAST_RESIDUE;
57  import static ffx.potential.bonded.AminoAcidUtils.ResiduePosition.MIDDLE_RESIDUE;
58  import static ffx.potential.bonded.BondedUtils.buildBond;
59  import static ffx.potential.bonded.BondedUtils.buildH;
60  import static ffx.potential.bonded.BondedUtils.buildHeavy;
61  import static ffx.potential.bonded.BondedUtils.intxyz;
62  import static java.lang.String.format;
63  import static org.apache.commons.math3.util.FastMath.toDegrees;
64  
65  /**
66   * Utilities for creating Nucleic Acid residues.
67   *
68   * @author Michael Schnieders
69   * @since 1.0
70   */
71  public class NucleicAcidUtils {
72  
73    private static final Logger logger = Logger.getLogger(NucleicAcidUtils.class.getName());
74  
75    /**
76     * The 4 RNA bases, 4 DNA bases, and mono- di- and triphosphate.
77     */
78    public enum NA {
79      ADENINE, CYTOSINE, GUANINE, URACIL, DEOXYADENINE, DEOXYCYTOSINE, DEOXYGUANINE, THYMINE, MONOPHOSPHATE, DIPHOSPHATE, TRIPHOSPHATE
80    }
81  
82    public enum NucleicAcid1 {
83      A, G, C, U, D, B, I, T, O, W, H, X
84    }
85  
86    /**
87     * Since enumeration values must start with a letter, an 'M' is added to modified bases whose IUPAC
88     * name starts with an integer.
89     */
90    public enum NucleicAcid3 {
91      ADE, GUA, CYT, URI, DAD, DGU, DCY, DTY, THY, MP1, DP2, TP3, UNK, M2MG, H2U, M2G, OMC, OMG, PSU, M5MC, M7MG, M5MU, M1MA, YYG;
92  
93      /**
94       * Best-guess parse of a String to an NA3.
95       *
96       * @param name Parse to NA3.
97       * @return Corresponding NA3.
98       * @throws IllegalArgumentException For 'DU', which has no implemented NA3.
99       */
100     public static NucleicAcid3 parse(String name) throws IllegalArgumentException {
101       // Only semi-abnormal cases: THY parses to DT instead of T, and DU throws an exception.
102       return switch (name.toUpperCase()) {
103         case "ADE", "A" -> NucleicAcid3.ADE;
104         case "CYT", "C" -> NucleicAcid3.CYT;
105         case "GUA", "G" -> NucleicAcid3.GUA;
106         case "URI", "U" -> NucleicAcid3.URI;
107         case "DAD", "DA" -> NucleicAcid3.DAD;
108         case "DCY", "DC" -> NucleicAcid3.DCY;
109         case "DGU", "DG" -> NucleicAcid3.DGU;
110         case "DTY", "THY", "DT" -> NucleicAcid3.DTY;
111         case "MPO" -> NucleicAcid3.MP1;
112         case "DPO" -> NucleicAcid3.DP2;
113         case "TPO" -> NucleicAcid3.TP3;
114         case "DU" -> throw new IllegalArgumentException(" No NA3 value exists for deoxy-uracil!");
115         default -> NucleicAcid3.UNK;
116       };
117     }
118   }
119 
120   /**
121    * Biotype keys for nucleic acid backbone atom types. These are consistent with parameter files
122    * from TINKER v. 6.1 (June 2012).
123    */
124   public static final int[] NA_O5 = {1001, 1031, 1062, 1090, 1117, 1146, 1176, 1203, 0, 0, 0, 0,
125       1300, 1334, 1363, 1400, 1431, 1465, 1492, 1523, 1559, 1589, 1624};
126   /** Constant <code>NA_C5</code> */
127   public static final int[] NA_C5 = {1002, 1032, 1063, 1091, 1118, 1147, 1177, 1204, 0, 0, 0, 0,
128       1301, 1335, 1364, 1401, 1432, 1466, 1493, 1524, 1560, 1590, 1625};
129   /** Constant <code>NA_H51</code> */
130   public static final int[] NA_H51 = {1003, 1033, 1064, 1092, 1119, 1148, 1178, 1205, 0, 0, 0, 0,
131       1302, 1336, 1365, 1402, 1433, 1467, 1494, 1525, 1561, 1591, 1626};
132   /** Constant <code>NA_H52</code> */
133   public static final int[] NA_H52 = {1004, 1034, 1065, 1093, 1120, 1149, 1179, 1206, 0, 0, 0, 0,
134       1303, 1337, 1366, 1403, 1434, 1468, 1495, 1526, 1562, 1592, 1627};
135   /** Constant <code>NA_C4</code> */
136   public static final int[] NA_C4 = {1005, 1035, 1066, 1094, 1121, 1150, 1180, 1207, 0, 0, 0, 0,
137       1304, 1338, 1367, 1404, 1435, 1469, 1496, 1527, 1563, 1593, 1628};
138   /** Constant <code>NA_H4</code> */
139   public static final int[] NA_H4 = {1006, 1036, 1067, 1095, 1122, 1151, 1181, 1208, 0, 0, 0, 0,
140       1305, 1339, 1368, 1405, 1436, 1470, 1497, 1528, 1564, 1594, 1629};
141   /** Constant <code>NA_O4</code> */
142   public static final int[] NA_O4 = {1007, 1037, 1068, 1096, 1123, 1152, 1182, 1209, 0, 0, 0, 0,
143       1306, 1340, 1369, 1406, 1437, 1471, 1498, 1529, 1565, 1595, 1630};
144   /** Constant <code>NA_C1</code> */
145   public static final int[] NA_C1 = {1008, 1038, 1069, 1097, 1124, 1153, 1183, 1210, 0, 0, 0, 0,
146       1307, 1341, 1370, 1407, 1438, 1472, 1499, 1530, 1566, 1596, 1631};
147   /** Constant <code>NA_H1</code> */
148   public static final int[] NA_H1 = {1009, 1039, 1070, 1098, 1125, 1154, 1184, 1211, 0, 0, 0, 0,
149       1308, 1342, 1371, 1408, 1439, 1473, 1500, 1531, 1567, 1597, 1632};
150   /** Constant <code>NA_C3</code> */
151   public static final int[] NA_C3 = {1010, 1040, 1071, 1099, 1126, 1155, 1185, 1212, 0, 0, 0, 0,
152       1309, 1343, 1372, 1409, 1440, 1474, 1501, 1532, 1568, 1598, 1633};
153   /** Constant <code>NA_H3</code> */
154   public static final int[] NA_H3 = {1011, 1041, 1072, 1100, 1127, 1156, 1186, 1213, 0, 0, 0, 0,
155       1310, 1344, 1373, 1410, 1441, 1475, 1502, 1533, 1569, 1599, 1634};
156   /** Constant <code>NA_C2</code> */
157   public static final int[] NA_C2 = {1012, 1042, 1073, 1101, 1128, 1157, 1187, 1214, 0, 0, 0, 0,
158       1311, 1345, 1374, 1411, 1442, 1476, 1503, 1534, 1570, 1600, 1635};
159   /** Constant <code>NA_H21</code> */
160   public static final int[] NA_H21 = {1013, 1043, 1074, 1102, 1129, 1158, 1188, 1215, 0, 0, 0, 0,
161       1312, 1346, 1375, 1412, 1443, 1477, 1504, 1535, 1571, 1601, 1636};
162   /** Constant <code>NA_O2</code> */
163   public static final int[] NA_O2 = {1014, 1044, 1075, 1103, 0, 0, 0, 0, 0, 0, 0, 0, 1313, 1347,
164       1376, 1413, 1444, 1478, 1505, 1536, 1572, 1602, 1637};
165   /** Constant <code>NA_H22</code> */
166   public static final int[] NA_H22 = {1015, 1045, 1076, 1104, 1130, 1159, 1189, 1216, 0, 0, 0, 0,
167       1314, 1348, 1377, 0, 0, 1479, 1506, 1537, 1573, 1603, 1638};
168   /** Constant <code>NA_O3</code> */
169   public static final int[] NA_O3 = {1016, 1046, 1077, 1105, 1131, 1160, 1190, 1217, 0, 0, 0, 0,
170       1315, 1349, 1378, 1414, 1445, 1480, 1507, 1538, 1574, 1604, 1639};
171   /** Constant <code>NA_P</code> */
172   public static final int[] NA_P = {1230, 1230, 1230, 1230, 1242, 1242, 1242, 1242, 0, 0, 0, 0, 1230,
173       1230, 1230, 1230, 1230, 1230, 1230, 1230, 1230, 1230, 1230};
174   /** Constant <code>NA_OP</code> */
175   public static final int[] NA_OP = {1231, 1231, 1231, 1231, 1243, 1243, 1243, 1243, 0, 0, 0, 0,
176       1231, 1231, 1231, 1231, 1231, 1231, 1231, 1231, 1231, 1231, 1231};
177   /**
178    * Should be NA_HO5' (' replaced by T in the name).
179    *
180    * <p>Constant <code>NA_HO5T</code>
181    */
182   public static final int[] NA_HO5T = {1233, 1233, 1233, 1233, 1245, 1245, 1245, 1245, 0, 0, 0, 0,
183       1233, 1233, 1233, 1233, 1233, 1233, 1233, 1233, 1233, 1233, 1233};
184   /**
185    * Should be NA_HO3' (' replaced by T in the name).
186    *
187    * <p>Constant <code>NA_HO3T</code>
188    */
189   public static final int[] NA_HO3T = {1238, 1238, 1238, 1238, 1250, 1250, 1250, 1250, 0, 0, 0, 0,
190       1238, 1238, 1238, 1238, 1238, 1238, 1238, 1238, 1238, 1238, 1238};
191 
192   /** Constant <code>nucleicAcidList</code> */
193   public static final List<NucleicAcid3> nucleicAcidList = Arrays.asList(NucleicAcid3.values());
194 
195   /** Constant <code>NA1toNA3</code> */
196   public static final HashMap<NucleicAcid1, NucleicAcid3> NA1toNA3 = new HashMap<>();
197 
198   /** Substitutes from old PDB nucleic acid atom names to new ones. */
199   private static final Map<String, String> newNucleicAcidAtomNames;
200 
201   /** Repeating atomic numbers of a nucleic acid chain. */
202   public static final int[] NAPATTERN = {8, 6, 6, 6, 8};
203 
204   static {
205     newNucleicAcidAtomNames = Map.of("HO'", "HO2'", "H3T", "HO3'", "H5'1", "H5'", "H5'2", "H5''",
206         "H5T", "HO5'", "H2'1", "H2'", "H2'2", "H2''");
207 
208     NucleicAcid1[] na1 = NucleicAcid1.values();
209     NucleicAcid3[] na3 = NucleicAcid3.values();
210     for (int i = 0; i < NucleicAcid1.values().length; i++) {
211       NA1toNA3.put(na1[i], na3[i]);
212     }
213   }
214 
215   /**
216    * Assign atom types for a nucleic acid polymer.
217    *
218    * @param residues A list of residues that form the nucleic acid polymer.
219    * @param forceField The ForceField in use.
220    * @param bondList A list of created bonds.
221    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
222    * @throws ffx.potential.bonded.BondedUtils.MissingAtomTypeException if any.
223    */
224   public static void assignNucleicAcidAtomTypes(List<Residue> residues, ForceField forceField,
225       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
226 
227     // A reference to the O3* atom of the previous base.
228     Atom pSugarO3 = null;
229 
230     // Loop over residues.
231     int numberOfResidues = residues.size();
232     for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
233 
234       // Match the residue name to a known nucleic acid residue.
235       Residue residue = residues.get(residueNumber);
236       String residueName = residue.getName().toUpperCase();
237       NucleicAcid3 nucleicAcid = null;
238       int naNumber = -1;
239       for (NucleicAcid3 nucleic : nucleicAcidList) {
240         naNumber++;
241         String nuc3 = nucleic.toString();
242         nuc3 = nuc3.substring(nuc3.length() - 3);
243         if (nuc3.equalsIgnoreCase(residueName)) {
244           nucleicAcid = nucleic;
245           break;
246         }
247       }
248 
249       // Do atom name conversions.
250       List<Atom> resAtoms = residue.getAtomList();
251       int nAtoms = resAtoms.size();
252       for (Atom atom : resAtoms) {
253         String name = atom.getName();
254         name = name.replace('*', '\'');
255         name = newNucleicAcidAtomNames.getOrDefault(name, name);
256         atom.setName(name);
257       }
258 
259       // Set a position flag.
260       ResiduePosition position = MIDDLE_RESIDUE;
261       boolean lastRes = false;
262       if (residueNumber == 0) {
263         position = FIRST_RESIDUE;
264       }
265 
266       if (residueNumber == numberOfResidues - 1) {
267         lastRes = true;
268         if (position != FIRST_RESIDUE) {
269           position = LAST_RESIDUE;
270         }
271       }
272 
273       // Check if this is a lone 3'-terminal phosphate cap; if so, will skip a lot of logic.
274       boolean threePrimeCap = isThreePrimeCap(nAtoms, resAtoms);
275 
276       // Check if the sugar is deoxyribose and change the residue name if necessary.
277       boolean isDNA = false;
278       Residue resToCheck = threePrimeCap ? residues.get(residueNumber - 1) : residue;
279       Atom sugarO2 = (Atom) resToCheck.getAtomNode("O2'");
280 
281       if (sugarO2 == null) {
282         // Assume deoxyribose (DNA) since there is an O2* atom.
283         isDNA = true;
284         if (!residueName.startsWith("D")) {
285           switch (nucleicAcid) {
286             case ADE -> {
287               nucleicAcid = NucleicAcid3.DAD;
288               residueName = "DAD";
289               residue.setName(residueName);
290             }
291             case CYT -> {
292               nucleicAcid = NucleicAcid3.DCY;
293               residueName = "DCY";
294               residue.setName(residueName);
295             }
296             case GUA -> {
297               nucleicAcid = NucleicAcid3.DGU;
298               residueName = "DGU";
299               residue.setName(residueName);
300             }
301             case THY -> {
302               nucleicAcid = NucleicAcid3.DTY;
303               residueName = "DTY";
304               residue.setName(residueName);
305             }
306             default -> {
307             }
308           }
309         }
310       } else {
311         // Assume ribose (RNA) since there is no O2* atom.
312         if (residueName.startsWith("D")) {
313           switch (nucleicAcid) {
314             case DAD -> {
315               nucleicAcid = NucleicAcid3.ADE;
316               residueName = "ADE";
317               residue.setName(residueName);
318             }
319             case DCY -> {
320               nucleicAcid = NucleicAcid3.CYT;
321               residueName = "CYT";
322               residue.setName(residueName);
323             }
324             case DGU -> {
325               nucleicAcid = NucleicAcid3.GUA;
326               residueName = "GUA";
327               residue.setName(residueName);
328             }
329             case DTY -> {
330               nucleicAcid = NucleicAcid3.THY;
331               residueName = "THY";
332               residue.setName(residueName);
333             }
334             default -> {
335             }
336           }
337         }
338       }
339 
340       // Build the phosphate atoms of the current residue.
341       Atom phosphate = null;
342       Atom sugarO5 = null;
343       Atom sugarO3 = null;
344 
345       if (threePrimeCap) {
346         if (logger.isLoggable(Level.FINE)) {
347           logger.fine(format(" EXPERIMENTAL: Adding 3'-terminal phosphate cap %s", residue));
348         }
349         Residue priorResidue = residues.get(residueNumber - 1);
350         int phosType = isDNA ? 1252 : 1240;
351         int oxygenType = isDNA ? 1253 : 1241;
352 
353         phosphate = buildHeavy(residue, "P", pSugarO3, phosType, forceField, bondList);
354         buildHeavy(residue, "OP1", phosphate, oxygenType, forceField, bondList);
355         buildHeavy(residue, "OP2", phosphate, oxygenType, forceField, bondList);
356         switch (nAtoms) {
357           case 3 ->
358               buildOP3(residue, phosphate, oxygenType, forceField, bondList, priorResidue, false);
359           case 4 -> {
360             MSNode bogusO5s = residue.getAtomNode("O5'");
361             if (bogusO5s != null) {
362               bogusO5s.setName("OP3");
363             }
364             buildHeavy(residue, "OP3", phosphate, oxygenType, forceField, bondList);
365           }
366           case 5, 6 -> logger.severe(" Currently, FFX does not support partially-protonated "
367               + "3'-terminal phosphate caps from PDB files!");
368         }
369       } else {
370         if (position == FIRST_RESIDUE) {
371           /*
372            * The 5' O5' oxygen of the nucleic acid is generally terminated by
373            * 1.) A phosphate group PO3 (-3).
374            * 2.) A hydrogen.
375            * If the base has phosphate atom we will assume a PO3 group.
376            */
377           phosphate = (Atom) residue.getAtomNode("P");
378           if (phosphate != null) {
379             Residue nextRes = lastRes ? null : residues.get(residueNumber + 1);
380             if (isDNA) {
381               phosphate = buildHeavy(residue, "P", null, 1247, forceField, bondList);
382               buildHeavy(residue, "OP1", phosphate, 1248, forceField, bondList);
383               buildHeavy(residue, "OP2", phosphate, 1248, forceField, bondList);
384               buildOP3(residue, phosphate, 1248, forceField, bondList, nextRes, true);
385               sugarO5 = buildHeavy(residue, "O5'", phosphate, 1246, forceField, bondList);
386             } else {
387               phosphate = buildHeavy(residue, "P", null, 1235, forceField, bondList);
388               buildHeavy(residue, "OP1", phosphate, 1236, forceField, bondList);
389               buildHeavy(residue, "OP2", phosphate, 1236, forceField, bondList);
390               buildOP3(residue, phosphate, 1236, forceField, bondList, nextRes, true);
391               sugarO5 = buildHeavy(residue, "O5'", phosphate, 1234, forceField, bondList);
392             }
393           } else if (isDNA) {
394             Atom O5 = residue.getAtomByName("O5'", true);
395             if (O5 == null) {
396               sugarO5 = buildSugarO5(residue, 1244, forceField);
397             } else {
398               sugarO5 = buildHeavy(residue, "O5'", null, 1244, forceField, bondList);
399             }
400           } else {
401             Atom O5 = residue.getAtomByName("O5'", true);
402             if (O5 == null) {
403               sugarO5 = buildSugarO5(residue, 1232, forceField);
404             } else {
405               sugarO5 = buildHeavy(residue, "O5'", null, 1232, forceField, bondList);
406             }
407           }
408         } else {
409           phosphate = buildHeavy(residue, "P", pSugarO3, NA_P[naNumber], forceField, bondList);
410           buildHeavy(residue, "OP1", phosphate, NA_OP[naNumber], forceField, bondList);
411           buildHeavy(residue, "OP2", phosphate, NA_OP[naNumber], forceField, bondList);
412           sugarO5 = buildHeavy(residue, "O5'", phosphate, NA_O5[naNumber], forceField, bondList);
413         }
414 
415         // Build the ribose sugar atoms of the current base.
416         Atom sugarC5 = buildHeavy(residue, "C5'", sugarO5, NA_C5[naNumber], forceField, bondList);
417         Atom sugarC4 = buildHeavy(residue, "C4'", sugarC5, NA_C4[naNumber], forceField, bondList);
418         Atom sugarO4 = buildHeavy(residue, "O4'", sugarC4, NA_O4[naNumber], forceField, bondList);
419         Atom sugarC1 = buildHeavy(residue, "C1'", sugarO4, NA_C1[naNumber], forceField, bondList);
420         Atom sugarC3 = buildHeavy(residue, "C3'", sugarC4, NA_C3[naNumber], forceField, bondList);
421         Atom sugarC2 = buildHeavy(residue, "C2'", sugarC3, NA_C2[naNumber], forceField, bondList);
422         buildBond(sugarC2, sugarC1, forceField, bondList);
423         if (position == LAST_RESIDUE || numberOfResidues == 1) {
424           if (isDNA) {
425             sugarO3 = buildHeavy(residue, "O3'", sugarC3, 1249, forceField, bondList);
426           } else {
427             sugarO3 = buildHeavy(residue, "O3'", sugarC3, 1237, forceField, bondList);
428           }
429         } else {
430           boolean nextResIsCap = (residues.get(residueNumber + 1).getAtomList().size() < 7);
431           int o3Type = NA_O3[naNumber];
432           if (nextResIsCap) {
433             logger.fine(" Applying a 3'-terminal-phos-cap O3' atom type to residue " + residue);
434             o3Type = isDNA ? 1251 : 1239;
435           }
436           sugarO3 = buildHeavy(residue, "O3'", sugarC3, o3Type, forceField, bondList);
437         }
438         if (!isDNA) {
439           sugarO2 = buildHeavy(residue, "O2'", sugarC2, NA_O2[naNumber], forceField, bondList);
440         }
441 
442         // Build the backbone hydrogen atoms.
443         if (position == FIRST_RESIDUE && phosphate == null) {
444           buildH(residue, "HO5'", sugarO5, 1.00e0, sugarC5, 109.5e0, sugarC4, 180.0e0, 0,
445               NA_HO5T[naNumber], forceField, bondList);
446         }
447         buildH(residue, "H5'", sugarC5, 1.09e0, sugarO5, 109.5e0, sugarC4, 109.5e0, 1,
448             NA_H51[naNumber], forceField, bondList);
449         buildH(residue, "H5''", sugarC5, 1.09e0, sugarO5, 109.5e0, sugarC4, 109.5e0, -1,
450             NA_H52[naNumber], forceField, bondList);
451         buildH(residue, "H4'", sugarC4, 1.09e0, sugarC5, 109.5e0, sugarC3, 109.5e0, -1,
452             NA_H4[naNumber], forceField, bondList);
453         buildH(residue, "H3'", sugarC3, 1.09e0, sugarC4, 109.5e0, sugarC2, 109.5e0, -1,
454             NA_H3[naNumber], forceField, bondList);
455         if (isDNA) {
456           buildH(residue, "H2'", sugarC2, 1.09e0, sugarC3, 109.5e0, sugarC1, 109.5e0, -1,
457               NA_H21[naNumber], forceField, bondList);
458           buildH(residue, "H2''", sugarC2, 1.09e0, sugarC3, 109.5e0, sugarC1, 109.5e0, 1,
459               NA_H22[naNumber], forceField, bondList);
460         } else {
461           buildH(residue, "H2'", sugarC2, 1.09e0, sugarC3, 109.5e0, sugarC1, 109.5e0, -1,
462               NA_H21[naNumber], forceField, bondList);
463           // Add the NA_O2' Methyl for OMC and OMG
464           // Note: may be completely out-of-date with the current Chemical Component Dictionary.
465           if (nucleicAcid == NucleicAcid3.OMC || nucleicAcid == NucleicAcid3.OMG) {
466             Atom CM2 = buildHeavy(residue, "CM2", sugarO2, 1427, forceField, bondList);
467             Atom HM21 = buildH(residue, "HM21", CM2, 1.08e0, sugarO2, 109.5e0, sugarC2, 0.0e0, 0,
468                 1428, forceField, bondList);
469             buildH(residue, "HM22", CM2, 1.08e0, sugarO2, 109.5e0, HM21, 109.5e0, 1, 1429,
470                 forceField, bondList);
471             buildH(residue, "HM23", CM2, 1.08e0, sugarO2, 109.5e0, HM21, 109.5e0, -1, 1430,
472                 forceField, bondList);
473           } else {
474             buildH(residue, "HO2'", sugarO2, 1.00e0, sugarC2, 109.5e0, sugarC3, 180.0e0, 0,
475                 NA_H22[naNumber], forceField, bondList);
476           }
477         }
478         buildH(residue, "H1'", sugarC1, 1.09e0, sugarO4, 109.5e0, sugarC2, 109.5e0, -1,
479             NA_H1[naNumber], forceField, bondList);
480         if (position == LAST_RESIDUE || numberOfResidues == 1) {
481           buildH(residue, "HO3'", sugarO3, 1.00e0, sugarC3, 109.5e0, sugarC4, 180.0e0, 0,
482               NA_HO3T[naNumber], forceField, bondList);
483           // Else, if it is terminated by a 3' phosphate cap:
484           // Will need to see how PDB would label a 3' phosphate cap.
485         }
486 
487         // Build the nucleic acid base.
488         assignNucleicAcidBaseAtomTypes(nucleicAcid, residue, sugarC1, sugarO4, sugarC2, forceField,
489             bondList);
490       }
491 
492       // Do some checks on the current base to make sure all atoms have been assigned an atom type.
493       resAtoms = residue.getAtomList();
494       for (Atom atom : resAtoms) {
495         AtomType atomType = atom.getAtomType();
496         if (atomType == null) {
497           throw new MissingAtomTypeException(residue, atom);
498         }
499         int numberOfBonds = atom.getNumBonds();
500         if (numberOfBonds != atomType.valence) {
501           if (atom == sugarO3 && numberOfBonds == atomType.valence - 1 && position != LAST_RESIDUE
502               && numberOfResidues != 1) {
503             continue;
504           }
505           logger.warning(
506               format(" An atom for residue %s has the wrong number of bonds:\n %s", residueName,
507                   atom));
508           logger.info(format(" Expected: %d Actual: %d.", atomType.valence, numberOfBonds));
509           for (Bond bond : atom.getBonds()) {
510             logger.info(" " + bond.toString());
511           }
512         }
513       }
514 
515       // Save a reference to the current O3* oxygen.
516       pSugarO3 = sugarO3;
517     }
518   }
519 
520   private static boolean isThreePrimeCap(int nAtoms, List<Atom> resAtoms) {
521     boolean threePrimeCap = false;
522     if (nAtoms < 7) {
523       boolean unrecognized = false;
524       for (Atom atom : resAtoms) {
525         String atomName = atom.getName();
526         // Recognized: P, OP[1-3], O5' (will be renamed to OP3), HOP[1-3], and DOP[1-3])
527         if (!(atomName.equals("P") || atomName.matches("[HD]?OP[1-3]") || atomName.matches(
528             "O5'"))) {
529           unrecognized = true;
530           break;
531         }
532       }
533       threePrimeCap = !unrecognized;
534     }
535     return threePrimeCap;
536   }
537 
538   /**
539    * Assign atom types to the nucleic acid base.
540    *
541    * @param nucleicAcid The nucleic acid base to use.
542    * @param residue The residue node.
543    * @param C1s The C1* attachment atom.
544    * @param O4s The O4* attachment atom.
545    * @param C2s The C2* attachment atom.
546    * @param forceField The ForceField in use.
547    * @param bondList List of created bonds.
548    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
549    * @since 1.0
550    */
551   private static void assignNucleicAcidBaseAtomTypes(NucleicAcid3 nucleicAcid, Residue residue,
552       Atom C1s, Atom O4s, Atom C2s, ForceField forceField, List<Bond> bondList)
553       throws MissingHeavyAtomException, MissingAtomTypeException {
554     double glyco = 0.0; // DEFAULT - should it be more like 240?
555 
556     // glycosyl torsion: pyrimidines:O4′-C1′-N1-C2 | purines:O4′-C1′-N9-C4
557     Atom o4p = residue.getAtomByName("O4'", true);
558     Atom c1p = residue.getAtomByName("C1'", true);
559     Atom baseN = null;
560     Atom baseC = null;
561     if (nucleicAcid == NucleicAcid3.DAD || nucleicAcid == NucleicAcid3.DGU) {
562       baseN = residue.getAtomByName("N9", true);
563       baseC = residue.getAtomByName("C4", true);
564     } else if (nucleicAcid == NucleicAcid3.DCY || nucleicAcid == NucleicAcid3.DTY) {
565       baseN = residue.getAtomByName("N1", true);
566       baseC = residue.getAtomByName("C2", true);
567     }
568 
569     if (o4p != null && c1p != null && baseN != null && baseC != null) {
570       double angle = DoubleMath.dihedralAngle(o4p.getXYZ(null), c1p.getXYZ(null), baseN.getXYZ(null), baseC.getXYZ(null));
571       glyco = toDegrees(angle);
572     }
573 
574     switch (nucleicAcid) {
575       case ADE -> buildADE(residue, C1s, O4s, C2s, glyco, forceField, bondList);
576       case M1MA -> buildM1MA(residue, C1s, forceField, bondList);
577       case CYT -> buildCYT(residue, C1s, O4s, C2s, glyco, forceField, bondList);
578       case OMC -> buildOMC(residue, C1s, O4s, C2s, glyco, forceField, bondList);
579       case M5MC -> buildM5MC(residue, C1s, forceField, bondList);
580       case GUA -> buildGUA(residue, C1s, O4s, C2s, glyco, forceField, bondList);
581       case OMG -> buildOMG(residue, C1s, O4s, C2s, glyco, forceField, bondList);
582       case YYG -> buildYYG(residue, C1s, forceField, bondList);
583       case M2MG -> buildM2MG(residue, C1s, forceField, bondList);
584       case M2G -> buildM2G(residue, C1s, forceField, bondList);
585       case M7MG -> buildM7MG(residue, C1s, forceField, bondList);
586       case URI -> buildURI(residue, C1s, O4s, C2s, glyco, forceField, bondList);
587       case PSU -> buildPSU(residue, C1s, forceField, bondList);
588       case H2U -> buildH2U(residue, C1s, forceField, bondList);
589       case M5MU -> buildM5MU(residue, C1s, forceField, bondList);
590       case DAD -> buildDAD(residue, C1s, O4s, C2s, glyco, forceField, bondList);
591       case DCY -> buildDCY(residue, C1s, O4s, C2s, glyco, forceField, bondList);
592       case DGU -> buildDGU(residue, C1s, O4s, C2s, glyco, forceField, bondList);
593       case DTY -> buildDTY(residue, C1s, O4s, C2s, glyco, forceField, bondList);
594     }
595   }
596 
597   /**
598    * buildADE.
599    *
600    * @param residue a {@link ffx.potential.bonded.Residue} object.
601    * @param C1s a {@link ffx.potential.bonded.Atom} object.
602    * @param O4s a {@link ffx.potential.bonded.Atom} object.
603    * @param C2s a {@link ffx.potential.bonded.Atom} object.
604    * @param glyco a double.
605    * @param forceField a {@link ForceField} object.
606    * @param bondList a {@link java.util.List} object.
607    * @return a {@link ffx.potential.bonded.Residue} object.
608    */
609   private static Residue buildADE(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
610       ForceField forceField, List<Bond> bondList) {
611     Atom N9 = buildHeavy(residue, "N9", C1s, 1.48, O4s, 108.1, C2s, 113.7, 1, 1017, forceField,
612         bondList);
613     Atom C8 = buildHeavy(residue, "C8", N9, 1.37, C1s, 128.4, O4s, glyco + 180.0, 0, 1021,
614         forceField, bondList);
615     Atom N7 = buildHeavy(residue, "N7", C8, 1.30, N9, 113.8, C1s, 180.0, 0, 1020, forceField,
616         bondList);
617     Atom C5 = buildHeavy(residue, "C5", N7, 1.39, C8, 104.0, N9, 0.0, 0, 1019, forceField, bondList);
618     Atom C6 = buildHeavy(residue, "C6", C5, 1.40, N7, 132.4, C8, 180.0, 0, 1025, forceField,
619         bondList);
620     Atom N6 = buildHeavy(residue, "N6", C6, 1.34, C5, 123.5, N7, 0.0, 0, 1027, forceField, bondList);
621     Atom N1 = buildHeavy(residue, "N1", C6, 1.35, C5, 117.4, N7, 180.0, 0, 1024, forceField,
622         bondList);
623     Atom C2 = buildHeavy(residue, "C2", N1, 1.33, C6, 118.8, C5, 0.0, 0, 1023, forceField, bondList);
624     Atom N3 = buildHeavy(residue, "N3", C2, 1.32, N1, 129.2, C6, 0.0, 0, 1022, forceField, bondList);
625     Atom C4 = buildHeavy(residue, "C4", N3, 1.35, C2, 110.9, N1, 0.0, 0, 1018, forceField, bondList);
626     buildBond(C4, C5, forceField, bondList);
627     buildBond(C4, N9, forceField, bondList);
628     buildH(residue, "H8", C8, 1.08e0, N7, 123.1e0, C5, 180.0e0, 0, 1030, forceField, bondList);
629     buildH(residue, "H61", N6, 1.00e0, C6, 120.0e0, N7, 180.0e0, 0, 1028, forceField, bondList);
630     buildH(residue, "H62", N6, 1.00e0, C6, 120.0e0, N7, 0.0e0, 0, 1029, forceField, bondList);
631     buildH(residue, "H2", C2, 1.08e0, N3, 115.4e0, C4, 180.0e0, 0, 1026, forceField, bondList);
632     return residue;
633   }
634 
635   /**
636    * buildCYT.
637    *
638    * @param residue a {@link ffx.potential.bonded.Residue} object.
639    * @param C1s a {@link ffx.potential.bonded.Atom} object.
640    * @param O4s a {@link ffx.potential.bonded.Atom} object.
641    * @param C2s a {@link ffx.potential.bonded.Atom} object.
642    * @param glyco a double.
643    * @param forceField a {@link ForceField} object.
644    * @param bondList a {@link java.util.List} object.
645    * @return a {@link ffx.potential.bonded.Residue} object.
646    */
647   private static Residue buildCYT(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
648       ForceField forceField, List<Bond> bondList) {
649     Atom N1 = buildHeavy(residue, "N1", C1s, 1.48, O4s, 108.1, C2s, 113.7, 1, 1078, forceField,
650         bondList);
651     Atom C2 = buildHeavy(residue, "C2", N1, 1.37, C1s, 117.8, O4s, glyco + 180, 0, 1079, forceField,
652         bondList);
653     Atom O2 = buildHeavy(residue, "O2", C2, 1.24, N1, 118.9, C1s, 0.0, 0, 1084, forceField,
654         bondList);
655     Atom N3 = buildHeavy(residue, "N3", C2, 1.38, N1, 118.7, C1s, 180.0, 0, 1080, forceField,
656         bondList);
657     Atom C4 = buildHeavy(residue, "C4", N3, 1.34, C2, 120.6, N1, 0.0, 0, 1081, forceField, bondList);
658     Atom N4 = buildHeavy(residue, "N4", C4, 1.32, N3, 118.3, O2, 180.0, 0, 1085, forceField,
659         bondList);
660     Atom C5 = buildHeavy(residue, "C5", C4, 1.43, N3, 121.6, C2, 0.0, 0, 1082, forceField, bondList);
661     Atom C6 = buildHeavy(residue, "C6", C5, 1.36, C4, 116.9, N3, 0.0, 0, 1083, forceField, bondList);
662     buildBond(C6, N1, forceField, bondList);
663     buildH(residue, "H41", N4, 1.00e0, C4, 120.0e0, N3, 0.0e0, 0, 1086, forceField, bondList);
664     buildH(residue, "H42", N4, 1.00e0, C4, 120.0e0, N3, 180.0e0, 0, 1087, forceField, bondList);
665     buildH(residue, "H5", C5, 1.08e0, C4, 121.6e0, N3, 180.0e0, 0, 1088, forceField, bondList);
666     buildH(residue, "H6", C6, 1.08e0, C5, 119.4e0, C4, 180.0e0, 0, 1089, forceField, bondList);
667     return residue;
668   }
669 
670   /**
671    * buildDAD.
672    *
673    * @param residue a {@link ffx.potential.bonded.Residue} object.
674    * @param C1s a {@link ffx.potential.bonded.Atom} object.
675    * @param O4s a {@link ffx.potential.bonded.Atom} object.
676    * @param C2s a {@link ffx.potential.bonded.Atom} object.
677    * @param glyco a double.
678    * @param forceField a {@link ForceField} object.
679    * @param bondList a {@link java.util.List} object.
680    * @return a {@link ffx.potential.bonded.Residue} object.
681    */
682   private static Residue buildDAD(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
683       ForceField forceField, List<Bond> bondList) {
684 
685     Atom N9 = buildHeavy(residue, "N9", C1s, 1.48, O4s, 108.1, C2s, 113.7, 1, 1132, forceField,
686         bondList);
687     Atom C8 = buildHeavy(residue, "C8", N9, 1.37, C1s, 128.4, O4s, glyco + 180, 0, 1136, forceField,
688         bondList);
689     Atom N7 = buildHeavy(residue, "N7", C8, 1.30, N9, 113.8, C1s, 180.0, 0, 1135, forceField,
690         bondList);
691     Atom C5 = buildHeavy(residue, "C5", N7, 1.39, C8, 104.0, N9, 0.0, 0, 1134, forceField, bondList);
692     Atom C6 = buildHeavy(residue, "C6", C5, 1.40, N7, 132.4, C8, 180.0, 0, 1140, forceField,
693         bondList);
694     Atom N6 = buildHeavy(residue, "N6", C6, 1.34, C5, 123.5, N7, 0.0, 0, 1142, forceField, bondList);
695     Atom N1 = buildHeavy(residue, "N1", C6, 1.35, C5, 117.4, N7, 180.0, 0, 1139, forceField,
696         bondList);
697     Atom C2 = buildHeavy(residue, "C2", N1, 1.33, C6, 118.8, C5, 0.0, 0, 1138, forceField, bondList);
698     Atom N3 = buildHeavy(residue, "N3", C2, 1.32, N1, 129.2, C6, 0.0, 0, 1137, forceField, bondList);
699     Atom C4 = buildHeavy(residue, "C4", N3, 1.35, C2, 110.9, N1, 0.0, 0, 1133, forceField, bondList);
700     buildBond(C4, C5, forceField, bondList);
701     buildBond(C4, N9, forceField, bondList);
702     buildH(residue, "H8", C8, 1.08e0, N7, 123.1e0, C5, 180.0e0, 0, 1145, forceField, bondList);
703     buildH(residue, "H61", N6, 1.00e0, C6, 120.0e0, N7, 180.0e0, 0, 1143, forceField, bondList);
704     buildH(residue, "H62", N6, 1.00e0, C6, 120.0e0, N7, 0.0e0, 0, 1144, forceField, bondList);
705     buildH(residue, "H2", C2, 1.08e0, N3, 115.4e0, C4, 180.0e0, 0, 1141, forceField, bondList);
706     return residue;
707   }
708 
709   /**
710    * buildDCY.
711    *
712    * @param residue a {@link ffx.potential.bonded.Residue} object.
713    * @param C1s a {@link ffx.potential.bonded.Atom} object.
714    * @param O4s a {@link ffx.potential.bonded.Atom} object.
715    * @param C2s a {@link ffx.potential.bonded.Atom} object.
716    * @param glyco a double.
717    * @param forceField a {@link ForceField} object.
718    * @param bondList a {@link java.util.List} object.
719    * @return a {@link ffx.potential.bonded.Residue} object.
720    */
721   private static Residue buildDCY(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
722       ForceField forceField, List<Bond> bondList) {
723 
724     Atom N1 = buildHeavy(residue, "N1", C1s, 1.48, O4s, 108.1, C2s, 113.7, 1, 1191, forceField,
725         bondList);
726     Atom C2 = buildHeavy(residue, "C2", N1, 1.37, C1s, 117.8, O4s, glyco, 0, 1192, forceField,
727         bondList);
728     Atom O2 = buildHeavy(residue, "O2", C2, 1.24, N1, 118.9, C1s, 0.0, 0, 1197, forceField,
729         bondList);
730     Atom N3 = buildHeavy(residue, "N3", C2, 1.38, N1, 118.7, C1s, 180, 0, 1193, forceField,
731         bondList);
732     Atom C4 = buildHeavy(residue, "C4", N3, 1.34, C2, 120.6, N1, 0.0, 0, 1194, forceField, bondList);
733     Atom N4 = buildHeavy(residue, "N4", C4, 1.32, N3, 118.3, C2, 180.0, 0, 1198, forceField,
734         bondList);
735     Atom C5 = buildHeavy(residue, "C5", C4, 1.43, N3, 121.6, C2, 0.0, 0, 1195, forceField, bondList);
736     Atom C6 = buildHeavy(residue, "C6", C5, 1.36, C4, 116.9, N3, 0.0, 0, 1196, forceField, bondList);
737     buildBond(C6, N1, forceField, bondList);
738     buildH(residue, "H41", N4, 1.00e0, C4, 120.0e0, N3, 0.0e0, 0, 1199, forceField, bondList);
739     buildH(residue, "H42", N4, 1.00e0, C4, 120.0e0, N3, 180.0e0, 0, 1200, forceField, bondList);
740     buildH(residue, "H5", C5, 1.08e0, C4, 121.6e0, N3, 180.0e0, 0, 1201, forceField, bondList);
741     buildH(residue, "H6", C6, 1.08e0, C5, 119.4e0, C4, 180.0e0, 0, 1202, forceField, bondList);
742     return residue;
743   }
744 
745   /**
746    * buildDGU.
747    *
748    * @param residue a {@link ffx.potential.bonded.Residue} object.
749    * @param C1s a {@link ffx.potential.bonded.Atom} object.
750    * @param O4s a {@link ffx.potential.bonded.Atom} object.
751    * @param C2s a {@link ffx.potential.bonded.Atom} object.
752    * @param glyco a double.
753    * @param forceField a {@link ForceField} object.
754    * @param bondList a {@link java.util.List} object.
755    * @return a {@link ffx.potential.bonded.Residue} object.
756    */
757   private static Residue buildDGU(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
758       ForceField forceField, List<Bond> bondList) {
759 
760     Atom N9 = buildHeavy(residue, "N9", C1s, 1.48, O4s, 108.1, C2s, 113.7, 1, 1161, forceField,
761         bondList);
762     Atom C8 = buildHeavy(residue, "C8", N9, 1.38, C1s, 128.4, O4s, glyco + 180, 0, 1165, forceField,
763         bondList);
764     Atom N7 = buildHeavy(residue, "N7", C8, 1.31, N9, 114.0, C1s, 180.0, 0, 1164, forceField,
765         bondList);
766     Atom C5 = buildHeavy(residue, "C5", N7, 1.39, C8, 103.8, N9, 0.0, 0, 1163, forceField, bondList);
767     Atom C6 = buildHeavy(residue, "C6", C5, 1.40, N7, 130.1, C8, 180.0, 0, 1169, forceField,
768         bondList);
769     Atom O6 = buildHeavy(residue, "O6", C6, 1.23, C5, 128.8, N7, 0.0, 0, 1174, forceField, bondList);
770     Atom N1 = buildHeavy(residue, "N1", C6, 1.40, C5, 111.4, N7, 180.0, 0, 1168, forceField,
771         bondList);
772     Atom C2 = buildHeavy(residue, "C2", N1, 1.38, C6, 125.2, C5, 0.0, 0, 1167, forceField, bondList);
773     Atom N2 = buildHeavy(residue, "N2", C2, 1.34, N1, 116.1, C6, 180.0, 0, 1171, forceField,
774         bondList);
775     Atom N3 = buildHeavy(residue, "N3", C2, 1.33, N1, 123.3, O6, 0.0, 0, 1166, forceField, bondList);
776     Atom C4 = buildHeavy(residue, "C4", N3, 1.36, C2, 112.3, N1, 0.0, 0, 1162, forceField, bondList);
777     buildBond(C4, C5, forceField, bondList);
778     buildBond(C4, N9, forceField, bondList);
779     buildH(residue, "H8", C8, 1.08e0, N7, 123.0e0, C5, 180.0e0, 0, 1175, forceField, bondList);
780     buildH(residue, "H1", N1, 1.00e0, C6, 117.4e0, C5, 180.0e0, 0, 1170, forceField, bondList);
781     buildH(residue, "H21", N2, 1.00e0, C2, 120.0e0, N1, 0.0e0, 0, 1172, forceField, bondList);
782     buildH(residue, "H22", N2, 1.00e0, C2, 120.0e0, N1, 180.0e0, 0, 1173, forceField, bondList);
783     return residue;
784   }
785 
786   /**
787    * buildDTY.
788    *
789    * @param residue a {@link ffx.potential.bonded.Residue} object.
790    * @param C1s a {@link ffx.potential.bonded.Atom} object.
791    * @param O4s a {@link ffx.potential.bonded.Atom} object.
792    * @param C2s a {@link ffx.potential.bonded.Atom} object.
793    * @param glyco a double.
794    * @param forceField a {@link ForceField} object.
795    * @param bondList a {@link java.util.List} object.
796    * @return a {@link ffx.potential.bonded.Residue} object.
797    */
798   private static Residue buildDTY(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
799       ForceField forceField, List<Bond> bondList) {
800     Atom N1 = buildHeavy(residue, "N1", C1s, 1.48, O4s, 108.1, C2s, 113.7, 1, 1218, forceField,
801         bondList);
802     Atom C2 = buildHeavy(residue, "C2", N1, 1.37, C1s, 117.1, O4s, glyco, 0, 1219, forceField,
803         bondList);
804     Atom O2 = buildHeavy(residue, "O2", C2, 1.22, N1, 122.9, C1s, 0.0, 0, 1224, forceField,
805         bondList);
806     Atom N3 = buildHeavy(residue, "N3", C2, 1.38, N1, 115.4, C1s, 180.0, 0, 1220, forceField,
807         bondList);
808     Atom C4 = buildHeavy(residue, "C4", N3, 1.38, C2, 126.4, N1, 0.0, 0, 1221, forceField, bondList);
809     Atom O4 = buildHeavy(residue, "O4", C4, 1.23, N3, 120.5, C2, 180.0, 0, 1226, forceField,
810         bondList);
811     Atom C5 = buildHeavy(residue, "C5", C4, 1.44, N3, 114.1, C2, 0.0, 0, 1222, forceField, bondList);
812     Atom C7 = buildHeavy(residue, "C7", C5, 1.50, C4, 117.5, N3, 180.0, 0, 1227, forceField,
813         bondList);
814     Atom C6 = buildHeavy(residue, "C6", C5, 1.34, C4, 120.8, N3, 0.0, 0, 1223, forceField, bondList);
815     buildBond(C6, N1, forceField, bondList);
816     buildH(residue, "H3", N3, 1.00e0, C2, 116.8e0, N1, 180.0e0, 0, 1225, forceField, bondList);
817     Atom H = buildH(residue, "H71", C7, 1.09e0, C5, 109.5e0, C4, 0.0e0, 0, 1228, forceField,
818         bondList);
819     buildH(residue, "H72", C7, 1.09e0, C5, 109.5e0, H, 109.5e0, 1, 1228, forceField, bondList);
820     buildH(residue, "H73", C7, 1.09e0, C5, 109.5e0, H, 109.5e0, -1, 1228, forceField, bondList);
821     buildH(residue, "H6", C6, 1.08e0, C5, 119.4e0, C4, 180.0e0, 0, 1229, forceField, bondList);
822     return residue;
823   }
824 
825   /**
826    * buildGUA.
827    *
828    * @param residue a {@link ffx.potential.bonded.Residue} object.
829    * @param C1s a {@link ffx.potential.bonded.Atom} object.
830    * @param O4s a {@link ffx.potential.bonded.Atom} object.
831    * @param C2s a {@link ffx.potential.bonded.Atom} object.
832    * @param glyco a double.
833    * @param forceField a {@link ForceField} object.
834    * @param bondList a {@link java.util.List} object.
835    * @return a {@link ffx.potential.bonded.Residue} object.
836    */
837   private static Residue buildGUA(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
838       ForceField forceField, List<Bond> bondList) {
839     Atom N9 = buildHeavy(residue, "N9", C1s, 1.48, O4s, 108.1, C2s, 113.7, 1, 1047, forceField,
840         bondList);
841     Atom C8 = buildHeavy(residue, "C8", N9, 1.38, C1s, 128.4, O4s, glyco + 180, 0, 1051, forceField,
842         bondList);
843     Atom N7 = buildHeavy(residue, "N7", C8, 1.31, N9, 114.0, C1s, 180.0, 0, 1050, forceField,
844         bondList);
845     Atom C5 = buildHeavy(residue, "C5", N7, 1.39, C8, 103.8, N9, 0.0, 0, 1049, forceField, bondList);
846     Atom C6 = buildHeavy(residue, "C6", C5, 1.40, N7, 130.1, C8, 180.0, 0, 1055, forceField,
847         bondList);
848     Atom O6 = buildHeavy(residue, "O6", C6, 1.23, C5, 128.8, N7, 0.0, 0, 1060, forceField, bondList);
849     Atom N1 = buildHeavy(residue, "N1", C6, 1.40, C5, 111.4, N7, 180.0, 0, 1054, forceField,
850         bondList);
851     Atom C2 = buildHeavy(residue, "C2", N1, 1.38, C6, 125.2, C5, 0.0, 0, 1053, forceField, bondList);
852     Atom N2 = buildHeavy(residue, "N2", C2, 1.34, N1, 116.1, C6, 180.0, 0, 1057, forceField,
853         bondList);
854     Atom N3 = buildHeavy(residue, "N3", C2, 1.33, N1, 123.3, O6, 0.0, 0, 1052, forceField, bondList);
855     Atom C4 = buildHeavy(residue, "C4", N3, 1.36, C2, 112.3, N1, 0.0, 0, 1048, forceField, bondList);
856     buildBond(C4, C5, forceField, bondList);
857     buildBond(C4, N9, forceField, bondList);
858     buildH(residue, "H8", C8, 1.08e0, N7, 123.0e0, C5, 180.0e0, 0, 1061, forceField, bondList);
859     buildH(residue, "H1", N1, 1.00e0, C6, 117.4e0, C5, 180.0e0, 0, 1056, forceField, bondList);
860     buildH(residue, "H21", N2, 1.00e0, C2, 120.0e0, N1, 0.0e0, 0, 1058, forceField, bondList);
861     buildH(residue, "H22", N2, 1.00e0, C2, 120.0e0, N1, 180.0e0, 0, 1059, forceField, bondList);
862     return residue;
863   }
864 
865   /**
866    * buildH2U.
867    *
868    * @param residue a {@link ffx.potential.bonded.Residue} object.
869    * @param C1s a {@link ffx.potential.bonded.Atom} object.
870    * @param forceField a {@link ForceField} object.
871    * @param bondList a {@link java.util.List} object.
872    * @return a {@link ffx.potential.bonded.Residue} object.
873    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
874    */
875   private static Residue buildH2U(Residue residue, Atom C1s, ForceField forceField,
876       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
877 
878     Atom N1 = buildHeavy(residue, "N1", C1s, 1350, forceField, bondList);
879     Atom C2 = buildHeavy(residue, "C2", N1, 1351, forceField, bondList);
880     Atom O2 = buildHeavy(residue, "O2", C2, 1356, forceField, bondList);
881     Atom N3 = buildHeavy(residue, "N3", C2, 1352, forceField, bondList);
882     Atom C4 = buildHeavy(residue, "C4", N3, 1353, forceField, bondList);
883     Atom O4 = buildHeavy(residue, "O4", C4, 1358, forceField, bondList);
884     Atom C5 = buildHeavy(residue, "C5", C4, 1354, forceField, bondList);
885     Atom C6 = buildHeavy(residue, "C6", C5, 1355, forceField, bondList);
886     buildBond(C6, N1, forceField, bondList);
887     buildH(residue, "H3", N3, 1.00e0, C2, 116.5e0, N1, 180.0e0, 0, 1357, forceField, bondList);
888     buildH(residue, "H51", C5, 1.08e0, C4, 109.5e0, C6, 109.5e0, 1, 1359, forceField, bondList);
889     buildH(residue, "H52", C5, 1.08e0, C4, 109.5e0, C6, 109.5e0, -1, 1360, forceField, bondList);
890     buildH(residue, "H61", C6, 1.08e0, C5, 109.5e0, N1, 109.5e0, 1, 1361, forceField, bondList);
891     buildH(residue, "H62", C6, 1.08e0, C5, 109.5e0, N1, 109.5e0, -1, 1362, forceField, bondList);
892     return residue;
893   }
894 
895   /**
896    * buildM1MA.
897    *
898    * @param residue a {@link ffx.potential.bonded.Residue} object.
899    * @param C1s a {@link ffx.potential.bonded.Atom} object.
900    * @param forceField a {@link ForceField} object.
901    * @param bondList a {@link java.util.List} object.
902    * @return a {@link ffx.potential.bonded.Residue} object.
903    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
904    */
905   private static Residue buildM1MA(Residue residue, Atom C1s, ForceField forceField,
906       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
907     Atom N9 = buildHeavy(residue, "N9", C1s, 1605, forceField, bondList);
908     Atom C8 = buildHeavy(residue, "C8", N9, 1609, forceField, bondList);
909     Atom N7 = buildHeavy(residue, "N7", C8, 1608, forceField, bondList);
910     Atom C5 = buildHeavy(residue, "C5", N7, 1607, forceField, bondList);
911     Atom C6 = buildHeavy(residue, "C6", C5, 1613, forceField, bondList);
912     Atom N6 = buildHeavy(residue, "N6", C6, 1615, forceField, bondList);
913     Atom N1 = buildHeavy(residue, "N1", C6, 1612, forceField, bondList);
914     Atom C2 = buildHeavy(residue, "C2", N1, 1611, forceField, bondList);
915     Atom N3 = buildHeavy(residue, "N3", C2, 1610, forceField, bondList);
916     Atom C4 = buildHeavy(residue, "C4", N3, 1606, forceField, bondList);
917     Atom CM1 = buildHeavy(residue, "CM1", N1, 1619, forceField, bondList);
918     buildBond(C4, C5, forceField, bondList);
919     buildBond(C4, N9, forceField, bondList);
920     buildH(residue, "H2", C2, 1.08e0, N3, 115.4e0, C4, 180.0e0, 0, 1614, forceField, bondList);
921     buildH(residue, "H6", C6, 1.08e0, C5, 109.5e0, C4, 180.0e0, 0, 1623, forceField, bondList);
922     buildH(residue, "H8", C8, 1.08e0, N7, 123.1e0, C5, 180.0e0, 0, 1618, forceField, bondList);
923     buildH(residue, "HN61", N6, 1.00e0, C6, 109.5e0, C5, 0.0e0, 0, 1616, forceField, bondList);
924     buildH(residue, "HN62", N6, 1.00e0, C6, 109.5e0, C5, 109.5e0, 0, 1617, forceField, bondList);
925     Atom HM11 = buildH(residue, "HM11", CM1, 1.08e0, N1, 109.5e0, C2, 0.0e0, 0, 1620, forceField,
926         bondList);
927     buildH(residue, "HM12", CM1, 1.08e0, N1, 109.5e0, HM11, 109.5e0, 1, 1621, forceField, bondList);
928     buildH(residue, "HM13", CM1, 1.08e0, N1, 109.5e0, HM11, 109.5e0, -1, 1622, forceField, bondList);
929     return residue;
930   }
931 
932   /**
933    * buildM2G.
934    *
935    * @param residue a {@link ffx.potential.bonded.Residue} object.
936    * @param C1s a {@link ffx.potential.bonded.Atom} object.
937    * @param forceField a {@link ForceField} object.
938    * @param bondList a {@link java.util.List} object.
939    * @return a {@link ffx.potential.bonded.Residue} object.
940    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
941    */
942   private static Residue buildM2G(Residue residue, Atom C1s, ForceField forceField,
943       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
944     Atom N9 = buildHeavy(residue, "N9", C1s, 1379, forceField, bondList);
945     Atom C8 = buildHeavy(residue, "C8", N9, 1383, forceField, bondList);
946     Atom N7 = buildHeavy(residue, "N7", C8, 1382, forceField, bondList);
947     Atom C5 = buildHeavy(residue, "C5", N7, 1381, forceField, bondList);
948     Atom C6 = buildHeavy(residue, "C6", C5, 1387, forceField, bondList);
949     Atom O6 = buildHeavy(residue, "O6", C6, 1390, forceField, bondList);
950     Atom N1 = buildHeavy(residue, "N1", C6, 1386, forceField, bondList);
951     Atom C2 = buildHeavy(residue, "C2", N1, 1385, forceField, bondList);
952     Atom N2 = buildHeavy(residue, "N2", C2, 1389, forceField, bondList);
953     Atom N3 = buildHeavy(residue, "N3", C2, 1384, forceField, bondList);
954     Atom C4 = buildHeavy(residue, "C4", N3, 1380, forceField, bondList);
955     Atom CM1 = buildHeavy(residue, "CM1", N2, 1392, forceField, bondList);
956     Atom CM2 = buildHeavy(residue, "CM2", N2, 1396, forceField, bondList);
957     buildBond(C4, C5, forceField, bondList);
958     buildBond(C4, N9, forceField, bondList);
959     buildH(residue, "H8", C8, 1.08e0, N7, 123.0e0, C5, 180.0e0, 0, 1391, forceField, bondList);
960     buildH(residue, "H1", N1, 1.00e0, C6, 117.4e0, C5, 180.0e0, 0, 1388, forceField, bondList);
961     Atom HM11 = buildH(residue, "HM11", CM1, 1.08e0, N2, 109.5e0, C2, 0.0e0, 0, 1393, forceField,
962         bondList);
963     buildH(residue, "HM12", CM1, 1.08e0, N2, 109.5e0, HM11, 109.5e0, 1, 1394, forceField, bondList);
964     buildH(residue, "HM13", CM1, 1.08e0, N2, 109.5e0, HM11, 109.5e0, -1, 1395, forceField, bondList);
965     Atom HM21 = buildH(residue, "HM21", CM2, 1.08e0, N2, 109.5e0, C2, 0.0e0, 0, 1397, forceField,
966         bondList);
967     buildH(residue, "HM22", CM2, 1.08e0, N2, 109.5e0, HM21, 109.5e0, 1, 1398, forceField, bondList);
968     buildH(residue, "HM23", CM2, 1.08e0, N2, 109.5e0, HM21, 109.5e0, -1, 1399, forceField, bondList);
969     return residue;
970   }
971 
972   /**
973    * buildM2MG.
974    *
975    * @param residue a {@link ffx.potential.bonded.Residue} object.
976    * @param C1s a {@link ffx.potential.bonded.Atom} object.
977    * @param forceField a {@link ForceField} object.
978    * @param bondList a {@link java.util.List} object.
979    * @return a {@link ffx.potential.bonded.Residue} object.
980    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
981    */
982   private static Residue buildM2MG(Residue residue, Atom C1s, ForceField forceField,
983       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
984     Atom N9 = buildHeavy(residue, "N9", C1s, 1316, forceField, bondList);
985     Atom C8 = buildHeavy(residue, "C8", N9, 1320, forceField, bondList);
986     Atom N7 = buildHeavy(residue, "N7", C8, 1319, forceField, bondList);
987     Atom C5 = buildHeavy(residue, "C5", N7, 1318, forceField, bondList);
988     Atom C6 = buildHeavy(residue, "C6", C5, 1324, forceField, bondList);
989     Atom O6 = buildHeavy(residue, "O6", C6, 1328, forceField, bondList);
990     Atom N1 = buildHeavy(residue, "N1", C6, 1323, forceField, bondList);
991     Atom C2 = buildHeavy(residue, "C2", N1, 1322, forceField, bondList);
992     Atom N2 = buildHeavy(residue, "N2", C2, 1326, forceField, bondList);
993     Atom N3 = buildHeavy(residue, "N3", C2, 1321, forceField, bondList);
994     Atom C4 = buildHeavy(residue, "C4", N3, 1317, forceField, bondList);
995     Atom CM2 = buildHeavy(residue, "CM2", N2, 1330, forceField, bondList);
996     buildBond(C4, C5, forceField, bondList);
997     buildBond(C4, N9, forceField, bondList);
998     buildH(residue, "H8", C8, 1.08e0, N7, 123.0e0, C5, 180.0e0, 0, 1329, forceField, bondList);
999     buildH(residue, "H1", N1, 1.00e0, C6, 117.4e0, C5, 180.0e0, 0, 1325, forceField, bondList);
1000     buildH(residue, "H2", N2, 1.00e0, C2, 120.0e0, N1, 0.0e0, 0, 1327, forceField, bondList);
1001     Atom HM21 = buildH(residue, "HM21", CM2, 1.08e0, N2, 109.5e0, C2, 0.0e0, 0, 1331, forceField,
1002         bondList);
1003     buildH(residue, "HM22", CM2, 1.08e0, N2, 109.5e0, HM21, 109.5e0, 1, 1332, forceField, bondList);
1004     buildH(residue, "HM23", CM2, 1.08e0, N2, 109.5e0, HM21, 109.5e0, -1, 1333, forceField, bondList);
1005     return residue;
1006   }
1007 
1008   /**
1009    * buildM5MC.
1010    *
1011    * @param residue a {@link ffx.potential.bonded.Residue} object.
1012    * @param C1s a {@link ffx.potential.bonded.Atom} object.
1013    * @param forceField a {@link ForceField} object.
1014    * @param bondList a {@link java.util.List} object.
1015    * @return a {@link ffx.potential.bonded.Residue} object.
1016    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
1017    */
1018   private static Residue buildM5MC(Residue residue, Atom C1s, ForceField forceField,
1019       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
1020     Atom N1 = buildHeavy(residue, "N1", C1s, 1508, forceField, bondList);
1021     Atom C2 = buildHeavy(residue, "C2", N1, 1509, forceField, bondList);
1022     Atom O2 = buildHeavy(residue, "O2", C2, 1514, forceField, bondList);
1023     Atom N3 = buildHeavy(residue, "N3", C2, 1510, forceField, bondList);
1024     Atom C4 = buildHeavy(residue, "C4", N3, 1511, forceField, bondList);
1025     Atom N4 = buildHeavy(residue, "N4", C4, 1515, forceField, bondList);
1026     Atom C5 = buildHeavy(residue, "C5", C4, 1512, forceField, bondList);
1027     Atom C6 = buildHeavy(residue, "C6", C5, 1513, forceField, bondList);
1028     Atom CM5 = buildHeavy(residue, "CM5", C5, 1519, forceField, bondList);
1029     buildBond(C6, N1, forceField, bondList);
1030     buildH(residue, "H41", N4, 1.00e0, C4, 120.0e0, N3, 0.0e0, 0, 1516, forceField, bondList);
1031     buildH(residue, "H42", N4, 1.00e0, C4, 120.0e0, C5, 0.0e0, 0, 1517, forceField, bondList);
1032     buildH(residue, "H6", C6, 1.08e0, C5, 119.4e0, C4, 180.0e0, 0, 1518, forceField, bondList);
1033     Atom HM51 = buildH(residue, "HM51", CM5, 1.08e0, C5, 109.5e0, C4, 0.0e0, 0, 1520, forceField,
1034         bondList);
1035     buildH(residue, "HM52", CM5, 1.08e0, C5, 109.5e0, HM51, 109.5e0, 1, 1521, forceField, bondList);
1036     buildH(residue, "HM53", CM5, 1.08e0, C5, 109.5e0, HM51, 109.5e0, -1, 1522, forceField, bondList);
1037     return residue;
1038   }
1039 
1040   /**
1041    * buildM5MU.
1042    *
1043    * @param residue a {@link ffx.potential.bonded.Residue} object.
1044    * @param C1s a {@link ffx.potential.bonded.Atom} object.
1045    * @param forceField a {@link ForceField} object.
1046    * @param bondList a {@link java.util.List} object.
1047    * @return a {@link ffx.potential.bonded.Residue} object.
1048    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
1049    */
1050   private static Residue buildM5MU(Residue residue, Atom C1s, ForceField forceField,
1051       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
1052 
1053     Atom N1 = buildHeavy(residue, "N1", C1s, 1575, forceField, bondList);
1054     Atom C2 = buildHeavy(residue, "C2", N1, 1576, forceField, bondList);
1055     Atom O2 = buildHeavy(residue, "O2", C2, 1581, forceField, bondList);
1056     Atom N3 = buildHeavy(residue, "N3", C2, 1577, forceField, bondList);
1057     Atom C4 = buildHeavy(residue, "C4", N3, 1578, forceField, bondList);
1058     Atom O4 = buildHeavy(residue, "O4", C4, 1583, forceField, bondList);
1059     Atom C5 = buildHeavy(residue, "C5", C4, 1579, forceField, bondList);
1060     Atom C6 = buildHeavy(residue, "C6", C5, 1580, forceField, bondList);
1061     Atom C5M = buildHeavy(residue, "C5M", C5, 1585, forceField, bondList);
1062     buildBond(C6, N1, forceField, bondList);
1063     buildH(residue, "H3", N3, 1.00e0, C2, 116.5e0, N1, 180.0e0, 0, 1582, forceField, bondList);
1064     buildH(residue, "H6", C6, 1.08e0, C5, 118.6e0, C4, 180.0e0, 0, 1584, forceField, bondList);
1065     Atom H5M1 = buildH(residue, "H5M1", C5M, 1.08e0, C5, 109.5e0, C6, 0.0e0, 0, 1586, forceField,
1066         bondList);
1067     buildH(residue, "H5M2", C5M, 1.08e0, C5, 109.5e0, H5M1, 109.5e0, 1, 1587, forceField, bondList);
1068     buildH(residue, "H5M3", C5M, 1.08e0, C5, 109.5e0, H5M1, 109.5e0, -1, 1588, forceField, bondList);
1069     return residue;
1070   }
1071 
1072   /**
1073    * buildM7MG.
1074    *
1075    * @param residue a {@link ffx.potential.bonded.Residue} object.
1076    * @param C1s a {@link ffx.potential.bonded.Atom} object.
1077    * @param forceField a {@link ForceField} object.
1078    * @param bondList a {@link java.util.List} object.
1079    * @return a {@link ffx.potential.bonded.Residue} object.
1080    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
1081    */
1082   private static Residue buildM7MG(Residue residue, Atom C1s, ForceField forceField,
1083       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
1084     Atom N9 = buildHeavy(residue, "N9", C1s, 1539, forceField, bondList);
1085     Atom C8 = buildHeavy(residue, "C8", N9, 1543, forceField, bondList);
1086     Atom N7 = buildHeavy(residue, "N7", C8, 1542, forceField, bondList);
1087     Atom C5 = buildHeavy(residue, "C5", N7, 1541, forceField, bondList);
1088     Atom C6 = buildHeavy(residue, "C6", C5, 1547, forceField, bondList);
1089     Atom O6 = buildHeavy(residue, "O6", C6, 1552, forceField, bondList);
1090     Atom N1 = buildHeavy(residue, "N1", C6, 1546, forceField, bondList);
1091     Atom C2 = buildHeavy(residue, "C2", N1, 1545, forceField, bondList);
1092     Atom N2 = buildHeavy(residue, "N2", C2, 1549, forceField, bondList);
1093     Atom N3 = buildHeavy(residue, "N3", C2, 1544, forceField, bondList);
1094     Atom C4 = buildHeavy(residue, "C4", N3, 1540, forceField, bondList);
1095     Atom CM7 = buildHeavy(residue, "CM7", N7, 1555, forceField, bondList);
1096     buildBond(C4, C5, forceField, bondList);
1097     buildBond(C4, N9, forceField, bondList);
1098     buildH(residue, "H81", C8, 1.08e0, N7, 109.5e0, N9, 109.5e0, 1, 1553, forceField, bondList);
1099     buildH(residue, "H82", C8, 1.08e0, N7, 109.5e0, N9, 109.5e0, -1, 1554, forceField, bondList);
1100     buildH(residue, "H1", N1, 1.00e0, C6, 117.4e0, C5, 180.0e0, 0, 1548, forceField, bondList);
1101     buildH(residue, "H21", N2, 1.00e0, C2, 120.0e0, N1, 0.0e0, 0, 1550, forceField, bondList);
1102     buildH(residue, "H22", N2, 1.00e0, C2, 120.0e0, N3, 0.0e0, 0, 1551, forceField, bondList);
1103     Atom HM71 = buildH(residue, "HM71", CM7, 1.08e0, N7, 109.5e0, C8, 0.0e0, 0, 1556, forceField,
1104         bondList);
1105     buildH(residue, "HM72", CM7, 1.08e0, N7, 109.5e0, HM71, 109.5e0, 1, 1557, forceField, bondList);
1106     buildH(residue, "HM73", CM7, 1.08e0, N7, 109.5e0, HM71, 109.5e0, -1, 1558, forceField, bondList);
1107     return residue;
1108   }
1109 
1110   /**
1111    * buildOMC.
1112    *
1113    * @param residue a {@link ffx.potential.bonded.Residue} object.
1114    * @param C1s a {@link ffx.potential.bonded.Atom} object.
1115    * @param O4s a {@link ffx.potential.bonded.Atom} object.
1116    * @param C2s a {@link ffx.potential.bonded.Atom} object.
1117    * @param glyco a double.
1118    * @param forceField a {@link ForceField} object.
1119    * @param bondList a {@link java.util.List} object.
1120    * @return a {@link ffx.potential.bonded.Residue} object.
1121    */
1122   private static Residue buildOMC(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
1123       ForceField forceField, List<Bond> bondList) {
1124     return buildCYT(residue, C1s, O4s, C2s, glyco, forceField, bondList);
1125   }
1126 
1127   /**
1128    * buildOMG.
1129    *
1130    * @param residue a {@link ffx.potential.bonded.Residue} object.
1131    * @param C1s a {@link ffx.potential.bonded.Atom} object.
1132    * @param O4s a {@link ffx.potential.bonded.Atom} object.
1133    * @param C2s a {@link ffx.potential.bonded.Atom} object.
1134    * @param glyco a double.
1135    * @param forceField a {@link ForceField} object.
1136    * @param bondList a {@link java.util.List} object.
1137    * @return a {@link ffx.potential.bonded.Residue} object.
1138    */
1139   private static Residue buildOMG(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
1140       ForceField forceField, List<Bond> bondList) {
1141     return buildGUA(residue, C1s, O4s, C2s, glyco, forceField, bondList);
1142   }
1143 
1144   /**
1145    * Builds a missing OP3 atom. Can be applied either to 5'-terminal or 3'-terminal phosphates.
1146    *
1147    * @param residue Residue to build an OP3 for.
1148    * @param phosphate The phosphate atom.
1149    * @param aType The atom type to use.
1150    * @param forceField Force field in use.
1151    * @param bondList List of bonds.
1152    * @param adjacentResidue A Residue either 5' or 3' of residue.
1153    * @param at5prime If this residue is at the 5' terminus (vs. 3').
1154    * @throws MissingHeavyAtomException Thrown if a needed heavy atom is missing.
1155    */
1156   private static void buildOP3(Residue residue, Atom phosphate, int aType, ForceField forceField,
1157       List<Bond> bondList, Residue adjacentResidue, boolean at5prime)
1158       throws MissingHeavyAtomException, MissingAtomTypeException {
1159     Atom P = null;
1160     Atom OP1 = null;
1161     Atom OP2 = null;
1162     Atom riboOxygen = null;
1163     Atom riboCarbon = null;
1164     boolean foundOP3 = false;
1165 
1166     List<Atom> residueAtoms = residue.getAtomList();
1167     // Label for a break statement: OP3 causes a break from a switch to the surrounding for loop.
1168     AtomLoop:
1169     for (Atom at : residueAtoms) {
1170       String atName = at.getName().toUpperCase();
1171       switch (atName) {
1172         case "OP3":
1173           buildHeavy(residue, "OP3", phosphate, aType, forceField, bondList);
1174           foundOP3 = true;
1175           break AtomLoop;
1176         case "OP1":
1177           OP1 = at;
1178           break;
1179         case "OP2":
1180           OP2 = at;
1181           break;
1182         case "P":
1183           P = at;
1184           break;
1185         case "O5'":
1186           riboOxygen = at;
1187           break;
1188         case "C5'":
1189           riboCarbon = at;
1190           break;
1191       }
1192     }
1193 
1194     if (!at5prime) {
1195       riboOxygen = (Atom) adjacentResidue.getAtomNode("O3'");
1196       riboCarbon = (Atom) adjacentResidue.getAtomNode("C3'");
1197       if (riboCarbon == null || riboOxygen == null) {
1198         logger.severe(format(
1199             " Could not find either O3' " + "or C3' in residue %s prior to presumed 3' "
1200                 + "phosphate cap %s", adjacentResidue, residue));
1201       }
1202     }
1203 
1204     if (!foundOP3) {
1205       logger.fine(
1206           format(" EXPERIMENTAL: OP3 of residue %s not found, being rebuilt based on ideal geometry",
1207               residue));
1208       if (P == null || OP1 == null || OP2 == null) {
1209         throw new IllegalArgumentException(format(
1210             " Attempted to build OP3 for residue %s, but one of P, OP1, OP2, O5', or C5' were null",
1211             residue));
1212       }
1213       if (at5prime) {
1214         if (riboOxygen == null) {
1215           logger.fine(" Attempting to find O5' of subsequent residue");
1216           riboOxygen = getNextResAtom(residue, adjacentResidue, "O5'");
1217         }
1218         if (riboCarbon == null) {
1219           logger.fine(" Attempting to find C5' of subsequent residue");
1220           riboCarbon = getNextResAtom(residue, adjacentResidue, "C5'");
1221         }
1222       }
1223 
1224       // Borrow the P-OP1 distance.
1225       double[] xyzOP1 = new double[3];
1226       double[] xyzP = new double[3];
1227       xyzOP1 = OP1.getXYZ(xyzOP1);
1228       xyzP = P.getXYZ(xyzP);
1229       double distance = dist(xyzP, xyzOP1);
1230 
1231       // Borrow the O5'-P-OP1 angle.
1232       double[] xyzRiboO = new double[3];
1233       xyzRiboO = riboOxygen.getXYZ(xyzRiboO);
1234       double angle = toDegrees(DoubleMath.bondAngle(xyzRiboO, xyzP, xyzOP1));
1235 
1236       // Borrow the C5'-O5'-P-OP1 dihedral (which will have +120 or +240 degrees added on).
1237       double[] xyzRiboC = new double[3];
1238       xyzRiboC = riboCarbon.getXYZ(xyzRiboC);
1239       double dihedralOP1 = toDegrees(DoubleMath.dihedralAngle(xyzRiboC, xyzRiboO, xyzP, xyzOP1));
1240 
1241       Atom OP3 = buildH(residue, "OP3", P, distance, riboOxygen, angle, riboCarbon,
1242           dihedralOP1 + 120, 0, aType, forceField, bondList);
1243 
1244       // Measure OP3-OP2 distance for test dihedral + 120 degrees.
1245       double[] xyzOP2 = new double[3];
1246       xyzOP2 = OP2.getXYZ(xyzOP2);
1247       double[] xyzChiral1 = new double[3];
1248       xyzChiral1 = OP3.getXYZ(xyzChiral1);
1249       double distChiral1 = dist(xyzChiral1, xyzOP2);
1250 
1251       // Measure OP3-OP2 distance for test dihedral + 240 degrees.
1252       intxyz(OP3, P, distance, riboOxygen, angle, riboCarbon, dihedralOP1 + 240, 0);
1253       double[] xyzChiral2 = new double[3];
1254       xyzChiral2 = OP3.getXYZ(xyzChiral2);
1255       double distChiral2 = dist(xyzChiral2, xyzOP2);
1256 
1257       if (logger.isLoggable(Level.FINE)) {
1258         logger.fine(
1259             format(" Bond: %10.5f Angle: %10.5f Dihedral: %10.5f", distance, angle, dihedralOP1));
1260         logger.fine(format(" OP2 position: %s", Arrays.toString(xyzOP2)));
1261         logger.fine(format(" Position 1: %10.6g %10.6g %10.6g with distance %10.6g", xyzChiral1[0],
1262             xyzChiral1[1], xyzChiral1[2], distChiral1));
1263         logger.fine(format(" Position 2: %10.6g %10.6g %10.6g with distance %10.6g", xyzChiral2[0],
1264             xyzChiral2[1], xyzChiral2[2], distChiral2));
1265       }
1266       if (distChiral1 > distChiral2) {
1267         logger.fine(" Picked dihedral +120");
1268         OP3.setXYZ(xyzChiral1);
1269       } else {
1270         logger.fine(" Picked dihedral +240");
1271         OP3.setXYZ(xyzChiral2);
1272       }
1273     }
1274   }
1275 
1276   /**
1277    * buildPSU.
1278    *
1279    * @param residue a {@link ffx.potential.bonded.Residue} object.
1280    * @param C1s a {@link ffx.potential.bonded.Atom} object.
1281    * @param forceField a {@link ForceField} object.
1282    * @param bondList a {@link java.util.List} object.
1283    * @return a {@link ffx.potential.bonded.Residue} object.
1284    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
1285    */
1286   private static Residue buildPSU(Residue residue, Atom C1s, ForceField forceField,
1287       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
1288     // C1s bonds to C5 in PsuedoUridine
1289     Atom C5 = buildHeavy(residue, "C5", C1s, 1485, forceField, bondList);
1290     Atom C6 = buildHeavy(residue, "C6", C5, 1486, forceField, bondList);
1291     Atom N1 = buildHeavy(residue, "N1", C6, 1481, forceField, bondList);
1292     Atom C2 = buildHeavy(residue, "C2", N1, 1482, forceField, bondList);
1293     Atom O2 = buildHeavy(residue, "O2", C2, 1487, forceField, bondList);
1294     Atom N3 = buildHeavy(residue, "N3", C2, 1483, forceField, bondList);
1295     Atom C4 = buildHeavy(residue, "C4", N3, 1484, forceField, bondList);
1296     Atom O4 = buildHeavy(residue, "O4", C4, 1489, forceField, bondList);
1297     buildBond(C4, C5, forceField, bondList);
1298     buildH(residue, "H1", N1, 1.00e0, C2, 120.0e0, O2, 0.0e0, 0, 1491, forceField, bondList);
1299     buildH(residue, "H3", N3, 1.00e0, C2, 120.0e0, O2, 0.0e0, 0, 1488, forceField, bondList);
1300     buildH(residue, "H6", C6, 1.08e0, C5, 120.0e0, C1s, 0.0e0, 0, 1490, forceField, bondList);
1301     return residue;
1302   }
1303 
1304   /**
1305    * Builds a sugar O5 atom. No bonds are made to it.
1306    *
1307    * @param residue The residue to operate on.
1308    * @param lookUp Biotype lookup.
1309    * @param forceField The ForceField to use.
1310    * @return The created sugar O5 atoms.
1311    */
1312   private static Atom buildSugarO5(Residue residue, int lookUp, ForceField forceField) {
1313     Atom C5 = residue.getAtomByName("C5'", true);
1314     Atom C4 = residue.getAtomByName("C4'", true);
1315     Atom C3 = residue.getAtomByName("C3'", true);
1316     return buildHeavy(residue, "O5'", C5, 1.43, C4, 109.5, C3, 180.0, 0, 1244, forceField);
1317   }
1318 
1319   /**
1320    * buildURI.
1321    *
1322    * @param residue a {@link ffx.potential.bonded.Residue} object.
1323    * @param C1s a {@link ffx.potential.bonded.Atom} object.
1324    * @param O4s a {@link ffx.potential.bonded.Atom} object.
1325    * @param C2s a {@link ffx.potential.bonded.Atom} object.
1326    * @param glyco a double.
1327    * @param forceField a {@link ForceField} object.
1328    * @param bondList a {@link java.util.List} object.
1329    * @return a {@link ffx.potential.bonded.Residue} object.
1330    */
1331   private static Residue buildURI(Residue residue, Atom C1s, Atom O4s, Atom C2s, double glyco,
1332       ForceField forceField, List<Bond> bondList) {
1333     Atom N1 = buildHeavy(residue, "N1", C1s, 1.48, O4s, 108.1, C2s, 113.7, 1, 1106, forceField,
1334         bondList);
1335     Atom C2 = buildHeavy(residue, "C2", N1, 1.38, C1s, 117.1, O4s, glyco, 0, 1107, forceField,
1336         bondList);
1337     Atom O2 = buildHeavy(residue, "O2", C2, 1.22, N1, 123.2, C1s, 0.0, 0, 1112, forceField,
1338         bondList);
1339     Atom N3 = buildHeavy(residue, "N3", C2, 1.37, N1, 114.8, C1s, 180.0, 0, 1108, forceField,
1340         bondList);
1341     Atom C4 = buildHeavy(residue, "C4", N3, 1.38, C2, 127.0, N1, 0.0, 0, 1109, forceField, bondList);
1342     Atom O4 = buildHeavy(residue, "O4", C4, 1.23, N3, 119.8, C2, 180.0, 0, 1114, forceField,
1343         bondList);
1344     Atom C5 = buildHeavy(residue, "C5", C4, 1.44, N3, 114.7, C2, 0.0, 0, 1110, forceField, bondList);
1345     Atom C6 = buildHeavy(residue, "C6", C5, 1.34, O4, 119.2, C4, 0.0, 0, 1111, forceField, bondList);
1346     buildBond(C6, N1, forceField, bondList);
1347     buildH(residue, "H3", N3, 1.00e0, C2, 116.5e0, N1, 180.0e0, 0, 1113, forceField, bondList);
1348     buildH(residue, "H5", C5, 1.08e0, C4, 120.4e0, N3, 180.0e0, 0, 1115, forceField, bondList);
1349     buildH(residue, "H6", C6, 1.08e0, C5, 118.6e0, C4, 180.0e0, 0, 1116, forceField, bondList);
1350 
1351     return residue;
1352   }
1353 
1354   /**
1355    * buildYYG.
1356    *
1357    * @param residue a {@link ffx.potential.bonded.Residue} object.
1358    * @param C1s a {@link ffx.potential.bonded.Atom} object.
1359    * @param forceField a {@link ForceField} object.
1360    * @param bondList a {@link java.util.List} object.
1361    * @return a {@link ffx.potential.bonded.Residue} object.
1362    * @throws ffx.potential.bonded.BondedUtils.MissingHeavyAtomException if any.
1363    */
1364   private static Residue buildYYG(Residue residue, Atom C1s, ForceField forceField,
1365       List<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
1366     Atom N9 = buildHeavy(residue, "N9", C1s, 1640, forceField, bondList);
1367     Atom C8 = buildHeavy(residue, "C8", N9, 1644, forceField, bondList);
1368     Atom N7 = buildHeavy(residue, "N7", C8, 1643, forceField, bondList);
1369     Atom C5 = buildHeavy(residue, "C5", N7, 1642, forceField, bondList);
1370     Atom C6 = buildHeavy(residue, "C6", C5, 1648, forceField, bondList);
1371     Atom O6 = buildHeavy(residue, "O6", C6, 1650, forceField, bondList);
1372     Atom N1 = buildHeavy(residue, "N1", C6, 1647, forceField, bondList);
1373     Atom C2 = buildHeavy(residue, "C2", N1, 1646, forceField, bondList);
1374     Atom N2 = buildHeavy(residue, "N2", C2, 1649, forceField, bondList);
1375     Atom N3 = buildHeavy(residue, "N3", C2, 1645, forceField, bondList);
1376     Atom C3 = buildHeavy(residue, "C3", N3, 1652, forceField, bondList);
1377     Atom C4 = buildHeavy(residue, "C4", N3, 1641, forceField, bondList);
1378     Atom C11 = buildHeavy(residue, "C11", N2, 1657, forceField, bondList);
1379     Atom C10 = buildHeavy(residue, "C10", C11, 1658, forceField, bondList);
1380     Atom C12 = buildHeavy(residue, "C12", C11, 1656, forceField, bondList);
1381     Atom C13 = buildHeavy(residue, "C13", C12, 1662, forceField, bondList);
1382     Atom C14 = buildHeavy(residue, "C14", C13, 1665, forceField, bondList);
1383     Atom C15 = buildHeavy(residue, "C15", C14, 1668, forceField, bondList);
1384     Atom C16 = buildHeavy(residue, "C16", C15, 1675, forceField, bondList);
1385     Atom O17 = buildHeavy(residue, "O17", C16, 1676, forceField, bondList);
1386     Atom O18 = buildHeavy(residue, "O18", C16, 1674, forceField, bondList);
1387     Atom C19 = buildHeavy(residue, "C19", O18, 1670, forceField, bondList);
1388     Atom N20 = buildHeavy(residue, "N20", C15, 1677, forceField, bondList);
1389     Atom C21 = buildHeavy(residue, "C21", N20, 1679, forceField, bondList);
1390     Atom O22 = buildHeavy(residue, "O22", C21, 1680, forceField, bondList);
1391     Atom O23 = buildHeavy(residue, "O23", C21, 1681, forceField, bondList);
1392     Atom C24 = buildHeavy(residue, "C24", O23, 1682, forceField, bondList);
1393     buildBond(C4, C5, forceField, bondList);
1394     buildBond(C4, N9, forceField, bondList);
1395     buildBond(N1, C12, forceField, bondList);
1396     buildH(residue, "H8", C8, 1.08e0, N7, 123.0e0, C5, 180.0e0, 0, 1651, forceField, bondList);
1397     Atom H31 = buildH(residue, "H31", C3, 1.08e0, N3, 109.5e0, C4, 0.0e0, 0, 1653, forceField,
1398         bondList);
1399     buildH(residue, "H32", C3, 1.08e0, N3, 109.5e0, H31, 109.5e0, 1, 1654, forceField, bondList);
1400     buildH(residue, "H33", C3, 1.08e0, N3, 109.5e0, H31, 109.5e0, -1, 1655, forceField, bondList);
1401     Atom H101 = buildH(residue, "H101", C10, 1.08e0, C11, 109.5e0, N2, 0.0e0, 0, 1659, forceField,
1402         bondList);
1403     buildH(residue, "H102", C10, 1.08e0, C11, 109.5e0, H101, 109.5e0, 1, 1660, forceField, bondList);
1404     buildH(residue, "H103", C10, 1.08e0, C11, 109.5e0, H101, 109.5e0, -1, 1661, forceField,
1405         bondList);
1406     buildH(residue, "H131", C13, 1.08e0, C12, 109.5e0, C14, 109.5e0, 1, 1663, forceField, bondList);
1407     buildH(residue, "H132", C13, 1.08e0, C12, 109.5e0, C14, 109.5e0, -1, 1664, forceField, bondList);
1408     buildH(residue, "H141", C14, 1.08e0, C13, 109.5e0, C15, 109.5e0, 1, 1666, forceField, bondList);
1409     buildH(residue, "H142", C14, 1.08e0, C13, 109.5e0, C15, 109.5e0, -1, 1667, forceField, bondList);
1410     buildH(residue, "H15", C15, 1.08e0, C14, 109.5e0, O18, 180.e0, 0, 1669, forceField, bondList);
1411     Atom H191 = buildH(residue, "H191", C19, 1.08e0, O18, 109.5e0, C16, 0.0e0, 0, 1671, forceField,
1412         bondList);
1413     buildH(residue, "H192", C19, 1.08e0, O18, 109.5e0, H191, 109.5e0, 1, 1672, forceField, bondList);
1414     buildH(residue, "H193", C19, 1.08e0, O18, 109.5e0, H191, 109.5e0, -1, 1673, forceField,
1415         bondList);
1416     buildH(residue, "HN2", N20, 1.00e0, C15, 109.5e0, O22, 180.0e0, 0, 1678, forceField, bondList);
1417     Atom H241 = buildH(residue, "H241", C24, 1.08e0, O23, 109.5e0, C21, 0.0e0, 0, 1683, forceField,
1418         bondList);
1419     buildH(residue, "H242", C24, 1.08e0, O23, 109.5e0, H241, 109.5e0, 1, 1684, forceField, bondList);
1420     buildH(residue, "H243", C24, 1.08e0, O23, 109.5e0, H241, 109.5e0, -1, 1685, forceField,
1421         bondList);
1422     return residue;
1423   }
1424 
1425   /**
1426    * Intended for 5'-terminal phosphate caps listed as their own residue: find a given atom name in
1427    * the next residue.
1428    *
1429    * @param residue Current residue.
1430    * @param nextResidue Next residue to search in.
1431    * @param atomName Atom name to look for.
1432    * @return The requested atom.
1433    */
1434   private static Atom getNextResAtom(Residue residue, Residue nextResidue, String atomName) {
1435     if (nextResidue == null) {
1436       throw new IllegalArgumentException(
1437           String.format(" Residue %s: No subsequent residue to find atom %s in!", residue,
1438               atomName));
1439     }
1440     List<Atom> nrAtoms = nextResidue.getAtomList();
1441     for (Atom atom : nrAtoms) {
1442       if (atom.getName().equalsIgnoreCase(atomName)) {
1443         return atom;
1444       }
1445     }
1446     throw new IllegalArgumentException(
1447         String.format(" Residue %s: Subsequent residue %s did not contain an atom %s", residue,
1448             nextResidue, atomName));
1449   }
1450 
1451   /**
1452    * getNucleicAcidNumber.
1453    *
1454    * @param residueName a {@link String} object.
1455    * @return The index of the nucleic acid in the nucleicAcidList.
1456    */
1457   public static int getNucleicAcidNumber(String residueName) {
1458     int nucleicAcidNumber = -1;
1459     for (NucleicAcid3 nucleicAcid : nucleicAcidList) {
1460       nucleicAcidNumber++;
1461       if (nucleicAcid.toString().equalsIgnoreCase(residueName)) {
1462         break;
1463       }
1464     }
1465     return nucleicAcidNumber;
1466   }
1467 }