Improve atom mapping between chiral centers (#1111)

* New method of mapping atoms bonded to chiral atoms when there are 4
bonded atoms and 2 are known; use the sign of the improper torsion
instead of a delta. Should be more robust to any distortions in
structures.

* Make chiral improper cutoff a class variable

* Separate out case where chiral centers only have one unmapped atom each

* Use map_atoms

* Remove old 1 to 1 chiral map code

* Warn about large deltas when mapping based on chiral sign

* Add function to set chiral improper cut. Use forward declare.

* Add chiralimpcut keyword

* Version 6.29.5. Revision bump for new chiral map code and addition of
chiralimpcut keyword to atommap

* Add chiralimpcut keyword to manual.

* Update manual pdf

doc/DocumentChecksums.txt

@@ -1,3 +1,3 @@
 b37726e7a841f6fc695ecd7fb040ffbf  CpptrajDevelopmentGuide.lyx
-151841b1d3283890c425aaae5a923e7f  cpptraj.lyx
+786c4aec62d60a3f9265dec52b767e02  cpptraj.lyx
 07c4039e732fc2eb1df0c1e0863cb949  CpptrajManual.lyx

doc/cpptraj.lyx

@@ -24177,7 +24177,7 @@ atommap
 \end_layout
 
 \begin_layout LyX-Code
-        [changenames]
+        [changenames] [chiralimpcut <cut>]
 \end_layout
 
 \begin_deeper
@@ -24246,6 +24246,14 @@ changenames If specified,
  change names of mapped atoms in <target> to match those in <reference>.
 \end_layout
 
+\begin_layout Description
+chiralimpcut
+\begin_inset space ~
+\end_inset
+
+<cut> sets the improper dihedral angle cutoff for determining mapping via chirality (default 10 deg.).
+\end_layout
+
 \end_deeper
 \begin_layout Standard
 Attempt to map the atoms of <target> to those of <reference> based on structural similarity.
@@ -24269,6 +24277,34 @@ rmsfit
 \series default
  keyword is useful in cases where the atom mapping will not be complete (e.g.
  two ligands with the same scaffold but different substituents).
+ If 
+\series bold
+changenames
+\series default
+ is specified,
+ in addition to remapping,
+ the target atom names will be changed to match the reference atom names.
+\end_layout
+
+\begin_layout Standard
+Part of the mapping process involves mapping atoms bonded to chiral centers,
+ which in turn involves calculating so-called 
+\begin_inset Quotes eld
+\end_inset
+
+improper
+\begin_inset Quotes erd
+\end_inset
+
+ dihedral angles to determine the orientation of the bonded atoms.
+ By default the code assumes the orientation of the atoms around the chiral centers are fairly close,
+ hence there is a small improper dihedral angle cutoff of 10 degrees.
+ However,
+ this can be increased via the 
+\series bold
+chiralimpcut
+\series default
+ keyword to handle cases where one structure is distorted with respect to the other.
 \end_layout
 
 \begin_layout Standard

src/Action_AtomMap.cpp

@@ -22,13 +22,16 @@ Action_AtomMap::Action_AtomMap() :
 void Action_AtomMap::Help() const {
   mprintf("\t<target> <reference> [mapout <filename>] [maponly]\n"
           "\t[rmsfit [ rmsout <rmsout> ]] [mode {all | byres}]\n"
-          "\t[changenames]\n"
+          "\t[changenames] [chiralimpcut <cut>]\n"
           "  Attempt to create a map from atoms in <target> to atoms in <reference> solely\n"
           "  based on how they are bonded; remap <target> so it matches <reference>.\n"
           "  If 'rmsfit' is specified, calculate the RMSD of remapped target to reference.\n"
           "  If 'changenames' is specified change target atom names to match reference.\n"
+          "  'chiralimpcut' sets the improper dihedral angle cutoff for determining\n"
+          "    mapping via chirality (default 10 deg.)\n"
           "  Modes: 'all'   : try to map all atoms at once (default).\n"
-          "         'byres' : map residue by residue (assumes 1 to 1 residue correspondence).\n");
+          "         'byres' : map residue by residue (assumes 1 to 1 residue correspondence).\n"
+         );
 }
 
 // DESTRUCTOR
@@ -49,6 +52,7 @@ Action::RetType Action_AtomMap::Init(ArgList& actionArgs, ActionInit& init, int
   // Prevent non-master from writing.
   if (!init.TrajComm().Master()) outputfile = 0;
 # endif
+  double chiral_improper_cut_deg = actionArgs.getKeyDouble("chiralimpcut", 10.0);
   maponly_ = actionArgs.hasKey("maponly");
   rmsfit_ = actionArgs.hasKey("rmsfit");
   renameAtoms_ = actionArgs.hasKey("changenames");
@@ -96,6 +100,8 @@ Action::RetType Action_AtomMap::Init(ArgList& actionArgs, ActionInit& init, int
     mprintf("\tMap will only be written, not used to remap input trajectories.\n");
   else
     mprintf("\tAtoms in input trajectories matching target will be remapped.\n");
+  mprintf("\tCutoff for determining mapping via chiral improper dihedral angle: %g deg.\n",
+          chiral_improper_cut_deg);
   if (!maponly_) {
     if (rmsfit_) {
       mprintf("\tWill RMS-fit mapped atoms in tgt to reference.\n");
@@ -119,6 +125,7 @@ Action::RetType Action_AtomMap::Init(ArgList& actionArgs, ActionInit& init, int
   // ---------------------------------------------
   int err = 0;
   StructureMapper Mapper;
+  Mapper.SetChiralImproperCut(chiral_improper_cut_deg);
   switch (mode_) {
     case BY_RES: err = Mapper.CreateMapByResidue( RefFrame_, TgtFrame_, debug_ ); break;
     case ALL   : err = Mapper.CreateMap( RefFrame_, TgtFrame_, debug_ ); break;

src/StructureMapper.cpp

@@ -2,6 +2,7 @@
 #include <list>
 #include "StructureMapper.h"
 #include "CpptrajStdio.h"
+#include "DataSet_Coords_REF.h"
 #include "TorsionRoutines.h"
 #include "Structure/Chirality.h"
 
@@ -12,7 +13,8 @@ StructureMapper::StructureMapper() :
   refTop_(0),
   tgtTop_(0),
   refFrame_(0),
-  tgtFrame_(0)
+  tgtFrame_(0),
+  chiral_improper_cut_(10.0*Constants::DEGRAD)
 {}
 
 /** DESTRUCTOR */
@@ -35,6 +37,16 @@ void StructureMapper::clearPseudoTopFrame() {
   tgtFrame_ = 0;
 }
 
+/** Set chiral improper angle cutoff (in degrees, convert to radians). */
+void StructureMapper::SetChiralImproperCut(double cutInDeg) {
+  if (cutInDeg < 0) {
+    mprintf("Warning: Negative value given for chiral improper angle cutoff.\n"
+            "Warning: Setting to 0.0.\n");
+    chiral_improper_cut_ = 0;
+  } else
+    chiral_improper_cut_ = cutInDeg * Constants::DEGRAD;
+}
+
 // StructureMapper::mapBondsToUnique()
 /** For each atom R in reference already mapped to unique atom T in 
   * target, try to match up non-mapped reference atoms r bonded to R to 
@@ -401,6 +413,48 @@ const
   return numMappedAtoms;
 }
 
+/** Map given ref atom to target atom. */
+void StructureMapper::map_atoms(int iR, int iT, AtomMap& Ref, AtomMap& Tgt, int& numMappedAtoms,
+                                const char* desc)
+{
+  if (debug_>0)
+    mprintf("    Mapping tgt atom %i:%s to ref atom %i:%s based on %s.\n",
+            iT+1, Tgt[iT].c_str(), iR+1, Ref[iR].c_str(), desc );
+  AMap_[ iR ] = iT;
+  ++numMappedAtoms;
+  // Once an atom has been mapped set its mapped flag
+  Ref[iR].SetMapped();
+  Tgt[iT].SetMapped();
+}
+
+/** \return true if the absolute value of the given value is less than
+  *         the chiral improper cutoff.
+  */
+bool StructureMapper::chiralImproperCutSatisfied(double deltaIn) const {
+  double delta;
+  if (deltaIn < 0.0)
+    delta = -deltaIn;
+  else
+    delta = deltaIn;
+  return (delta < chiral_improper_cut_);
+}
+
+/** Warning for when atoms will be mapped but delta seems large. */
+void StructureMapper::check_large_delta(int iR, int iT, AtomMap const& Ref, AtomMap const& Tgt,
+                                        int ratom, int tatom, double delta)
+const
+{
+  if (!chiralImproperCutSatisfied(delta))
+    mprintf("Warning: Ref %i:%s and Tgt %i:%s bonded to chiral centers\n"
+            "Warning: %i:%s | %i:%s have the same sign but a large delta\n"
+            "Warning: (%.4f deg). This can indicate structural problems in either\n"
+            "Warning: the target or reference. Mapping atoms, but it is recommended\n"
+            "Warning: the structures be visually inspected for problems.\n",
+            iR+1, Ref[iR].c_str(), iT+1, Tgt[iT].c_str(),
+            ratom+1, Ref[ratom].c_str(), tatom+1, Tgt[tatom].c_str(),
+            delta*Constants::RADDEG);
+}
+
 // StructureMapper::mapChiral()
 /** Given two atommaps and a map relating the two, find chiral centers for
   * which at least 3 of the atoms have been mapped. Assign the remaining
@@ -532,42 +586,79 @@ int StructureMapper::mapChiral(AtomMap& Ref, AtomMap& Tgt) {
       if (debug_>1) mprintf("    Tgt Improper %i [%3i,%3i,%3i,%3i]= %lf\n",i,
                            uT[0]+1, uT[1]+1, uT[2]+1, nT[i]+1, dT[i]*Constants::RADDEG);
     }
-    // Match impropers to each other using a cutoff. Note that all torsions
-    // are in radians.
-    // NOTE: 10.0 degrees seems reasonable? Also there is currently no 
-    //       check for repeated deltas.
-    for (int i=0; i<notunique_r; i++) {
-      for (int j=0; j<notunique_t; j++) {
-        double delta = dR[i] - dT[j];
-        if (delta<0.0) delta=-delta;
-        if (delta<0.17453292519943295769236907684886) {
-          if (debug_>0)
-            mprintf("    Mapping tgt atom %i:%s to ref atom %i:%s based on chirality.\n",
-                    nT[j]+1, Tgt[nT[j]].c_str(), nR[i]+1, Ref[nR[i]].c_str() );
-          AMap_[ nR[i] ] = nT[j];
-          ++numMappedAtoms;
-          // Once an atom has been mapped set its mapped flag
-          Ref[nR[i]].SetMapped();
-          Tgt[nT[j]].SetMapped();
-        } else if (notunique_r == 1 && notunique_t == 1) {
-          // This is the only non-mapped atom of the chiral center but for
-          // some reason the improper dihedral doesnt match. Map it but warn
-          // the user.
-          mprintf("Warning: Ref %i:%s and Tgt %i:%s are the only unmapped atoms of chiral\n"
-                  "Warning: centers %i:%s | %i:%s, but the improper dihedral angles do not\n"
-                  "Warning: match (%.4f rad != %.4f rad). This can indicate structural problems\n"
-                  "Warning: in either the target or reference. Mapping atoms, but it is\n"
-                  "Warning: recommended the structures be visually inspected for problems.\n",
-                  nR[i]+1, Ref[nR[i]].c_str(), nT[j]+1, Tgt[nT[j]].c_str(),
-                  ratom+1, Ref[ratom].c_str(), tatom+1, Tgt[tatom].c_str(),
-                  dR[i], dT[j]);
-          AMap_[ nR[i] ] = nT[j];
-          ++numMappedAtoms;
-          Ref[nR[i]].SetMapped();
-          Tgt[nT[j]].SetMapped();
-        }
+    // Try different mapping methods
+    bool successful_map = false;
+    if (notunique_r == 1 && notunique_t == 1) {
+      // Only 1 non-mapped atom of each chiral center.
+      // If for some reason the improper dihedral doesnt match, map it
+      // but warn the user.
+      if (!chiralImproperCutSatisfied(dR[0] - dT[0])) {
+        mprintf("Warning: Ref %i:%s and Tgt %i:%s are the only unmapped atoms of chiral\n"
+                "Warning: centers %i:%s | %i:%s, but the improper dihedral angles do not\n"
+                "Warning: match (%.4f deg != %.4f deg). This can indicate structural problems\n"
+                "Warning: in either the target or reference. Mapping atoms, but it is\n"
+                "Warning: recommended the structures be visually inspected for problems.\n",
+                nR[0]+1, Ref[nR[0]].c_str(), nT[0]+1, Tgt[nT[0]].c_str(),
+                ratom+1, Ref[ratom].c_str(), tatom+1, Tgt[tatom].c_str(),
+                dR[0]*Constants::RADDEG, dT[0]*Constants::RADDEG);
+      }
+      map_atoms( nR[0], nT[0], Ref, Tgt, numMappedAtoms, "only unmapped atom" );
+      successful_map = true;
+    } else if (notunique_r == 2 && notunique_t == 2) {
+      // Chiral center atom, 2 mapped atoms, 2 unmapped atoms.
+      // Expect one atom to have positive dihedral, one to have negative.
+      // Map positive to positive and negative to negative.
+      int ref_0, ref_1, tgt_0, tgt_1;
+      if (dR[0] < 0)
+        ref_0 = -1;
+      else
+        ref_0 = 1;
+      if (dR[1] < 0)
+        ref_1 = -1;
+      else
+        tgt_1 = 1;
+      if (dT[0] < 0)
+        tgt_0 = -1;
+      else
+        tgt_0 = 1;
+      if (dT[1] < 0)
+        tgt_1 = -1;
+      else
+        tgt_1 = 1;
+      if ( (ref_0 + ref_1 == 0) && (tgt_0 + tgt_1 == 0) ) {
+        // Each has a positive and negative angle.
+        if (ref_0 == tgt_0) {
+          check_large_delta(nR[0], nT[0], Ref, Tgt, ratom, tatom, dR[0]-dT[0]);
+          check_large_delta(nR[1], nT[1], Ref, Tgt, ratom, tatom, dR[1]-dT[1]);
+
+          map_atoms( nR[0], nT[0], Ref, Tgt, numMappedAtoms, "chirality sign" );
+          map_atoms( nR[1], nT[1], Ref, Tgt, numMappedAtoms, "chirality sign" );
+          successful_map = true;
+        } else if (ref_0 == tgt_1) {
+          check_large_delta(nR[0], nT[1], Ref, Tgt, ratom, tatom, dR[0]-dT[1]);
+          check_large_delta(nR[1], nT[0], Ref, Tgt, ratom, tatom, dR[1]-dT[0]);
+
+          map_atoms( nR[0], nT[1], Ref, Tgt, numMappedAtoms, "chirality sign" );
+          map_atoms( nR[1], nT[0], Ref, Tgt, numMappedAtoms, "chirality sign" );
+          successful_map = true;
+        } 
       }
     }
+    if (!successful_map) {
+      // Chiral map general (previous default) method.
+      // Match impropers to each other using a cutoff. Note that all torsions
+      // are in radians.
+      // NOTE: 10.0 degrees seems reasonable? Also there is currently no 
+      //       check for repeated deltas.
+      for (int i=0; i < notunique_r; i++) {
+        for (int j=0; j < notunique_t; j++) {
+          double delta = dR[i] - dT[j];
+          if (chiralImproperCutSatisfied(delta)) {
+            map_atoms( nR[i], nT[j], Ref, Tgt, numMappedAtoms, "chirality angle" );
+          }
+        }
+      }
+    } // END chiral map general method 
     // Check if ref atom or tgt atom is now completely mapped
     Ref.MarkAtomComplete(ratom,false);
     Tgt.MarkAtomComplete(tatom,false);

src/StructureMapper.h

@@ -1,7 +1,7 @@
 #ifndef INC_STRUCTUREMAPPER_H
 #define INC_STRUCTUREMAPPER_H
 #include "AtomMap.h"
-#include "DataSet_Coords_REF.h"
+class DataSet_Coords_REF;
 /// Attempt to create a one-to-one mapping between atoms in two structures.
 class StructureMapper {
   public:
@@ -11,6 +11,9 @@ class StructureMapper {
     /// DESTRUCTOR
     ~StructureMapper();
 
+    /// Set the improper angle cutoff (in deg.) for mapping via chirality.
+    void SetChiralImproperCut(double);
+
     int CreateMap(DataSet_Coords_REF*, DataSet_Coords_REF*, int);
     int CreateMapByResidue(DataSet_Coords_REF*, DataSet_Coords_REF*, int);
     int operator[](int idx)         const { return AMap_[idx];   }
@@ -31,6 +34,10 @@ class StructureMapper {
     /// Map unmapped atoms bonded to mapped chiral centers via priority
     int mapChiral_viaPriority(MapType&, AtomMap&, AtomMap&, int, int);
 
+    void map_atoms(int, int, AtomMap&, AtomMap&, int&, const char*);
+    bool chiralImproperCutSatisfied(double) const;
+    void check_large_delta(int, int, AtomMap const&, AtomMap const&, int, int, double) const;
+
     int mapChiral(AtomMap&, AtomMap&);
     int mapByIndex(AtomMap&, AtomMap&);
     int mapUniqueRefToTgt(AtomMap&, AtomMap&, int);
@@ -50,5 +57,6 @@ class StructureMapper {
     Topology* tgtTop_; ///< Pseudo topology for current tgt atoms being mapped
     Frame* refFrame_;  ///< Pseudo frame for current ref atoms being mapped to
     Frame* tgtFrame_;  ///< Pseudo frame for current tgt atoms being mapped
+    double chiral_improper_cut_; ///< Cutoff for mapping chiral based on improper angle (radians)
 };
 #endif

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