#include <omp.h>
#include <thread>
#include <string>
+#include <math.h>
+#include <algorithm>
#include <gromacs/trajectoryanalysis.h>
#include <gromacs/pbcutil/pbc.h>
#include <gromacs/math/vec.h>
#include <gromacs/math/do_fit.h>
-#define MAX_NTRICVEC 12
-
using namespace gmx;
using gmx::RVec;
-struct node {
- short int n;
- RVec r;
- bool yep;
-};
+double F (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
+ return sqrt( pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
+ pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2) +
+ pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) );
+}
-void make_graph(int mgwi_natoms, rvec *mgwi_x, std::vector< std::vector< node > > &mgwi_graph)
-{
- mgwi_graph.resize(mgwi_natoms);
- for (int i = 0; i < mgwi_natoms; i++) {
- mgwi_graph[i].resize(mgwi_natoms);
- }
- for (int i = 0; i < mgwi_natoms; i++) {
- for (int j = 0; j < mgwi_natoms; j++) {
- rvec_sub(mgwi_x[i], mgwi_x[j], mgwi_graph[i][j].r);
- mgwi_graph[i][j].n = 0;
- }
- }
+double Fx (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
+ return -(2 * cos(B) * cos(C) * (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) -
+ 2 * sin(B) * (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y)) +
+ 2 * cos(B) * sin(C) * (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x))) /
+ (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
+ pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
+ pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
}
-void update_graph(std::vector< std::vector< node > > &ugwi_graph, rvec *ugwi_x, long double ugwi_epsi) {
- rvec ugwi_temp;
- int ugwi_for = ugwi_graph.size();
- for (int i = 0; i < ugwi_for; i++) {
- for (int j = i; j < ugwi_for; j++) {
- rvec_sub(ugwi_x[i], ugwi_x[j], ugwi_temp);
- rvec_dec(ugwi_temp, ugwi_graph[i][j].r.as_vec());
- if (norm(ugwi_temp) <= ugwi_epsi) {
- if (i == j) {
- ugwi_graph[i][j].n++;
- }
- else {
- ugwi_graph[i][j].n++;
- ugwi_graph[j][i].n++;
- }
- }
- }
- }
+double Fy (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
+ return -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) *
+ (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x)) -
+ 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) *
+ (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) +
+ 2 * cos(B) * sin(A) * (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y))) /
+ (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
+ pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
+ pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
}
-void check_domains(long double cd_delta, int cd_frames, std::vector< std::vector< std::vector< node > > > &cd_graph) {
- int cd_for1 = cd_graph.size(), cd_for2 = cd_graph[1].size();
- for (int k = 0; k < cd_for1; k++) {
- for (int i = 0; i < cd_for2; i++) {
- for (int j = 0; j < cd_for2; j++) {
- if (cd_graph[k][i][j].n >= cd_frames * cd_delta) {
- cd_graph[k][i][j].yep = true;
- }
- else {
- cd_graph[k][i][j].yep = false;
- }
- }
- }
- }
+double Fz (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
+ return -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) *
+ (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) -
+ 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) *
+ (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x)) +
+ 2 * cos(A) * cos(B) * (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y))) /
+ (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
+ pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
+ pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
}
-void find_domain_sizes(std::vector< std::vector< std::vector< node > > > fds_graph, std::vector< std::vector< int > > &fds_domsizes) {
- fds_domsizes.resize(fds_graph.size());
- int fds_for1 = fds_graph.size(), fds_for2 = fds_graph[1].size();
- for (int i = 0; i < fds_for1; i++) {
- fds_domsizes[i].resize(fds_for2, 0);
- for (int j = 0; j < fds_for2; j++) {
- for (int k = 0; k < fds_for2; k++) {
- if (fds_graph[i][j][k].yep) {
- fds_domsizes[i][j]++;
- }
- }
- }
- }
+double FA (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
+ return -(2 * (cos(A) * cos(B) * (biy + y) - cos(B) * sin(A) * (biz + z)) * (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y)) -
+ 2 * ((biy + y) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + (biz + z) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) *
+ (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x)) +
+ 2 * ((biy + y) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + (biz + z) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) *
+ (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x))) /
+ (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
+ pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
+ pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
}
-void get_maxsized_domain(std::vector< int > &gmd_max_d, std::vector< std::vector< std::vector< node > > > gmd_graph, std::vector< std::vector< int > > gmd_domsizes) {
- int gmd_number1 = 0, gmd_number2 = 0;
- int gmd_for1 = gmd_domsizes.size(), gmd_for2 = gmd_domsizes[0].size(), gmd_for3 = gmd_graph[1][1].size();
- for (int i = 0; i < gmd_for1; i++) {
- for (int j = 0; j < gmd_for2; j++) {
- if (gmd_domsizes[i][j] >= gmd_domsizes[gmd_number1][gmd_number2]) {
- gmd_number1 = i;
- gmd_number2 = j;
- }
- }
- }
- gmd_max_d.resize(0);
- for (int i = 0; i < gmd_for3; i++) {
- if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
- gmd_max_d.push_back(i);
- }
- }
+double FB (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
+ return -(2 * (cos(A) * cos(B) * sin(C) * (biz + z) - sin(B) * sin(C) * (bix + x) + cos(B) * sin(A) * sin(C) * (biy + y)) *
+ (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x)) +
+ 2 * (cos(A) * cos(B) * cos(C) * (biz + z) - cos(C) * sin(B) * (bix + x) + cos(B) * cos(C) * sin(A) * (biy + y)) *
+ (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) -
+ 2 * (cos(B) * (bix + x) + sin(A) * sin(B) * (biy + y) + cos(A) * sin(B) * (biz + z)) *
+ (aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y))) /
+ (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
+ pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
+ pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
}
-void delete_domain_from_graph(std::vector< std::vector< std::vector< node > > > &ddf_graph, std::vector< int > ddf_domain) {
- int ddfg_for1 = ddf_domain.size(), ddfg_for2 = ddf_graph.size(), ddfg_for3 = ddf_graph[1].size();
- for (int i = 0; i < ddfg_for1; i++) {
- for (int k = 0; k < ddfg_for2; k++) {
- for (int j = 0; j < ddfg_for3; j++) {
- if (ddf_graph[k][ddf_domain[i]][j].yep) {
- ddf_graph[k][ddf_domain[i]][j].yep = false;
- }
- if (ddf_graph[k][j][ddf_domain[i]].yep) {
- ddf_graph[k][j][ddf_domain[i]].yep = false;
- }
- }
- }
+double FC (double aix, double aiy, double aiz, double bix, double biy, double biz, double x, double y, double z, double A, double B, double C) {
+ return (2 * ((biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (bix + x)) *
+ (aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x)) -
+ 2 * ((biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (bix + x)) *
+ (aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x))) /
+ (2 * sqrt(pow(aiz + sin(B) * (bix + x) - cos(A) * cos(B) * (biz + z) - cos(B) * sin(A) * (biy + y), 2) +
+ pow(aix + (biy + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - (biz + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * (bix + x), 2) +
+ pow(aiy - (biy + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + (biz + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * (bix + x), 2)));
+}
+
+void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
+ RVec temp;
+ 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);
}
}
-bool check_domsizes(std::vector< std::vector< int > > cd_domsizes, int cd_domain_min_size) {
- int cd_for1 = cd_domsizes.size(), cd_for2 = cd_domsizes[0].size();
- for (int i = 0; i < cd_for1; i++) {
- for (int j = 0; j < cd_for2; j++) {
- if (cd_domsizes[i][j] >= cd_domain_min_size) {
- return true;
- }
- }
+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]);
}
- return false;
+ 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 print_domains(std::vector< std::vector< int > > pd_domains, std::vector< int > index, std::string fnNdx_) {
- FILE *fpNdx_;
- fpNdx_ = std::fopen(fnNdx_.c_str(), "w+");
- int write_count;
- for (int i = 0; i < pd_domains.size(); i++) {
- std::fprintf(fpNdx_, "[domain_%d]\n", i + 1);
- write_count = 0;
- for (int j = 0; j < pd_domains[i].size(); j++) {
- write_count++;
- if (write_count > 20) {
- write_count -= 20;
- std::fprintf(fpNdx_, "\n");
+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;
+ double f21 = 0, f22 = 0, f23 = 0, f24 = 0, f25 = 0, f26 = 0;
+ double l1 = 1, l2 = 1, l3 = 1, l4 = 1, l5 = 1, l6 = 1, 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;
+
+ int count;
+ while (true) {
+ 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], b[pairs[i].second][1], b[pairs[i].second][2], x, y, z, A, B, C);
+ fx += Fx(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);
+ fy += Fy(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);
+ fz += Fz(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);
+ fa += FA(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);
+ fb += FB(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);
+ fc += FC(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);
+ }
+ 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], b[pairs[i].second][1], b[pairs[i].second][2], x - l * fx, y - l * fy, z - l * fz, A - l * fa, B - l * fb, C - l * fc);
+ /*
+ f2 += F(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 - l1 * fx, y - l2 * fy, z - l3 * fz, A - l4 * fa, B - l5 * fb, C - l6 * fc);
+ f21 += F(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 - l1 * fx, y, z, A, B, C);
+ f22 += F(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 - l2 * fy, z, A, B, C);
+ f23 += F(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 - l3 * fz, A, B, C);
+ f24 += F(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 - l4 * fa, B, C);
+ f25 += F(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 - l5 * fb, C);
+ f26 += F(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 - l6 * fc);
+ */
}
- std::fprintf(fpNdx_, "%5d ", index[pd_domains[i][j]] + 1);
+ std::cout << f1 << " " << f2 << "\n";
+ if (f2 < f1) {
+ x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
+ //x -= l1 * fx; y -= l2 * fy; z -= l3 * fz; A -= l4 * fa; B -= l5 * fb; C -= l6 * fc;
+ fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
+ break;
+ } else {
+ l *= 0.85;
+ /*
+ if (f21 > f1) {
+ l1 /= 2;
+ }
+ if (f22 > f1) {
+ l2 /= 2;
+ }
+ if (f23 > f1) {
+ l3 /= 2;
+ }
+ if (f24 > f1) {
+ l4 /= 2;
+ }
+ if (f25 > f1) {
+ l5 /= 2;
+ }
+ if (f26 > f1) {
+ l6 /= 2;
+ }
+ */
+ }
+ }
+ if (f1 - f2 > 0 && f1 - f2 < 0.00001) {
+ break;
}
- std::fprintf(fpNdx_,"\n\n");
+ f1 = 0; f2 = 0;
}
- std::fprintf(fpNdx_,"\n");
- std::fclose(fpNdx_);
+ ApplyFit(b, x, y, z, A, B, C);
}
-
-/*! \brief
- * Class used to compute free volume in a simulations box.
- *
- * Inherits TrajectoryAnalysisModule and all functions from there.
- * Does not implement any new functionality.
- *
- * \ingroup module_trajectoryanalysis
- */
-class Domains : public TrajectoryAnalysisModule
+class Fitng : public TrajectoryAnalysisModule
{
public:
- Domains();
- virtual ~Domains();
+ Fitng();
+ virtual ~Fitng();
//! Set the options and setting
virtual void initOptions(IOptionsContainer *options,
std::string fnNdx_;
- std::vector< std::vector< std::vector< node > > > graph;
-
- std::vector< std::vector< int > > domains;
- std::vector< std::vector< int > > domsizes;
-
std::vector< int > index;
- std::vector< int > numbers;
std::vector< std::vector < RVec > > trajectory;
Selection selec;
int frames = 0;
- int domain_min_size = 5; // should be selectable
real **w_rls;
- int bone;
- double delta = 0.90; //0.95 // should be selectable
- double epsi = 0.15; //0.3 колебания внутри домена // should be selectable
- int domains_ePBC;
+ int Fitng_ePBC;
// Copy and assign disallowed by base.
};
-Domains::Domains(): TrajectoryAnalysisModule()
+Fitng::Fitng(): TrajectoryAnalysisModule()
{
}
-Domains::~Domains()
+Fitng::~Fitng()
{
}
void
-Domains::initOptions(IOptionsContainer *options,
+Fitng::initOptions(IOptionsContainer *options,
TrajectoryAnalysisSettings *settings)
{
static const char *const desc[] = {
.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("domains")
- .description("Index file from the domains"));
+ .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"));
+ //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"));
+ //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(DoubleOption("delta")
+ // .store(&delta)
+ // .description("domain membership probability"));
// Control input settings
settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
settings->setPBC(true);
}
void
-Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
+Fitng::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top)
{
- domains_ePBC = top.ePBC();
+ Fitng_ePBC = top.ePBC();
}
void
-Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
+Fitng::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
const t_trxframe &fr)
{
t_pbc pbc;
matrix boxx;
copy_mat(fr.box, boxx);
if (ppbc != NULL) {
- set_pbc(ppbc, domains_ePBC, boxx);
+ 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++) {
index.push_back(*ai);
}
- trajectory.resize(2);
- trajectory[0].resize(selec.atomCount());
+ trajectory.resize(1);
+ trajectory.back().resize(selec.atomCount());
for (int i = 0; i < selec.atomCount(); i++) {
- trajectory[0][i] = fr.x[index[i]];
- }
-
- bone = index.size() - domain_min_size + 1;
- graph.resize(bone);
- snew(w_rls, bone);
- for (int i = 0; i < bone; i++) {
- snew(w_rls[i], index.size());
- for (int j = 0; j < index.size(); j++) {
- if (j >= i && j <= i + domain_min_size - 1) {
- w_rls[i][j] = 1;
- } else {
- w_rls[i][j] = 0;
- }
- }
- rvec *etalon;
- snew(etalon, index.size());
- for (int j = 0; j < index.size(); j++) {
- copy_rvec(trajectory[0][j].as_vec(), etalon[j]);
- }
- reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[i]);
- make_graph(index.size(), etalon, graph[i]);
- sfree(etalon);
+ trajectory.back()[i] = fr.x[index[i]];
}
- trajectory[1].resize(index.size());
}
void
-Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+Fitng::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
TrajectoryAnalysisModuleData *pdata)
{
+ trajectory.resize(trajectory.size() + 1);
+ trajectory.back().resize(index.size());
for (int i = 0; i < index.size(); i++) {
- trajectory[1][i] = fr.x[index[i]];
- }
- frames++;
-
- #pragma omp parallel
- {
- #pragma omp for schedule(dynamic)
- for (int j = 0; j < bone; j++) {
- rvec *etalon, *traj;
- snew(etalon, index.size());
- for (int k = 0; k < index.size(); k++) {
- copy_rvec(trajectory[0][k].as_vec(), etalon[k]);
- }
- snew(traj, index.size());
- for (int k = 0; k < index.size(); k++) {
- copy_rvec(trajectory[1][k].as_vec(), traj[k]);
- }
- reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[j]);
- reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[j]);
- do_fit(index.size(), w_rls[j], etalon, traj);
- update_graph(graph[j], traj, epsi);
- sfree(etalon);
- sfree(traj);
- }
+ trajectory.back()[i] = fr.x[index[i]];
}
- std::cout << "frame: " << frames << " analyzed\n";
}
-//domains -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'
-//domains -s '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.tpr' -f '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.10th.xtc' -select 'name CA'
-//domains -s '/home/toluk/Data/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Data/reca_rd/reca_rd.mono.xtc' - select 'name CA'
+//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'
void
-Domains::finishAnalysis(int nframes)
+Fitng::finishAnalysis(int nframes)
{
- frames -= 1;
-
- std::cout << "final cheking\n";
- check_domains(delta, frames, graph);
-
- std::cout << "finding domains' sizes\n";
- find_domain_sizes(graph, domsizes);
-
- std::cout << "finding domains\n";
- std::vector< int > a;
- a.resize(0);
- while (check_domsizes(domsizes, domain_min_size)) {
- domains.push_back(a);
- get_maxsized_domain(domains.back(), graph, domsizes);
- delete_domain_from_graph(graph, domains.back());
- domsizes.resize(0);
- find_domain_sizes(graph, domsizes);
+ 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 i = 0; i < bone; i++) {
- sfree(w_rls[i]);
+
+ double dist1 = 0, dist2 = 0;
+
+ for (int i = 0; i < index.size(); i++) {
+ dist1 += F(trajectory[0][i][0], trajectory[0][i][1], trajectory[0][i][2], trajectory.back()[i][0], trajectory.back()[i][1], trajectory.back()[i][2], 0, 0, 0, 0, 0, 0);
}
- sfree(w_rls);
+
+ std::cout << "\n\n\n";
+
+ MyFitNew(trajectory[0], trajectory.back(), pairs);
+
+ for (int i = 0; i < index.size(); i++) {
+ dist2 += F(trajectory[0][i][0], trajectory[0][i][1], trajectory[0][i][2], trajectory.back()[i][0], trajectory.back()[i][1], trajectory.back()[i][2], 0, 0, 0, 0, 0, 0);
+ }
+
+ std::cout << "\n\n\n" << "old dist = " << dist1 << "\nnew dist = " << dist2 << "\n\n";
+
}
void
-Domains::writeOutput()
+Fitng::writeOutput()
{
- std::cout << "making output file\n";
- print_domains(domains, index, fnNdx_); // see function for details | numbers from index
- std::cout << "\n END \n";
+
}
int
main(int argc, char *argv[])
{
- return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
+ return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Fitng>(argc, argv);
}
biod.pnpi.spb.ru / alexxy / gromacs-fitng.git / commitdiff