View Javadoc
1   // ******************************************************************************
2   //
3   // Title:       Force Field X.
4   // Description: Force Field X - Software for Molecular Biophysics.
5   // Copyright:   Copyright (c) Michael J. Schnieders 2001-2024.
6   //
7   // This file is part of Force Field X.
8   //
9   // Force Field X is free software; you can redistribute it and/or modify it
10  // under the terms of the GNU General Public License version 3 as published by
11  // the Free Software Foundation.
12  //
13  // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14  // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15  // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16  // details.
17  //
18  // You should have received a copy of the GNU General Public License along with
19  // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20  // Place, Suite 330, Boston, MA 02111-1307 USA
21  //
22  // Linking this library statically or dynamically with other modules is making a
23  // combined work based on this library. Thus, the terms and conditions of the
24  // GNU General Public License cover the whole combination.
25  //
26  // As a special exception, the copyright holders of this library give you
27  // permission to link this library with independent modules to produce an
28  // executable, regardless of the license terms of these independent modules, and
29  // to copy and distribute the resulting executable under terms of your choice,
30  // provided that you also meet, for each linked independent module, the terms
31  // and conditions of the license of that module. An independent module is a
32  // module which is not derived from or based on this library. If you modify this
33  // library, you may extend this exception to your version of the library, but
34  // you are not obligated to do so. If you do not wish to do so, delete this
35  // exception statement from your version.
36  //
37  // ******************************************************************************
38  package ffx.potential.parsers;
39  
40  import com.github.javaparser.resolution.declarations.ResolvedInterfaceDeclaration;
41  import ffx.crystal.Crystal;
42  import ffx.crystal.SymOp;
43  import ffx.potential.MolecularAssembly;
44  import ffx.potential.Utilities.FileType;
45  import ffx.potential.bonded.AminoAcidUtils;
46  import ffx.potential.bonded.Atom;
47  import ffx.potential.bonded.Bond;
48  import ffx.potential.bonded.Residue;
49  import ffx.potential.extended.ExtendedSystem;
50  import ffx.potential.parameters.AtomType;
51  import ffx.potential.parameters.BondType;
52  import ffx.potential.parameters.ForceField;
53  import org.apache.commons.configuration2.CompositeConfiguration;
54  import org.jogamp.vecmath.Vector3d;
55  
56  import java.io.*;
57  import java.lang.reflect.Array;
58  import java.util.ArrayList;
59  import java.util.HashMap;
60  import java.util.List;
61  import java.util.OptionalDouble;
62  import java.util.logging.Level;
63  import java.util.logging.Logger;
64  import java.util.regex.Matcher;
65  import java.util.regex.Pattern;
66  
67  import static ffx.potential.bonded.Bond.logNoBondType;
68  import static java.lang.Double.parseDouble;
69  import static java.lang.Integer.parseInt;
70  import static java.lang.String.format;
71  
72  /**
73   * The XYZFilter class parses TINKER Cartesian coordinate (*.XYZ) files.
74   *
75   * @author Michael J. Schnieders
76   * @since 1.0
77   */
78  public class XPHFilter extends SystemFilter {
79  
80    private static final Logger logger = Logger.getLogger(XPHFilter.class.getName());
81    private BufferedReader bufferedReader = null;
82    private int snapShot;
83    private String remarkLine;
84    private ExtendedSystem extendedSystem;
85    private MolecularAssembly system = null;
86  
87    /**
88     * Constructor for XPHFilter.
89     *
90     * @param file a {@link File} object.
91     * @param system a {@link MolecularAssembly} object.
92     * @param forceField a {@link ForceField} object.
93     * @param properties a {@link CompositeConfiguration} object.
94     */
95    public XPHFilter(File file, MolecularAssembly system, ForceField forceField,
96                     CompositeConfiguration properties, ExtendedSystem esvSystem) {
97      super(file, system, forceField, properties);
98      this.fileType = FileType.XPH;
99      extendedSystem = esvSystem;
100     this.system = system;
101   }
102 
103   public XPHFilter(SystemFilter systemFilter, ExtendedSystem esvSystem) {
104     super(systemFilter.getFile(), systemFilter.getActiveMolecularSystem(),
105             systemFilter.getActiveMolecularSystem().getForceField(),
106             systemFilter.getActiveMolecularSystem().getProperties());
107 
108     extendedSystem = esvSystem;
109   }
110 
111   /**
112    * readOnto
113    *
114    * @param newFile a {@link File} object.
115    * @param oldSystem a {@link MolecularAssembly} object.
116    * @return a boolean.
117    */
118   public static boolean readOnto(File newFile, MolecularAssembly oldSystem) {
119     if (newFile == null || !newFile.exists() || oldSystem == null) {
120       return false;
121     }
122     try (BufferedReader br = new BufferedReader(new FileReader(newFile))) {
123       String data = br.readLine();
124       if (data == null) {
125         return false;
126       }
127       String[] tokens = data.trim().split(" +");
128       int num_atoms = parseInt(tokens[0]);
129       if (num_atoms != oldSystem.getAtomList().size()) {
130         return false;
131       }
132 
133       br.mark(10000);
134       data = br.readLine();
135       if (!readPBC(data, oldSystem)) {
136         br.reset();
137       }
138 
139       double[][] d = new double[num_atoms][3];
140       for (int i = 0; i < num_atoms; i++) {
141         if (!br.ready()) {
142           return false;
143         }
144         data = br.readLine();
145         if (data == null) {
146           logger.warning(format(" Check atom %d.", (i + 1)));
147           return false;
148         }
149         tokens = data.trim().split(" +");
150         if (tokens.length < 6) {
151           logger.warning(format(" Check atom %d.", (i + 1)));
152           return false;
153         }
154         d[i][0] = parseDouble(tokens[2]);
155         d[i][1] = parseDouble(tokens[3]);
156         d[i][2] = parseDouble(tokens[4]);
157       }
158       List<Atom> atoms = oldSystem.getAtomList();
159       for (Atom a : atoms) {
160         int index = a.getIndex() - 1;
161         a.setXYZ(d[index]);
162       }
163       oldSystem.center();
164       oldSystem.setFile(newFile);
165       return true;
166     } catch (Exception e) {
167       return false;
168     }
169   }
170 
171   private static boolean firstTokenIsInteger(String data) {
172     if (data == null) {
173       return false;
174     }
175 
176     // Check for a blank line.
177     data = data.trim();
178     if (data.equals("")) {
179       return false;
180     }
181 
182     // Check if the first token in an integer.
183     try {
184       String[] tokens = data.split(" +");
185       parseInt(tokens[0]);
186       return true;
187     } catch (NumberFormatException e) {
188       return false;
189     }
190   }
191 
192   /**
193    * Attempt to parse the String as unit cell parameters.
194    *
195    * @param data The String to parse.
196    * @return false if the first token in the String is an integer and true otherwise.
197    */
198   private static boolean readPBC(String data, MolecularAssembly activeMolecularAssembly) {
199     if (firstTokenIsInteger(data)) {
200       return false;
201     }
202 
203     String[] tokens = data.trim().split(" +");
204     if (tokens.length == 6) {
205       CompositeConfiguration config = activeMolecularAssembly.getProperties();
206       double a = parseDouble(tokens[0]);
207       double b = parseDouble(tokens[1]);
208       double c = parseDouble(tokens[2]);
209       double alpha = parseDouble(tokens[3]);
210       double beta = parseDouble(tokens[4]);
211       double gamma = parseDouble(tokens[5]);
212       config.setProperty("a-axis", a);
213       config.setProperty("b-axis", b);
214       config.setProperty("c-axis", c);
215       config.setProperty("alpha", alpha);
216       config.setProperty("beta", beta);
217       config.setProperty("gamma", gamma);
218 
219       Crystal crystal = activeMolecularAssembly.getCrystal();
220       if (crystal != null) {
221         crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
222       }
223     }
224     return true;
225   }
226 
227   /** {@inheritDoc} */
228   @Override
229   public void closeReader() {
230     if (bufferedReader != null) {
231       try {
232         bufferedReader.close();
233       } catch (IOException ex) {
234         logger.warning(format(" Exception in closing XYZ filter: %s", ex));
235       }
236     }
237   }
238 
239   @Override
240   public int countNumModels(){
241     File xphFile = activeMolecularAssembly.getFile();
242     int nAtoms = activeMolecularAssembly.getAtomArray().length;
243     Pattern crystInfoPattern = Pattern.compile(
244             "^ *(?:[0-9]+\\.[0-9]+ +){3}(?:-?[0-9]+\\.[0-9]+ +){2}(?:-?[0-9]+\\.[0-9]+) *$");
245 
246     try (BufferedReader br = new BufferedReader(new FileReader(xphFile))) {
247       String line = br.readLine();
248       int nSnaps = 0;
249       // For each header line, will read either X or X+1 lines, where X is the number of atoms.
250       while (line != null) {
251         assert parseInt(line.trim().split("\\s+")[0]) == nAtoms;
252         // Read either the crystal information *or* the first line of the snapshot.
253         line = br.readLine();
254         Matcher m = crystInfoPattern.matcher(line);
255         if (m.matches()) {
256           // If this is crystal information, move onto the first line of the snapshot.
257           br.readLine();
258         }
259         // Read lines 2-X of the XYZ.
260         for (int i = 1; i < nAtoms; i++) {
261           br.readLine();
262         }
263 
264         String data = br.readLine();
265 
266         //Read past blanklines
267         while (data != null && data.trim().equals("")) {
268           data = br.readLine();
269         }
270 
271         // Read Past ESV
272 
273         assert data != null;
274         if(data.contains("ESV")) {
275 
276           while (data != null && !data.trim().equals("")) {
277             data = br.readLine();
278           }
279         }
280 
281         ++nSnaps;
282         line = br.readLine();
283       }
284       return nSnaps;
285     } catch (Exception ex) {
286       logger.log(Level.WARNING,
287               String.format(" Exception reading trajectory file %s: %s", xphFile, ex));
288       return 1;
289     }
290   }
291 
292   /** {@inheritDoc} */
293   @Override
294   public OptionalDouble getLastReadLambda() {
295     String[] toks = remarkLine.split("\\s+");
296     int nToks = toks.length;
297     for (int i = 0; i < (nToks - 1); i++) {
298       if (toks[i].equals("Lambda:")) {
299         return OptionalDouble.of(Double.parseDouble(toks[i + 1]));
300       }
301     }
302     return OptionalDouble.empty();
303   }
304 
305   @Override
306   public String[] getRemarkLines() {
307     return new String[] {remarkLine};
308   }
309 
310   @Override
311   public int getSnapshot() {
312     return snapShot;
313   }
314 
315   /**
316    * {@inheritDoc}
317    *
318    * <p>Parse the XYZ File
319    */
320   @Override
321   public boolean readFile() {
322     File xphFile = activeMolecularAssembly.getFile();
323 
324     if (forceField == null) {
325       logger.warning(format(" No force field is associated with %s.", xphFile.toString()));
326       return false;
327     }
328     try (BufferedReader br = new BufferedReader(new FileReader(xphFile))) {
329       String data = br.readLine();
330       // Read blank lines at the top of the file
331       while (data != null && data.trim().equals("")) {
332         data = br.readLine();
333       }
334       if (data == null) {
335         return false;
336       }
337       String[] tokens = data.trim().split(" +", 2);
338       int numberOfAtoms = parseInt(tokens[0]);
339       if (numberOfAtoms < 1) {
340         return false;
341       }
342       if (tokens.length == 2) {
343         getActiveMolecularSystem().setName(tokens[1]);
344       }
345       logger.info(format(" Opening %s with %d atoms\n", xphFile.getName(), numberOfAtoms));
346       remarkLine = data.trim();
347 
348       // The header line is reasonable. Check for periodic box dimensions.
349       br.mark(10000);
350       data = br.readLine();
351       if (!readPBC(data, activeMolecularAssembly)) {
352         br.reset();
353       }
354 
355       // Prepare to parse atom lines.
356       HashMap<Integer, Integer> labelHash = new HashMap<>();
357       int[] label = new int[numberOfAtoms];
358       int[][] bonds = new int[numberOfAtoms][8];
359       double[][] d = new double[numberOfAtoms][3];
360       boolean renumber = false;
361       atomList = new ArrayList<>();
362       // Loop over the expected number of atoms.
363       for (int i = 0; i < numberOfAtoms; i++) {
364         if (!br.ready()) {
365           return false;
366         }
367         data = br.readLine();
368         if (data == null) {
369           logger.warning(
370                   format(" Check atom %d in %s.", (i + 1), activeMolecularAssembly.getFile().getName()));
371           return false;
372         }
373         tokens = data.trim().split(" +");
374         if (tokens.length < 6) {
375           logger.warning(
376                   format(" Check atom %d in %s.", (i + 1), activeMolecularAssembly.getFile().getName()));
377           return false;
378         }
379         // Valid number of tokens, so try to parse this line.
380         label[i] = parseInt(tokens[0]);
381         // Check for valid atom numbering, or flag for re-numbering.
382         if (label[i] != i + 1) {
383           renumber = true;
384         }
385         String atomName = tokens[1];
386         d[i][0] = parseDouble(tokens[2]);
387         d[i][1] = parseDouble(tokens[3]);
388         d[i][2] = parseDouble(tokens[4]);
389         int type = parseInt(tokens[5]);
390         AtomType atomType = forceField.getAtomType(Integer.toString(type));
391         if (atomType == null) {
392           StringBuilder message = new StringBuilder("Check atom type ");
393           message.append(type).append(" for Atom ").append(i + 1);
394           message.append(" in ").append(activeMolecularAssembly.getFile().getName());
395           logger.warning(message.toString());
396           return false;
397         }
398         Atom a = new Atom(i + 1, atomName, atomType, d[i]);
399         atomList.add(a);
400         // Bond Data
401         int numberOfBonds = tokens.length - 6;
402         for (int b = 0; b < 8; b++) {
403           if (b < numberOfBonds) {
404             int bond = parseInt(tokens[6 + b]);
405             bonds[i][b] = bond;
406           } else {
407             bonds[i][b] = 0;
408           }
409         }
410       }
411 
412       // Check if this is an archive.
413       if (br.ready()) {
414         // Read past blank lines between archive files
415         data = br.readLine().trim();
416         while (data.equals("") && br.ready()) {
417           data = br.readLine().trim();
418         }
419         tokens = data.split(" +", 2);
420         if (tokens.length > 0) {
421           try {
422             int archiveNumberOfAtoms = parseInt(tokens[0]);
423             if (archiveNumberOfAtoms == numberOfAtoms) {
424               setType(FileType.ARC);
425             }
426           } catch (NumberFormatException e) {
427             //
428           }
429         }
430       }
431       // Try to renumber
432       if (renumber) {
433         for (int i = 0; i < numberOfAtoms; i++) {
434           if (labelHash.containsKey(label[i])) {
435             logger.warning(format(" Two atoms have the same index: %d.", label[i]));
436             return false;
437           }
438           labelHash.put(label[i], i + 1);
439         }
440         for (int i = 0; i < numberOfAtoms; i++) {
441           int j = -1;
442           while (j < 3 && bonds[i][++j] > 0) {
443             bonds[i][j] = labelHash.get(bonds[i][j]);
444           }
445         }
446       }
447       bondList = new ArrayList<>();
448       int[] c = new int[2];
449       for (int a1 = 1; a1 <= numberOfAtoms; a1++) {
450         int j = -1;
451         while (j < 7 && bonds[a1 - 1][++j] > 0) {
452           int a2 = bonds[a1 - 1][j];
453           if (a1 < a2) {
454             if (a2 > numberOfAtoms) {
455               logger.warning(format(" Check the bond between %d and %d in %s.", a1, a2,
456                       activeMolecularAssembly.getFile().getName()));
457               return false;
458             }
459             // Check for bidirectional connection
460             boolean bidirectional = false;
461             int k = -1;
462             while (k < 7 && bonds[a2 - 1][++k] > 0) {
463               int a3 = bonds[a2 - 1][k];
464               if (a3 == a1) {
465                 bidirectional = true;
466                 break;
467               }
468             }
469             if (!bidirectional) {
470               logger.warning(format(" Check the bond between %d and %d in %s.", a1, a2,
471                       activeMolecularAssembly.getFile().getName()));
472               return false;
473             }
474             Atom atom1 = atomList.get(a1 - 1);
475             Atom atom2 = atomList.get(a2 - 1);
476             if (atom1 == null || atom2 == null) {
477               logger.warning(format(" Check the bond between %d and %d in %s.", a1, a2,
478                       activeMolecularAssembly.getFile().getName()));
479               return false;
480             }
481             Bond bond = new Bond(atom1, atom2);
482             BondType bondType = forceField.getBondType(atom1.getAtomType(), atom2.getAtomType());
483             if (bondType == null) {
484               logNoBondType(atom1, atom2, forceField);
485             } else {
486               bond.setBondType(bondType);
487             }
488             bondList.add(bond);
489           }
490         }
491       }
492       // Read ESVs
493       while (data != null && data.equals("") && br.ready()) {
494         data = br.readLine().trim();
495       }
496 
497       if (data != null) {
498         tokens = data.split(" +", 2);
499 
500         if (tokens[0].equalsIgnoreCase("ESV")) {
501           int numOfESVs = parseInt(tokens[1]);
502           data = br.readLine().trim();
503 
504           List<Residue> residueList = extendedSystem.getExtendedResidueList();
505 
506           if (numOfESVs == residueList.size()) {
507             int switchIndex = extendedSystem.getTitratingResidueList().size();
508             for (int i = 0; i < residueList.size(); i++) {
509               tokens = data.split(" +");
510 
511               if (i < switchIndex) {
512                 extendedSystem.setTitrationLambda(residueList.get(i), parseDouble(tokens[2]));
513               } else {
514                 extendedSystem.setTautomerLambda(residueList.get(i), parseDouble(tokens[2]));
515               }
516 
517               data = br.readLine().trim();
518             }
519 
520             for (Atom atom : system.getAtomList()) {
521               int atomIndex = atom.getIndex() - 1;
522               atom.setOccupancy(extendedSystem.getTitrationLambda(atomIndex));
523               atom.setTempFactor(extendedSystem.getTautomerLambda(atomIndex));
524             }
525 
526           } else {
527             logger.severe(
528                     " Number of ESVs in archive doesn't match extended system residue list size.");
529             return false;
530           }
531         }
532       }
533       return true;
534     } catch (IOException e) {
535       logger.severe(e.toString());
536     }
537     return false;
538   }
539 
540   /** {@inheritDoc} */
541   @Override
542   public boolean readNext() {
543     return readNext(false);
544   }
545 
546   /**
547    * {@inheritDoc}
548    *
549    * <p>Reads the next snap-shot of an archive into the activeMolecularAssembly. After calling this
550    * function, a BufferedReader will remain open until the <code>close</code> method is called.
551    */
552   @Override
553   public boolean readNext(boolean resetPosition) {
554     return readNext(resetPosition, true);
555   }
556 
557   /**
558    * {@inheritDoc}
559    *
560    * <p>Reads the next snap-shot of an archive into the activeMolecularAssembly. After calling this
561    * function, a BufferedReader will remain open until the <code>close</code> method is called.
562    */
563   @Override
564   public boolean readNext(boolean resetPosition, boolean print) {
565     return readNext(resetPosition, print, true);
566   }
567 
568   /**
569    * Reads the next snap-shot of an archive into the activeMolecularAssembly. After calling this
570    * function, a BufferedReader will remain open until the <code>close</code> method is called.
571    */
572   public boolean readNext(boolean resetPosition, boolean print, boolean parse) {
573     try {
574       String data;
575       Atom[] atoms = activeMolecularAssembly.getAtomArray();
576       int nSystem = atoms.length;
577 
578       if (bufferedReader == null && !resetPosition) {
579         bufferedReader = new BufferedReader(new FileReader(currentFile));
580         // Read past the first N + 1 lines that begin with an integer.
581         for (int i = 0; i < nSystem + 1; i++) {
582           data = bufferedReader.readLine();
583           while (!firstTokenIsInteger(data)) {
584             data = bufferedReader.readLine();
585           }
586         }
587         snapShot = 1;
588       } else if (resetPosition) {
589         // Reset the reader to the beginning of the file. Do not skip reading the first entry if resetPostion is true.
590         bufferedReader = new BufferedReader(new FileReader(currentFile));
591         snapShot = 0;
592       }
593 
594       snapShot++;
595 
596       data = bufferedReader.readLine();
597 
598       // Read past blank lines
599       while (data != null && data.trim().equals("")) {
600         data = bufferedReader.readLine();
601       }
602       if (data == null) {
603         return false;
604       }
605 
606       // Read Past ESV
607       if (data.contains("ESV")) {
608         while (data != null && !data.trim().equals("")) {
609           data = bufferedReader.readLine();
610         }
611 
612         data = bufferedReader.readLine();
613 
614         if (data == null) {
615           return false;
616         }
617       }
618 
619       if (print) {
620         logger.info(format("\n Attempting to read snapshot %d.", snapShot));
621       }
622       try {
623         int nArchive = parseInt(data.trim().split(" +")[0]);
624         if (nArchive != nSystem) {
625           String message = format("Number of atoms mismatch (Archive: %d, System: %d).", nArchive,
626                   nSystem);
627           if (dieOnMissingAtom) {
628             logger.severe(message);
629           }
630           logger.warning(message);
631           return false;
632         }
633       } catch (NumberFormatException e) {
634         logger.warning(e.toString());
635         return false;
636       }
637 
638       remarkLine = data;
639 
640       // The header line is reasonable. Check for periodic box dimensions.
641       bufferedReader.mark(10000);
642       data = bufferedReader.readLine();
643       if (!readPBC(data, activeMolecularAssembly)) {
644         bufferedReader.reset();
645       }
646 
647       String[] tokens;
648       for (int i = 0; i < nSystem; i++) {
649         data = bufferedReader.readLine();
650         // Read past blank lines
651         while (data != null && data.trim().equals("")) {
652           data = bufferedReader.readLine();
653         }
654         tokens = data.trim().split(" +");
655         if (tokens.length < 6) {
656           String message = format("Check atom %d in %s.", (i + 1), currentFile.getName());
657           logger.warning(message);
658           return false;
659         }
660         double x = parseDouble(tokens[2]);
661         double y = parseDouble(tokens[3]);
662         double z = parseDouble(tokens[4]);
663         int xyzIndex = atoms[i].getIndex();
664         if (xyzIndex != i + 1) {
665           String message = format("Archive atom index %d being read onto system atom index %d.",
666                   i + 1, xyzIndex);
667           logger.warning(message);
668         }
669         atoms[i].moveTo(x, y, z);
670       }
671 
672       // Read ESVs
673       int counter = 0;
674       while (data != null && !data.contains("ESV") && bufferedReader.ready()) {
675         data = bufferedReader.readLine().trim();
676         counter++;
677         if (counter > 5) {
678           logger.severe(" Read through too many lines");
679         }
680       }
681 
682       if (data != null) {
683         tokens = data.split(" +", 2);
684         if (tokens[0].equalsIgnoreCase("ESV")) {
685           int numOfESVs = parseInt(tokens[1]);
686           data = bufferedReader.readLine().trim();
687 
688           List<Residue> residueList = extendedSystem.getExtendedResidueList();
689 
690           if (numOfESVs == residueList.size()) {
691             int switchIndex = extendedSystem.getTitratingResidueList().size();
692             for (int i = 0; i < residueList.size(); i++) {
693               tokens = data.split(" +");
694 
695               if (i < switchIndex) {
696                 extendedSystem.setTitrationLambda(residueList.get(i), parseDouble(tokens[2]));
697               } else {
698                 extendedSystem.setTautomerLambda(residueList.get(i), parseDouble(tokens[2]));
699               }
700 
701               data = bufferedReader.readLine().trim();
702             }
703           } else {
704             logger.severe(
705                     " Number of ESVs in archive doesn't match extended system residue list size.");
706             return false;
707           }
708         }
709       }
710       return true;
711 
712     } catch (FileNotFoundException e) {
713       String message = format("Exception opening file %s.", currentFile);
714       logger.log(Level.WARNING, message, e);
715     } catch (IOException e) {
716       String message = format("Exception reading from file %s.", currentFile);
717       logger.log(Level.WARNING, message, e);
718     }
719     return false;
720   }
721 
722   public ExtendedSystem getExtendedSystem() {
723     return extendedSystem;
724   }
725 
726   public void setExtendedSystem(ExtendedSystem esvSystem) {
727     this.extendedSystem = esvSystem;
728   }
729 
730   /** {@inheritDoc} */
731   @Override
732   public boolean writeFile(File saveFile, boolean append, String[] extraLines) {
733     if (saveFile == null) {
734       return false;
735     }
736 
737     File newFile = saveFile;
738     if (!append) {
739       newFile = version(saveFile);
740     }
741     activeMolecularAssembly.setFile(newFile);
742     if (activeMolecularAssembly.getName() == null) {
743       activeMolecularAssembly.setName(newFile.getName());
744     }
745 
746     try (FileWriter fw = new FileWriter(newFile,
747             append && newFile.exists()); BufferedWriter bw = new BufferedWriter(fw)) {
748       // XYZ File First Line
749       int numberOfAtoms = activeMolecularAssembly.getAtomList().size();
750       StringBuilder sb = new StringBuilder(
751               format("%7d  %s", numberOfAtoms, activeMolecularAssembly.getName()));
752       if (extraLines != null) {
753         for (String line : extraLines) {
754           line = line.replaceAll("\n", " ");
755           sb.append(" ").append(line);
756         }
757       }
758       String output = sb.append("\n").toString();
759       bw.write(output);
760 
761       Crystal crystal = activeMolecularAssembly.getCrystal();
762       if (crystal != null && !crystal.aperiodic()) {
763         Crystal uc = crystal.getUnitCell();
764         String params = format("%14.8f%14.8f%14.8f%14.8f%14.8f%14.8f\n", uc.a, uc.b, uc.c, uc.alpha,
765                 uc.beta, uc.gamma);
766         bw.write(params);
767       }
768 
769       Atom a2;
770       StringBuilder line;
771       StringBuilder[] lines = new StringBuilder[numberOfAtoms];
772       // XYZ File Atom Lines
773       List<Atom> atoms = activeMolecularAssembly.getAtomList();
774       Vector3d offset = activeMolecularAssembly.getOffset();
775       for (Atom a : atoms) {
776         if (vdwH) {
777           line = new StringBuilder(
778                   format("%7d %3s%14.8f%14.8f%14.8f%6d", a.getIndex(), a.getAtomType().name,
779                           a.getRedX() - offset.x, a.getRedY() - offset.y, a.getRedZ() - offset.z,
780                           a.getType()));
781         } else {
782           line = new StringBuilder(
783                   format("%7d %3s%14.8f%14.8f%14.8f%6d", a.getIndex(), a.getAtomType().name,
784                           a.getX() - offset.x, a.getY() - offset.y, a.getZ() - offset.z, a.getType()));
785         }
786         for (Bond b : a.getBonds()) {
787           a2 = b.get1_2(a);
788           line.append(format("%8d", a2.getIndex()));
789         }
790         lines[a.getIndex() - 1] = line.append("\n");
791       }
792 
793       StringBuilder[] esvLines = null;
794       if (extendedSystem != null) {
795         List<Residue> residueList = extendedSystem.getTitratingResidueList();
796         esvLines = new StringBuilder[extendedSystem.getExtendedResidueList().size()];
797         double ESV;
798 
799         List<Atom> residueAtoms;
800         int CAIndex = -1;
801 
802         // Titrating residues
803         for (int i = 0; i < residueList.size(); i++) {
804           ESV = extendedSystem.getTitrationLambda(residueList.get(i));
805 
806           residueAtoms = residueList.get(i).getAtomList(true);
807           for (Atom atom : residueAtoms) {
808             if (atom.getAtomType().name.equalsIgnoreCase("CA")) {
809               CAIndex = atom.getIndex();
810             }
811           }
812 
813           line = new StringBuilder(
814                   format("%7d%7d%14.8f%7s\n", i, CAIndex, ESV, residueList.get(i).getAminoAcid3()));
815 
816           esvLines[i] = line;
817 
818           CAIndex = -1;
819         }
820 
821         int offsetIndex = residueList.size();
822         residueList = extendedSystem.getTautomerizingResidueList();
823         // Tautomerizing residues
824         for (int i = 0; i < residueList.size(); i++) {
825           ESV = extendedSystem.getTautomerLambda(residueList.get(i));
826 
827           residueAtoms = residueList.get(i).getAtomList(true);
828           for (Atom atom : residueAtoms) {
829             if (atom.getAtomType().name.equalsIgnoreCase("CA")) {
830               CAIndex = atom.getIndex();
831             }
832           }
833 
834           line = new StringBuilder(format("%7d%7d%14.8f%7s\n", offsetIndex + i, CAIndex, ESV,
835                   residueList.get(i).getAminoAcid3()));
836 
837           esvLines[offsetIndex + i] = line;
838 
839           CAIndex = -1;
840         }
841 
842       }
843       try {
844         for (int i = 0; i < numberOfAtoms; i++) {
845           bw.write(lines[i].toString());
846         }
847 
848         if (esvLines != null) {
849           String esvHead = format("\n%7s%7d\n", "ESV", esvLines.length);
850           bw.write(esvHead);
851 
852           for (int i = 0; i < esvLines.length; i++) {
853             bw.write(esvLines[i].toString());
854           }
855           bw.write("\n");
856         }
857       } catch (IOException e) {
858         String message = format(" There was an unexpected error writing to %s.",
859                 getActiveMolecularSystem().toString());
860         logger.log(Level.WARNING, message, e);
861         return false;
862       }
863     } catch (IOException e) {
864       String message = format(" There was an unexpected error writing to %s.",
865               getActiveMolecularSystem().toString());
866       logger.log(Level.WARNING, message, e);
867       return false;
868     }
869     return true;
870   }
871 
872   /**
873    * writeFileAsP1
874    *
875    * @param saveFile a {@link File} object.
876    * @param append a boolean.
877    * @param crystal a {@link Crystal} object.
878    * @return a boolean.
879    */
880   public boolean writeFileAsP1(File saveFile, boolean append, Crystal crystal) {
881     return writeFileAsP1(saveFile, append, crystal, null);
882   }
883 
884   /**
885    * writeFileAsP1
886    *
887    * @param saveFile a {@link File} object.
888    * @param append a boolean.
889    * @param crystal a {@link Crystal} object.
890    * @param extraLines Additional lines to print in the header.
891    * @return a boolean.
892    */
893   public boolean writeFileAsP1(File saveFile, boolean append, Crystal crystal, String[] extraLines) {
894     if (saveFile == null) {
895       return false;
896     }
897 
898     File newFile = saveFile;
899     if (!append) {
900       newFile = version(saveFile);
901     }
902     activeMolecularAssembly.setFile(newFile);
903     activeMolecularAssembly.setName(newFile.getName());
904 
905     try (FileWriter fw = new FileWriter(newFile,
906             append && newFile.exists()); BufferedWriter bw = new BufferedWriter(fw)) {
907       int nSymm = crystal.spaceGroup.symOps.size();
908       // XYZ File First Line
909       int numberOfAtoms = activeMolecularAssembly.getAtomList().size() * nSymm;
910       StringBuilder sb = new StringBuilder(
911               format("%7d  %s", numberOfAtoms, activeMolecularAssembly.toString()));
912       if (extraLines != null) {
913         for (String line : extraLines) {
914           line = line.replaceAll("\n", " ");
915           sb.append(" ").append(line);
916         }
917       }
918       String output = sb.append("\n").toString();
919       bw.write(output);
920 
921       if (!crystal.aperiodic()) {
922         Crystal uc = crystal.getUnitCell();
923         String params = format("%14.8f%14.8f%14.8f%14.8f%14.8f%14.8f\n", uc.a, uc.b, uc.c, uc.alpha,
924                 uc.beta, uc.gamma);
925         bw.write(params);
926       }
927 
928       Atom a2;
929       StringBuilder line;
930       StringBuilder[] lines = new StringBuilder[numberOfAtoms];
931       // XYZ File Atom Lines
932       Atom[] atoms = activeMolecularAssembly.getAtomArray();
933       double[] xyz = new double[3];
934       for (int iSym = 0; iSym < nSymm; iSym++) {
935         SymOp symOp = crystal.spaceGroup.getSymOp(iSym);
936         int indexOffset = iSym * atoms.length;
937         for (Atom a : atoms) {
938           int index = a.getIndex() + indexOffset;
939           String id = a.getAtomType().name;
940           if (vdwH) {
941             a.getRedXYZ(xyz);
942           } else {
943             xyz[0] = a.getX();
944             xyz[1] = a.getY();
945             xyz[2] = a.getZ();
946           }
947           crystal.applySymOp(xyz, xyz, symOp);
948           int type = a.getType();
949           line = new StringBuilder(
950                   format("%7d %3s%14.8f%14.8f%14.8f%6d", index, id, xyz[0], xyz[1], xyz[2], type));
951           for (Bond b : a.getBonds()) {
952             a2 = b.get1_2(a);
953             line.append(format("%8d", a2.getIndex() + indexOffset));
954           }
955           lines[index - 1] = line.append("\n");
956         }
957       }
958       try {
959         for (int i = 0; i < numberOfAtoms; i++) {
960           bw.write(lines[i].toString());
961         }
962       } catch (IOException e) {
963         String message = format(" There was an unexpected error writing to %s.",
964                 getActiveMolecularSystem().toString());
965         logger.log(Level.WARNING, message, e);
966         return false;
967       }
968     } catch (IOException e) {
969       String message = format(" There was an unexpected error writing to %s.",
970               getActiveMolecularSystem().toString());
971       logger.log(Level.WARNING, message, e);
972       return false;
973     }
974     return true;
975   }
976 }