*/
#include <iostream>
-#include <chrono>
-#include <cstdint>
+//#include <chrono>
+//#include <cstdint>
#include <cfloat>
#include <omp.h>
#include <gromacs/trajectoryanalysis.h>
-#include <gromacs/pbcutil/pbc.h>
-#include <gromacs/utility/smalloc.h>
-#include <gromacs/math/do_fit.h>
+//#include <gromacs/pbcutil/pbc.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"
+//#include "gromacs/topology/topology.h"
+//#include "gromacs/topology/index.h"
-//#include "/home/titov_ai/Develop/fittingLib/fittinglib.h"
-
-#define MAX_NTRICVEC 12
+// #define MAX_NTRICVEC 12
using namespace gmx;
for (unsigned int i = 0; i < pairs.size(); i++) {
midA += a[pairs[i].first];
midB += b[pairs[i].second];
- //rvec_inc(midA, a[pairs[i].first]);
- //rvec_inc(midB, b[pairs[i].second]);
}
midA[0] /= pairs.size();
midA[1] /= pairs.size();
RVec ma, mb;
CalcMid(a, b, ma, mb, pairs);
ma -= mb;
- //rvec_dec(ma, mb);
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]),
const int domainSearchAlgorythm, const int timeStepBetweenWindows,
const unsigned int sliceNum, const double epsilon, const double delta,
const std::string outPutFileName) {
- graph.resize(1);
- graph.front().resize(sliceNum);
- for (unsigned i = 0; i < sliceNum; i++) {
- setGraph(graph.front()[i], reference);
- }
+ graph.resize(0);
structIndex = index;
ref = reference;
if (epsilon != -1) { eps = epsilon; }
if (domainSearchAlgorythm != -1) { dsa = domainSearchAlgorythm; }
if (timeStepBetweenWindows != -1) { ts = timeStepBetweenWindows; }
if (outPutFileName != "") { outPut = outPutFileName; }
+ domsizes.resize(0);
+ domsizes.resize(sliceNum);
}
void update(const std::vector< std::vector< RVec > > frame, const int frameNumber) {
if (frameNumber % ts == 0) {
graph.resize(graph.size() + 1);
- graph.back().resize(graph.front().size());
+ graph.back().resize(structIndex.size() - static_cast< unsigned long >(dms) + 1);
for (unsigned i = 0; i < graph.front().size(); i++) {
setGraph(graph.back()[i], ref);
}
}
bool searchDomainSizes() {
- // it goes infinite without these 2 lines, curious even if it was resized before
- // resizing inside vectors with 0 is not enough
- domsizes.resize(0);
- domsizes.resize(graph.front().size());
bool flag = false;
for (unsigned int i = 0; i < graph.front().size(); i++) {
+ domsizes[i].resize(0);
domsizes[i].resize(graph.front()[i].size(), 0);
for (unsigned int j = 0; j < graph.front()[i].size(); j++) {
for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
if (domains.size() == 0) {
return;
}
- FILE *ndxFile, *slFile;
+ FILE *ndxFile, *slFile;
ndxFile = std::fopen(outPut.c_str(), "a+");
- slFile = std::fopen(("SelectionList-" + outPut.substr(0, outPut.size() - 4)).c_str(), "a+");
+ slFile = std::fopen(("SelectionList-" + outPut.substr(0, outPut.size() - 4)).c_str(), "a+");
int writeCount;
for (unsigned int i = 0; i < domains.size(); i++) {
std::fprintf(ndxFile, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", currentFrame - window, i + 1, dms, eps, dlt);
}
};
-struct node {
- short int n;
- RVec r;
- bool yep;
-};
-
-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 (unsigned int i = 0; i < mgwi_natoms; i++) {
- mgwi_graph[i].resize(mgwi_natoms);
- }
- 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) {
- RVec ugwi_temp;
- 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 (static_cast< double >(norm(ugwi_temp)) <= ugwi_epsi) {
- if (i == j) {
- ugwi_graph[i][j].n++;
- }
- else {
- ugwi_graph[i][j].n++;
- ugwi_graph[j][i].n++;
- }
- }
- }
- }
-}
-
-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;
- }
- else {
- cd_graph[k][i][j].yep = false;
- }
- }
- }
- }
-}
-
-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());
- 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< 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);
- 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< 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;
- }
- if (gmd_domsizes[i][j] <= gmd_domsizes[gmd_number1][gmd_number2] && gmd_domsizes[i][j] >= gmd_min_dom_size) {
- gmd_number1 = i;
- gmd_number2 = j;
- }
- }
- }
- gmd_min_d.resize(0);
- 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< 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;
- }
- if (ddf_graph[k][j][ddf_domain[i]].yep) {
- ddf_graph[k][j][ddf_domain[i]].yep = false;
- }
- }
- }
- }
-}
-
-bool check_domsizes(std::vector< std::vector< int > > cd_domsizes, int cd_domain_min_size) {
- 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< unsigned int > > pd_domains, std::vector< int > index, std::string fnNdx_, int dms, double epsi, double delta, int window, int st_pos) {
- FILE *fpNdx_, *fpNdx2_;
- fpNdx_ = std::fopen(fnNdx_.c_str(), "a+");
- fpNdx2_ = std::fopen(("SelectionList-" + fnNdx_.substr(0, fnNdx_.size() - 4)).c_str(), "a+");
- int write_count;
- for (unsigned int i1 = 0; i1 < pd_domains.size(); i1++) {
- std::fprintf(fpNdx_, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", static_cast< int >(st_pos) * window / 10, i1 + 1, dms, epsi, delta);
- std::fprintf(fpNdx2_, "group %cdomain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', static_cast< int >(st_pos) * window / 10, i1 + 1, dms, epsi, delta, '"');
- write_count = 0;
- for (unsigned int j = 0; j < pd_domains[i1].size(); j++) {
- write_count++;
- if (write_count > 20) {
- write_count -= 20;
- std::fprintf(fpNdx_, "\n");
- }
- std::fprintf(fpNdx_, "%5d ", index[pd_domains[i1][j]] + 1);
- std::printf("%5d ", index[pd_domains[i1][j]] + 1);
- }
- std::printf("\n\n");
- std::fprintf(fpNdx_,"\n\n");
- }
- std::fprintf(fpNdx_,"\n");
- std::fclose(fpNdx_);
- //std::fprintf(fpNdx2_,"\n");
- std::fclose(fpNdx2_);
-}
-
/*! \brief
* \ingroup module_trajectoryanalysis
*/
private:
std::string fnNdx_;
-
- std::vector< std::vector< std::vector< std::vector< node > > > > graph;
-
domainsType testSubject;
-
-
- std::vector< std::vector< unsigned int > > domains;
- std::vector< std::vector< int > > domsizes;
-
std::vector< int > index;
std::vector< RVec > trajectory;
std::vector< RVec > reference;
Selection selec;
- int frames = 0;
- int window = -1; // selectable
- int domain_min_size = -1; // selectable
-
- std::vector< std::vector< std::pair< unsigned int, unsigned int > > > w_rls2;
+ std::vector< std::vector< std::pair< unsigned int, unsigned int > > > fitPairs;
unsigned int bone;
- double delta = -1; // selectable
- double epsi = -1; // selectable
- int twStep = -1;
-
- int DomainSearchingAlgorythm = -1; // selectable
+ int window = -1; // selectable
+ int domain_min_size = -1; // selectable
+ double delta = -1; // selectable
+ double epsi = -1; // selectable
+ int twStep = -1; // selectable
+ int DomainSearchingAlgorythm = -1; // selectable
// Copy and assign disallowed by base.
};
}
bone = static_cast< unsigned int >(index.size() - static_cast< unsigned int >(domain_min_size) + 1);
- graph.resize(1);
- graph.front().resize(bone);
-
- w_rls2.resize(bone);
+ fitPairs.resize(bone);
for (unsigned int i = 0; i < bone; i++) {
- w_rls2[i].resize(0);
+ fitPairs[i].resize(0);
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));
+ fitPairs[i].push_back(std::make_pair(j + i, j + i));
}
- //make_graph(index.size(), reference, graph.front()[i]);
}
trajectory.resize(index.size());
void
Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
{
- /*if (frnr % twStep == 0) {
- graph.resize(graph.size() + 1);
- graph.back().resize(bone);
- for (unsigned int i = 0; i < bone; i++) {
- make_graph(index.size(), reference, graph.back()[i]);
- }
- }
-
- int temp = graph.size() - window / twStep;
-
- // most probably == 0
- if (temp > 0) {
- domains.resize(0);
- std::vector< unsigned int > a;
- a.resize(0);
- domsizes.resize(0);
- check_domains(delta, window, graph[0]);
- std::cout << "finding domains' sizes\n";
- find_domain_sizes(graph[0], domsizes);
- while (check_domsizes(domsizes, domain_min_size)) {
- domains.push_back(a);
- if (DomainSearchingAlgorythm == 0) {
- get_maxsized_domain(domains.back(), graph[0], domsizes);
- } else if (DomainSearchingAlgorythm == 1) {
- 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[0], domains.back());
- domsizes.resize(0);
- find_domain_sizes(graph[0], domsizes);
- }
- graph.erase(graph.begin());
- std::cout << (frnr + 1) / twStep - (window / twStep)<< " window's domains have been evaluated\n";
- print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta, window, (frnr + 1) / twStep - (window / twStep)); // see function for details | numbers from index
- }*/
-
for (unsigned int i = 0; i < index.size(); i++) {
trajectory[i] = fr.x[index[i]];
}
- frames++;
-
- /*#pragma omp parallel for ordered schedule(dynamic) shared(w_rls2, graph) firstprivate(trajectory, reference, epsi, frnr, window)
- for (unsigned int j = 0; j < bone; j++) {
- std::vector < RVec > temp = trajectory;
- MyFitNew(reference, temp, w_rls2[j], 0);
- for (unsigned int i = 0; i < graph.size(); i++) {
- update_graph(graph[i][j], temp, epsi);
- }
- }
- #pragma omp barrier*/
std::vector< std::vector< RVec > > tsTemp;
+ tsTemp.resize(0);
tsTemp.resize(bone, trajectory);
- #pragma omp parallel for ordered schedule(dynamic) shared(w_rls2) firstprivate(trajectory, reference)
+ #pragma omp parallel for ordered schedule(dynamic) shared(fitPairs) firstprivate(reference)
for (unsigned int i = 0; i < bone; i++) {
- MyFitNew(reference, tsTemp[i], w_rls2[i], 0);
+ MyFitNew(reference, tsTemp[i], fitPairs[i], 0);
}
#pragma omp barrier
+
testSubject.update(tsTemp, frnr);
- //std::cout << "frame: " << frnr << " analyzed\n";
}
void
biod.pnpi.spb.ru / alexxy / gromacs-domains.git / commitdiff