std::vector< double > betaAngles;
};
+struct colorAnglesPairs
+{
+ std::vector< std::pair<size_t, size_t> > a1, a2;
+ std::vector< std::pair<size_t, size_t> > 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 {
}
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;
// функция выделения индексов для бэталистов
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];
}
// вычисление внутренних углов в краске
-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);
#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) {
}
}
-//// определение близко ли к белку находится краска
+// определение близко ли к белку находится краска
//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) {
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;
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";
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();
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_;
// 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);
}
// считывание красок
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);
// разбор топологии для индексации бета листов
}
// подсчёт углов в красках
std::vector< colorLocalAngles > colorFormation;
- colorsAnglesEvaluation(trajectoryFrame, colorsNames, colorFormation);
+ colorsAnglesEvaluation(trajectoryFrame, colorsVectorsIndex, colorFormation);
// рассчёт положения относительно белка
/*
* формально можно совместить определение близости с белком и определение ближайших бэта-листов
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);
}
// расчёт угла и среднего угла с ближайшими бета листами
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;
}
biod.pnpi.spb.ru / alexxy / gromacs-colorvec.git / commitdiff