From 3d469f7c686e23cdd792571f9b0eea953949c530 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 10:09:41 -0400 Subject: [PATCH 01/39] Always print a warning when NA residue not recognized. Add some debug info --- src/Action_NAstruct.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index dc3183b817..055eb47d56 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -1503,11 +1503,10 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { NA_Reference::RetType err = refBases_.SetupBaseRef( currentBase, setup.Top(), *resnum, *masterDSL_, dataname_ ); if (err == NA_Reference::NOT_FOUND) { - // Residue not recognized. Print a warning if the user specified this range. - if (!resRange_.Empty()) { - mprintf("Warning: Residue %i:%s not recognized as NA residue.\n", - *resnum+1, setup.Top().Res(*resnum).c_str()); - } + // Residue not recognized. Print a warning. + mprintf("Warning: Residue %i:%s not recognized as a nucleic acid residue.\n", + *resnum+1, setup.Top().Res(*resnum).c_str()); + mprintf("Warning: For non-standard names, use the 'resmap' keyword.\n"); continue; } else if (err == NA_Reference::BASE_ERROR) { mprinterr("Error: Could not set up residue %s for NA structure analysis.\n", @@ -1564,6 +1563,8 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { if (c5neighbor == Bases_[idx].ResNum()) base->SetC5Idx( idx ); if (c3neighbor == Bases_[idx].ResNum()) base->SetC3Idx( idx ); } + if (debug_ > 0) + mprintf("DEBUG: Base res %i C3res=%i C5res=%i\n", base->ResNum()+1, base->C3resIdx(), base->C5resIdx()); } // END loop over bases setting 5' and 3' neighbors // For each base, find 5' terminii, trace them to the 3' terminii, set up strands From 6f0ca0a5fb9b5dff092607fa332f986760675449 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 11:22:27 -0400 Subject: [PATCH 02/39] Start the axesout keyword --- src/Action_NAstruct.cpp | 20 ++++++++++++++++---- src/Action_NAstruct.h | 5 +++++ 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 055eb47d56..445ef4ab5d 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -5,14 +5,15 @@ #include "StringRoutines.h" // integerToString #include "DistRoutines.h" #include "Constants.h" // RADDEG +#include "Trajout_Single.h" #ifdef NASTRUCTDEBUG -#include "PDBfile.h" +#include "PDBfile.h" //FIXME remove #endif #ifdef MPI # include "DataSet_float.h" // internal pointer needed for sync #endif -// CONSTRUCTOR +/** CONSTRUCTOR */ Action_NAstruct::Action_NAstruct() : puckerMethod_(NA_Base::ALTONA), HBdistCut2_(12.25), // Hydrogen Bond distance cutoff^2: 3.5^2 @@ -37,12 +38,23 @@ Action_NAstruct::Action_NAstruct() : ssout_(0), stepout_(0), helixout_(0), - masterDSL_(0) + masterDSL_(0), # ifdef NASTRUCTDEBUG - ,calcparam_(true) + calcparam_(true), # endif + axesOut_(0), + axesParm_(0) {} +/** DESTRUCTOR */ +Action_NAstruct::~Action_NAstruct() { + if (axesOut_ != 0) { + axesOut_->EndTraj(); + delete axesOut_; + } +} + +/** Print help text. */ void Action_NAstruct::Help() const { mprintf("\t[] [resrange ] [sscalc] [naout ]\n" "\t[noheader] [resmap :{A,C,G,T,U} ...] [calcnohb]\n" diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 14ade00b44..ef9a3347ea 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -7,6 +7,7 @@ #include "AxisType.h" #include "Range.h" #include "DataSet_1D.h" +class Trajout_Single; /// Basic Nucleic acid structure analysis. /** Calculate nucleic acid base/base pair structural parameters. * Algorithms for calculation of base/base pair structural parameters @@ -31,6 +32,8 @@ class Action_NAstruct: public Action { public: Action_NAstruct(); + /// DESTRUCTOR - Needed in case axes trajectories are being written + ~Action_NAstruct(); DispatchObject* Alloc() const { return (DispatchObject*)new Action_NAstruct(); } void Help() const; private: @@ -199,5 +202,7 @@ class Action_NAstruct: public Action { // DEBUG - used to trigger AxisPDBwriter for first call of calculateParameters bool calcparam_; # endif + Trajout_Single* axesOut_; ///< Output trajectory for base axes + Topology* axesParm_; ///< Pseudo-topology for base axes }; #endif From 980954d606846cde3b8b6d559a9a1aad4a43c7f4 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 12:31:04 -0400 Subject: [PATCH 03/39] Start initializing base axes traj --- src/Action_NAstruct.cpp | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 445ef4ab5d..6e9acceae5 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -52,6 +52,9 @@ Action_NAstruct::~Action_NAstruct() { axesOut_->EndTraj(); delete axesOut_; } + if (axesParm_ != 0) { + delete axesParm_; + } } /** Print help text. */ @@ -137,6 +140,32 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int grooveCalcType_ = PP_OO; if (actionArgs.hasKey("altona")) puckerMethod_=NA_Base::ALTONA; else if (actionArgs.hasKey("cremer")) puckerMethod_=NA_Base::CREMER; + // See if we want axes pseudo-trajectories + std::string axesout = actionArgs.GetStringKey("axesout"); + if (!axesout.empty()) { + axesOut_ = new Trajout_Single(); + axesOut_->SetDebug( debug_ ); +# ifdef MPI + axesOut_->SetTrajComm( trajComm_ ); +# endif + std::string axesoutargStr; + std::string axesoutarg = actionArgs.GetStringKey("axesoutarg"); + while (!axesoutarg.empty()) { + axesoutargStr.append(" " + axesoutarg); + axesoutarg = actionArgs.GetStringKey("axesoutarg"); + } + ArgList axesOutArglist( axesoutargStr ); + if (axesOut_->InitEnsembleTrajWrite( axesout, axesOutArglist, init.DSL(), + TrajectoryFile::UNKNOWN_TRAJ, init.DSL().EnsembleNum()) ) + { + mprinterr("Error: Could not init base axes trajectory '%s'\n", axesout.c_str()); + return Action::ERR; + } + std::string axesparmout = actionArgs.GetStringKey("axesparmout"); + axesParm_ = new Topology(); + axesParm_->SetDebug( debug_ ); + axesParm_->SetParmName( "BaseAxes", axesparmout ); + } // Get residue range resRange_.SetRange(actionArgs.GetStringKey("resrange")); if (!resRange_.Empty()) From b2a7fc06214f82e2c8a717bda456a693c0f1d2f8 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 12:45:21 -0400 Subject: [PATCH 04/39] Put pseudo traj init into separate function --- src/Action_NAstruct.cpp | 52 ++++++++++++++++++++++++++++++++++++++++- src/Action_NAstruct.h | 5 ++++ 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 6e9acceae5..6534ebe54d 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -57,6 +57,46 @@ Action_NAstruct::~Action_NAstruct() { } } +/** Set up axes pseudo trajectory. */ +int Action_NAstruct::setup_axes_pseudoTraj(const char* description, + const char* trajKeyword, + const char* argKeyword, + const char* parmKeyword, + const char* topName, + DataSetList const& DSL, + ArgList& actionArgs, + Trajout_Single** outputTraj, + Topology** outputParm) +const +{ + std::string axesout = actionArgs.GetStringKey(trajKeyword); + if (!axesout.empty()) { + *outputTraj = new Trajout_Single(); + (*outputTraj)->SetDebug( debug_ ); +# ifdef MPI + (*outputTraj)->SetTrajComm( trajComm_ ); +# endif + std::string axesoutargStr; + std::string axesoutarg = actionArgs.GetStringKey(argKeyword); + while (!axesoutarg.empty()) { + axesoutargStr.append(" " + axesoutarg); + axesoutarg = actionArgs.GetStringKey(argKeyword); + } + ArgList axesOutArglist( axesoutargStr ); + if ((*outputTraj)->InitEnsembleTrajWrite( axesout, axesOutArglist, DSL, + TrajectoryFile::UNKNOWN_TRAJ, DSL.EnsembleNum()) ) + { + mprinterr("Error: Could not init %s trajectory '%s'\n", description, axesout.c_str()); + return 1; + } + std::string axesparmout = actionArgs.GetStringKey(parmKeyword); + *outputParm = new Topology(); + (*outputParm)->SetDebug( debug_ ); + (*outputParm)->SetParmName( topName, axesparmout ); + } + return 0; +} + /** Print help text. */ void Action_NAstruct::Help() const { mprintf("\t[] [resrange ] [sscalc] [naout ]\n" @@ -141,6 +181,11 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int if (actionArgs.hasKey("altona")) puckerMethod_=NA_Base::ALTONA; else if (actionArgs.hasKey("cremer")) puckerMethod_=NA_Base::CREMER; // See if we want axes pseudo-trajectories + if (setup_axes_pseudoTraj("base axes", "axesout", "axesoutarg", "axesparmout", + "BaseAxes", init.DSL(), actionArgs, + &axesOut_, &axesParm_)) + return Action::ERR; +/* std::string axesout = actionArgs.GetStringKey("axesout"); if (!axesout.empty()) { axesOut_ = new Trajout_Single(); @@ -165,7 +210,7 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int axesParm_ = new Topology(); axesParm_->SetDebug( debug_ ); axesParm_->SetParmName( "BaseAxes", axesparmout ); - } + }*/ // Get residue range resRange_.SetRange(actionArgs.GetStringKey("resrange")); if (!resRange_.Empty()) @@ -1693,6 +1738,11 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { } // END loop over bases in strand } // END loop over strands } // END if sscalc + + // Set up axes pseudo-topologies + if (axesParm_ != 0) { + // 1 pseudo bond type + } return Action::OK; } diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index ef9a3347ea..db15bf37cc 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -128,6 +128,11 @@ class Action_NAstruct: public Action { typedef std::vector StrandArray; ///< Hold indices into Bases_ for strand beg/end typedef std::map RefMapType; ///< Map custom res names to target types // ----- Functions --------------------------- + /// Set up an axes pseudo-trajectory + int setup_axes_pseudoTraj(const char*, const char*, const char*, + const char*, const char*, + DataSetList const&, ArgList&, + Trajout_Single**, Topology**) const; /// Recursively travel to 3' terminal base static int follow_base_to_3prime(Barray&, unsigned int, std::vector&, int); /// Recursively travel sugar-phosphate backbone to find the next residue in a strand. From b4b894d145eecd1ae71e706c3873d7c301bfaa4f Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 12:47:03 -0400 Subject: [PATCH 05/39] Rename function --- src/Action_NAstruct.cpp | 24 ++++++++++++------------ src/Action_NAstruct.h | 10 +++++----- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 6534ebe54d..272db1977b 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -58,15 +58,15 @@ Action_NAstruct::~Action_NAstruct() { } /** Set up axes pseudo trajectory. */ -int Action_NAstruct::setup_axes_pseudoTraj(const char* description, - const char* trajKeyword, - const char* argKeyword, - const char* parmKeyword, - const char* topName, - DataSetList const& DSL, - ArgList& actionArgs, - Trajout_Single** outputTraj, - Topology** outputParm) +int Action_NAstruct::init_axes_pseudoTraj(const char* description, + const char* trajKeyword, + const char* argKeyword, + const char* parmKeyword, + const char* topName, + DataSetList const& DSL, + ArgList& actionArgs, + Trajout_Single** outputTraj, + Topology** outputParm) const { std::string axesout = actionArgs.GetStringKey(trajKeyword); @@ -181,9 +181,9 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int if (actionArgs.hasKey("altona")) puckerMethod_=NA_Base::ALTONA; else if (actionArgs.hasKey("cremer")) puckerMethod_=NA_Base::CREMER; // See if we want axes pseudo-trajectories - if (setup_axes_pseudoTraj("base axes", "axesout", "axesoutarg", "axesparmout", - "BaseAxes", init.DSL(), actionArgs, - &axesOut_, &axesParm_)) + if (init_axes_pseudoTraj("base axes", "axesout", "axesoutarg", "axesparmout", + "BaseAxes", init.DSL(), actionArgs, + &axesOut_, &axesParm_)) return Action::ERR; /* std::string axesout = actionArgs.GetStringKey("axesout"); diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index db15bf37cc..6172156096 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -128,11 +128,11 @@ class Action_NAstruct: public Action { typedef std::vector StrandArray; ///< Hold indices into Bases_ for strand beg/end typedef std::map RefMapType; ///< Map custom res names to target types // ----- Functions --------------------------- - /// Set up an axes pseudo-trajectory - int setup_axes_pseudoTraj(const char*, const char*, const char*, - const char*, const char*, - DataSetList const&, ArgList&, - Trajout_Single**, Topology**) const; + /// Initialize an axes pseudo-trajectory + int init_axes_pseudoTraj(const char*, const char*, const char*, + const char*, const char*, + DataSetList const&, ArgList&, + Trajout_Single**, Topology**) const; /// Recursively travel to 3' terminal base static int follow_base_to_3prime(Barray&, unsigned int, std::vector&, int); /// Recursively travel sugar-phosphate backbone to find the next residue in a strand. From e7e7ff3e2f77d17ab8f26b611bf265fd854aa887 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 13:16:02 -0400 Subject: [PATCH 06/39] Set up topology for axes pseudo trajectory --- src/Action_NAstruct.cpp | 52 ++++++++++++++++++++++++++++++++++++++++- src/Action_NAstruct.h | 2 ++ 2 files changed, 53 insertions(+), 1 deletion(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 272db1977b..b2f536b6be 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -6,6 +6,7 @@ #include "DistRoutines.h" #include "Constants.h" // RADDEG #include "Trajout_Single.h" +#include "ParmFile.h" #ifdef NASTRUCTDEBUG #include "PDBfile.h" //FIXME remove #endif @@ -1741,11 +1742,60 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { // Set up axes pseudo-topologies if (axesParm_ != 0) { - // 1 pseudo bond type + // Create residue information for each base axes + std::vector axesResidues; + axesResidues.reserve( Bases_.size() ); + for (Barray::iterator base = Bases_.begin(); base != Bases_.end(); ++base) { + Residue const& res = setup.Top().Res( base->ResNum() ); + axesResidues.push_back( res ); + } + if (setup_axes_pseudoTraj( *axesParm_, axesResidues )) + return Action::ERR; } return Action::OK; } +/** Set up topology for an axes pseudo-trajectory. */ +int Action_NAstruct::setup_axes_pseudoTraj(Topology& pseudo, + std::vector const& axesResidues) +const +{ + if (pseudo.Natom() > 0) { + mprintf("\tAxes pseudo-topology '%s' is already set up.\n", pseudo.c_str()); + return 0; + } + // 1 pseudo bond type, Rk = 0.0, Req = 1.0 Ang., will be bond parm index 0 + BondArray bonds; + pseudo.AddBondParm( BondParmType(0.0, 1.0) ); + int natom = 0; + for (std::vector::const_iterator res = axesResidues.begin(); + res != axesResidues.end(); + ++res) + { + // Order is origin, x, y, z + pseudo.AddTopAtom(Atom("Orig", 0), *res); + pseudo.AddTopAtom(Atom("X", 0), *res); + pseudo.AddTopAtom(Atom("Y", 0), *res); + pseudo.AddTopAtom(Atom("Z", 0), *res); + // Bond x y and z to origin + pseudo.AddBond(natom, natom+1, 0); + pseudo.AddBond(natom, natom+2, 0); + pseudo.AddBond(natom, natom+3, 0); + natom += 4; + } + + pseudo.CommonSetup(); + if (!pseudo.OriginalFilename().empty()) { + mprintf("\tWriting axes pseudo-topology to '%s'\n", pseudo.OriginalFilename().full()); + ParmFile pfile; + if (pfile.WriteTopology(pseudo, pseudo.OriginalFilename(), ParmFile::UNKNOWN_PARM, debug_)) { + mprinterr("Error: Could not write axes pseudo-topology '%s'\n", pseudo.c_str()); + return 1; + } + } + return 0; +} + // Action_NAstruct::CalculateHbonds() void Action_NAstruct::CalculateHbonds() { for (BPmap::iterator BP = BasePairs_.begin(); BP != BasePairs_.end(); ++BP) diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 6172156096..f2ac8abe3a 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -133,6 +133,8 @@ class Action_NAstruct: public Action { const char*, const char*, DataSetList const&, ArgList&, Trajout_Single**, Topology**) const; + /// Set up an axes pseudo-trajectory + int setup_axes_pseudoTraj(Topology&, std::vector const&) const; /// Recursively travel to 3' terminal base static int follow_base_to_3prime(Barray&, unsigned int, std::vector&, int); /// Recursively travel sugar-phosphate backbone to find the next residue in a strand. From d4e8efd8a6cc832c10f9199c670c656381588029 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 13:21:27 -0400 Subject: [PATCH 07/39] Set up the trajectory as well --- src/Action_NAstruct.cpp | 13 +++++++++++-- src/Action_NAstruct.h | 2 +- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index b2f536b6be..7b41644d20 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -1749,7 +1749,7 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { Residue const& res = setup.Top().Res( base->ResNum() ); axesResidues.push_back( res ); } - if (setup_axes_pseudoTraj( *axesParm_, axesResidues )) + if (setup_axes_pseudoTraj( *axesParm_, *axesOut_, axesResidues, setup.Nframes() )) return Action::ERR; } return Action::OK; @@ -1757,7 +1757,9 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { /** Set up topology for an axes pseudo-trajectory. */ int Action_NAstruct::setup_axes_pseudoTraj(Topology& pseudo, - std::vector const& axesResidues) + Trajout_Single& outtraj, + std::vector const& axesResidues, + int Nframes) const { if (pseudo.Natom() > 0) { @@ -1793,6 +1795,13 @@ const return 1; } } + + if (outtraj.SetupTrajWrite( &pseudo, CoordinateInfo(), Nframes)) { + mprinterr("Error: Could not set up axes output trajectory.\n"); + return 1; + } + mprintf(" "); //TODO this is a kludge; PrintInfo should be a string. + outtraj.PrintInfo(0); return 0; } diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index f2ac8abe3a..e0344aad31 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -134,7 +134,7 @@ class Action_NAstruct: public Action { DataSetList const&, ArgList&, Trajout_Single**, Topology**) const; /// Set up an axes pseudo-trajectory - int setup_axes_pseudoTraj(Topology&, std::vector const&) const; + int setup_axes_pseudoTraj(Topology&, Trajout_Single&, std::vector const&, int) const; /// Recursively travel to 3' terminal base static int follow_base_to_3prime(Barray&, unsigned int, std::vector&, int); /// Recursively travel sugar-phosphate backbone to find the next residue in a strand. From 99239725a1e3745b993e97055d18effae04bd1b1 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 13:40:33 -0400 Subject: [PATCH 08/39] Write psuedo axes frame --- src/Action_NAstruct.cpp | 26 +++++++++++++++++++++++++- src/Action_NAstruct.h | 4 +++- 2 files changed, 28 insertions(+), 2 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 7b41644d20..745d5d12d4 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -390,6 +390,19 @@ static void WriteAxes(PDBfile& outfile, int resnum, const char* resname, NA_Axis // ----------------------------------------------------------------------------- #endif +/// Add NA_Axis to given Frame +static void axesToFrame(Frame& frame, NA_Axis const& axis) { + // Origin + Vec3 oxyz = axis.Oxyz(); + frame.AddVec3( oxyz ); + // X vector + frame.AddVec3( axis.Rx() + oxyz ); + // Y vector + frame.AddVec3( axis.Ry() + oxyz ); + // Z vector + frame.AddVec3( axis.Rz() + oxyz ); +} + // Action_NAstruct::SetupBaseAxes() /** For each residue in Bases (set up in Setup()), get the corresponding input * coords and fit the reference coords on top of input coords. This sets up @@ -1749,7 +1762,7 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { Residue const& res = setup.Top().Res( base->ResNum() ); axesResidues.push_back( res ); } - if (setup_axes_pseudoTraj( *axesParm_, *axesOut_, axesResidues, setup.Nframes() )) + if (setup_axes_pseudoTraj( *axesParm_, *axesOut_, axesFrame_, axesResidues, setup.Nframes() )) return Action::ERR; } return Action::OK; @@ -1758,6 +1771,7 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { /** Set up topology for an axes pseudo-trajectory. */ int Action_NAstruct::setup_axes_pseudoTraj(Topology& pseudo, Trajout_Single& outtraj, + Frame& frame, std::vector const& axesResidues, int Nframes) const @@ -1802,6 +1816,8 @@ const } mprintf(" "); //TODO this is a kludge; PrintInfo should be a string. outtraj.PrintInfo(0); + + frame = Frame(natom); return 0; } @@ -1868,6 +1884,14 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { // Determine base pair step parameters DetermineStepParameters(frameNum); + + // Output base axes if needed + if (axesOut_ != 0) { + for (std::vector::const_iterator base = Bases_.begin(); + base != Bases_.end(); ++base) + axesToFrame( axesFrame_, base->Axis() ); + } + nframes_++; return Action::OK; } diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index e0344aad31..f7a94acf1b 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -134,7 +134,8 @@ class Action_NAstruct: public Action { DataSetList const&, ArgList&, Trajout_Single**, Topology**) const; /// Set up an axes pseudo-trajectory - int setup_axes_pseudoTraj(Topology&, Trajout_Single&, std::vector const&, int) const; + int setup_axes_pseudoTraj(Topology&, Trajout_Single&, Frame&, + std::vector const&, int) const; /// Recursively travel to 3' terminal base static int follow_base_to_3prime(Barray&, unsigned int, std::vector&, int); /// Recursively travel sugar-phosphate backbone to find the next residue in a strand. @@ -211,5 +212,6 @@ class Action_NAstruct: public Action { # endif Trajout_Single* axesOut_; ///< Output trajectory for base axes Topology* axesParm_; ///< Pseudo-topology for base axes + Frame axesFrame_; }; #endif From ddfbcd98fb0e68640aac1164be208943ec62b205 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 13:42:07 -0400 Subject: [PATCH 09/39] Update help and dependencies --- src/Action_NAstruct.cpp | 1 + src/cpptrajdepend | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 745d5d12d4..b56262c8a9 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -106,6 +106,7 @@ void Action_NAstruct::Help() const { "\t[bpmode {3dna|babcock}]\n" "\t[hbcut ] [origincut ] [altona | cremer]\n" "\t[zcut ] [zanglecut ] [groovecalc {simple | 3dna}]\n" + "\t[axesout [axesoutarg ...] [axesparmout ]]\n" "\t[{ %s | allframes}]\n", DataSetList::RefArgs); mprintf(" Perform nucleic acid structure analysis. Base pairing can be determined\n" " in multiple ways:\n" diff --git a/src/cpptrajdepend b/src/cpptrajdepend index c18848d7cb..abc6e5dcd8 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -54,7 +54,7 @@ Action_Molsurf.o : Action_Molsurf.cpp Action.h ActionState.h Action_Molsurf.h Ar Action_MultiDihedral.o : Action_MultiDihedral.cpp Action.h ActionState.h Action_MultiDihedral.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h DihedralSearch.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TypeNameHolder.h Unit.h Vec3.h Action_MultiPucker.o : Action_MultiPucker.cpp Action.h ActionState.h Action_MultiPucker.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.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 Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Pucker.h Pucker_PuckerMask.h Pucker_PuckerSearch.h Pucker_PuckerToken.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h TypeNameHolder.h Unit.h Vec3.h Action_MultiVector.o : Action_MultiVector.cpp Action.h ActionState.h Action_MultiVector.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Vector.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h -Action_NAstruct.o : Action_NAstruct.cpp Action.h ActionState.h Action_NAstruct.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h AxisType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_float.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h PDBfile.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h +Action_NAstruct.o : Action_NAstruct.cpp Action.h ActionFrameCounter.h ActionState.h Action_NAstruct.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h AxisType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_float.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h PDBfile.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ParmFile.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajectoryFile.h Trajout_Single.h TypeNameHolder.h Unit.h Vec3.h Action_NMRrst.o : Action_NMRrst.cpp Action.h ActionState.h Action_NMRrst.h ArgList.h AssociatedData.h Atom.h AtomMap.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_double.h DataSet_float.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MapAtom.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h ViewRst.h Action_NativeContacts.o : Action_NativeContacts.cpp Action.h ActionState.h Action_NativeContacts.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_MatrixDbl.h DataSet_integer.h Dimension.h DispatchObject.h DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h PDBfile.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Action_OrderParameter.o : Action_OrderParameter.cpp Action.h ActionState.h Action_OrderParameter.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.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 DistRoutines.h FileIO.h FileName.h FileTypes.h Frame.h ImageOption.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OnlineVarT.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h From 3cec5ab7c0e642d7ba09f93b811fd4a02b2e294f Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 13:53:50 -0400 Subject: [PATCH 10/39] Actually write the base axes frames --- src/Action_NAstruct.cpp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index b56262c8a9..db536db979 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -353,6 +353,11 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int mprintf("\tUsing simple groove width calculation (P-P and O-O base pair distances).\n"); else if (grooveCalcType_ == HASSAN_CALLADINE) mprintf("\tUsing groove width calculation of El Hassan & Calladine.\n"); + if (axesOut_ != 0) { + mprintf("\tWriting base axes pseudo trajectory to '%s'\n", axesOut_->Traj().Filename().full()); + if (!axesParm_->OriginalFilename().empty()) + mprintf("\tWriting base axes pseudo topology to '%s'\n", axesParm_->OriginalFilename().full()); + } mprintf("# Citations: Babcock MS; Pednault EPD; Olson WK; \"Nucleic Acid Structure\n" "# Analysis: Mathematics for Local Cartesian and Helical Structure\n" "# Parameters That Are Truly Comparable Between Structures\",\n" @@ -1790,10 +1795,10 @@ const ++res) { // Order is origin, x, y, z - pseudo.AddTopAtom(Atom("Orig", 0), *res); - pseudo.AddTopAtom(Atom("X", 0), *res); - pseudo.AddTopAtom(Atom("Y", 0), *res); - pseudo.AddTopAtom(Atom("Z", 0), *res); + pseudo.AddTopAtom(Atom("Orig", "C"), *res); + pseudo.AddTopAtom(Atom("X", "H"), *res); + pseudo.AddTopAtom(Atom("Y", "H"), *res); + pseudo.AddTopAtom(Atom("Z", "H"), *res); // Bond x y and z to origin pseudo.AddBond(natom, natom+1, 0); pseudo.AddBond(natom, natom+2, 0); @@ -1888,9 +1893,11 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { // Output base axes if needed if (axesOut_ != 0) { + axesFrame_.ClearAtoms(); for (std::vector::const_iterator base = Bases_.begin(); base != Bases_.end(); ++base) axesToFrame( axesFrame_, base->Axis() ); + axesOut_->WriteSingle(frm.TrajoutNum(), axesFrame_); } nframes_++; From 3ba0b756aed29ddf49adb8837e9fdd6b611ee291 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Mon, 10 Jul 2023 15:19:42 -0400 Subject: [PATCH 11/39] Start base pair axes pseudo traj --- src/Action_NAstruct.cpp | 46 +++++++++++++++++------------------------ src/Action_NAstruct.h | 5 ++++- 2 files changed, 23 insertions(+), 28 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index db536db979..3c554d51aa 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -44,7 +44,9 @@ Action_NAstruct::Action_NAstruct() : calcparam_(true), # endif axesOut_(0), - axesParm_(0) + axesParm_(0), + bpAxesOut_(0), + bpAxesParm_(0) {} /** DESTRUCTOR */ @@ -56,6 +58,13 @@ Action_NAstruct::~Action_NAstruct() { if (axesParm_ != 0) { delete axesParm_; } + if (bpAxesOut_ != 0) { + bpAxesOut_->EndTraj(); + delete bpAxesOut_; + } + if (bpAxesParm_ != 0) { + delete bpAxesParm_; + } } /** Set up axes pseudo trajectory. */ @@ -187,32 +196,10 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int "BaseAxes", init.DSL(), actionArgs, &axesOut_, &axesParm_)) return Action::ERR; -/* - std::string axesout = actionArgs.GetStringKey("axesout"); - if (!axesout.empty()) { - axesOut_ = new Trajout_Single(); - axesOut_->SetDebug( debug_ ); -# ifdef MPI - axesOut_->SetTrajComm( trajComm_ ); -# endif - std::string axesoutargStr; - std::string axesoutarg = actionArgs.GetStringKey("axesoutarg"); - while (!axesoutarg.empty()) { - axesoutargStr.append(" " + axesoutarg); - axesoutarg = actionArgs.GetStringKey("axesoutarg"); - } - ArgList axesOutArglist( axesoutargStr ); - if (axesOut_->InitEnsembleTrajWrite( axesout, axesOutArglist, init.DSL(), - TrajectoryFile::UNKNOWN_TRAJ, init.DSL().EnsembleNum()) ) - { - mprinterr("Error: Could not init base axes trajectory '%s'\n", axesout.c_str()); - return Action::ERR; - } - std::string axesparmout = actionArgs.GetStringKey("axesparmout"); - axesParm_ = new Topology(); - axesParm_->SetDebug( debug_ ); - axesParm_->SetParmName( "BaseAxes", axesparmout ); - }*/ + if (init_axes_pseudoTraj("basepair axes", "bpaxesout", "bpaxesoutarg", "bpaxesparmout", + "BasePairAxes", init.DSL(), actionArgs, + &bpAxesOut_, &bpAxesParm_)) + return Action::ERR; // Get residue range resRange_.SetRange(actionArgs.GetStringKey("resrange")); if (!resRange_.Empty()) @@ -358,6 +345,11 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int if (!axesParm_->OriginalFilename().empty()) mprintf("\tWriting base axes pseudo topology to '%s'\n", axesParm_->OriginalFilename().full()); } + if (bpAxesOut_ != 0) { + mprintf("\tWriting base pair axes pseudo trajectory to '%s'\n", bpAxesOut_->Traj().Filename().full()); + if (!bpAxesParm_->OriginalFilename().empty()) + mprintf("\tWriting base pair axes pseudo topology to '%s'\n", bpAxesParm_->OriginalFilename().full()); + } mprintf("# Citations: Babcock MS; Pednault EPD; Olson WK; \"Nucleic Acid Structure\n" "# Analysis: Mathematics for Local Cartesian and Helical Structure\n" "# Parameters That Are Truly Comparable Between Structures\",\n" diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index f7a94acf1b..a000ce89cb 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -212,6 +212,9 @@ class Action_NAstruct: public Action { # endif Trajout_Single* axesOut_; ///< Output trajectory for base axes Topology* axesParm_; ///< Pseudo-topology for base axes - Frame axesFrame_; + Frame axesFrame_; ///< Frame for base axes pseudo traj + Trajout_Single* bpAxesOut_; ///< Output trajectory for base pair axes + Topology* bpAxesParm_; ///< Pseudo-topology for base pair axes + Frame bpAxesFrame_; ///< Frame for base pair axes pseudo traj }; #endif From aa90a534f5d826693233377c10729669c997e759 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 07:41:24 -0400 Subject: [PATCH 12/39] Only set up base axes topology once --- src/Action_NAstruct.cpp | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 3c554d51aa..96dfc1c19f 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -1751,17 +1751,21 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { } // END loop over strands } // END if sscalc - // Set up axes pseudo-topologies + // Set up base axes pseudo-topology if (axesParm_ != 0) { - // Create residue information for each base axes - std::vector axesResidues; - axesResidues.reserve( Bases_.size() ); - for (Barray::iterator base = Bases_.begin(); base != Bases_.end(); ++base) { - Residue const& res = setup.Top().Res( base->ResNum() ); - axesResidues.push_back( res ); + if (axesParm_->Natom() > 0) { + mprintf("\tBase axes pseudo-topology is already set up.\n"); + } else { + // Create residue information for each base axes + std::vector axesResidues; + axesResidues.reserve( Bases_.size() ); + for (Barray::iterator base = Bases_.begin(); base != Bases_.end(); ++base) { + Residue const& res = setup.Top().Res( base->ResNum() ); + axesResidues.push_back( res ); + } + if (setup_axes_pseudoTraj( *axesParm_, *axesOut_, axesFrame_, axesResidues, setup.Nframes() )) + return Action::ERR; } - if (setup_axes_pseudoTraj( *axesParm_, *axesOut_, axesFrame_, axesResidues, setup.Nframes() )) - return Action::ERR; } return Action::OK; } From 03fbb4264666fcf744a1450a46fb1aaf8ef37d34 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 07:44:51 -0400 Subject: [PATCH 13/39] Add check that number of bases has not changed --- src/Action_NAstruct.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 96dfc1c19f..5a947453d3 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -1755,6 +1755,13 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { if (axesParm_ != 0) { if (axesParm_->Natom() > 0) { mprintf("\tBase axes pseudo-topology is already set up.\n"); + // Check that number of bases has not changed + if ((unsigned int)axesParm_->Nres() != Bases_.size()) { + mprinterr("Error: Number of bases has changed from %i to %zu.\n" + "Error: Base axes pseudo-topology is already set up for %i bases.\n", + axesParm_->Nres(), Bases_.size(), axesParm_->Nres()); + return Action::ERR; + } } else { // Create residue information for each base axes std::vector axesResidues; From 6ed7a8a369a6ea5f6cbf117ff918e9e7f49156a5 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 08:05:40 -0400 Subject: [PATCH 14/39] Check for incompatibilities --- src/Action_NAstruct.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 5a947453d3..bb26450fa4 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -228,6 +228,14 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int return Action::ERR; } # endif + // Check for base pair mode incompatibilities + if (findBPmode_ == ALL) { + if (bpAxesOut_ != 0) { + mprinterr("Error: Cannot use 'allframes' mode with 'bpaxesout' since the # of\n" + "Error: base pairs can change each frame.\n"); + return Action::ERR; + } + } // For guess/specify modes, get base pairing type std::string bptype = actionArgs.GetStringKey("bptype"); while (!bptype.empty()) { From 985c465fe3f0f1221dab69f76e80fd3dd1cb2449 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 08:29:25 -0400 Subject: [PATCH 15/39] Do base pair axes out --- src/Action_NAstruct.cpp | 54 +++++++++++++++++++++++++++++++++++++---- src/Action_NAstruct.h | 4 ++- 2 files changed, 52 insertions(+), 6 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index bb26450fa4..05b5131f7f 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -46,7 +46,9 @@ Action_NAstruct::Action_NAstruct() : axesOut_(0), axesParm_(0), bpAxesOut_(0), - bpAxesParm_(0) + bpAxesParm_(0), + setupNframes_(0), + setupTop_(0) {} /** DESTRUCTOR */ @@ -760,6 +762,7 @@ int Action_NAstruct::DetermineBasePairing() { } // END if base to base origin distance < cut } // END base2 loop } // END base1 loop + return 0; } @@ -1760,6 +1763,8 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { } // END if sscalc // Set up base axes pseudo-topology + setupNframes_ = setup.Nframes(); + setupTop_ = setup.TopAddress(); if (axesParm_ != 0) { if (axesParm_->Natom() > 0) { mprintf("\tBase axes pseudo-topology is already set up.\n"); @@ -1778,7 +1783,7 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { Residue const& res = setup.Top().Res( base->ResNum() ); axesResidues.push_back( res ); } - if (setup_axes_pseudoTraj( *axesParm_, *axesOut_, axesFrame_, axesResidues, setup.Nframes() )) + if (setup_axes_pseudoTraj( *axesParm_, *axesOut_, axesFrame_, axesResidues )) return Action::ERR; } } @@ -1789,8 +1794,7 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { int Action_NAstruct::setup_axes_pseudoTraj(Topology& pseudo, Trajout_Single& outtraj, Frame& frame, - std::vector const& axesResidues, - int Nframes) + std::vector const& axesResidues) const { if (pseudo.Natom() > 0) { @@ -1827,7 +1831,7 @@ const } } - if (outtraj.SetupTrajWrite( &pseudo, CoordinateInfo(), Nframes)) { + if (outtraj.SetupTrajWrite( &pseudo, CoordinateInfo(), setupNframes_)) { mprinterr("Error: Could not set up axes output trajectory.\n"); return 1; } @@ -1892,6 +1896,37 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { # endif findBPmode_ = REFERENCE; } + // Set up base pair axes pseudo-topology. This is done inside the action + // because we may not know about base pairing until the first frame. + if (bpAxesParm_ != 0) { + if (bpAxesParm_->Natom() > 0) { + //mprintf("\tBase pair axes pseudo-topology is already set up.\n"); + // Check that number of base pairs has not changed + if ((unsigned int)bpAxesParm_->Nres() != BasePairs_.size()) { + mprinterr("Error: Number of base pairs has changed from %i to %zu.\n" + "Error: Base pair axes pseudo-topology is already set up for %i bases.\n", + bpAxesParm_->Nres(), BasePairs_.size(), bpAxesParm_->Nres()); + return Action::ERR; + } + } else { + // Create residue information for each base pair axes + std::vector bpAxesResidues; + bpAxesResidues.reserve( BasePairs_.size() ); + for (BPmap::const_iterator it = BasePairs_.begin(); it != BasePairs_.end(); ++it) + { + //if (it->second.nhb_ < 1 && skipIfNoHB_) continue; + BPtype const& BP = it->second; + int b1 = BP.base1idx_; + //int b2 = BP.base2idx_; + NA_Base const& base1 = Bases_[b1]; + //NA_Base const& base2 = Bases_[b2]; + Residue const& res = setupTop_->Res( base1.ResNum() ); + bpAxesResidues.push_back( res ); + } + if (setup_axes_pseudoTraj( *bpAxesParm_, *bpAxesOut_, bpAxesFrame_, bpAxesResidues )) + return Action::ERR; + } + } // Determine strand parameters if desired if (sscalc_) DetermineStrandParameters(frameNum); @@ -1910,6 +1945,15 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { axesToFrame( axesFrame_, base->Axis() ); axesOut_->WriteSingle(frm.TrajoutNum(), axesFrame_); } + // Output base pair axes if needed + if (bpAxesOut_ != 0) { + bpAxesFrame_.ClearAtoms(); + for (BPmap::const_iterator it = BasePairs_.begin(); it != BasePairs_.end(); ++it) { + //if (it->second.nhb_ < 1 && skipIfNoHB_) continue; + axesToFrame( bpAxesFrame_, it->second.bpaxis_ ); + } + bpAxesOut_->WriteSingle(frm.TrajoutNum(), bpAxesFrame_); + } nframes_++; return Action::OK; diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index a000ce89cb..9fb8bbb719 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -135,7 +135,7 @@ class Action_NAstruct: public Action { Trajout_Single**, Topology**) const; /// Set up an axes pseudo-trajectory int setup_axes_pseudoTraj(Topology&, Trajout_Single&, Frame&, - std::vector const&, int) const; + std::vector const&) const; /// Recursively travel to 3' terminal base static int follow_base_to_3prime(Barray&, unsigned int, std::vector&, int); /// Recursively travel sugar-phosphate backbone to find the next residue in a strand. @@ -216,5 +216,7 @@ class Action_NAstruct: public Action { Trajout_Single* bpAxesOut_; ///< Output trajectory for base pair axes Topology* bpAxesParm_; ///< Pseudo-topology for base pair axes Frame bpAxesFrame_; ///< Frame for base pair axes pseudo traj + int setupNframes_; ///< Set in Setup(); number of expected frames to write (pseudo-traj) + Topology* setupTop_; ///< Set in Setup(); current topology }; #endif From 9a5bcbaa576598abae5639be8b8942be512b084c Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 08:39:12 -0400 Subject: [PATCH 16/39] Save base pair step axis for output --- src/Action_NAstruct.cpp | 4 ++-- src/Action_NAstruct.h | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 05b5131f7f..17d0efc682 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -1403,8 +1403,8 @@ int Action_NAstruct::DetermineStepParameters(int frameNum) { } StepType& currentStep = entry->second; // Calc step parameters - NA_Axis midFrame; - calculateParameters(BP1.bpaxis_, BP2.bpaxis_, &midFrame, Param); + calculateParameters(BP1.bpaxis_, BP2.bpaxis_, &(currentStep.stepaxis_), Param); + NA_Axis& midFrame = currentStep.stepaxis_; // Calculate zP: difference in step phosphate atoms along the Z axis // of the step middle frame. float Zp = 0.0; diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 9fb8bbb719..0e5f7e0b7c 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -92,6 +92,7 @@ class Action_NAstruct: public Action { }; /// Hold a base pair step. struct StepType { + NA_Axis stepaxis_; ///< Base pair step axis DataSet_1D* shift_; DataSet_1D* slide_; DataSet_1D* rise_; From a806645c9c89ddd084534fee5946bf053510f7c6 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 08:56:12 -0400 Subject: [PATCH 17/39] Add base pair step axes write --- src/Action_NAstruct.cpp | 62 +++++++++++++++++++++++++++++++++++++++++ src/Action_NAstruct.h | 19 +++++++------ 2 files changed, 73 insertions(+), 8 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 17d0efc682..43bb672cde 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -47,6 +47,8 @@ Action_NAstruct::Action_NAstruct() : axesParm_(0), bpAxesOut_(0), bpAxesParm_(0), + stepAxesOut_(0), + stepAxesParm_(0), setupNframes_(0), setupTop_(0) {} @@ -67,6 +69,13 @@ Action_NAstruct::~Action_NAstruct() { if (bpAxesParm_ != 0) { delete bpAxesParm_; } + if (stepAxesOut_ != 0) { + stepAxesOut_->EndTraj(); + delete stepAxesOut_; + } + if (stepAxesParm_ != 0) { + delete stepAxesParm_; + } } /** Set up axes pseudo trajectory. */ @@ -118,6 +127,8 @@ void Action_NAstruct::Help() const { "\t[hbcut ] [origincut ] [altona | cremer]\n" "\t[zcut ] [zanglecut ] [groovecalc {simple | 3dna}]\n" "\t[axesout [axesoutarg ...] [axesparmout ]]\n" + "\t[bpaxesout [bpaxesoutarg ...] [bpaxesparmout ]]\n" + "\t[stepaxesout [stepaxesoutarg ...] [stepaxesparmout ]]\n" "\t[{ %s | allframes}]\n", DataSetList::RefArgs); mprintf(" Perform nucleic acid structure analysis. Base pairing can be determined\n" " in multiple ways:\n" @@ -202,6 +213,10 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int "BasePairAxes", init.DSL(), actionArgs, &bpAxesOut_, &bpAxesParm_)) return Action::ERR; + if (init_axes_pseudoTraj("step axes", "stepaxesout", "stepaxesoutarg", "stepaxesparmout", + "StepAxes", init.DSL(), actionArgs, + &stepAxesOut_, &stepAxesParm_)) + return Action::ERR; // Get residue range resRange_.SetRange(actionArgs.GetStringKey("resrange")); if (!resRange_.Empty()) @@ -237,6 +252,11 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int "Error: base pairs can change each frame.\n"); return Action::ERR; } + if (stepAxesOut_ != 0) { + mprinterr("Error: Cannot use 'allframes' mode with 'stepaxesout' since the # of\n" + "Error: base pair steps can change each frame.\n"); + return Action::ERR; + } } // For guess/specify modes, get base pairing type std::string bptype = actionArgs.GetStringKey("bptype"); @@ -360,6 +380,11 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int if (!bpAxesParm_->OriginalFilename().empty()) mprintf("\tWriting base pair axes pseudo topology to '%s'\n", bpAxesParm_->OriginalFilename().full()); } + if (stepAxesOut_ != 0) { + mprintf("\tWriting base pair step axes pseudo trajectory to '%s'\n", stepAxesOut_->Traj().Filename().full()); + if (!stepAxesParm_->OriginalFilename().empty()) + mprintf("\tWriting base pair step axes pseudo topology to '%s'\n", stepAxesParm_->OriginalFilename().full()); + } mprintf("# Citations: Babcock MS; Pednault EPD; Olson WK; \"Nucleic Acid Structure\n" "# Analysis: Mathematics for Local Cartesian and Helical Structure\n" "# Parameters That Are Truly Comparable Between Structures\",\n" @@ -1896,6 +1921,7 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { # endif findBPmode_ = REFERENCE; } + // Set up base pair axes pseudo-topology. This is done inside the action // because we may not know about base pairing until the first frame. if (bpAxesParm_ != 0) { @@ -1927,6 +1953,7 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { return Action::ERR; } } + // Determine strand parameters if desired if (sscalc_) DetermineStrandParameters(frameNum); @@ -1937,6 +1964,34 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { // Determine base pair step parameters DetermineStepParameters(frameNum); + // Set up base pair step axes pseudo-topology. This is done inside the action + // because we may not know about base pairing until the first frame. + if (stepAxesParm_ != 0) { + if (stepAxesParm_->Natom() > 0) { + //mprintf("\tBase pair step axes pseudo-topology is already set up.\n"); + // Check that number of base pair steps has not changed + if ((unsigned int)stepAxesParm_->Nres() != Steps_.size()) { + mprinterr("Error: Number of base pair steps has changed from %i to %zu.\n" + "Error: Base pair step axes pseudo-topology is already set up for %i bases.\n", + stepAxesParm_->Nres(), Steps_.size(), stepAxesParm_->Nres()); + return Action::ERR; + } + } else { + // Create residue information for each base pair step axes + std::vector stepAxesResidues; + stepAxesResidues.reserve( Steps_.size() ); + for (StepMap::const_iterator it = Steps_.begin(); it != Steps_.end(); ++it) + { + StepType const& BS = it->second; + NA_Base const& base1 = Bases_[BS.b1idx_]; + Residue const& res = setupTop_->Res( base1.ResNum() ); + stepAxesResidues.push_back( res ); + } + if (setup_axes_pseudoTraj( *stepAxesParm_, *stepAxesOut_, stepAxesFrame_, stepAxesResidues )) + return Action::ERR; + } + } + // Output base axes if needed if (axesOut_ != 0) { axesFrame_.ClearAtoms(); @@ -1954,6 +2009,13 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { } bpAxesOut_->WriteSingle(frm.TrajoutNum(), bpAxesFrame_); } + // Output base pair step axes if needed + if (stepAxesOut_ != 0) { + stepAxesFrame_.ClearAtoms(); + for (StepMap::const_iterator it = Steps_.begin(); it != Steps_.end(); ++it) + axesToFrame( stepAxesFrame_, it->second.stepaxis_ ); + stepAxesOut_->WriteSingle(frm.TrajoutNum(), stepAxesFrame_); + } nframes_++; return Action::OK; diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 0e5f7e0b7c..d1b32a5930 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -211,13 +211,16 @@ class Action_NAstruct: public Action { // DEBUG - used to trigger AxisPDBwriter for first call of calculateParameters bool calcparam_; # endif - Trajout_Single* axesOut_; ///< Output trajectory for base axes - Topology* axesParm_; ///< Pseudo-topology for base axes - Frame axesFrame_; ///< Frame for base axes pseudo traj - Trajout_Single* bpAxesOut_; ///< Output trajectory for base pair axes - Topology* bpAxesParm_; ///< Pseudo-topology for base pair axes - Frame bpAxesFrame_; ///< Frame for base pair axes pseudo traj - int setupNframes_; ///< Set in Setup(); number of expected frames to write (pseudo-traj) - Topology* setupTop_; ///< Set in Setup(); current topology + Trajout_Single* axesOut_; ///< Output trajectory for base axes + Topology* axesParm_; ///< Pseudo-topology for base axes + Frame axesFrame_; ///< Frame for base axes pseudo traj + Trajout_Single* bpAxesOut_; ///< Output trajectory for base pair axes + Topology* bpAxesParm_; ///< Pseudo-topology for base pair axes + Frame bpAxesFrame_; ///< Frame for base pair axes pseudo traj + Trajout_Single* stepAxesOut_; ///< Output trajectory for base pair step axes + Topology* stepAxesParm_; ///< Pseudo-topology for base pair step axes + Frame stepAxesFrame_; ///< Frame for base pair step axes pseudo traj + int setupNframes_; ///< Set in Setup(); number of expected frames to write (pseudo-traj) + Topology* setupTop_; ///< Set in Setup(); current topology }; #endif From 23abaab4b98bfb7e978af1ee20a38883918a2f88 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 08:57:34 -0400 Subject: [PATCH 18/39] Fix comments --- src/Action_NAstruct.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 43bb672cde..3f2520f2fa 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -1922,7 +1922,7 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { findBPmode_ = REFERENCE; } - // Set up base pair axes pseudo-topology. This is done inside the action + // Set up base pair axes pseudo-topology. This is done inside DoAction // because we may not know about base pairing until the first frame. if (bpAxesParm_ != 0) { if (bpAxesParm_->Natom() > 0) { @@ -1964,7 +1964,7 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { // Determine base pair step parameters DetermineStepParameters(frameNum); - // Set up base pair step axes pseudo-topology. This is done inside the action + // Set up base pair step axes pseudo-topology. This is done inside DoAction // because we may not know about base pairing until the first frame. if (stepAxesParm_ != 0) { if (stepAxesParm_->Natom() > 0) { From 380693151ff1e68f3c0c132ce809f71110b7f145 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 09:20:33 -0400 Subject: [PATCH 19/39] Move the base axes write to before base pair parameters are determined since that can flip base axes --- src/Action_NAstruct.cpp | 22 ++++++++++++---------- test/Test_NAstruct/RunTest.sh | 7 +++++-- 2 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 3f2520f2fa..a54cd74b6b 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -1121,8 +1121,8 @@ int Action_NAstruct::DetermineStrandParameters(int frameNum) { //mprintf("DEBUG: Strand %u, bases %u to %u\n", SP.strandidx_, b1idx, b2idx); // Get bases - NA_Base& base1 = Bases_[b1idx]; - NA_Base& base2 = Bases_[b2idx]; //TODO copy? + NA_Base const& base1 = Bases_[b1idx]; + NA_Base const& base2 = Bases_[b2idx]; // Calc parameters between bases in the strand calculateParameters(base2.Axis(), base1.Axis(), &commonAxis, Param); // Store data @@ -1958,6 +1958,16 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { if (sscalc_) DetermineStrandParameters(frameNum); + // Output base axes if needed. Do it here because DeterminePairParameters() + // can flip base axes. + if (axesOut_ != 0) { + axesFrame_.ClearAtoms(); + for (std::vector::const_iterator base = Bases_.begin(); + base != Bases_.end(); ++base) + axesToFrame( axesFrame_, base->Axis() ); + axesOut_->WriteSingle(frm.TrajoutNum(), axesFrame_); + } + // Determine base parameters DeterminePairParameters(frameNum); @@ -1992,14 +2002,6 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { } } - // Output base axes if needed - if (axesOut_ != 0) { - axesFrame_.ClearAtoms(); - for (std::vector::const_iterator base = Bases_.begin(); - base != Bases_.end(); ++base) - axesToFrame( axesFrame_, base->Axis() ); - axesOut_->WriteSingle(frm.TrajoutNum(), axesFrame_); - } // Output base pair axes if needed if (bpAxesOut_ != 0) { bpAxesFrame_.ClearAtoms(); diff --git a/test/Test_NAstruct/RunTest.sh b/test/Test_NAstruct/RunTest.sh index 2ab2c09890..fb71e3a562 100755 --- a/test/Test_NAstruct/RunTest.sh +++ b/test/Test_NAstruct/RunTest.sh @@ -4,7 +4,8 @@ # Clean CleanFiles nastruct.in BP.*.dat BPstep.*.dat bases.pdb baseaxes.pdb basepairaxes.pdb \ - Helix.*.dat Param.pdb SS.mol1.dat SS.mol1.selected.dat + Helix.*.dat Param.pdb SS.mol1.dat SS.mol1.selected.dat \ + axes.bases.pdb axes.bp.mol2 # Test 2 TESTNAME='NAstruct tests' @@ -14,7 +15,9 @@ INPUT="-i nastruct.in" cat > nastruct.in < Date: Tue, 11 Jul 2023 09:25:32 -0400 Subject: [PATCH 20/39] Add some axes write tests --- test/Test_NAstruct/RunTest.sh | 8 +- test/Test_NAstruct/axes.bases.pdb.save | 247 ++++++++++++++++++++++++ test/Test_NAstruct/axes.bp.mol2.save | 222 +++++++++++++++++++++ test/Test_NAstruct/axes.step.parm7.save | 134 +++++++++++++ 4 files changed, 609 insertions(+), 2 deletions(-) create mode 100644 test/Test_NAstruct/axes.bases.pdb.save create mode 100644 test/Test_NAstruct/axes.bp.mol2.save create mode 100644 test/Test_NAstruct/axes.step.parm7.save diff --git a/test/Test_NAstruct/RunTest.sh b/test/Test_NAstruct/RunTest.sh index fb71e3a562..5fcdee9e7b 100755 --- a/test/Test_NAstruct/RunTest.sh +++ b/test/Test_NAstruct/RunTest.sh @@ -5,7 +5,7 @@ # Clean CleanFiles nastruct.in BP.*.dat BPstep.*.dat bases.pdb baseaxes.pdb basepairaxes.pdb \ Helix.*.dat Param.pdb SS.mol1.dat SS.mol1.selected.dat \ - axes.bases.pdb axes.bp.mol2 + axes.bases.pdb axes.bp.mol2 axes.step.dcd axes.step.parm7 # Test 2 TESTNAME='NAstruct tests' @@ -17,7 +17,8 @@ parm ../adh026.3.pdb trajin ../adh026.3.pdb nastruct naout adh026.dat \ axesout axes.bases.pdb \ - bpaxesout axes.bp.mol2 + bpaxesout axes.bp.mol2 \ + stepaxesout axes.step.dcd stepaxesparmout axes.step.parm7 nastruct naout baseref.dat baseref Atomic_G.pdb.nastruct nastruct naout groove.dat groovecalc 3dna nastruct naout GuessBP.dat guessbp @@ -33,6 +34,9 @@ DoTest BPstep.groove.dat.save BPstep.groove.dat DoTest BP.adh026.dat.save BP.GuessBP.dat DoTest BPstep.adh026.dat.save BPstep.GuessBP.dat DoTest Helix.adh026.dat.save Helix.GuessBP.dat +DoTest axes.bases.pdb.save axes.bases.pdb +DoTest axes.bp.mol2.save axes.bp.mol2 +DoTest axes.step.parm7.save axes.step.parm7 -I %VERSION # Single strand cat > nastruct.in <MOLECULE +Cpptraj Generated mol2 file. + 32 24 8 0 0 +SMALL +NO_CHARGES + + +@ATOM + 1 Orig 16.1006 0.6576 -26.9263 Orig 1 G5 0.000000 + 2 X 16.3905 1.2916 -26.2094 X 1 G5 0.000000 + 3 Y 15.4695 1.3474 -27.2811 Y 1 G5 0.000000 + 4 Z 15.3811 0.3080 -26.3262 Z 1 G5 0.000000 + 5 Orig 14.9712 -1.5411 -24.3127 Orig 2 G 0.000000 + 6 X 14.9559 -0.5839 -24.0239 X 2 G 0.000000 + 7 Y 14.2470 -1.3525 -24.9761 Y 2 G 0.000000 + 8 Z 14.2816 -1.7604 -23.6224 Z 2 G 0.000000 + 9 Orig 13.5140 -2.0311 -21.2086 Orig 3 G 0.000000 + 10 X 13.1353 -1.1212 -21.3782 X 3 G 0.000000 + 11 Y 12.9133 -2.4121 -21.9115 Y 3 G 0.000000 + 12 Z 12.8100 -2.1954 -20.5178 Z 3 G 0.000000 + 13 Orig 12.2467 -1.4109 -17.6502 Orig 4 C 0.000000 + 14 X 11.6896 -0.8029 -18.2158 X 4 C 0.000000 + 15 Y 11.9148 -2.1983 -18.1697 Y 4 C 0.000000 + 16 Z 11.4854 -1.5126 -17.0098 Z 4 C 0.000000 + 17 Orig 10.2509 0.5164 -14.5640 Orig 5 G 0.000000 + 18 X 9.7225 0.8166 -15.3582 X 5 G 0.000000 + 19 Y 10.2890 -0.4097 -14.9394 Y 5 G 0.000000 + 20 Z 9.4027 0.2878 -14.0861 Z 5 G 0.000000 + 21 Orig 7.0428 1.9286 -12.1872 Orig 6 C 0.000000 + 22 X 6.6519 1.6464 -13.0633 X 6 C 0.000000 + 23 Y 7.3893 1.0017 -12.0432 Y 6 C 0.000000 + 24 Z 6.1901 1.6815 -11.7271 Z 6 C 0.000000 + 25 Orig 3.7050 3.0335 -10.5734 Orig 7 C 0.000000 + 26 X 3.5253 2.2660 -11.1886 X 7 C 0.000000 + 27 Y 4.3081 2.4535 -10.0258 Y 7 C 0.000000 + 28 Z 2.9278 2.7609 -10.0062 Z 7 C 0.000000 + 29 Orig -0.3242 2.4056 -10.1776 Orig 8 C3 0.000000 + 30 X -0.0571 1.4817 -10.4515 X 8 C3 0.000000 + 31 Y 0.2128 2.3123 -9.3392 Y 8 C3 0.000000 + 32 Z -1.1243 2.0346 -9.7064 Z 8 C3 0.000000 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 1 4 1 + 4 5 6 1 + 5 5 7 1 + 6 5 8 1 + 7 9 10 1 + 8 9 11 1 + 9 9 12 1 + 10 13 14 1 + 11 13 15 1 + 12 13 16 1 + 13 17 18 1 + 14 17 19 1 + 15 17 20 1 + 16 21 22 1 + 17 21 23 1 + 18 21 24 1 + 19 25 26 1 + 20 25 27 1 + 21 25 28 1 + 22 29 30 1 + 23 29 31 1 + 24 29 32 1 +@SUBSTRUCTURE + 1 G5 1 **** 0 **** **** + 2 G 5 **** 0 **** **** + 3 G 9 **** 0 **** **** + 4 C 13 **** 0 **** **** + 5 G 17 **** 0 **** **** + 6 C 21 **** 0 **** **** + 7 C 25 **** 0 **** **** + 8 C3 29 **** 0 **** **** +@MOLECULE +Cpptraj Generated mol2 file. + 32 24 8 0 0 +SMALL +NO_CHARGES + + +@ATOM + 1 Orig 16.0119 0.2594 -26.8514 Orig 1 G5 0.000000 + 2 X 16.2690 0.9081 -26.1351 X 1 G5 0.000000 + 3 Y 15.2651 0.8633 -27.1303 Y 1 G5 0.000000 + 4 Z 15.3984 -0.2038 -26.2118 Z 1 G5 0.000000 + 5 Orig 15.1578 -1.6816 -24.4587 Orig 2 G 0.000000 + 6 X 14.9827 -0.7932 -24.0345 X 2 G 0.000000 + 7 Y 14.4431 -1.5000 -25.1341 Y 2 G 0.000000 + 8 Z 14.4808 -2.1032 -23.8555 Z 2 G 0.000000 + 9 Orig 13.8510 -2.3965 -20.6422 Orig 3 G 0.000000 + 10 X 13.4383 -1.4932 -20.7592 X 3 G 0.000000 + 11 Y 13.2860 -2.7510 -21.3873 Y 3 G 0.000000 + 12 Z 13.1364 -2.6379 -19.9856 Z 3 G 0.000000 + 13 Orig 11.9209 -0.9002 -17.1342 Orig 4 C 0.000000 + 14 X 11.3209 -0.2633 -17.6183 X 4 C 0.000000 + 15 Y 11.6516 -1.6308 -17.7617 Y 4 C 0.000000 + 16 Z 11.1675 -1.1463 -16.5243 Z 4 C 0.000000 + 17 Orig 10.1037 0.8873 -13.8721 Orig 5 G 0.000000 + 18 X 9.5283 1.1627 -14.6423 X 5 G 0.000000 + 19 Y 10.1097 -0.0528 -14.2129 Y 5 G 0.000000 + 20 Z 9.2858 0.6866 -13.3329 Z 5 G 0.000000 + 21 Orig 7.0514 1.7403 -12.1675 Orig 6 C 0.000000 + 22 X 6.5167 1.3830 -12.9333 X 6 C 0.000000 + 23 Y 7.4021 0.8219 -11.9840 Y 6 C 0.000000 + 24 Z 6.2826 1.5698 -11.5512 Z 6 C 0.000000 + 25 Orig 3.8075 2.9895 -10.4554 Orig 7 C 0.000000 + 26 X 3.6634 2.2589 -11.1229 X 7 C 0.000000 + 27 Y 4.3526 2.3679 -9.8928 Y 7 C 0.000000 + 28 Z 2.9816 2.7066 -9.9676 Z 7 C 0.000000 + 29 Orig -0.1535 2.4052 -10.4691 Orig 8 C3 0.000000 + 30 X -0.0012 1.4882 -10.8378 X 8 C3 0.000000 + 31 Y 0.3973 2.1742 -9.6671 Y 8 C3 0.000000 + 32 Z -0.9742 2.0799 -9.9992 Z 8 C3 0.000000 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 1 4 1 + 4 5 6 1 + 5 5 7 1 + 6 5 8 1 + 7 9 10 1 + 8 9 11 1 + 9 9 12 1 + 10 13 14 1 + 11 13 15 1 + 12 13 16 1 + 13 17 18 1 + 14 17 19 1 + 15 17 20 1 + 16 21 22 1 + 17 21 23 1 + 18 21 24 1 + 19 25 26 1 + 20 25 27 1 + 21 25 28 1 + 22 29 30 1 + 23 29 31 1 + 24 29 32 1 +@SUBSTRUCTURE + 1 G5 1 **** 0 **** **** + 2 G 5 **** 0 **** **** + 3 G 9 **** 0 **** **** + 4 C 13 **** 0 **** **** + 5 G 17 **** 0 **** **** + 6 C 21 **** 0 **** **** + 7 C 25 **** 0 **** **** + 8 C3 29 **** 0 **** **** +@MOLECULE +Cpptraj Generated mol2 file. + 32 24 8 0 0 +SMALL +NO_CHARGES + + +@ATOM + 1 Orig 14.9252 0.4732 -27.2106 Orig 1 G5 0.000000 + 2 X 15.2002 1.1520 -26.5298 X 1 G5 0.000000 + 3 Y 14.1964 1.0822 -27.5235 Y 1 G5 0.000000 + 4 Z 14.2982 0.0630 -26.5484 Z 1 G5 0.000000 + 5 Orig 14.2491 -1.9434 -24.0972 Orig 2 G 0.000000 + 6 X 14.0678 -1.0191 -23.7613 X 2 G 0.000000 + 7 Y 13.5037 -1.8498 -24.7572 Y 2 G 0.000000 + 8 Z 13.6077 -2.3135 -23.4252 Z 2 G 0.000000 + 9 Orig 13.5563 -2.1715 -20.8287 Orig 3 G 0.000000 + 10 X 13.0907 -1.3020 -20.9937 X 3 G 0.000000 + 11 Y 12.9838 -2.6096 -21.5218 Y 3 G 0.000000 + 12 Z 12.8814 -2.3997 -20.1270 Z 3 G 0.000000 + 13 Orig 12.1978 -1.0008 -17.0450 Orig 4 C 0.000000 + 14 X 11.5728 -0.4185 -17.5649 X 4 C 0.000000 + 15 Y 11.8445 -1.8057 -17.5217 Y 4 C 0.000000 + 16 Z 11.5018 -1.1150 -16.3362 Z 4 C 0.000000 + 17 Orig 10.6600 0.6594 -13.5992 Orig 5 G 0.000000 + 18 X 10.0282 0.8970 -14.3370 X 5 G 0.000000 + 19 Y 10.6695 -0.2900 -13.9130 Y 5 G 0.000000 + 20 Z 9.8849 0.4542 -13.0016 Z 5 G 0.000000 + 21 Orig 7.3054 1.9585 -11.8952 Orig 6 C 0.000000 + 22 X 6.8706 1.6565 -12.7436 X 6 C 0.000000 + 23 Y 7.6425 1.0303 -11.7377 Y 6 C 0.000000 + 24 Z 6.4704 1.7410 -11.3898 Z 6 C 0.000000 + 25 Orig 3.9524 3.3271 -10.2654 Orig 7 C 0.000000 + 26 X 3.7554 2.5924 -10.9144 X 7 C 0.000000 + 27 Y 4.5627 2.7171 -9.7600 Y 7 C 0.000000 + 28 Z 3.1851 3.0306 -9.6968 Z 7 C 0.000000 + 29 Orig -0.2495 3.0294 -10.2703 Orig 8 C3 0.000000 + 30 X 0.0856 2.1591 -10.6312 X 8 C3 0.000000 + 31 Y 0.2817 2.8876 -9.4350 Y 8 C3 0.000000 + 32 Z -1.0276 2.5578 -9.8555 Z 8 C3 0.000000 +@BOND + 1 1 2 1 + 2 1 3 1 + 3 1 4 1 + 4 5 6 1 + 5 5 7 1 + 6 5 8 1 + 7 9 10 1 + 8 9 11 1 + 9 9 12 1 + 10 13 14 1 + 11 13 15 1 + 12 13 16 1 + 13 17 18 1 + 14 17 19 1 + 15 17 20 1 + 16 21 22 1 + 17 21 23 1 + 18 21 24 1 + 19 25 26 1 + 20 25 27 1 + 21 25 28 1 + 22 29 30 1 + 23 29 31 1 + 24 29 32 1 +@SUBSTRUCTURE + 1 G5 1 **** 0 **** **** + 2 G 5 **** 0 **** **** + 3 G 9 **** 0 **** **** + 4 C 13 **** 0 **** **** + 5 G 17 **** 0 **** **** + 6 C 21 **** 0 **** **** + 7 C 25 **** 0 **** **** + 8 C3 29 **** 0 **** **** diff --git a/test/Test_NAstruct/axes.step.parm7.save b/test/Test_NAstruct/axes.step.parm7.save new file mode 100644 index 0000000000..d37d7d0a24 --- /dev/null +++ b/test/Test_NAstruct/axes.step.parm7.save @@ -0,0 +1,134 @@ +%VERSION VERSION_STAMP = V0001.000 DATE = 07/11/23 09:21:41 +%FLAG TITLE +%FORMAT(20a4) +StepAxes +%FLAG POINTERS +%FORMAT(10I8) + 28 0 21 0 0 0 0 0 0 0 + 49 7 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 0 4 0 + 0 +%FLAG ATOM_NAME +%FORMAT(20a4) +OrigX Y Z OrigX Y Z OrigX Y Z OrigX Y Z OrigX Y Z +OrigX Y Z OrigX Y Z +%FLAG CHARGE +%FORMAT(5E16.8) + 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 + 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 + 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 + 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 + 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 + 0.00000000E+00 0.00000000E+00 0.00000000E+00 +%FLAG ATOMIC_NUMBER +%FORMAT(10I8) + 6 1 1 1 6 1 1 1 6 1 + 1 1 6 1 1 1 6 1 1 1 + 6 1 1 1 6 1 1 1 +%FLAG MASS +%FORMAT(5E16.8) + 1.20107000E+01 1.00794000E+00 1.00794000E+00 1.00794000E+00 1.20107000E+01 + 1.00794000E+00 1.00794000E+00 1.00794000E+00 1.20107000E+01 1.00794000E+00 + 1.00794000E+00 1.00794000E+00 1.20107000E+01 1.00794000E+00 1.00794000E+00 + 1.00794000E+00 1.20107000E+01 1.00794000E+00 1.00794000E+00 1.00794000E+00 + 1.20107000E+01 1.00794000E+00 1.00794000E+00 1.00794000E+00 1.20107000E+01 + 1.00794000E+00 1.00794000E+00 1.00794000E+00 +%FLAG ATOM_TYPE_INDEX +%FORMAT(10I8) + 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 1 1 + 1 1 1 1 1 1 1 1 +%FLAG NUMBER_EXCLUDED_ATOMS +%FORMAT(10I8) + 3 2 1 1 3 2 1 1 3 2 + 1 1 3 2 1 1 3 2 1 1 + 3 2 1 1 3 2 1 1 +%FLAG NONBONDED_PARM_INDEX +%FORMAT(10I8) + +%FLAG RESIDUE_LABEL +%FORMAT(20a4) +G5 G G C G C C +%FLAG RESIDUE_POINTER +%FORMAT(10I8) + 1 5 9 13 17 21 25 +%FLAG BOND_FORCE_CONSTANT +%FORMAT(5E16.8) + 0.00000000E+00 +%FLAG BOND_EQUIL_VALUE +%FORMAT(5E16.8) + 1.00000000E+00 +%FLAG ANGLE_FORCE_CONSTANT +%FORMAT(5E16.8) + +%FLAG ANGLE_EQUIL_VALUE +%FORMAT(5E16.8) + +%FLAG DIHEDRAL_FORCE_CONSTANT +%FORMAT(5E16.8) + +%FLAG DIHEDRAL_PERIODICITY +%FORMAT(5E16.8) + +%FLAG DIHEDRAL_PHASE +%FORMAT(5E16.8) + +%FLAG SCEE_SCALE_FACTOR +%FORMAT(5E16.8) + +%FLAG SCNB_SCALE_FACTOR +%FORMAT(5E16.8) + +%FLAG SOLTY +%FORMAT(5E16.8) + +%FLAG LENNARD_JONES_ACOEF +%FORMAT(5E16.8) + +%FLAG LENNARD_JONES_BCOEF +%FORMAT(5E16.8) + +%FLAG BONDS_INC_HYDROGEN +%FORMAT(10I8) + 0 3 1 0 6 1 0 9 1 12 + 15 1 12 18 1 12 21 1 24 27 + 1 24 30 1 24 33 1 36 39 1 + 36 42 1 36 45 1 48 51 1 48 + 54 1 48 57 1 60 63 1 60 66 + 1 60 69 1 72 75 1 72 78 1 + 72 81 1 +%FLAG BONDS_WITHOUT_HYDROGEN +%FORMAT(10I8) + +%FLAG ANGLES_INC_HYDROGEN +%FORMAT(10I8) + +%FLAG ANGLES_WITHOUT_HYDROGEN +%FORMAT(10I8) + +%FLAG DIHEDRALS_INC_HYDROGEN +%FORMAT(10I8) + +%FLAG DIHEDRALS_WITHOUT_HYDROGEN +%FORMAT(10I8) + +%FLAG EXCLUDED_ATOMS_LIST +%FORMAT(10I8) + 2 3 4 3 4 4 0 6 7 8 + 7 8 8 0 10 11 12 11 12 12 + 0 14 15 16 15 16 16 0 18 19 + 20 19 20 20 0 22 23 24 23 24 + 24 0 26 27 28 27 28 28 0 +%FLAG HBOND_ACOEF +%FORMAT(5E16.8) + +%FLAG HBOND_BCOEF +%FORMAT(5E16.8) + +%FLAG HBCUT +%FORMAT(5E16.8) + +%FLAG AMBER_ATOM_TYPE +%FORMAT(20a4) + + From 59fc90d192397b6fffb5ca6250664bc2ddb8d968 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 09:27:12 -0400 Subject: [PATCH 21/39] Add test for step axes coords --- test/Test_NAstruct/RunTest.sh | 5 +++-- test/Test_NAstruct/axes.step.crd.save | 28 +++++++++++++++++++++++++++ 2 files changed, 31 insertions(+), 2 deletions(-) create mode 100644 test/Test_NAstruct/axes.step.crd.save diff --git a/test/Test_NAstruct/RunTest.sh b/test/Test_NAstruct/RunTest.sh index 5fcdee9e7b..d881fd2ff5 100755 --- a/test/Test_NAstruct/RunTest.sh +++ b/test/Test_NAstruct/RunTest.sh @@ -5,7 +5,7 @@ # Clean CleanFiles nastruct.in BP.*.dat BPstep.*.dat bases.pdb baseaxes.pdb basepairaxes.pdb \ Helix.*.dat Param.pdb SS.mol1.dat SS.mol1.selected.dat \ - axes.bases.pdb axes.bp.mol2 axes.step.dcd axes.step.parm7 + axes.bases.pdb axes.bp.mol2 axes.step.crd axes.step.parm7 # Test 2 TESTNAME='NAstruct tests' @@ -18,7 +18,7 @@ trajin ../adh026.3.pdb nastruct naout adh026.dat \ axesout axes.bases.pdb \ bpaxesout axes.bp.mol2 \ - stepaxesout axes.step.dcd stepaxesparmout axes.step.parm7 + stepaxesout axes.step.crd stepaxesparmout axes.step.parm7 nastruct naout baseref.dat baseref Atomic_G.pdb.nastruct nastruct naout groove.dat groovecalc 3dna nastruct naout GuessBP.dat guessbp @@ -36,6 +36,7 @@ DoTest BPstep.adh026.dat.save BPstep.GuessBP.dat DoTest Helix.adh026.dat.save Helix.GuessBP.dat DoTest axes.bases.pdb.save axes.bases.pdb DoTest axes.bp.mol2.save axes.bp.mol2 +DoTest axes.step.crd.save axes.step.crd DoTest axes.step.parm7.save axes.step.parm7 -I %VERSION # Single strand diff --git a/test/Test_NAstruct/axes.step.crd.save b/test/Test_NAstruct/axes.step.crd.save new file mode 100644 index 0000000000..314759889a --- /dev/null +++ b/test/Test_NAstruct/axes.step.crd.save @@ -0,0 +1,28 @@ +Cpptraj Generated trajectory + 15.536 -0.442 -25.620 15.681 0.395 -25.092 14.844 0.026 -26.169 14.829 + -0.727 -24.972 14.243 -1.786 -22.761 14.036 -0.810 -22.698 13.556 -1.885 + -23.481 13.546 -1.978 -22.070 12.880 -1.721 -19.429 12.393 -0.935 -19.809 + 12.407 -2.324 -20.071 12.147 -1.854 -18.763 11.249 -0.447 -16.107 10.680 + 0.010 -16.791 11.105 -1.321 -16.572 10.439 -0.613 -15.544 8.647 1.223 + -13.376 8.162 1.231 -14.250 8.851 0.251 -13.498 7.796 0.985 -12.907 + 5.374 2.481 -11.380 5.063 1.929 -12.154 5.860 1.689 -11.010 4.557 + 2.221 -10.865 1.690 2.720 -10.375 1.741 1.838 -10.845 2.301 2.375 + -9.662 0.900 2.397 -9.855 + 15.585 -0.711 -25.655 15.619 0.086 -25.052 14.822 -0.300 -26.154 14.939 + -1.154 -25.033 14.504 -2.039 -22.550 14.200 -1.100 -22.392 13.857 -2.121 + -23.308 13.805 -2.372 -21.918 12.886 -1.648 -18.888 12.359 -0.856 -19.194 + 12.458 -2.207 -19.599 12.152 -1.892 -18.255 11.012 -0.006 -15.503 10.408 + 0.462 -16.147 10.887 -0.861 -16.007 10.226 -0.230 -14.928 8.578 1.314 + -13.020 7.996 1.271 -13.832 8.753 0.332 -13.094 7.783 1.128 -12.441 + 5.429 2.365 -11.311 5.075 1.799 -12.056 5.913 1.573 -10.939 4.629 + 2.137 -10.757 1.827 2.697 -10.462 1.830 1.850 -10.994 2.394 2.261 + -9.764 1.003 2.393 -9.983 + 14.587 -0.735 -25.654 14.634 0.107 -25.117 13.816 -0.363 -26.170 13.953 + -1.125 -24.987 13.903 -2.057 -22.463 13.567 -1.119 -22.376 13.230 -2.231 + -23.183 13.243 -2.357 -21.774 12.877 -1.586 -18.937 12.321 -0.836 -19.295 + 12.409 -2.225 -19.547 12.190 -1.758 -18.230 11.429 -0.171 -15.322 10.776 + 0.244 -15.956 11.258 -1.066 -15.733 10.691 -0.331 -14.667 8.983 1.309 + -12.747 8.424 1.275 -13.576 9.177 0.332 -12.838 8.176 1.097 -12.195 + 5.629 2.643 -11.080 5.288 2.101 -11.849 6.118 1.843 -10.733 4.826 + 2.385 -10.542 1.851 3.178 -10.268 1.929 2.337 -10.803 2.475 2.800 + -9.584 1.073 2.792 -9.773 From edb74620bf7cbe5fa9f2e6431430f76d6d96a954 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 09:37:30 -0400 Subject: [PATCH 22/39] Add ability to customize atom axis names --- src/Action_NAstruct.cpp | 16 ++++++++++++---- src/Action_NAstruct.h | 4 ++++ 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index a54cd74b6b..66202a5f1a 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -217,6 +217,10 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int "StepAxes", init.DSL(), actionArgs, &stepAxesOut_, &stepAxesParm_)) return Action::ERR; + axesNameO_ = actionArgs.GetStringKey("axesnameo", "Orig"); + axesNameX_ = actionArgs.GetStringKey("axesnamex", "X"); + axesNameY_ = actionArgs.GetStringKey("axesnamey", "Y"); + axesNameZ_ = actionArgs.GetStringKey("axesnamez", "Z"); // Get residue range resRange_.SetRange(actionArgs.GetStringKey("resrange")); if (!resRange_.Empty()) @@ -385,6 +389,10 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int if (!stepAxesParm_->OriginalFilename().empty()) mprintf("\tWriting base pair step axes pseudo topology to '%s'\n", stepAxesParm_->OriginalFilename().full()); } + if (axesOut_ != 0 || bpAxesOut_ != 0 || stepAxesOut_ != 0) { + mprintf("\tAxes pseudo atom names: origin='%s' X='%s' Y='%s' Z='%s'\n", + *axesNameO_, *axesNameX_, *axesNameY_, axesNameZ_); + } mprintf("# Citations: Babcock MS; Pednault EPD; Olson WK; \"Nucleic Acid Structure\n" "# Analysis: Mathematics for Local Cartesian and Helical Structure\n" "# Parameters That Are Truly Comparable Between Structures\",\n" @@ -1835,10 +1843,10 @@ const ++res) { // Order is origin, x, y, z - pseudo.AddTopAtom(Atom("Orig", "C"), *res); - pseudo.AddTopAtom(Atom("X", "H"), *res); - pseudo.AddTopAtom(Atom("Y", "H"), *res); - pseudo.AddTopAtom(Atom("Z", "H"), *res); + pseudo.AddTopAtom(Atom(axesNameO_, "C"), *res); + pseudo.AddTopAtom(Atom(axesNameX_, "H"), *res); + pseudo.AddTopAtom(Atom(axesNameY_, "H"), *res); + pseudo.AddTopAtom(Atom(axesNameZ_, "H"), *res); // Bond x y and z to origin pseudo.AddBond(natom, natom+1, 0); pseudo.AddBond(natom, natom+2, 0); diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index d1b32a5930..d413d06764 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -222,5 +222,9 @@ class Action_NAstruct: public Action { Frame stepAxesFrame_; ///< Frame for base pair step axes pseudo traj int setupNframes_; ///< Set in Setup(); number of expected frames to write (pseudo-traj) Topology* setupTop_; ///< Set in Setup(); current topology + NameType axesNameO_; + NameType axesNameX_; + NameType axesNameY_; + NameType axesNameZ_; }; #endif From 45e38f1eee5ced1101f213ce98b88d023d09ef1a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 09:40:41 -0400 Subject: [PATCH 23/39] Update help and fix var names --- src/Action_NAstruct.cpp | 19 ++++++++++--------- src/Action_NAstruct.h | 8 ++++---- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 66202a5f1a..e62c3d076a 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -129,6 +129,7 @@ void Action_NAstruct::Help() const { "\t[axesout [axesoutarg ...] [axesparmout ]]\n" "\t[bpaxesout [bpaxesoutarg ...] [bpaxesparmout ]]\n" "\t[stepaxesout [stepaxesoutarg ...] [stepaxesparmout ]]\n" + "\t[axisnameo ] [axisnamex ] [axisnamey ] [axisnamez ]\n" "\t[{ %s | allframes}]\n", DataSetList::RefArgs); mprintf(" Perform nucleic acid structure analysis. Base pairing can be determined\n" " in multiple ways:\n" @@ -217,10 +218,10 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int "StepAxes", init.DSL(), actionArgs, &stepAxesOut_, &stepAxesParm_)) return Action::ERR; - axesNameO_ = actionArgs.GetStringKey("axesnameo", "Orig"); - axesNameX_ = actionArgs.GetStringKey("axesnamex", "X"); - axesNameY_ = actionArgs.GetStringKey("axesnamey", "Y"); - axesNameZ_ = actionArgs.GetStringKey("axesnamez", "Z"); + axisNameO_ = actionArgs.GetStringKey("axisnameo", "Orig"); + axisNameX_ = actionArgs.GetStringKey("axisnamex", "X"); + axisNameY_ = actionArgs.GetStringKey("axisnamey", "Y"); + axisNameZ_ = actionArgs.GetStringKey("axisnamez", "Z"); // Get residue range resRange_.SetRange(actionArgs.GetStringKey("resrange")); if (!resRange_.Empty()) @@ -391,7 +392,7 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int } if (axesOut_ != 0 || bpAxesOut_ != 0 || stepAxesOut_ != 0) { mprintf("\tAxes pseudo atom names: origin='%s' X='%s' Y='%s' Z='%s'\n", - *axesNameO_, *axesNameX_, *axesNameY_, axesNameZ_); + *axisNameO_, *axisNameX_, *axisNameY_, axisNameZ_); } mprintf("# Citations: Babcock MS; Pednault EPD; Olson WK; \"Nucleic Acid Structure\n" "# Analysis: Mathematics for Local Cartesian and Helical Structure\n" @@ -1843,10 +1844,10 @@ const ++res) { // Order is origin, x, y, z - pseudo.AddTopAtom(Atom(axesNameO_, "C"), *res); - pseudo.AddTopAtom(Atom(axesNameX_, "H"), *res); - pseudo.AddTopAtom(Atom(axesNameY_, "H"), *res); - pseudo.AddTopAtom(Atom(axesNameZ_, "H"), *res); + pseudo.AddTopAtom(Atom(axisNameO_, "C"), *res); + pseudo.AddTopAtom(Atom(axisNameX_, "H"), *res); + pseudo.AddTopAtom(Atom(axisNameY_, "H"), *res); + pseudo.AddTopAtom(Atom(axisNameZ_, "H"), *res); // Bond x y and z to origin pseudo.AddBond(natom, natom+1, 0); pseudo.AddBond(natom, natom+2, 0); diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index d413d06764..877a5f8194 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -222,9 +222,9 @@ class Action_NAstruct: public Action { Frame stepAxesFrame_; ///< Frame for base pair step axes pseudo traj int setupNframes_; ///< Set in Setup(); number of expected frames to write (pseudo-traj) Topology* setupTop_; ///< Set in Setup(); current topology - NameType axesNameO_; - NameType axesNameX_; - NameType axesNameY_; - NameType axesNameZ_; + NameType axisNameO_; + NameType axisNameX_; + NameType axisNameY_; + NameType axisNameZ_; }; #endif From 8f5899914473471fd0f0f357139b4f98d0dfbc52 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 09:42:50 -0400 Subject: [PATCH 24/39] Remove old code --- src/Action_NAstruct.cpp | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index e62c3d076a..d4eff900c6 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -755,28 +755,6 @@ int Action_NAstruct::DetermineBasePairing() { # endif // Stagger (vertical separation) must be less than a cutoff. if ( fabs(Param[2]) < staggerCut_ ) { -/* // Figure out if z vectors point in same (<90 deg) or opposite (>90 deg) direction - bool AntiParallel; - double theta = base1->Axis().Rz().Angle( base2->Axis().Rz() ); - double t_delta; // Deviation from linear - if (theta > Constants::PIOVER2) { // If theta(Z) > 90 deg. -# ifdef NASTRUCTDEBUG - mprintf("\t%s is anti-parallel to %s (%g deg)\n", base1->ResName(), base2->ResName(), - theta * Constants::RADDEG); -# endif - AntiParallel = true; - t_delta = Constants::PI - theta; - } else { -# ifdef NASTRUCTDEBUG - mprintf("\t%s is parallel to %s (%g deg)\n", base1->ResName(), base2->ResName(), - theta * Constants::RADDEG); -# endif - AntiParallel = false; - t_delta = theta; - } -# ifdef NASTRUCTDEBUG - mprintf("\tDeviation from linear: %g deg.\n", t_delta * Constants::RADDEG); -# endif*/ // Deviation from linear must be less than cutoff if (z_deviation_from_linear < z_angle_cut_) { int NHB = CalcNumHB(*base1, *base2, n_wc_hb); From 8cc4547d632f44c2eccb7f122224d301e4836f24 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 10:00:32 -0400 Subject: [PATCH 25/39] Start adding user-specified base pairing --- src/Action_NAstruct.cpp | 95 +++++++++++++++++++++++++++++++++++++++++ src/Action_NAstruct.h | 4 ++ 2 files changed, 99 insertions(+) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index d4eff900c6..539ab72f76 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -665,6 +665,101 @@ Action_NAstruct::BPmap::iterator return entry; } +/** User-specified base pairing. */ +int Action_NAstruct::SpecifiedBasePairing() { + int n_wc_hb; +# ifdef NASTRUCTDEBUG + mprintf("\n=================== Specified Base Pairing ===================\n"); +# endif + for (PairArray::const_iterator it = specifiedPairs_.begin(); + it != specifiedPairs_.end(); + ++it) + { + // User-specified base pair #s start from 1 + int b1idx = it->first - 1; + int b2idx = it->second - 1; + if (b1idx < 0 || (unsigned int)b1idx >= Bases_.size()) { + mprinterr("Error: Specified base # %u is out of range.\n", it->first); + return 1; + } + NA_Base& base1 = Bases_[b1idx]; + if (b2idx < 0 || (unsigned int)b2idx >= Bases_.size()) { + mprinterr("Error: Specified base # %u is out of range.\n", it->second); + return 1; + } + NA_Base& base2 = Bases_[b2idx]; +/*# ifdef NASTRUCTDEBUG + // Glycosidic N-N distance + if (base1.HasNXatom() && base2.HasNXatom()) { + double n_n_dist2 = DIST2_NoImage(base1.NXxyz(), base2.NXxyz()); + mprintf("DEBUG: NX-NX distance= %f\n", sqrt(n_n_dist2)); + } +# endif*/ + NA_Axis b1Axis = base1.Axis(); + NA_Axis b2Axis = base2.Axis(); + // Determine if base Z axis vectors are aligned with strand direction + bool is_z = false; + int b1_5to3 = axis_points_5p_to_3p( base1 ); + int b2_5to3 = axis_points_5p_to_3p( base2 ); + // TODO trap errors here + // If antiparallel and both bases are aligned 3' to 5', may be ZDNA + if (b1_5to3 == 0 && b2_5to3 == 0) { +# ifdef NASTRUCTDEBUG + mprintf("Both bases aligned 3' to 5', ZDNA\n"); +# endif + b1Axis.FlipXZ(); + b2Axis.FlipXZ(); + is_z = true; + } + // Determine if base Z axis vectors point in same (theta <= 90) or + // opposite (theta > 90) directions. + bool is_antiparallel; + double z_theta = b1Axis.Rz().Angle( b2Axis.Rz() ); + double z_deviation_from_linear; + if (z_theta > Constants::PIOVER2) { // If theta(Z) > 90 deg. +# ifdef NASTRUCTDEBUG + mprintf("\t%s is anti-parallel to %s (%g deg)\n", base1.ResName(), base2.ResName(), + z_theta * Constants::RADDEG); +# endif + is_antiparallel = true; + z_deviation_from_linear = Constants::PI - z_theta; + // Antiparallel - flip Y and Z axes of complimentary base + b2Axis.FlipYZ(); + } else { +# ifdef NASTRUCTDEBUG + mprintf("\t%s is parallel to %s (%g deg)\n", base1.ResName(), base2.ResName(), + z_theta * Constants::RADDEG); +# endif + is_antiparallel = false; + z_deviation_from_linear = z_theta; + // Parallel - no flip needed if 3dna. + // If using Babcock convention, flip X and Y axes. + if (bpConvention_ == BP_BABCOCK) + b2Axis.FlipXY(); + } +# ifdef NASTRUCTDEBUG + mprintf("\tDeviation from linear: %g deg.\n", z_deviation_from_linear * Constants::RADDEG); + // Calculate parameters between axes. + double Param[6]; + calculateParameters(b2Axis, b1Axis, 0, Param); + mprintf(" Shear= %6.2f Stretch= %6.2f Stagger= %6.2f Buck= %6.2f Prop= %6.2f Open= %6.2f\n", + Param[0], Param[1], Param[2], Param[5]*Constants::RADDEG, Param[4]*Constants::RADDEG, Param[3]*Constants::RADDEG); +# endif + int NHB = CalcNumHB(base1, base2, n_wc_hb); + BPmap::iterator entry = AddBasePair(b1idx, base1, b2idx, base2); +# ifdef NASTRUCTDEBUG + mprintf(", %i hbonds\n", NHB); +# endif + entry->second.nhb_ = NHB; + entry->second.n_wc_hb_ = n_wc_hb; + entry->second.isAnti_ = is_antiparallel; + entry->second.isZ_ = is_z; + } // END loop over pairs + + return 0; + +} + // Action_NAstruct::DetermineBasePairing() /** Determine which bases are paired from the individual base axes and set up * entry in BasePairs_ if one not already present. diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 877a5f8194..ea6f2023da 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -128,6 +128,7 @@ class Action_NAstruct: public Action { typedef std::map StepMap; ///< Map of BP indices to Steps typedef std::vector StrandArray; ///< Hold indices into Bases_ for strand beg/end typedef std::map RefMapType; ///< Map custom res names to target types + typedef std::vector< std::pair > PairArray; ///< Specified base pair #s // ----- Functions --------------------------- /// Initialize an axes pseudo-trajectory int init_axes_pseudoTraj(const char*, const char*, const char*, @@ -155,6 +156,8 @@ class Action_NAstruct: public Action { BPmap::iterator AddBasePair(int, NA_Base const&, int, NA_Base const&); /// Determine which bases are paired geometrically, set base pair data. int DetermineBasePairing(); + /// Set up base pairs based on user specification + int SpecifiedBasePairing(); /// Calculate translational/rotational parameters between two axes. int calculateParameters(NA_Axis const&, NA_Axis const&, NA_Axis*, double*); /// Calculate helical parameters between two axes. @@ -183,6 +186,7 @@ class Action_NAstruct: public Action { BPmap BasePairs_; ///< Hold base pairs StepMap Steps_; ///< Hold base pair steps. StrandArray Strands_; ///< Hold strand info + PairArray specifiedPairs_; ///< User-specified base pairing NA_Base::PmethodType puckerMethod_; ///< Pucker calculation method. double HBdistCut2_; ///< distance Cutoff^2 for determining hydrogen bonds double originCut2_; ///< Cutoff^2 for determining base-pairing vi origins From 6b8257418ba84acd7a5197a4a8b50d729bc86f01 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 10:30:01 -0400 Subject: [PATCH 26/39] Enable specifiedbp mode --- src/Action_NAstruct.cpp | 51 ++++++++++++++++++++++++++++++++++++++++- src/Action_NAstruct.h | 2 +- 2 files changed, 51 insertions(+), 2 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 539ab72f76..de5645db69 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -130,7 +130,9 @@ void Action_NAstruct::Help() const { "\t[bpaxesout [bpaxesoutarg ...] [bpaxesparmout ]]\n" "\t[stepaxesout [stepaxesoutarg ...] [stepaxesparmout ]]\n" "\t[axisnameo ] [axisnamex ] [axisnamey ] [axisnamez ]\n" - "\t[{ %s | allframes}]\n", DataSetList::RefArgs); + "\t[{ %s |\n" + "\t allframes |\n" + "\t specifiedbp pairs -,... }]\n", DataSetList::RefArgs); mprintf(" Perform nucleic acid structure analysis. Base pairing can be determined\n" " in multiple ways:\n" " - If 'first' (default) or a reference is specified, determine base\n" @@ -236,6 +238,8 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int findBPmode_ = REFERENCE; else if (actionArgs.hasKey("allframes")) findBPmode_ = ALL; + else if (actionArgs.hasKey("specifiedbp")) + findBPmode_ = SPECIFIED; else if (actionArgs.hasKey("guessbp")) { mprintf("Warning: 'guessbp' is deprecated. Defaulting to 'first'.\n"); findBPmode_ = FIRST; @@ -250,6 +254,41 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int return Action::ERR; } # endif + // Check for user-specified base pairs + if (findBPmode_ == SPECIFIED) { + std::string pairsarg = actionArgs.GetStringKey("pairs"); + // Format is b1-b2,... + while (!pairsarg.empty()) { + ArgList pairslist( pairsarg, "," ); + for (int iarg = 0; iarg < pairslist.Nargs(); iarg++) { + ArgList bpair( pairslist[iarg], "-" ); + if (bpair.Nargs() != 2) { + mprinterr("Error: Malformed base pair argument: %s\n", pairslist[iarg].c_str()); + return Action::ERR; + } + if (!validInteger(bpair[0])) { + mprinterr("Error: Expected an integer, got '%s'\n", bpair[0].c_str()); + return Action::ERR; + } + if (!validInteger(bpair[1])) { + mprinterr("Error: Expected an integer, got '%s'\n", bpair[1].c_str()); + return Action::ERR; + } + int b1idx = convertToInteger(bpair[0]); + int b2idx = convertToInteger(bpair[1]); + if (b1idx < 1 || b2idx < 1) { + mprinterr("Error: Base pair arg '%s', base #s must be > 0.\n", pairslist[iarg].c_str()); + return Action::ERR; + } + specifiedPairs_.push_back( std::pair( b1idx, b2idx ) ); + } // END loop over specified pairs + pairsarg = actionArgs.GetStringKey("pairs"); + } // END checking for 'pairs' keywords + if (specifiedPairs_.empty()) { + mprinterr("Error: No 'pairs' arguments for 'specifiedbp'\n"); + return Action::ERR; + } + } // Check for base pair mode incompatibilities if (findBPmode_ == ALL) { if (bpAxesOut_ != 0) { @@ -359,6 +398,11 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int // Determine Base Pairing if ( DetermineBasePairing() ) return Action::ERR; mprintf("\tSet up %zu base pairs.\n", BasePairs_.size() ); + } else if (findBPmode_ == SPECIFIED) { + mprintf("\tUser specified base pairs:"); + for (PairArray::const_iterator it = specifiedPairs_.begin(); it != specifiedPairs_.end(); ++it) + mprintf(" %u-%u", it->first, it->second); + mprintf("\n"); } else if (findBPmode_ == ALL) mprintf("\tBase pairs will be determined for each frame.\n"); else if (findBPmode_ == FIRST) @@ -1971,6 +2015,11 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { // calculated as part of determining base pairing. if ( SetupBaseAxes(frm.Frm()) ) return Action::ERR; if ( DetermineBasePairing() ) return Action::ERR; + } else if ( findBPmode_ == SPECIFIED ) { + // User-specified base pairing. + if ( SetupBaseAxes(frm.Frm()) ) return Action::ERR; + if ( SpecifiedBasePairing() ) return Action::ERR; + findBPmode_ = REFERENCE; } else if ( findBPmode_ == FIRST) { // Base pairs need to be determined from first frame. # ifdef MPI diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index ea6f2023da..575221b84b 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -51,7 +51,7 @@ class Action_NAstruct: public Action { enum GrooveType { PP_OO = 0, HASSAN_CALLADINE }; enum BP_ConventionType { BP_3DNA = 0, BP_BABCOCK }; /// How to find base pairs: first frame, reference structure, all frames. - enum FindType { FIRST = 0, REFERENCE, ALL }; + enum FindType { FIRST = 0, REFERENCE, ALL, SPECIFIED }; // ----- Data Structures --------------------- /// Hold consecutive bases struct Stype { From 22e355e17375d21339a994444de0dbf8f621b178 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 10:30:14 -0400 Subject: [PATCH 27/39] Add test for user-specified base pairing --- test/Test_NAstruct/RunTest.sh | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/test/Test_NAstruct/RunTest.sh b/test/Test_NAstruct/RunTest.sh index d881fd2ff5..8e4c9b4a19 100755 --- a/test/Test_NAstruct/RunTest.sh +++ b/test/Test_NAstruct/RunTest.sh @@ -5,7 +5,8 @@ # Clean CleanFiles nastruct.in BP.*.dat BPstep.*.dat bases.pdb baseaxes.pdb basepairaxes.pdb \ Helix.*.dat Param.pdb SS.mol1.dat SS.mol1.selected.dat \ - axes.bases.pdb axes.bp.mol2 axes.step.crd axes.step.parm7 + axes.bases.pdb axes.bp.mol2 axes.step.crd axes.step.parm7 \ + BP.specified.dat BPstep.specified.dat Helix.specified.dat # Test 2 TESTNAME='NAstruct tests' @@ -39,6 +40,18 @@ DoTest axes.bp.mol2.save axes.bp.mol2 DoTest axes.step.crd.save axes.step.crd DoTest axes.step.parm7.save axes.step.parm7 -I %VERSION +# User-specified base pairing +cat > nastruct.in < nastruct.in < Date: Tue, 11 Jul 2023 10:32:12 -0400 Subject: [PATCH 28/39] Formally deprecate guessbp in favor of specifiedbp --- src/Action_NAstruct.cpp | 4 ++-- test/Test_NAstruct/RunTest.sh | 4 ---- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index de5645db69..5f094e9b3f 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -241,8 +241,8 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int else if (actionArgs.hasKey("specifiedbp")) findBPmode_ = SPECIFIED; else if (actionArgs.hasKey("guessbp")) { - mprintf("Warning: 'guessbp' is deprecated. Defaulting to 'first'.\n"); - findBPmode_ = FIRST; + mprinterr("Error: 'guessbp' is deprecated. Consider using 'specifiedbp' instead.\n"); + return Action::ERR; } else if (actionArgs.hasKey("first")) findBPmode_ = FIRST; else diff --git a/test/Test_NAstruct/RunTest.sh b/test/Test_NAstruct/RunTest.sh index 8e4c9b4a19..348401efc7 100755 --- a/test/Test_NAstruct/RunTest.sh +++ b/test/Test_NAstruct/RunTest.sh @@ -22,7 +22,6 @@ nastruct naout adh026.dat \ stepaxesout axes.step.crd stepaxesparmout axes.step.parm7 nastruct naout baseref.dat baseref Atomic_G.pdb.nastruct nastruct naout groove.dat groovecalc 3dna -nastruct naout GuessBP.dat guessbp EOF RunCpptraj "NAstruct command test." DoTest BP.adh026.dat.save BP.adh026.dat @@ -32,9 +31,6 @@ DoTest BP.adh026.dat.save BP.baseref.dat DoTest BPstep.adh026.dat.save BPstep.baseref.dat DoTest Helix.adh026.dat.save Helix.baseref.dat DoTest BPstep.groove.dat.save BPstep.groove.dat -DoTest BP.adh026.dat.save BP.GuessBP.dat -DoTest BPstep.adh026.dat.save BPstep.GuessBP.dat -DoTest Helix.adh026.dat.save Helix.GuessBP.dat DoTest axes.bases.pdb.save axes.bases.pdb DoTest axes.bp.mol2.save axes.bp.mol2 DoTest axes.step.crd.save axes.step.crd From 7f36c6ce9828c6547ba84d28488bd5325e95889b Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 10:34:50 -0400 Subject: [PATCH 29/39] Add specifiedbp example to help text --- src/Action_NAstruct.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 5f094e9b3f..b7cadc8ce2 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -139,6 +139,9 @@ void Action_NAstruct::Help() const { " pairing using geometric criteria in a manner similar to 3DNA.\n" " - If 'allframes' is specified, base pairing will be determined\n" " using geometric criteria for every single frame.\n" + " - If 'specifiedbp' is specified, base pairing is given by subsequent\n" + " 'pairs -,...' arguments, where and are the residue\n" + " numbers of bases in the base pair, e.g. 'pairs 1-16,2-15,3-14,4-13'.\n" " If 'calcnohb' is specified NA parameters will be calculated even if no\n" " hydrogen bonds present between base pairs.\n" " Base pair parameters are written to 'BP.', base pair step parameters\n" From 6eb6de9a45653b9b8ad8bcd400170a67df87d593 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 10:39:40 -0400 Subject: [PATCH 30/39] Test passing traj args to axes output --- test/Test_NAstruct/RunTest.sh | 2 +- test/Test_NAstruct/axes.bases.pdb.save | 48 -------------------------- 2 files changed, 1 insertion(+), 49 deletions(-) diff --git a/test/Test_NAstruct/RunTest.sh b/test/Test_NAstruct/RunTest.sh index 348401efc7..6512775936 100755 --- a/test/Test_NAstruct/RunTest.sh +++ b/test/Test_NAstruct/RunTest.sh @@ -17,7 +17,7 @@ cat > nastruct.in < Date: Tue, 11 Jul 2023 10:54:17 -0400 Subject: [PATCH 31/39] Fix up help, add new keywords. --- doc/cpptraj.lyx | 205 +++++++++++++++++++++++++++++++++++++++- src/Action_NAstruct.cpp | 6 +- 2 files changed, 205 insertions(+), 6 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index eae887b22f..e688b8acff 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -33250,15 +33250,33 @@ nastruct \end_layout \begin_layout LyX-Code - [{ first | reference | ref | refindex <#> | allframes | guessbp}] -\begin_inset Separator latexpar -\end_inset + [axesout [axesoutarg ...] [axesparmout ]] +\end_layout + +\begin_layout LyX-Code + [bpaxesout [bpaxesoutarg ...] [bpaxesparmout ]] +\end_layout + +\begin_layout LyX-Code + [stepaxesout [stepaxesoutarg ...] [stepaxesparmout ]] +\end_layout +\begin_layout LyX-Code + [axisnameo ] [axisnamex ] [axisnamey ] [axisnamez + ] +\end_layout +\begin_layout LyX-Code + [{ first | reference | ref | refindex <#> | \end_layout \begin_layout LyX-Code - [bptype {anti | para} ...] + allframes | +\end_layout + +\begin_layout LyX-Code + specifiedbp pairs -,... + }] \end_layout \begin_deeper @@ -33445,6 +33463,123 @@ literal "true" \end_layout \end_deeper +\begin_layout Description +[axesout +\begin_inset space ~ +\end_inset + +] Trajectory file to write base axes to. +\end_layout + +\begin_deeper +\begin_layout Description +[axesoutarg +\begin_inset space ~ +\end_inset + +] Trajectory argument to pass to base axes trajectory file (can specify + more than once). +\end_layout + +\begin_layout Description +[axesparmout +\begin_inset space ~ +\end_inset + +] Topology file to write base axes pseudo topology to. +\end_layout + +\end_deeper +\begin_layout Description +[bpaxesout +\begin_inset space ~ +\end_inset + +] Trajectory file to write base pair axes to. +\end_layout + +\begin_deeper +\begin_layout Description +[bpaxesoutarg +\begin_inset space ~ +\end_inset + +] Trajectory argument to pass to base pair axes trajectory file (can + specify more than once). +\end_layout + +\begin_layout Description +[bpaxesparmout +\begin_inset space ~ +\end_inset + +] Topology file to write base pair axes pseudo topology to. +\end_layout + +\end_deeper +\begin_layout Description +[stepaxesout +\begin_inset space ~ +\end_inset + +] Trajectory file to write base pair step axes to. +\end_layout + +\begin_deeper +\begin_layout Description +[stepaxesoutarg +\begin_inset space ~ +\end_inset + +] Trajectory argument to pass to base pair step axes trajectory file + (can specify more than once). +\end_layout + +\begin_layout Description +[stepaxesparmout +\begin_inset space ~ +\end_inset + +] Topology file to write base pair step axes pseudo topology to. +\end_layout + +\end_deeper +\begin_layout Description +[axisnameo +\begin_inset space ~ +\end_inset + +] Change name of axis origin pseudo atom (default 'Orig'). +\end_layout + +\begin_layout Description +[axisnamex +\begin_inset space ~ +\end_inset + +] Change name of axis origin pseudo atom (default 'X'). +\end_layout + +\begin_layout Description +[axisnamey +\begin_inset space ~ +\end_inset + +] Change name of axis origin pseudo atom (default 'Y'). +\end_layout + +\begin_layout Description +[axisnamez +\begin_inset space ~ +\end_inset + +] Change name of axis origin pseudo atom (default 'Z'). +\end_layout + +\begin_layout Standard +How to determine base pairing: +\end_layout + \begin_layout Description [first] Use first frame to determine base pairing (default). \end_layout @@ -33481,6 +33616,23 @@ ref [allframes] If specified determine base pairing each frame. \end_layout +\begin_layout Description +[specifiedbp +\begin_inset space ~ +\end_inset + +pairs +\begin_inset space ~ +\end_inset + +-,...] User specified base pairing. + Base pairs are specified in a comma-separated list after the 'pairs' keyword + as -, where and are the residue numbers of bases in the + base pair, e.g. + 'pairs 1-16,2-15,3-14,4-13'. + Can specify 'pairs' multiple times. +\end_layout + \begin_layout Standard DataSets Created: \end_layout @@ -33687,6 +33839,15 @@ Base pairs are determined either once from the first frame or from a reference allframes \series default is specified. + Base pairing can also be specified via the +\series bold +specifiedbp +\series default + and +\series bold +pairs +\series default + keywords. Base pairing is determined first by base reference axis origin distance, then by stagger, then by angle between base Z axes, then finally by hydrogen bonding (at least one hydrogen bond must be present). @@ -33784,6 +33945,42 @@ sscalc . \end_layout +\begin_layout Standard +Base axes, base pair axes, and base pair step axes can be written to trajectory + files using the +\series bold +axesout +\series default +, +\series bold +bpaxesout +\series default +, and +\series bold +stepaxesout +\series default + and related keywords. + The axes are written using 4 points: an origin, and X Y and Z which are + bonded to the origin. + The names of these pseudo atoms can be changed using the +\series bold +axisnameo +\series default +, +\series bold +axisnamex +\series default +, +\series bold +axisnamey +\series default +, and +\series bold +axisnamez +\series default + keywords. +\end_layout + \begin_layout Subsubsection* Custom Nucleic Acid Base References \end_layout diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index b7cadc8ce2..85d70a528b 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -131,8 +131,8 @@ void Action_NAstruct::Help() const { "\t[stepaxesout [stepaxesoutarg ...] [stepaxesparmout ]]\n" "\t[axisnameo ] [axisnamex ] [axisnamey ] [axisnamez ]\n" "\t[{ %s |\n" - "\t allframes |\n" - "\t specifiedbp pairs -,... }]\n", DataSetList::RefArgs); + "\t allframes |\n" + "\t specifiedbp pairs -,... }]\n", DataSetList::RefArgs); mprintf(" Perform nucleic acid structure analysis. Base pairing can be determined\n" " in multiple ways:\n" " - If 'first' (default) or a reference is specified, determine base\n" @@ -142,6 +142,8 @@ void Action_NAstruct::Help() const { " - If 'specifiedbp' is specified, base pairing is given by subsequent\n" " 'pairs -,...' arguments, where and are the residue\n" " numbers of bases in the base pair, e.g. 'pairs 1-16,2-15,3-14,4-13'.\n" + " - If 'reference', 'ref', or 'refindex' is specified, use a reference\n" + " structure to determine base pairing.\n" " If 'calcnohb' is specified NA parameters will be calculated even if no\n" " hydrogen bonds present between base pairs.\n" " Base pair parameters are written to 'BP.', base pair step parameters\n" From 6a7b94b618048f5afbb425479425c51752f27e3f Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 11:26:31 -0400 Subject: [PATCH 32/39] Add wchbonly keyword --- src/Action_NAstruct.cpp | 23 ++++++++++++++++++----- src/Action_NAstruct.h | 1 + 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 85d70a528b..2472edc9e0 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -35,6 +35,7 @@ Action_NAstruct::Action_NAstruct() : skipIfNoHB_(true), spaceBetweenFrames_(true), sscalc_(false), + wc_hb_only_(false), bpout_(0), ssout_(0), stepout_(0), @@ -123,7 +124,7 @@ void Action_NAstruct::Help() const { mprintf("\t[] [resrange ] [sscalc] [naout ]\n" "\t[noheader] [resmap :{A,C,G,T,U} ...] [calcnohb]\n" "\t[noframespaces] [baseref ] ...\n" - "\t[bpmode {3dna|babcock}]\n" + "\t[bpmode {3dna|babcock}] [wchbonly]\n" "\t[hbcut ] [origincut ] [altona | cremer]\n" "\t[zcut ] [zanglecut ] [groovecalc {simple | 3dna}]\n" "\t[axesout [axesoutarg ...] [axesparmout ]]\n" @@ -146,6 +147,8 @@ void Action_NAstruct::Help() const { " structure to determine base pairing.\n" " If 'calcnohb' is specified NA parameters will be calculated even if no\n" " hydrogen bonds present between base pairs.\n" + " If 'wchbonly' is specified only report the number of Watson-Crick-Franklin\n" + " hydrogen bonds in output.\n" " Base pair parameters are written to 'BP.', base pair step parameters\n" " are written to 'BPstep.', and helix parameters are written to\n" " Helix.'.\n" @@ -187,6 +190,7 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int mprinterr("Error: Unrecognized keyword '%s' for 'bpmode'.\n", bpmode.c_str()); return Action::ERR; } + wc_hb_only_ = actionArgs.hasKey("wchbonly"); double hbcut = actionArgs.getKeyDouble("hbcut", -1); if (hbcut > 0) HBdistCut2_ = hbcut * hbcut; @@ -380,6 +384,10 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int } mprintf("\tHydrogen bond cutoff for determining base pairs is %.2f Angstroms.\n", sqrt( HBdistCut2_ ) ); + if (wc_hb_only_) + mprintf("\tOnly reporting total # of Watson-Crick-Franklin hydrogen bonds.\n"); + else + mprintf("\tReporting total # of all hydrogen bonds.\n"); mprintf("\tBase reference axes origin cutoff for determining base pairs is %.2f Angstroms.\n", sqrt( originCut2_ ) ); mprintf("\tBase Z height cutoff (stagger) for determining base pairs is %.2f Angstroms.\n", @@ -618,7 +626,7 @@ Action_NAstruct::HbondType Action_NAstruct::ID_HBtype(NA_Base const& base1, int /** Given two NA_Bases for which IDs have been given and input coords set, * calculate the number of hydrogen bonds between them. */ -// TODO Identify type of base pairing (WC, Hoog., etc) +// TODO Identify type of base pairing (WCF, Hoog., etc) int Action_NAstruct::CalcNumHB(NA_Base const& base1, NA_Base const& base2, int& n_WC) { int Nhbonds = 0; n_WC = 0; @@ -764,23 +772,25 @@ int Action_NAstruct::SpecifiedBasePairing() { // opposite (theta > 90) directions. bool is_antiparallel; double z_theta = b1Axis.Rz().Angle( b2Axis.Rz() ); +# ifdef NASTRUCTDEBUG double z_deviation_from_linear; +# endif if (z_theta > Constants::PIOVER2) { // If theta(Z) > 90 deg. # ifdef NASTRUCTDEBUG mprintf("\t%s is anti-parallel to %s (%g deg)\n", base1.ResName(), base2.ResName(), z_theta * Constants::RADDEG); + z_deviation_from_linear = Constants::PI - z_theta; # endif is_antiparallel = true; - z_deviation_from_linear = Constants::PI - z_theta; // Antiparallel - flip Y and Z axes of complimentary base b2Axis.FlipYZ(); } else { # ifdef NASTRUCTDEBUG mprintf("\t%s is parallel to %s (%g deg)\n", base1.ResName(), base2.ResName(), z_theta * Constants::RADDEG); + z_deviation_from_linear = z_theta; # endif is_antiparallel = false; - z_deviation_from_linear = z_theta; // Parallel - no flip needed if 3dna. // If using Babcock convention, flip X and Y axes. if (bpConvention_ == BP_BABCOCK) @@ -1433,7 +1443,10 @@ int Action_NAstruct::DeterminePairParameters(int frameNum) { BP.opening_->Add(frameNum, &opening); BP.prop_->Add(frameNum, &prop); BP.buckle_->Add(frameNum, &buckle); - BP.hbonds_->Add(frameNum, &(BP.n_wc_hb_)); + if (wc_hb_only_) + BP.hbonds_->Add(frameNum, &(BP.n_wc_hb_)); + else + BP.hbonds_->Add(frameNum, &(BP.nhb_)); static const int ONE = 1; if (BP.nhb_ > 0) BP.isBP_->Add(frameNum, &ONE); diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 575221b84b..6c1fb1b69b 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -204,6 +204,7 @@ class Action_NAstruct: public Action { bool skipIfNoHB_; ///< When true, do not calc parameters when BP not present bool spaceBetweenFrames_; ///< If false do not print spaces between frames in naout bool sscalc_; ///< If true determine params for consecutive bases in strands + bool wc_hb_only_; ///< If true, only report # of WC hydrogen bonds. CpptrajFile* bpout_; ///< Base pair out (BP.). CpptrajFile* ssout_; ///< Single strand out (SS.). CpptrajFile* stepout_; ///< Base pair step out (BPstep.). From f210d190409271641e539f22f4b4e9ede29e44c8 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 11:31:29 -0400 Subject: [PATCH 33/39] Change the keyword to allhb --- src/Action_NAstruct.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 2472edc9e0..e1f4459bae 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -124,7 +124,7 @@ void Action_NAstruct::Help() const { mprintf("\t[] [resrange ] [sscalc] [naout ]\n" "\t[noheader] [resmap :{A,C,G,T,U} ...] [calcnohb]\n" "\t[noframespaces] [baseref ] ...\n" - "\t[bpmode {3dna|babcock}] [wchbonly]\n" + "\t[bpmode {3dna|babcock}] [allhb]\n" "\t[hbcut ] [origincut ] [altona | cremer]\n" "\t[zcut ] [zanglecut ] [groovecalc {simple | 3dna}]\n" "\t[axesout [axesoutarg ...] [axesparmout ]]\n" @@ -147,7 +147,8 @@ void Action_NAstruct::Help() const { " structure to determine base pairing.\n" " If 'calcnohb' is specified NA parameters will be calculated even if no\n" " hydrogen bonds present between base pairs.\n" - " If 'wchbonly' is specified only report the number of Watson-Crick-Franklin\n" + " If 'allhb' is specified report the total number of hydrogen bonds detected\n" + " instead of just the number of Watson-Crick-Franklin hydrogen bonds.\n" " hydrogen bonds in output.\n" " Base pair parameters are written to 'BP.', base pair step parameters\n" " are written to 'BPstep.', and helix parameters are written to\n" @@ -190,7 +191,7 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int mprinterr("Error: Unrecognized keyword '%s' for 'bpmode'.\n", bpmode.c_str()); return Action::ERR; } - wc_hb_only_ = actionArgs.hasKey("wchbonly"); + wc_hb_only_ = !actionArgs.hasKey("allhb"); double hbcut = actionArgs.getKeyDouble("hbcut", -1); if (hbcut > 0) HBdistCut2_ = hbcut * hbcut; From bf499134043741b3ab916c92898cd89b2694ef3f Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 11:32:58 -0400 Subject: [PATCH 34/39] Document allhb keyword. --- doc/cpptraj.lyx | 7 ++++++- src/Action_NAstruct.cpp | 1 - 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index e688b8acff..9843c2e03c 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -33238,7 +33238,7 @@ nastruct \end_layout \begin_layout LyX-Code - [bpmode {3dna|babcock}] + [bpmode {3dna|babcock}] [allhb] \end_layout \begin_layout LyX-Code @@ -33385,6 +33385,11 @@ literal "true" parallel. \end_layout +\begin_layout Description +[allhb] Report the total number of hydrogen bonds detected instead of just + the number of Watson-Crick-Franklin hydrogen bonds. +\end_layout + \begin_layout Description [hbcut \begin_inset space ~ diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index e1f4459bae..9d16c8b9c5 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -149,7 +149,6 @@ void Action_NAstruct::Help() const { " hydrogen bonds present between base pairs.\n" " If 'allhb' is specified report the total number of hydrogen bonds detected\n" " instead of just the number of Watson-Crick-Franklin hydrogen bonds.\n" - " hydrogen bonds in output.\n" " Base pair parameters are written to 'BP.', base pair step parameters\n" " are written to 'BPstep.', and helix parameters are written to\n" " Helix.'.\n" From 372b50ce20b3decfeafde57a0fe9b7fc136cde43 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 11:33:32 -0400 Subject: [PATCH 35/39] 6.19.6. Revision bump for adding 'axesout' etc, 'allhb', and 'specifiedhb' etc keywords to the 'nastruct' action. --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index 183804cea8..45cab90c8e 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.19.5" +#define CPPTRAJ_INTERNAL_VERSION "V6.19.6" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif From 15837498bfa5f7392755e2afe636c4ec6afa871e Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 11:43:02 -0400 Subject: [PATCH 36/39] Protect tests in parallel --- test/Test_NAstruct/RunTest.sh | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/test/Test_NAstruct/RunTest.sh b/test/Test_NAstruct/RunTest.sh index 6512775936..f1905d9d14 100755 --- a/test/Test_NAstruct/RunTest.sh +++ b/test/Test_NAstruct/RunTest.sh @@ -17,8 +17,6 @@ cat > nastruct.in < nastruct.in < nastruct.in < Date: Tue, 11 Jul 2023 11:57:26 -0400 Subject: [PATCH 37/39] Do proper search for residue # in Bases instead of using user specified number as internal index into Bases --- src/Action_NAstruct.cpp | 22 ++++++++++++++++------ src/Action_NAstruct.h | 2 ++ 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 9d16c8b9c5..30da4aad90 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -722,6 +722,16 @@ Action_NAstruct::BPmap::iterator return entry; } +/** Search Bases_ for the specified residue number, return index into Bases_. */ +int Action_NAstruct::find_index_in_bases(int resnum) const { + for (Barray::const_iterator it = Bases_.begin(); it != Bases_.end(); ++it) { + if (it->ResNum() == resnum) { + return (int)(it - Bases_.begin()); + } + } + return -1; +} + /** User-specified base pairing. */ int Action_NAstruct::SpecifiedBasePairing() { int n_wc_hb; @@ -733,15 +743,15 @@ int Action_NAstruct::SpecifiedBasePairing() { ++it) { // User-specified base pair #s start from 1 - int b1idx = it->first - 1; - int b2idx = it->second - 1; - if (b1idx < 0 || (unsigned int)b1idx >= Bases_.size()) { - mprinterr("Error: Specified base # %u is out of range.\n", it->first); + int b1idx = find_index_in_bases(it->first - 1); + int b2idx = find_index_in_bases(it->second - 1); + if (b1idx < 0) { + mprinterr("Error: Specified base residue # %u not found in set up bases.\n", it->first); return 1; } NA_Base& base1 = Bases_[b1idx]; - if (b2idx < 0 || (unsigned int)b2idx >= Bases_.size()) { - mprinterr("Error: Specified base # %u is out of range.\n", it->second); + if (b2idx < 0) { + mprinterr("Error: Specified base residue # %u not found in set up bases.\n", it->second); return 1; } NA_Base& base2 = Bases_[b2idx]; diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 6c1fb1b69b..49c511bd26 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -156,6 +156,8 @@ class Action_NAstruct: public Action { BPmap::iterator AddBasePair(int, NA_Base const&, int, NA_Base const&); /// Determine which bases are paired geometrically, set base pair data. int DetermineBasePairing(); + /// Find index in bases for given internal residue # + int find_index_in_bases(int) const; /// Set up base pairs based on user specification int SpecifiedBasePairing(); /// Calculate translational/rotational parameters between two axes. From f23698cfef53cd668ade98b2e7fc42a45cf415bc Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 12:00:52 -0400 Subject: [PATCH 38/39] Remove old debug print of axes, superseded by new functionality. --- src/Action_NAstruct.cpp | 21 +-------------------- 1 file changed, 1 insertion(+), 20 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 30da4aad90..76cdbdeb64 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -511,8 +511,6 @@ int Action_NAstruct::SetupBaseAxes(Frame const& InputFrame) { Frame refFrame(maxResSize_); // Hold copy of base reference coords for RMS fit Frame inpFrame(maxResSize_); // Hold copy of input base coords for RMS fit # ifdef NASTRUCTDEBUG - PDBfile baseaxesfile; - baseaxesfile.OpenWrite("baseaxes.pdb"); PDBfile basesfile; basesfile.OpenWrite("bases.pdb"); mprintf("\n=================== Setup Base Axes ===================\n"); @@ -578,12 +576,7 @@ int Action_NAstruct::SetupBaseAxes(Frame const& InputFrame) { } # endif } // END loop over bases -# ifdef NASTRUCTDEBUG - // DEBUG - Write base axis to file - for (std::vector::iterator base = Bases_.begin(); - base != Bases_.end(); ++base) - WriteAxes(baseaxesfile, base->ResNum()+1, base->ResName(), base->Axis()); -# endif + return 0; } @@ -1368,8 +1361,6 @@ int Action_NAstruct::axis_points_5p_to_3p(NA_Base const& base1) const { int Action_NAstruct::DeterminePairParameters(int frameNum) { double Param[6]; # ifdef NASTRUCTDEBUG - PDBfile basepairaxesfile; - basepairaxesfile.OpenWrite("basepairaxes.pdb"); mprintf("\n=================== Determine BP Parameters ===================\n"); # endif // NOTE: iterator cannot be const because bpaxis_ needs to be updated @@ -1469,10 +1460,6 @@ int Action_NAstruct::DeterminePairParameters(int frameNum) { bp_axes_vec[4] = BP.bpaxis_.Oxyz()[1]; bp_axes_vec[5] = BP.bpaxis_.Oxyz()[2]; BP.axes_nxyz_->Add(frameNum, bp_axes_vec); -# endif -# ifdef NASTRUCTDEBUG - // DEBUG - write base pair axes - WriteAxes(basepairaxesfile, b1+1, base1.ResName(), BP.bpaxis_); # endif } // Calculate base parameters. @@ -1489,8 +1476,6 @@ int Action_NAstruct::DeterminePairParameters(int frameNum) { int Action_NAstruct::DetermineStepParameters(int frameNum) { double Param[6]; # ifdef NASTRUCTDEBUG - PDBfile stepaxesfile; - stepaxesfile.OpenWrite("stepaxes.pdb"); mprintf("\n=================== Determine BPstep Parameters ===================\n"); # endif if (BasePairs_.size() < 2) return 0; @@ -1674,10 +1659,6 @@ int Action_NAstruct::DetermineStepParameters(int frameNum) { currentStep.incl_->Add(frameNum, &incl); currentStep.tip_->Add(frameNum, &tip); currentStep.htwist_->Add(frameNum, &htwist); -# ifdef NASTRUCTDEBUG - // DEBUG - write base pair step axes - WriteAxes(stepaxesfile, base1.ResNum()+1, base1.ResName(), midFrame); -# endif } // END second base pair found } // END second base pair valid } // END loop over base pairs From 14ba457dbdb6eb747aee04721cb8aeaeddc10069 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Jul 2023 12:29:40 -0400 Subject: [PATCH 39/39] Fix error in printf statement --- src/Action_NAstruct.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 76cdbdeb64..881b7e1529 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -449,7 +449,7 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int } if (axesOut_ != 0 || bpAxesOut_ != 0 || stepAxesOut_ != 0) { mprintf("\tAxes pseudo atom names: origin='%s' X='%s' Y='%s' Z='%s'\n", - *axisNameO_, *axisNameX_, *axisNameY_, axisNameZ_); + *axisNameO_, *axisNameX_, *axisNameY_, *axisNameZ_); } mprintf("# Citations: Babcock MS; Pednault EPD; Olson WK; \"Nucleic Acid Structure\n" "# Analysis: Mathematics for Local Cartesian and Helical Structure\n"