#include <iostream>
#include <chrono>
+#include <cstdint>
+#include <cfloat>
#include <omp.h>
#include <gromacs/trajectoryanalysis.h>
#include <gromacs/utility/smalloc.h>
#include <gromacs/math/do_fit.h>
+#include <gromacs/trajectoryanalysis/topologyinformation.h>
+#include "gromacs/topology/topology.h"
+#include "gromacs/topology/index.h"
+
+
#define MAX_NTRICVEC 12
using namespace gmx;
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);
- //}
+ F += t04;
+ Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
+ 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);
+ 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);
+ 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);
+ 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);
+ 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 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);
- }
- //}
+ double t0 = 0, t1 = 0, t2 = 0;
+ for (unsigned int i = 0; i < b.size(); i++) {
+ t0 = static_cast< double >(b[i][0]);
+ t1 = static_cast< double >(b[i][1]);
+ t2 = static_cast< double >(b[i][2]);
+ b[i][0] = static_cast< float >((t2 + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (t1 + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (t0 + x));
+ b[i][1] = static_cast< float >((t1 + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (t2 + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (t0 + x));
+ b[i][2] = static_cast< float >(cos(A) * cos(B) * (t2 + z) - sin(B) * (t0 + x) + cos(B) * sin(A) * (t1 + y));
+ }
}
-void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< int, int > > pairs) {
+void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< unsigned int, unsigned int > > pairs) {
midA[0] = 0;
midA[1] = 0;
midA[2] = 0;
midB[1] = 0;
midB[2] = 0;
- for (int i = 0; i < pairs.size(); i++) {
+ for (unsigned int i = 0; i < pairs.size(); i++) {
rvec_inc(midA, a[pairs[i].first]);
rvec_inc(midB, b[pairs[i].second]);
}
midB[2] /= pairs.size();
}
-void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< int, int > > pairs, double FtCnst) {
+void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned 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);
+ double x = static_cast< double >(ma[0]), y = static_cast< double >(ma[1]), z = static_cast< double >(ma[2]), A = 0, B = 0, C = 0;
+ for (unsigned int i = 0; i < pairs.size(); i++) {
+ f1 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
+ static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y, static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
}
if (FtCnst == 0) {
FtCnst = 0.00001;
} 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);
+ f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0; l = 1;
+ for (unsigned int i = 0; i < pairs.size(); i++) {
+ searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc,
+ static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
+ static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y,
+ static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
}
- while (true) {
+ while (f2 != f1) {
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);
+ for (unsigned int i = 0; i < pairs.size(); i++) {
+ f2 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
+ static_cast< double >(b[pairs[i].second][0]) + (x - l * fx), static_cast< double >(b[pairs[i].second][1]) + (y - l * fy),
+ static_cast< double >(b[pairs[i].second][2]) + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
}
count++;
if (f2 < f1) {
fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
break;
} else {
- l *= 0.85;
+ //l *= 0.85;
+ l *= 0.50;
+ /*if (count >= 14) {
+ std::cout << count << " " << l << "\n";
+ }*/
+ if (l == DBL_MIN) { /*DBL_TRUE_MIN*/
+ break;
+ }
}
}
count++;
+ if (count % 100000 == 0) {
+ std::cout << "+100k\n";
+ }
}
ApplyFit(b, x, y, z, A, B, C);
}
bool yep;
};
-void make_graph(int mgwi_natoms, std::vector < RVec > mgwi_x, std::vector< std::vector< node > > &mgwi_graph)
+void make_graph(unsigned long mgwi_natoms, std::vector < RVec > mgwi_x, std::vector< std::vector< node > > &mgwi_graph)
{
mgwi_graph.resize(mgwi_natoms);
- for (int i = 0; i < mgwi_natoms; i++) {
+ for (unsigned 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++) {
+ for (unsigned int i = 0; i < mgwi_natoms; i++) {
+ for (unsigned int j = 0; j < mgwi_natoms; j++) {
rvec_sub(mgwi_x[i].as_vec(), mgwi_x[j].as_vec(), mgwi_graph[i][j].r.as_vec());
mgwi_graph[i][j].n = 0;
}
}
}
-void update_graph(std::vector< std::vector< node > > &ugwi_graph, std::vector < RVec > ugwi_x, long double ugwi_epsi) {
+void update_graph(std::vector< std::vector< node > > &ugwi_graph, std::vector < 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++) {
+ for (unsigned int i = 0; i < ugwi_graph.size(); i++) {
+ for (unsigned int j = i; j < ugwi_graph.size(); j++) {
rvec_sub(ugwi_x[i].as_vec(), ugwi_x[j].as_vec(), ugwi_temp.as_vec());
rvec_dec(ugwi_temp.as_vec(), ugwi_graph[i][j].r.as_vec());
- if (norm(ugwi_temp) <= ugwi_epsi) {
+ if (static_cast< double >(norm(ugwi_temp)) <= ugwi_epsi) {
if (i == j) {
ugwi_graph[i][j].n++;
}
}
}
-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++) {
+void check_domains(double cd_delta, int cd_frames, std::vector< std::vector< std::vector< node > > > &cd_graph) {
+ for (unsigned int k = 0; k < cd_graph.size(); k++) {
+ for (unsigned int i = 0; i < cd_graph[1].size(); i++) {
+ for (unsigned int j = 0; j < cd_graph[1].size(); j++) {
if (cd_graph[k][i][j].n >= cd_frames * cd_delta) {
cd_graph[k][i][j].yep = true;
}
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++) {
+ for (unsigned int i = 0; i < fds_graph.size(); i++) {
+ fds_domsizes[i].resize(fds_graph[1].size(), 0);
+ for (unsigned int j = 0; j < fds_graph[1].size(); j++) {
+ for (unsigned int k = 0; k < fds_graph[1].size(); k++) {
if (fds_graph[i][j][k].yep) {
fds_domsizes[i][j]++;
}
}
}
-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();
- for (int i = 0; i < gmd_for1; i++) {
- for (int j = 0; j < gmd_for2; j++) {
+void get_maxsized_domain(std::vector< unsigned int > &gmd_max_d, std::vector< std::vector< std::vector< node > > > gmd_graph, std::vector< std::vector< int > > gmd_domsizes) {
+ unsigned int gmd_number1 = 0, gmd_number2 = 0;
+ for (unsigned int i = 0; i < gmd_domsizes.size(); i++) {
+ for (unsigned int j = 0; j < gmd_domsizes[0].size(); j++) {
if (gmd_domsizes[i][j] >= gmd_domsizes[gmd_number1][gmd_number2]) {
gmd_number1 = i;
gmd_number2 = j;
}
}
gmd_max_d.resize(0);
- int gmd_for3 = gmd_graph[gmd_number1][gmd_number2].size();
- for (int i = 0; i < gmd_for3; i++) {
+ for (unsigned int i = 0; i < gmd_graph[gmd_number1][gmd_number2].size(); i++) {
if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
gmd_max_d.push_back(i);
}
}
}
-void get_min_domain(std::vector< int > &gmd_min_d, std::vector< std::vector< std::vector< node > > > gmd_graph, std::vector< std::vector< int > > gmd_domsizes, int gmd_min_dom_size) {
- int gmd_number1 = 0, gmd_number2 = 0;
- int gmd_for1 = gmd_domsizes.size(), gmd_for2 = gmd_domsizes[0].size();
- for (int i = 0; i < gmd_for1; i++) {
- for (int j = 0; j < gmd_for2; j++) {
+void get_min_domain(std::vector< unsigned int > &gmd_min_d, std::vector< std::vector< std::vector< node > > > gmd_graph, std::vector< std::vector< int > > gmd_domsizes, int gmd_min_dom_size) {
+ unsigned int gmd_number1 = 0, gmd_number2 = 0;
+ for (unsigned int i = 0; i < gmd_domsizes.size(); i++) {
+ for (unsigned int j = 0; j < gmd_domsizes[0].size(); j++) {
if (gmd_domsizes[gmd_number1][gmd_number2] < gmd_min_dom_size && gmd_domsizes[i][j] >= gmd_min_dom_size) {
gmd_number1 = i;
gmd_number2 = j;
}
}
gmd_min_d.resize(0);
- int gmd_for3 = gmd_graph[gmd_number1][gmd_number2].size();
- for (int i = 0; i < gmd_for3; i++) {
+ for (unsigned int i = 0; i < gmd_graph[gmd_number1][gmd_number2].size(); i++) {
if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
gmd_min_d.push_back(i);
}
}
}
-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++) {
+void delete_domain_from_graph(std::vector< std::vector< std::vector< node > > > &ddf_graph, std::vector< unsigned int > ddf_domain) {
+ for (unsigned int i = 0; i < ddf_domain.size(); i++) {
+ for (unsigned int k = 0; k < ddf_graph.size(); k++) {
+ for (unsigned int j = 0; j < ddf_graph[1].size(); j++) {
if (ddf_graph[k][ddf_domain[i]][j].yep) {
ddf_graph[k][ddf_domain[i]][j].yep = false;
}
}
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++) {
+ for (unsigned int i = 0; i < cd_domsizes.size(); i++) {
+ for (unsigned int j = 0; j < cd_domsizes[0].size(); j++) {
if (cd_domsizes[i][j] >= cd_domain_min_size) {
return true;
}
return false;
}
-void print_domains(std::vector< std::vector< int > > pd_domains, std::vector< int > index, std::string fnNdx_, int dms, double epsi, double delta) {
+void print_domains(std::vector< std::vector< unsigned int > > pd_domains, std::vector< int > index, std::string fnNdx_, int dms, double epsi, double delta) {
FILE *fpNdx_;
fpNdx_ = std::fopen(fnNdx_.c_str(), "w+");
int write_count;
- for (int i = 0; i < pd_domains.size(); i++) {
+ for (unsigned int i = 0; i < pd_domains.size(); i++) {
std::fprintf(fpNdx_, "[domain_№_%3d_%3d_%4.3f_%4.3f]\n", i + 1, dms, epsi, delta);
write_count = 0;
- for (int j = 0; j < pd_domains[i].size(); j++) {
+ for (unsigned int j = 0; j < pd_domains[i].size(); j++) {
write_count++;
if (write_count > 20) {
write_count -= 20;
virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top);
- virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
- const t_trxframe &fr);
-
//! Call for each frame of the trajectory
// virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
std::string fnNdx_;
- std::vector< std::vector< std::vector< node > > > graph;
+ //std::vector< std::vector< std::vector< node > > > graph;
+ std::vector< std::vector< std::vector< std::vector< node > > > > graph;
- std::vector< std::vector< int > > domains;
+
+ std::vector< std::vector< unsigned int > > domains;
std::vector< std::vector< int > > domsizes;
std::vector< int > index;
std::vector< int > numbers;
- std::vector< std::vector < RVec > > trajectory;
+ std::vector< RVec > trajectory;
+ std::vector< RVec > reference;
Selection selec;
int frames = 0;
- int domain_min_size = 5; // selectable
+ int domain_min_size = 4; // selectable
- std::vector< std::vector< std::pair< int, int > > > w_rls2;
- int bone;
- double delta = 0.90; // selectable
- double epsi = 0.15; // selectable
+ std::vector< std::vector< std::pair< unsigned int, unsigned int > > > w_rls2;
+ unsigned int bone;
+ double delta = 0.95; // selectable
+ double epsi = 0.10; // selectable
- int domains_ePBC;
int DomainSearchingAlgorythm = 0; // selectable
// Copy and assign disallowed by base.
};
// Add option for domain min size constant
options->addOption(gmx::IntegerOption("dms")
.store(&domain_min_size)
- .description("minimum domain size"));
+ .description("minimum domain size, should be >= 4"));
// Add option for Domain's Searching Algorythm
options->addOption(gmx::IntegerOption("DSA")
.store(&DomainSearchingAlgorythm)
.description("domain membership probability"));
// Control input settings
settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
+ settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
settings->setPBC(true);
}
void
Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top)
-{
-}
-
-void
-Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
- const t_trxframe &fr)
{
index.resize(0);
for (ArrayRef< const int >::iterator ai = selec.atomIndices().begin(); (ai < selec.atomIndices().end()); ai++) {
index.push_back(*ai);
}
- trajectory.resize(2);
- trajectory[0].resize(selec.atomCount());
- for (int i = 0; i < selec.atomCount(); i++) {
- trajectory[0][i] = fr.x[index[i]];
+ reference.resize(0);
+ if (top.hasFullTopology()) {
+ for (unsigned int i = 0; i < index.size(); i++) {
+ reference.push_back(top.x().at(index[i]));
+ }
}
- bone = index.size() - domain_min_size + 1;
- graph.resize(bone);
+ graph.resize(0);
+
+ bone = static_cast< unsigned int >(index.size() - static_cast< unsigned int >(domain_min_size) + 1);
+ graph[0].resize(bone);
w_rls2.resize(bone);
- for (int i = 0; i < bone; i++) {
+ for (unsigned int i = 0; i < bone; i++) {
w_rls2[i].resize(0);
- for (int j = 0; j < domain_min_size; j++) {
+ for (unsigned int j = 0; j < static_cast< unsigned int >(domain_min_size); j++) {
w_rls2[i].push_back(std::make_pair(j + i, j + i));
}
- make_graph(index.size(), trajectory[0], graph[i]);
+ make_graph(index.size(), reference, graph[0][i]);
}
- trajectory[1].resize(index.size());
+
+ trajectory.resize(index.size());
}
void
-Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
- TrajectoryAnalysisModuleData *pdata)
+Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
{
- for (int i = 0; i < index.size(); i++) {
- trajectory[1][i] = fr.x[index[i]];
+ for (unsigned int i = 0; i < index.size(); i++) {
+ trajectory[i] = fr.x[index[i]];
}
frames++;
- //#pragma omp parallel
- //{
- //#pragma omp for schedule(dynamic)
- for (int j = 0; j < bone; j++) {
- std::vector < RVec > temp = trajectory[1];
- MyFitNew(trajectory[0], temp, w_rls2[j], 0);
- update_graph(graph[j], temp, epsi);
- }
- //}
+ #pragma omp parallel for ordered schedule(dynamic) shared(w_rls2, graph) firstprivate(trajectory, reference, epsi)
+ for (unsigned int j = 0; j < bone; j++) {
+ std::vector < RVec > temp = trajectory;
+ MyFitNew(reference, temp, w_rls2[j], 0);
+ update_graph(graph[0][j], temp, epsi);
+ }
+ #pragma omp barrier
std::cout << "frame: " << frames << " analyzed\n";
}
frames -= 1;
std::cout << "final cheking\n";
- check_domains(delta, frames, graph);
+ check_domains(delta, frames, graph[0]);
std::cout << "finding domains' sizes\n";
- find_domain_sizes(graph, domsizes);
+ find_domain_sizes(graph[0], domsizes);
std::cout << "finding domains\n";
- std::vector< int > a;
+ std::vector< unsigned int > a;
a.resize(0);
while (check_domsizes(domsizes, domain_min_size)) {
domains.push_back(a);
if (DomainSearchingAlgorythm == 0) {
- get_maxsized_domain(domains.back(), graph, domsizes);
+ get_maxsized_domain(domains.back(), graph[0], domsizes);
} else if (DomainSearchingAlgorythm == 1) {
- get_min_domain(domains.back(), graph, domsizes, domain_min_size);
+ get_min_domain(domains.back(), graph[0], domsizes, domain_min_size);
}
std::cout << "new domain: " << domains.back().size() << " atoms\n";
- delete_domain_from_graph(graph, domains.back());
+ delete_domain_from_graph(graph[0], domains.back());
domsizes.resize(0);
- find_domain_sizes(graph, domsizes);
+ find_domain_sizes(graph[0], domsizes);
}
}
biod.pnpi.spb.ru / alexxy / gromacs-domains.git / commitdiff