From 5ee6328ff4740e7e4af653ac159fde5a30453a66 Mon Sep 17 00:00:00 2001 From: Anatoly Titov Date: Wed, 10 May 2017 17:52:34 +0300 Subject: [PATCH] fixes --- src/spacetimecorr.cpp | 132 ++++++++++++++++++++++++++++-------------- 1 file changed, 87 insertions(+), 45 deletions(-) diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp index ce95293..6e5d624 100644 --- a/src/spacetimecorr.cpp +++ b/src/spacetimecorr.cpp @@ -77,14 +77,12 @@ using gmx::RVec; void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const char* file_name) { - FILE *file; + FILE *file; file = std::fopen(file_name, "w+"); for (int i = 1; i < trajectory.size(); i++) { std::fprintf(file, "MODEL%9d\n", i); for (int j = 0; j < trajectory[i].size(); j++) { std::fprintf(file, "ATOM %5d CA LYS 1 %8.3f%8.3f%8.3f 1.00 0.00\n", j, trajectory[i][j][0] * 10, trajectory[i][j][1] * 10, trajectory[i][j][2] * 10); - // 123456 123456789012345678901 8 6 4567890123456 - //std::fprintf(file, "1234567890123456789012345678901234567890123456789012345678901234567890\n"); } std::fprintf(file, "ENDMDL\n"); } @@ -108,10 +106,9 @@ void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const ch 79 - 80 LString(2) Charge on the atom. */ -void make_correlation_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start) +void make_correlation_file(std::vector< std::vector< std::vector< double > > > correlations, const char* file_name, int start) { - - FILE *file; + FILE *file; file = std::fopen(file_name, "w+"); std::cout << "start printing\n"; for (int i = start; i < correlations.size(); i++) { @@ -132,8 +129,7 @@ void make_correlation_file(std::vector< std::vector< std::vector< float > > > co void make_rout_file(std::vector< int > dist, std::vector< std::vector< int > > rout, const char* file_name) { - - FILE *file; + FILE *file; file = std::fopen(file_name, "w+"); std::cout << "start printing\n"; for (int i = 0; i < rout.size(); i++) { @@ -419,39 +415,65 @@ void Domains::finishAnalysis(int nframes) { //make_pdb_trajectory(frankenstein_trajectory, "/home/toluk/Develop/samples/reca_rd/frank_with_duo_test.pdb"); - std::vector< std::vector< std::vector< float > > > crltns; - int k = 1000, m = 1; - crltns.resize(k); + std::vector< std::vector< std::vector< double > > > crltns; + int k = 950, m = 949; + crltns.resize(k + 1); for (int i = 0; i < crltns.size(); i++) { crltns[i].resize(index.size()); for (int j = 0; j < crltns[i].size(); j++) { crltns[i][j].resize(index.size(), 0); } } - #pragma omp parallel - { - #pragma omp for schedule(dynamic) - for (int i = m; i < k; i += 1) { + + //#pragma omp parallel + //{ + //#pragma omp for schedule(dynamic) + for (int i = m; i <= k; i += 1) { std::cout << " \n " << " k = " << i << " \n"; - std::vector< float > medx, medy; - medx.resize(index.size(), 0); - medy.resize(index.size(), 0); + //std::vector< float > medx, medy; + //medx.resize(index.size(), 0); + //medy.resize(index.size(), 0); + rvec *medx, *medy, temp_zero; + snew(medx, index.size()); + snew(medy, index.size()); + temp_zero[0] = 0; + temp_zero[1] = 0; + temp_zero[2] = 0; + for (int i = 0; i < index.size(); i++) { + copy_rvec(temp_zero, medx[i]); + copy_rvec(temp_zero, medy[i]); + } for (int j = 0; j < frankenstein_trajectory.size() - k; j++) { for (int f = 0; f < index.size(); f++) { rvec temp; - rvec_sub(frankenstein_trajectory[j][f], frankenstein_trajectory[j + 1][f], temp); - //rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + 1][f], temp); - medx[f] += norm(temp); - rvec_sub(frankenstein_trajectory[j + k - 1][f], frankenstein_trajectory[j + k][f], temp); - //rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + k][f], temp); - medy[f] += norm(temp); + //rvec_sub(frankenstein_trajectory[j][f], frankenstein_trajectory[j + 1][f], temp); + rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j][f], temp); + rvec_inc(medx[f], temp); + //medx[f] += norm(temp); + //rvec_sub(frankenstein_trajectory[j + k - 1][f], frankenstein_trajectory[j + k][f], temp); + rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + k][f], temp); + rvec_inc(medy[f], temp); + //medy[f] += norm(temp); } } + for (int j = 0; j < index.size(); j++) { - medx[j] /= frankenstein_trajectory.size(); - medy[j] /= frankenstein_trajectory.size(); + //medx[j] /= frankenstein_trajectory.size(); + rvec temp; + real temp2 = frankenstein_trajectory.size() - k; + temp2 = 1 / temp2; + copy_rvec(medx[j], temp); + svmul(temp2, temp, medx[j]); + //medy[j] /= frankenstein_trajectory.size(); + copy_rvec(medy[j], temp); + svmul(temp2, temp, medy[j]); } - std::vector< std::vector< float > > a, b, c; + + /*for (int j = 0; j < index.size(); j++) { + std::cout << medx[j][0] << " " << medx[j][1] << " " << medx[j][2] << " " << medy[j][0] << " " << medy[j][1] << " " << medy[j][2] << "\n"; + }*/ + + std::vector< std::vector< double > > a, b, c; a.resize(index.size()); b.resize(index.size()); c.resize(index.size()); @@ -460,32 +482,52 @@ Domains::finishAnalysis(int nframes) b[j].resize(index.size(), 0); c[j].resize(index.size(), 0); } - for (int j = 1; j < frankenstein_trajectory.size() - k; j++) { + for (int j = 0; j < frankenstein_trajectory.size() - k; j++) { for (int f1 = 0; f1 < index.size(); f1++) { - for (int f2 = f1; f2 < index.size(); f2++) { + for (int f2 = 0; f2 < index.size(); f2++) { rvec temp1, temp2; - rvec_sub(frankenstein_trajectory[j][f1], frankenstein_trajectory[j + 1][f1], temp1); - //rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j + 1][f1], temp1); - rvec_sub(frankenstein_trajectory[j + k - 1][f2], frankenstein_trajectory[j + k][f2], temp2); - //rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + k][f2], temp2); - a[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp2) - medy[f2]); - b[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp1) - medx[f1]); - c[f1][f2] += (norm(temp2) - medy[f2]) * (norm(temp2) - medy[f2]); + //rvec_sub(frankenstein_trajectory[j][f1], frankenstein_trajectory[j + 1][f1], temp1); + rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j][f1], temp1); + rvec_dec(temp1, medx[f1]); + //rvec_sub(frankenstein_trajectory[j + k - 1][f2], frankenstein_trajectory[j + k][f2], temp2); + rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + k][f2], temp2); + rvec_dec(temp2, medy[f2]); + //a[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp2) - medy[f2]); + //b[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp1) - medx[f1]); + //c[f1][f2] += (norm(temp2) - medy[f2]) * (norm(temp2) - medy[f2]); + a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]); + b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]); + c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]); + //std::cout << a[f1][f2] << " " << b[f1][f2] << " " << c[f1][f2] << "\n"; + } } } for (int j = 0; j < index.size(); j++) { - for (int f = j; f < index.size(); f++) { - crltns[i][j][f] = a[j][f] / (std::sqrt(b[j][f]) * std::sqrt(c[j][f])); + for (int f = 0; f < index.size(); f++) { + crltns[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f])); } } + sfree(medx); + sfree(medy); } - } + //} // std::cout << "\noutput file\n"; -// make_correlation_file(crltns, "/home/toluk/Develop/samples/reca_rd/corrs2.txt", m); + make_correlation_file(crltns, "/home/toluk/Develop/samples/reca_rd/corrs2.txt", m); // std::cout << "\n\n file done \n\n"; - std::cout << "\n\n\n\n"; + + + + + + + + + + + + /*std::cout << "\n\n\n\n"; for (int i2 = 0; i2 < 75; i2++) { long int count_here = 0; for (int i = 0; i < k; i++) { @@ -499,7 +541,7 @@ Domains::finishAnalysis(int nframes) } std::cout << i2 << " " << count_here << "\n"; } - std::cout << "\n\n\n\n"; + std::cout << "\n\n\n\n";*/ @@ -535,11 +577,11 @@ Domains::finishAnalysis(int nframes) std::cout << "building routs\n" ; - +/* for (int i = 0; i < index.size(); i++) { - std::cout << dive_in(graph, i, 0, frames, 0); + std::cout << i << " " << dive_in(graph, i, 0, frames, 0) << " "; } - +*/ /*while (flag) { flag = false; for (long int i = 0; i < rout_old.size(); i++) { -- 2.22.0