From a5a8fc99e129168d685d0198b8f70338809ffb7a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Feb 2022 11:06:08 -0500 Subject: [PATCH 01/34] Start adding summary by parts --- src/Action_HydrogenBond.cpp | 16 ++++++++-------- src/Action_HydrogenBond.h | 35 +++++++++++++++++++++++++++-------- 2 files changed, 35 insertions(+), 16 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 2d8f57f0db..442a53760d 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -593,21 +593,21 @@ void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,i ds = (DataSet_integer*) masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"solventhb",hbidx)); if (UVseriesout_ != 0) UVseriesout_->AddDataSet( ds ); - ds->AddVal( fnum, 1 ); + //ds->AddVal( fnum, 1 ); } Hbond hb; if (udonor) { // Do not care about which solvent acceptor if (ds != 0) ds->SetLegend( CurrentParm_->TruncResAtomName(h_atom) + "-V" ); - hb = Hbond(dist,angle,ds,-1,h_atom,d_atom); + hb = Hbond(ds, -1, h_atom, d_atom, splitFrames_); } else { // Do not care about which solvent donor if (ds != 0) ds->SetLegend( CurrentParm_->TruncResAtomName(a_atom) + "-V" ); - hb = Hbond(dist,angle,ds,a_atom,-1,-1); + hb = Hbond(ds, a_atom, -1, -1, splitFrames_); } - UV_Map_.insert(it, std::pair(hbidx,hb)); + it = UV_Map_.insert(it, std::pair(hbidx,hb)); } else { // mprintf("DBG1: OLD hbond : %8i .. %8i - %8i\n", a_atom+1,h_atom+1,d_atom+1); - it->second.Update(dist,angle,fnum); } + it->second.Update(dist, angle, fnum, splitFrames_); } // Action_HydrogenBond::CalcSolvHbonds() @@ -677,13 +677,13 @@ void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, DataSet_integer* ds = 0; if (series_) { ds = UUset(a_atom, h_atom, d_atom); - ds->AddVal( fnum, 1 ); + //ds->AddVal( fnum, 1 ); } - UU_Map_.insert(it, std::pair(hbidx,Hbond(dist,angle,ds,a_atom,h_atom,d_atom))); + it = UU_Map_.insert(it, std::pair(hbidx,Hbond(ds, a_atom, h_atom, d_atom, splitFrames_))); } else { // mprintf("DBG1: OLD hbond : %8i .. %8i - %8i\n", a_atom+1,h_atom+1,d_atom+1); - it->second.Update(dist,angle,fnum); } + it->second.Update(dist, angle, fnum, splitFrames_); } // Action_HydrogenBond::CalcSiteHbonds() diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 5fc105764d..ebc3675c7e 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -4,6 +4,7 @@ #include "Action.h" #include "ImageOption.h" #include "DataSet_integer.h" +#include "OnlineVarT.h" #ifdef TIMER # include "Timer.h" #endif @@ -77,6 +78,7 @@ class Action_HydrogenBond : public Action { AtomMask SolventAcceptorMask_; AtomMask Mask_; ImageOption imageOpt_; ///< Used to determine if imaging should be performed + Iarray splitFrames_; ///< For calculating hydrogen bonds by parts # ifdef TIMER Timer t_action_; Timer t_uu_; @@ -147,14 +149,22 @@ class Action_HydrogenBond::Site { Iarray hlist_; ///< List of hydrogen indices int idx_; ///< Heavy atom index }; - +// ----------------------------------------------------------------------------- /// Track specific hydrogen bond. class Action_HydrogenBond::Hbond { public: Hbond() : dist_(0.0), angle_(0.0), data_(0), A_(-1), H_(-1), D_(-1), frames_(0) {} /// New hydrogen bond - Hbond(double d, double a, DataSet_integer* s, int ia, int ih, int id) : - dist_(d), angle_(a), data_(s), A_(ia), H_(ih), D_(id), frames_(1) {} + Hbond(DataSet_integer* s, int ia, int ih, int id, Iarray const& splits) : + dist_(0), angle_(0), data_(s), A_(ia), H_(ih), D_(id), frames_(0) + { + if (!splits.empty()) { + partsDist_.resize(splits.size()+1); + partsAng_.resize(splits.size()+1); + } + } +// Hbond(double d, double a, DataSet_integer* s, int ia, int ih, int id) : +// dist_(d), angle_(a), data_(s), A_(ia), H_(ih), D_(id), frames_(1) {} # ifdef _OPENMP /// Just record that hbond exists Hbond(double d, double a, int ia, int ih, int id) : @@ -189,11 +199,18 @@ class Action_HydrogenBond::Hbond { return (frames_ > rhs.frames_); } /// Update distance/angle/time series - void Update(double d, double a, int f) { - dist_ += d; - angle_ += a; + void Update(double distIn, double angIn, int fnum, Iarray const& splitFrames) { + dist_ += distIn; + angle_ += angIn; ++frames_; - if (data_ != 0) data_->AddVal(f, 1); + if (data_ != 0) data_->AddVal(fnum, 1); + if (!splitFrames.empty()) { + // Find the correct part NOTE assumes fnum never out of range + int part = 0; + while (fnum >= splitFrames[part]) part++; + partsDist_[part].accumulate( distIn ); + partsAng_[part].accumulate( angIn ); + } } void CalcAvg(); void FinishSeries(unsigned int); @@ -206,8 +223,10 @@ class Action_HydrogenBond::Hbond { int H_; ///< Hydrogen atom index int D_; ///< Donor atom index int frames_; ///< # frames this hydrogen bond has been present + std::vector> partsDist_; ///< Hold avg. distance, split by parts + std::vector> partsAng_; ///< Hold avg. angle, split by parts }; - +// ----------------------------------------------------------------------------- /// Track solvent bridge between 2 or more solute residues. class Action_HydrogenBond::Bridge { public: From 20b8aa676fa30b517cde5363c96f49aabba906ff Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Feb 2022 11:36:37 -0500 Subject: [PATCH 02/34] Add splitframe keyword and split analysis to avg. solute hbond printout --- src/Action_HydrogenBond.cpp | 34 +++++++++++++++++++++++++++++++++- src/Action_HydrogenBond.h | 12 ++++++++++-- 2 files changed, 43 insertions(+), 3 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 442a53760d..914e60c772 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -89,6 +89,25 @@ Action::RetType Action_HydrogenBond::Init(ArgList& actionArgs, ActionInit& init, acut_ = actionArgs.getKeyDouble("angle",135.0); noIntramol_ = actionArgs.hasKey("nointramol"); bridgeByAtom_ = actionArgs.hasKey("bridgebyatom"); + // Hbond split analysis + std::string splitarg = actionArgs.GetStringKey("splitframe"); + if (!splitarg.empty()) { + ArgList splits( splitarg, "," ); + if (splits.Nargs() < 1) { + mprinterr("Error: Invalid argument for 'splitframe': %s\n", splitarg.c_str()); + return Action::ERR; + } + int sf = splits.getNextInteger(-1); // User frame #s start at 1 + while (sf > 0) { + splitFrames_.push_back( sf - 1 ); + sf = splits.getNextInteger(-1); + } + if ((int)splitFrames_.size() < splits.Nargs()) { + mprinterr("Error: Invalid split frame arguments.\n"); + splits.CheckForMoreArgs(); + return Action::ERR; + } + } // Convert angle cutoff to radians acut_ *= Constants::DEGRAD; double dcut = actionArgs.getKeyDouble("dist",3.0); @@ -194,6 +213,12 @@ Action::RetType Action_HydrogenBond::Init(ArgList& actionArgs, ActionInit& init, mprintf("\tWriting # Hbond v time results to %s\n", DF->DataFilename().full()); if (avgout_ != 0) mprintf("\tWriting Hbond avgs to %s\n",avgout_->Filename().full()); + if (!splitFrames_.empty()) { + mprintf("\tWill split analysis at frames:"); + for (Iarray::const_iterator it = splitFrames_.begin(); it != splitFrames_.end(); ++it) + mprintf(" %i", *it + 1); + mprintf("\n"); + } if (calcSolvent_) { if (solvout_ != 0) mprintf("\tWriting solute-solvent hbond avgs to %s\n", solvout_->Filename().full()); @@ -1298,9 +1323,16 @@ void Action_HydrogenBond::Print() { Hname.append("_" + integerToString(hbond->H()+1)); Dname.append("_" + integerToString(hbond->D()+1)); } - avgout_->Printf("%-*s %*s %*s %8i %12.4f %12.4f %12.4f\n", + avgout_->Printf("%-*s %*s %*s %8i %12.4f %12.4f %12.4f", NUM, Aname.c_str(), NUM, Hname.c_str(), NUM, Dname.c_str(), hbond->Frames(), avg, hbond->Dist(), hbond->Angle()); + if (!splitFrames_.empty()) { + for (unsigned int idx = 0; idx != hbond->Nparts(); idx++) + avgout_->Printf(" %8i %12.4f %12.4f %12.4f", + hbond->PartFrames(idx), hbond->PartFrac(idx, Nframes_), + hbond->PartDist(idx), hbond->PartAngle(idx)); + } + avgout_->Printf("\n"); } } diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index ebc3675c7e..022192db06 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -5,6 +5,7 @@ #include "ImageOption.h" #include "DataSet_integer.h" #include "OnlineVarT.h" +//#inc lude "CpptrajStdio.h" // DEBUG #ifdef TIMER # include "Timer.h" #endif @@ -179,6 +180,12 @@ class Action_HydrogenBond::Hbond { int A() const { return A_; } int H() const { return H_; } int D() const { return D_; } + // Summary by parts + unsigned int Nparts() const { return partsDist_.size(); } + unsigned int PartFrames(unsigned int idx) const { return (unsigned int)partsDist_[idx].nData(); } + double PartFrac(unsigned int idx, unsigned int Nframes) const { return partsDist_[idx].nData() / (double)Nframes; } + double PartDist(unsigned int idx) const { return partsDist_[idx].mean(); } + double PartAngle(unsigned int idx) const { return partsAng_[idx].mean(); } # ifdef MPI DataSet_integer* Data() const { return data_; } /// New hydrogen bond with given # frames @@ -205,9 +212,10 @@ class Action_HydrogenBond::Hbond { ++frames_; if (data_ != 0) data_->AddVal(fnum, 1); if (!splitFrames.empty()) { + //mprintf("DEBUG: SPLIT: fnum= %i partsDistSize= %zu splitFramesSize= %zu\n", fnum, partsDist_.size(), splitFrames.size()); // Find the correct part NOTE assumes fnum never out of range - int part = 0; - while (fnum >= splitFrames[part]) part++; + unsigned int part = 0; + while (part < splitFrames.size() && fnum >= splitFrames[part]) part++; partsDist_[part].accumulate( distIn ); partsAng_[part].accumulate( angIn ); } From d13625b7559750e3c12d198d0478bb6315e7605c Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Feb 2022 11:41:47 -0500 Subject: [PATCH 03/34] Ensure angle is converted to degrees --- src/Action_HydrogenBond.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 914e60c772..69079e655f 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1330,7 +1330,7 @@ void Action_HydrogenBond::Print() { for (unsigned int idx = 0; idx != hbond->Nparts(); idx++) avgout_->Printf(" %8i %12.4f %12.4f %12.4f", hbond->PartFrames(idx), hbond->PartFrac(idx, Nframes_), - hbond->PartDist(idx), hbond->PartAngle(idx)); + hbond->PartDist(idx), hbond->PartAngle(idx)*Constants::RADDEG); } avgout_->Printf("\n"); } From 4304b022ab4ee23f7e5ee8fa9eca495f25894685 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Feb 2022 11:43:26 -0500 Subject: [PATCH 04/34] Add first test for hbond splitframe --- test/Test_Hbond_Split/RunTest.sh | 32 ++++++++++++++++++ test/Test_Hbond_Split/avghb.dat.save | 49 ++++++++++++++++++++++++++++ 2 files changed, 81 insertions(+) create mode 100755 test/Test_Hbond_Split/RunTest.sh create mode 100644 test/Test_Hbond_Split/avghb.dat.save diff --git a/test/Test_Hbond_Split/RunTest.sh b/test/Test_Hbond_Split/RunTest.sh new file mode 100755 index 0000000000..f0deed0a6e --- /dev/null +++ b/test/Test_Hbond_Split/RunTest.sh @@ -0,0 +1,32 @@ +#!/bin/bash + +. ../MasterTest.sh + +# Clean +CleanFiles hbond.in nhb.dat avghb.dat + +INPUT="-i hbond.in" + +# Solute-solute, split by parts +TestUUsplit() { + UNITNAME='Solute Hbond test' + CheckFor netcdf + if [ $? -eq 0 ] ; then + cat > hbond.in < Date: Tue, 15 Feb 2022 14:11:24 -0500 Subject: [PATCH 05/34] Make the split by parts use overall trajout num --- src/Action_HydrogenBond.cpp | 16 ++++++++-------- src/Action_HydrogenBond.h | 12 ++++++------ 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 69079e655f..14396e0dba 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -587,7 +587,7 @@ double Action_HydrogenBond::Angle(const double* XA, const double* XH, const doub } // Action_HydrogenBond::AddUV() -void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,int h_atom,int d_atom,bool udonor) +void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,int h_atom,int d_atom,bool udonor, int onum) { int hbidx, solventres, soluteres; // TODO: Option to use solvent mol num? @@ -632,7 +632,7 @@ void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,i } else { // mprintf("DBG1: OLD hbond : %8i .. %8i - %8i\n", a_atom+1,h_atom+1,d_atom+1); } - it->second.Update(dist, angle, fnum, splitFrames_); + it->second.Update(dist, angle, fnum, splitFrames_, onum); } // Action_HydrogenBond::CalcSolvHbonds() @@ -660,7 +660,7 @@ void Action_HydrogenBond::CalcSolvHbonds(int frameNum, double dist2, thread_HBs_[numHB].push_back( Hbond(sqrt(dist2), angle, a_atom, *h_atom, d_atom, (int)soluteDonor) ); # else ++numHB; - AddUV(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom, soluteDonor); + AddUV(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom, soluteDonor, frmIn.TrajoutNum()); # endif } } @@ -691,7 +691,7 @@ DataSet_integer* Action_HydrogenBond::UUset(int a_atom, int h_atom, int d_atom) } // Action_HydrogenBond::AddUU() -void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, int h_atom, int d_atom) +void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, int h_atom, int d_atom, int onum) { // Index UU hydrogen bonds by DonorH-Acceptor Hpair hbidx(h_atom, a_atom); @@ -708,7 +708,7 @@ void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, } else { // mprintf("DBG1: OLD hbond : %8i .. %8i - %8i\n", a_atom+1,h_atom+1,d_atom+1); } - it->second.Update(dist, angle, fnum, splitFrames_); + it->second.Update(dist, angle, fnum, splitFrames_, onum); } // Action_HydrogenBond::CalcSiteHbonds() @@ -730,7 +730,7 @@ void Action_HydrogenBond::CalcSiteHbonds(int frameNum, double dist2, thread_HBs_[numHB].push_back( Hbond(sqrt(dist2), angle, a_atom, *h_atom, d_atom) ); # else ++numHB; - AddUU(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom); + AddUU(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom, frmIn.TrajoutNum()); # endif } } @@ -844,7 +844,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { for (std::vector::iterator it = thread_HBs_.begin(); it != thread_HBs_.end(); ++it) { numHB += (int)it->size(); for (Harray::const_iterator hb = it->begin(); hb != it->end(); ++hb) - AddUU(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D()); + AddUU(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D(), frm.TrajoutNum()); it->clear(); } # endif @@ -911,7 +911,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { for (std::vector::iterator it = thread_HBs_.begin(); it != thread_HBs_.end(); ++it) { numHB += (int)it->size(); for (Harray::const_iterator hb = it->begin(); hb != it->end(); ++hb) - AddUV(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D(), (bool)hb->Frames()); + AddUV(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D(), (bool)hb->Frames(), frm.TrajoutNum()); it->clear(); } # endif diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 022192db06..d64e4b54ba 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -33,8 +33,8 @@ class Action_HydrogenBond : public Action { inline double Angle(const double*, const double*, const double*, Box const&) const; inline int UU_Set_Idx(int,int) const; inline DataSet_integer* UUset(int,int,int); - void AddUU(double,double,int,int,int,int); - void AddUV(double,double,int,int,int,int,bool); + void AddUU(double,double,int,int,int,int,int); + void AddUV(double,double,int,int,int,int,bool,int); void CalcSiteHbonds(int,double,Site const&,const double*,int,const double*, Frame const&, int&); void CalcSolvHbonds(int,double,Site const&,const double*,int,const double*, @@ -206,16 +206,16 @@ class Action_HydrogenBond::Hbond { return (frames_ > rhs.frames_); } /// Update distance/angle/time series - void Update(double distIn, double angIn, int fnum, Iarray const& splitFrames) { + void Update(double distIn, double angIn, int fnum, Iarray const& splitFrames, int onum) { dist_ += distIn; angle_ += angIn; ++frames_; if (data_ != 0) data_->AddVal(fnum, 1); if (!splitFrames.empty()) { - //mprintf("DEBUG: SPLIT: fnum= %i partsDistSize= %zu splitFramesSize= %zu\n", fnum, partsDist_.size(), splitFrames.size()); - // Find the correct part NOTE assumes fnum never out of range + //mprintf("DEBUG: SPLIT: onum= %i partsDistSize= %zu splitFramesSize= %zu\n", onum, partsDist_.size(), splitFrames.size()); + // Find the correct part NOTE assumes onum never out of range unsigned int part = 0; - while (part < splitFrames.size() && fnum >= splitFrames[part]) part++; + while (part < splitFrames.size() && onum >= splitFrames[part]) part++; partsDist_[part].accumulate( distIn ); partsAng_[part].accumulate( angIn ); } From d8064d177cd32c14b89971f077c235780462053d Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Feb 2022 14:53:02 -0500 Subject: [PATCH 06/34] Consolidate code for converting hbond to flat arrays for communication. Ensure part data is communicated. --- src/Action_HydrogenBond.cpp | 52 ++++++++++++++++++++++++++++--------- src/Action_HydrogenBond.h | 5 ++-- 2 files changed, 43 insertions(+), 14 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 14396e0dba..6ebe3a65c7 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1012,6 +1012,25 @@ std::vector Action_HydrogenBond::GetRankNhbonds( int num_hb, Parallel::Comm return nhb_on_rank; } +/** Add Hbond class to flat arrays. */ +void Action_HydrogenBond::HbondToArray(std::vector& Dvals, std::vector& Ivals, Hbond const& hb) +{ + Dvals.push_back( hb.Dist() ); + Dvals.push_back( hb.Angle() ); + for (unsigned int idx = 0; idx != hb.Nparts(); idx++) { + Dvals.push_back( hb.PartDist(idx).nData() ); + Dvals.push_back( hb.PartDist(idx).mean() ); + Dvals.push_back( hb.PartDist(idx).M2() ); + Dvals.push_back( hb.PartAngle(idx).nData() ); + Dvals.push_back( hb.PartAngle(idx).mean() ); + Dvals.push_back( hb.PartAngle(idx).M2() ); + } + Ivals.push_back( hb.A() ); + Ivals.push_back( hb.H() ); + Ivals.push_back( hb.D() ); + Ivals.push_back( hb.Frames() ); +} + /** PARALLEL NOTES: * The following tags are used for MPI send/receive: * 1300 : Array containing hbond double info on rank. @@ -1039,12 +1058,18 @@ int Action_HydrogenBond::SyncAction() { } // Need to send hbond data from all ranks to master. - std::vector Dvals; // Hold dist_ and angle_ for each hbond + std::vector Dvals; // Hold dist_ and angle_ for each hbond, as well as n_/mean_/M2_ for dist/angle each part std::vector Ivals; // Hold A_, H_, D_, and frames_ for each hbond + unsigned int dvalsPerHbond; + if (!splitFrames_.empty()) + dvalsPerHbond = 2 + ((splitFrames_.size()+1) * 6); + else + dvalsPerHbond = 2; // Need to know how many hbonds on each thread. std::vector nuu_on_rank = GetRankNhbonds( UU_Map_.size(), trajComm_ ); std::vector nuv_on_rank = GetRankNhbonds( UV_Map_.size(), trajComm_ ); if (trajComm_.Master()) { + // MASTER RANK for (int rank = 1; rank < trajComm_.Size(); rank++) { int n_total_on_rank = nuu_on_rank[rank] + nuv_on_rank[rank]; @@ -1052,8 +1077,8 @@ int Action_HydrogenBond::SyncAction() { { std::vector Svals; if (series_) Svals.reserve( n_total_on_rank ); - Dvals.resize( 2 * n_total_on_rank ); - Ivals.resize( 4 * n_total_on_rank ); + Dvals.resize( dvalsPerHbond * n_total_on_rank ); + Ivals.resize( 4 * n_total_on_rank ); trajComm_.Recv( &(Dvals[0]), Dvals.size(), MPI_DOUBLE, rank, 1300 ); trajComm_.Recv( &(Ivals[0]), Ivals.size(), MPI_INT, rank, 1301 ); // Loop over all received hydrogen bonds @@ -1062,7 +1087,7 @@ int Action_HydrogenBond::SyncAction() { const int* IV = &Ivals[0]; const double* DV = &Dvals[0]; // UU Hbonds - for (int in = 0; in != nuu_on_rank[rank]; in++, IV += 4, DV += 2) + for (int in = 0; in != nuu_on_rank[rank]; in++, IV += 4, DV += dvalsPerHbond) { Hpair hbidx(IV[1], IV[0]); UUmapType::iterator it = UU_Map_.lower_bound( hbidx ); @@ -1081,7 +1106,7 @@ int Action_HydrogenBond::SyncAction() { Svals.push_back( ds ); } // UV Hbonds - for (int in = 0; in != nuv_on_rank[rank]; in++, IV += 4, DV += 2) + for (int in = 0; in != nuv_on_rank[rank]; in++, IV += 4, DV += dvalsPerHbond) { int hbidx; if (IV[1] != -1) @@ -1135,27 +1160,30 @@ int Action_HydrogenBond::SyncAction() { hb->second.FinishSeries( Nframes_ ); } } else { + // NON-MASTER RANK if (!UU_Map_.empty()) { unsigned int ntotal = UU_Map_.size() + UV_Map_.size(); - Dvals.reserve( ntotal * 2 ); + Dvals.reserve( ntotal * dvalsPerHbond ); Ivals.reserve( ntotal * 4 ); // Store UU bonds in flat arrays. for (UUmapType::const_iterator it = UU_Map_.begin(); it != UU_Map_.end(); ++it) { - Dvals.push_back( it->second.Dist() ); + /*Dvals.push_back( it->second.Dist() ); Dvals.push_back( it->second.Angle() ); Ivals.push_back( it->second.A() ); Ivals.push_back( it->second.H() ); Ivals.push_back( it->second.D() ); - Ivals.push_back( it->second.Frames() ); + Ivals.push_back( it->second.Frames() );*/ + HbondToArray(Dvals, Ivals, it->second); } // Store UV bonds in flat arrays for (UVmapType::const_iterator it = UV_Map_.begin(); it != UV_Map_.end(); ++it) { - Dvals.push_back( it->second.Dist() ); + /*Dvals.push_back( it->second.Dist() ); Dvals.push_back( it->second.Angle() ); Ivals.push_back( it->second.A() ); Ivals.push_back( it->second.H() ); Ivals.push_back( it->second.D() ); - Ivals.push_back( it->second.Frames() ); + Ivals.push_back( it->second.Frames() );*/ + HbondToArray(Dvals, Ivals, it->second); } trajComm_.Send( &(Dvals[0]), Dvals.size(), MPI_DOUBLE, 0, 1300 ); trajComm_.Send( &(Ivals[0]), Ivals.size(), MPI_INT, 0, 1301 ); @@ -1176,7 +1204,7 @@ int Action_HydrogenBond::SyncAction() { } } } - } + } // END COMMUNICATING HBOND DATA TO MASTER if (calcSolvent_) { // Sync bridging data @@ -1330,7 +1358,7 @@ void Action_HydrogenBond::Print() { for (unsigned int idx = 0; idx != hbond->Nparts(); idx++) avgout_->Printf(" %8i %12.4f %12.4f %12.4f", hbond->PartFrames(idx), hbond->PartFrac(idx, Nframes_), - hbond->PartDist(idx), hbond->PartAngle(idx)*Constants::RADDEG); + hbond->PartDist(idx).mean(), hbond->PartAngle(idx).mean()*Constants::RADDEG); } avgout_->Printf("\n"); } diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index d64e4b54ba..0ea8236ecd 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -45,6 +45,7 @@ class Action_HydrogenBond : public Action { std::string MemoryUsage(size_t, size_t, size_t) const; # ifdef MPI static std::vector GetRankNhbonds(int,Parallel::Comm const&); + static void HbondToArray(std::vector&, std::vector&, Hbond const&); # endif typedef std::vector Sarray; @@ -184,8 +185,8 @@ class Action_HydrogenBond::Hbond { unsigned int Nparts() const { return partsDist_.size(); } unsigned int PartFrames(unsigned int idx) const { return (unsigned int)partsDist_[idx].nData(); } double PartFrac(unsigned int idx, unsigned int Nframes) const { return partsDist_[idx].nData() / (double)Nframes; } - double PartDist(unsigned int idx) const { return partsDist_[idx].mean(); } - double PartAngle(unsigned int idx) const { return partsAng_[idx].mean(); } + Stats const& PartDist(unsigned int idx) const { return partsDist_[idx]; } + Stats const& PartAngle(unsigned int idx) const { return partsAng_[idx]; } # ifdef MPI DataSet_integer* Data() const { return data_; } /// New hydrogen bond with given # frames From 6b4d42fd2aef837857f30729bcbfdf32f10352eb Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Feb 2022 18:54:23 -0500 Subject: [PATCH 07/34] Try to do combination for parts. Does not appear to work yet --- src/Action_HydrogenBond.cpp | 14 ++++++++++---- src/Action_HydrogenBond.h | 23 ++++++++++++++++++++++- 2 files changed, 32 insertions(+), 5 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 6ebe3a65c7..bfb1e640ae 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1061,10 +1061,14 @@ int Action_HydrogenBond::SyncAction() { std::vector Dvals; // Hold dist_ and angle_ for each hbond, as well as n_/mean_/M2_ for dist/angle each part std::vector Ivals; // Hold A_, H_, D_, and frames_ for each hbond unsigned int dvalsPerHbond; - if (!splitFrames_.empty()) - dvalsPerHbond = 2 + ((splitFrames_.size()+1) * 6); - else + unsigned int nParts; + if (!splitFrames_.empty()) { + nParts = splitFrames_.size() + 1; + dvalsPerHbond = 2 + (nParts * 6); + } else { + nParts = 0; dvalsPerHbond = 2; + } // Need to know how many hbonds on each thread. std::vector nuu_on_rank = GetRankNhbonds( UU_Map_.size(), trajComm_ ); std::vector nuv_on_rank = GetRankNhbonds( UV_Map_.size(), trajComm_ ); @@ -1097,10 +1101,12 @@ int Action_HydrogenBond::SyncAction() { // Hbond on rank that has not been found on master if (series_) ds = UUset(IV[0], IV[1], IV[2]); - UU_Map_.insert(it, std::pair(hbidx,Hbond(DV[0],DV[1],ds,IV[0],IV[1],IV[2],IV[3]))); + it = UU_Map_.insert(it, std::pair(hbidx,Hbond(DV[0],DV[1],ds,IV[0],IV[1],IV[2],IV[3]))); + it->second.SetupParts(nParts, DV+2); } else { // Hbond on rank and master. Update on master. it->second.Combine(DV[0], DV[1], IV[3]); + it->second.CombineParts(nParts, DV+2); ds = it->second.Data(); } Svals.push_back( ds ); diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 0ea8236ecd..d567fa763e 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -189,7 +189,7 @@ class Action_HydrogenBond::Hbond { Stats const& PartAngle(unsigned int idx) const { return partsAng_[idx]; } # ifdef MPI DataSet_integer* Data() const { return data_; } - /// New hydrogen bond with given # frames + /// CONSTRUCTOR - New hydrogen bond with given # frames Hbond(double d, double a, DataSet_integer* s, int ia, int ih, int id, int n) : dist_(d), angle_(a), data_(s), A_(ia), H_(ih), D_(id), frames_(n) {} /// Update distance/angle/number frames @@ -198,6 +198,27 @@ class Action_HydrogenBond::Hbond { angle_ += a; frames_ += n; } + /// Set up parts with given double array containing N, mean, and M2 for each part + void SetupParts(unsigned int nparts, const double* dvals) { + if (nparts == 0) return; + partsDist_.clear(); + partsAng_.clear(); + partsDist_.reserve( nparts ); + partsAng_.reserve( nparts ); + const double* dptr = dvals; + for (unsigned int idx = 0; idx != nparts; idx++, dptr += 6) { + partsDist_.push_back( Stats(dptr[0], dptr[1], dptr[2]) ); + partsAng_.push_back( Stats(dptr[3], dptr[4], dptr[5]) ); + } + } + /// Update parts with given double array containing N, mean, and M2 for each part + void CombineParts(unsigned int nparts, const double* dvals) { + const double* dptr = dvals; + for (unsigned int idx = 0; idx != nparts; idx++, dptr += 6) { + partsDist_[idx].Combine( Stats(dptr[0], dptr[1], dptr[2]) ); + partsAng_[idx].Combine( Stats(dptr[3], dptr[4], dptr[5]) ); + } + } # endif /// First sort by frames (descending), then distance (ascending). bool operator<(const Hbond& rhs) const { From 348fcc5b9a6291087576b1f0e752afc68e6b4346 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 08:44:14 -0500 Subject: [PATCH 08/34] When combining, handle case where one (or both) Stats have nothing in them. --- src/OnlineVarT.h | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/src/OnlineVarT.h b/src/OnlineVarT.h index fb43b6edba..183d961e18 100644 --- a/src/OnlineVarT.h +++ b/src/OnlineVarT.h @@ -16,7 +16,7 @@ template class Stats { public: /// Default CONSTRUCTOR - Stats() : n_(0.0), mean_(0.0), M2_(0.0) {} + Stats() : n_(0), mean_(0), M2_(0) {} # ifdef MPI /// CONSTRUCTOR taking N, mean, and M2 Stats(Float n, Float mean, Float m2) : n_(n), mean_(mean), M2_(m2) {} @@ -37,20 +37,26 @@ class Stats { Float mean() const { return mean_; }; /// \return Current variance Float variance() const { - if (n_ < 2) return 0.0; + if (n_ < 2) return 0; return M2_ / (n_ - 1.0); }; /// \return Current number of data points in mean/variance Float nData() const { return n_; }; /// Combine given Stats with this Stats void Combine(Stats const& rhs) { - Float nX = n_ + rhs.n_; - Float delta = rhs.mean_ - mean_; - // Combined mean - mean_ = mean_ + delta * (rhs.n_ / nX); - // Combined variance - M2_ = M2_ + rhs.M2_ + ((delta*delta) * ((n_*rhs.n_) / nX)); - n_ = nX; + if (((unsigned long)n_) == 0) { + n_ = rhs.n_; + mean_ = rhs.mean_; + M2_ = rhs.M2_; + } else if (((unsigned long)rhs.n_) != 0) { + Float nX = n_ + rhs.n_; + Float delta = rhs.mean_ - mean_; + // Combined mean + mean_ = mean_ + delta * (rhs.n_ / nX); + // Combined variance + M2_ = M2_ + rhs.M2_ + ((delta*delta) * ((n_*rhs.n_) / nX)); + n_ = nX; + } } private: Float n_; @@ -68,7 +74,7 @@ template class StatsMap { public: StatsMap() : - n_(0.0), min_(0), max_(0) + n_(0), min_(0), max_(0) {} # ifdef MPI /// CONSTRUCTOR - intended to create StatsMap from data transfered via MPI. @@ -130,7 +136,7 @@ class StatsMap { Value mean(Key i) { return mean_[i]; }; Value variance(Key i) { - if (n_ < 2) return 0.0; + if (n_ < 2) return 0; return M2_[i] / (n_ - 1.0); }; From 2a0ed8db8934f9c0bf473c9efd78771a8242492d Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 09:07:45 -0500 Subject: [PATCH 09/34] Enable split summary for solvent solute hbonds --- src/Action_HydrogenBond.cpp | 25 +++++++++++++++++-------- src/Action_HydrogenBond.h | 2 ++ 2 files changed, 19 insertions(+), 8 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index bfb1e640ae..86e7f6f76e 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1132,10 +1132,12 @@ int Action_HydrogenBond::SyncAction() { ds->SetLegend( CreateHBlegend(*CurrentParm_, IV[0], IV[1], IV[2]) ); if (UVseriesout_ != 0) UVseriesout_->AddDataSet( ds ); } - UV_Map_.insert(it, std::pair(hbidx,Hbond(DV[0],DV[1],ds,IV[0],IV[1],IV[2],IV[3]))); + it = UV_Map_.insert(it, std::pair(hbidx,Hbond(DV[0],DV[1],ds,IV[0],IV[1],IV[2],IV[3]))); + it->second.SetupParts(nParts, DV+2); } else { // Hbond on rank and master. Update on master. it->second.Combine(DV[0], DV[1], IV[3]); + it->second.CombineParts(nParts, DV+2); ds = it->second.Data(); } Svals.push_back( ds ); @@ -1301,6 +1303,14 @@ std::string Action_HydrogenBond::MemoryUsage(size_t n_uu_pairs, size_t n_uv_pair return ByteString( memTotal, BYTE_DECIMAL ); } +/** Print summary by parts for given hbond. */ +void Action_HydrogenBond::summary_Parts(CpptrajFile* avgout, Hbond const& hb) const { + for (unsigned int idx = 0; idx != hb.Nparts(); idx++) + avgout->Printf(" %8i %12.4f %12.4f %12.4f", + hb.PartFrames(idx), hb.PartFrac(idx, Nframes_), + hb.PartDist(idx).mean(), hb.PartAngle(idx).mean()*Constants::RADDEG); +} + // Action_HydrogenBond::Print() /** Print average occupancies over all frames for all detected Hbonds. */ void Action_HydrogenBond::Print() { @@ -1360,12 +1370,8 @@ void Action_HydrogenBond::Print() { avgout_->Printf("%-*s %*s %*s %8i %12.4f %12.4f %12.4f", NUM, Aname.c_str(), NUM, Hname.c_str(), NUM, Dname.c_str(), hbond->Frames(), avg, hbond->Dist(), hbond->Angle()); - if (!splitFrames_.empty()) { - for (unsigned int idx = 0; idx != hbond->Nparts(); idx++) - avgout_->Printf(" %8i %12.4f %12.4f %12.4f", - hbond->PartFrames(idx), hbond->PartFrac(idx, Nframes_), - hbond->PartDist(idx).mean(), hbond->PartAngle(idx).mean()*Constants::RADDEG); - } + if (!splitFrames_.empty()) + summary_Parts(avgout_, *hbond); avgout_->Printf("\n"); } } @@ -1406,9 +1412,12 @@ void Action_HydrogenBond::Print() { Hname.append("_" + integerToString(hbond->H()+1)); } } - solvout_->Printf("%-*s %*s %*s %8i %12.4f %12.4f %12.4f\n", + solvout_->Printf("%-*s %*s %*s %8i %12.4f %12.4f %12.4f", NUM, Aname.c_str(), NUM, Hname.c_str(), NUM, Dname.c_str(), hbond->Frames(), avg, hbond->Dist(), hbond->Angle()); + if (!splitFrames_.empty()) + summary_Parts(solvout_, *hbond); + solvout_->Printf("\n"); } HbondList.clear(); } diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index d567fa763e..f73bc76698 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -43,6 +43,8 @@ class Action_HydrogenBond : public Action { void UpdateSeries(); /// Determine memory usage from # hbonds and time series std::string MemoryUsage(size_t, size_t, size_t) const; + /// Print summary of parts for hbond + void summary_Parts(CpptrajFile*, Hbond const&) const; # ifdef MPI static std::vector GetRankNhbonds(int,Parallel::Comm const&); static void HbondToArray(std::vector&, std::vector&, Hbond const&); From 77e999164d55cf1cb207a8e5c60a5dcf7b321ac5 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 09:08:08 -0500 Subject: [PATCH 10/34] Add test for solute solvent hbonds with split --- test/Test_Hbond_Split/RunTest.sh | 18 ++++++++++- test/Test_Hbond_Split/solvavg.dat.save | 44 ++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 test/Test_Hbond_Split/solvavg.dat.save diff --git a/test/Test_Hbond_Split/RunTest.sh b/test/Test_Hbond_Split/RunTest.sh index f0deed0a6e..97d10dc985 100755 --- a/test/Test_Hbond_Split/RunTest.sh +++ b/test/Test_Hbond_Split/RunTest.sh @@ -9,7 +9,7 @@ INPUT="-i hbond.in" # Solute-solute, split by parts TestUUsplit() { - UNITNAME='Solute Hbond test' + UNITNAME='Solute Hbond test with split' CheckFor netcdf if [ $? -eq 0 ] ; then cat > hbond.in < hbond.in < Date: Wed, 16 Feb 2022 09:18:21 -0500 Subject: [PATCH 11/34] Write a proper header for split --- src/Action_HydrogenBond.cpp | 28 ++++++++++++++++++++++---- src/Action_HydrogenBond.h | 2 ++ test/Test_Hbond_Split/avghb.dat.save | 2 +- test/Test_Hbond_Split/solvavg.dat.save | 2 +- 4 files changed, 28 insertions(+), 6 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 86e7f6f76e..fc86b27db1 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1303,6 +1303,20 @@ std::string Action_HydrogenBond::MemoryUsage(size_t n_uu_pairs, size_t n_uv_pair return ByteString( memTotal, BYTE_DECIMAL ); } +/** Print header for summary by parts. */ +void Action_HydrogenBond::summary_Parts_header(CpptrajFile* avgout, unsigned int nParts) +{ + if (nParts < 1) return; + for (unsigned int idx = 0; idx != nParts; idx++) { + std::string spart(integerToString(idx+1)); + std::string frames( "Frames" + spart); + std::string frac( "Frac" + spart); + std::string avgdist("AvgDist" + spart); + std::string avgang( "AvgAng" + spart); + avgout->Printf(" %8s %12s %12s %12s", frames.c_str(), frac.c_str(), avgdist.c_str(), avgang.c_str()); + } +} + /** Print summary by parts for given hbond. */ void Action_HydrogenBond::summary_Parts(CpptrajFile* avgout, Hbond const& hb) const { for (unsigned int idx = 0; idx != hb.Nparts(); idx++) @@ -1354,8 +1368,11 @@ void Action_HydrogenBond::Print() { UU_Map_.clear(); // Sort and Print sort( HbondList.begin(), HbondList.end() ); - avgout_->Printf("%-*s %*s %*s %8s %12s %12s %12s\n", NUM, "#Acceptor", - NUM, "DonorH", NUM, "Donor", "Frames", "Frac", "AvgDist", "AvgAng"); + avgout_->Printf("%-*s %*s %*s %8s %12s %12s %12s", NUM, "#Acceptor", + NUM, "DonorH", NUM, "Donor", "Frames", "Frac", "AvgDist", "AvgAng"); + if (!splitFrames_.empty()) + summary_Parts_header(avgout_, splitFrames_.size()+1); + avgout_->Printf("\n"); for (Harray::const_iterator hbond = HbondList.begin(); hbond != HbondList.end(); ++hbond ) { double avg = ((double)hbond->Frames()) / ((double) Nframes_); @@ -1388,8 +1405,11 @@ void Action_HydrogenBond::Print() { sort( HbondList.begin(), HbondList.end() ); // Calc averages and print solvout_->Printf("#Solute-Solvent Hbonds:\n"); - solvout_->Printf("%-*s %*s %*s %8s %12s %12s %12s\n", NUM, "#Acceptor", - NUM, "DonorH", NUM, "Donor", "Count", "Frac", "AvgDist", "AvgAng"); + solvout_->Printf("%-*s %*s %*s %8s %12s %12s %12s", NUM, "#Acceptor", + NUM, "DonorH", NUM, "Donor", "Count", "Frac", "AvgDist", "AvgAng"); + if (!splitFrames_.empty()) + summary_Parts_header(solvout_, splitFrames_.size()+1); + solvout_->Printf("\n"); for (Harray::const_iterator hbond = HbondList.begin(); hbond != HbondList.end(); ++hbond ) { // Average has slightly diff meaning since for any given frame multiple diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index f73bc76698..75cb851192 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -43,6 +43,8 @@ class Action_HydrogenBond : public Action { void UpdateSeries(); /// Determine memory usage from # hbonds and time series std::string MemoryUsage(size_t, size_t, size_t) const; + /// Print header for summary of parts + static void summary_Parts_header(CpptrajFile*, unsigned int); /// Print summary of parts for hbond void summary_Parts(CpptrajFile*, Hbond const&) const; # ifdef MPI diff --git a/test/Test_Hbond_Split/avghb.dat.save b/test/Test_Hbond_Split/avghb.dat.save index e9d3e55370..8b0bab2b03 100644 --- a/test/Test_Hbond_Split/avghb.dat.save +++ b/test/Test_Hbond_Split/avghb.dat.save @@ -1,4 +1,4 @@ -#Acceptor DonorH Donor Frames Frac AvgDist AvgAng +#Acceptor DonorH Donor Frames Frac AvgDist AvgAng Frames1 Frac1 AvgDist1 AvgAng1 Frames2 Frac2 AvgDist2 AvgAng2 ILE_19@O THR_12@H THR_12@N 69 0.6900 2.8626 160.5248 37 0.3700 2.8638 160.2253 32 0.3200 2.8612 160.8711 VAL_2@O GLU_13@H GLU_13@N 59 0.5900 2.8517 159.8143 36 0.3600 2.8510 158.5372 23 0.2300 2.8527 161.8133 TYR_11@O ILE_4@H ILE_4@N 53 0.5300 2.8682 160.4581 31 0.3100 2.8728 162.3399 22 0.2200 2.8616 157.8065 diff --git a/test/Test_Hbond_Split/solvavg.dat.save b/test/Test_Hbond_Split/solvavg.dat.save index 0e00cdf0f7..3a309b7d69 100644 --- a/test/Test_Hbond_Split/solvavg.dat.save +++ b/test/Test_Hbond_Split/solvavg.dat.save @@ -1,5 +1,5 @@ #Solute-Solvent Hbonds: -#Acceptor DonorH Donor Count Frac AvgDist AvgAng +#Acceptor DonorH Donor Count Frac AvgDist AvgAng Frames1 Frac1 AvgDist1 AvgAng1 Frames2 Frac2 AvgDist2 AvgAng2 GLU_5@OE1 SolventH SolventDnr 24 2.4000 2.7018 162.8196 8 0.8000 2.7293 160.2363 16 1.6000 2.6881 164.1113 GLU_5@OE2 SolventH SolventDnr 17 1.7000 2.6891 158.7293 7 0.7000 2.6849 158.4141 10 1.0000 2.6920 158.9500 GLY_7@O SolventH SolventDnr 15 1.5000 2.7616 157.9700 5 0.5000 2.7877 156.6524 10 1.0000 2.7486 158.6288 From f9fac8ed506e791cba43444f62e9f1450111a8ef Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 09:21:33 -0500 Subject: [PATCH 12/34] Fix clean rules --- test/Test_Hbond_Split/RunTest.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/Test_Hbond_Split/RunTest.sh b/test/Test_Hbond_Split/RunTest.sh index 97d10dc985..d30690a8d1 100755 --- a/test/Test_Hbond_Split/RunTest.sh +++ b/test/Test_Hbond_Split/RunTest.sh @@ -3,7 +3,7 @@ . ../MasterTest.sh # Clean -CleanFiles hbond.in nhb.dat avghb.dat +CleanFiles hbond.in nhb.dat avghb.dat solvhb.dat solvavg.dat INPUT="-i hbond.in" From 13ffeb5ea989db375082529dab0066be489fd7b4 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 09:31:43 -0500 Subject: [PATCH 13/34] Add some error checking and help text --- src/Action_HydrogenBond.cpp | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index fc86b27db1..f2eb887099 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -49,6 +49,7 @@ void Action_HydrogenBond::Help() const { "\t[solventdonor ] [solventacceptor ]\n" "\t[solvout ] [bridgeout ] [bridgebyatom]\n" "\t[series [uuseries ] [uvseries ]]\n" + "\t[splitframe \n" " Hydrogen bond is defined as A-HD, where A is acceptor heavy atom, H is\n" " hydrogen, D is donor heavy atom. Hydrogen bond is formed when\n" " A to D distance < dcut and A-H-D angle > acut; if acut < 0 it is ignored.\n" @@ -60,7 +61,10 @@ void Action_HydrogenBond::Help() const { " searched for in .\n" " If both donormask and acceptor mask are specified no automatic searching will occur.\n" " If donorhmask is specified atoms in that mask will be paired with atoms in\n" - " donormask instead of automatically searching for hydrogen atoms.\n"); + " donormask instead of automatically searching for hydrogen atoms.\n" + " The 'splitframe' keyword can be used to divide the average hydrogen bond\n" + " analysis into parts, e.g. 'splitframe 250,500,1000' will divide analysis\n" + " into 1-249,250-499,500-999, etc.\n"); } // Action_HydrogenBond::Init() @@ -99,7 +103,16 @@ Action::RetType Action_HydrogenBond::Init(ArgList& actionArgs, ActionInit& init, } int sf = splits.getNextInteger(-1); // User frame #s start at 1 while (sf > 0) { - splitFrames_.push_back( sf - 1 ); + // Since user frame #s start at 1, subtract 1 for actual frame index. + sf--; + // Check that frame arg is valid + if (!splitFrames_.empty()) { + if (sf <= splitFrames_.back()) { + mprinterr("Error: 'splitframe' #s must be in increasing order.\n"); + return Action::ERR; + } + } + splitFrames_.push_back( sf ); sf = splits.getNextInteger(-1); } if ((int)splitFrames_.size() < splits.Nargs()) { From c682c13141108ab45f3636c1200d300fd53790dc Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 09:44:28 -0500 Subject: [PATCH 14/34] Start adding split for bridge --- src/Action_HydrogenBond.cpp | 18 +++++++++++------- src/Action_HydrogenBond.h | 22 ++++++++++++++++------ 2 files changed, 27 insertions(+), 13 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index f2eb887099..bbdd0e6768 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -941,7 +941,10 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { for (RmapType::const_iterator bridge = solvent2solute_.begin(); bridge != solvent2solute_.end(); ++bridge) { - // If solvent molecule is bound to 2 or more different residues, + // bridge->first is solvent residue number. + // bridge->second is a set of solute residue numbers the solvent + // residue is bound to. + // If solvent molecule is bound to 2 or more different solute residues, // it is bridging. if ( bridge->second.size() > 1) { bool isBridge = true; @@ -966,25 +969,26 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { isBridge = (nequal < bridge->second.size()); } if (isBridge) { + // numHB is used to track the number of bridges ++numHB; // Bridging Solvent residue number bridgeID.append(integerToString( bridge->first+1 ) + "("); + // Loop over solute residues this solvent is bound to. for (std::set::const_iterator res = bridge->second.begin(); res != bridge->second.end(); ++res) // Solute residue number being bridged bridgeID.append( integerToString( *res+1 ) + "+" ); bridgeID.append("),"); - // Find bridge in map based on this combo of residues (bridge.second) + // Find bridge in map based on this combo of residues (bridge->second) BmapType::iterator b_it = BridgeMap_.lower_bound( bridge->second ); if (b_it == BridgeMap_.end() || b_it->first != bridge->second) // New Bridge - BridgeMap_.insert( b_it, std::pair,Bridge>(bridge->second, Bridge()) ); - else - // Increment bridge #frames - b_it->second.Update(frameNum); + b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(bridge->second, Bridge(splitFrames_)) ); + // Increment bridge #frames + b_it->second.Update(frameNum, splitFrames_, frm.TrajoutNum()); } } - } + } // END LOOP OVER solvent2solute_ if (bridgeID.empty()) bridgeID.assign("None"); NumBridge_->Add(frameNum, &numHB); diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 75cb851192..a3ef9ca458 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -264,8 +264,11 @@ class Action_HydrogenBond::Hbond { /// Track solvent bridge between 2 or more solute residues. class Action_HydrogenBond::Bridge { public: - /// Constructor - new bridge - Bridge() : data_(0), frames_(1) {} + /// CONSTRUCTOR - new bridge + Bridge(Iarray const& splits) : data_(0), frames_(0) { + if (!splits.empty()) + partsFrames_.assign(splits.size()+1, 0); + } # ifdef MPI /// Constructor - new bridge with given # frames Bridge(int f) : data_(0), frames_(f) {} @@ -274,9 +277,15 @@ class Action_HydrogenBond::Bridge { # endif int Frames() const { return frames_; } /// Update frames/time series - void Update(int f) { + void Update(int fnum, Iarray const& splitFrames, int onum) { ++frames_; - if (data_ != 0) data_->AddVal(f, 1); + if (data_ != 0) data_->AddVal(fnum, 1); + if (!splitFrames.empty()) { + // Find the correct part NOTE assumes onum never out of range + unsigned int part = 0; + while (part < splitFrames.size() && onum >= splitFrames[part]) part++; + partsFrames_[part]++; + } } /// \return true if bridge has more frames than rhs. bool operator()(Bridge const& rhs) const { @@ -288,7 +297,8 @@ class Action_HydrogenBond::Bridge { return false; } private: - DataSet_integer* data_; ///< Hold time series data - int frames_; ///< # frames this bridge has been present. + DataSet_integer* data_; ///< Hold time series data TODO + int frames_; ///< # frames this bridge has been present. + Iarray partsFrames_; ///< Hold # frames bridge present for each part. }; #endif From 91efb77e45e6ea9d41ee67cdc1cb388a8ae68888 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 10:13:07 -0500 Subject: [PATCH 15/34] Add parts to bridge printout --- src/Action_HydrogenBond.cpp | 46 ++++++++++++++++++++++++++++++------- src/Action_HydrogenBond.h | 16 ++++--------- 2 files changed, 42 insertions(+), 20 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index bbdd0e6768..b246fd5a45 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1342,6 +1342,28 @@ void Action_HydrogenBond::summary_Parts(CpptrajFile* avgout, Hbond const& hb) co hb.PartDist(idx).mean(), hb.PartAngle(idx).mean()*Constants::RADDEG); } +/** Used to associate Bridge with solute atoms/residues and sort Bridges by frames. */ +class Action_HydrogenBond::bridgeSorter { + public: + /// CONSTRUCTOR - list of solute atoms/residues, Bridge info + bridgeSorter(std::set const& uIdx, Bridge const& bridge) : + uIdx_(uIdx), bridge_(bridge) {} + /// Used to sort by bridge # frames + bool operator<(bridgeSorter const& rhs) const { + if (bridge_.Frames() == rhs.bridge_.Frames()) + return (uIdx_ < rhs.uIdx_); + else + return (bridge_.Frames() > rhs.bridge_.Frames()); + } + /// \return List of solute atom/residue #s + std::set const& Uidx() const { return uIdx_; } + /// \return Bridging info + Bridge const& Binfo() const { return bridge_; } + private: + std::set uIdx_; ///< Hold solute atom/residue #s + Bridge bridge_; ///< Hold bridging info. +}; + // Action_HydrogenBond::Print() /** Print average occupancies over all frames for all detected Hbonds. */ void Action_HydrogenBond::Print() { @@ -1458,6 +1480,7 @@ void Action_HydrogenBond::Print() { } HbondList.clear(); } + // BRIDGING INFO if (bridgeout_ != 0 && calcSolvent_) { if (bridgeByAtom_) @@ -1465,26 +1488,33 @@ void Action_HydrogenBond::Print() { else bridgeout_->Printf("#Bridging Solute Residues:\n"); // Place bridging values in a vector for sorting - typedef std::vector Bvec; + typedef std::vector Bvec; Bvec bridgevector; + bridgevector.reserve( BridgeMap_.size() ); for (BmapType::const_iterator it = BridgeMap_.begin(); it != BridgeMap_.end(); ++it) - bridgevector.push_back( Bpair(it->first, it->second.Frames()) ); - std::sort( bridgevector.begin(), bridgevector.end(), bridge_cmp() ); + bridgevector.push_back( bridgeSorter(it->first, it->second) ); + std::sort( bridgevector.begin(), bridgevector.end() ); for (Bvec::const_iterator bv = bridgevector.begin(); bv != bridgevector.end(); ++bv) { if (bridgeByAtom_) { bridgeout_->Printf("Bridge Atm"); - for (std::set::const_iterator atm = bv->first.begin(); - atm != bv->first.end(); ++atm) + for (std::set::const_iterator atm = bv->Uidx().begin(); + atm != bv->Uidx().end(); ++atm) bridgeout_->Printf(" %s", CurrentParm_->TruncAtomNameNum(*atm).c_str()); } else { bridgeout_->Printf("Bridge Res"); - for (std::set::const_iterator res = bv->first.begin(); - res != bv->first.end(); ++res) + for (std::set::const_iterator res = bv->Uidx().begin(); + res != bv->Uidx().end(); ++res) bridgeout_->Printf(" %i:%s", *res+1, CurrentParm_->Res( *res ).Name().Formatted(4).c_str()); } - bridgeout_->Printf(", %i frames.\n", bv->second); + bridgeout_->Printf(", %i frames.", bv->Binfo().Frames()); + if (!splitFrames_.empty()) { + bridgeout_->Printf(" Parts:"); + for (unsigned int idx = 0; idx != bv->Binfo().Nparts(); idx++) + bridgeout_->Printf(" %i", bv->Binfo().PartFrames(idx)); + } + bridgeout_->Printf("\n"); } } } diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index a3ef9ca458..f08164a446 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -29,6 +29,7 @@ class Action_HydrogenBond : public Action { class Site; class Hbond; class Bridge; + class bridgeSorter; inline double Angle(const double*, const double*, const double*, Box const&) const; inline int UU_Set_Idx(int,int) const; @@ -118,18 +119,6 @@ class Action_HydrogenBond : public Action { bool hasSolventAcceptor_; ///< If true a solvent acceptor mask was specified. bool calcSolvent_; ///< If true solute-solvent hbonds and bridges will be calcd. bool bridgeByAtom_; ///< If true determine bridging by atom. - // TODO replace with class - typedef std::pair< std::set,int > Bpair; - /// \return true if 1) p0 frames > p1 frames, 2) p0 res < p1 res - struct bridge_cmp { - inline bool operator()(Bpair const& p0, Bpair const& p1) const - { - if (p0.second == p1.second) - return (p0.first < p1.first); - else - return (p0.second > p1.second); - } - }; }; // ----- CLASSES --------------------------------------------------------------- @@ -296,6 +285,9 @@ class Action_HydrogenBond::Bridge { else return false; } + // Summary by parts + unsigned int Nparts() const { return partsFrames_.size(); } + unsigned int PartFrames(unsigned int idx) const { return (unsigned int)partsFrames_[idx]; } private: DataSet_integer* data_; ///< Hold time series data TODO int frames_; ///< # frames this bridge has been present. From 88d9566843b33ba4765a80aaa920bf17e336ffd1 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 10:13:42 -0500 Subject: [PATCH 16/34] Add parts to bridge printout --- test/Test_Hbond_Split/solvavg.dat.save | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/Test_Hbond_Split/solvavg.dat.save b/test/Test_Hbond_Split/solvavg.dat.save index 3a309b7d69..dd2db3f157 100644 --- a/test/Test_Hbond_Split/solvavg.dat.save +++ b/test/Test_Hbond_Split/solvavg.dat.save @@ -37,8 +37,8 @@ SolventAcc TRP_4@HE1 TRP_4@NE1 1 0.1000 2.7128 SolventAcc SER_1@H1 SER_1@N 1 0.1000 2.9279 139.1256 0 0.0000 0.0000 0.0000 1 0.1000 2.9279 139.1256 SolventAcc TRP_11@H TRP_11@N 1 0.1000 2.9678 166.9306 0 0.0000 0.0000 0.0000 1 0.1000 2.9678 166.9306 #Bridging Solute Residues: -Bridge Res 10:THR 11:TRP , 5 frames. -Bridge Res 5:GLU 12:LYS , 4 frames. -Bridge Res 2:TRP 3:THR , 2 frames. -Bridge Res 3:THR 5:GLU , 2 frames. -Bridge Res 5:GLU 10:THR , 2 frames. +Bridge Res 10:THR 11:TRP , 5 frames. Parts: 4 1 +Bridge Res 5:GLU 12:LYS , 4 frames. Parts: 2 2 +Bridge Res 2:TRP 3:THR , 2 frames. Parts: 2 0 +Bridge Res 3:THR 5:GLU , 2 frames. Parts: 2 0 +Bridge Res 5:GLU 10:THR , 2 frames. Parts: 0 2 From 9ff2b506a16f6b68e356c62b4fc00685b637ea0e Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 10:28:57 -0500 Subject: [PATCH 17/34] Fix in parallel --- src/Action_HydrogenBond.cpp | 29 ++++++++++++----------------- src/Action_HydrogenBond.h | 13 +++++++++++++ 2 files changed, 25 insertions(+), 17 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index b246fd5a45..87fee30743 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1192,22 +1192,10 @@ int Action_HydrogenBond::SyncAction() { Ivals.reserve( ntotal * 4 ); // Store UU bonds in flat arrays. for (UUmapType::const_iterator it = UU_Map_.begin(); it != UU_Map_.end(); ++it) { - /*Dvals.push_back( it->second.Dist() ); - Dvals.push_back( it->second.Angle() ); - Ivals.push_back( it->second.A() ); - Ivals.push_back( it->second.H() ); - Ivals.push_back( it->second.D() ); - Ivals.push_back( it->second.Frames() );*/ HbondToArray(Dvals, Ivals, it->second); } // Store UV bonds in flat arrays for (UVmapType::const_iterator it = UV_Map_.begin(); it != UV_Map_.end(); ++it) { - /*Dvals.push_back( it->second.Dist() ); - Dvals.push_back( it->second.Angle() ); - Ivals.push_back( it->second.A() ); - Ivals.push_back( it->second.H() ); - Ivals.push_back( it->second.D() ); - Ivals.push_back( it->second.Frames() );*/ HbondToArray(Dvals, Ivals, it->second); } trajComm_.Send( &(Dvals[0]), Dvals.size(), MPI_DOUBLE, 0, 1300 ); @@ -1233,7 +1221,7 @@ int Action_HydrogenBond::SyncAction() { if (calcSolvent_) { // Sync bridging data - // iArray will contain for each bridge: Nres, res1, ..., resN, Frames + // iArray will contain for each bridge: Nres, res1, ..., resN, Frames, Npart1, ..., NpartN std::vector iArray; int iSize; if (trajComm_.Master()) { @@ -1251,13 +1239,16 @@ int Action_HydrogenBond::SyncAction() { for (int ir = 0; ir != iArray[idx]; ir++, i2++) residues.insert( iArray[i2] ); BmapType::iterator b_it = BridgeMap_.lower_bound( residues ); - if (b_it == BridgeMap_.end() || b_it->first != residues ) + if (b_it == BridgeMap_.end() || b_it->first != residues ) { // New Bridge - BridgeMap_.insert( b_it, std::pair,Bridge>(residues, Bridge(iArray[i2])) ); - else + b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(residues, Bridge(iArray[i2])) ); + b_it->second.SetupParts(nParts, &iArray[0] + i2 + 1); + } else { // Increment bridge #frames b_it->second.Combine( iArray[i2] ); - idx = i2 + 1; + b_it->second.CombineParts(nParts, &iArray[0] + i2 + 1); + } + idx = i2 + 1 + nParts; } } } else { @@ -1268,6 +1259,10 @@ int Action_HydrogenBond::SyncAction() { for ( std::set::const_iterator r = b->first.begin(); r != b->first.end(); ++r) iArray.push_back( *r ); // Bridging res iArray.push_back( b->second.Frames() ); // # frames + if (nParts > 0) { + for (unsigned int part = 0; part != nParts; part++) + iArray.push_back( b->second.PartFrames(part) ); + } } // Since the size of each bridge can be different (i.e. differing #s of // residues may be bridged), first send size of the transport array. diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index f08164a446..f4fbea7518 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -263,6 +263,19 @@ class Action_HydrogenBond::Bridge { Bridge(int f) : data_(0), frames_(f) {} /// Increment number of frames void Combine(int n) { frames_ += n; } + /// Set up parts with given int array containing # frames for each part + void SetupParts(unsigned int nparts, const int* ivals) { + if (nparts == 0) return; + partsFrames_.clear(); + partsFrames_.reserve( nparts ); + for (unsigned int idx = 0; idx != nparts; idx++) + partsFrames_.push_back( ivals[idx] ); + } + /// Update parts with given int array containing # frames for each part + void CombineParts(unsigned int nparts, const int* ivals) { + for (unsigned int idx = 0; idx != nparts; idx++) + partsFrames_[idx] += ivals[idx]; + } # endif int Frames() const { return frames_; } /// Update frames/time series From 06414b5e105d764f7ebc1fbbbd616bc7a13f5b3e Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 10:34:21 -0500 Subject: [PATCH 18/34] Clean up code doc --- src/Action_HydrogenBond.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 87fee30743..8a68fc1a0b 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1075,7 +1075,7 @@ int Action_HydrogenBond::SyncAction() { } // Need to send hbond data from all ranks to master. - std::vector Dvals; // Hold dist_ and angle_ for each hbond, as well as n_/mean_/M2_ for dist/angle each part + std::vector Dvals; // Hold dist_ and angle_ for each hbond, (as well as n_/mean_/M2_ for dist/angle each part) std::vector Ivals; // Hold A_, H_, D_, and frames_ for each hbond unsigned int dvalsPerHbond; unsigned int nParts; @@ -1221,10 +1221,11 @@ int Action_HydrogenBond::SyncAction() { if (calcSolvent_) { // Sync bridging data - // iArray will contain for each bridge: Nres, res1, ..., resN, Frames, Npart1, ..., NpartN + // iArray will contain for each bridge: Nres, res1, ..., resN, Frames[, Npart1, ..., NpartN] std::vector iArray; int iSize; if (trajComm_.Master()) { + // MASTER RANK for (int rank = 1; rank < trajComm_.Size(); rank++) { // Receive size of iArray @@ -1252,6 +1253,7 @@ int Action_HydrogenBond::SyncAction() { } } } else { + // NON-MASTER // Construct bridge info array. for (BmapType::const_iterator b = BridgeMap_.begin(); b != BridgeMap_.end(); ++b) { @@ -1270,7 +1272,7 @@ int Action_HydrogenBond::SyncAction() { trajComm_.Send( &iSize, 1, MPI_INT, 0, 1302 ); trajComm_.Send( &(iArray[0]), iSize, MPI_INT, 0, 1303 ); } - } + } // END COMMUNICATING BRIDGE DATA TO MASTER return 0; } #endif From 005bf5f55ba602b594a9ef8f7bf4c57aa6f98a04 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 10:52:35 -0500 Subject: [PATCH 19/34] Improve code docs --- src/Action_HydrogenBond.cpp | 19 +++++++++++++++++-- src/Action_HydrogenBond.h | 10 +++++++++- 2 files changed, 26 insertions(+), 3 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 8a68fc1a0b..69b11b9825 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -577,6 +577,7 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { } // Action_HydrogenBond::Angle() +/** Calculate angle between 3 atoms, optionally with imaging. */ double Action_HydrogenBond::Angle(const double* XA, const double* XH, const double* XD, Box const& boxIn) const { if (imageOpt_.ImagingType() == ImageOption::NO_IMAGE) @@ -649,6 +650,9 @@ void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,i } // Action_HydrogenBond::CalcSolvHbonds() +/** Calculate hydrogen bonds between solute site and solvent acceptor, + * or solvent site and solute acceptor. + */ void Action_HydrogenBond::CalcSolvHbonds(int frameNum, double dist2, Site const& SiteD, const double* XYZD, int a_atom, const double* XYZA, @@ -692,6 +696,7 @@ int Action_HydrogenBond::UU_Set_Idx(int a_atom, int h_atom) const { } // Action_HydrogenBond::UUset() +/** \return solute-solute hydrogen bond time series set with legend set. */ DataSet_integer* Action_HydrogenBond::UUset(int a_atom, int h_atom, int d_atom) { std::string hblegend = CurrentParm_->TruncResAtomName(a_atom) + "-" + CurrentParm_->TruncResAtomName(d_atom) + "-" + @@ -704,6 +709,7 @@ DataSet_integer* Action_HydrogenBond::UUset(int a_atom, int h_atom, int d_atom) } // Action_HydrogenBond::AddUU() +/** Add or update a solute-solute hydrogen bond with given angle/distance. */ void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, int h_atom, int d_atom, int onum) { // Index UU hydrogen bonds by DonorH-Acceptor @@ -725,6 +731,11 @@ void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, } // Action_HydrogenBond::CalcSiteHbonds() +/** Calculate hydrogen bonds between given solute donor site and + * solute acceptor atom. + * The distance cutoff should already be satisfied between donor and + * acceptor heavy atoms. + */ void Action_HydrogenBond::CalcSiteHbonds(int frameNum, double dist2, Site const& SiteD, const double* XYZD, int a_atom, const double* XYZA, @@ -1020,6 +1031,7 @@ static inline std::string CreateHBlegend(Topology const& topIn, int a_atom, int } // Action_HydrogenBond::GetRankNhbonds() +/** Determine how many hydrogen bonds are on each rank. */ std::vector Action_HydrogenBond::GetRankNhbonds( int num_hb, Parallel::Comm const& commIn ) { std::vector nhb_on_rank; @@ -1275,9 +1287,12 @@ int Action_HydrogenBond::SyncAction() { } // END COMMUNICATING BRIDGE DATA TO MASTER return 0; } -#endif +#endif /* MPI */ // Action_HydrogenBond::UpdateSeries() +/** Ensure all time series data is up-to-date with Nframes. + * Should only be called once. + */ void Action_HydrogenBond::UpdateSeries() { if (seriesUpdated_) return; if (series_ && Nframes_ > 0) { @@ -1286,11 +1301,11 @@ void Action_HydrogenBond::UpdateSeries() { for (UVmapType::iterator hb = UV_Map_.begin(); hb != UV_Map_.end(); ++hb) hb->second.FinishSeries(Nframes_); } - // Should only be called once. seriesUpdated_ = true; } // Action_Hbond::MemoryUsage() +/** Estimate the memory usage of the hbond command. */ std::string Action_HydrogenBond::MemoryUsage(size_t n_uu_pairs, size_t n_uv_pairs, size_t nFrames) const { static const size_t sizeHbond = sizeof(Hbond); diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index f4fbea7518..b25f4ef2d8 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -30,14 +30,20 @@ class Action_HydrogenBond : public Action { class Hbond; class Bridge; class bridgeSorter; - + /// Calcuate angle with optional imaging. inline double Angle(const double*, const double*, const double*, Box const&) const; + /// \return solute-solute hydrogen bond index based on donor/acceptor atoms inline int UU_Set_Idx(int,int) const; + /// \return solute-solute hydrogen bond time series DataSet. inline DataSet_integer* UUset(int,int,int); + /// Add/update solute-solute hydrogen bond void AddUU(double,double,int,int,int,int,int); + /// Add/update solute-solvent hydrogen bond void AddUV(double,double,int,int,int,int,bool,int); + /// Calculate UU hydrogen bonds between site and acceptor atom. void CalcSiteHbonds(int,double,Site const&,const double*,int,const double*, Frame const&, int&); + /// Calculate UV hydrogen bonds between solute/solvent site and solute/solvent acceptor atom. void CalcSolvHbonds(int,double,Site const&,const double*,int,const double*, Frame const&, int&, bool); /// Update all hydrogen bond time series @@ -49,7 +55,9 @@ class Action_HydrogenBond : public Action { /// Print summary of parts for hbond void summary_Parts(CpptrajFile*, Hbond const&) const; # ifdef MPI + /// \return Array containing # hbonds on each rank static std::vector GetRankNhbonds(int,Parallel::Comm const&); + /// Flatten given Hbond into arrays static void HbondToArray(std::vector&, std::vector&, Hbond const&); # endif From 87517fa2b2673dd5135e268d557e9ede5f7c5ca2 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 10:57:36 -0500 Subject: [PATCH 20/34] Use CreateHBlegend function to create all hb legends --- src/Action_HydrogenBond.cpp | 35 ++++++++++++++++------------------- 1 file changed, 16 insertions(+), 19 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 69b11b9825..fee4da085f 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -266,6 +266,19 @@ inline bool IsFON(Atom const& atm) { atm.Element() == Atom::NITROGEN); } +/** Create legend for hydrogen bond based on given atoms. */ +static inline std::string CreateHBlegend(Topology const& topIn, int a_atom, int h_atom, int d_atom) +{ + if (a_atom == -1) + return (topIn.TruncResAtomName(h_atom) + "-V"); + else if (d_atom == -1) + return (topIn.TruncResAtomName(a_atom) + "-V"); + else + return (topIn.TruncResAtomName(a_atom) + "-" + + topIn.TruncResAtomName(d_atom) + "-" + + topIn[h_atom].Name().Truncated()); +} + // Action_HydrogenBond::Setup() Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { CurrentParm_ = setup.TopAddress(); @@ -636,10 +649,10 @@ void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,i } Hbond hb; if (udonor) { // Do not care about which solvent acceptor - if (ds != 0) ds->SetLegend( CurrentParm_->TruncResAtomName(h_atom) + "-V" ); + if (ds != 0) ds->SetLegend( CreateHBlegend(*CurrentParm_, -1, h_atom, d_atom) ); hb = Hbond(ds, -1, h_atom, d_atom, splitFrames_); } else { // Do not care about which solvent donor - if (ds != 0) ds->SetLegend( CurrentParm_->TruncResAtomName(a_atom) + "-V" ); + if (ds != 0) ds->SetLegend( CreateHBlegend(*CurrentParm_, a_atom, -1, -1) ); hb = Hbond(ds, a_atom, -1, -1, splitFrames_); } it = UV_Map_.insert(it, std::pair(hbidx,hb)); @@ -698,13 +711,10 @@ int Action_HydrogenBond::UU_Set_Idx(int a_atom, int h_atom) const { // Action_HydrogenBond::UUset() /** \return solute-solute hydrogen bond time series set with legend set. */ DataSet_integer* Action_HydrogenBond::UUset(int a_atom, int h_atom, int d_atom) { - std::string hblegend = CurrentParm_->TruncResAtomName(a_atom) + "-" + - CurrentParm_->TruncResAtomName(d_atom) + "-" + - (*CurrentParm_)[h_atom].Name().Truncated(); DataSet_integer* ds = (DataSet_integer*) masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"solutehb",UU_Set_Idx(a_atom,h_atom))); if (UUseriesout_ != 0) UUseriesout_->AddDataSet( ds ); - ds->SetLegend( hblegend ); + ds->SetLegend( CreateHBlegend(*CurrentParm_, a_atom, h_atom, d_atom) ); return ds; } @@ -1017,19 +1027,6 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { } #ifdef MPI -// TODO Use in other hbond functions -static inline std::string CreateHBlegend(Topology const& topIn, int a_atom, int h_atom, int d_atom) -{ - if (a_atom == -1) - return (topIn.TruncResAtomName(h_atom) + "-V"); - else if (d_atom == -1) - return (topIn.TruncResAtomName(a_atom) + "-V"); - else - return (topIn.TruncResAtomName(a_atom) + "-" + - topIn.TruncResAtomName(d_atom) + "-" + - topIn[h_atom].Name().Truncated()); -} - // Action_HydrogenBond::GetRankNhbonds() /** Determine how many hydrogen bonds are on each rank. */ std::vector Action_HydrogenBond::GetRankNhbonds( int num_hb, Parallel::Comm const& commIn ) From 3105115178e595fc652dc67c9f9f74d0082f3368 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 11:05:07 -0500 Subject: [PATCH 21/34] Start adding time series for bridges --- src/Action_HydrogenBond.cpp | 14 +++++++++++--- src/Action_HydrogenBond.h | 3 ++- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index fee4da085f..160b86949d 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -29,6 +29,7 @@ Action_HydrogenBond::Action_HydrogenBond() : Nframes_(0), debug_(0), series_(false), + Bseries_(false), seriesUpdated_(false), useAtomNum_(false), noIntramol_(false), @@ -1002,9 +1003,16 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { bridgeID.append("),"); // Find bridge in map based on this combo of residues (bridge->second) BmapType::iterator b_it = BridgeMap_.lower_bound( bridge->second ); - if (b_it == BridgeMap_.end() || b_it->first != bridge->second) - // New Bridge - b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(bridge->second, Bridge(splitFrames_)) ); + if (b_it == BridgeMap_.end() || b_it->first != bridge->second) { + // New Bridge + DataSet_integer* bds = 0; + if (Bseries_) { + bds = (DataSet_integer*) + masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"bridge",BridgeMap_.size())); + //if (Bseriesout_ != 0) Bseriesout_->AddDataSet( bds ); + } + b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(bridge->second, Bridge(bds, splitFrames_)) ); + } // Increment bridge #frames b_it->second.Update(frameNum, splitFrames_, frm.TrajoutNum()); } diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index b25f4ef2d8..ebf9a92c18 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -117,6 +117,7 @@ class Action_HydrogenBond : public Action { int Nframes_; ///< Number of frames action has been active int debug_; bool series_; ///< If true track hbond time series. + bool Bseries_; ///< If true track bridge time series. bool seriesUpdated_; ///< If false hbond time series need to be finished. bool useAtomNum_; ///< If true include atom numbers in labels/legends bool noIntramol_; ///< If true ignore intramolecular hydrogen bonds/bridges. @@ -262,7 +263,7 @@ class Action_HydrogenBond::Hbond { class Action_HydrogenBond::Bridge { public: /// CONSTRUCTOR - new bridge - Bridge(Iarray const& splits) : data_(0), frames_(0) { + Bridge(DataSet_integer* bds, Iarray const& splits) : data_(bds), frames_(0) { if (!splits.empty()) partsFrames_.assign(splits.size()+1, 0); } From b443c156256dd6ccc5c50659026a7d7b2f80f915 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 11:14:54 -0500 Subject: [PATCH 22/34] Add bridging time series --- src/Action_HydrogenBond.cpp | 19 ++++++++++++++++++- src/Action_HydrogenBond.h | 3 ++- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 160b86949d..f3e7247902 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -20,6 +20,7 @@ Action_HydrogenBond::Action_HydrogenBond() : BridgeID_(0), UUseriesout_(0), UVseriesout_(0), + Bseriesout_(0), avgout_(0), solvout_(0), bridgeout_(0), @@ -50,6 +51,7 @@ void Action_HydrogenBond::Help() const { "\t[solventdonor ] [solventacceptor ]\n" "\t[solvout ] [bridgeout ] [bridgebyatom]\n" "\t[series [uuseries ] [uvseries ]]\n" + "\t[bseries [bseriesfile ]]\n" "\t[splitframe \n" " Hydrogen bond is defined as A-HD, where A is acceptor heavy atom, H is\n" " hydrogen, D is donor heavy atom. Hydrogen bond is formed when\n" @@ -84,6 +86,11 @@ Action::RetType Action_HydrogenBond::Init(ArgList& actionArgs, ActionInit& init, UVseriesout_ = init.DFL().AddDataFile(actionArgs.GetStringKey("uvseries"), actionArgs); init.DSL().SetDataSetsPending(true); } + Bseries_ = actionArgs.hasKey("bseries"); + if (Bseries_) { + Bseriesout_ = init.DFL().AddDataFile(actionArgs.GetStringKey("bseriesfile"), actionArgs); + init.DSL().SetDataSetsPending(true); + } std::string avgname = actionArgs.GetStringKey("avgout"); std::string solvname = actionArgs.GetStringKey("solvout"); if (solvname.empty()) solvname = avgname; @@ -252,6 +259,11 @@ Action::RetType Action_HydrogenBond::Init(ArgList& actionArgs, ActionInit& init, if (UVseriesout_ != 0) mprintf("\tWriting solute-solvent time series to %s\n", UVseriesout_->DataFilename().full()); } + if (Bseries_) { + mprintf("\tTime series data for each bridge will be saved for analysis.\n"); + if (Bseriesout_ != 0) + mprintf("\tWriting bridge time series to '%s'\n", Bseriesout_->DataFilename().full()); + } if (imageOpt_.UseImage()) mprintf("\tImaging enabled.\n"); masterDSL_ = init.DslPtr(); @@ -1009,7 +1021,12 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { if (Bseries_) { bds = (DataSet_integer*) masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"bridge",BridgeMap_.size())); - //if (Bseriesout_ != 0) Bseriesout_->AddDataSet( bds ); + // Create a legend from the indices. + std::string blegend("B"); + for (std::set::const_iterator brs = bridge->second.begin(); brs != bridge->second.end(); ++brs) + blegend.append("_" + integerToString(*brs + 1)); + bds->SetLegend( blegend ); + if (Bseriesout_ != 0) Bseriesout_->AddDataSet( bds ); } b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(bridge->second, Bridge(bds, splitFrames_)) ); } diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index ebf9a92c18..e563799749 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -107,7 +107,8 @@ class Action_HydrogenBond : public Action { DataSet* NumBridge_; ///< Hold # solute-solvent bridges per frame. DataSet* BridgeID_; ///< Hold info on each bridge per frame. DataFile* UUseriesout_; ///< File to write UU time series to. - DataFile* UVseriesout_; ///< File to write UN time series to. + DataFile* UVseriesout_; ///< File to write UV time series to. + DataFile* Bseriesout_; ///< File to write bridge time series to. CpptrajFile* avgout_; ///< File to write UU averages to. CpptrajFile* solvout_; ///< File to write UV averages to. CpptrajFile* bridgeout_; ///< File to write bridge totals to. From 6ed20645c11f15344fab898d2910587efe4eb314 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 11:15:30 -0500 Subject: [PATCH 23/34] Add bridges time series test --- test/Test_Hbond_Split/RunTest.sh | 6 ++++-- test/Test_Hbond_Split/bridges.dat.save | 11 +++++++++++ 2 files changed, 15 insertions(+), 2 deletions(-) create mode 100644 test/Test_Hbond_Split/bridges.dat.save diff --git a/test/Test_Hbond_Split/RunTest.sh b/test/Test_Hbond_Split/RunTest.sh index d30690a8d1..196e2ce779 100755 --- a/test/Test_Hbond_Split/RunTest.sh +++ b/test/Test_Hbond_Split/RunTest.sh @@ -3,7 +3,7 @@ . ../MasterTest.sh # Clean -CleanFiles hbond.in nhb.dat avghb.dat solvhb.dat solvavg.dat +CleanFiles hbond.in nhb.dat avghb.dat solvhb.dat solvavg.dat bridges.dat INPUT="-i hbond.in" @@ -33,10 +33,12 @@ TestUVsplit() { parm ../tz2.ortho.parm7 trajin ../tz2.ortho.nc hbond hb out solvhb.dat :1-13 solventacceptor :WAT@O solventdonor :WAT \ - solvout solvavg.dat bridgeout solvavg.dat splitframe 5 + solvout solvavg.dat bridgeout solvavg.dat splitframe 5 \ + bseries bseriesfile bridges.dat EOF RunCpptraj "$UNITNAME" DoTest solvavg.dat.save solvavg.dat + DoTest bridges.dat.save bridges.dat fi } diff --git a/test/Test_Hbond_Split/bridges.dat.save b/test/Test_Hbond_Split/bridges.dat.save new file mode 100644 index 0000000000..d2eab49f45 --- /dev/null +++ b/test/Test_Hbond_Split/bridges.dat.save @@ -0,0 +1,11 @@ +#Frame B_10_11 B_2_3 B_3_5 B_5_12 B_5_10 + 1 1 1 1 0 0 + 2 1 1 1 0 0 + 3 1 0 0 1 0 + 4 1 0 0 1 0 + 5 0 0 0 1 0 + 6 0 0 0 1 0 + 7 0 0 0 0 0 + 8 0 0 0 0 0 + 9 0 0 0 0 1 + 10 1 0 0 0 1 From b468d74c4e812a495ec0f1578f24c67cb6d8ae30 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 11:59:54 -0500 Subject: [PATCH 24/34] Add bridge time series sync --- src/Action_HydrogenBond.cpp | 103 +++++++++++++++++++++++-------- src/Action_HydrogenBond.h | 26 ++++---- test/Test_Hbond_Split/RunTest.sh | 1 + 3 files changed, 92 insertions(+), 38 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index f3e7247902..3b348fd8c6 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -292,6 +292,15 @@ static inline std::string CreateHBlegend(Topology const& topIn, int a_atom, int topIn[h_atom].Name().Truncated()); } +/** Create legend for bridge based on given indices. */ +static inline std::string CreateBridgeLegend(std::set indices) +{ + std::string blegend("B"); + for (std::set::const_iterator brs = indices.begin(); brs != indices.end(); ++brs) + blegend.append("_" + integerToString(*brs + 1)); + return blegend; +} + // Action_HydrogenBond::Setup() Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { CurrentParm_ = setup.TopAddress(); @@ -1022,10 +1031,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { bds = (DataSet_integer*) masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"bridge",BridgeMap_.size())); // Create a legend from the indices. - std::string blegend("B"); - for (std::set::const_iterator brs = bridge->second.begin(); brs != bridge->second.end(); ++brs) - blegend.append("_" + integerToString(*brs + 1)); - bds->SetLegend( blegend ); + bds->SetLegend( CreateBridgeLegend( bridge->second ) ); if (Bseriesout_ != 0) Bseriesout_->AddDataSet( bds ); } b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(bridge->second, Bridge(bds, splitFrames_)) ); @@ -1051,6 +1057,19 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { return Action::OK; } +/** Ensure DataSet time series has at least N frames, fill out if necessary. */ +void Action_HydrogenBond::FinishSeries(DataSet_integer* data, unsigned int N) { + static const int ZERO = 0; + if (data != 0 && N > 0) { + if ( data->Size() < N ) { + data->Add( N-1, &ZERO ); +# ifdef MPI + data->SetNeedsSync( false ); +# endif + } + } +} + #ifdef MPI // Action_HydrogenBond::GetRankNhbonds() /** Determine how many hydrogen bonds are on each rank. */ @@ -1214,9 +1233,9 @@ int Action_HydrogenBond::SyncAction() { // updated to reflect overall # frames. if (series_) { for (UUmapType::iterator hb = UU_Map_.begin(); hb != UU_Map_.end(); ++hb) - hb->second.FinishSeries( Nframes_ ); + FinishSeries( hb->second.Data(), Nframes_ ); for (UVmapType::iterator hb = UV_Map_.begin(); hb != UV_Map_.end(); ++hb) - hb->second.FinishSeries( Nframes_ ); + FinishSeries( hb->second.Data(), Nframes_ ); } } else { // NON-MASTER RANK @@ -1262,6 +1281,7 @@ int Action_HydrogenBond::SyncAction() { // MASTER RANK for (int rank = 1; rank < trajComm_.Size(); rank++) { + std::vector Svals; // Receive size of iArray trajComm_.Recv( &iSize, 1, MPI_INT, rank, 1302 ); //mprintf("DEBUG: Receiving %i bridges from rank %i\n", iSize, rank); @@ -1274,17 +1294,50 @@ int Action_HydrogenBond::SyncAction() { for (int ir = 0; ir != iArray[idx]; ir++, i2++) residues.insert( iArray[i2] ); BmapType::iterator b_it = BridgeMap_.lower_bound( residues ); + DataSet_integer* bds = 0; if (b_it == BridgeMap_.end() || b_it->first != residues ) { - // New Bridge - b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(residues, Bridge(iArray[i2])) ); + // Bridge not found on master. Create new Bridge. + if (Bseries_) { + bds = (DataSet_integer*) + masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"bridge",BridgeMap_.size())); + // Create a legend from the indices. + bds->SetLegend( CreateBridgeLegend( residues ) ); + } + b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(residues, Bridge(bds, iArray[i2])) ); b_it->second.SetupParts(nParts, &iArray[0] + i2 + 1); + if (bds != 0) mprintf("DEBUG: '%s' was not on master.\n", bds->legend()); + if (Bseriesout_ != 0) Bseriesout_->AddDataSet( bds ); } else { - // Increment bridge #frames + // Bridge on master and rank. Increment bridge #frames. + bds = b_it->second.Data(); b_it->second.Combine( iArray[i2] ); b_it->second.CombineParts(nParts, &iArray[0] + i2 + 1); + if (bds != 0) mprintf("DEBUG: '%s' was already on master.\n", bds->legend()); } + Svals.push_back( bds ); idx = i2 + 1 + nParts; } + // Update all time series + if (Bseries_) { + for (unsigned int in = 0; in != Svals.size(); in++) { + DataSet_integer* ds = Svals[in]; + //ds->Resize( Nframes_ ); + //int* d_beg = ds->Ptr() + rank_offsets[ rank ]; + rprintf("Receiving %i frames of bridge series data for %s, starting frame %i, # frames %i from rank %i (%i)\n", + Nframes_, ds->legend(), rank_offsets[rank], rank_frames[rank], rank, 1304 + in); + ds->Recv(Nframes_, rank_offsets[ rank ], rank_frames[ rank ], + rank, 1304 + in, trajComm_); + //trajComm_.Recv( d_beg, rank_frames[ rank ], MPI_INT, rank, 1304 + in ); + ds->SetNeedsSync( false ); + } + } + } // END LOOP OVER MASTER RANKS + // At this point we have all bridges from all ranks. Mark all bridge sets + // smaller than Nframes_ as synced and ensure the time series has been + // updated to reflect overall # frames. + if (Bseries_) { + for (BmapType::iterator b = BridgeMap_.begin(); b != BridgeMap_.end(); ++b) + FinishSeries( b->second.Data(), Nframes_ ); } } else { // NON-MASTER @@ -1305,6 +1358,16 @@ int Action_HydrogenBond::SyncAction() { iSize = (int)iArray.size(); trajComm_.Send( &iSize, 1, MPI_INT, 0, 1302 ); trajComm_.Send( &(iArray[0]), iSize, MPI_INT, 0, 1303 ); + // Send series data to master + if (Bseries_) { + int in = 0; // For tag + for (BmapType::const_iterator b = BridgeMap_.begin(); b != BridgeMap_.end(); ++b, in++) { + rprintf("Sending %zu frames of %s to master (%i).\n", b->second.Data()->Size(), b->second.Data()->legend(), 1304+in); + b->second.Data()->Send( 0, 1304 + in, trajComm_ ); + //trajComm_.Send( hb->second.Data()->Ptr(), hb->second.Data()->Size(), MPI_INT, 0, 1304 + in ); + b->second.Data()->SetNeedsSync( false ); + } + } } } // END COMMUNICATING BRIDGE DATA TO MASTER return 0; @@ -1319,9 +1382,13 @@ void Action_HydrogenBond::UpdateSeries() { if (seriesUpdated_) return; if (series_ && Nframes_ > 0) { for (UUmapType::iterator hb = UU_Map_.begin(); hb != UU_Map_.end(); ++hb) - hb->second.FinishSeries(Nframes_); + FinishSeries(hb->second.Data(), Nframes_); for (UVmapType::iterator hb = UV_Map_.begin(); hb != UV_Map_.end(); ++hb) - hb->second.FinishSeries(Nframes_); + FinishSeries(hb->second.Data(), Nframes_); + } + if (Bseries_ && Nframes_ > 0) { + for (BmapType::iterator b = BridgeMap_.begin(); b != BridgeMap_.end(); ++b) + FinishSeries( b->second.Data(), Nframes_ ); } seriesUpdated_ = true; } @@ -1562,17 +1629,3 @@ void Action_HydrogenBond::Hbond::CalcAvg() { angle_ *= Constants::RADDEG; } -/** For filling in series data. */ -const int Action_HydrogenBond::Hbond::ZERO = 0; - -/** Calculate average angle/distance. */ -void Action_HydrogenBond::Hbond::FinishSeries(unsigned int N) { - if (data_ != 0 && N > 0) { - if ( data_->Size() < N ) { - data_->Add( N-1, &ZERO ); -# ifdef MPI - data_->SetNeedsSync( false ); -# endif - } - } -} diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index e563799749..575e3cae32 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -60,6 +60,8 @@ class Action_HydrogenBond : public Action { /// Flatten given Hbond into arrays static void HbondToArray(std::vector&, std::vector&, Hbond const&); # endif + /// Ensure DataSet time series has at least N frames, fill out if necessary. + static inline void FinishSeries(DataSet_integer*, unsigned int); typedef std::vector Sarray; typedef std::pair Hpair; @@ -168,8 +170,6 @@ class Action_HydrogenBond::Hbond { partsAng_.resize(splits.size()+1); } } -// Hbond(double d, double a, DataSet_integer* s, int ia, int ih, int id) : -// dist_(d), angle_(a), data_(s), A_(ia), H_(ih), D_(id), frames_(1) {} # ifdef _OPENMP /// Just record that hbond exists Hbond(double d, double a, int ia, int ih, int id) : @@ -178,12 +178,13 @@ class Action_HydrogenBond::Hbond { Hbond(double d, double a, int ia, int ih, int id, int sd) : dist_(d), angle_(a), data_(0), A_(ia), H_(ih), D_(id), frames_(sd) {} # endif - double Dist() const { return dist_; } - double Angle() const { return angle_; } - int Frames() const { return frames_; } - int A() const { return A_; } - int H() const { return H_; } - int D() const { return D_; } + double Dist() const { return dist_; } + double Angle() const { return angle_; } + int Frames() const { return frames_; } + int A() const { return A_; } + int H() const { return H_; } + int D() const { return D_; } + DataSet_integer* Data() const { return data_; } // Summary by parts unsigned int Nparts() const { return partsDist_.size(); } unsigned int PartFrames(unsigned int idx) const { return (unsigned int)partsDist_[idx].nData(); } @@ -191,7 +192,6 @@ class Action_HydrogenBond::Hbond { Stats const& PartDist(unsigned int idx) const { return partsDist_[idx]; } Stats const& PartAngle(unsigned int idx) const { return partsAng_[idx]; } # ifdef MPI - DataSet_integer* Data() const { return data_; } /// CONSTRUCTOR - New hydrogen bond with given # frames Hbond(double d, double a, DataSet_integer* s, int ia, int ih, int id, int n) : dist_(d), angle_(a), data_(s), A_(ia), H_(ih), D_(id), frames_(n) {} @@ -246,9 +246,7 @@ class Action_HydrogenBond::Hbond { } } void CalcAvg(); - void FinishSeries(unsigned int); private: - static const int ZERO; double dist_; ///< Used to calculate average distance of hydrogen bond double angle_; ///< Used to calculate average angle of hydrogen bond DataSet_integer* data_; ///< Hold time series data @@ -270,7 +268,7 @@ class Action_HydrogenBond::Bridge { } # ifdef MPI /// Constructor - new bridge with given # frames - Bridge(int f) : data_(0), frames_(f) {} + Bridge(DataSet_integer* bds, int f) : data_(bds), frames_(f) {} /// Increment number of frames void Combine(int n) { frames_ += n; } /// Set up parts with given int array containing # frames for each part @@ -287,7 +285,9 @@ class Action_HydrogenBond::Bridge { partsFrames_[idx] += ivals[idx]; } # endif - int Frames() const { return frames_; } + /// \return internal data set + DataSet_integer* Data() const { return data_; } + int Frames() const { return frames_; } /// Update frames/time series void Update(int fnum, Iarray const& splitFrames, int onum) { ++frames_; diff --git a/test/Test_Hbond_Split/RunTest.sh b/test/Test_Hbond_Split/RunTest.sh index 196e2ce779..2025130e62 100755 --- a/test/Test_Hbond_Split/RunTest.sh +++ b/test/Test_Hbond_Split/RunTest.sh @@ -35,6 +35,7 @@ trajin ../tz2.ortho.nc hbond hb out solvhb.dat :1-13 solventacceptor :WAT@O solventdonor :WAT \ solvout solvavg.dat bridgeout solvavg.dat splitframe 5 \ bseries bseriesfile bridges.dat +run EOF RunCpptraj "$UNITNAME" DoTest solvavg.dat.save solvavg.dat From c538fb04a5384125d1acda4689a475d2cff4ac80 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 13:04:37 -0500 Subject: [PATCH 25/34] Comment out some debug info --- src/Action_HydrogenBond.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 3b348fd8c6..f134487ae5 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1305,14 +1305,14 @@ int Action_HydrogenBond::SyncAction() { } b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(residues, Bridge(bds, iArray[i2])) ); b_it->second.SetupParts(nParts, &iArray[0] + i2 + 1); - if (bds != 0) mprintf("DEBUG: '%s' was not on master.\n", bds->legend()); + //if (bds != 0) mprintf("DEBUG: '%s' was not on master.\n", bds->legend()); if (Bseriesout_ != 0) Bseriesout_->AddDataSet( bds ); } else { // Bridge on master and rank. Increment bridge #frames. bds = b_it->second.Data(); b_it->second.Combine( iArray[i2] ); b_it->second.CombineParts(nParts, &iArray[0] + i2 + 1); - if (bds != 0) mprintf("DEBUG: '%s' was already on master.\n", bds->legend()); + //if (bds != 0) mprintf("DEBUG: '%s' was already on master.\n", bds->legend()); } Svals.push_back( bds ); idx = i2 + 1 + nParts; @@ -1323,8 +1323,8 @@ int Action_HydrogenBond::SyncAction() { DataSet_integer* ds = Svals[in]; //ds->Resize( Nframes_ ); //int* d_beg = ds->Ptr() + rank_offsets[ rank ]; - rprintf("Receiving %i frames of bridge series data for %s, starting frame %i, # frames %i from rank %i (%i)\n", - Nframes_, ds->legend(), rank_offsets[rank], rank_frames[rank], rank, 1304 + in); + //rprintf("Receiving %i frames of bridge series data for %s, starting frame %i, # frames %i from rank %i (%i)\n", + // Nframes_, ds->legend(), rank_offsets[rank], rank_frames[rank], rank, 1304 + in); ds->Recv(Nframes_, rank_offsets[ rank ], rank_frames[ rank ], rank, 1304 + in, trajComm_); //trajComm_.Recv( d_beg, rank_frames[ rank ], MPI_INT, rank, 1304 + in ); @@ -1362,7 +1362,7 @@ int Action_HydrogenBond::SyncAction() { if (Bseries_) { int in = 0; // For tag for (BmapType::const_iterator b = BridgeMap_.begin(); b != BridgeMap_.end(); ++b, in++) { - rprintf("Sending %zu frames of %s to master (%i).\n", b->second.Data()->Size(), b->second.Data()->legend(), 1304+in); + //rprintf("Sending %zu frames of %s to master (%i).\n", b->second.Data()->Size(), b->second.Data()->legend(), 1304+in); b->second.Data()->Send( 0, 1304 + in, trajComm_ ); //trajComm_.Send( hb->second.Data()->Ptr(), hb->second.Data()->Size(), MPI_INT, 0, 1304 + in ); b->second.Data()->SetNeedsSync( false ); From e9aac3b3a2e97b42bdcf8e343b3447c551b60132 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 13:13:15 -0500 Subject: [PATCH 26/34] Make bridge aspect bridge_ to make it easier to select --- src/Action_HydrogenBond.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index f134487ae5..97f31f3fc6 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -293,9 +293,9 @@ static inline std::string CreateHBlegend(Topology const& topIn, int a_atom, int } /** Create legend for bridge based on given indices. */ -static inline std::string CreateBridgeLegend(std::set indices) +static inline std::string CreateBridgeLegend(std::string const& prefix, std::set indices) { - std::string blegend("B"); + std::string blegend(prefix); for (std::set::const_iterator brs = indices.begin(); brs != indices.end(); ++brs) blegend.append("_" + integerToString(*brs + 1)); return blegend; @@ -1029,9 +1029,9 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { DataSet_integer* bds = 0; if (Bseries_) { bds = (DataSet_integer*) - masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"bridge",BridgeMap_.size())); + masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,CreateBridgeLegend("bridge",bridge->second),BridgeMap_.size())); // Create a legend from the indices. - bds->SetLegend( CreateBridgeLegend( bridge->second ) ); + bds->SetLegend( CreateBridgeLegend( "B", bridge->second ) ); if (Bseriesout_ != 0) Bseriesout_->AddDataSet( bds ); } b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(bridge->second, Bridge(bds, splitFrames_)) ); @@ -1299,9 +1299,9 @@ int Action_HydrogenBond::SyncAction() { // Bridge not found on master. Create new Bridge. if (Bseries_) { bds = (DataSet_integer*) - masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"bridge",BridgeMap_.size())); + masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,CreateBridgeLegend("bridge",residues),BridgeMap_.size())); // Create a legend from the indices. - bds->SetLegend( CreateBridgeLegend( residues ) ); + bds->SetLegend( CreateBridgeLegend( "B", residues ) ); } b_it = BridgeMap_.insert( b_it, std::pair,Bridge>(residues, Bridge(bds, iArray[i2])) ); b_it->second.SetupParts(nParts, &iArray[0] + i2 + 1); From f745dd24823876db5790318fd9e360aac6546219 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 13:27:49 -0500 Subject: [PATCH 27/34] Fix serial compile --- src/Action_HydrogenBond.cpp | 37 +++++++++++++++++++++++++++++++++---- src/Action_HydrogenBond.h | 4 ++++ 2 files changed, 37 insertions(+), 4 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 97f31f3fc6..bb4314dbed 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -635,7 +635,13 @@ double Action_HydrogenBond::Angle(const double* XA, const double* XH, const doub } } +/** This version is for when fnum == onum (non-MPI) */ +void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,int h_atom,int d_atom,bool udonor) { + AddUV(dist, angle, fnum, a_atom, h_atom, d_atom, udonor, fnum); +} + // Action_HydrogenBond::AddUV() +/** Add or update a solute-solvent hydrogen bond with given angle/distance. */ void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,int h_atom,int d_atom,bool udonor, int onum) { int hbidx, solventres, soluteres; @@ -712,7 +718,11 @@ void Action_HydrogenBond::CalcSolvHbonds(int frameNum, double dist2, thread_HBs_[numHB].push_back( Hbond(sqrt(dist2), angle, a_atom, *h_atom, d_atom, (int)soluteDonor) ); # else ++numHB; - AddUV(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom, soluteDonor, frmIn.TrajoutNum()); + AddUV(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom, soluteDonor +# ifdef MPI + ,frmIn.TrajoutNum() +# endif + ); # endif } } @@ -740,6 +750,13 @@ DataSet_integer* Action_HydrogenBond::UUset(int a_atom, int h_atom, int d_atom) return ds; } +/** Add or update a solute-solute hydrogen bond with given angle/distance, + * input frame == output frame (non-MPI). + */ +void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, int h_atom, int d_atom) { + AddUU(dist, angle, fnum, a_atom, h_atom, d_atom, fnum); +} + // Action_HydrogenBond::AddUU() /** Add or update a solute-solute hydrogen bond with given angle/distance. */ void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, int h_atom, int d_atom, int onum) @@ -786,7 +803,11 @@ void Action_HydrogenBond::CalcSiteHbonds(int frameNum, double dist2, thread_HBs_[numHB].push_back( Hbond(sqrt(dist2), angle, a_atom, *h_atom, d_atom) ); # else ++numHB; - AddUU(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom, frmIn.TrajoutNum()); + AddUU(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom +# ifdef MPI + ,frmIn.TrajoutNum() +# endif + ); # endif } } @@ -900,7 +921,11 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { for (std::vector::iterator it = thread_HBs_.begin(); it != thread_HBs_.end(); ++it) { numHB += (int)it->size(); for (Harray::const_iterator hb = it->begin(); hb != it->end(); ++hb) - AddUU(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D(), frm.TrajoutNum()); + AddUU(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D() +# ifdef MPI + ,frm.TrajoutNum() +# endif + ); it->clear(); } # endif @@ -967,7 +992,11 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { for (std::vector::iterator it = thread_HBs_.begin(); it != thread_HBs_.end(); ++it) { numHB += (int)it->size(); for (Harray::const_iterator hb = it->begin(); hb != it->end(); ++hb) - AddUV(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D(), (bool)hb->Frames(), frm.TrajoutNum()); + AddUV(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D(), (bool)hb->Frames() +# ifdef MPI + ,frm.TrajoutNum() +# endif + ); it->clear(); } # endif diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 575e3cae32..bb86312ad3 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -36,8 +36,12 @@ class Action_HydrogenBond : public Action { inline int UU_Set_Idx(int,int) const; /// \return solute-solute hydrogen bond time series DataSet. inline DataSet_integer* UUset(int,int,int); + /// Add/update solute-solute hydrogen bond, input frame == output frame (non-MPI) + void AddUU(double,double,int,int,int,int); /// Add/update solute-solute hydrogen bond void AddUU(double,double,int,int,int,int,int); + /// Add/update solute-solvent hydrogen bond, input frame == output frame (non-MPI) + void AddUV(double,double,int,int,int,int,bool); /// Add/update solute-solvent hydrogen bond void AddUV(double,double,int,int,int,int,bool,int); /// Calculate UU hydrogen bonds between site and acceptor atom. From e28af8a71d320ae9ff96491cb27e4f1d1e04df95 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 13:40:17 -0500 Subject: [PATCH 28/34] Actually fix the serial/mpi compile --- src/Action_HydrogenBond.cpp | 63 +++++++++++-------------------------- src/Action_HydrogenBond.h | 8 ++--- 2 files changed, 20 insertions(+), 51 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index bb4314dbed..ab30002255 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -635,13 +635,7 @@ double Action_HydrogenBond::Angle(const double* XA, const double* XH, const doub } } -/** This version is for when fnum == onum (non-MPI) */ -void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,int h_atom,int d_atom,bool udonor) { - AddUV(dist, angle, fnum, a_atom, h_atom, d_atom, udonor, fnum); -} - // Action_HydrogenBond::AddUV() -/** Add or update a solute-solvent hydrogen bond with given angle/distance. */ void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,int h_atom,int d_atom,bool udonor, int onum) { int hbidx, solventres, soluteres; @@ -697,7 +691,8 @@ void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,i void Action_HydrogenBond::CalcSolvHbonds(int frameNum, double dist2, Site const& SiteD, const double* XYZD, int a_atom, const double* XYZA, - Frame const& frmIn, int& numHB, bool soluteDonor) + Frame const& frmIn, int& numHB, bool soluteDonor, + int trajoutNum) { // The two sites are close enough to hydrogen bond. int d_atom = SiteD.Idx(); @@ -718,11 +713,7 @@ void Action_HydrogenBond::CalcSolvHbonds(int frameNum, double dist2, thread_HBs_[numHB].push_back( Hbond(sqrt(dist2), angle, a_atom, *h_atom, d_atom, (int)soluteDonor) ); # else ++numHB; - AddUV(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom, soluteDonor -# ifdef MPI - ,frmIn.TrajoutNum() -# endif - ); + AddUV(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom, soluteDonor, trajoutNum); # endif } } @@ -750,13 +741,6 @@ DataSet_integer* Action_HydrogenBond::UUset(int a_atom, int h_atom, int d_atom) return ds; } -/** Add or update a solute-solute hydrogen bond with given angle/distance, - * input frame == output frame (non-MPI). - */ -void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, int h_atom, int d_atom) { - AddUU(dist, angle, fnum, a_atom, h_atom, d_atom, fnum); -} - // Action_HydrogenBond::AddUU() /** Add or update a solute-solute hydrogen bond with given angle/distance. */ void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, int h_atom, int d_atom, int onum) @@ -788,7 +772,8 @@ void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, void Action_HydrogenBond::CalcSiteHbonds(int frameNum, double dist2, Site const& SiteD, const double* XYZD, int a_atom, const double* XYZA, - Frame const& frmIn, int& numHB) + Frame const& frmIn, int& numHB, + int trajoutNum) { // The two sites are close enough to hydrogen bond. int d_atom = SiteD.Idx(); @@ -803,11 +788,7 @@ void Action_HydrogenBond::CalcSiteHbonds(int frameNum, double dist2, thread_HBs_[numHB].push_back( Hbond(sqrt(dist2), angle, a_atom, *h_atom, d_atom) ); # else ++numHB; - AddUU(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom -# ifdef MPI - ,frmIn.TrajoutNum() -# endif - ); + AddUU(sqrt(dist2), angle, frameNum, a_atom, *h_atom, d_atom, trajoutNum); # endif } } @@ -854,9 +835,9 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { if ( !(dist2 > dcut2_) ) { // Site 0 donor, Site 1 acceptor - CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, Site1.Idx(), XYZ1, frm.Frm(), numHB); + CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, Site1.Idx(), XYZ1, frm.Frm(), numHB, frm.TrajoutNum()); // Site 1 donor, Site 0 acceptor - CalcSiteHbonds(frameNum, dist2, Site1, XYZ1, Site0.Idx(), XYZ0, frm.Frm(), numHB); + CalcSiteHbonds(frameNum, dist2, Site1, XYZ1, Site0.Idx(), XYZ0, frm.Frm(), numHB, frm.TrajoutNum()); } } } @@ -867,7 +848,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { const double* XYZ1 = frm.Frm().XYZ( *a_atom ); double dist2 = DIST2( imageOpt_.ImagingType(), XYZ0, XYZ1, frm.Frm().BoxCrd() ); if ( !(dist2 > dcut2_) ) - CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, *a_atom, XYZ1, frm.Frm(), numHB); + CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, *a_atom, XYZ1, frm.Frm(), numHB, frm.TrajoutNum()); } } } @@ -896,9 +877,9 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { if ( !(dist2 > dcut2_) ) { // Site 0 donor, Site 1 acceptor - CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, Site1.Idx(), XYZ1, frm.Frm(), numHB); + CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, Site1.Idx(), XYZ1, frm.Frm(), numHB, frm.TrajoutNum()); // Site 1 donor, Site 0 acceptor - CalcSiteHbonds(frameNum, dist2, Site1, XYZ1, Site0.Idx(), XYZ0, frm.Frm(), numHB); + CalcSiteHbonds(frameNum, dist2, Site1, XYZ1, Site0.Idx(), XYZ0, frm.Frm(), numHB, frm.TrajoutNum()); } } // Loop over solute acceptor-only @@ -907,7 +888,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { const double* XYZ1 = frm.Frm().XYZ( *a_atom ); double dist2 = DIST2( imageOpt_.ImagingType(), XYZ0, XYZ1, frm.Frm().BoxCrd() ); if ( !(dist2 > dcut2_) ) - CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, *a_atom, XYZ1, frm.Frm(), numHB); + CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, *a_atom, XYZ1, frm.Frm(), numHB, frm.TrajoutNum()); } } # ifdef _OPENMP @@ -921,11 +902,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { for (std::vector::iterator it = thread_HBs_.begin(); it != thread_HBs_.end(); ++it) { numHB += (int)it->size(); for (Harray::const_iterator hb = it->begin(); hb != it->end(); ++hb) - AddUU(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D() -# ifdef MPI - ,frm.TrajoutNum() -# endif - ); + AddUU(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D(), frm.TrajoutNum()); it->clear(); } # endif @@ -961,9 +938,9 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { if ( !(dist2 > dcut2_) ) { // Solvent site donor, solute site acceptor - CalcSolvHbonds(frameNum, dist2, Vsite, VXYZ, Both_[sidx].Idx(), UXYZ, frm.Frm(), numHB, false); + CalcSolvHbonds(frameNum, dist2, Vsite, VXYZ, Both_[sidx].Idx(), UXYZ, frm.Frm(), numHB, false, frm.TrajoutNum()); // Solvent site acceptor, solute site donor - CalcSolvHbonds(frameNum, dist2, Both_[sidx], UXYZ, Vsite.Idx(), VXYZ, frm.Frm(), numHB, true); + CalcSolvHbonds(frameNum, dist2, Both_[sidx], UXYZ, Vsite.Idx(), VXYZ, frm.Frm(), numHB, true, frm.TrajoutNum()); } } // Loop over solute sites that are donor only @@ -973,7 +950,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { double dist2 = DIST2( imageOpt_.ImagingType(), VXYZ, UXYZ, frm.Frm().BoxCrd() ); if ( !(dist2 > dcut2_) ) // Solvent site acceptor, solute site donor - CalcSolvHbonds(frameNum, dist2, Both_[sidx], UXYZ, Vsite.Idx(), VXYZ, frm.Frm(), numHB, true); + CalcSolvHbonds(frameNum, dist2, Both_[sidx], UXYZ, Vsite.Idx(), VXYZ, frm.Frm(), numHB, true, frm.TrajoutNum()); } // Loop over solute sites that are acceptor only for (Iarray::const_iterator a_atom = Acceptor_.begin(); a_atom != Acceptor_.end(); ++a_atom) @@ -982,7 +959,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { double dist2 = DIST2( imageOpt_.ImagingType(), VXYZ, UXYZ, frm.Frm().BoxCrd() ); if ( !(dist2 > dcut2_) ) // Solvent site donor, solute site acceptor - CalcSolvHbonds(frameNum, dist2, Vsite, VXYZ, *a_atom, UXYZ, frm.Frm(), numHB, false); + CalcSolvHbonds(frameNum, dist2, Vsite, VXYZ, *a_atom, UXYZ, frm.Frm(), numHB, false, frm.TrajoutNum()); } } // END loop over solvent sites # ifdef _OPENMP @@ -992,11 +969,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { for (std::vector::iterator it = thread_HBs_.begin(); it != thread_HBs_.end(); ++it) { numHB += (int)it->size(); for (Harray::const_iterator hb = it->begin(); hb != it->end(); ++hb) - AddUV(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D(), (bool)hb->Frames() -# ifdef MPI - ,frm.TrajoutNum() -# endif - ); + AddUV(hb->Dist(), hb->Angle(), frameNum, hb->A(), hb->H(), hb->D(), (bool)hb->Frames(), frm.TrajoutNum()); it->clear(); } # endif diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index bb86312ad3..298208fea7 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -36,20 +36,16 @@ class Action_HydrogenBond : public Action { inline int UU_Set_Idx(int,int) const; /// \return solute-solute hydrogen bond time series DataSet. inline DataSet_integer* UUset(int,int,int); - /// Add/update solute-solute hydrogen bond, input frame == output frame (non-MPI) - void AddUU(double,double,int,int,int,int); /// Add/update solute-solute hydrogen bond void AddUU(double,double,int,int,int,int,int); - /// Add/update solute-solvent hydrogen bond, input frame == output frame (non-MPI) - void AddUV(double,double,int,int,int,int,bool); /// Add/update solute-solvent hydrogen bond void AddUV(double,double,int,int,int,int,bool,int); /// Calculate UU hydrogen bonds between site and acceptor atom. void CalcSiteHbonds(int,double,Site const&,const double*,int,const double*, - Frame const&, int&); + Frame const&, int&, int); /// Calculate UV hydrogen bonds between solute/solvent site and solute/solvent acceptor atom. void CalcSolvHbonds(int,double,Site const&,const double*,int,const double*, - Frame const&, int&, bool); + Frame const&, int&, bool, int); /// Update all hydrogen bond time series void UpdateSeries(); /// Determine memory usage from # hbonds and time series From 7f27b807ad370e84304cdb77990a9155c5a7357e Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 13:43:33 -0500 Subject: [PATCH 29/34] Add missing bracket --- src/Action_HydrogenBond.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index ab30002255..e275f62673 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -52,7 +52,7 @@ void Action_HydrogenBond::Help() const { "\t[solvout ] [bridgeout ] [bridgebyatom]\n" "\t[series [uuseries ] [uvseries ]]\n" "\t[bseries [bseriesfile ]]\n" - "\t[splitframe \n" + "\t[splitframe ]\n" " Hydrogen bond is defined as A-HD, where A is acceptor heavy atom, H is\n" " hydrogen, D is donor heavy atom. Hydrogen bond is formed when\n" " A to D distance < dcut and A-H-D angle > acut; if acut < 0 it is ignored.\n" From 82e5010a50dd45564a7524e1e8c836c52d990d01 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 13:45:33 -0500 Subject: [PATCH 30/34] Remove old code --- src/Action_HydrogenBond.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index e275f62673..12f9b9f79b 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -667,7 +667,6 @@ void Action_HydrogenBond::AddUV(double dist, double angle,int fnum, int a_atom,i ds = (DataSet_integer*) masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"solventhb",hbidx)); if (UVseriesout_ != 0) UVseriesout_->AddDataSet( ds ); - //ds->AddVal( fnum, 1 ); } Hbond hb; if (udonor) { // Do not care about which solvent acceptor @@ -754,7 +753,6 @@ void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, DataSet_integer* ds = 0; if (series_) { ds = UUset(a_atom, h_atom, d_atom); - //ds->AddVal( fnum, 1 ); } it = UU_Map_.insert(it, std::pair(hbidx,Hbond(ds, a_atom, h_atom, d_atom, splitFrames_))); } else { From 592047e4a409aae31cc502d4eaaa86b6622cd500 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 13:58:25 -0500 Subject: [PATCH 31/34] Enable the hbond split test. --- test/Makefile | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Makefile b/test/Makefile index 08e10fb0d8..ab45a436e1 100644 --- a/test/Makefile +++ b/test/Makefile @@ -31,6 +31,7 @@ test.center: test.hbond: @-cd Test_Hbond && ./RunTest.sh $(OPT) + @-cd Test_Hbond_Split && ./RunTest.sh $(OPT) test.image: @-cd Test_Image && ./RunTest.sh $(OPT) From 788a5a0593b68313a65fd96eed80735e7f4b94c8 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 14:26:12 -0500 Subject: [PATCH 32/34] Update documentation for new hbond keywords. --- doc/cpptraj.lyx | 63 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 43744f52dc..86ccfba770 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -27297,6 +27297,14 @@ hbond \end_layout +\begin_layout LyX-Code + [bseries [bseriesfile ]] +\end_layout + +\begin_layout LyX-Code + [splitframe ] +\end_layout + \begin_deeper \begin_layout Description [] Data set name. @@ -27457,6 +27465,50 @@ hbond \end_layout \end_deeper +\begin_layout Description +[bseries] Save bridge formed (1.0) or not formed (0.0) per frame for any detected + bridge. + Bridges are saved with aspect [bridge_], where is + an underscore ('_') delimited list of bridged atom/residue numbers (depending + on +\series bold +bridgebyatom +\series default +). +\end_layout + +\begin_deeper +\begin_layout Description +[bseriesfile +\begin_inset space ~ +\end_inset + +] File to write bridge time series data to. +\end_layout + +\end_deeper +\begin_layout Description +[splitframe +\begin_inset space ~ +\end_inset + +] If specified, aplit the average hydrogen bond ( +\series bold +avgout +\series default +, +\series bold +solvout +\series default +, +\series bold +bridgeout +\series default +) analysis into sections delimited by the frame numbers. + For example, 'splitframe 250,500,1000' will divide analysis into frames + 1-249, 250-499, 500-999, and 1000 to end. +\end_layout + \begin_layout Standard Data Sets Created: \end_layout @@ -27520,6 +27572,17 @@ series not present. \end_layout +\begin_layout Description +[bridge_] (bseries only) Time series for bridge; 1 for + present, 0 for not present. + The is an underscore ('_') delimited list of bridged atom/residue + numbers (depending on +\series bold +bridgebyatom +\series default +). +\end_layout + \end_deeper \begin_layout Standard From 77f643b1c3cf99b537484161f16d282db1855b9b Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 14:27:14 -0500 Subject: [PATCH 33/34] 6.2.9. Revision bump for hbond bseries and splitframe keywords. --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index 7d44ad8ba7..ccd9e63a83 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.2.8" +#define CPPTRAJ_INTERNAL_VERSION "V6.2.9" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif From 911c67fc07695dc69696f5b9a8e64a77520d17bd Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 16 Feb 2022 14:55:54 -0500 Subject: [PATCH 34/34] Update dependencies --- src/cpptrajdepend | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 5eff776bfb..72f04e8806 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -37,7 +37,7 @@ Action_FixImagedBonds.o : Action_FixImagedBonds.cpp Action.h ActionState.h Actio Action_GIST.o : Action_GIST.cpp Action.h ActionState.h Action_GIST.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_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridDbl.h DataSet_GridFlt.h DataSet_MatrixFlt.h Dimension.h DispatchObject.h DistRoutines.h Ewald.h EwaldOptions.h Ewald_ParticleMesh.h ExclusionArray.h FileIO.h FileName.h FileTypes.h Frame.h GIST_PME.h Grid.h GridBin.h ImageOption.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h PairList.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ProgressBar.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SplineFxnTable.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h cuda_kernels/GistCudaSetup.cuh helpme_standalone.h Action_Grid.o : Action_Grid.cpp Action.h ActionState.h Action_Grid.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_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h Grid.h GridAction.h GridBin.h MaskArray.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 SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Action_GridFreeEnergy.o : Action_GridFreeEnergy.cpp Action.h ActionState.h Action_GridFreeEnergy.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_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h Grid.h GridAction.h GridBin.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 -Action_HydrogenBond.o : Action_HydrogenBond.cpp Action.h ActionState.h Action_HydrogenBond.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_integer.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 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_HydrogenBond.o : Action_HydrogenBond.cpp Action.h ActionState.h Action_HydrogenBond.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_integer.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 TorsionRoutines.h TypeNameHolder.h Unit.h Vec3.h Action_Image.o : Action_Image.cpp Action.h ActionState.h Action_Image.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 ImageRoutines.h ImageTypes.h Image_List.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_InfraredSpectrum.o : Action_InfraredSpectrum.cpp Action.h ActionState.h Action_InfraredSpectrum.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h ComplexArray.h Constants.h CoordinateInfo.h Corr.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Vector.h DataSet_double.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 ProgressBar.h PubFFT.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_Jcoupling.o : Action_Jcoupling.cpp Action.h ActionState.h Action_Jcoupling.h ArgList.h AssociatedData.h Atom.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_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 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 @@ -92,7 +92,7 @@ AnalysisList.o : AnalysisList.cpp ActionState.h Analysis.h AnalysisList.h Analys Analysis_AmdBias.o : Analysis_AmdBias.cpp ActionState.h Analysis.h AnalysisState.h Analysis_AmdBias.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 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 Analysis_AutoCorr.o : Analysis_AutoCorr.cpp ActionState.h Analysis.h AnalysisState.h Analysis_AutoCorr.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_1D.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 Analysis_Average.o : Analysis_Average.cpp ActionState.h Analysis.h AnalysisState.h Analysis_Average.h ArgList.h Array1D.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 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 -Analysis_Clustering.o : Analysis_Clustering.cpp ActionState.h Analysis.h AnalysisState.h Analysis_Clustering.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cluster/Algorithm.h Cluster/BestReps.h Cluster/CentroidArray.h Cluster/Cframes.h Cluster/Control.h Cluster/DrawGraph.h Cluster/List.h Cluster/Metric.h Cluster/MetricArray.h Cluster/Node.h Cluster/Sieve.h Cluster/Silhouette.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 Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h +Analysis_Clustering.o : Analysis_Clustering.cpp ActionState.h Analysis.h AnalysisState.h Analysis_Clustering.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cluster/Algorithm.h Cluster/BestReps.h Cluster/CentroidArray.h Cluster/Cframes.h Cluster/Control.h Cluster/DrawGraph.h Cluster/List.h Cluster/Metric.h Cluster/MetricArray.h Cluster/Node.h Cluster/Sieve.h Cluster/Silhouette.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 OnlineVarT.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 Analysis_ConstantPHStats.o : Analysis_ConstantPHStats.cpp ActionState.h Analysis.h AnalysisState.h Analysis_ConstantPHStats.h ArgList.h AssociatedData.h Atom.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h Cph.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_1D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_Mesh.h DataSet_pH.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 Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h Analysis_Corr.o : Analysis_Corr.cpp ActionState.h Analysis.h AnalysisState.h Analysis_Corr.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_1D.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 Analysis_CrankShaft.o : Analysis_CrankShaft.cpp ActionState.h Analysis.h AnalysisState.h Analysis_CrankShaft.h Analysis_Statistics.h ArgList.h Array1D.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_float.h DataSet_integer.h DataSet_string.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