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 mtn.origResName = resName;
833 resName = mtn.resName;
834 String atomName = name.toUpperCase();
835
836 int isAA = AminoAcidUtils.getAminoAcidNumber(resName);
837 int isNA = NucleicAcidUtils.getNucleicAcidNumber(resName);
838
839 if ((isNA != -1 && naBackboneNames.contains(atomName)) || (isAA != -1 && backboneNames.contains(atomName))) {
840 printAtom = true;
841 } else {
842 // grab pur-pur or pyr-pyr alchem. atoms
843 ArrayList<String> alchAtoms = mtn.getAlchemicalAtoms(false);
844 if (alchAtoms == null) {
845 // test to see if atom is involved in glycosyl torsion and if needs renaming (pyr-pur/pur-pyr)
846 String newName = mtn.isNonAlchemicalAtom(atomName);
847 if (newName != null) { // if not null -- use name it
848 printAtom = true;
849 if (newName.startsWith("~")) { // switch from purine to pyrmidine or v.v.
850 // switch name and include it as an alchemical atom
851 name = newName.substring(1);
852 logger.info(format(" DELETING atom %d %s of %s %d in chain %s", serial, atomName, resName, resSeq, chainID));
853 } else {
854 // replace name but do not include as an alchemical atom
855 name = newName;
856 }
857 doBreak = false;
858 } else if (!atomName.contains("'")) {
859 logger.info(format(" DELETING atom %d %s of %s %d in chain %s", serial, atomName, resName, resSeq, chainID));
860 doBreak = true;
861 } else {
862 printAtom = true;
863 doBreak = false;
864 }
865 } else {
866 if (alchAtoms.contains(atomName) && !atomName.contains("'")) {
867 logger.info(format(" DELETING atom %d %s of %s %d in chain %s", serial, atomName, resName, resSeq, chainID));
868 doBreak = true;
869 } else {
870 printAtom = true;
871 doBreak = false;
872 }
873 }
874 break;
875 }
876 }
877 }
878 if (doBreak) {
879 break;
880 }
881 }
882
883 if (constantPH) {
884 AminoAcid3 aa3 = AminoAcidUtils.getAminoAcid(resName.toUpperCase());
885 if (constantPHResidueMap.containsKey(aa3)) {
886 String atomName = name.toUpperCase();
887 AminoAcid3 aa3PH = constantPHResidueMap.get(aa3);
888 resName = aa3PH.name();
889 if (constantPhBackboneNames.contains(atomName)) {
890 logger.info(format(" %s-%d %s", resName, resSeq, atomName));
891 } else if (!atomName.startsWith("H")) {
892 logger.info(format(" %s-%d %s", resName, resSeq, atomName));
893 } else {
894 logger.info(format(" %s-%d %s skipped", resName, resSeq, atomName));
895 break;
896 }
897 }
898 } else if (rotamerTitration) {
899 AminoAcid3 aa3 = AminoAcidUtils.getAminoAcid(resName.toUpperCase());
900 if (rotamerResidueMap.containsKey(aa3) && resNumberList.contains(resSeq)) {
901 AminoAcid3 aa3rotamer = rotamerResidueMap.get(aa3);
902 resName = aa3rotamer.name();
903 }
904 }
905 d = new double[3];
906 d[0] = parseDouble(line.substring(30, 38).trim());
907 d[1] = parseDouble(line.substring(38, 46).trim());
908 d[2] = parseDouble(line.substring(46, 54).trim());
909 occupancy = 1.0;
910 tempFactor = 1.0;
911 try {
912 occupancy = parseDouble(line.substring(54, 60).trim());
913 tempFactor = parseDouble(line.substring(60, 66).trim());
914 } catch (NumberFormatException | StringIndexOutOfBoundsException e) {
915 // Use default values.
916 if (printMissingFields) {
917 logger.info(" Missing occupancy and b-factors set to 1.0.");
918 printMissingFields = false;
919 } else if (logger.isLoggable(Level.FINE)) {
920 logger.fine(" Missing occupancy and b-factors set to 1.0.");
921 }
922 }
923
924 double bfactor = properties.getDouble("reset-bfactors", -1.0);
925 if (bfactor >= 0.0) {
926 tempFactor = bfactor;
927 }
928
929 newAtom = new Atom(0, name, altLoc, d, resName, resSeq, chainID, occupancy, tempFactor, segID);
930
931 // Check if this is a modified residue.
932 if (modRes.containsKey(resName.toUpperCase())) {
933 newAtom.setModRes(true);
934 }
935 returnedAtom = (Atom) activeMolecularAssembly.addMSNode(newAtom);
936 if (returnedAtom != newAtom) {
937 // A previously added atom has been retained.
938 atoms.put(serial, returnedAtom);
939 if (logger.isLoggable(Level.FINE)) {
940 logger.fine(returnedAtom + " has been retained over\n" + newAtom);
941 }
942 } else {
943 // The new atom has been added.
944 atoms.put(serial, newAtom);
945 // Check if the newAtom took the xyzIndex of a previous alternate conformer.
946 if (newAtom.getIndex() == 0) {
947 newAtom.setXyzIndex(xyzIndex++);
948 }
949 if (printAtom) {
950 logger.info(newAtom.toString());
951 }
952 }
953 break;
954 }
955 break;
956 case HETATM:
957 // =============================================================================
958 // 1 - 6 Record name "HETATM"
959 // 7 - 11 Integer serial Atom serial number.
960 // 13 - 16 Atom name Atom name.
961 // 17 Character altLoc Alternate location indicator.
962 // 18 - 20 Residue name resName Residue name.
963 // 22 Character chainID Chain identifier.
964 // 23 - 26 Integer resSeq Residue sequence number.
965 // 27 AChar iCode Code for insertion of residues.
966 // 31 - 38 Real(8.3) x Orthogonal coordinates for X.
967 // 39 - 46 Real(8.3) y Orthogonal coordinates for Y.
968 // 47 - 54 Real(8.3) z Orthogonal coordinates for Z.
969 // 55 - 60 Real(6.2) occupancy Occupancy.
970 // 61 - 66 Real(6.2) tempFactor Temperature factor.
971 // 77 - 78 LString(2) element Element symbol; right-justified.
972 // 79 - 80 LString(2) charge Charge on the atom.
973 // =============================================================================
974 serial = Hybrid36.decode(5, line.substring(6, 11));
975 name = line.substring(12, 16).trim();
976 altLoc = line.substring(16, 17).toUpperCase().charAt(0);
977 if (!altLocs.contains(altLoc)) {
978 altLocs.add(altLoc);
979 }
980 if (!altLoc.equals(' ') && !altLoc.equals(currentAltLoc)) {
981 break;
982 }
983 // if (!altLoc.equals(' ') && !altLoc.equals('A') && !altLoc.equals(currentAltLoc)) {
984 // break;
985 //}
986 resName = line.substring(17, 20).trim();
987 chainID = line.substring(21, 22).charAt(0);
988 segID = getSegID(chainID);
989 resSeq = Hybrid36.decode(4, line.substring(22, 26));
990
991 char insertionCode = line.charAt(26);
992 if (insertionCode != ' ' && !containsInsCode) {
993 containsInsCode = true;
994 logger.warning(" FFX support for files with " + "insertion codes is experimental. "
995 + "Residues will be renumbered to " + "eliminate insertion codes (52A "
996 + "becomes 53, 53 becomes 54, etc)");
997 }
998
999 int offset = insertionCodeCount.getOrDefault(chainID, 0);
1000 String pdbResNum = format("%c%d%c", chainID, resSeq, insertionCode);
1001 if (!pdbToNewResMap.containsKey(pdbResNum)) {
1002 if (insertionCode != ' ') {
1003 ++offset;
1004 insertionCodeCount.put(chainID, offset);
1005 }
1006 resSeq += offset;
1007 if (offset != 0) {
1008 logger.info(
1009 format(" Chain %c " + "molecule %s-%s renumbered to %c %s-%d", chainID,
1010 pdbResNum.substring(1).trim(), resName, chainID, resName, resSeq));
1011 }
1012 String newNum = format("%c%d", chainID, resSeq);
1013 pdbToNewResMap.put(pdbResNum, newNum);
1014 } else {
1015 resSeq += offset;
1016 }
1017
1018 d = new double[3];
1019 d[0] = parseDouble(line.substring(30, 38).trim());
1020 d[1] = parseDouble(line.substring(38, 46).trim());
1021 d[2] = parseDouble(line.substring(46, 54).trim());
1022 occupancy = 1.0;
1023 tempFactor = 1.0;
1024 try {
1025 occupancy = parseDouble(line.substring(54, 60).trim());
1026 tempFactor = parseDouble(line.substring(60, 66).trim());
1027 } catch (NumberFormatException | StringIndexOutOfBoundsException e) {
1028 // Use default values.
1029 if (printMissingFields) {
1030 logger.info(" Missing occupancy and b-factors set to 1.0.");
1031 printMissingFields = false;
1032 } else if (logger.isLoggable(Level.FINE)) {
1033 logger.fine(" Missing occupancy and b-factors set to 1.0.");
1034 }
1035 }
1036
1037 double bfactor = properties.getDouble("reset-bfactors", -1.0);
1038 if (bfactor >= 0.0) {
1039 tempFactor = bfactor;
1040 }
1041
1042 newAtom = new Atom(0, name, altLoc, d, resName, resSeq, chainID, occupancy, tempFactor, segID);
1043 newAtom.setHetero(true);
1044 // Check if this is a modified residue.
1045 if (modRes.containsKey(resName.toUpperCase())) {
1046 newAtom.setModRes(true);
1047 }
1048 returnedAtom = (Atom) activeMolecularAssembly.addMSNode(newAtom);
1049 if (returnedAtom != newAtom) {
1050 // A previously added atom has been retained.
1051 atoms.put(serial, returnedAtom);
1052 if (logger.isLoggable(Level.FINE)) {
1053 logger.fine(returnedAtom + " has been retained over\n" + newAtom);
1054 }
1055 } else {
1056 // The new atom has been added.
1057 atoms.put(serial, newAtom);
1058 newAtom.setXyzIndex(xyzIndex++);
1059 }
1060 break;
1061 case CRYST1:
1062 // =============================================================================
1063 // The CRYST1 record presents the unit cell parameters, space group, and Z
1064 // value. If the structure was not determined by crystallographic means, CRYST1
1065 // simply provides the unitary values, with an appropriate REMARK.
1066 //
1067 // 7 - 15 Real(9.3) a a (Angstroms).
1068 // 16 - 24 Real(9.3) b b (Angstroms).
1069 // 25 - 33 Real(9.3) c c (Angstroms).
1070 // 34 - 40 Real(7.2) alpha alpha (degrees).
1071 // 41 - 47 Real(7.2) beta beta (degrees).
1072 // 48 - 54 Real(7.2) gamma gamma (degrees).
1073 // 56 - 66 LString sGroup Space group.
1074 // 67 - 70 Integer z Z value.
1075 // =============================================================================
1076 if (line.length() < 55) {
1077 logger.severe(" CRYST1 record is improperly formatted.");
1078 }
1079 double aaxis = parseDouble(line.substring(6, 15).trim());
1080 double baxis = parseDouble(line.substring(15, 24).trim());
1081 double caxis = parseDouble(line.substring(24, 33).trim());
1082 double alpha = parseDouble(line.substring(33, 40).trim());
1083 double beta = parseDouble(line.substring(40, 47).trim());
1084 double gamma = parseDouble(line.substring(47, 54).trim());
1085 int limit = min(line.length(), 66);
1086 String sg = line.substring(55, limit).trim();
1087 properties.addProperty("a-axis", aaxis);
1088 properties.addProperty("b-axis", baxis);
1089 properties.addProperty("c-axis", caxis);
1090 properties.addProperty("alpha", alpha);
1091 properties.addProperty("beta", beta);
1092 properties.addProperty("gamma", gamma);
1093 properties.addProperty("spacegroup", SpaceGroupInfo.pdb2ShortName(sg));
1094 break;
1095 case CONECT:
1096 // =============================================================================
1097 // 7 - 11 Integer serial Atom serial number
1098 // 12 - 16 Integer serial Serial number of bonded atom
1099 // 17 - 21 Integer serial Serial number of bonded atom
1100 // 22 - 26 Integer serial Serial number of bonded atom
1101 // 27 - 31 Integer serial Serial number of bonded atom
1102 //
1103 // CONECT records involving atoms for which the coordinates are not present
1104 // in the entry (e.g., symmetry-generated) are not given.
1105 // CONECT records involving atoms for which the coordinates are missing due
1106 // to disorder, are also not provided.
1107 // =============================================================================
1108 conects.add(line);
1109 break;
1110 case LINK:
1111 // =============================================================================
1112 // The LINK records specify connectivity between residues that is not implied by
1113 // the primary structure. Connectivity is expressed in terms of the atom names.
1114 // They also include the distance associated with each linkage following the
1115 // symmetry operations at the end of each record.
1116 // 13 - 16 Atom name1 Atom name.
1117 // 17 Character altLoc1 Alternate location indicator.
1118 // 18 - 20 Residue name resName1 Residue name.
1119 // 22 Character chainID1 Chain identifier.
1120 // 23 - 26 Integer resSeq1 Residue sequence number.
1121 // 27 AChar iCode1 Insertion code.
1122 // 43 - 46 Atom name2 Atom name.
1123 // 47 Character altLoc2 Alternate location indicator.
1124 // 48 - 50 Residue name resName2 Residue name.
1125 // 52 Character chainID2 Chain identifier.
1126 // 53 - 56 Integer resSeq2 Residue sequence number.
1127 // 57 AChar iCode2 Insertion code.
1128 // 60 - 65 SymOP sym1 Symmetry operator atom 1.
1129 // 67 - 72 SymOP sym2 Symmetry operator atom 2.
1130 // 74 – 78 Real(5.2) Length Link distance
1131 // =============================================================================
1132 char a1 = line.charAt(16);
1133 char a2 = line.charAt(46);
1134 if (a1 != a2) {
1135 // logger.info(format(" Ignoring LINK record as alternate locations do not match\n
1136 // %s.", line));
1137 break;
1138 }
1139 if (currentAltLoc == 'A') {
1140 if ((a1 == ' ' || a1 == 'A') && (a2 == ' ' || a2 == 'A')) {
1141 links.add(line);
1142 }
1143 } else if (a1 == currentAltLoc && a2 == currentAltLoc) {
1144 links.add(line);
1145 }
1146 break;
1147 case SSBOND:
1148 // =============================================================================
1149 // The SSBOND record identifies each disulfide bond in protein and polypeptide
1150 // structures by identifying the two residues involved in the bond.
1151 // The disulfide bond distance is included after the symmetry operations at
1152 // the end of the SSBOND record.
1153 //
1154 // 8 - 10 Integer serNum Serial number.
1155 // 12 - 14 LString(3) "CYS" Residue name.
1156 // 16 Character chainID1 Chain identifier.
1157 // 18 - 21 Integer seqNum1 Residue sequence number.
1158 // 22 AChar icode1 Insertion code.
1159 // 26 - 28 LString(3) "CYS" Residue name.
1160 // 30 Character chainID2 Chain identifier.
1161 // 32 - 35 Integer seqNum2 Residue sequence number.
1162 // 36 AChar icode2 Insertion code.
1163 // 60 - 65 SymOP sym1 Symmetry oper for 1st resid
1164 // 67 - 72 SymOP sym2 Symmetry oper for 2nd resid
1165 // 74 – 78 Real(5.2) Length Disulfide bond distance
1166 //
1167 // If SG of cysteine is disordered then there are possible alternate linkages.
1168 // wwPDB practice is to put together all possible SSBOND records. This is
1169 // problematic because the alternate location identifier is not specified in
1170 // the SSBOND record.
1171 //
1172 // Notes:
1173 // SSBOND records may be invalid if chain IDs are reused.
1174 // SSBOND records are applied by FFX to all conformers.
1175 // =============================================================================
1176 ssbonds.add(line);
1177 break;
1178 case HELIX:
1179 // =============================================================================
1180 // HELIX records are used to identify the position of helices in the molecule.
1181 // Helices are named, numbered, and classified by type. The residues where the
1182 // helix begins and ends are noted, as well as the total length.
1183 //
1184 // 8 - 10 Integer serNum Serial number of the helix. This starts
1185 // at 1 and increases incrementally.
1186 // 12 - 14 LString(3) helixID Helix identifier. In addition to a serial
1187 // number, each helix is given an
1188 // alphanumeric character helix identifier.
1189 // 16 - 18 Residue name initResName Name of the initial residue.
1190 // 20 Character initChainID Chain identifier for the chain containing
1191 // this helix.
1192 // 22 - 25 Integer initSeqNum Sequence number of the initial residue.
1193 // 26 AChar initICode Insertion code of the initial residue.
1194 // 28 - 30 Residue name endResName Name of the terminal residue of the helix.
1195 // 32 Character endChainID Chain identifier for the chain containing
1196 // this helix.
1197 // 34 - 37 Integer endSeqNum Sequence number of the terminal residue.
1198 // 38 AChar endICode Insertion code of the terminal residue.
1199 // 39 - 40 Integer helixClass Helix class (see below).
1200 // 41 - 70 String comment Comment about this helix.
1201 // 72 - 76 Integer length Length of this helix.
1202 //
1203 // CLASS NUMBER
1204 // TYPE OF HELIX (COLUMNS 39 - 40)
1205 // --------------------------------------------------------------
1206 // Right-handed alpha (default) 1
1207 // Right-handed omega 2
1208 // Right-handed pi 3
1209 // Right-handed gamma 4
1210 // Right-handed 3 - 10 5
1211 // Left-handed alpha 6
1212 // Left-handed omega 7
1213 // Left-handed gamma 8
1214 // 2 - 7 ribbon/helix 9
1215 // Polyproline 10
1216 // =============================================================================
1217 case SHEET:
1218 // =============================================================================
1219 // SHEET records are used to identify the position of sheets in the molecule.
1220 // Sheets are both named and numbered. The residues where the sheet begins and
1221 // ends are noted.
1222 //
1223 // 8 - 10 Integer strand Strand number which starts at 1 for each
1224 // strand within a sheet and increases by one.
1225 // 12 - 14 LString(3) sheetID Sheet identifier.
1226 // 15 - 16 Integer numStrands Number of strands in sheet.
1227 // 18 - 20 Residue name initResName Residue name of initial residue.
1228 // 22 Character initChainID Chain identifier of initial residue in strand.
1229 // 23 - 26 Integer initSeqNum Sequence number of initial residue in strand.
1230 // 27 AChar initICode Insertion code of initial residue in strand.
1231 // 29 - 31 Residue name endResName Residue name of terminal residue.
1232 // 33 Character endChainID Chain identifier of terminal residue.
1233 // 34 - 37 Integer endSeqNum Sequence number of terminal residue.
1234 // 38 AChar endICode Insertion code of terminal residue.
1235 // 39 - 40 Integer sense Sense of strand with respect to previous
1236 // strand in the sheet. 0 if first strand,
1237 // 1 if parallel,and -1 if anti-parallel.
1238 // 42 - 45 Atom curAtom Registration. Atom name in current strand.
1239 // 46 - 48 Residue name curResName Registration. Residue name in current strand
1240 // 50 Character curChainId Registration. Chain identifier in current strand.
1241 // 51 - 54 Integer curResSeq Registration. Residue sequence number
1242 // in current strand.
1243 // 55 AChar curICode Registration. Insertion code in current strand.
1244 // 57 - 60 Atom prevAtom Registration. Atom name in previous strand.
1245 // 61 - 63 Residue name prevResName Registration. Residue name in previous strand.
1246 // 65 Character prevChainId Registration. Chain identifier in previous strand.
1247 // 66 - 69 Integer prevResSeq Registration. Residue sequence number
1248 // in previous strand.
1249 // 70 AChar prevICode Registration. Insertion code in
1250 // previous strand.
1251 // =============================================================================
1252 structs.add(line);
1253 break;
1254 case MODEL: // Currently, no handling in initial read.
1255 break;
1256 case MTRIX1:
1257 // ================================================================================
1258 // MTRIXn (n = 1, 2, or 3) records present transformations expressing
1259 // non-crystallographic symmetry.
1260 // MTRIXn will appear only when such transformations are required to generate an
1261 // entire asymmetric unit,
1262 // such as a large viral structure.
1263 //
1264 // 8 - 10 Integer serial Serial number.
1265 // 11 - 20 Real(10.6) m[n][1] Mn1
1266 // 21 - 30 Real(10.6) m[n][2] Mn2
1267 // 31 - 40 Real(10.6) m[n][3] Mn3
1268 // 21 - 30 Real(10.6) v[n] Vn
1269 // 60 Integer iGiven 1 if coordinates for the representations which are
1270 // approximately related by the transformations of the
1271 // molecule are contained in the entry. Otherwise, blank.
1272 // =================================================================================
1273 StringBuilder MTRX1 = new StringBuilder(line.substring(11, 55));
1274 properties.addProperty("MTRIX1", MTRX1);
1275 break;
1276 case MTRIX2:
1277 StringBuilder MTRX2 = new StringBuilder(line.substring(11, 55));
1278 properties.addProperty("MTRIX2", MTRX2);
1279 break;
1280 case MTRIX3:
1281 StringBuilder MTRX3 = new StringBuilder(line.substring(11, 55));
1282 properties.addProperty("MTRIX3", MTRX3);
1283 break;
1284 case REMARK:
1285 remarkLines.add(line.trim());
1286 if (line.contains("Lambda:")) {
1287 Matcher m = lambdaPattern.matcher(line);
1288 if (m.find()) {
1289 lastReadLambda = Double.parseDouble(m.group(1));
1290 }
1291 }
1292 // =================================================================================
1293 // REMARK 350: presents all transformations, both crystallographic and non-crystallographic,
1294 // needed to generate the biomolecule. These transformations operate on the coordinates in the
1295 // entry. Both author and computational descriptions of assemblies are provided, if applicable.
1296 // For strict ncs case where more than one assembly presents in asymmetric unit, only one
1297 // chain with unit matrix will reported in REMARK 350, the other chain will be generated
1298 // by rotation and translation.
1299 //
1300 // 20 - 23 Integer serial Serial number.
1301 // 24 - 33 Real(10.6) m[n][1] Mn1
1302 // 34 - 43 Real(10.6) m[n][2] Mn2
1303 // 44 - 53 Real(10.6) m[n][3] Mn3
1304 // 59 - 68 Real(10.6) v[n] Vn
1305 // =================================================================================
1306 if (line.length() >= 68) {
1307 String remarkType = line.substring(7, 10).trim();
1308 if (remarkType.matches("\\d+") && parseInt(remarkType) == 350 && line.substring(13,
1309 18).equalsIgnoreCase("BIOMT")) {
1310 properties.addProperty("BIOMTn", new StringBuilder(line.substring(24, 68)));
1311 }
1312 }
1313 break;
1314 default:
1315 break;
1316 }
1317 line = br.readLine();
1318 }
1319
1320 } catch (FileNotFoundException fileNotFoundException) {
1321 logger.log(Level.SEVERE, " PDB file not found", fileNotFoundException);
1322 }
1323 }
1324 xyzIndex--;
1325 setFileRead(true);
1326 } catch (IOException e) {
1327 logger.exiting(PDBFilter.class.getName(), "readFile", e);
1328 return false;
1329 }
1330
1331 // Locate disulfide bonds; bond parameters are assigned below.
1332 List<Bond> ssBondList = locateDisulfideBonds(ssbonds, activeMolecularAssembly, pdbToNewResMap);
1333
1334 // Record the number of atoms read in from the PDB file before applying
1335 // algorithms that may build new atoms.
1336 int pdbAtoms = activeMolecularAssembly.getAtomArray().length;
1337
1338 // Build missing backbone atoms in loops.
1339 buildMissingResidues(xyzIndex, activeMolecularAssembly, seqRes, dbRef);
1340
1341 // Assign atom types. Missing side-chains atoms and missing hydrogen will be built in.
1342 bondList = assignAtomTypes(activeMolecularAssembly, fileStandard);
1343
1344 // Assign disulfide bonds parameters and log their creation.
1345 buildDisulfideBonds(ssBondList, activeMolecularAssembly, bondList);
1346
1347 // Finally, re-number the atoms if missing atoms were created.
1348 int currentN = activeMolecularAssembly.getAtomArray().length;
1349 boolean renumber = forceField.getBoolean("renumber-pdb", false);
1350 if (pdbAtoms != currentN) {
1351 logger.info(format(" Renumbering PDB file due to built atoms (%d vs %d)", currentN, pdbAtoms));
1352 numberAtoms(activeMolecularAssembly);
1353 } else if (renumber) {
1354 logger.info(" Renumbering PDB file due to renumber-pdb flag.");
1355 numberAtoms(activeMolecularAssembly);
1356 }
1357 return true;
1358 }
1359
1360 /** {@inheritDoc} */
1361 @Override
1362 public boolean readNext() {
1363 return readNext(false);
1364 }
1365
1366 /** {@inheritDoc} */
1367 @Override
1368 public boolean readNext(boolean resetPosition) {
1369 return readNext(resetPosition, false);
1370 }
1371
1372 /** {@inheritDoc} */
1373 @Override
1374 public boolean readNext(boolean resetPosition, boolean print) {
1375 return readNext(resetPosition, print, true);
1376 }
1377
1378 /** {@inheritDoc} */
1379 @Override
1380 public boolean readNext(boolean resetPosition, boolean print, boolean parse) {
1381 modelsRead = resetPosition ? 1 : modelsRead + 1;
1382 if (!parse) {
1383 if (print) {
1384 logger.info(format(" Skipped Model %d.", modelsRead));
1385 }
1386 return true;
1387 }
1388 remarkLines = new ArrayList<>(remarkLines.size());
1389 // ^ is beginning of line, \\s+ means "one or more whitespace", (\\d+) means match and capture
1390 // one or more digits.
1391 Pattern modelPatt = Pattern.compile("^MODEL\\s+(\\d+)");
1392 boolean eof = true;
1393 for (MolecularAssembly system : systems) {
1394 try {
1395 BufferedReader currentReader;
1396 if (readers.containsKey(system)) {
1397 currentReader = readers.get(system);
1398 try {
1399 if (!currentReader.ready()) {
1400 currentReader = new BufferedReader(new FileReader(readFile));
1401 // Mark the start of the file.
1402 currentReader.mark(0);
1403 readers.remove(system);
1404 readers.put(system, currentReader);
1405 } else if (resetPosition) {
1406 // If the BufferedReader has been opened, and reset is requested, reset the position.
1407 currentReader.reset();
1408 }
1409 } catch (Exception exception) {
1410 // If all structures in the PDB file have been read, the currentReader may have closed.
1411 // The try block will catch this case and reset to the beginning of the file.
1412 currentReader = new BufferedReader(new FileReader(readFile));
1413 // Mark the start of the file.
1414 currentReader.mark(0);
1415 readers.remove(system);
1416 readers.put(system, currentReader);
1417 }
1418 } else {
1419 currentReader = new BufferedReader(new FileReader(readFile));
1420 // Mark the start of the file.
1421 currentReader.mark(0);
1422 readers.put(system, currentReader);
1423 }
1424
1425 // Skip to appropriate model.
1426 String line = currentReader.readLine();
1427 while (line != null) {
1428 line = line.trim();
1429 Matcher m = modelPatt.matcher(line);
1430 if (m.find()) {
1431 int modelNum = parseInt(m.group(1));
1432 if (modelNum == modelsRead) {
1433 if (print) {
1434 logger.log(Level.INFO, format(" Reading model %d for %s", modelNum, currentFile));
1435 }
1436 eof = false;
1437 break;
1438 }
1439 }
1440 line = currentReader.readLine();
1441 }
1442 if (eof) {
1443 if (logger.isLoggable(Level.FINEST)) {
1444 logger.log(Level.FINEST, format("\n End of file reached for %s", readFile));
1445 }
1446 currentReader.close();
1447 return false;
1448 }
1449
1450 // Begin parsing the model.
1451 currentChainID = null;
1452 currentSegID = null;
1453 boolean modelDone = false;
1454 line = currentReader.readLine();
1455 while (line != null) {
1456 line = line.trim();
1457 String recID = line.substring(0, Math.min(6, line.length())).trim();
1458 try {
1459 Record record = Record.valueOf(recID);
1460 boolean hetatm = true;
1461 switch (record) {
1462 // =============================================================================
1463 //
1464 // 7 - 11 Integer serial Atom serial number.
1465 // 13 - 16 Atom name Atom name.
1466 // 17 Character altLoc Alternate location indicator.
1467 // 18 - 20 Residue name resName Residue name.
1468 // 22 Character chainID Chain identifier.
1469 // 23 - 26 Integer resSeq Residue sequence number.
1470 // 27 AChar iCode Code for insertion of residues.
1471 // 31 - 38 Real(8.3) x Orthogonal coordinates for X in
1472 // Angstroms.
1473 // 39 - 46 Real(8.3) y Orthogonal coordinates for Y in
1474 // Angstroms.
1475 // 47 - 54 Real(8.3) z Orthogonal coordinates for Z in
1476 // Angstroms.
1477 // 55 - 60 Real(6.2) occupancy Occupancy.
1478 // 61 - 66 Real(6.2) tempFactor Temperature factor.
1479 // 77 - 78 LString(2) element Element symbol, right-justified.
1480 // 79 - 80 LString(2) charge Charge on the atom.
1481 // =============================================================================
1482 // 1 2 3 4 5 6 7
1483 // 123456789012345678901234567890123456789012345678901234567890123456789012345678
1484 // ATOM 1 N ILE A 16 60.614 71.140 -10.592 1.00 7.38 N
1485 // ATOM 2 CA ILE A 16 60.793 72.149 -9.511 1.00 6.91 C
1486 case ATOM:
1487 hetatm = false;
1488 case HETATM:
1489 String name = line.substring(12, 16).trim();
1490 if (name.toUpperCase().contains("1H") || name.toUpperCase().contains("2H")
1491 || name.toUpperCase().contains("3H")) {
1492 // VERSION3_2 is presently just a placeholder for "anything non-standard".
1493 fileStandard = VERSION3_2;
1494 }
1495 Character altLoc = line.substring(16, 17).toUpperCase().charAt(0);
1496 if (!altLoc.equals(' ') && !altLoc.equals(currentAltLoc)) {
1497 break;
1498 }
1499 // if (!altLoc.equals(' ') && !altLoc.equals('A') && !altLoc.equals(currentAltLoc)) {
1500 // break;
1501 // }
1502 String resName = line.substring(17, 20).trim();
1503 Character chainID = line.substring(21, 22).charAt(0);
1504 String segID = getExistingSegID(chainID);
1505 int resSeq = Hybrid36.decode(4, line.substring(22, 26));
1506 double[] d = new double[3];
1507 d[0] = parseDouble(line.substring(30, 38).trim());
1508 d[1] = parseDouble(line.substring(38, 46).trim());
1509 d[2] = parseDouble(line.substring(46, 54).trim());
1510 double occupancy = 1.0;
1511 double tempFactor = 1.0;
1512 Atom newAtom = new Atom(0, name, altLoc, d, resName, resSeq, chainID, occupancy,
1513 tempFactor, segID);
1514 newAtom.setHetero(hetatm);
1515 // Check if this is a modified residue.
1516 if (modRes.containsKey(resName.toUpperCase())) {
1517 newAtom.setModRes(true);
1518 }
1519
1520 Atom returnedAtom = activeMolecularAssembly.findAtom(newAtom);
1521 if (returnedAtom != null) {
1522 returnedAtom.setXYZ(d);
1523 double[] retXYZ = new double[3];
1524 returnedAtom.getXYZ(retXYZ);
1525 } else {
1526 String message = format(" Could not find atom %s in assembly", newAtom);
1527 if (dieOnMissingAtom) {
1528 logger.severe(message);
1529 } else {
1530 logger.warning(message);
1531 }
1532 }
1533 break;
1534 case CRYST1:
1535 // =============================================================================
1536 // The CRYST1 record presents the unit cell parameters, space group, and Z
1537 // value. If the structure was not determined by crystallographic means, CRYST1
1538 // simply provides the unitary values, with an appropriate REMARK.
1539 //
1540 // 7 - 15 Real(9.3) a a (Angstroms).
1541 // 16 - 24 Real(9.3) b b (Angstroms).
1542 // 25 - 33 Real(9.3) c c (Angstroms).
1543 // 34 - 40 Real(7.2) alpha alpha (degrees).
1544 // 41 - 47 Real(7.2) beta beta (degrees).
1545 // 48 - 54 Real(7.2) gamma gamma (degrees).
1546 // 56 - 66 LString sGroup Space group.
1547 // 67 - 70 Integer z Z value.
1548 // =============================================================================
1549 logger.fine(" Crystal record found.");
1550 if (line.length() < 55) {
1551 logger.severe(" CRYST1 record is improperly formatted.");
1552 }
1553 double aaxis = parseDouble(line.substring(6, 15).trim());
1554 double baxis = parseDouble(line.substring(15, 24).trim());
1555 double caxis = parseDouble(line.substring(24, 33).trim());
1556 double alpha = parseDouble(line.substring(33, 40).trim());
1557 double beta = parseDouble(line.substring(40, 47).trim());
1558 double gamma = parseDouble(line.substring(47, 54).trim());
1559 int limit = min(line.length(), 66);
1560 String sg = line.substring(55, limit).trim();
1561 // properties.clearProperty("a-axis");
1562 // properties.clearProperty("b-axis");
1563 // properties.clearProperty("c-axis");
1564 // properties.clearProperty("alpha");
1565 // properties.clearProperty("beta");
1566 // properties.clearProperty("gamma");
1567 // properties.clearProperty("spacegroup");
1568 //
1569 // properties.addProperty("a-axis", aaxis);
1570 // properties.addProperty("b-axis", baxis);
1571 // properties.addProperty("c-axis", caxis);
1572 // properties.addProperty("alpha", alpha);
1573 // properties.addProperty("beta", beta);
1574 // properties.addProperty("gamma", gamma);
1575 // properties.addProperty("spacegroup", SpaceGroupInfo.pdb2ShortName(sg));
1576 Crystal crystal = activeMolecularAssembly.getCrystal();
1577 SpaceGroup spaceGroup = SpaceGroupDefinitions.spaceGroupFactory(sg);
1578 if (Objects.equals(crystal.spaceGroup.shortName, spaceGroup.shortName)) {
1579 crystal.changeUnitCellParameters(aaxis, baxis, caxis, alpha, beta, gamma);
1580 } else {
1581 // TODO: Handle changes in space groups... Means recalculating force field terms.
1582 logger.warning(format(" Original space group %s could not be changed to %s",
1583 crystal.spaceGroup.shortName, spaceGroup.shortName));
1584 }
1585 break;
1586 case ENDMDL:
1587 case END: // Technically speaking, END should be at the end of the file, not end of
1588 // the model.
1589 logger.log(Level.FINE, format(" Model %d successfully read", modelsRead));
1590 modelDone = true;
1591 break;
1592 case REMARK:
1593 remarkLines.add(line.trim());
1594 if (line.contains("Lambda:")) {
1595 Matcher m = lambdaPattern.matcher(line);
1596 if (m.find()) {
1597 lastReadLambda = Double.parseDouble(m.group(1));
1598 }
1599 }
1600 break;
1601 default:
1602 break;
1603 }
1604 } catch (Exception ex) {
1605 // Do nothing; it's not an ATOM/HETATM line.
1606 }
1607 if (modelDone) {
1608 break;
1609 }
1610 line = currentReader.readLine();
1611 }
1612 return true;
1613 } catch (IOException ex) {
1614 logger.info(
1615 format(" Exception in parsing frame %d of %s:" + " %s", modelsRead, system.toString(),
1616 ex));
1617 }
1618 }
1619 return false;
1620 }
1621
1622 /**
1623 * Specify the alternate location.
1624 *
1625 * @param molecularAssembly The MolecularAssembly to populate.
1626 * @param altLoc The alternate location to use.
1627 */
1628 public void setAltID(MolecularAssembly molecularAssembly, Character altLoc) {
1629 setMolecularSystem(molecularAssembly);
1630 currentAltLoc = altLoc;
1631 }
1632
1633 /**
1634 * Sets whether this PDBFilter should log each time it saves to a file.
1635 *
1636 * @param logWrites a boolean.
1637 */
1638 public void setLogWrites(boolean logWrites) {
1639 this.logWrites = logWrites;
1640 }
1641
1642 /**
1643 * setModelNumbering.
1644 *
1645 * @param modelsWritten the number of models written.
1646 */
1647 public void setModelNumbering(int modelsWritten) {
1648 this.modelsWritten = modelsWritten;
1649 }
1650
1651 public void setLMN(int[] lmn) {
1652 if(lmn[0] >= 1 && lmn[1] >= 1 && lmn[2] >= 1){
1653 this.lmn = lmn;
1654 }else{
1655 // Provided dimensions are not handled. Revert to P1.
1656 this.lmn = new int[]{1,1,1};
1657 }
1658 }
1659
1660 /**
1661 * setSymOp.
1662 *
1663 * @param symOp a int.
1664 */
1665 public void setSymOp(int symOp) {
1666 this.nSymOp = symOp;
1667 }
1668
1669 /**
1670 * Expand the current system to P1 during the save operation.
1671 *
1672 * @param file The file to write.
1673 * @return Return true on a successful write.
1674 */
1675 public boolean writeFileAsP1(File file) {
1676 // XYZ File First Line
1677 final int l = lmn[0];
1678 final int m = lmn[1];
1679 final int n = lmn[2];
1680 final int numReplicates = l * m * n;
1681 Crystal crystal = activeMolecularAssembly.getCrystal();
1682 int nSymOps = crystal.getUnitCell().spaceGroup.getNumberOfSymOps();
1683
1684 if (nSymOps == 1 && l <= 1 && m <= 1 && n <= 1) {
1685 // This is a P1 system.
1686 if (!writeFile(file, false)) {
1687 logger.info(format(" Save failed for %s", activeMolecularAssembly));
1688 return false;
1689 } else {
1690 return true;
1691 }
1692 } else {
1693 Polymer[] polymers = activeMolecularAssembly.getChains();
1694 int chainCount = 0;
1695 for (Polymer polymer : polymers) {
1696 Character chainID = Polymer.CHAIN_IDS.charAt(chainCount++);
1697 polymer.setChainID(chainID);
1698 polymer.setSegID(chainID.toString());
1699 }
1700 nSymOp = 0;
1701 logWrites = false;
1702 boolean writeEnd = false;
1703 if (!writeFile(file, false, false, writeEnd)) {
1704 logger.info(format(" Save failed for %s", activeMolecularAssembly));
1705 return false;
1706 } else {
1707 for (int i = 0; i < l; i++) {
1708 for (int j = 0; j < m; j++) {
1709 for (int k = 0; k < n; k++) {
1710 lValue = i;
1711 mValue = j;
1712 nValue = k;
1713 for (int iSym = 0; iSym < nSymOps; iSym++) {
1714 nSymOp = iSym;
1715 for (Polymer polymer : polymers) {
1716 Character chainID = Polymer.CHAIN_IDS.charAt(chainCount++);
1717 polymer.setChainID(chainID);
1718 polymer.setSegID(chainID.toString());
1719 }
1720 // If the last sym op to be written.
1721 writeEnd = iSym == nSymOps - 1 && i == l - 1 && j == m - 1 && k == n - 1;
1722 if (!writeFile(file, true, false, writeEnd)) {
1723 logger.info(format(" Save failed for %s", activeMolecularAssembly));
1724 return false;
1725 }
1726 }
1727 }
1728 }
1729 }
1730 }
1731
1732 // Reset the chainIDs.
1733 chainCount = 0;
1734 for (Polymer polymer : polymers) {
1735 Character chainID = Polymer.CHAIN_IDS.charAt(chainCount++);
1736 polymer.setChainID(chainID);
1737 polymer.setSegID(chainID.toString());
1738 }
1739
1740 }
1741
1742 return true;
1743 }
1744
1745 /**
1746 * writeFile
1747 *
1748 * @param saveFile a {@link java.io.File} object.
1749 * @param append Whether to append to saveFile (vs over-write).
1750 * @param printLinear Ignored (remains to present a different method signature).
1751 * @param writeEnd True if this is the final model.
1752 * @return Success of writing.
1753 */
1754 public boolean writeFile(File saveFile, boolean append, boolean printLinear, boolean writeEnd) {
1755 return writeFile(saveFile, append, Collections.emptySet(), writeEnd, true);
1756 }
1757
1758 /**
1759 * writeFile
1760 *
1761 * @param saveFile a {@link java.io.File} object to save to.
1762 * @param append Whether to append to saveFile (vs over-write).
1763 * @param toExclude A {@link java.util.Set} of {@link ffx.potential.bonded.Atom}s to exclude
1764 * from writing.
1765 * @param writeEnd True if this is the final model.
1766 * @param versioning True if the file being saved to should be versioned. False if the file
1767 * being saved to should be overwritten.
1768 * @return Success of writing.
1769 */
1770 public boolean writeFile(File saveFile, boolean append, Set<Atom> toExclude, boolean writeEnd,
1771 boolean versioning) {
1772 return writeFile(saveFile, append, toExclude, writeEnd, versioning, null);
1773 }
1774
1775 /**
1776 * writeFile
1777 *
1778 * @param saveFile a {@link java.io.File} object to save to.
1779 * @param append Whether to append to saveFile (vs over-write).
1780 * @param toExclude A {@link java.util.Set} of {@link ffx.potential.bonded.Atom}s to exclude
1781 * from writing.
1782 * @param writeEnd True if this is the final model.
1783 * @param versioning True if the file being saved to should be versioned. False if the file
1784 * being saved to should be overwritten.
1785 * @param extraLines Extra comment/header lines to write.
1786 * @return Success of writing.
1787 */
1788 public boolean writeFile(File saveFile, boolean append, Set<Atom> toExclude, boolean writeEnd,
1789 boolean versioning, String[] extraLines) {
1790 if (standardizeAtomNames) {
1791 logger.info(" Setting atom names to PDB standard.");
1792 renameAtomsToPDBStandard(activeMolecularAssembly);
1793 }
1794 final Set<Atom> atomExclusions = toExclude == null ? Collections.emptySet() : toExclude;
1795 if (saveFile == null) {
1796 return false;
1797 }
1798 if (vdwH) {
1799 logger.info(" Saving hydrogen to van der Waals centers instead of nuclear locations.");
1800 }
1801 if (nSymOp > -1) {
1802 logger.info(format(" Saving atoms using the symmetry operator:\n%s\n",
1803 activeMolecularAssembly.getCrystal().getUnitCell().spaceGroup.getSymOp(nSymOp)
1804 .toString()));
1805 }
1806
1807 // Create StringBuilders for ATOM, ANISOU and TER records that can be reused.
1808 StringBuilder sb = new StringBuilder("ATOM ");
1809 StringBuilder anisouSB = new StringBuilder("ANISOU");
1810 StringBuilder terSB = new StringBuilder("TER ");
1811 StringBuilder model = null;
1812 for (int i = 6; i < 80; i++) {
1813 sb.append(' ');
1814 anisouSB.append(' ');
1815 terSB.append(' ');
1816 }
1817
1818 File newFile = saveFile;
1819 if (!append) {
1820 if (versioning) {
1821 newFile = version(saveFile);
1822 }
1823 } else if (modelsWritten >= 0) {
1824 model = new StringBuilder(format("MODEL %-4d", ++modelsWritten));
1825 model.append(repeat(" ", 65));
1826 }
1827 activeMolecularAssembly.setFile(newFile);
1828 if (activeMolecularAssembly.getName() == null) {
1829 activeMolecularAssembly.setName(newFile.getName());
1830 }
1831 if (logWrites) {
1832 logger.log(Level.INFO, " Saving {0}", newFile.getName());
1833 }
1834
1835 try (FileWriter fw = new FileWriter(newFile, append);
1836 BufferedWriter bw = new BufferedWriter(fw)) {
1837 /*
1838 Will come before CRYST1 and ATOM records, but after anything
1839 written by writeFileWithHeader (particularly X-ray refinement statistics).
1840 */
1841 String[] headerLines = activeMolecularAssembly.getHeaderLines();
1842 for (String line : headerLines) {
1843 bw.write(format("%s\n", line));
1844 }
1845 if (extraLines != null) {
1846 if (rotamerTitration && extraLines[0].contains("REMARK")) {
1847 for (String line : extraLines) {
1848 bw.write(line + "\n");
1849 }
1850 } else {
1851 for (String line : extraLines) {
1852 bw.write(format("REMARK 999 %s\n", line));
1853 }
1854 }
1855 }
1856 if (model != null) {
1857 bw.write(model.toString());
1858 bw.newLine();
1859 }
1860 // =============================================================================
1861 // The CRYST1 record presents the unit cell parameters, space group, and Z
1862 // value. If the structure was not determined by crystallographic means, CRYST1
1863 // simply provides the unitary values, with an appropriate REMARK.
1864 //
1865 // 7 - 15 Real(9.3) a a (Angstroms).
1866 // 16 - 24 Real(9.3) b b (Angstroms).
1867 // 25 - 33 Real(9.3) c c (Angstroms).
1868 // 34 - 40 Real(7.2) alpha alpha (degrees).
1869 // 41 - 47 Real(7.2) beta beta (degrees).
1870 // 48 - 54 Real(7.2) gamma gamma (degrees).
1871 // 56 - 66 LString sGroup Space group.
1872 // 67 - 70 Integer z Z value.
1873 // =============================================================================
1874 if (nSymOp < 0) {
1875 // Write out the unit cell.
1876 Crystal crystal = activeMolecularAssembly.getCrystal();
1877 if (crystal != null && !crystal.aperiodic()) {
1878 Crystal c = crystal.getUnitCell();
1879 if (lmn[0] > 0 || lmn[1] > 0 || lmn[2] > 0) {
1880 c.a = c.a * lmn[0];
1881 c.b = c.b * lmn[1];
1882 c.c = c.c * lmn[2];
1883 }
1884 bw.write(c.toCRYST1());
1885 }
1886 } else if (nSymOp == 0) {
1887 // Write a P1 cell.
1888 Crystal crystal = activeMolecularAssembly.getCrystal();
1889 if (crystal != null && !crystal.aperiodic()) {
1890 Crystal c = crystal.getUnitCell();
1891 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");
1892 bw.write(p1.toCRYST1());
1893 }
1894 }
1895 // =============================================================================
1896 // The SSBOND record identifies each disulfide bond in protein and polypeptide
1897 // structures by identifying the two residues involved in the bond.
1898 // The disulfide bond distance is included after the symmetry operations at
1899 // the end of the SSBOND record.
1900 //
1901 // 8 - 10 Integer serNum Serial number.
1902 // 12 - 14 LString(3) "CYS" Residue name.
1903 // 16 Character chainID1 Chain identifier.
1904 // 18 - 21 Integer seqNum1 Residue sequence number.
1905 // 22 AChar icode1 Insertion code.
1906 // 26 - 28 LString(3) "CYS" Residue name.
1907 // 30 Character chainID2 Chain identifier.
1908 // 32 - 35 Integer seqNum2 Residue sequence number.
1909 // 36 AChar icode2 Insertion code.
1910 // 60 - 65 SymOP sym1 Symmetry oper for 1st resid
1911 // 67 - 72 SymOP sym2 Symmetry oper for 2nd resid
1912 // 74 – 78 Real(5.2) Length Disulfide bond distance
1913 //
1914 // If SG of cysteine is disordered then there are possible alternate linkages.
1915 // wwPDB practice is to put together all possible SSBOND records. This is
1916 // problematic because the alternate location identifier is not specified in
1917 // the SSBOND record.
1918 // =============================================================================
1919 int serNum = 1;
1920 Polymer[] polymers = activeMolecularAssembly.getChains();
1921 if (polymers != null) {
1922 for (Polymer polymer : polymers) {
1923 List<Residue> residues = polymer.getResidues();
1924 for (Residue residue : residues) {
1925 if (residue.getName().equalsIgnoreCase("CYS")) {
1926 List<Atom> cysAtoms = residue.getAtomList().stream()
1927 .filter(a -> !atomExclusions.contains(a)).toList();
1928 Atom SG1 = null;
1929 for (Atom atom : cysAtoms) {
1930 String atName = atom.getName().toUpperCase();
1931 if (atName.equals("SG") || atName.equals("SH")
1932 || atom.getAtomType().atomicNumber == 16) {
1933 SG1 = atom;
1934 break;
1935 }
1936 }
1937 List<Bond> bonds = SG1.getBonds();
1938 for (Bond bond : bonds) {
1939 Atom SG2 = bond.get1_2(SG1);
1940 if (SG2.getAtomType().atomicNumber == 16 && !atomExclusions.contains(SG2)) {
1941 if (SG1.getIndex() < SG2.getIndex()) {
1942 bond.energy(false);
1943 bw.write(format("SSBOND %3d CYS %1s %4s CYS %1s %4s %36s %5.2f\n", serNum++,
1944 SG1.getChainID().toString(), Hybrid36.encode(4, SG1.getResidueNumber()),
1945 SG2.getChainID().toString(), Hybrid36.encode(4, SG2.getResidueNumber()), "",
1946 bond.getValue()));
1947 }
1948 }
1949 }
1950 }
1951 }
1952 }
1953 }
1954 // =============================================================================
1955 //
1956 // 7 - 11 Integer serial Atom serial number.
1957 // 13 - 16 Atom name Atom name.
1958 // 17 Character altLoc Alternate location indicator.
1959 // 18 - 20 Residue name resName Residue name.
1960 // 22 Character chainID Chain identifier.
1961 // 23 - 26 Integer resSeq Residue sequence number.
1962 // 27 AChar iCode Code for insertion of residues.
1963 // 31 - 38 Real(8.3) x Orthogonal coordinates for X in Angstroms.
1964 // 39 - 46 Real(8.3) y Orthogonal coordinates for Y in Angstroms.
1965 // 47 - 54 Real(8.3) z Orthogonal coordinates for Z in Angstroms.
1966 // 55 - 60 Real(6.2) occupancy Occupancy.
1967 // 61 - 66 Real(6.2) tempFactor Temperature factor.
1968 // 77 - 78 LString(2) element Element symbol, right-justified.
1969 // 79 - 80 LString(2) charge Charge on the atom.
1970 // =============================================================================
1971 // 1 2 3 4 5 6 7
1972 // 123456789012345678901234567890123456789012345678901234567890123456789012345678
1973 // ATOM 1 N ILE A 16 60.614 71.140 -10.592 1.00 7.38 N
1974 // ATOM 2 CA ILE A 16 60.793 72.149 -9.511 1.00 6.91 C
1975 MolecularAssembly[] molecularAssemblies = this.getMolecularAssemblyArray();
1976 int serial = 1;
1977 if (nSymOp > 0) {
1978 serial = serialP1;
1979 }
1980
1981 // Loop over biomolecular chains
1982 if (polymers != null) {
1983 for (Polymer polymer : polymers) {
1984 currentSegID = polymer.getName();
1985 currentChainID = polymer.getChainID();
1986 sb.setCharAt(21, currentChainID);
1987 // Loop over residues
1988 List<Residue> residues = polymer.getResidues();
1989 for (Residue residue : residues) {
1990 String resName = residue.getName();
1991 if (resName.length() > 3) {
1992 resName = resName.substring(0, 3);
1993 }
1994 int resID = residue.getResidueNumber();
1995 sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
1996 sb.replace(22, 26, format("%4s", Hybrid36.encode(4, resID)));
1997 // Loop over atoms
1998 List<Atom> residueAtoms = residue.getAtomList().stream()
1999 .filter(a -> !atomExclusions.contains(a)).collect(Collectors.toList());
2000 boolean altLocFound = false;
2001 for (Atom atom : residueAtoms) {
2002 if (mutate) {
2003 for (Mutation mtn : mutations) {
2004 if (resID == mtn.resID) {
2005 ArrayList<String> alchAtoms = mtn.getAlchemicalAtoms(true);
2006 if (alchAtoms != null) {
2007 if (residue.getBackboneAtoms().contains(atom) && alchAtoms.contains(atom.getName())) {
2008 logger.info(format(" MUTATION atom is %d chain %s",serial, currentChainID));
2009 }
2010 } else {
2011 // treating pur-pyr or pyr-pur N9/N1 & C2/C4 as alchemical
2012 if (residue.getBackboneAtoms().contains(atom)) {
2013 logger.info(format(" MUTATION atom is %d chain %s",serial, currentChainID));
2014 }
2015 }
2016 }
2017 }
2018 }
2019 writeAtom(atom, serial++, sb, anisouSB, bw);
2020 Character altLoc = atom.getAltLoc();
2021 if (altLoc != null && !altLoc.equals(' ')) {
2022 altLocFound = true;
2023 }
2024 }
2025 // Write out alternate conformers
2026 if (altLocFound) {
2027 for (int ma = 1; ma < molecularAssemblies.length; ma++) {
2028 MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
2029 Polymer altPolymer = altMolecularAssembly.getPolymer(currentChainID, currentSegID, false);
2030 Residue altResidue = altPolymer.getResidue(resName, resID, false, Residue.ResidueType.AA);
2031 if (altResidue == null) {
2032 resName = AminoAcid3.UNK.name();
2033 altResidue = altPolymer.getResidue(resName, resID, false, Residue.ResidueType.AA);
2034 }
2035 residueAtoms = altResidue.getAtomList().stream()
2036 .filter(a -> !atomExclusions.contains(a)).collect(Collectors.toList());
2037 for (Atom atom : residueAtoms) {
2038 if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc()
2039 .equals('A')) {
2040 sb.replace(17, 20, padLeft(atom.getResidueName().toUpperCase(), 3));
2041 writeAtom(atom, serial++, sb, anisouSB, bw);
2042 }
2043 }
2044 }
2045 }
2046 }
2047 terSB.replace(6, 11, format("%5s", Hybrid36.encode(5, serial++)));
2048 terSB.replace(12, 16, " ");
2049 terSB.replace(16, 26, sb.substring(16, 26));
2050 bw.write(terSB.toString());
2051 bw.newLine();
2052 }
2053 }
2054 sb.replace(0, 6, "HETATM");
2055 sb.setCharAt(21, 'A');
2056
2057 Character chainID = 'A';
2058 if (polymers != null) {
2059 chainID = polymers[0].getChainID();
2060 }
2061 activeMolecularAssembly.setChainIDAndRenumberMolecules(chainID);
2062
2063 // Loop over molecules, ions and then water.
2064 List<MSNode> molecules = activeMolecularAssembly.getMolecules();
2065 int numMolecules = molecules.size();
2066 for (int i = 0; i < numMolecules; i++) {
2067 Molecule molecule = (Molecule) molecules.get(i);
2068 chainID = molecule.getChainID();
2069 sb.setCharAt(21, chainID);
2070 String resName = molecule.getResidueName();
2071 int resID = molecule.getResidueNumber();
2072 if (resName.length() > 3) {
2073 resName = resName.substring(0, 3);
2074 }
2075 sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
2076 sb.replace(22, 26, format("%4s", Hybrid36.encode(4, resID)));
2077 // List<Atom> moleculeAtoms = molecule.getAtomList();
2078 List<Atom> moleculeAtoms = molecule.getAtomList().stream()
2079 .filter(a -> !atomExclusions.contains(a)).collect(Collectors.toList());
2080 boolean altLocFound = false;
2081 for (Atom atom : moleculeAtoms) {
2082 writeAtom(atom, serial++, sb, anisouSB, bw);
2083 Character altLoc = atom.getAltLoc();
2084 if (altLoc != null && !altLoc.equals(' ')) {
2085 altLocFound = true;
2086 }
2087 }
2088 // Write out alternate conformers
2089 if (altLocFound) {
2090 for (int ma = 1; ma < molecularAssemblies.length; ma++) {
2091 MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
2092 MSNode altmolecule = altMolecularAssembly.getMolecules().get(i);
2093 moleculeAtoms = altmolecule.getAtomList();
2094 for (Atom atom : moleculeAtoms) {
2095 if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc()
2096 .equals('A')) {
2097 writeAtom(atom, serial++, sb, anisouSB, bw);
2098 }
2099 }
2100 }
2101 }
2102 }
2103
2104 List<MSNode> ions = activeMolecularAssembly.getIons();
2105 for (int i = 0; i < ions.size(); i++) {
2106 Molecule ion = (Molecule) ions.get(i);
2107 chainID = ion.getChainID();
2108 sb.setCharAt(21, chainID);
2109 String resName = ion.getResidueName();
2110 int resID = ion.getResidueNumber();
2111 if (resName.length() > 3) {
2112 resName = resName.substring(0, 3);
2113 }
2114 sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
2115 sb.replace(22, 26, format("%4s", Hybrid36.encode(4, resID)));
2116 // List<Atom> ionAtoms = ion.getAtomList();
2117 List<Atom> ionAtoms = ion.getAtomList().stream().filter(a -> !atomExclusions.contains(a))
2118 .collect(Collectors.toList());
2119 boolean altLocFound = false;
2120 for (Atom atom : ionAtoms) {
2121 writeAtom(atom, serial++, sb, anisouSB, bw);
2122 Character altLoc = atom.getAltLoc();
2123 if (altLoc != null && !altLoc.equals(' ')) {
2124 altLocFound = true;
2125 }
2126 }
2127 // Write out alternate conformers
2128 if (altLocFound) {
2129 for (int ma = 1; ma < molecularAssemblies.length; ma++) {
2130 MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
2131 MSNode altion = altMolecularAssembly.getIons().get(i);
2132 ionAtoms = altion.getAtomList();
2133 for (Atom atom : ionAtoms) {
2134 if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc()
2135 .equals('A')) {
2136 writeAtom(atom, serial++, sb, anisouSB, bw);
2137 }
2138 }
2139 }
2140 }
2141 }
2142
2143 List<MSNode> water = activeMolecularAssembly.getWater();
2144 for (int i = 0; i < water.size(); i++) {
2145 Molecule wat = (Molecule) water.get(i);
2146 chainID = wat.getChainID();
2147 sb.setCharAt(21, chainID);
2148 String resName = wat.getResidueName();
2149 int resID = wat.getResidueNumber();
2150 if (resName.length() > 3) {
2151 resName = resName.substring(0, 3);
2152 }
2153 sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
2154 sb.replace(22, 26, format("%4s", Hybrid36.encode(4, resID)));
2155 List<Atom> waterAtoms = wat.getAtomList().stream().filter(a -> !atomExclusions.contains(a))
2156 .collect(Collectors.toList());
2157 boolean altLocFound = false;
2158 for (Atom atom : waterAtoms) {
2159 writeAtom(atom, serial++, sb, anisouSB, bw);
2160 Character altLoc = atom.getAltLoc();
2161 if (altLoc != null && !altLoc.equals(' ')) {
2162 altLocFound = true;
2163 }
2164 }
2165 // Write out alternate conformers
2166 if (altLocFound) {
2167 for (int ma = 1; ma < molecularAssemblies.length; ma++) {
2168 MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
2169 MSNode altwater = altMolecularAssembly.getWater().get(i);
2170 waterAtoms = altwater.getAtomList();
2171 for (Atom atom : waterAtoms) {
2172 if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc()
2173 .equals('A')) {
2174 writeAtom(atom, serial++, sb, anisouSB, bw);
2175 }
2176 }
2177 }
2178 }
2179 }
2180
2181 if (model != null) {
2182 bw.write("ENDMDL");
2183 bw.newLine();
2184 }
2185
2186 if (writeEnd) {
2187 bw.write("END");
2188 bw.newLine();
2189 }
2190
2191 if (nSymOp >= 0) {
2192 serialP1 = serial;
2193 }
2194
2195 } catch (Exception e) {
2196 String message = "Exception writing to file: " + saveFile;
2197 logger.log(Level.WARNING, message, e);
2198 return false;
2199 }
2200 return true;
2201 }
2202
2203 /**
2204 * {@inheritDoc}
2205 *
2206 * <p>Write out the Atomic information in PDB format.
2207 */
2208 @Override
2209 public boolean writeFile(File saveFile, boolean append) {
2210 return writeFile(saveFile, append, false, true);
2211 }
2212
2213 public boolean writeFile(File saveFile, boolean append, String[] extraLines) {
2214 return writeFile(saveFile, append, Collections.emptySet(), false, !append, extraLines);
2215 }
2216
2217 /**
2218 * Writes out the atomic information in PDB format.
2219 *
2220 * @param saveFile The file to save information to.
2221 * @param append True if the current data should be appended to the saveFile (as in arc
2222 * files).
2223 * @param versioning True if the saveFile should be versioned. False if the saveFile should be
2224 * overwritten.
2225 * @return Success of writing.
2226 */
2227 public boolean writeFile(File saveFile, boolean append, boolean versioning) {
2228 return writeFile(saveFile, append, Collections.emptySet(), true, versioning);
2229 }
2230
2231 /**
2232 * writeFileWithHeader.
2233 *
2234 * @param saveFile a {@link java.io.File} object.
2235 * @param header a {@link java.lang.String} object.
2236 * @param append a boolean.
2237 * @return a boolean.
2238 */
2239 public boolean writeFileWithHeader(File saveFile, String header, boolean append) {
2240 if (standardizeAtomNames) {
2241 logger.info(" Setting atom names to PDB standard.");
2242 renameAtomsToPDBStandard(activeMolecularAssembly);
2243 }
2244 activeMolecularAssembly.setFile(saveFile);
2245 activeMolecularAssembly.setName(saveFile.getName());
2246
2247 try (FileWriter fw = new FileWriter(saveFile, append); BufferedWriter bw = new BufferedWriter(
2248 fw)) {
2249 bw.write(header);
2250 bw.newLine();
2251 } catch (Exception e) {
2252 String message = " Exception writing to file: " + saveFile;
2253 logger.log(Level.WARNING, message, e);
2254 return false;
2255 }
2256 if (writeFile(saveFile, true)) {
2257 logger.log(Level.INFO, " Wrote PDB to file {0}", saveFile.getPath());
2258 return true;
2259 } else {
2260 logger.log(Level.INFO, " Error writing to file {0}", saveFile.getPath());
2261 return false;
2262 }
2263 }
2264
2265 /**
2266 * writeFileWithHeader.
2267 *
2268 * @param saveFile a {@link java.io.File} object.
2269 * @param header a {@link java.lang.String} object.
2270 * @return a boolean.
2271 */
2272 public boolean writeFileWithHeader(File saveFile, String header) {
2273 return writeFileWithHeader(saveFile, header, true);
2274 }
2275
2276 /**
2277 * writeFileWithHeader.
2278 *
2279 * @param saveFile a {@link java.io.File} object.
2280 * @param header a {@link java.lang.StringBuilder} object.
2281 * @return a boolean.
2282 */
2283 public boolean writeFileWithHeader(File saveFile, StringBuilder header) {
2284 return writeFileWithHeader(saveFile, header.toString());
2285 }
2286
2287 /**
2288 * Get unique SegID for possibly duplicate chain IDs.
2289 *
2290 * @param c chain ID just read.
2291 * @return a unique segID.
2292 */
2293 private String getExistingSegID(Character c) {
2294 if (c.equals(' ')) {
2295 c = 'A';
2296 }
2297
2298 // If the chain ID has not changed, return the existing segID.
2299 if (c.equals(currentChainID)) {
2300 return currentSegID;
2301 } else {
2302 currentChainID = null;
2303 }
2304
2305 List<String> segIDs = segidMap.get(c);
2306 if (segIDs != null) {
2307 if (segIDs.size() > 1) {
2308 if (currentSegID == null) {
2309 currentChainID = c;
2310 currentSegID = segIDs.get(0);
2311 return segIDs.get(0);
2312 } else if (currentSegID.length() == 1) {
2313 currentChainID = c;
2314 currentSegID = segIDs.get(1);
2315 return segIDs.get(1);
2316 } else if (currentSegID.length() == 2) {
2317 String s = currentSegID.substring(0,1);
2318 int num = -2;
2319 try {
2320 num = Integer.parseInt(s);
2321 } catch (NumberFormatException e) {
2322 logger.severe(" SegID of length 2 does not start with an integer.");
2323 }
2324 currentChainID = c;
2325 currentSegID = segIDs.get(num+1);
2326 return segIDs.get(num+1);
2327 } else {
2328 logger.info(" Too many repeated chains. Using single letter for segID.");
2329 }
2330 }
2331 return segIDs.get(0);
2332 } else {
2333 logger.log(Level.INFO, format(" Creating SegID for to chain %s", c));
2334 return getSegID(c);
2335 }
2336 }
2337
2338 /**
2339 * Convert possibly duplicate chain IDs into unique segIDs.
2340 *
2341 * @param c chain ID just read.
2342 * @return a unique segID.
2343 */
2344 private String getSegID(Character c) {
2345 if (c.equals(' ')) {
2346 c = 'A';
2347 }
2348
2349 // If the chain ID has not changed, return the existing segID.
2350 if (c.equals(currentChainID)) {
2351 return currentSegID;
2352 }
2353
2354 // Loop through existing segIDs to find the first one that is unused.
2355 int count = 0;
2356 for (String segID : segIDs) {
2357 if (segID.endsWith(c.toString())) {
2358 count++;
2359 }
2360 }
2361
2362 // If the count is greater than 0, then append it.
2363 String newSegID;
2364 if (count == 0) {
2365 newSegID = c.toString();
2366 } else {
2367 newSegID = count + c.toString();
2368 }
2369 segIDs.add(newSegID);
2370 currentChainID = c;
2371 currentSegID = newSegID;
2372
2373 if (segidMap.containsKey(c)) {
2374 segidMap.get(c).add(newSegID);
2375 } else {
2376 List<String> newChainList = new ArrayList<>();
2377 newChainList.add(newSegID);
2378 segidMap.put(c, newChainList);
2379 }
2380
2381 return newSegID;
2382 }
2383
2384 /**
2385 * writeAtom
2386 *
2387 * @param atom a {@link ffx.potential.bonded.Atom} object.
2388 * @param serial a int.
2389 * @param sb a {@link java.lang.StringBuilder} object.
2390 * @param anisouSB a {@link java.lang.StringBuilder} object.
2391 * @param bw a {@link java.io.BufferedWriter} object.
2392 * @throws java.io.IOException if any.
2393 */
2394 private void writeAtom(Atom atom, int serial, StringBuilder sb, StringBuilder anisouSB,
2395 BufferedWriter bw) throws IOException {
2396 String name = atom.getName();
2397 int nameLength = name.length();
2398 if (nameLength > 4) {
2399 name = name.substring(0, 4);
2400 } else if (nameLength == 1) {
2401 name = name + " ";
2402 } else if (nameLength == 2) {
2403 if (atom.getAtomType().valence == 0) {
2404 name = name + " ";
2405 } else {
2406 name = name + " ";
2407 }
2408 }
2409 double[] xyz = vdwH ? atom.getRedXYZ() : atom.getXYZ(null);
2410 if (nSymOp >= 0) {
2411 Crystal crystal = activeMolecularAssembly.getCrystal().getUnitCell();
2412 SymOp symOp = crystal.spaceGroup.getSymOp(nSymOp);
2413 double[] newXYZ = new double[xyz.length];
2414 crystal.applySymOp(xyz, newXYZ, symOp);
2415 if (lValue > 0 || mValue > 0 || nValue > 0) {
2416 double[] translation = new double[] {lValue, mValue, nValue};
2417 crystal.getUnitCell().toCartesianCoordinates(translation, translation);
2418 newXYZ[0] += translation[0];
2419 newXYZ[1] += translation[1];
2420 newXYZ[2] += translation[2];
2421 }
2422 xyz = newXYZ;
2423 }
2424 sb.replace(6, 16, format("%5s " + padLeft(name.toUpperCase(), 4), Hybrid36.encode(5, serial)));
2425 Character altLoc = atom.getAltLoc();
2426 sb.setCharAt(16, Objects.requireNonNullElse(altLoc, ' '));
2427
2428
2429 /*
2430 * On the following code:
2431 * #1: StringBuilder.replace will allow for longer strings, expanding the StringBuilder's length if necessary.
2432 *
2433 * #2: sb was never re-initialized, so if there was overflow,
2434 * sb would continue to be > 80 characters long, resulting in broken PDB files
2435 *
2436 * #3: It may be wiser to have XYZ coordinates result in shutdown, not
2437 * truncation of coordinates. #4: Excessive B-factors aren't much of an
2438 * issue; if the B-factor is past 999.99, that's the difference between
2439 * "density extends to Venus" and "density extends to Pluto".
2440 */
2441 StringBuilder decimals = new StringBuilder();
2442 for (int i = 0; i < 3; i++) {
2443 try {
2444 decimals.append(StringUtils.fwFpDec(xyz[i], 8, 3));
2445 } catch (IllegalArgumentException ex) {
2446 String newValue = StringUtils.fwFpTrunc(xyz[i], 8, 3);
2447 logger.info(format(" XYZ %d coordinate %8.3f for atom %s "
2448 + "overflowed bounds of 8.3f string specified by PDB "
2449 + "format; truncating value to %s", i, xyz[i], atom, newValue));
2450 decimals.append(newValue);
2451 }
2452 }
2453 try {
2454 decimals.append(StringUtils.fwFpDec(atom.getOccupancy(), 6, 2));
2455 } catch (IllegalArgumentException ex) {
2456 logger.severe(
2457 format(" Occupancy %f for atom %s is impossible; " + "value must be between 0 and 1",
2458 atom.getOccupancy(), atom));
2459 }
2460 try {
2461 decimals.append(StringUtils.fwFpDec(atom.getTempFactor(), 6, 2));
2462 } catch (IllegalArgumentException ex) {
2463 String newValue = StringUtils.fwFpTrunc(atom.getTempFactor(), 6, 2);
2464 logger.info(format(" Atom temp factor %6.2f for atom %s overflowed "
2465 + "bounds of 6.2f string specified by PDB format; truncating " + "value to %s",
2466 atom.getTempFactor(), atom, newValue));
2467 decimals.append(newValue);
2468 }
2469 sb.replace(30, 66, decimals.toString());
2470
2471 name = Atom.ElementSymbol.values()[atom.getAtomicNumber() - 1].toString();
2472 name = name.toUpperCase();
2473 if (atom.isDeuterium()) {
2474 name = "D";
2475 }
2476 sb.replace(76, 78, padLeft(name, 2));
2477 sb.replace(78, 80, format("%2d", 0));
2478 bw.write(sb.toString());
2479 bw.newLine();
2480 // =============================================================================
2481 // 1 - 6 Record name "ANISOU"
2482 // 7 - 11 Integer serial Atom serial number.
2483 // 13 - 16 Atom name Atom name.
2484 // 17 Character altLoc Alternate location indicator
2485 // 18 - 20 Residue name resName Residue name.
2486 // 22 Character chainID Chain identifier.
2487 // 23 - 26 Integer resSeq Residue sequence number.
2488 // 27 AChar iCode Insertion code.
2489 // 29 - 35 Integer u[0][0] U(1,1)
2490 // 36 - 42 Integer u[1][1] U(2,2)
2491 // 43 - 49 Integer u[2][2] U(3,3)
2492 // 50 - 56 Integer u[0][1] U(1,2)
2493 // 57 - 63 Integer u[0][2] U(1,3)
2494 // 64 - 70 Integer u[1][2] U(2,3)
2495 // 77 - 78 LString(2) element Element symbol, right-justified.
2496 // 79 - 80 LString(2) charge Charge on the atom.
2497 // =============================================================================
2498 double[] anisou = atom.getAnisou(null);
2499 if (anisou != null) {
2500 anisouSB.replace(6, 80, sb.substring(6, 80));
2501 anisouSB.replace(28, 70,
2502 format("%7d%7d%7d%7d%7d%7d", (int) (anisou[0] * 1e4), (int) (anisou[1] * 1e4),
2503 (int) (anisou[2] * 1e4), (int) (anisou[3] * 1e4), (int) (anisou[4] * 1e4),
2504 (int) (anisou[5] * 1e4)));
2505 bw.write(anisouSB.toString());
2506 bw.newLine();
2507 }
2508 }
2509
2510 /** PDB records that are recognized. */
2511 private enum Record {
2512 ANISOU, ATOM, CONECT, CRYST1, DBREF, END, MODEL, ENDMDL, HELIX, HETATM, LINK, MTRIX1, MTRIX2, MTRIX3, MODRES, SEQRES, SHEET, SSBOND, REMARK
2513 }
2514
2515 /** Presently, VERSION3_3 is default, and VERSION3_2 is anything non-standard. */
2516 public enum PDBFileStandard {
2517 VERSION2_3, VERSION3_0, VERSION3_1, VERSION3_2, VERSION3_3
2518 }
2519
2520 public static class Mutation {
2521
2522 /** Residue ID of the residue to mutate. */
2523 final int resID;
2524 /** Residue name after mutation. */
2525 final String resName;
2526 /** Character for the chain ID of the residue that will be mutated. */
2527 final char chainChar;
2528 /** Residue name before mutation. */
2529 String origResName;
2530
2531 public Mutation(int resID, char chainChar, String newResName) {
2532 newResName = newResName.toUpperCase();
2533 if (newResName.length() != 3) {
2534 logger.log(Level.WARNING, format("Invalid mutation target: %s.", newResName));
2535 }
2536 int isAA = AminoAcidUtils.getAminoAcidNumber(newResName);
2537 int isNA = NucleicAcidUtils.getNucleicAcidNumber(newResName);
2538 if (isAA == -1 && isNA == -1) {
2539 logger.log(Level.WARNING, format("Invalid mutation target: %s.", newResName));
2540 }
2541 this.resID = resID;
2542 this.chainChar = chainChar;
2543 this.resName = newResName;
2544 }
2545
2546 /**
2547 * Check to see if an atom is involved in the mutated base's glycosyl torsion. If the mutation is a switch from
2548 * purine to pyrimidine or vice versa, it will return '~name', meaning the name should be replaced but to include it
2549 * as an alchemical atom.
2550 * @param atomName atom name to check
2551 * @return new name to use if it is involved in glycosyl torsion, null otherwise
2552 */
2553 public String isNonAlchemicalAtom(String atomName) {
2554 if (isWtPurine()) {
2555 if (atomName.equals("N9")) {
2556 if (isMtnPyrimidine()) {
2557 return "~N1";
2558 }
2559 return atomName;
2560 } else if (atomName.equals("C4")) {
2561 if (isMtnPyrimidine()) {
2562 return "~C2";
2563 }
2564 return atomName;
2565 }
2566 return null;
2567 }
2568
2569 if (isWtPyrimidine()) {
2570 if (atomName.equals("N1")) {
2571 if (isMtnPurine()) {
2572 // here
2573 return "~N9";
2574 }
2575 return atomName;
2576 } else if (atomName.equals("C2")) {
2577 if (isMtnPurine()) {
2578 // here
2579 return "~C4";
2580 }
2581 return atomName;
2582 }
2583 return null;
2584 }
2585
2586 return null;
2587 }
2588
2589 /**
2590 * Determines what atoms should be alchemical for a purine to purine or pyrimidine to pyrimidine mutation.
2591 * @param isWriting true if writing the pdb, false if reading the pdb
2592 * @return ArrayList of alchemical atoms, null if not a pur-pur or pyr-pyr mutation
2593 */
2594 public ArrayList<String> getAlchemicalAtoms(boolean isWriting) {
2595 // Log warning that the mutation input is the same residue and return nothing so prev. functionality is not changed
2596 if (resName.equals(origResName)) {
2597 logger.severe("Desired Mutation residue is the same as the original.");
2598 return null;
2599 }
2600
2601 boolean purpur;
2602 if (isMtnPurine() && isWtPurine()) {
2603 purpur = true;
2604 } else if (isMtnPyrimidine() && isWtPyrimidine()) {
2605 purpur = false;
2606 } else {
2607 // Return nothing so previous functionality is not changed
2608 return null;
2609 }
2610
2611 String res;
2612
2613 // look at the mutation residue if writing or the original (wild type) residue if reading
2614 if (isWriting) {
2615 res = resName;
2616 } else {
2617 res = origResName;
2618 }
2619
2620 ArrayList<String> list = new ArrayList<>();
2621 if (purpur) { // purine: either A to G or G to A
2622 if (res.equals("DAD") || res.equals("DA")) { // MTN is A
2623 list.add("N6");
2624 list.add("H61");
2625 list.add("H62");
2626 list.add("H2");
2627 } else { // MTN is G
2628 list.add("H1");
2629 list.add("N2");
2630 list.add("H21");
2631 list.add("H22");
2632 list.add("O6");
2633 }
2634 } else { // pyrimidine: either T to C or C to T
2635 if (res.equals("DTY") || res.equals("DT")) { // MTN is T
2636 list.add("H3");
2637 list.add("O4");
2638 list.add("C7");
2639 list.add("H71");
2640 list.add("H72");
2641 list.add("H73");
2642 } else { // MTN is C
2643 list.add("N4");
2644 list.add("H41");
2645 list.add("H42");
2646 list.add("H5");
2647 }
2648 }
2649 return list;
2650 }
2651
2652 /**
2653 * Determine if the mutation residue is purine.
2654 * @return true if mutation is purine
2655 */
2656 public boolean isMtnPurine() {
2657 return resName.equals("DA") || resName.equals("DG") || resName.equals("DAD") || resName.equals("DGU");
2658 }
2659
2660 /**
2661 * Determine if the mutation residue is pyrimidine.
2662 * @return true if mutation is pyrimidine
2663 */
2664 public boolean isMtnPyrimidine() {
2665 return resName.equals("DC") || resName.equals("DT") || resName.equals("DCY") || resName.equals("DTY");
2666 }
2667
2668 /**
2669 * Determine if original (wild type) residue is purine.
2670 * @return true if original residue is purine
2671 */
2672 public boolean isWtPurine() {
2673 return origResName.equals("DA") || origResName.equals("DG") || origResName.equals("DAD") || origResName.equals("DGU");
2674 }
2675
2676 /**
2677 * Determine if original (wild type) residue is pyrimidine.
2678 * @return true if original residue is pyrimidine
2679 */
2680 public boolean isWtPyrimidine() {
2681 return origResName.equals("DT") || origResName.equals("DC") || origResName.equals("DTY") || origResName.equals("DCY");
2682 }
2683 }
2684 }