#include <gromacs/utility/smalloc.h>
#include <gromacs/math/do_fit.h>
+#include "fittingn.cpp"
+
#define MAX_NTRICVEC 12
using namespace gmx;
}
}
+void update_graph_v2(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++) {
+ rvec_sub(ugwi_x[i].as_vec(), ugwi_x[j].as_vec(), 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++) {
std::fprintf(fpNdx_, "\n");
}
std::fprintf(fpNdx_, "%5d ", index[pd_domains[i][j]] + 1);
+ std::printf("%5d ", index[pd_domains[i][j]] + 1);
}
+ std::printf("\n\n");
std::fprintf(fpNdx_,"\n\n");
}
std::fprintf(fpNdx_,"\n");
int domain_min_size = 5; // selectable
real **w_rls;
+ std::vector< std::vector< std::pair< int, int > > > w_rls2;
int bone;
double delta = 0.90; // selectable
double epsi = 0.15; // selectable
Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top)
{
- domains_ePBC = top.ePBC();
}
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<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
+ for (ArrayRef< const int >::iterator ai = selec.atomIndices().begin(); (ai < selec.atomIndices().end()); ai++) {
index.push_back(*ai);
}
trajectory.resize(2);
trajectory[0][i] = fr.x[index[i]];
}
+
+
+
bone = index.size() - domain_min_size + 1;
graph.resize(bone);
snew(w_rls, bone);
+ w_rls2.resize(bone);
for (int i = 0; i < bone; i++) {
snew(w_rls[i], index.size());
for (int j = 0; j < index.size(); j++) {
w_rls[i][j] = 0;
}
}
+ w_rls2[i].resize(0);
+ for (int j = 0; j < domain_min_size; j++) {
+ w_rls2[i].push_back(std::make_pair(j + i, j + i));
+ }
rvec *etalon;
snew(etalon, index.size());
for (int j = 0; j < index.size(); j++) {
for (int i = 0; i < index.size(); i++) {
trajectory[1][i] = fr.x[index[i]];
}
+ std::vector < RVec > temp;
frames++;
#pragma omp parallel
{
#pragma omp for schedule(dynamic)
for (int j = 0; j < bone; j++) {
- rvec *etalon, *traj;
+ /*rvec *etalon, *traj;
snew(etalon, index.size());
for (int k = 0; k < index.size(); k++) {
copy_rvec(trajectory[0][k].as_vec(), etalon[k]);
do_fit(index.size(), w_rls[j], etalon, traj);
update_graph(graph[j], traj, epsi);
sfree(etalon);
- sfree(traj);
+ sfree(traj);*/
+
+ std::vector < RVec > temp = trajectory[1];
+ MyFitNew(trajectory[0], temp, w_rls2[j]);
+ update_graph_v2(graph[j], temp, epsi);
}
}
std::cout << "frame: " << frames << " analyzed\n";
{
return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+