diff --git a/README.md b/README.md index 9605ff09d8..864166cf91 100644 --- a/README.md +++ b/README.md @@ -180,9 +180,6 @@ Original code for the 'xtalsymm' Action. Alrun N. Koller (Heinrich-Heine-University, Düsseldorf, Germany) Original implementation of matrix/vector functionality in PTRAJ, including matrix diagonalization, IRED analysis, eigenmode analysis, and vector time correlations. -* Holger Gohlke (Heinrich-Heine-University, Düsseldorf, Germany) -Original code for DSSP (secstruct). - * Michael Crowley (University of Southern California, Los Angeles, CA, USA) Original code for dealing with truncated octahedral unit cells. diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf index 25bd9be940..f4a3bc090c 100644 Binary files a/doc/CpptrajManual.pdf and b/doc/CpptrajManual.pdf differ diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 35887a3aff..18cc1c59e1 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -30504,6 +30504,10 @@ secstruct [assignout ] [totalout [ptrajformat] \end_layout +\begin_layout LyX-Code + [betadetail] +\end_layout + \begin_layout LyX-Code [namen ] [nameh ] [nameca ] \end_layout @@ -30556,6 +30560,12 @@ secstruct frame, similar to ptraj output. \end_layout +\begin_layout Description +[betadetail] Record anti-parallel beta and parallel beta in place of extended + and bridge secondary structure. + If a residue could be both only anti-parallel is reported. +\end_layout + \begin_layout Description [namen \begin_inset space ~ @@ -30680,12 +30690,20 @@ None \end_deeper \begin_layout Standard -\shape italic -Note that when not using \series bold -ptrajformat -\series default -, data sets are not generated until +As of version 4.18.0, this command now produces output that better conforms + with the original definitions in Kabsch and Sander 1983; namely that Extended + beta (i.e. + 2 or more consecutive beta bridges of the same type) and beta Bridge (i.e. + an isolated beta bridge) are now reported instead of anti-parallel and + parallel beta. + To restore the original behavior the 'betadetail' keyword must be specified. +\end_layout + +\begin_layout Standard + +\shape italic +Note that the residue and [avgss] data sets are not generated until \series bold run \series default @@ -30802,9 +30820,11 @@ where STRING is a string of characters (one for each residue) where each \shape italic ptraj \shape default - outputs). + had outputed and is retained for backwards compatibility). The various secondary structure types and their corresponding integer/character - are listed below: + are listed below. + If 'betadetail' is specified what is reported and the characters used change + slightly. \begin_inset Separator latexpar \end_inset @@ -30825,7 +30845,7 @@ ptraj \begin_inset Text \begin_layout Plain Layout -Character +STRING (betadetail) \end_layout \end_inset @@ -30843,7 +30863,7 @@ Integer \begin_inset Text \begin_layout Plain Layout -DSSP Char +DSSP \end_layout \end_inset @@ -30852,7 +30872,7 @@ DSSP Char \begin_inset Text \begin_layout Plain Layout -SS type +SS type (betadetail) \end_layout \end_inset @@ -30901,7 +30921,7 @@ None \begin_inset Text \begin_layout Plain Layout -b +E (b) \end_layout \end_inset @@ -30928,7 +30948,7 @@ b \begin_inset Text \begin_layout Plain Layout -Parallel Beta-sheet +Extended beta (parallel beta) \end_layout \end_inset @@ -30966,7 +30986,7 @@ B \begin_inset Text \begin_layout Plain Layout -Anti-parallel Beta-sheet +Isolated beta (anti-parallel beta) \end_layout \end_inset @@ -31210,7 +31230,7 @@ totalout \end_layout \begin_layout LyX-Code - SSS TH HHHTTSBBBB TTTBBBB SS S + SSS TH HHHTTSEEEE TTTEEEE SS S \end_layout \begin_layout Standard diff --git a/src/Action_DSSP.cpp b/src/Action_DSSP.cpp index 25fdd7e161..9075abcc06 100644 --- a/src/Action_DSSP.cpp +++ b/src/Action_DSSP.cpp @@ -1,47 +1,302 @@ #include // sqrt +#include // std::fill, std::min, std::max +#include // SET #include "Action_DSSP.h" #include "CpptrajStdio.h" +#include "DataSet.h" #include "DistRoutines.h" -/// Hbond energy calc prefactor for kcal/mol: q1*q2*E, 0.42*0.20*332 -const double Action_DSSP::DSSP_fac = 27.888; +#ifdef DSSPDEBUG +#include "Constants.h" +#endif +#ifdef _OPENMP +# include +#endif + +/** From the original Kabsch & Sander 1983 paper. Obtained via + * q1 = 0.42e + * q2 = 0.20e + * f = 332 (approximate conversion factor to get energies in kcal/mol) + * fac = q1*q2*f + */ +const double Action_DSSP::DSSP_fac_ = 27.888; + +/** DSSP default cutoff for determining hbonds, from 1983 Kabsch & Sander paper. */ +const double Action_DSSP::DSSP_cut_ = -0.5; + +/** DSSP 1 character SS assignment. */ +const char Action_DSSP::DSSP_char_[] = { ' ', 'E', 'B', 'G', 'H', 'I', 'T', 'S' }; + +/** Used for output to STRING data set. */ +const char* Action_DSSP::SSchar_[] = { "0", "E", "B", "G", "H", "I", "T", "S" }; + +/** Full SS names. */ +const char* Action_DSSP::DSSP_name_[]={"None", "Extended", "Bridge", "3-10", "Alpha", "Pi", "Turn", "Bend"}; + + +// ----- SSres ----------------------------------------------------------------- +Action_DSSP::SSres::SSres() : + resDataSet_(0), + chirality_(0), + bend_(0), + sstype_(NONE), + num_(-1), + C_(-1), + O_(-1), + N_(-1), + H_(-1), + CA_(-1), + prevIdx_(-1), + nextIdx_(-1), + bridge1idx_(-1), + b1type_(NO_BRIDGE), + bridge2idx_(-1), + b2type_(NO_BRIDGE), + resChar_(' '), + isSelected_(false) +{ + std::fill(SScount_, SScount_ + NSSTYPE_, 0); + std::fill(Bcount_, Bcount_ + NBRIDGETYPE_, 0); + std::fill(turnChar_, turnChar_ + NTURNTYPE_, ' '); +} + +Action_DSSP::SSres::SSres(SSres const& rhs) : + resDataSet_(rhs.resDataSet_), + chirality_(rhs.chirality_), + bend_(rhs.bend_), + sstype_(rhs.sstype_), + num_(rhs.num_), + C_(rhs.C_), + O_(rhs.O_), + N_(rhs.N_), + H_(rhs.H_), + CA_(rhs.CA_), + prevIdx_(rhs.prevIdx_), + nextIdx_(rhs.nextIdx_), + bridge1idx_(rhs.bridge1idx_), + b1type_(rhs.b1type_), + bridge2idx_(rhs.bridge2idx_), + b2type_(rhs.b2type_), + resChar_(rhs.resChar_), + isSelected_(rhs.isSelected_) +{ + std::copy(rhs.SScount_, rhs.SScount_ + NSSTYPE_, SScount_); + std::copy(rhs.Bcount_, rhs.Bcount_ + NBRIDGETYPE_, Bcount_); + std::copy(rhs.turnChar_, rhs.turnChar_ + NTURNTYPE_, turnChar_); +} + +Action_DSSP::SSres& Action_DSSP::SSres::operator=(SSres const& rhs) { + if (this == &rhs) return *this; + resDataSet_ = rhs.resDataSet_; + chirality_ = rhs.chirality_; + bend_ = rhs.bend_; + sstype_ = rhs.sstype_; + num_ = rhs.num_; + C_ = rhs.C_; + O_ = rhs.O_; + N_ = rhs.N_; + H_ = rhs.H_; + CA_ = rhs.CA_; + prevIdx_ = rhs.prevIdx_; + nextIdx_ = rhs.nextIdx_; + bridge1idx_ = rhs.bridge1idx_; + b1type_ = rhs.b1type_; + bridge2idx_ = rhs.bridge2idx_; + b2type_ = rhs.b2type_; + resChar_ = rhs.resChar_; + isSelected_ = rhs.isSelected_; + std::copy(rhs.SScount_, rhs.SScount_ + NSSTYPE_, SScount_); + std::copy(rhs.Bcount_, rhs.Bcount_ + NBRIDGETYPE_, Bcount_); + std::copy(rhs.turnChar_, rhs.turnChar_ + NTURNTYPE_, turnChar_); + return *this; +} + + +void Action_DSSP::SSres::Deselect() { + isSelected_ = false; + C_ = -1; + H_ = -1; + N_ = -1; + O_ = -1; + CA_ = -1; +} + +void Action_DSSP::SSres::Unassign() { + sstype_ = NONE; + bridge1idx_ = -1; + b1type_ = NO_BRIDGE; + bridge2idx_ = -1; + b2type_ = NO_BRIDGE; + std::fill(turnChar_, turnChar_ + NTURNTYPE_, ' '); +} + +/** Accumulate SS data. */ +void Action_DSSP::SSres::AccumulateData(int frameNum, bool useString, bool betaDetail) +{ + SScount_[sstype_]++; + int idata = (int)sstype_; + const char* sdata = SSchar_[sstype_]; + if (sstype_ == EXTENDED || sstype_ == BRIDGE) { + if (b1type_ == PARALLEL || b2type_ == PARALLEL) { + Bcount_[PARALLEL]++; + if (betaDetail) { + idata = (int)EXTENDED; + sdata = "b"; + } + } + if (b1type_ == ANTIPARALLEL || b2type_ == ANTIPARALLEL) { + Bcount_[ANTIPARALLEL]++; + if (betaDetail) { + idata = (int)BRIDGE; + sdata = "B"; + } + } + // TODO bulge? + } else + Bcount_[NO_BRIDGE]++; + if (useString) + resDataSet_->Add(frameNum, sdata); + else + resDataSet_->Add(frameNum, &idata); +} + +/** Set turn beginning. */ +void Action_DSSP::SSres::SetTurnBegin(TurnType typeIn) { + if (turnChar_[typeIn] == '<') + turnChar_[typeIn] = 'X'; + else + turnChar_[typeIn] = '>'; +} + +void Action_DSSP::SSres::SetTurnEnd(TurnType typeIn) { + if (turnChar_[typeIn] == '>') + turnChar_[typeIn] = 'X'; + else + turnChar_[typeIn] = '<'; +} + +void Action_DSSP::SSres::SetTurn(TurnType typeIn) { + // Do not overwrite an existing end character + if (turnChar_[typeIn] == ' ') { + if (typeIn == T3) turnChar_[typeIn] = '3'; + else if (typeIn == T4) turnChar_[typeIn] = '4'; + else if (typeIn == T5) turnChar_[typeIn] = '5'; + } +} + +int Action_DSSP::SSres::ssPriority(SStype typeIn) { + switch (typeIn) { + case ALPHA : return 8; + case BRIDGE : return 7; + case EXTENDED : return 6; + case H3_10 : return 5; + case HPI : return 4; + case TURN : return 3; + case BEND : return 2; + case NONE : return 1; + } + return 0; +} + +int Action_DSSP::SSres::SSpriority() const { + return ssPriority(sstype_); +} + +void Action_DSSP::SSres::SetSS(SStype typeIn) { + // TODO check if the priority check is necessary + //if (ssPriority(typeIn) > ssPriority(sstype_)) + sstype_ = typeIn; +} + +bool Action_DSSP::SSres::HasTurnStart(TurnType typeIn) const { + if (turnChar_[typeIn] == '>' || + turnChar_[typeIn] == 'X') + return true; + return false; +} + +void Action_DSSP::SSres::SetBridge(int idx, BridgeType btypeIn) { + if (bridge1idx_ == -1) { + bridge1idx_ = idx; + b1type_ = btypeIn; + } else if (bridge2idx_ == -1) { + bridge2idx_ = idx; + b2type_ = btypeIn; + } else + mprinterr("Error: Too many bridges for %i (to %i)\n", Num(), idx+1); +} + +bool Action_DSSP::SSres::HasBridge() const { + if (bridge1idx_ != -1) return true; + return false; +} + +bool Action_DSSP::SSres::IsBridgedWith(int idx2) const { + if (bridge1idx_ == idx2) return true; + if (bridge2idx_ == idx2) return true; + return false; +} + +/*char Action_DSSP::SSres::StrandChar() const { + // TODO ever b2? + return ssChar_[B1]; +}*/ + +void Action_DSSP::SSres::PrintSSchar() const { + static const char btypeChar[] = { ' ', 'p', 'A' }; + mprintf("\t%6i %c %c %c %c %c(%6i) %c(%6i) %6i %6i %6i %c\n", num_+1, resChar_, + turnChar_[T3], turnChar_[T4], turnChar_[T5], + btypeChar[b1type_], bridge1idx_+1, btypeChar[b2type_], bridge2idx_+1, + Bcount_[NO_BRIDGE], Bcount_[PARALLEL], Bcount_[ANTIPARALLEL], + DSSP_char_[sstype_]); +} + +#ifdef MPI +/** Sync residue SS and beta counts to master. */ +void Action_DSSP::SSres::SyncToMaster(Parallel::Comm const& commIn) { + int SSprob[NSSTYPE_ + NBRIDGETYPE_]; + int Buffer[NSSTYPE_ + NBRIDGETYPE_]; + std::copy( SScount_, SScount_ + NSSTYPE_, SSprob ); + std::copy( Bcount_, Bcount_ + NBRIDGETYPE_, SSprob + NSSTYPE_ ); + commIn.ReduceMaster( Buffer, SSprob, NSSTYPE_ + NBRIDGETYPE_, MPI_INT, MPI_SUM ); + if (commIn.Master()) { + std::copy( Buffer, Buffer + NSSTYPE_, SScount_ ); + std::copy( Buffer + NSSTYPE_, Buffer + NSSTYPE_ + NBRIDGETYPE_, Bcount_ ); + } +} +#endif -const int Action_DSSP::NSSTYPE = 8; +// ----- Action_DSSP ----------------------------------------------------------- -// CONSTRUCTOR Action_DSSP::Action_DSSP() : debug_(0), - outfile_(0), - dsspFile_(0), - assignout_(0), - Nres_(0), - Nselected_(0), - Nframe_(0), - printString_(false), BB_N_("N"), BB_H_("H"), BB_C_("C"), BB_O_("O"), - BB_CA_("CA") + BB_CA_("CA"), + outfile_(0), + dsspFile_(0), + assignout_(0), + printString_(false), + betaDetail_(false) {} +// Action_DSSP::Help() void Action_DSSP::Help() const { mprintf("\t[] [out ] [] [sumout ]\n" "\t[assignout ] [totalout ] [ptrajformat]\n" + "\t[betadetail]\n" "\t[namen ] [nameh ] [nameca ]\n" "\t[namec ] [nameo ]\n" " Calculate secondary structure content for residues in .\n" " If sumout not specified, the filename specified by out is used with .sum suffix.\n"); } -const char Action_DSSP::dssp_char[] = { ' ', 'E', 'B', 'G', 'H', 'I', 'T', 'S' }; -const char* Action_DSSP::SSchar[] = { "0", "b", "B", "G", "H", "I", "T", "S" }; -const char* Action_DSSP::SSname[]={"None", "Para", "Anti", "3-10", "Alpha", "Pi", "Turn", "Bend"}; - // Action_DSSP::Init() Action::RetType Action_DSSP::Init(ArgList& actionArgs, ActionInit& init, int debugIn) { debug_ = debugIn; - Nframe_ = 0; + Nframes_ = 0; // Get keywords outfile_ = init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs ); std::string temp = actionArgs.GetStringKey("sumout"); @@ -51,6 +306,7 @@ Action::RetType Action_DSSP::Init(ArgList& actionArgs, ActionInit& init, int deb DataFile* totalout = init.DFL().AddDataFile( actionArgs.GetStringKey("totalout"), actionArgs ); assignout_ = init.DFL().AddCpptrajFile(actionArgs.GetStringKey("assignout"), "SS assignment"); printString_ = actionArgs.hasKey("ptrajformat"); + betaDetail_ = actionArgs.hasKey("betadetail"); temp = actionArgs.GetStringKey("namen"); if (!temp.empty()) BB_N_ = temp; temp = actionArgs.GetStringKey("nameh"); @@ -70,11 +326,21 @@ Action::RetType Action_DSSP::Init(ArgList& actionArgs, ActionInit& init, int deb if (dsetname_.empty()) dsetname_ = init.DSL().GenerateDefaultName("DSSP"); // Set up Z labels - if (outfile_ != 0) - outfile_->ProcessArgs("zlabels None,Para,Anti,3-10,Alpha,Pi,Turn,Bend"); + if (outfile_ != 0) { + if (betaDetail_) + outfile_->ProcessArgs("zlabels None,Para,Anti,3-10,Alpha,Pi,Turn,Bend"); + else + outfile_->ProcessArgs("zlabels None,Ext,Bridge,3-10,Alpha,Pi,Turn,Bend"); + } // Create data sets for total fraction SS vs time. - for (int i = 0; i < NSSTYPE; i++) { - totalDS_[i] = init.DSL().AddSet(DataSet::FLOAT, MetaData(dsetname_, SSname[i])); + for (int i = 0; i < NSSTYPE_; i++) { + const char* aspect = DSSP_name_[i]; + if (betaDetail_ && (SStype)i == EXTENDED) + aspect = "Para"; + else if (betaDetail_ && (SStype)i == BRIDGE) + aspect = "Anti"; + SSname_.push_back( std::string(aspect) ); + totalDS_[i] = init.DSL().AddSet(DataSet::FLOAT, MetaData(dsetname_, aspect)); if (totalDS_[i] == 0) { mprinterr("Error: Could not create DSSP total frac v time data set.\n"); return Action::ERR; @@ -82,20 +348,36 @@ Action::RetType Action_DSSP::Init(ArgList& actionArgs, ActionInit& init, int deb // For now dont add 'None' so colors match up. if (totalout != 0 && i > 0) totalout->AddDataSet( totalDS_[i] ); } +# ifdef _OPENMP + // Each thread needs temp. space to store found hbonds every frame + // to avoid memory clashes when adding/updating in map. +# pragma omp parallel + { +# pragma omp master + { + CO_NH_bondsArray_.resize( omp_get_num_threads() ); + } + } +# endif mprintf( " SECSTRUCT: Calculating secondary structure using mask [%s]\n",Mask_.MaskString()); +# ifdef _OPENMP + mprintf("\tParallelizing calculation with %zu threads.\n", CO_NH_bondsArray_.size()); +# endif if (outfile_ != 0) mprintf("\tDumping results to %s\n", outfile_->DataFilename().full()); if (dsspFile_ != 0) mprintf("\tSum results to %s\n", dsspFile_->DataFilename().full()); + if (betaDetail_) + mprintf("\tWill print parallel/anti-parallel beta in place of extended/bridge\n"); if (printString_) { mprintf("\tSS data for each residue will be stored as a string.\n"); - for (int i = 0; i < NSSTYPE; i++) - mprintf("\t\t%s = %s\n", SSchar[i], SSname[i]); + for (int i = 0; i < NSSTYPE_; i++) + mprintf("\t\t%s = %s\n", SSchar_[i], SSname_[i].c_str()); } else { mprintf("\tSS data for each residue will be stored as integers.\n"); - for (int i = 0; i < NSSTYPE; i++) - mprintf("\t\t%i = %s\n", i, SSname[i]); + for (int i = 0; i < NSSTYPE_; i++) + mprintf("\t\t%i = %s\n", i, SSname_[i].c_str()); } if (assignout_ != 0) mprintf("\tOverall assigned SS will be written to %s\n", assignout_->Filename().full()); @@ -109,426 +391,646 @@ Action::RetType Action_DSSP::Init(ArgList& actionArgs, ActionInit& init, int deb return Action::OK; } +static inline void PrintAtom(Topology const& top, NameType const& name, int idx) { + if (idx > -1) + mprintf(" '%s'=%-12s", *name, top.AtomMaskName(idx/3).c_str()); + else + mprintf(" '%s'=%-12s", *name, "NONE"); +} + // Action_DSSP::Setup() /** Set up secondary structure calculation for all residues selected by the - * mask expression. A residue is selected if at least one of the following - * atoms named "C ", "O ", "N ", or "H " (i.e. standard atom protein - * BB atom names) is selected. The residue is only initialized if it has not - * been previously selected and set up by a prior call to setup. + * mask. A residue is selected if at least one atom in the residue is + * selected by the mask. The coordinate indices (i.e. atom # * 3) for + * the C, O, N, H, and CA atoms are set up if those atoms are present. + * The residue is only initialized if it has not been previously selected + * and set up by a prior call to Setup(). */ -// NOTE: Currently relatively memory-intensive. Eventually set up so that SecStruct and -// CO_HN_Hbond members exist only for selected residues? Use Map? -Action::RetType Action_DSSP::Setup(ActionSetup& setup) { +Action::RetType Action_DSSP::Setup(ActionSetup& setup) +{ // Set up mask for this parm - if ( setup.Top().SetupIntegerMask( Mask_ ) ) return Action::ERR; + if ( setup.Top().SetupCharMask( Mask_ ) ) return Action::ERR; if ( Mask_.None() ) { mprintf("Warning: DSSP: Mask has no atoms.\n"); return Action::SKIP; } - // Initially mark all residues already set up as not selected and - // reset all atom coordinate indices. - for (unsigned int res = 0; res != SecStruct_.size(); res++) { - SecStruct_[res].isSelected = false; - SecStruct_[res].C = -1; - SecStruct_[res].H = -1; - SecStruct_[res].N = -1; - SecStruct_[res].O = -1; - SecStruct_[res].CA = -1; - } - - // Set up SecStruct for each solute residue + // Deselect any existing residues + for (SSarrayType::iterator it = Residues_.begin(); it != Residues_.end(); ++it) + it->Deselect(); + // Set up for each solute residue Range soluteRes = setup.Top().SoluteResidues(); - Nres_ = soluteRes.Back() + 1; - if (debug_>0) mprintf("\tDSSP: Setting up for %i residues.\n", Nres_); - - // Set up for each residue of the current Parm if not already set-up. - SSres RES; - RES.sstype = NONE; - RES.isSelected = false; - RES.C = -1; - RES.O = -1; - RES.N = -1; - RES.H = -1; - RES.CA = -1; - RES.CO_HN_Hbond.assign( Nres_, 0 ); - std::fill( RES.SSprob, RES.SSprob + NSSTYPE, 0 ); - RES.resDataSet = 0; - // Only resize SecStruct if current # residues > previous # residues - if (Nres_ > (int)SecStruct_.size()) - SecStruct_.resize(Nres_, RES); - - // Go through all atoms in mask. Determine which residues have their C, - // O, N, or H atoms selected. Store the actual coordinate index - // (i.e. atom# * 3) instead of atom # for slight speed gain. - for (AtomMask::const_iterator atom = Mask_.begin(); atom!=Mask_.end(); ++atom) { - int res = setup.Top()[*atom].ResNum(); - // If residue is out of bounds skip it - if ( res < Nres_ ) { - //fprintf(stdout,"DEBUG: Atom %i Res %i [%s]\n",*atom,res,P->names[*atom]); - SecStruct_[res].isSelected = true; - if ( setup.Top()[*atom].Name() == BB_C_) - SecStruct_[res].C = (*atom) * 3; - else if ( setup.Top()[*atom].Name() == BB_O_) - SecStruct_[res].O = (*atom) * 3; - else if ( setup.Top()[*atom].Name() == BB_N_) - SecStruct_[res].N = (*atom) * 3; - else if ( setup.Top()[*atom].Name() == BB_H_) - SecStruct_[res].H = (*atom) * 3; - else if ( setup.Top()[*atom].Name() == BB_CA_) - SecStruct_[res].CA = (*atom) * 3; - } - } - - // For each residue selected in the mask, check if residue is missing atoms. - // Set up DataSet if necessary. - Nselected_ = 0; - std::vector missingResidues; + mprintf("\tSetting up for %i solute residues.\n", soluteRes.Size()); + if ((unsigned int)soluteRes.Size() > Residues_.size()) + Residues_.resize( soluteRes.Size() ); MetaData md(dsetname_, "res"); - DataSet::DataType dt = DataSet::INTEGER; - if (printString_) dt = DataSet::STRING; - for (int res = 0; res < Nres_; ++res) { - if (SecStruct_[res].isSelected) { - // Check if C-O/N-H selected. - SecStruct_[res].hasCO = (SecStruct_[res].C != -1 && - SecStruct_[res].O != -1); - SecStruct_[res].hasNH = (SecStruct_[res].N != -1 && - SecStruct_[res].H != -1); - if (!SecStruct_[res].hasCO || !SecStruct_[res].hasNH || SecStruct_[res].CA == -1) - { - missingResidues.push_back( setup.Top().TruncResNameNum( res ) ); - if (debug_ > 0) { - mprintf("Warning: Not all BB atoms found for res %s:", missingResidues.back().c_str()); - if (SecStruct_[res].N==-1) mprintf(" N"); - if (SecStruct_[res].H==-1) mprintf(" H"); - if (SecStruct_[res].C==-1) mprintf(" C"); - if (SecStruct_[res].O==-1) mprintf(" O"); - if (SecStruct_[res].CA==-1) mprintf(" CA"); - mprintf("\n"); + DataSet::DataType dt; + if (printString_) + dt = DataSet::STRING; + else + dt = DataSet::INTEGER; + unsigned int nResSelected = 0; + SSarrayType::iterator Res = Residues_.begin(); + for (Range::const_iterator ridx = soluteRes.begin(); ridx != soluteRes.end(); ++ridx, ++Res) + { + Residue const& thisRes = setup.Top().Res( *ridx ); + if (Res->Num() != -1) { + // Residue has previously been set up. Check that indices match. + if (Res->Num() != *ridx) { + mprinterr("Error: Solute residue index %i does not match previously setup\n" + "Error: index %i\n", *ridx, Res->Num()); + return Action::ERR; + } + } else { + // Set up Residue. TODO also molecule index? + Res->SetNum( *ridx ); + Res->SetResChar( thisRes.SingleCharName() ); + // Determine the previous and next residues + int prevresnum = -1; + int nextresnum = -1; + for (int at = thisRes.FirstAtom(); at != thisRes.LastAtom(); at++) { + if ( setup.Top()[at].Element() != Atom::HYDROGEN ) { + for (Atom::bond_iterator ib = setup.Top()[at].bondbegin(); + ib != setup.Top()[at].bondend(); ++ib) + { + if ( setup.Top()[*ib].ResNum() < *ridx ) { + if (prevresnum != -1) + mprintf("Warning: Multiple previous residues for res %i\n", *ridx+1); + else + prevresnum = setup.Top()[*ib].ResNum(); + } else if ( setup.Top()[*ib].ResNum() > *ridx ) { + if (nextresnum != -1) + mprintf("Warning: Multiple next residues for res %i\n", *ridx+1); + else + nextresnum = setup.Top()[*ib].ResNum(); + } + } } } - // Set up dataset if necessary - if (SecStruct_[res].resDataSet == 0) { - md.SetIdx( res+1 ); - md.SetLegend( setup.Top().TruncResNameNum(res) ); - // Setup dataset for this residue - SecStruct_[res].resDataSet = Init_.DSL().AddSet( dt, md ); - if (SecStruct_[res].resDataSet == 0) { - mprinterr("Error: Could not allocate DSSP data set for residue %i\n", res+1); +# ifdef DSSPDEBUG + mprintf("\t %8i < %8i < %8i\n", prevresnum+1, *ridx+1, nextresnum+1); +# endif + // Here we assume that residues are sequential! + if (prevresnum > -1) Res->SetPrevIdx( Res-Residues_.begin()-1 ); + if (nextresnum > -1) Res->SetNextIdx( Res-Residues_.begin()+1 ); + } + // Determine if this residue is selected + if (Mask_.AtomsInCharMask(thisRes.FirstAtom(), thisRes.LastAtom())) { + Res->SetSelected( true ); + ++nResSelected; + // Determine atom indices + for (int at = thisRes.FirstAtom(); at != thisRes.LastAtom(); at++) + { + if ( setup.Top()[at].Name() == BB_C_ ) Res->SetC( at*3 ); + else if ( setup.Top()[at].Name() == BB_O_ ) Res->SetO( at*3 ); + else if ( setup.Top()[at].Name() == BB_N_ ) Res->SetN( at*3 ); + else if ( setup.Top()[at].Name() == BB_H_ ) Res->SetH( at*3 ); + else if ( setup.Top()[at].Name() == BB_CA_ ) Res->SetCA( at*3 ); + } + // Check if residue is missing atoms + if (Res->IsMissingAtoms()) { + mprintf("Warning: Res %s is missing atoms", setup.Top().TruncResNameNum( *ridx ).c_str()); + if (Res->C() == -1) mprintf(" %s", *BB_C_); + if (Res->O() == -1) mprintf(" %s", *BB_N_); + if (Res->N() == -1) mprintf(" %s", *BB_O_); + if (Res->H() == -1) mprintf(" %s", *BB_H_); + if (Res->CA() == -1) mprintf(" %s", *BB_CA_); + mprintf("\n"); + } + // Set up DataSet if necessary + if (Res->Dset() == 0) { + md.SetIdx( *ridx+1 ); + md.SetLegend( setup.Top().TruncResNameNum( *ridx ) ); + // Setup DataSet for this residue + Res->SetDset( Init_.DSL().AddSet( dt, md ) ); + if (Res->Dset() == 0) { + mprinterr("Error: Could not allocate DSSP data set for residue %i\n", *ridx+1); return Action::ERR; } - if (outfile_ != 0) outfile_->AddDataSet(SecStruct_[res].resDataSet); + if (outfile_ != 0) outfile_->AddDataSet( Res->Dset() ); } - ++Nselected_; - } - } - if (!missingResidues.empty()) { - mprintf("Warning: Not all BB atoms found for %zu residues:", missingResidues.size()); - for (std::vector::const_iterator mr = missingResidues.begin(); - mr != missingResidues.end(); ++mr) - mprintf(" %s", mr->c_str()); - mprintf("\nInfo: This is expected for Proline and terminal/non-standard residues.\n" - "Info: Expected BB atom names: N=[%s] H=[%s] C=[%s] O=[%s] CA=[%s]\n", - *BB_N_, *BB_H_, *BB_C_, *BB_O_, *BB_CA_ ); - mprintf("Info: Re-run with action debug level >= 1 to see which residues are missing atoms.\n"); + } // END residue is selected } + mprintf("\t%u of %i solute residues selected.\n", nResSelected, soluteRes.Size()); - // Count number of selected residues - mprintf("\tMask [%s] corresponds to %u residues.\n", Mask_.MaskString(), Nselected_); -# ifdef DSSPDEBUG - // DEBUG - Print atom nums for each residue set up - for (int res=0; res < Nres_; res++) { - if (SecStruct_[res].isSelected) { - mprintf("DEBUG: Res %i", res + 1); - if (SecStruct_[res].hasCO) - mprintf(" C=%s O=%s",setup.Top().AtomMaskName(SecStruct_[res].C/3).c_str(), - setup.Top().AtomMaskName(SecStruct_[res].O/3).c_str()); - if (SecStruct_[res].hasNH) - mprintf(" N=%s H=%s",setup.Top().AtomMaskName(SecStruct_[res].N/3).c_str(), - setup.Top().AtomMaskName(SecStruct_[res].H/3).c_str()); - if (SecStruct_[res].CA != -1) - mprintf(" CA=%s",setup.Top().AtomMaskName(SecStruct_[res].CA/3).c_str()); + // DEBUG - print each residue set up. + if (debug_ > 0) { + for (SSarrayType::const_iterator it = Residues_.begin(); it != Residues_.end(); ++it) + { + mprintf(" %8i", it->Num() + 1); + PrintAtom(setup.Top(), BB_C_, it->C()); + PrintAtom(setup.Top(), BB_O_, it->O()); + PrintAtom(setup.Top(), BB_N_, it->N()); + PrintAtom(setup.Top(), BB_H_, it->H()); + PrintAtom(setup.Top(), BB_CA_, it->CA()); + mprintf(" Prev=%8i Next=%8i", it->PrevIdx(), it->NextIdx()); mprintf("\n"); } } -# endif + return Action::OK; } -// Action_DSSP::isBonded() -/** Return 1 if residue 1 CO bonded to residue 2 NH. - * Ensure residue numbers are valid and residues are selected. - */ -int Action_DSSP::isBonded(int res1, int res2) { - if (res1<0 || res2<0 || res1>=Nres_ || res2>=Nres_) return 0; - return SecStruct_[res1].CO_HN_Hbond[res2]; +void Action_DSSP::AssignBridge(int idx1in, int idx2in, BridgeType btypeIn) { + // By convention, always make idx1 the lower one + int idx1, idx2; + if (idx1in < idx2in) { + idx1 = idx1in; + idx2 = idx2in; + } else { + idx1 = idx2in; + idx2 = idx1in; + } + + SSres& Resi = Residues_[idx1]; + SSres& Resj = Residues_[idx2]; +# ifdef DSSPDEBUG + if (btypeIn == ANTIPARALLEL) { + mprintf("\t\tAssignBridge %i to %i, Antiparallel\n", idx1+1, idx2+1); + } else { + mprintf("\t\tAssignBridge %i to %i, Parallel\n", idx1+1, idx2+1); + } +# endif + // Do not duplicate bridges + if (Resi.IsBridgedWith(idx2)) { +# ifdef DSSPDEBUG + mprintf("\t\tAlready present.\n"); +# endif + return; + } + + Resi.SetBridge( idx2, btypeIn ); + Resj.SetBridge( idx1, btypeIn ); } -/// \return true if type1 has priority over type2. -bool Action_DSSP::HasPriority( SStype type1, SStype type2 ) { - switch (type1) { - case ALPHA: return true; - case ANTI: - if (type2 != ALPHA) return true; - break; - case PARA: - if (type2 != ALPHA && type2 != ANTI) return true; - break; - case H3_10: - if (type2 != ALPHA && type2 != ANTI && type2 != PARA) return true; - break; - case HPI: - if (type2 != ALPHA && type2 != ANTI && type2 != PARA && - type2 != H3_10) return true; - break; - case TURN: - if (type2 != ALPHA && type2 != ANTI && type2 != PARA && - type2 != H3_10 && type2 != HPI) return true; - break; - case BEND: - if (type2 != ALPHA && type2 != ANTI && type2 != PARA && - type2 != H3_10 && type2 != HPI && type2 != TURN) return true; - break; - case NONE: break; +/* +void Action_DSSP::AssignBridge(int idx1in, int idx2in, BridgeType btypeIn, char& currentStrandChar) { + // By convention, always make idx1 the lower one + int idx1, idx2; + if (idx1in < idx2in) { + idx1 = idx1in; + idx2 = idx2in; + } else { + idx1 = idx2in; + idx2 = idx1in; } - return false; + if (btypeIn == ANTIPARALLEL) + mprintf("\t\tAssignBridge %i to %i, Antiparallel\n", idx1+1, idx2+1); + else + mprintf("\t\tAssignBridge %i to %i, Parallel\n", idx1+1, idx2+1); + SSres& Resi = Residues_[idx1]; + SSres& Resj = Residues_[idx2]; + // Do not duplicate bridges + if (Resi.IsBridgedWith(idx2)) { + mprintf("\t\tAlready present.\n"); + return; + } + // Determine if we are already part of a ladder. + char ladderChar = ' '; + char resiLadderChar = Residues_[Resi.PrevIdx()].StrandChar(); + if (resiLadderChar == ' ') + resiLadderChar = Residues_[Resi.NextIdx()].StrandChar(); + char resjLadderChar = Residues_[Resj.PrevIdx()].StrandChar(); + if (resjLadderChar == ' ') + resjLadderChar = Residues_[Resj.NextIdx()].StrandChar(); + if (resiLadderChar == ' ' && resjLadderChar == ' ') { + // If both are blank, new ladder. + ladderChar = currentStrandChar; + ++currentStrandChar; + } else if (resiLadderChar != resjLadderChar) { + // They do not match. New ladder. + ladderChar = currentStrandChar; + ++currentStrandChar; + } else + ladderChar = resiLadderChar; + //else if (resiLadderChar != ' ') + // ladderChar = resiLadderChar; + //else + // ladderChar = resjLadderChar; + // Set the bridge; adjust character case if needed + if (btypeIn == ANTIPARALLEL) + ladderChar = toupper( ladderChar ); + mprintf("\t\tResi strand %c, resj strand %c, ladder char %c\n", resiLadderChar, resjLadderChar, ladderChar); + Resi.SetBridge( idx2, ladderChar ); + Resj.SetBridge( idx1, ladderChar ); } +*/ -// Action_DSSP::SSassign() -/** Assign all residues from res1 to res2-1 the secondary structure sstype - * only if it has priority. - * Assumes given residue range is valid. - */ -void Action_DSSP::SSassign(int res1, int res2, SStype typeIn, bool force) { -# ifdef DSSPDEBUG - mprintf("DEBUG:\tCalling SSassign from %i to %i, %s:", res1+1, res2, SSname[typeIn]); -# endif - for (int res = res1; res < res2; res++) { - if (res==Nres_) break; - if ( HasPriority(typeIn, SecStruct_[res].sstype) || force ) { -# ifdef DSSPDEBUG - mprintf(" %i", res+1); // DEBUG -# endif - SecStruct_[res].sstype = typeIn; - } +// TODO use Num()? +static inline int AbsResDelta(int idx1, int idx2) { + int resDelta = idx1 - idx2; + if (resDelta < 0) resDelta = -resDelta; + return resDelta; +} + +static inline void SetMin(int& resGapSize, int& sres0, int& sres1, int Nres, int Pres) +{ + int itmp = AbsResDelta(Nres, Pres); + if (itmp < resGapSize) { + resGapSize = itmp; + sres0 = std::min(Pres, Nres); + sres1 = std::max(Pres, Nres); } -# ifdef DSSPDEBUG - mprintf("\n"); // DEBUG -# endif } - -// Action_DSSP::DoAction() -/** Determine secondary structure by hydrogen bonding pattern. */ -Action::RetType Action_DSSP::DoAction(int frameNum, ActionFrame& frm) { - int resi, resj; - const double *C, *O, *H, *N; - double rON, rCH, rOH, rCN, E; - // Determine C=O to H-N hydrogen bonds for each residue to each other residue + + +int Action_DSSP::OverHbonds(int frameNum, ActionFrame& frm) +{ + t_total_.Start(); + t_calchb_.Start(); + // ----- Determine hydrogen bonding ------------ + int resi; + int Nres = (int)Residues_.size(); #ifdef _OPENMP -#pragma omp parallel private(resi,resj,C,O,H,N,rON, rCH, rOH, rCN, E) + int mythread; +#pragma omp parallel private(resi, mythread) { + mythread = omp_get_thread_num(); + CO_NH_bondsArray_[mythread].clear(); #pragma omp for -#endif - for (resi = 0; resi < Nres_; resi++) { - if (SecStruct_[resi].isSelected) { - // Reset previous SS assignment/Hbond status - SecStruct_[resi].sstype = NONE; - SecStruct_[resi].CO_HN_Hbond.assign( Nres_, 0 ); - if (SecStruct_[resi].hasCO) { - C = frm.Frm().CRD(SecStruct_[resi].C); - O = frm.Frm().CRD(SecStruct_[resi].O); - for (resj = 0; resj < Nres_; resj++) { - if (SecStruct_[resj].isSelected && resi != resj && SecStruct_[resj].hasNH) - { - N = frm.Frm().CRD(SecStruct_[resj].N); - H = frm.Frm().CRD(SecStruct_[resj].H); - - rON = 1.0/sqrt(DIST2_NoImage(O, N)); - rCH = 1.0/sqrt(DIST2_NoImage(C, H)); - rOH = 1.0/sqrt(DIST2_NoImage(O, H)); - rCN = 1.0/sqrt(DIST2_NoImage(C, N)); - - E = DSSP_fac * (rON + rCH - rOH - rCN); - if (E < -0.5) { -# ifdef DSSPDEBUG - mprintf("DEBUG: %i-CO --> %i-NH E= %g\n", resi+1, resj+1, E); -# endif - SecStruct_[resi].CO_HN_Hbond[resj] = 1; - } - } - } - } - } - } -#ifdef _OPENMP -} // END pragma omp parallel -#endif +#else /* _OPENMP */ + CO_NH_bonds_.clear(); +#endif /* _OPENMP */ + for (resi = 0; resi < Nres; resi++) + { + Residues_[resi].Unassign(); + if (Residues_[resi].IsSelected()) + { + SSres& ResCO = Residues_[resi]; + if (ResCO.HasCO()) + { + const double* Cxyz = frm.Frm().CRD( ResCO.C() ); + const double* Oxyz = frm.Frm().CRD( ResCO.O() ); + for (int resj = 0; resj < Nres; resj++) + { + // Need to consider adjacent residues since delta (2 res) turns possible + if (resi != resj) { + SSres& ResNH = Residues_[resj]; + if (ResNH.IsSelected() && ResNH.HasNH()) + { + const double* Nxyz = frm.Frm().CRD( ResNH.N() ); + const double* Hxyz = frm.Frm().CRD( ResNH.H() ); - // Determine Secondary Structure based on Hbonding pattern. - // In case of structural overlap, priority is given to the structure first - // in this list (see p. 2587 & 2595 in the Kabsch & Sander paper): - // H = 4-helix (alpha helix) - // B = residue in isolated beta bridge (anti) - // E = extended strand, participates in beta-ladder (para) - // G = 3-helix (310 helix) - // I = 5-helix (pi helix) - // T = H-bonded turn - // S = bend - for (resi=0; resi < Nres_; resi++) { - if (SecStruct_[resi].isSelected) { -# ifdef DSSPDEBUG - mprintf("DEBUG: Residue %i -----\n", resi+1); -# endif - // Alpha helices - if ( isBonded( resi - 1, resi+3 ) && isBonded( resi, resi + 4) ) - SSassign(resi, resi+4, ALPHA, false); - - // Beta sheets - only needed if SS not already assigned to alpha - if ( SecStruct_[resi].sstype != ALPHA ) { - for (resj=0; resj < Nres_; resj++) { - if (SecStruct_[resj].isSelected) { - // Only consider residues spaced more than 2 apart - int abs_resi_resj = resi - resj; - if (abs_resi_resj<0) abs_resi_resj = -abs_resi_resj; - if (abs_resi_resj > 2) { - if ( (isBonded(resi-1, resj) && isBonded(resj, resi+1)) || - (isBonded(resj-1, resi) && isBonded(resi, resj+1)) ) - { - // Parallel - // NOTE: Not checking if ANTI since not geometrically possible - // to be anti-parallel and parallel at the same time. + double rON = 1.0/sqrt(DIST2_NoImage(Oxyz, Nxyz)); + double rCH = 1.0/sqrt(DIST2_NoImage(Cxyz, Hxyz)); + double rOH = 1.0/sqrt(DIST2_NoImage(Oxyz, Hxyz)); + double rCN = 1.0/sqrt(DIST2_NoImage(Cxyz, Nxyz)); + + double E = DSSP_fac_ * (rON + rCH - rOH - rCN); + if (E < DSSP_cut_) { # ifdef DSSPDEBUG - mprintf("DEBUG:\tAssigning %i to parallel beta.\n", resi+1); + mprintf("DBG: %i-CO --> %i-NH E= %g\n", resi+1, resj+1, E); # endif - SecStruct_[resi].sstype = PARA; - break; - } else if ( (isBonded(resi-1, resj+1) && isBonded(resj-1, resi+1)) || - (isBonded(resi, resj ) && isBonded(resj, resi )) ) - { - // Anti-parallel -# ifdef DSSPDEBUG - mprintf("DEBUG:\tAssigning %i to anti-parallel beta.\n", resi+1); +# ifdef _OPENMP + CO_NH_bondsArray_[mythread].insert( HbondPairType(resi, resj) ); +# else + CO_NH_bonds_.insert( HbondPairType(resi, resj) ); # endif - SecStruct_[resi].sstype = ANTI; - break; } - } - } - } - } - - // 3-10 helix - if ( isBonded( resi - 1, resi+2 ) && isBonded( resi, resi + 3) ) - SSassign(resi, resi+3, H3_10, false); +//# ifdef DSSPDEBUG +// else if (resDelta < 6) +// mprintf("DBG: No hbond %i-CO --> %i-NH E= %g\n", resi+1, resj+1, E); +//# endif + } // END ResNH selected + } // END residues spaced > 2 apart + } // END inner loop over residues + } // END has CO + } // END ResCO selected + } // END outer loop over residues +#ifdef _OPENMP +} // END pragma omp parallel + for (unsigned int thread = 1; thread < CO_NH_bondsArray_.size(); thread++) + for (HbondMapType::const_iterator hb = CO_NH_bondsArray_[thread].begin(); + hb != CO_NH_bondsArray_[thread].end(); ++hb) + CO_NH_bondsArray_[0].insert( *hb ); + HbondMapType const& CO_NH_bonds_ = CO_NH_bondsArray_[0]; +#endif + t_calchb_.Stop(); - // Pi helix - if ( isBonded( resi - 1, resi+4 ) && isBonded( resi, resi + 5) ) - SSassign(resi, resi+5, HPI, false); - - // n-Turn, n=3,4,5 - for (int n=5; n > 2; n--) { - if ( isBonded(resi, resi + n) ) { - SSassign(resi+1, resi + n, TURN, false); // FIXME: Should this be resi? - break; - } + t_assign_.Start(); + // ----- Do basic assignment ------------------- + for (HbondMapType::const_iterator hb0 = CO_NH_bonds_.begin(); hb0 != CO_NH_bonds_.end(); ++hb0) + { + int riidx = hb0->first; + int rjidx = hb0->second; + SSres const& Resi = Residues_[riidx]; + SSres const& Resj = Residues_[rjidx]; +# ifdef DSSPDEBUG + mprintf("Res %8i %c -- %8i %c", Resi.Num()+1, Resi.ResChar(), + Resj.Num()+1, Resj.ResChar()); // DBG +# endif + // Spacing between residues i and j + int resDelta = Resj.Num() - Resi.Num(); +# ifdef DSSPDEBUG + mprintf("(%4i)\n", resDelta); +# endif + // Check for H bond from CO i to NH i+n + if (resDelta == 3) { + // 3-TURN + Residues_[riidx ].SetTurnBegin(T3); + Residues_[riidx+1].SetTurn(T3); + Residues_[riidx+2].SetTurn(T3); + Residues_[riidx+3].SetTurnEnd(T3); + Residues_[riidx+1].SetSS( TURN ); + Residues_[riidx+2].SetSS( TURN ); + } else if (resDelta == 4) { + // 4-TURN + Residues_[riidx ].SetTurnBegin(T4); + Residues_[riidx+1].SetTurn(T4); + Residues_[riidx+2].SetTurn(T4); + Residues_[riidx+3].SetTurn(T4); + Residues_[riidx+4].SetTurnEnd(T4); + Residues_[riidx+1].SetSS( TURN ); + Residues_[riidx+2].SetSS( TURN ); + Residues_[riidx+3].SetSS( TURN ); + } else if (resDelta == 5) { + // 5-TURN + Residues_[riidx ].SetTurnBegin(T5); + Residues_[riidx+1].SetTurn(T5); + Residues_[riidx+2].SetTurn(T5); + Residues_[riidx+3].SetTurn(T5); + Residues_[riidx+4].SetTurn(T5); + Residues_[riidx+5].SetTurnEnd(T5); + Residues_[riidx+1].SetSS( TURN ); + Residues_[riidx+2].SetSS( TURN ); + Residues_[riidx+3].SetSS( TURN ); + Residues_[riidx+4].SetSS( TURN ); + } + // Look for bridge. Start with the premise that this bond is part of one + // of the 4 potential bridge patterns, then check if the compliment exists. + HbondMapType::iterator hb; + // Here we want absolute value of spacing. + if (resDelta < 0) resDelta = -resDelta; + if (resDelta > 2) { + // Assume (i,j). Look for (j,i) + hb = CO_NH_bonds_.find( HbondPairType(hb0->second, hb0->first) ); + if (hb != CO_NH_bonds_.end()) { +# ifdef DSSPDEBUG + mprintf("\t\t%i ANTI-PARALLELa with %i (%i)\n", hb0->first+1, hb0->second+1, hb->first+1); +# endif + AssignBridge(hb0->first, hb0->second, ANTIPARALLEL); } + } + resDelta = AbsResDelta(hb0->first+1, hb0->second-1); + if (resDelta > 2) { + // Assume (i-1, j+1). Look for (j-1, i+1) + hb = CO_NH_bonds_.find( HbondPairType(hb0->second-2, hb0->first+2) ); + if (hb != CO_NH_bonds_.end()) { +# ifdef DSSPDEBUG + mprintf("\t\t%i ANTI-PARALLELb with %i (%i)\n", hb0->first+2, hb0->second, hb->first+1); +# endif + AssignBridge(hb0->first+1, hb0->second-1, ANTIPARALLEL); + } + } + resDelta = AbsResDelta(hb0->first+1, hb0->second); + if (resDelta > 2) { + // Assume (i-1, j). Check for (j, i+1) PARALLEL + hb = CO_NH_bonds_.find( HbondPairType(hb0->second, hb0->first+2) ); + if (hb != CO_NH_bonds_.end()) { +# ifdef DSSPDEBUG + mprintf("\t\t%i PARALLELa with %i (%i)\n", hb0->first+2, hb0->second+1, hb->first+1); +# endif + AssignBridge(hb0->first+1, hb0->second, PARALLEL); + } + } + resDelta = AbsResDelta(hb0->second-1, hb0->first); + if (resDelta > 2) { + // Assume (j, i+1). Check for (i-1, j) + hb = CO_NH_bonds_.find( HbondPairType(hb0->second-2, hb0->first) ); + if (hb != CO_NH_bonds_.end()) { +# ifdef DSSPDEBUG + mprintf("\t\t%i PARALLELb with %i (%i)\n", hb0->second, hb0->first+1, hb->first+1); +# endif + AssignBridge(hb0->second-1, hb0->first, PARALLEL); + } + } + } // END loop over Hbonds - // Bend - has lowest priority, so only do if no assignment. - if (SecStruct_[resi].sstype == NONE && resi > 1 && resi < Nres_ - 2) - { - if (SecStruct_[resi-2].CA != -1 && - SecStruct_[resi ].CA != -1 && - SecStruct_[resi+2].CA != -1) - { - const double* CAm2 = frm.Frm().CRD(SecStruct_[resi-2].CA); - const double* CA0 = frm.Frm().CRD(SecStruct_[resi ].CA); - const double* CAp2 = frm.Frm().CRD(SecStruct_[resi+2].CA); - Vec3 CA1( CA0[0]-CAm2[0], CA0[1]-CAm2[1], CA0[2]-CAm2[2] ); - Vec3 CA2( CAp2[0]-CA0[0], CAp2[1]-CA0[1], CAp2[2]-CA0[2] ); - CA1.Normalize(); - CA2.Normalize(); - // 1.221730476 rad = 70 degrees - if (CA1.Angle(CA2) > 1.221730476) { + // ----- Do SS assignment ---------------------- + // Priority is 'H', 'B', 'E', 'G', 'I', 'T', 'S' None + // 8 7 6 5 4 3 2 1 + for (resi = 0; resi < Nres; resi++) + { +# ifdef DSSPDEBUG + mprintf("Residue %i\n", resi+1); +# endif + SSres& Resi = Residues_[resi]; + int prevRes = resi - 1; + int nextRes = resi + 1; + int priority = Resi.SSpriority(); + if ( Resi.HasTurnStart(T4) && prevRes > -1 && Residues_[prevRes].HasTurnStart(T4) ) + { + // Alpha helix. +# ifdef DSSPDEBUG + mprintf("ALPHA helix starting at %i\n", resi+1); +# endif + Residues_[resi ].SetSS( ALPHA ); + Residues_[resi+1].SetSS( ALPHA ); + Residues_[resi+2].SetSS( ALPHA ); + Residues_[resi+3].SetSS( ALPHA ); + } else if (Resi.SS() != ALPHA) { + if (priority < 6) { + // Priority < 6 means not alpha or beta assigned yet. + // Check for Beta structure + bool prevHasBridge = (prevRes > -1 && Residues_[prevRes].HasBridge()); + bool nextHasBridge = (nextRes < Nres && Residues_[nextRes].HasBridge()); + if (Resi.HasBridge()) { + // Regular Beta. Check if previous is assigned EXTENDED in case it + // was assigned via a Beta bulge. + if ( prevHasBridge || nextHasBridge || Residues_[prevRes].SS() == EXTENDED ) + { # ifdef DSSPDEBUG - mprintf("DEBUG: Bend calc %i-%i-%i: %g rad.\n", resi-1, resi+1, resi+3, CA1.Angle(CA2)); + mprintf("Extended BETA bridge at %i\n", resi+1); # endif - SecStruct_[resi].sstype = BEND; + Resi.SetSS( EXTENDED ); + } else { +# ifdef DSSPDEBUG + mprintf("Isolated BETA bridge at %i.\n", resi+1); +# endif + Resi.SetSS( BRIDGE ); } - } - } - } - } // End Initial SS assignment over all residues - - // Change 3-10 and Pi helices that are less than minimal size to Turn - SStype lastType = NONE; - int resStart = -1; - for (resi = 0; resi < Nres_; resi++) { - if (SecStruct_[resi].isSelected) { - if (lastType != SecStruct_[resi].sstype) { - // Secondary structure type has changed. - if (lastType == H3_10) { + } else if (prevHasBridge && nextHasBridge) { + // Potential Beta bulge. Check other strand. + int presb1 = Residues_[prevRes].Bridge1Idx(); + int presb2 = Residues_[prevRes].Bridge2Idx(); + int nresb1 = Residues_[nextRes].Bridge1Idx(); + int nresb2 = Residues_[nextRes].Bridge2Idx(); + int prest1 = Residues_[prevRes].Bridge1Type(); + int prest2 = Residues_[prevRes].Bridge2Type(); + int nrest1 = Residues_[nextRes].Bridge1Type(); + int nrest2 = Residues_[nextRes].Bridge2Type(); # ifdef DSSPDEBUG - mprintf("DEBUG: 3-10 helix length is %i\n", resi - resStart); + mprintf("Potential bulge? Prev res bridge res: %i %i Next res bridge res: %i %i\n", presb1, presb2, nresb1, nresb2); # endif - if (resi - resStart < 3) - SSassign(resStart, resi, TURN, true); - } else if (lastType == HPI) { + // The largest allowed gap in the other strand is 5 residues. + // Since we know that next res and previous res both have at + // least 1 bridge, Next B1 - Prev B1 can be the gap to beat. + // Need to also make sure the bridge types match. + int resGapSize = -1; + int sres0 = -1; + int sres1 = -1; + if (nrest1 == prest1) { + resGapSize = AbsResDelta(nresb1, presb1); + sres0 = std::min(presb1, nresb1); + sres1 = std::max(presb1, nresb1); + } + if (presb2 != -1 && nrest1 == prest2) + SetMin( resGapSize, sres0, sres1, nresb1, presb2 ); + if (nresb2 != -1) { + if (nrest2 == prest1) + SetMin( resGapSize, sres0, sres1, nresb2, presb1 ); + if (presb2 != -1 && nrest2 == prest2) + SetMin( resGapSize, sres0, sres1, nresb2, presb2 ); + } # ifdef DSSPDEBUG - mprintf("DEBUG: PI helix length is %i\n", resi - resStart); + mprintf("Min res gap size on other strand = %i (%i to %i)\n", resGapSize, sres0, sres1); # endif - if (resi - resStart < 5) - SSassign(resStart, resi, TURN, true); + // Minimum allowed gap is 4 residues in between, so 5 residues total. + if (resGapSize > -1 && resGapSize < 6) { +# ifdef DSSPDEBUG + mprintf("Beta bulge.\n"); +# endif + if (Residues_[prevRes].SS() != ALPHA) + Residues_[prevRes].SetSS( EXTENDED ); + Resi.SetSS( EXTENDED ); + if (Residues_[nextRes].SS() != ALPHA) + Residues_[nextRes].SetSS( EXTENDED ); + // Set extended on other strand as well + for (int sres = sres0; sres != sres1; sres++) + if (Residues_[sres].SS() != ALPHA) + Residues_[sres].SetSS( EXTENDED ); + } } - resStart = resi; -# ifdef DSSPDEBUG - mprintf("DEBUG: ResStart=%i for type %s\n", resi+1, SSname[SecStruct_[resi].sstype]); -# endif - } - lastType = SecStruct_[resi].sstype; + } // END check for Beta structure + } // END not alpha + } // END loop over residues + + // Check for 3-10 helices. Do this separately so we dont assign regions + // that are too small because other residues have already been assigned. + for (resi = 1; resi < Nres-2; resi++) { + if (Residues_[resi ].SSpriority() < 6 && + Residues_[resi+1].SSpriority() < 6 && + Residues_[resi+2].SSpriority() < 6 && + Residues_[resi ].HasTurnStart(T3) && + Residues_[resi-1].HasTurnStart(T3)) + { + // 3-10 helix +# ifdef DSSPDEBUG + mprintf("3-10 helix starting at %i\n", resi+1); +# endif + Residues_[resi ].SetSS( H3_10 ); + Residues_[resi+1].SetSS( H3_10 ); + Residues_[resi+2].SetSS( H3_10 ); } } - - // Store data for each residue. Get statistics. - int totalSS[NSSTYPE]; - std::fill( totalSS, totalSS + NSSTYPE, 0 ); - for (resi=0; resi < Nres_; resi++) { - if (SecStruct_[resi].isSelected) { - totalSS[SecStruct_[resi].sstype]++; - SecStruct_[resi].SSprob[SecStruct_[resi].sstype]++; - if (printString_) - SecStruct_[resi].resDataSet->Add(frameNum, SSchar[SecStruct_[resi].sstype]); - else - SecStruct_[resi].resDataSet->Add(frameNum, &(SecStruct_[resi].sstype)); + // Check for PI helices, similar to 3-10 helices. + for (resi = 1; resi < Nres-4; resi++) { + if (Residues_[resi ].SSpriority() < 5 && + Residues_[resi+1].SSpriority() < 5 && + Residues_[resi+2].SSpriority() < 5 && + Residues_[resi+3].SSpriority() < 5 && + Residues_[resi+4].SSpriority() < 5 && + Residues_[resi ].HasTurnStart(T5) && + Residues_[resi-1].HasTurnStart(T5)) + { + // PI helix +# ifdef DSSPDEBUG + mprintf("PI helix starting at %i\n", resi+1); +# endif + Residues_[resi ].SetSS( HPI ); + Residues_[resi+1].SetSS( HPI ); + Residues_[resi+2].SetSS( HPI ); + Residues_[resi+3].SetSS( HPI ); + Residues_[resi+4].SetSS( HPI ); + } + } + // Check for bends. Only do if no other assignment. + for (resi = 0; resi < Nres; resi++) { + if (Residues_[resi].IsSelected() && Residues_[resi].SS() == NONE) + { + int im2 = resi - 2; + if (im2 > -1) { + int ip2 = resi + 2; + if (ip2 < Nres) { + SSres& Resi = Residues_[resi]; + if (Residues_[im2].CA() != -1 && Resi.CA() != -1 && Residues_[ip2].CA() != -1) { + const double* CAm2 = frm.Frm().CRD(Residues_[im2].CA()); + const double* CA0 = frm.Frm().CRD(Resi.CA()); + const double* CAp2 = frm.Frm().CRD(Residues_[ip2].CA()); + Vec3 CA1( CA0[0]-CAm2[0], CA0[1]-CAm2[1], CA0[2]-CAm2[2] ); + Vec3 CA2( CAp2[0]-CA0[0], CAp2[1]-CA0[1], CAp2[2]-CA0[2] ); + CA1.Normalize(); + CA2.Normalize(); + double bAngle = CA1.Angle(CA2); +# ifdef DSSPDEBUG + mprintf("DEBUG: Bend calc %i-%i-%i: %g deg.\n", resi-1, resi+1, resi+3, bAngle*Constants::RADDEG); +# endif + // 1.221730476 rad = 70 degrees + if (bAngle > 1.221730476) { + Resi.SetSS( BEND ); + } + } + } + } + } // END selected and no assignment + } + + // ----- Store data for each res. Get stats ---- + int totalSS[NSSTYPE_]; + std::fill( totalSS, totalSS + NSSTYPE_, 0 ); + int Nselected = 0; + for (resi=0; resi < Nres; resi++) { + SSres& Resi = Residues_[resi]; + if (Resi.IsSelected()) { + if (betaDetail_ && (Resi.SS() == EXTENDED || Resi.SS() == BRIDGE)) + { + if (Resi.Bridge1Type() == ANTIPARALLEL || + Resi.Bridge2Type() == ANTIPARALLEL) + totalSS[BRIDGE]++; + else if (Resi.Bridge1Type() == PARALLEL || + Resi.Bridge2Type() == PARALLEL) + totalSS[EXTENDED]++; + } else + totalSS[Residues_[resi].SS()]++; + Residues_[resi].AccumulateData(frameNum, printString_, betaDetail_); + Nselected++; } } - for (int i = 0; i < NSSTYPE; i++) { + for (int i = 0; i < NSSTYPE_; i++) { float fvar = (float)totalSS[i]; - fvar /= (float)Nselected_; + fvar /= (float)Nselected; totalDS_[i]->Add(frameNum, &fvar); } - ++Nframe_; + + t_assign_.Stop(); + t_total_.Stop(); + return 0; +} + +Action::RetType Action_DSSP::DoAction(int frameNum, ActionFrame& frm) +{ + OverHbonds(frameNum, frm); + Nframes_++; + // DEBUG - Print basic assignment + if (debug_ > 1) { + for (SSarrayType::const_iterator it = Residues_.begin(); it != Residues_.end(); ++it) + it->PrintSSchar(); + } return Action::OK; } #ifdef MPI int Action_DSSP::SyncAction() { - // Consolidate SSprob data to master. - int SSprob[8]; - for (std::vector::iterator res = SecStruct_.begin(); res != SecStruct_.end(); ++res) - { - std::fill(SSprob, SSprob+8, 0); - Init_.TrajComm().ReduceMaster( SSprob, res->SSprob, 8, MPI_INT, MPI_SUM ); - if (Init_.TrajComm().Master()) { - for (int idx = 0; idx != 8; idx++) - res->SSprob[idx] = SSprob[idx]; - } - } + // Consolidate SScount data to master. + for (SSarrayType::iterator res = Residues_.begin(); res != Residues_.end(); ++res) + res->SyncToMaster( Init_.TrajComm() ); + // Calc total number of frames. int total_frames = 0; - Init_.TrajComm().ReduceMaster( &total_frames, &Nframe_, 1, MPI_INT, MPI_SUM ); + Init_.TrajComm().ReduceMaster( &total_frames, &Nframes_, 1, MPI_INT, MPI_SUM ); if (Init_.TrajComm().Master()) - Nframe_ = total_frames; + Nframes_ = total_frames; return 0; } #endif @@ -536,12 +1038,16 @@ int Action_DSSP::SyncAction() { // Action_DSSP::Print() void Action_DSSP::Print() { if (dsetname_.empty()) return; + t_total_.WriteTiming(1,"DSSP total"); + t_calchb_.WriteTiming(2, "Calc Hbonds", t_total_.Total()); + t_assign_.WriteTiming(2, "Assignment ", t_total_.Total()); + // Try not to print empty residues. Find the minimum and maximum residue // for which there is data. Output res nums start from 1. int min_res = -1; int max_res = -1; - for (int resi = 0; resi != (int)SecStruct_.size(); resi++) { - if (SecStruct_[resi].resDataSet != 0) { + for (int resi = 0; resi != (int)Residues_.size(); resi++) { + if (Residues_[resi].Dset() != 0) { if (min_res < 0) min_res = resi; if (resi > max_res) max_res = resi; } @@ -552,25 +1058,38 @@ void Action_DSSP::Print() { } // Calculate average of each SS type across all residues. if (dsspFile_ != 0) { - std::vector dsspData_(NSSTYPE); + std::vector dsspData_(NSSTYPE_); Dimension Xdim( min_res + 1, 1, "Residue" ); MetaData md(dsetname_, "avgss", MetaData::NOT_TS); // Set up a dataset for each SS type. TODO: NONE type? - for (int ss = 1; ss < NSSTYPE; ss++) { + for (int ss = 1; ss < NSSTYPE_; ss++) { md.SetIdx(ss); - md.SetLegend( SSname[ss] ); + md.SetLegend( SSname_[ss] ); dsspData_[ss] = Init_.DSL().AddSet(DataSet::DOUBLE, md); dsspData_[ss]->SetDim(Dimension::X, Xdim); dsspFile_->AddDataSet( dsspData_[ss] ); } // Calc the avg SS type for each residue that has data. - int idx = 0; + int idx = 0; + unsigned int norm = Nframes_; for (int resi = min_res; resi < max_res+1; resi++) { - if (SecStruct_[resi].resDataSet != 0) { - for (int ss = 1; ss < NSSTYPE; ss++) { - double avg = (double)SecStruct_[resi].SSprob[ss]; - avg /= (double)Nframe_; + SSres& Resi = Residues_[resi]; + if (Resi.Dset() != 0) { + //int Nframe = 0; + //for (int ss = 0; ss < NSSTYPE_; ss++) + // Nframe += Resi.SScount((SStype)ss); + //mprintf("DEBUG: Total frames for residue %i is %i\n", Resi.Num()+1, Nframe); + for (int ss = 1; ss < NSSTYPE_; ss++) { + double avg; + if (betaDetail_ && (SStype)ss == EXTENDED) + avg = (double)Resi.Bcount(PARALLEL); + else if (betaDetail_ && (SStype)ss == BRIDGE) + avg = (double)Resi.Bcount(ANTIPARALLEL); + else + avg = (double)Resi.SScount((SStype)ss); + //mprintf("DEBUG:\t\tCount for type %i is %f\n", ss, avg); + avg /= (double)norm; dsspData_[ss]->Add(idx, &avg); } ++idx; @@ -585,18 +1104,26 @@ void Action_DSSP::Print() { for (int resi = min_res; resi < max_res+1; resi++) { if (startRes == -1) startRes = resi; // Convert residue name. - resLine += Residue::ConvertResName( SecStruct_[resi].resDataSet->Meta().Legend() ); + SSres& Resi = Residues_[resi]; + resLine += Resi.ResChar(); // Figure out which SS element is dominant for res if selected - if (SecStruct_[resi].resDataSet != 0) { + if (Resi.Dset() != 0) { int dominantType = 0; int ssmax = 0; - for (int ss = 0; ss < NSSTYPE; ss++) { - if ( SecStruct_[resi].SSprob[ss] > ssmax ) { - ssmax = SecStruct_[resi].SSprob[ss]; + for (int ss = 0; ss < NSSTYPE_; ss++) { + int sscount; + if (betaDetail_ && (SStype)ss == EXTENDED) + sscount = Resi.Bcount(PARALLEL); + else if (betaDetail_ && (SStype)ss == BRIDGE) + sscount = Resi.Bcount(ANTIPARALLEL); + else + sscount = Resi.SScount((SStype)ss); + if ( sscount > ssmax ) { + ssmax = sscount; dominantType = ss; } } - ssLine += dssp_char[dominantType]; + ssLine += DSSP_char_[dominantType]; } else ssLine += '-'; total++; diff --git a/src/Action_DSSP.h b/src/Action_DSSP.h index ba68a46a18..ca552d5c0d 100644 --- a/src/Action_DSSP.h +++ b/src/Action_DSSP.h @@ -1,7 +1,31 @@ #ifndef INC_ACTION_DSSP_H #define INC_ACTION_DSSP_H +#include +#include +#include #include "Action.h" -/// Calculate protein secondary structure using DSSP algorithm. +#include "NameType.h" +#include "CharMask.h" +#include "Timer.h" +class DataSet; +/// Do protein secondary structure assignment. +/** Based on protein secondary structure definitions given in: + * Kabsch, W.; Sander, C.; \"Dictionary of Protein Secondary Structure: + * Pattern Recognition of Hydrogen-Bonded and Geometrical Features. + * Biopolymers (1983), V.22, pp.2577-2637. + * This is the second major version of this algorithm in CPPTRAJ; it + * has been completely rewritten by DRR. It correctly implements + * detection of Extended and Bridge regions, including beta bulges. + * The various secondary structure types that can be assigned (in order + * of decreasing priority) are: + * H - Alpha helix (4) + * B - Single bridge (ladder of length 1) + * E - Extended beta (ladder length > 1) + * G - 3-10 helix (3) + * I - Pi helix (5) + * T - Turn (3, 4, 5) + * S - Bend (Angle CA i-2, i, i+2 > 70 deg.) + */ class Action_DSSP : public Action { public: Action_DSSP(); @@ -15,51 +39,155 @@ class Action_DSSP : public Action { int SyncAction(); # endif void Print(); - // Enum and static vars - //enum SStype { ALPHA=0, ANTI, PARA, H3_10, HPI, TURN, BEND, NONE }; - enum SStype { NONE=0, PARA, ANTI, H3_10, ALPHA, HPI, TURN, BEND }; - static const int NSSTYPE; ///< # of secondary structure types. - static const double DSSP_fac; ///< DSSP factor for calc. Hbond energy. - static const char dssp_char[]; ///< DSSP 1 character SS names - static const char* SSchar[]; ///< PTRAJ 1 character SS names - static const char* SSname[]; ///< Full SS names - /// Hold SS-related data for each residue - struct SSres { - std::vector CO_HN_Hbond; ///< 1 if this res CO bonded to res X NH. - DataSet* resDataSet; ///< DataSet for SS assignment each frame. - int SSprob[8]; ///< Hold count for each SS type - SStype sstype; ///< Assigned secondary structure - int C; ///< Coord idx of BB carbon - int O; ///< Coord idx of BB oxygen - int N; ///< Coord idx of BB nitrogen - int H; ///< Coord idx of BB hydrogen - int CA; ///< Coord idx of BB alpha carbon - bool isSelected; ///< True if calculating SS for this residue. - bool hasCO; ///< True if both C and O atoms selected. - bool hasNH; ///< True if both N and H atoms selected. - }; - std::vector SecStruct_; ///< Hold SS-related data for all residues - // Class variables - int debug_; - DataFile* outfile_; ///< Output Data file - DataFile* dsspFile_; ///< Sum output file - std::string dsetname_; ///< DSSP data set name - CpptrajFile* assignout_; ///< Assignment output file. - AtomMask Mask_; ///< Mask used to determine selected residues - int Nres_; ///< Current total # of residues - unsigned int Nselected_; ///< Current # residues selected. - int Nframe_; ///< # of frames, for calculating SS avg. - bool printString_; ///< If true print 1 char per residue indicating ss type - ActionInit Init_; ///< Hold pointers to master DSL/DFL - DataSet* totalDS_[8]; - NameType BB_N_; - NameType BB_H_; - NameType BB_C_; - NameType BB_O_; - NameType BB_CA_; - // Private fns - inline int isBonded(int, int); - inline void SSassign(int, int, SStype, bool); - static inline bool HasPriority(SStype, SStype); + + /// Class that will hold SS info for each residue + class SSres; + /// Secondary structure types + enum SStype { NONE=0, EXTENDED, BRIDGE, H3_10, ALPHA, HPI, TURN, BEND }; + static const int NSSTYPE_ = 8; ///< # of secondary structure types. + static const char DSSP_char_[]; ///< DSSP 1 char names corresponding to SStype + static const char* DSSP_name_[]; ///< Full secondary structure names corresponding to SStype + static const char* SSchar_[]; ///< PTRAJ 1 character names corresponding to SStype + /// Turn types + enum TurnType { T3 = 0, T4, T5 }; + static const int NTURNTYPE_ = 3; + /// Beta types + enum BetaType { B1 = 0, B2, S }; + static const int NBETATYPE_ = 3; + /// Bridge direction types + enum BridgeType { NO_BRIDGE = 0, PARALLEL, ANTIPARALLEL }; + static const int NBRIDGETYPE_ = 3; + + static const double DSSP_fac_; ///< Original DSSP factor for calc. H-bond "energy" + static const double DSSP_cut_; ///< Original DSSP H-bond energy cutoff in kcal/mol + + typedef std::vector SSarrayType; + typedef std::vector Sarray; + typedef std::pair HbondPairType; + typedef std::set HbondMapType; + + int OverHbonds(int, ActionFrame&); + + void AssignBridge(int, int, BridgeType); + + /// Map resi (CO) to resj (NH) hydrogen bonds +# ifdef _OPENMP + std::vector CO_NH_bondsArray_; +# else + HbondMapType CO_NH_bonds_; +# endif + int debug_; ///< Action debug level + unsigned int Nframes_; ///< Number of frames processed, for total SS normalization + SSarrayType Residues_; ///< Hold SS data for all residues. + Sarray SSname_; ///< Hold full secondary structure names + NameType BB_N_; ///< Protein N atom name ('N') + NameType BB_H_; ///< Protein N-H atom name ('H') + NameType BB_C_; ///< Protein C atom name ('C') + NameType BB_O_; ///< Protein C-O atom name ('O') + NameType BB_CA_; ///< Protein alpha C name ('CA') + CharMask Mask_; ///< Mask used to determine selected residues. + DataFile* outfile_; ///< Output Data file + DataFile* dsspFile_; ///< Sum output file + CpptrajFile* assignout_; ///< Assignment output file. + std::string dsetname_; ///< DSSP data set name + DataSet* totalDS_[NSSTYPE_]; ///< Hold total SS each frame for each SS type + ActionInit Init_; ///< Hold pointers to master DSL/DFL + bool printString_; ///< If true print 1 char per residue indicating ss type + bool betaDetail_; ///< If true use para/anti in place of extended/bridge + Timer t_total_; + Timer t_calchb_; + Timer t_assign_; +}; + +// ============================================================================= +class Action_DSSP::SSres { + public: + SSres(); + SSres(SSres const&); + SSres& operator=(SSres const&); + int SScount(SStype t) const { return SScount_[t]; } + int Bcount(BridgeType t) const { return Bcount_[t]; } + SStype SS() const { return sstype_; } + int Num() const { return num_; } + char ResChar() const { return resChar_; } + bool IsSelected() const { return isSelected_; } + int C() const { return C_; } + int O() const { return O_; } + int N() const { return N_; } + int H() const { return H_; } + int CA() const { return CA_; } + int PrevIdx() const { return prevIdx_; } + int NextIdx() const { return nextIdx_; } + int Bridge1Idx() const { return bridge1idx_; } + BridgeType Bridge1Type() const { return b1type_;} + int Bridge2Idx() const { return bridge2idx_; } + BridgeType Bridge2Type() const { return b2type_;} + DataSet* Dset() const { return resDataSet_; } + + bool IsMissingAtoms() const { return (C_==-1 || O_==-1 || N_==-1 || H_==-1 || CA_==-1); } + bool HasCO() const { return (C_!=-1 && O_!=-1); } + bool HasNH() const { return (N_!=-1 && H_!=-1); } + bool HasCA() const { return (CA_!=-1); } + + bool HasTurnStart(TurnType) const; + bool HasBridge() const; + bool IsBridgedWith(int) const; +// char StrandChar() const; + void PrintSSchar() const; + + void SetTurnBegin(TurnType); + void SetTurn(TurnType); + void SetTurnEnd(TurnType); + /// Set a bridge between this res and other res index into Residues_ + void SetBridge(int, BridgeType); + + void AccumulateData(int, bool, bool); + + void SetSS(SStype); + void SetNum(int i) { num_ = i; } + void SetResChar(char c) { resChar_ = c; } + void SetSelected(bool b) { isSelected_ = b; } + void SetC(int i) { C_ = i; } + void SetO(int i) { O_ = i; } + void SetN(int i) { N_ = i; } + void SetH(int i) { H_ = i; } + void SetCA(int i) { CA_ = i; } + void SetPrevIdx(int i) { prevIdx_ = i; } + void SetNextIdx(int i) { nextIdx_ = i; } + void SetDset(DataSet* d) { resDataSet_ = d; } + /// Deselect this residue and reset coordinate indices. + void Deselect(); + /// Reset hbonds, pattern and SS type assignments + void Unassign(); + /// \return Relative priority of currently assigned SS type + int SSpriority() const; +# ifdef MPI + void SyncToMaster(Parallel::Comm const&); +# endif + private: + /// \return Relative priority of given SS type + static inline int ssPriority(SStype); + + DataSet* resDataSet_; ///< DataSet for SS assignment each frame for this res. + double chirality_; ///< Dihedral CA[i-1, i, i+1, i+2] + double bend_; ///< Angle CA[i-2, i, i+2] + int SScount_[NSSTYPE_]; ///< Hold count for each SS type + int Bcount_[NBRIDGETYPE_]; ///< Hold count for Beta types + SStype sstype_; ///< SS assignment for this frame + int num_; ///< Residue index in Topology + int C_; ///< Coord idx of BB carbon + int O_; ///< Coord idx of BB oxygen + int N_; ///< Coord idx of BB nitrogen + int H_; ///< Coord idx of BB hydrogen + int CA_; ///< Coord idx of BB alpha carbon + int prevIdx_; ///< Index in Residues_ of previous residue + int nextIdx_; ///< Index in Residues_ of next residue + int bridge1idx_; ///< Index in Residues_ of res this is bridged to + BridgeType b1type_; ///< Type of bridge1 + int bridge2idx_; ///< Index in Residues_ of res this is bridged to + BridgeType b2type_; ///< Type of bridge2 + char resChar_; ///< Single char residue ID + char turnChar_[NTURNTYPE_]; ///< Character if part of N turn + bool isSelected_; ///< True if calculating SS for this residue. }; #endif diff --git a/src/Version.h b/src/Version.h index d1d8d986ce..10d1a0d9f6 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V4.17.0" +#define CPPTRAJ_INTERNAL_VERSION "V4.18.0" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 94b9602c77..cce298300f 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -21,7 +21,7 @@ Action_Contacts.o : Action_Contacts.cpp Action.h ActionState.h Action_Contacts.h Action_CreateCrd.o : Action_CreateCrd.cpp Action.h ActionState.h Action_CreateCrd.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_CRD.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 StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h Action_CreateReservoir.o : Action_CreateReservoir.cpp Action.h ActionState.h Action_CreateReservoir.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_1D.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h NetcdfFile.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h RemdReservoirNC.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h Action_DNAionTracker.o : Action_DNAionTracker.cpp Action.h ActionState.h Action_DNAionTracker.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 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 SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h -Action_DSSP.o : Action_DSSP.cpp Action.h ActionState.h Action_DSSP.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 DistRoutines.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 SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h +Action_DSSP.o : Action_DSSP.cpp Action.h ActionState.h Action_DSSP.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 DistRoutines.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 SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h Action_Density.o : Action_Density.cpp Action.h ActionState.h Action_Density.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 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 OnlineVarT.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h Action_Diffusion.o : Action_Diffusion.cpp Action.h ActionState.h Action_Diffusion.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_1D.h DataSet_Coords.h DataSet_Coords_REF.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 SymbolExporting.h TextFormat.h Timer.h Topology.h Vec3.h Action_Dihedral.o : Action_Dihedral.cpp Action.h ActionState.h Action_Dihedral.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_double.h Dimension.h DispatchObject.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TorsionRoutines.h Vec3.h diff --git a/test/Test_DSSP/RunTest.sh b/test/Test_DSSP/RunTest.sh index 264748dd5d..71413a1e2a 100755 --- a/test/Test_DSSP/RunTest.sh +++ b/test/Test_DSSP/RunTest.sh @@ -5,7 +5,7 @@ # Clean CleanFiles cpptraj.in dssp.dat dssp.dat.sum dssp.sum.agr dssp.gnu \ dssp2.gnu dssp2.sum.agr total.agr assign.4.dat \ - DSSP.assign.dat + DSSP.assign.dat ftufabi.assign.dat INPUT="-i cpptraj.in" TOP="../DPDP.parm7" @@ -18,7 +18,7 @@ if [ $? -eq 0 ] ; then noprogress trajin ../DPDP.nc secstruct out dssp.gnu sumout dssp.sum.agr totalout total.agr \ - assignout DSSP.assign.dat + assignout DSSP.assign.dat betadetail EOF RunCpptraj "Secstruct (DSSP) command test." DoTest dssp.gnu.save dssp.gnu @@ -29,7 +29,7 @@ EOF # Test 2 cat > cpptraj.in < cpptraj.in < cpptraj.in <