From 9cb93836fef2e1fad0d3649c8a536d564cda2b07 Mon Sep 17 00:00:00 2001 From: Anatoly Date: Tue, 27 Oct 2020 16:07:33 +0300 Subject: [PATCH] - implemented an optimization of pre-search names->index for the vectors - added color stacking search and its output --- src/colorvec.cpp | 192 +++++++++++++++++++++++++++++++++++------------ 1 file changed, 145 insertions(+), 47 deletions(-) diff --git a/src/colorvec.cpp b/src/colorvec.cpp index f883911..dc9b765 100644 --- a/src/colorvec.cpp +++ b/src/colorvec.cpp @@ -63,17 +63,23 @@ struct colorLocalAngles { std::vector< double > betaAngles; }; +struct colorAnglesPairs +{ + std::vector< std::pair > a1, a2; + std::vector< std::pair > b1, b2; +}; + // функция парсинга одной строки void parseBetaListDATLine(const std::string ¤tLine, std::vector< std::vector< unsigned int > > &localInputBL) { - size_t equalCount = 0; + size_t equalCount {0}; std::vector< unsigned int > a; a.resize(0); - for (size_t i = 0; i < currentLine.size(); ++i) { + for (size_t i {0}; i < currentLine.size(); ++i) { if (currentLine[i] == '=') { // подсчитываем число "пустых" символов ++equalCount; } else { size_t temp {i}; - for (size_t j = i + 1; j < currentLine.size(); ++j) { + for (size_t j {i + 1}; j < currentLine.size(); ++j) { if (currentLine[j] == 'B' || currentLine[j] == 'E') { ++temp; } else { @@ -82,7 +88,7 @@ void parseBetaListDATLine(const std::string ¤tLine, std::vector< std::vect } if (temp - i > 3) { localInputBL.push_back(a); - for (size_t j = i; j <= temp; ++j) { + for (size_t j {i}; j <= temp; ++j) { localInputBL.back().push_back(j - equalCount); } i = temp; @@ -114,22 +120,12 @@ void betaListDigestion(const std::string &inputFile, std::vector< std::vector< s // функция выделения индексов для бэталистов inline void aminoacidsIndexation(const std::vector< size_t > &inputIndex, const gmx::TopologyInformation &top, std::vector< std::vector< size_t > > &aminoacidsIndex) { aminoacidsIndex.resize(0); - for (size_t i = 0; i < inputIndex.size(); ++i) { + for (size_t i {0}; i < inputIndex.size(); ++i) { aminoacidsIndex.resize(std::max(aminoacidsIndex.size(), static_cast< size_t >(top.atoms()->atom[inputIndex[i]].resind + 1))); aminoacidsIndex[top.atoms()->atom[inputIndex[i]].resind].push_back(inputIndex[i]); } } -// функция поиска RVec в кадре по имени->индексу -gmx::RVec returnRVec(const std::vector< gmx::RVec > &frame, const std::vector< std::pair< std::string, size_t > > &colorsIndex, const std::string &toFind) { - for (auto &i : colorsIndex) { - if (i.first == toFind) { - return frame[i.second]; - } - } - throw "WRONG COLOR ATOM NAME TO EVALUATE VECTORS"; -} - // функция векторного произведения двух RVec'ов inline void RVecVecMultiply(gmx::RVec &n, const gmx::RVec &a, const gmx::RVec &b) { n[0] = a[1] * b[2] - a[2] * b[1]; @@ -143,27 +139,29 @@ inline double RVecAngle(const gmx::RVec &a, const gmx::RVec &b) { } // вычисление внутренних углов в краске -void colorsAnglesEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex, +void colorsAnglesEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< colorAnglesPairs > &inputPairs, std::vector< colorLocalAngles > &colorStruct) { - colorStruct.resize(colorsIndex.size()); + colorStruct.resize(inputPairs.size()); #pragma omp parallel for ordered schedule(dynamic) - for (size_t i = 0; i < colorsIndex.size(); ++i) { + for (size_t i = 0; i < inputPairs.size(); ++i) { colorStruct[i].betaAngles.resize(0); // "левый" блок - colorStruct[i].a1 = returnRVec(frame, colorsIndex[i], "CAF") - returnRVec(frame, colorsIndex[i], "CAJ") + returnRVec(frame, colorsIndex[i], "CAK") - returnRVec(frame, colorsIndex[i], "CAO"); - colorStruct[i].a1 /= colorStruct[i].a1.norm(); - colorStruct[i].b1 = returnRVec(frame, colorsIndex[i], "CAO") - returnRVec(frame, colorsIndex[i], "CAJ") + returnRVec(frame, colorsIndex[i], "CAM") - returnRVec(frame, colorsIndex[i], "CAH") + returnRVec(frame, colorsIndex[i], "CAK") - returnRVec(frame, colorsIndex[i], "CAF"); - colorStruct[i].b1 /= colorStruct[i].b1.norm(); + for (size_t j {0}; j < inputPairs[i].a1.size(); ++j) { + colorStruct[i].a1 += frame[inputPairs[i].a1[j].first] - frame[inputPairs[i].a1[j].second]; + } + for (size_t j {0}; j < inputPairs[i].b1.size(); ++j) { + colorStruct[i].b1 += frame[inputPairs[i].b1[j].first] - frame[inputPairs[i].b1[j].second]; + } RVecVecMultiply(colorStruct[i].n1, colorStruct[i].a1, colorStruct[i].b1); - colorStruct[i].n1 /= colorStruct[i].n1.norm(); // "правый" блок - colorStruct[i].a2 = returnRVec(frame, colorsIndex[i], "CAB") - returnRVec(frame, colorsIndex[i], "CBQ") + returnRVec(frame, colorsIndex[i], "CBP") - returnRVec(frame, colorsIndex[i], "CBL"); - colorStruct[i].a2 /= colorStruct[i].a2.norm(); - colorStruct[i].b2 = returnRVec(frame, colorsIndex[i], "CBL") - returnRVec(frame, colorsIndex[i], "CBQ") + returnRVec(frame, colorsIndex[i], "CBN") - returnRVec(frame, colorsIndex[i], "CBS") + returnRVec(frame, colorsIndex[i], "CBP") - returnRVec(frame, colorsIndex[i], "CAB"); - colorStruct[i].b2 /= colorStruct[i].b2.norm(); + for (size_t j {0}; j < inputPairs[i].a2.size(); ++j) { + colorStruct[i].a2 += frame[inputPairs[i].a2[j].first] - frame[inputPairs[i].a2[j].second]; + } + for (size_t j {0}; j < inputPairs[i].b2.size(); ++j) { + colorStruct[i].b2 += frame[inputPairs[i].b2[j].first] - frame[inputPairs[i].b2[j].second]; + } RVecVecMultiply(colorStruct[i].n2, colorStruct[i].a2, colorStruct[i].b2); colorStruct[i].n2 *= -1; // для "сонаправленности" векторов нормали - colorStruct[i].n2 /= colorStruct[i].n2.norm(); colorStruct[i].a12 = RVecAngle(colorStruct[i].a1, colorStruct[i].a2); colorStruct[i].b12 = RVecAngle(colorStruct[i].b1, colorStruct[i].b2); colorStruct[i].n12 = RVecAngle(colorStruct[i].n1, colorStruct[i].n2); @@ -171,6 +169,39 @@ void colorsAnglesEvaluation(const std::vector< gmx::RVec > &frame, const std::ve #pragma omp barrier } +// функция поиска индекса в краске по имени->индексу +size_t returnRVecIndex(const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex, const std::string &toFind) { + for (auto &i : colorsIndex) { + for (auto &j : i) { + if (j.first == toFind) { + return j.second; + } + } + } + throw "WRONG COLOR ATOM NAME TO EVALUATE VECTORS"; +} + +// заполнение пресетов для поиска средних векторов в краске +void colorAnglesPairsEvaluation(const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex, + std::vector< colorAnglesPairs > &output) { + for (size_t i {0}; i < output.size(); ++i) { + output[i].a1.resize(0); + output[i].a1.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CAF"), returnRVecIndex(colorsIndex, "CAJ"))); + output[i].a1.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CAK"), returnRVecIndex(colorsIndex, "CAO"))); + output[i].b1.resize(0); + output[i].b1.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CAO"), returnRVecIndex(colorsIndex, "CAJ"))); + output[i].b1.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CAM"), returnRVecIndex(colorsIndex, "CAH"))); + output[i].b1.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CAK"), returnRVecIndex(colorsIndex, "CAF"))); + output[i].a2.resize(0); + output[i].a2.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CAB"), returnRVecIndex(colorsIndex, "CBQ"))); + output[i].a2.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CBP"), returnRVecIndex(colorsIndex, "CBL"))); + output[i].b2.resize(0); + output[i].b2.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CBL"), returnRVecIndex(colorsIndex, "CBQ"))); + output[i].b2.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CBN"), returnRVecIndex(colorsIndex, "CBS"))); + output[i].b2.push_back(std::make_pair(returnRVecIndex(colorsIndex, "CBP"), returnRVecIndex(colorsIndex, "CAB"))); + } +} + // вычисление направляющего вектора в бэта-листе inline void betaListsRVecsEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< std::vector< unsigned int > > &inputBetaLists, std::vector< gmx::RVec > &temp, const std::vector< size_t > &inputCA) { @@ -183,7 +214,7 @@ inline void betaListsRVecsEvaluation(const std::vector< gmx::RVec > &frame, cons } } -//// определение близко ли к белку находится краска +// определение близко ли к белку находится краска //bool isNearPeptide(const t_trxframe &fr, const t_pbc *pbc, const std::vector< gmx::RVec > &inputFrame, // /*gmx::AnalysisNeighborhood &nbhood, */const std::vector< size_t > &inputIndex, // const std::vector< std::pair< std::string, size_t > > &inputColor) { @@ -260,9 +291,9 @@ inline void searchNearBetaLists(const t_pbc *inputPBC, outputList.resize(inputBLists.size(), false); gmx::RVec tempRVec; for (size_t i {0}; i < inputColor.size(); ++i) { - for (size_t j = 0; j < inputBLists.size(); ++j) { - for (size_t k = 0; k < inputBLists[j].size(); ++k) { - for (size_t m = 0; m < inputAminoacids[inputBLists[j][k]].size(); ++m) { + for (size_t j {0}; j < inputBLists.size(); ++j) { + for (size_t k {0}; k < inputBLists[j].size(); ++k) { + for (size_t m {0}; m < inputAminoacids[inputBLists[j][k]].size(); ++m) { pbc_dx(inputPBC, inputFrame[inputAminoacids[inputBLists[j][k]][m]], inputFrame[inputColor[i].second], tempRVec); if (tempRVec.norm() <= cutOff) { outputList[j] = true; @@ -284,11 +315,54 @@ inline void computeAnglesColorsVsBeta(const gmx::RVec &inputBetaRVec, colorLocal colorFormation.betaAngles.push_back(RVecAngle(inputBetaRVec, colorFormation.n2)); } +void colorsStackingSearch(const t_pbc *inputPBC, const std::vector< gmx::RVec > &inputFrame, + const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex, + const float radius, const float epsi, + std::vector< std::vector< int > > &output) { + std::vector< std::vector< int > > outputTemp; + outputTemp.resize(colorsIndex.size()); + for (size_t i {0}; i < colorsIndex.size(); ++i) { + outputTemp[i].resize(0); + outputTemp[i].resize(colorsIndex.size(), 0); + } + gmx::RVec temp; + #pragma omp parallel for ordered schedule(dynamic) + for (size_t i = 0; i < colorsIndex.size(); ++i) { + for (size_t j {i + 1}; j < colorsIndex.size(); ++j) { + for (size_t k1 {0}; k1 < colorsIndex[i].size(); ++k1) { + for (size_t k2 {0}; k2 < colorsIndex[j].size(); ++k2) { + pbc_dx(inputPBC, inputFrame[colorsIndex[i][k1].second], inputFrame[colorsIndex[j][k2].second], temp); + if (temp.norm() <= radius) { + ++outputTemp[i][j]; + ++outputTemp[j][i]; + break; + } + } + } + } + } + #pragma omp barrier + output.resize(0); + output.resize(colorsIndex.size()); + for (size_t i {0}; i < colorsIndex.size(); ++i) { + output[i].resize(0); + } + for (size_t i {0}; i < outputTemp.size(); ++i) { + for (size_t j {0}; j < outputTemp[i].size(); ++j) { + if (static_cast< float >(outputTemp[i][j]) >= static_cast< float >(colorsIndex.size()) * epsi) { + output[i].push_back(j); + output[j].push_back(i); + } + } + } +} + // функция записи в файл значений углов для кадра -void anglesFileDump(const int frameNum, const std::string &output, const std::vector< bool > &toPeptide, const std::vector< colorLocalAngles > &colorFormation) { +void anglesFileDump(const int frameNum, const std::string &output, const std::vector< bool > &toPeptide, + const std::vector< colorLocalAngles > &colorFormation, const std::vector< std::vector< int > > inputStack) { std::ofstream file(output, std::ofstream::app); file << "frame =" << std::setw(8) << frameNum << std::endl; - for (size_t i = 0; i < colorFormation.size(); ++i) { + for (size_t i {0}; i < colorFormation.size(); ++i) { file << "color #" << std::setw(3) << i; if (toPeptide[i]) { file << std::setw(4) << "yes"; @@ -306,20 +380,29 @@ void anglesFileDump(const int frameNum, const std::string &output, const std::ve std::vector< double > betaTemp; betaTemp.resize(0); betaTemp.resize(6, 0.); // magic number, meh - for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) { + for (size_t j {0}; j < colorFormation[i].betaAngles.size(); ++j) { betaTemp[j % 6] += colorFormation[i].betaAngles[j]; } - for (size_t j = 0; j < betaTemp.size(); ++j) { + for (size_t j {0}; j < betaTemp.size(); ++j) { betaTemp[j] /= static_cast< double >(temp); } file << std::setw(4) << temp; - for (size_t j = 0; j < betaTemp.size(); ++j) { + for (size_t j {0}; j < betaTemp.size(); ++j) { file << std::setw(8) << std::setprecision(2) << betaTemp[j]; } - for (size_t j = 0; j < colorFormation[i].betaAngles.size(); ++j) { + for (size_t j {0}; j < colorFormation[i].betaAngles.size(); ++j) { file << std::setw(8) << std::setprecision(2) << colorFormation[i].betaAngles[j]; } } + if (inputStack[i].size() == 0) { + file << std::setw(4) << "no"; + } else { + file << std::setw(4) << "yes"; + file << std::setw(4) << inputStack[i].size(); + for (size_t j {0}; j < inputStack[i].size(); ++j) { + file << std::setw(4) << inputStack[i][j]; + } + } file << std::endl; } file.close(); @@ -361,11 +444,14 @@ class colorVec : public gmx::TrajectoryAnalysisModule gmx::SelectionList sel_; std::string fnOut {"name"}; // selectable std::string fnBetaListsDat; // selectable - float effRad {0.8}; // selectable + float effRad {0.4}; // selectable + float stackRad {0.5}; // selectable + float stackCoef {0.66}; // selectable std::vector< size_t > index; std::vector< size_t > indexCA; std::vector< std::vector< size_t > > aminoacidsIndex; std::vector< std::vector< std::pair< std::string, size_t > > > colorsNames; + std::vector< colorAnglesPairs > colorsVectorsIndex; std::vector< std::vector< std::vector< unsigned int > > > betaLists; gmx::AnalysisNeighborhood nb_; @@ -413,7 +499,15 @@ colorVec::initOptions( gmx::IOptionsContainer *options, // Add option for effRad constant options->addOption(gmx::FloatOption("efRad") .store(&effRad) - .description("max distance from colors to peptide in nm to consider to be \"near\"")); + .description("max distance from colors to peptide in nm to consider to be \"near\" - default == 0.4")); + // Add option for effRad constant + options->addOption(gmx::FloatOption("stackRad") + .store(&stackRad) + .description("max distance from color to color in nm to consider to be \"stacked\" - default == 0.5")); + // Add option for effRad constant + options->addOption(gmx::FloatOption("stackCoef") + .store(&stackCoef) + .description("min % (1 == 100%) stacked mass of a colorto be considered \"stacked\" - default == 0.66")); // Control input settings settings->setFlags(gmx::TrajectoryAnalysisSettings::efNoUserPBC); settings->setFlag(gmx::TrajectoryAnalysisSettings::efUseTopX); @@ -437,11 +531,13 @@ colorVec::initAnalysis( const gmx::TrajectoryAnalysisSettings &settings, } // считывание красок colorsNames.resize(sel_.size() - 2); - for (size_t i = 2; i < sel_.size(); ++i) { + for (size_t i {2}; i < sel_.size(); ++i) { for (gmx::ArrayRef< const int >::iterator ai {sel_[i].atomIndices().begin()}; (ai < sel_[i].atomIndices().end()); ai++) { colorsNames[i - 2].push_back(std::make_pair(*(top.atoms()->atomname[*ai]), static_cast< size_t >(*ai))); } } + colorsVectorsIndex.resize(colorsNames.size()); + colorAnglesPairsEvaluation(colorsNames, colorsVectorsIndex); // разбор dat файла для создания бэта листов betaListDigestion(fnBetaListsDat, betaLists); // разбор топологии для индексации бета листов @@ -481,7 +577,7 @@ colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::Trajecto } // подсчёт углов в красках std::vector< colorLocalAngles > colorFormation; - colorsAnglesEvaluation(trajectoryFrame, colorsNames, colorFormation); + colorsAnglesEvaluation(trajectoryFrame, colorsVectorsIndex, colorFormation); // рассчёт положения относительно белка /* * формально можно совместить определение близости с белком и определение ближайших бэта-листов @@ -494,7 +590,7 @@ colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::Trajecto std::vector< bool > colorsToPeptide; colorsToPeptide.resize(0); colorsToPeptide.resize(colorsNames.size(), false); - for (size_t i = 0; i < colorsNames.size(); ++i) { + for (size_t i {0}; i < colorsNames.size(); ++i) { colorsToPeptide[i] = isNearPeptide(pbc, trajectoryFrame, index, colorsNames[i], effRad); } // расчёт угла и среднего угла с ближайшими бета листами @@ -505,20 +601,22 @@ colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::Trajecto colorsToBetaAngles.resize(0); std::vector< gmx::RVec > betaListsRVecs; betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs, indexCA); - for (size_t i = 0; i < colorsToPeptide.size(); ++i) { + for (size_t i {0}; i < colorsToPeptide.size(); ++i) { if (colorsToPeptide[i]) { searchNearBetaLists(pbc, trajectoryFrame, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex, effRad); } } - for (size_t i = 0; i < colorsToBeta.size(); ++i) { - for (size_t j = 0; j < colorsToBeta[i].size(); ++j) { + for (size_t i {0}; i < colorsToBeta.size(); ++i) { + for (size_t j {0}; j < colorsToBeta[i].size(); ++j) { if (colorsToBeta[i][j]) { computeAnglesColorsVsBeta(betaListsRVecs[j], colorFormation[i]); } } } + std::vector< std::vector< int > > colorStack; + colorsStackingSearch(pbc, trajectoryFrame, colorsNames, stackRad, stackCoef, colorStack); // вывод в файл(ы) информацию для фрейма - anglesFileDump(frnr, fnOut, colorsToPeptide, colorFormation); + anglesFileDump(frnr, fnOut, colorsToPeptide, colorFormation, colorStack); std::cout << "frame number = " << std::setw(8) << frnr << " trajectory time = " << std::setw(8) << fr.time << std::endl; } -- 2.22.0