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 static ffx.potential.bonded.BondedUtils.numberAtoms;
41 import static ffx.potential.bonded.NamingUtils.renameAtomsToPDBStandard;
42 import static ffx.potential.bonded.PolymerUtils.assignAtomTypes;
43 import static ffx.potential.bonded.PolymerUtils.buildDisulfideBonds;
44 import static ffx.potential.bonded.PolymerUtils.buildMissingResidues;
45 import static ffx.potential.bonded.PolymerUtils.locateDisulfideBonds;
46 import static ffx.potential.parsers.PDBFilter.PDBFileStandard.VERSION3_2;
47 import static ffx.potential.parsers.PDBFilter.PDBFileStandard.VERSION3_3;
48 import static ffx.utilities.StringUtils.padLeft;
49 import static java.lang.Double.parseDouble;
50 import static java.lang.Integer.parseInt;
51 import static java.lang.String.format;
52 import static org.apache.commons.lang3.StringUtils.repeat;
53 import static org.apache.commons.math3.util.FastMath.min;
54
55 import ffx.crystal.Crystal;
56 import ffx.crystal.SpaceGroup;
57 import ffx.crystal.SpaceGroupDefinitions;
58 import ffx.crystal.SpaceGroupInfo;
59 import ffx.crystal.SymOp;
60 import ffx.potential.MolecularAssembly;
61 import ffx.potential.Utilities.FileType;
62 import ffx.potential.bonded.*;
63 import ffx.potential.bonded.AminoAcidUtils.AminoAcid3;
64 import ffx.potential.parameters.ForceField;
65 import ffx.utilities.Hybrid36;
66 import ffx.utilities.StringUtils;
67 import java.io.BufferedReader;
68 import java.io.BufferedWriter;
69 import java.io.File;
70 import java.io.FileNotFoundException;
71 import java.io.FileReader;
72 import java.io.FileWriter;
73 import java.io.IOException;
74 import java.util.ArrayList;
75 import java.util.Collections;
76 import java.util.HashMap;
77 import java.util.List;
78 import java.util.Map;
79 import java.util.Objects;
80 import java.util.OptionalDouble;
81 import java.util.Set;
82 import java.util.logging.Level;
83 import java.util.logging.Logger;
84 import java.util.regex.Matcher;
85 import java.util.regex.Pattern;
86 import java.util.stream.Collectors;
87 import org.apache.commons.configuration2.CompositeConfiguration;
88
89 /**
90 * The PDBFilter class parses data from a Protein DataBank (*.PDB) file. The following records are
91 * recognized: ANISOU, ATOM, CONECT, CRYST1, END, HELIX, HETATM, LINK, SHEET, SSBOND, REMARK. The
92 * rest are currently ignored.
93 *
94 * @author Michael J. Schnieders
95 * @see <a href="http://www.wwpdb.org/documentation/format32/v3.2.html">PDB format 3.2</a>
96 * @since 1.0
97 */
98 public final class PDBFilter extends SystemFilter {
99
100 private static final Logger logger = Logger.getLogger(PDBFilter.class.getName());
101 private static final Set<String> backboneNames;
102 private static final Set<String> constantPhBackboneNames;
103 private static final Set<String> naBackboneNames;
104
105 static {
106 String[] names = {"C", "CA", "N", "O", "OXT", "OT2"};
107 backboneNames = Set.of(names);
108
109 String[] constantPhNames = {"C", "CA", "N", "O", "OXT", "OT2", "H", "HA", "H1", "H2", "H3"};
110 constantPhBackboneNames = Set.of(constantPhNames);
111
112 String[] naNames = {"P", "OP1", "OP2", "O5'", "C5'", "C4'", "O4'", "C3'", "O3'", "C2'", "C1'"};
113 naBackboneNames = Set.of(naNames);
114 }
115
116 /** Map of SEQRES entries. */
117 private final Map<Character, String[]> seqRes = new HashMap<>();
118 /** Map of DBREF entries. */
119 private final Map<Character, int[]> dbRef = new HashMap<>();
120 /** List of altLoc characters seen in the PDB file. */
121 private final List<Character> altLocs = new ArrayList<>();
122 /**
123 * List of segIDs defined for the PDB file.
124 *
125 * <p>The expectation is for chain naming from A-Z, then from 0-9. For large systems, chain names
126 * are sometimes reused due to limitations in the PDB format.
127 *
128 * <p>However, we define segIDs to always be unique. For the first A-Z,0-9 series chainID ==
129 * segID. Then, for second A-Z,0-9 series, the segID = 1A-1Z,10-19, and for the third series segID
130 * = 2A-2Z,20-29, and so on.
131 */
132 private final List<String> segIDs = new ArrayList<>();
133
134 private final Map<Character, List<String>> segidMap = new HashMap<>();
135 /** Maps a chain to the number of insertion codes encountered in that chain. */
136 private final Map<Character, Integer> insertionCodeCount = new HashMap<>();
137 /**
138 * Maps chainIDResNumInsCode to renumbered chainIDResNum. For example, residue 52A in chain C might
139 * be renumbered to residue 53, and mapped as "C52A" to "C53".
140 */
141 private final Map<String, String> pdbToNewResMap = new HashMap<>();
142 /** List of modified residues * */
143 private final Map<String, String> modRes = new HashMap<>();
144 /** Keep track of ATOM record serial numbers to match them with ANISOU records. */
145 private final HashMap<Integer, Atom> atoms = new HashMap<>();
146
147 private final Map<MolecularAssembly, BufferedReader> readers = new HashMap<>();
148 /** The current altLoc - i.e., the one we are defining a chemical system for. */
149 private Character currentAltLoc = 'A';
150 /** Character for the current chain ID. */
151 private Character currentChainID = null;
152 /** String for the current SegID. */
153 private String currentSegID = null;
154 /** Flag to indicate a mutation is requested. */
155 private boolean mutate = false;
156 private List<Mutation> mutations = null;
157 private List<Integer> resNumberList = null;
158 /** Flag to indicate if missing fields should be printed (i.e. missing B-factors). */
159 private boolean printMissingFields = true;
160 /** Number of symmetry operators when expanding to a P1 unit cell (-1 saves as current spacegroup). */
161 private int nSymOp = -1;
162 /** Number of replicates in A lattice direction (-1 defaults to unit cell). */
163 private int lValue = -1;
164 /** Number of replicates in B lattice direction (-1 defaults to unit cell). */
165 private int mValue = -1;
166 /** Number of replicates in C lattice direction (-1 defaults to unit cell). */
167 private int nValue = -1;
168 /**
169 * The serial field continues from the previous asymmetric unit when expanding to P1. This is not
170 * used when saving as the current spacegroup.
171 */
172 private int serialP1 = 0;
173 /** Assume current standard. */
174 private PDBFileStandard fileStandard = VERSION3_3;
175 /** If false, skip logging "Saving file". */
176 private boolean logWrites = true;
177 /** Keep track of the current MODEL in the file. */
178 private int modelsRead = 1;
179 /** Tracks output MODEL numbers. Unused if below zero. */
180 private int modelsWritten = -1;
181 /** Replicates vector dimensions if saving as expanded. */
182 private int[] lmn = new int[]{1,1,1};
183 private String versionFileName;
184
185 private final File readFile;
186 private List<String> remarkLines = Collections.emptyList();
187 private double lastReadLambda = Double.NaN;
188
189 /**
190 * If true, read in titratable residues in their fully protonated form.
191 */
192 private boolean constantPH = false;
193 /**
194 * If true, read in titratable residues in their protonated form.
195 */
196 private boolean rotamerTitration = false;
197 /**
198 * List of residue to rename for constantPH simulations.
199 */
200 private static final HashMap<AminoAcid3, AminoAcid3> constantPHResidueMap = new HashMap<>();
201
202 static {
203 // Lysine
204 constantPHResidueMap.put(AminoAcidUtils.AminoAcid3.LYD, AminoAcidUtils.AminoAcid3.LYS);
205 // Cysteine
206 constantPHResidueMap.put(AminoAcidUtils.AminoAcid3.CYD, AminoAcidUtils.AminoAcid3.CYS);
207 // Histidine
208 constantPHResidueMap.put(AminoAcidUtils.AminoAcid3.HID, AminoAcidUtils.AminoAcid3.HIS);
209 constantPHResidueMap.put(AminoAcidUtils.AminoAcid3.HIE, AminoAcidUtils.AminoAcid3.HIS);
210 // Aspartate
211 constantPHResidueMap.put(AminoAcidUtils.AminoAcid3.ASP, AminoAcidUtils.AminoAcid3.ASD);
212 constantPHResidueMap.put(AminoAcidUtils.AminoAcid3.ASH, AminoAcidUtils.AminoAcid3.ASD);
213 // Glutamate
214 constantPHResidueMap.put(AminoAcidUtils.AminoAcid3.GLU, AminoAcidUtils.AminoAcid3.GLD);
215 constantPHResidueMap.put(AminoAcidUtils.AminoAcid3.GLH, AminoAcidUtils.AminoAcid3.GLD);
216 }
217
218 /**
219 * List of residue to rename for rotamer titration simulations.
220 */
221 private static final HashMap<AminoAcid3, AminoAcid3> rotamerResidueMap = new HashMap<>();
222
223 static {
224 // Lysine
225 rotamerResidueMap.put(AminoAcidUtils.AminoAcid3.LYD, AminoAcidUtils.AminoAcid3.LYS);
226 // Histidine
227 rotamerResidueMap.put(AminoAcidUtils.AminoAcid3.HID, AminoAcidUtils.AminoAcid3.HIS);
228 rotamerResidueMap.put(AminoAcidUtils.AminoAcid3.HIE, AminoAcidUtils.AminoAcid3.HIS);
229 // Aspartate
230 rotamerResidueMap.put(AminoAcidUtils.AminoAcid3.ASP, AminoAcidUtils.AminoAcid3.ASH);
231 // Glutamate
232 rotamerResidueMap.put(AminoAcidUtils.AminoAcid3.GLU, AminoAcidUtils.AminoAcid3.GLH);
233 //Cysteine
234 rotamerResidueMap.put(AminoAcidUtils.AminoAcid3.CYD, AminoAcidUtils.AminoAcid3.CYS);
235 }
236
237 /**
238 * Constructor for PDBFilter.
239 *
240 * @param files a {@link java.util.List} object.
241 * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
242 * @param forceField a {@link ffx.potential.parameters.ForceField} object.
243 * @param properties a {@link org.apache.commons.configuration2.CompositeConfiguration}
244 * object.
245 */
246 public PDBFilter(List<File> files, MolecularAssembly molecularAssembly, ForceField forceField,
247 CompositeConfiguration properties) {
248 super(files, molecularAssembly, forceField, properties);
249 bondList = new ArrayList<>();
250 this.fileType = FileType.PDB;
251 readFile = files.get(0);
252 }
253
254 /**
255 * Parse the PDB File from a URL.
256 *
257 * @param file a {@link java.io.File} object.
258 * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
259 * @param forceField a {@link ffx.potential.parameters.ForceField} object.
260 * @param properties a {@link org.apache.commons.configuration2.CompositeConfiguration}
261 * object.
262 */
263 public PDBFilter(File file, MolecularAssembly molecularAssembly, ForceField forceField,
264 CompositeConfiguration properties) {
265 super(file, molecularAssembly, forceField, properties);
266 bondList = new ArrayList<>();
267 this.fileType = FileType.PDB;
268 readFile = file;
269 }
270
271 /**
272 * Parse the PDB File from a URL.
273 *
274 * @param file a {@link java.io.File} object.
275 * @param molecularAssemblies a {@link java.util.List} object.
276 * @param forceField a {@link ffx.potential.parameters.ForceField} object.
277 * @param properties a {@link org.apache.commons.configuration2.CompositeConfiguration}
278 * object.
279 */
280 public PDBFilter(File file, List<MolecularAssembly> molecularAssemblies, ForceField forceField,
281 CompositeConfiguration properties) {
282 super(file, molecularAssemblies, forceField, properties);
283 bondList = new ArrayList<>();
284 this.fileType = FileType.PDB;
285 readFile = file;
286 }
287
288 /**
289 * Constructor for PDBFilter with residue numbers.
290 *
291 * @param file a {@link java.util.List} object.
292 * @param molecularAssembly a {@link ffx.potential.MolecularAssembly} object.
293 * @param forceField a {@link ffx.potential.parameters.ForceField} object.
294 * @param properties a {@link org.apache.commons.configuration2.CompositeConfiguration}
295 * object.
296 * @param resNumberList a List of integer residue numbers for constant pH rotamer
297 * optimization.
298 */
299 public PDBFilter(File file, MolecularAssembly molecularAssembly, ForceField forceField,
300 CompositeConfiguration properties, List<Integer> resNumberList) {
301 super(file, molecularAssembly, forceField, properties);
302 bondList = new ArrayList<>();
303 this.fileType = FileType.PDB;
304 this.readFile = file;
305 this.resNumberList = resNumberList;
306 //this.chainList = chainList;
307 }
308
309
310 /**
311 * Simple method useful for converting files to PDB format.
312 *
313 * @param atom a {@link ffx.potential.bonded.Atom} object.
314 * @return Returns a PDB ATOM or HETATM record for the passed Atom.
315 */
316 public static String toPDBAtomLine(Atom atom) {
317 StringBuilder sb;
318 if (atom.isHetero()) {
319 sb = new StringBuilder("HETATM");
320 } else {
321 sb = new StringBuilder("ATOM ");
322 }
323 sb.append(repeat(" ", 74));
324
325 String name = atom.getName();
326 int nameLength = name.length();
327 if (nameLength > 4) {
328 name = name.substring(0, 4);
329 } else if (nameLength == 1) {
330 name = name + " ";
331 } else if (nameLength == 2) {
332 name = name + " ";
333 }
334 int serial = atom.getXyzIndex();
335 sb.replace(6, 16, format("%5s " + padLeft(name.toUpperCase(), 4), Hybrid36.encode(5, serial)));
336
337 Character altLoc = atom.getAltLoc();
338 if (altLoc != null) {
339 sb.setCharAt(16, altLoc);
340 } else {
341 char blankChar = ' ';
342 sb.setCharAt(16, blankChar);
343 }
344
345 String resName = atom.getResidueName();
346 sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
347
348 char chain = atom.getChainID();
349 sb.setCharAt(21, chain);
350
351 int resID = atom.getResidueNumber();
352 sb.replace(22, 26, format("%4s", Hybrid36.encode(4, resID)));
353
354 double[] xyz = atom.getXYZ(null);
355 StringBuilder decimals = new StringBuilder();
356 for (int i = 0; i < 3; i++) {
357 try {
358 decimals.append(StringUtils.fwFpDec(xyz[i], 8, 3));
359 } catch (IllegalArgumentException ex) {
360 String newValue = StringUtils.fwFpTrunc(xyz[i], 8, 3);
361 logger.info(
362 format(" XYZ coordinate %8.3f for atom %s overflowed PDB format and is truncated to %s.",
363 xyz[i], atom, newValue));
364 decimals.append(newValue);
365 }
366 }
367 try {
368 decimals.append(StringUtils.fwFpDec(atom.getOccupancy(), 6, 2));
369 } catch (IllegalArgumentException ex) {
370 logger.severe(
371 format(" Occupancy %6.2f for atom %s must be between 0 and 1.", atom.getOccupancy(),
372 atom));
373 }
374 try {
375 decimals.append(StringUtils.fwFpDec(atom.getTempFactor(), 6, 2));
376 } catch (IllegalArgumentException ex) {
377 String newValue = StringUtils.fwFpTrunc(atom.getTempFactor(), 6, 2);
378 logger.info(
379 format(" B-factor %6.2f for atom %s overflowed the PDB format and is truncated to %s.",
380 atom.getTempFactor(), atom, newValue));
381 decimals.append(newValue);
382 }
383
384 sb.replace(30, 66, decimals.toString());
385 sb.replace(78, 80, format("%2d", 0));
386 sb.append("\n");
387 return sb.toString();
388 }
389
390 public void setConstantPH(boolean constantPH) {
391 this.constantPH = constantPH;
392 }
393
394 public void setRotamerTitration(boolean rotamerTitration) {
395 this.rotamerTitration = rotamerTitration;
396 }
397
398 /** clearSegIDs */
399 public void clearSegIDs() {
400 segIDs.clear();
401 }
402
403 /** {@inheritDoc} */
404 @Override
405 public void closeReader() {
406 for (MolecularAssembly system : systems) {
407 BufferedReader br = readers.get(system);
408 if (br != null) {
409 try {
410 br.close();
411 } catch (IOException ex) {
412 logger.warning(format(" Exception in closing system %s: %s", system.toString(), ex));
413 }
414 }
415 }
416 }
417
418 @Override
419 public int countNumModels() {
420 Set<File> files = systems.stream().map(MolecularAssembly::getFile).map(File::toString).distinct()
421 .map(File::new).collect(Collectors.toSet());
422
423 // Dangers of parallelism are minimized by: unique files/filenames, read-only access.
424 return files.parallelStream().mapToInt((File fi) -> {
425 int nModelsLocal = 0;
426 try (BufferedReader br = new BufferedReader(new FileReader(fi))) {
427 String line = br.readLine();
428 while (line != null) {
429 if (line.startsWith("MODEL")) {
430 ++nModelsLocal;
431 }
432 line = br.readLine();
433 }
434 nModelsLocal = Math.max(1, nModelsLocal);
435 } catch (IOException ex) {
436 logger.info(format(" Exception in parsing file %s: %s", fi, ex));
437 }
438 return nModelsLocal;
439 }).sum();
440 }
441
442 /**
443 * Get the list of alternate locations encountered.
444 *
445 * @return the alternate location list.
446 */
447 public List<Character> getAltLocs() {
448 return altLocs;
449 }
450
451 /** {@inheritDoc} */
452 @Override
453 public OptionalDouble getLastReadLambda() {
454 return Double.isNaN(lastReadLambda) ? OptionalDouble.empty() : OptionalDouble.of(lastReadLambda);
455 }
456
457 /**
458 * Returns all the remark lines found by the last readFile call.
459 *
460 * @return Remark lines from the last readFile call.
461 */
462 @Override
463 public String[] getRemarkLines() {
464 int nRemarks = remarkLines.size();
465 return remarkLines.toArray(new String[nRemarks]);
466 }
467
468 @Override
469 public int getSnapshot() {
470 return modelsRead;
471 }
472
473 /**
474 * Mutate residue(s) as the PDB file is being parsed.
475 *
476 * @param mutations a {@link java.util.List} object.
477 */
478 public void mutate(List<Mutation> mutations) {
479 mutate = true;
480 if (this.mutations == null) {
481 this.mutations = new ArrayList<>();
482 }
483 this.mutations.addAll(mutations);
484 }
485
486 /** Parse the PDB File */
487 @Override
488 public boolean readFile() {
489 remarkLines = new ArrayList<>();
490 // First atom is #1, to match xyz file format
491 int xyzIndex = 1;
492 setFileRead(false);
493 systems.add(activeMolecularAssembly);
494
495 List<String> conects = new ArrayList<>();
496 List<String> links = new ArrayList<>();
497 List<String> ssbonds = new ArrayList<>();
498 List<String> structs = new ArrayList<>();
499 try {
500 for (File file : files) {
501 currentFile = file;
502 if (mutate) {
503 List<Character> chainIDs = new ArrayList<>();
504 try (BufferedReader br = new BufferedReader(new FileReader(file))) {
505 String line = br.readLine();
506 while (line != null) {
507 // Replace all tabs w/ 4x spaces
508 line = line.replaceAll("\t", " ");
509 String identity = line;
510 if (line.length() > 6) {
511 identity = line.substring(0, 6);
512 }
513 identity = identity.trim().toUpperCase();
514 Record record;
515 try {
516 record = Record.valueOf(identity);
517 } catch (Exception e) {
518 // Continue until the record is recognized.
519 line = br.readLine();
520 continue;
521 }
522 switch (record) {
523 case ANISOU, HETATM, ATOM -> {
524 char c22 = line.charAt(21);
525 boolean idFound = false;
526 for (Character chainID : chainIDs) {
527 if (c22 == chainID) {
528 idFound = true;
529 break;
530 }
531 }
532 if (!idFound) {
533 chainIDs.add(c22);
534 }
535 }
536 }
537 line = br.readLine();
538 }
539 for (Mutation mtn : mutations) {
540 if (!chainIDs.contains(mtn.chainChar)) {
541 if (chainIDs.size() == 1) {
542 logger.warning(
543 format(" Chain ID %c for mutation not found: only one chain %c found.",
544 mtn.chainChar, chainIDs.get(0)));
545 } else {
546 logger.warning(
547 format(" Chain ID %c for mutation not found: mutation will not proceed.",
548 mtn.chainChar));
549 }
550 }
551 }
552 } catch (IOException ioException) {
553 logger.fine(format(" Exception %s in parsing file to find chain IDs", ioException));
554 }
555 }
556
557 // Check that the current file exists and that we can read it.
558 if (currentFile == null || !currentFile.exists() || !currentFile.canRead()) {
559 return false;
560 }
561
562 // Open the current file for parsing.
563 try (BufferedReader br = new BufferedReader(new FileReader(currentFile))) {
564 // Echo the alternate location being parsed.
565 if (currentAltLoc == 'A') {
566 logger.info(format(" Reading %s", currentFile.getName()));
567 } else {
568 logger.info(format(" Reading %s alternate location %s", currentFile.getName(), currentAltLoc));
569
570 }
571 activeMolecularAssembly.setAlternateLocation(currentAltLoc);
572
573 // Reset the current chain and segID.
574 currentChainID = null;
575 currentSegID = null;
576 boolean containsInsCode = false;
577
578 // Read the first line of the file.
579 String line = br.readLine();
580
581 // Parse until END or ENDMDL is found, or to the end of the file.
582 while (line != null) {
583 // Replace all tabs w/ 4x spaces
584 line = line.replaceAll("\t", " ");
585 String identity = line;
586 if (line.length() > 6) {
587 identity = line.substring(0, 6);
588 }
589 identity = identity.trim().toUpperCase();
590 Record record;
591 try {
592 record = Record.valueOf(identity);
593 } catch (Exception e) {
594 // Continue until the record is recognized.
595 line = br.readLine();
596 continue;
597 }
598
599 // Switch on the known record.
600 switch (record) {
601 case ENDMDL:
602 case END:
603 // Setting "line" to null will exit the loop.
604 line = null;
605 continue;
606 case DBREF:
607 // =============================================================================
608 // 1 - 6 Record name "DBREF "
609 // 8 - 11 IDcode idCode ID code of this entry.
610 // 13 Character chainID Chain identifier.
611 // 15 - 18 Integer seqBegin Initial sequence number of the
612 // PDB sequence segment.
613 // 19 AChar insertBegin Initial insertion code of the
614 // PDB sequence segment.
615 // 21 - 24 Integer seqEnd Ending sequence number of the
616 // PDB sequence segment.
617 // 25 AChar insertEnd Ending insertion code of the
618 // PDB sequence segment.
619 // 27 - 32 LString database Sequence database name.
620 // 34 - 41 LString dbAccession Sequence database accession code.
621 // 43 - 54 LString dbIdCode Sequence database identification code.
622 // 56 - 60 Integer dbseqBegin Initial sequence number of the
623 // database seqment.
624 // 61 AChar idbnsBeg Insertion code of initial residue of the
625 // segment, if PDB is the reference.
626 // 63 - 67 Integer dbseqEnd Ending sequence number of the
627 // database segment.
628 // 68 AChar dbinsEnd Insertion code of the ending residue of
629 // the segment, if PDB is the reference.
630 // =============================================================================
631 Character chainID = line.substring(12, 13).toUpperCase().charAt(0);
632 int seqBegin = parseInt(line.substring(14, 18).trim());
633 int seqEnd = parseInt(line.substring(20, 24).trim());
634 int[] seqRange = dbRef.computeIfAbsent(chainID, k -> new int[2]);
635 seqRange[0] = seqBegin;
636 seqRange[1] = seqEnd;
637 break;
638 case SEQRES:
639 // =============================================================================
640 // 1 - 6 Record name "SEQRES"
641 // 8 - 10 Integer serNum Serial number of the SEQRES record for the
642 // current chain. Starts at 1 and increments
643 // by one each line. Reset to 1 for each chain.
644 // 12 Character chainID Chain identifier. This may be any single
645 // legal character, including a blank which is
646 // is used if there is only one chain.
647 // 14 - 17 Integer numRes Number of residues in the chain.
648 // This value is repeated on every record.
649 // 20 - 22 Residue name resName Residue name.
650 // 24 - 26 Residue name resName Residue name.
651 // 28 - 30 Residue name resName Residue name.
652 // 32 - 34 Residue name resName Residue name.
653 // 36 - 38 Residue name resName Residue name.
654 // 40 - 42 Residue name resName Residue name.
655 // 44 - 46 Residue name resName Residue name.
656 // 48 - 50 Residue name resName Residue name.
657 // 52 - 54 Residue name resName Residue name.
658 // 56 - 58 Residue name resName Residue name.
659 // 60 - 62 Residue name resName Residue name.
660 // 64 - 66 Residue name resName Residue name.
661 // 68 - 70 Residue name resName Residue name.
662 // =============================================================================
663 activeMolecularAssembly.addHeaderLine(line);
664 chainID = line.substring(11, 12).toUpperCase().charAt(0);
665 int serNum = parseInt(line.substring(7, 10).trim());
666 String[] chain = seqRes.get(chainID);
667 int numRes = parseInt(line.substring(13, 17).trim());
668 if (chain == null) {
669 chain = new String[numRes];
670 seqRes.put(chainID, chain);
671 }
672 int resID = (serNum - 1) * 13;
673 int end = line.length();
674 for (int start = 19; start + 3 <= end; start += 4) {
675 String res = line.substring(start, start + 3).trim();
676 if (res.isEmpty()) {
677 break;
678 }
679 chain[resID++] = res;
680 }
681 break;
682 case MODRES:
683 String modResName = line.substring(12, 15).trim();
684 String stdName = line.substring(24, 27).trim();
685 modRes.put(modResName.toUpperCase(), stdName.toUpperCase());
686 activeMolecularAssembly.addHeaderLine(line);
687 // =============================================================================
688 // 1 - 6 Record name "MODRES"
689 // 8 - 11 IDcode idCode ID code of this entry.
690 // 13 - 15 Residue name resName Residue name used in this entry.
691 // 17 Character chainID Chain identifier.
692 // 19 - 22 Integer seqNum Sequence number.
693 // 23 AChar iCode Insertion code.
694 // 25 - 27 Residue name stdRes Standard residue name.
695 // 30 - 70 String comment Description of the residue modification.
696 // =============================================================================
697 break;
698 case ANISOU:
699 // =============================================================================
700 // 1 - 6 Record name "ANISOU"
701 // 7 - 11 Integer serial Atom serial number.
702 // 13 - 16 Atom name Atom name.
703 // 17 Character altLoc Alternate location indicator
704 // 18 - 20 Residue name resName Residue name.
705 // 22 Character chainID Chain identifier.
706 // 23 - 26 Integer resSeq Residue sequence number.
707 // 27 AChar iCode Insertion code.
708 // 29 - 35 Integer u[0][0] U(1,1)
709 // 36 - 42 Integer u[1][1] U(2,2)
710 // 43 - 49 Integer u[2][2] U(3,3)
711 // 50 - 56 Integer u[0][1] U(1,2)
712 // 57 - 63 Integer u[0][2] U(1,3)
713 // 64 - 70 Integer u[1][2] U(2,3)
714 // 77 - 78 LString(2) element Element symbol, right-justified.
715 // 79 - 80 LString(2) charge Charge on the atom.
716 // =============================================================================
717 boolean deleteAnisou = properties.getBoolean("delete-anisou", false);
718 double resetBfactors = properties.getDouble("reset-bfactors", -1.0);
719 if (deleteAnisou || resetBfactors >= 0.0) {
720 break;
721 }
722 Integer serial = Hybrid36.decode(5, line.substring(6, 11));
723 Character altLoc = line.substring(16, 17).toUpperCase().charAt(0);
724 if (!altLocs.contains(altLoc)) {
725 altLocs.add(altLoc);
726 }
727 if (!altLoc.equals(' ') && !altLoc.equals('A') && !altLoc.equals(currentAltLoc)) {
728 break;
729 }
730 double[] adp = new double[6];
731 adp[0] = parseInt(line.substring(28, 35).trim()) * 1.0e-4;
732 adp[1] = parseInt(line.substring(35, 42).trim()) * 1.0e-4;
733 adp[2] = parseInt(line.substring(42, 49).trim()) * 1.0e-4;
734 adp[3] = parseInt(line.substring(49, 56).trim()) * 1.0e-4;
735 adp[4] = parseInt(line.substring(56, 63).trim()) * 1.0e-4;
736 adp[5] = parseInt(line.substring(63, 70).trim()) * 1.0e-4;
737 if (atoms.containsKey(serial)) {
738 Atom a = atoms.get(serial);
739 a.setAltLoc(altLoc);
740 a.setAnisou(adp);
741 } else {
742 logger.info(
743 format(" No ATOM record for ANISOU serial number %d has been found.", serial));
744 logger.info(format(" This ANISOU record will be ignored:\n %s", line));
745 }
746 break;
747 case ATOM:
748 // =============================================================================
749 // 1 - 6 Record name "ATOM "
750 // 7 - 11 Integer serial Atom serial number.
751 // 13 - 16 Atom name Atom name.
752 // 17 Character altLoc Alternate location indicator.
753 // 18 - 20 Residue name resName Residue name.
754 // 22 Character chainID Chain identifier.
755 // 23 - 26 Integer resSeq Residue sequence number.
756 // 27 AChar iCode Code for insertion of residues.
757 // 31 - 38 Real(8.3) x Orthogonal coordinates for X in Angstroms.
758 // 39 - 46 Real(8.3) y Orthogonal coordinates for Y in Angstroms.
759 // 47 - 54 Real(8.3) z Orthogonal coordinates for Z in Angstroms.
760 // 55 - 60 Real(6.2) occupancy Occupancy.
761 // 61 - 66 Real(6.2) tempFactor Temperature factor.
762 // 77 - 78 LString(2) element Element symbol, right-justified.
763 // 79 - 80 LString(2) charge Charge on the atom.
764 // =============================================================================
765 String name;
766 String resName;
767 String segID;
768 int resSeq;
769 boolean printAtom;
770 double[] d;
771 double occupancy;
772 double tempFactor;
773 Atom newAtom;
774 Atom returnedAtom;
775 // If it's a misnamed water, it will fall through to HETATM.
776 if (!line.substring(17, 20).trim().equals("HOH")) {
777 serial = Hybrid36.decode(5, line.substring(6, 11));
778 name = line.substring(12, 16).trim();
779 if (name.toUpperCase().contains("1H") || name.toUpperCase().contains("2H")
780 || name.toUpperCase().contains("3H")) {
781 // VERSION3_2 is presently just a placeholder for "anything non-standard".
782 fileStandard = VERSION3_2;
783 }
784 altLoc = line.substring(16, 17).toUpperCase().charAt(0);
785 if (!altLocs.contains(altLoc)) {
786 altLocs.add(altLoc);
787 }
788 if (!altLoc.equals(' ') && !altLoc.equals(currentAltLoc)) {
789 break;
790 }
791 // if (!altLoc.equals(' ') && !altLoc.equals('A') && !altLoc.equals(currentAltLoc)) {
792 // break;
793 // }
794 resName = line.substring(17, 20).trim();
795 chainID = line.substring(21, 22).charAt(0);
796 segID = getSegID(chainID);
797 resSeq = Hybrid36.decode(4, line.substring(22, 26));
798
799 char insertionCode = line.charAt(26);
800 if (insertionCode != ' ' && !containsInsCode) {
801 containsInsCode = true;
802 logger.warning(
803 " FFX support for files with " + "insertion codes is experimental. "
804 + "Residues will be renumbered to " + "eliminate insertion codes (52A "
805 + "becomes 53, 53 becomes 54, etc)");
806 }
807
808 int offset = insertionCodeCount.getOrDefault(chainID, 0);
809 String pdbResNum = format("%c%d%c", chainID, resSeq, insertionCode);
810 if (!pdbToNewResMap.containsKey(pdbResNum)) {
811 if (insertionCode != ' ') {
812 ++offset;
813 insertionCodeCount.put(chainID, offset);
814 }
815 resSeq += offset;
816 if (offset != 0) {
817 logger.info(
818 format(" Chain %c " + "residue %s-%s renumbered to %c %s-%d", chainID,
819 pdbResNum.substring(1).trim(), resName, chainID, resName, resSeq));
820 }
821 String newNum = format("%c%d", chainID, resSeq);
822 pdbToNewResMap.put(pdbResNum, newNum);
823 } else {
824 resSeq += offset;
825 }
826
827 printAtom = false;
828 if (mutate) {
829 boolean doBreak = false;
830 for (Mutation mtn : mutations) {
831 if (chainID == mtn.chainChar && resSeq == mtn.resID) {
832 String atomName = name.toUpperCase();
833
834 int isAA = AminoAcidUtils.getAminoAcidNumber(resName);
835 int isNA = NucleicAcidUtils.getNucleicAcidNumber(resName);
836
837 if ((isNA != -1 && naBackboneNames.contains(atomName)) || (isAA != -1 && backboneNames.contains(atomName))) {
838 printAtom = true;
839 resName = mtn.resName;
840 } else {
841 // pyrimidines: need N1 & C2 | purines: need N9 & C4
842 if (resName.equals("DA") || resName.equals("DG") || resName.equals("DAD") || resName.equals("DGU")) {
843 boolean isMtnPyrimidine = mtn.resName.equals("DCY") || mtn.resName.equals("DTY");
844 // log the deletion to get alchemical atoms from WT (don't include H primes)
845 if (!atomName.contains("'") || !atomName.startsWith("H")) {
846 logger.info(format(" DELETING atom %d %s of %s %d in chain %s", serial, atomName, resName, resSeq, chainID));
847 }
848 if (atomName.equals("N9")) {
849 printAtom = true;
850 resName = mtn.resName;
851 if (isMtnPyrimidine) {
852 name = "N1"; // change N9 to N1
853 }
854 } else if (atomName.equals("C4")) {
855 printAtom = true;
856 resName = mtn.resName;
857 if (isMtnPyrimidine) {
858 name = "C2"; // change C4 to C2
859 }
860 } else {
861 doBreak = true;
862 break;
863 }
864 } else if (resName.equals("DC") || resName.equals("DT") || resName.equals("DCY") || resName.equals("DTY")) {
865 boolean isMtnPurine = mtn.resName.equals("DAD") || mtn.resName.equals("DGU");
866 if (!atomName.contains("'") || !atomName.startsWith("H")) {
867 logger.info(format(" DELETING atom %d %s of %s %d in chain %s", serial, atomName, resName, resSeq, chainID));
868 }
869 if (atomName.equals("N1")) {
870 printAtom = true;
871 resName = mtn.resName;
872 if (isMtnPurine) {
873 name = "N9"; // change N1 to N9
874 }
875 } else if (atomName.equals("C2")) {
876 printAtom = true;
877 resName = mtn.resName;
878 if (isMtnPurine) {
879 name = "C4"; // change C2 to C4
880 }
881 } else {
882 doBreak = true;
883 break;
884 }
885 } else {
886 logger.info(format(" Deleting atom %s of %s %d", atomName, resName, resSeq));
887 // don't have alchemical atom logging because this would be for AA
888 doBreak = true;
889 break;
890 }
891 }
892 }
893 }
894 if (doBreak) {
895 break;
896 }
897 }
898
899 if (constantPH) {
900 AminoAcid3 aa3 = AminoAcidUtils.getAminoAcid(resName.toUpperCase());
901 if (constantPHResidueMap.containsKey(aa3)) {
902 String atomName = name.toUpperCase();
903 AminoAcid3 aa3PH = constantPHResidueMap.get(aa3);
904 resName = aa3PH.name();
905 if (constantPhBackboneNames.contains(atomName)) {
906 logger.info(format(" %s-%d %s", resName, resSeq, atomName));
907 } else if (!atomName.startsWith("H")) {
908 logger.info(format(" %s-%d %s", resName, resSeq, atomName));
909 } else {
910 logger.info(format(" %s-%d %s skipped", resName, resSeq, atomName));
911 break;
912 }
913 }
914 } else if (rotamerTitration) {
915 AminoAcid3 aa3 = AminoAcidUtils.getAminoAcid(resName.toUpperCase());
916 if (rotamerResidueMap.containsKey(aa3) && resNumberList.contains(resSeq)) {
917 AminoAcid3 aa3rotamer = rotamerResidueMap.get(aa3);
918 resName = aa3rotamer.name();
919 }
920 }
921 d = new double[3];
922 d[0] = parseDouble(line.substring(30, 38).trim());
923 d[1] = parseDouble(line.substring(38, 46).trim());
924 d[2] = parseDouble(line.substring(46, 54).trim());
925 occupancy = 1.0;
926 tempFactor = 1.0;
927 try {
928 occupancy = parseDouble(line.substring(54, 60).trim());
929 tempFactor = parseDouble(line.substring(60, 66).trim());
930 } catch (NumberFormatException | StringIndexOutOfBoundsException e) {
931 // Use default values.
932 if (printMissingFields) {
933 logger.info(" Missing occupancy and b-factors set to 1.0.");
934 printMissingFields = false;
935 } else if (logger.isLoggable(Level.FINE)) {
936 logger.fine(" Missing occupancy and b-factors set to 1.0.");
937 }
938 }
939
940 double bfactor = properties.getDouble("reset-bfactors", -1.0);
941 if (bfactor >= 0.0) {
942 tempFactor = bfactor;
943 }
944
945 newAtom = new Atom(0, name, altLoc, d, resName, resSeq, chainID, occupancy, tempFactor, segID);
946
947 // Check if this is a modified residue.
948 if (modRes.containsKey(resName.toUpperCase())) {
949 newAtom.setModRes(true);
950 }
951 returnedAtom = (Atom) activeMolecularAssembly.addMSNode(newAtom);
952 if (returnedAtom != newAtom) {
953 // A previously added atom has been retained.
954 atoms.put(serial, returnedAtom);
955 if (logger.isLoggable(Level.FINE)) {
956 logger.fine(returnedAtom + " has been retained over\n" + newAtom);
957 }
958 } else {
959 // The new atom has been added.
960 atoms.put(serial, newAtom);
961 // Check if the newAtom took the xyzIndex of a previous alternate conformer.
962 if (newAtom.getIndex() == 0) {
963 newAtom.setXyzIndex(xyzIndex++);
964 }
965 if (printAtom) {
966 logger.info(newAtom.toString());
967 }
968 }
969 break;
970 }
971 break;
972 case HETATM:
973 // =============================================================================
974 // 1 - 6 Record name "HETATM"
975 // 7 - 11 Integer serial Atom serial number.
976 // 13 - 16 Atom name Atom name.
977 // 17 Character altLoc Alternate location indicator.
978 // 18 - 20 Residue name resName Residue name.
979 // 22 Character chainID Chain identifier.
980 // 23 - 26 Integer resSeq Residue sequence number.
981 // 27 AChar iCode Code for insertion of residues.
982 // 31 - 38 Real(8.3) x Orthogonal coordinates for X.
983 // 39 - 46 Real(8.3) y Orthogonal coordinates for Y.
984 // 47 - 54 Real(8.3) z Orthogonal coordinates for Z.
985 // 55 - 60 Real(6.2) occupancy Occupancy.
986 // 61 - 66 Real(6.2) tempFactor Temperature factor.
987 // 77 - 78 LString(2) element Element symbol; right-justified.
988 // 79 - 80 LString(2) charge Charge on the atom.
989 // =============================================================================
990 serial = Hybrid36.decode(5, line.substring(6, 11));
991 name = line.substring(12, 16).trim();
992 altLoc = line.substring(16, 17).toUpperCase().charAt(0);
993 if (!altLocs.contains(altLoc)) {
994 altLocs.add(altLoc);
995 }
996 if (!altLoc.equals(' ') && !altLoc.equals(currentAltLoc)) {
997 break;
998 }
999 // if (!altLoc.equals(' ') && !altLoc.equals('A') && !altLoc.equals(currentAltLoc)) {
1000 // break;
1001 //}
1002 resName = line.substring(17, 20).trim();
1003 chainID = line.substring(21, 22).charAt(0);
1004 segID = getSegID(chainID);
1005 resSeq = Hybrid36.decode(4, line.substring(22, 26));
1006
1007 char insertionCode = line.charAt(26);
1008 if (insertionCode != ' ' && !containsInsCode) {
1009 containsInsCode = true;
1010 logger.warning(" FFX support for files with " + "insertion codes is experimental. "
1011 + "Residues will be renumbered to " + "eliminate insertion codes (52A "
1012 + "becomes 53, 53 becomes 54, etc)");
1013 }
1014
1015 int offset = insertionCodeCount.getOrDefault(chainID, 0);
1016 String pdbResNum = format("%c%d%c", chainID, resSeq, insertionCode);
1017 if (!pdbToNewResMap.containsKey(pdbResNum)) {
1018 if (insertionCode != ' ') {
1019 ++offset;
1020 insertionCodeCount.put(chainID, offset);
1021 }
1022 resSeq += offset;
1023 if (offset != 0) {
1024 logger.info(
1025 format(" Chain %c " + "molecule %s-%s renumbered to %c %s-%d", chainID,
1026 pdbResNum.substring(1).trim(), resName, chainID, resName, resSeq));
1027 }
1028 String newNum = format("%c%d", chainID, resSeq);
1029 pdbToNewResMap.put(pdbResNum, newNum);
1030 } else {
1031 resSeq += offset;
1032 }
1033
1034 d = new double[3];
1035 d[0] = parseDouble(line.substring(30, 38).trim());
1036 d[1] = parseDouble(line.substring(38, 46).trim());
1037 d[2] = parseDouble(line.substring(46, 54).trim());
1038 occupancy = 1.0;
1039 tempFactor = 1.0;
1040 try {
1041 occupancy = parseDouble(line.substring(54, 60).trim());
1042 tempFactor = parseDouble(line.substring(60, 66).trim());
1043 } catch (NumberFormatException | StringIndexOutOfBoundsException e) {
1044 // Use default values.
1045 if (printMissingFields) {
1046 logger.info(" Missing occupancy and b-factors set to 1.0.");
1047 printMissingFields = false;
1048 } else if (logger.isLoggable(Level.FINE)) {
1049 logger.fine(" Missing occupancy and b-factors set to 1.0.");
1050 }
1051 }
1052
1053 double bfactor = properties.getDouble("reset-bfactors", -1.0);
1054 if (bfactor >= 0.0) {
1055 tempFactor = bfactor;
1056 }
1057
1058 newAtom = new Atom(0, name, altLoc, d, resName, resSeq, chainID, occupancy, tempFactor, segID);
1059 newAtom.setHetero(true);
1060 // Check if this is a modified residue.
1061 if (modRes.containsKey(resName.toUpperCase())) {
1062 newAtom.setModRes(true);
1063 }
1064 returnedAtom = (Atom) activeMolecularAssembly.addMSNode(newAtom);
1065 if (returnedAtom != newAtom) {
1066 // A previously added atom has been retained.
1067 atoms.put(serial, returnedAtom);
1068 if (logger.isLoggable(Level.FINE)) {
1069 logger.fine(returnedAtom + " has been retained over\n" + newAtom);
1070 }
1071 } else {
1072 // The new atom has been added.
1073 atoms.put(serial, newAtom);
1074 newAtom.setXyzIndex(xyzIndex++);
1075 }
1076 break;
1077 case CRYST1:
1078 // =============================================================================
1079 // The CRYST1 record presents the unit cell parameters, space group, and Z
1080 // value. If the structure was not determined by crystallographic means, CRYST1
1081 // simply provides the unitary values, with an appropriate REMARK.
1082 //
1083 // 7 - 15 Real(9.3) a a (Angstroms).
1084 // 16 - 24 Real(9.3) b b (Angstroms).
1085 // 25 - 33 Real(9.3) c c (Angstroms).
1086 // 34 - 40 Real(7.2) alpha alpha (degrees).
1087 // 41 - 47 Real(7.2) beta beta (degrees).
1088 // 48 - 54 Real(7.2) gamma gamma (degrees).
1089 // 56 - 66 LString sGroup Space group.
1090 // 67 - 70 Integer z Z value.
1091 // =============================================================================
1092 if (line.length() < 55) {
1093 logger.severe(" CRYST1 record is improperly formatted.");
1094 }
1095 double aaxis = parseDouble(line.substring(6, 15).trim());
1096 double baxis = parseDouble(line.substring(15, 24).trim());
1097 double caxis = parseDouble(line.substring(24, 33).trim());
1098 double alpha = parseDouble(line.substring(33, 40).trim());
1099 double beta = parseDouble(line.substring(40, 47).trim());
1100 double gamma = parseDouble(line.substring(47, 54).trim());
1101 int limit = min(line.length(), 66);
1102 String sg = line.substring(55, limit).trim();
1103 properties.addProperty("a-axis", aaxis);
1104 properties.addProperty("b-axis", baxis);
1105 properties.addProperty("c-axis", caxis);
1106 properties.addProperty("alpha", alpha);
1107 properties.addProperty("beta", beta);
1108 properties.addProperty("gamma", gamma);
1109 properties.addProperty("spacegroup", SpaceGroupInfo.pdb2ShortName(sg));
1110 break;
1111 case CONECT:
1112 // =============================================================================
1113 // 7 - 11 Integer serial Atom serial number
1114 // 12 - 16 Integer serial Serial number of bonded atom
1115 // 17 - 21 Integer serial Serial number of bonded atom
1116 // 22 - 26 Integer serial Serial number of bonded atom
1117 // 27 - 31 Integer serial Serial number of bonded atom
1118 //
1119 // CONECT records involving atoms for which the coordinates are not present
1120 // in the entry (e.g., symmetry-generated) are not given.
1121 // CONECT records involving atoms for which the coordinates are missing due
1122 // to disorder, are also not provided.
1123 // =============================================================================
1124 conects.add(line);
1125 break;
1126 case LINK:
1127 // =============================================================================
1128 // The LINK records specify connectivity between residues that is not implied by
1129 // the primary structure. Connectivity is expressed in terms of the atom names.
1130 // They also include the distance associated with each linkage following the
1131 // symmetry operations at the end of each record.
1132 // 13 - 16 Atom name1 Atom name.
1133 // 17 Character altLoc1 Alternate location indicator.
1134 // 18 - 20 Residue name resName1 Residue name.
1135 // 22 Character chainID1 Chain identifier.
1136 // 23 - 26 Integer resSeq1 Residue sequence number.
1137 // 27 AChar iCode1 Insertion code.
1138 // 43 - 46 Atom name2 Atom name.
1139 // 47 Character altLoc2 Alternate location indicator.
1140 // 48 - 50 Residue name resName2 Residue name.
1141 // 52 Character chainID2 Chain identifier.
1142 // 53 - 56 Integer resSeq2 Residue sequence number.
1143 // 57 AChar iCode2 Insertion code.
1144 // 60 - 65 SymOP sym1 Symmetry operator atom 1.
1145 // 67 - 72 SymOP sym2 Symmetry operator atom 2.
1146 // 74 – 78 Real(5.2) Length Link distance
1147 // =============================================================================
1148 char a1 = line.charAt(16);
1149 char a2 = line.charAt(46);
1150 if (a1 != a2) {
1151 // logger.info(format(" Ignoring LINK record as alternate locations do not match\n
1152 // %s.", line));
1153 break;
1154 }
1155 if (currentAltLoc == 'A') {
1156 if ((a1 == ' ' || a1 == 'A') && (a2 == ' ' || a2 == 'A')) {
1157 links.add(line);
1158 }
1159 } else if (a1 == currentAltLoc && a2 == currentAltLoc) {
1160 links.add(line);
1161 }
1162 break;
1163 case SSBOND:
1164 // =============================================================================
1165 // The SSBOND record identifies each disulfide bond in protein and polypeptide
1166 // structures by identifying the two residues involved in the bond.
1167 // The disulfide bond distance is included after the symmetry operations at
1168 // the end of the SSBOND record.
1169 //
1170 // 8 - 10 Integer serNum Serial number.
1171 // 12 - 14 LString(3) "CYS" Residue name.
1172 // 16 Character chainID1 Chain identifier.
1173 // 18 - 21 Integer seqNum1 Residue sequence number.
1174 // 22 AChar icode1 Insertion code.
1175 // 26 - 28 LString(3) "CYS" Residue name.
1176 // 30 Character chainID2 Chain identifier.
1177 // 32 - 35 Integer seqNum2 Residue sequence number.
1178 // 36 AChar icode2 Insertion code.
1179 // 60 - 65 SymOP sym1 Symmetry oper for 1st resid
1180 // 67 - 72 SymOP sym2 Symmetry oper for 2nd resid
1181 // 74 – 78 Real(5.2) Length Disulfide bond distance
1182 //
1183 // If SG of cysteine is disordered then there are possible alternate linkages.
1184 // wwPDB practice is to put together all possible SSBOND records. This is
1185 // problematic because the alternate location identifier is not specified in
1186 // the SSBOND record.
1187 //
1188 // Notes:
1189 // SSBOND records may be invalid if chain IDs are reused.
1190 // SSBOND records are applied by FFX to all conformers.
1191 // =============================================================================
1192 ssbonds.add(line);
1193 break;
1194 case HELIX:
1195 // =============================================================================
1196 // HELIX records are used to identify the position of helices in the molecule.
1197 // Helices are named, numbered, and classified by type. The residues where the
1198 // helix begins and ends are noted, as well as the total length.
1199 //
1200 // 8 - 10 Integer serNum Serial number of the helix. This starts
1201 // at 1 and increases incrementally.
1202 // 12 - 14 LString(3) helixID Helix identifier. In addition to a serial
1203 // number, each helix is given an
1204 // alphanumeric character helix identifier.
1205 // 16 - 18 Residue name initResName Name of the initial residue.
1206 // 20 Character initChainID Chain identifier for the chain containing
1207 // this helix.
1208 // 22 - 25 Integer initSeqNum Sequence number of the initial residue.
1209 // 26 AChar initICode Insertion code of the initial residue.
1210 // 28 - 30 Residue name endResName Name of the terminal residue of the helix.
1211 // 32 Character endChainID Chain identifier for the chain containing
1212 // this helix.
1213 // 34 - 37 Integer endSeqNum Sequence number of the terminal residue.
1214 // 38 AChar endICode Insertion code of the terminal residue.
1215 // 39 - 40 Integer helixClass Helix class (see below).
1216 // 41 - 70 String comment Comment about this helix.
1217 // 72 - 76 Integer length Length of this helix.
1218 //
1219 // CLASS NUMBER
1220 // TYPE OF HELIX (COLUMNS 39 - 40)
1221 // --------------------------------------------------------------
1222 // Right-handed alpha (default) 1
1223 // Right-handed omega 2
1224 // Right-handed pi 3
1225 // Right-handed gamma 4
1226 // Right-handed 3 - 10 5
1227 // Left-handed alpha 6
1228 // Left-handed omega 7
1229 // Left-handed gamma 8
1230 // 2 - 7 ribbon/helix 9
1231 // Polyproline 10
1232 // =============================================================================
1233 case SHEET:
1234 // =============================================================================
1235 // SHEET records are used to identify the position of sheets in the molecule.
1236 // Sheets are both named and numbered. The residues where the sheet begins and
1237 // ends are noted.
1238 //
1239 // 8 - 10 Integer strand Strand number which starts at 1 for each
1240 // strand within a sheet and increases by one.
1241 // 12 - 14 LString(3) sheetID Sheet identifier.
1242 // 15 - 16 Integer numStrands Number of strands in sheet.
1243 // 18 - 20 Residue name initResName Residue name of initial residue.
1244 // 22 Character initChainID Chain identifier of initial residue in strand.
1245 // 23 - 26 Integer initSeqNum Sequence number of initial residue in strand.
1246 // 27 AChar initICode Insertion code of initial residue in strand.
1247 // 29 - 31 Residue name endResName Residue name of terminal residue.
1248 // 33 Character endChainID Chain identifier of terminal residue.
1249 // 34 - 37 Integer endSeqNum Sequence number of terminal residue.
1250 // 38 AChar endICode Insertion code of terminal residue.
1251 // 39 - 40 Integer sense Sense of strand with respect to previous
1252 // strand in the sheet. 0 if first strand,
1253 // 1 if parallel,and -1 if anti-parallel.
1254 // 42 - 45 Atom curAtom Registration. Atom name in current strand.
1255 // 46 - 48 Residue name curResName Registration. Residue name in current strand
1256 // 50 Character curChainId Registration. Chain identifier in current strand.
1257 // 51 - 54 Integer curResSeq Registration. Residue sequence number
1258 // in current strand.
1259 // 55 AChar curICode Registration. Insertion code in current strand.
1260 // 57 - 60 Atom prevAtom Registration. Atom name in previous strand.
1261 // 61 - 63 Residue name prevResName Registration. Residue name in previous strand.
1262 // 65 Character prevChainId Registration. Chain identifier in previous strand.
1263 // 66 - 69 Integer prevResSeq Registration. Residue sequence number
1264 // in previous strand.
1265 // 70 AChar prevICode Registration. Insertion code in
1266 // previous strand.
1267 // =============================================================================
1268 structs.add(line);
1269 break;
1270 case MODEL: // Currently, no handling in initial read.
1271 break;
1272 case MTRIX1:
1273 // ================================================================================
1274 // MTRIXn (n = 1, 2, or 3) records present transformations expressing
1275 // non-crystallographic symmetry.
1276 // MTRIXn will appear only when such transformations are required to generate an
1277 // entire asymmetric unit,
1278 // such as a large viral structure.
1279 //
1280 // 8 - 10 Integer serial Serial number.
1281 // 11 - 20 Real(10.6) m[n][1] Mn1
1282 // 21 - 30 Real(10.6) m[n][2] Mn2
1283 // 31 - 40 Real(10.6) m[n][3] Mn3
1284 // 21 - 30 Real(10.6) v[n] Vn
1285 // 60 Integer iGiven 1 if coordinates for the representations which are
1286 // approximately related by the transformations of the
1287 // molecule are contained in the entry. Otherwise, blank.
1288 // =================================================================================
1289 StringBuilder MTRX1 = new StringBuilder(line.substring(11, 55));
1290 properties.addProperty("MTRIX1", MTRX1);
1291 break;
1292 case MTRIX2:
1293 StringBuilder MTRX2 = new StringBuilder(line.substring(11, 55));
1294 properties.addProperty("MTRIX2", MTRX2);
1295 break;
1296 case MTRIX3:
1297 StringBuilder MTRX3 = new StringBuilder(line.substring(11, 55));
1298 properties.addProperty("MTRIX3", MTRX3);
1299 break;
1300 case REMARK:
1301 remarkLines.add(line.trim());
1302 if (line.contains("Lambda:")) {
1303 Matcher m = lambdaPattern.matcher(line);
1304 if (m.find()) {
1305 lastReadLambda = Double.parseDouble(m.group(1));
1306 }
1307 }
1308 // =================================================================================
1309 // REMARK 350: presents all transformations, both crystallographic and non-crystallographic,
1310 // needed to generate the biomolecule. These transformations operate on the coordinates in the
1311 // entry. Both author and computational descriptions of assemblies are provided, if applicable.
1312 // For strict ncs case where more than one assembly presents in asymmetric unit, only one
1313 // chain with unit matrix will reported in REMARK 350, the other chain will be generated
1314 // by rotation and translation.
1315 //
1316 // 20 - 23 Integer serial Serial number.
1317 // 24 - 33 Real(10.6) m[n][1] Mn1
1318 // 34 - 43 Real(10.6) m[n][2] Mn2
1319 // 44 - 53 Real(10.6) m[n][3] Mn3
1320 // 59 - 68 Real(10.6) v[n] Vn
1321 // =================================================================================
1322 if (line.length() >= 68) {
1323 String remarkType = line.substring(7, 10).trim();
1324 if (remarkType.matches("\\d+") && parseInt(remarkType) == 350 && line.substring(13,
1325 18).equalsIgnoreCase("BIOMT")) {
1326 properties.addProperty("BIOMTn", new StringBuilder(line.substring(24, 68)));
1327 }
1328 }
1329 break;
1330 default:
1331 break;
1332 }
1333 line = br.readLine();
1334 }
1335
1336 } catch (FileNotFoundException fileNotFoundException) {
1337 logger.log(Level.SEVERE, " PDB file not found", fileNotFoundException);
1338 }
1339 }
1340 xyzIndex--;
1341 setFileRead(true);
1342 } catch (IOException e) {
1343 logger.exiting(PDBFilter.class.getName(), "readFile", e);
1344 return false;
1345 }
1346
1347 // Locate disulfide bonds; bond parameters are assigned below.
1348 List<Bond> ssBondList = locateDisulfideBonds(ssbonds, activeMolecularAssembly, pdbToNewResMap);
1349
1350 // Record the number of atoms read in from the PDB file before applying
1351 // algorithms that may build new atoms.
1352 int pdbAtoms = activeMolecularAssembly.getAtomArray().length;
1353
1354 // Build missing backbone atoms in loops.
1355 buildMissingResidues(xyzIndex, activeMolecularAssembly, seqRes, dbRef);
1356
1357 // Assign atom types. Missing side-chains atoms and missing hydrogen will be built in.
1358 bondList = assignAtomTypes(activeMolecularAssembly, fileStandard);
1359
1360 // Assign disulfide bonds parameters and log their creation.
1361 buildDisulfideBonds(ssBondList, activeMolecularAssembly, bondList);
1362
1363 // Finally, re-number the atoms if missing atoms were created.
1364 int currentN = activeMolecularAssembly.getAtomArray().length;
1365 boolean renumber = forceField.getBoolean("renumber-pdb", false);
1366 if (pdbAtoms != currentN) {
1367 logger.info(format(" Renumbering PDB file due to built atoms (%d vs %d)", currentN, pdbAtoms));
1368 numberAtoms(activeMolecularAssembly);
1369 } else if (renumber) {
1370 logger.info(" Renumbering PDB file due to renumber-pdb flag.");
1371 numberAtoms(activeMolecularAssembly);
1372 }
1373 return true;
1374 }
1375
1376 /** {@inheritDoc} */
1377 @Override
1378 public boolean readNext() {
1379 return readNext(false);
1380 }
1381
1382 /** {@inheritDoc} */
1383 @Override
1384 public boolean readNext(boolean resetPosition) {
1385 return readNext(resetPosition, false);
1386 }
1387
1388 /** {@inheritDoc} */
1389 @Override
1390 public boolean readNext(boolean resetPosition, boolean print) {
1391 return readNext(resetPosition, print, true);
1392 }
1393
1394 /** {@inheritDoc} */
1395 @Override
1396 public boolean readNext(boolean resetPosition, boolean print, boolean parse) {
1397 modelsRead = resetPosition ? 1 : modelsRead + 1;
1398 if (!parse) {
1399 if (print) {
1400 logger.info(format(" Skipped Model %d.", modelsRead));
1401 }
1402 return true;
1403 }
1404 remarkLines = new ArrayList<>(remarkLines.size());
1405 // ^ is beginning of line, \\s+ means "one or more whitespace", (\\d+) means match and capture
1406 // one or more digits.
1407 Pattern modelPatt = Pattern.compile("^MODEL\\s+(\\d+)");
1408 boolean eof = true;
1409 for (MolecularAssembly system : systems) {
1410 try {
1411 BufferedReader currentReader;
1412 if (readers.containsKey(system)) {
1413 currentReader = readers.get(system);
1414 try {
1415 if (!currentReader.ready()) {
1416 currentReader = new BufferedReader(new FileReader(readFile));
1417 // Mark the start of the file.
1418 currentReader.mark(0);
1419 readers.remove(system);
1420 readers.put(system, currentReader);
1421 } else if (resetPosition) {
1422 // If the BufferedReader has been opened, and reset is requested, reset the position.
1423 currentReader.reset();
1424 }
1425 } catch (Exception exception) {
1426 // If all structures in the PDB file have been read, the currentReader may have closed.
1427 // The try block will catch this case and reset to the beginning of the file.
1428 currentReader = new BufferedReader(new FileReader(readFile));
1429 // Mark the start of the file.
1430 currentReader.mark(0);
1431 readers.remove(system);
1432 readers.put(system, currentReader);
1433 }
1434 } else {
1435 currentReader = new BufferedReader(new FileReader(readFile));
1436 // Mark the start of the file.
1437 currentReader.mark(0);
1438 readers.put(system, currentReader);
1439 }
1440
1441 // Skip to appropriate model.
1442 String line = currentReader.readLine();
1443 while (line != null) {
1444 line = line.trim();
1445 Matcher m = modelPatt.matcher(line);
1446 if (m.find()) {
1447 int modelNum = parseInt(m.group(1));
1448 if (modelNum == modelsRead) {
1449 if (print) {
1450 logger.log(Level.INFO, format(" Reading model %d for %s", modelNum, currentFile));
1451 }
1452 eof = false;
1453 break;
1454 }
1455 }
1456 line = currentReader.readLine();
1457 }
1458 if (eof) {
1459 if (logger.isLoggable(Level.FINEST)) {
1460 logger.log(Level.FINEST, format("\n End of file reached for %s", readFile));
1461 }
1462 currentReader.close();
1463 return false;
1464 }
1465
1466 // Begin parsing the model.
1467 boolean modelDone = false;
1468 line = currentReader.readLine();
1469 while (line != null) {
1470 line = line.trim();
1471 String recID = line.substring(0, Math.min(6, line.length())).trim();
1472 try {
1473 Record record = Record.valueOf(recID);
1474 boolean hetatm = true;
1475 switch (record) {
1476 // =============================================================================
1477 //
1478 // 7 - 11 Integer serial Atom serial number.
1479 // 13 - 16 Atom name Atom name.
1480 // 17 Character altLoc Alternate location indicator.
1481 // 18 - 20 Residue name resName Residue name.
1482 // 22 Character chainID Chain identifier.
1483 // 23 - 26 Integer resSeq Residue sequence number.
1484 // 27 AChar iCode Code for insertion of residues.
1485 // 31 - 38 Real(8.3) x Orthogonal coordinates for X in
1486 // Angstroms.
1487 // 39 - 46 Real(8.3) y Orthogonal coordinates for Y in
1488 // Angstroms.
1489 // 47 - 54 Real(8.3) z Orthogonal coordinates for Z in
1490 // Angstroms.
1491 // 55 - 60 Real(6.2) occupancy Occupancy.
1492 // 61 - 66 Real(6.2) tempFactor Temperature factor.
1493 // 77 - 78 LString(2) element Element symbol, right-justified.
1494 // 79 - 80 LString(2) charge Charge on the atom.
1495 // =============================================================================
1496 // 1 2 3 4 5 6 7
1497 // 123456789012345678901234567890123456789012345678901234567890123456789012345678
1498 // ATOM 1 N ILE A 16 60.614 71.140 -10.592 1.00 7.38 N
1499 // ATOM 2 CA ILE A 16 60.793 72.149 -9.511 1.00 6.91 C
1500 case ATOM:
1501 hetatm = false;
1502 case HETATM:
1503 String name = line.substring(12, 16).trim();
1504 if (name.toUpperCase().contains("1H") || name.toUpperCase().contains("2H")
1505 || name.toUpperCase().contains("3H")) {
1506 // VERSION3_2 is presently just a placeholder for "anything non-standard".
1507 fileStandard = VERSION3_2;
1508 }
1509 Character altLoc = line.substring(16, 17).toUpperCase().charAt(0);
1510 if (!altLoc.equals(' ') && !altLoc.equals(currentAltLoc)) {
1511 break;
1512 }
1513 // if (!altLoc.equals(' ') && !altLoc.equals('A') && !altLoc.equals(currentAltLoc)) {
1514 // break;
1515 // }
1516 String resName = line.substring(17, 20).trim();
1517 Character chainID = line.substring(21, 22).charAt(0);
1518 String segID = getExistingSegID(chainID);
1519 int resSeq = Hybrid36.decode(4, line.substring(22, 26));
1520 double[] d = new double[3];
1521 d[0] = parseDouble(line.substring(30, 38).trim());
1522 d[1] = parseDouble(line.substring(38, 46).trim());
1523 d[2] = parseDouble(line.substring(46, 54).trim());
1524 double occupancy = 1.0;
1525 double tempFactor = 1.0;
1526 Atom newAtom = new Atom(0, name, altLoc, d, resName, resSeq, chainID, occupancy,
1527 tempFactor, segID);
1528 newAtom.setHetero(hetatm);
1529 // Check if this is a modified residue.
1530 if (modRes.containsKey(resName.toUpperCase())) {
1531 newAtom.setModRes(true);
1532 }
1533
1534 Atom returnedAtom = activeMolecularAssembly.findAtom(newAtom);
1535 if (returnedAtom != null) {
1536 returnedAtom.setXYZ(d);
1537 double[] retXYZ = new double[3];
1538 returnedAtom.getXYZ(retXYZ);
1539 } else {
1540 String message = format(" Could not find atom %s in assembly", newAtom);
1541 if (dieOnMissingAtom) {
1542 logger.severe(message);
1543 } else {
1544 logger.warning(message);
1545 }
1546 }
1547 break;
1548 case CRYST1:
1549 // =============================================================================
1550 // The CRYST1 record presents the unit cell parameters, space group, and Z
1551 // value. If the structure was not determined by crystallographic means, CRYST1
1552 // simply provides the unitary values, with an appropriate REMARK.
1553 //
1554 // 7 - 15 Real(9.3) a a (Angstroms).
1555 // 16 - 24 Real(9.3) b b (Angstroms).
1556 // 25 - 33 Real(9.3) c c (Angstroms).
1557 // 34 - 40 Real(7.2) alpha alpha (degrees).
1558 // 41 - 47 Real(7.2) beta beta (degrees).
1559 // 48 - 54 Real(7.2) gamma gamma (degrees).
1560 // 56 - 66 LString sGroup Space group.
1561 // 67 - 70 Integer z Z value.
1562 // =============================================================================
1563 logger.fine(" Crystal record found.");
1564 if (line.length() < 55) {
1565 logger.severe(" CRYST1 record is improperly formatted.");
1566 }
1567 double aaxis = parseDouble(line.substring(6, 15).trim());
1568 double baxis = parseDouble(line.substring(15, 24).trim());
1569 double caxis = parseDouble(line.substring(24, 33).trim());
1570 double alpha = parseDouble(line.substring(33, 40).trim());
1571 double beta = parseDouble(line.substring(40, 47).trim());
1572 double gamma = parseDouble(line.substring(47, 54).trim());
1573 int limit = min(line.length(), 66);
1574 String sg = line.substring(55, limit).trim();
1575 // properties.clearProperty("a-axis");
1576 // properties.clearProperty("b-axis");
1577 // properties.clearProperty("c-axis");
1578 // properties.clearProperty("alpha");
1579 // properties.clearProperty("beta");
1580 // properties.clearProperty("gamma");
1581 // properties.clearProperty("spacegroup");
1582 //
1583 // properties.addProperty("a-axis", aaxis);
1584 // properties.addProperty("b-axis", baxis);
1585 // properties.addProperty("c-axis", caxis);
1586 // properties.addProperty("alpha", alpha);
1587 // properties.addProperty("beta", beta);
1588 // properties.addProperty("gamma", gamma);
1589 // properties.addProperty("spacegroup", SpaceGroupInfo.pdb2ShortName(sg));
1590 Crystal crystal = activeMolecularAssembly.getCrystal();
1591 SpaceGroup spaceGroup = SpaceGroupDefinitions.spaceGroupFactory(sg);
1592 if (Objects.equals(crystal.spaceGroup.shortName, spaceGroup.shortName)) {
1593 crystal.changeUnitCellParameters(aaxis, baxis, caxis, alpha, beta, gamma);
1594 } else {
1595 // TODO: Handle changes in space groups... Means recalculating force field terms.
1596 logger.warning(format(" Original space group %s could not be changed to %s",
1597 crystal.spaceGroup.shortName, spaceGroup.shortName));
1598 }
1599 break;
1600 case ENDMDL:
1601 case END: // Technically speaking, END should be at the end of the file, not end of
1602 // the model.
1603 logger.log(Level.FINE, format(" Model %d successfully read", modelsRead));
1604 modelDone = true;
1605 break;
1606 case REMARK:
1607 remarkLines.add(line.trim());
1608 if (line.contains("Lambda:")) {
1609 Matcher m = lambdaPattern.matcher(line);
1610 if (m.find()) {
1611 lastReadLambda = Double.parseDouble(m.group(1));
1612 }
1613 }
1614 break;
1615 default:
1616 break;
1617 }
1618 } catch (Exception ex) {
1619 // Do nothing; it's not an ATOM/HETATM line.
1620 }
1621 if (modelDone) {
1622 break;
1623 }
1624 line = currentReader.readLine();
1625 }
1626 return true;
1627 } catch (IOException ex) {
1628 logger.info(
1629 format(" Exception in parsing frame %d of %s:" + " %s", modelsRead, system.toString(),
1630 ex));
1631 }
1632 }
1633 return false;
1634 }
1635
1636 /**
1637 * Specify the alternate location.
1638 *
1639 * @param molecularAssembly The MolecularAssembly to populate.
1640 * @param altLoc The alternate location to use.
1641 */
1642 public void setAltID(MolecularAssembly molecularAssembly, Character altLoc) {
1643 setMolecularSystem(molecularAssembly);
1644 currentAltLoc = altLoc;
1645 }
1646
1647 /**
1648 * Sets whether this PDBFilter should log each time it saves to a file.
1649 *
1650 * @param logWrites a boolean.
1651 */
1652 public void setLogWrites(boolean logWrites) {
1653 this.logWrites = logWrites;
1654 }
1655
1656 /**
1657 * setModelNumbering.
1658 *
1659 * @param modelsWritten the number of models written.
1660 */
1661 public void setModelNumbering(int modelsWritten) {
1662 this.modelsWritten = modelsWritten;
1663 }
1664
1665 public void setLMN(int[] lmn) {
1666 if(lmn[0] >= 1 && lmn[1] >= 1 && lmn[2] >= 1){
1667 this.lmn = lmn;
1668 }else{
1669 // Provided dimensions are not handled. Revert to P1.
1670 this.lmn = new int[]{1,1,1};
1671 }
1672 }
1673
1674 /**
1675 * setSymOp.
1676 *
1677 * @param symOp a int.
1678 */
1679 public void setSymOp(int symOp) {
1680 this.nSymOp = symOp;
1681 }
1682
1683 /**
1684 * Expand the current system to P1 during the save operation.
1685 *
1686 * @param file The file to write.
1687 * @return Return true on a successful write.
1688 */
1689 public boolean writeFileAsP1(File file) {
1690 // XYZ File First Line
1691 final int l = lmn[0];
1692 final int m = lmn[1];
1693 final int n = lmn[2];
1694 final int numReplicates = l * m * n;
1695 Crystal crystal = activeMolecularAssembly.getCrystal();
1696 int nSymOps = crystal.getUnitCell().spaceGroup.getNumberOfSymOps();
1697
1698 if (nSymOps == 1 && l <= 1 && m <= 1 && n <= 1) {
1699 // This is a P1 system.
1700 if (!writeFile(file, false)) {
1701 logger.info(format(" Save failed for %s", activeMolecularAssembly));
1702 return false;
1703 } else {
1704 return true;
1705 }
1706 } else {
1707 Polymer[] polymers = activeMolecularAssembly.getChains();
1708 int chainCount = 0;
1709 for (Polymer polymer : polymers) {
1710 Character chainID = Polymer.CHAIN_IDS.charAt(chainCount++);
1711 polymer.setChainID(chainID);
1712 polymer.setSegID(chainID.toString());
1713 }
1714 nSymOp = 0;
1715 logWrites = false;
1716 boolean writeEnd = false;
1717 if (!writeFile(file, false, false, writeEnd)) {
1718 logger.info(format(" Save failed for %s", activeMolecularAssembly));
1719 return false;
1720 } else {
1721 for (int i = 0; i < l; i++) {
1722 for (int j = 0; j < m; j++) {
1723 for (int k = 0; k < n; k++) {
1724 lValue = i;
1725 mValue = j;
1726 nValue = k;
1727 for (int iSym = 0; iSym < nSymOps; iSym++) {
1728 nSymOp = iSym;
1729 for (Polymer polymer : polymers) {
1730 Character chainID = Polymer.CHAIN_IDS.charAt(chainCount++);
1731 polymer.setChainID(chainID);
1732 polymer.setSegID(chainID.toString());
1733 }
1734 // If the last sym op to be written.
1735 writeEnd = iSym == nSymOps - 1 && i == l - 1 && j == m - 1 && k == n - 1;
1736 if (!writeFile(file, true, false, writeEnd)) {
1737 logger.info(format(" Save failed for %s", activeMolecularAssembly));
1738 return false;
1739 }
1740 }
1741 }
1742 }
1743 }
1744 }
1745
1746 // Reset the chainIDs.
1747 chainCount = 0;
1748 for (Polymer polymer : polymers) {
1749 Character chainID = Polymer.CHAIN_IDS.charAt(chainCount++);
1750 polymer.setChainID(chainID);
1751 polymer.setSegID(chainID.toString());
1752 }
1753
1754 }
1755
1756 return true;
1757 }
1758
1759 /**
1760 * writeFile
1761 *
1762 * @param saveFile a {@link java.io.File} object.
1763 * @param append Whether to append to saveFile (vs over-write).
1764 * @param printLinear Ignored (remains to present a different method signature).
1765 * @param writeEnd True if this is the final model.
1766 * @return Success of writing.
1767 */
1768 public boolean writeFile(File saveFile, boolean append, boolean printLinear, boolean writeEnd) {
1769 return writeFile(saveFile, append, Collections.emptySet(), writeEnd, true);
1770 }
1771
1772 /**
1773 * writeFile
1774 *
1775 * @param saveFile a {@link java.io.File} object to save to.
1776 * @param append Whether to append to saveFile (vs over-write).
1777 * @param toExclude A {@link java.util.Set} of {@link ffx.potential.bonded.Atom}s to exclude
1778 * from writing.
1779 * @param writeEnd True if this is the final model.
1780 * @param versioning True if the file being saved to should be versioned. False if the file
1781 * being saved to should be overwritten.
1782 * @return Success of writing.
1783 */
1784 public boolean writeFile(File saveFile, boolean append, Set<Atom> toExclude, boolean writeEnd,
1785 boolean versioning) {
1786 return writeFile(saveFile, append, toExclude, writeEnd, versioning, null);
1787 }
1788
1789 /**
1790 * writeFile
1791 *
1792 * @param saveFile a {@link java.io.File} object to save to.
1793 * @param append Whether to append to saveFile (vs over-write).
1794 * @param toExclude A {@link java.util.Set} of {@link ffx.potential.bonded.Atom}s to exclude
1795 * from writing.
1796 * @param writeEnd True if this is the final model.
1797 * @param versioning True if the file being saved to should be versioned. False if the file
1798 * being saved to should be overwritten.
1799 * @param extraLines Extra comment/header lines to write.
1800 * @return Success of writing.
1801 */
1802 public boolean writeFile(File saveFile, boolean append, Set<Atom> toExclude, boolean writeEnd,
1803 boolean versioning, String[] extraLines) {
1804 if (standardizeAtomNames) {
1805 logger.info(" Setting atom names to PDB standard.");
1806 renameAtomsToPDBStandard(activeMolecularAssembly);
1807 }
1808 final Set<Atom> atomExclusions = toExclude == null ? Collections.emptySet() : toExclude;
1809 if (saveFile == null) {
1810 return false;
1811 }
1812 if (vdwH) {
1813 logger.info(" Saving hydrogen to van der Waals centers instead of nuclear locations.");
1814 }
1815 if (nSymOp > -1) {
1816 logger.info(format(" Saving atoms using the symmetry operator:\n%s\n",
1817 activeMolecularAssembly.getCrystal().getUnitCell().spaceGroup.getSymOp(nSymOp)
1818 .toString()));
1819 }
1820
1821 // Create StringBuilders for ATOM, ANISOU and TER records that can be reused.
1822 StringBuilder sb = new StringBuilder("ATOM ");
1823 StringBuilder anisouSB = new StringBuilder("ANISOU");
1824 StringBuilder terSB = new StringBuilder("TER ");
1825 StringBuilder model = null;
1826 for (int i = 6; i < 80; i++) {
1827 sb.append(' ');
1828 anisouSB.append(' ');
1829 terSB.append(' ');
1830 }
1831
1832 File newFile = saveFile;
1833 if (!append) {
1834 if (versioning) {
1835 newFile = version(saveFile);
1836 }
1837 } else if (modelsWritten >= 0) {
1838 model = new StringBuilder(format("MODEL %-4d", ++modelsWritten));
1839 model.append(repeat(" ", 65));
1840 }
1841 activeMolecularAssembly.setFile(newFile);
1842 if (activeMolecularAssembly.getName() == null) {
1843 activeMolecularAssembly.setName(newFile.getName());
1844 }
1845 if (logWrites) {
1846 logger.log(Level.INFO, " Saving {0}", newFile.getName());
1847 }
1848
1849 try (FileWriter fw = new FileWriter(newFile, append);
1850 BufferedWriter bw = new BufferedWriter(fw)) {
1851 /*
1852 Will come before CRYST1 and ATOM records, but after anything
1853 written by writeFileWithHeader (particularly X-ray refinement statistics).
1854 */
1855 String[] headerLines = activeMolecularAssembly.getHeaderLines();
1856 for (String line : headerLines) {
1857 bw.write(format("%s\n", line));
1858 }
1859 if (extraLines != null) {
1860 if (rotamerTitration && extraLines[0].contains("REMARK")) {
1861 for (String line : extraLines) {
1862 bw.write(line + "\n");
1863 }
1864 } else {
1865 for (String line : extraLines) {
1866 bw.write(format("REMARK 999 %s\n", line));
1867 }
1868 }
1869 }
1870 if (model != null) {
1871 bw.write(model.toString());
1872 bw.newLine();
1873 }
1874 // =============================================================================
1875 // The CRYST1 record presents the unit cell parameters, space group, and Z
1876 // value. If the structure was not determined by crystallographic means, CRYST1
1877 // simply provides the unitary values, with an appropriate REMARK.
1878 //
1879 // 7 - 15 Real(9.3) a a (Angstroms).
1880 // 16 - 24 Real(9.3) b b (Angstroms).
1881 // 25 - 33 Real(9.3) c c (Angstroms).
1882 // 34 - 40 Real(7.2) alpha alpha (degrees).
1883 // 41 - 47 Real(7.2) beta beta (degrees).
1884 // 48 - 54 Real(7.2) gamma gamma (degrees).
1885 // 56 - 66 LString sGroup Space group.
1886 // 67 - 70 Integer z Z value.
1887 // =============================================================================
1888 if (nSymOp < 0) {
1889 // Write out the unit cell.
1890 Crystal crystal = activeMolecularAssembly.getCrystal();
1891 if (crystal != null && !crystal.aperiodic()) {
1892 Crystal c = crystal.getUnitCell();
1893 if (lmn[0] > 0 || lmn[1] > 0 || lmn[2] > 0) {
1894 c.a = c.a * lmn[0];
1895 c.b = c.b * lmn[1];
1896 c.c = c.c * lmn[2];
1897 }
1898 bw.write(c.toCRYST1());
1899 }
1900 } else if (nSymOp == 0) {
1901 // Write a P1 cell.
1902 Crystal crystal = activeMolecularAssembly.getCrystal();
1903 if (crystal != null && !crystal.aperiodic()) {
1904 Crystal c = crystal.getUnitCell();
1905 Crystal p1 = new Crystal((lmn[0]>0)? c.a * lmn[0] : c.a, (lmn[1]>0)? c.b * lmn[1] : c.b, (lmn[2]>0)? c.c * lmn[2] : c.c, c.alpha, c.beta, c.gamma, "P1");
1906 bw.write(p1.toCRYST1());
1907 }
1908 }
1909 // =============================================================================
1910 // The SSBOND record identifies each disulfide bond in protein and polypeptide
1911 // structures by identifying the two residues involved in the bond.
1912 // The disulfide bond distance is included after the symmetry operations at
1913 // the end of the SSBOND record.
1914 //
1915 // 8 - 10 Integer serNum Serial number.
1916 // 12 - 14 LString(3) "CYS" Residue name.
1917 // 16 Character chainID1 Chain identifier.
1918 // 18 - 21 Integer seqNum1 Residue sequence number.
1919 // 22 AChar icode1 Insertion code.
1920 // 26 - 28 LString(3) "CYS" Residue name.
1921 // 30 Character chainID2 Chain identifier.
1922 // 32 - 35 Integer seqNum2 Residue sequence number.
1923 // 36 AChar icode2 Insertion code.
1924 // 60 - 65 SymOP sym1 Symmetry oper for 1st resid
1925 // 67 - 72 SymOP sym2 Symmetry oper for 2nd resid
1926 // 74 – 78 Real(5.2) Length Disulfide bond distance
1927 //
1928 // If SG of cysteine is disordered then there are possible alternate linkages.
1929 // wwPDB practice is to put together all possible SSBOND records. This is
1930 // problematic because the alternate location identifier is not specified in
1931 // the SSBOND record.
1932 // =============================================================================
1933 int serNum = 1;
1934 Polymer[] polymers = activeMolecularAssembly.getChains();
1935 if (polymers != null) {
1936 for (Polymer polymer : polymers) {
1937 List<Residue> residues = polymer.getResidues();
1938 for (Residue residue : residues) {
1939 if (residue.getName().equalsIgnoreCase("CYS")) {
1940 List<Atom> cysAtoms = residue.getAtomList().stream()
1941 .filter(a -> !atomExclusions.contains(a)).toList();
1942 Atom SG1 = null;
1943 for (Atom atom : cysAtoms) {
1944 String atName = atom.getName().toUpperCase();
1945 if (atName.equals("SG") || atName.equals("SH")
1946 || atom.getAtomType().atomicNumber == 16) {
1947 SG1 = atom;
1948 break;
1949 }
1950 }
1951 List<Bond> bonds = SG1.getBonds();
1952 for (Bond bond : bonds) {
1953 Atom SG2 = bond.get1_2(SG1);
1954 if (SG2.getAtomType().atomicNumber == 16 && !atomExclusions.contains(SG2)) {
1955 if (SG1.getIndex() < SG2.getIndex()) {
1956 bond.energy(false);
1957 bw.write(format("SSBOND %3d CYS %1s %4s CYS %1s %4s %36s %5.2f\n", serNum++,
1958 SG1.getChainID().toString(), Hybrid36.encode(4, SG1.getResidueNumber()),
1959 SG2.getChainID().toString(), Hybrid36.encode(4, SG2.getResidueNumber()), "",
1960 bond.getValue()));
1961 }
1962 }
1963 }
1964 }
1965 }
1966 }
1967 }
1968 // =============================================================================
1969 //
1970 // 7 - 11 Integer serial Atom serial number.
1971 // 13 - 16 Atom name Atom name.
1972 // 17 Character altLoc Alternate location indicator.
1973 // 18 - 20 Residue name resName Residue name.
1974 // 22 Character chainID Chain identifier.
1975 // 23 - 26 Integer resSeq Residue sequence number.
1976 // 27 AChar iCode Code for insertion of residues.
1977 // 31 - 38 Real(8.3) x Orthogonal coordinates for X in Angstroms.
1978 // 39 - 46 Real(8.3) y Orthogonal coordinates for Y in Angstroms.
1979 // 47 - 54 Real(8.3) z Orthogonal coordinates for Z in Angstroms.
1980 // 55 - 60 Real(6.2) occupancy Occupancy.
1981 // 61 - 66 Real(6.2) tempFactor Temperature factor.
1982 // 77 - 78 LString(2) element Element symbol, right-justified.
1983 // 79 - 80 LString(2) charge Charge on the atom.
1984 // =============================================================================
1985 // 1 2 3 4 5 6 7
1986 // 123456789012345678901234567890123456789012345678901234567890123456789012345678
1987 // ATOM 1 N ILE A 16 60.614 71.140 -10.592 1.00 7.38 N
1988 // ATOM 2 CA ILE A 16 60.793 72.149 -9.511 1.00 6.91 C
1989 MolecularAssembly[] molecularAssemblies = this.getMolecularAssemblyArray();
1990 int serial = 1;
1991 if (nSymOp > 0) {
1992 serial = serialP1;
1993 }
1994
1995 // Loop over biomolecular chains
1996 if (polymers != null) {
1997 for (Polymer polymer : polymers) {
1998 currentSegID = polymer.getName();
1999 currentChainID = polymer.getChainID();
2000 sb.setCharAt(21, currentChainID);
2001 // Loop over residues
2002 List<Residue> residues = polymer.getResidues();
2003 for (Residue residue : residues) {
2004 String resName = residue.getName();
2005 if (resName.length() > 3) {
2006 resName = resName.substring(0, 3);
2007 }
2008 int resID = residue.getResidueNumber();
2009 sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
2010 sb.replace(22, 26, format("%4s", Hybrid36.encode(4, resID)));
2011 // Loop over atoms
2012 List<Atom> residueAtoms = residue.getAtomList().stream()
2013 .filter(a -> !atomExclusions.contains(a)).collect(Collectors.toList());
2014 boolean altLocFound = false;
2015 for (Atom atom : residueAtoms) {
2016 if (mutate) {
2017 for (Mutation mtn : mutations) {
2018 if (resID == mtn.resID) {
2019 if (residue.getBackboneAtoms().contains(atom)) {
2020 logger.info(format(" MUTATION atom is %d chain %s",serial, currentChainID));
2021 }
2022 }
2023 }
2024 }
2025 writeAtom(atom, serial++, sb, anisouSB, bw);
2026 Character altLoc = atom.getAltLoc();
2027 if (altLoc != null && !altLoc.equals(' ')) {
2028 altLocFound = true;
2029 }
2030 }
2031 // Write out alternate conformers
2032 if (altLocFound) {
2033 for (int ma = 1; ma < molecularAssemblies.length; ma++) {
2034 MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
2035 Polymer altPolymer = altMolecularAssembly.getPolymer(currentChainID, currentSegID, false);
2036 Residue altResidue = altPolymer.getResidue(resName, resID, false, Residue.ResidueType.AA);
2037 if (altResidue == null) {
2038 resName = AminoAcid3.UNK.name();
2039 altResidue = altPolymer.getResidue(resName, resID, false, Residue.ResidueType.AA);
2040 }
2041 residueAtoms = altResidue.getAtomList().stream()
2042 .filter(a -> !atomExclusions.contains(a)).collect(Collectors.toList());
2043 for (Atom atom : residueAtoms) {
2044 if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc()
2045 .equals('A')) {
2046 sb.replace(17, 20, padLeft(atom.getResidueName().toUpperCase(), 3));
2047 writeAtom(atom, serial++, sb, anisouSB, bw);
2048 }
2049 }
2050 }
2051 }
2052 }
2053 terSB.replace(6, 11, format("%5s", Hybrid36.encode(5, serial++)));
2054 terSB.replace(12, 16, " ");
2055 terSB.replace(16, 26, sb.substring(16, 26));
2056 bw.write(terSB.toString());
2057 bw.newLine();
2058 }
2059 }
2060 sb.replace(0, 6, "HETATM");
2061 sb.setCharAt(21, 'A');
2062
2063 Character chainID = 'A';
2064 if (polymers != null) {
2065 chainID = polymers[0].getChainID();
2066 }
2067 activeMolecularAssembly.setChainIDAndRenumberMolecules(chainID);
2068
2069 // Loop over molecules, ions and then water.
2070 List<MSNode> molecules = activeMolecularAssembly.getMolecules();
2071 int numMolecules = molecules.size();
2072 for (int i = 0; i < numMolecules; i++) {
2073 Molecule molecule = (Molecule) molecules.get(i);
2074 chainID = molecule.getChainID();
2075 sb.setCharAt(21, chainID);
2076 String resName = molecule.getResidueName();
2077 int resID = molecule.getResidueNumber();
2078 if (resName.length() > 3) {
2079 resName = resName.substring(0, 3);
2080 }
2081 sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
2082 sb.replace(22, 26, format("%4s", Hybrid36.encode(4, resID)));
2083 // List<Atom> moleculeAtoms = molecule.getAtomList();
2084 List<Atom> moleculeAtoms = molecule.getAtomList().stream()
2085 .filter(a -> !atomExclusions.contains(a)).collect(Collectors.toList());
2086 boolean altLocFound = false;
2087 for (Atom atom : moleculeAtoms) {
2088 writeAtom(atom, serial++, sb, anisouSB, bw);
2089 Character altLoc = atom.getAltLoc();
2090 if (altLoc != null && !altLoc.equals(' ')) {
2091 altLocFound = true;
2092 }
2093 }
2094 // Write out alternate conformers
2095 if (altLocFound) {
2096 for (int ma = 1; ma < molecularAssemblies.length; ma++) {
2097 MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
2098 MSNode altmolecule = altMolecularAssembly.getMolecules().get(i);
2099 moleculeAtoms = altmolecule.getAtomList();
2100 for (Atom atom : moleculeAtoms) {
2101 if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc()
2102 .equals('A')) {
2103 writeAtom(atom, serial++, sb, anisouSB, bw);
2104 }
2105 }
2106 }
2107 }
2108 }
2109
2110 List<MSNode> ions = activeMolecularAssembly.getIons();
2111 for (int i = 0; i < ions.size(); i++) {
2112 Molecule ion = (Molecule) ions.get(i);
2113 chainID = ion.getChainID();
2114 sb.setCharAt(21, chainID);
2115 String resName = ion.getResidueName();
2116 int resID = ion.getResidueNumber();
2117 if (resName.length() > 3) {
2118 resName = resName.substring(0, 3);
2119 }
2120 sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
2121 sb.replace(22, 26, format("%4s", Hybrid36.encode(4, resID)));
2122 // List<Atom> ionAtoms = ion.getAtomList();
2123 List<Atom> ionAtoms = ion.getAtomList().stream().filter(a -> !atomExclusions.contains(a))
2124 .collect(Collectors.toList());
2125 boolean altLocFound = false;
2126 for (Atom atom : ionAtoms) {
2127 writeAtom(atom, serial++, sb, anisouSB, bw);
2128 Character altLoc = atom.getAltLoc();
2129 if (altLoc != null && !altLoc.equals(' ')) {
2130 altLocFound = true;
2131 }
2132 }
2133 // Write out alternate conformers
2134 if (altLocFound) {
2135 for (int ma = 1; ma < molecularAssemblies.length; ma++) {
2136 MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
2137 MSNode altion = altMolecularAssembly.getIons().get(i);
2138 ionAtoms = altion.getAtomList();
2139 for (Atom atom : ionAtoms) {
2140 if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc()
2141 .equals('A')) {
2142 writeAtom(atom, serial++, sb, anisouSB, bw);
2143 }
2144 }
2145 }
2146 }
2147 }
2148
2149 List<MSNode> water = activeMolecularAssembly.getWater();
2150 for (int i = 0; i < water.size(); i++) {
2151 Molecule wat = (Molecule) water.get(i);
2152 chainID = wat.getChainID();
2153 sb.setCharAt(21, chainID);
2154 String resName = wat.getResidueName();
2155 int resID = wat.getResidueNumber();
2156 if (resName.length() > 3) {
2157 resName = resName.substring(0, 3);
2158 }
2159 sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
2160 sb.replace(22, 26, format("%4s", Hybrid36.encode(4, resID)));
2161 List<Atom> waterAtoms = wat.getAtomList().stream().filter(a -> !atomExclusions.contains(a))
2162 .collect(Collectors.toList());
2163 boolean altLocFound = false;
2164 for (Atom atom : waterAtoms) {
2165 writeAtom(atom, serial++, sb, anisouSB, bw);
2166 Character altLoc = atom.getAltLoc();
2167 if (altLoc != null && !altLoc.equals(' ')) {
2168 altLocFound = true;
2169 }
2170 }
2171 // Write out alternate conformers
2172 if (altLocFound) {
2173 for (int ma = 1; ma < molecularAssemblies.length; ma++) {
2174 MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
2175 MSNode altwater = altMolecularAssembly.getWater().get(i);
2176 waterAtoms = altwater.getAtomList();
2177 for (Atom atom : waterAtoms) {
2178 if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc()
2179 .equals('A')) {
2180 writeAtom(atom, serial++, sb, anisouSB, bw);
2181 }
2182 }
2183 }
2184 }
2185 }
2186
2187 if (model != null) {
2188 bw.write("ENDMDL");
2189 bw.newLine();
2190 }
2191
2192 if (writeEnd) {
2193 bw.write("END");
2194 bw.newLine();
2195 }
2196
2197 if (nSymOp >= 0) {
2198 serialP1 = serial;
2199 }
2200
2201 } catch (Exception e) {
2202 String message = "Exception writing to file: " + saveFile;
2203 logger.log(Level.WARNING, message, e);
2204 return false;
2205 }
2206 return true;
2207 }
2208
2209 /**
2210 * {@inheritDoc}
2211 *
2212 * <p>Write out the Atomic information in PDB format.
2213 */
2214 @Override
2215 public boolean writeFile(File saveFile, boolean append) {
2216 return writeFile(saveFile, append, false, true);
2217 }
2218
2219 public boolean writeFile(File saveFile, boolean append, String[] extraLines) {
2220 return writeFile(saveFile, append, Collections.emptySet(), false, !append, extraLines);
2221 }
2222
2223 /**
2224 * Writes out the atomic information in PDB format.
2225 *
2226 * @param saveFile The file to save information to.
2227 * @param append True if the current data should be appended to the saveFile (as in arc
2228 * files).
2229 * @param versioning True if the saveFile should be versioned. False if the saveFile should be
2230 * overwritten.
2231 * @return Success of writing.
2232 */
2233 public boolean writeFile(File saveFile, boolean append, boolean versioning) {
2234 return writeFile(saveFile, append, Collections.emptySet(), true, versioning);
2235 }
2236
2237 /**
2238 * writeFileWithHeader.
2239 *
2240 * @param saveFile a {@link java.io.File} object.
2241 * @param header a {@link java.lang.String} object.
2242 * @param append a boolean.
2243 * @return a boolean.
2244 */
2245 public boolean writeFileWithHeader(File saveFile, String header, boolean append) {
2246 if (standardizeAtomNames) {
2247 logger.info(" Setting atom names to PDB standard.");
2248 renameAtomsToPDBStandard(activeMolecularAssembly);
2249 }
2250 activeMolecularAssembly.setFile(saveFile);
2251 activeMolecularAssembly.setName(saveFile.getName());
2252
2253 try (FileWriter fw = new FileWriter(saveFile, append); BufferedWriter bw = new BufferedWriter(
2254 fw)) {
2255 bw.write(header);
2256 bw.newLine();
2257 } catch (Exception e) {
2258 String message = " Exception writing to file: " + saveFile;
2259 logger.log(Level.WARNING, message, e);
2260 return false;
2261 }
2262 if (writeFile(saveFile, true)) {
2263 logger.log(Level.INFO, " Wrote PDB to file {0}", saveFile.getPath());
2264 return true;
2265 } else {
2266 logger.log(Level.INFO, " Error writing to file {0}", saveFile.getPath());
2267 return false;
2268 }
2269 }
2270
2271 /**
2272 * writeFileWithHeader.
2273 *
2274 * @param saveFile a {@link java.io.File} object.
2275 * @param header a {@link java.lang.String} object.
2276 * @return a boolean.
2277 */
2278 public boolean writeFileWithHeader(File saveFile, String header) {
2279 return writeFileWithHeader(saveFile, header, true);
2280 }
2281
2282 /**
2283 * writeFileWithHeader.
2284 *
2285 * @param saveFile a {@link java.io.File} object.
2286 * @param header a {@link java.lang.StringBuilder} object.
2287 * @return a boolean.
2288 */
2289 public boolean writeFileWithHeader(File saveFile, StringBuilder header) {
2290 return writeFileWithHeader(saveFile, header.toString());
2291 }
2292
2293 /**
2294 * Get unique SegID for possibly duplicate chain IDs.
2295 *
2296 * @param c chain ID just read.
2297 * @return a unique segID.
2298 */
2299 private String getExistingSegID(Character c) {
2300 if (c.equals(' ')) {
2301 c = 'A';
2302 }
2303
2304 List<String> segIDs = segidMap.get(c);
2305 if (segIDs != null) {
2306 String segID = segIDs.get(0);
2307 if (segIDs.size() > 1) {
2308 logger.log(Level.INFO, format(" Multiple SegIDs for to chain %s; using %s.", c, segID));
2309 }
2310 return segID;
2311 } else {
2312 logger.log(Level.INFO, format(" Creating SegID for to chain %s", c));
2313 return getSegID(c);
2314 }
2315 }
2316
2317 /**
2318 * Convert possibly duplicate chain IDs into unique segIDs.
2319 *
2320 * @param c chain ID just read.
2321 * @return a unique segID.
2322 */
2323 private String getSegID(Character c) {
2324 if (c.equals(' ')) {
2325 c = 'A';
2326 }
2327
2328 // If the chain ID has not changed, return the existing segID.
2329 if (c.equals(currentChainID)) {
2330 return currentSegID;
2331 }
2332
2333 // Loop through existing segIDs to find the first one that is unused.
2334 int count = 0;
2335 for (String segID : segIDs) {
2336 if (segID.endsWith(c.toString())) {
2337 count++;
2338 }
2339 }
2340
2341 // If the count is greater than 0, then append it.
2342 String newSegID;
2343 if (count == 0) {
2344 newSegID = c.toString();
2345 } else {
2346 newSegID = count + c.toString();
2347 }
2348 segIDs.add(newSegID);
2349 currentChainID = c;
2350 currentSegID = newSegID;
2351
2352 if (segidMap.containsKey(c)) {
2353 segidMap.get(c).add(newSegID);
2354 } else {
2355 List<String> newChainList = new ArrayList<>();
2356 newChainList.add(newSegID);
2357 segidMap.put(c, newChainList);
2358 }
2359
2360 return newSegID;
2361 }
2362
2363 /**
2364 * writeAtom
2365 *
2366 * @param atom a {@link ffx.potential.bonded.Atom} object.
2367 * @param serial a int.
2368 * @param sb a {@link java.lang.StringBuilder} object.
2369 * @param anisouSB a {@link java.lang.StringBuilder} object.
2370 * @param bw a {@link java.io.BufferedWriter} object.
2371 * @throws java.io.IOException if any.
2372 */
2373 private void writeAtom(Atom atom, int serial, StringBuilder sb, StringBuilder anisouSB,
2374 BufferedWriter bw) throws IOException {
2375 String name = atom.getName();
2376 int nameLength = name.length();
2377 if (nameLength > 4) {
2378 name = name.substring(0, 4);
2379 } else if (nameLength == 1) {
2380 name = name + " ";
2381 } else if (nameLength == 2) {
2382 if (atom.getAtomType().valence == 0) {
2383 name = name + " ";
2384 } else {
2385 name = name + " ";
2386 }
2387 }
2388 double[] xyz = vdwH ? atom.getRedXYZ() : atom.getXYZ(null);
2389 if (nSymOp >= 0) {
2390 Crystal crystal = activeMolecularAssembly.getCrystal().getUnitCell();
2391 SymOp symOp = crystal.spaceGroup.getSymOp(nSymOp);
2392 double[] newXYZ = new double[xyz.length];
2393 crystal.applySymOp(xyz, newXYZ, symOp);
2394 if (lValue > 0 || mValue > 0 || nValue > 0) {
2395 double[] translation = new double[] {lValue, mValue, nValue};
2396 crystal.getUnitCell().toCartesianCoordinates(translation, translation);
2397 newXYZ[0] += translation[0];
2398 newXYZ[1] += translation[1];
2399 newXYZ[2] += translation[2];
2400 }
2401 xyz = newXYZ;
2402 }
2403 sb.replace(6, 16, format("%5s " + padLeft(name.toUpperCase(), 4), Hybrid36.encode(5, serial)));
2404 Character altLoc = atom.getAltLoc();
2405 sb.setCharAt(16, Objects.requireNonNullElse(altLoc, ' '));
2406
2407
2408 /*
2409 * On the following code:
2410 * #1: StringBuilder.replace will allow for longer strings, expanding the StringBuilder's length if necessary.
2411 *
2412 * #2: sb was never re-initialized, so if there was overflow,
2413 * sb would continue to be > 80 characters long, resulting in broken PDB files
2414 *
2415 * #3: It may be wiser to have XYZ coordinates result in shutdown, not
2416 * truncation of coordinates. #4: Excessive B-factors aren't much of an
2417 * issue; if the B-factor is past 999.99, that's the difference between
2418 * "density extends to Venus" and "density extends to Pluto".
2419 */
2420 StringBuilder decimals = new StringBuilder();
2421 for (int i = 0; i < 3; i++) {
2422 try {
2423 decimals.append(StringUtils.fwFpDec(xyz[i], 8, 3));
2424 } catch (IllegalArgumentException ex) {
2425 String newValue = StringUtils.fwFpTrunc(xyz[i], 8, 3);
2426 logger.info(format(" XYZ %d coordinate %8.3f for atom %s "
2427 + "overflowed bounds of 8.3f string specified by PDB "
2428 + "format; truncating value to %s", i, xyz[i], atom, newValue));
2429 decimals.append(newValue);
2430 }
2431 }
2432 try {
2433 decimals.append(StringUtils.fwFpDec(atom.getOccupancy(), 6, 2));
2434 } catch (IllegalArgumentException ex) {
2435 logger.severe(
2436 format(" Occupancy %f for atom %s is impossible; " + "value must be between 0 and 1",
2437 atom.getOccupancy(), atom));
2438 }
2439 try {
2440 decimals.append(StringUtils.fwFpDec(atom.getTempFactor(), 6, 2));
2441 } catch (IllegalArgumentException ex) {
2442 String newValue = StringUtils.fwFpTrunc(atom.getTempFactor(), 6, 2);
2443 logger.info(format(" Atom temp factor %6.2f for atom %s overflowed "
2444 + "bounds of 6.2f string specified by PDB format; truncating " + "value to %s",
2445 atom.getTempFactor(), atom, newValue));
2446 decimals.append(newValue);
2447 }
2448 sb.replace(30, 66, decimals.toString());
2449
2450 name = Atom.ElementSymbol.values()[atom.getAtomicNumber() - 1].toString();
2451 name = name.toUpperCase();
2452 if (atom.isDeuterium()) {
2453 name = "D";
2454 }
2455 sb.replace(76, 78, padLeft(name, 2));
2456 sb.replace(78, 80, format("%2d", 0));
2457 bw.write(sb.toString());
2458 bw.newLine();
2459 // =============================================================================
2460 // 1 - 6 Record name "ANISOU"
2461 // 7 - 11 Integer serial Atom serial number.
2462 // 13 - 16 Atom name Atom name.
2463 // 17 Character altLoc Alternate location indicator
2464 // 18 - 20 Residue name resName Residue name.
2465 // 22 Character chainID Chain identifier.
2466 // 23 - 26 Integer resSeq Residue sequence number.
2467 // 27 AChar iCode Insertion code.
2468 // 29 - 35 Integer u[0][0] U(1,1)
2469 // 36 - 42 Integer u[1][1] U(2,2)
2470 // 43 - 49 Integer u[2][2] U(3,3)
2471 // 50 - 56 Integer u[0][1] U(1,2)
2472 // 57 - 63 Integer u[0][2] U(1,3)
2473 // 64 - 70 Integer u[1][2] U(2,3)
2474 // 77 - 78 LString(2) element Element symbol, right-justified.
2475 // 79 - 80 LString(2) charge Charge on the atom.
2476 // =============================================================================
2477 double[] anisou = atom.getAnisou(null);
2478 if (anisou != null) {
2479 anisouSB.replace(6, 80, sb.substring(6, 80));
2480 anisouSB.replace(28, 70,
2481 format("%7d%7d%7d%7d%7d%7d", (int) (anisou[0] * 1e4), (int) (anisou[1] * 1e4),
2482 (int) (anisou[2] * 1e4), (int) (anisou[3] * 1e4), (int) (anisou[4] * 1e4),
2483 (int) (anisou[5] * 1e4)));
2484 bw.write(anisouSB.toString());
2485 bw.newLine();
2486 }
2487 }
2488
2489 /** PDB records that are recognized. */
2490 private enum Record {
2491 ANISOU, ATOM, CONECT, CRYST1, DBREF, END, MODEL, ENDMDL, HELIX, HETATM, LINK, MTRIX1, MTRIX2, MTRIX3, MODRES, SEQRES, SHEET, SSBOND, REMARK
2492 }
2493
2494 /** Presently, VERSION3_3 is default, and VERSION3_2 is anything non-standard. */
2495 public enum PDBFileStandard {
2496 VERSION2_3, VERSION3_0, VERSION3_1, VERSION3_2, VERSION3_3
2497 }
2498
2499 public static class Mutation {
2500
2501 /** Residue ID of the residue to mutate. */
2502 final int resID;
2503 /** Residue name after mutation. */
2504 final String resName;
2505 /** Character for the chain ID of the residue that will be mutated. */
2506 final char chainChar;
2507
2508 public Mutation(int resID, char chainChar, String newResName) {
2509 newResName = newResName.toUpperCase();
2510 if (newResName.length() != 3) {
2511 logger.log(Level.WARNING, format("Invalid mutation target: %s.", newResName));
2512 }
2513 int isAA = AminoAcidUtils.getAminoAcidNumber(newResName);
2514 int isNA = NucleicAcidUtils.getNucleicAcidNumber(newResName);
2515 if (isAA == -1 && isNA == -1) {
2516 logger.log(Level.WARNING, format("Invalid mutation target: %s.", newResName));
2517 }
2518 this.resID = resID;
2519 this.chainChar = chainChar;
2520 this.resName = newResName;
2521 }
2522 }
2523 }