From 2f58addb3a20c3d05b72fcbadd114405db9cf546 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 09:49:57 -0500 Subject: [PATCH 01/13] Test adding sscalc to nastruct --- src/Action_NAstruct.cpp | 54 ++++++++++++++++++++++++++++++++++++++++- src/Action_NAstruct.h | 3 +++ 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 5e4af0a7bb..9c575ac19f 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -37,7 +37,7 @@ Action_NAstruct::Action_NAstruct() : {} void Action_NAstruct::Help() const { - mprintf("\t[] [resrange ] [naout ]\n" + mprintf("\t[] [resrange ] [naout ] [sscalc]\n" "\t[noheader] [resmap :{A,C,G,T,U} ...] [calcnohb]\n" "\t[noframespaces] [baseref ] ...\n" "\t[hbcut ] [origincut ] [altona | cremer]\n" @@ -127,6 +127,7 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int findBPmode_ = FIRST; else findBPmode_ = FIRST; + sscalc_ = actionArgs.hasKey("sscalc"); # ifdef MPI if (findBPmode_ == ALL && trajComm_.Size() > 1) { mprinterr("Error: Currently 'allframes' does not work with > 1 process per trajectory" @@ -195,6 +196,8 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int mprintf("Scanning all NA residues\n"); else mprintf("Scanning residues %s\n",resRange_.RangeArg()); + if (sscalc_) + mprintf("\tWill determine parameters for consecutive bases in strands.\n"); if (bpout_ != 0) { mprintf("\tBase pair parameters written to %s\n", bpout_->Filename().full()); mprintf("\tBase pair step parameters written to %s\n", stepout_->Filename().full()); @@ -964,6 +967,52 @@ int Action_NAstruct::GetBaseIdxStep(int idx, int Nsteps) const { return GetBaseIdxStep( base.C5resIdx(), Nsteps + 1); } +/** For bases in a single strand, get the values of buckle, propeller twist, + * opening, shear, stretch, and stagger. + */ +int Action_NAstruct::DetermineStrandParameters(int frameNum) { + double Param[6]; + NA_Axis commonAxis; + // Loop over strands + for (StrandArray::const_iterator strand = Strands_.begin(); + strand != Strands_.end(); ++strand) + { + mprintf("DEBUG: Strand %li, %i-%i\n", strand-Strands_.begin(), strand->first, strand->second); + for (int b1idx = strand->first; b1idx < strand->second; b1idx++) + { + int b2idx = b1idx + 1; + mprintf("DEBUG:\tStrand pair: %i to %i\n", b1idx, b2idx); + NA_Base& base1 = Bases_[b1idx]; + NA_Base& base2 = Bases_[b2idx]; //TODO copy? + // Calc parameters between bases in the strand + calculateParameters(base2.Axis(), base1.Axis(), &commonAxis, Param); + // Store data + Param[3] *= Constants::RADDEG; + Param[4] *= Constants::RADDEG; + Param[5] *= Constants::RADDEG; + //mprintf("DBG: BP %i # hbonds = %i\n", nbasepair+1, NumberOfHbonds_[nbasepair]); + // Convert everything to float to save space + float shear = (float)Param[0]; + float stretch = (float)Param[1]; + float stagger = (float)Param[2]; + float opening = (float)Param[3]; + float prop = (float)Param[4]; + float buckle = (float)Param[5]; + // Add to DataSets + mprintf("DEBUG:\tShear= %f stretch= %f stagger= %f\n", shear, stretch, stagger); + mprintf("DEBUG:\tOpeni= %f propell= %f buckle_= %f\n", opening, prop, buckle); + //BP.shear_->Add(frameNum, &shear); + //BP.stretch_->Add(frameNum, &stretch); + //BP.stagger_->Add(frameNum, &stagger); + //BP.opening_->Add(frameNum, &opening); + //BP.prop_->Add(frameNum, &prop); + //BP.buckle_->Add(frameNum, &buckle); + //BP.hbonds_->Add(frameNum, &(BP.n_wc_hb_)); + } + } + return 0; +} + // Action_NAstruct::DeterminePairParameters() /** For each base pair, get the values of buckle, propeller twist, * opening, shear, stretch, and stagger. Also determine the origin and @@ -1476,6 +1525,9 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { # endif findBPmode_ = REFERENCE; } + // Determine strand parameters if desired + if (sscalc_) + DetermineStrandParameters(frameNum); // Determine base parameters DeterminePairParameters(frameNum); diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 0dc800f4b4..00c6d4bb47 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -124,6 +124,8 @@ class Action_NAstruct: public Action { int helicalParameters(NA_Axis const&, NA_Axis const&, double *); /// \return index of base that is N steps away in specified direction from another base. int GetBaseIdxStep(int, int) const; + /// Determine individual base parameters in single strands. + int DetermineStrandParameters(int); /// Determine individual base and base pair parameters. int DeterminePairParameters(int); /// Determine base pair steps and step parameters, including HC groove calc. @@ -156,6 +158,7 @@ class Action_NAstruct: public Action { bool seriesUpdated_; ///< If false, check that time series data is nframes long 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 CpptrajFile* bpout_; ///< Base pair out (BP.). CpptrajFile* stepout_; ///< Base pair step out (BPstep.). CpptrajFile* helixout_; ///< Helical parameters out (Helix.). From 20ba49132a14511bc676491c0ff1055216f8d00f Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 10:33:48 -0500 Subject: [PATCH 02/13] Add data sets --- src/Action_NAstruct.cpp | 42 ++++++++++++++++++++++++++++++++++------- src/Action_NAstruct.h | 14 ++++++++++++++ 2 files changed, 49 insertions(+), 7 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 9c575ac19f..59579f01c4 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -982,6 +982,35 @@ int Action_NAstruct::DetermineStrandParameters(int frameNum) { { int b2idx = b1idx + 1; mprintf("DEBUG:\tStrand pair: %i to %i\n", b1idx, b2idx); + + // Get/add data set for strandpair + Rpair strandpair(b1idx, b2idx); + Smap::iterator entry = StrandPairs_.find( strandpair ); + if (entry == StrandPairs_.end()) { + // New strand pair + Stype SP; + // MetaData md = NewStepType(BS, BP1.base1idx_, BP1.base2idx_, + // BP2.base1idx_, BP2.base2idx_, Steps_.size()+1); + MetaData md(dataname_, StrandPairs_.size()+1); + md.SetLegend( Bases_[b1idx].BaseName() + "-" + Bases_[b2idx].BaseName() ); + md.SetAspect("dx"); + SP.dx_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("dy"); + SP.dy_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("dz"); + SP.dz_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("rx"); + SP.rx_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("ry"); + SP.ry_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("rz"); + SP.rz_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + + entry = StrandPairs_.insert( entry, std::pair(strandpair, SP) ); // FIXME does entry make more efficient? + } + Stype& currentSpair = entry->second; + + // Get bases NA_Base& base1 = Bases_[b1idx]; NA_Base& base2 = Bases_[b2idx]; //TODO copy? // Calc parameters between bases in the strand @@ -1001,13 +1030,12 @@ int Action_NAstruct::DetermineStrandParameters(int frameNum) { // Add to DataSets mprintf("DEBUG:\tShear= %f stretch= %f stagger= %f\n", shear, stretch, stagger); mprintf("DEBUG:\tOpeni= %f propell= %f buckle_= %f\n", opening, prop, buckle); - //BP.shear_->Add(frameNum, &shear); - //BP.stretch_->Add(frameNum, &stretch); - //BP.stagger_->Add(frameNum, &stagger); - //BP.opening_->Add(frameNum, &opening); - //BP.prop_->Add(frameNum, &prop); - //BP.buckle_->Add(frameNum, &buckle); - //BP.hbonds_->Add(frameNum, &(BP.n_wc_hb_)); + currentSpair.dx_->Add(frameNum, &shear); + currentSpair.dy_->Add(frameNum, &stretch); + currentSpair.dz_->Add(frameNum, &stagger); + currentSpair.rx_->Add(frameNum, &opening); + currentSpair.ry_->Add(frameNum, &prop); + currentSpair.rz_->Add(frameNum, &buckle); } } return 0; diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 00c6d4bb47..ca72ced77c 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -43,6 +43,18 @@ class Action_NAstruct: public Action { /// How to find base pairs: first frame, reference structure, all frames, guess. enum FindType { FIRST = 0, REFERENCE, ALL, GUESS }; // ----- Data Structures --------------------- + /// Hold consecutive bases + struct Stype { + DataSet_1D* dx_; + DataSet_1D* dy_; + DataSet_1D* dz_; + DataSet_1D* rx_; + DataSet_1D* ry_; + DataSet_1D* rz_; + unsigned int strandidx_; ///< Index in Strands_ + unsigned int base1idx_; ///< Index of first base in Bases_ + unsigned int base2idx_; ///< Index of second base in Bases_ + }; /// Hold a base pair. struct BPtype { NA_Axis bpaxis_; ///< Base pair axis. @@ -95,6 +107,7 @@ class Action_NAstruct: public Action { // ----- Type Definitions -------------------- typedef std::vector Barray; ///< Array of NA bases typedef std::pair Rpair; ///< Pair of residue numbers / BP indices + typedef std::map Smap; ///< Map of residue numbers to strand pairs typedef std::map BPmap; ///< Map of residue numbers to BP typedef std::map StepMap; ///< Map of BP indices to Steps typedef std::vector StrandArray; ///< Hold indices into Bases_ for strand beg/end @@ -140,6 +153,7 @@ class Action_NAstruct: public Action { NA_Reference refBases_; ///< Hold reference bases RefMapType nameToRef_; ///< Map residue names to custom reference Barray Bases_; ///< Hold nucleobases + Smap StrandPairs_; ///< Hold consecutive bases in strands BPmap BasePairs_; ///< Hold base pairs StepMap Steps_; ///< Hold base pair steps. StrandArray Strands_; ///< Hold strand info From 7f602e95e60364bd782d46f3f708e5a4bdcbac13 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 10:48:17 -0500 Subject: [PATCH 03/13] Strand sets can be allocated during setup --- src/Action_NAstruct.cpp | 139 ++++++++++++++++++++++------------------ src/Action_NAstruct.h | 4 +- 2 files changed, 78 insertions(+), 65 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 59579f01c4..cf08d615ff 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -973,70 +973,40 @@ int Action_NAstruct::GetBaseIdxStep(int idx, int Nsteps) const { int Action_NAstruct::DetermineStrandParameters(int frameNum) { double Param[6]; NA_Axis commonAxis; - // Loop over strands - for (StrandArray::const_iterator strand = Strands_.begin(); - strand != Strands_.end(); ++strand) + // Loop over strand pairs + for (Smap::const_iterator it = StrandPairs_.begin(); it != StrandPairs_.end(); ++it) { - mprintf("DEBUG: Strand %li, %i-%i\n", strand-Strands_.begin(), strand->first, strand->second); - for (int b1idx = strand->first; b1idx < strand->second; b1idx++) - { - int b2idx = b1idx + 1; - mprintf("DEBUG:\tStrand pair: %i to %i\n", b1idx, b2idx); - - // Get/add data set for strandpair - Rpair strandpair(b1idx, b2idx); - Smap::iterator entry = StrandPairs_.find( strandpair ); - if (entry == StrandPairs_.end()) { - // New strand pair - Stype SP; - // MetaData md = NewStepType(BS, BP1.base1idx_, BP1.base2idx_, - // BP2.base1idx_, BP2.base2idx_, Steps_.size()+1); - MetaData md(dataname_, StrandPairs_.size()+1); - md.SetLegend( Bases_[b1idx].BaseName() + "-" + Bases_[b2idx].BaseName() ); - md.SetAspect("dx"); - SP.dx_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); - md.SetAspect("dy"); - SP.dy_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); - md.SetAspect("dz"); - SP.dz_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); - md.SetAspect("rx"); - SP.rx_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); - md.SetAspect("ry"); - SP.ry_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); - md.SetAspect("rz"); - SP.rz_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); - - entry = StrandPairs_.insert( entry, std::pair(strandpair, SP) ); // FIXME does entry make more efficient? - } - Stype& currentSpair = entry->second; - - // Get bases - NA_Base& base1 = Bases_[b1idx]; - NA_Base& base2 = Bases_[b2idx]; //TODO copy? - // Calc parameters between bases in the strand - calculateParameters(base2.Axis(), base1.Axis(), &commonAxis, Param); - // Store data - Param[3] *= Constants::RADDEG; - Param[4] *= Constants::RADDEG; - Param[5] *= Constants::RADDEG; - //mprintf("DBG: BP %i # hbonds = %i\n", nbasepair+1, NumberOfHbonds_[nbasepair]); - // Convert everything to float to save space - float shear = (float)Param[0]; - float stretch = (float)Param[1]; - float stagger = (float)Param[2]; - float opening = (float)Param[3]; - float prop = (float)Param[4]; - float buckle = (float)Param[5]; - // Add to DataSets - mprintf("DEBUG:\tShear= %f stretch= %f stagger= %f\n", shear, stretch, stagger); - mprintf("DEBUG:\tOpeni= %f propell= %f buckle_= %f\n", opening, prop, buckle); - currentSpair.dx_->Add(frameNum, &shear); - currentSpair.dy_->Add(frameNum, &stretch); - currentSpair.dz_->Add(frameNum, &stagger); - currentSpair.rx_->Add(frameNum, &opening); - currentSpair.ry_->Add(frameNum, &prop); - currentSpair.rz_->Add(frameNum, &buckle); - } + int b1idx = it->first.first; + int b2idx = it->first.second; + Stype const& SP = it->second; + 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? + // Calc parameters between bases in the strand + calculateParameters(base2.Axis(), base1.Axis(), &commonAxis, Param); + // Store data + Param[3] *= Constants::RADDEG; + Param[4] *= Constants::RADDEG; + Param[5] *= Constants::RADDEG; + //mprintf("DBG: BP %i # hbonds = %i\n", nbasepair+1, NumberOfHbonds_[nbasepair]); + // Convert everything to float to save space + float shear = (float)Param[0]; + float stretch = (float)Param[1]; + float stagger = (float)Param[2]; + float opening = (float)Param[3]; + float prop = (float)Param[4]; + float buckle = (float)Param[5]; + // Add to DataSets + mprintf("DEBUG:\tShear= %f stretch= %f stagger= %f\n", shear, stretch, stagger); + mprintf("DEBUG:\tOpeni= %f propell= %f buckle_= %f\n", opening, prop, buckle); + SP.dx_->Add(frameNum, &shear); + SP.dy_->Add(frameNum, &stretch); + SP.dz_->Add(frameNum, &stagger); + SP.rx_->Add(frameNum, &opening); + SP.ry_->Add(frameNum, &prop); + SP.rz_->Add(frameNum, &buckle); } return 0; } @@ -1496,6 +1466,49 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { if (GuessBasePairing(setup.Top())) return Action::ERR; findBPmode_ = REFERENCE; } + // Determine strand data + if (sscalc_) { + for (StrandArray::const_iterator strand = Strands_.begin(); + strand != Strands_.end(); ++strand) + { + mprintf("DEBUG: Strand %li, %i-%i\n", strand-Strands_.begin(), strand->first, strand->second); + for (int b1idx = strand->first; b1idx < strand->second; b1idx++) + { + int b2idx = b1idx + 1; + mprintf("DEBUG:\tStrand pair: %i to %i\n", b1idx, b2idx); + + // Get/add data set for strandpair + Rpair strandpair(b1idx, b2idx); + Smap::iterator entry = StrandPairs_.find( strandpair ); + if (entry == StrandPairs_.end()) { + // New strand pair + Stype SP; + // MetaData md = NewStepType(BS, BP1.base1idx_, BP1.base2idx_, + // BP2.base1idx_, BP2.base2idx_, Steps_.size()+1); + MetaData md(dataname_, StrandPairs_.size()+1); + md.SetLegend( Bases_[b1idx].BaseName() + "-" + Bases_[b2idx].BaseName() ); + md.SetAspect("dx"); + SP.dx_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("dy"); + SP.dy_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("dz"); + SP.dz_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("rx"); + SP.rx_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("ry"); + SP.ry_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + md.SetAspect("rz"); + SP.rz_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); + + SP.strandidx_ = (unsigned int)(strand - Strands_.begin()); + //SP.base1idx_ = b1idx_; + //SP.base2idx_ = b2idx_; + + entry = StrandPairs_.insert( entry, std::pair(strandpair, SP) ); // FIXME does entry make more efficient? + } + } // END loop over bases in strand + } // END loop over strands + } // END if sscalc return Action::OK; } diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index ca72ced77c..032705c173 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -52,8 +52,8 @@ class Action_NAstruct: public Action { DataSet_1D* ry_; DataSet_1D* rz_; unsigned int strandidx_; ///< Index in Strands_ - unsigned int base1idx_; ///< Index of first base in Bases_ - unsigned int base2idx_; ///< Index of second base in Bases_ + //unsigned int base1idx_; ///< Index of first base in Bases_ + //unsigned int base2idx_; ///< Index of second base in Bases_ }; /// Hold a base pair. struct BPtype { From 73aa8cc77d67147f58d06998bf864faae7eb99f8 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 10:52:10 -0500 Subject: [PATCH 04/13] Ensure strand time series are updated --- src/Action_NAstruct.cpp | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index cf08d615ff..e2ff0989b5 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -1829,8 +1829,8 @@ void Action_NAstruct::UpdateSeries() { if (seriesUpdated_) return; if (nframes_ > 0) { // Base pair data - for (BPmap::iterator it = BasePairs_.begin(); it != BasePairs_.end(); ++it) { - BPtype& BP = it->second; + for (BPmap::const_iterator it = BasePairs_.begin(); it != BasePairs_.end(); ++it) { + BPtype const& BP = it->second; UpdateTimeSeries( nframes_, BP.shear_ ); UpdateTimeSeries( nframes_, BP.stretch_ ); UpdateTimeSeries( nframes_, BP.stagger_ ); @@ -1845,8 +1845,8 @@ void Action_NAstruct::UpdateSeries() { UpdateTimeSeries( nframes_, Bases_[BP.base2idx_].Pucker() ); } // Step and helix data - for (StepMap::iterator it = Steps_.begin(); it != Steps_.end(); ++it) { - StepType& BS = it->second; + for (StepMap::const_iterator it = Steps_.begin(); it != Steps_.end(); ++it) { + StepType const& BS = it->second; UpdateTimeSeries( nframes_, BS.shift_ ); UpdateTimeSeries( nframes_, BS.slide_ ); UpdateTimeSeries( nframes_, BS.rise_ ); @@ -1863,6 +1863,16 @@ void Action_NAstruct::UpdateSeries() { UpdateTimeSeries( nframes_, BS.majGroove_ ); UpdateTimeSeries( nframes_, BS.minGroove_ ); } + // Strand data + for (Smap::const_iterator it = StrandPairs_.begin(); it != StrandPairs_.end(); ++it) { + Stype const& SP = it->second; + UpdateTimeSeries( nframes_, SP.dx_ ); + UpdateTimeSeries( nframes_, SP.dy_ ); + UpdateTimeSeries( nframes_, SP.dz_ ); + UpdateTimeSeries( nframes_, SP.rx_ ); + UpdateTimeSeries( nframes_, SP.ry_ ); + UpdateTimeSeries( nframes_, SP.rz_ ); + } } // Should only be called once. seriesUpdated_ = true; From 6dcd1ed6b08836440bfb4c13f9c662eceef1c558 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 11:06:32 -0500 Subject: [PATCH 05/13] Add file output for sscalc --- src/Action_NAstruct.cpp | 52 +++++++++++++++++++++++++++++++++++++---- src/Action_NAstruct.h | 1 + 2 files changed, 49 insertions(+), 4 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index e2ff0989b5..2e88c88a37 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -29,7 +29,11 @@ Action_NAstruct::Action_NAstruct() : seriesUpdated_(false), skipIfNoHB_(true), spaceBetweenFrames_(true), - bpout_(0), stepout_(0), helixout_(0), + sscalc_(false), + bpout_(0), + ssout_(0), + stepout_(0), + helixout_(0), masterDSL_(0) # ifdef NASTRUCTDEBUG ,calcparam_(true) @@ -73,6 +77,7 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int debug_ = debugIn; masterDSL_ = init.DslPtr(); // Get keywords + sscalc_ = actionArgs.hasKey("sscalc"); std::string outputsuffix = actionArgs.GetStringKey("naout"); if (!outputsuffix.empty()) { // Set up output files. @@ -81,6 +86,10 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int stepout_ = init.DFL().AddCpptrajFile(FName.PrependFileName("BPstep."), "Base Pair Step"); helixout_ = init.DFL().AddCpptrajFile(FName.PrependFileName("Helix."), "Helix"); if (bpout_ == 0 || stepout_ == 0 || helixout_ == 0) return Action::ERR; + if (sscalc_) { + ssout_ = init.DFL().AddCpptrajFile(FName.PrependFileName("SS."), "Single Strand"); + if (ssout_ == 0) return Action::ERR; + } } double hbcut = actionArgs.getKeyDouble("hbcut", -1); if (hbcut > 0) @@ -127,7 +136,6 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int findBPmode_ = FIRST; else findBPmode_ = FIRST; - sscalc_ = actionArgs.hasKey("sscalc"); # ifdef MPI if (findBPmode_ == ALL && trajComm_.Size() > 1) { mprinterr("Error: Currently 'allframes' does not work with > 1 process per trajectory" @@ -196,8 +204,6 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int mprintf("Scanning all NA residues\n"); else mprintf("Scanning residues %s\n",resRange_.RangeArg()); - if (sscalc_) - mprintf("\tWill determine parameters for consecutive bases in strands.\n"); if (bpout_ != 0) { mprintf("\tBase pair parameters written to %s\n", bpout_->Filename().full()); mprintf("\tBase pair step parameters written to %s\n", stepout_->Filename().full()); @@ -205,6 +211,11 @@ Action::RetType Action_NAstruct::Init(ArgList& actionArgs, ActionInit& init, int if (!printheader_) mprintf("\tHeader line will not be written.\n"); if (!spaceBetweenFrames_) mprintf("\tNo spaces will be written between frames.\n"); } + if (sscalc_) { + mprintf("\tWill determine parameters for consecutive bases in strands.\n"); + if (ssout_ != 0) + mprintf("\tSingle strand parameters written to %s\n", ssout_->Filename().full()); + } mprintf("\tHydrogen bond cutoff for determining base pairs is %.2f Angstroms.\n", sqrt( HBdistCut2_ ) ); if (findBPmode_ != GUESS) { @@ -1023,6 +1034,7 @@ int Action_NAstruct::DeterminePairParameters(int frameNum) { basepairaxesfile.OpenWrite("basepairaxes.pdb"); mprintf("\n=================== Determine BP Parameters ===================\n"); # endif + // NOTE: iterator cannot be const because bpaxis_ needs to be updated for (BPmap::iterator it = BasePairs_.begin(); it != BasePairs_.end(); ++it) { if (it->second.nhb_ < 1 && skipIfNoHB_) continue; @@ -1880,6 +1892,7 @@ void Action_NAstruct::UpdateSeries() { // Output Format Strings static const char* BP_OUTPUT_FMT = "%8i %8i %8i %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %2.0f %2.0f"; +static const char* SS_OUTPUT_FMT = "%8i %8i %8i %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f"; static const char* GROOVE_FMT = " %10.4f %10.4f"; static const char* HELIX_OUTPUT_FMT = "%8i %4i-%-4i %4i-%-4i %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f"; static const char* STEP_OUTPUT_FMT = "%8i %4i-%-4i %4i-%-4i %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f"; @@ -1927,6 +1940,37 @@ void Action_NAstruct::Print() { } } + // ---------- Single strand parameters --------- + if (sscalc_ && ssout_ != 0) { + // Check that there is actually data + if ( StrandPairs_.empty() || nframes_ < 1) + mprintf("Warning: Could not write single strand file %s: No SS data.\n", + ssout_->Filename().full()); + else { + mprintf("\tSingle strand output file %s; %i frames, %zu base pairs.\n", + ssout_->Filename().full(), nframes_, StrandPairs_.size()); + // File header + if (printheader_) { + bpout_->Printf("%-8s %8s %8s %10s %10s %10s %10s %10s %10s\n", + "#Frame","Base1","Base2", "DX","DY","DZ", "RX","RY","RZ"); + } + // Loop over all frames + for (int frame = 0; frame < nframes_; ++frame) { + for (Smap::const_iterator it = StrandPairs_.begin(); + it != StrandPairs_.end(); ++it) + { + Stype const& SP = it->second; + ssout_->Printf(SS_OUTPUT_FMT, frame+1, + Bases_[it->first.first].ResNum()+1, Bases_[it->first.second].ResNum()+1, + SP.dx_->Dval(frame), SP.dy_->Dval(frame), SP.dz_->Dval(frame), + SP.rx_->Dval(frame), SP.ry_->Dval(frame), SP.rz_->Dval(frame)); + ssout_->Printf("\n"); + } + if (spaceBetweenFrames_) ssout_->Printf("\n"); + } + } + } + // ---------- Base pair step parameters -------- // Check that there is actually data if ( Steps_.empty() || nframes_ < 1 ) { diff --git a/src/Action_NAstruct.h b/src/Action_NAstruct.h index 032705c173..468bfb1c91 100644 --- a/src/Action_NAstruct.h +++ b/src/Action_NAstruct.h @@ -174,6 +174,7 @@ class Action_NAstruct: public Action { bool spaceBetweenFrames_; ///< If false do not print spaces between frames in naout bool sscalc_; ///< If true determine params for consecutive bases in strands CpptrajFile* bpout_; ///< Base pair out (BP.). + CpptrajFile* ssout_; ///< Single strand out (SS.). CpptrajFile* stepout_; ///< Base pair step out (BPstep.). CpptrajFile* helixout_; ///< Helical parameters out (Helix.). std::string dataname_; ///< NA DataSet name (default NA). From a12c730997cd9aa13bff492295ccfa3663f23e4e Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 11:08:46 -0500 Subject: [PATCH 06/13] Hide some debug info --- src/Action_NAstruct.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index 2e88c88a37..107d049725 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -1010,8 +1010,8 @@ int Action_NAstruct::DetermineStrandParameters(int frameNum) { float prop = (float)Param[4]; float buckle = (float)Param[5]; // Add to DataSets - mprintf("DEBUG:\tShear= %f stretch= %f stagger= %f\n", shear, stretch, stagger); - mprintf("DEBUG:\tOpeni= %f propell= %f buckle_= %f\n", opening, prop, buckle); + //mprintf("DEBUG:\tShear= %f stretch= %f stagger= %f\n", shear, stretch, stagger); + //mprintf("DEBUG:\tOpeni= %f propell= %f buckle_= %f\n", opening, prop, buckle); SP.dx_->Add(frameNum, &shear); SP.dy_->Add(frameNum, &stretch); SP.dz_->Add(frameNum, &stagger); @@ -1483,11 +1483,13 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { for (StrandArray::const_iterator strand = Strands_.begin(); strand != Strands_.end(); ++strand) { - mprintf("DEBUG: Strand %li, %i-%i\n", strand-Strands_.begin(), strand->first, strand->second); + if (debug_ > 0) + mprintf("DEBUG: Strand %li, %i-%i\n", strand-Strands_.begin(), strand->first, strand->second); for (int b1idx = strand->first; b1idx < strand->second; b1idx++) { int b2idx = b1idx + 1; - mprintf("DEBUG:\tStrand pair: %i to %i\n", b1idx, b2idx); + if (debug_ > 0) + mprintf("DEBUG:\tStrand pair: %i to %i\n", b1idx, b2idx); // Get/add data set for strandpair Rpair strandpair(b1idx, b2idx); From 4686005ac16b3e2602f5a34f0f35461b1f3402e7 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 11:10:08 -0500 Subject: [PATCH 07/13] Hide more debug info --- 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 107d049725..901cdfcf2a 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -990,7 +990,7 @@ int Action_NAstruct::DetermineStrandParameters(int frameNum) { int b1idx = it->first.first; int b2idx = it->first.second; Stype const& SP = it->second; - mprintf("DEBUG: Strand %u, bases %u to %u\n", SP.strandidx_, b1idx, b2idx); + //mprintf("DEBUG: Strand %u, bases %u to %u\n", SP.strandidx_, b1idx, b2idx); // Get bases NA_Base& base1 = Bases_[b1idx]; From 81e839ae14bbfb594aac5b0e3f7208f4045f5e38 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 11:13:05 -0500 Subject: [PATCH 08/13] Add single strand test --- test/Test_NAstruct/RunTest.sh | 12 +++++++++++- test/Test_NAstruct/SS.mol1.dat.save | 24 ++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 1 deletion(-) create mode 100644 test/Test_NAstruct/SS.mol1.dat.save diff --git a/test/Test_NAstruct/RunTest.sh b/test/Test_NAstruct/RunTest.sh index e21bef5a97..f15912b847 100755 --- a/test/Test_NAstruct/RunTest.sh +++ b/test/Test_NAstruct/RunTest.sh @@ -4,7 +4,7 @@ # Clean CleanFiles nastruct.in BP.*.dat BPstep.*.dat bases.pdb baseaxes.pdb basepairaxes.pdb \ - Helix.*.dat Param.pdb + Helix.*.dat Param.pdb SS.mol1.dat # Test 2 TESTNAME='NAstruct tests' @@ -31,6 +31,16 @@ DoTest BP.adh026.dat.save BP.GuessBP.dat DoTest BPstep.adh026.dat.save BPstep.GuessBP.dat DoTest Helix.adh026.dat.save Helix.GuessBP.dat +# Single strand +cat > nastruct.in < Date: Thu, 11 Feb 2021 11:23:12 -0500 Subject: [PATCH 09/13] Add manual entry and update help --- doc/cpptraj.lyx | 51 +++++++++++++++++++++++++++++++++++++++-- src/Action_NAstruct.cpp | 2 +- 2 files changed, 50 insertions(+), 3 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 13446a031e..28a0a4e455 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -29579,7 +29579,7 @@ nastruct \end_inset - [] [resrange ] [naout ] + [] [resrange ] [sscalc] [naout ] \end_layout \begin_layout LyX-Code @@ -29627,6 +29627,10 @@ name>] Output data set name. ] Residue range to search for nucleic acids in (default all). \end_layout +\begin_layout Description +[sscalc] Calculate parameters between consectuive bases in strands. +\end_layout + \begin_layout Description [naout \begin_inset space ~ @@ -29635,6 +29639,11 @@ name>] Output data set name. ] File name suffix for output files; BP. for base pair parameters , BPstep. for base pair step parameters, and Helix. for base pair step helical parameters. + If +\series bold +sscalc +\series default + is specified, also SS. for parameters of consecutive bases in strands. \end_layout \begin_layout Description @@ -29929,11 +29938,39 @@ Helical steps: [htwist]:X Helical twist. \end_layout +\begin_layout Standard +Strands (sscalc only): +\end_layout + +\begin_layout Description +[dx]:X Strand pair X (starting from 1) X displacement. +\end_layout + +\begin_layout Description +[dy]:X Y displacement. +\end_layout + +\begin_layout Description +[dz]:X Z displacement. +\end_layout + +\begin_layout Description +[rx]:X Relative rotation around X axis. +\end_layout + +\begin_layout Description +[ry]:X Relative rotation around Y axis. +\end_layout + +\begin_layout Description +[rz]:X Relateive rotation around Z axis. +\end_layout + \end_deeper \begin_layout Standard \shape italic -Note that data sets are not created until base pairing is determined. +Note that base pair data sets are not created until base pairing is determined. \end_layout \begin_layout Standard @@ -30066,6 +30103,16 @@ run writedata NApucker.dat NA[pucker] \end_layout +\begin_layout Standard +Note that while the underlying procedure is geared towards calculating parameter +s for base pairs, the code can be made to calculate parameters between consecuti +ve bases in single strands by specifying +\series bold +sscalc +\series default +. +\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 901cdfcf2a..cfdd2015fb 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -41,7 +41,7 @@ Action_NAstruct::Action_NAstruct() : {} void Action_NAstruct::Help() const { - mprintf("\t[] [resrange ] [naout ] [sscalc]\n" + mprintf("\t[] [resrange ] [sscalc] [naout ]\n" "\t[noheader] [resmap :{A,C,G,T,U} ...] [calcnohb]\n" "\t[noframespaces] [baseref ] ...\n" "\t[hbcut ] [origincut ] [altona | cremer]\n" From 0cba11eef05630dcb6d4ee8f61abb923205e1556 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 11:32:27 -0500 Subject: [PATCH 10/13] Revision bump for sscalc keyword --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index fcf9b83978..1fe194b5c4 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 "V5.0.4" +#define CPPTRAJ_INTERNAL_VERSION "V5.0.5" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif From 52eb74ed791610e879c9003c82caea63dd6c0ac7 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 13:22:48 -0500 Subject: [PATCH 11/13] Make map insertion a bit more efficient in nastruct. Make SendFrame const --- src/Action_NAstruct.cpp | 26 +++++++++++++------------- src/Box.cpp | 2 +- src/Box.h | 2 +- src/Frame.cpp | 2 +- src/Frame.h | 2 +- src/Matrix_3x3.cpp | 2 +- src/Matrix_3x3.h | 2 +- src/Parallel.cpp | 2 +- src/Parallel.h | 2 +- 9 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/Action_NAstruct.cpp b/src/Action_NAstruct.cpp index cfdd2015fb..dadb3edd02 100644 --- a/src/Action_NAstruct.cpp +++ b/src/Action_NAstruct.cpp @@ -408,9 +408,9 @@ Action_NAstruct::HbondType Action_NAstruct::ID_HBtype(NA_Base const& base1, int return (ATpair(base1, b1, base2, b2)); else if ( base1.Type() == NA_Base::THY && base2.Type() == NA_Base::ADE ) return (ATpair(base2, b2, base1, b1)); - else if ( base1.Type() == NA_Base::ADE && base2.Type() == NA_Base::URA ) // FIXME: OK for A-U? + else if ( base1.Type() == NA_Base::ADE && base2.Type() == NA_Base::URA ) // A-U has same WC pattern as A-T return (ATpair(base1, b1, base2, b2)); - else if ( base1.Type() == NA_Base::URA && base2.Type() == NA_Base::ADE ) // FIXME: OK for A-U? + else if ( base1.Type() == NA_Base::URA && base2.Type() == NA_Base::ADE ) // A-U has same WC pattern as A-T return (ATpair(base2, b2, base1, b1)); return OTHER; } @@ -517,7 +517,7 @@ int Action_NAstruct::DetermineBasePairing() { # ifdef NASTRUCTDEBUG mprintf("\n=================== Setup Base Pairing ===================\n"); # endif - // FIXME does data set name gen belong in Init()? + // Loop over all possible pairs of bases for (Barray::const_iterator base1 = Bases_.begin(); base1 != Bases_.end(); ++base1) { for (Barray::const_iterator base2 = base1 + 1; base2 != Bases_.end(); ++base2) @@ -1164,8 +1164,8 @@ int Action_NAstruct::DetermineStepParameters(int frameNum) { // pair steps are indexed by base pair indices. Rpair steppair(BP1.bpidx_, BP2.bpidx_); // Base pair step. Try to find existing base pair step. - StepMap::iterator entry = Steps_.find( steppair ); - if (entry == Steps_.end()) { + StepMap::iterator entry = Steps_.lower_bound( steppair ); + if (entry == Steps_.end() || entry->first != steppair) { // New base pair step StepType BS; MetaData md = NewStepType(BS, BP1.base1idx_, BP1.base2idx_, @@ -1211,7 +1211,7 @@ int Action_NAstruct::DetermineStepParameters(int frameNum) { // BS.minGroove_->legend(), BS.P_p1_+1, BS.P_p2_+1, BS.p_m1_+1, BS.p_m2_+1); } } - entry = Steps_.insert( entry, std::pair(steppair, BS) ); // FIXME does entry make more efficient? + entry = Steps_.insert( entry, std::pair(steppair, BS) ); # ifdef NASTRUCTDEBUG mprintf(" New base pair step: %s\n", md.Legend().c_str()); # endif @@ -1493,8 +1493,8 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { // Get/add data set for strandpair Rpair strandpair(b1idx, b2idx); - Smap::iterator entry = StrandPairs_.find( strandpair ); - if (entry == StrandPairs_.end()) { + Smap::iterator entry = StrandPairs_.lower_bound( strandpair ); + if (entry == StrandPairs_.end() || entry->first != strandpair) { // New strand pair Stype SP; // MetaData md = NewStepType(BS, BP1.base1idx_, BP1.base2idx_, @@ -1518,7 +1518,7 @@ Action::RetType Action_NAstruct::Setup(ActionSetup& setup) { //SP.base1idx_ = b1idx_; //SP.base2idx_ = b2idx_; - entry = StrandPairs_.insert( entry, std::pair(strandpair, SP) ); // FIXME does entry make more efficient? + entry = StrandPairs_.insert( entry, std::pair(strandpair, SP) ); } } // END loop over bases in strand } // END loop over strands @@ -1556,7 +1556,7 @@ Action::RetType Action_NAstruct::DoAction(int frameNum, ActionFrame& frm) { int err = 0; if (trajComm_.Master()) { // TODO MasterBcast? for (int rank = 1; rank < trajComm_.Size(); rank++) - frm.ModifyFrm().SendFrame( rank, trajComm_); // FIXME make SendFrame const + frm.Frm().SendFrame( rank, trajComm_); err = SetupBaseAxes(frm.Frm()); if (err==0) err = DetermineBasePairing(); } else { @@ -1710,8 +1710,8 @@ int Action_NAstruct::SyncAction() { int ii = 0; for (int in = 0; in != nsteps_on_rank[rank]; in++, ii += iSize) { Rpair steppair( iArray[ii], iArray[ii+1] ); - StepMap::iterator entry = Steps_.find( steppair ); - if (entry == Steps_.end()) { + StepMap::iterator entry = Steps_.lower_bound( steppair ); + if (entry == Steps_.end() || entry->first != steppair) { // New base pair step //mprintf("NEW BASE PAIR STEP: %i %i\n", iArray[0], iArray[1]); StepType BS; @@ -1729,7 +1729,7 @@ int Action_NAstruct::SyncAction() { BS.minGroove_ = (DataSet_1D*)masterDSL_->AddSet(DataSet::FLOAT, md); } } - entry = Steps_.insert( entry, std::pair(steppair, BS) ); // FIXME does entry make more efficient? + entry = Steps_.insert( entry, std::pair(steppair, BS) ); } //else mprintf("EXISTING BASE PAIR STEP: %i %i\n", entry->first.first, entry->first.second); // Synchronize all step data sets from rank. diff --git a/src/Box.cpp b/src/Box.cpp index c378106592..242acf342f 100644 --- a/src/Box.cpp +++ b/src/Box.cpp @@ -78,7 +78,7 @@ int Box::SyncBox(Parallel::Comm const& commIn) { } /** Send box info to recvrank. */ -int Box::SendBox(int recvrank, Parallel::Comm const& commIn) { +int Box::SendBox(int recvrank, Parallel::Comm const& commIn) const { commIn.Send( box_, 6, MPI_DOUBLE, recvrank, 1801 ); unitCell_.SendMatrix(recvrank, commIn); fracCell_.SendMatrix(recvrank, commIn); diff --git a/src/Box.h b/src/Box.h index 051bd26380..3edcf96083 100644 --- a/src/Box.h +++ b/src/Box.h @@ -25,7 +25,7 @@ class Box { void swap(Box&); # ifdef MPI int SyncBox(Parallel::Comm const&); - int SendBox(int, Parallel::Comm const&); + int SendBox(int, Parallel::Comm const&) const; int RecvBox(int, Parallel::Comm const&); # endif /// Remove all box information diff --git a/src/Frame.cpp b/src/Frame.cpp index 40fa5adab6..32dbd0aa50 100644 --- a/src/Frame.cpp +++ b/src/Frame.cpp @@ -1435,7 +1435,7 @@ void Frame::SetOrthoBoundingBox(std::vector const& Radii, double offset) #ifdef MPI // TODO: Change Frame class so everything can be sent in one MPI call. /** Send contents of this Frame to recvrank. */ -int Frame::SendFrame(int recvrank, Parallel::Comm const& commIn) { +int Frame::SendFrame(int recvrank, Parallel::Comm const& commIn) const { //rprintf("SENDING TO %i\n", recvrank); // DEBUG commIn.Send( X_, ncoord_, MPI_DOUBLE, recvrank, 1212 ); if (V_ != 0) diff --git a/src/Frame.h b/src/Frame.h index 610ab36618..8a388482a0 100644 --- a/src/Frame.h +++ b/src/Frame.h @@ -270,7 +270,7 @@ class Frame { void SetOrthoBoundingBox(std::vector const& Radii, double); # ifdef MPI // ----- Parallel Routines ------------------- - int SendFrame(int, Parallel::Comm const&); + int SendFrame(int, Parallel::Comm const&) const; int RecvFrame(int, Parallel::Comm const&); int SumToMaster(Parallel::Comm const&); # endif diff --git a/src/Matrix_3x3.cpp b/src/Matrix_3x3.cpp index 3960939479..f65cb4e4b4 100644 --- a/src/Matrix_3x3.cpp +++ b/src/Matrix_3x3.cpp @@ -537,7 +537,7 @@ void Matrix_3x3::SyncMatrix(Parallel::Comm const& commIn) { commIn.MasterBcast( M_, 9, MPI_DOUBLE ); } -int Matrix_3x3::SendMatrix(int recvrank, Parallel::Comm const& commIn) { +int Matrix_3x3::SendMatrix(int recvrank, Parallel::Comm const& commIn) const { commIn.Send( M_, 9, MPI_DOUBLE, recvrank, 1900 ); return 0; } diff --git a/src/Matrix_3x3.h b/src/Matrix_3x3.h index 689b392e0b..2302a8b5d7 100644 --- a/src/Matrix_3x3.h +++ b/src/Matrix_3x3.h @@ -97,7 +97,7 @@ class Matrix_3x3 { double* Dptr() { return M_; } # ifdef MPI void SyncMatrix(Parallel::Comm const&); - int SendMatrix(int, Parallel::Comm const&); + int SendMatrix(int, Parallel::Comm const&) const; int RecvMatrix(int, Parallel::Comm const&); # endif private: diff --git a/src/Parallel.cpp b/src/Parallel.cpp index d5e77af9e0..52f3bf63df 100644 --- a/src/Parallel.cpp +++ b/src/Parallel.cpp @@ -439,7 +439,7 @@ const } /** Send data to specified rank. */ -int Parallel::Comm::Send(void* sendbuffer, int sendcount, MPI_Datatype sendtype, int dest, int tag) +int Parallel::Comm::Send(const void* sendbuffer, int sendcount, MPI_Datatype sendtype, int dest, int tag) const { int err = MPI_Send( sendbuffer, sendcount, sendtype, dest, tag, comm_ ); diff --git a/src/Parallel.h b/src/Parallel.h index 532315c0f9..38ecd579a4 100644 --- a/src/Parallel.h +++ b/src/Parallel.h @@ -142,7 +142,7 @@ class Parallel::Comm { /// SendBuffer, Count, DataType, RecvBuffer int AllGather(void*, int, MPI_Datatype, void*) const; /// Buffer, Count, DataType, Destination Rank, Tag - int Send(void*, int, MPI_Datatype, int, int) const; + int Send(const void*, int, MPI_Datatype, int, int) const; /// Buffer, Count, DataType, Source Rank, Tag int Recv(void*, int, MPI_Datatype, int, int) const; /// Buffer, Count, DataType From 25454bb3de45d541c4655bad4f6c49b013574447 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 13:29:32 -0500 Subject: [PATCH 12/13] SendFrame is now const --- src/Energy_Sander.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Energy_Sander.cpp b/src/Energy_Sander.cpp index 7e135d61af..8c06e2a8a5 100644 --- a/src/Energy_Sander.cpp +++ b/src/Energy_Sander.cpp @@ -216,7 +216,7 @@ int Energy_Sander::Initialize(Topology const& topIn, Frame& fIn, Parallel::Comm err = CommonInit( topIn, fIn ); // Master sends reference frame. for (int rank = 1; rank < commIn.Size(); rank++) - fIn.SendFrame( rank, commIn ); // FIXME make SendFrame const + fIn.SendFrame( rank, commIn ); } else { // Check if master was able to write temporary top file. if (commIn.CheckError( err )) return 1; From 2e0ecb490d3e135e2b56afbfdfd9305a6fdb4267 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 11 Feb 2021 13:40:18 -0500 Subject: [PATCH 13/13] Resolve some FIXME comments --- src/Action_MakeStructure.h | 2 +- src/Action_OrderParameter.h | 1 - src/Analysis_Spline.cpp | 9 ++++----- src/Analysis_Spline.h | 1 - src/AtomMap.h | 2 +- src/ClusterNode.h | 1 - src/DataIO_RemLog.h | 2 +- src/DataSet.h | 2 +- src/DataSetList.h | 2 +- src/DataSet_2D.h | 2 +- src/DataSet_3D.h | 2 +- src/DataSet_Coords_TRJ.h | 2 +- src/DataSet_GridDbl.h | 2 +- src/DataSet_GridFlt.h | 2 +- src/DataSet_MatrixDbl.h | 2 +- src/DataSet_MatrixFlt.h | 2 +- src/Ewald_Regular.h | 2 +- src/Frame.h | 2 +- src/MetaData.h | 2 +- src/OutputTrajCommon.h | 4 ++-- src/Parallel.h | 2 +- src/ParameterTypes.h | 2 +- src/TrajoutList.h | 2 +- src/Trajout_Single.h | 4 ++-- 24 files changed, 26 insertions(+), 30 deletions(-) diff --git a/src/Action_MakeStructure.h b/src/Action_MakeStructure.h index 36f1891a44..c54000b22c 100644 --- a/src/Action_MakeStructure.h +++ b/src/Action_MakeStructure.h @@ -20,7 +20,7 @@ class Action_MakeStructure : public Action { SS_TYPE(double ph,double ps,double ph2,double ps2,int t,std::string const& n) : phi(ph), psi(ps), phi2(ph2), psi2(ps2), isTurn(t), type_arg(n) {} bool empty() { return (isTurn == -1); } - double phi, psi, phi2, psi2; // Angle(s) FIXME: Should this be an array? + double phi, psi, phi2, psi2; // Angle(s) TODO: Should this be an array? int isTurn; // 0=phi/psi, 1=Turn (phi/psi/phi2/psi2), 2=Single Dihedral(phi) std::string type_arg; }; diff --git a/src/Action_OrderParameter.h b/src/Action_OrderParameter.h index 1a87845d6d..d1a4c5adf2 100644 --- a/src/Action_OrderParameter.h +++ b/src/Action_OrderParameter.h @@ -41,7 +41,6 @@ class Action_OrderParameter : public Action { AtomMask tailend_mask_; AtomMask unsat_mask_; - // FIXME: // std::vector > masks_; std::vector masks_; std::vector > dbonds_; diff --git a/src/Analysis_Spline.cpp b/src/Analysis_Spline.cpp index 53807d9fe4..4f5f081bb9 100644 --- a/src/Analysis_Spline.cpp +++ b/src/Analysis_Spline.cpp @@ -3,7 +3,6 @@ #include "Constants.h" // SMALL Analysis_Spline::Analysis_Spline() : - outfile_(0), meshsize_(0), meshmin_(0.0), meshmax_(0.0), @@ -21,7 +20,7 @@ void Analysis_Spline::Help() const { Analysis::RetType Analysis_Spline::Setup(ArgList& analyzeArgs, AnalysisSetup& setup, int debugIn) { std::string setname = analyzeArgs.GetStringKey("name"); - outfile_ = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs); + DataFile* outfile = setup.DFL().AddDataFile(analyzeArgs.GetStringKey("out"), analyzeArgs); meshsize_ = analyzeArgs.getKeyInt("meshsize", 0); meshfactor_ = -1.0; if (meshsize_ < 3) { @@ -66,7 +65,7 @@ Analysis::RetType Analysis_Spline::Setup(ArgList& analyzeArgs, AnalysisSetup& se ds->SetLegend( "Spline(" + (*dsIn)->Meta().Legend() + ")" ); // TODO: Set individually based on input_dsets_ ds->SetDim(Dimension::X, Xdim); - if (outfile_ != 0) outfile_->AddDataSet( ds ); + if (outfile != 0) outfile->AddDataSet( ds ); output_dsets_.push_back( (DataSet_Mesh*)ds ); } @@ -83,10 +82,10 @@ Analysis::RetType Analysis_Spline::Setup(ArgList& analyzeArgs, AnalysisSetup& se mprintf(" Mesh max= %f\n", meshmax_); else mprintf(" Mesh max will be input set max.\n"); - if (outfile_ != 0) { + if (outfile != 0) { if (!setname.empty()) mprintf("\tOutput set name: %s\n", setname.c_str()); - mprintf("\tOutfile name: %s\n", outfile_->DataFilename().base()); + mprintf("\tOutfile name: %s\n", outfile->DataFilename().base()); } //for (Array1D::const_iterator set = input_dsets_.begin(); set != input_dsets_.end(); ++set) // mprintf("\t%s\n", (*set)->legend()); diff --git a/src/Analysis_Spline.h b/src/Analysis_Spline.h index 4902fa5c35..f5dc2e117a 100644 --- a/src/Analysis_Spline.h +++ b/src/Analysis_Spline.h @@ -12,7 +12,6 @@ class Analysis_Spline : public Analysis { Analysis::RetType Setup(ArgList&, AnalysisSetup&, int); Analysis::RetType Analyze(); private: - DataFile* outfile_; // FIXME: May not need to be class var Array1D input_dsets_; std::vector output_dsets_; int meshsize_; ///< Mesh size will be DataSet.Size * meshfactor diff --git a/src/AtomMap.h b/src/AtomMap.h index b3013e9fac..4d61aee317 100644 --- a/src/AtomMap.h +++ b/src/AtomMap.h @@ -8,7 +8,7 @@ class AtomMap { AtomMap(); MapAtom& operator[](int); - const MapAtom& operator[](int i) const { return mapatoms_[i]; } // FIXME: bounds + const MapAtom& operator[](int i) const { return mapatoms_[i]; } /// \return the number of atoms in the AtomMap. int Natom() const { return (int)mapatoms_.size(); } /// Set the debug level of the AtomMap. diff --git a/src/ClusterNode.h b/src/ClusterNode.h index c86a6c6993..313b054a77 100644 --- a/src/ClusterNode.h +++ b/src/ClusterNode.h @@ -24,7 +24,6 @@ class ClusterNode { void CalcEccentricity(DataSet_Cmatrix const&); /// Calculate centroid of members of this cluster. void CalculateCentroid(ClusterDist* Cdist) { - // FIXME: Could potentially get rid of this branch. if (centroid_ == 0) centroid_ = Cdist->NewCentroid( frameList_ ); else diff --git a/src/DataIO_RemLog.h b/src/DataIO_RemLog.h index 5df9ba3837..c6e3989676 100644 --- a/src/DataIO_RemLog.h +++ b/src/DataIO_RemLog.h @@ -20,7 +20,7 @@ class DataIO_RemLog : public DataIO { enum LogType { UNKNOWN = 0, TREMD, HREMD, MREMD, RXSGLD, PHREMD }; static const char* LogDescription[]; typedef std::vector Sarray; // TODO FileName array? - typedef std::map TmapType; // FIXME: Use ReplicaMap + typedef std::map TmapType; // TODO: Use ReplicaMap typedef DataSet_RemLog::IdxArray IdxArray; /// Read remlog header diff --git a/src/DataSet.h b/src/DataSet.h index 76a8ec7c0f..27b5aa57e2 100644 --- a/src/DataSet.h +++ b/src/DataSet.h @@ -49,7 +49,7 @@ class DataSet { virtual SizeArray DimSizes() const { return SizeArray(); } /// Print DataSet information //TODO return string instead? virtual void Info() const = 0; - /// Write data to file given start indices. FIXME Buffer? Should this function take number of elements as well? + /// Write data to file given start indices. virtual void WriteBuffer(CpptrajFile&, SizeArray const&) const = 0; /// \return value of coordinate for specified dimension d and position p. /** NOTE: It is assumed this can ALWAYS be represented as double precision. */ diff --git a/src/DataSetList.h b/src/DataSetList.h index c04d096194..d3756d05f9 100644 --- a/src/DataSetList.h +++ b/src/DataSetList.h @@ -30,7 +30,7 @@ class DataSetList { DataSetList& operator+=(DataSetList const&); /// \return DataSet at didx. - DataSet* operator[](int didx) const { return DataList_[didx]; } // FIXME: No bounds check + DataSet* operator[](int didx) const { return DataList_[didx]; } /// DataSetList default iterator typedef DataListType::const_iterator const_iterator; /// Iterator to beginning of dataset list diff --git a/src/DataSet_2D.h b/src/DataSet_2D.h index 5bd3be4cc5..5224c40d40 100644 --- a/src/DataSet_2D.h +++ b/src/DataSet_2D.h @@ -12,7 +12,7 @@ class DataSet_2D : public DataSet { DataSet(tIn, MATRIX_2D, fIn, 2) {} // TODO enable Append? int Append(DataSet*) { return 1; } - int Allocate(SizeArray const& s) { return Allocate2D(s[0], s[1]); } // FIXME bounds check + int Allocate(SizeArray const& s) { return Allocate2D(s[0], s[1]); } // TODO: 1 allocate using MatrixKind? /// Set up matrix for given # columns and rows. // TODO replace with Allocate virtual int Allocate2D(size_t, size_t) = 0; diff --git a/src/DataSet_3D.h b/src/DataSet_3D.h index 7967f48772..e280f9a89f 100644 --- a/src/DataSet_3D.h +++ b/src/DataSet_3D.h @@ -4,7 +4,7 @@ #include "GridBin.h" class Box; /// Interface for 3D DataSets. -// FIXME: Use DataSet Dims? +// TODO: Use DataSet Dims? class DataSet_3D : public DataSet { public: DataSet_3D() : gridBin_(0) {} diff --git a/src/DataSet_Coords_TRJ.h b/src/DataSet_Coords_TRJ.h index 47464ff786..d2b3d01da2 100644 --- a/src/DataSet_Coords_TRJ.h +++ b/src/DataSet_Coords_TRJ.h @@ -22,7 +22,7 @@ class DataSet_Coords_TRJ : public DataSet_Coords { void Info() const; void Add( size_t, const void* ) { return; } int Allocate(SizeArray const&) { return 0; } - size_t MemUsageInBytes() const { return IDX_.DataSize() + readFrame_.DataSize(); } // FIXME + size_t MemUsageInBytes() const { return IDX_.DataSize() + readFrame_.DataSize(); } // ----- DataSet_Coords functions ------------ /// DISABLED: Add a frame. void AddFrame(Frame const& fIn) { } diff --git a/src/DataSet_GridDbl.h b/src/DataSet_GridDbl.h index 4ebe9614d3..1025535e53 100644 --- a/src/DataSet_GridDbl.h +++ b/src/DataSet_GridDbl.h @@ -26,7 +26,7 @@ class DataSet_GridDbl : public DataSet_3D { // ----- DataSet functions ------------------- size_t Size() const { return grid_.size(); } # ifdef MPI - // FIXME: Currently just sums up. Should this be a separate Sync function? + // TODO: Currently just sums up. Should this be a separate Sync function? int Sync(size_t, std::vector const&, Parallel::Comm const&); # endif void Info() const { return; } diff --git a/src/DataSet_GridFlt.h b/src/DataSet_GridFlt.h index cb7169b364..12740dabae 100644 --- a/src/DataSet_GridFlt.h +++ b/src/DataSet_GridFlt.h @@ -26,7 +26,7 @@ class DataSet_GridFlt : public DataSet_3D { // ----- DataSet functions ------------------- size_t Size() const { return grid_.size(); } # ifdef MPI - // FIXME: Currently just sums up. Should this be a separate Sync function? + // TODO: Currently just sums up. Should this be a separate Sync function? int Sync(size_t, std::vector const&, Parallel::Comm const&); # endif void Info() const { return; } diff --git a/src/DataSet_MatrixDbl.h b/src/DataSet_MatrixDbl.h index e2166f87a0..c135c78bf1 100644 --- a/src/DataSet_MatrixDbl.h +++ b/src/DataSet_MatrixDbl.h @@ -14,7 +14,7 @@ class DataSet_MatrixDbl : public DataSet_2D { // ----- DataSet functions ------------------- size_t Size() const { return mat_.size(); } # ifdef MPI - // FIXME: Currently just sums up. Should this be a separate Sync function? + // TODO: Currently just sums up. Should this be a separate Sync function? int Sync(size_t, std::vector const&, Parallel::Comm const&); # endif void Info() const { return; } diff --git a/src/DataSet_MatrixFlt.h b/src/DataSet_MatrixFlt.h index f0c71ad752..d60c80ad4f 100644 --- a/src/DataSet_MatrixFlt.h +++ b/src/DataSet_MatrixFlt.h @@ -11,7 +11,7 @@ class DataSet_MatrixFlt : public DataSet_2D { // ----- DataSet functions ------------------- size_t Size() const { return mat_.size(); } # ifdef MPI - // FIXME: Currently just sums up. Should this be a separate Sync function? + // TODO: Currently just sums up. Should this be a separate Sync function? int Sync(size_t, std::vector const&, Parallel::Comm const&); int SendSet(int, Parallel::Comm const&); int RecvSet(int, Parallel::Comm const&); diff --git a/src/Ewald_Regular.h b/src/Ewald_Regular.h index f8bd958e44..0957f6a83b 100644 --- a/src/Ewald_Regular.h +++ b/src/Ewald_Regular.h @@ -39,7 +39,7 @@ class Ewald_Regular : public Ewald { Iarray mlim2_; ///< Hold m2 reciprocal indices int multCut_; ///< Hold index after which multiplier should be 2.0. # endif - double maxexp_; ///< Determines how far out recip vectors go? FIXME check! + double maxexp_; ///< Determines how far out recip vectors go? TODO check! double rsumTol_; ///< Reciprocal space sum tolerance. int mlimit_[3]; ///< Number of units in each direction to calc recip. sum. / nfft int maxmlim_; ///< The max of the three mlimit_ values. / pme spline order diff --git a/src/Frame.h b/src/Frame.h index 8a388482a0..a6ed3e06e2 100644 --- a/src/Frame.h +++ b/src/Frame.h @@ -286,7 +286,7 @@ class Frame { double T_; ///< Temperature double pH_; ///< pH double redox_; ///< RedOx potential - double time_; ///< Time FIXME Should this be float? + double time_; ///< Time double* X_; ///< Coord array, X0 Y0 Z0 X1 Y1 Z1 ... double* V_; ///< Velocities (same arrangement as Coords). double* F_; ///< Frame (same arrangement as Coords). diff --git a/src/MetaData.h b/src/MetaData.h index 818e66fba5..fb56b41ac4 100644 --- a/src/MetaData.h +++ b/src/MetaData.h @@ -25,7 +25,7 @@ class MetaData { CHI2, CHI3, CHI4, CHI5, PUCKER, NOE, - DIST, COVAR, MWCOVAR, //FIXME: May need more descriptive names + DIST, COVAR, MWCOVAR, CORREL, DISTCOVAR, IDEA, IREDMAT, DIHCOVAR, IREDVEC, diff --git a/src/OutputTrajCommon.h b/src/OutputTrajCommon.h index 575b87e961..3ad6a4478d 100644 --- a/src/OutputTrajCommon.h +++ b/src/OutputTrajCommon.h @@ -44,8 +44,8 @@ class OutputTrajCommon { /// \return 'true' if output frame range has been set up inline bool HasRange() const { return hasRange_ || !frameCount_.DefaultSettings(); } private: - FileName trajName_; // FIXME: Save this here? - Topology* trajParm_;// FIXME: Save this here? + FileName trajName_; ///< Output file name + Topology* trajParm_; ///< Pointer to associated Topology CoordinateInfo cInfo_; ///< Coordinate info for trajectory int NframesToWrite_; ///< Expected number of frames to be written. // Track frame numbers diff --git a/src/Parallel.h b/src/Parallel.h index 38ecd579a4..6beff59e1b 100644 --- a/src/Parallel.h +++ b/src/Parallel.h @@ -8,7 +8,7 @@ # undef MPI # define CPPTRAJ_MPI # include -# include // off_t FIXME necessary? +# include // off_t # ifdef PARALLEL_DEBUG_VERBOSE # include // for FILE # endif diff --git a/src/ParameterTypes.h b/src/ParameterTypes.h index 74cefe2f58..42cfb325b3 100644 --- a/src/ParameterTypes.h +++ b/src/ParameterTypes.h @@ -399,7 +399,7 @@ class NonbondParmType { void SetNHBterms(int n) { hbarray_.assign( n, HB_ParmType() ); } /// Set specified HB term HB_ParmType& SetHB(int i) { return hbarray_[i]; } - /// Set specified nbindex location to given value. FIXME no bounds check + /// Set specified nbindex location to given value. void SetNbIdx(int idx, int nbidx) { nbindex_[idx] = nbidx; } /// Add given LJ term to nonbond array and update nonbond index array. /** Certain routines in sander (like the 1-4 calcs) do NOT use the diff --git a/src/TrajoutList.h b/src/TrajoutList.h index 2dba9a6fad..292e3ba4d2 100644 --- a/src/TrajoutList.h +++ b/src/TrajoutList.h @@ -19,7 +19,7 @@ class TrajoutList { ~TrajoutList(); void SetDebug(int); void Clear(); - /// FIXME legacy function to maintain pytraj compatibility + /// NOTE: legacy function to maintain pytraj compatibility int AddTrajout(std::string const&, ArgList const&, Topology*); /// Add output trajectory to the list and associate with given topology. int AddTrajout(std::string const&, ArgList const&, DataSetList const&, Topology*); diff --git a/src/Trajout_Single.h b/src/Trajout_Single.h index 2e245774e0..e1988f2d76 100644 --- a/src/Trajout_Single.h +++ b/src/Trajout_Single.h @@ -12,10 +12,10 @@ class DataSetList; * however in case we do want to inherit in the future in a manner * similar to Trajin. */ -// FIXME: Should open state be tracked and error given if setup called +// TODO: Should open state be tracked and error given if setup called // twice for already open traj? Should append be default if // traj is opened twice? -// FIXME: Should take ArgList&, not const, so we can determine if there are unknown args later. +// TODO: Should take ArgList&, not const, so we can determine if there are unknown args later. class Trajout_Single { public: Trajout_Single() : trajio_(0), debug_(0) {}