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