+float DsspTool::CalculateDihedralAngle(const int &A, const int &B, const int &C, const int &D, const t_trxframe &fr, const t_pbc *pbc){
+ float result{360}, u{}, v{};
+ gmx::RVec vBA{}, vDC{}, vBD{}, p{}, x{}, y{};
+ pbc_dx(pbc, fr.x[A], fr.x[B], vBA.as_vec());
+ pbc_dx(pbc, fr.x[C], fr.x[D], vDC.as_vec());
+ pbc_dx(pbc, fr.x[D], fr.x[B], vBD.as_vec());
+
+ for(std::size_t i{XX}, j{j + 1}, k{i + 2}; i <= ZZ; ++i, ++j, ++k){
+ if (j > 2){
+ j -= 3;
+ }
+ if (k > 2){
+ k -= 3;
+ }
+ p[i] = (vBD[j] * vBA[k]) - (vBD[k] * vBA[j]);
+ x[i] = (vBD[j] * vDC[k]) - (vBD[k] * vDC[j]);
+
+ }
+ for(std::size_t i{XX}, j{j + 1}, k{i + 2}; i <= ZZ; ++i, ++j, ++k){
+ if (j > 2){
+ j -= 3;
+ }
+ if (k > 2){
+ k -= 3;
+ }
+ y[i] = (vBD[j] * x[k]) - (vBD[k] * x[j]);
+ }
+
+ u = (x[XX] * x[XX]) + (x[YY] * x[YY]) + (x[ZZ] * x[ZZ]);
+ v = (y[XX] * y[XX]) + (y[YY] * y[YY]) + (y[ZZ] * y[ZZ]);
+
+ if (u > 0 and v > 0){
+ u = ((p[XX] * x[XX]) + (p[YY] * x[YY]) + (p[ZZ] * x[ZZ])) / std::sqrt(u);
+ v = ((p[XX] * y[XX]) + (p[YY] * y[YY]) + (p[ZZ] * y[ZZ])) / std::sqrt(v);
+ if (u != 0 or v != 0){
+ result = std::atan2(v, u) * gmx::c_rad2Deg;
+ }
+ }
+ return result;
+}
+
+void DsspTool::calculateDihedrals(const t_trxframe &fr, const t_pbc *pbc){
+ const float epsilon = 29;
+ const float phi_min = -75 - epsilon;
+ const float phi_max = -75 + epsilon;
+ const float psi_min = 145 - epsilon;
+ const float psi_max = 145 + epsilon;
+ std::vector<float> phi, psi;
+ phi.resize(0);
+ psi.resize(0);
+ phi.resize(IndexMap.size(), 360);
+ psi.resize(IndexMap.size(), 360);
+
+ for (std::size_t i = 1; i + 1 < IndexMap.size(); ++i){ // TODO add index verifictaion (check if those atom indexes exist)
+ phi[i] = CalculateDihedralAngle(static_cast<int>(IndexMap[i - 1].getIndex(backboneAtomTypes::AtomC)),
+ static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomN)),
+ static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomCA)),
+ static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
+ fr,
+ pbc);
+ psi[i] = CalculateDihedralAngle(static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomN)),
+ static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomCA)),
+ static_cast<int>(IndexMap[i].getIndex(backboneAtomTypes::AtomC)),
+ static_cast<int>(IndexMap[i + 1].getIndex(backboneAtomTypes::AtomN)),
+ fr,
+ pbc);
+ }
+
+ for (std::size_t i = 1; i + 3 < IndexMap.size(); ++i){
+ switch (initParams.pp_stretch){
+ case 2: {
+ if (phi_min > phi[i] or phi[i] > phi_max or
+ phi_min > phi[i + 1] or phi[i + 1]> phi_max){
+ continue;
+ }
+
+ if (psi_min > psi[i] or psi[i] > psi_max or
+ psi_min > psi[i + 1] or psi[i + 1] > psi_max){
+ continue;
+ }
+
+ switch (PatternSearch.SecondaryStructuresStatusMap[i].getStatus(turnsTypes::Turn_PP)){
+ case HelixPositions::None:
+ PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start, turnsTypes::Turn_PP);
+ break;
+
+ case HelixPositions::End:
+ PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start_AND_End, turnsTypes::Turn_PP);
+ break;
+
+ default:
+ break;
+ }
+
+ PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(HelixPositions::End, turnsTypes::Turn_PP);
+ /* Пропустил проверку того, что заменяемая ак - петля */
+ PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Helix_PP);
+ PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(secondaryStructureTypes::Helix_PP);
+ break;
+ }
+ case 3:{
+ if (phi_min > phi[i] or phi[i] > phi_max or
+ phi_min > phi[i + 1] or phi[i + 1]> phi_max or
+ phi_min > phi[i + 2] or phi[i + 2]> phi_max){
+ continue;
+ }
+
+ if (psi_min > psi[i] or psi[i] > psi_max or
+ psi_min > psi[i + 1] or psi[i + 1] > psi_max or
+ psi_min > psi[i + 2] or psi[i + 2] > psi_max){
+ continue;
+ }
+
+ switch (PatternSearch.SecondaryStructuresStatusMap[i].getStatus(turnsTypes::Turn_PP)){
+ case HelixPositions::None:
+ PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start, turnsTypes::Turn_PP);
+ break;
+
+ case HelixPositions::End:
+ PatternSearch.SecondaryStructuresStatusMap[i].setStatus(HelixPositions::Start_AND_End, turnsTypes::Turn_PP);
+ break;
+
+ default:
+ break;
+ }
+
+ PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(HelixPositions::Middle, turnsTypes::Turn_PP);
+ PatternSearch.SecondaryStructuresStatusMap[i + 2].setStatus(HelixPositions::End, turnsTypes::Turn_PP);
+ /* Пропустил проверку того, что заменяемая ак - петля */
+ PatternSearch.SecondaryStructuresStatusMap[i].setStatus(secondaryStructureTypes::Helix_PP);
+ PatternSearch.SecondaryStructuresStatusMap[i + 1].setStatus(secondaryStructureTypes::Helix_PP);
+ PatternSearch.SecondaryStructuresStatusMap[i + 2].setStatus(secondaryStructureTypes::Helix_PP);
+
+ break;
+ }
+ default:
+ throw std::runtime_error("Unsupported stretch length");
+ }
+ }
+
+}
+