View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2025.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.potential.parsers;
39  
40  import ffx.crystal.Crystal;
41  import ffx.crystal.LatticeSystem;
42  import ffx.crystal.SpaceGroup;
43  import ffx.crystal.SpaceGroupDefinitions;
44  import ffx.potential.MolecularAssembly;
45  import ffx.potential.Utilities;
46  import ffx.potential.bonded.Angle;
47  import ffx.potential.bonded.Atom;
48  import ffx.potential.bonded.Bond;
49  import ffx.potential.bonded.MSGroup;
50  import ffx.potential.bonded.MSNode;
51  import ffx.potential.bonded.Molecule;
52  import ffx.potential.bonded.Polymer;
53  import ffx.potential.bonded.Residue;
54  import ffx.potential.parameters.ForceField;
55  import org.apache.commons.configuration2.CompositeConfiguration;
56  import org.openscience.cdk.AtomContainer;
57  import org.openscience.cdk.config.AtomTypeFactory;
58  import org.openscience.cdk.graph.rebond.RebondTool;
59  import org.openscience.cdk.interfaces.IAtom;
60  import org.openscience.cdk.interfaces.IAtomType;
61  import org.openscience.cdk.interfaces.IBond;
62  import org.openscience.cdk.isomorphism.AtomMatcher;
63  import org.openscience.cdk.isomorphism.BondMatcher;
64  import org.openscience.cdk.isomorphism.Pattern;
65  import org.openscience.cdk.isomorphism.VentoFoggia;
66  import org.rcsb.cif.CifIO;
67  import org.rcsb.cif.model.Column;
68  import org.rcsb.cif.model.FloatColumn;
69  import org.rcsb.cif.schema.StandardSchemata;
70  import org.rcsb.cif.schema.core.AtomSite;
71  import org.rcsb.cif.schema.core.Cell;
72  import org.rcsb.cif.schema.core.Chemical;
73  import org.rcsb.cif.schema.core.CifCoreBlock;
74  import org.rcsb.cif.schema.core.CifCoreFile;
75  import org.rcsb.cif.schema.core.Symmetry;
76  
77  import javax.vecmath.Point3d;
78  import java.io.BufferedReader;
79  import java.io.BufferedWriter;
80  import java.io.File;
81  import java.io.FileWriter;
82  import java.io.IOException;
83  import java.nio.file.Path;
84  import java.nio.file.Paths;
85  import java.util.ArrayList;
86  import java.util.Arrays;
87  import java.util.List;
88  import java.util.Objects;
89  import java.util.logging.Level;
90  import java.util.logging.Logger;
91  
92  import static ffx.crystal.SpaceGroupConversions.hrConversion;
93  import static ffx.numerics.math.DoubleMath.dihedralAngle;
94  import static ffx.numerics.math.DoubleMath.dist;
95  import static ffx.potential.bonded.BondedUtils.intxyz;
96  import static java.lang.String.format;
97  import static org.apache.commons.io.FilenameUtils.getExtension;
98  import static org.apache.commons.io.FilenameUtils.getFullPath;
99  import static org.apache.commons.io.FilenameUtils.getName;
100 import static org.apache.commons.io.FilenameUtils.removeExtension;
101 import static org.apache.commons.math3.util.FastMath.sqrt;
102 import static org.openscience.cdk.tools.periodictable.PeriodicTable.getCovalentRadius;
103 import static org.openscience.cdk.tools.periodictable.PeriodicTable.getSymbol;
104 
105 /**
106  * The CIFFilter class parses CIF coordinate (*.CIF) files.
107  *
108  * @author Aaron J. Nessler
109  * @since 1.0
110  */
111 public class CIFFilter extends SystemFilter {
112   /**
113    * A logger for this class.
114    */
115   private static final Logger logger = Logger.getLogger(CIFFilter.class.getName());
116   /**
117    * Set bond tolerance (distance to determine if bond is made).
118    */
119   private double bondTolerance = 0.2;
120   /**
121    * The CIF file reader.
122    */
123   private final BufferedReader bufferedReader = null;
124   /**
125    * The CIF file object.
126    */
127   private CifCoreFile cifFile;
128   /**
129    * List of output files created from converted CIF file.
130    */
131   private final List<String> createdFiles = new ArrayList<>();
132   /**
133    * The current snapshot.
134    */
135   private int snapShot;
136   /**
137    * The space group name.
138    */
139   private String sgName = null;
140   /**
141    * The space group number.
142    */
143   private int sgNum = -1;
144   /**
145    * The crystallographic information for this system.
146    */
147   private Crystal crystal;
148   /**
149    * Whether to fix lattice parameters.
150    */
151   private boolean fixLattice = false;
152   /**
153    * Whether to use PDB format.
154    */
155   private boolean usePDB = false;
156   /**
157    * The number of entities in the asymmetric unit.
158    */
159   private int zPrime = -1;
160   /**
161    * Maximum atomic covalent radius for CDK Rebonder Tool
162    */
163   public static final double MAX_COVALENT_RADIUS = 2.0;
164 
165   /**
166    * Minimum bond distance for CDK Rebonder Tool
167    */
168   public static final double MIN_BOND_DISTANCE = 0.5;
169 
170   /**
171    * Constructor for CIFFilter on a single file and a single assembly.
172    *
173    * @param file Input file
174    * @param molecularAssembly Active assembly
175    * @param forceField Force field for save file.
176    * @param properties Properties for save file.
177    */
178   public CIFFilter(File file, MolecularAssembly molecularAssembly, ForceField forceField,
179       CompositeConfiguration properties, boolean saveCIF) {
180     super(file, molecularAssembly, forceField, properties);
181     this.fileType = Utilities.FileType.CIF;
182 
183     if (!saveCIF) {
184       String dir = getFullPath(molecularAssembly.getFile().getAbsolutePath()) + File.separator;
185       String cifName = removeExtension(file.getName()) + ".cif";
186       Path path = Paths.get(dir + cifName);
187       try {
188         cifFile = CifIO.readFromPath(path).as(StandardSchemata.CIF_CORE);
189       } catch (Exception ex) {
190         logger.info(" Failed to create CIF file object.\n " + ex);
191       }
192       String ext = getExtension(file.getName());
193       currentFile = new File(removeExtension(currentFile.getAbsolutePath()) + ext);
194       if(ext.equalsIgnoreCase("pdb")){
195         usePDB = true;
196       }
197     }
198   }
199 
200   /**
201    * Constructor for CIFFilter on a single file and multiple assemblies.
202    *
203    * @param file Input file
204    * @param molecularAssemblies Active assemblies
205    * @param forceField Force field for save file.
206    * @param properties Properties for save file.
207    */
208   public CIFFilter(File file, List<MolecularAssembly> molecularAssemblies, ForceField forceField,
209       CompositeConfiguration properties, boolean saveCIF) {
210     super(file, molecularAssemblies, forceField, properties);
211     this.fileType = Utilities.FileType.CIF;
212 
213     if (!saveCIF) {
214       String dir =
215           getFullPath(molecularAssemblies.get(0).getFile().getAbsolutePath()) + File.separator;
216       String cifName = removeExtension(file.getName()) + ".cif";
217       Path path = Paths.get(dir + cifName);
218       try {
219         cifFile = CifIO.readFromPath(path).as(StandardSchemata.CIF_CORE);
220       } catch (Exception ex) {
221         logger.info(" Failed to create CIF file object.\n" + ex);
222       }
223       currentFile = new File(removeExtension(currentFile.getAbsolutePath()) + ".xyz");
224     }
225   }
226 
227   /**
228    * Constructor for CIFFilter on a multiple files and a single assembly.
229    *
230    * @param files Input files
231    * @param molecularAssembly Active assembly
232    * @param forceField Force field for save file.
233    * @param properties Properties for save file.
234    */
235   public CIFFilter(List<File> files, MolecularAssembly molecularAssembly, ForceField forceField,
236       CompositeConfiguration properties, boolean saveCIF) {
237     super(files, molecularAssembly, forceField, properties);
238     this.fileType = Utilities.FileType.CIF;
239 
240     if (!saveCIF) {
241       String dir = getFullPath(molecularAssembly.getFile().getAbsolutePath()) + File.separator;
242       String cifName = removeExtension(files.get(0).getName()) + ".cif";
243       Path path = Paths.get(dir + cifName);
244       try {
245         cifFile = CifIO.readFromPath(path).as(StandardSchemata.CIF_CORE);
246       } catch (Exception ex) {
247         logger.info(" Failed to create CIF file object.\n" + ex);
248       }
249       currentFile = new File(removeExtension(currentFile.getAbsolutePath()) + ".xyz");
250     }
251   }
252 
253   /**
254    * Add bonds between atoms.
255    *
256    * @param atoms To potentially be bonded together.
257    */
258   public static int bondAtoms(Atom[] atoms, double bondTolerance) {
259     int bondCount = 0;
260     int nAtoms = atoms.length;
261     for (int i = 0; i < nAtoms; i++) {
262       Atom atom1 = atoms[i];
263       String atomIelement = getAtomElement(atom1);
264       double radiusI =
265           (getCovalentRadius(atomIelement) == null) ? 0 : getCovalentRadius(atomIelement);
266       if (radiusI == 0) {
267         logger.warning(format(" Atom %d (%s) element (%s) not determined.", i + 1, atom1.getName(),
268             atomIelement));
269         return -1;
270       }
271       double[] xyz = {atom1.getX(), atom1.getY(), atom1.getZ()};
272       for (int j = i + 1; j < nAtoms; j++) {
273         Atom atom2 = atoms[j];
274         String atomJelement = getAtomElement(atom2);
275         double radiusJ =
276             (getCovalentRadius(atomJelement) == null) ? 0 : getCovalentRadius(atomJelement);
277         if (radiusJ == 0) {
278           logger.warning(format(" Atom %d element (%s) not determined.", j + 1, atomJelement));
279           return -1;
280         }
281         double[] xyz2 = {atom2.getX(), atom2.getY(), atom2.getZ()};
282         double length = dist(xyz, xyz2);
283         double bondLength = radiusI + radiusJ + bondTolerance;
284         if (length < bondLength) {
285           bondCount++;
286           Bond bond = new Bond(atom1, atom2);
287           atom1.setBond(bond);
288           atom2.setBond(bond);
289           logger.finest(format(
290               " Bonded atom %d (%s) with atom %d (%s): bond length (%4.4f Å) < tolerance (%4.4f Å)",
291               i + 1, atomIelement, j + 1, atomJelement, length, bondLength));
292         }
293       }
294     }
295     return bondCount;
296   }
297 
298   /** {@inheritDoc} */
299   @Override
300   public void closeReader() {
301     if (bufferedReader != null) {
302       try {
303         bufferedReader.close();
304       } catch (IOException ex) {
305         logger.warning(format(" Exception in closing CIF filter: %s", ex));
306       }
307     }
308   }
309 
310   /**
311    * Finds all atoms that are bonded to one another. See potential/Utilities/collectAtoms
312    *
313    * @param seed Starting atom.
314    * @param atoms that are bonded.
315    */
316   public static void collectAtoms(Atom seed, ArrayList<Atom> atoms) {
317     logger.finest(format(" Atom: %s", seed.getName()));
318     atoms.add(seed);
319     for (Bond b : seed.getBonds()) {
320       Atom nextAtom = b.get1_2(seed);
321       if (!atoms.contains(nextAtom)) {
322         collectAtoms(nextAtom, atoms);
323       }
324     }
325   }
326 
327   /**
328    * Parse atom name to determine atomic element.
329    *
330    * @param atom Atom whose element we desire
331    * @return String specifying atom element.
332    */
333   public static String getAtomElement(Atom atom) {
334     return getAtomElement(atom.getName());
335   }
336 
337   /**
338    * Parse atom name to determine atomic element.
339    *
340    * @param name Name of atom whose element we desire
341    * @return String specifying atom element.
342    */
343   private static String getAtomElement(String name) {
344     String value = "";
345     try {
346       value = name.replaceAll("[()]", "").replaceAll("_", "").replaceAll("-", "")
347           .replaceAll(" +", "").split("[0-9]")[0];
348     } catch (Exception e) {
349       logger.severe(" Error extracting atom element. Please ensure the CIF is formatted correctly.\n" + e);
350     }
351     return value;
352   }
353 
354   /**
355    * Obtain a list of output files written from the conversion. Used in file not committed to Git...
356    *
357    * @return String[] containing output file names.
358    */
359   public String[] getCreatedFileNames() {
360     String[] files = new String[createdFiles.size()];
361     return createdFiles.toArray(files);
362   }
363 
364   /**
365    * Parse CIF block to determine crystal and atom information.
366    * @param block CifCoreBlock object to be interpreted.
367    * @return Atoms from CIF. (Crystal updated as global variable)
368    */
369   private Atom[] parseBlock(CifCoreBlock block){
370     // Parse CIF file "blocks"
371     logger.info("\n Block ID: " + block.getBlockHeader());
372     Chemical chemical = block.getChemical();
373     Column<String[]> nameCommon = chemical.getNameCommon();
374     Column<String[]> nameSystematic = chemical.getNameSystematic();
375     int rowCount = nameCommon.getRowCount();
376     if (rowCount > 1) {
377       logger.info(" Chemical components");
378       for (int i = 0; i < rowCount; i++) {
379         logger.info(format("  %s", nameCommon.getColumnName()));
380       }
381     } else if (rowCount > 0) {
382       logger.info(format(" Chemical component: %s", nameCommon.getColumnName()));
383     }
384 
385     // Determine the space group from CIF.
386     Symmetry symmetry = block.getSymmetry();
387     if (sgNum == -1 && sgName == null || sgName.isEmpty()) {
388       if (symmetry.getIntTablesNumber().getRowCount() > 0) {
389         sgNum = symmetry.getIntTablesNumber().get(0);
390         logger.info(format(" CIF International Tables Number: %d", sgNum));
391       }
392       if (symmetry.getSpaceGroupNameH_M().getRowCount() > 0) {
393         sgName = symmetry.getSpaceGroupNameH_M().get(0);
394         logger.info(format(" CIF Hermann–Mauguin Space Group: %s", sgName));
395       } else if (block.getSpaceGroup().getNameH_mFull().getRowCount() > 0) {
396         sgName = block.getSpaceGroup().getNameH_mFull().get(0);
397         logger.info(format(" CIF Hermann–Mauguin Space Group: %s", sgName));
398       } else if (block.getSpaceGroup().getNameH_mAlt().getRowCount() > 0) {
399         sgName = block.getSpaceGroup().getNameH_mAlt().get(0);
400         logger.info(format(" CIF Hermann–Mauguin Space Group: %s", sgName));
401       }
402     } else {
403       // Parse user defined space group.
404       if (sgNum != -1) {
405         logger.info(format(" Command line International Tables Number: %d", sgNum));
406       } else {
407         logger.info(format(" Command line space group name: %s", sgName));
408       }
409     }
410 
411     // Set space group.
412     SpaceGroup sg = null;
413     if (sgName != null) {
414       sg = SpaceGroupDefinitions.spaceGroupFactory(sgName);
415     }
416     if (sg == null) {
417       logger.finer(" Space group name not found. Attempting to use space group number.");
418       sg = SpaceGroupDefinitions.spaceGroupFactory(sgNum);
419     }
420 
421     // Fall back to P1.
422     if (sg == null) {
423       logger.warning(" The space group could not be determined from the CIF file (using P1).");
424       sg = SpaceGroupDefinitions.spaceGroupFactory(1);
425     }
426 
427     // Generate FFX Crystal for CIF file.
428     Cell cell = block.getCell();
429     if(cell.isDefined()) {
430       double a = cell.getLengthA().get(0);
431       double b = cell.getLengthB().get(0);
432       double c = cell.getLengthC().get(0);
433       double alpha = cell.getAngleAlpha().get(0);
434       double beta = cell.getAngleBeta().get(0);
435       double gamma = cell.getAngleGamma().get(0);
436       assert sg != null;
437       LatticeSystem latticeSystem = sg.latticeSystem;
438       double[] latticeParameters = {a, b, c, alpha, beta, gamma};
439       if (latticeSystem.validParameters(a, b, c, alpha, beta, gamma)) {
440         crystal = new Crystal(a, b, c, alpha, beta, gamma, sg.pdbName);
441       } else {
442         // Fix lattice parameters if they disobey space group.
443         if (fixLattice) {
444           logger.info(
445                   " Attempting to patch disagreement between lattice system and lattice parameters.");
446           boolean fixed = false;
447           // Check if Rhombohedral lattice has been named Hexagonal
448           if (latticeSystem == LatticeSystem.HEXAGONAL_LATTICE
449                   && LatticeSystem.RHOMBOHEDRAL_LATTICE.validParameters(a, b, c, alpha, beta, gamma)) {
450             crystal = hrConversion(a, b, c, alpha, beta, gamma, sg);
451             latticeSystem = crystal.spaceGroup.latticeSystem;
452             if (latticeSystem.validParameters(crystal.a, crystal.b, crystal.c, crystal.alpha,
453                     crystal.beta, crystal.gamma)) {
454               fixed = true;
455             }
456             // Check if Hexagonal lattice has been named Rhombohedral
457           } else if (latticeSystem == LatticeSystem.RHOMBOHEDRAL_LATTICE
458                   && LatticeSystem.HEXAGONAL_LATTICE.validParameters(a, b, c, alpha, beta, gamma)) {
459             crystal = hrConversion(a, b, c, alpha, beta, gamma, sg);
460             latticeSystem = crystal.spaceGroup.latticeSystem;
461             if (latticeSystem.validParameters(crystal.a, crystal.b, crystal.c, crystal.alpha,
462                     crystal.beta, crystal.gamma)) {
463               fixed = true;
464             }
465           }
466           if (!fixed) {
467             double[] newLatticeParameters = latticeSystem.fixParameters(a, b, c, alpha, beta, gamma);
468             if (newLatticeParameters == latticeParameters) {
469               logger.warning(" Conversion Failed: The proposed lattice parameters for " + sg.pdbName
470                       + " do not satisfy the " + latticeSystem + ".");
471               return null;
472             } else {
473               a = newLatticeParameters[0];
474               b = newLatticeParameters[1];
475               c = newLatticeParameters[2];
476               alpha = newLatticeParameters[3];
477               beta = newLatticeParameters[4];
478               gamma = newLatticeParameters[5];
479               crystal = new Crystal(a, b, c, alpha, beta, gamma, sg.pdbName);
480             }
481           }
482         } else {
483           logger.warning(" Conversion Failed: The proposed lattice parameters for " + sg.pdbName
484                   + " do not satisfy the " + latticeSystem + ".");
485           logger.info(" Use \"--fixLattice\" or \"--fl\" flag to attempt to fix automatically.");
486           return null;
487         }
488       }
489       activeMolecularAssembly.setName(block.getBlockHeader());
490       logger.info(" New Crystal: " + crystal);
491       activeMolecularAssembly.setCrystal(crystal);
492     }
493     // Parse CIF Atoms
494     AtomSite atomSite = block.getAtomSite();
495     Column<String[]> label = atomSite.getLabel();
496     Column<String[]> typeSymbol = atomSite.getTypeSymbol();
497     FloatColumn fractX = atomSite.getFractX();
498     FloatColumn fractY = atomSite.getFractY();
499     FloatColumn fractZ = atomSite.getFractZ();
500 
501     FloatColumn cartX = atomSite.getCartnX();
502     FloatColumn cartY = atomSite.getCartnY();
503     FloatColumn cartZ = atomSite.getCartnZ();
504 
505     int nAtoms = label.getRowCount();
506     if (nAtoms < 1) {
507       logger.warning(" CIF file did not contain coordinates.");
508       return null;
509     }
510     if (logger.isLoggable(Level.FINE)) {
511       logger.fine(format("\n Number of Atoms in CIF: %d", nAtoms));
512     }
513     Atom[] atoms = new Atom[nAtoms];
514 
515     // Define per atom information from the CIF.
516     String resName = "CIF";
517     double bfactor = 1.0;
518     char altLoc = ' ';
519     char chain = 'A';
520 //    String segID = "A";
521     String[] symbols = new String[nAtoms];
522 
523     // Loop over atoms in CIF and generate FFX Atom(s).
524     for (int i = 0; i < nAtoms; i++) {
525       double occupancy = 1.0;
526       if(atomSite.getOccupancy().isDefined()){
527         occupancy = atomSite.getOccupancy().get(i);
528       }
529       // Assigning each atom their own resID prevents comparator from treating them as the same atom.
530       if (typeSymbol.getRowCount() > 0) {
531         symbols[i] = typeSymbol.getStringData(i);
532       } else {
533         symbols[i] = getAtomElement(label.getStringData(i));
534       }
535       double x = (fractX.isDefined()) ? fractX.get(i) : cartX.get(i);
536       double y = (fractY.isDefined()) ? fractY.get(i) : cartY.get(i);
537       double z = (fractZ.isDefined()) ? fractZ.get(i) : cartZ.get(i);
538       double[] xyz = {x, y, z};
539       if (fractX.isDefined() && crystal != null) {
540         crystal.toCartesianCoordinates(xyz, xyz);
541       }
542       atoms[i] = new Atom(i + 1, label.getStringData(i), altLoc, xyz, resName, i, chain, occupancy,
543               bfactor, symbols[i]);
544       atoms[i].setHetero(true);
545       if (logger.isLoggable(Level.FINE)) {
546         logger.fine(format(
547                 " Atom (%2d) Name: " + atoms[i].getName() + " Label: " + label.getStringData(i)
548                         + " Symbol: " + symbols[i], i));
549       }
550     }
551     return atoms;
552   }
553 
554   /**
555    * {@inheritDoc}
556    *
557    * <p>Parse the CIF File
558    */
559   @Override
560   public boolean readFile() {
561     // Open the input file (with electrostatics turned off).
562     System.setProperty("mpoleterm", "false");
563     Atom[] inputAtoms = activeMolecularAssembly.getAtomArray();
564     int nInputAtoms = inputAtoms.length;
565 
566     int numHydrogen = 0;
567     for (Atom atom : inputAtoms) {
568       if (atom.isHydrogen()) {
569         numHydrogen++;
570       }
571     }
572 
573     List<MSNode> entitiesInput = activeMolecularAssembly.getNodeList();
574     int numEntitiesInput = entitiesInput.size();
575     if (logger.isLoggable(Level.FINE)) {
576       logger.fine(format(" Number of entities in input: %d", numEntitiesInput));
577       for (MSNode entity : entitiesInput) {
578         logger.fine(format(" Entity: " + entity.getName()));
579         int size = entity.getAtomList().size();
580         logger.fine(format("   Entity Size: %3d", size));
581         if (size > 0) {
582           logger.fine(format("   Entity First Atom: " + entity.getAtomList().get(0).toString()));
583         } else {
584           logger.warning(" Entity did not contain atoms.");
585         }
586       }
587     }
588 
589     for (CifCoreBlock block : cifFile.getBlocks()) {
590       Atom[] atoms = parseBlock(block);
591       logger.info(crystal.toString());
592       if(atoms == null || atoms.length == 0){
593         logger.warning(" Atoms could not be obtained from CIF.");
594         return false;
595       }
596 
597       int nAtoms = atoms.length;
598 
599       // Define assembly object to contain CIF info.
600       MolecularAssembly outputAssembly = new MolecularAssembly(block.getBlockHeader());
601       int zPrime;
602       if (this.zPrime > 0) {
603         zPrime = this.zPrime;
604       } else if (nAtoms % nInputAtoms == 0) {
605         zPrime = nAtoms / nInputAtoms;
606       } else if (nAtoms % (nInputAtoms - numHydrogen) == 0) {
607         zPrime = nAtoms / (nInputAtoms - numHydrogen);
608       } else {
609         zPrime = 1;
610       }
611 
612       if (zPrime > 1) {
613         logger.info(format(" Detected more than one copy in asymmetric unit of CIF file (Z'=%d)."
614             + " -- attempting to separate.", zPrime));
615       }
616       ArrayList<ArrayList<Atom>> xyzatoms = new ArrayList<>();
617       int atomIndex = 1;
618       int[] moleculeNumbers = activeMolecularAssembly.getMoleculeNumbers();
619       if(logger.isLoggable(Level.FINER)) {
620         for (int i = 0; i < moleculeNumbers.length; i++) {
621           logger.finer(format(" Molecule Number (atom %2d): %2d", i, moleculeNumbers[i]));
622         }
623       }
624       // For each entity (protein, molecule, ion, etc) in the XYZ
625       for (int i = 0; i < numEntitiesInput; i++) {
626         MSNode mol = entitiesInput.get(i);
627         int shiftIndex = mol.getAtomList().get(0).getXyzIndex() - 1;
628         logger.info(format(" Size of Entity: %2d", mol.getAtomList().size()));
629         xyzatoms.add((ArrayList<Atom>) mol.getAtomList());
630         ArrayList<Atom> currentAtoms = xyzatoms.get(i);
631         // Set up input (XYZ) file contents as CDK variable
632         AtomContainer xyzCDKAtoms = new AtomContainer();
633         // Number of hydrogen in an entity from the XYZ.
634         int numMolHydrogen = 0;
635         int atomCount = 0;
636         for (Atom atom : currentAtoms) {
637           if (atom.isHydrogen()) {
638             numMolHydrogen++;
639           }
640           String atomName = getSymbol(atom.getAtomType().atomicNumber);
641           xyzCDKAtoms.addAtom(
642               new org.openscience.cdk.Atom(atomName, new Point3d(atom.getXYZ(null))));
643           atomCount++;
644         }
645         logger.info(format(" Number of atoms in XYZ CDK: %2d", atomCount));
646         int numInputMolAtoms = currentAtoms.size();
647         if (logger.isLoggable(Level.FINE)) {
648           logger.fine(format(" Current entity number of atoms: %d (%d + %dH)", numInputMolAtoms,
649                   numInputMolAtoms - numMolHydrogen, numMolHydrogen));
650         }
651 
652         // Add known input bonds; a limitation is that bonds all are given a Bond order of 1.
653         List<Bond> bonds = mol.getBondList();
654         IBond.Order order = IBond.Order.SINGLE;
655         int xyzBonds = bonds.size();
656         if (xyzBonds == 0) {
657           logger.warning(" No bonds detected in input structure. Please check input.\n " +
658               "If correct, separate non-bonded entities into multiple CIFs.");
659         }
660         for (Bond xyzBond : bonds) {
661           int atom0Index = xyzBond.getAtom(0).getXyzIndex();
662           int atom1Index = xyzBond.getAtom(1).getXyzIndex();
663           if (logger.isLoggable(Level.FINER)) {
664             logger.finer(
665                 format(" Bonded atom 1: %d, Bonded atom 2: %d", atom0Index,
666                     atom1Index));
667           }
668           logger.fine(format(" Atom 0 index: %2d Atom 1 Index: %2d shiftIndex: %2d Value 0: %2d Value 1: %2d",
669                   atom0Index, atom1Index, shiftIndex, atom0Index - shiftIndex - 1, atom1Index - shiftIndex - 1));
670           xyzCDKAtoms.addBond(atom0Index - shiftIndex - 1, atom1Index - shiftIndex - 1, order);
671         }
672 
673         // Assign CDK atom types for the input molecule.
674         AtomTypeFactory factory = AtomTypeFactory.getInstance(
675             "org/openscience/cdk/config/data/jmol_atomtypes.txt", xyzCDKAtoms.getBuilder());
676 
677         // Assign atom types to CDK object.
678         for (IAtom atom : xyzCDKAtoms.atoms()) {
679           setAtomTypes(factory, atom);
680           try {
681             factory.configure(atom);
682           } catch (Exception ex) {
683             logger.info(" Failed to configure atoms from CIF.\n" + ex + "\n" + Arrays.toString(
684                 ex.getStackTrace()));
685           }
686         }
687 
688         ArrayList<ArrayList<Integer>> zindices = new ArrayList<>();
689         int counter = 0;
690         // Bond atoms from CIF.
691         int cifBonds = bondAtoms(atoms, bondTolerance);
692         if (logger.isLoggable(Level.FINE)) {
693           logger.fine(
694               format(" Created %d bonds between CIF atoms (%d in input).", cifBonds, xyzBonds));
695         }
696         List<Atom> atomPool = new ArrayList<>(Arrays.asList(atoms));
697 
698         try {
699           while (!atomPool.isEmpty()) {
700             ArrayList<Atom> molecule = new ArrayList<>();
701             collectAtoms(atomPool.get(0), molecule);
702             if (logger.isLoggable(Level.FINER)) {
703               logger.finer(format(" Molecule (%d) Size: %d", counter, molecule.size()));
704             }
705             ArrayList<Integer> indices = new ArrayList<>();
706             while (!molecule.isEmpty()) {
707               Atom atom = molecule.remove(0);
708               indices.add(atom.getIndex());
709               atomPool.remove(atom);
710             }
711 
712             if (logger.isLoggable(Level.FINER)) {
713               logger.finer(format(
714                   " Molecule %d: %d atoms are ready and %d remain (%d atoms in input, %d atoms in CIF). ",
715                   counter + 1, indices.size(), atomPool.size(), nInputAtoms, nAtoms));
716             }
717             zindices.add(indices);
718             counter++;
719           }
720         } catch (Exception e) {
721           logger.info(e + "\n" + Arrays.toString(e.getStackTrace()));
722           logger.warning(" Failed to separate copies within the asymmetric unit.");
723           return false;
724         }
725         zPrime = zindices.size();
726         // Set up CIF contents as CDK variable
727         AtomContainer[] cifCDKAtomsArr = new AtomContainer[zPrime];
728         AtomContainer cifCDKAtoms = new AtomContainer();
729         for (int j = 0; j < zPrime; j++) {
730           if (zPrime > 1) {
731             logger.info(format("\n Attempting entity %d of %d", j + 1, zPrime));
732           }
733           ArrayList<Integer> currentList = zindices.get(j);
734           int cifMolAtoms = currentList.size();
735           if (logger.isLoggable(Level.FINE)) {
736             logger.fine(format(" CIF atoms in current: %d", cifMolAtoms));
737           }
738           // Detect if CIF contains multiple copies (Z'>1)
739           if (cifMolAtoms % numInputMolAtoms == 0
740               || cifMolAtoms % (numInputMolAtoms - numMolHydrogen) == 0) {
741             cifCDKAtomsArr[j] = new AtomContainer();
742             for (Integer integer : currentList) {
743               cifCDKAtomsArr[j].addAtom(new org.openscience.cdk.Atom(atoms[integer-1].getSegID(),
744                   new Point3d(atoms[integer - 1].getXYZ(null))));
745             }
746             AtomContainer nullCDKAtoms = new AtomContainer();
747 
748             for (IAtom atom : cifCDKAtomsArr[j].atoms()) {
749               if (atom.toString() == null) {
750                 nullCDKAtoms.addAtom(atom);
751               }
752             }
753             cifCDKAtomsArr[j].remove(nullCDKAtoms);
754 
755             // Compute bonds for CIF molecule.
756             factory = AtomTypeFactory.getInstance(
757                 "org/openscience/cdk/config/data/jmol_atomtypes.txt",
758                 cifCDKAtomsArr[j].getBuilder());
759             for (IAtom atom : cifCDKAtomsArr[j].atoms()) {
760               setAtomTypes(factory, atom);
761               try {
762                 factory.configure(atom);
763               } catch (Exception ex) {
764                 logger.info(" Failed to configure CIF atoms.\n" + ex + "\n" + Arrays.toString(
765                     ex.getStackTrace()));
766               }
767             }
768 
769             RebondTool rebonder = new RebondTool(MAX_COVALENT_RADIUS, MIN_BOND_DISTANCE,
770                 bondTolerance);
771             try {
772               rebonder.rebond(cifCDKAtomsArr[j]);
773             } catch (Exception ex) {
774               logger.info(
775                   "Failed to rebond CIF atoms.\n" + ex + "\n" + Arrays.toString(ex.getStackTrace()));
776             }
777 
778             int cifMolBonds = cifCDKAtomsArr[j].getBondCount();
779             if (logger.isLoggable(Level.FINE)) {
780               logger.fine(format(" Number of CIF bonds: %d (%d in input)", cifMolBonds, xyzBonds));
781             }
782             // Number of bonds matches.
783             // If cifMolBonds == 0 then ion or atom with implicit hydrogen (e.g. water, methane, etc.)
784             if (cifMolBonds != 0 && cifMolBonds % xyzBonds == 0) {
785               Pattern pattern = VentoFoggia.findIdentical(
786                   xyzCDKAtoms, AtomMatcher.forElement(), BondMatcher.forAny());
787               int[] p = pattern.match(cifCDKAtomsArr[j]);
788               int pLength = p.length;
789               if (p != null && pLength == numInputMolAtoms) {
790                 // Used matched atoms to update the positions of the input file atoms.
791                 for (int k = 0; k < pLength; k++) {
792                   if (logger.isLoggable(Level.FINEST)) {
793                     logger.finest(
794                         format(" %d input %s -> CIF %s", k, xyzCDKAtoms.getAtom(k).getSymbol(),
795                             cifCDKAtomsArr[j].getAtom(p[k]).getSymbol()));
796                   }
797                   Point3d point3d = cifCDKAtomsArr[j].getAtom(p[k]).getPoint3d();
798                   currentAtoms.get(k).setXYZ(new double[]{point3d.x, point3d.y, point3d.z});
799                 }
800               } else {
801                 if (logger.isLoggable(Level.FINE)) {
802                   logger.fine(
803                       format(" Atoms from CIF (%d) and input (%d) structures don't match.", p.length,
804                           nAtoms));
805                 }
806                 continue;
807               }
808             } else if ((xyzBonds - numMolHydrogen) == 0
809                 || cifMolBonds % ((xyzBonds - numMolHydrogen)) == 0) {
810               // Hydrogen most likely missing from file. If zero, then potentially water/methane (implicit hydrogen atoms).
811               if (logger.isLoggable(Level.FINE)) {
812                 logger.info(" CIF may contain implicit hydrogen -- attempting to patch.");
813               }
814               // Match heavy atoms between CIF and input
815               Pattern pattern = VentoFoggia.findSubstructure(
816                   cifCDKAtomsArr[j], AtomMatcher.forElement(), BondMatcher.forAny());
817               int[] p = pattern.match(xyzCDKAtoms);
818               int pLength = p.length;
819               if (p != null && pLength == numInputMolAtoms - numMolHydrogen) {
820                 // Used matched atoms to update the positions of the input file atoms.
821                 for (int k = 0; k < pLength; k++) {
822                   Point3d point3d = cifCDKAtomsArr[j].getAtom(k).getPoint3d();
823                   currentAtoms.get(p[k]).setXYZ(new double[]{point3d.x, point3d.y, point3d.z});
824                 }
825                 // Add in hydrogen atoms
826                 Atom lastKnownAtom1 = null;
827                 int knownCount = 0;
828                 for (Atom hydrogen : currentAtoms) {
829                   if (hydrogen.isHydrogen()) {
830                     Bond bond0 = hydrogen.getBonds().get(0);
831                     Atom atom1 = bond0.get1_2(hydrogen);
832                     double angle0_2 = 0;
833                     List<Angle> anglesList = hydrogen.getAngles(); // Same as bond
834                     Atom atom2 = null;
835                     switch (anglesList.size()) {
836                       case 0 ->
837                         // H-Cl No angles
838                         // Place hydrogen slightly off center of bonded atom (~1Å away).
839                           hydrogen.moveTo(new double[]{atom1.getX() - 0.6, atom1.getY() - 0.6,
840                               atom1.getZ() - 0.6});
841                       case 1 -> {
842                         // H-O=C Need torsion
843                         // H-O-H
844                         for (Angle angle : anglesList) {
845                           atom2 = angle.get1_3(hydrogen);
846                           if (atom2 != null) {
847                             angle0_2 = angle.angleType.angle[0];
848                           }
849                           if (angle0_2 != 0) {
850                             //if angle0_2 is found then no need to look for another atom.
851                             break;
852                           }
853                         }
854                         assert atom2 != null;
855                         List<Bond> bonds2 = atom2.getBonds();
856                         Atom atom3 =
857                             (bonds2.size() > 1 && atom1 == bonds2.get(0).get1_2(atom2)) ? bonds2.get(
858                                 1).get1_2(atom2) : bonds2.get(0).get1_2(atom2);
859                         double diAng = dihedralAngle(hydrogen.getXYZ(null), atom1.getXYZ(null),
860                             atom2.getXYZ(null), atom3.getXYZ(null));
861                         if (atom1 != atom3) {
862                           intxyz(hydrogen, atom1, bond0.bondType.distance, atom2, angle0_2, atom3,
863                               Math.toDegrees(diAng), 0);
864                         } else {
865                           // Likely water as Atom 3 is not unique. Since no hydrogen atoms
866                           //  are present, there isn't a third atom...
867                           double[] coord = new double[]{atom2.getX(), atom2.getY(), atom3.getZ()};
868                           double mag = sqrt(
869                               coord[0] * coord[0] + coord[1] * coord[1] + coord[2] * coord[2]);
870                           coord[0] /= mag;
871                           coord[1] /= mag;
872                           coord[2] /= mag;
873                           hydrogen.moveTo(atom1.getX() - coord[0], atom1.getY() - coord[1],
874                               atom1.getZ() - coord[2]);
875                         }
876                       }
877                       default -> {
878                         // H-C(-C)(-C)-C
879                         Atom atom2B = null;
880                         double angle0_2B = 0;
881                         Atom proposedAtom;
882                         double proposedAngle = 0;
883                         int chiral = 1;
884                         for (Angle angle : anglesList) {
885                           proposedAtom = angle.get1_3(hydrogen);
886                           if (proposedAtom != null && !proposedAtom.isHydrogen()) {
887                             proposedAngle = angle.angleType.angle[0];
888                           }
889                           if (proposedAngle != 0) {
890                             // If angle1_3 is found then no need to look for another atom.
891                             if (angle0_2 != 0) {
892                               atom2B = proposedAtom;
893                               angle0_2B = proposedAngle;
894                             } else {
895                               atom2 = proposedAtom;
896                               angle0_2 = proposedAngle;
897                             }
898                             proposedAngle = 0.0;
899                           }
900                           if (angle0_2 != 0 && angle0_2B != 0) {
901                             break;
902                           }
903                         }
904                         if (lastKnownAtom1 == null || lastKnownAtom1 != atom1) {
905                           lastKnownAtom1 = atom1;
906                           knownCount = 0;
907                         } else {
908                           knownCount++;
909                         }
910                         if (angle0_2B == 0.0) {
911                           // Hydrogen position depends on other hydrogen, use generic location
912                           chiral = 0;
913                           assert atom2 != null;
914                           List<Bond> bonds2 = atom2.getBonds();
915                           atom2B =
916                               (atom1 == bonds2.get(0).get1_2(atom2)) ? bonds2.get(1).get1_2(atom2)
917                                   : bonds2.get(0).get1_2(atom2);
918                           // Evenly space out hydrogen
919                           angle0_2B = 180.0 - 120.0 * knownCount;
920                         } else if (anglesList.size() == 2) {
921                           // Trigonal hydrogen use equipartition between selected atoms
922                           chiral = 3;
923                         } else if (knownCount == 1) {
924                           // Already created one hydrogen (chiral = 1), use other.
925                           chiral = -1;
926                         }
927                         //TODO discern whether to use chiral = 1 or -1 when angles are to three heavy atoms.
928                         intxyz(hydrogen, atom1, bond0.bondType.distance, atom2, angle0_2, atom2B,
929                             angle0_2B, chiral);
930                       }
931                     }
932                   }
933                 }
934               } else {
935                 if (logger.isLoggable(Level.FINE)) {
936                   logger.fine(" Could not match heavy atoms between CIF and input.");
937                 }
938                 if (p != null && logger.isLoggable(Level.FINE)) {
939                   logger.fine(
940                       format(" Matched %d atoms out of %d in CIF (%d in input)", pLength, nAtoms,
941                           nInputAtoms - numMolHydrogen));
942                 }
943                 continue;
944               }
945             } else {
946               logger.info(format(" CIF (%d) and input ([%d+%dH=]%d) have a different number of bonds.",
947                   cifMolBonds, xyzBonds - numMolHydrogen, numMolHydrogen, xyzBonds));
948               continue;
949             }
950             cifCDKAtoms.add(cifCDKAtomsArr[j]);
951             MSGroup molecule;
952             if (mol instanceof Polymer) {
953               molecule = new Polymer(((Polymer) mol).getChainID(), mol.getName(), true);
954             } else {
955               molecule = new Molecule(mol.getName());
956             }
957             List<Atom> atomList = new ArrayList<>();
958             for (Atom atom : currentAtoms) {
959               if(logger.isLoggable(Level.FINER)) {
960                 logger.finer(format(" Atom Residue: %s %2d", atom.getResidueName(), atom.getResidueNumber()));
961               }
962               Atom molAtom;
963               if (molecule instanceof Polymer) {
964                 molAtom = new Atom(atomIndex++, atom.getName(), atom.getAltLoc(), atom.getXYZ(null),
965                     atom.getResidueName(), atom.getResidueNumber(), atom.getChainID(), atom.getOccupancy(),
966                     atom.getTempFactor(), atom.getSegID());
967                 molAtom.setAtomType(atom.getAtomType());
968               } else {
969                 molAtom = new Atom(atomIndex++, atom.getName(), atom.getAtomType(), atom.getXYZ(null));
970                 if(atom.getResidueName() != null){
971                   molAtom.setResName(atom.getResidueName());
972                 }
973               }
974               atomList.add(molAtom);
975             }
976             List<Bond> bondList = mol.getBondList();
977             for (Bond bond : bondList) {
978               Atom a1 = bond.getAtom(0);
979               Atom a2 = bond.getAtom(1);
980               if(logger.isLoggable(Level.FINE)) {
981                 logger.fine(format(" Bonded atom 1: %d, Bonded atom 2: %d", a1.getXyzIndex(),
982                     a2.getXyzIndex()));
983               }
984               Atom newA1 = atomList.get(a1.getIndex() - shiftIndex - 1);
985               Atom newA2 = atomList.get(a2.getIndex() - shiftIndex - 1);
986               Bond bond2 = new Bond(newA1, newA2);
987               bond2.setBondType(bond.getBondType());
988             }
989             Residue res = null;
990             for (Atom atom : atomList) {
991               if(molecule instanceof Polymer){
992                 if(res == null){
993                   res = new Residue(atom.getResidueName(), atom.getResidueNumber(), ((Polymer) mol).getResidue(atom.getResidueNumber()).getResidueType());
994                 } else if(res.getResidueNumber() != atom.getResidueNumber()){
995                   if(logger.isLoggable(Level.FINER)) {
996                     logger.finer(" Added Residue " + res.getResidueNumber() + " " + res.getName());
997                   }
998                   molecule.addMSNode(res);
999                   res = new Residue(atom.getResidueName(), atom.getResidueNumber(), ((Polymer) mol).getResidue(atom.getResidueNumber()).getResidueType());
1000                 }
1001                 res.addMSNode(atom);
1002               }else {
1003                 molecule.addMSNode(atom);
1004               }
1005             }
1006             if(molecule instanceof Polymer && res != null){
1007               if(logger.isLoggable(Level.FINER)) {
1008                 logger.finer(" Added Final Residue " + res.getResidueNumber() + " " + res.getName());
1009               }
1010               molecule.addMSNode(res);
1011             }
1012             outputAssembly.addMSNode(molecule);
1013           } else {
1014             if (logger.isLoggable(Level.INFO)) {
1015               logger.info(
1016                   format(" Number of atoms in CIF (%d) molecule do not match input (%d + %dH = %d).",
1017                       cifMolAtoms, nInputAtoms - numMolHydrogen, numMolHydrogen, nInputAtoms));
1018             }
1019           }
1020         }
1021       }
1022       // If no atoms, then conversion has failed... use active assembly
1023       if (logger.isLoggable(Level.FINE)) {
1024         logger.fine(format("\n Output Assembly Atoms: %d", outputAssembly.getAtomList().size()));
1025       }
1026       outputAssembly.setPotential(activeMolecularAssembly.getPotentialEnergy());
1027       outputAssembly.setCrystal(crystal);
1028       outputAssembly.setForceField(activeMolecularAssembly.getForceField());
1029       outputAssembly.setFile(activeMolecularAssembly.getFile());
1030       outputAssembly.setName(activeMolecularAssembly.getName());
1031       setMolecularSystem(outputAssembly);
1032 
1033       if (outputAssembly.getAtomList().isEmpty()) {
1034         logger.info(" Atom types could not be matched. File could not be written.");
1035       } else if (!writeOutputFile()) {
1036         logger.info(" Input File could not be written.");
1037       }
1038       // Finsished with this block. Clean up/reset variables for next block.
1039       outputAssembly.destroy();
1040       sgNum = -1;
1041       sgName = null;
1042     }
1043     return true;
1044   }
1045 
1046   /**
1047    * Reads the next snapshot of an archive into the activeMolecularAssembly. After calling this
1048    * function, a BufferedReader will remain open until the <code>close</code> method is called.
1049    */
1050   @Override
1051   public boolean readNext() {
1052     return readNext(false, true);
1053   }
1054 
1055   /**
1056    * Reads the next snapshot of an archive into the activeMolecularAssembly. After calling this
1057    * function, a BufferedReader will remain open until the <code>close</code> method is called.
1058    */
1059   @Override
1060   public boolean readNext(boolean resetPosition) {
1061     return readNext(resetPosition, true);
1062   }
1063 
1064   /**
1065    * Reads the next snapshot of an archive into the activeMolecularAssembly. After calling this
1066    * function, a BufferedReader will remain open until the <code>close</code> method is called.
1067    */
1068   @Override
1069   public boolean readNext(boolean resetPosition, boolean print) {
1070     return readNext(resetPosition, print, true);
1071   }
1072 
1073   /**
1074    * Reads the next snapshot of an archive into the activeMolecularAssembly. After calling this
1075    * function, a BufferedReader will remain open until the <code>close</code> method is called.
1076    */
1077   @Override
1078   public boolean readNext(boolean resetPosition, boolean print, boolean parse) {
1079     List<CifCoreBlock> blocks = cifFile.getBlocks();
1080     CifCoreBlock currentBlock;
1081     if (!parse) {
1082       snapShot++;
1083       if (print) {
1084         logger.info(format(" Skipped Block: %d", snapShot));
1085       }
1086       return true;
1087     } else if (resetPosition) {
1088       currentBlock = blocks.get(0);
1089       snapShot = 0;
1090     } else if (++snapShot < blocks.size()) {
1091       currentBlock = blocks.get(snapShot);
1092     } else {
1093       if (print) {
1094         logger.info(" Reached end of available blocks in CIF file.");
1095       }
1096       return false;
1097     }
1098     if (print) {
1099       logger.info(" Current Block: " + currentBlock.getBlockHeader());
1100     }
1101     return true;
1102   }
1103 
1104   /**
1105    * Specify the atom types for atoms created by factory.
1106    *
1107    * @param factory Atom generator
1108    * @param atom Atom of interest
1109    */
1110   public static void setAtomTypes(AtomTypeFactory factory, IAtom atom) {
1111     String atomTypeName = atom.getAtomTypeName();
1112     if (atomTypeName == null || atomTypeName.isEmpty()) {
1113       IAtomType[] types = factory.getAtomTypes(atom.getSymbol());
1114       if (types.length > 0) {
1115         IAtomType atomType = types[0];
1116         atom.setAtomTypeName(atomType.getAtomTypeName());
1117       } else {
1118         logger.info(" No atom type found for " + atom);
1119       }
1120     }
1121   }
1122 
1123   /**
1124    * Set buffer value when bonding atoms together.
1125    *
1126    * @param bondTolerance Value added to VdW radii when determining bonding.
1127    */
1128   public void setBondTolerance(double bondTolerance) {
1129     this.bondTolerance = bondTolerance;
1130   }
1131 
1132   /**
1133    * Determine whether lattice parameters can be manipulated to follow lattice system constraints.
1134    *
1135    * @param fixLattice True if lattices should be fixed.
1136    */
1137   public void setFixLattice(boolean fixLattice) {
1138     this.fixLattice = fixLattice;
1139   }
1140 
1141   /**
1142    * Override the space group of a CIF conversion based on space group name.
1143    *
1144    * @param sgName Name of the desired space group
1145    */
1146   public void setSgName(String sgName) {
1147     this.sgName = sgName;
1148   }
1149 
1150   /**
1151    * Override the space group of a CIF conversion based on space group number.
1152    *
1153    * @param sgNum Number of the desired space group
1154    */
1155   public void setSgNum(int sgNum) {
1156     this.sgNum = sgNum;
1157   }
1158 
1159   /**
1160    * Override the number of copies in the asymmetric unit (Z').
1161    *
1162    * @param zPrime Number of the desired space group
1163    */
1164   public void setZprime(int zPrime) {
1165     this.zPrime = zPrime;
1166   }
1167 
1168   /**
1169    * Write CIF files for multiple File objects.
1170    *
1171    * @return whether conversion was successful.
1172    */
1173   public boolean writeFiles() {
1174     for (File file : files) {
1175       File saveFile = new File(removeExtension(file.getAbsolutePath()) + ".cif");
1176       if (!writeFile(saveFile, true, null)) {
1177         return false;
1178       }
1179     }
1180     return true;
1181   }
1182 
1183   /**
1184    * Save a CIF file for the given molecular assembly.
1185    *
1186    * @return whether conversion was successful.
1187    */
1188   @Override
1189   public boolean writeFile(File saveFile, boolean append) {
1190     return writeFile(saveFile, append, null);
1191   }
1192 
1193   /**
1194    * Save a CIF file for the given molecular assembly.
1195    *
1196    * @return whether conversion was successful.
1197    */
1198   @Override
1199   public boolean writeFile(File saveFile, boolean append, String[] extraLines) {
1200     try {
1201       if (!append) {
1202         SystemFilter.setVersioning(Versioning.PREFIX);
1203         saveFile = version(saveFile);
1204       }
1205       BufferedWriter bw = new BufferedWriter(new FileWriter(saveFile, append));
1206       if (extraLines != null) {
1207         for (String line : extraLines) {
1208           line = line.replaceAll("\n", " ");
1209           bw.write("### " + line + "\n");
1210         }
1211       }
1212       bw.write("\ndata_" + activeMolecularAssembly.getName());
1213       Crystal xtal = activeMolecularAssembly.getCrystal().getUnitCell();
1214       bw.write("\n_symmetry_cell_setting\t" + xtal.spaceGroup.latticeSystem.name().toLowerCase()
1215           .replaceAll("_lattice", ""));
1216       bw.write("\n_symmetry_space_group_name_H-M\t" + "'" + xtal.spaceGroup.shortName + "'");
1217       bw.write("\n_symmetry_Int_Tables_number\t" + xtal.spaceGroup.number);
1218       bw.write("\nloop_\n_symmetry_equiv_pos_site_id\n_symmetry_equiv_pos_as_xyz");
1219       int numSymOps = xtal.spaceGroup.getNumberOfSymOps();
1220       for (int i = 0; i < numSymOps; i++) {
1221         bw.write(format(
1222             "\n%d " + xtal.spaceGroup.getSymOp(i).toXYZString().toLowerCase().replaceAll(" +", ""),
1223             i));
1224       }
1225       bw.write(format("\n_cell_length_a\t%4.4f", xtal.a));
1226       bw.write(format("\n_cell_length_b\t%4.4f", xtal.b));
1227       bw.write(format("\n_cell_length_c\t%4.4f", xtal.c));
1228       bw.write(format("\n_cell_angle_alpha\t%4.4f", xtal.alpha));
1229       bw.write(format("\n_cell_angle_beta\t%4.4f", xtal.beta));
1230       bw.write(format("\n_cell_angle_gamma\t%4.4f", xtal.gamma));
1231       bw.write(format("\n_cell_volume\t%4.4f", xtal.volume));
1232       int numEntities = (zPrime < 1) ? activeMolecularAssembly.getMolecules().size() : zPrime;
1233       if (numEntities > 1) {
1234         if (zPrime < 1) {
1235           logger.info(format(" Molecules detected, guessing a Z' of %d. Set manually using --zp.",
1236               numEntities));
1237         }
1238         bw.write(format("\n_cell_formula_units_Z\t%3d", numEntities));
1239       }
1240       bw.write("\nloop_");
1241       bw.write("\n_atom_site_label");
1242       bw.write("\n_atom_site_type_symbol");
1243       bw.write("\n_atom_site_fract_x");
1244       bw.write("\n_atom_site_fract_y");
1245       bw.write("\n_atom_site_fract_z");
1246       Atom[] atoms = activeMolecularAssembly.getAtomArray();
1247       int count = 1;
1248       for (Atom atom : atoms) {
1249         String name = atom.getName();
1250         String symbol = getSymbol(atom.getAtomicNumber());
1251         if (Objects.equals(name, symbol)) {
1252           name += count++;
1253         }
1254         double[] xyzC = {atom.getX(), atom.getY(), atom.getZ()};
1255         double[] xyzF = new double[3];
1256         xtal.toFractionalCoordinates(xyzC, xyzF);
1257         bw.write(format("\n%-3s %2s %8.6f %8.6f %8.6f", name, symbol, xyzF[0], xyzF[1], xyzF[2]));
1258       }
1259       bw.write("\n#END\n");
1260       bw.close();
1261       logger.info(format("\n Wrote CIF file: %s \n", saveFile.getAbsolutePath()));
1262     } catch (Exception ex) {
1263       logger.info(format("\n Failed to write out CIF file: %s \n" + ex, saveFile.getAbsolutePath()));
1264       return false;
1265     }
1266     if (!createdFiles.contains(saveFile.getAbsolutePath())) {
1267       createdFiles.add(saveFile.getAbsolutePath());
1268     }
1269     return true;
1270   }
1271 
1272   /**
1273    * Write an output file.
1274    *
1275    * @return Whether file was written successfully.
1276    */
1277   private boolean writeOutputFile() {
1278     String dir = getFullPath(files.get(0).getAbsolutePath()) + File.separator;
1279     String fileName = removeExtension(getName(files.get(0).getAbsolutePath()));
1280     String spacegroup;
1281     if(activeMolecularAssembly.getCrystal() != null) {
1282       spacegroup = activeMolecularAssembly.getCrystal().getUnitCell().spaceGroup.shortName;
1283     }else{
1284       spacegroup = null;
1285     }
1286     List<MSNode> entities = activeMolecularAssembly.getAllBondedEntities();
1287 
1288     File saveFile;
1289     if(usePDB){
1290       // Change name for different space groups. Input cannot handle multiple space groups in same file.
1291       if (entities.size() > 1) {
1292         fileName += "_z" + entities.size();
1293       }
1294       saveFile = new File(dir + fileName + ".pdb");
1295     } else if (cifFile.getBlocks().size() > 1) {
1296       // Change name for different space groups. Input cannot handle multiple space groups in same file.
1297       if (entities.size() > 1) {
1298         fileName += "_z" + entities.size();
1299       }
1300       fileName += "_" + spacegroup.replaceAll("\\/", "");
1301       saveFile = new File(dir + fileName + ".arc");
1302     } else {
1303       // If only structure then create XYZ.
1304       saveFile = XYZFilter.version(new File(dir + fileName + ".xyz"));
1305     }
1306     ForceField ff = activeMolecularAssembly.getForceField();
1307     CompositeConfiguration properties = activeMolecularAssembly.getProperties();
1308     if(usePDB){
1309       PDBFilter pdbFilter = new PDBFilter(saveFile, activeMolecularAssembly, ff, properties);
1310       pdbFilter.writeFile(saveFile, true);
1311     }else {
1312       XYZFilter xyzFilter = new XYZFilter(saveFile, activeMolecularAssembly, ff, properties);
1313       xyzFilter.writeFile(saveFile, true);
1314     }
1315     logger.info("\n Saved output file:        " + saveFile.getAbsolutePath());
1316     writePropertiesFile(dir, fileName, spacegroup, entities);
1317     if (!createdFiles.contains(saveFile.getAbsolutePath())) {
1318       createdFiles.add(saveFile.getAbsolutePath());
1319     }
1320     return true;
1321   }
1322 
1323   /**
1324    * Write an output file.
1325    *
1326    * @return Whether file was written successfully.
1327    */
1328   public boolean writeOutputFile(MolecularAssembly ma) {
1329     setMolecularSystem(ma);
1330     return writeOutputFile();
1331   }
1332 
1333   /**
1334    * Write a properties file.
1335    * @param dir String for directory location.
1336    * @param fileName String for file name.
1337    * @param spacegroup String for space group.
1338    * @param entities List of MSNodes in assembly.
1339    */
1340   public void writePropertiesFile(String dir, String fileName, String spacegroup, List<MSNode> entities) {
1341     // Write out a property file.
1342     File propertyFile = new File(dir + fileName + ".properties");
1343     if (!propertyFile.exists()) {
1344       try {
1345         FileWriter fw = new FileWriter(propertyFile, false);
1346         BufferedWriter bw = new BufferedWriter(fw);
1347         ForceField forceField = activeMolecularAssembly.getForceField();
1348         String parameters = forceField.getString("parameters", "none");
1349         if (parameters != null && !parameters.equalsIgnoreCase("none")) {
1350           bw.write(format("parameters %s\n", parameters));
1351         }
1352         String forceFieldProperty = forceField.getString("forcefield", "none");
1353         if (forceFieldProperty != null && !forceFieldProperty.equalsIgnoreCase("none")) {
1354           bw.write(format("forcefield %s\n", forceFieldProperty));
1355         }
1356         String patch = forceField.getString("patch", "none");
1357         if (patch != null && !patch.equalsIgnoreCase("none")) {
1358           bw.write(format("patch %s\n", patch));
1359         }
1360         if(spacegroup != null) {
1361           bw.write(format("spacegroup %s\n", spacegroup));
1362         }
1363         if (entities.size() > 1) {
1364           bw.write("intermolecular-softcore true\n");
1365         }
1366         logger.info("\n Saved properties file: " + propertyFile.getAbsolutePath() + "\n");
1367         bw.close();
1368       } catch (Exception ex) {
1369         logger.info("Failed to write files.\n" + ex);
1370       }
1371     } else {
1372       logger.info("\n Property file already exists:  " + propertyFile.getAbsolutePath() + "\n");
1373     }
1374   }
1375 }