2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Implements gmx::analysismodules::Freevolume.
39 * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
40 * \ingroup module_trajectoryanalysis
53 #include <gromacs/trajectoryanalysis.h>
54 #include <gromacs/utility/smalloc.h>
55 #include <gromacs/math/do_fit.h>
56 #include <gromacs/trajectoryanalysis/topologyinformation.h>
57 #include "gromacs/topology/topology.h"
58 #include "gromacs/topology/index.h"
64 const std::vector< std::vector< std::vector< double > > > read_correlation_matrix_file(const char* file_name, unsigned long size, unsigned int length)
67 std::vector< std::vector< std::vector< double > > > crrlts;
68 file = std::fopen(file_name, "r+");
69 std::vector< double > c;
71 std::vector< std::vector< double > > b;
73 crrlts.resize(length + 1, b);
75 std::cout << "reading correlations from a file:\n";
76 for (unsigned int i = 0; i < length + 1; i++) {
77 int t0 = 0, t1 = std::fscanf(file, "%d\n", &t0);
78 std::cout << t0 << " | ";
82 for (unsigned int j = 0; j < size; j++) {
83 for (unsigned int f = 0; f < size; f++) {
84 int t2 = std::fscanf(file, "%lf ", &crrlts[i][j][f]);
93 template< typename T >
94 unsigned int posSrch(std::vector< T > a, T b) {
95 for (unsigned i01 = 0; i01 < a.size(); i01++) {
104 void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std::string > resNames,
105 std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > rout,
106 const char* file_name01, const unsigned int tauD, const unsigned int window_size)
109 file01 = std::fopen(file_name01, "w+");
111 unsigned int f, pre = 0;
112 std::vector< std::tuple< int, int, std::vector< int > > > table;
114 for (unsigned int i01 = 0; i01 < rout.size(); i01++) {
116 for (unsigned int i02 = i01 + 1; i02 <rout.size(); i02++) {
117 if (rout[i02] == rout[i01]) {
124 std::fprintf(file01, "\n Starting time point = %d | correlations >= %0.2f | tau = %d | window = %d\n\n", i01, crl_border, tauD, window_size);
126 std::fprintf(file01, "\n Starting time point = [%d ; %d] | correlations >= %0.2f | tau = %d | window = %d\n\n", i01, f, crl_border, tauD, window_size);
128 for (unsigned int i02 = 0; i02 < rout[i01].size(); i02++) {
129 for (unsigned int i03 = 0; i03 < rout[i01][i02].size(); i03++) {
130 std::fprintf(file01, "cgo_arrow (id %3d), (id %3d), radius=0.075\n", indx[rout[i01][i02][i03].first] + 1, indx[rout[i01][i02][i03].second] + 1);
132 std::fprintf(file01, "\n");
135 for (unsigned int i02 = 0; i02 < rout[i01].size(); i02++) {
136 for (unsigned int i03 = 0; i03 < rout[i01][i02].size(); i03++) {
138 for (unsigned int i04 = 0; i04 < table.size(); i04++) {
139 if (std::get<0>(table[i04]) == indx[rout[i01][i02][i03].first]) {
140 std::get<1>(table[i04])++;
141 std::get<2>(table[i04]).push_back(indx[rout[i01][i02][i03].second]);
143 } else if (std::get<0>(table[i04]) == indx[rout[i01][i02][i03].second]) {
144 std::get<1>(table[i04])++;
145 std::get<2>(table[i04]).push_back(indx[rout[i01][i02][i03].first]);
150 std::vector< int > a;
152 a.push_back(indx[rout[i01][i02][i03].second]);
153 table.push_back(std::make_tuple(indx[rout[i01][i02][i03].first], 1, a));
155 a.push_back(indx[rout[i01][i02][i03].first]);
156 table.push_back(std::make_tuple(indx[rout[i01][i02][i03].second], 1, a));
161 for (unsigned int i02 = 0; i02 < table.size(); i02++) {
162 std::fprintf(file01, "residue %s connections %d | ", (resNames[posSrch(indx, std::get<0>(table[i02]))]).c_str(), std::get<1>(table[i02]));
163 for (unsigned int i03 = 0; i03 < std::get<2>(table[i02]).size(); i03++) {
164 std::fprintf(file01, "%s ", (resNames[posSrch(indx, std::get<2>(table[i02])[i03])]).c_str());
166 std::fprintf(file01, "\n");
169 std::vector< std::vector< std::pair< unsigned int, unsigned int > > > temp01, temp02;
170 std::vector< std::pair< unsigned int, unsigned int > > temp03, temp04, temp05;
175 for (unsigned int i02 = 0; i02 < temp01.size(); i02++) {
176 for (unsigned int i03 = 0; i03 < temp01[i02].size(); i03++) {
177 temp03.push_back(temp01[i02][i03]);
180 for (unsigned int i02 = 0; i02 < temp02.size(); i02++) {
181 for (unsigned int i03 = 0; i03 < temp02[i02].size(); i03++) {
182 temp04.push_back(temp02[i02][i03]);
185 std::sort(temp03.begin(), temp03.end());
186 std::sort(temp04.begin(), temp04.end());
188 std::set_difference(temp03.begin(), temp03.end(), temp04.begin(), temp04.end(), std::inserter(temp05, temp05.begin()));
189 std::fprintf(file01, "minus:\n");
190 for (unsigned int i02 = 0; i02 < temp05.size(); i02++) {
191 std::fprintf(file01, "cgo_arrow (id %3d), (id %3d), radius=0.075\n", indx[temp05[i02].first] + 1, indx[temp05[i02].second] + 1);
194 std::set_difference(temp04.begin(), temp04.end(), temp03.begin(), temp03.end(), std::inserter(temp05, temp05.begin()));
195 std::fprintf(file01, "plus:\n");
196 for (unsigned int i02 = 0; i02 < temp05.size(); i02++) {
197 std::fprintf(file01, "cgo_arrow (id %3d), (id %3d), radius=0.075\n", indx[temp05[i02].first] + 1, indx[temp05[i02].second] + 1);
206 bool isitsubset (std::vector< int > a, std::vector< int > b) {
210 std::sort(a.begin(), a.end());
211 std::sort(b.begin(), b.end());
213 for (unsigned int i = 0; i < a.size(); i++) {
225 void correlation_evaluation(const std::vector< RVec > ref, const std::vector< std::vector< RVec > > trajectory, const unsigned int tauE, const unsigned int window_size, const char* file_name) {
226 std::vector< std::vector< std::vector< double > > > crl;
227 std::vector< std::vector< RVec > > trj;
229 for (unsigned int istart = 0; istart < trajectory.size() - window_size - tauE; istart++) {
231 trj.resize(window_size + tauE);
232 for (unsigned int i = 0; i < window_size + tauE; i++) {
233 for (unsigned int j = 0; j < trajectory[i].size(); j++) {
234 trj[i].push_back(trajectory[i + istart][j]);
237 crl.resize(tauE + 1);
238 for (unsigned int i = 0; i < crl.size(); i++) {
239 crl[i].resize(trajectory.front().size());
240 for (unsigned int j = 0; j < trajectory.front().size(); j++) {
241 crl[i][j].resize(trajectory.front().size());
242 for (unsigned int k = 0; k < trajectory.front().size(); k++) {
247 #pragma omp parallel for ordered schedule(dynamic) shared(crl) firstprivate(trj, ref)
248 for (unsigned int i = 0; i <= tauE; i++) {
249 std::vector< std::vector< double > > a, b, c;
250 std::vector< double > d;
252 for (unsigned int j = 0; j < trajectory.front().size(); j++) {
256 for (unsigned int k = 0; k < trajectory.front().size(); k++) {
262 for (unsigned int j = 0; j < trj.size() - tauE; j++) {
263 for (unsigned int f1 = 0; f1 < trj.front().size(); f1++) {
264 for (unsigned int f2 = 0; f2 < trj.front().size(); f2++) {
266 temp1 = trj[j][f1] - ref[f1];
267 temp2 = trj[j + i][f2] - ref[f2];
268 a[f1][f2] += (static_cast<double>(temp1[0]) * static_cast<double>(temp2[0]) +
269 static_cast<double>(temp1[1]) * static_cast<double>(temp2[1]) +
270 static_cast<double>(temp1[2]) * static_cast<double>(temp2[2]));
271 b[f1][f2] += (static_cast<double>(temp1[0]) * static_cast<double>(temp1[0]) +
272 static_cast<double>(temp1[1]) * static_cast<double>(temp1[1]) +
273 static_cast<double>(temp1[2]) * static_cast<double>(temp1[2]));
274 c[f1][f2] += (static_cast<double>(temp2[0]) * static_cast<double>(temp2[0]) +
275 static_cast<double>(temp2[1]) * static_cast<double>(temp2[1]) +
276 static_cast<double>(temp2[2]) * static_cast<double>(temp2[2]));
280 for (unsigned int j = 0; j < trj.front().size(); j++) {
281 for (unsigned int f = 0; f < trj.front().size(); f++) {
282 crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
288 for (unsigned int i1 = 0; i1 < crl.size(); i1++) {
289 for (unsigned int i2 = 0; i2 < crl[i1].size(); i2++) {
290 for (unsigned int i3 = 0; i3 < crl[i1][i2].size(); i3++) {
291 crl[i1][i2][i3] = std::round(crl[i1][i2][i3] * 10000) / 10000;
297 file = std::fopen(file_name, "a");
298 for (unsigned int i = 0; i < crl.size(); i++) {
299 std::fprintf(file, "%d %d\n", istart, i);
300 for (unsigned int j = 0; j < crl[i].size(); j++) {
301 for (unsigned int f = 0; f < crl[i][j].size(); f++) {
302 std::fprintf(file, "%.6f ", crl[i][j][f]); //~16
304 std::fprintf(file, "\n");
308 if (istart % 10 == 0 && istart > 0) {
311 std::cout << " | done " << istart << " | ";
316 void graph_calculation( std::vector< std::vector< std::pair< double, long int > > > &graph, std::vector< std::vector< unsigned int > > &s_graph,
317 std::vector< std::vector< std::pair< unsigned int, unsigned int > > > &s_graph_rbr,
318 const std::vector< RVec > ref,
319 const std::vector< std::vector< std::vector<double > > > crl,
320 const double crl_b, const double e_rad, const unsigned int tauB, const unsigned int tauE) {
321 graph.resize(ref.size());
322 for (unsigned int i = 0; i < ref.size(); i++) {
323 graph[i].resize(ref.size(), std::make_pair(0, -1));
326 for (unsigned int i = tauB; i <= tauE; i++) {
327 for (unsigned int j = 0; j < ref.size(); j++) {
328 for (unsigned int f = j; f < ref.size(); f++) {
329 temp = ref[i] - ref[f];
330 if (std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j])) >= crl_b &&
331 static_cast<double>(norm(temp)) <= e_rad && std::abs(graph[j][f].first) < std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j]))) {
332 if (std::abs(crl[i][j][f]) > std::abs(crl[i][f][j])) {
333 graph[j][f].first = crl[i][j][f];
335 graph[j][f].first = crl[i][f][j];
337 graph[j][f].second = i;
342 std::vector< bool > graph_flags;
343 graph_flags.resize(ref.size(), true);
344 std::vector< unsigned int > a;
346 std::vector< std::pair< unsigned int, unsigned int > > b;
348 std::vector< unsigned int > que1, que2, que3;
349 for (unsigned int i = 0; i < ref.size(); i++) {
350 if (graph_flags[i]) {
351 s_graph.push_back(a);
352 s_graph_rbr.push_back(b);
358 graph_flags[i] = false;
359 while(que1.size() > 0) {
361 for (unsigned int k = 0; k < que1.size(); k++) {
362 for (unsigned int j = 0; j < ref.size(); j++) {
363 if (graph[que1[k]][j].second > -1 && graph_flags[j]) {
365 graph_flags[j] = false;
370 for (unsigned int j = 0; j < que2.size(); j++) {
371 que3.push_back(que2[j]);
374 s_graph.back() = que3;
375 for (unsigned int j = 0; j < que3.size(); j++) {
376 for (unsigned int k = 0; k < ref.size(); k++) {
377 if (graph[que3[j]][k].second > -1) {
378 s_graph_rbr.back().push_back(std::make_pair(que3[j], k));
386 bool myfunction (const std::pair< int, double > i, const std::pair< int, double > j) {
387 return i.second < j.second;
390 void graph_back_bone_evaluation(std::vector< std::vector < std::pair< unsigned int, unsigned int > > > &rout_n, const unsigned long indxSize,
391 const std::vector< std::vector< std::pair< double, long > > > graph, const std::vector< std::vector< unsigned int > > s_graph,
392 const std::vector< std::vector< std::pair< unsigned int, unsigned int > > > s_graph_rbr) {
393 std::vector< double > key;
394 std::vector< long > path;
395 std::vector< std::pair< unsigned int, double > > que;
396 std::vector< std::pair< unsigned int, unsigned int > > a;
398 for (unsigned int i = 0; i < s_graph.size(); i++) {
403 if (s_graph[i].size() > 2) {
404 key.resize(indxSize, 2);
405 path.resize(indxSize, -1);
406 key[s_graph[i][0]] = 0;
407 for (unsigned int j = 0; j < s_graph[i].size(); j++) {
408 que.push_back(std::make_pair(s_graph[i][j], key[s_graph[i][j]]));
410 std::sort(que.begin(), que.end(), myfunction);
411 while (!que.empty()) {
413 que.erase(que.begin());
414 for (unsigned int j = 0; j < s_graph_rbr[i].size(); j++) {
416 if (s_graph_rbr[i][j].first == v) {
417 u = s_graph_rbr[i][j].second;
418 } else if (s_graph_rbr[i][j].second == v) {
419 u = s_graph_rbr[i][j].first;
422 unsigned long pos = 0;
423 for (unsigned long k = 0; k < que.size(); k++) {
424 if (que[k].first == u) {
430 if (flag && key[static_cast< unsigned long >(u)] > 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first)) {
431 path[static_cast< unsigned long >(u)] = v;
432 key[static_cast< unsigned long >(u)] = 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first);
433 que[pos].second = key[static_cast< unsigned long >(u)];
434 sort(que.begin(), que.end(), myfunction);
440 for (unsigned int j = 0; j < indxSize; j++) {
442 rout_n.back().push_back(std::make_pair(j, path[j]));
449 gmx::RVec evaluate_com(std::vector< RVec > frame) {
454 for (unsigned int i = 0; i < frame.size(); i++) {
457 temp[0] /= frame.size();
458 temp[1] /= frame.size();
459 temp[2] /= frame.size();
464 * \ingroup module_trajectoryanalysis
467 class SpaceTimeCorr : public TrajectoryAnalysisModule
472 virtual ~SpaceTimeCorr();
474 //! Set the options and setting
475 virtual void initOptions(IOptionsContainer *options,
476 TrajectoryAnalysisSettings *settings);
478 //! First routine called by the analysis framework
479 // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
480 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
481 const TopologyInformation &top);
483 virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings &settings,
484 const t_trxframe &fr);
486 //! Call for each frame of the trajectory
487 // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
488 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
489 TrajectoryAnalysisModuleData *pdata);
491 //! Last routine called by the analysis framework
492 // virtual void finishAnalysis(t_pbc *pbc);
493 virtual void finishAnalysis(int nframes);
495 //! Routine to write output, that is additional over the built-in
496 virtual void writeOutput();
502 std::vector< std::vector< RVec > > trajectory;
503 std::vector< RVec > reference;
505 std::vector< int > index;
506 std::vector< std::string > resNms;
509 int tau = 0; // selectable
511 float crl_border = 0; // selectable
512 float eff_rad = 1.5; // selectable
513 std::string OutPutName; // selectable
514 int mode = 0; // selectable
515 // Copy and assign disallowed by base.
518 SpaceTimeCorr::SpaceTimeCorr(): TrajectoryAnalysisModule()
522 SpaceTimeCorr::~SpaceTimeCorr()
527 SpaceTimeCorr::initOptions(IOptionsContainer *options,
528 TrajectoryAnalysisSettings *settings)
530 static const char *const desc[] = {
531 "[THISMODULE] to be done"
533 // Add the descriptive text (program help text) to the options
534 settings->setHelpText(desc);
535 // Add option for output file name
536 //options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
537 // .store(&fnNdx_).defaultBasename("domains")
538 // .description("Index file from the domains"));
539 // Add option for working mode
540 options->addOption(gmx::IntegerOption("mode")
542 .description("default 0 - creation of sigle matrixes' file | 1 - reading prepared correlation file to make graphs and etc."));
543 // Add option for tau constant
544 options->addOption(gmx::IntegerOption("tau")
546 .description("max number of frames to shift to evaluate correlations | all starts from array [0 ; tau]"));
547 // Add option for window constant
548 options->addOption(gmx::IntegerOption("window")
550 .description("size of window in frames that correlation evaluation based on"));
551 // Add option for crl_border constant
552 options->addOption(FloatOption("crl")
554 .description("make graphs based on corrs > crl const"));
555 // Add option for eff_rad constant
556 options->addOption(FloatOption("ef_rad")
558 .description("effective radius for atoms to see each other to evaluate corrs"));
559 // Add option for selection list
560 options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
561 .required().dynamicMask().multiValue()
562 .description("Domains to form rigid skeleton"));
563 // Add option for output file names
564 options->addOption(StringOption("out_put")
566 .description("<your name here> + <local file tag>.txt"));
567 // Control input settings
568 settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
569 settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
570 settings->setPBC(true);
574 SpaceTimeCorr::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation &top)
576 ArrayRef< const int > atomind = sel_[0].atomIndices();
578 for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
579 index.push_back(*ai);
580 resNms.push_back(*(top.atoms()->resinfo[top.atoms()->atom[*ai].resind].name) + std::to_string((top.atoms()->atom[*ai].resind + 1)));
582 trajectory.resize(0);
585 if (top.hasFullTopology()) {
586 for (unsigned int i = 0; i < index.size(); i++) {
587 reference.push_back(top.x().at(index[i]));
593 SpaceTimeCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, const t_trxframe &fr)
597 // -s *.tpr -f *.xtc -n *.ndx -sf 'name' -mode <> -tau <> -window <> -crl <> -ef_rad <> -out_put <>
600 SpaceTimeCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
602 trajectory.resize(trajectory.size() + 1);
603 trajectory.back().resize(index.size());
604 for (unsigned int i = 0; i < index.size(); i++) {
605 trajectory.back()[i] = fr.x[index[i]];
611 SpaceTimeCorr::finishAnalysis(int nframes)
613 std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > rout_new;
614 unsigned int k = 1000, m = 0;
616 k = static_cast< unsigned int>(tau);
623 std::cout << "\nCorrelation's evaluation - start\n";
625 correlation_evaluation(reference, trajectory, k, static_cast< unsigned int >(window), (OutPutName + "-DumpData.txt").c_str());
627 std::cout << "Correlation's evaluation - end\n";
628 } else if (mode == 1) {
629 rout_new.resize(trajectory.size() - static_cast< unsigned long>(window) - k);
631 file = std::fopen((OutPutName + "-DumpData.txt").c_str(), "r+");
632 #pragma omp parallel for ordered schedule(dynamic) shared(rout_new) firstprivate(reference, crl_border, eff_rad, m, k)
633 for(unsigned long i = 0; i < trajectory.size() - static_cast< unsigned int>(window) - k; i++) {
634 std::vector< std::vector< std::vector< double > > > crltns;
635 std::vector< std::vector< std::pair< double, long int > > > graph;
636 std::vector< std::vector< unsigned int > > sub_graph;
637 std::vector< std::vector< std::pair< unsigned int, unsigned int > > > sub_graph_rbr;
641 crltns.resize(k + 1);
642 for (unsigned int i1 = 0; i1 < crltns.size(); i1++) {
643 crltns[i1].resize(trajectory.front().size());
644 for (unsigned int i2 = 0; i2 < trajectory.front().size(); i2++) {
645 crltns[i1][i2].resize(trajectory.front().size());
646 for (unsigned int i3 = 0; i3 < trajectory.front().size(); i3++) {
647 crltns[i1][i2][i3] = 0.;
651 for (unsigned int i1 = 0; i1 < k + 1; i1++) {
652 int t0, t1, t2 = std::fscanf(file, "%d %d\n", &t0, &t1);
653 for (unsigned int i2 = 0; i2 < trajectory.front().size(); i2++) {
654 for (unsigned int i3 = 0; i3 < trajectory.front().size(); i3++) {
655 t2 = std::fscanf(file, "%lf ", &crltns[i1][i2][i3]);
659 if (i % 10 == 0 && i > 0) {
662 std::cout << " | mtrxs " << i << " red | ";
665 graph_calculation(graph, sub_graph, sub_graph_rbr, reference, crltns, static_cast< double >(crl_border), static_cast< double >(eff_rad), m, k);
666 graph_back_bone_evaluation(rout_new[i], reference.size()/*index.size()*/, graph, sub_graph, sub_graph_rbr);
670 std::cout << "\n Almost - end\n";
671 make_rout_file(static_cast< double >(crl_border), index, resNms, rout_new, (OutPutName + "-routs.txt").c_str(), static_cast< unsigned int >(tau), static_cast< unsigned int >(window));
673 std::cout << "Finish Analysis - end\n\n";
677 SpaceTimeCorr::writeOutput()
683 * The main function for the analysis template.
686 main(int argc, char *argv[])
688 return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SpaceTimeCorr>(argc, argv);