std::vector< std::pair<size_t, size_t> > b1, b2;
};
+struct massChargePair {
+ double m{0.}, ch{0.};
+};
+
// функция парсинга одной строки
-void parseBetaListDATLine(const std::string ¤tLine, std::vector< std::vector< unsigned int > > &localInputBL) {
+void parseBetaListDATLine(const std::string ¤tLine,
+ std::vector< std::vector< unsigned int > > &localInputBL) {
size_t equalCount {0};
std::vector< unsigned int > a;
a.resize(0);
}
// функция нахождения бета-листов в структуре по файлу ДССП
-void betaListDigestion(const std::string &inputFile, std::vector< std::vector< std::vector< unsigned int > > > &inputBL) {
+void betaListDigestion(const std::string &inputFile,
+ std::vector< std::vector< std::vector< unsigned int > > > &inputBL) {
inputBL.resize(0);
std::ifstream file(inputFile);
std::string line;
}
// функция выделения индексов для бэталистов
-inline void aminoacidsIndexation(const std::vector< size_t > &inputIndex, const gmx::TopologyInformation &top, std::vector< std::vector< size_t > > &aminoacidsIndex) {
+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) {
aminoacidsIndex.resize(std::max(aminoacidsIndex.size(), static_cast< size_t >(top.atoms()->atom[inputIndex[i]].resind + 1)));
}
// вычисление внутренних углов в краске
-void colorsAnglesEvaluation(const std::vector< gmx::RVec > &frame, const std::vector< colorAnglesPairs > &inputPairs,
- std::vector< colorLocalAngles > &colorStruct) {
+void colorsAnglesEvaluation(const std::vector< gmx::RVec > &frame,
+ const std::vector< colorAnglesPairs > &inputPairs,
+ std::vector< colorLocalAngles > &colorStruct) {
+ colorStruct.resize(0);
colorStruct.resize(inputPairs.size());
#pragma omp parallel for ordered schedule(dynamic)
for (size_t i = 0; i < inputPairs.size(); ++i) {
}
// функция поиска индекса в краске по имени->индексу
-size_t returnRVecIndex(const std::vector< std::vector< std::pair< std::string, size_t > > > &colorsIndex, const std::string &toFind) {
+size_t returnRVecIndex(const 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;
- }
+ if (i.first == toFind) {
+ return i.second;
}
}
throw "WRONG COLOR ATOM NAME TO EVALUATE VECTORS";
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].a1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAF"), returnRVecIndex(colorsIndex[i], "CAJ")));
+ output[i].a1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAK"), returnRVecIndex(colorsIndex[i], "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].b1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAO"), returnRVecIndex(colorsIndex[i], "CAJ")));
+ output[i].b1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAM"), returnRVecIndex(colorsIndex[i], "CAH")));
+ output[i].b1.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAK"), returnRVecIndex(colorsIndex[i], "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].a2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CAB"), returnRVecIndex(colorsIndex[i], "CBQ")));
+ output[i].a2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CBP"), returnRVecIndex(colorsIndex[i], "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")));
+ output[i].b2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CBL"), returnRVecIndex(colorsIndex[i], "CBQ")));
+ output[i].b2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CBN"), returnRVecIndex(colorsIndex[i], "CBS")));
+ output[i].b2.push_back(std::make_pair(returnRVecIndex(colorsIndex[i], "CBP"), returnRVecIndex(colorsIndex[i], "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) {
+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) {
temp.resize(0);
gmx::RVec tempA;
for (const auto &i : inputBetaLists) {
}
}
-// определение близко ли к белку находится краска
-//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) {
-// /*gmx::AnalysisNeighborhoodSearch nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
-// gmx::AnalysisNeighborhoodPair pair;*/
-
-// for (size_t i = 0; i < inputColor.size(); ++i) {
-// std::cout << inputColor[i].first << " " << i << " / " << inputColor.size() << std::endl;
-// gmx::AnalysisNeighborhood nbhood;
-// nbhood.setCutoff(0.8);
-// gmx::AnalysisNeighborhoodSearch nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
-// gmx::AnalysisNeighborhoodPair pair;
-// gmx::AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(inputFrame[inputColor[i].second].as_vec());
-// int count1 {0};
-// while (pairSearch.findNextPair(&pair)) {
-// std::cout << ++count1 << " ";
-// for (size_t j = 0; j < inputIndex.size(); ++j) {
-// if (pair.refIndex() == inputIndex[j]) {
-// return true;
-// }
-// }
-// }
-// std::cout << inputColor[i].first << " atom done" << std::endl;
-// }
-// return false;
-//}
-
// определение близко ли к белку находится краска
bool isNearPeptide( const t_pbc *inputPBC,
- const std::vector< gmx::RVec > &inputFrame, const std::vector< size_t > &inputIndex,
+ const std::vector< gmx::RVec > &inputFrame,
+ const std::vector< size_t > &inputIndex,
const std::vector< std::pair< std::string, size_t > > &inputColor, const double cutOff) {
gmx::RVec tempRVec;
for (size_t i {0}; i < inputColor.size(); ++i) {
}
// поиск ближайших к краске бэта-листов
-//inline void searchNearBetaLists(const t_trxframe &fr, const t_pbc *pbc, const std::vector< gmx::RVec > &inputFrame,
-// gmx::AnalysisNeighborhood &nbhood,
-// const std::vector< std::vector< unsigned int > > &inputBLists,
-// const std::vector< std::pair< std::string, size_t > > &inputColor,
-// std::vector< bool > &outputList, const std::vector< std::vector< size_t > > &inputAminoacids) {
-// outputList.resize(0);
-// outputList.resize(inputBLists.size(), false);
-// gmx::AnalysisNeighborhoodSearch nbsearch = nbhood.initSearch(pbc, gmx::AnalysisNeighborhoodPositions(fr.x, fr.natoms));
-// gmx::AnalysisNeighborhoodPair pair;
-// for (const auto &i : inputColor) {
-// while (nbsearch.startPairSearch(inputFrame[i.second].as_vec()).findNextPair(&pair)) {
-// 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) {
-// if (pair.testIndex() == inputAminoacids[inputBLists[j][k]][m]) {
-// outputList[j] = true;
-// }
-// }
-// }
-// }
-// }
-// }
-//}
-
inline void searchNearBetaLists(const t_pbc *inputPBC,
const std::vector< gmx::RVec > &inputFrame, const std::vector< std::vector< unsigned int > > &inputBLists,
const std::vector< std::pair< std::string, size_t > > &inputColor, std::vector< bool > &outputList,
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(0);
outputTemp.resize(colorsIndex.size());
for (size_t i {0}; i < colorsIndex.size(); ++i) {
outputTemp[i].resize(0);
if (temp.norm() <= radius) {
++outputTemp[i][j];
++outputTemp[j][i];
- break;
}
}
}
}
}
+// запоминаем массу/заряды атомов красок
+void colorMassChargesEval( const std::vector< std::pair< std::string, size_t > > &color,
+ const gmx::TopologyInformation &top,
+ std::vector< std::pair< size_t, massChargePair > > &outputIndexMassChargePair) {
+ outputIndexMassChargePair.resize(0);
+ massChargePair temp;
+ for (size_t i {0}; i < color.size(); ++i) {
+ temp.m = top.atoms()->atom[color[i].second].m;
+ temp.ch = top.atoms()->atom[color[i].second].q;
+ outputIndexMassChargePair.push_back(std::make_pair(color[i].second, temp));
+ }
+}
+
+// вычисляем дипольный вектор краски
+gmx::RVec colorDipoleEval(const std::vector< std::pair< size_t, massChargePair > > &inputIndexMassChargePair,
+ const std::vector< gmx::RVec > &inputFrame) {
+ gmx::RVec temp {0., 0., 0.};
+ double temp2 {0.};
+ for (auto &i : inputIndexMassChargePair) {
+ temp += inputFrame[i.first] * static_cast< float >(i.second.m);
+ temp2 += i.second.m;
+ }
+ temp /= static_cast< float >(temp2);
+ return temp;
+}
+
// функция записи в файл значений углов для кадра
-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) {
+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,
+ const std::vector< std::vector< double > > &dipoleBetaAngles) {
std::ofstream file(output, std::ofstream::app);
file << "frame =" << std::setw(8) << frameNum << std::endl;
for (size_t i {0}; i < colorFormation.size(); ++i) {
file << std::setw(8) << std::setprecision(2) << colorFormation[i].b12;
file << std::setw(8) << std::setprecision(2) << colorFormation[i].n12;
if (colorFormation[i].betaAngles.size() == 0) {
- file << std::setw(4) << 0;
+ file << std::setw(4) << 0 << " [] ";
} else {
int temp = colorFormation[i].betaAngles.size() / 6;
std::vector< double > betaTemp;
for (size_t j {0}; j < betaTemp.size(); ++j) {
betaTemp[j] /= static_cast< double >(temp);
}
- file << std::setw(4) << temp;
+ file << std::setw(4) << temp << " [ ";
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) {
file << std::setw(8) << std::setprecision(2) << colorFormation[i].betaAngles[j];
}
+ file << " ] ";
}
if (inputStack[i].size() == 0) {
- file << std::setw(4) << "no";
+ file << std::setw(4) << "no []";
} else {
- file << std::setw(4) << "yes";
+ 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 << " ] ";
+ }
+ if (dipoleBetaAngles[i].size() == 0) {
+ file << "no []";
+ } else {
+ double tempDBangle {0.};
+ for (size_t j {0}; j < dipoleBetaAngles[i].size(); ++j) {
+ tempDBangle += dipoleBetaAngles[i][j];
+ }
+ file << std::setw(8) << std::setprecision(2) << tempDBangle / dipoleBetaAngles[i].size();
+ for (size_t j {0}; j < dipoleBetaAngles[i].size(); ++j) {
+ file << std::setw(8) << std::setprecision(2) << dipoleBetaAngles[i][j];
+ }
}
file << std::endl;
}
private:
- gmx::SelectionList sel_;
- std::string fnOut {"name"}; // selectable
- std::string fnBetaListsDat; // 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_;
+ gmx::SelectionList sel_;
+ std::string fnOut {"name"}; // selectable
+ std::string fnBetaListsDat; // 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< colorLocalAngles > colorsAngles;
+ std::vector< std::vector< std::pair< size_t, massChargePair > > > colorsMassCharges;
+ std::vector< std::vector< std::vector< unsigned int > > > betaLists;
+
+ std::vector< gmx::RVec > trajectoryFrame;
// Copy and assign disallowed by base.
};
{
// считывание индекса
index.resize(0);
- for (gmx::ArrayRef< const int >::iterator ai {sel_.front().atomIndices().begin()}; (ai < sel_.front().atomIndices().end()); ai++) {
+ for (gmx::ArrayRef< const int >::iterator ai {sel_.front().atomIndices().begin()}; (ai < sel_.front().atomIndices().end()); ++ai) {
index.push_back(static_cast< size_t >(*ai));
}
// считывание индекса
indexCA.resize(0);
- for (gmx::ArrayRef< const int >::iterator ai {sel_[1].atomIndices().begin()}; (ai < sel_[1].atomIndices().end()); ai++) {
+ for (gmx::ArrayRef< const int >::iterator ai {sel_[1].atomIndices().begin()}; (ai < sel_[1].atomIndices().end()); ++ai) {
indexCA.push_back(static_cast< size_t >(*ai));
}
// считывание красок
+ colorsNames.resize(0);
colorsNames.resize(sel_.size() - 2);
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++) {
+ 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(0);
colorsVectorsIndex.resize(colorsNames.size());
colorAnglesPairsEvaluation(colorsNames, colorsVectorsIndex);
// разбор dat файла для создания бэта листов
betaListDigestion(fnBetaListsDat, betaLists);
// разбор топологии для индексации бета листов
aminoacidsIndexation(index, top, aminoacidsIndex);
- // задание радиуса рассматриваемых соседей
- nb_.setCutoff(effRad);
+
+ // запоминание распределения зарядов в красках
+ colorsMassCharges.resize(0);
+ colorsMassCharges.resize(colorsNames.size());
+ for (size_t i {0}; i < colorsNames.size(); ++i) {
+ colorMassChargesEval(colorsNames[i], top, colorsMassCharges[i]);
+ }
+
// вывод базовых значений
std::cout << std::endl;
std::cout << "effective radius = " << std::setw(30) << effRad << std::endl;
std::cout << "Beta frames = " << std::setw(30) << betaLists.size() << std::endl;
std::cout << "residue # = " << std::setw(30) << aminoacidsIndex.size() << std::endl;
std::cout << "colors # = " << std::setw(30) << colorsNames.size() << std::endl;
-
- /*
- * формально можно сделать:
- * найти соотношение между именами для углов и индексными номерами
- * заполнить std::vector< std::vector< size_t > > nameToIndex;
- * и в последствии считать:
- * a1 = frame[nameToIndex[0][1] + ... - ... - ...
- * b1 a2 b2 по аналогии только для [@<-[1, 2, 3]]
- *
- */
}
void
colorVec::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, gmx::TrajectoryAnalysisModuleData *pdata)
{
- std::vector< gmx::RVec > trajectoryFrame;
trajectoryFrame.resize(0);
// считывания текущего фрейма траектории
for (size_t i {0}; i < fr.natoms; ++i) {
trajectoryFrame.push_back(fr.x[i]);
}
// подсчёт углов в красках
- std::vector< colorLocalAngles > colorFormation;
- colorsAnglesEvaluation(trajectoryFrame, colorsVectorsIndex, colorFormation);
+ colorsAngles.resize(0);
+ colorsAnglesEvaluation(trajectoryFrame, colorsVectorsIndex, colorsAngles);
// рассчёт положения относительно белка
- /*
- * формально можно совместить определение близости с белком и определение ближайших бэта-листов
- * для этого думаю нужно как-то соотнести или сделать соотношение между индексом и бэта-листами (т.е. хранить что бы за О(1) делать)
- *
- */
-
-
-
- std::vector< bool > colorsToPeptide;
+ std::vector< bool > colorsToPeptide;
colorsToPeptide.resize(0);
colorsToPeptide.resize(colorsNames.size(), false);
for (size_t i {0}; i < colorsNames.size(); ++i) {
std::vector< std::vector< double > > colorsToBetaAngles;
colorsToBetaAngles.resize(0);
std::vector< gmx::RVec > betaListsRVecs;
+ // подсчёт направляющей beta листа
betaListsRVecsEvaluation(trajectoryFrame, betaLists[frnr], betaListsRVecs, indexCA);
+ // поиск ближайших beta листов
for (size_t i {0}; i < colorsToPeptide.size(); ++i) {
if (colorsToPeptide[i]) {
searchNearBetaLists(pbc, trajectoryFrame, betaLists[frnr], colorsNames[i], colorsToBeta[i], aminoacidsIndex, effRad);
}
- }
+ }
+ // подсчёт углов между векторами краски и beta листами
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]);
+ computeAnglesColorsVsBeta(betaListsRVecs[j], colorsAngles[i]);
}
}
}
+ // поиск стакированных красок
std::vector< std::vector< int > > colorStack;
colorsStackingSearch(pbc, trajectoryFrame, colorsNames, stackRad, stackCoef, colorStack);
+ // рассчёт дипольных моментов для красок
+ std::vector< gmx::RVec > colorsDipole;
+ std::vector< std::vector< double > > dipoleBetaAngles;
+ colorsDipole.resize(0);
+ dipoleBetaAngles.resize(0);
+ dipoleBetaAngles.resize(colorsNames.size());
+ for (auto &i : colorsMassCharges) {
+ colorsDipole.push_back(colorDipoleEval(i, trajectoryFrame));
+ }
+ for (size_t i {0}; i < colorsToBeta.size(); ++i) {
+ dipoleBetaAngles[i].resize(0);
+ for (size_t j {0}; j < colorsToBeta[i].size(); ++j) {
+ if (colorsToBeta[i][j]) {
+ dipoleBetaAngles[i].push_back(RVecAngle(colorsDipole[i], betaListsRVecs[j]));
+ }
+ }
+ }
// вывод в файл(ы) информацию для фрейма
- anglesFileDump(frnr, fnOut, colorsToPeptide, colorFormation, colorStack);
+ anglesFileDump(frnr, fnOut, colorsToPeptide, colorsAngles, colorStack, dipoleBetaAngles);
std::cout << "frame number = " << std::setw(8) << frnr << " trajectory time = " << std::setw(8) << fr.time << std::endl;
}
biod.pnpi.spb.ru / alexxy / gromacs-colorvec.git / commitdiff