Несовпадают углы. Добавил coutы чтобы найти причину
authorГорелов Сергей Васильевич <gorelov_sv@pnpi.nrcki.ru>
Mon, 3 Oct 2022 09:37:45 +0000 (12:37 +0300)
committerГорелов Сергей Васильевич <gorelov_sv@pnpi.nrcki.ru>
Mon, 3 Oct 2022 09:37:45 +0000 (12:37 +0300)
src/dssptools.cpp

index ef2e921c6810919368a3f3cc64b89c5882f24d04..77b7e32d4f1f4d2b95405ceacf48e62a7fa8e590 100644 (file)
@@ -769,9 +769,14 @@ float DsspTool::CalculateDihedralAngle(const int &A, const int &B, const int &C,
     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]);
 
+    std::cout << "u = " << u << std::endl;
+    std::cout << "v = " << v << std::endl;
+
     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);
+        std::cout << "new u = " << u << std::endl;
+        std::cout << "new v = " << v << std::endl;
         if (u != 0 or v != 0){
             result = std::atan2(v, u) * gmx::c_rad2Deg;
         }
@@ -785,13 +790,14 @@ void DsspTool::calculateDihedrals(const t_trxframe &fr, const t_pbc *pbc){
     const float phi_max = -75 + epsilon; // -46
     const float psi_min = 145 - epsilon; // 116
     const float psi_max = 145 + epsilon; // 176
-    std::vector<float>          phi, psi;
-    phi.resize(0);
-    psi.resize(0);
-    phi.resize(IndexMap.size(), 360);
-    psi.resize(IndexMap.size(), 360);
+    std::vector<float>          phi(IndexMap.size(), 360), psi(IndexMap.size(), 360);
+//    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)
+        std::cout << "For resi " << i << " :" << std::endl;
         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)),