#include <iostream>
#include <chrono>
+#include <string.h>
+#include <algorithm>
+
#include <omp.h>
#include <gromacs/trajectoryanalysis.h>
#include <gromacs/utility/smalloc.h>
#include <gromacs/math/do_fit.h>
-#include "fittingn.h"
-
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 FtCnst) {
+ 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 (FtCnst == 0) {
+ FtCnst = 0.00001;
+ }
+ if (f1 == 0) {
+ return;
+ } else {
+ int count = 0;
+ while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
+ 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);
+ }
+ 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++;
+ }
+ ApplyFit(b, x, y, z, A, B, C);
+ }
+}
+
+void printPDBtraj(const char* Fname, std::vector< std::vector < RVec > > trj, std::vector< std::vector< std::pair< int, int > > > prs) {
+ FILE* a = freopen (Fname, "w", stdout);
+ for (int i = 1; i < trj.size(); i++) {
+ for (int j = 0; j < prs.size(); j++) {
+ for (int k = 0; k < prs[j].size(); k++) {
+ printf("ATOM %5d CA CA %5d %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s \n", prs[j][k].first, prs[j][k].first, (trj[i][prs[j][k].first][0] * 10), (trj[i][prs[j][k].first][1] * 10), (trj[i][prs[j][k].first][2] * 10), 1.0, 20.0, " ", " ");
+ }
+ }
+ printf("ENDMDL\n");
+ }
+ printf("END\n");
+ fclose(a);
+}
+
+bool myfunction (const std::vector< std::pair< int, int > > i, const std::vector< std::pair< int, int > > j) {
+ return i.size() < j.size();
+}
+
+void domain_reading(SelectionList sList, std::vector< int > ind, std::vector< std::vector < std::pair< int, int > > > &d_ind_num) {
+ d_ind_num.resize(sList.size() - 1);
+ for (int i = 1; i < sList.size(); i++) {
+ d_ind_num[i - 1].resize(0);
+ for (ArrayRef< const int >::iterator ai = sList[i].atomIndices().begin(); (ai < sList[i].atomIndices().end()); ai++) {
+ d_ind_num[i - 1].push_back(std::make_pair(*ai, 0));
+ }
+ }
+
+ for (int i = 0; i < d_ind_num.size(); i++) {
+ int k;
+ for (int j = 0; j < d_ind_num[i].size(); j++) {
+ k = 0;
+ while (ind[k] != d_ind_num[i][j].first) {
+ k++;
+ }
+ d_ind_num[i][j].second = k;
+ }
+ }
+
+ // its a point of interest as we can either take biggest domains of the given set and work from those or we take them in given order
+ //std::sort(d_ind_num.begin(), d_ind_num.end(), myfunction);
+
+ for (int i = 0; i < d_ind_num.size(); i++) {
+ if (d_ind_num[i].size() >= 4) {
+ for (int k = 0; k < d_ind_num[i].size(); k++) {
+ for (int j = i + 1; j < d_ind_num.size(); j++) {
+ for (int x = 0; x < d_ind_num[j].size(); x++) {
+ if (d_ind_num[j][x].first == d_ind_num[i][k].first) {
+ d_ind_num[j].erase(d_ind_num[j].begin() + x);
+ x--;
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+void domain_expansion(std::vector< int > indx, std::vector< std::vector < std::pair< int, int > > > &d_ind_num, std::vector< RVec > frame) {
+ std::vector< bool > flag;
+ flag.resize(indx.size(), 1);
+ for (int i = 0; i < d_ind_num.size(); i++) {
+ for (int j = 0; j < d_ind_num[i].size(); j++) {
+ flag[d_ind_num[i][j].second] = 0;
+ }
+ }
+ float dist;
+ int x3 = -1;
+ for (int i = 0; i < indx.size(); i++) {
+ if (flag[i]) {
+ dist = 0;
+ for (int x1 = 0; x1 < d_ind_num.size(); x1++) {
+ for (int x2 = 0; x2 < d_ind_num[x1].size(); x2++) {
+ if (dist == 0) {
+ dist = (frame[i] - frame[d_ind_num[0][0].first]).norm();
+ x3 = 0;
+ }
+ if ((frame[i] - frame[d_ind_num[x1][x2].first]).norm() < dist) {
+ x3 = x1;
+ dist = (frame[i] - frame[d_ind_num[x1][x2].first]).norm();
+ }
+ }
+ }
+ d_ind_num[x3].push_back(std::make_pair(indx[i], i));
+ }
+ }
+}
+
+void MakeDomainFitPairs(std::vector< std::vector< std::pair< int, int > > > &p, std::vector< std::vector< std::pair< int, int > > > d) {
+ p.resize(d.size());
+ for (int i = 0; i < d.size(); i++) {
+ p[i].resize(0);
+ for (int j = 0; j < d[i].size(); j++) {
+ p[i].push_back(std::make_pair(d[i][j].second, d[i][j].second));
+ }
+ }
+}
+
+void TrajFitting(std::vector< std::vector < RVec > > &trj, double FC, std::vector< std::vector< std::pair< int, int > > > prs) {
+ std::vector< std::vector < RVec > > clone;
+ clone.resize(0);
+ // its needed to be thought through on whats shared and whats private
+ //#pragma omp parallel for shared(trajectory, pairs, FitConst) private(clone)
+ for (int i = 1; i < trj.size(); i++) {
+ clone.resize(prs.size(), trj[i]);
+ for (int j = 0; j < prs.size(); j++) {
+ MyFitNew(trj[0], clone[j], prs[j], FC);
+ for (int k = 0; k < prs[j].size(); k++) {
+ trj[i][prs[j][k].first] = clone[j][prs[j][k].first];
+ }
+ }
+ clone.resize(0);
+ std::cout << "frame " << i << " fitted\n";
+ }
+}
+
class Fitng : public TrajectoryAnalysisModule
{
public:
virtual void writeOutput();
private:
-
- std::string fnNdx_;
+ SelectionList sel_;
std::vector< int > index;
+ std::vector< std::vector< std::pair< int, int > > > domains_index_number;
std::vector< std::vector < RVec > > trajectory;
- Selection selec;
- int frames = 0;
-
- real **w_rls;
-
- int Fitng_ePBC;
- // Copy and assign disallowed by base.
+ //Selection selec;
+ int frames = 0;
+ int method = 1;
+ double FitConst = 0;
+ std::string OutPutTrjName;
};
Fitng::Fitng(): TrajectoryAnalysisModule()
// Add the descriptive text (program help text) to the options
settings->setHelpText(desc);
// Add option for selecting a subset of atoms
- options->addOption(SelectionOption("select")
- .store(&selec).required()
- .description("Atoms that are considered as part of the excluded volume"));
- // Add option for output file name
- options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
- .store(&fnNdx_).defaultBasename("Fitng")
- .description("Index file from the Fitng"));
- // Add option for etalon_frame constant
- //options->addOption(gmx::IntegerOption("dms")
- // .store(&domain_min_size)
- // .description("minimum domain size"));
- // Add option for epsi constant
- //options->addOption(DoubleOption("epsilon")
- // .store(&epsi)
- // .description("thermal vibrations' constant"));
- // Add option for delta constant
- //options->addOption(DoubleOption("delta")
- // .store(&delta)
- // .description("domain membership probability"));
+ //options->addOption(SelectionOption("select")
+ // .store(&selec).required()
+ // .description("Atoms that are considered as part of the excluded volume"));
+ // Add option for selection list
+ options->addOption(SelectionOption("select_residue_and_domains").storeVector(&sel_)
+ .required().dynamicMask().multiValue()
+ .description("Selected group and domains based on it"));
+ // Add option for Fit constant
+ options->addOption(DoubleOption("FitConst")
+ .store(&FitConst)
+ .description("Fitting untill diff <= FitConst, by default == 0.00001"));
+ // Add option for output pdb traj (meh edition)
+ options->addOption(StringOption("pdb_out")
+ .store(&OutPutTrjName)
+ .description("transformed trajectory"));
+ // Add option for trajectory Fit transformation
+ options->addOption(IntegerOption("method")
+ .store(&method)
+ .description("1 - fitting all -> all | 2 - keep only domains, all -> all | 3 - keep only domains, domain -> domain | 4 - extend domains to full coverage, domain -> domain"));
// Control input settings
settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
settings->setPBC(true);
}
-void
-Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings,
- const TopologyInformation &top)
-{
- Fitng_ePBC = top.ePBC();
-}
+// Fitng -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA'
+// Fitng -s '/home/toluk/Develop/FitSamples/cube_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/cube.pdb' -select 'name NA'
+// Fitng -s '/home/toluk/Develop/FitSamples/cube_with_termal_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_termal.pdb' -select 'name NA'
+// Fitng -s '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_1st_fr_D100.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_D100.pdb' -select 'name NA'
+// Fitng -s '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_1st_fr_D10.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_D10.pdb' -select 'name NA'
+// Fitng -s '/home/toluk/Develop/FitSamples/garmoshka_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/garmoshka.pdb' -select 'name NA'
+// Fitng -s '/home/toluk/Develop/FitSamples/garmoshka_2nd_fr.pdb' -f '/home/toluk/Develop/FitSamples/garmoshka.pdb' -select 'name NA'
+// Fitng -s '/home/toluk/Develop/FitSamples/rubber_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/rubber.pdb' -select 'name NA'
+// rubber_1st_fr
+
+// /home/toluk/Develop/samples/reca_test_ground/
+// -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL0' -n '/home/toluk/Develop/samples/reca_test_ground/t0.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/FullCAFit.pdb' -e 1000
+// -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsFit.pdb' -e 1000 -method 2
+// -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsOnlyFit.pdb' -e 1000 -method 3
+// -s '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_test_ground/reca_rd.mono.xtc' -sf '/home/toluk/Develop/samples/reca_test_ground/SL2' -n '/home/toluk/Develop/samples/reca_test_ground/t2.ndx' -pdb_out '/home/toluk/Develop/samples/reca_test_ground/CAdomainsWExtraOnlyFit.pdb' -e 1000 -method 4
void
-Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
- const t_trxframe &fr)
+Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
{
- t_pbc pbc;
- t_pbc *ppbc = settings.hasPBC() ? &pbc : NULL;
- matrix boxx;
- copy_mat(fr.box, boxx);
- if (ppbc != NULL) {
- set_pbc(ppbc, Fitng_ePBC, boxx);
- }
- ConstArrayRef< int > atomind = selec.atomIndices();
index.resize(0);
- for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
+ for (ArrayRef< const int >::iterator ai = sel_[0].atomIndices().begin(); (ai < sel_[0].atomIndices().end()); ai++) {
index.push_back(*ai);
}
- trajectory.resize(1);
- trajectory.back().resize(selec.atomCount());
+ trajectory.resize(0);
- for (int i = 0; i < selec.atomCount(); i++) {
- trajectory.back()[i] = fr.x[index[i]];
+ if (sel_.size() > 1 && method >= 2) {
+ domain_reading(sel_, index, domains_index_number);
}
}
+void
+Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr)
+{
+
+}
+
void
Fitng::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
TrajectoryAnalysisModuleData *pdata)
}
}
-// Fitng -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA'
-// Fitng -s '/home/toluk/Develop/FitSamples/cube_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/cube.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/cube_with_termal_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_termal.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_1st_fr_D100.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_D100.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_1st_fr_D10.pdb' -f '/home/toluk/Develop/FitSamples/cube_with_small_cube_floating_D10.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/garmoshka_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/garmoshka.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/garmoshka_2nd_fr.pdb' -f '/home/toluk/Develop/FitSamples/garmoshka.pdb' -select 'name NA'
-// Fitng -s '/home/toluk/Develop/FitSamples/rubber_1st_fr.pdb' -f '/home/toluk/Develop/FitSamples/rubber.pdb' -select 'name NA'
-// rubber_1st_fr
void
Fitng::finishAnalysis(int nframes)
{
- std::vector< std::pair< int, int > > pairs;
+ std::vector< std::vector< std::pair< int, int > > > pairs;
+
+ switch (method) {
+ case 1:
+ pairs.resize(1);
+ for (int i = 0; i < index.size(); i++) {
+ pairs[0].push_back(std::make_pair(i, i));
+ }
+ TrajFitting(trajectory, FitConst, pairs);
+ break;
+ case 2:
+ pairs.resize(1);
+ for (int i = 0; i < domains_index_number.size(); i++) {
+ for (int j = 0; j < domains_index_number[i].size(); j++) {
+ pairs[0].push_back(std::make_pair(domains_index_number[i][j].second, domains_index_number[i][j].second));
+ }
+ }
+ TrajFitting(trajectory, FitConst, pairs);
+ break;
+ case 3:
+ MakeDomainFitPairs(pairs, domains_index_number);
+ TrajFitting(trajectory, FitConst, pairs);
+ break;
+ case 4:
+ domain_expansion(index, domains_index_number, trajectory[0]);
+ MakeDomainFitPairs(pairs, domains_index_number);
+ TrajFitting(trajectory, FitConst, pairs);
+ break;
+ }
+
+ printPDBtraj(OutPutTrjName.c_str(), trajectory, pairs);
+}
+
+void
+Fitng::writeOutput()
+{
+
+}
+
+
+/*! \brief
+ * The main function for the analysis template.
+ */
+int
+main(int argc, char *argv[])
+{
+ return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Fitng>(argc, argv);
+}
+
+
+
+/* std::vector< std::pair< int, int > > pairs;
for (int i = 0; i < index.size(); i++) {
pairs.push_back(std::make_pair(i, i));
}
for (int j = 1; j < trajectory.size(); j++) {
for (int i = 0; i < index.size(); i++) {
- dist1 += ftn::F( static_cast< double >(trajectory[0][i][0]), static_cast< double >(trajectory[0][i][1]), static_cast< double >(trajectory[0][i][2]),
+ dist1 += ( static_cast< double >(trajectory[0][i][0]), static_cast< double >(trajectory[0][i][1]), static_cast< double >(trajectory[0][i][2]),
static_cast< double >(trajectory[j][i][0] + 0), static_cast< double >(trajectory[j][i][1] + 0), static_cast< double >(trajectory[j][i][2] + 0),
static_cast< double >(0), static_cast< double >(0), static_cast< double >(0) );
}
- std::cout << "\nbasic dist = " << dist1 /*/ index.size() */<< "\n";
+ std::cout << "\nbasic dist = " << dist1 /* index.size() << "\n";
//
// My Fit
//
start = std::chrono::system_clock::now();
- ftn::MyFitNew(trajectory[0], trajectory[j], pairs);
+ /*ftn:: MyFitNew(trajectory[0], trajectory[j], pairs, FitConst);
end = std::chrono::system_clock::now();
std::cout << std::chrono::duration_cast<std::chrono::microseconds> (end-start).count() << " ";
for (int i = 0; i < index.size(); i++) {
- dist2 += ftn::F( static_cast< double >(trajectory[0][i][0]), static_cast< double >(trajectory[0][i][1]), static_cast< double >(trajectory[0][i][2]),
+ dist2 += /*ftn::F ( static_cast< double >(trajectory[0][i][0]), static_cast< double >(trajectory[0][i][1]), static_cast< double >(trajectory[0][i][2]),
static_cast< double >(trajectory[j][i][0] + 0), static_cast< double >(trajectory[j][i][1] + 0), static_cast< double >(trajectory[j][i][2] + 0),
static_cast< double >(0), static_cast< double >(0), static_cast< double >(0));
}
- std::cout << "my fit dist = " << dist2 /*/ index.size()*/<< "\n";
+ //std::cout << "my fit dist = " << dist2 /* index.size()<< "\n";
//
// Old Fit
//end = std::chrono::system_clock::now();
//std::cout << std::chrono::duration_cast<std::chrono::microseconds> (end-start).count() << "\n";
for (int i = 0; i < index.size(); i++) {
- dist2 += ftn::F( static_cast< double >(traj1[i][0]), static_cast< double >(traj1[i][1]), static_cast< double >(traj1[i][2]),
+ dist2 += /*ftn:: F( static_cast< double >(traj1[i][0]), static_cast< double >(traj1[i][1]), static_cast< double >(traj1[i][2]),
static_cast< double >(traj2[j][i][0] + 0), static_cast< double >(traj2[j][i][1] + 0), static_cast< double >(traj2[j][i][2] + 0),
static_cast< double >(0), static_cast< double >(0), static_cast< double >(0));
}
- std::cout << "old fit dist = " << dist2 /*/ index.size()*/<< "\n";
+ std::cout << "old fit dist = " << dist2 / index.size()<< "\n";
dist1 = 0;
dist2 = 0;
}
- /*int co = 0;
- freopen ("/home/toluk/Develop/FitSamples/old_fit_result.pdb","w",stdout);
+ printPDBtraj(OutPutTrjName.c_str(), trajectory);
+
+ int co = 0;
+ freopen ("/home/toluk/Develop/FitSamples/old_fit_result.pdb", "w", stdout);
for (int l = 1; l < trajectory.size(); l++) {
for (int i = 0; i < index.size(); i++) {
//cout << D*i + (float)l * 6;
printf("ENDMDL\n");
}
printf("END\n");
- freopen ("/home/toluk/Develop/FitSamples/new_fit_result.pdb","w",stdout);
+ freopen ("/home/toluk/Develop/FitSamples/new_fit_result.pdb", "w", stdout);
for (int l = 1; l < trajectory.size(); l++) {
for (int i = 0; i < index.size(); i++) {
//cout << D*i + (float)l * 6;
printf("ENDMDL\n");
}
printf("END\n");
- fclose(stdout);*/
+ fclose(stdout);
sfree(traj1);
for (int i = 0; i < trajectory.size(); i++) {
sfree(traj2[i]);
}
sfree(traj2);
-}
-
-void
-Fitng::writeOutput()
-{
-
-}
-
-
-/*! \brief
- * The main function for the analysis template.
- */
-int
-main(int argc, char *argv[])
-{
- return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Fitng>(argc, argv);
-}
+*/
biod.pnpi.spb.ru / alexxy / gromacs-fitng.git / commitdiff