Ensure all PDB-related Amber Topology sections are read and written

doc/cpptraj.lyx

@@ -12521,7 +12521,7 @@ parm
 \end_layout
 
 \begin_layout LyX-Code
-      }]
+      }] [nomolsearch] [renumresidues]
 \end_layout
 
 \begin_deeper
@@ -12560,16 +12560,6 @@ reference "subsec:Parameter/Reference-Tagging"
 
 \begin_layout Description
 
-\series bold
-[nobondsearch]
-\series default
- Optional; if specified do not search for bonds via geometry if Topology
- does not include bond information.
- May cause some Actions to fail.
-\end_layout
-
-\begin_layout Description
-
 \series bold
 [bondsearch <offset>]
 \series default
@@ -12601,6 +12591,34 @@ pairlist Uses a pair list to search for bonds between atoms.
 \end_layout
 
 \end_deeper
+\begin_layout Standard
+
+\series bold
+Advanced Options - Not recommended for general use
+\end_layout
+
+\begin_layout Description
+
+\series bold
+[nobondsearch]
+\series default
+ If specified do not search for bonds via geometry if Topology does not
+ include bond information.
+ May cause some Actions to fail.
+\end_layout
+
+\begin_layout Description
+[nomolsearch] If specified do not search for molecule information.
+ May cause some Actions to fail.
+\end_layout
+
+\begin_layout Description
+[renumresidues] If specified, ensure that any residue cannot be part of
+ more than 1 molecule (can occur with e.g.
+ alternate sites).
+ Residues will be renumbered according to molecule information in that case.
+\end_layout
+
 \end_deeper
 \begin_layout Standard
 Read in parameter file.

src/BondSearch.cpp

@@ -44,7 +44,7 @@ Box CreateBoundingBox(Frame const& frameIn) {
 /** Add bonds within residues to top using coords in frameIn. */
 void BondsWithinResidues(Topology& top, Frame const& frameIn, double offset) {
   // ----- STEP 1: Determine bonds within residues
-  std::vector<AtomExtra> const& Extra = top.Extra();
+  std::vector<char> const& AtomAltLoc = top.AtomAltLoc();
   for (Topology::res_iterator res = top.ResStart(); res != top.ResEnd(); ++res)
   {
     int stopatom = res->LastAtom();
@@ -55,15 +55,15 @@ void BondsWithinResidues(Topology& top, Frame const& frameIn, double offset) {
       if (a1Elt==Atom::HYDROGEN && top[atom1].Nbonds() > 0 )
         continue;
       // Determine if atom1 has an alternate location
-      bool hasAltLoc1 = (!Extra.empty() && Extra[atom1].AtomAltLoc() != ' ');
+      bool hasAltLoc1 = (!AtomAltLoc.empty() && AtomAltLoc[atom1] != ' ');
       // Loop over all other atoms in the residue
       for (int atom2 = atom1 + 1; atom2 != stopatom; ++atom2) {
         // Determine if atom2 has an alternate location. Only matters if
         // atom1 has an alternate location.
-        if (hasAltLoc1 && Extra[atom2].AtomAltLoc() != ' ') {
+        if (hasAltLoc1 && AtomAltLoc[atom2] != ' ') {
           // If alternate location info present for both atoms, make sure the
           // location IDs match.
-          if (Extra[atom1].AtomAltLoc() != Extra[atom2].AtomAltLoc()) continue;
+          if (AtomAltLoc[atom1] != AtomAltLoc[atom2]) continue;
         }
         Atom::AtomicElementType a2Elt = top[atom2].Element();
         double D2 = DIST2_NoImage(frameIn.XYZ(atom1), frameIn.XYZ(atom2) );

src/Exec_AddMissingRes.cpp

@@ -765,7 +765,7 @@ const
   mprintf("DEBUG: LJ radius= %g\n", CAtop.GetVDWradius(0));
   // Final setup
   CAtop.SetParmName("caseq", "seq.ca.mol2");
-  CAtop.CommonSetup(true, 2); // molecule search, exclude bonds
+  CAtop.CommonSetup(true, 2, false); // molecule search, exclude bonds, do not renumber residues
   CAtop.Summary();
 
   // Write CA top
@@ -964,7 +964,7 @@ const
   // Finish new top
   //newTop.SetParmName("seqpdb", "seq.pdb");
   newTop.SetParmName("seqmol2", "seq.mol2");
-  newTop.CommonSetup( false ); // No molecule search
+  newTop.CommonSetup( false, false ); // No molecule search, no renumbering residues
   newTop.Summary();
   // DEBUG print new topology
   mprintf("New topology:\n");
@@ -1213,7 +1213,7 @@ const
                     originalAtToNew[bnd->A2()] );
   // Finish new top and write
   newTop.SetParmName("newpdb", "temp.pdb");
-  newTop.CommonSetup( false ); // No molecule search
+  newTop.CommonSetup( false, false ); // No molecule search, no renumbering residues
   newTop.Summary();
   if (WriteStructure("temp.pdb", &newTop, newFrame, TrajectoryFile::PDBFILE)) {
     mprinterr("Error: Write of temp.pdb failed.\n");
@@ -1260,7 +1260,7 @@ const
   mprintf("DEBUG: LJ radius= %g\n", CAtop.GetVDWradius(0));
   // Final setup
   CAtop.SetParmName("capdb", "temp.ca.mol2");
-  CAtop.CommonSetup(true, 2); // molecule search, exclude bonds
+  CAtop.CommonSetup(true, 2, false); // molecule search, exclude bonds, no renumbering residues
   CAtop.Summary();
 
   // Write CA top

src/Exec_Top.cpp

@@ -1,20 +1,12 @@
 #include "Exec_Top.h"
 #include "CpptrajStdio.h"
 #include "TopInfo.h"
+#include "ParmFile.h"
 
 void Exec_LoadParm::Help() const {
-  mprintf("\t<filename> [{[TAG] | name <setname>}]\n"
-          "\t [{ nobondsearch |\n"
-          "\t    [bondsearch <offset>] [searchtype {grid|pairlist}]\n"
-          "\t  }]\n"
-          "  Add <filename> to the topology list.\n"
-          "  For topologies that may not have bond information, 'bondsearch <offset>'\n"
-          "  controls the offset that will be added to atom-atom distances when\n"
-          "  searching for bonds (default 0.2 Ang), and 'searchtype' specifies\n"
-          "  different algorithms that can be used for searching bonds (still\n"
-          "  experimental. Bond searching can be skipped via 'nobondsearch' (not\n"
-          "  recommended).\n"
-          "  Use 'help Formats parm' for format-specific options.\n");
+  mprintf("\t<filename> [{[TAG] | name <setname>}]\n%s", ParmFile::ReadTopologyKeywords());
+  mprintf("  Add <filename> to the topology list.\n%s", ParmFile::ReadTopologyHelp());
+  mprintf("  Use 'help Formats parm' for format-specific options.\n");
 }
 // -----------------------------------------------------------------------------
 void Exec_ParmInfo::Help() const {

src/ParmFile.cpp

@@ -78,6 +78,29 @@ int ParmFile::ReadTopology(Topology& t, FileName const& n, int d) {
   return ReadTopology(t, n, ArgList(), d);
 }
 
+/** Keywords for ReadTopology() */
+const char* ParmFile::ReadTopologyKeywords() {
+  return
+  "\t [{ nobondsearch |\n"
+  "\t    [bondsearch <offset>] [searchtype {grid|pairlist}]\n"
+  "\t  }] [nomolsearch] [renumresidues]\n";
+}
+
+/** More detailed help for ReadTopology() */
+const char* ParmFile::ReadTopologyHelp() {
+  return
+  "  For topologies that may not have bond information, 'bondsearch <offset>'\n"
+  "  controls the offset that will be added to atom-atom distances when\n"
+  "  searching for bonds (default 0.2 Ang), and 'searchtype' specifies\n"
+  "  alternative (and still experimental) algorithms that can be used in\n"
+  "  place of the usual bond search algorithm.\n"
+  "  ** ADVANCED OPTIONS - NOT RECOMMENDED FOR GENERAL USE. **\n"
+  "  The 'nobondsearch' keyword can be specified to skip searching for bonds.\n"
+  "  The 'nomolsearch' keyword can be specified to skip molecule determintation via bonds.\n"
+  "  The 'renumresidues' keyword can be specified to ensure that any residue cannot\n"
+  "  be part of more than 1 molecule (can occur with e.g. alternate sites).\n";
+}
+
 // ParmFile::ReadTopology()
 int ParmFile::ReadTopology(Topology& Top, FileName const& fnameIn, 
                            ArgList const& argListIn, int debugIn) 
@@ -116,8 +139,12 @@ int ParmFile::ReadTopology(Topology& Top, FileName const& fnameIn,
   }
   double bondoffset = argIn.getKeyDouble("bondsearch", -1.0);
   bool molsearch = !argIn.hasKey("nomolsearch");
+  bool renumberResidues = argIn.hasKey("renumresidues");
   if (!molsearch)
     mprintf("\tDisabling molecule search. Topology will have no molecule info.\n");
+  if (renumberResidues)
+    mprintf("\tIf any residue corresponds to more than 1 molecule, residues will be renumbered\n"
+            "\t  according to molecule information.\n");
   // Only force bond search when 'bondsearch' is specified.
 //  bool bondsearch = false;
 //  if (argIn.Contains("bondsearch")) {
@@ -148,7 +175,7 @@ int ParmFile::ReadTopology(Topology& Top, FileName const& fnameIn,
   int err = parmio->ReadParm( parmName_.Full(), Top);
   // Perform setup common to all parm files.
   if (err == 0) 
-    err = Top.CommonSetup( molsearch );
+    err = Top.CommonSetup( molsearch, renumberResidues );
   else
     mprinterr("Error reading topology file '%s'\n", parmName_.full());
   delete parmio;

src/ParmFile.h

@@ -23,6 +23,10 @@ class ParmFile {
     static void WriteOptions(std::string const& fkey){ FileTypes::Options(PF_WriteKeyArray,PF_AllocArray,UNKNOWN_PARM,fkey,FileTypes::WRITEOPT);}
     /// CONSTRUCTOR
     ParmFile();
+    /// ReadTopology() keywords
+    static const char* ReadTopologyKeywords();
+    /// ReadTopology() help
+    static const char* ReadTopologyHelp();
     /// Read topology file with optional arguments and debug level
     int ReadTopology(Topology&, FileName const&, ArgList const&,int);
     /// Read topology file