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