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;
39  
40  import edu.rit.pj.ParallelTeam;
41  import ffx.crystal.Crystal;
42  import ffx.crystal.SymOp;
43  import ffx.potential.bonded.Atom;
44  import ffx.potential.bonded.Bond;
45  import ffx.potential.bonded.MSGroup;
46  import ffx.potential.bonded.MSNode;
47  import ffx.potential.bonded.Molecule;
48  import ffx.potential.bonded.Polymer;
49  import ffx.potential.bonded.RendererCache;
50  import ffx.potential.bonded.Residue;
51  import ffx.potential.bonded.Residue.ResidueType;
52  import ffx.potential.parameters.ForceField;
53  import ffx.utilities.StringUtils;
54  import org.apache.commons.configuration2.CompositeConfiguration;
55  import org.jogamp.java3d.Appearance;
56  import org.jogamp.java3d.BoundingSphere;
57  import org.jogamp.java3d.BranchGroup;
58  import org.jogamp.java3d.ColoringAttributes;
59  import org.jogamp.java3d.GeometryArray;
60  import org.jogamp.java3d.Group;
61  import org.jogamp.java3d.LineArray;
62  import org.jogamp.java3d.LineAttributes;
63  import org.jogamp.java3d.Link;
64  import org.jogamp.java3d.Material;
65  import org.jogamp.java3d.Node;
66  import org.jogamp.java3d.RenderingAttributes;
67  import org.jogamp.java3d.Shape3D;
68  import org.jogamp.java3d.SharedGroup;
69  import org.jogamp.java3d.Switch;
70  import org.jogamp.java3d.Transform3D;
71  import org.jogamp.java3d.TransformGroup;
72  import org.jogamp.java3d.utils.picking.PickTool;
73  import org.jogamp.vecmath.Color3f;
74  import org.jogamp.vecmath.Matrix3d;
75  import org.jogamp.vecmath.Point3d;
76  import org.jogamp.vecmath.Vector3d;
77  
78  import java.awt.GraphicsEnvironment;
79  import java.io.File;
80  import java.io.Serial;
81  import java.util.ArrayList;
82  import java.util.Arrays;
83  import java.util.HashMap;
84  import java.util.HashSet;
85  import java.util.Iterator;
86  import java.util.List;
87  import java.util.ListIterator;
88  import java.util.Set;
89  import java.util.logging.Level;
90  import java.util.logging.Logger;
91  import java.util.stream.Collectors;
92  import java.util.stream.Stream;
93  
94  import static ffx.crystal.SymOp.applyCartesianSymOp;
95  import static ffx.numerics.math.DoubleMath.length;
96  import static ffx.numerics.math.DoubleMath.sub;
97  import static java.lang.String.format;
98  
99  /**
100  * The MolecularAssembly class is a collection of Polymers, Hetero Molecules, Ions and Water
101  *
102  * @author Michael J. Schnieders
103  * @since 1.0
104  */
105 public class MolecularAssembly extends MSGroup {
106 
107   @Serial
108   private static final long serialVersionUID = 1L;
109 
110   private static final Logger logger = Logger.getLogger(MolecularAssembly.class.getName());
111   private static final double[] a = new double[3];
112   private Character alternateLocation = ' ';
113   private final List<String> headerLines = new ArrayList<>();
114   // Data Nodes
115   private final MSNode ions = new MSNode("Ions");
116   private final HashMap<String, Molecule> ionHashMap = new HashMap<>();
117   private final MSNode water = new MSNode("Water");
118   private final HashMap<String, Molecule> waterHashMap = new HashMap<>();
119   private final MSNode molecules = new MSNode("Hetero Molecules");
120   private final HashMap<String, Molecule> moleculeHashMap = new HashMap<>();
121   private final List<BranchGroup> myNewShapes = new ArrayList<>();
122   // MolecularAssembly member variables
123   protected ForceField forceField;
124   // The File that was read to create this MolecularAssembly.
125   private File file;
126   // The File to use for writing archives.
127   private File archiveFile;
128   private ForceFieldEnergy potentialEnergy;
129   private CompositeConfiguration properties;
130   private Vector3d offset;
131   private int cycles = 1;
132   private int currentCycle = 1;
133   // 3D Graphics Nodes - There is a diagram explaining the MolecularAssembly
134   // Scenegraph below
135   private BranchGroup branchGroup;
136   private TransformGroup originToRot;
137   private Transform3D originToRotT3D;
138   private Vector3d originToRotV3D;
139   private TransformGroup rotToCOM;
140   private Transform3D rotToCOMT3D;
141   private Vector3d rotToCOMV3D;
142   private BranchGroup base;
143   private Switch switchGroup;
144   private Shape3D wire;
145   private BranchGroup vrml;
146   private TransformGroup vrmlTG;
147   private Transform3D vrmlTd;
148   private BranchGroup childNodes;
149   private Atom[] atomLookUp;
150   private LineAttributes lineAttributes;
151   private boolean visible = false;
152   private FractionalMode fractionalMode = FractionalMode.MOLECULE;
153   private double[][] fractionalCoordinates;
154   private boolean titrateConformer = false;
155   private Atom atomInitial;
156 
157   /**
158    * Constructor for MolecularAssembly.
159    *
160    * @param name a {@link java.lang.String} object.
161    */
162   public MolecularAssembly(String name) {
163     super(name);
164     getAtomNode().setName("MacroMolecules");
165     add(molecules);
166     add(ions);
167     add(water);
168   }
169 
170   /**
171    * Constructor for MolecularAssembly.
172    *
173    * @param name a {@link java.lang.String} object.
174    * @param Polymers a {@link ffx.potential.bonded.MSNode} object.
175    */
176   public MolecularAssembly(String name, MSNode Polymers) {
177     super(name, Polymers);
178   }
179 
180   /**
181    * Constructor for MolecularAssembly.
182    *
183    * @param name a {@link java.lang.String} object.
184    * @param Polymers a {@link ffx.potential.bonded.MSNode} object.
185    * @param properties a {@link org.apache.commons.configuration2.CompositeConfiguration}
186    *     object.
187    */
188   public MolecularAssembly(String name, MSNode Polymers, CompositeConfiguration properties) {
189     this(name, Polymers);
190     this.properties = properties;
191   }
192 
193   /**
194    * Set the alternate location.
195    * @param alternateLocation The alternate location.
196    */
197   public void setAlternateLocation(Character alternateLocation) {
198     this.alternateLocation = alternateLocation;
199   }
200 
201   /**
202    * Get the alternate location.
203    * @return The alternate location.
204    */
205   public Character getAlternateLocation() {
206     return alternateLocation;
207   }
208 
209   /**
210    * Adds a header line to this MolecularAssembly (particularly for PDB formats)
211    *
212    * @param line Line to add.
213    */
214   public void addHeaderLine(String line) {
215     headerLines.add(line);
216   }
217 
218   /** {@inheritDoc} */
219   @Override
220   public MSNode addMSNode(MSNode o) {
221     List<MSNode> Polymers = getAtomNodeList();
222     if (o instanceof Atom) {
223       Atom atom = (Atom) o;
224       if (atom.isModRes()) {
225         return getResidue(atom, true, Residue.ResidueType.AA);
226       } else if (!atom.isHetero()) {
227         String resName = atom.getResidueName();
228         switch (resName){
229           case "HIS":
230           case "HIE":
231           case "HID":
232           case "ASP":
233           case "ASH":
234           case "GLU":
235           case "GLH":
236           case "LYS":
237           case "LYD":
238             for (Residue residue: getResidueList()){
239               if(residue.getResidueNumber() == atom.getResidueNumber()){
240                 for (Atom currentAtom : residue.getAtomList()) {
241                   if (atom.getResidueNumber() == currentAtom.getResidueNumber() && atom.getChainID() == currentAtom.getChainID()
242                           && atom.getAltLoc() != currentAtom.getAltLoc() && !atom.getResidueName().equals(currentAtom.getResidueName()) &&
243                           atom.getName().equals(currentAtom.getName())) {
244                     titrateConformer = true;
245                     atomInitial = currentAtom;
246                     return getResidue(atom, true);
247                   } else {
248                     titrateConformer = false;
249                     atomInitial = null;
250                   }
251                 }
252               }
253             }
254             break;
255           default:
256             break;
257         }
258         return getResidue(atom, true);
259       } else {
260         return getMolecule(atom, true);
261       }
262     } else if (o instanceof Residue) {
263       Residue residue = (Residue) o;
264       Character chainID = residue.getChainID();
265       String segID = residue.getSegID();
266       int index = Polymers.indexOf(new Polymer(chainID, segID));
267       // See if the polymer already exists.
268       if (index != -1) {
269         Polymer c = (Polymer) Polymers.get(index);
270         setFinalized(false);
271         return c.addMSNode(residue);
272       } else {
273         Polymer newc = new Polymer(chainID, segID);
274         getAtomNode().add(newc);
275         setFinalized(false);
276         return newc.addMSNode(residue);
277       }
278     } else if (o instanceof Polymer) {
279       Polymer c = (Polymer) o;
280       int index = Polymers.indexOf(c);
281       if (index == -1) {
282         getAtomNode().add(c);
283         setFinalized(false);
284         return c;
285       } else {
286         return Polymers.get(index);
287       }
288     } else if (o instanceof Molecule) {
289       Molecule m = (Molecule) o;
290       String key = m.getKey();
291       if (m.getAtomNode().getChildCount() == 1) {
292         ions.add(m);
293         if (ionHashMap.containsKey(key)) {
294           logger.info(" Ion map already contains " + m);
295         } else {
296           ionHashMap.put(m.getKey(), m);
297         }
298         return m;
299       } else if (Utilities.isWaterOxygen((Atom) m.getAtomNode().getChildAt(0))) {
300         water.add(m);
301         if (waterHashMap.containsKey(key)) {
302           logger.info(" Water map already contains " + m);
303         } else {
304           waterHashMap.put(m.getKey(), m);
305         }
306         return m;
307       } else {
308         molecules.add(m);
309         if (moleculeHashMap.containsKey(key)) {
310           logger.info(" Molecule map already contains " + m);
311         } else {
312           moleculeHashMap.put(m.getKey(), m);
313         }
314         return m;
315       }
316     } else {
317       String message = "Programming error in MolecularAssembly addNode";
318       logger.log(Level.SEVERE, message);
319       return o;
320     }
321   }
322 
323   /**
324    * Applies a randomly drawn density to a molecular system's crystal.
325    *
326    * @param ucDensity Target density.
327    */
328   public void applyRandomDensity(double ucDensity) {
329     if (ucDensity > 0) {
330       logger.info(format("\n Applying random unit cell axes with target density %6.3f (g/cc).", ucDensity));
331       Crystal crystal = getCrystal();
332       if (!crystal.aperiodic()) {
333         double mass = getMass();
334         crystal.randomParameters(ucDensity, mass);
335         logger.info(crystal.toString());
336         potentialEnergy.setCrystal(crystal);
337       } else {
338         logger.fine(String.format(" Potential %s is an aperiodic system!", potentialEnergy));
339       }
340     }
341   }
342 
343   /**
344    * Randomizes position in the unit cell of each molecule by applying a Cartesian SymOp with a
345    * random translation.
346    *
347    * @param x SymOp with translation range -x/2 .. x/2 (0 for random placement in the unit cell,
348    *     negative for no SymOp)
349    */
350   public void applyRandomSymOp(double x) {
351     // The Unit Cell is needed here (not the ReplicatesCrystal).
352     Crystal crystal = getCrystal().getUnitCell();
353     if (crystal.aperiodic() || x < 0) {
354       return;
355     }
356     double[] xyz = new double[3];
357     List<MSNode> molecules = getMolecules();
358     int moleculeNum = 1;
359     for (MSNode msNode : molecules) {
360       List<Atom> atoms = msNode.getAtomList();
361       SymOp symOp;
362       if (x == 0) {
363         double[] translation = crystal.getRandomCartTranslation();
364         symOp = SymOp.randomSymOpFactory(translation);
365       } else {
366         symOp = SymOp.randomSymOpFactory(x);
367       }
368       logger.info(
369           format(
370               "\n Applying random Cartesian SymOp to molecule %d:\n%s",
371               moleculeNum, symOp));
372       for (Atom atom : atoms) {
373         atom.getXYZ(xyz);
374         applyCartesianSymOp(xyz, xyz, symOp);
375         atom.setXYZ(xyz);
376       }
377       moleculeNum++;
378     }
379   }
380 
381   /** center */
382   public void center() {
383     double[] center = getMultiScaleCenter(false);
384     offset = new Vector3d(center);
385     if (vrml != null) {
386       vrmlTd.set(offset);
387       vrmlTG.setTransform(vrmlTd);
388     }
389     offset.negate();
390     originToRotV3D.set(offset);
391     originToRotT3D.setTranslation(originToRotV3D);
392     originToRot.setTransform(originToRotT3D);
393     rotToCOMT3D.setIdentity();
394     rotToCOM.setTransform(rotToCOMT3D);
395     offset.negate();
396     rotateAbout(offset);
397     originToRotT3D.get(offset);
398   }
399 
400   /**
401    * centerAt
402    *
403    * @param d an array of double.
404    */
405   public void centerAt(double[] d) {
406     double[] Rc = {0, 0, 0};
407     double[] c = new double[3];
408     ListIterator<Atom> li;
409     int i, num = getAtomList().size();
410     for (li = getAtomList().listIterator(); li.hasNext(); ) {
411       (li.next()).getXYZ(a);
412       Rc[0] += a[0];
413       Rc[1] += a[1];
414       Rc[2] += a[2];
415     }
416     for (i = 0; i < 3; i++) {
417       Rc[i] /= num;
418     }
419     sub(d, Rc, c);
420     for (li = getAtomList().listIterator(); li.hasNext(); ) {
421       (li.next()).move(c);
422     }
423   }
424 
425   /**
426    * centerView
427    *
428    * @param rot a boolean.
429    * @param trans a boolean.
430    */
431   public void centerView(boolean rot, boolean trans) {
432     originToRot.getTransform(originToRotT3D);
433     if (rot) {
434       Matrix3d m3d = new Matrix3d();
435       m3d.setIdentity();
436       originToRotT3D.setRotation(m3d);
437       // rotToCOMT3D.setRotation(m3d);
438     }
439     if (trans) {
440       originToRotV3D.set(offset);
441       originToRotT3D.set(originToRotV3D);
442     }
443     originToRot.setTransform(originToRotT3D);
444     // rotToCOM.setTransform(rotToCOMT3D);
445   }
446 
447   /** Compute fractional coordinates. */
448   public void computeFractionalCoordinates() {
449 
450     // Count up the number of fractional coordinate entities.
451     fractionalCount();
452 
453     Crystal unitCell = getCrystal().getUnitCell();
454     double[] com = new double[3];
455 
456     switch (fractionalMode) {
457       case MOLECULE:
458         int iMolecule = 0;
459         Polymer[] polymers = getChains();
460         if (polymers != null && polymers.length > 0) {
461           // Find the center of mass
462           for (Polymer polymer : polymers) {
463             List<Atom> list = polymer.getAtomList();
464             com[0] = 0.0;
465             com[1] = 0.0;
466             com[2] = 0.0;
467             double totalMass = 0.0;
468             for (Atom atom : list) {
469               double m = atom.getMass();
470               com[0] += atom.getX() * m;
471               com[1] += atom.getY() * m;
472               com[2] += atom.getZ() * m;
473               totalMass += m;
474             }
475             com[0] /= totalMass;
476             com[1] /= totalMass;
477             com[2] /= totalMass;
478             unitCell.toFractionalCoordinates(com, fractionalCoordinates[iMolecule++]);
479           }
480         }
481 
482         // Loop over each molecule
483         List<MSNode> molecules = getMolecules();
484         for (MSNode molecule : molecules) {
485           List<Atom> list = molecule.getAtomList();
486           // Find the center of mass
487           com[0] = 0.0;
488           com[1] = 0.0;
489           com[2] = 0.0;
490           double totalMass = 0.0;
491           for (Atom atom : list) {
492             double m = atom.getMass();
493             com[0] += atom.getX() * m;
494             com[1] += atom.getY() * m;
495             com[2] += atom.getZ() * m;
496             totalMass += m;
497           }
498           com[0] /= totalMass;
499           com[1] /= totalMass;
500           com[2] /= totalMass;
501           unitCell.toFractionalCoordinates(com, fractionalCoordinates[iMolecule++]);
502         }
503 
504         // Loop over each water
505         List<MSNode> water = getWater();
506         for (MSNode wat : water) {
507           List<Atom> list = wat.getAtomList();
508           // Find the center of mass
509           com[0] = 0.0;
510           com[1] = 0.0;
511           com[2] = 0.0;
512           double totalMass = 0.0;
513           for (Atom atom : list) {
514             double m = atom.getMass();
515             com[0] += atom.getX() * m;
516             com[1] += atom.getY() * m;
517             com[2] += atom.getZ() * m;
518             totalMass += m;
519           }
520           com[0] /= totalMass;
521           com[1] /= totalMass;
522           com[2] /= totalMass;
523           unitCell.toFractionalCoordinates(com, fractionalCoordinates[iMolecule++]);
524         }
525 
526         // Loop over each ion
527         List<MSNode> ions = getIons();
528         for (MSNode ion : ions) {
529           List<Atom> list = ion.getAtomList();
530           // Find the center of mass
531           com[0] = 0.0;
532           com[1] = 0.0;
533           com[2] = 0.0;
534           double totalMass = 0.0;
535           for (Atom atom : list) {
536             double m = atom.getMass();
537             com[0] += atom.getX() * m;
538             com[1] += atom.getY() * m;
539             com[2] += atom.getZ() * m;
540             totalMass += m;
541           }
542           com[0] /= totalMass;
543           com[1] /= totalMass;
544           com[2] /= totalMass;
545           unitCell.toFractionalCoordinates(com, fractionalCoordinates[iMolecule++]);
546         }
547         break;
548       case ATOM:
549         Atom[] atoms = getAtomArray();
550         int nAtoms = atoms.length;
551         for (int i = 0; i < nAtoms; i++) {
552           atoms[i].getXYZ(com);
553           unitCell.toFractionalCoordinates(com, fractionalCoordinates[i]);
554         }
555         break;
556       case OFF:
557         break;
558     }
559   }
560 
561   /** createBox */
562   public void createBox() {
563     int vertices = 8;
564     LineArray la =
565         new LineArray(
566             4 * vertices,
567             GeometryArray.COORDINATES | GeometryArray.COLOR_4 | GeometryArray.NORMALS);
568     la.setCapability(LineArray.ALLOW_COORDINATE_WRITE);
569     la.setCapability(LineArray.ALLOW_COORDINATE_READ);
570     la.setCapability(LineArray.ALLOW_COLOR_WRITE);
571     la.setCapability(LineArray.ALLOW_COUNT_READ);
572     la.setCapability(LineArray.ALLOW_INTERSECT);
573     la.setCapability(LineArray.ALLOW_FORMAT_READ);
574 
575     ColoringAttributes cola =
576         new ColoringAttributes(new Color3f(), ColoringAttributes.SHADE_GOURAUD);
577     Appearance app = new Appearance();
578     lineAttributes = new LineAttributes();
579     lineAttributes.setLineWidth(RendererCache.bondwidth);
580     lineAttributes.setCapability(LineAttributes.ALLOW_WIDTH_WRITE);
581     lineAttributes.setLineAntialiasingEnable(true);
582     app.setLineAttributes(lineAttributes);
583     app.setCapability(Appearance.ALLOW_LINE_ATTRIBUTES_READ);
584     app.setCapability(Appearance.ALLOW_LINE_ATTRIBUTES_WRITE);
585     RenderingAttributes ra = new RenderingAttributes();
586     ra.setAlphaTestValue(0.1f);
587     ra.setAlphaTestFunction(RenderingAttributes.GREATER);
588     ra.setDepthBufferEnable(true);
589     ra.setDepthBufferWriteEnable(true);
590     app.setRenderingAttributes(ra);
591     app.setColoringAttributes(cola);
592     Shape3D wireframe = new Shape3D(la, app);
593     // PickTool.setCapabilities(wire, PickTool.INTERSECT_COORD);
594     wireframe.setUserData(this);
595     wireframe.setBounds(new BoundingSphere(new Point3d(0, 0, 0), 10.0));
596     try {
597       wireframe.setBoundsAutoCompute(false);
598     } catch (Exception e) {
599       logger.warning(" Exception in setting bounds for wireframe:\n" + e);
600     }
601     wireframe.setCapability(Shape3D.ALLOW_GEOMETRY_READ);
602     wireframe.setCapability(Shape3D.ALLOW_APPEARANCE_READ);
603     // return wire;
604   }
605 
606   /**
607    * deleteMolecule
608    *
609    * @param molecule a {@link ffx.potential.bonded.Molecule} object.
610    */
611   public void deleteMolecule(Molecule molecule) {
612     List<MSNode> list = ions.getChildList();
613     for (MSNode node : list) {
614       Molecule m = (Molecule) node;
615       if (molecule == m) {
616         ions.remove(m);
617         ionHashMap.remove(m.getKey());
618         return;
619       }
620     }
621     list = water.getChildList();
622     for (MSNode node : list) {
623       Molecule m = (Molecule) node;
624       if (molecule == m) {
625         water.remove(m);
626         waterHashMap.remove(m.getKey());
627         return;
628       }
629     }
630     list = molecules.getChildList();
631     for (MSNode node : list) {
632       Molecule m = (Molecule) node;
633       if (molecule == m) {
634         molecules.remove(m);
635         moleculeHashMap.remove(m.getKey());
636         return;
637       }
638     }
639   }
640 
641   /** {@inheritDoc} */
642   @Override
643   public boolean destroy() {
644     try {
645       if (potentialEnergy == null) {
646         finishDestruction();
647         return true;
648       } else {
649         // potentialEnergy.destroy() will loop around to call finishDestruction().
650         // This is a poor construction, but is necessary due to prior decisions.
651         return potentialEnergy.destroy();
652       }
653     } catch (Exception ex) {
654       logger.warning(String.format(" Exception in destroying a MolecularAssembly: %s", ex));
655       logger.info(Utilities.stackTraceToString(ex));
656       return false;
657     }
658   }
659 
660   /** detach */
661   public void detach() {
662     synchronized (this) {
663       if (branchGroup != null && branchGroup.isLive()) {
664         branchGroup.detach();
665       }
666     }
667   }
668 
669   /** {@inheritDoc} */
670   @Override
671   public void finalize(boolean finalizeGroups, ForceField forceField) {
672     setFinalized(false);
673     if (finalizeGroups) {
674       bondTime = 0;
675       angleTime = 0;
676       stretchBendTime = 0;
677       ureyBradleyTime = 0;
678       outOfPlaneBendTime = 0;
679       torsionTime = 0;
680       piOrbitalTorsionTime = 0;
681       torsionTorsionTime = 0;
682       List<MSNode> Polymers = getAtomNodeList();
683       for (MSNode msNode : Polymers) {
684         MSGroup group = (MSGroup) msNode;
685         if (logger.isLoggable(Level.FINE)) {
686           logger.fine(" Finalizing bonded terms for polymer " + group.toString());
687         }
688         try {
689           group.finalize(true, forceField);
690         } catch (Exception e) {
691           String message = "Fatal exception finalizing " + group.toString();
692           logger.log(Level.SEVERE, message, e);
693           System.exit(-1);
694         }
695         if (logger.isLoggable(Level.FINE)) {
696           Runtime runtime = Runtime.getRuntime();
697           long occupiedMemory = runtime.totalMemory() - runtime.freeMemory();
698           long MB = 1024 * 1024;
699           logger.fine(
700               "\n In-Use Memory   (Mb): "
701                   + occupiedMemory / MB
702                   + "\n Free Memory     (Mb): "
703                   + runtime.freeMemory() / MB
704                   + "\n Total Memory    (Mb): "
705                   + runtime.totalMemory() / MB);
706         }
707       }
708       for (MSNode m : molecules.getChildList()) {
709         Molecule molecule = (Molecule) m;
710         molecule.finalize(true, forceField);
711       }
712       for (MSNode m : water.getChildList()) {
713         Molecule molecule = (Molecule) m;
714         molecule.finalize(true, forceField);
715       }
716       for (MSNode m : ions.getChildList()) {
717         Molecule molecule = (Molecule) m;
718         molecule.finalize(true, forceField);
719       }
720       if (logger.isLoggable(Level.FINE)) {
721         StringBuilder sb = new StringBuilder("\n Time to create bonded energy terms\n\n");
722         sb.append(format(" Bond Streching     %10.3f\n", bondTime * 1.0e-9));
723         sb.append(format(" Angle Bending      %10.3f\n", angleTime * 1.0e-9));
724         sb.append(format(" Stretch-Bend       %10.3f\n", stretchBendTime * 1.0e-9));
725         sb.append(format(" Urey-Bradley       %10.3f\n", ureyBradleyTime * 1.0e-9));
726         sb.append(format(" Out-of-Plane Bend  %10.3f\n", outOfPlaneBendTime * 1.0e-9));
727         sb.append(format(" Torsionanl Angle   %10.3f\n", torsionTime * 1.0e-9));
728         sb.append(format(" Pi-Orbital Torsion %10.3f\n", piOrbitalTorsionTime * 1.0e-9));
729         sb.append(format(" Torsion-Torsion    %10.3f\n", torsionTorsionTime * 1.0e-9));
730         logger.fine(sb.toString());
731       }
732     }
733     if (!GraphicsEnvironment.isHeadless()) {
734       createScene(!finalizeGroups);
735       center();
736     }
737     removeLeaves();
738 
739     // Apply the Heavy Hydrogen flag.
740     boolean heavyHydrogen = forceField.getBoolean("HEAVY_HYDROGENS", false);
741     heavyHydrogen = forceField.getBoolean("HEAVY_HYDROGEN", heavyHydrogen);
742     if (heavyHydrogen) {
743       applyHeavyHydrogen();
744     }
745 
746     setFinalized(true);
747   }
748 
749   /**
750    * findAtom
751    *
752    * @param atom a {@link ffx.potential.bonded.Atom} object.
753    * @return a {@link ffx.potential.bonded.Atom} object.
754    */
755   public Atom findAtom(Atom atom) {
756     if (!atom.isHetero() || atom.isModRes()) {
757       Polymer polymer = getPolymer(atom.getChainID(), atom.getSegID(), false);
758       if (polymer != null) {
759         Residue res = polymer.getResidue(atom.getResidueName(), atom.getResidueNumber(), false);
760         if (res != null) {
761           MSNode node = res.getAtomNode();
762           return (Atom) node.contains(atom);
763         }
764       }
765       return null;
766     } else {
767       String resName = atom.getResidueName();
768       int resNum = atom.getResidueNumber();
769       String segID = atom.getSegID();
770       String key = resNum + resName + segID;
771       Molecule m = ionHashMap.get(key);
772       if (m == null) {
773         m = waterHashMap.get(key);
774         if (m == null) {
775           m = moleculeHashMap.get(key);
776           if (m == null) {
777             return null;
778           }
779         }
780       }
781       return (Atom) m.contains(atom);
782     }
783   }
784 
785   /**
786    * Count the number of fractional coordinate entities in the system. If fractionalMode is MOLECULE,
787    * then the count is equal to the number of molecules. If fractionalMode is ATOM, then the count is
788    * the number of atoms. Otherwise the count is zero.
789    *
790    * @return The number of fractional coordinate entities.
791    */
792   public int fractionalCount() {
793     int count = 0;
794     switch (fractionalMode) {
795       case MOLECULE:
796         // Move polymers togethers.
797         Polymer[] polymers = getChains();
798         if (polymers != null && polymers.length > 0) {
799           count += polymers.length;
800         }
801         List<MSNode> molecules = getMolecules();
802         if (molecules != null) {
803           count += molecules.size();
804         }
805 
806         List<MSNode> water = getWater();
807         if (water != null) {
808           count += water.size();
809         }
810 
811         List<MSNode> ions = getIons();
812         if (ions != null) {
813           count += ions.size();
814         }
815         break;
816       case ATOM:
817         count = getAtomArray().length;
818         break;
819       case OFF:
820         count = 0;
821         break;
822     }
823 
824     if (fractionalCoordinates == null || fractionalCoordinates.length != count) {
825       fractionalCoordinates = new double[count][3];
826     }
827 
828     return count;
829   }
830 
831   /**
832    * getActiveAtomArray
833    *
834    * @return an array of active {@link ffx.potential.bonded.Atom} objects.
835    */
836   public Atom[] getActiveAtomArray() {
837     List<Atom> atoms = getAtomList();
838     List<Atom> activeAtoms = new ArrayList<>();
839     for (Atom a : atoms) {
840       if (a.isActive()) {
841         activeAtoms.add(a);
842       }
843     }
844     Atom[] atomArray = activeAtoms.toArray(new Atom[0]);
845     Arrays.sort(atomArray);
846     return atomArray;
847   }
848 
849   /**
850    * Gets all bonded entities in this MolecularAssembly, where an entity can be a polymer, molecule,
851    * monoatomic ion, or monoatomic gas (i.e. noble gas atoms).
852    *
853    * @return All bonded groups of atoms and all singleton atoms.
854    */
855   public List<MSNode> getAllBondedEntities() {
856     // Initial construction via HashSet to eliminate duplicates.
857     Set<MSNode> allBondedNodes = new HashSet<>(getIons());
858     allBondedNodes.addAll(getMolecules());
859     allBondedNodes.addAll(getWater());
860     Polymer[] polys = getChains();
861     if (polys != null && polys.length > 0) {
862       allBondedNodes.addAll(Arrays.asList(polys));
863     }
864     return new ArrayList<>(allBondedNodes);
865   }
866 
867   /**
868    * Return an Array of all atoms in the System.
869    *
870    * @return an array of {@link ffx.potential.bonded.Atom} objects.
871    */
872   public Atom[] getAtomArray() {
873     List<Atom> atoms = getAtomList();
874     Atom[] atomArray = atoms.toArray(new Atom[0]);
875     for (int i = 0; i < atoms.size(); i++) {
876       atomArray[i].setXyzIndex(i + 1);
877     }
878 
879     return atomArray;
880   }
881 
882   /**
883    * getAtomFromWireVertex
884    *
885    * @param i a int.
886    * @return a {@link ffx.potential.bonded.Atom} object.
887    */
888   public Atom getAtomFromWireVertex(int i) {
889     if (atomLookUp != null && atomLookUp.length > i) {
890       return atomLookUp[i];
891     }
892     return null;
893   }
894 
895   /**
896    * getBackBoneAtoms
897    *
898    * @return a {@link java.util.List} object.
899    */
900   public List<Atom> getBackBoneAtoms() {
901     List<Atom> backbone = new ArrayList<>();
902     List<Residue> residues = getResidueList();
903     for (Residue residue : residues) {
904       backbone.addAll(residue.getBackboneAtoms());
905     }
906     return backbone;
907   }
908 
909   /**
910    * Getter for the field <code>branchGroup</code>.
911    *
912    * @return a {@link org.jogamp.java3d.BranchGroup} object.
913    */
914   public BranchGroup getBranchGroup() {
915     return branchGroup;
916   }
917 
918   /**
919    * getChain
920    *
921    * @param name a {@link java.lang.String} object.
922    * @return a {@link ffx.potential.bonded.Polymer} object.
923    */
924   public Polymer getChain(String name) {
925     for (MSNode node : getAtomNodeList()) {
926       if (node instanceof Polymer) {
927         String chainName = node.getName();
928         if (chainName.equalsIgnoreCase(name)) {
929           return (Polymer) node;
930         } else if (name.contentEquals(" ")) {
931           return (Polymer) node;
932         }
933         /* TODO: Right now if a molecular assembly has no chain ID, the first node is returned.
934         This functionality will not work for cases where multiple chains exist and one of those chains has no name.*/
935       }
936     }
937     return null;
938   }
939 
940   /**
941    * getChainNames
942    *
943    * @return an array of {@link java.lang.String} objects.
944    */
945   public String[] getChainNames() {
946     List<String> temp = new ArrayList<>();
947     for (MSNode node : getAtomNodeList()) {
948       if (node instanceof Polymer) {
949         temp.add(node.getName());
950       }
951     }
952     if (temp.isEmpty()) {
953       return null;
954     }
955 
956     String[] names = new String[temp.size()];
957     for (int i = 0; i < temp.size(); i++) {
958       names[i] = temp.get(i);
959     }
960 
961     return names;
962   }
963 
964   /**
965    * getChains
966    *
967    * @return an array of {@link ffx.potential.bonded.Polymer} objects.
968    */
969   public Polymer[] getChains() {
970     List<Polymer> polymers = new ArrayList<>();
971     for (MSNode node : getAtomNodeList()) {
972       if (node instanceof Polymer) {
973         polymers.add((Polymer) node);
974       }
975     }
976     if (polymers.isEmpty()) {
977       return null;
978     }
979     return polymers.toArray(new Polymer[0]);
980   }
981 
982   /**
983    * Sums up charge of the system, checking nonstandard residues for non-unitary charges.
984    *
985    * @param alwaysLog Log non-unitary charge warnings for all nodes
986    * @return System charge
987    */
988   public double getCharge(boolean alwaysLog) {
989     double totalCharge = 0;
990     for (MSNode node : getNodeList()) {
991       double charge = 0;
992       boolean isNonstandard = false;
993       for (Atom atom : node.getAtomList()) {
994         charge += atom.getCharge(forceField);
995         if (atom.isModRes()) {
996           isNonstandard = true;
997         }
998       }
999       if ((alwaysLog || isNonstandard) && (Math.abs(Math.round(charge) - charge) > 1.0E-5)) {
1000         logger.warning(
1001             String.format(" Node %s has non-unitary charge %12.8f", node, charge));
1002       }
1003       totalCharge += charge;
1004     }
1005     return totalCharge;
1006   }
1007 
1008   /**
1009    * getCrystal
1010    *
1011    * @return a {@link ffx.crystal.Crystal} object.
1012    */
1013   public Crystal getCrystal() {
1014     if (potentialEnergy == null) {
1015       return null;
1016     }
1017     return potentialEnergy.getCrystal();
1018   }
1019 
1020   /**
1021    * Set the Crystal for the Potential of this MolecularAssembly.
1022    *
1023    * @param crystal Crystal instance.
1024    */
1025   public void setCrystal(Crystal crystal) {
1026     if (potentialEnergy != null) {
1027       potentialEnergy.setCrystal(crystal);
1028     }
1029   }
1030 
1031   /**
1032    * Getter for the field <code>currentCycle</code>.
1033    *
1034    * @return a int.
1035    */
1036   public int getCurrentCycle() {
1037     return currentCycle;
1038   }
1039 
1040   /**
1041    * Setter for the field <code>currentCycle</code>.
1042    *
1043    * @param c a int.
1044    */
1045   public void setCurrentCycle(int c) {
1046     if (c <= cycles && c > 0) {
1047       currentCycle = c;
1048       for (Atom atom : getAtomList()) {
1049         atom.setCurrentCycle(currentCycle);
1050       }
1051     }
1052   }
1053 
1054   /**
1055    * Getter for the field <code>cycles</code>.
1056    *
1057    * @return a int.
1058    */
1059   public int getCycles() {
1060     return cycles;
1061   }
1062 
1063   /**
1064    * Setter for the field <code>cycles</code>.
1065    *
1066    * @param c a int.
1067    */
1068   public void setCycles(int c) {
1069     cycles = c;
1070   }
1071 
1072   /** {@inheritDoc} */
1073   @Override
1074   public double getExtent() {
1075     double[] Rc = {0, 0, 0};
1076     int num = getAtomList().size();
1077     for (Atom atom : getAtomList()) {
1078       atom.getXYZ(a);
1079       Rc[0] += a[0];
1080       Rc[1] += a[1];
1081       Rc[2] += a[2];
1082     }
1083 
1084     for (int i = 0; i < 3; i++) {
1085       Rc[i] /= num;
1086     }
1087 
1088     double r, d = 0;
1089     double[] xyz = new double[3];
1090     for (Atom atom : getAtomList()) {
1091       atom.getXYZ(xyz);
1092       sub(xyz, Rc, xyz);
1093       r = length(xyz);
1094       if (d < r) {
1095         d = r;
1096       }
1097     }
1098     return d;
1099   }
1100 
1101   /**
1102    * Getter for the field <code>file</code>.
1103    *
1104    * @return a {@link java.io.File} object.
1105    */
1106   public File getFile() {
1107     return file;
1108   }
1109 
1110   /**
1111    * Setter for the field <code>file</code>.
1112    *
1113    * @param file a {@link java.io.File} object.
1114    */
1115   public void setFile(File file) {
1116     if (file == null) {
1117       // Keep the current file.
1118       return;
1119     }
1120     this.file = file;
1121   }
1122 
1123   /**
1124    * Getter for the field <code>archiveFile</code>.
1125    *
1126    * @return a {@link java.io.File} object.
1127    */
1128   public File getArchiveFile() {
1129     return archiveFile;
1130   }
1131 
1132   /**
1133    * Set the File for writing out an archive.
1134    *
1135    * @param archiveFile The archive file.
1136    */
1137   public void setArchiveFile(File archiveFile) {
1138     if (archiveFile == null) {
1139       // Keep the current archive file.
1140       return;
1141     }
1142     this.archiveFile = archiveFile;
1143   }
1144 
1145   /**
1146    * Getter for the field <code>forceField</code>.
1147    *
1148    * @return a {@link ffx.potential.parameters.ForceField} object.
1149    */
1150   public ForceField getForceField() {
1151     return forceField;
1152   }
1153 
1154   /**
1155    * Setter for the field <code>forceField</code>.
1156    *
1157    * @param forceField a {@link ffx.potential.parameters.ForceField} object.
1158    */
1159   public void setForceField(ForceField forceField) {
1160     this.forceField = forceField;
1161   }
1162 
1163   /**
1164    * Getter for the field <code>fractionalMode</code>.
1165    *
1166    * @return a {@link ffx.potential.MolecularAssembly.FractionalMode} object.
1167    */
1168   public FractionalMode getFractionalMode() {
1169     return fractionalMode;
1170   }
1171 
1172   /**
1173    * Setter for the field <code>fractionalMode</code>.
1174    *
1175    * @param mode a {@link ffx.potential.MolecularAssembly.FractionalMode} object.
1176    */
1177   public void setFractionalMode(FractionalMode mode) {
1178     fractionalMode = mode;
1179   }
1180 
1181   /**
1182    * Gets the header lines associated with this MolecularAssembly (particularly for PDB)
1183    *
1184    * @return Header lines.
1185    */
1186   public String[] getHeaderLines() {
1187     String[] ret = new String[headerLines.size()];
1188     headerLines.toArray(ret);
1189     return ret;
1190   }
1191 
1192   /**
1193    * Getter for the field <code>ions</code>.
1194    *
1195    * @return a {@link java.util.List} object.
1196    */
1197   public List<MSNode> getIons() {
1198     return ions.getChildList();
1199   }
1200 
1201   /**
1202    * getMass.
1203    *
1204    * @return a double.
1205    */
1206   public double getMass() {
1207     double mass = 0;
1208     for (Atom atom : getAtomArray()) {
1209       mass += atom.getMass();
1210     }
1211     return mass;
1212   }
1213 
1214   /**
1215    * Construct a boolean orray of flags to indicate atoms treated by a neural network.
1216    *
1217    * @return The boolean orray of flags to indicate atoms treated by a neural network.
1218    */
1219   public boolean[] getNeuralNetworkIdentity() {
1220     Atom[] atomArray = getAtomArray();
1221     int nAtoms = atomArray.length;
1222     boolean[] neuralNetwork = new boolean[nAtoms];
1223     for (int i = 0; i < nAtoms; i++) {
1224       neuralNetwork[i] = atomArray[i].isNeuralNetwork();
1225     }
1226     return neuralNetwork;
1227   }
1228 
1229   /**
1230    * This method assigns a unique integer to every molecule in the MolecularAssembly beginning at 0.
1231    * An integer array with these values for each atom is returned.
1232    *
1233    * @return an array of molecule numbers for each atom.
1234    */
1235   public int[] getMoleculeNumbers() {
1236     int[] moleculeNumber = new int[getAtomList().size()];
1237     int current = 0;
1238     // Loop over polymers together
1239     Polymer[] polymers = getChains();
1240     if (polymers != null) {
1241       for (Polymer polymer : polymers) {
1242         List<Atom> atomList = polymer.getAtomList();
1243         for (Atom atom : atomList) {
1244           moleculeNumber[atom.getXyzIndex() - 1] = current;
1245           atom.setMoleculeNumber(current);
1246         }
1247         current++;
1248       }
1249     }
1250 
1251     // Loop over each molecule
1252     for (MSNode molecule : molecules.getChildList()) {
1253       List<Atom> atomList = molecule.getAtomList();
1254       for (Atom atom : atomList) {
1255         moleculeNumber[atom.getXyzIndex() - 1] = current;
1256         atom.setMoleculeNumber(current);
1257       }
1258       current++;
1259     }
1260 
1261     // Loop over each water
1262     for (MSNode wat : water.getChildList()) {
1263       List<Atom> atomList = wat.getAtomList();
1264       for (Atom atom : atomList) {
1265         moleculeNumber[atom.getXyzIndex() - 1] = current;
1266         atom.setMoleculeNumber(current);
1267       }
1268       current++;
1269     }
1270 
1271     // Loop over each ion
1272     for (MSNode ion : ions.getChildList()) {
1273       List<Atom> atomList = ion.getAtomList();
1274       for (Atom atom : atomList) {
1275         moleculeNumber[atom.getXyzIndex() - 1] = current;
1276         atom.setMoleculeNumber(current);
1277       }
1278       current++;
1279     }
1280 
1281     return moleculeNumber;
1282   }
1283 
1284   /**
1285    * Getter for the field <code>molecules</code>.
1286    *
1287    * @return a {@link java.util.List} object.
1288    */
1289   public List<MSNode> getMolecules() {
1290     return molecules.getChildList();
1291   }
1292 
1293   public Molecule[] getMoleculeArray() {
1294     Molecule[] m = new Molecule[molecules.getChildCount()];
1295     int i = 0;
1296     for (MSNode msNode : molecules.getChildList()) {
1297       m[i++] = (Molecule) msNode;
1298     }
1299     return m;
1300   }
1301 
1302   /**
1303    * This method sets all HETATM molecules, including water and ions, to use the given chainID and
1304    * then renumbers the molecules.
1305    * <p>
1306    * If there is a polymer with the given chainID, then numbering begins after the final ResID.
1307    * Otherwise numbering begins at 1.
1308    * <p>
1309    * This is useful for producing consistent PDB file output.
1310    *
1311    * @param chainID The Character chainID to use for all HETATM molecules.
1312    */
1313   public void setChainIDAndRenumberMolecules(Character chainID) {
1314     int resID = 1;
1315 
1316     Polymer polymer = getPolymer(chainID, chainID.toString(), false);
1317     if (polymer != null) {
1318       List<Residue> residues = polymer.getResidues();
1319       for (Residue residue : residues) {
1320         int resID2 = residue.getResidueNumber();
1321         if (resID2 >= resID) {
1322           resID = resID2 + 1;
1323         }
1324       }
1325     }
1326     for (MSNode m : getMolecules()) {
1327       Molecule molecule = (Molecule) m;
1328       molecule.setChainID(chainID);
1329       molecule.setResidueNum(resID++);
1330     }
1331     for (MSNode ion : getIons()) {
1332       Molecule m = (Molecule) ion;
1333       m.setChainID(chainID);
1334       m.setResidueNum(resID++);
1335     }
1336     for (MSNode wat : getWater()) {
1337       Molecule water = (Molecule) wat;
1338       water.setChainID(chainID);
1339       water.setResidueNum(resID++);
1340     }
1341   }
1342 
1343   /**
1344    * getNodeList
1345    *
1346    * @return a {@link java.util.List} object.
1347    */
1348   public List<MSNode> getNodeList() {
1349     List<MSNode> residues = new ArrayList<>();
1350     ListIterator<MSNode> li, lj;
1351     MSNode o;
1352     Polymer c;
1353     for (li = getAtomNodeList().listIterator(); li.hasNext(); ) {
1354       o = li.next();
1355       if (o instanceof Polymer) {
1356         c = (Polymer) o;
1357         for (lj = c.getAtomNodeList().listIterator(); lj.hasNext(); ) {
1358           o = lj.next();
1359           if (o instanceof Residue) {
1360             residues.add(o);
1361           }
1362         }
1363       }
1364     }
1365 
1366     residues.addAll(ions.getChildList());
1367     residues.addAll(water.getChildList());
1368     residues.addAll(molecules.getChildList());
1369 
1370     return residues;
1371   }
1372 
1373   /**
1374    * Getter for the field <code>offset</code>.
1375    *
1376    * @return a {@link org.jogamp.vecmath.Vector3d} object.
1377    */
1378   public Vector3d getOffset() {
1379     if (offset == null) {
1380       offset = new Vector3d(0.0, 0.0, 0.0);
1381     }
1382 
1383     return offset;
1384   }
1385 
1386   /**
1387    * Setter for the field <code>offset</code>.
1388    *
1389    * @param o a {@link org.jogamp.vecmath.Vector3d} object.
1390    */
1391   public void setOffset(Vector3d o) {
1392     offset = o;
1393   }
1394 
1395   /**
1396    * Getter for the field <code>originToRot</code>.
1397    *
1398    * @return a {@link org.jogamp.java3d.TransformGroup} object.
1399    */
1400   public TransformGroup getOriginToRot() {
1401     return originToRot;
1402   }
1403 
1404   /**
1405    * getParallelTeam.
1406    *
1407    * @return a {@link edu.rit.pj.ParallelTeam} object.
1408    */
1409   public ParallelTeam getParallelTeam() {
1410     if (potentialEnergy != null) {
1411       return potentialEnergy.getParallelTeam();
1412     } else {
1413       return null;
1414     }
1415   }
1416 
1417   /**
1418    * getPolymer
1419    *
1420    * @param chainID a {@link java.lang.Character} object.
1421    * @param segID a {@link java.lang.String} object.
1422    * @param create a boolean.
1423    * @return a {@link ffx.potential.bonded.Polymer} object.
1424    */
1425   public Polymer getPolymer(Character chainID, String segID, boolean create) {
1426     for (MSNode node : getAtomNodeList()) {
1427       if (node instanceof Polymer) {
1428         Polymer polymer = (Polymer) node;
1429         if (polymer.getName().equals(segID) && polymer.getChainID().equals(chainID)) {
1430           return (Polymer) node;
1431         }
1432       }
1433     }
1434     if (create) {
1435       Polymer polymer = new Polymer(chainID, segID, true);
1436       addMSNode(polymer);
1437       return polymer;
1438     }
1439 
1440     return null;
1441   }
1442 
1443   /**
1444    * Getter for the field <code>potentialEnergy</code>.
1445    *
1446    * @return a {@link ffx.potential.ForceFieldEnergy} object.
1447    */
1448   public ForceFieldEnergy getPotentialEnergy() {
1449     return potentialEnergy;
1450   }
1451 
1452   /**
1453    * Getter for the field <code>properties</code>.
1454    *
1455    * @return a {@link org.apache.commons.configuration2.CompositeConfiguration} object.
1456    */
1457   public CompositeConfiguration getProperties() {
1458     return properties == null ? forceField.getProperties() : properties;
1459   }
1460 
1461   /**
1462    * getResidueList
1463    *
1464    * @return a {@link java.util.List} object.
1465    */
1466   public List<Residue> getResidueList() {
1467     List<Residue> residues = new ArrayList<>();
1468     ListIterator<MSNode> li, lj;
1469     MSNode o;
1470     Polymer c;
1471     for (li = getAtomNodeList().listIterator(); li.hasNext(); ) {
1472       o = li.next();
1473       if (o instanceof Polymer) {
1474         c = (Polymer) o;
1475         for (lj = c.getAtomNodeList().listIterator(); lj.hasNext(); ) {
1476           o = lj.next();
1477           if (o instanceof Residue) {
1478             residues.add((Residue) o);
1479           }
1480         }
1481       }
1482     }
1483     return residues;
1484   }
1485 
1486   /**
1487    * getTransformGroup
1488    *
1489    * @return a {@link org.jogamp.java3d.TransformGroup} object.
1490    */
1491   public TransformGroup getTransformGroup() {
1492     return originToRot;
1493   }
1494 
1495   /**
1496    * getWater
1497    *
1498    * @return a {@link java.util.List} object.
1499    */
1500   public List<MSNode> getWater() {
1501     return water.getChildList();
1502   }
1503 
1504   /**
1505    * getWireFrame
1506    *
1507    * @return a {@link org.jogamp.java3d.Node} object.
1508    */
1509   public Node getWireFrame() {
1510     return wire;
1511   }
1512 
1513   /**
1514    * isVisible
1515    *
1516    * @return a boolean.
1517    */
1518   public boolean isVisible() {
1519     return visible;
1520   }
1521 
1522   /**
1523    * loadVRML
1524    *
1525    * @return a {@link org.jogamp.java3d.BranchGroup} object.
1526    */
1527   public BranchGroup loadVRML() {
1528     //        try {
1529     //            VrmlLoader loader = new VrmlLoader();
1530     //            VrmlScene scene = null;
1531     //            if (vrmlFile != null && vrmlFile.exists()) {
1532     //                scene = (VrmlScene) loader.load(vrmlFile.getAbsolutePath());
1533     //            } else if (vrmlURL != null) {
1534     //                scene = (VrmlScene) loader.load(vrmlURL);
1535     //            } else {
1536     //                return null;
1537     //            }
1538     //            BranchGroup bg = scene.getSceneGroup();
1539     //            recurseVRML(bg);
1540     //            bg.setCapability(BranchGroup.ALLOW_DETACH);
1541     //            bg.setCapability(BranchGroup.ALLOW_BOUNDS_READ);
1542     //            bg.compile();
1543     //            return bg;
1544     //        } catch (Exception e) {
1545     //            String message = "Fatal exception loading VRML.\n";
1546     //            logger.log(Level.SEVERE, message, e);
1547     //            System.exit(-1);
1548     //            return null;
1549     //        }
1550     return null;
1551   }
1552 
1553   /**
1554    * Moves the center of all chemical entities into the primary unit cell. Somewhat experimental
1555    * feature; use with caution.
1556    */
1557   public void moveAllIntoUnitCell() {
1558     moveIntoUnitCell(getChains());
1559     moveIntoUnitCell(getWater());
1560     moveIntoUnitCell(getIons());
1561     moveIntoUnitCell(getMolecules());
1562   }
1563 
1564   /**
1565    * moveCenter
1566    *
1567    * @param d an array of double.
1568    */
1569   public void moveCenter(double[] d) {
1570     for (Atom atom : getAtomList()) {
1571       atom.move(d);
1572     }
1573   }
1574 
1575   /** Move to fractional coordinates. */
1576   public void moveToFractionalCoordinates() {
1577 
1578     if (fractionalCoordinates == null) {
1579       return;
1580     }
1581 
1582     Crystal unitCell = getCrystal().getUnitCell();
1583     double[] com = new double[3];
1584 
1585     switch (fractionalMode) {
1586       case MOLECULE:
1587         int iMolecule = 0;
1588         Polymer[] polymers = getChains();
1589         if (polymers != null && polymers.length > 0) {
1590           // Find the center of mass
1591           for (Polymer polymer : polymers) {
1592             List<Atom> list = polymer.getAtomList();
1593             double totalMass = 0.9;
1594             com[0] = 0.0;
1595             com[1] = 0.0;
1596             com[2] = 0.0;
1597             for (Atom atom : list) {
1598               double m = atom.getMass();
1599               com[0] += atom.getX() * m;
1600               com[1] += atom.getY() * m;
1601               com[2] += atom.getZ() * m;
1602               totalMass += m;
1603             }
1604             com[0] /= totalMass;
1605             com[1] /= totalMass;
1606             com[2] /= totalMass;
1607             // Find the new center of mass in fractional coordinates.
1608             unitCell.toFractionalCoordinates(com, com);
1609             // Find the reciprocal translation vector.
1610             double[] frac = fractionalCoordinates[iMolecule++];
1611             com[0] = frac[0] - com[0];
1612             com[1] = frac[1] - com[1];
1613             com[2] = frac[2] - com[2];
1614             // Convert the fractional translation vector to Cartesian coordinates.
1615             unitCell.toCartesianCoordinates(com, com);
1616             // Move all atoms.
1617             for (Atom atom : list) {
1618               atom.move(com);
1619             }
1620           }
1621         }
1622 
1623         // Loop over each molecule
1624         List<MSNode> molecules = getMolecules();
1625         for (MSNode molecule : molecules) {
1626           List<Atom> list = molecule.getAtomList();
1627           // Find the center of mass
1628           com[0] = 0.0;
1629           com[1] = 0.0;
1630           com[2] = 0.0;
1631           double totalMass = 0.0;
1632           for (Atom atom : list) {
1633             double m = atom.getMass();
1634             com[0] += atom.getX() * m;
1635             com[1] += atom.getY() * m;
1636             com[2] += atom.getZ() * m;
1637             totalMass += m;
1638           }
1639           com[0] /= totalMass;
1640           com[1] /= totalMass;
1641           com[2] /= totalMass;
1642           // Find the new center of mass in fractional coordinates.
1643           unitCell.toFractionalCoordinates(com, com);
1644           // Find the reciprocal translation vector to the previous COM.
1645           double[] frac = fractionalCoordinates[iMolecule++];
1646           com[0] = frac[0] - com[0];
1647           com[1] = frac[1] - com[1];
1648           com[2] = frac[2] - com[2];
1649           // Convert the fractional translation vector to Cartesian coordinates.
1650           unitCell.toCartesianCoordinates(com, com);
1651           // Move all atoms.
1652           for (Atom atom : list) {
1653             atom.move(com);
1654           }
1655         }
1656 
1657         // Loop over each water
1658         List<MSNode> water = getWater();
1659         for (MSNode wat : water) {
1660           List<Atom> list = wat.getAtomList();
1661           // Find the center of mass
1662           com[0] = 0.0;
1663           com[1] = 0.0;
1664           com[2] = 0.0;
1665           double totalMass = 0.0;
1666           for (Atom atom : list) {
1667             double m = atom.getMass();
1668             com[0] += atom.getX() * m;
1669             com[1] += atom.getY() * m;
1670             com[2] += atom.getZ() * m;
1671             totalMass += m;
1672           }
1673           com[0] /= totalMass;
1674           com[1] /= totalMass;
1675           com[2] /= totalMass;
1676           // Find the new center of mass in fractional coordinates.
1677           unitCell.toFractionalCoordinates(com, com);
1678           // Find the reciprocal translation vector to the previous COM.
1679           double[] frac = fractionalCoordinates[iMolecule++];
1680           com[0] = frac[0] - com[0];
1681           com[1] = frac[1] - com[1];
1682           com[2] = frac[2] - com[2];
1683           // Convert the fractional translation vector to Cartesian coordinates.
1684           unitCell.toCartesianCoordinates(com, com);
1685 
1686           double r = length(com);
1687 
1688           // Warn if an atom is moved more than 1 Angstrom.
1689           if (r > 1.0) {
1690             int i = iMolecule - 1;
1691             logger.info(String.format(" %d R: %16.8f", i, r));
1692             logger.info(
1693                 String.format(" %d FRAC %16.8f %16.8f %16.8f", i, frac[0], frac[1], frac[2]));
1694             logger.info(String.format(" %d COM  %16.8f %16.8f %16.8f", i, com[0], com[1], com[2]));
1695           }
1696 
1697           // Move all atoms.
1698           for (Atom atom : list) {
1699             atom.move(com);
1700           }
1701         }
1702 
1703         // Loop over each ion
1704         List<MSNode> ions = getIons();
1705         for (MSNode ion : ions) {
1706           List<Atom> list = ion.getAtomList();
1707           // Find the center of mass
1708           com[0] = 0.0;
1709           com[1] = 0.0;
1710           com[2] = 0.0;
1711           double totalMass = 0.0;
1712           for (Atom atom : list) {
1713             double m = atom.getMass();
1714             com[0] += atom.getX() * m;
1715             com[1] += atom.getY() * m;
1716             com[2] += atom.getZ() * m;
1717             totalMass += m;
1718           }
1719           com[0] /= totalMass;
1720           com[1] /= totalMass;
1721           com[2] /= totalMass;
1722           // Find the new center of mass in fractional coordinates.
1723           unitCell.toFractionalCoordinates(com, com);
1724           // Find the reciprocal translation vector to the previous COM.
1725           double[] frac = fractionalCoordinates[iMolecule++];
1726           com[0] = frac[0] - com[0];
1727           com[1] = frac[1] - com[1];
1728           com[2] = frac[2] - com[2];
1729           // Convert the fractional translation vector to Cartesian coordinates.
1730           unitCell.toCartesianCoordinates(com, com);
1731           // Move all atoms.
1732           for (Atom atom : list) {
1733             atom.move(com);
1734           }
1735         }
1736         break;
1737       case ATOM:
1738         Atom[] atoms = getAtomArray();
1739         int nAtoms = atoms.length;
1740         for (int i = 0; i < nAtoms; i++) {
1741           // Convert the stored factional coordinates to Cartesian coordinates in the current
1742           // unitcell.
1743           unitCell.toCartesianCoordinates(fractionalCoordinates[i], com);
1744           atoms[i].moveTo(com);
1745         }
1746         break;
1747       case OFF:
1748         break;
1749     }
1750   }
1751 
1752   /**
1753    * Rotate about a point in given in the System's Local Coordinates
1754    *
1755    * @param v Vector3d
1756    */
1757   public void rotateAbout(Vector3d v) {
1758     Vector3d newRotPoint = new Vector3d(v);
1759     originToRot.getTransform(originToRotT3D);
1760     originToRotT3D.get(originToRotV3D);
1761     originToRotT3D.setTranslation(new Vector3d(0, 0, 0));
1762     rotToCOM.getTransform(rotToCOMT3D);
1763     rotToCOMT3D.get(rotToCOMV3D);
1764     newRotPoint.add(rotToCOMV3D);
1765     originToRotT3D.transform(newRotPoint);
1766     newRotPoint.add(originToRotV3D);
1767     originToRotT3D.setTranslation(newRotPoint);
1768     rotToCOMV3D.set(v);
1769     rotToCOMV3D.negate();
1770     rotToCOMT3D.setTranslation(rotToCOMV3D);
1771     originToRot.setTransform(originToRotT3D);
1772     rotToCOM.setTransform(rotToCOMT3D);
1773   }
1774 
1775   /**
1776    * sceneGraphChange
1777    *
1778    * @param newShapes a {@link java.util.List} object.
1779    */
1780   public void sceneGraphChange(List<BranchGroup> newShapes) {
1781     if (newShapes == null) {
1782       newShapes = myNewShapes;
1783     }
1784 
1785     if (newShapes.isEmpty()) {
1786       return;
1787     }
1788 
1789     boolean reCompile = false;
1790     // Check for nodes (new and/or recycled) being added to this
1791     // MolecularAssembly
1792     for (ListIterator<BranchGroup> li = newShapes.listIterator(); li.hasNext(); ) {
1793       BranchGroup group = li.next();
1794       li.remove();
1795       // This is code for cycling between two MolecularAssemblies
1796       if (group.getUserData() != null) {
1797         logger.info(format("%s %s", group, group.getUserData().toString()));
1798         /*
1799          * Object userData = group.getUserData(); if (userData!=this) {
1800          * // The appearance has already been set during a recursive
1801          * call to setView, // although we need to turn back on
1802          * Pickablility. TransformGroup tg = (TransformGroup)
1803          * group.getChild(0); Shape3D shape = (Shape3D) tg.getChild(0);
1804          * shape.setPickable(true); group.setUserData(this); if
1805          * (!reCompile) { if (childNodes.isLive()) {
1806          * childNodes.detach(); } reCompile = true; }
1807          * childNodes.moveTo(group);
1808          */
1809       } else {
1810         // This is a new group since it has no userData.
1811         // We can not query for the identity of its parent later, so
1812         // we will store it as a userData reference.
1813         group.setUserData(this);
1814         if (!reCompile) {
1815           if (childNodes.isLive()) {
1816             childNodes.detach();
1817           }
1818 
1819           reCompile = true;
1820         }
1821 
1822         childNodes.addChild(group);
1823       }
1824     }
1825     if (reCompile) {
1826       childNodes.compile();
1827       base.addChild(childNodes);
1828     }
1829   }
1830 
1831   /** {@inheritDoc} */
1832   @Override
1833   public void setColor(RendererCache.ColorModel newColorModel, Color3f color, Material mat) {
1834     for (MSNode msNode : getAtomNodeList()) {
1835       MSGroup group = (MSGroup) msNode;
1836       group.setColor(newColorModel, color, mat);
1837     }
1838 
1839     for (MSNode m : molecules.getChildList()) {
1840       m.setColor(newColorModel, color, mat);
1841     }
1842 
1843     for (MSNode m : water.getChildList()) {
1844       m.setColor(newColorModel, color, mat);
1845     }
1846 
1847     for (MSNode m : ions.getChildList()) {
1848       m.setColor(newColorModel, color, mat);
1849     }
1850   }
1851 
1852   /**
1853    * setPotential
1854    *
1855    * @param potentialEnergy a {@link ffx.potential.ForceFieldEnergy} object.
1856    */
1857   public void setPotential(ForceFieldEnergy potentialEnergy) {
1858     this.potentialEnergy = potentialEnergy;
1859   }
1860 
1861   /** {@inheritDoc} */
1862   @Override
1863   public void setView(RendererCache.ViewModel newViewModel, List<BranchGroup> newShapes) {
1864     // Just Detach the whole system branch group
1865     if (newViewModel == RendererCache.ViewModel.DESTROY) {
1866       if (switchGroup != null) {
1867         switchGroup.setWhichChild(Switch.CHILD_NONE);
1868       }
1869       visible = false;
1870     } else if (newViewModel == RendererCache.ViewModel.SHOWVRML) {
1871       switchGroup.setWhichChild(Switch.CHILD_ALL);
1872     } else if (newViewModel == RendererCache.ViewModel.HIDEVRML) {
1873       switchGroup.setWhichChild(0);
1874     } else {
1875       setWireWidth(RendererCache.bondwidth);
1876       if (newViewModel == RendererCache.ViewModel.DETAIL && childNodes.isLive()) {
1877         childNodes.detach();
1878       }
1879       /*
1880        We'll collect new Scenegraph Shapes in our newShapeNode This is to avoid the case where
1881        setView is called from the root node and all new shapes for every MolecularAssembly would
1882        then be put into the same List.
1883        */
1884       super.setView(newViewModel, myNewShapes);
1885       List<Molecule> moleculeList = getList(Molecule.class, new ArrayList<>());
1886       for (Molecule m : moleculeList) {
1887         m.setView(newViewModel, myNewShapes);
1888       }
1889       for (MSNode m : molecules.getChildList()) {
1890         m.setView(newViewModel, myNewShapes);
1891       }
1892       for (MSNode m : water.getChildList()) {
1893         m.setView(newViewModel, myNewShapes);
1894       }
1895       for (MSNode m : ions.getChildList()) {
1896         m.setView(newViewModel, myNewShapes);
1897       }
1898       if (newViewModel == RendererCache.ViewModel.INVISIBLE) {
1899         switchGroup.setWhichChild(0);
1900       }
1901       if (newViewModel == RendererCache.ViewModel.DETAIL) {
1902         childNodes.compile();
1903         base.addChild(childNodes);
1904       }
1905     }
1906   }
1907 
1908   /**
1909    * setWireWidth
1910    *
1911    * @param f a float.
1912    */
1913   public void setWireWidth(float f) {
1914     if (wire == null) {
1915       return;
1916     }
1917 
1918     lineAttributes.setLineWidth(f);
1919   }
1920 
1921   /**
1922    * The MolecularAssembly BranchGroup has two TransformGroups between it and the "base" node where
1923    * geometry is attached. If the point between the two transformations is where user rotation
1924    * occurs. For example, if rotating about the center of mass of the system, the RotToCOM
1925    * transformation will be an identity transformation (ie. none). If rotation is about some atom or
1926    * group of atoms within the system, then the RotToCOM transformation will be a translation from
1927    * that point to the COM.
1928    *
1929    * @param zero boolean
1930    * @return BranchGroup
1931    */
1932   private BranchGroup createScene(boolean zero) {
1933     originToRotT3D = new Transform3D();
1934     originToRotV3D = new Vector3d();
1935     originToRot = new TransformGroup(originToRotT3D);
1936     branchGroup = new BranchGroup();
1937     rotToCOM = new TransformGroup();
1938     rotToCOMT3D = new Transform3D();
1939     rotToCOMV3D = new Vector3d();
1940     // Set capabilities needed for picking and moving the MolecularAssembly
1941     branchGroup.setCapability(BranchGroup.ALLOW_DETACH);
1942     originToRot.setCapability(TransformGroup.ALLOW_TRANSFORM_WRITE);
1943     originToRot.setCapability(TransformGroup.ALLOW_TRANSFORM_READ);
1944     originToRot.setCapability(TransformGroup.ENABLE_PICK_REPORTING);
1945     rotToCOM.setCapability(TransformGroup.ALLOW_TRANSFORM_READ);
1946     rotToCOM.setCapability(TransformGroup.ALLOW_TRANSFORM_WRITE);
1947     // Put the MolecularAssembly in the middle of the scene
1948     if (zero) {
1949       originToRotV3D.set(0.0, 0.0, 0.0);
1950       originToRotT3D.set(originToRotV3D);
1951       originToRot.setTransform(originToRotT3D);
1952     }
1953     wire = renderWire();
1954     switchGroup = new Switch(Switch.CHILD_NONE);
1955     switchGroup.setCapability(Switch.ALLOW_SWITCH_WRITE);
1956     base = new BranchGroup();
1957     base.setCapability(BranchGroup.ALLOW_CHILDREN_EXTEND);
1958     base.setCapability(BranchGroup.ALLOW_CHILDREN_WRITE);
1959     childNodes = new BranchGroup();
1960     childNodes.setCapability(BranchGroup.ALLOW_DETACH);
1961     childNodes.setCapability(BranchGroup.ALLOW_CHILDREN_EXTEND);
1962     childNodes.setCapability(BranchGroup.ALLOW_CHILDREN_WRITE);
1963     switchGroup.addChild(base);
1964     if (wire != null) {
1965       base.addChild(wire);
1966     }
1967     vrml = loadVRML();
1968     if (vrml != null) {
1969       vrmlTG = new TransformGroup();
1970       vrmlTd = new Transform3D();
1971       vrmlTG.setTransform(vrmlTd);
1972       vrmlTG.setCapability(TransformGroup.ALLOW_TRANSFORM_WRITE);
1973       vrmlTG.addChild(vrml);
1974       switchGroup.addChild(vrmlTG);
1975       setView(RendererCache.ViewModel.INVISIBLE, null);
1976     }
1977     switchGroup.setWhichChild(Switch.CHILD_ALL);
1978     rotToCOM.addChild(switchGroup);
1979     originToRot.addChild(rotToCOM);
1980     branchGroup.addChild(originToRot);
1981     branchGroup.compile();
1982     return branchGroup;
1983   }
1984 
1985   void finishDestruction() {
1986     detach();
1987     super.destroy();
1988   }
1989 
1990   /**
1991    * Move mass from heavy atoms to their attached hydrogens.
1992    *
1993    * <p>The mass of each hydrogen is scaled by a factor of 3. The mass of the heavy atom is reduced
1994    * by 2 AMU.
1995    */
1996   private void applyHeavyHydrogen() {
1997     List<Bond> bonds = getBondList();
1998     for (Bond bond : bonds) {
1999       Atom a1 = bond.getAtom(0);
2000       Atom a2 = bond.getAtom(1);
2001       if (a1.isHydrogen() && a2.isHeavy()) {
2002         double hydrogenMass = a1.getMass();
2003         double heavyAtomMass = a2.getMass();
2004         if (hydrogenMass < 1.1) {
2005           a2.setMass(heavyAtomMass - 2.0 * hydrogenMass);
2006           a1.setMass(3.0 * hydrogenMass);
2007         }
2008       } else if (a1.isHeavy() && a2.isHydrogen()) {
2009         double heavyAtomMass = a1.getMass();
2010         double hydrogenMass = a2.getMass();
2011         if (hydrogenMass < 1.1) {
2012           a1.setMass(heavyAtomMass - 2.0 * hydrogenMass);
2013           a2.setMass(3.0 * hydrogenMass);
2014         }
2015       }
2016     }
2017   }
2018 
2019   private Atom getResidue(Atom atom, boolean create) {
2020     return getResidue(atom, create, Residue.ResidueType.UNK);
2021   }
2022 
2023   private Atom getResidue(Atom atom, boolean create, ResidueType defaultRT) {
2024     Character chainID = atom.getChainID();
2025     String resName = atom.getResidueName();
2026     int resNum = atom.getResidueNumber();
2027     String segID = atom.getSegID();
2028     // Find/Create the chain
2029     Polymer polymer = getPolymer(chainID, segID, create);
2030     if (polymer == null) {
2031       return null;
2032     }
2033 
2034     Residue res;
2035     if(titrateConformer){
2036       res = polymer.getResidue(atomInitial.getResidueName(), atomInitial.getResidueNumber(), create, defaultRT);
2037       res.setName(atom.getResidueName());
2038     } else {
2039       res = polymer.getResidue(resName, resNum, create, defaultRT);
2040     }
2041     if (create && res != null) {
2042       res.setTitrateConformers(titrateConformer);
2043       if(titrateConformer){
2044         res.setAtomInitial(atomInitial);
2045       }
2046       return (Atom) res.addMSNode(atom);
2047     }
2048     return null;
2049   }
2050 
2051   private Atom getMolecule(Atom atom, boolean create) {
2052     String resName = atom.getResidueName();
2053     int resNum = atom.getResidueNumber();
2054     Character chainID = atom.getChainID();
2055     String segID = atom.getSegID();
2056 
2057     String key = resNum + resName + segID;
2058     Molecule m = ionHashMap.get(key);
2059     if (m == null) {
2060       m = waterHashMap.get(key);
2061       if (m == null) {
2062         m = moleculeHashMap.get(key);
2063       }
2064     }
2065     if (m != null) {
2066       return (Atom) m.addMSNode(atom);
2067     }
2068 
2069     if (create) {
2070       boolean isWater = false;
2071       boolean isIon = false;
2072       if (StringUtils.looksLikeWater(resName)) {
2073         if (!resName.equalsIgnoreCase("DOD")) {
2074           resName = StringUtils.STANDARD_WATER_NAME;
2075         }
2076         isWater = true;
2077       } else if (StringUtils.looksLikeIon(resName)) {
2078         resName = StringUtils.tryParseIon(resName);
2079         isIon = true;
2080       }
2081       atom.setResName(resName);
2082       m = new Molecule(resName, resNum, chainID, segID);
2083       m.addMSNode(atom);
2084       if (resName == null) {
2085         logger.warning(format(" Attempting to create a molecule %s with a null name on atom %s! Defaulting to creating a generic Molecule.",
2086                 m, atom));
2087         molecules.add(m);
2088         moleculeHashMap.put(key, m);
2089       } else if (isWater) {
2090         water.add(m);
2091         waterHashMap.put(key, m);
2092       } else if (isIon) {
2093         ions.add(m);
2094         ionHashMap.put(key, m);
2095       } else {
2096         molecules.add(m);
2097         moleculeHashMap.put(key, m);
2098       }
2099       return atom;
2100     } else {
2101       return null;
2102     }
2103   }
2104 
2105   private void recurseVRML(Node node) {
2106     if (node instanceof Shape3D) {
2107       Shape3D s3d = (Shape3D) node;
2108       PickTool.setCapabilities(s3d, PickTool.INTERSECT_COORD);
2109     } else if (node instanceof SharedGroup) {
2110       SharedGroup sg = (SharedGroup) node;
2111       for (Iterator<Node> e = sg.getAllChildren(); e.hasNext(); ) {
2112         recurseVRML(e.next());
2113       }
2114     } else if (node instanceof BranchGroup) {
2115       BranchGroup bg = (BranchGroup) node;
2116       for (Iterator<Node> e = bg.getAllChildren(); e.hasNext(); ) {
2117         recurseVRML(e.next());
2118       }
2119     } else if (node instanceof TransformGroup) {
2120       TransformGroup vrmlTG1 = (TransformGroup) node;
2121       for (Iterator<Node> e = vrmlTG1.getAllChildren(); e.hasNext(); ) {
2122         node = e.next();
2123         recurseVRML(node);
2124       }
2125     } else if (node instanceof Link) {
2126       Link link = (Link) node;
2127       recurseVRML(link.getSharedGroup());
2128     } else if (node instanceof Group) {
2129       Group group = (Group) node;
2130       for (Iterator<Node> e = group.getAllChildren(); e.hasNext(); ) {
2131         recurseVRML(e.next());
2132       }
2133     }
2134   }
2135 
2136   /** {@inheritDoc} */
2137   @Override
2138   protected void removeLeaves() {
2139     super.removeLeaves();
2140     MSNode macroNode = getAtomNode();
2141     if (macroNode != null) {
2142       if (macroNode.getChildCount() > 0) {
2143         getAtomNode().setName("Macromolecules " + "(" + macroNode.getChildCount() + ")");
2144       } else if (macroNode.getParent() == this) {
2145         removeChild(macroNode);
2146       }
2147     }
2148 
2149     if (molecules.getChildCount() == 0) {
2150       removeChild(molecules);
2151     } else {
2152       molecules.setName("Hetero Molecules " + "(" + molecules.getChildCount() + ")");
2153     }
2154 
2155     if (ions.getChildCount() == 0) {
2156       removeChild(ions);
2157     } else {
2158       ions.setName("Ions " + "(" + ions.getChildCount() + ")");
2159     }
2160 
2161     if (water.getChildCount() == 0) {
2162       removeChild(water);
2163     } else {
2164       water.setName("Water " + "(" + water.getChildCount() + ")");
2165     }
2166   }
2167 
2168   private Shape3D renderWire() {
2169     List<Bond> bonds = getBondList();
2170     int numbonds = bonds.size();
2171     if (numbonds < 1) {
2172       return null;
2173     }
2174 
2175     Vector3d bondmidpoint = new Vector3d();
2176     double[] mid = {0, 0, 0};
2177     Vector3d v1 = new Vector3d();
2178     Vector3d v2 = new Vector3d();
2179     float[] a1 = {0, 0, 0};
2180     float[] a2 = {0, 0, 0};
2181     float[] col = new float[4];
2182 
2183     Atom atom1, atom2;
2184     LineArray la =
2185         new LineArray(
2186             4 * numbonds,
2187             GeometryArray.COORDINATES | GeometryArray.COLOR_4 | GeometryArray.NORMALS);
2188     la.setCapability(LineArray.ALLOW_COORDINATE_WRITE);
2189     la.setCapability(LineArray.ALLOW_COORDINATE_READ);
2190     la.setCapability(LineArray.ALLOW_COLOR_WRITE);
2191     la.setCapability(LineArray.ALLOW_COUNT_READ);
2192     la.setCapability(LineArray.ALLOW_INTERSECT);
2193     la.setCapability(LineArray.ALLOW_FORMAT_READ);
2194     atomLookUp = new Atom[4 * numbonds];
2195     int i = 0;
2196     col[3] = 0.9f;
2197     for (Bond bond : bonds) {
2198       bond.setWire(la, i);
2199       atom1 = bond.getAtom(0);
2200       atom2 = bond.getAtom(1);
2201       atom1.getV3D(v1);
2202       atom2.getV3D(v2);
2203       a1[0] = (float) v1.x;
2204       a1[1] = (float) v1.y;
2205       a1[2] = (float) v1.z;
2206       a2[0] = (float) v2.x;
2207       a2[1] = (float) v2.y;
2208       a2[2] = (float) v2.z;
2209       // Find the bond center
2210       bondmidpoint.add(v1, v2);
2211       bondmidpoint.scale(0.5d);
2212       bondmidpoint.get(mid);
2213 
2214       // Atom #1
2215       Atom.AtomColor.get(atom1.getAtomicNumber()).get(col);
2216       atomLookUp[i] = atom1;
2217       la.setCoordinate(i, a1);
2218       la.setColor(i, col);
2219       la.setNormal(i, a2);
2220       i++;
2221 
2222       atomLookUp[i] = atom1;
2223       la.setCoordinate(i, mid);
2224       la.setColor(i, col);
2225       la.setNormal(i, a2);
2226       i++;
2227 
2228       // Atom #2
2229       Atom.AtomColor.get(atom2.getAtomicNumber()).get(col);
2230       atomLookUp[i] = atom2;
2231       la.setCoordinate(i, a2);
2232       la.setColor(i, col);
2233       la.setNormal(i, a1);
2234       i++;
2235 
2236       atomLookUp[i] = atom2;
2237       la.setCoordinate(i, mid);
2238       la.setColor(i, col);
2239       la.setNormal(i, a1);
2240       i++;
2241     }
2242 
2243     ColoringAttributes cola =
2244         new ColoringAttributes(new Color3f(), ColoringAttributes.SHADE_GOURAUD);
2245     Appearance app = new Appearance();
2246     lineAttributes = new LineAttributes();
2247     lineAttributes.setLineWidth(RendererCache.bondwidth);
2248     lineAttributes.setCapability(LineAttributes.ALLOW_WIDTH_WRITE);
2249     lineAttributes.setLineAntialiasingEnable(true);
2250     app.setLineAttributes(lineAttributes);
2251     app.setCapability(Appearance.ALLOW_LINE_ATTRIBUTES_READ);
2252     app.setCapability(Appearance.ALLOW_LINE_ATTRIBUTES_WRITE);
2253     RenderingAttributes ra = new RenderingAttributes();
2254     ra.setAlphaTestValue(0.1f);
2255     ra.setAlphaTestFunction(RenderingAttributes.GREATER);
2256     ra.setDepthBufferEnable(true);
2257     ra.setDepthBufferWriteEnable(true);
2258     app.setRenderingAttributes(ra);
2259     app.setColoringAttributes(cola);
2260     Shape3D wireframe = new Shape3D(la, app);
2261     // PickTool.setCapabilities(wire, PickTool.INTERSECT_COORD);
2262     wireframe.setUserData(this);
2263     wireframe.setBounds(new BoundingSphere(new Point3d(0, 0, 0), 1000.0));
2264     try {
2265       wireframe.setBoundsAutoCompute(false);
2266     } catch (Exception e) {
2267       logger.warning("Unable to set boundsAutoCompute to false.\n" + e);
2268     }
2269 
2270     wireframe.setCapability(Shape3D.ALLOW_GEOMETRY_READ);
2271     wireframe.setCapability(Shape3D.ALLOW_APPEARANCE_READ);
2272     wireframe.setCapability(Shape3D.ALLOW_LOCAL_TO_VWORLD_READ);
2273     return wireframe;
2274   }
2275 
2276   /**
2277    * Returns a list of all water molecules in the system, including both those properly under the
2278    * water node and those under the molecules node.
2279    *
2280    * @return All water in the system
2281    */
2282   private List<Molecule> getAllWaterInclMistyped() {
2283     Stream<Molecule> mistyped =
2284         getMolecules()
2285             .parallelStream()
2286             .map((MSNode m) -> (Molecule) m)
2287             .filter(
2288                 (Molecule m) -> {
2289                   List<Atom> atoms = m.getAtomList();
2290                   if (atoms.size() != 3) {
2291                     return false;
2292                   }
2293                   int nO = 0;
2294                   int nH = 0;
2295                   for (Atom atom : atoms) {
2296                     int el = atom.getAtomicNumber();
2297                     if (el == 1) {
2298                       ++nH;
2299                     } else if (nH == 8) {
2300                       ++nO;
2301                     }
2302                   }
2303                   return nO == 1 && nH == 2;
2304                 });
2305     return Stream.concat(mistyped, getWater().stream().map((MSNode m) -> (Molecule) m))
2306         .distinct()
2307         .collect(Collectors.toList());
2308   }
2309 
2310   /** Renames water protons to H1 and H2. */
2311   void renameWaterProtons() {
2312     for (Molecule water : getAllWaterInclMistyped()) {
2313       Atom H1 = water.getAtomByName("H1", false);
2314       Atom H2 = water.getAtomByName("H2", false);
2315       if (H1 != null && H2 != null) {
2316         continue;
2317       }
2318       for (Atom a : water.getAtomList()) {
2319         if (a.getAtomicNumber() == 1) {
2320           if (H1 == null) {
2321             H1 = a;
2322             H1.setName("H1");
2323           } else if (H2 == null) {
2324             H2 = a;
2325             H2.setName("H2");
2326           }
2327         }
2328       }
2329     }
2330   }
2331 
2332   private void moveIntoUnitCell(MSNode[] groups) {
2333     if (groups != null && groups.length > 0) {
2334       moveIntoUnitCell(Arrays.asList(groups));
2335     }
2336   }
2337 
2338   /**
2339    * Move the center of each listed chemical entity into the primary unit cell.
2340    *
2341    * @param groups Move each MSNode into the primary unit cell.
2342    */
2343   private void moveIntoUnitCell(List<MSNode> groups) {
2344     Crystal cryst = getCrystal().getUnitCell();
2345     if (cryst.aperiodic()) {
2346       return;
2347     }
2348     for (MSNode group : groups) {
2349       double[] com = new double[3];
2350       double[] xyz = new double[3];
2351       double[] translate = new double[3];
2352       List<Atom> atoms = group.getAtomList();
2353       double totMass = 0;
2354       for (Atom atom : atoms) {
2355         double mass = atom.getMass();
2356         totMass += mass;
2357         xyz = atom.getXYZ(xyz);
2358         com[0] += mass * xyz[0];
2359         com[1] += mass * xyz[1];
2360         com[2] += mass * xyz[2];
2361       }
2362       com[0] /= totMass;
2363       com[1] /= totMass;
2364       com[2] /= totMass;
2365 
2366       // Move the COM to the primary unit cell
2367       cryst.toPrimaryCell(com, translate);
2368 
2369       // The translation vector is difference between the new location and the current COM.
2370       translate[0] -= com[0];
2371       translate[1] -= com[1];
2372       translate[2] -= com[2];
2373 
2374       for (Atom atom : atoms) {
2375         atom.move(translate);
2376       }
2377     }
2378   }
2379 
2380   public enum FractionalMode {
2381     OFF,
2382     MOLECULE,
2383     ATOM
2384   }
2385 }