secondaryStructures::secondaryStructures(){
}
-void secondaryStructures::initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez){
+void secondaryStructures::initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez, const int _pp_stretch){
SecondaryStructuresStatusMap.resize(0);
SecondaryStructuresStringLine.resize(0);
std::vector<std::size_t> temp; temp.resize(0),
PiHelixPreference = PiHelicesPreferencez;
ResInfoMap = &ResInfoMatrix;
+ pp_stretch = _pp_stretch;
SecondaryStructuresStatusMap.resize(ResInfoMatrix.size());
SecondaryStructuresStringLine.resize(ResInfoMatrix.size(), '~');
}
-void secondaryStructures::secondaryStructuresData::setStatus(const secondaryStructureTypes secondaryStructureTypeName){
- SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)] = true;
- SecondaryStructuresStatus = secondaryStructureTypeName;
+void secondaryStructures::secondaryStructuresData::setStatus(const secondaryStructureTypes secondaryStructureTypeName, bool status){
+ SecondaryStructuresStatusArray[static_cast<std::size_t>(secondaryStructureTypeName)] = status;
+ if (status){
+ SecondaryStructuresStatus = secondaryStructureTypeName;
+ }
+ else{
+ SecondaryStructuresStatus = secondaryStructureTypes::Loop;
+ }
}
void secondaryStructures::secondaryStructuresData::setStatus(const HelixPositions helixPosition, const turnsTypes turn){
for(const bridgeTypes &bridgeType : {bridgeTypes::ParallelBridge, bridgeTypes::AntiParallelBridge}){
if (SecondaryStructuresStatusMap[i].hasBridges(bridgeType) && SecondaryStructuresStatusMap[i + j].hasBridges(bridgeType) ){
std::size_t i_partner{SecondaryStructuresStatusMap[i].getBridgePartnerIndex(bridgeType)}, j_partner{SecondaryStructuresStatusMap[i + j].getBridgePartnerIndex(bridgeType)}, second_strand{};
-// std::cout << "i = " << i << std::endl;
-// std::cout << "j = " << j << std::endl;
-// std::cout << "i_partner = " << i_partner << std::endl;
-// std::cout << "j_partner = " << j_partner << std::endl;
if ( abs(i_partner - j_partner) < 6){
if (i_partner < j_partner){
second_strand = i_partner;
second_strand = j_partner;
}
for(int k{abs(i_partner - j_partner)}; k >= 0; --k){
+ if (SecondaryStructuresStatusMap[second_strand + k].getStatus(secondaryStructureTypes::Bridge)){
+ SecondaryStructuresStatusMap[second_strand + k].setStatus(secondaryStructureTypes::Bridge, false);
+ }
+
SecondaryStructuresStatusMap[second_strand + k].setStatus(secondaryStructureTypes::Strand);
- std::cout << second_strand + k << " is strand" << std::endl;
}
is_Estrand = true;
}
-
}
}
if (is_Estrand){
for(std::size_t k{0}; k <= j; ++k){
+ if (SecondaryStructuresStatusMap[i + k].getStatus(secondaryStructureTypes::Bridge)){
+ SecondaryStructuresStatusMap[i + k].setStatus(secondaryStructureTypes::Bridge, false);
+ }
SecondaryStructuresStatusMap[i + k].setStatus(secondaryStructureTypes::Strand);
- std::cout << i + k << " is strand" << std::endl;
}
break;
}
void secondaryStructures::analyzeTurnsAndHelicesPatterns(){
for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
std::size_t stride {static_cast<std::size_t>(i) + 3};
-// std::cout << "Testing Helix_" << stride << std::endl;
for(std::size_t j {0}; j + stride < SecondaryStructuresStatusMap.size(); ++j){
-// std::cout << "Testing " << j << " and " << j + stride << std::endl;
if ( hasHBondBetween(j + stride, j) && NoChainBreaksBetween(j, j + stride) ){
-// std::cout << j << " and " << j + stride << " has hbond!" << std::endl;
SecondaryStructuresStatusMap[j + stride].setStatus(HelixPositions::End, i);
for (std::size_t k {1}; k < stride; ++k){
}
-void secondaryStructures::analyzePPHelicesPatterns(){}
-
std::string secondaryStructures::patternSearch(){
analyzeBridgesAndStrandsPatterns();
analyzeTurnsAndHelicesPatterns();
-// analyzePPHelicesPatterns();
// for(std::size_t i {0}; i < ResInfoMap->size(); ++i){
// std::cout << (*ResInfoMap)[i].info->nr << " " << *((*ResInfoMap)[i].info->name) << std::endl;
}
}
+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");
+ }
+ }
+
+}
+
void DsspTool::calculateHBondEnergy(ResInfo& Donor,
ResInfo& Acceptor,
const t_trxframe& fr,
// std::cout << IndexMap[i].info->nr << " " << *(IndexMap[i].info->name) << std::endl;
// }
- PatternSearch.initiateSearch(IndexMap, initParams.PPHelices);
+ PatternSearch.initiateSearch(IndexMap, initParams.PPHelices, initParams.pp_stretch);
calculateBends(fr, pbc);
Storage.storageData(frnr, PatternSearch.patternSearch());
biod.pnpi.spb.ru / alexxy / gromacs-dssp.git / blobdiff