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 diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 2d8f57f0db..12f9b9f79b 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), @@ -29,6 +30,7 @@ Action_HydrogenBond::Action_HydrogenBond() : Nframes_(0), debug_(0), series_(false), + Bseries_(false), seriesUpdated_(false), useAtomNum_(false), noIntramol_(false), @@ -49,6 +51,8 @@ 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" " A to D distance < dcut and A-H-D angle > acut; if acut < 0 it is ignored.\n" @@ -60,7 +64,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() @@ -79,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; @@ -89,6 +101,34 @@ 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) { + // 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()) { + 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 +234,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()); @@ -213,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(); @@ -228,6 +279,28 @@ 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()); +} + +/** Create legend for bridge based on given indices. */ +static inline std::string CreateBridgeLegend(std::string const& prefix, std::set indices) +{ + std::string blegend(prefix); + 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(); @@ -539,6 +612,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) @@ -562,7 +636,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? @@ -593,28 +667,31 @@ 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 - if (ds != 0) ds->SetLegend( CurrentParm_->TruncResAtomName(h_atom) + "-V" ); - hb = Hbond(dist,angle,ds,-1,h_atom,d_atom); + 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" ); - hb = Hbond(dist,angle,ds,a_atom,-1,-1); + if (ds != 0) ds->SetLegend( CreateHBlegend(*CurrentParm_, 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_, onum); } // 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, - 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(); @@ -635,7 +712,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, trajoutNum); # endif } } @@ -654,19 +731,18 @@ 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; } // Action_HydrogenBond::AddUU() -void Action_HydrogenBond::AddUU(double dist, double angle, int fnum, int a_atom, int h_atom, int d_atom) +/** 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 Hpair hbidx(h_atom, a_atom); @@ -677,20 +753,25 @@ 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 ); } - 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_, onum); } // 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, - 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(); @@ -705,7 +786,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, trajoutNum); # endif } } @@ -752,9 +833,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()); } } } @@ -765,7 +846,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()); } } } @@ -794,9 +875,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 @@ -805,7 +886,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 @@ -819,7 +900,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 @@ -855,9 +936,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 @@ -867,7 +948,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) @@ -876,7 +957,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 @@ -886,7 +967,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 @@ -903,7 +984,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; @@ -928,25 +1012,35 @@ 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); + 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_,CreateBridgeLegend("bridge",bridge->second),BridgeMap_.size())); + // Create a legend from the indices. + 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_)) ); + } + // 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); @@ -963,21 +1057,22 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { return Action::OK; } -#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()); +/** 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. */ std::vector Action_HydrogenBond::GetRankNhbonds( int num_hb, Parallel::Comm const& commIn ) { std::vector nhb_on_rank; @@ -987,6 +1082,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. @@ -1014,12 +1128,22 @@ 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; + 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_ ); 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]; @@ -1027,8 +1151,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 @@ -1037,7 +1161,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 ); @@ -1047,16 +1171,18 @@ 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 ); } // 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) @@ -1076,10 +1202,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 ); @@ -1105,32 +1233,23 @@ 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 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.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 ); trajComm_.Send( &(Ivals[0]), Ivals.size(), MPI_INT, 0, 1301 ); @@ -1151,16 +1270,18 @@ int Action_HydrogenBond::SyncAction() { } } } - } + } // END COMMUNICATING HBOND DATA TO MASTER 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()) { + // 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); @@ -1173,16 +1294,53 @@ 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 ) - // New Bridge - BridgeMap_.insert( b_it, std::pair,Bridge>(residues, Bridge(iArray[i2])) ); - else - // Increment bridge #frames + DataSet_integer* bds = 0; + if (b_it == BridgeMap_.end() || b_it->first != residues ) { + // Bridge not found on master. Create new Bridge. + if (Bseries_) { + bds = (DataSet_integer*) + masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,CreateBridgeLegend("bridge",residues),BridgeMap_.size())); + // Create a legend from the indices. + 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); + //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] ); - idx = i2 + 1; + 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 // Construct bridge info array. for (BmapType::const_iterator b = BridgeMap_.begin(); b != BridgeMap_.end(); ++b) { @@ -1190,32 +1348,53 @@ 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. 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; } -#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) { 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_ ); } - // 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); @@ -1242,6 +1421,50 @@ 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++) + 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); +} + +/** 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() { @@ -1285,8 +1508,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_); @@ -1298,9 +1524,12 @@ 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()) + summary_Parts(avgout_, *hbond); + avgout_->Printf("\n"); } } @@ -1316,8 +1545,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 @@ -1340,12 +1572,16 @@ 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(); } + // BRIDGING INFO if (bridgeout_ != 0 && calcSolvent_) { if (bridgeByAtom_) @@ -1353,26 +1589,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"); } } } @@ -1386,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 5fc105764d..298208fea7 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -4,6 +4,8 @@ #include "Action.h" #include "ImageOption.h" #include "DataSet_integer.h" +#include "OnlineVarT.h" +//#inc lude "CpptrajStdio.h" // DEBUG #ifdef TIMER # include "Timer.h" #endif @@ -27,23 +29,39 @@ class Action_HydrogenBond : public Action { class Site; 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); - void AddUU(double,double,int,int,int,int); - void AddUV(double,double,int,int,int,int,bool); + /// 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&); + 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 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 + /// \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 + /// 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; @@ -77,6 +95,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_; @@ -90,7 +109,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. @@ -100,6 +120,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. @@ -110,18 +131,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 --------------------------------------------------------------- @@ -147,14 +156,20 @@ 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); + } + } # ifdef _OPENMP /// Just record that hbond exists Hbond(double d, double a, int ia, int ih, int id) : @@ -163,15 +178,21 @@ 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_; } -# ifdef MPI + 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_; } - /// New hydrogen bond with given # frames + // 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; } + Stats const& PartDist(unsigned int idx) const { return partsDist_[idx]; } + Stats const& PartAngle(unsigned int idx) const { return partsAng_[idx]; } +# ifdef MPI + /// 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 @@ -180,6 +201,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 { @@ -189,16 +231,22 @@ 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, int onum) { + dist_ += distIn; + angle_ += angIn; ++frames_; - if (data_ != 0) data_->AddVal(f, 1); + if (data_ != 0) data_->AddVal(fnum, 1); + if (!splitFrames.empty()) { + //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() && onum >= splitFrames[part]) part++; + partsDist_[part].accumulate( distIn ); + partsAng_[part].accumulate( angIn ); + } } 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 @@ -206,24 +254,50 @@ 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: - /// Constructor - new bridge - Bridge() : data_(0), frames_(1) {} + /// CONSTRUCTOR - new bridge + Bridge(DataSet_integer* bds, Iarray const& splits) : data_(bds), 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) {} + 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 + 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_; } + /// \return internal data set + DataSet_integer* Data() const { return data_; } + 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 { @@ -234,8 +308,12 @@ 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 - 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 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); }; 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 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 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) diff --git a/test/Test_Hbond_Split/RunTest.sh b/test/Test_Hbond_Split/RunTest.sh new file mode 100755 index 0000000000..2025130e62 --- /dev/null +++ b/test/Test_Hbond_Split/RunTest.sh @@ -0,0 +1,51 @@ +#!/bin/bash + +. ../MasterTest.sh + +# Clean +CleanFiles hbond.in nhb.dat avghb.dat solvhb.dat solvavg.dat bridges.dat + +INPUT="-i hbond.in" + +# Solute-solute, split by parts +TestUUsplit() { + UNITNAME='Solute Hbond test with split' + CheckFor netcdf + if [ $? -eq 0 ] ; then + cat > hbond.in < hbond.in <