Add 'rescut' and 'bondoffset' keywords to 'prepareforleap'

doc/DocumentChecksums.txt

@@ -0,0 +1,3 @@
+b37726e7a841f6fc695ecd7fb040ffbf  CpptrajDevelopmentGuide.lyx
+c9b19c2bade122b63363f364f86a0a69  cpptraj.lyx
+07c4039e732fc2eb1df0c1e0863cb949  CpptrajManual.lyx

doc/GeneratePDFs.sh

@@ -0,0 +1,61 @@
+#!/bin/bash
+
+LYX=`which lyx`
+if [ -z "$LYX" ] ; then
+  echo "No lyx present."
+  exit 1
+fi
+MD5SUM=`which md5sum`
+if [ -z "$MD5SUM" ] ; then
+  echo "No md5sum present."
+  exit 1
+fi
+
+GEN_MANUAL=0
+C0=`grep CpptrajManual.lyx DocumentChecksums.txt`
+M0=`md5sum CpptrajManual.lyx`
+echo $C0
+echo $M0
+if [ "$C0" != "$M0" ] ; then
+  GEN_MANUAL=1
+fi
+C1=`grep cpptraj.lyx DocumentChecksums.txt`
+M1=`md5sum cpptraj.lyx`
+echo $C1
+echo $M1
+if [ "$C1" != "$M1" ] ; then
+  GEN_MANUAL=1
+fi
+if [ $GEN_MANUAL -eq 1 ] ; then
+  lyx -batch --export pdf2 CpptrajManual.lyx
+  if [ $? -ne 0 ] ; then
+    echo "Generation of manual failed."
+    exit 1
+  fi
+  ls -l CpptrajManual.pdf
+fi
+
+GEN_DEV=0
+C2=`grep CpptrajDevelopmentGuide.lyx DocumentChecksums.txt`
+M2=`md5sum CpptrajDevelopmentGuide.lyx`
+echo $C2
+echo $M2
+if [ "$C2" != "$M2" ] ; then
+  GEN_DEV=1
+fi
+if [ $GEN_DEV -eq 1 ] ; then
+  lyx -batch --export pdf2 CpptrajDevelopmentGuide.lyx
+  if [ $? -ne 0 ] ; then
+    echo "Generation of dev guide failed."
+    exit 1
+  fi
+  ls -l CpptrajDevelopmentGuide.pdf
+fi
+
+# Update checksums if needed
+if [ $GEN_MANUAL -ne 0 -o $GEN_DEV -ne 0 ] ; then
+  echo "Updating document checksums."
+  md5sum *.lyx > DocumentChecksums.txt
+  ls -l DocumentChecksums.txt
+fi
+exit 0

doc/cpptraj.lyx

@@ -8500,6 +8500,10 @@ prepareforleap crdset <coords set> [frame <#>] name <out coords set>
                   	  sugarmask <sugarmask> [noc1search] [nosplitres]
 \end_layout
 
+\begin_layout LyX-Code
+                        [rescut <residue cutoff>] [bondoffset <offset>]
+\end_layout
+
 \begin_layout LyX-Code
                         [resmapfile <file>]
 \end_layout
@@ -8828,11 +8832,42 @@ resmapfile
 \end_layout
 
 \begin_layout Description
-[noc1search] If specified disable search for missing sugar C1 atom bonds.
+[noc1search] If specified,
+ disable search for missing linkages to sugar C1 atom bonds.
+\end_layout
+
+\begin_layout Description
+[nosplitres] If specified,
+ do not attempt to split off functional groups from sugars into separate residues.
+\end_layout
+
+\begin_layout Description
+[rescut
+\begin_inset space ~
+\end_inset
+
+<residue
+\begin_inset space ~
+\end_inset
+
+cutoff>] Initial distance cutoff (default 8 Ang.) for residue center to residue center distance when looking for missing sugar linkages.
 \end_layout
 
 \begin_layout Description
-[nosplitres] If specified do not attempt to split off functional groups from sugars into separate residues.
+[bondoffset
+\begin_inset space ~
+\end_inset
+
+<offset>] Offset (default 0.2 Ang.) to add to 
+\begin_inset Quotes eld
+\end_inset
+
+ideal
+\begin_inset Quotes erd
+\end_inset
+
+ bond distances when looking for missing sugar linkages.
+ Can be increased to accommodate distorted structures.
 \end_layout
 
 \begin_layout Description

src/Exec_PrepareForLeap.cpp

@@ -653,6 +653,7 @@ void Exec_PrepareForLeap::Help() const
           "\t  [cysmask <cysmask>] [disulfidecut <cut>] [newcysname <name>]}]\n"
           "\t[{nosugars |\n"
           "\t  sugarmask <sugarmask> [noc1search] [nosplitres]\n"
+          "\t  [rescut <residue cutoff>] [bondoffset <offset>]\n"
           "\t  [resmapfile <file>]\n"
           "\t  [hasglycam] [determinesugarsby {geom|name}]\n"
           "\t }]\n"
@@ -785,6 +786,8 @@ Exec::RetType Exec_PrepareForLeap::Execute(CpptrajState& State, ArgList& argIn)
   if (prepare_sugars) {
     // Init options
     if (sugarBuilder.InitOptions( argIn.hasKey("hasglycam"),
+                                  argIn.getKeyDouble("rescut", 8.0),
+                                  argIn.getKeyDouble("bondoffset", 0.2),
                                   argIn.GetStringKey("sugarmask"),
                                   argIn.GetStringKey("determinesugarsby", "geometry"),
                                   argIn.GetStringKey("resmapfile") ))

src/Structure/SugarBuilder.cpp

@@ -16,13 +16,17 @@ using namespace Cpptraj::Structure;
 
 /** CONSTRUCTOR */
 SugarBuilder::SugarBuilder(int debugIn) :
+  rescut2_(64.0), // default initial cutoff 8 Ang
+  offset_(0.2),
   hasGlycam_(false),
   useSugarName_(false),
   debug_(debugIn)
 {}
 
 /** Initialize options. */
 int SugarBuilder::InitOptions(bool hasGlycamIn,
+                              double rescutIn,
+                              double offsetIn,
                               std::string const& sugarmaskstrIn,
                               std::string const& determineSugarsBy,
                               std::string const& resmapfile)
@@ -32,6 +36,11 @@ int SugarBuilder::InitOptions(bool hasGlycamIn,
     return 1;
   }
   hasGlycam_ = hasGlycamIn;
+  if (rescutIn <= 0)
+    rescut2_ = 64.0;
+  else
+    rescut2_ = rescutIn * rescutIn;
+  offset_ = offsetIn;
   sugarmaskstr_ = sugarmaskstrIn;
 
   // Check options
@@ -64,6 +73,8 @@ int SugarBuilder::InitOptions(bool hasGlycamIn,
   // Write options to stdout
   if (hasGlycam_)
     mprintf("\tAssuming sugars already have glycam residue names.\n");
+  mprintf("\tInitial residue-residue cutoff for determining linkages: %g Ang.\n", sqrt(rescut2_));
+  mprintf("\tOffset added to bond cutoffs for determining linkages: %g Ang.\n", offset_);
   return 0;
 }
 
@@ -1685,13 +1696,23 @@ int SugarBuilder::IdentifySugar(Sugar& sugarIn, Topology& topIn,
 }
 
 // -----------------------------------------------------------------------------
+/** Cache residue centers. */
+void SugarBuilder::CacheResidueCenters(Topology const& topIn, Frame const& frameIn)
+{
+  residueCenters_.clear();
+  mprintf("\tCaching residue-residue center distances.\n");
+  residueCenters_.reserve( topIn.Nres() );
+  for (int rr = 0; rr < topIn.Nres(); rr++)
+    residueCenters_.push_back( frameIn.VGeometricCenter( topIn.Res(rr).FirstAtom(), topIn.Res(rr).LastAtom() ) );
+}
+
 /** Attempt to find any missing linkages to the anomeric carbon in sugar. */
 int SugarBuilder::FindSugarC1Linkages(int rnum1, int c_beg,
                                       Topology& topIn, Frame const& frameIn,
                                       NameType const& solventResName)
 const
 {
-  Residue const& res1 = topIn.SetRes(rnum1);
+  //Residue const& res1 = topIn.SetRes(rnum1);
   // If the anomeric atom is already bonded to another residue, skip this.
   for (Atom::bond_iterator bat = topIn[c_beg].bondbegin();
                            bat != topIn[c_beg].bondend(); ++bat)
@@ -1705,9 +1726,9 @@ const
   }
   // TODO pass these in
   // residue first atom to residue first atom cutoff^2
-  const double rescut2 = 64.0;
+  //const double rescut2 = 64.0;
   // bond cutoff offset
-  const double offset = 0.2;
+  //const double offset = 0.2;
   // index of atom to be bonded to c_beg
   int closest_at = -1;
   // distance^2 of atom to be bonded to c_beg
@@ -1716,29 +1737,38 @@ const
   Atom::AtomicElementType a1Elt = topIn[c_beg].Element(); // Should always be C
   if (debug_ > 0)
     mprintf("DEBUG: Anomeric ring carbon: %s\n", topIn.ResNameNumAtomNameNum(c_beg).c_str());
+  // Calculate residue centers
+  if (residueCenters_.empty()) {
+    mprinterr("Internal Error: Residue center distances were not cached.\n");
+    return 1;
+  }
   // Loop over other residues
   for (int rnum2 = 0; rnum2 < topIn.Nres(); rnum2++)
   {
     if (rnum2 != rnum1) {
       Residue const& res2 = topIn.Res(rnum2);
       // Ignore solvent residues
       if (res2.Name() != solventResName) {
-        int at1 = res1.FirstAtom();
+        //int at1 = res1.FirstAtom();
         int at2 = res2.FirstAtom();
         // Initial residue-residue distance based on first atoms in each residue
-        double dist2_1 = DIST2_NoImage( frameIn.XYZ(at1), frameIn.XYZ(at2) );
-        if (dist2_1 < rescut2) {
+        //double dist2_1 = DIST2_NoImage( frameIn.XYZ(at1), frameIn.XYZ(at2) );
+        // Initial residue-residue center distance
+        double dist2_1 = DIST2_NoImage( residueCenters_[rnum1].Dptr(), residueCenters_[rnum2].Dptr() );
+        if (dist2_1 < rescut2_) {
           if (debug_ > 1)
             mprintf("DEBUG: Residue %s to %s = %f\n",
                     topIn.TruncResNameOnumId(rnum1).c_str(), topIn.TruncResNameOnumId(rnum2).c_str(),
                     sqrt(dist2_1));
           // Do the rest of the atoms in res2 to the anomeric carbon
           for (; at2 != res2.LastAtom(); ++at2)
           {
-            if (!topIn[c_beg].IsBondedTo(at2)) {
+            if (topIn[at2].Element() != Atom::HYDROGEN && !topIn[c_beg].IsBondedTo(at2)) {
               double D2 = DIST2_NoImage( frameIn.XYZ(c_beg), frameIn.XYZ(at2) );
               Atom::AtomicElementType a2Elt = topIn[at2].Element();
-              double cutoff2 = Atom::GetBondLength(a1Elt, a2Elt) + offset;
+              double cutoff2 = Atom::GetBondLength(a1Elt, a2Elt) + offset_;
+              //    mprintf("DEBUG: Atom %s to %s = %f cut %f\n",
+              //            topIn.AtomMaskName(c_beg).c_str(), topIn.AtomMaskName(at2).c_str(), sqrt(D2), cutoff2);
               cutoff2 *= cutoff2;
               if (D2 < cutoff2) {
                 if (debug_ > 1)
@@ -1798,6 +1828,9 @@ int SugarBuilder::FixSugarsStructure(Topology& topIn, Frame& frameIn,
     mprintf("Warning: No sugar atoms selected by %s\n", sugarMask.MaskString());
     return 0;
   }
+  // Cache residue distances if we will be looking for linkages
+  if (c1bondsearch)
+    CacheResidueCenters(topIn, frameIn);
   //CharMask cmask( sugarMask.ConvertToCharMask(), sugarMask.Nselected() );
   // Get sugar residue numbers.
   Iarray sugarResNums = topIn.ResnumsSelectedBy( sugarMask );
@@ -1828,6 +1861,7 @@ int SugarBuilder::FixSugarsStructure(Topology& topIn, Frame& frameIn,
   }
 
   if (c1bondsearch) {
+
     // Loop over sugar indices to see if anomeric C is missing bonds
     for (std::vector<Sugar>::const_iterator sugar = Sugars_.begin();
                                             sugar != Sugars_.end(); ++sugar)

src/Structure/SugarBuilder.h

@@ -17,8 +17,8 @@ class SugarBuilder {
     typedef std::vector<Sugar> Array;
     /// CONSTRUCTOR - Take debug level
     SugarBuilder(int);
-    /// Init options: hasGlycam, sugar mask str, determineSugarsBy, resmapfile
-    int InitOptions(bool, std::string const&, std::string const&, std::string const&);
+    /// Init options: hasGlycam, residue cutoff, bond offset, sugar mask str, determineSugarsBy, resmapfile
+    int InitOptions(bool, double, double, std::string const&, std::string const&, std::string const&);
     /// \return true if given res name is a recognized PDB sugar
     bool IsRecognizedPdbSugar(NameType const&) const;
     /// ID sugar rings, find missing C1 links, split off functional groups
@@ -83,6 +83,8 @@ class SugarBuilder {
                       Cpptraj::Structure::ResStatArray&,
                       std::set<BondType>&,
                       std::set<BondType>&);
+    /// Cache residue centers for determining linkages
+    void CacheResidueCenters(Topology const&, Frame const&);
     /// Try to find missing linkages to anomeric carbon in sugar.
     int FindSugarC1Linkages(int, int, Topology&, Frame const&, NameType const&) const;
 
@@ -107,6 +109,9 @@ class SugarBuilder {
 
     Array Sugars_;             ///< Array of found sugars
     std::string sugarmaskstr_; ///< Mask string for selecting sugars
+    double rescut2_;           ///< Cutoff for determining if residue centers are close enough to bond
+    double offset_;            ///< Bond cutoff offset for determining if atoms are bonded
+    std::vector<Vec3> residueCenters_; ///< Residue centers for determining bonding.
     bool hasGlycam_;           ///< If true, assume sugars already have glycam names
     bool useSugarName_;        ///< If true, base form/chirality on name instead of geometry
     AtomMap myMap_;            ///< Used to determine unique atoms for chirality

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 "V6.29.1"
+#define CPPTRAJ_INTERNAL_VERSION "V6.29.2"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif