From dea5d065c62ba56b58e3a7d4a805083badf11375 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Apr 2017 12:02:28 -0400 Subject: [PATCH 01/56] DRR - Cpptraj: Add hbond timing info --- src/Action_Hbond.cpp | 18 ++++++++++++++++++ src/Action_Hbond.h | 6 ++++++ 2 files changed, 24 insertions(+) diff --git a/src/Action_Hbond.cpp b/src/Action_Hbond.cpp index ab22cf323d..aebacf9a59 100644 --- a/src/Action_Hbond.cpp +++ b/src/Action_Hbond.cpp @@ -576,11 +576,13 @@ int Action_Hbond::AtomsAreHbonded(Frame const& currentFrame, int frameNum, /** Calculate distance between all donors and acceptors. Store Hbond info. */ Action::RetType Action_Hbond::DoAction(int frameNum, ActionFrame& frm) { + t_action_.Start(); int D, H; // accept ... H-D if (Image_.ImagingEnabled()) frm.Frm().BoxCrd().ToRecip(ucell_, recip_); // SOLUTE-SOLUTE HBONDS + t_uu_.Start(); # ifdef HBOND_OPENMP int dAt, hAt, aAt, didx, aidx, hbidx, mol1 = -1, numHB = 0; double dist, dist2, angle; @@ -674,7 +676,9 @@ Action::RetType Action_Hbond::DoAction(int frameNum, ActionFrame& frm) { NumHbonds_->Add(frameNum, &numHB); //mprintf("HBOND: Scanned %i hbonds.\n",hbidx); # endif + t_uu_.Stop(); if (calcSolvent_) { + t_ud_va_.Start(); // Contains info about which residue(s) a Hbonding solvent mol is // Hbonded to. std::map< int, std::set > solvent2solute; @@ -704,6 +708,8 @@ Action::RetType Action_Hbond::DoAction(int frameNum, ActionFrame& frm) { //mprintf("DEBUG: # Solvent Acceptor to Solute Donor Hbonds is %i\n", numHB); solventHbonds += numHB; } + t_ud_va_.Stop(); + t_vd_ua_.Start(); // SOLVENT DONOR-SOLUTE ACCEPTOR // Index by solute acceptor atom if (hasSolventDonor_) { @@ -729,9 +735,11 @@ Action::RetType Action_Hbond::DoAction(int frameNum, ActionFrame& frm) { //mprintf("DEBUG: # Solvent Donor to Solute Acceptor Hbonds is %i\n", numHB); solventHbonds += numHB; } + t_vd_ua_.Stop(); NumSolvent_->Add(frameNum, &solventHbonds); // Determine number of bridging waters. + t_bridge_.Start(); numHB = 0; std::string bridgeID; for (std::map< int, std::set >::const_iterator bridge = solvent2solute.begin(); @@ -759,10 +767,12 @@ Action::RetType Action_Hbond::DoAction(int frameNum, ActionFrame& frm) { bridgeID.assign("None"); NumBridge_->Add(frameNum, &numHB); BridgeID_->Add(frameNum, bridgeID.c_str()); + t_bridge_.Stop(); } ++Nframes_; + t_action_.Stop(); return Action::OK; } @@ -1005,6 +1015,14 @@ void Action_Hbond::Print() { mprintf("\t%zu unique solute-solvent bridging interactions.\n", BridgeMap_.size()); } + t_uu_.WriteTiming( 2,"Solute-Solute :",t_action_.Total()); + if (calcSolvent_) { + t_ud_va_.WriteTiming( 2,"Solute Donor-Solvent Acceptor :",t_action_.Total()); + t_vd_ua_.WriteTiming( 2,"Solvent Donor-Solute Acceptor :",t_action_.Total()); + t_bridge_.WriteTiming(2,"Bridging waters :",t_action_.Total()); + } + t_action_.WriteTiming(1,"Total:"); + // Ensure all series have been updated for all frames. UpdateSeries(); diff --git a/src/Action_Hbond.h b/src/Action_Hbond.h index 37803eb510..679b9a7bf1 100644 --- a/src/Action_Hbond.h +++ b/src/Action_Hbond.h @@ -6,6 +6,7 @@ #include "Action.h" #include "ImagedAction.h" #include "DataSet_integer.h" +#include "Timer.h" /// Action to calculate the Hbonds present in each frame. class Action_Hbond : public Action { public: @@ -117,5 +118,10 @@ class Action_Hbond : public Action { void SyncMap(HBmapType&, std::vector const&, std::vector const&, const char*, Parallel::Comm const&) const; # endif + Timer t_action_; + Timer t_uu_; + Timer t_ud_va_; + Timer t_vd_ua_; + Timer t_bridge_; }; #endif From 4d65a8376e81aa6d99542dda4735357ae66c1fc0 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 11 Apr 2017 16:12:01 -0400 Subject: [PATCH 02/56] DRR - Cpptraj: Test new hydrogen bond action that arranges things by site to avoid doubled calculation and hopefully make parallelization simpler. --- src/Action_HydrogenBond.cpp | 341 ++++++++++++++++++++++++++++++++++++++++++++ src/Action_HydrogenBond.h | 98 +++++++++++++ 2 files changed, 439 insertions(+) create mode 100644 src/Action_HydrogenBond.cpp create mode 100644 src/Action_HydrogenBond.h diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp new file mode 100644 index 0000000000..882de56d1f --- /dev/null +++ b/src/Action_HydrogenBond.cpp @@ -0,0 +1,341 @@ +#include "Action_HydrogenBond.h" +#include "CpptrajStdio.h" +#include "Constants.h" + +// CONSTRUCTOR +Action_HydrogenBond::Action_HydrogenBond() : + CurrentParm_(0), + masterDSL_(0), + NumHbonds_(0), + NumSolvent_(0), + NumBridge_(0), + BridgeID_(0), + UUseriesout_(0), + UVseriesout_(0), + avgout_(0), + solvout_(0), + bridgeout_(0), + dcut2_(0.0), + acut_(0.0), + debug_(0), + series_(0), + useAtomNum_(false), + noIntramol_(false), + hasDonorMask_(false), + hasDonorHmask_(false), + hasAcceptorMask_(false), + hasSolventDonor_(false), + calcSolvent_(false), + hasSolventAcceptor_(false) +{} + +// void Action_HydrogenBond::Help() +void Action_HydrogenBond::Help() const { + mprintf("\t[] [out ] [] [angle ] [dist ]\n" + "\t[donormask [donorhmask ]] [acceptormask ]\n" + "\t[avgout ] [printatomnum] [nointramol] [image]\n" + "\t[solventdonor ] [solventacceptor ]\n" + "\t[solvout ] [bridgeout ]\n" + "\t[series [uuseries ] [uvseries ]]\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" + " Search for hydrogen bonds using atoms in the region specified by mask.\n" + " If just specified donors and acceptors will be automatically searched for.\n" + " If donormask is specified but not acceptormask, acceptors will be\n" + " automatically searched for in .\n" + " If acceptormask is specified but not donormask, donors will be automatically\n" + " 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"); +} + +// Action_HydrogenBond::Init() +Action::RetType Action_HydrogenBond::Init(ArgList& actionArgs, ActionInit& init, int debugIn) +{ +# ifdef MPI + trajComm_ = init.TrajComm(); +# endif + debug_ = debugIn; + // Get keywords + Image_.InitImaging( (actionArgs.hasKey("image")) ); + DataFile* DF = init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs ); + series_ = actionArgs.hasKey("series"); + if (series_) { + UUseriesout_ = init.DFL().AddDataFile(actionArgs.GetStringKey("uuseries"), actionArgs); + UVseriesout_ = init.DFL().AddDataFile(actionArgs.GetStringKey("uvseries"), actionArgs); + init.DSL().SetDataSetsPending(true); + } + std::string avgname = actionArgs.GetStringKey("avgout"); + std::string solvname = actionArgs.GetStringKey("solvout"); + if (solvname.empty()) solvname = avgname; + std::string bridgename = actionArgs.GetStringKey("bridgeout"); + if (bridgename.empty()) bridgename = solvname; + + useAtomNum_ = actionArgs.hasKey("printatomnum"); + acut_ = actionArgs.getKeyDouble("angle",135.0); + noIntramol_ = actionArgs.hasKey("nointramol"); + // Convert angle cutoff to radians + acut_ *= Constants::DEGRAD; + double dcut = actionArgs.getKeyDouble("dist",3.0); + dcut = actionArgs.getKeyDouble("distance", dcut); // for PTRAJ compat. + dcut2_ = dcut * dcut; + // Get donor mask + std::string mask = actionArgs.GetStringKey("donormask"); + if (!mask.empty()) { + DonorMask_.SetMaskString(mask); + hasDonorMask_=true; + // Get donorH mask (if specified) + mask = actionArgs.GetStringKey("donorhmask"); + if (!mask.empty()) { + DonorHmask_.SetMaskString(mask); + hasDonorHmask_=true; + } + } + // Get acceptor mask + mask = actionArgs.GetStringKey("acceptormask"); + if (!mask.empty()) { + AcceptorMask_.SetMaskString(mask); + hasAcceptorMask_=true; + } + // Get solvent donor mask + mask = actionArgs.GetStringKey("solventdonor"); + if (!mask.empty()) { + SolventDonorMask_.SetMaskString(mask); + hasSolventDonor_ = true; + calcSolvent_ = true; + } + // Get solvent acceptor mask + mask = actionArgs.GetStringKey("solventacceptor"); + if (!mask.empty()) { + SolventAcceptorMask_.SetMaskString(mask); + hasSolventAcceptor_ = true; + calcSolvent_ = true; + } + // Get generic mask + Mask_.SetMaskString(actionArgs.GetMaskNext()); + + // Setup datasets + hbsetname_ = actionArgs.GetStringNext(); + if (hbsetname_.empty()) + hbsetname_ = init.DSL().GenerateDefaultName("HB"); + NumHbonds_ = init.DSL().AddSet(DataSet::INTEGER, MetaData(hbsetname_, "UU")); + if (NumHbonds_==0) return Action::ERR; + if (DF != 0) DF->AddDataSet( NumHbonds_ ); + avgout_ = init.DFL().AddCpptrajFile(avgname, "Avg. solute-solute HBonds"); + if (calcSolvent_) { + NumSolvent_ = init.DSL().AddSet(DataSet::INTEGER, MetaData(hbsetname_, "UV")); + if (NumSolvent_ == 0) return Action::ERR; + if (DF != 0) DF->AddDataSet( NumSolvent_ ); + NumBridge_ = init.DSL().AddSet(DataSet::INTEGER, MetaData(hbsetname_, "Bridge")); + if (NumBridge_ == 0) return Action::ERR; + if (DF != 0) DF->AddDataSet( NumBridge_ ); + BridgeID_ = init.DSL().AddSet(DataSet::STRING, MetaData(hbsetname_, "ID")); + if (BridgeID_ == 0) return Action::ERR; + if (DF != 0) DF->AddDataSet( BridgeID_ ); + solvout_ = init.DFL().AddCpptrajFile(solvname,"Avg. solute-solvent HBonds"); + bridgeout_ = init.DFL().AddCpptrajFile(bridgename,"Solvent bridging info"); + } + + mprintf( " HBOND: "); + if (!hasDonorMask_ && !hasAcceptorMask_) + mprintf("Searching for Hbond donors/acceptors in region specified by %s\n", + Mask_.MaskString()); + else if (hasDonorMask_ && !hasAcceptorMask_) + mprintf("Donor mask is %s, acceptors will be searched for in region specified by %s\n", + DonorMask_.MaskString(), Mask_.MaskString()); + else if (hasAcceptorMask_ && !hasDonorMask_) + mprintf("Acceptor mask is %s, donors will be searched for in a region specified by %s\n", + AcceptorMask_.MaskString(), Mask_.MaskString()); + else + mprintf("Donor mask is %s, Acceptor mask is %s\n", + DonorMask_.MaskString(), AcceptorMask_.MaskString()); + if (hasDonorHmask_) + mprintf("\tSeparate donor H mask is %s\n", DonorHmask_.MaskString() ); + if (noIntramol_) + mprintf("\tOnly looking for intermolecular hydrogen bonds.\n"); + if (hasSolventDonor_) + mprintf("\tWill search for hbonds between solute and solvent donors in [%s]\n", + SolventDonorMask_.MaskString()); + if (hasSolventAcceptor_) + mprintf("\tWill search for hbonds between solute and solvent acceptors in [%s]\n", + SolventAcceptorMask_.MaskString()); + mprintf("\tDistance cutoff = %.3f, Angle Cutoff = %.3f\n",dcut,acut_*Constants::RADDEG); + if (DF != 0) + 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 (calcSolvent_ && solvout_ != 0) + mprintf("\tWriting solute-solvent hbond avgs to %s\n", solvout_->Filename().full()); + if (calcSolvent_ && bridgeout_ != 0) + mprintf("\tWriting solvent bridging info to %s\n", bridgeout_->Filename().full()); + if (useAtomNum_) + mprintf("\tAtom numbers will be written to output.\n"); + if (series_) { + mprintf("\tTime series data for each hbond will be saved for analysis.\n"); + if (UUseriesout_ != 0) mprintf("\tWriting solute-solute time series to %s\n", + UUseriesout_->DataFilename().full()); + if (UVseriesout_ != 0) mprintf("\tWriting solute-solvent time series to %s\n", + UVseriesout_->DataFilename().full()); + } + if (Image_.UseImage()) + mprintf("\tImaging enabled.\n"); + masterDSL_ = init.DslPtr(); + + return Action::OK; +} + +inline bool IsFON(Atom const& atm) { + return (atm.Element() == Atom::FLUORINE || + atm.Element() == Atom::OXYGEN || + atm.Element() == Atom::NITROGEN); +} + +// Action_HydrogenBond::Setup() +Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { + CurrentParm_ = setup.TopAddress(); + Image_.SetupImaging( setup.CoordInfo().TrajBox().Type() ); + + // Set up generic mask + if (!hasDonorMask_ || !hasAcceptorMask_) { + if ( setup.Top().SetupIntegerMask( Mask_ ) ) return Action::ERR; + if ( Mask_.None() ) { + mprintf("Warning: Mask has no atoms.\n"); + return Action::SKIP; + } + } + + // ACCEPTOR MASK SETUP + CharMask Amask; + if (hasAcceptorMask_) { + // Acceptor mask specified + if ( setup.Top().SetupIntegerMask( AcceptorMask_ ) ) return Action::ERR; + if (AcceptorMask_.None()) { + mprintf("Warning: AcceptorMask has no atoms.\n"); + return Action::SKIP; + } + Acceptor_ = AcceptorMask_.Selected(); + // Set up the char mask + Amask = CharMask(AcceptorMask_.ConvertToCharMask(), AcceptorMask_.Nselected()); + } else { + // No specified acceptor mask; search generic mask. + // Setup the char mask + Amask = CharMask( setup.Top().Natom() ); + for (AtomMask::const_iterator at = Mask_.begin(); at != Mask_.end(); ++at) { + if (IsFON( setup.Top()[*at] )) { + Amask.Select(*at); + Acceptor_.push_back( *at ); + } + } + } + + // DONOR/ACCEPTOR SETUP + if (hasDonorMask_) { + // Donor heavy atom mask specified + if ( setup.Top().SetupIntegerMask( DonorMask_ ) ) return Action::ERR; + if (DonorMask_.None()) { + mprintf("Warning: DonorMask has no atoms.\n"); + return Action::SKIP; + } + if ( hasDonorHmask_ ) { + // Donor hydrogen mask also specified + if ( setup.Top().SetupIntegerMask( DonorHmask_ ) ) return Action::ERR; + if ( DonorHmask_.None() ) { + mprintf("Warning: Donor H mask has no atoms.\n"); + return Action::SKIP; + } + if ( DonorHmask_.Nselected() != DonorMask_.Nselected() ) { + mprinterr("Error: There is not a 1 to 1 correspondance between donor and donorH masks.\n"); + mprinterr("Error: donor (%i atoms), donorH (%i atoms).\n",DonorMask_.Nselected(), + DonorHmask_.Nselected()); + return Action::ERR; + } + for (int midx = 0; midx != DonorMask_.Nselected(); midx++) + if ( Amask.AtomInCharMask(DonorMask_[midx]) ) + Both_.push_back( Site(DonorMask_[midx], DonorHmask_[midx]) ); + else + Donor_.push_back( Site(DonorMask_[midx], DonorHmask_[midx]) ); + } else { + // No donor hydrogen mask; use any hydrogens bonded to donor heavy atoms. + for (AtomMask::const_iterator at = DonorMask_.begin(); at != DonorMask_.end(); ++at) { + Iarray Hatoms; + for (Atom::bond_iterator batom = setup.Top()[*at].bondbegin(); + batom != setup.Top()[*at].bondend(); ++batom) + if (setup.Top()[*batom].Element() == Atom::HYDROGEN) + Hatoms.push_back( *batom ); + if (!Hatoms.empty()) { + if ( Amask.AtomInCharMask(*at) ) + Both_.push_back( Site(*at, Hatoms) ); + else + Donor_.push_back( Site(*at, Hatoms) ); + } + } + } + } else { + // No specified donor mask; search generic mask. + for (AtomMask::const_iterator at = Mask_.begin(); at != Mask_.end(); ++at) { + if (IsFON( setup.Top()[*at] )) { + Iarray Hatoms; + for (Atom::bond_iterator batom = setup.Top()[*at].bondbegin(); + batom != setup.Top()[*at].bondend(); ++batom) + if (setup.Top()[*batom].Element() == Atom::HYDROGEN) + Hatoms.push_back( *batom ); + if (!Hatoms.empty()) { + if ( Amask.AtomInCharMask(*at) ) + Both_.push_back( Site(*at, Hatoms) ); + else + Donor_.push_back( Site(*at, Hatoms) ); + } + } + } + } + if (calcSolvent_) { + // Set up solvent donor/acceptor masks + if (hasSolventDonor_) { + if (setup.Top().SetupIntegerMask( SolventDonorMask_ )) return Action::ERR; + if (SolventDonorMask_.None()) { + mprintf("Warning: SolventDonorMask has no atoms.\n"); + return Action::SKIP; + } + } + if (hasSolventAcceptor_) { + if (setup.Top().SetupIntegerMask( SolventAcceptorMask_ )) return Action::ERR; + if (SolventAcceptorMask_.None()) { + mprintf("Warning: SolventAcceptorMask has no atoms.\n"); + return Action::SKIP; + } + } + } + + mprintf("Acceptor atoms (%zu):\n", Acceptor_.size()); + for (Iarray::const_iterator at = Acceptor_.begin(); at != Acceptor_.end(); ++at) + mprintf("\t%20s %8i\n", setup.Top().TruncResAtomName(*at).c_str(), *at+1); + mprintf("Donor/acceptor sites (%zu):\n", Both_.size()); + for (Sarray::const_iterator si = Both_.begin(); si != Both_.end(); ++si) { + mprintf("\t%20s %8i", setup.Top().TruncResAtomName(si->Idx()).c_str(), si->Idx()+1); + for (Iarray::const_iterator at = si->Hbegin(); at != si->Hend(); ++at) + mprintf(" %s", setup.Top()[*at].c_str()); + mprintf("\n"); + } + mprintf("Donor sites (%zu):\n", Donor_.size()); + for (Sarray::const_iterator si = Donor_.begin(); si != Donor_.end(); ++si) { + mprintf("\t%20s %8i", setup.Top().TruncResAtomName(si->Idx()).c_str(), si->Idx()+1); + for (Iarray::const_iterator at = si->Hbegin(); at != si->Hend(); ++at) + mprintf(" %s", setup.Top()[*at].c_str()); + mprintf("\n"); + } + + + return Action::OK; +} + +// Action_HydrogenBond::DoAction() +Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { + + return Action::OK; +} + +// Action_HydrogenBond::Print() +void Action_HydrogenBond::Print() { +} diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h new file mode 100644 index 0000000000..af1fbdacc4 --- /dev/null +++ b/src/Action_HydrogenBond.h @@ -0,0 +1,98 @@ +#ifndef INC_ACTION_HYDROGENBOND_H +#define INC_ACTION_HYDROGENBOND_H +#include "Action.h" +#include "ImagedAction.h" +#include "DataSet_integer.h" +#include "Timer.h" +class Action_HydrogenBond : public Action { + public: + Action_HydrogenBond(); + DispatchObject* Alloc() const { return (DispatchObject*)new Action_HydrogenBond(); } + void Help() const; + private: + Action::RetType Init(ArgList&, ActionInit&, int); + Action::RetType Setup(ActionSetup&); + Action::RetType DoAction(int, ActionFrame&); +# ifdef MPI + int SyncAction(); + Parallel::Comm trajComm_; +# endif + void Print(); + + typedef std::vector Iarray; + + class Site; + class Hbond; + + typedef std::vector Sarray; + + ImagedAction Image_; ///< Hold imaging info. + Sarray Donor_; ///< Array of sites that are just donor. + Sarray Both_; ///< Array of sites that are donor and acceptor + Iarray Acceptor_; ///< Array of acceptor-only atom indices + + std::string hbsetname_; + AtomMask DonorMask_; + AtomMask DonorHmask_; + AtomMask AcceptorMask_; + AtomMask SolventDonorMask_; + AtomMask SolventAcceptorMask_; + AtomMask Mask_; + Topology* CurrentParm_; ///< Used to set atom/residue labels + DataSetList* masterDSL_; + DataSet* NumHbonds_; + DataSet* NumSolvent_; + DataSet* NumBridge_; + DataSet* BridgeID_; + DataFile* UUseriesout_; + DataFile* UVseriesout_; + CpptrajFile* avgout_; + CpptrajFile* solvout_; + CpptrajFile* bridgeout_; + double dcut2_; + double acut_; + int debug_; + bool series_; + bool useAtomNum_; + bool noIntramol_; + bool hasDonorMask_; + bool hasDonorHmask_; + bool hasAcceptorMask_; + bool hasSolventDonor_; + bool calcSolvent_; + bool hasSolventAcceptor_; + +}; +// ----- CLASSES --------------------------------------------------------------- +class Action_HydrogenBond::Site { + public: + Site() : idx_(-1), isV_(false) {} + /// Solute site - heavy atom, hydrogen atom + Site(int d, int h) : hlist_(1,h), idx_(d), isV_(false) {} + /// Solute site - heavy atom, list of hydrogen atoms + Site(int d, Iarray const& H) : hlist_(H), idx_(d), isV_(false) {} + /// \return heavy atom index + int Idx() const { return idx_; } + /// \return iterator to beginning of hydrogen indices + Iarray::const_iterator Hbegin() const { return hlist_.begin(); } + /// \return iterator to end of hydrogen indices + Iarray::const_iterator Hend() const { return hlist_.end(); } + private: + Iarray hlist_; ///< List of hydrogen indices + int idx_; ///< Heavy atom index + bool isV_; ///< True if site is solvent +}; + +class Action_HydrogenBond::Hbond { + public: + Hbond() : dist_(0.0), angle_(0.0), data_(0), A_(-1), H_(-1), D_(-1), frames_(0) {} + private: + 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 + int A_; ///< Acceptor atom index + int H_; ///< Hydrogen atom index + int D_; ///< Donor atom index + int frames_; ///< # frames this hydrogen bond has been present +}; +#endif From df323da8d82d5e235360613df78388352c7be2c6 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 12 Apr 2017 08:29:31 -0400 Subject: [PATCH 03/56] DRR - Cpptraj: Add test hbond2 action --- src/Command.cpp | 2 ++ src/cpptrajdepend | 3 ++- src/cpptrajfiles | 1 + 3 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/Command.cpp b/src/Command.cpp index f854143a32..3d22b43c88 100644 --- a/src/Command.cpp +++ b/src/Command.cpp @@ -124,6 +124,7 @@ #include "Action_Volume.h" #include "Action_Align.h" #include "Action_Remap.h" +#include "Action_HydrogenBond.h" // ----- ANALYSIS -------------------------------------------------------------- #include "Analysis_Hist.h" #include "Analysis_Corr.h" @@ -273,6 +274,7 @@ void Command::Init() { Command::AddCmd( new Action_GridFreeEnergy(),Cmd::ACT, 1, "gfe" ); // HIDDEN Command::AddCmd( new Action_Grid(), Cmd::ACT, 1, "grid" ); Command::AddCmd( new Action_Hbond(), Cmd::ACT, 1, "hbond" ); + Command::AddCmd( new Action_HydrogenBond(), Cmd::ACT, 1, "hbond2" ); Command::AddCmd( new Action_Image(), Cmd::ACT, 1, "image" ); Command::AddCmd( new Action_Jcoupling(), Cmd::ACT, 1, "jcoupling" ); Command::AddCmd( new Action_LESsplit(), Cmd::ACT, 1, "lessplit" ); diff --git a/src/cpptrajdepend b/src/cpptrajdepend index a88aaa2d65..76ab31487c 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -34,6 +34,7 @@ Action_GIST.o : Action_GIST.cpp Action.h ActionState.h Action_GIST.h ArgList.h A Action_Grid.o : Action_Grid.cpp Action.h ActionState.h Action_Grid.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 PDBfile.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_GridFreeEnergy.o : Action_GridFreeEnergy.cpp Action.h ActionState.h Action_GridFreeEnergy.h ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h TextFormat.h Timer.h Topology.h Vec3.h Action_Hbond.o : Action_Hbond.cpp Action.h ActionState.h Action_Hbond.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h TextFormat.h Timer.h Topology.h TorsionRoutines.h Vec3.h +Action_HydrogenBond.o : Action_HydrogenBond.cpp Action.h ActionState.h Action_HydrogenBond.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ImagedAction.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_Image.o : Action_Image.cpp Action.h ActionState.h Action_Image.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h Action_Jcoupling.o : Action_Jcoupling.cpp Action.h ActionState.h Action_Jcoupling.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.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 ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h TorsionRoutines.h Vec3.h Action_LESsplit.o : Action_LESsplit.cpp Action.h ActionFrameCounter.h ActionState.h Action_LESsplit.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h TrajectoryFile.h TrajectoryIO.h Trajout_Single.h Vec3.h @@ -326,7 +327,7 @@ TrajoutList.o : TrajoutList.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra. Vec3.o : Vec3.cpp Constants.h CpptrajStdio.h Vec3.h ViewRst.o : ViewRst.cpp ActionFrameCounter.h ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h ParmFile.h ParmIO.h Range.h ReplicaDimArray.h Residue.h Topology.h TrajectoryFile.h TrajectoryIO.h Trajout_Single.h Vec3.h ViewRst.h Action_Esander.o : Action_Esander.cpp Action.h ActionState.h Action_Esander.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h Energy_Sander.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h TextFormat.h Timer.h Topology.h Vec3.h -Command.o : Command.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Action_Align.h Action_Angle.h Action_AreaPerMol.h Action_AtomMap.h Action_AtomicCorr.h Action_AtomicFluct.h Action_AutoImage.h Action_Average.h Action_Bounds.h Action_Box.h Action_Center.h Action_Channel.h Action_CheckChirality.h Action_CheckStructure.h Action_Closest.h Action_ClusterDihedral.h Action_Contacts.h Action_CreateCrd.h Action_CreateReservoir.h Action_DNAionTracker.h Action_DSSP.h Action_Density.h Action_Diffusion.h Action_Dihedral.h Action_Dipole.h Action_DistRmsd.h Action_Distance.h Action_Energy.h Action_Esander.h Action_FilterByData.h Action_FixAtomOrder.h Action_GIST.h Action_Grid.h Action_GridFreeEnergy.h Action_Hbond.h Action_Image.h Action_Jcoupling.h Action_LESsplit.h Action_LIE.h Action_MakeStructure.h Action_Mask.h Action_Matrix.h Action_MinImage.h Action_Molsurf.h Action_MultiDihedral.h Action_MultiVector.h Action_NAstruct.h Action_NMRrst.h Action_NativeContacts.h Action_OrderParameter.h Action_Outtraj.h Action_PairDist.h Action_Pairwise.h Action_Principal.h Action_Projection.h Action_Pucker.h Action_Radgyr.h Action_Radial.h Action_RandomizeIons.h Action_Remap.h Action_ReplicateCell.h Action_Rmsd.h Action_Rotate.h Action_RunningAvg.h Action_STFC_Diffusion.h Action_Scale.h Action_SetVelocity.h Action_Spam.h Action_Strip.h Action_Surf.h Action_SymmetricRmsd.h Action_Temperature.h Action_Translate.h Action_Unstrip.h Action_Unwrap.h Action_Vector.h Action_VelocityAutoCorr.h Action_Volmap.h Action_Volume.h Action_Watershell.h Analysis.h AnalysisList.h AnalysisState.h Analysis_AmdBias.h Analysis_AutoCorr.h Analysis_Average.h Analysis_Clustering.h Analysis_Corr.h Analysis_CrankShaft.h Analysis_CrdFluct.h Analysis_CrossCorr.h Analysis_CurveFit.h Analysis_Divergence.h Analysis_FFT.h Analysis_Hist.h Analysis_IRED.h Analysis_Integrate.h Analysis_KDE.h Analysis_Lifetime.h Analysis_LowestCurve.h Analysis_Matrix.h Analysis_MeltCurve.h Analysis_Modes.h Analysis_MultiHist.h Analysis_Multicurve.h Analysis_Overlap.h Analysis_PhiPsi.h Analysis_Regression.h Analysis_RemLog.h Analysis_Rms2d.h Analysis_RmsAvgCorr.h Analysis_Rotdif.h Analysis_RunningAvg.h Analysis_Spline.h Analysis_State.h Analysis_Statistics.h Analysis_TI.h Analysis_Timecorr.h Analysis_VectorMath.h Analysis_Wavelet.h ArgList.h Array1D.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AxisType.h BaseIOtype.h Box.h BufferedLine.h CharMask.h ClusterDist.h ClusterList.h ClusterMap.h ClusterNode.h ClusterSieve.h Cmd.h CmdInput.h CmdList.h Command.h ComplexArray.h CoordinateInfo.h Corr.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_RemLog.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_string.h Deprecated.h DihedralSearch.h Dimension.h DispatchObject.h DistRoutines.h Energy.h Energy_Sander.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Ewald.h Exec.h Exec_Analyze.h Exec_Calc.h Exec_ClusterMap.h Exec_CombineCoords.h Exec_Commands.h Exec_CompareTop.h Exec_CrdAction.h Exec_CrdOut.h Exec_DataFile.h Exec_DataFilter.h Exec_DataSetCmd.h Exec_GenerateAmberRst.h Exec_Help.h Exec_LoadCrd.h Exec_LoadTraj.h Exec_ParmBox.h Exec_ParmSolvent.h Exec_ParmStrip.h Exec_ParmWrite.h Exec_PermuteDihedrals.h Exec_Precision.h Exec_PrintData.h Exec_ReadData.h Exec_ReadInput.h Exec_RotateDihedral.h Exec_RunAnalysis.h Exec_ScaleDihedralK.h Exec_SequenceAlign.h Exec_System.h Exec_Top.h Exec_Traj.h Exec_ViewRst.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h Grid.h GridAction.h GridBin.h HistBin.h Hungarian.h ImageTypes.h ImagedAction.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h NetcdfFile.h OnlineVarT.h OutputTrajCommon.h PDBfile.h PairList.h Parallel.h ParameterTypes.h PubFFT.h RPNcalc.h Random.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Spline.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Traj_AmberNetcdf.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h molsurf.h +Command.o : Command.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Action_Align.h Action_Angle.h Action_AreaPerMol.h Action_AtomMap.h Action_AtomicCorr.h Action_AtomicFluct.h Action_AutoImage.h Action_Average.h Action_Bounds.h Action_Box.h Action_Center.h Action_Channel.h Action_CheckChirality.h Action_CheckStructure.h Action_Closest.h Action_ClusterDihedral.h Action_Contacts.h Action_CreateCrd.h Action_CreateReservoir.h Action_DNAionTracker.h Action_DSSP.h Action_Density.h Action_Diffusion.h Action_Dihedral.h Action_Dipole.h Action_DistRmsd.h Action_Distance.h Action_Energy.h Action_Esander.h Action_FilterByData.h Action_FixAtomOrder.h Action_GIST.h Action_Grid.h Action_GridFreeEnergy.h Action_Hbond.h Action_HydrogenBond.h Action_Image.h Action_Jcoupling.h Action_LESsplit.h Action_LIE.h Action_MakeStructure.h Action_Mask.h Action_Matrix.h Action_MinImage.h Action_Molsurf.h Action_MultiDihedral.h Action_MultiVector.h Action_NAstruct.h Action_NMRrst.h Action_NativeContacts.h Action_OrderParameter.h Action_Outtraj.h Action_PairDist.h Action_Pairwise.h Action_Principal.h Action_Projection.h Action_Pucker.h Action_Radgyr.h Action_Radial.h Action_RandomizeIons.h Action_Remap.h Action_ReplicateCell.h Action_Rmsd.h Action_Rotate.h Action_RunningAvg.h Action_STFC_Diffusion.h Action_Scale.h Action_SetVelocity.h Action_Spam.h Action_Strip.h Action_Surf.h Action_SymmetricRmsd.h Action_Temperature.h Action_Translate.h Action_Unstrip.h Action_Unwrap.h Action_Vector.h Action_VelocityAutoCorr.h Action_Volmap.h Action_Volume.h Action_Watershell.h Analysis.h AnalysisList.h AnalysisState.h Analysis_AmdBias.h Analysis_AutoCorr.h Analysis_Average.h Analysis_Clustering.h Analysis_Corr.h Analysis_CrankShaft.h Analysis_CrdFluct.h Analysis_CrossCorr.h Analysis_CurveFit.h Analysis_Divergence.h Analysis_FFT.h Analysis_Hist.h Analysis_IRED.h Analysis_Integrate.h Analysis_KDE.h Analysis_Lifetime.h Analysis_LowestCurve.h Analysis_Matrix.h Analysis_MeltCurve.h Analysis_Modes.h Analysis_MultiHist.h Analysis_Multicurve.h Analysis_Overlap.h Analysis_PhiPsi.h Analysis_Regression.h Analysis_RemLog.h Analysis_Rms2d.h Analysis_RmsAvgCorr.h Analysis_Rotdif.h Analysis_RunningAvg.h Analysis_Spline.h Analysis_State.h Analysis_Statistics.h Analysis_TI.h Analysis_Timecorr.h Analysis_VectorMath.h Analysis_Wavelet.h ArgList.h Array1D.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMap.h AtomMask.h AxisType.h BaseIOtype.h Box.h BufferedLine.h CharMask.h ClusterDist.h ClusterList.h ClusterMap.h ClusterNode.h ClusterSieve.h Cmd.h CmdInput.h CmdList.h Command.h ComplexArray.h CoordinateInfo.h Corr.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Cmatrix.h DataSet_Coords.h DataSet_Coords_CRD.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_Mat3x3.h DataSet_MatrixDbl.h DataSet_MatrixFlt.h DataSet_Mesh.h DataSet_Modes.h DataSet_RemLog.h DataSet_Vector.h DataSet_double.h DataSet_float.h DataSet_integer.h DataSet_string.h Deprecated.h DihedralSearch.h Dimension.h DispatchObject.h DistRoutines.h Energy.h Energy_Sander.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h Ewald.h Exec.h Exec_Analyze.h Exec_Calc.h Exec_ClusterMap.h Exec_CombineCoords.h Exec_Commands.h Exec_CompareTop.h Exec_CrdAction.h Exec_CrdOut.h Exec_DataFile.h Exec_DataFilter.h Exec_DataSetCmd.h Exec_GenerateAmberRst.h Exec_Help.h Exec_LoadCrd.h Exec_LoadTraj.h Exec_ParmBox.h Exec_ParmSolvent.h Exec_ParmStrip.h Exec_ParmWrite.h Exec_PermuteDihedrals.h Exec_Precision.h Exec_PrintData.h Exec_ReadData.h Exec_ReadInput.h Exec_RotateDihedral.h Exec_RunAnalysis.h Exec_ScaleDihedralK.h Exec_SequenceAlign.h Exec_System.h Exec_Top.h Exec_Traj.h Exec_ViewRst.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h Grid.h GridAction.h GridBin.h HistBin.h Hungarian.h ImageTypes.h ImagedAction.h InputTrajCommon.h MapAtom.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h NetcdfFile.h OnlineVarT.h OutputTrajCommon.h PDBfile.h PairList.h Parallel.h ParameterTypes.h PubFFT.h RPNcalc.h Random.h Range.h ReferenceAction.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Spline.h SymmetricRmsdCalc.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Traj_AmberNetcdf.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h molsurf.h Cpptraj.o : Cpptraj.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Cmd.h CmdInput.h CmdList.h Command.h CoordinateInfo.h Cpptraj.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h ParmFile.h ParmIO.h Range.h ReadLine.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h StringRoutines.h TextFormat.h Timer.h TopInfo.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h Version.h Energy_Sander.o : Energy_Sander.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h Energy_Sander.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterTypes.h ParmFile.h ParmIO.h Range.h ReplicaDimArray.h Residue.h Topology.h Vec3.h ReadLine.o : ReadLine.cpp Action.h ActionFrameCounter.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h BaseIOtype.h Box.h CharMask.h Cmd.h CmdInput.h CmdList.h Command.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataIO.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOut.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FrameArray.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterTypes.h Range.h ReadLine.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h TrajoutList.h Trajout_Single.h Vec3.h diff --git a/src/cpptrajfiles b/src/cpptrajfiles index f308a545f9..592b424b36 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -36,6 +36,7 @@ COMMON_SOURCES=ActionFrameCounter.cpp \ Action_Grid.cpp \ Action_GridFreeEnergy.cpp \ Action_Hbond.cpp \ + Action_HydrogenBond.cpp \ Action_Image.cpp \ Action_Jcoupling.cpp \ Action_LESsplit.cpp \ From 35705f9f28847edc556672fc9f268f5ce8ee5626 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 12 Apr 2017 10:30:23 -0400 Subject: [PATCH 04/56] DRR - Cpptraj: Finish initial solute-solute hbond calculation. Seems to be decent speedup so far but requires much more testing. --- src/Action_HydrogenBond.cpp | 222 ++++++++++++++++++++++++++++++++++++-------- src/Action_HydrogenBond.h | 21 ++++- src/AtomMask.h | 2 + 3 files changed, 206 insertions(+), 39 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 882de56d1f..2ca7d17e2a 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1,6 +1,8 @@ +#include // sqrt #include "Action_HydrogenBond.h" #include "CpptrajStdio.h" #include "Constants.h" +#include "TorsionRoutines.h" // CONSTRUCTOR Action_HydrogenBond::Action_HydrogenBond() : @@ -186,6 +188,7 @@ Action::RetType Action_HydrogenBond::Init(ArgList& actionArgs, ActionInit& init, return Action::OK; } +// IsFON() inline bool IsFON(Atom const& atm) { return (atm.Element() == Atom::FLUORINE || atm.Element() == Atom::OXYGEN || @@ -207,7 +210,6 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { } // ACCEPTOR MASK SETUP - CharMask Amask; if (hasAcceptorMask_) { // Acceptor mask specified if ( setup.Top().SetupIntegerMask( AcceptorMask_ ) ) return Action::ERR; @@ -215,19 +217,14 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { mprintf("Warning: AcceptorMask has no atoms.\n"); return Action::SKIP; } - Acceptor_ = AcceptorMask_.Selected(); - // Set up the char mask - Amask = CharMask(AcceptorMask_.ConvertToCharMask(), AcceptorMask_.Nselected()); } else { // No specified acceptor mask; search generic mask. - // Setup the char mask - Amask = CharMask( setup.Top().Natom() ); + AcceptorMask_.ResetMask(); for (AtomMask::const_iterator at = Mask_.begin(); at != Mask_.end(); ++at) { - if (IsFON( setup.Top()[*at] )) { - Amask.Select(*at); - Acceptor_.push_back( *at ); - } + if (IsFON( setup.Top()[*at] )) + AcceptorMask_.AddSelectedAtom( *at ); } + AcceptorMask_.SetNatoms( Mask_.NmaskAtoms() ); } // DONOR/ACCEPTOR SETUP @@ -251,45 +248,97 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { DonorHmask_.Nselected()); return Action::ERR; } - for (int midx = 0; midx != DonorMask_.Nselected(); midx++) - if ( Amask.AtomInCharMask(DonorMask_[midx]) ) - Both_.push_back( Site(DonorMask_[midx], DonorHmask_[midx]) ); - else - Donor_.push_back( Site(DonorMask_[midx], DonorHmask_[midx]) ); + int maxAtom = std::max( AcceptorMask_.back(), DonorMask_.back() ) + 1; + AtomMask::const_iterator a_atom = AcceptorMask_.begin(); + AtomMask::const_iterator d_atom = DonorMask_.begin(); + AtomMask::const_iterator h_atom = DonorHmask_.begin(); + bool isDonor, isAcceptor; + int H_at = -1; + for (int at = 0; at != maxAtom; at++) + { + isDonor = false; + isAcceptor = false; + if ( d_atom != DonorMask_.end() && *d_atom == at ) { + isDonor = true; + ++d_atom; + H_at = *(h_atom++); + } + if ( a_atom != AcceptorMask_.end() && *a_atom == at ) { + isAcceptor = true; + ++a_atom; + } + if (isDonor && isAcceptor) + Both_.push_back( Site(at, H_at) ); + else if (isDonor) + Donor_.push_back( Site(at, H_at) ); + else if (isAcceptor) + Acceptor_.push_back( at ); + } } else { // No donor hydrogen mask; use any hydrogens bonded to donor heavy atoms. - for (AtomMask::const_iterator at = DonorMask_.begin(); at != DonorMask_.end(); ++at) { + int maxAtom = std::max( AcceptorMask_.back(), DonorMask_.back() ) + 1; + AtomMask::const_iterator a_atom = AcceptorMask_.begin(); + AtomMask::const_iterator d_atom = DonorMask_.begin(); + bool isDonor, isAcceptor; + for (int at = 0; at != maxAtom; at++) + { Iarray Hatoms; - for (Atom::bond_iterator batom = setup.Top()[*at].bondbegin(); - batom != setup.Top()[*at].bondend(); ++batom) - if (setup.Top()[*batom].Element() == Atom::HYDROGEN) - Hatoms.push_back( *batom ); - if (!Hatoms.empty()) { - if ( Amask.AtomInCharMask(*at) ) - Both_.push_back( Site(*at, Hatoms) ); - else - Donor_.push_back( Site(*at, Hatoms) ); - } + isDonor = false; + isAcceptor = false; + if ( d_atom != DonorMask_.end() && *d_atom == at ) { + ++d_atom; + for (Atom::bond_iterator H_at = setup.Top()[at].bondbegin(); + H_at != setup.Top()[at].bondend(); ++H_at) + if (setup.Top()[*H_at].Element() == Atom::HYDROGEN) + Hatoms.push_back( *H_at ); + isDonor = !Hatoms.empty(); + } + if ( a_atom != AcceptorMask_.end() && *a_atom == at ) { + isAcceptor = true; + ++a_atom; + } + if (isDonor && isAcceptor) + Both_.push_back( Site(at, Hatoms) ); + else if (isDonor) + Donor_.push_back( Site(at, Hatoms) ); + else if (isAcceptor) + Acceptor_.push_back( at ); } } } else { // No specified donor mask; search generic mask. - for (AtomMask::const_iterator at = Mask_.begin(); at != Mask_.end(); ++at) { - if (IsFON( setup.Top()[*at] )) { - Iarray Hatoms; - for (Atom::bond_iterator batom = setup.Top()[*at].bondbegin(); - batom != setup.Top()[*at].bondend(); ++batom) - if (setup.Top()[*batom].Element() == Atom::HYDROGEN) - Hatoms.push_back( *batom ); - if (!Hatoms.empty()) { - if ( Amask.AtomInCharMask(*at) ) - Both_.push_back( Site(*at, Hatoms) ); - else - Donor_.push_back( Site(*at, Hatoms) ); + int maxAtom = std::max( AcceptorMask_.back(), Mask_.back() ) + 1; + AtomMask::const_iterator a_atom = AcceptorMask_.begin(); + AtomMask::const_iterator d_atom = Mask_.begin(); + bool isDonor, isAcceptor; + for (int at = 0; at != maxAtom; at++) + { + Iarray Hatoms; + isDonor = false; + isAcceptor = false; + if ( d_atom != Mask_.end() && *d_atom == at) { + ++d_atom; + if ( IsFON( setup.Top()[at] ) ) { + for (Atom::bond_iterator H_at = setup.Top()[at].bondbegin(); + H_at != setup.Top()[at].bondend(); ++H_at) + if (setup.Top()[*H_at].Element() == Atom::HYDROGEN) + Hatoms.push_back( *H_at ); + isDonor = !Hatoms.empty(); } } + if ( a_atom != AcceptorMask_.end() && *a_atom == at ) { + isAcceptor = true; + ++a_atom; + } + if (isDonor && isAcceptor) + Both_.push_back( Site(at, Hatoms) ); + else if (isDonor) + Donor_.push_back( Site(at, Hatoms) ); + else if (isAcceptor) + Acceptor_.push_back( at ); } } + if (calcSolvent_) { // Set up solvent donor/acceptor masks if (hasSolventDonor_) { @@ -330,12 +379,109 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { return Action::OK; } +// Action_HydrogenBond::Angle() +double Action_HydrogenBond::Angle(const double* XA, const double* XH, const double* XD) const +{ + if (Image_.ImageType() == NOIMAGE) + return (CalcAngle(XA, XH, XD)); + else { + double angle; + Vec3 VH = Vec3(XH); + Vec3 H_A = MinImagedVec(VH, Vec3(XA), ucell_, recip_); + Vec3 H_D = Vec3(XD) - VH; + double rha = H_A.Magnitude2(); + double rhd = H_D.Magnitude2(); + if (rha > Constants::SMALL && rhd > Constants::SMALL) { + angle = (H_A * H_D) / sqrt(rha * rhd); + if (angle > 1.0) angle = 1.0; + else if (angle < -1.0) angle = -1.0; + angle = acos(angle); + } else + angle = 0.0; + return angle; + } +} + +// Action_HydrogenBond::CalcSiteHbonds() +void Action_HydrogenBond::CalcSiteHbonds(int hbidx, int frameNum, double dist2, + Site const& SiteD, const double* XYZD, + int a_atom, const double* XYZA, + Frame const& frmIn, int& numHB) +{ + // The two sites are close enough to hydrogen bond. + int d_atom = SiteD.Idx(); + // Determine if angle cutoff is satisfied + for (Iarray::const_iterator h_atom = SiteD.Hbegin(); h_atom != SiteD.Hend(); ++h_atom) + { + double angle = Angle(XYZA, frmIn.XYZ(*h_atom), XYZD); + if ( !(angle < acut_) ) + { + ++numHB; + double dist = sqrt(dist2); + // Index UU hydrogen bonds by DonorH-Acceptor + HBmapType::iterator it = UU_Map_.find( hbidx ); + if (it == UU_Map_.end()) + { + DataSet_integer* ds = 0; + if (series_) { + std::string hblegend = CurrentParm_->TruncResAtomName(a_atom) + "-" + + CurrentParm_->TruncResAtomName(d_atom) + "-" + + (*CurrentParm_)[*h_atom].Name().Truncated(); + ds = (DataSet_integer*) + masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"solutehb",hbidx)); + if (UUseriesout_ != 0) UUseriesout_->AddDataSet( ds ); + ds->SetLegend( hblegend ); + ds->AddVal( frameNum, 1 ); + } + UU_Map_.insert(it, std::pair(hbidx,Hbond(dist,angle,ds,a_atom,*h_atom,d_atom))); + } else + it->second.Update(dist,angle,frameNum); + } + } +} + // Action_HydrogenBond::DoAction() Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { + t_action_.Start(); + if (Image_.ImagingEnabled()) + frm.Frm().BoxCrd().ToRecip(ucell_, recip_); + + int numHB = 0; + int hbidx = 0; + for (unsigned int sidx0 = 0; sidx0 != Both_.size(); sidx0++) + { + const double* XYZ0 = frm.Frm().XYZ( Both_[sidx0].Idx() ); + // Loop over sites that can be both donor and acceptor + for (unsigned int sidx1 = sidx0 + 1; sidx1 != Both_.size(); sidx1++, hbidx++) + { + const double* XYZ1 = frm.Frm().XYZ( Both_[sidx1].Idx() ); + double dist2 = DIST2( XYZ0, XYZ1, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); + if ( !(dist2 > dcut2_) ) + { + // Site 0 donor, Site 1 acceptor + CalcSiteHbonds(hbidx, frameNum, dist2, Both_[sidx0], XYZ0, Both_[sidx1].Idx(), XYZ1, frm.Frm(), numHB); + // Site 1 donor, Site 0 acceptor + CalcSiteHbonds(hbidx, frameNum, dist2, Both_[sidx1], XYZ1, Both_[sidx0].Idx(), XYZ0, frm.Frm(), numHB); + } + } + // Loop over acceptor-only + for (Iarray::const_iterator a_atom = Acceptor_.begin(); a_atom != Acceptor_.end(); ++a_atom, hbidx++) + { + const double* XYZ1 = frm.Frm().XYZ( *a_atom ); + double dist2 = DIST2( XYZ0, XYZ1, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); + if ( !(dist2 > dcut2_) ) + CalcSiteHbonds(hbidx, frameNum, dist2, Both_[sidx0], XYZ0, *a_atom, XYZ1, frm.Frm(), numHB); + } + } + + NumHbonds_->Add(frameNum, &numHB); + t_action_.Stop(); return Action::OK; } // Action_HydrogenBond::Print() void Action_HydrogenBond::Print() { + mprintf(" HBOND:\n"); + t_action_.WriteTiming(1,"Total:"); } diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index af1fbdacc4..002fdf3b07 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -1,5 +1,6 @@ #ifndef INC_ACTION_HYDROGENBOND_H #define INC_ACTION_HYDROGENBOND_H +#include #include "Action.h" #include "ImagedAction.h" #include "DataSet_integer.h" @@ -24,13 +25,21 @@ class Action_HydrogenBond : public Action { class Site; class Hbond; + inline double Angle(const double*, const double*, const double*) const; + void CalcSiteHbonds(int,int,double,Site const&,const double*,int,const double*, + Frame const&, int&); + typedef std::vector Sarray; + typedef std::map HBmapType; ImagedAction Image_; ///< Hold imaging info. Sarray Donor_; ///< Array of sites that are just donor. Sarray Both_; ///< Array of sites that are donor and acceptor Iarray Acceptor_; ///< Array of acceptor-only atom indices + HBmapType UU_Map_; + HBmapType UV_Map_; + std::string hbsetname_; AtomMask DonorMask_; AtomMask DonorHmask_; @@ -38,6 +47,8 @@ class Action_HydrogenBond : public Action { AtomMask SolventDonorMask_; AtomMask SolventAcceptorMask_; AtomMask Mask_; + Matrix_3x3 ucell_, recip_; + Timer t_action_; Topology* CurrentParm_; ///< Used to set atom/residue labels DataSetList* masterDSL_; DataSet* NumHbonds_; @@ -61,8 +72,8 @@ class Action_HydrogenBond : public Action { bool hasSolventDonor_; bool calcSolvent_; bool hasSolventAcceptor_; - }; + // ----- CLASSES --------------------------------------------------------------- class Action_HydrogenBond::Site { public: @@ -86,6 +97,14 @@ class Action_HydrogenBond::Site { 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) {} + void Update(double d, double a, int f) { + dist_ += d; + angle_ += a; + if (data_ != 0) data_->AddVal(f, 1); + } private: double dist_; ///< Used to calculate average distance of hydrogen bond double angle_; ///< Used to calculate average angle of hydrogen bond diff --git a/src/AtomMask.h b/src/AtomMask.h index 54577c61d8..e92dd6d876 100644 --- a/src/AtomMask.h +++ b/src/AtomMask.h @@ -37,6 +37,8 @@ class AtomMask : public MaskTokenArray { int back() const { return Selected_.back(); } /// \return selected atom at idx const int& operator[](int idx) const { return Selected_[idx]; } + /// \return Number of atoms from corresponding topology. + int NmaskAtoms() const { return Natom_; } /// Flip current mask expression. void InvertMaskExpression(); /// Invert current mask From 452629942f7037fa51fdcb3678bbd7c7f3cbfc0f Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 12 Apr 2017 11:07:41 -0400 Subject: [PATCH 05/56] DRR - Cpptraj: Ensure donor-only sites will be calculated. --- src/Action_HydrogenBond.cpp | 23 +++++++++++++++-------- src/Action_HydrogenBond.h | 4 ++-- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 2ca7d17e2a..28ebb56158 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -227,6 +227,7 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { AcceptorMask_.SetNatoms( Mask_.NmaskAtoms() ); } + Sarray donorOnly; // DONOR/ACCEPTOR SETUP if (hasDonorMask_) { // Donor heavy atom mask specified @@ -270,7 +271,7 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { if (isDonor && isAcceptor) Both_.push_back( Site(at, H_at) ); else if (isDonor) - Donor_.push_back( Site(at, H_at) ); + donorOnly.push_back( Site(at, H_at) ); else if (isAcceptor) Acceptor_.push_back( at ); } @@ -300,7 +301,7 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { if (isDonor && isAcceptor) Both_.push_back( Site(at, Hatoms) ); else if (isDonor) - Donor_.push_back( Site(at, Hatoms) ); + donorOnly.push_back( Site(at, Hatoms) ); else if (isAcceptor) Acceptor_.push_back( at ); } @@ -333,11 +334,15 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { if (isDonor && isAcceptor) Both_.push_back( Site(at, Hatoms) ); else if (isDonor) - Donor_.push_back( Site(at, Hatoms) ); + donorOnly.push_back( Site(at, Hatoms) ); else if (isAcceptor) Acceptor_.push_back( at ); } } + // Place donor-only sites at the end of Both_ + bothEnd_ = Both_.size(); + for (Sarray::const_iterator site = donorOnly.begin(); site != donorOnly.end(); ++site) + Both_.push_back( *site ); if (calcSolvent_) { // Set up solvent donor/acceptor masks @@ -360,15 +365,16 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { mprintf("Acceptor atoms (%zu):\n", Acceptor_.size()); for (Iarray::const_iterator at = Acceptor_.begin(); at != Acceptor_.end(); ++at) mprintf("\t%20s %8i\n", setup.Top().TruncResAtomName(*at).c_str(), *at+1); - mprintf("Donor/acceptor sites (%zu):\n", Both_.size()); - for (Sarray::const_iterator si = Both_.begin(); si != Both_.end(); ++si) { + mprintf("Donor/acceptor sites (%u):\n", bothEnd_); + Sarray::const_iterator END = Both_.begin() + bothEnd_; + for (Sarray::const_iterator si = Both_.begin(); si != END; ++si) { mprintf("\t%20s %8i", setup.Top().TruncResAtomName(si->Idx()).c_str(), si->Idx()+1); for (Iarray::const_iterator at = si->Hbegin(); at != si->Hend(); ++at) mprintf(" %s", setup.Top()[*at].c_str()); mprintf("\n"); } - mprintf("Donor sites (%zu):\n", Donor_.size()); - for (Sarray::const_iterator si = Donor_.begin(); si != Donor_.end(); ++si) { + mprintf("Donor sites (%zu):\n", Both_.size() - bothEnd_); + for (Sarray::const_iterator si = END; si != Both_.end(); ++si) { mprintf("\t%20s %8i", setup.Top().TruncResAtomName(si->Idx()).c_str(), si->Idx()+1); for (Iarray::const_iterator at = si->Hbegin(); at != si->Hend(); ++at) mprintf(" %s", setup.Top()[*at].c_str()); @@ -448,11 +454,12 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { int numHB = 0; int hbidx = 0; + // Loop over all donor sites for (unsigned int sidx0 = 0; sidx0 != Both_.size(); sidx0++) { const double* XYZ0 = frm.Frm().XYZ( Both_[sidx0].Idx() ); // Loop over sites that can be both donor and acceptor - for (unsigned int sidx1 = sidx0 + 1; sidx1 != Both_.size(); sidx1++, hbidx++) + for (unsigned int sidx1 = sidx0 + 1; sidx1 < bothEnd_; sidx1++, hbidx++) { const double* XYZ1 = frm.Frm().XYZ( Both_[sidx1].Idx() ); double dist2 = DIST2( XYZ0, XYZ1, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 002fdf3b07..c070a26306 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -33,8 +33,7 @@ class Action_HydrogenBond : public Action { typedef std::map HBmapType; ImagedAction Image_; ///< Hold imaging info. - Sarray Donor_; ///< Array of sites that are just donor. - Sarray Both_; ///< Array of sites that are donor and acceptor + Sarray Both_; ///< Array of donor sites that can also be acceptors Iarray Acceptor_; ///< Array of acceptor-only atom indices HBmapType UU_Map_; @@ -62,6 +61,7 @@ class Action_HydrogenBond : public Action { CpptrajFile* bridgeout_; double dcut2_; double acut_; + unsigned int bothEnd_; ///< Index in Both_ where donor-only sites begin int debug_; bool series_; bool useAtomNum_; From 78f12eb0bf866bc50a84efaea469331bc26f3a00 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 12 Apr 2017 11:55:15 -0400 Subject: [PATCH 06/56] DRR - Cpptraj: Add print code. MPI code added too but not yet working. --- src/Action_HydrogenBond.cpp | 381 +++++++++++++++++++++++++++++++++++++++++++- src/Action_HydrogenBond.h | 34 +++- 2 files changed, 412 insertions(+), 3 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 28ebb56158..1867085cd8 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -1,8 +1,10 @@ #include // sqrt +#include // sort TODO may be able to do without this #include "Action_HydrogenBond.h" #include "CpptrajStdio.h" #include "Constants.h" #include "TorsionRoutines.h" +#include "StringRoutines.h" // ByteString // CONSTRUCTOR Action_HydrogenBond::Action_HydrogenBond() : @@ -19,6 +21,7 @@ Action_HydrogenBond::Action_HydrogenBond() : bridgeout_(0), dcut2_(0.0), acut_(0.0), + Nframes_(0), debug_(0), series_(0), useAtomNum_(false), @@ -455,6 +458,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { int numHB = 0; int hbidx = 0; // Loop over all donor sites + t_uu_.Start(); for (unsigned int sidx0 = 0; sidx0 != Both_.size(); sidx0++) { const double* XYZ0 = frm.Frm().XYZ( Both_[sidx0].Idx() ); @@ -480,15 +484,390 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { CalcSiteHbonds(hbidx, frameNum, dist2, Both_[sidx0], XYZ0, *a_atom, XYZ1, frm.Frm(), numHB); } } + t_uu_.Stop(); NumHbonds_->Add(frameNum, &numHB); + Nframes_++; t_action_.Stop(); 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()); +} + +// Action_HydrogenBond::SyncMap +void Action_HydrogenBond::SyncMap(HBmapType& mapIn, std::vector const& rank_frames, + std::vector const& rank_offsets, + const char* aspect, Parallel::Comm const& commIn) +const +{ + // Need to know how many hbonds on each thread. + int num_hb = (int)mapIn.size(); + std::vector nhb_on_rank; + if (commIn.Master()) + nhb_on_rank.resize( commIn.Size() ); + commIn.GatherMaster( &num_hb, 1, MPI_INT, &nhb_on_rank[0] ); + std::vector dArray; + std::vector iArray; + if (commIn.Master()) { + for (int rank = 1; rank < commIn.Size(); rank++) { + if (nhb_on_rank[rank] > 0) { + //mprintf("DEBUG:\tReceiving %i hbonds from rank %i.\n", nhb_on_rank[rank], rank); + dArray.resize( 2 * nhb_on_rank[rank] ); + iArray.resize( 5 * nhb_on_rank[rank] ); + commIn.Recv( &(dArray[0]), dArray.size(), MPI_DOUBLE, rank, 1300 ); + commIn.Recv( &(iArray[0]), iArray.size(), MPI_INT, rank, 1301 ); + HbondType HB; + int ii = 0, id = 0; + for (int in = 0; in != nhb_on_rank[rank]; in++, ii += 5, id += 2) { + HBmapType::iterator it = mapIn.find( iArray[ii] ); // hbidx + if (it == mapIn.end() ) { + // Hbond on rank that has not been found on master + HB.dist = dArray[id ]; + HB.angle = dArray[id+1]; + HB.A = iArray[ii+1]; + HB.H = iArray[ii+2]; + HB.D = iArray[ii+3]; + HB.Frames = iArray[ii+4]; + HB.data_ = 0; + //mprintf("\tNEW Hbond %i: %i-%i-%i D=%g A=%g %i frames", iArray[ii], + // HB.A+1, HB.H+1, HB.D+1, HB.dist, + // HB.angle, HB.Frames); + if (series_) { + HB.data_ = (DataSet_integer*) + masterDSL_->AddSet( DataSet::INTEGER, + MetaData(hbsetname_, aspect, iArray[ii]) ); + // FIXME: This may be incorrect if CurrentParm_ has changed + HB.data_->SetLegend( CreateHBlegend(*CurrentParm_, HB.A, HB.H, HB.D) ); + // mprintf(" \"%s\"", HB.data_->legend()); + } + //mprintf("\n"); + mapIn.insert( it, std::pair(iArray[ii], HB) ); + } else { + // Hbond on rank and master. Update on master. + //mprintf("\tAPPENDING Hbond %i: %i-%i-%i D=%g A=%g %i frames\n", iArray[ii], + // it->second.A+1, it->second.H+1, it->second.D+1, dArray[id], + // dArray[id+1], iArray[ii+4]); + it->second.dist += dArray[id ]; + it->second.angle += dArray[id+1]; + it->second.Frames += iArray[ii+4]; + if (series_) + HB.data_ = it->second.data_; + } + if (series_) { + HB.data_->Resize( Nframes_ ); + int* d_beg = HB.data_->Ptr() + rank_offsets[ rank ]; + //mprintf("\tResizing hbond series data to %i, starting frame %i, # frames %i\n", + // Nframes_, rank_offsets[rank], rank_frames[rank]); + commIn.Recv( d_beg, rank_frames[ rank ], MPI_INT, rank, 1304 + in ); + HB.data_->SetNeedsSync( false ); + } + } // END master loop over hbonds from rank + } + } // END master loop over ranks + // At this point we have all hbond sets from all ranks. Mark all HB sets + // smaller than Nframes_ as synced and ensure the time series has been + // updated to reflect overall # frames. + if (series_) { + const int ZERO = 0; + for (HBmapType::iterator hb = mapIn.begin(); hb != mapIn.end(); ++hb) + if ((int)hb->second.data_->Size() < Nframes_) { + hb->second.data_->SetNeedsSync( false ); + hb->second.data_->Add( Nframes_-1, &ZERO ); + } + } + } else { + if (mapIn.size() > 0) { + dArray.reserve( 2 * mapIn.size() ); + iArray.reserve( 5 * mapIn.size() ); + for (HBmapType::const_iterator hb = mapIn.begin(); hb != mapIn.end(); ++hb) { + dArray.push_back( hb->second.dist ); + dArray.push_back( hb->second.angle ); + iArray.push_back( hb->first ); + iArray.push_back( hb->second.A ); + iArray.push_back( hb->second.H ); + iArray.push_back( hb->second.D ); + iArray.push_back( hb->second.Frames ); + } + commIn.Send( &(dArray[0]), dArray.size(), MPI_DOUBLE, 0, 1300 ); + commIn.Send( &(iArray[0]), iArray.size(), MPI_INT, 0, 1301 ); + // Send series data to master + if (series_) { + int in = 0; // For tag + for (HBmapType::const_iterator hb = mapIn.begin(); hb != mapIn.end(); ++hb, in++) { + commIn.Send( hb->second.data_->Ptr(), hb->second.data_->Size(), + MPI_INT, 0, 1304 + in ); + hb->second.data_->SetNeedsSync( false ); + } + } + } + } +} + +/** PARALLEL NOTES: + * The following tags are used for MPI send/receive: + * 1300 : Array containing hbond double info on rank. + * 1301 : Array containing hbond integer info on rank. + * 1302 : Number of bridges to expect from rank. + * 1303 : Array containing bridge integer info on rank. + * 1304+X: Array of hbond X series info from rank. + */ +int Action_HydrogenBond::SyncAction() { + // Make sure all time series are updated at this point. + UpdateSeries(); + // TODO consolidate # frames / offset calc code with Action_NAstruct + // Get total number of frames + std::vector rank_frames( trajComm_.Size() ); + trajComm_.GatherMaster( &Nframes_, 1, MPI_INT, &rank_frames[0] ); + if (trajComm_.Master()) { + for (int rank = 1; rank < trajComm_.Size(); rank++) + Nframes_ += rank_frames[rank]; + } + // Convert rank frames to offsets. + std::vector rank_offsets( trajComm_.Size(), 0 ); + if (trajComm_.Master()) { + for (int rank = 1; rank < trajComm_.Size(); rank++) + rank_offsets[rank] = rank_offsets[rank-1] + rank_frames[rank-1]; + } + // Need to send hbond data from all ranks to master. + SyncMap( UU_Map_, rank_frames, rank_offsets, "solutehb", trajComm_ ); + if (calcSolvent_) { + SyncMap( UV_Map_, rank_frames, rank_offsets, "solventhb", trajComm_ ); + // iArray will contain for each bridge: Nres, res1, ..., resN, Frames + std::vector iArray; + int iSize; + if (trajComm_.Master()) { + for (int rank = 1; rank < trajComm_.Size(); rank++) + { + // Receive size of iArray + trajComm_.Recv( &iSize, 1, MPI_INT, rank, 1302 ); + //mprintf("DEBUG: Receiving %i bridges from rank %i\n", iSize, rank); + iArray.resize( iSize ); + trajComm_.Recv( &(iArray[0]), iSize, MPI_INT, rank, 1303 ); + unsigned int idx = 0; + while (idx < iArray.size()) { + std::set residues; + unsigned int i2 = idx + 1; + for (int ir = 0; ir != iArray[idx]; ir++, i2++) + residues.insert( iArray[i2] ); + BridgeType::iterator b_it = BridgeMap_.find( residues ); + if (b_it == BridgeMap_.end() ) // New Bridge + BridgeMap_.insert( b_it, std::pair,int>(residues, iArray[i2]) ); + else // Increment bridge #frames + b_it->second += iArray[i2]; + idx = i2 + 1; + } + } + } else { + // Construct bridge info array. + for (BridgeType::const_iterator b = BridgeMap_.begin(); b != BridgeMap_.end(); ++b) + { + iArray.push_back( b->first.size() ); // # of bridging res + 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 + } + // 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 ); + } + } + return 0; +} +#endif + +// Action_HydrogenBond::UpdateSeries() +void Action_HydrogenBond::UpdateSeries() { + if (seriesUpdated_) return; + if (series_ && Nframes_ > 0) { + for (HBmapType::iterator hb = UU_Map_.begin(); hb != UU_Map_.end(); ++hb) + hb->second.FinishSeries(Nframes_); + for (HBmapType::iterator hb = UV_Map_.begin(); hb != UV_Map_.end(); ++hb) + hb->second.FinishSeries(Nframes_); + } + // Should only be called once. + seriesUpdated_ = true; +} + +// Action_Hbond::MemoryUsage() +std::string Action_HydrogenBond::MemoryUsage(size_t nPairs, size_t nFrames) const { + static const size_t HBmapTypeElt = 32 + sizeof(int) + + (2*sizeof(double) + sizeof(DataSet_integer*) + 4*sizeof(int)); + size_t memTotal = nPairs * HBmapTypeElt; + // If calculating series every hbond will have time series. + // NOTE: This does not include memory used by DataSet. + if (series_ && nFrames > 0) { + size_t seriesSet = (nFrames * sizeof(int)) + sizeof(std::vector); + memTotal += (seriesSet * nPairs); + } +/* + // Current memory used by bridging solvent + static const size_t BridgeTypeElt = 32 + sizeof(std::set) + sizeof(int); + for (BridgeType::const_iterator it = BridgeMap_.begin(); it != BridgeMap_.end(); ++it) + memTotal += (it->first.size() * sizeof(int)); + memTotal += (BridgeMap_.size() * BridgeTypeElt); +*/ + return ByteString( memTotal, BYTE_DECIMAL ); +} + // Action_HydrogenBond::Print() +/** Print average occupancies over all frames for all detected Hbonds. */ void Action_HydrogenBond::Print() { - mprintf(" HBOND:\n"); + typedef std::vector Harray; + Harray HbondList; // For sorting + std::string Aname, Hname, Dname; + + // Final memory usage + mprintf(" HBOND: Actual memory usage is %s\n", + MemoryUsage(UU_Map_.size()+UV_Map_.size(), Nframes_).c_str()); + mprintf("\t%zu solute-solute hydrogen bonds.\n", UU_Map_.size()); + if (calcSolvent_) { + mprintf("\t%zu solute-solvent hydrogen bonds.\n", UV_Map_.size()); +// mprintf("\t%zu unique solute-solvent bridging interactions.\n", BridgeMap_.size()); + } + + t_uu_.WriteTiming( 2,"Solute-Solute :",t_action_.Total()); + if (calcSolvent_) { + t_ud_va_.WriteTiming( 2,"Solute Donor-Solvent Acceptor :",t_action_.Total()); + t_vd_ua_.WriteTiming( 2,"Solvent Donor-Solute Acceptor :",t_action_.Total()); + t_bridge_.WriteTiming(2,"Bridging waters :",t_action_.Total()); + } t_action_.WriteTiming(1,"Total:"); + + // Ensure all series have been updated for all frames. + UpdateSeries(); + + if (CurrentParm_ == 0) return; + // Calculate necessary column width for strings based on how many residues. + // ResName+'_'+ResNum+'@'+AtomName | NUM = 4+1+R+1+4 = R+10 + int NUM = DigitWidth( CurrentParm_->Nres() ) + 10; + // If useAtomNum_ +'_'+AtomNum += 1+A + if (useAtomNum_) NUM += ( DigitWidth( CurrentParm_->Natom() ) + 1 ); + + // Solute Hbonds + if (avgout_ != 0) { + // Place all detected Hbonds in a list and sort. + for (HBmapType::const_iterator it = UU_Map_.begin(); it != UU_Map_.end(); ++it) { + HbondList.push_back( it->second ); + // Calculate average distance and angle for this hbond. + HbondList.back().CalcAvg(); + } + 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"); + for (Harray::const_iterator hbond = HbondList.begin(); hbond != HbondList.end(); ++hbond ) + { + double avg = ((double)hbond->Frames()) / ((double) Nframes_); + Aname = CurrentParm_->TruncResAtomName(hbond->A()); + Hname = CurrentParm_->TruncResAtomName(hbond->H()); + Dname = CurrentParm_->TruncResAtomName(hbond->D()); + if (useAtomNum_) { + Aname.append("_" + integerToString(hbond->A()+1)); + 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", + NUM, Aname.c_str(), NUM, Hname.c_str(), NUM, Dname.c_str(), + hbond->Frames(), avg, hbond->Dist(), hbond->Angle()); + } + } + + // Solute-solvent Hbonds + if (solvout_ != 0 && calcSolvent_) { + HbondList.clear(); + for (HBmapType::const_iterator it = UV_Map_.begin(); it != UV_Map_.end(); ++it) { + HbondList.push_back( it->second ); + // Calculate average distance and angle for this hbond. + HbondList.back().CalcAvg(); + } + UV_Map_.clear(); + 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"); + for (Harray::const_iterator hbond = HbondList.begin(); hbond != HbondList.end(); ++hbond ) + { + // Average has slightly diff meaning since for any given frame multiple + // solvent can bond to the same solute. + double avg = ((double)hbond->Frames()) / ((double) Nframes_); + if (hbond->A()==-1) // Solvent acceptor + Aname = "SolventAcc"; + else { + Aname = CurrentParm_->TruncResAtomName(hbond->A()); + if (useAtomNum_) Aname.append("_" + integerToString(hbond->A()+1)); + } + if (hbond->D()==-1) { // Solvent donor + Dname = "SolventDnr"; + Hname = "SolventH"; + } else { + Dname = CurrentParm_->TruncResAtomName(hbond->D()); + Hname = CurrentParm_->TruncResAtomName(hbond->H()); + if (useAtomNum_) { + Dname.append("_" + integerToString(hbond->D()+1)); + Hname.append("_" + integerToString(hbond->H()+1)); + } + } + solvout_->Printf("%-*s %*s %*s %8i %12.4f %12.4f %12.4f\n", + NUM, Aname.c_str(), NUM, Hname.c_str(), NUM, Dname.c_str(), + hbond->Frames(), avg, hbond->Dist(), hbond->Angle()); + } + } +/* + // BRIDGING INFO + if (bridgeout_ != 0 && calcSolvent_) { + bridgeout_->Printf("#Bridging Solute Residues:\n"); + // Place bridging values in a vector for sorting + std::vector, int> > bridgevector; + for (BridgeType::const_iterator it = BridgeMap_.begin(); + it != BridgeMap_.end(); ++it) + bridgevector.push_back( *it ); + std::sort( bridgevector.begin(), bridgevector.end(), bridge_cmp() ); + for (std::vector, int> >::const_iterator bv = bridgevector.begin(); + bv != bridgevector.end(); + ++bv) + { + bridgeout_->Printf("Bridge Res"); + for (std::set::const_iterator res = bv->first.begin(); + res != bv->first.end(); ++res) + bridgeout_->Printf(" %i:%s", *res+1, CurrentParm_->Res( *res ).c_str()); + bridgeout_->Printf(", %i frames.\n", bv->second); + } + } +*/ +} + +// ===== Action_HydrogenBond::Hbond ============================================ +/** Calculate average distance and angle for hbond. */ +void Action_HydrogenBond::Hbond::CalcAvg() { + double dFrames = (double)frames_; + dist_ /= dFrames; + angle_ /= dFrames; + 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 ); + } } diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index c070a26306..1f81be3d87 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -28,6 +28,10 @@ class Action_HydrogenBond : public Action { inline double Angle(const double*, const double*, const double*) const; void CalcSiteHbonds(int,int,double,Site const&,const double*,int,const double*, Frame const&, int&); + /// Update all hydrogen bond time series + void UpdateSeries(); + /// Determine memory usage from # hbonds and time series + std::string MemoryUsage(size_t, size_t) const; typedef std::vector Sarray; typedef std::map HBmapType; @@ -48,6 +52,10 @@ class Action_HydrogenBond : public Action { AtomMask Mask_; Matrix_3x3 ucell_, recip_; Timer t_action_; + Timer t_uu_; + Timer t_ud_va_; + Timer t_vd_ua_; + Timer t_bridge_; Topology* CurrentParm_; ///< Used to set atom/residue labels DataSetList* masterDSL_; DataSet* NumHbonds_; @@ -62,9 +70,11 @@ class Action_HydrogenBond : public Action { double dcut2_; double acut_; unsigned int bothEnd_; ///< Index in Both_ where donor-only sites begin + int Nframes_; ///< Number of frames action has been active int debug_; - bool series_; - bool useAtomNum_; + bool series_; ///< If true track hbond 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_; bool hasDonorMask_; bool hasDonorHmask_; @@ -75,6 +85,7 @@ class Action_HydrogenBond : public Action { }; // ----- CLASSES --------------------------------------------------------------- +/// Potential hydrogen bond site. Can be either donor or donor/acceptor. class Action_HydrogenBond::Site { public: Site() : idx_(-1), isV_(false) {} @@ -94,18 +105,37 @@ class Action_HydrogenBond::Site { bool isV_; ///< True if site is solvent }; +/// 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) {} + 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_; } + /// First sort by frames (descending), then distance (ascending). + bool operator<(const Hbond& rhs) const { + if (frames_ == rhs.frames_) + return (dist_ < rhs.dist_); + else + return (frames_ > rhs.frames_); + } + /// Update distance/angle/time series void Update(double d, double a, int f) { dist_ += d; angle_ += a; + ++frames_; if (data_ != 0) data_->AddVal(f, 1); } + 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 From 023e553c1ba32d82c70f46b2a9b611459490333c Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 12 Apr 2017 12:39:20 -0400 Subject: [PATCH 07/56] DRR - Cpptraj: Ensure UU hydrogen bonds in map are truely unique by using donor H atom - acceptor atom indices as key. --- src/Action_HydrogenBond.cpp | 25 ++++++++++++++----------- src/Action_HydrogenBond.h | 5 +++-- 2 files changed, 17 insertions(+), 13 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 1867085cd8..3c402cf6e2 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -412,7 +412,7 @@ double Action_HydrogenBond::Angle(const double* XA, const double* XH, const doub } // Action_HydrogenBond::CalcSiteHbonds() -void Action_HydrogenBond::CalcSiteHbonds(int hbidx, int frameNum, double dist2, +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) @@ -428,23 +428,27 @@ void Action_HydrogenBond::CalcSiteHbonds(int hbidx, int frameNum, double dist2, ++numHB; double dist = sqrt(dist2); // Index UU hydrogen bonds by DonorH-Acceptor + Hpair hbidx(*h_atom, a_atom); HBmapType::iterator it = UU_Map_.find( hbidx ); if (it == UU_Map_.end()) { +// mprintf("DBG1: NEW hbond : %8i .. %8i - %8i\n", a_atom+1,*h_atom+1,d_atom+1); DataSet_integer* ds = 0; if (series_) { std::string hblegend = CurrentParm_->TruncResAtomName(a_atom) + "-" + CurrentParm_->TruncResAtomName(d_atom) + "-" + (*CurrentParm_)[*h_atom].Name().Truncated(); ds = (DataSet_integer*) - masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"solutehb",hbidx)); + masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"solutehb",UU_Map_.size())); if (UUseriesout_ != 0) UUseriesout_->AddDataSet( ds ); ds->SetLegend( hblegend ); ds->AddVal( frameNum, 1 ); } - UU_Map_.insert(it, std::pair(hbidx,Hbond(dist,angle,ds,a_atom,*h_atom,d_atom))); - } else + UU_Map_.insert(it, std::pair(hbidx,Hbond(dist,angle,ds,a_atom,*h_atom,d_atom))); + } else { +// mprintf("DBG1: OLD hbond : %8i .. %8i - %8i\n", a_atom+1,*h_atom+1,d_atom+1); it->second.Update(dist,angle,frameNum); + } } } } @@ -455,33 +459,32 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { if (Image_.ImagingEnabled()) frm.Frm().BoxCrd().ToRecip(ucell_, recip_); - int numHB = 0; - int hbidx = 0; // Loop over all donor sites + int numHB = 0; t_uu_.Start(); for (unsigned int sidx0 = 0; sidx0 != Both_.size(); sidx0++) { const double* XYZ0 = frm.Frm().XYZ( Both_[sidx0].Idx() ); // Loop over sites that can be both donor and acceptor - for (unsigned int sidx1 = sidx0 + 1; sidx1 < bothEnd_; sidx1++, hbidx++) + for (unsigned int sidx1 = sidx0 + 1; sidx1 < bothEnd_; sidx1++) { const double* XYZ1 = frm.Frm().XYZ( Both_[sidx1].Idx() ); double dist2 = DIST2( XYZ0, XYZ1, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); if ( !(dist2 > dcut2_) ) { // Site 0 donor, Site 1 acceptor - CalcSiteHbonds(hbidx, frameNum, dist2, Both_[sidx0], XYZ0, Both_[sidx1].Idx(), XYZ1, frm.Frm(), numHB); + CalcSiteHbonds(frameNum, dist2, Both_[sidx0], XYZ0, Both_[sidx1].Idx(), XYZ1, frm.Frm(), numHB); // Site 1 donor, Site 0 acceptor - CalcSiteHbonds(hbidx, frameNum, dist2, Both_[sidx1], XYZ1, Both_[sidx0].Idx(), XYZ0, frm.Frm(), numHB); + CalcSiteHbonds(frameNum, dist2, Both_[sidx1], XYZ1, Both_[sidx0].Idx(), XYZ0, frm.Frm(), numHB); } } // Loop over acceptor-only - for (Iarray::const_iterator a_atom = Acceptor_.begin(); a_atom != Acceptor_.end(); ++a_atom, hbidx++) + for (Iarray::const_iterator a_atom = Acceptor_.begin(); a_atom != Acceptor_.end(); ++a_atom) { const double* XYZ1 = frm.Frm().XYZ( *a_atom ); double dist2 = DIST2( XYZ0, XYZ1, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); if ( !(dist2 > dcut2_) ) - CalcSiteHbonds(hbidx, frameNum, dist2, Both_[sidx0], XYZ0, *a_atom, XYZ1, frm.Frm(), numHB); + CalcSiteHbonds(frameNum, dist2, Both_[sidx0], XYZ0, *a_atom, XYZ1, frm.Frm(), numHB); } } t_uu_.Stop(); diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 1f81be3d87..7e668fc8a4 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -26,7 +26,7 @@ class Action_HydrogenBond : public Action { class Hbond; inline double Angle(const double*, const double*, const double*) const; - void CalcSiteHbonds(int,int,double,Site const&,const double*,int,const double*, + void CalcSiteHbonds(int,double,Site const&,const double*,int,const double*, Frame const&, int&); /// Update all hydrogen bond time series void UpdateSeries(); @@ -34,7 +34,8 @@ class Action_HydrogenBond : public Action { std::string MemoryUsage(size_t, size_t) const; typedef std::vector Sarray; - typedef std::map HBmapType; + typedef std::pair Hpair; + typedef std::map HBmapType; ImagedAction Image_; ///< Hold imaging info. Sarray Both_; ///< Array of donor sites that can also be acceptors From 3005b652e2e997e26be80e970923c57a61e421b7 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Wed, 12 Apr 2017 12:51:43 -0400 Subject: [PATCH 08/56] DRR - Cpptraj: Try to improve insertion speed into hbond map --- src/Action_HydrogenBond.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 3c402cf6e2..83ee9b7651 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -429,8 +429,8 @@ void Action_HydrogenBond::CalcSiteHbonds(int frameNum, double dist2, double dist = sqrt(dist2); // Index UU hydrogen bonds by DonorH-Acceptor Hpair hbidx(*h_atom, a_atom); - HBmapType::iterator it = UU_Map_.find( hbidx ); - if (it == UU_Map_.end()) + HBmapType::iterator it = UU_Map_.lower_bound( hbidx ); + if (it == UU_Map_.end() || it->first != hbidx) { // mprintf("DBG1: NEW hbond : %8i .. %8i - %8i\n", a_atom+1,*h_atom+1,d_atom+1); DataSet_integer* ds = 0; From 2a613e99f6c3ea6bb53edf7a85a846300f975b5e Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 08:10:23 -0400 Subject: [PATCH 09/56] DRR - Cpptraj: Add solvent Hbond calc. Speedup in initial testing is about 3.4x so far. --- src/Action_HydrogenBond.cpp | 202 +++++++++++++++++++++++++++++++++++++------- src/Action_HydrogenBond.h | 9 +- 2 files changed, 179 insertions(+), 32 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 83ee9b7651..96fd1118ce 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -347,24 +347,6 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { for (Sarray::const_iterator site = donorOnly.begin(); site != donorOnly.end(); ++site) Both_.push_back( *site ); - if (calcSolvent_) { - // Set up solvent donor/acceptor masks - if (hasSolventDonor_) { - if (setup.Top().SetupIntegerMask( SolventDonorMask_ )) return Action::ERR; - if (SolventDonorMask_.None()) { - mprintf("Warning: SolventDonorMask has no atoms.\n"); - return Action::SKIP; - } - } - if (hasSolventAcceptor_) { - if (setup.Top().SetupIntegerMask( SolventAcceptorMask_ )) return Action::ERR; - if (SolventAcceptorMask_.None()) { - mprintf("Warning: SolventAcceptorMask has no atoms.\n"); - return Action::SKIP; - } - } - } - mprintf("Acceptor atoms (%zu):\n", Acceptor_.size()); for (Iarray::const_iterator at = Acceptor_.begin(); at != Acceptor_.end(); ++at) mprintf("\t%20s %8i\n", setup.Top().TruncResAtomName(*at).c_str(), *at+1); @@ -384,6 +366,67 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { mprintf("\n"); } + if (calcSolvent_) { + int at_beg = 0; + int at_end = 0; + // Set up solvent donor/acceptor masks + if (hasSolventDonor_) { + if (setup.Top().SetupIntegerMask( SolventDonorMask_ )) return Action::ERR; + if (SolventDonorMask_.None()) { + mprintf("Warning: SolventDonorMask has no atoms.\n"); + return Action::SKIP; + } + at_beg = SolventDonorMask_[0]; + at_end = SolventDonorMask_.back() + 1; + } + if (hasSolventAcceptor_) { + if (setup.Top().SetupIntegerMask( SolventAcceptorMask_ )) return Action::ERR; + if (SolventAcceptorMask_.None()) { + mprintf("Warning: SolventAcceptorMask has no atoms.\n"); + return Action::SKIP; + } + if (!hasSolventDonor_) { + at_beg = SolventAcceptorMask_[0]; + at_end = SolventAcceptorMask_.back() + 1; + } else { + at_beg = std::min( SolventDonorMask_[0], SolventAcceptorMask_[0] ); + at_end = std::max( SolventDonorMask_.back(), SolventAcceptorMask_.back() ) + 1; + } + } + AtomMask::const_iterator a_atom = SolventAcceptorMask_.begin(); + AtomMask::const_iterator d_atom = SolventDonorMask_.begin(); + bool isDonor, isAcceptor; + for (int at = at_beg; at != at_end; at++) + { + Iarray Hatoms; + isDonor = false; + isAcceptor = false; + if ( d_atom != SolventDonorMask_.end() && *d_atom == at ) { + ++d_atom; + if ( IsFON( setup.Top()[at] ) ) { + for (Atom::bond_iterator H_at = setup.Top()[at].bondbegin(); + H_at != setup.Top()[at].bondend(); ++H_at) + if (setup.Top()[*H_at].Element() == Atom::HYDROGEN) + Hatoms.push_back( *H_at ); + } + isDonor = !Hatoms.empty(); + } + if ( a_atom != SolventAcceptorMask_.end() && *a_atom == at ) { + isAcceptor = true; + ++a_atom; + } + if (isDonor || isAcceptor) + SolventSites_.push_back( Site(at, Hatoms) ); + } + + mprintf("Solvent sites (%zu):\n", SolventSites_.size()); + for (Sarray::const_iterator si = SolventSites_.begin(); si != SolventSites_.end(); ++si) { + mprintf("\t%20s %8i", setup.Top().TruncResAtomName(si->Idx()).c_str(), si->Idx()+1); + for (Iarray::const_iterator at = si->Hbegin(); at != si->Hend(); ++at) + mprintf(" %s", setup.Top()[*at].c_str()); + mprintf("\n"); + } + } return Action::OK; } @@ -391,6 +434,8 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { // Action_HydrogenBond::Angle() double Action_HydrogenBond::Angle(const double* XA, const double* XH, const double* XD) const { + if (acut_ < 0.0) // Indicates skip angle calc + return 0.0; if (Image_.ImageType() == NOIMAGE) return (CalcAngle(XA, XH, XD)); else { @@ -411,6 +456,63 @@ double Action_HydrogenBond::Angle(const double* XA, const double* XH, const doub } } +// Action_HydrogenBond::CalcSolvHbonds() +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) +{ + // The two sites are close enough to hydrogen bond. + int d_atom = SiteD.Idx(); + // Determine if angle cutoff is satisfied + for (Iarray::const_iterator h_atom = SiteD.Hbegin(); h_atom != SiteD.Hend(); ++h_atom) + { + double angle = 0.0; + bool angleSatisfied = true; + // For ions, donor atom will be same as h atom so no angle needed. + if (d_atom != *h_atom) { + angle = Angle(XYZA, frmIn.XYZ(*h_atom), XYZD); + angleSatisfied = !(angle < acut_); + } + if (angleSatisfied) + { + ++numHB; + double dist = sqrt(dist2); + // Index U-H .. V hydrogen bonds by solute H atom. + // Index U .. H-V hydrogen bonds by solute A atom. + int hbidx; + if (soluteDonor) + hbidx = *h_atom; + else + hbidx = a_atom; + UVmapType::iterator it = UV_Map_.lower_bound( hbidx ); + if (it == UV_Map_.end() || it->first != hbidx) + { +// mprintf("DBG1: NEW hbond : %8i .. %8i - %8i\n", a_atom+1,*h_atom+1,d_atom+1); + DataSet_integer* ds = 0; + if (series_) { + ds = (DataSet_integer*) + masterDSL_->AddSet(DataSet::INTEGER,MetaData(hbsetname_,"solventhb",hbidx)); + if (UVseriesout_ != 0) UVseriesout_->AddDataSet( ds ); + ds->AddVal( frameNum, 1 ); + } + Hbond hb; + if (soluteDonor) { // 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); + } 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); + } + 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,frameNum); + } + } + } +} + // Action_HydrogenBond::CalcSiteHbonds() void Action_HydrogenBond::CalcSiteHbonds(int frameNum, double dist2, Site const& SiteD, const double* XYZD, @@ -459,13 +561,13 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { if (Image_.ImagingEnabled()) frm.Frm().BoxCrd().ToRecip(ucell_, recip_); - // Loop over all donor sites - int numHB = 0; + // Loop over all solute donor sites t_uu_.Start(); + int numHB = 0; for (unsigned int sidx0 = 0; sidx0 != Both_.size(); sidx0++) { const double* XYZ0 = frm.Frm().XYZ( Both_[sidx0].Idx() ); - // Loop over sites that can be both donor and acceptor + // Loop over solute sites that can be both donor and acceptor for (unsigned int sidx1 = sidx0 + 1; sidx1 < bothEnd_; sidx1++) { const double* XYZ1 = frm.Frm().XYZ( Both_[sidx1].Idx() ); @@ -478,7 +580,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { CalcSiteHbonds(frameNum, dist2, Both_[sidx1], XYZ1, Both_[sidx0].Idx(), XYZ0, frm.Frm(), numHB); } } - // Loop over acceptor-only + // Loop over solute acceptor-only for (Iarray::const_iterator a_atom = Acceptor_.begin(); a_atom != Acceptor_.end(); ++a_atom) { const double* XYZ1 = frm.Frm().XYZ( *a_atom ); @@ -487,9 +589,52 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { CalcSiteHbonds(frameNum, dist2, Both_[sidx0], XYZ0, *a_atom, XYZ1, frm.Frm(), numHB); } } + NumHbonds_->Add(frameNum, &numHB); t_uu_.Stop(); - NumHbonds_->Add(frameNum, &numHB); + // Loop over all solvent sites + if (calcSolvent_) { + t_uv_.Start(); + numHB = 0; + for (unsigned int vidx = 0; vidx != SolventSites_.size(); vidx++) + { + Site const& Vsite = SolventSites_[vidx]; + const double* VXYZ = frm.Frm().XYZ( Vsite.Idx() ); + // Loop over solute sites that can be both donor and acceptor + for (unsigned int sidx = 0; sidx < bothEnd_; sidx++) + { + const double* UXYZ = frm.Frm().XYZ( Both_[sidx].Idx() ); + double dist2 = DIST2( VXYZ, UXYZ, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); + if ( !(dist2 > dcut2_) ) + { + // Solvent site donor, solute site acceptor + CalcSolvHbonds(frameNum, dist2, Vsite, VXYZ, Both_[sidx].Idx(), UXYZ, frm.Frm(), numHB, false); + // Solvent site acceptor, solute site donor + CalcSolvHbonds(frameNum, dist2, Both_[sidx], UXYZ, Vsite.Idx(), VXYZ, frm.Frm(), numHB, true); + } + } + // Loop over solute sites that are donor only + for (unsigned int sidx = bothEnd_; sidx < Both_.size(); sidx++) + { + const double* UXYZ = frm.Frm().XYZ( Both_[sidx].Idx() ); + double dist2 = DIST2( VXYZ, UXYZ, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); + if ( !(dist2 > dcut2_) ) + // Solvent site acceptor, solute site donor + CalcSolvHbonds(frameNum, dist2, Both_[sidx], UXYZ, Vsite.Idx(), VXYZ, frm.Frm(), numHB, true); + } + // Loop over solute sites that are acceptor only + for (Iarray::const_iterator a_atom = Acceptor_.begin(); a_atom != Acceptor_.end(); ++a_atom) + { + const double* UXYZ = frm.Frm().XYZ( *a_atom ); + double dist2 = DIST2( VXYZ, UXYZ, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); + if ( !(dist2 > dcut2_) ) + // Solvent site donor, solute site acceptor + CalcSolvHbonds(frameNum, dist2, Vsite, VXYZ, *a_atom, UXYZ, frm.Frm(), numHB, false); + } + } // END loop over solvent sites + NumSolvent_->Add(frameNum, &numHB); + t_uv_.Stop(); + } Nframes_++; t_action_.Stop(); @@ -699,7 +844,7 @@ void Action_HydrogenBond::UpdateSeries() { if (series_ && Nframes_ > 0) { for (HBmapType::iterator hb = UU_Map_.begin(); hb != UU_Map_.end(); ++hb) hb->second.FinishSeries(Nframes_); - for (HBmapType::iterator hb = UV_Map_.begin(); hb != UV_Map_.end(); ++hb) + for (UVmapType::iterator hb = UV_Map_.begin(); hb != UV_Map_.end(); ++hb) hb->second.FinishSeries(Nframes_); } // Should only be called once. @@ -743,11 +888,10 @@ void Action_HydrogenBond::Print() { // mprintf("\t%zu unique solute-solvent bridging interactions.\n", BridgeMap_.size()); } - t_uu_.WriteTiming( 2,"Solute-Solute :",t_action_.Total()); + t_uu_.WriteTiming( 2, "Solute-Solute :",t_action_.Total()); if (calcSolvent_) { - t_ud_va_.WriteTiming( 2,"Solute Donor-Solvent Acceptor :",t_action_.Total()); - t_vd_ua_.WriteTiming( 2,"Solvent Donor-Solute Acceptor :",t_action_.Total()); - t_bridge_.WriteTiming(2,"Bridging waters :",t_action_.Total()); + t_uv_.WriteTiming( 2, "Solute-Solvent :",t_uv_.Total()); + t_bridge_.WriteTiming(2,"Bridging waters :",t_action_.Total()); } t_action_.WriteTiming(1,"Total:"); @@ -794,7 +938,7 @@ void Action_HydrogenBond::Print() { // Solute-solvent Hbonds if (solvout_ != 0 && calcSolvent_) { HbondList.clear(); - for (HBmapType::const_iterator it = UV_Map_.begin(); it != UV_Map_.end(); ++it) { + for (UVmapType::const_iterator it = UV_Map_.begin(); it != UV_Map_.end(); ++it) { HbondList.push_back( it->second ); // Calculate average distance and angle for this hbond. HbondList.back().CalcAvg(); diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 7e668fc8a4..3661e270c0 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -28,6 +28,8 @@ class Action_HydrogenBond : public Action { inline double Angle(const double*, const double*, const double*) const; 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*, + Frame const&, int&, bool); /// Update all hydrogen bond time series void UpdateSeries(); /// Determine memory usage from # hbonds and time series @@ -36,13 +38,15 @@ class Action_HydrogenBond : public Action { typedef std::vector Sarray; typedef std::pair Hpair; typedef std::map HBmapType; + typedef std::map UVmapType; ImagedAction Image_; ///< Hold imaging info. Sarray Both_; ///< Array of donor sites that can also be acceptors Iarray Acceptor_; ///< Array of acceptor-only atom indices + Sarray SolventSites_; HBmapType UU_Map_; - HBmapType UV_Map_; + UVmapType UV_Map_; std::string hbsetname_; AtomMask DonorMask_; @@ -54,8 +58,7 @@ class Action_HydrogenBond : public Action { Matrix_3x3 ucell_, recip_; Timer t_action_; Timer t_uu_; - Timer t_ud_va_; - Timer t_vd_ua_; + Timer t_uv_; Timer t_bridge_; Topology* CurrentParm_; ///< Used to set atom/residue labels DataSetList* masterDSL_; From f91a99365936750dff5fb21b47b88ed94cbd7bf7 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 08:16:51 -0400 Subject: [PATCH 10/56] DRR - Cpptraj: HBmapType to UUmapType --- src/Action_HydrogenBond.cpp | 16 ++++++++-------- src/Action_HydrogenBond.h | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 96fd1118ce..7a54f86df5 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -531,7 +531,7 @@ void Action_HydrogenBond::CalcSiteHbonds(int frameNum, double dist2, double dist = sqrt(dist2); // Index UU hydrogen bonds by DonorH-Acceptor Hpair hbidx(*h_atom, a_atom); - HBmapType::iterator it = UU_Map_.lower_bound( hbidx ); + UUmapType::iterator it = UU_Map_.lower_bound( hbidx ); if (it == UU_Map_.end() || it->first != hbidx) { // mprintf("DBG1: NEW hbond : %8i .. %8i - %8i\n", a_atom+1,*h_atom+1,d_atom+1); @@ -656,7 +656,7 @@ static inline std::string CreateHBlegend(Topology const& topIn, int a_atom, int } // Action_HydrogenBond::SyncMap -void Action_HydrogenBond::SyncMap(HBmapType& mapIn, std::vector const& rank_frames, +void Action_HydrogenBond::SyncMap(UUmapType& mapIn, std::vector const& rank_frames, std::vector const& rank_offsets, const char* aspect, Parallel::Comm const& commIn) const @@ -680,7 +680,7 @@ const HbondType HB; int ii = 0, id = 0; for (int in = 0; in != nhb_on_rank[rank]; in++, ii += 5, id += 2) { - HBmapType::iterator it = mapIn.find( iArray[ii] ); // hbidx + UUmapType::iterator it = mapIn.find( iArray[ii] ); // hbidx if (it == mapIn.end() ) { // Hbond on rank that has not been found on master HB.dist = dArray[id ]; @@ -730,7 +730,7 @@ const // updated to reflect overall # frames. if (series_) { const int ZERO = 0; - for (HBmapType::iterator hb = mapIn.begin(); hb != mapIn.end(); ++hb) + for (UUmapType::iterator hb = mapIn.begin(); hb != mapIn.end(); ++hb) if ((int)hb->second.data_->Size() < Nframes_) { hb->second.data_->SetNeedsSync( false ); hb->second.data_->Add( Nframes_-1, &ZERO ); @@ -740,7 +740,7 @@ const if (mapIn.size() > 0) { dArray.reserve( 2 * mapIn.size() ); iArray.reserve( 5 * mapIn.size() ); - for (HBmapType::const_iterator hb = mapIn.begin(); hb != mapIn.end(); ++hb) { + for (UUmapType::const_iterator hb = mapIn.begin(); hb != mapIn.end(); ++hb) { dArray.push_back( hb->second.dist ); dArray.push_back( hb->second.angle ); iArray.push_back( hb->first ); @@ -754,7 +754,7 @@ const // Send series data to master if (series_) { int in = 0; // For tag - for (HBmapType::const_iterator hb = mapIn.begin(); hb != mapIn.end(); ++hb, in++) { + for (UUmapType::const_iterator hb = mapIn.begin(); hb != mapIn.end(); ++hb, in++) { commIn.Send( hb->second.data_->Ptr(), hb->second.data_->Size(), MPI_INT, 0, 1304 + in ); hb->second.data_->SetNeedsSync( false ); @@ -842,7 +842,7 @@ int Action_HydrogenBond::SyncAction() { void Action_HydrogenBond::UpdateSeries() { if (seriesUpdated_) return; if (series_ && Nframes_ > 0) { - for (HBmapType::iterator hb = UU_Map_.begin(); hb != UU_Map_.end(); ++hb) + for (UUmapType::iterator hb = UU_Map_.begin(); hb != UU_Map_.end(); ++hb) hb->second.FinishSeries(Nframes_); for (UVmapType::iterator hb = UV_Map_.begin(); hb != UV_Map_.end(); ++hb) hb->second.FinishSeries(Nframes_); @@ -908,7 +908,7 @@ void Action_HydrogenBond::Print() { // Solute Hbonds if (avgout_ != 0) { // Place all detected Hbonds in a list and sort. - for (HBmapType::const_iterator it = UU_Map_.begin(); it != UU_Map_.end(); ++it) { + for (UUmapType::const_iterator it = UU_Map_.begin(); it != UU_Map_.end(); ++it) { HbondList.push_back( it->second ); // Calculate average distance and angle for this hbond. HbondList.back().CalcAvg(); diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 3661e270c0..ad80a1c16e 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -37,7 +37,7 @@ class Action_HydrogenBond : public Action { typedef std::vector Sarray; typedef std::pair Hpair; - typedef std::map HBmapType; + typedef std::map UUmapType; typedef std::map UVmapType; ImagedAction Image_; ///< Hold imaging info. @@ -45,7 +45,7 @@ class Action_HydrogenBond : public Action { Iarray Acceptor_; ///< Array of acceptor-only atom indices Sarray SolventSites_; - HBmapType UU_Map_; + UUmapType UU_Map_; UVmapType UV_Map_; std::string hbsetname_; From f12ebf29521289b9c85d563b49f861eb7b5dfac6 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 08:35:38 -0400 Subject: [PATCH 11/56] DRR - Cpptraj: Ensure solvent is ignored when solute donor or acceptor mask not explicitly specified. --- src/Action_HydrogenBond.cpp | 54 ++++++++++++++++++++++++++------------------- src/Action_HydrogenBond.h | 2 +- 2 files changed, 32 insertions(+), 24 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 7a54f86df5..ae10ab143e 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -224,14 +224,16 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { // No specified acceptor mask; search generic mask. AcceptorMask_.ResetMask(); for (AtomMask::const_iterator at = Mask_.begin(); at != Mask_.end(); ++at) { - if (IsFON( setup.Top()[*at] )) + // Since an acceptor mask was not specified ignore solvent. + int molnum = setup.Top()[*at].MolNum(); + if (!setup.Top().Mol(molnum).IsSolvent() && IsFON( setup.Top()[*at] )) AcceptorMask_.AddSelectedAtom( *at ); } AcceptorMask_.SetNatoms( Mask_.NmaskAtoms() ); } + // SOLUTE DONOR/ACCEPTOR SITE SETUP Sarray donorOnly; - // DONOR/ACCEPTOR SETUP if (hasDonorMask_) { // Donor heavy atom mask specified if ( setup.Top().SetupIntegerMask( DonorMask_ ) ) return Action::ERR; @@ -317,29 +319,34 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { bool isDonor, isAcceptor; for (int at = 0; at != maxAtom; at++) { - Iarray Hatoms; - isDonor = false; - isAcceptor = false; - if ( d_atom != Mask_.end() && *d_atom == at) { - ++d_atom; - if ( IsFON( setup.Top()[at] ) ) { - for (Atom::bond_iterator H_at = setup.Top()[at].bondbegin(); - H_at != setup.Top()[at].bondend(); ++H_at) - if (setup.Top()[*H_at].Element() == Atom::HYDROGEN) - Hatoms.push_back( *H_at ); - isDonor = !Hatoms.empty(); + // Since an acceptor mask was not specified ignore solvent. + int molnum = setup.Top()[at].MolNum(); + if (!setup.Top().Mol(molnum).IsSolvent()) + { + Iarray Hatoms; + isDonor = false; + isAcceptor = false; + if ( d_atom != Mask_.end() && *d_atom == at) { + ++d_atom; + if ( IsFON( setup.Top()[at] ) ) { + for (Atom::bond_iterator H_at = setup.Top()[at].bondbegin(); + H_at != setup.Top()[at].bondend(); ++H_at) + if (setup.Top()[*H_at].Element() == Atom::HYDROGEN) + Hatoms.push_back( *H_at ); + isDonor = !Hatoms.empty(); + } } + if ( a_atom != AcceptorMask_.end() && *a_atom == at ) { + isAcceptor = true; + ++a_atom; + } + if (isDonor && isAcceptor) + Both_.push_back( Site(at, Hatoms) ); + else if (isDonor) + donorOnly.push_back( Site(at, Hatoms) ); + else if (isAcceptor) + Acceptor_.push_back( at ); } - if ( a_atom != AcceptorMask_.end() && *a_atom == at ) { - isAcceptor = true; - ++a_atom; - } - if (isDonor && isAcceptor) - Both_.push_back( Site(at, Hatoms) ); - else if (isDonor) - donorOnly.push_back( Site(at, Hatoms) ); - else if (isAcceptor) - Acceptor_.push_back( at ); } } // Place donor-only sites at the end of Both_ @@ -366,6 +373,7 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { mprintf("\n"); } + // SOLVENT SITE SETUP if (calcSolvent_) { int at_beg = 0; int at_end = 0; diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index ad80a1c16e..ae974b5a5a 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -43,7 +43,7 @@ class Action_HydrogenBond : public Action { ImagedAction Image_; ///< Hold imaging info. Sarray Both_; ///< Array of donor sites that can also be acceptors Iarray Acceptor_; ///< Array of acceptor-only atom indices - Sarray SolventSites_; + Sarray SolventSites_; ///< Array of solvent donor/acceptor sites UUmapType UU_Map_; UVmapType UV_Map_; From 902943e4199042501bddcc720aef4cbc7640b323 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 08:54:20 -0400 Subject: [PATCH 12/56] DRR - Cpptraj: Recognize ions --- src/Action_HydrogenBond.cpp | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index ae10ab143e..3d712cf15d 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -320,7 +320,8 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { for (int at = 0; at != maxAtom; at++) { // Since an acceptor mask was not specified ignore solvent. - int molnum = setup.Top()[at].MolNum(); + Atom const& atom = setup.Top()[at]; + int molnum = atom.MolNum(); if (!setup.Top().Mol(molnum).IsSolvent()) { Iarray Hatoms; @@ -328,9 +329,9 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { isAcceptor = false; if ( d_atom != Mask_.end() && *d_atom == at) { ++d_atom; - if ( IsFON( setup.Top()[at] ) ) { - for (Atom::bond_iterator H_at = setup.Top()[at].bondbegin(); - H_at != setup.Top()[at].bondend(); ++H_at) + if ( IsFON( atom ) ) { + for (Atom::bond_iterator H_at = atom.bondbegin(); + H_at != atom.bondend(); ++H_at) if (setup.Top()[*H_at].Element() == Atom::HYDROGEN) Hatoms.push_back( *H_at ); isDonor = !Hatoms.empty(); @@ -411,11 +412,16 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { isAcceptor = false; if ( d_atom != SolventDonorMask_.end() && *d_atom == at ) { ++d_atom; - if ( IsFON( setup.Top()[at] ) ) { - for (Atom::bond_iterator H_at = setup.Top()[at].bondbegin(); - H_at != setup.Top()[at].bondend(); ++H_at) + Atom const& atom = setup.Top()[at]; + if ( IsFON( atom ) ) { + for (Atom::bond_iterator H_at = atom.bondbegin(); + H_at != atom.bondend(); ++H_at) if (setup.Top()[*H_at].Element() == Atom::HYDROGEN) Hatoms.push_back( *H_at ); + } else if ( atom.Nbonds() == 0 ) { + // If no bonds to this atom assume it is an ion. Set the H atom + // to be the same as D atom; this will skip the angle calc. + Hatoms.push_back( at ); } isDonor = !Hatoms.empty(); } From 2248a30821f91d372d53868cd3cb1afaa13a7dee Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 09:06:53 -0400 Subject: [PATCH 13/56] DRR - Cpptraj: Add solute/solvent hydrogen and ion count --- src/Action_HydrogenBond.cpp | 19 +++++++++++++++---- src/Action_HydrogenBond.h | 4 ++++ 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 3d712cf15d..edf52bffb6 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -355,24 +355,28 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { for (Sarray::const_iterator site = donorOnly.begin(); site != donorOnly.end(); ++site) Both_.push_back( *site ); - mprintf("Acceptor atoms (%zu):\n", Acceptor_.size()); + mprintf(" Acceptor atoms (%zu):\n", Acceptor_.size()); for (Iarray::const_iterator at = Acceptor_.begin(); at != Acceptor_.end(); ++at) mprintf("\t%20s %8i\n", setup.Top().TruncResAtomName(*at).c_str(), *at+1); - mprintf("Donor/acceptor sites (%u):\n", bothEnd_); + unsigned int hcount = 0; + mprintf(" Donor/acceptor sites (%u):\n", bothEnd_); Sarray::const_iterator END = Both_.begin() + bothEnd_; for (Sarray::const_iterator si = Both_.begin(); si != END; ++si) { + hcount += si->n_hydrogens(); mprintf("\t%20s %8i", setup.Top().TruncResAtomName(si->Idx()).c_str(), si->Idx()+1); for (Iarray::const_iterator at = si->Hbegin(); at != si->Hend(); ++at) mprintf(" %s", setup.Top()[*at].c_str()); mprintf("\n"); } - mprintf("Donor sites (%zu):\n", Both_.size() - bothEnd_); + mprintf(" Donor sites (%zu):\n", Both_.size() - bothEnd_); for (Sarray::const_iterator si = END; si != Both_.end(); ++si) { + hcount += si->n_hydrogens(); mprintf("\t%20s %8i", setup.Top().TruncResAtomName(si->Idx()).c_str(), si->Idx()+1); for (Iarray::const_iterator at = si->Hbegin(); at != si->Hend(); ++at) mprintf(" %s", setup.Top()[*at].c_str()); mprintf("\n"); } + mprintf(" %u solute hydrogens.\n", hcount); // SOLVENT SITE SETUP if (calcSolvent_) { @@ -433,13 +437,20 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { SolventSites_.push_back( Site(at, Hatoms) ); } - mprintf("Solvent sites (%zu):\n", SolventSites_.size()); + hcount = 0; + unsigned int icount = 0; + mprintf(" Solvent sites (%zu):\n", SolventSites_.size()); for (Sarray::const_iterator si = SolventSites_.begin(); si != SolventSites_.end(); ++si) { + if (si->IsIon()) + icount++; + else + hcount += si->n_hydrogens(); mprintf("\t%20s %8i", setup.Top().TruncResAtomName(si->Idx()).c_str(), si->Idx()+1); for (Iarray::const_iterator at = si->Hbegin(); at != si->Hend(); ++at) mprintf(" %s", setup.Top()[*at].c_str()); mprintf("\n"); } + mprintf(" %u solvent hydrogens, %u ions.\n", hcount, icount); } return Action::OK; diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index ae974b5a5a..6f216f2816 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -99,6 +99,10 @@ class Action_HydrogenBond::Site { Site(int d, Iarray const& H) : hlist_(H), idx_(d), isV_(false) {} /// \return heavy atom index int Idx() const { return idx_; } + /// \return number of hydrogen indices + unsigned int n_hydrogens() const { return hlist_.size(); } + /// \return true if site is an ion (D atom == H atom) + bool IsIon() const { return (hlist_.size()==1 && hlist_[0] == idx_); } /// \return iterator to beginning of hydrogen indices Iarray::const_iterator Hbegin() const { return hlist_.begin(); } /// \return iterator to end of hydrogen indices From 86e1b0f6aaffbd01c59c20519d28033537b91c15 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 09:33:21 -0400 Subject: [PATCH 14/56] DRR - Cpptraj: Add bridging calc --- src/Action_HydrogenBond.cpp | 83 ++++++++++++++++++++++++++++++++------------- src/Action_HydrogenBond.h | 24 +++++++++++-- 2 files changed, 81 insertions(+), 26 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index edf52bffb6..0cfef7b13a 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -503,13 +503,20 @@ void Action_HydrogenBond::CalcSolvHbonds(int frameNum, double dist2, { ++numHB; double dist = sqrt(dist2); - // Index U-H .. V hydrogen bonds by solute H atom. - // Index U .. H-V hydrogen bonds by solute A atom. - int hbidx; - if (soluteDonor) + int hbidx, solventmol, soluteres; + // TODO: Option to use solvent mol num? + if (soluteDonor) { + // Index U-H .. V hydrogen bonds by solute H atom. hbidx = *h_atom; - else + soluteres = (*CurrentParm_)[d_atom].ResNum(); + solventmol = (*CurrentParm_)[a_atom].ResNum(); + } else { + // Index U .. H-V hydrogen bonds by solute A atom. hbidx = a_atom; + soluteres = (*CurrentParm_)[a_atom].ResNum(); + solventmol = (*CurrentParm_)[d_atom].ResNum(); + } + solvent2solute_[solventmol].insert( soluteres ); UVmapType::iterator it = UV_Map_.lower_bound( hbidx ); if (it == UV_Map_.end() || it->first != hbidx) { @@ -620,6 +627,7 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { // Loop over all solvent sites if (calcSolvent_) { t_uv_.Start(); + solvent2solute_.clear(); numHB = 0; for (unsigned int vidx = 0; vidx != SolventSites_.size(); vidx++) { @@ -659,6 +667,40 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { } // END loop over solvent sites NumSolvent_->Add(frameNum, &numHB); t_uv_.Stop(); + + // Determine number of bridging waters + t_bridge_.Start(); + numHB = 0; + std::string bridgeID; + for (RmapType::const_iterator bridge = solvent2solute_.begin(); + bridge != solvent2solute_.end(); ++bridge) + { + // If solvent molecule is bound to 2 or more different residues, + // it is bridging. + if ( bridge->second.size() > 1) { + ++numHB; + // Bridging Solvent residue number + bridgeID.append(integerToString( bridge->first+1 ) + "("); + 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) + 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,int>(bridge->second, 1) ); + else + // Increment bridge #frames + b_it->second++; + } + } + if (bridgeID.empty()) + bridgeID.assign("None"); + NumBridge_->Add(frameNum, &numHB); + BridgeID_->Add(frameNum, bridgeID.c_str()); + t_bridge_.Stop(); } Nframes_++; @@ -835,7 +877,7 @@ int Action_HydrogenBond::SyncAction() { unsigned int i2 = idx + 1; for (int ir = 0; ir != iArray[idx]; ir++, i2++) residues.insert( iArray[i2] ); - BridgeType::iterator b_it = BridgeMap_.find( residues ); + BmapType::iterator b_it = BridgeMap_.find( residues ); if (b_it == BridgeMap_.end() ) // New Bridge BridgeMap_.insert( b_it, std::pair,int>(residues, iArray[i2]) ); else // Increment bridge #frames @@ -845,7 +887,7 @@ int Action_HydrogenBond::SyncAction() { } } else { // Construct bridge info array. - for (BridgeType::const_iterator b = BridgeMap_.begin(); b != BridgeMap_.end(); ++b) + for (BmapType::const_iterator b = BridgeMap_.begin(); b != BridgeMap_.end(); ++b) { iArray.push_back( b->first.size() ); // # of bridging res for ( std::set::const_iterator r = b->first.begin(); r != b->first.end(); ++r) @@ -878,22 +920,20 @@ void Action_HydrogenBond::UpdateSeries() { // Action_Hbond::MemoryUsage() std::string Action_HydrogenBond::MemoryUsage(size_t nPairs, size_t nFrames) const { - static const size_t HBmapTypeElt = 32 + sizeof(int) + - (2*sizeof(double) + sizeof(DataSet_integer*) + 4*sizeof(int)); - size_t memTotal = nPairs * HBmapTypeElt; + static const size_t sizeUUmap = 32 + sizeof(int) + + (2*sizeof(double) + sizeof(DataSet_integer*) + 4*sizeof(int)); + size_t memTotal = nPairs * sizeUUmap; // If calculating series every hbond will have time series. // NOTE: This does not include memory used by DataSet. if (series_ && nFrames > 0) { size_t seriesSet = (nFrames * sizeof(int)) + sizeof(std::vector); memTotal += (seriesSet * nPairs); } -/* // Current memory used by bridging solvent - static const size_t BridgeTypeElt = 32 + sizeof(std::set) + sizeof(int); - for (BridgeType::const_iterator it = BridgeMap_.begin(); it != BridgeMap_.end(); ++it) + static const size_t sizeBmap = 32 + sizeof(std::set) + sizeof(int); + for (BmapType::const_iterator it = BridgeMap_.begin(); it != BridgeMap_.end(); ++it) memTotal += (it->first.size() * sizeof(int)); - memTotal += (BridgeMap_.size() * BridgeTypeElt); -*/ + memTotal += (BridgeMap_.size() * sizeBmap); return ByteString( memTotal, BYTE_DECIMAL ); } @@ -1001,19 +1041,17 @@ void Action_HydrogenBond::Print() { hbond->Frames(), avg, hbond->Dist(), hbond->Angle()); } } -/* // BRIDGING INFO if (bridgeout_ != 0 && calcSolvent_) { bridgeout_->Printf("#Bridging Solute Residues:\n"); // Place bridging values in a vector for sorting - std::vector, int> > bridgevector; - for (BridgeType::const_iterator it = BridgeMap_.begin(); - it != BridgeMap_.end(); ++it) + typedef std::vector, int> > Bvec; + Bvec bridgevector; + for (BmapType::const_iterator it = BridgeMap_.begin(); + it != BridgeMap_.end(); ++it) bridgevector.push_back( *it ); std::sort( bridgevector.begin(), bridgevector.end(), bridge_cmp() ); - for (std::vector, int> >::const_iterator bv = bridgevector.begin(); - bv != bridgevector.end(); - ++bv) + for (Bvec::const_iterator bv = bridgevector.begin(); bv != bridgevector.end(); ++bv) { bridgeout_->Printf("Bridge Res"); for (std::set::const_iterator res = bv->first.begin(); @@ -1022,7 +1060,6 @@ void Action_HydrogenBond::Print() { bridgeout_->Printf(", %i frames.\n", bv->second); } } -*/ } // ===== Action_HydrogenBond::Hbond ============================================ diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index 6f216f2816..e95df911ce 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -39,14 +39,18 @@ class Action_HydrogenBond : public Action { typedef std::pair Hpair; typedef std::map UUmapType; typedef std::map UVmapType; + typedef std::map< int,std::set > RmapType; + typedef std::map< std::set,int > BmapType; - ImagedAction Image_; ///< Hold imaging info. - Sarray Both_; ///< Array of donor sites that can also be acceptors - Iarray Acceptor_; ///< Array of acceptor-only atom indices + ImagedAction Image_; ///< Hold imaging info. + Sarray Both_; ///< Array of donor sites that can also be acceptors + Iarray Acceptor_; ///< Array of acceptor-only atom indices Sarray SolventSites_; ///< Array of solvent donor/acceptor sites UUmapType UU_Map_; UVmapType UV_Map_; + RmapType solvent2solute_; ///< Map solvent mol # to residues it is bound to each frame + BmapType BridgeMap_; ///< Map residues involved in bridging to # frames bridge present std::string hbsetname_; AtomMask DonorMask_; @@ -86,6 +90,20 @@ class Action_HydrogenBond : public Action { bool hasSolventDonor_; bool calcSolvent_; bool hasSolventAcceptor_; + // TODO replace with class + /// Return true if first bridge has more frames than second. + struct bridge_cmp { + inline bool operator()(std::pair< std::set, int> const& first, + std::pair< std::set, int> const& second) const + { + if (first.second > second.second) + return true; + else if (first.second < second.second) + return false; + else + return (first.second < second.second); + } + }; }; // ----- CLASSES --------------------------------------------------------------- From 3932a7a3a17b94553d84d6914ee4e7897650305c Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 09:44:56 -0400 Subject: [PATCH 15/56] DRR - Cpptraj: Properly init some variables. Fix spacing. --- 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 0cfef7b13a..64c4a184ef 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -21,9 +21,11 @@ Action_HydrogenBond::Action_HydrogenBond() : bridgeout_(0), dcut2_(0.0), acut_(0.0), + bothEnd_(0), Nframes_(0), debug_(0), - series_(0), + series_(false), + seriesUpdated_(false), useAtomNum_(false), noIntramol_(false), hasDonorMask_(false), @@ -953,9 +955,9 @@ void Action_HydrogenBond::Print() { // mprintf("\t%zu unique solute-solvent bridging interactions.\n", BridgeMap_.size()); } - t_uu_.WriteTiming( 2, "Solute-Solute :",t_action_.Total()); + t_uu_.WriteTiming( 2,"Solute-Solute :",t_action_.Total()); if (calcSolvent_) { - t_uv_.WriteTiming( 2, "Solute-Solvent :",t_uv_.Total()); + t_uv_.WriteTiming( 2,"Solute-Solvent :",t_uv_.Total()); t_bridge_.WriteTiming(2,"Bridging waters :",t_action_.Total()); } t_action_.WriteTiming(1,"Total:"); From cde9c408719c436d82a90a0252fa2c72f50dc9ec Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 09:51:58 -0400 Subject: [PATCH 16/56] DRR - Cpptraj: Remove unused class variable. --- src/Action_HydrogenBond.h | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index e95df911ce..b0fe69e815 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -110,11 +110,11 @@ class Action_HydrogenBond : public Action { /// Potential hydrogen bond site. Can be either donor or donor/acceptor. class Action_HydrogenBond::Site { public: - Site() : idx_(-1), isV_(false) {} + Site() : idx_(-1) {} /// Solute site - heavy atom, hydrogen atom - Site(int d, int h) : hlist_(1,h), idx_(d), isV_(false) {} + Site(int d, int h) : hlist_(1,h), idx_(d) {} /// Solute site - heavy atom, list of hydrogen atoms - Site(int d, Iarray const& H) : hlist_(H), idx_(d), isV_(false) {} + Site(int d, Iarray const& H) : hlist_(H), idx_(d) {} /// \return heavy atom index int Idx() const { return idx_; } /// \return number of hydrogen indices @@ -128,7 +128,6 @@ class Action_HydrogenBond::Site { private: Iarray hlist_; ///< List of hydrogen indices int idx_; ///< Heavy atom index - bool isV_; ///< True if site is solvent }; /// Track specific hydrogen bond. From 34c36b27d206562cf8c97bed95ad237cf6c04fae Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 10:04:24 -0400 Subject: [PATCH 17/56] DRR - Cpptraj: Add nointramol functionality. Slows down the code a little since it introduces an if statement in some inner loops but is still faster than original. --- src/Action_HydrogenBond.cpp | 43 ++++++++++++++++++++++++++----------------- 1 file changed, 26 insertions(+), 17 deletions(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 64c4a184ef..0385441da3 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -357,11 +357,11 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { for (Sarray::const_iterator site = donorOnly.begin(); site != donorOnly.end(); ++site) Both_.push_back( *site ); - mprintf(" Acceptor atoms (%zu):\n", Acceptor_.size()); + mprintf(" Acceptor-only atoms (%zu)\n", Acceptor_.size()); for (Iarray::const_iterator at = Acceptor_.begin(); at != Acceptor_.end(); ++at) mprintf("\t%20s %8i\n", setup.Top().TruncResAtomName(*at).c_str(), *at+1); unsigned int hcount = 0; - mprintf(" Donor/acceptor sites (%u):\n", bothEnd_); + mprintf(" Donor/acceptor sites (%u)\n", bothEnd_); Sarray::const_iterator END = Both_.begin() + bothEnd_; for (Sarray::const_iterator si = Both_.begin(); si != END; ++si) { hcount += si->n_hydrogens(); @@ -370,7 +370,7 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { mprintf(" %s", setup.Top()[*at].c_str()); mprintf("\n"); } - mprintf(" Donor sites (%zu):\n", Both_.size() - bothEnd_); + mprintf(" Donor-only sites (%zu)\n", Both_.size() - bothEnd_); for (Sarray::const_iterator si = END; si != Both_.end(); ++si) { hcount += si->n_hydrogens(); mprintf("\t%20s %8i", setup.Top().TruncResAtomName(si->Idx()).c_str(), si->Idx()+1); @@ -441,7 +441,7 @@ Action::RetType Action_HydrogenBond::Setup(ActionSetup& setup) { hcount = 0; unsigned int icount = 0; - mprintf(" Solvent sites (%zu):\n", SolventSites_.size()); + mprintf(" Solvent sites (%zu)\n", SolventSites_.size()); for (Sarray::const_iterator si = SolventSites_.begin(); si != SolventSites_.end(); ++si) { if (si->IsIon()) icount++; @@ -598,29 +598,38 @@ Action::RetType Action_HydrogenBond::DoAction(int frameNum, ActionFrame& frm) { // Loop over all solute donor sites t_uu_.Start(); int numHB = 0; + int mol0 = -1; for (unsigned int sidx0 = 0; sidx0 != Both_.size(); sidx0++) { - const double* XYZ0 = frm.Frm().XYZ( Both_[sidx0].Idx() ); + Site const& Site0 = Both_[sidx0]; + const double* XYZ0 = frm.Frm().XYZ( Site0.Idx() ); + if (noIntramol_) + mol0 = (*CurrentParm_)[Site0.Idx()].MolNum(); // Loop over solute sites that can be both donor and acceptor for (unsigned int sidx1 = sidx0 + 1; sidx1 < bothEnd_; sidx1++) { - const double* XYZ1 = frm.Frm().XYZ( Both_[sidx1].Idx() ); - double dist2 = DIST2( XYZ0, XYZ1, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); - if ( !(dist2 > dcut2_) ) - { - // Site 0 donor, Site 1 acceptor - CalcSiteHbonds(frameNum, dist2, Both_[sidx0], XYZ0, Both_[sidx1].Idx(), XYZ1, frm.Frm(), numHB); - // Site 1 donor, Site 0 acceptor - CalcSiteHbonds(frameNum, dist2, Both_[sidx1], XYZ1, Both_[sidx0].Idx(), XYZ0, frm.Frm(), numHB); + Site const& Site1 = Both_[sidx1]; + if (mol0 != (*CurrentParm_)[Site1.Idx()].MolNum()) { + const double* XYZ1 = frm.Frm().XYZ( Site1.Idx() ); + double dist2 = DIST2( XYZ0, XYZ1, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); + if ( !(dist2 > dcut2_) ) + { + // Site 0 donor, Site 1 acceptor + CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, Site1.Idx(), XYZ1, frm.Frm(), numHB); + // Site 1 donor, Site 0 acceptor + CalcSiteHbonds(frameNum, dist2, Site1, XYZ1, Site0.Idx(), XYZ0, frm.Frm(), numHB); + } } } // Loop over solute acceptor-only for (Iarray::const_iterator a_atom = Acceptor_.begin(); a_atom != Acceptor_.end(); ++a_atom) { - const double* XYZ1 = frm.Frm().XYZ( *a_atom ); - double dist2 = DIST2( XYZ0, XYZ1, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); - if ( !(dist2 > dcut2_) ) - CalcSiteHbonds(frameNum, dist2, Both_[sidx0], XYZ0, *a_atom, XYZ1, frm.Frm(), numHB); + if (mol0 != (*CurrentParm_)[*a_atom].MolNum()) { + const double* XYZ1 = frm.Frm().XYZ( *a_atom ); + double dist2 = DIST2( XYZ0, XYZ1, Image_.ImageType(), frm.Frm().BoxCrd(), ucell_, recip_ ); + if ( !(dist2 > dcut2_) ) + CalcSiteHbonds(frameNum, dist2, Site0, XYZ0, *a_atom, XYZ1, frm.Frm(), numHB); + } } } NumHbonds_->Add(frameNum, &numHB); From 256bb2f7442277c6d81dd54d369a3f8d2938fefe Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 11:36:36 -0400 Subject: [PATCH 18/56] DRR - Cpptraj: Start fixing the MPI consolidation --- src/Action_HydrogenBond.cpp | 93 ++++++++++++++++++++++++++++++++++++++++++++- src/Action_HydrogenBond.h | 15 ++++++++ 2 files changed, 107 insertions(+), 1 deletion(-) diff --git a/src/Action_HydrogenBond.cpp b/src/Action_HydrogenBond.cpp index 0385441da3..5ac0bee39a 100644 --- a/src/Action_HydrogenBond.cpp +++ b/src/Action_HydrogenBond.cpp @@ -732,7 +732,7 @@ static inline std::string CreateHBlegend(Topology const& topIn, int a_atom, int topIn.TruncResAtomName(d_atom) + "-" + topIn[h_atom].Name().Truncated()); } - +/* // Action_HydrogenBond::SyncMap void Action_HydrogenBond::SyncMap(UUmapType& mapIn, std::vector const& rank_frames, std::vector const& rank_offsets, @@ -841,6 +841,17 @@ const } } } +*/ + +// Action_HydrogenBond::GetRankNhbonds() +std::vector Action_HydrogenBond::GetRankNhbonds( int num_hb, Parallel::Comm const& commIn ) +{ + std::vector nhb_on_rank; + if (commIn.Master()) + nhb_on_rank.resize( commIn.Size() ); + commIn.GatherMaster( &num_hb, 1, MPI_INT, &nhb_on_rank[0] ); + return nhb_on_rank; +} /** PARALLEL NOTES: * The following tags are used for MPI send/receive: @@ -867,10 +878,90 @@ int Action_HydrogenBond::SyncAction() { for (int rank = 1; rank < trajComm_.Size(); rank++) rank_offsets[rank] = rank_offsets[rank-1] + rank_frames[rank-1]; } + // Need to send hbond data from all ranks to master. + std::vector Dvals; // Hold dist_ and angle_ for each hbond + std::vector Ivals; // Hold A_, H_, D_, and frames_ for each hbond + // Need to know how many hbonds on each thread. + std::vector nhb_on_rank = GetRankNhbonds( UU_Map_.size(), trajComm_ ); + if (trajComm_.Master()) { + for (int rank = 1; rank < trajComm_.Size(); rank++) + { + if (nhb_on_rank[rank] > 0) + { + Dvals.resize( 2 * nhb_on_rank[rank] ); + Ivals.resize( 4 * nhb_on_rank[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 + // Dvals = dist, angle + // Avals = A, H, D, frames + const int* IV = &Ivals[0]; + const double* DV = &Dvals[0]; + for (int in = 0; in != nhb_on_rank[rank]; in++, IV += 4, DV += 2) + { + Hpair hbidx(IV[1], IV[0]); + UUmapType::iterator it = UU_Map_.lower_bound( hbidx ); + DataSet_integer* ds = 0; + if (it == UU_Map_.end() || it->first != hbidx) + { + // Hbond on rank that has not been found on master + if (series_) { + ds = (DataSet_integer*) + masterDSL_->AddSet(DataSet::INTEGER, MetaData(hbsetname_,"solutehb",UU_Map_.size())); + ds->SetLegend( CreateHBlegend(*CurrentParm_, 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]))); + } else { + // Hbond on rank and master. Update on master. + it->second.Combine(DV[0], DV[1], IV[3]); + ds = it->second.Data(); + } + // Update time series + if (series_) { + ds->Resize( Nframes_ ); + int* d_beg = ds->Ptr() + rank_offsets[ rank ]; + //mprintf("\tResizing hbond series data to %i, starting frame %i, # frames %i\n", + // Nframes_, rank_offsets[rank], rank_frames[rank]); + trajComm_.Recv( d_beg, rank_frames[ rank ], MPI_INT, rank, 1304 + in ); + ds->SetNeedsSync( false ); + } + } + } + } + } else { + if (!UU_Map_.empty()) { + // Store UU bonds in flat arrays. + Dvals.reserve( UU_Map_.size() * 2 ); + Ivals.reserve( UU_Map_.size() * 4 ); + 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() ); + } + trajComm_.Send( &(Dvals[0]), Dvals.size(), MPI_DOUBLE, 0, 1300 ); + trajComm_.Send( &(Ivals[0]), Ivals.size(), MPI_INT, 0, 1301 ); + // Send series data to master + if (series_) { + int in = 0; // For tag + for (UUmapType::const_iterator hb = UU_Map_.begin(); hb != UU_Map_.end(); ++hb, in++) { + trajComm_.Send( hb->second.Data()->Ptr(), hb->second.Data()->Size(), MPI_INT, 0, 1304 + in ); + hb->second.Data()->SetNeedsSync( false ); + } + } + } + } + +/* SyncMap( UU_Map_, rank_frames, rank_offsets, "solutehb", trajComm_ ); +*/ if (calcSolvent_) { +/* SyncMap( UV_Map_, rank_frames, rank_offsets, "solventhb", trajComm_ ); +*/ // iArray will contain for each bridge: Nres, res1, ..., resN, Frames std::vector iArray; int iSize; diff --git a/src/Action_HydrogenBond.h b/src/Action_HydrogenBond.h index b0fe69e815..973d09fc66 100644 --- a/src/Action_HydrogenBond.h +++ b/src/Action_HydrogenBond.h @@ -34,6 +34,9 @@ class Action_HydrogenBond : public Action { void UpdateSeries(); /// Determine memory usage from # hbonds and time series std::string MemoryUsage(size_t, size_t) const; +# ifdef MPI + static std::vector GetRankNhbonds(int,Parallel::Comm const&); +# endif typedef std::vector Sarray; typedef std::pair Hpair; @@ -143,6 +146,18 @@ class Action_HydrogenBond::Hbond { int A() const { return A_; } int H() const { return H_; } int D() const { return D_; } +# ifdef MPI + DataSet_integer* Data() const { return data_; } + /// 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 + void Combine(double d, double a, int n) { + dist_ += d; + angle_ += a; + frames_ += n; + } +# endif /// First sort by frames (descending), then distance (ascending). bool operator<(const Hbond& rhs) const { if (frames_ == rhs.frames_) From 75cf54688865d8f0418e6a8e4dc1c4cbf1a84554 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Thu, 13 Apr 2017 12:59:29 -0400 Subject: [PATCH 19/56] DRR - Cpptraj: Add sort option for data file write. This allows comparison between series data of the old and new versions of the hbond code since the order in which hbonds are searched for is slightly different --- src/DataFile.cpp | 5 ++++- src/DataFile.h | 1 + src/DataSetList.cpp | 4 ++++ src/DataSetList.h | 2 ++ 4 files changed, 11 insertions(+), 1 deletion(-) diff --git a/src/DataFile.cpp b/src/DataFile.cpp index 065ca7ab8a..63f04bdf55 100644 --- a/src/DataFile.cpp +++ b/src/DataFile.cpp @@ -29,6 +29,7 @@ DataFile::DataFile() : dfType_(DATAFILE), dflWrite_(true), setDataSetPrecision_(false), //TODO: Just use default_width_ > -1? + sortSets_(false), default_width_(-1), default_precision_(-1), dataio_(0), @@ -83,7 +84,7 @@ const FileTypes::KeyToken DataFile::DF_KeyArray[] = { void DataFile::WriteHelp() { mprintf("\t[]\n" - "\t[{xlabel|ylabel|zlabel}