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