From 3a57da19729b4c676b4329f9be26b915736fe7ca Mon Sep 17 00:00:00 2001 From: Anatoly Titov Date: Wed, 21 Nov 2018 14:32:19 +0300 Subject: [PATCH] - merged with FitNew - switched full corr output place in code --- src/spacetimecorr.cpp | 371 ++++++++++++++++++++++++++++++++++-------- 1 file changed, 307 insertions(+), 64 deletions(-) diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp index f930fb4..6fa1cbe 100644 --- a/src/spacetimecorr.cpp +++ b/src/spacetimecorr.cpp @@ -49,12 +49,255 @@ #include #include -#include "fittingn.cpp" +#include +#include + +//#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; @@ -68,7 +311,7 @@ void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const ch } } -void make_correlation_file(std::vector< std::vector< std::vector< double > > > correlations, const char* file_name, int start) +void make_correlation_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start) { FILE *file; file = std::fopen(file_name, "w+"); @@ -89,7 +332,7 @@ void make_correlation_file(std::vector< std::vector< std::vector< double > > > c std::fclose(file); } -void make_rout_file2(double crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< int, int > > > rout, const char* file_name) +void make_rout_file2(float crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< int, int > > > rout, const char* file_name) { FILE *file; file = std::fopen(file_name, "w+"); @@ -103,7 +346,7 @@ void make_rout_file2(double crl_border, std::vector< int > indx, std::vector< st std::fclose(file); } -void make_best_corrs_graphics(std::vector< std::vector< std::vector< double > > > correlations, +void make_best_corrs_graphics(std::vector< std::vector< std::vector< float > > > correlations, std::vector< std::vector< std::pair< int, int > > > rout_pairs, std::vector< int > indx, const char* file_name) @@ -143,7 +386,7 @@ bool isitsubset (std::vector< int > a, std::vector< int > b) { } } -void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< double > > > &crl, int tauS, int tauE) { +void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< float > > > &crl, int tauS, int tauE) { crl.resize(tauE + 1); for (int i = 0; i < crl.size(); i++) { crl[i].resize(traj.front().size()); @@ -155,12 +398,16 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame temp_zero[0] = 0; temp_zero[1] = 0; temp_zero[2] = 0; - #pragma omp parallel shared(crl) + + std::vector< float > d; + d.resize(traj.front().size(), 0); + + #pragma omp parallel shared(crl, temp_zero, d) { #pragma omp for schedule(dynamic) for (int i = tauS; i <= tauE; i += 1) { - real tmp2; - rvec *medx, *medy, temp, temp1, temp2; + rvec *medx, *medy, temp1, temp2; + real tmp; snew(medx, traj.front().size()); snew(medy, traj.front().size()); for (int i = 0; i < traj.front().size(); i++) { @@ -169,39 +416,30 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame } for (int j = 0; j < traj.size() - i - 1; j++) { for (int f = 0; f < traj.front().size(); f++) { - rvec_sub(traj[b_frame][f], traj[j][f], temp); - rvec_inc(medx[f], temp); + rvec_sub(traj[b_frame][f], traj[j][f], temp1); + rvec_inc(medx[f], temp1); - rvec_sub(traj[b_frame][f], traj[j + i][f], temp); - rvec_inc(medy[f], temp); + rvec_sub(traj[b_frame][f], traj[j + i][f], temp1); + rvec_inc(medy[f], temp1); } } - - for (int j = 0; j < traj.front().size(); j++) { - tmp2 = 1 / (traj.size() - 1 - i); - - /*temp2 = traj.size() - 1; - temp2 -= i; - temp2 = 1 / temp2;*/ - - copy_rvec(medx[j], temp); - svmul(tmp2, temp, medx[j]); - copy_rvec(medy[j], temp); - svmul(tmp2, temp, medy[j]); - } - std::vector< std::vector< double > > a, b, c; - a.resize(traj.front().size()); - b.resize(traj.front().size()); - c.resize(traj.front().size()); for (int j = 0; j < traj.front().size(); j++) { - a[j].resize(traj.front().size(), 0); - b[j].resize(traj.front().size(), 0); - c[j].resize(traj.front().size(), 0); + tmp = traj.size() - 1; + tmp -= i; + tmp = 1 / tmp; + //tmp = 1 / (traj.size() - 1 - i); + copy_rvec(medx[j], temp1); + svmul(tmp, temp1, medx[j]); + copy_rvec(medy[j], temp1); + svmul(tmp, temp1, medy[j]); } + std::vector< std::vector< float > > a, b, c; + a.resize(traj.front().size(), d); + b.resize(traj.front().size(), d); + c.resize(traj.front().size(), d); for (int j = 0; j < traj.size() - i - 1; j++) { for (int f1 = 0; f1 < traj.front().size(); f1++) { for (int f2 = 0; f2 < traj.front().size(); f2++) { - rvec_sub(traj[b_frame][f1], traj[j][f1], temp1); rvec_dec(temp1, medx[f1]); @@ -221,7 +459,7 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame } sfree(medx); sfree(medy); - std::cout << i << "\n"; + std::cout << i << " corr done\n"; } } #pragma omp barrier @@ -292,9 +530,9 @@ void domain_chopping(SelectionList seList, std::vector< int > indx, std::vector< } } -void graph_calculation(std::vector< std::vector< std::pair< double, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr, +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< double > > > crl, double crl_b, double e_rad, int tauE) { + std::vector< std::vector< std::vector< float > > > crl, float crl_b, float e_rad, int tauE) { graph.resize(traj.front().size()); for (int i = 0; i < traj.front().size(); i++) { graph[i].resize(traj.front().size(), std::make_pair(0, -1)); @@ -316,6 +554,7 @@ void graph_calculation(std::vector< std::vector< std::pair< double, int > > > &g } } } + std::cout << "crl analysed\n"; std::vector< bool > graph_flags; graph_flags.resize(traj.front().size(), true); std::vector< int > a; @@ -361,15 +600,15 @@ void graph_calculation(std::vector< std::vector< std::pair< double, int > > > &g } } -bool myfunction (const std::pair< int, double > i, const std::pair< int, double > j) { +bool myfunction (const std::pair< int, float > i, const std::pair< int, float > j) { return i.second < j.second; } void graph_back_bone_evaluation(std::vector< std::vector < std::pair< int, int > > > &rout_n, int indxSize, - std::vector< std::vector< std::pair< double, int > > > graph, std::vector< std::vector< int > > s_graph, std::vector< std::vector< std::pair< int, int > > > s_graph_rbr) { - std::vector< double > key; + 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< float > key; std::vector< int > path; - std::vector< std::pair< int, double > > que; + std::vector< std::pair< int, float > > que; std::vector< std::pair< int, int > > a; int v; for (int i = 0; i < s_graph.size(); i++) { @@ -473,9 +712,8 @@ class Domains : public TrajectoryAnalysisModule int frames = 0; int basic_frame = 0; int tau = 0; - int graph_tau = 0; - double crl_border = 0; - double 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; @@ -504,20 +742,16 @@ Domains::initOptions(IOptionsContainer *options, options->addOption(FileNameOption("on").filetype(eftIndex).outputFile() .store(&fnNdx_).defaultBasename("domains") .description("Index file from the domains")); - // Add option for graph_tau constant - options->addOption(gmx::IntegerOption("Gtau") - .store(&graph_tau) - .description("number of frames for graph to see for time travel")); // Add option for tau constant options->addOption(gmx::IntegerOption("tau") .store(&tau) .description("number of frames for time travel")); // Add option for crl_border constant - options->addOption(DoubleOption("crl") + options->addOption(FloatOption("crl") .store(&crl_border) .description("make graph based on corrs > constant")); // Add option for eff_rad constant - options->addOption(DoubleOption("ef_rad") + options->addOption(FloatOption("ef_rad") .store(&eff_rad) .description("effective radius for atoms to evaluate corrs")); // Add option for selection list @@ -558,6 +792,8 @@ Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, 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++) { @@ -574,27 +810,30 @@ Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, } // -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 void Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata) { - trajectory.back().resize(0); - 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()); - for (int i = 0; i < domains_local.size(); i++) { - std::vector< RVec > basic, traj; - basic = trajectory[basic_frame]; - traj = trajectory.back(); - MyFitNew(basic, 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 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++; } @@ -602,7 +841,7 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, void Domains::finishAnalysis(int nframes) { - std::vector< std::vector< std::vector< double > > > crltns; + std::vector< std::vector< std::vector< float > > > crltns; int k = 1000, m = 0; if (tau > 0) { k = tau; @@ -616,6 +855,7 @@ Domains::finishAnalysis(int nframes) 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"; /* @@ -630,7 +870,7 @@ Domains::finishAnalysis(int nframes) 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]) >= (double)i1 / 100 && i3 != i4) { + if (std::abs(crltns[i2][i3][i4]) >= (float)i1 / 100 && i3 != i4) { i5++; } } @@ -644,12 +884,12 @@ Domains::finishAnalysis(int nframes) std::cout << "graph evaluation: start\n"; - std::vector< std::vector< std::pair< double, int > > > graph; + 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; - graph_calculation(graph, sub_graph, sub_graph_rbr, frankenstein_trajectory, basic_frame, crltns, crl_border, eff_rad, graph_tau/*k*/); + 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; @@ -661,8 +901,11 @@ Domains::finishAnalysis(int nframes) std::cout << "routs evaluation: end\n"; - make_rout_file2(crl_border, index, rout_new, "Routs_DomainsNFit_5000_0.70_1_t.txt"); - make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_DomainsNFit_5000_0.70_1_t.txt"); + + // 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"); + std::cout << "Finish Analysis - end\n\n"; } -- 2.22.0