* \ingroup module_trajectoryanalysis
*/
+
#include <iostream>
#include <algorithm>
+#include <vector>
+#include <math.h>
#include <omp.h>
-#include <set>
#include <gromacs/trajectoryanalysis.h>
#include <gromacs/utility/smalloc.h>
#include <gromacs/math/do_fit.h>
-
-#include <vector>
-#include <math.h>
-
-//#include "fittingn.cpp"
-
using namespace gmx;
using gmx::RVec;
-double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
- return sqrt( pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
- pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
- pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2) );
-}
-
-void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
- double t01, t02, t03, t04, t05, t06, t07;
- t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
- t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
- t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
- t04 = sqrt(t01 + t02 + t03);
- t05 = (aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x);
- t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
- t07 = (aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x);
-
- //#pragma omp parallel sections
- //{
- //#pragma omp section
- F += t04;
- //#pragma omp section
- Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
- //#pragma omp section
- Fy += -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) * t07 - 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) * t05 + 2 * cos(B) * sin(A) * t06) / (2 * t04);
- //#pragma omp section
- Fz += -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) * t05 - 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) * t07 + 2 * cos(A) * cos(B) * t06) / (2 * t04);
- //#pragma omp section
- Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
- 2 * (biy_plus_y * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + biz_plus_z * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) * t07 +
- 2 * (biy_plus_y * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + biz_plus_z * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) * t05) / (2 * t04);
- //#pragma omp section
- Fb += -(2 * (cos(A) * cos(B) * sin(C) * biz_plus_z - sin(B) * sin(C) * bix_plus_x + cos(B) * sin(A) * sin(C) * biy_plus_y) * t07 +
- 2 * (cos(A) * cos(B) * cos(C) * biz_plus_z - cos(C) * sin(B) * bix_plus_x + cos(B) * cos(C) * sin(A) * biy_plus_y) * t05 -
- 2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
- //#pragma omp section
- Fc += (2 * (biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * bix_plus_x) * t05 -
- 2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
- //}
-}
-
-void searchL (double &fLd, double &fLdd, double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C, double L, double fx, double fy, double fz, double fa, double fb, double fc) {
- double AmLfa, BmLfb, CmLfc, BIXpXmLfx, BIYpYmLfy, BIZpZmLfz, t1, t2, t3, t4, t5, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16, t17, t18, t19, t20;
- AmLfa = A - L * fa;
- BmLfb = B - L * fb;
- CmLfc = C - L * fc;
- BIXpXmLfx = bix + x - L * fx;
- BIYpYmLfy = biy + y - L * fy;
- BIZpZmLfz = biz + z - L * fz;
- t1 = (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc));
- t2 = (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc));
- t3 = (cos(AmLfa) * sin(CmLfc) - cos(CmLfc) * sin(AmLfa) * sin(BmLfb));
- t4 = (sin(AmLfa) * sin(CmLfc) + cos(AmLfa) * cos(CmLfc) * sin(BmLfb));
- t5 = pow(aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy, 2);
- t6 = pow(aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx, 2);
- t7 = pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2);
-
- t8 = fc * sin(AmLfa) * sin(CmLfc);
- t9 = fa * cos(AmLfa) * cos(CmLfc);
- t10 = fb * cos(AmLfa) * cos(BmLfb) * sin(CmLfc);
- t11 = fc * cos(AmLfa) * cos(CmLfc) * sin(BmLfb);
- t12 = fa * sin(AmLfa) * sin(BmLfb) * sin(CmLfc);
-
- t13 = fa * cos(AmLfa) * sin(BmLfb) * sin(CmLfc);
- t14 = fc * cos(AmLfa) * sin(CmLfc);
- t15 = fa * cos(CmLfc) * sin(AmLfa);
- t16 = fb * cos(BmLfb) * sin(AmLfa) * sin(CmLfc);
- t17 = fc * cos(CmLfc) * sin(AmLfa) * sin(BmLfb);
-
- t18 = fx * cos(BmLfb) * sin(CmLfc);
- t19 = fc * cos(BmLfb) * cos(CmLfc);
- t20 = fb * sin(BmLfb) * sin(CmLfc);
-
- fLd += (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- ( BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
- 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) -
- fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy + fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
- 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
- fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * t3 +
- fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) /
- (2 * sqrt( t5 + t6 + pow(aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx, 2)));
-
-
- fLdd += ( (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (2 * fb * fx * cos(BmLfb) - pow(fb, 2) * sin(BmLfb) * BIXpXmLfx - 2 * fa * fy * cos(AmLfa) * cos(BmLfb) + 2 * fa * fz * cos(BmLfb) * sin(AmLfa) + 2 * fb * fz * cos(AmLfa) * sin(BmLfb) +
- 2 * fb * fy * sin(AmLfa) * sin(BmLfb) + pow(fa, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz + pow(fb, 2) * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz +
- pow(fa, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + pow(fb, 2) * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy + 2 * fa * fb * cos(AmLfa) * sin(BmLfb) * BIYpYmLfy -
- 2 * fa * fb * sin(AmLfa) * sin(BmLfb) * BIZpZmLfz) +
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
- fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) +
- fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fy * (cos(AmLfa) * sin(CmLfc) -
- cos(CmLfc) * sin(AmLfa) * sin(BmLfb)) + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) -
- fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (2 * BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) +
- fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) -
- fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) +
- t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- 2 * fy * t3 + 2 * fz * t4 +
- 2 * fx * cos(BmLfb) * cos(CmLfc) - 2 * fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - 2 * fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
- (2 * aix + 2 * t3 * BIYpYmLfy - 2 * (sin(AmLfa) * sin(CmLfc) +
- cos(AmLfa) * cos(CmLfc) * sin(BmLfb)) * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (pow(fa, 2) * sin(AmLfa) * sin(CmLfc) +
- pow(fc, 2) * sin(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(AmLfa) * cos(CmLfc) + pow(fa, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) +
- pow(fb, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) * sin(BmLfb) + 2 * fa * fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) +
- 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * sin(CmLfc) - 2 * fa * fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + BIYpYmLfy * (pow(fa, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) -
- pow(fc, 2) * cos(AmLfa) * sin(CmLfc) - 2 * fa * fc * cos(CmLfc) * sin(AmLfa) - pow(fa, 2) * cos(AmLfa) * sin(CmLfc) + pow(fb, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) +
- pow(fc, 2) * cos(CmLfc) * sin(AmLfa) * sin(BmLfb) - 2 * fa * fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
- 2 * fb * fc * cos(BmLfb) * sin(AmLfa) * sin(CmLfc)) - 2 * fz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) +
- fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- 2 * fy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) -
- fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) + 2 * fb * fx * cos(CmLfc) * sin(BmLfb) + 2 * fc * t18 + pow(fb, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx +
- pow(fc, 2) * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx - 2 * fb * fc * sin(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
- (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) *
- (2 * BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + 2 * fy * t1 - 2 * fz * t2 + 2 * BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + 2 * t18 + 2 * t19 * BIXpXmLfx - 2 * t20 * BIXpXmLfx) +
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) *
- (2 * fx * sin(BmLfb) - 2 * fy * cos(BmLfb) * sin(AmLfa) + 2 * fb * cos(BmLfb) * BIXpXmLfx - 2 * fz * cos(AmLfa) * cos(BmLfb) - 2 * fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- 2 * fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + 2 * fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + 2 * fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
- (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * (cos(CmLfc) * sin(AmLfa) - cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (pow(fa, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) -
- pow(fc, 2) * cos(CmLfc) * sin(AmLfa) - 2 * fa * t14 - pow(fa, 2) * cos(CmLfc) * sin(AmLfa) + pow(fb, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) +
- pow(fc, 2) * cos(AmLfa) * sin(BmLfb) * sin(CmLfc) - 2 * fb * fc * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) + 2 * fa * t16 +
- 2 * fa * t17) + BIYpYmLfy * (pow(fa, 2) * cos(AmLfa) * cos(CmLfc) + pow(fc, 2) * cos(AmLfa) * cos(CmLfc) - 2 * fa * t8 +
- pow(fa, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fb, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) + pow(fc, 2) * sin(AmLfa) * sin(BmLfb) * sin(CmLfc) -
- 2 * fa * t10 - 2 * fa * t11 - 2 * fb * t19 * sin(AmLfa)) - 2 * fz * (t8 - t9 + t10 + t11 - t12) - 2 * fy * (t13 - t14 - t15 + t16 + t17) -
- 2 * fc * fx * cos(BmLfb) * cos(CmLfc) + 2 * fb * fx * sin(BmLfb) * sin(CmLfc) +
- pow(fb, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + pow(fc, 2) * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx + 2 * fb * fc * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx)) /
- (2 * sqrt( t5 + t6 + t7)) -
- ( (2 * (aiy - t1 * BIYpYmLfy + t2 * BIZpZmLfz - cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * t1 - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) + t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) -
- 2 * (aiz + sin(BmLfb) * BIXpXmLfx - cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy) +
- 2 * (aix + t3 * BIYpYmLfy - t4 * BIZpZmLfz - cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx)) *
- ( (2 * aix + 2 * t3 * BIYpYmLfy - 2 * t4 * BIZpZmLfz - 2 * cos(BmLfb) * cos(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (fa * cos(AmLfa) * sin(CmLfc) + fc * cos(CmLfc) * sin(AmLfa) + fb * cos(AmLfa) * cos(BmLfb) * cos(CmLfc) - t15 * sin(BmLfb) - fc * cos(AmLfa) * sin(BmLfb) * sin(CmLfc)) +
- BIYpYmLfy * (fa * sin(AmLfa) * sin(CmLfc) - fc * cos(AmLfa) * cos(CmLfc) + t9 * sin(BmLfb) + fb * cos(BmLfb) * cos(CmLfc) * sin(AmLfa) - fc * sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) -
- fy * t3 + fz * t4 + fx * cos(BmLfb) * cos(CmLfc) - fb * cos(CmLfc) * sin(BmLfb) * BIXpXmLfx - fc * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) +
- (2 * aiy - 2 * t1 * BIYpYmLfy + 2 * t2 * BIZpZmLfz - 2 * cos(BmLfb) * sin(CmLfc) * BIXpXmLfx) *
- (BIZpZmLfz * (t8 - t9 + t10 + t11 - t12) + fy * (cos(AmLfa) * cos(CmLfc) + sin(AmLfa) * sin(BmLfb) * sin(CmLfc)) - fz * t2 + BIYpYmLfy * (t13 - t14 - t15 + t16 + t17) +
- t18 + t19 * BIXpXmLfx - t20 * BIXpXmLfx) - (2 * aiz + 2 * sin(BmLfb) * BIXpXmLfx - 2 * cos(AmLfa) * cos(BmLfb) * BIZpZmLfz - 2 * cos(BmLfb) * sin(AmLfa) * BIYpYmLfy) *
- (fx * sin(BmLfb) - fy * cos(BmLfb) * sin(AmLfa) + fb * cos(BmLfb) * BIXpXmLfx - fz * cos(AmLfa) * cos(BmLfb) - fa * cos(AmLfa) * cos(BmLfb) * BIYpYmLfy +
- fa * cos(BmLfb) * sin(AmLfa) * BIZpZmLfz + fb * cos(AmLfa) * sin(BmLfb) * BIZpZmLfz + fb * sin(AmLfa) * sin(BmLfb) * BIYpYmLfy))) /
- (4 * pow ( t5 + t6 + t7, 3 / 2));
-
-}
-
-
-void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
- RVec temp;
- //#pragma omp parallel
- //{
- //#pragma omp for schedule(dynamic)
- for (int i = 0; i < b.size(); i++) {
- temp = b[i];
- b[i][0] = (temp[2] + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (temp[1] + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (temp[0] + x);
- b[i][1] = (temp[1] + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (temp[2] + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (temp[0] + x);
- b[i][2] = cos(A) * cos(B) * (temp[2] + z) - sin(B) * (temp[0] + x) + cos(B) * sin(A) * (temp[1] + y);
- }
- //}
-}
-
-void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs) {
- midA[0] = 0;
- midA[1] = 0;
- midA[2] = 0;
-
- midB[0] = 0;
- midB[1] = 0;
- midB[2] = 0;
-
- for (int i = 0; i < pairs.size(); i++) {
- rvec_inc(midA, a[pairs[i].first]);
- rvec_inc(midB, b[pairs[i].second]);
- }
- midA[0] /= pairs.size();
- midA[1] /= pairs.size();
- midA[2] /= pairs.size();
-
- midB[0] /= pairs.size();
- midB[1] /= pairs.size();
- midB[2] /= pairs.size();
-}
-
-void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs) {
- double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
- RVec ma, mb;
- CalcMid(a, b, ma, mb, pairs);
- rvec_dec(ma, mb);
- double x = ma[0], y = ma[1], z = ma[2], A = 0, B = 0, C = 0;
- for (int i = 0; i < pairs.size(); i++) {
- f1 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
- }
- if (f1 == 0) {
- return;
- } else {
- double fLd, fLdd;
- //l = 0;
- int count = 0;
- while (f1 - f2 < 0 || f1 - f2 > 0.00001) {
- f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
- for (int i = 0; i < pairs.size(); i++) {
- searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + x, b[pairs[i].second][1] + y, b[pairs[i].second][2] + z, A, B, C);
- }
- /*do {
- fLd = 0;
- fLdd = 0;
- for (int i = 0; i < pairs.size(); i++) {
- searchL(fLd, fLdd, a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0], b[pairs[i].second][1], b[pairs[i].second][2], x, y, z, A, B, C, l, fx, fy, fz, fa, fb, fc);
- }
- l -= fLd / fLdd;
- } while (std::abs(fLd) > 1);*/
- //std::cout << " Ltempo= " << Ltempo << " iter= " << Lco << "\n";
- while (true) {
- f2 = 0;
- for (int i = 0; i < pairs.size(); i++) {
- f2 += F(a[pairs[i].first][0], a[pairs[i].first][1], a[pairs[i].first][2], b[pairs[i].second][0] + (x - l * fx), b[pairs[i].second][1] + (y - l * fy), b[pairs[i].second][2] + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
- }
- count++;
- if (f2 < f1) {
- x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
- fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
- break;
- } else {
- l *= 0.85;
- }
- }
- count++;
- }
- //std::cout << "iteration count: " << count << "\n";
- ApplyFit(b, x, y, z, A, B, C);
- }
-}
-
-void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const char* file_name)
-{
- FILE *file;
- file = std::fopen(file_name, "w+");
- for (int i = 1; i < trajectory.size(); i++) {
- std::fprintf(file, "MODEL%9d\n", i);
- for (int j = 0; j < trajectory[i].size(); j++) {
- std::fprintf(file, "ATOM %5d CA LYS 1 %8.3f%8.3f%8.3f 1.00 0.00\n", j, trajectory[i][j][0] * 10, trajectory[i][j][1] * 10, trajectory[i][j][2] * 10);
- }
- std::fprintf(file, "ENDMDL\n");
- }
-}
-
void make_correlation_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start)
{
FILE *file;
}
for (int j = 0; j < correlations[i].size(); j++) {
for (int f = 0; f < correlations[i][j].size(); f++) {
- std::fprintf(file, "%3.2f ", correlations[i][j][f]);
+ std::fprintf(file, "%3.4f ", correlations[i][j][f]);
}
std::fprintf(file, "\n");
}
std::fprintf(file, "correlations >= %0.2f\n\n", crl_border);
for (int i = 0; i < rout.size(); i++) {
for (int j = 0; j < rout[i].size(); j++) {
- std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.15\n", indx[rout[i][j].first] + 1, indx[rout[i][j].second] + 1);
+ std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.15\n", indx[rout[i][j].first]/* + 1*/, indx[rout[i][j].second]/* + 1*/);
}
std::fprintf(file, "\n\n");
}
file = std::fopen(file_name, "w+");
for (int i = 0; i < rout_pairs.size(); i++) {
for (int j = 0; j < rout_pairs[i].size(); j++) {
- std::fprintf(file, "%3d %3d\n", indx[rout_pairs[i][j].first] + 1, indx[rout_pairs[i][j].second] + 1);
+ std::fprintf(file, "%3d %3d\n", indx[rout_pairs[i][j].first]/* + 1*/, indx[rout_pairs[i][j].second]/* + 1*/);
for (int k = 0; k < correlations.size(); k++) {
std::fprintf(file, "%3.5f ", correlations[k][rout_pairs[i][j].first][rout_pairs[i][j].second]);
}
std::fclose(file);
}
-bool mysortfunc (std::vector< int > a, std::vector< int > b) { return (a.size() > b.size()); }
+bool mysortfunc (std::vector< int > a, std::vector< int > b) {
+ return (a.size() > b.size());
+}
bool isitsubset (std::vector< int > a, std::vector< int > b) {
if (b.size() == 0) {
#pragma omp barrier
}
-void domain_chopping(SelectionList seList, std::vector< int > indx, std::vector< std::vector< int > > &filled_domains, std::vector< RVec > frame) {
- //ConstArrayRef< int > atomInd = seList[0].atomIndices();
- ArrayRef< const int > atomInd = seList[0].atomIndices();
- std::vector< std::vector< int > > init_domains;
- init_domains.resize(seList.size());
- for (int i = 0; i < seList.size(); i++) {
- if (seList.size() != 1 && i == 0) {
- continue;
- }
- atomInd = seList[i].atomIndices();
- /*for (ConstArrayRef<int>::iterator ai = atomInd.begin(); (ai < atomInd.end()); ai++) {
- init_domains[i].push_back(*ai);
- }*/
- for (ArrayRef< const int >::iterator ai = atomInd.begin(); (ai < atomInd.end()); ai++) {
- init_domains[i].push_back(*ai);
- }
- }
- for (int i = 0; i < init_domains.size(); i++) {
- for (int j = 0; j < init_domains[i].size(); j++) {
- int k = 0;
- while (indx[k] != init_domains[i][j]) {
- k++;
- }
- init_domains[i][j] = k;
- }
- }
- for (int i = 0; i < init_domains.size(); i++) {
- if (init_domains[i].size() >= 2) {
- filled_domains.push_back(init_domains[i]);
- for (int k = 0; k < init_domains[i].size(); k++) {
- for (int j = i + 1; j < init_domains.size(); j++) {
- for (int x = 0; x < init_domains[j].size(); x++) {
- if (init_domains[j][x] == init_domains[i][k]) {
- init_domains[j].erase(init_domains[j].begin() + x);
- }
- }
- }
- }
- }
- }
- init_domains.resize(0);
- init_domains = filled_domains;
- std::vector< bool > flags;
- flags.resize(indx.size(), true);
- for (int i = 0; i < init_domains.size(); i++) {
- for (int j = 0; j < init_domains[i].size(); j++) {
- flags[init_domains[i][j]] = false;
- }
- }
- int a;
- rvec temp;
- for (int i = 0; i < indx.size(); i++) {
- if (flags[i]) {
- float dist = 90000001;
- for (int j = 0; j < init_domains.size(); j++) {
- for (int k = 0; k < init_domains[j].size(); k++) {
- rvec_sub(frame[i], frame[init_domains[j][k]], temp);
- if (norm(temp) <= dist) {
- dist = norm(temp);
- a = j;
- }
- }
- }
- filled_domains[a].push_back(i);
- flags[i] = false;
- }
- }
-}
-
void graph_calculation(std::vector< std::vector< std::pair< float, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,
std::vector< std::vector< RVec > > traj, int b_frame,
std::vector< std::vector< std::vector< float > > > crl, float crl_b, float e_rad, int tauE) {
int frames = 0;
int basic_frame = 0;
int tau = 0;
- float crl_border = 0;
- float eff_rad = 1.5;
+ float crl_border = 0;
+ float eff_rad = 1.5;
real **w_rls;
std::vector< std::vector< std::pair< int, int > > > w_rls2;
-
- int domains_ePBC;
// Copy and assign disallowed by base.
};
Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top)
{
- domains_ePBC = top.ePBC();
}
void
Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
const t_trxframe &fr)
{
- //std::vector< std::vector< int > > domains;
- //ConstArrayRef< int > atomind = sel_[0].atomIndices();
ArrayRef< const int > atomind = sel_[0].atomIndices();
index.resize(0);
- /*for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
- index.push_back(*ai);
- }*/
for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
index.push_back(*ai);
}
for (int i = 0; i < index.size(); i++) {
trajectory.back()[i] = fr.x[index[i]];
}
-
- domain_chopping(sel_, index, domains_local, trajectory.back());
-
- frankenstein_trajectory.resize(trajectory.size());
- frankenstein_trajectory.back() = trajectory.back();
- trajectory.resize(trajectory.size() + 1);
- trajectory.back().resize(0);
- trajectory.back().resize(index.size());
-
- w_rls2.resize(domains_local.size() + 1);
- for (int i = 0; i < domains_local.size(); i++) {
- w_rls2[i].resize(domains_local[i].size());
- for (int j = 0; j < domains_local[i].size(); j++) {
- w_rls2[i][j] = std::make_pair(domains_local[i][j], domains_local[i][j]);
- }
- }
- w_rls2[domains_local.size()].resize(index.size());
- for (int i = 0; i < index.size(); i++) {
- w_rls2.back()[i] = std::make_pair(i, i);
- }
frames++;
}
// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test6.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionList5' -tau 5 -crl 0.10
// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/CorrsTestDomainsNfit.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionListDomainsNFit' -tau 5000 -crl 0.75 -ef_rad 1
// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/TestCa.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelListCa' -tau 100 -crl 0.75 -ef_rad 1
+// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/testCa.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelListCa' -tau 9000 -crl 0.60 -ef_rad 1
+
+// -s '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest1.pdb' -f '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest.pdb' -n '/home/toluk/Develop/samples/CubeTermal/testCa.ndx' -sf '/home/toluk/Develop/samples/CubeTermal/SelListCa' -tau 900 -crl 0.20 -ef_rad 1
void
Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
TrajectoryAnalysisModuleData *pdata)
{
+ trajectory.resize(trajectory.size() + 1);
+ trajectory.back().resize(index.size());
for (int i = 0; i < index.size(); i++) {
trajectory.back()[i] = fr.x[index[i]];
}
- frankenstein_trajectory.resize(frankenstein_trajectory.size() + 1);
- frankenstein_trajectory.back().resize(index.size());
- #pragma omp parallel shared(frankenstein_trajectory, trajectory)
- {
- #pragma omp for schedule(dynamic)
- for (int i = 0; i < domains_local.size(); i++) {
- std::vector< RVec > traj = trajectory.back();
- MyFitNew(trajectory[basic_frame], traj, w_rls2[i]);
- for (int j = 0; j < domains_local[i].size(); j++) {
- frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]];
- }
- }
- }
- #pragma omp barrier
- std::cout << "frame № " << frames + 1 <<" analysed\n";
frames++;
}
Domains::finishAnalysis(int nframes)
{
std::vector< std::vector< std::vector< float > > > crltns;
+ std::vector< std::vector< std::pair< float, int > > > graph;
+ std::vector< std::vector< int > > sub_graph;
+ std::vector< std::vector< std::pair< int, int > > > sub_graph_rbr, rout_new;
int k = 1000, m = 0;
if (tau > 0) {
k = tau;
}
- /*
- *
- *
- *
- */
+ std::cout << "\nCorrelation's evaluation - start\n";
- std::cout << "Correlation's evaluation - start\n\n";
correlation_evaluation(frankenstein_trajectory, basic_frame, crltns, m, k);
- make_correlation_file(crltns, "full crltns file Ca.txt", 0);
- std::cout << "Correlation's evaluation - end\n\n";
-
- /*
- *
- *
- *
- */
-
- //number of corelations in the matrixes
- /*for (int i1 = 0; i1 < 100; i1++) {
- int i5 = 0;
- for (int i2 = 1; i2 < crltns.size(); i2++) {
- for (int i3 = 0; i3 < crltns[i2].size(); i3++) {
- for (int i4 = 0; i4 < crltns[i2][i3].size(); i4++) {
- if (std::abs(crltns[i2][i3][i4]) >= (float)i1 / 100 && i3 != i4) {
- i5++;
- }
- }
- }
- }
- std::cout << i1 << " - " << i5 << " | ";
- if (i1 % 10 == 0) {
- std::cout << "\n" ;
- }
- }*/
-
- std::cout << "graph evaluation: start\n";
-
- std::vector< std::vector< std::pair< float, int > > > graph;
- std::vector< std::vector< int > > sub_graph;
- std::vector< std::vector< std::pair< int, int > > > sub_graph_rbr;
+ std::cout << "Correlation's evaluation - end\n" << "graph evaluation: start\n";
graph_calculation(graph, sub_graph, sub_graph_rbr, frankenstein_trajectory, basic_frame, crltns, crl_border, eff_rad, k);
-
- std::vector< std::vector < std::pair< int, int > > > rout_new;
-
- std::cout << "graph evaluation: end\n";
- std::cout << "routs evaluation: start\n";
+ std::cout << "graph evaluation: end\n" << "routs evaluation: start\n";
graph_back_bone_evaluation(rout_new, index.size(), graph, sub_graph, sub_graph_rbr);
std::cout << "routs evaluation: end\n";
-
- // Routs_DomainsNFit_5000_0.75_1_t.txt
- make_rout_file2(crl_border, index, rout_new, "Routs_Ca.txt");
- make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_Ca.txt");
+ make_correlation_file(crltns, "CubeT.txt", 0);
+ make_rout_file2(crl_border, index, rout_new, "Routs_Cube.txt");
+ make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_Cube.txt");
std::cout << "Finish Analysis - end\n\n";
}
biod.pnpi.spb.ru / alexxy / gromacs-testing.git / commitdiff