#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;
-};
-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 <string>
-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.
*
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.
// 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);
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<int>::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";
+
}
biod.pnpi.spb.ru / alexxy / gromacs-testing.git / commitdiff