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