From a0d5e59ba2d0ac98bbce9e56569a3848a24418b5 Mon Sep 17 00:00:00 2001 From: Anatoly Titov Date: Tue, 25 Apr 2017 11:40:45 +0300 Subject: [PATCH] fresh beginning --- src/spacetimecorr.cpp | 347 ++++++++++++++---------------------------- 1 file changed, 113 insertions(+), 234 deletions(-) diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp index ac6e380..3242d27 100644 --- a/src/spacetimecorr.cpp +++ b/src/spacetimecorr.cpp @@ -54,152 +54,40 @@ #include #include -#define MAX_NTRICVEC 12 -using namespace gmx; -using gmx::RVec; -struct node { - short int n; - RVec r; - bool yep; -}; -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; - } - } -} -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++; - } - } - } - } -} -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; - } - } - } - } -} -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]++; - } - } - } - } -} +#include -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); - } - } -} +#include "gromacs/analysisdata/analysisdata.h" +#include "gromacs/analysisdata/modules/average.h" +#include "gromacs/analysisdata/modules/histogram.h" +#include "gromacs/analysisdata/modules/plot.h" +#include "gromacs/math/vec.h" +#include "gromacs/options/basicoptions.h" +#include "gromacs/options/filenameoption.h" +#include "gromacs/options/ioptionscontainer.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/selection/selection.h" +#include "gromacs/selection/selectionoption.h" +#include "gromacs/trajectory/trajectoryframe.h" +#include "gromacs/trajectoryanalysis/analysissettings.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/exceptions.h" +#include "gromacs/utility/stringutil.h" -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; - } - } - } - } -} -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; - } - } - } - return false; -} -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"); - } - std::fprintf(fpNdx_, "%5d ", index[pd_domains[i][j]] + 1); - } - std::fprintf(fpNdx_,"\n\n"); - } - std::fprintf(fpNdx_,"\n"); - std::fclose(fpNdx_); -} +using namespace gmx; + +using gmx::RVec; + /*! \brief * Class used to compute free volume in a simulations box. * @@ -242,23 +130,19 @@ class Domains : public TrajectoryAnalysisModule private: std::string fnNdx_; - - std::vector< std::vector< std::vector< node > > > graph; + SelectionList sel_; + Selection selec; std::vector< std::vector< int > > domains; - std::vector< std::vector< int > > domsizes; + std::vector< std::vector< RVec > > trajectory; + std::vector< std::vector< RVec > > frankenstein_trajectory; + 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 - + int basic_frame = 0; 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; // Copy and assign disallowed by base. @@ -282,25 +166,16 @@ Domains::initOptions(IOptionsContainer *options, // 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") + options->addOption(SelectionOption("select_type") .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("domains") .description("Index file from the domains")); - // 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_domains").storeVector(&sel_) + .required().dynamicMask().multiValue() + .description("Domains to form rigid skeleton")); // Control input settings settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC); settings->setPBC(true); @@ -317,119 +192,123 @@ void Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr) { - t_pbc pbc; - t_pbc *ppbc = settings.hasPBC() ? &pbc : NULL; - matrix boxx; - copy_mat(fr.box, boxx); - if (ppbc != NULL) { - set_pbc(ppbc, domains_ePBC, boxx); - } ConstArrayRef< int > atomind = selec.atomIndices(); index.resize(0); for (ConstArrayRef::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) { index.push_back(*ai); } - trajectory.resize(2); + + //checkSelections(sel_); + domains.resize(sel_.size()); + for (int i = 0; i < sel_.size(); i++) { + for (int j = 0; j < sel_[i].size(); j++) { //distance / modules / trajectoryanalyziz + domains[i].push_back(sel_[i][j]); + } + } + + trajectory.resize(1); trajectory[0].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[1].resize(index.size()); + frames++; } void Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata) { + std::vector< std::vector< int > > domains_local; + std::vector< bool > flags; + + flags.resize(index.size(), true); + domains_local = domains; // технически оно за О(1) должно делаться, но не точно + + 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]]; + trajectory.back()[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]); + for (int i = 0; i < domains.size(); i++) { + for (int j = 0; j < domains[i].size(); j++) { + flags[domains[i][j]] = false; + } + } + + for (int i = 0; i < index.size(); i++) { + if (flags[i]) { + int a, b, dist = 9001; + rvec temp; + for (int j = 0; j < domains.size(); j++) { + for (int k = 0; k < domains[j].size(); k++) { + rvec_sub(trajectory.back()[i], trajectory.back()[domains[j][k]], temp); + if (norm(temp) <= dist) { + dist = norm(temp); + a = j; + b = 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); + domains_local[a].push_back(i); + flags[i] = false; + } + } + + snew(w_rls, domains_local.size()); + for (int i = 0; i < domains_local.size(); i++) { + snew(w_rls[i], index.size()); + for (int j = 0; j < index.size(); j++) { + w_rls[i][j] = 0; + } + for (int j = 0; j < domains_local[i].size(); j++) { + w_rls[i][domains_local[i][j]] = 1; + } + } + + frankenstein_trajectory.resize(trajectory.size()); + frankenstein_trajectory.back().resize(index.size()); + + for (int i = 0; i < domains_local.size(); i++) { + rvec *basic, *traj; + snew(basic, index.size()); + for (int k = 0; k < index.size(); k++) { + copy_rvec(trajectory[basic_frame][k].as_vec(), basic[k]); + } + snew(traj, index.size()); + for (int k = 0; k < index.size(); k++) { + copy_rvec(trajectory.back()[k].as_vec(), traj[k]); } + reset_x(index.size(), NULL, index.size(), NULL, basic, w_rls[i]); + reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[i]); + do_fit(index.size(), w_rls[i], basic, traj); + + for (int j = 0; j < domains_local[i].size(); j++) { + frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]]; + } + + sfree(basic); + sfree(traj); } - std::cout << "frame: " << frames << " analyzed\n"; + + frames++; } //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' void Domains::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); - } - for (int i = 0; i < bone; i++) { - sfree(w_rls[i]); - } - sfree(w_rls); + } void Domains::writeOutput() { - std::cout << "making output file\n"; - print_domains(domains, index, fnNdx_); // see function for details | numbers from index - std::cout << "\n END \n"; + } -- 2.22.0