include/newfit.h

   1 #ifndef NEWFIT_H
   2 #define NEWFIT_H
   3
   4 #include "gromacs/math/vectypes.h"
   5 #include "vector"
   6
   7 using gmx::RVec;
   8
   9 namespace NewFit {
  10
  11     void center_frame(std::vector< RVec > &cf_frame, std::vector< int > cf_atoms) {
  12         if (cf_atoms.empty()) {
  13             cf_atoms.resize(0);
  14             for (int i = 0; i < cf_frame.size(); i++) {
  15                 cf_atoms.push_back(i);
  16             }
  17         }
  18         RVec cf_temp;
  19         cf_temp[0] = 0;
  20         cf_temp[1] = 0;
  21         cf_temp[2] = 0;
  22         for (int i = 0; i < cf_atoms.size(); i++) {
  23             cf_temp += cf_frame[cf_atoms[i]];
  24         }
  25         cf_temp[0] /= cf_atoms.size();
  26         cf_temp[1] /= cf_atoms.size();
  27         cf_temp[2] /= cf_atoms.size();
  28         for (int i = 0; i < cf_atoms.size(); i++) {
  29             cf_frame[cf_atoms[i]] -= cf_temp;
  30         }
  31     }
  32
  33     void implement_fit(std::vector< RVec > &if_frame, RVec if_step, double if_A, double if_B, double if_C) {
  34         RVec temp;
  35         for (int i = 0; i < if_frame.size(); i++) {
  36             if_frame[i] += if_step;
  37             temp = if_frame[i];
  38             if_frame[i][0] = cos(if_C) * cos(if_B) * temp[0] + (-sin(if_C)* cos(if_A) + sin(if_B) * cos(if_C) * sin(if_A)) * temp[1] + (sin(if_C) * sin(if_A) + sin(if_B) * cos(if_C) * cos(if_A)) * temp[2];
  39             if_frame[i][1] = sin(if_C) * cos(if_B) * temp[0] + (cos(if_C) * cos(if_A) + sin(if_C) * sin(if_B) * sin(if_A)) * temp[1] + (-cos(if_C) * sin(if_A) + sin(if_C) * sin(if_B) * cos(if_A)) * temp[2];
  40             if_frame[i][2] = -sin(if_B) * temp[0] + cos(if_B) * sin(if_A) * temp[1] + cos(if_B) * cos(if_A) * temp[2];
  41         }
  42     }
  43
  44     double F(double aix, double aiy, double aiz, double bix, double biy, double biz, double cA, double sA, double cB, double sB, double cC, double sC) {
  45         return  gmx::square(aix + biy*(cA*sC - cC*sA*sB) - biz*(sA*sC + cA*cC*sB) - cB*cC*bix) +
  46                 gmx::square(aiy - biy*(cA*cC + sA*sB*sC) + biz*(cC*sA - cA*sB*sC) - cB*sC*bix) +
  47                 gmx::square(aiz + sB*bix - cA*cB*biz - cB*sA*biy) ;
  48     }
  49
  50     double Fx(double temp1, double temp2, double temp3, double cB, double sB, double cC, double sC) {
  51         return 2*sB*temp3 -
  52                2*cB*cC*temp1 -
  53                2*cB*sC*temp2 ;
  54     }
  55
  56     double Fy(double temp1, double temp2, double temp3, double cA, double sA, double cB, double sB, double cC, double sC) {
  57         return 2*(cA*sC - cC*sA*sB)*temp1 -
  58                2*(cA*cC + sA*sB*sC)*temp2 -
  59                2*cB*sA*temp3 ;
  60     }
  61
  62     double Fz(double temp1, double temp2, double temp3, double cA, double sA, double cB, double sB, double cC, double sC) {
  63         return 2*(cC*sA - cA*sB*sC)*temp2 -
  64                2*(sA*sC + cA*cC*sB)*temp1 -
  65                2*cA*cB*temp3 ;
  66     }
  67
  68     double FA(double temp1, double temp2, double temp3, double biy_pl_y, double biz_pl_z, double cA, double sA, double cB, double sB, double cC, double sC) {
  69         return 2*(biy_pl_y*(cC*sA - cA*sB*sC) + biz_pl_z*(cA*cC + sA*sB*sC))*temp2 -
  70                2*(cA*cB*biy_pl_y - cB*sA*biz_pl_z)*temp3 -
  71                2*(biy_pl_y*(sA*sC + cA*cC*sB) + biz_pl_z*(cA*sC - cC*sA*sB))*temp1 ;
  72     }
  73
  74     double FB(double temp1, double temp2, double temp3, double bix_pl_x, double biy_pl_y, double biz_pl_z, double cA, double sA, double cB, double sB, double cC, double sC) {
  75         return 2*(cB*bix_pl_x + sA*sB*biy_pl_y + cA*sB*biz_pl_z)*temp3 -
  76                2*(cA*cB*cC*biz_pl_z - cC*sB*bix_pl_x + cB*cC*sA*biy_pl_y)*temp1 -
  77                2*(cA*cB*sC*biz_pl_z - sB*sC*bix_pl_x + cB*sA*sC*biy_pl_y)*temp2 ;
  78     }
  79
  80     double FC(double temp1, double temp2, double bix_pl_x, double biy_pl_y, double biz_pl_z, double cA, double sA, double cB, double sB, double cC, double sC) {
  81         return 2*(biy_pl_y*(cA*cC + sA*sB*sC) - biz_pl_z*(cC*sA - cA*sB*sC) + cB*sC*bix_pl_x)*temp1 -
  82                2*(biz_pl_z*(sA*sC + cA*cC*sB) - biy_pl_y*(cA*sC - cC*sA*sB) + cB*cC*bix_pl_x)*temp2 ;
  83     }
  84
  85     double dist(std::vector< RVec > alfa, std::vector< RVec > beta, std::vector< std::pair< int, int > > cros) {
  86         double distance = 0;
  87         for (int i = 0; i < cros.size(); i++) {
  88             distance += F(alfa[cros[i].first][0], alfa[cros[i].first][1], alfa[cros[i].first][2], beta[cros[i].second][0], beta[cros[i].second][1], beta[cros[i].second][2], 1, 0, 1, 0, 1, 0);
  89         }
  90         return distance;
  91     }
  92
  93     double dist_old_format(rvec *alfa, rvec *beta, std::vector< std::pair< int, int > > cros) {
  94         double distance = 0;
  95         for (int i = 0; i < cros.size(); i++) {
  96             distance += F(alfa[cros[i].first][0], alfa[cros[i].first][1], alfa[cros[i].first][2], beta[cros[i].second][0], beta[cros[i].second][1], beta[cros[i].second][2], 1, 0, 1, 0, 1, 0);
  97         }
  98         return distance;
  99     }
 100
 101     void new_fit(std::vector< RVec > alfa, std::vector< RVec > &beta, std::vector< std::pair< int, int > > cros, double prec) {
 102         double bx, by, bz;
 103         RVec step;
 104         double x1 = 0, y1 = 0, z1 = 0, A1 = 0, B1 = 0, C1 = 0;
 105         double x2 = 0, y2 = 0, z2 = 0, A2 = 0, B2 = 0, C2 = 0;
 106         double cA, sA, cB, sB, cC, sC, temp1, temp2, temp3;
 107         long double H1 = 0.0001, H2 = 0.0001, F1 = 0, F2 = 0, H1_temp, H2_temp, H1_const = 1, H2_const = M_PI/180;
 108         F1 = dist(alfa, beta, cros) / cros.size();
 109         bool flag = true;
 110         do {
 111             F2 = 0;
 112             for (int i = 0; i < cros.size(); i++) {
 113                 bx = beta[cros[i].second][0] + x1;
 114                 by = beta[cros[i].second][1] + y1;
 115                 bz = beta[cros[i].second][2] + z1;
 116                 cA = cos(A1);
 117                 sA = sin(A1);
 118                 cB = cos(B1);
 119                 sB = sin(B1);
 120                 cC = cos(C1);
 121                 sC = sin(C1);
 122                 temp1 = alfa[cros[i].first][0] + by*(cA*sC - cC*sA*sB) - bz*(sA*sC + cA*cC*sB) - cB*cC*bx;
 123                 temp2 = alfa[cros[i].first][1] - by*(cA*cC + sA*sB*sC) + bz*(cC*sA - cA*sB*sC) - cB*sC*bx;
 124                 temp3 = alfa[cros[i].first][2] + sB*bx - cA*cB*bz - cB*sA*by;
 125                 x2 += Fx(temp1, temp2, temp3,                     cB, sB, cC, sC);
 126                 y2 += Fy(temp1, temp2, temp3,             cA, sA, cB, sB, cC, sC);
 127                 z2 += Fz(temp1, temp2, temp3,             cA, sA, cB, sB, cC, sC);
 128                 A2 += FA(temp1, temp2, temp3,     by, bz, cA, sA, cB, sB, cC, sC);
 129                 B2 += FB(temp1, temp2, temp3, bx, by, bz, cA, sA, cB, sB, cC, sC);
 130                 C2 += FC(temp1, temp2,        bx, by, bz, cA, sA, cB, sB, cC, sC);
 131             }
 132             H1_temp = std::max(std::fabs(H1*x2), std::max(std::fabs(H1*y2), std::fabs(H1*z2)));
 133             if (H1_temp > H1_const) {
 134                 H1 /= 10;
 135             } else {
 136                 x2 = x1 - H1 * x2;
 137                 y2 = y1 - H1 * y2;
 138                 z2 = z1 - H1 * z2;
 139             }
 140             H2_temp = std::max(std::fabs(H2*A2), std::max(std::fabs(H2*B2), std::fabs(H2*C2)));
 141             if (H2_temp > H2_const) {
 142                 H2 /= 10;
 143             } else {
 144                 A2 = A1 - H2 * A2;
 145                 B2 = B1 - H2 * B2;
 146                 C2 = C1 - H2 * C2;
 147             }
 148             for (int i = 0; i < cros.size(); i++) {
 149                 cA = cos(A2);
 150                 sA = sin(A2);
 151                 cB = cos(B2);
 152                 sB = sin(B2);
 153                 cC = cos(C2);
 154                 sC = sin(C2);
 155                 F2 += F(alfa[cros[i].first][0], alfa[cros[i].first][1], alfa[cros[i].first][2], beta[cros[i].second][0] + x2, beta[cros[i].second][1] + y2, beta[cros[i].second][2] + z2, cA, sA, cB, sB, cC, sC);
 156             }
 157             F2 /= cros.size();
 158             if (F2 <= F1) {
 159                 x1 = x2;
 160                 y1 = y2;
 161                 z1 = z2;
 162                 A1 = A2;
 163                 B1 = B2;
 164                 C1 = C2;
 165                 x2 = 0;
 166                 y2 = 0;
 167                 z2 = 0;
 168                 A2 = 0;
 169                 B2 = 0;
 170                 C2 = 0;
 171                 if (std::fabs(F1 - F2) <= prec) {
 172                     flag = false;
 173                 }
 174                 F1 = F2;
 175             }
 176         } while (flag);
 177         step[0] = x1;
 178         step[1] = y1;
 179         step[2] = z1;
 180         implement_fit(beta, step, A1, B1, C1);
 181     }
 182 }
 183 #endif // NEWFIT_H