diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf index b6c344a580..432281bdc6 100644 Binary files a/doc/CpptrajManual.pdf and b/doc/CpptrajManual.pdf differ diff --git a/doc/DocumentChecksums.txt b/doc/DocumentChecksums.txt index 5bc24381c3..7cd22a8665 100644 --- a/doc/DocumentChecksums.txt +++ b/doc/DocumentChecksums.txt @@ -1,3 +1,3 @@ b37726e7a841f6fc695ecd7fb040ffbf CpptrajDevelopmentGuide.lyx -151841b1d3283890c425aaae5a923e7f cpptraj.lyx +786c4aec62d60a3f9265dec52b767e02 cpptraj.lyx 07c4039e732fc2eb1df0c1e0863cb949 CpptrajManual.lyx diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 60732e1e34..f5e0287a3f 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -24177,7 +24177,7 @@ atommap \end_layout \begin_layout LyX-Code - [changenames] + [changenames] [chiralimpcut ] \end_layout \begin_deeper @@ -24246,6 +24246,14 @@ changenames If specified, change names of mapped atoms in to match those in . \end_layout +\begin_layout Description +chiralimpcut +\begin_inset space ~ +\end_inset + + 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 to those of 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 diff --git a/src/Action_AtomMap.cpp b/src/Action_AtomMap.cpp index e1441c03b2..ea747d9756 100644 --- a/src/Action_AtomMap.cpp +++ b/src/Action_AtomMap.cpp @@ -22,13 +22,16 @@ Action_AtomMap::Action_AtomMap() : void Action_AtomMap::Help() const { mprintf("\t [mapout ] [maponly]\n" "\t[rmsfit [ rmsout ]] [mode {all | byres}]\n" - "\t[changenames]\n" + "\t[changenames] [chiralimpcut ]\n" " Attempt to create a map from atoms in to atoms in solely\n" " based on how they are bonded; remap so it matches .\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; diff --git a/src/StructureMapper.cpp b/src/StructureMapper.cpp index a502f173c4..a03553a470 100644 --- a/src/StructureMapper.cpp +++ b/src/StructureMapper.cpp @@ -2,6 +2,7 @@ #include #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; i0) - 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); diff --git a/src/StructureMapper.h b/src/StructureMapper.h index 9315975907..bf4cc24799 100644 --- a/src/StructureMapper.h +++ b/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 diff --git a/src/Version.h b/src/Version.h index 7957c39291..355123d3f4 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes 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