Fix ModXNA compatibility

doc/ChangeLog.v7.md

@@ -85,3 +85,17 @@ Not Yet Ready
 - Introduce the `byatom` keyword to `checkchirality` action to check the chirality of any atoms in specified mask that are chiral centers.
 
 - Cannot yet build topology files for polarizable force fields (IPOL > 0).
+
+ModXNA Compatibility
+====================
+The old `sequence` command functionality and some new features have been added to ensure backwards compatibility with ModXNA (V1.9.2 as of right now).
+
+The following files/classes/functions have been retained from CPPTRAJ V6.31.0 with very few changes: `Structure/BuildAtom.h`, `Structure/Model.h`, `Structure/Model.cpp`
+
+The `OldBuilder` class in `Structure/OldBuilder.cpp` and `Structure/OldBuilder.h` is the `Builder` class from CPPTRAJ V6.31.0.
+
+The `Zmatrix::SetupICsAroundBond()` function has been retained from CPPTRAJ V6.31.0 `Zmatrix` class so that `sequence` can be made to work the way it used to.
+
+Default bond Req parameters are set when reading in mol2 files. This is needed so that the RK/REQ columns are printed for the `bonds` command, which is expected/needed for modxna.sh to determine the name of the H atom bonded to O3' (`--3cap` flag).
+
+ModXNA metadata (head/tail/anchor atom information) will be read from the title of Mol2 files if present and retained in the Topology. This metadata will be written back out to Mol2 files if another title has not been specified. This allows the `sequence` command to automatically recognized when units have ModXNA information and switch back to the old `sequence` mode automatically. 

src/Exec_Sequence.cpp

@@ -5,23 +5,21 @@
 #include "Parm/AssignParams.h"
 #include "Parm/GB_Params.h"
 #include "Structure/Builder.h"
+#include "Structure/OldBuilder.h" // TODO needed for modXNA
 #include "Structure/Creator.h"
 
-/** Generate and build the specified sequence. */
-int Exec_Sequence::generate_sequence(DataSet_Coords* OUT,
-                                     DataSetList const& DSL,
-                                     Sarray const& main_sequence,
-                                     Cpptraj::Structure::Creator const& creator)
-const
+/** CONSTRUCTOR */
+Exec_Sequence::Exec_Sequence() : Exec(COORDS), debug_(0), mode_(UNSPECIFIED) {}
+
+/** Get units in the sequence. */
+int Exec_Sequence::get_units(Uarray& Units,
+                             Sarray const& main_sequence,
+                             Cpptraj::Structure::Creator const& creator)
 {
-  std::string title;
   // First, get all units in order.
-  typedef std::vector<DataSet_Coords*> Uarray;
-  Uarray Units;
+  Units.clear();
   Units.reserve( main_sequence.size() );
 
-  int total_natom = 0;
-
   for (Sarray::const_iterator it = main_sequence.begin(); it != main_sequence.end(); ++it)
   {
     DataSet_Coords* unit = creator.IdTemplateFromName( *it );
@@ -37,10 +35,119 @@ const
     if (unit->Size() > 1) {
       mprintf("Warning: Unit '%s' has more than 1 frame. Only using first frame.\n", unit->legend());
     }
+    if (unit->Top().Modxna().HasModxna()) {
+      mprintf("\tUnit '%s' has ModXNA info.\n", unit->legend());
+      if (mode_ == UNSPECIFIED) {
+        mprintf("Warning: No sequence mode specified and ModXNA info detected; using old sequence mode.\n");
+        mode_ = OLD;
+      }
+    }
     Units.push_back( unit );
-    // Use the first unit name as title to match leap behavior
-    if (title.empty())
-      title = *it;
+  } // END loop over sequence
+  mprintf("\tFound %zu units.\n", Units.size());
+  if (Units.empty()) {
+    mprinterr("Error: No units.\n");
+    return 1;
+  }
+  return 0;
+}
+
+/** Old routine for building the specified sequence. TODO needed for modXNA */
+int Exec_Sequence::old_generate_sequence(DataSet_Coords* OUT, Uarray const& Units)
+const
+{
+  mprintf("\tUsing the old 'sequence' routine compatible with modXNA.\n");
+  // First, get all units in order.
+  typedef std::vector<int> Iarray;
+  Iarray connectAt0, connectAt1;
+  connectAt0.reserve( Units.size() );
+  connectAt1.reserve( Units.size() );
+  int total_natom = 0;
+
+  for (Uarray::const_iterator it = Units.begin(); it != Units.end(); ++it)
+  {
+    DataSet_Coords* unit = *it;
+
+    // Needs to have connect associated data
+    AssociatedData* ad = unit->GetAssociatedData(AssociatedData::CONNECT);
+    if (ad == 0) {
+      mprinterr("Error: Unit '%s' does not have CONNECT data.\n", unit->legend());
+      return 1;
+    }
+    AssociatedData_Connect const& C = static_cast<AssociatedData_Connect const&>( *ad );
+    if (C.NconnectAtoms() < 2) {
+      mprinterr("Error: Not enough connect atoms in unit '%s'\n", unit->legend());
+      return 1;
+    }
+    // Update connect atom 1 indices based on their position in the sequence.
+    // Do not update connect atom 0 indices since they will not yet be connected.
+    connectAt0.push_back( C.Connect()[0] );
+    connectAt1.push_back( C.Connect()[1] + total_natom );
+    total_natom += unit->Top().Natom();
+  } // END loop over sequence
+  mprintf("\tFound %zu units.\n", Units.size());
+  if (Units.empty()) {
+    mprinterr("Error: No units.\n");
+    return 1;
+  }
+  for (unsigned int idx = 0; idx < Units.size(); idx++)
+    mprintf("\tUnit %s HEAD %i TAIL %i\n", Units[idx]->legend(), connectAt0[idx]+1, connectAt1[idx]+1);
+
+  Topology combinedTop;
+  combinedTop.SetDebug( debug_ );
+  combinedTop.SetParmName( OUT->Meta().Name(), FileName() );
+  combinedTop.AppendTop( Units.front()->Top() );
+  //combinedTop.SetParmBox( Units // TODO
+  combinedTop.Brief("Sequence topology:");
+
+  Frame CombinedFrame = Units.front()->AllocateFrame();
+  Units.front()->GetFrame(0, CombinedFrame);
+
+
+  using namespace Cpptraj::Structure;
+  OldBuilder builder;
+  for (unsigned int idx = 1; idx < Units.size(); idx++) {
+    mprintf("\tConnect %s atom %i to %s atom %i\n",
+            Units[idx-1]->legend(), connectAt1[idx-1]+1,
+            Units[idx]->legend(),   connectAt0[idx]  +1);
+    Frame mol1frm = Units[idx]->AllocateFrame();
+    Units[idx]->GetFrame(0, mol1frm);
+    int bondat0 = connectAt1[idx-1];
+    int bondat1 = connectAt0[idx];
+    if (bondat0 < 0 || bondat1 < 0) {
+      mprinterr("Error: Invalid connect atom(s) between %s atom %i to %s atom %i\n",
+                Units[idx-1]->legend(), bondat0+1, Units[idx]->legend(), bondat1+1);
+      return 1;
+    }
+    if (builder.Combine( combinedTop, CombinedFrame, Units[idx]->Top(), mol1frm,
+                         connectAt1[idx-1], connectAt0[idx] )) {
+      mprinterr("Error: Sequence combine between units %u %s and %u %s failed.\n",
+                idx, Units[idx-1]->legend(), idx+1, Units[idx]->legend());
+      return 1;
+    }
+  }
+
+  // FIXME TODO common setup?
+
+  OUT->CoordsSetup(combinedTop, CombinedFrame.CoordsInfo());
+  OUT->AddFrame( CombinedFrame );
+
+  return 0;
+}
+
+
+/** Generate and build the specified sequence. */
+int Exec_Sequence::generate_sequence(DataSet_Coords* OUT,
+                                     std::string const& title,
+                                     Uarray const& Units,
+                                     Cpptraj::Structure::Creator const& creator)
+const
+{
+  // First, get total atom count. 
+  int total_natom = 0;
+  for (Uarray::const_iterator it = Units.begin(); it != Units.end(); ++it)
+  {
+    DataSet_Coords* unit = *it;
     total_natom += unit->Top().Natom();
   } // END loop over sequence
   mprintf("\tFound %zu units.\n", Units.size());
@@ -247,6 +354,19 @@ Exec::RetType Exec_Sequence::Execute(CpptrajState& State, ArgList& argIn)
 
   // Args
   int verbose = argIn.getKeyInt("verbose", 0);
+  mode_ = UNSPECIFIED;
+  std::string modeArg = argIn.GetStringKey("mode");
+  if (!modeArg.empty()) {
+    if (modeArg == "old")
+      mode_ = OLD;
+    else if (modeArg == "new")
+      mode_ = NEW;
+    else {
+      mprinterr("Error: Unrecognized 'mode' argument: %s\n", modeArg.c_str());
+      return CpptrajState::ERR;
+    }
+  }
+  std::string title = argIn.GetStringKey("title");
   // GB radii set. Default to mbondi
   Cpptraj::Parm::GB_Params gbradii;
   if (gbradii.Init_GB_Radii(argIn, Cpptraj::Parm::UNKNOWN_GB)) return CpptrajState::ERR;
@@ -286,9 +406,29 @@ Exec::RetType Exec_Sequence::Execute(CpptrajState& State, ArgList& argIn)
     mprintf(" %s", it->c_str());
   mprintf("\n");
   mprintf("\tOutput set name : %s\n", OUT->legend());
+  // If no title, use first unit name to be consistent with leap
+  if (title.empty())
+    title.assign( main_sequence.front() );
+
+  // Get units
+  Uarray Units;
+  if ( get_units(Units, main_sequence, creator) ) {
+    mprinterr("Error: Could not get units.\n");
+    return CpptrajState::ERR;
+  }
 
   // Execute
-  if (generate_sequence(OUT, State.DSL(), main_sequence, creator)) {
+  int err = 0;
+  switch (mode_) {
+    case UNSPECIFIED :
+    case NEW:
+      err = generate_sequence(OUT, title, Units, creator);
+      break;
+    case OLD:
+      err = old_generate_sequence(OUT, Units);
+      break;
+  }
+  if (err != 0) {
     mprinterr("Error: Could not generate sequence.\n");
     return CpptrajState::ERR;
   }

src/Exec_Sequence.h

@@ -9,16 +9,22 @@ class Creator;
 /// Create a molecule from a sequence of units 
 class Exec_Sequence : public Exec {
   public:
-    Exec_Sequence() : Exec(COORDS), debug_(0) {}
+    Exec_Sequence();
     void Help() const;
     DispatchObject* Alloc() const { return (DispatchObject*)new Exec_Sequence(); }
     RetType Execute(CpptrajState&, ArgList&);
   private:
     typedef std::vector<std::string> Sarray;
+    typedef std::vector<DataSet_Coords*> Uarray;
+    enum ModeType { UNSPECIFIED = 0, NEW, OLD };
 
-    int generate_sequence(DataSet_Coords*, DataSetList const&,
-                          Sarray const&, Cpptraj::Structure::Creator const&) const;
+    int get_units(Uarray&, Sarray const&, Cpptraj::Structure::Creator const&);
+
+    int old_generate_sequence(DataSet_Coords*, Uarray const&) const;
+
+    int generate_sequence(DataSet_Coords*, std::string const& title, Uarray const&, Cpptraj::Structure::Creator const&) const;
 
     int debug_;
+    ModeType mode_;
 };
 #endif

src/ModXNA_Info.cpp

@@ -44,3 +44,24 @@ void ModXNA_Info::PrintModxna() const {
   printmxna("\tAnchor      ", anchor_);
   printmxna("\tAnchorStrip ", anchorStrip_);
 }
+
+/** Generate ModXNA metadata string */
+std::string ModXNA_Info::GenMetadataString() const {
+  if (fragmentName_.empty()) return std::string("");
+  std::string out("modXNA-fragment:" + fragmentName_);
+  // Base
+  out.append(":HEAD01:" + head_);
+  out.append(":HEAD01STRIP:" + headStrip_);
+  if (!tail_.empty()) {
+    // BB
+    out.append(":TAIL01:" + tail_);
+    out.append(":TAIL01STRIP:" + tailStrip_);
+    if (!anchor_.empty()) {
+      // Sugar
+      out.append(":HEAD03:" + anchor_);
+      out.append(":HEAD03STRIP:" + anchorStrip_);
+    }
+  }
+  return out;
+}
+

src/ModXNA_Info.h

@@ -17,6 +17,8 @@ class ModXNA_Info {
     bool HasModxna() const { return (!fragmentName_.empty()); }
     /// Print to stdout
     void PrintModxna() const;
+    /// Generate ModXNA metadata string for e.g. writing ModXNA mol2 files
+    std::string GenMetadataString() const;
 
     std::string const& FragmentName() const { return fragmentName_; }
     std::string const& Head()         const { return head_; }

src/Parm_Mol2.cpp

@@ -3,6 +3,7 @@
 #include "BondSearch.h"
 #include "CpptrajStdio.h"
 #include "ModXNA_Info.h"
+#include "Parm/Bond_Params.h"
 
 // Parm_Mol2::ID_ParmFormat() 
 bool Parm_Mol2::ID_ParmFormat(CpptrajFile& fileIn) {
@@ -51,6 +52,9 @@ int Parm_Mol2::ReadParm(FileName const& fname, Topology &parmOut) {
     bondSearch.FindBonds( parmOut, searchType_, Coords, Offset_, debug_ );
   }
 
+  // For modXNA backwards-compatibility, generate bond parameters TODO remove
+  Cpptraj::Parm::GenerateBondParameters(parmOut);
+
   // No box
   parmOut.SetParmBox( Box() );
 

src/Structure/BuildAtom.h

@@ -0,0 +1,29 @@
+#ifndef INC_STRUCTURE_BUILDATOM_H
+#define INC_STRUCTURE_BUILDATOM_H
+#include "StructureEnum.h"
+namespace Cpptraj {
+namespace Structure {
+/// Old (<v7) class for holding information for an atom used when building/modelling new coordinates. TODO for modXNA
+class BuildAtom {
+  public:
+    BuildAtom() : ctype_(IS_UNKNOWN_CHIRALITY) {}
+    /// Set chirality
+    void SetChirality(ChiralType c) { ctype_ = c; }
+    /// Set size of priority array based on number of bonds
+    //void SetNbonds(int n) { priority_.assign(n, -1); }
+    /// \return Pointer to priority array
+    //int* PriorityPtr() { return &(priority_[0]); }
+    /// Used to modify the priority array
+    std::vector<int>& ModifyPriority() { return priority_; }
+
+    /// \return Atom chirality
+    ChiralType Chirality() const { return ctype_; }
+    /// \return Priority array
+    std::vector<int> const& Priority() const { return priority_; }
+  private:
+    ChiralType ctype_;          ///< Chirality around atom
+    std::vector<int> priority_; ///< Indices of bonded atoms, sorted by priority
+};
+}
+}
+#endif

src/Structure/CMakeLists.txt

@@ -12,6 +12,8 @@ target_sources(cpptraj_common_obj PRIVATE
   ${CMAKE_CURRENT_LIST_DIR}/InternalCoords.cpp
   ${CMAKE_CURRENT_LIST_DIR}/LeastSquaresPlane.cpp
   ${CMAKE_CURRENT_LIST_DIR}/MetalCenterFinder.cpp
+  ${CMAKE_CURRENT_LIST_DIR}/Model.cpp
+  ${CMAKE_CURRENT_LIST_DIR}/OldBuilder.cpp
   ${CMAKE_CURRENT_LIST_DIR}/PdbCleaner.cpp
   ${CMAKE_CURRENT_LIST_DIR}/RingFinder.cpp
   ${CMAKE_CURRENT_LIST_DIR}/Solvate.cpp