diff --git a/doc/ChangeLog.v7.md b/doc/ChangeLog.v7.md
index 2665ef855b..21382198e7 100644
--- a/doc/ChangeLog.v7.md
+++ b/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. 
diff --git a/src/Exec_Sequence.cpp b/src/Exec_Sequence.cpp
index 775dbe0819..6d6683f8da 100644
--- a/src/Exec_Sequence.cpp
+++ b/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;
   }
diff --git a/src/Exec_Sequence.h b/src/Exec_Sequence.h
index 12c71db23b..e2e45d8ae5 100644
--- a/src/Exec_Sequence.h
+++ b/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
diff --git a/src/ModXNA_Info.cpp b/src/ModXNA_Info.cpp
index d08b065f29..526ad2a8c5 100644
--- a/src/ModXNA_Info.cpp
+++ b/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;
+}
+    
diff --git a/src/ModXNA_Info.h b/src/ModXNA_Info.h
index 3b7714163f..73f23eba48 100644
--- a/src/ModXNA_Info.h
+++ b/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_; }
diff --git a/src/Parm_Mol2.cpp b/src/Parm_Mol2.cpp
index c32b105a16..b5e90e4403 100644
--- a/src/Parm_Mol2.cpp
+++ b/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() );
 
diff --git a/src/Structure/BuildAtom.h b/src/Structure/BuildAtom.h
new file mode 100644
index 0000000000..3fa2b56a0a
--- /dev/null
+++ b/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
diff --git a/src/Structure/CMakeLists.txt b/src/Structure/CMakeLists.txt
index a98547d991..68d133db81 100644
--- a/src/Structure/CMakeLists.txt
+++ b/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
diff --git a/src/Structure/Model.cpp b/src/Structure/Model.cpp
new file mode 100644
index 0000000000..affbf21f06
--- /dev/null
+++ b/src/Structure/Model.cpp
@@ -0,0 +1,258 @@
+#include "Model.h"
+#include "BuildAtom.h"
+#include "../CpptrajStdio.h"
+#include "../Topology.h"
+#include "../TorsionRoutines.h"
+#include "../Constants.h"
+
+/** Attempt to assign a reasonable value for theta internal coordinate for
+  * atom i given that atoms j and k have known positions.
+  */
+int Cpptraj::Structure::Model::AssignTheta(double& theta, int ai, int aj, int ak, Topology const& topIn, Frame const& frameIn, std::vector<bool> const& atomPositionKnown, int debug)
+{
+  // Figure out hybridization and chirality of atom j.
+  if (debug > 0)
+    mprintf("DEBUG: AssignTheta for atom j : %s\n", topIn.AtomMaskName(aj).c_str());
+
+  enum HybridizationType { SP = 0, SP2, SP3, UNKNOWN_HYBRIDIZATION };
+
+  Atom const& AJ = topIn[aj];
+  if (debug > 0) {
+    mprintf("DEBUG:\t\tI %s Nbonds: %i\n", topIn[ai].ElementName(), topIn[ai].Nbonds());
+    mprintf("DEBUG:\t\tJ %s Nbonds: %i\n", AJ.ElementName(), AJ.Nbonds());
+    mprintf("DEBUG:\t\tK %s Nbonds: %i\n", topIn[ak].ElementName(), topIn[ak].Nbonds());
+  }
+  // Sanity check
+  if (AJ.Nbonds() < 2) {
+    mprinterr("Internal Error: AssignTheta() called for atom J %s with fewer than 2 bonds.\n", topIn.AtomMaskName(aj).c_str());
+    return 1;
+  }
+
+  HybridizationType hybrid = UNKNOWN_HYBRIDIZATION;
+  // Handle specific elements
+  switch (AJ.Element()) {
+    case Atom::CARBON :
+      switch (AJ.Nbonds()) {
+        case 2 : hybrid = SP; break;
+        case 3 : hybrid = SP2; break;
+        case 4 : hybrid = SP3; break;
+      }
+      break;
+    case Atom::NITROGEN :
+      switch (AJ.Nbonds()) {
+        case 2 : hybrid = SP2; break;
+        case 3 :
+          // Check for potential SP2. If only 1 of the bonded atoms is
+          // hydrogen, assume SP2. TODO actually check for aromaticity.
+          int n_hydrogens = 0;
+          for (Atom::bond_iterator bat = AJ.bondbegin(); bat != AJ.bondend(); ++bat)
+            if (topIn[*bat].Element() == Atom::HYDROGEN)
+              n_hydrogens++;
+          if (n_hydrogens == 1)
+            hybrid = SP2;
+          else
+            hybrid = SP3;
+          break;
+      }
+      break;
+    case Atom::OXYGEN :
+      switch (AJ.Nbonds()) {
+        case 2 : hybrid = SP3; break;
+      }
+      break;
+    case Atom::SULFUR :
+      switch (AJ.Nbonds()) {
+        case 2 : hybrid = SP3; break;
+      }
+      break;
+    default: hybrid = UNKNOWN_HYBRIDIZATION; break;
+  }
+  // Fill in what values we can for known atoms
+/*  std::vector<double> knownTheta( AJ.Nbonds() );
+  int knownIdx = -1;
+  for (int idx = 0; idx < AJ.Nbonds(); idx++) {
+    int atnum = AJ.Bond(idx);
+    if (atnum != ak && atomPositionKnown[atnum]) {
+      knownTheta[idx] = CalcAngle(frameIn.XYZ(atnum),
+                                  frameIn.XYZ(aj),
+                                  frameIn.XYZ(ak));
+      mprintf("DEBUG:\t\tKnown theta for %s = %g\n", topIn.AtomMaskName(atnum).c_str(), knownTheta[idx]*Constants::RADDEG);
+      if (knownIdx == -1) knownIdx = idx; // FIXME handle more than 1 known
+    }
+  }
+  if (knownIdx == -1) {*/
+    //mprintf("DEBUG:\t\tNo known theta.\n");
+  if (hybrid == UNKNOWN_HYBRIDIZATION) {
+    // Assign a theta based on number of bonds 
+    switch (AJ.Nbonds()) {
+      case 4 : hybrid = SP3; break;
+      case 3 : hybrid = SP2; break;
+      case 2 : hybrid = SP; break;
+      default : mprinterr("Internal Error: AssignTheta(): Unhandled # bonds for %s (%i)\n", topIn.AtomMaskName(aj).c_str(), AJ.Nbonds()); return 1;
+    }
+  }
+  // Assign a theta based on hybridization
+  switch (hybrid) {
+    case SP3 : theta = 109.5 * Constants::DEGRAD; break;
+    case SP2 : theta = 120.0 * Constants::DEGRAD; break;
+    case SP  : theta = 180.0 * Constants::DEGRAD; break;
+    default : mprinterr("Internal Error: AssignTheta(): Unhandled hybridization for %s (%i)\n", topIn.AtomMaskName(aj).c_str(), AJ.Nbonds()); return 1;
+  }/*
+  } else {
+    theta = knownTheta[knownIdx]; // TODO just use above guess via hybrid?
+  }*/
+
+  return 0;
+}
+
+/// Recursive function to return depth from an atom along bonds
+static int atom_depth(int& depth,
+                      int at, Topology const& topIn, std::vector<bool>& visited, int maxdepth)
+{
+  if (depth == maxdepth) return 0;
+  depth++;
+  visited[at] = true;
+  int depthFromHere = 1;
+  for (Atom::bond_iterator bat = topIn[at].bondbegin(); bat != topIn[at].bondend(); ++bat)
+  {
+    if (!visited[*bat])
+      depthFromHere += atom_depth( depth, *bat, topIn, visited, maxdepth );
+  }
+  return depthFromHere;
+}
+
+static inline double wrap360(double phi) {
+  if (phi > Constants::PI)
+    return phi - Constants::TWOPI;
+  else if (phi < -Constants::PI)
+    return phi + Constants::TWOPI;
+  else
+    return phi;
+}
+
+/** Attempt to assign a reasonable value for phi internal coordinate for atom i
+  * given that atoms j k and l have known positions.
+  *   j - k
+  *  /     \
+  * i       l
+  */
+int Cpptraj::Structure::Model::AssignPhi(double& phi, int ai, int aj, int ak, int al,
+                                         Topology const& topIn, Frame const& frameIn,
+                                         std::vector<bool> const& atomPositionKnown,
+                                         BuildAtom const& AtomJ, int debug)
+{
+  // Figure out hybridization and chirality of atom j.
+  if (debug > 0)
+    mprintf("DEBUG: AssignPhi for atom j : %s\n", topIn.AtomMaskName(aj).c_str());
+
+  Atom const& AJ = topIn[aj];
+  if (debug > 0) mprintf("DEBUG:\t\tNbonds: %i\n", AJ.Nbonds());
+  // If atom J only has 2 bonds, ai-aj-ak-al is the only possibility.
+  if (AJ.Nbonds() < 3) {
+    if (debug > 0)
+      mprintf("DEBUG:\t\tFewer than 3 bonds. Setting phi to -180.\n");
+    phi = -180 * Constants::DEGRAD;
+    return 0;
+  }
+
+  // TODO check that atom i actually ends up on the list?
+  std::vector<int> const& priority = AtomJ.Priority();
+  if (debug > 0) {
+    mprintf("DEBUG: Original chirality around J %s is %s\n", topIn.AtomMaskName(aj).c_str(), chiralStr(AtomJ.Chirality()));
+    mprintf("DEBUG:\t\tPriority around J %s(%i):", 
+            topIn.AtomMaskName(aj).c_str(), (int)atomPositionKnown[aj]);
+    for (int idx = 0; idx < AJ.Nbonds(); idx++)
+      mprintf(" %s(%i)", topIn.AtomMaskName(priority[idx]).c_str(), (int)atomPositionKnown[priority[idx]]);
+    mprintf("\n");
+  }
+
+  // Fill in what values we can for known atoms
+  std::vector<double> knownPhi( AJ.Nbonds() );
+  int knownIdx = -1;
+  for (int idx = 0; idx < AJ.Nbonds(); idx++) {
+    int atnum = priority[idx];
+    if (atnum != ak && atomPositionKnown[atnum] &&
+                       atomPositionKnown[aj] &&
+                       atomPositionKnown[ak] &&
+                       atomPositionKnown[al])
+    {
+      knownPhi[idx] = Torsion(frameIn.XYZ(atnum),
+                              frameIn.XYZ(aj),
+                              frameIn.XYZ(ak),
+                              frameIn.XYZ(al));
+      if (debug > 0)
+        mprintf("DEBUG:\t\tKnown phi for %s = %g\n", topIn.AtomMaskName(atnum).c_str(), knownPhi[idx]*Constants::RADDEG);
+      if (knownIdx == -1) knownIdx = idx; // FIXME handle more than 1 known
+    }
+  }
+
+  // If we have to assign an initial phi, make trans the longer branch
+  if (knownIdx == -1) {
+    std::vector<bool> visited = atomPositionKnown;
+    // TODO: Ensure bonded atoms are not yet visited?
+    visited[aj] = true;
+    visited[ak] = true;
+    std::vector<int> depth( AJ.Nbonds() );
+    for (int idx = 0; idx < AJ.Nbonds(); idx++) {
+      int atnum = priority[idx];
+      if (atnum != ak) {
+        int currentDepth = 0;
+        depth[idx] = atom_depth(currentDepth, atnum, topIn, visited, 10);
+        if (debug > 0)
+          mprintf("DEBUG:\t\tAJ %s depth from %s is %i\n",
+                  topIn.AtomMaskName(aj).c_str(), topIn.AtomMaskName(atnum).c_str(), depth[idx]);
+        if (knownIdx == -1 && depth[idx] < 3) {
+          knownIdx = idx;
+          knownPhi[idx] = 0;
+        }
+      }
+    }
+  }
+
+  // The interval will be 360 / (number of bonds - 1)
+  double interval = Constants::TWOPI / (AJ.Nbonds() - 1);
+
+  double startPhi;
+  if (knownIdx == -1) {
+    startPhi = -180*Constants::DEGRAD;
+    if (debug > 0) mprintf("DEBUG:\t\tNo known phi. Setting to %g.\n", startPhi*Constants::RADDEG);
+    knownIdx = 0;
+  } else
+    startPhi = knownPhi[knownIdx];
+
+  if (AtomJ.Chirality() == IS_R) {
+    startPhi = -startPhi;
+    interval = -interval;
+  }
+  if (debug > 0) {
+    mprintf("DEBUG:\t\tStart phi is %g degrees\n", startPhi*Constants::RADDEG);
+    mprintf("DEBUG:\t\tInterval is %g degrees\n", interval * Constants::RADDEG);
+  }
+    
+  // Forward direction
+  double currentPhi = startPhi;
+  for (int idx = knownIdx; idx < AJ.Nbonds(); idx++) {
+    int atnum = priority[idx];
+    if (atnum != ak) {
+      if (atnum == ai) phi = currentPhi;
+      if (debug > 0)
+        mprintf("DEBUG:\t\t\t%s phi= %g\n", topIn.AtomMaskName(atnum).c_str(), currentPhi*Constants::RADDEG);
+      currentPhi += interval;
+      currentPhi = wrap360(currentPhi);
+    }
+  }
+  // Reverse direction
+  currentPhi = startPhi - interval;
+  for (int idx = knownIdx - 1; idx > -1; idx--) {
+    int atnum = priority[idx];
+    if (atnum != ak) {
+      if (atnum == ai) phi = currentPhi;
+      if (debug > 0)
+        mprintf("DEBUG:\t\t\t%s phi= %g\n", topIn.AtomMaskName(atnum).c_str(), currentPhi*Constants::RADDEG);
+      currentPhi -= interval;
+      currentPhi = wrap360(currentPhi);
+    }
+  }
+
+  return 0;
+}
diff --git a/src/Structure/Model.h b/src/Structure/Model.h
new file mode 100644
index 0000000000..3b1f4fafde
--- /dev/null
+++ b/src/Structure/Model.h
@@ -0,0 +1,20 @@
+#ifndef INC_STRUCTURE_MODEL_H
+#define INC_STRUCTURE_MODEL_H
+#include <vector>
+#include "StructureEnum.h"
+class Topology;
+class Frame;
+namespace Cpptraj {
+namespace Structure {
+class BuildAtom;
+/// Routines to generate model parameters TODO needed for modXNA
+namespace Model {
+/// Given atoms J K and L, attempt to assign a reasonable value for phi for atom I
+int AssignPhi(double&, int, int, int, int, Topology const&, Frame const&, std::vector<bool> const&, BuildAtom const&, int);
+/// Given atoms J and K, attempt to assign a reasonable value for theta for atom I
+int AssignTheta(double&, int, int, int, Topology const&, Frame const&, std::vector<bool> const&, int);
+
+} // END namespace Model
+} // END namespace Structure
+} // END namespace Cpptraj
+#endif
diff --git a/src/Structure/OldBuilder.cpp b/src/Structure/OldBuilder.cpp
new file mode 100644
index 0000000000..32f1c22412
--- /dev/null
+++ b/src/Structure/OldBuilder.cpp
@@ -0,0 +1,109 @@
+#include "OldBuilder.h"
+#include "BuildAtom.h"
+#include "Chirality.h"
+#include "Zmatrix.h"
+#include "../CpptrajStdio.h"
+#include "../Frame.h"
+#include "../Topology.h"
+#include <algorithm> // std::copy
+
+using namespace Cpptraj::Structure;
+
+/** CONSTRUCTOR */
+OldBuilder::OldBuilder() :
+  debug_(0)
+{}
+
+/** Combine two units. Fragment 1 will be merged into Fragment 0 and bonded. */
+int OldBuilder::Combine(Topology&       frag0Top, Frame&       frag0frm,
+                     Topology const& frag1Top, Frame const& frag1frm,
+                     int bondAt0, int bondAt1)
+{
+  int natom0 = frag0Top.Natom();
+  int newNatom = natom0 + frag1Top.Natom();
+
+  // Determine which "direction" we will be combining the fragments.
+  // Make atA belong to the smaller fragment. atB fragment will be "known".
+  // Ensure atB index is what it will be after fragments are combined.
+  Barray posKnown( newNatom, false );
+  int atA, atB;
+  if (frag0Top.HeavyAtomCount() < frag1Top.HeavyAtomCount()) {
+    // Fragment 1 is larger
+    atA = bondAt0;
+    atB = bondAt1 + natom0;
+    for (int at = frag0Top.Natom(); at != newNatom; at++)
+      posKnown[at] = true;
+  } else {
+    // Fragment 0 is larger or equal
+    atA = bondAt1 + natom0;
+    atB = bondAt0;
+    for (int at = 0; at != natom0; at++)
+     posKnown[at] = true;
+  }
+
+  // Combine fragment1 into fragment 0 topology
+  Topology& combinedTop = frag0Top;
+  combinedTop.AppendTop( frag1Top );
+  // Combined fragment1 into fragment 0 coords.
+  // Need to save the original coords in frame0 since SetupFrameV does not preserve.
+  double* tmpcrd0 = new double[natom0*3];
+  std::copy( frag0frm.xAddress(), frag0frm.xAddress()+frag0frm.size(), tmpcrd0 );
+  frag0frm.SetupFrameV( combinedTop.Atoms(), CoordinateInfo(frag0frm.BoxCrd(), false, false, false));
+  std::copy( tmpcrd0, tmpcrd0+natom0*3, frag0frm.xAddress() );
+  std::copy( frag1frm.xAddress(), frag1frm.xAddress()+frag1frm.size(), frag0frm.xAddress()+natom0*3 );
+  Frame& CombinedFrame = frag0frm;
+  delete[] tmpcrd0;
+
+  int chiralityDebug;
+  if (debug_ < 1)
+    chiralityDebug = 0;
+  else
+    chiralityDebug = debug_ - 1;
+  // Get the chirality around each atom before the bond is added.
+  BuildAtom AtomA;
+  if (combinedTop[atA].Nbonds() > 2)
+    AtomA.SetChirality( DetermineChirality(atA, combinedTop, CombinedFrame, chiralityDebug) );
+  if (debug_ > 0)
+    mprintf("DEBUG:\tAtom %4s chirality %6s\n", combinedTop.AtomMaskName(atA).c_str(), chiralStr(AtomA.Chirality()));
+  BuildAtom AtomB;
+  if (combinedTop[atB].Nbonds() > 2)
+    AtomB.SetChirality( DetermineChirality(atB, combinedTop, CombinedFrame, chiralityDebug) );
+  if (debug_ > 0)
+    mprintf("DEBUG:\tAtom %4s chirality %6s\n", combinedTop.AtomMaskName(atB).c_str(), chiralStr(AtomB.Chirality()));
+
+  // Create the bond
+  if (debug_ > 0)
+    mprintf("DEBUG: Bonding atom %s to %s\n", combinedTop.AtomMaskName(atA).c_str(), combinedTop.AtomMaskName(atB).c_str());
+  combinedTop.AddBond( atA, atB ); // TODO pseudo-parameter?
+  // // Regenerate the molecule info FIXME should Topology just do this?
+  if (combinedTop.DetermineMolecules()) return 1;
+
+  // Determine priorities
+  if (combinedTop[atA].Nbonds() > 2) {
+    //AtomA.SetNbonds(combinedTop[atA].Nbonds());
+    SetPriority(AtomA.ModifyPriority(), atA, combinedTop, CombinedFrame, chiralityDebug);
+  }
+  if (combinedTop[atB].Nbonds() > 2) {
+    //AtomB.SetNbonds(combinedTop[atB].Nbonds());
+    SetPriority(AtomB.ModifyPriority(), atB, combinedTop, CombinedFrame, chiralityDebug);
+  }
+
+  // Generate Zmatrix only for ICs involving bonded atoms
+  Zmatrix bondZmatrix;
+
+  bondZmatrix.SetDebug( debug_ );
+  if (bondZmatrix.SetupICsAroundBond(atA, atB, CombinedFrame, combinedTop, posKnown, AtomA, AtomB)) {
+    mprinterr("Error: Zmatrix setup for ICs around %s and %s failed.\n",
+              combinedTop.AtomMaskName(atA).c_str(),
+              combinedTop.AtomMaskName(atB).c_str());
+    return 1;
+  }
+  if (debug_ > 0)
+    bondZmatrix.print(&combinedTop);
+  if (bondZmatrix.SetToFrame( CombinedFrame, posKnown )) {
+    mprinterr("Error: Conversion from bondZmatrix to Cartesian coords failed.\n");
+    return 1;
+  }
+
+  return 0;
+}
diff --git a/src/Structure/OldBuilder.h b/src/Structure/OldBuilder.h
new file mode 100644
index 0000000000..b37339376e
--- /dev/null
+++ b/src/Structure/OldBuilder.h
@@ -0,0 +1,27 @@
+#ifndef INC_STRUCTURE_OLDBUILDER_H
+#define INC_STRUCTURE_OLDBUILDER_H
+#include <vector>
+class Topology;
+class Frame;
+namespace Cpptraj {
+namespace Structure {
+/// This is the old (<v7) builder used to attach different topology/frame combos using internal coordinates.
+/** Retaining this class for backwards compatibility with modXNA TODO. */
+class OldBuilder {
+  public:
+    /// CONSTRUCTOR
+    OldBuilder();
+    /// Combine second fragment into first fragment and bond
+    int Combine(Topology&, Frame&, Topology const&, Frame const&, int, int);
+    /// Set debug
+    void SetDebug(int d) { debug_ = d; }
+  private:
+    typedef std::vector<bool> Barray;
+
+    int debug_;
+};
+
+
+}
+}
+#endif
diff --git a/src/Structure/Zmatrix.cpp b/src/Structure/Zmatrix.cpp
index 2d2dab8ad3..fbca3e7493 100644
--- a/src/Structure/Zmatrix.cpp
+++ b/src/Structure/Zmatrix.cpp
@@ -10,6 +10,8 @@
 #include "../DistRoutines.h"
 #include "../Topology.h"
 #include "../TorsionRoutines.h"
+#include "BuildAtom.h" // TODO needed for modXNA
+#include "Model.h" // TODO needed for modXNA
 
 using namespace Cpptraj::Structure;
 
@@ -1411,3 +1413,208 @@ int Zmatrix::SetToFrame(Frame& frameOut, Barray& hasPosition) const {
 */
   return 0;
 }
+
+// =============================================================================
+/** Given two bonded atoms A and B, where B has a depth of at least 2
+  * (i.e., it is possible to have B be atom J where we can form J-K-L),
+  * set up a complete set of internal coordinates involving A and B in the
+  * direction of atom A. This means all internal coordinates with A and B
+  * as I and J (should be only 1), as J and K, and as K and L.
+  * TODO This routine is included for backwards compatibility with modXNA.
+  */
+int Zmatrix::SetupICsAroundBond(int atA, int atB, Frame const& frameIn, Topology const& topIn,
+                                std::vector<bool> const& atomPositionKnown,
+                                BuildAtom const& AtomA, BuildAtom const& AtomB)
+{
+  if (debug_ > 0)
+    mprintf("DEBUG: SetupICsAroundBond: atA= %s (%i)  atB= %s (%i)\n",
+            topIn.AtomMaskName(atA).c_str(), atA+1,
+            topIn.AtomMaskName(atB).c_str(), atB+1);
+  IC_.clear();
+
+  Barray hasIC( topIn.Natom(), false );
+  unsigned int nHasIC = 0;
+  // Mark known atoms as already having IC
+  for (std::vector<bool>::const_iterator it = atomPositionKnown.begin();
+                                         it != atomPositionKnown.end(); ++it)
+  {
+    if (*it) MARK( it - atomPositionKnown.begin(), hasIC, nHasIC );
+  }
+
+  // First, make sure atom B as a bond depth of at least 2.
+  // Choose K and L atoms given atA is I and atB is J.
+  Atom const& AJ = topIn[atB];
+  typedef std::pair<int,int> Apair;
+  std::vector<Apair> KLpairs;
+  for (Atom::bond_iterator kat = AJ.bondbegin(); kat != AJ.bondend(); ++kat)
+  {
+    if (*kat != atA) {
+      //mprintf("DEBUG: kat= %s\n", topIn.AtomMaskName(*kat).c_str());
+      Atom const& AK = topIn[*kat];
+      for (Atom::bond_iterator lat = AK.bondbegin(); lat != AK.bondend(); ++lat)
+      {
+        if (*lat != atB && *lat != atA) {
+          //mprintf("DEBUG: lat= %s\n", topIn.AtomMaskName(*lat).c_str());
+          KLpairs.push_back( Apair(*kat, *lat) );
+        }
+      }
+    }
+  }
+  if (debug_ > 0) {
+    for (std::vector<Apair>::const_iterator it = KLpairs.begin();
+                                            it != KLpairs.end(); ++it)
+      mprintf("DEBUG:\t\tKL pair %s - %s\n", topIn.AtomMaskName(it->first).c_str(),
+              topIn.AtomMaskName(it->second).c_str());
+  }
+  if (KLpairs.empty()) {
+    mprinterr("Error: SetFromFrameAroundBond(): Could not find an atom pair bonded to atom %s\n",
+              topIn.AtomMaskName(atB).c_str());
+    return 1;
+  }
+  // TODO be smarter about how K and L are selected?
+  double maxMass = topIn[KLpairs[0].first].Mass() + topIn[KLpairs[0].second].Mass();
+  unsigned int maxIdx = 0;
+  for (unsigned int idx = 1; idx < KLpairs.size(); idx++) {
+    double sumMass = topIn[KLpairs[idx].first].Mass() + topIn[KLpairs[idx].second].Mass();
+    if (sumMass > maxMass) {
+      maxMass = sumMass;
+      maxIdx = idx;
+    }
+  }
+  int atk0 = KLpairs[maxIdx].first;
+  int atl0 = KLpairs[maxIdx].second;
+  int modelDebug = 0;
+  if (debug_ > 0) {
+    mprintf("DEBUG: Chosen KL pair: %s - %s\n",topIn.AtomMaskName(atk0).c_str(),
+              topIn.AtomMaskName(atl0).c_str());
+    modelDebug = debug_ - 1;
+  }
+  // ---- I J: Set dist, theta, phi for atA atB K L internal coord ---
+  if (debug_ > 0)
+    mprintf("DEBUG: IC (i j) %i - %i - %i - %i\n", atA+1, atB+1, atk0+1, atl0+1);
+  double newDist = Atom::GetBondLength( topIn[atA].Element(), topIn[atB].Element() );
+  if (debug_ > 0) mprintf("DEBUG:\t\tnewDist= %g\n", newDist);
+  double newTheta = 0;
+  if (Cpptraj::Structure::Model::AssignTheta(newTheta, atA, atB, atk0, topIn, frameIn, atomPositionKnown, modelDebug)) {
+    mprinterr("Error: theta (i j) assignment failed.\n");
+    return 1;
+  }
+  if (debug_ > 0) mprintf("DEBUG:\t\tnewTheta = %g\n", newTheta*Constants::RADDEG);
+  double newPhi = 0;
+  if (Cpptraj::Structure::Model::AssignPhi(newPhi, atA, atB, atk0, atl0, topIn, frameIn,
+                                           atomPositionKnown, AtomB, modelDebug))
+  {
+    mprinterr("Error: phi (i j) assignment failed.\n");
+    return 1;
+  }
+  if (debug_ > 0) mprintf("DEBUG:\t\tnewPhi = %g\n", newPhi*Constants::RADDEG);
+  IC_.push_back(InternalCoords( atA, atB, atk0, atl0, newDist, newTheta*Constants::RADDEG, newPhi*Constants::RADDEG ));
+  if (debug_ > 0) {
+    mprintf("DEBUG: MODEL I J IC: ");
+    IC_.back().printIC(topIn);
+  }
+  MARK( atA, hasIC, nHasIC );
+  // ----- J K: Set up ICs for X atA atB K ---------------------------
+  Atom const& AJ1 = topIn[atA];
+  int ati = -1;
+  //Atom const& AK1 = topIn[atB];
+  //Atom const& AL1 = topIn[atk0];
+  for (Atom::bond_iterator iat = AJ1.bondbegin(); iat != AJ1.bondend(); ++iat)
+  {
+    if (*iat != atB) {
+      if (ati == -1) ati = *iat;
+      // Use existing dist
+      newDist = sqrt( DIST2_NoImage(frameIn.XYZ(*iat), frameIn.XYZ(atA)) );
+      // Set theta for I atA atB
+      newTheta = 0;
+      if (Cpptraj::Structure::Model::AssignTheta(newTheta, *iat, atA, atB, topIn, frameIn, atomPositionKnown, modelDebug)) {
+        mprinterr("Error: theta (j k) assignment failed.\n");
+        return 1;
+      }
+      if (debug_ > 0)
+        mprintf("DEBUG:\t\tnewTheta = %g\n", newTheta*Constants::RADDEG);
+      // Set phi for I atA atB K
+      newPhi = 0;
+      if (Cpptraj::Structure::Model::AssignPhi(newPhi, *iat, atA, atB, atk0, topIn, frameIn,
+                                               atomPositionKnown, AtomA, modelDebug))
+      {
+        mprinterr("Error: phi (j k) assignment failed.\n");
+        return 1;
+      }
+      if (debug_ > 0)
+        mprintf("DEBUG:\t\tnewPhi = %g\n", newPhi*Constants::RADDEG);
+      IC_.push_back(InternalCoords( *iat, atA, atB, atk0, newDist, newTheta*Constants::RADDEG, newPhi*Constants::RADDEG ));
+      if (debug_ > 0) {
+        mprintf("DEBUG: MODEL J K IC: ");
+        IC_.back().printIC(topIn);
+      }
+      MARK( *iat, hasIC, nHasIC );
+      // ----- K L: Set up ICs for X iat atA atB ---------------------
+      /*Atom const& AJ2 = topIn[*iat];
+      for (Atom::bond_iterator i2at = AJ2.bondbegin(); i2at != AJ2.bondend(); ++i2at)
+      {
+        if (*i2at != atA && *i2at != atB && !hasIC[*i2at]) {
+          // Use existing dist and theta
+          newDist = sqrt( DIST2_NoImage(frameIn.XYZ(*i2at), frameIn.XYZ(*iat)) );
+          newTheta = CalcAngle( frameIn.XYZ(*i2at), frameIn.XYZ(*iat), frameIn.XYZ(atA) );
+          // Set phi for X iat atA atB
+          newPhi = 0;
+          if (Cpptraj::Structure::Model::AssignPhi(newPhi, *i2at, *iat, atA, atB, topIn, frameIn, atomPositionKnown, atomChirality)) {
+            mprinterr("Error: phi (k l) assignment failed.\n");
+            return 1;
+          }
+          mprintf("DEBUG:\t\tnewPhi = %g\n", newPhi*Constants::RADDEG);
+          IC_.push_back(InternalCoords( *i2at, *iat, atA, atB, newDist, newTheta*Constants::RADDEG, newPhi*Constants::RADDEG ));
+          mprintf("DEBUG: MODEL K L IC: ");
+          IC_.back().printIC(topIn);
+          MARK( *i2at, hasIC, nHasIC );
+          // Trace from atA *iat *i2at outwards
+          ToTrace.push_back(atA);
+          ToTrace.push_back(*iat);
+          ToTrace.push_back(*i2at);
+          //if (traceMol(atA, *iat, *i2at, frameIn, topIn, topIn.Natom(), nHasIC, hasIC)) return 1;
+        }
+      } */
+    }
+  }
+  // Handle remaining atoms.
+  if (AJ1.Nbonds() > 1) {
+    if (AJ1.Nbonds() == 2) {
+      if (debug_ > 0) mprintf("DEBUG: 2 bonds to %s.\n", topIn.AtomMaskName(atA).c_str());
+      if (traceMol(atB, atA, ati, frameIn, topIn, topIn.Natom(), nHasIC, hasIC))
+        return 1;
+    } else {
+      // 3 or more bonds
+      std::vector<int> const& priority = AtomA.Priority();
+      int at0 = -1;
+      int at1 = -1;
+      std::vector<int> remainingAtoms;
+      if (debug_ > 0) mprintf("DEBUG: %i bonds to %s\n", AJ1.Nbonds(), topIn.AtomMaskName(atA).c_str());
+      for (std::vector<int>::const_iterator it = priority.begin(); it != priority.end(); ++it) {
+        if (*it != atB) {
+          if (debug_ > 0) mprintf("DEBUG:\t\t%s\n", topIn.AtomMaskName(*it).c_str());
+          if (at0 == -1)
+            at0 = *it;
+          else if (at1 == -1)
+            at1 = *it;
+          else
+            remainingAtoms.push_back( *it );
+        }
+      }
+      // at0 atA at1
+      if (traceMol(at1, atA, at0, frameIn, topIn, topIn.Natom(), nHasIC, hasIC))
+        return 1;
+      // at1 atA, at0
+      if (traceMol(at0, atA, at1, frameIn, topIn, topIn.Natom(), nHasIC, hasIC))
+        return 1;
+      // Remaining atoms.
+      for (std::vector<int>::const_iterator it = remainingAtoms.begin(); it != remainingAtoms.end(); ++it) {
+        if (traceMol(at0, atA, *it, frameIn, topIn, topIn.Natom(), nHasIC, hasIC))
+          return 1;
+      }
+    }
+  }
+
+  return 0;
+}
+
diff --git a/src/Structure/Zmatrix.h b/src/Structure/Zmatrix.h
index cd11068403..c33b7cf3f9 100644
--- a/src/Structure/Zmatrix.h
+++ b/src/Structure/Zmatrix.h
@@ -7,6 +7,7 @@ class Topology;
 class Molecule;
 namespace Cpptraj {
 namespace Structure {
+class BuildAtom; // TODO needed for modXNA
 /// Hold internal coordinates for a system.
 class Zmatrix {
     typedef std::vector<InternalCoords> ICarray;
@@ -84,6 +85,10 @@ class Zmatrix {
     static Vec3 AtomIposition(Vec3 const&, Vec3 const&, Vec3 const&, double, double, double);
     /// \return XYZ position of atom I for given internal coordinate
     static Vec3 AtomIposition(InternalCoords const&, Frame const&);
+
+    /// Get internal coordinates around bond in one direction. TODO for modXNA
+    int SetupICsAroundBond(int, int, Frame const&, Topology const&,
+                           std::vector<bool> const&, BuildAtom const&, BuildAtom const&);
   private:
     typedef std::vector<int> Iarray;
     /// Set seeds as 3 consecutive atoms with known positions
diff --git a/src/Structure/structuredepend b/src/Structure/structuredepend
index 2bdebfb6b5..219597caff 100644
--- a/src/Structure/structuredepend
+++ b/src/Structure/structuredepend
@@ -10,6 +10,8 @@ HisProt.o : HisProt.cpp ../ArgList.h ../Atom.h ../AtomMask.h ../Box.h ../Constan
 InternalCoords.o : InternalCoords.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h InternalCoords.h
 LeastSquaresPlane.o : LeastSquaresPlane.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../Molecule.h ../NameType.h ../Parallel.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Unit.h ../Vec3.h LeastSquaresPlane.h
 MetalCenterFinder.o : MetalCenterFinder.cpp ../ArgList.h ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../DistRoutines.h ../FileName.h ../Frame.h ../ImageOption.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h MetalCenterFinder.h
+Model.o : Model.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../TorsionRoutines.h ../Unit.h ../Vec3.h BuildAtom.h Model.h StructureEnum.h
+OldBuilder.o : OldBuilder.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h BuildAtom.h Chirality.h InternalCoords.h OldBuilder.h StructureEnum.h Zmatrix.h
 PdbCleaner.o : PdbCleaner.cpp ../ArgList.h ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h PdbCleaner.h
 RingFinder.o : RingFinder.cpp ../Atom.h ../AtomMask.h ../Box.h ../CharMask.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h RingFinder.h
 Solvate.o : Solvate.cpp ../Action.h ../ActionState.h ../Action_AutoImage.h ../ArgList.h ../AssociatedData.h ../Atom.h ../AtomMask.h ../BaseIOtype.h ../Box.h ../CharMask.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DataFile.h ../DataFileList.h ../DataSet.h ../DataSetList.h ../DataSet_Coords.h ../DataSet_Coords_REF.h ../Dimension.h ../DispatchObject.h ../FileIO.h ../FileName.h ../FileTypes.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../MetaData.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Parm/../AtomType.h ../Parm/../CmapParmHolder.h ../Parm/../Constants.h ../Parm/../FileName.h ../Parm/../NameType.h ../Parm/../ParameterTypes.h ../Parm/../Parm/ParmEnum.h ../Parm/../TypeNameHolder.h ../Parm/DihedralParmHolder.h ../Parm/ImproperParmHolder.h ../Parm/ParameterSet.h ../Parm/ParmEnum.h ../Parm/ParmHolder.h ../Range.h ../ReferenceFrame.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../TextFormat.h ../Timer.h ../Topology.h ../Unit.h ../Vec3.h Solvate.h
@@ -19,4 +21,4 @@ Sugar.o : Sugar.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../Coordinat
 SugarBuilder.o : SugarBuilder.cpp ../ArgList.h ../Atom.h ../AtomMap.h ../AtomMask.h ../Box.h ../CharMask.h ../Constants.h ../CoordinateInfo.h ../CpptrajFile.h ../CpptrajStdio.h ../DistRoutines.h ../FileIO.h ../FileName.h ../Frame.h ../ImageOption.h ../MapAtom.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../TorsionRoutines.h ../Unit.h ../Vec3.h Chirality.h FunctionalGroup.h FxnGroupBuilder.h LinkAtom.h ResStatArray.h StructureEnum.h StructureRoutines.h Sugar.h SugarBuilder.h SugarLinkAtoms.h SugarToken.h
 SugarLinkAtoms.o : SugarLinkAtoms.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../FileName.h ../Frame.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../StringRoutines.h ../SymbolExporting.h ../Topology.h ../Unit.h ../Vec3.h LinkAtom.h SugarLinkAtoms.h
 SugarToken.o : SugarToken.cpp ../ArgList.h ../CpptrajStdio.h SugarToken.h
-Zmatrix.o : Zmatrix.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../DistRoutines.h ../FileName.h ../Frame.h ../ImageOption.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../TorsionRoutines.h ../Unit.h ../Vec3.h GenerateConnectivityArrays.h InternalCoords.h Zmatrix.h
+Zmatrix.o : Zmatrix.cpp ../Atom.h ../AtomMask.h ../Box.h ../Constants.h ../CoordinateInfo.h ../CpptrajStdio.h ../DistRoutines.h ../FileName.h ../Frame.h ../ImageOption.h ../MaskToken.h ../Matrix_3x3.h ../ModXNA_Info.h ../Molecule.h ../NameType.h ../Parallel.h ../ParameterTypes.h ../Range.h ../ReplicaDimArray.h ../Residue.h ../Segment.h ../SymbolExporting.h ../Topology.h ../TorsionRoutines.h ../Unit.h ../Vec3.h BuildAtom.h GenerateConnectivityArrays.h InternalCoords.h Model.h StructureEnum.h Zmatrix.h
diff --git a/src/Structure/structurefiles b/src/Structure/structurefiles
index 3f53f820b1..ea2949bf19 100644
--- a/src/Structure/structurefiles
+++ b/src/Structure/structurefiles
@@ -12,6 +12,8 @@ STRUCTURE_SOURCES= \
   InternalCoords.cpp \
   LeastSquaresPlane.cpp \
   MetalCenterFinder.cpp \
+  Model.cpp \
+  OldBuilder.cpp \
   PdbCleaner.cpp \
   RingFinder.cpp \
   Solvate.cpp \
diff --git a/src/Topology.cpp b/src/Topology.cpp
index e3142ef857..76c75ef1f1 100644
--- a/src/Topology.cpp
+++ b/src/Topology.cpp
@@ -1932,6 +1932,9 @@ Topology* Topology::ModifyByMap(std::vector<int> const& MapIn, bool setupFullPar
   // Determine number of extra points
   newParm->DetermineNumExtraPoints();
 
+  // modXNA
+  newParm->modxna_ = modxna_;
+
   return newParm;
 }
 
diff --git a/src/Traj_Mol2File.cpp b/src/Traj_Mol2File.cpp
index 98bfbffcd5..4786b82f04 100644
--- a/src/Traj_Mol2File.cpp
+++ b/src/Traj_Mol2File.cpp
@@ -203,8 +203,13 @@ int Traj_Mol2File::setupTrajout(FileName const& fname, Topology* trajParm,
     }
   }
   // Set Title
-  if (Title().empty())
-    SetTitle("Cpptraj Generated mol2 file.");
+  if (Title().empty()) {
+    if (mol2Top_->Modxna().HasModxna()) {
+      mprintf("\tModXNA data present and no title set; writing ModXNA metadata in title.\n");
+      SetTitle(mol2Top_->Modxna().GenMetadataString());
+    } else
+      SetTitle("Cpptraj Generated mol2 file.");
+  }
   file_.SetMol2Title( Title() );
   // Set up number of bonds
   file_.SetMol2Nbonds( mol2Top_->Bonds().size() + mol2Top_->BondsH().size() );
diff --git a/src/Version.h b/src/Version.h
index 1fe64fc8a6..09fd3ebeda 100644
--- a/src/Version.h
+++ b/src/Version.h
@@ -12,7 +12,7 @@
  * Whenever a number that precedes <revision> is incremented, all subsequent
  * numbers should be reset to 0.
  */
-#define CPPTRAJ_INTERNAL_VERSION "V7.00.0"
+#define CPPTRAJ_INTERNAL_VERSION "V7.1.0"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif
diff --git a/src/cpptrajdepend b/src/cpptrajdepend
index c633777ac6..053e6367f1 100644
--- a/src/cpptrajdepend
+++ b/src/cpptrajdepend
@@ -364,7 +364,7 @@ Exec_ReadInput.o : Exec_ReadInput.cpp Action.h ActionList.h ActionState.h Analys
 Exec_RotateDihedral.o : Exec_RotateDihedral.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DihedralSearch.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_RotateDihedral.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h Unit.h Vec3.h
 Exec_RunAnalysis.o : Exec_RunAnalysis.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h Cmd.h CmdList.h Command.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_RunAnalysis.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h Unit.h Vec3.h
 Exec_ScaleDihedralK.o : Exec_ScaleDihedralK.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_ScaleDihedralK.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h Unit.h Vec3.h
-Exec_Sequence.o : Exec_Sequence.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h AssociatedData_Connect.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CmapParmHolder.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Parameters.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Sequence.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Parm/AssignParams.h Parm/DihedralParmHolder.h Parm/GB_Params.h Parm/ImproperParmHolder.h Parm/ParameterSet.h Parm/ParmEnum.h Parm/ParmHolder.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h Structure/Builder.h Structure/Creator.h Structure/StructureEnum.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h
+Exec_Sequence.o : Exec_Sequence.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h AssociatedData_Connect.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h CmapParmHolder.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Parameters.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Sequence.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Parm/AssignParams.h Parm/DihedralParmHolder.h Parm/GB_Params.h Parm/ImproperParmHolder.h Parm/ParameterSet.h Parm/ParmEnum.h Parm/ParmHolder.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h Structure/Builder.h Structure/Creator.h Structure/OldBuilder.h Structure/StructureEnum.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h
 Exec_SequenceAlign.o : Exec_SequenceAlign.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_SequenceAlign.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Unit.h Vec3.h
 Exec_Set.o : Exec_Set.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_StringVar.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Set.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h Unit.h Vec3.h
 Exec_Show.o : Exec_Show.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_StringVar.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h Exec.h Exec_Show.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h Unit.h Vec3.h
@@ -448,7 +448,7 @@ Parm_Amber.o : Parm_Amber.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h BondSearc
 Parm_CIF.o : Parm_CIF.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h BondSearch.h Box.h BufferedLine.h CIFfile.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h ParmIO.h Parm_CIF.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h Topology.h Unit.h Vec3.h
 Parm_CharmmPsf.o : Parm_CharmmPsf.cpp ArgList.h Atom.h AtomMask.h AtomType.h BaseIOtype.h BondSearch.h Box.h BufferedLine.h CharmmParamFile.h CmapParmHolder.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Mol.h Molecule.h NameType.h Parallel.h ParameterTypes.h Parm/AssignParams.h Parm/DihedralParmHolder.h Parm/ImproperParmHolder.h Parm/ParameterSet.h Parm/ParmEnum.h Parm/ParmHolder.h ParmIO.h Parm_CharmmPsf.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 Parm_Gromacs.o : Parm_Gromacs.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h BondSearch.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h ParmIO.h Parm_Gromacs.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h Topology.h Unit.h Vec3.h
-Parm_Mol2.o : Parm_Mol2.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h BondSearch.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Mol2File.h Molecule.h NameType.h Parallel.h ParameterTypes.h ParmIO.h Parm_Mol2.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h Unit.h Vec3.h
+Parm_Mol2.o : Parm_Mol2.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h BondSearch.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Mol2File.h Molecule.h NameType.h Parallel.h ParameterTypes.h Parm/Bond_Params.h ParmIO.h Parm_Mol2.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Topology.h Unit.h Vec3.h
 Parm_PDB.o : Parm_PDB.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h BondSearch.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h PDBfile.h Parallel.h ParameterTypes.h ParmIO.h Parm_PDB.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h Timer.h Topology.h Unit.h Vec3.h
 Parm_SDF.o : Parm_SDF.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h BondSearch.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h ParmIO.h Parm_SDF.h Range.h ReplicaDimArray.h Residue.h SDFfile.h Segment.h SymbolExporting.h Topology.h Unit.h Vec3.h
 Parm_Tinker.o : Parm_Tinker.cpp ArgList.h Atom.h AtomMask.h BaseIOtype.h BondSearch.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h ParmIO.h Parm_Tinker.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TinkerFile.h Topology.h Unit.h Vec3.h
@@ -499,6 +499,8 @@ Structure/HisProt.o : Structure/HisProt.cpp ArgList.h Atom.h AtomMask.h Box.h Co
 Structure/InternalCoords.o : Structure/InternalCoords.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/InternalCoords.h SymbolExporting.h Topology.h Unit.h Vec3.h
 Structure/LeastSquaresPlane.o : Structure/LeastSquaresPlane.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ReplicaDimArray.h Residue.h Segment.h Structure/LeastSquaresPlane.h SymbolExporting.h Unit.h Vec3.h
 Structure/MetalCenterFinder.o : Structure/MetalCenterFinder.cpp ArgList.h Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h FileName.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/MetalCenterFinder.h SymbolExporting.h Topology.h Unit.h Vec3.h
+Structure/Model.o : Structure/Model.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/BuildAtom.h Structure/Model.h Structure/StructureEnum.h SymbolExporting.h Topology.h TorsionRoutines.h Unit.h Vec3.h
+Structure/OldBuilder.o : Structure/OldBuilder.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/BuildAtom.h Structure/Chirality.h Structure/InternalCoords.h Structure/OldBuilder.h Structure/StructureEnum.h Structure/Zmatrix.h SymbolExporting.h Topology.h Unit.h Vec3.h
 Structure/PdbCleaner.o : Structure/PdbCleaner.cpp ArgList.h Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/PdbCleaner.h SymbolExporting.h Topology.h Unit.h Vec3.h
 Structure/RingFinder.o : Structure/RingFinder.cpp Atom.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/RingFinder.h SymbolExporting.h Topology.h Unit.h Vec3.h
 Structure/Solvate.o : Structure/Solvate.cpp Action.h ActionState.h Action_AutoImage.h ArgList.h AssociatedData.h Atom.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h Structure/../Parm/../AtomType.h Structure/../Parm/../CmapParmHolder.h Structure/../Parm/../Constants.h Structure/../Parm/../FileName.h Structure/../Parm/../NameType.h Structure/../Parm/../ParameterTypes.h Structure/../Parm/../Parm/ParmEnum.h Structure/../Parm/../TypeNameHolder.h Structure/../Parm/DihedralParmHolder.h Structure/../Parm/ImproperParmHolder.h Structure/../Parm/ParameterSet.h Structure/../Parm/ParmEnum.h Structure/../Parm/ParmHolder.h Structure/Solvate.h SymbolExporting.h TextFormat.h Timer.h Topology.h Unit.h Vec3.h
@@ -508,7 +510,7 @@ Structure/Sugar.o : Structure/Sugar.cpp Atom.h AtomMask.h Box.h Constants.h Coor
 Structure/SugarBuilder.o : Structure/SugarBuilder.cpp ArgList.h Atom.h AtomMap.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DistRoutines.h FileIO.h FileName.h Frame.h ImageOption.h MapAtom.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/Chirality.h Structure/FunctionalGroup.h Structure/FxnGroupBuilder.h Structure/LinkAtom.h Structure/ResStatArray.h Structure/StructureEnum.h Structure/StructureRoutines.h Structure/Sugar.h Structure/SugarBuilder.h Structure/SugarLinkAtoms.h Structure/SugarToken.h SymbolExporting.h Topology.h TorsionRoutines.h Unit.h Vec3.h
 Structure/SugarLinkAtoms.o : Structure/SugarLinkAtoms.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h Structure/LinkAtom.h Structure/SugarLinkAtoms.h SymbolExporting.h Topology.h Unit.h Vec3.h
 Structure/SugarToken.o : Structure/SugarToken.cpp ArgList.h CpptrajStdio.h Structure/SugarToken.h
-Structure/Zmatrix.o : Structure/Zmatrix.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h FileName.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/GenerateConnectivityArrays.h Structure/InternalCoords.h Structure/Zmatrix.h SymbolExporting.h Topology.h TorsionRoutines.h Unit.h Vec3.h
+Structure/Zmatrix.o : Structure/Zmatrix.cpp Atom.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h FileName.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/BuildAtom.h Structure/GenerateConnectivityArrays.h Structure/InternalCoords.h Structure/Model.h Structure/StructureEnum.h Structure/Zmatrix.h SymbolExporting.h Topology.h TorsionRoutines.h Unit.h Vec3.h
 StructureCheck.o : StructureCheck.cpp Atom.h AtomMask.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DistRoutines.h ExclusionArray.h FileIO.h FileName.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h PairList.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/LeastSquaresPlane.h Structure/RingFinder.h StructureCheck.h SymbolExporting.h Timer.h Topology.h Unit.h Vec3.h
 StructureMapper.o : StructureMapper.cpp AssociatedData.h Atom.h AtomMap.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h FileIO.h FileName.h Frame.h MapAtom.h MaskToken.h Matrix_3x3.h MetaData.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h Structure/Chirality.h Structure/StructureEnum.h StructureMapper.h SymbolExporting.h TextFormat.h Topology.h TorsionRoutines.h Unit.h Vec3.h
 SymmetricRmsdCalc.o : SymmetricRmsdCalc.cpp ArrayIterator.h Atom.h AtomMap.h AtomMask.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h DistRoutines.h FileName.h Frame.h Hungarian.h ImageOption.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h ModXNA_Info.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h SymmetricRmsdCalc.h Topology.h Unit.h Vec3.h
diff --git a/test/Makefile b/test/Makefile
index b781b9f8f0..d0859c944d 100644
--- a/test/Makefile
+++ b/test/Makefile
@@ -567,6 +567,9 @@ test.mutate:
 test.build:
 	@-cd Test_Build && ./RunTest.sh $(OPT)
 
+test.modxna:
+	@-cd Test_ModXNA && ./RunTest.sh $(OPT)
+
 # Only GPU-specific tests should go here
 GPUTESTS=test.closest \
          test.solvent \
@@ -753,7 +756,8 @@ COMPLETETESTS=test.general \
               test.minmaxdist \
               test.watershellclosest \
               test.mutate \
-              test.build
+              test.build \
+              test.modxna
 
 test.all:
 	$(MAKE) test.complete summary
diff --git a/test/Test_ModXNA/RunTest.sh b/test/Test_ModXNA/RunTest.sh
new file mode 100755
index 0000000000..a076d2389c
--- /dev/null
+++ b/test/Test_ModXNA/RunTest.sh
@@ -0,0 +1,44 @@
+#!/bin/bash
+
+. ../MasterTest.sh
+
+TESTNAME='ModXNA tests'
+
+# These tests will help ensure compatability with the modxna.sh driver script.
+
+CleanFiles cpptraj.in tmp.bonds tmp.Nucleotide.mol2
+
+INPUT='-i cpptraj.in'
+
+UNITNAME="Ensure 'sequence' works for ModXNA"
+CheckFor maxthreads 1
+if [ $? -eq 0 ] ; then
+  cat > cpptraj.in <<EOF
+parm tmp.BackboneSugar.mol2
+# Check that the RK/REQ column that ModXNA expects is written
+bonds @C1' @/H out tmp.bonds
+loadcrd tmp.BackboneSugar.mol2 name BackboneSugar parm tmp.BackboneSugar.mol2
+parm tmp.base-striped.mol2
+loadcrd tmp.base-striped.mol2 name Base parm tmp.base-striped.mol2
+dataset connect BackboneSugar headmask @C1'
+dataset connect Base tailmask @N9
+sequence Base BackboneSugar name Nucleotide
+change crdset Nucleotide mergeres firstres 1 lastres 2
+change crdset Nucleotide resname from * to AD3
+change crdset Nucleotide oresnums of :1 min 1 max 1
+# Charge on all atoms but O3' and HO3'
+set Q1 = crdset Nucleotide charge inmask !@O3',HO3'
+# Charge on all atoms plus O3'
+Q2 = \$Q1 + -0.4
+# What does charge on HO3' need to be to get target total 3cap charge?
+Q3 = -0.679652 - \$Q2
+change crdset Nucleotide charge of @HO3' to \$Q3
+change crdset Nucleotide charge of @O3' to -0.4
+crdout Nucleotide tmp.Nucleotide.mol2
+EOF
+  RunCpptraj "$UNITNAME"
+  DoTest tmp.Nucleotide.mol2.save tmp.Nucleotide.mol2
+  DoTest tmp.bonds.save tmp.bonds
+fi
+
+EndTest
diff --git a/test/Test_ModXNA/tmp.BackboneSugar.mol2 b/test/Test_ModXNA/tmp.BackboneSugar.mol2
new file mode 100644
index 0000000000..00ef0ecff5
--- /dev/null
+++ b/test/Test_ModXNA/tmp.BackboneSugar.mol2
@@ -0,0 +1,49 @@
+@<TRIPOS>MOLECULE
+Cpptraj Generated mol2 file.
+   19    19     1     0     0
+SMALL
+USER_CHARGES
+
+
+@<TRIPOS>ATOM
+      1 P          -2.7515    2.5361    1.1902 P          1 AD3      1.165356
+      2 OP1        -3.2333    3.1586   -0.0665 O2         1 AD3     -0.776462
+      3 OP2        -2.2870    3.4684    2.2360 O2         1 AD3     -0.776462
+      4 O5'        -1.6064    1.4790    0.8528 OS         1 AD3     -0.495631
+      5 C5'        -1.8840    0.8600   -0.4060 CJ         1 AD3     -0.008126
+      6 H5'        -1.8990    1.6040   -1.1900 H1         1 AD3      0.062001
+      7 H5''       -2.8510    0.3600   -0.4020 H1         1 AD3      0.062001
+      8 C4'        -0.7910   -0.1490   -0.7280 CT         1 AD3      0.133952
+      9 H4'        -1.1020   -0.7030   -1.6130 H1         1 AD3      0.096702
+     10 O4'         0.4020    0.5340   -1.0090 OS         1 AD3     -0.434692
+     11 C1'         1.4780    0.0560   -0.2260 CT         1 AD3      0.035441
+     12 C3'        -0.4660   -1.1380    0.4090 C7         1 AD3      0.058629
+     13 H3'        -1.2760   -1.2130    1.1320 H1         1 AD3      0.080996
+     14 C2'         0.8030   -0.5470    1.0040 CT         1 AD3     -0.100576
+     15 H2'         0.5540    0.2380    1.7080 HC         1 AD3      0.059041
+     16 O3'        -0.1480   -2.4180   -0.0780 OH         1 AD3     -0.616176
+     17 HO3'       -0.8860   -2.7790   -0.5500 HO         1 AD3      0.348720
+     18 H1'         1.9910   -0.7280   -0.7770 H2         1 AD3      0.151137
+     19 H2''        1.4000   -1.3010    1.5000 HC         1 AD3      0.059041
+@<TRIPOS>BOND
+    1     1     4 1
+    2     1     2 1
+    3     1     3 1
+    4     5     8 1
+    5     8    10 1
+    6     8    12 1
+    7    10    11 1
+    8    11    14 1
+    9    12    14 1
+   10    12    16 1
+   11     4     5 1
+   12     5     6 1
+   13     5     7 1
+   14     8     9 1
+   15    11    18 1
+   16    12    13 1
+   17    14    15 1
+   18    14    19 1
+   19    16    17 1
+@<TRIPOS>SUBSTRUCTURE
+      1 AD3               1 ****               0 ****  **** 
diff --git a/test/Test_ModXNA/tmp.Nucleotide.mol2.save b/test/Test_ModXNA/tmp.Nucleotide.mol2.save
new file mode 100644
index 0000000000..92b484d804
--- /dev/null
+++ b/test/Test_ModXNA/tmp.Nucleotide.mol2.save
@@ -0,0 +1,79 @@
+@<TRIPOS>MOLECULE
+Cpptraj Generated mol2 file.
+   33    35     1     0     0
+SMALL
+USER_CHARGES
+
+
+@<TRIPOS>ATOM
+      1 N9          2.7158    0.0933   -1.0181 N*         1 AD3     -0.026798
+      2 C8          2.4424    0.5815   -2.2741 C2         1 AD3      0.160711
+      3 H8          1.4463    0.8516   -2.5490 H5         1 AD3      0.187713
+      4 N7          4.5155    0.1983   -2.3242 NB         1 AD3     -0.617458
+      5 C5          4.8580   -0.2966   -1.0811 CB         1 AD3      0.072505
+      6 C6          6.0644   -0.7292   -0.5214 CA         1 AD3      0.689747
+      7 N6          7.2280   -0.6908   -1.1973 N2         1 AD3     -0.912238
+      8 H61         8.0090   -1.1650   -0.8043 H          1 AD3      0.416728
+      9 H62         7.2128   -0.5332   -2.1787 H          1 AD3      0.416728
+     10 N1          6.0630   -1.1730    0.7300 NC         1 AD3     -0.762348
+     11 C2          4.9178   -1.1767    1.4002 CQ         1 AD3      0.571639
+     12 H2          4.9793   -1.5432    2.4089 H5         1 AD3      0.059804
+     13 N3          3.7288   -0.7952    0.9951 NC         1 AD3     -0.741649
+     14 C4          3.7527   -0.3673   -0.2635 CB         1 AD3      0.380026
+     15 P          -2.7515    2.5361    1.1902 P          1 AD3      1.165356
+     16 OP1        -3.2333    3.1586   -0.0665 O2         1 AD3     -0.776462
+     17 OP2        -2.2870    3.4684    2.2360 O2         1 AD3     -0.776462
+     18 O5'        -1.6064    1.4790    0.8528 OS         1 AD3     -0.495631
+     19 C5'        -1.8840    0.8600   -0.4060 CJ         1 AD3     -0.008126
+     20 H5'        -1.8990    1.6040   -1.1900 H1         1 AD3      0.062001
+     21 H5''       -2.8510    0.3600   -0.4020 H1         1 AD3      0.062001
+     22 C4'        -0.7910   -0.1490   -0.7280 CT         1 AD3      0.133952
+     23 H4'        -1.1020   -0.7030   -1.6130 H1         1 AD3      0.096702
+     24 O4'         0.4020    0.5340   -1.0090 OS         1 AD3     -0.434692
+     25 C1'         1.4780    0.0560   -0.2260 CT         1 AD3      0.035441
+     26 C3'        -0.4660   -1.1380    0.4090 C7         1 AD3      0.058629
+     27 H3'        -1.2760   -1.2130    1.1320 H1         1 AD3      0.080996
+     28 C2'         0.8030   -0.5470    1.0040 CT         1 AD3     -0.100576
+     29 H2'         0.5540    0.2380    1.7080 HC         1 AD3      0.059041
+     30 O3'        -0.1480   -2.4180   -0.0780 OH         1 AD3     -0.400000
+     31 HO3'       -0.8860   -2.7790   -0.5500 HO         1 AD3      0.452888
+     32 H1'         1.9910   -0.7280   -0.7770 H2         1 AD3      0.151137
+     33 H2''        1.4000   -1.3010    1.5000 HC         1 AD3      0.059041
+@<TRIPOS>BOND
+    1     1     2 1
+    2     1    14 1
+    3     2     4 1
+    4     4     5 1
+    5     5     6 1
+    6     5    14 1
+    7     6     7 1
+    8     6    10 1
+    9    10    11 1
+   10    11    13 1
+   11    13    14 1
+   12    15    18 1
+   13    15    16 1
+   14    15    17 1
+   15    19    22 1
+   16    22    24 1
+   17    22    26 1
+   18    24    25 1
+   19    25    28 1
+   20    26    28 1
+   21    26    30 1
+   22    18    19 1
+   23     1    25 1
+   24     2     3 1
+   25     7     8 1
+   26     7     9 1
+   27    11    12 1
+   28    19    20 1
+   29    19    21 1
+   30    22    23 1
+   31    25    32 1
+   32    26    27 1
+   33    28    29 1
+   34    28    33 1
+   35    30    31 1
+@<TRIPOS>SUBSTRUCTURE
+      1 AD3               1 ****               0 ****  **** 
diff --git a/test/Test_ModXNA/tmp.base-striped.mol2 b/test/Test_ModXNA/tmp.base-striped.mol2
new file mode 100644
index 0000000000..b057e0eadd
--- /dev/null
+++ b/test/Test_ModXNA/tmp.base-striped.mol2
@@ -0,0 +1,40 @@
+@<TRIPOS>MOLECULE
+modXNA-fragment:DAA:HEAD01:@N9:HEAD01STRIP:@C1',H1',O4',HHO4,C2',H2'1,H2'3,H2'2
+   14    15     1     0     0
+SMALL
+USER_CHARGES
+
+
+@<TRIPOS>ATOM
+      1 N9         -1.1080   -0.2670    0.1170 N*         1 DAA     -0.026798
+      2 C8         -0.8260   -1.6030   -0.0450 C2         1 DAA      0.160711
+      3 H8         -1.6110   -2.3270   -0.0670 H5         1 DAA      0.187713
+      4 N7          0.4250   -1.8610   -0.1380 NB         1 DAA     -0.617458
+      5 C5          1.0350   -0.6260   -0.0370 CB         1 DAA      0.072505
+      6 C6          2.3710   -0.2130   -0.0570 CA         1 DAA      0.689747
+      7 N6          3.3900   -1.0750   -0.2320 N2         1 DAA     -0.912238
+      8 H61         4.3130   -0.7420   -0.0690 H          1 DAA      0.416728
+      9 H62         3.2170   -2.0510   -0.1560 H          1 DAA      0.416728
+     10 N1          2.6400    1.0800    0.0800 NC         1 DAA     -0.762348
+     11 C2          1.6280    1.9270    0.2190 CQ         1 DAA      0.571639
+     12 H2          1.9070    2.9600    0.3220 H5         1 DAA      0.059804
+     13 N3          0.3410    1.6700    0.2480 NC         1 DAA     -0.741649
+     14 C4          0.0970    0.3690    0.1220 CB         1 DAA      0.380026
+@<TRIPOS>BOND
+    1     1     2 1
+    2     1    14 1
+    3     2     4 1
+    4     4     5 1
+    5     5     6 1
+    6     5    14 1
+    7     6     7 1
+    8     6    10 1
+    9    10    11 1
+   10    11    13 1
+   11    13    14 1
+   12     2     3 1
+   13     7     8 1
+   14     7     9 1
+   15    11    12 1
+@<TRIPOS>SUBSTRUCTURE
+      1 DAA               1 ****               0 ****  **** 
diff --git a/test/Test_ModXNA/tmp.bonds.save b/test/Test_ModXNA/tmp.bonds.save
new file mode 100644
index 0000000000..36707c6afb
--- /dev/null
+++ b/test/Test_ModXNA/tmp.bonds.save
@@ -0,0 +1,2 @@
+#Bnd     RK    REQ Atom1   Atom2   A1 A2 T1 T2
+   4   0.00  1.090 :1@C1'  :1@H1'  11 18 CT H2