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