#include <gromacs/utility/smalloc.h>
#include <gromacs/math/do_fit.h>
-#include "fittingn.cpp"
+#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;
}
}
-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+");
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+");
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)
}
}
-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());
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++) {
}
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]);
}
sfree(medx);
sfree(medy);
- std::cout << i << "\n";
+ std::cout << i << " corr done\n";
}
}
#pragma omp barrier
}
}
-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));
}
}
}
+ std::cout << "crl analysed\n";
std::vector< bool > graph_flags;
graph_flags.resize(traj.front().size(), true);
std::vector< int > a;
}
}
-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++) {
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;
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
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++) {
}
// -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++;
}
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;
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";
/*
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++;
}
}
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;
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";
}
biod.pnpi.spb.ru / alexxy / gromacs-testing.git / commitdiff