Add ability to filter out alternate atom locations from PDB

doc/cpptraj.lyx

@@ -12399,10 +12399,17 @@ reference "tab:cpptraj-Topology-formats"
 
 \emph on
 PDB format:
+\begin_inset CommandInset label
+LatexCommand label
+name "subsec:cpptraj-parm-pdb"
+
+\end_inset
+
+
 \end_layout
 
 \begin_layout LyX-Code
-[pqr] [readbox] [conect] [noconect] [link] [nolink]
+[pqr] [readbox] [conect] [noconect] [link] [nolink] [keepaltloc <char>]
 \begin_inset Separator latexpar
 \end_inset
 
@@ -12436,6 +12443,19 @@ PDB format:
 [nolink] Do not read LINK records if present (default).
 \end_layout
 
+\begin_layout Description
+[keepaltloc
+\begin_inset space ~
+\end_inset
+
+<char>] If specified, only keep alternate atom location IDs matching the
+ specified character 
+\series bold
+<char>
+\series default
+.
+\end_layout
+
 \end_deeper
 \begin_layout Paragraph*
 
@@ -12469,10 +12489,34 @@ If this is the case
 \shape italic
 cpptraj
 \shape default
- will print a warning about duplicate atom names but will take no other
- action.
+ will print a warning about alternate location IDs being present but will
+ take no other action.
  Both residues are considered 'CYS' and the mask ':CYS@CA' would select
- both atom 806 and 809.
+ both atom 806 and 808.
+ If desired, a specific location ID can be kept via the 
+\series bold
+keepaltloc
+\series default
+ keyword.
+ If 
+\series bold
+keepaltloc
+\series default
+ is specified, it should also be specified for any 
+\series bold
+trajin
+\series default
+ commands (see 
+\begin_inset CommandInset ref
+LatexCommand vref
+reference "subsec:cpptraj-trajin-pdb"
+plural "false"
+caps "false"
+noprefix "false"
+
+\end_inset
+
+).
  Residue insertion codes are read in but also not used by the mask parser.
 \end_layout
 
@@ -15002,6 +15046,64 @@ shape Force reading of box information as shape matrix.
 \end_layout
 
 \end_deeper
+\begin_layout Subsubsection
+Options for PDB files:
+\begin_inset CommandInset label
+LatexCommand label
+name "subsec:cpptraj-trajin-pdb"
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout LyX-Code
+[keepaltloc <char>]
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+[keepaltloc
+\begin_inset space ~
+\end_inset
+
+<char>] If specified, only keep alternate atom location IDs matching the
+ specified character 
+\series bold
+<char>
+\series default
+.
+\end_layout
+
+\end_deeper
+\begin_layout Standard
+Note that if 
+\series bold
+keepaltloc
+\series default
+ is specified, the associated topology should not have alternate location
+ IDs, i.e.
+ if the topology is from a PDB the 
+\series bold
+keepaltloc
+\series default
+ keyword may need to be used with the 
+\series bold
+parm
+\series default
+ command (see 
+\begin_inset CommandInset ref
+LatexCommand vref
+reference "subsec:cpptraj-parm-pdb"
+plural "false"
+caps "false"
+noprefix "false"
+
+\end_inset
+
+).
+\end_layout
+
 \begin_layout Subsection
 trajout
 \begin_inset CommandInset label

src/PDBfile.cpp

@@ -127,6 +127,12 @@ PDBfile::PDB_RECTYPE PDBfile::NextRecord() {
   return recType_;
 }
 
+/** \return Atom alt. loc. code from PDB ATOM/HETATM line. */
+char PDBfile::pdb_AltLoc() const {
+  return linebuffer_[16];
+}
+
+/** \return Atom containing information from PDB ATOM/HETATM line. */
 Atom PDBfile::pdb_Atom(char& altLoc, int& atnum) {
   // ATOM or HETATM keyword.
   // Check line length before any modification.
@@ -214,7 +220,8 @@ void PDBfile::pdb_ChargeAndRadius(float& charge, float& radius) {
   sscanf(linebuffer_+54, "%f %f", &charge, &radius);
 }
 
-void PDBfile::pdb_Box(double* box) {
+/** Set box[0-5] with A B C ALPHA BETA GAMMA from CRYST1 record. */
+void PDBfile::readCRYST1(double* box) {
   // CRYST1 keyword. RECORD A B C ALPHA BETA GAMMA SGROUP Z
   unsigned int lb_size = strlen(linebuffer_);
   if (lb_size < 54) {
@@ -238,14 +245,30 @@ void PDBfile::pdb_Box(double* box) {
     box[ib] = atof( linebuffer_ + lb );
     linebuffer_[end] = savechar;
   }
-  mprintf("\tRead CRYST1 info from PDB: a=%g b=%g c=%g alpha=%g beta=%g gamma=%g\n",
-          box[0], box[1], box[2], box[3], box[4], box[5]);
-  // Warn if the box looks strange.
+}
+
+/** Print a warning to STDOUT if unit cell lengths are strange. */
+static inline void box_warning(const double* box) {
   if (box[0] == 1.0 && box[1] == 1.0 && box[2] == 1.0)
     mprintf("Warning: PDB cell lengths are all 1.0 Ang.;"
             " this usually indicates an invalid box.\n");
 }
 
+/** Read box info from CRYST1, verbose. */
+void PDBfile::pdb_Box_verbose(double* box) {
+  readCRYST1(box);
+  box_warning(box);
+  mprintf("\tRead CRYST1 info from PDB: a=%g b=%g c=%g alpha=%g beta=%g gamma=%g\n",
+          box[0], box[1], box[2], box[3], box[4], box[5]);
+}
+
+/** Read box info from CRYST1, warn only if box is strange looking. */
+void PDBfile::pdb_Box_terse(double* box) {
+  readCRYST1(box);
+  box_warning(box);
+}
+
+/** Read serial #s of atoms from a CONECT record. */
 int PDBfile::pdb_Bonds(int* bnd) {
   unsigned int lb_size = strlen(linebuffer_);
   int Nscan = 0;

src/PDBfile.h

@@ -17,6 +17,8 @@ class PDBfile : public CpptrajFile {
     static bool ID_PDB(CpptrajFile&);
     /// \return the type of the next PDB record read.
     PDB_RECTYPE NextRecord();
+    /// \return ATOM/HETATM alt. loc ID
+    char pdb_AltLoc() const;
     /// \return Atom info with name and element for ATOM/HETATM; set altLoc and #.
     Atom pdb_Atom(char&, int&);
     /// \return Atom info with name and element for ATOM/HETATM.
@@ -29,8 +31,10 @@ class PDBfile : public CpptrajFile {
     void pdb_OccupancyAndBfactor(float&, float&);
     /// Get charge and radius from PQR ATOM/HETATM record.
     void pdb_ChargeAndRadius(float&, float&);
-    /// Set given XYZ array with A/B/C/alpha/beta/gamma from CRYST1 record.
-    void pdb_Box(double*);
+    /// Set given XYZ array with A/B/C/alpha/beta/gamma from CRYST1 record, verbose.
+    void pdb_Box_verbose(double*);
+    /// Set given XYZ array with A/B/C/alpha/beta/gamma from CRYST1 record, terse.
+    void pdb_Box_terse(double*);
     /// Set given array with atom and #s of bonded atoms from CONECT record.
     int pdb_Bonds(int*);
     /// \return Link record.
@@ -84,6 +88,8 @@ class PDBfile : public CpptrajFile {
   private:
     /// \return true if the first 6 chars of buffer match a PDB keyword
     static bool IsPDBkeyword(std::string const&);
+    /// Read box info from CRYST1 record
+    void readCRYST1(double*);
 
     int anum_;               ///< Atom number for writing.
     PDB_RECTYPE recType_;    ///< Current record type.

src/Parm_PDB.cpp

@@ -6,15 +6,27 @@
 # include "Timer.h"
 #endif
 
+/** CONSTRUCTOR */
+Parm_PDB::Parm_PDB() :
+  ConectMode_(UNSPECIFIED),
+  LinkMode_(UNSPECIFIED),
+  keepAltLoc_(' '),
+  readAsPQR_(false),
+  readBox_(false) {}
+
+// Parm_PDB::ReadHelp()
 void Parm_PDB::ReadHelp() {
-  mprintf("\tpqr      : Read atomic charge/radius from occupancy/B-factor columns (PQR).\n"
-          "\treadbox  : Read unit cell information from CRYST1 record if present.\n"
-          "\tconect   : Read CONECT records if present (default).\n"
-          "\tnoconect : Do not read CONECT records if present.\n"
-          "\tlink     : Read LINK records if present.\n"
-          "\tnolink   : Do not read LINK records if present (default).\n");
+  mprintf("\tpqr               : Read atomic charge/radius from occupancy/B-factor columns (PQR).\n"
+          "\treadbox           : Read unit cell information from CRYST1 record if present.\n"
+          "\tconect            : Read CONECT records if present (default).\n"
+          "\tnoconect          : Do not read CONECT records if present.\n"
+          "\tlink              : Read LINK records if present.\n"
+          "\tnolink            : Do not read LINK records if present (default).\n"
+          "\tkeepaltloc <char> : If specified, alternate location ID to keep.\n"
+         );
 }
 
+// Parm_PDB::processReadArgs()
 int Parm_PDB::processReadArgs(ArgList& argIn) {
   readAsPQR_ = argIn.hasKey("pqr");
   readBox_ = argIn.hasKey("readbox");
@@ -26,9 +38,13 @@ int Parm_PDB::processReadArgs(ArgList& argIn) {
     LinkMode_ = READ;
   else if (argIn.hasKey("nolink"))
     LinkMode_ = SKIP;
+  std::string keepAltChar = argIn.GetStringKey("keepaltloc");
+  if (!keepAltChar.empty())
+    keepAltLoc_ = keepAltChar[0];
   return 0;
 }
 
+// Parm_PDB::ReadParm()
 int Parm_PDB::ReadParm(FileName const& fname, Topology &TopIn) {
   typedef std::vector<PDBfile::Link> Larray;
   PDBfile infile;
@@ -64,16 +80,19 @@ int Parm_PDB::ReadParm(FileName const& fname, Topology &TopIn) {
     mprintf("\tReading bond info from LINK records.\n");
   else
     mprintf("\tNot reading bond info from LINK records.\n");
+  if (keepAltLoc_ != ' ')
+    mprintf("\tWhen present, only reading alternate location ID %c\n", keepAltLoc_);
 # ifdef TIMER
   Timer time_total, time_atom;
   time_total.Start();
 # endif
   bool missingResidues = false;
+  int nAltLocSkipped = 0;
   // Loop over PDB records
   while ( infile.NextRecord() != PDBfile::END_OF_FILE ) {
     if (readBox_ && infile.RecType() == PDBfile::CRYST1) {
       // Box info from CRYST1 record.
-      infile.pdb_Box( XYZ );
+      infile.pdb_Box_verbose( XYZ );
       TopIn.SetParmBox( XYZ );
     } else if (infile.RecType() == PDBfile::CONECT && readConect) {
       // BOND - first element will be atom, next few are bonded atoms.
@@ -101,6 +120,11 @@ int Parm_PDB::ReadParm(FileName const& fname, Topology &TopIn) {
       // If this is an ATOM / HETATM keyword, add to topology.
       infile.pdb_XYZ( XYZ );
       Atom pdbAtom = infile.pdb_Atom(altLoc, atnum);
+      // Check if we are filtering alt loc IDs
+      if (keepAltLoc_ != ' ' && altLoc != ' ' && altLoc != keepAltLoc_) {
+        nAltLocSkipped++;
+        continue;
+      }
       if (atnum >= (int)serial.size())
         serial.resize( atnum+1, -1 );
       serial[atnum] = TopIn.Natom();
@@ -113,6 +137,12 @@ int Parm_PDB::ReadParm(FileName const& fname, Topology &TopIn) {
         TopIn.AddExtraAtomInfo( AtomExtra(occupancy, bfactor, altLoc) );
       }
       TopIn.AddTopAtom(pdbAtom, infile.pdb_Residue());
+      if (altLoc != ' ' && keepAltLoc_ == ' ') {
+        Residue const& lastRes = TopIn.Res(TopIn.Nres()-1);
+        mprintf("Warning: Atom %i %s in res %s %i %c has alternate location specifier %c\n",
+                atnum, *pdbAtom.Name(), *lastRes.Name(), lastRes.OriginalResNum(),
+                lastRes.ChainId(), altLoc);
+      }
       Coords.AddXYZ( XYZ );
 #     ifdef TIMER
       time_atom.Stop();
@@ -131,7 +161,9 @@ int Parm_PDB::ReadParm(FileName const& fname, Topology &TopIn) {
       if (readLink)
         mprintf("Warning: If molecule determination fails try not specifying 'link' instead.\n");
     }
-  }
+  } // END loop over PDB records
+  if (nAltLocSkipped > 0)
+    mprintf("\tSkipped %i alternate atom locations.\n", nAltLocSkipped);
   // Sanity check
   if (TopIn.Natom() < 1) {
     mprinterr("Error: No atoms present in PDB.\n");

src/Parm_PDB.h

@@ -3,7 +3,7 @@
 #include "ParmIO.h"
 class Parm_PDB : public ParmIO {
   public :
-    Parm_PDB() : ConectMode_(UNSPECIFIED), LinkMode_(UNSPECIFIED), readAsPQR_(false), readBox_(false) {}
+    Parm_PDB();
     static BaseIOtype* Alloc() { return (BaseIOtype*)new Parm_PDB(); }
     static void ReadHelp();
     bool ID_ParmFormat(CpptrajFile&);
@@ -15,6 +15,7 @@ class Parm_PDB : public ParmIO {
     enum ReadType { UNSPECIFIED = 0, READ, SKIP };
     ReadType ConectMode_; ///< Specify how to handle CONECT records.
     ReadType LinkMode_;   ///< Specify how to handle LINK records.
+    char keepAltLoc_;     ///< Alternate location to keep
     bool readAsPQR_;      ///< If true get charge and radius from occ/b factor cols
     bool readBox_;        ///< If true try to read CRYST1 record as box info.
 };

src/ReferenceAction.cpp

@@ -131,6 +131,7 @@ int ReferenceAction::InitRef(ArgList& argIn, DataSetList const& DSLin,
 
 // ReferenceAction::SetupRefMask()
 int ReferenceAction::SetupRefMask(Topology const& topIn) {
+  mprintf("\tReference topology: %s\n", topIn.c_str());
   if (refMask_.MaskStringSet()) {
     if (topIn.SetupIntegerMask( refMask_ )) return 1;
     mprintf("\tReference mask:");