commit for diploma

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index 47cfacb913d3e8a63a690152d39d06237e54b716..c47b4d62682e9d8130c710c72f0526c29909f87a 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -88,24 +88,6 @@ void make_pdb_trajectory(std::vector< std::vector< RVec > > trajectory, const ch
         std::fprintf(file, "ENDMDL\n");
     }
 }
-/*
-1 -  6        Record name     "ATOM  "
-7 - 11        Integer         Atom serial number.
-13 - 16        Atom            Atom name.
-17             Character       Alternate location indicator.
-18 - 20        Residue name    Residue name.
-22             Character       Chain identifier.
-23 - 26        Integer         Residue sequence number.
-27             AChar           Code for insertion of residues.
-31 - 38        Real(8.3)       Orthogonal coordinates for X in Angstroms.
-39 - 46        Real(8.3)       Orthogonal coordinates for Y in Angstroms.
-47 - 54        Real(8.3)       Orthogonal coordinates for Z in Angstroms.
-55 - 60        Real(6.2)       Occupancy.
-61 - 66        Real(6.2)       Temperature factor (Default = 0.0).
-73 - 76        LString(4)      Segment identifier, left-justified.
-77 - 78        LString(2)      Element symbol, right-justified.
-79 - 80        LString(2)      Charge on the atom.
-*/
 
 void make_correlation_file(std::vector< std::vector< std::vector< double > > > correlations, const char* file_name, int start)
 {
@@ -142,19 +124,12 @@ void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std
     std::fclose(file);
 }
 
-//template < typename T >
 bool mysortfunc (std::vector< int > a, std::vector< int > b) { return (a.size() > b.size()); }
 
 /*! \brief
- * Class used to compute free volume in a simulations box.
- *
- * Inherits TrajectoryAnalysisModule and all functions from there.
- * Does not implement any new functionality.
- *
  * \ingroup module_trajectoryanalysis
  */
 
-
 class Domains : public TrajectoryAnalysisModule
 {
     public:
@@ -196,7 +171,6 @@ class Domains : public TrajectoryAnalysisModule
         std::vector< std::vector< RVec > >                          trajectory;
         std::vector< std::vector< RVec > >                          frankenstein_trajectory;
 
-
         std::vector< int >                                          index;
         std::vector< int >                                          numbers;
         int                                                         frames              = 0;
@@ -402,31 +376,21 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
     frames++;
 }
 
-//spacetimecorr -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test.ndx'
-
 void
 Domains::finishAnalysis(int nframes)
 {
-    //make_pdb_trajectory(frankenstein_trajectory, "/home/toluk/Develop/samples/reca_rd/10k_frames.pdb");
-
-    //frankenstein_trajectory = trajectory;
+    make_pdb_trajectory(frankenstein_trajectory, "/home/toluk/Develop/samples/reca_rd/10k_frames.pdb");
 
-    std::vector< std::vector< std::vector< double > > > crltns, crltns_scalars, crltns_mini;
+    std::vector< std::vector< std::vector< double > > > crltns;
     int k = 1000, m = 0;
     if (tau > 0) {
         k = tau;
     }
     crltns.resize(k + 1);
-    crltns_scalars.resize(k + 1);
-    crltns_mini.resize(k + 1);
     for (int i = 0; i < crltns.size(); i++) {
         crltns[i].resize(index.size());
-        crltns_scalars[i].resize(index.size());
-        crltns_mini[i].resize(index.size());
         for (int j = 0; j < crltns[i].size(); j++) {
             crltns[i][j].resize(index.size(), 0);
-            crltns_scalars[i][j].resize(index.size(), 0);
-            crltns_mini[i][j].resize(index.size(), 0);
         }
     }
 
@@ -437,22 +401,15 @@ Domains::finishAnalysis(int nframes)
         for (int i = m; i <= k; i += 1) {
             std::cout << " k = " << i << " ";
 
-            rvec *medx, *medy, temp_zero, *medx_mini, *medy_mini;
-            std::vector< double > medx_scalar, medy_scalar;
+            rvec *medx, *medy, temp_zero;
             snew(medx, index.size());
             snew(medy, index.size());
-            medx_scalar.resize(index.size(), 0);
-            medy_scalar.resize(index.size(), 0);
-            snew(medx_mini, index.size());
-            snew(medy_mini, 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]);
-                copy_rvec(temp_zero, medx_mini[i]);
-                copy_rvec(temp_zero, medy_mini[i]);
             }
 
             for (int j = 0; j < frankenstein_trajectory.size() - i - 1; j++) {
@@ -462,18 +419,8 @@ Domains::finishAnalysis(int nframes)
                     rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j][f], temp);
                     rvec_inc(medx[f], temp);
 
-                    medx_scalar[f] += norm(temp);
-
                     rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + i][f], temp);
                     rvec_inc(medy[f], temp);
-
-                    medy_scalar[f] += norm(temp);
-
-                    rvec_sub(frankenstein_trajectory[j][f], frankenstein_trajectory[j + 1][f], temp);
-                    rvec_inc(medx_mini[f], temp);
-                    rvec_sub(frankenstein_trajectory[j + i][f], frankenstein_trajectory[j + i + 1][f], temp);
-                    rvec_inc(medy_mini[f], temp);
-
                 }
             }
 
@@ -487,37 +434,16 @@ Domains::finishAnalysis(int nframes)
                 svmul(temp2, temp, medx[j]);
                 copy_rvec(medy[j], temp);
                 svmul(temp2, temp, medy[j]);
-
-                medx_scalar[j] *= temp2;
-                medy_scalar[j] *= temp2;
-
-                copy_rvec(medx_mini[j], temp);
-                svmul(temp2, temp, medx_mini[j]);
-                copy_rvec(medy_mini[j], temp);
-                svmul(temp2, temp, medy_mini[j]);
-
             }
 
-            std::vector< std::vector< double > > a, b, c, as, bs, cs, am, bm, cm;
+            std::vector< std::vector< double > > a, b, c;
             a.resize(index.size());
             b.resize(index.size());
             c.resize(index.size());
-            as.resize(index.size());
-            bs.resize(index.size());
-            cs.resize(index.size());
-            am.resize(index.size());
-            bm.resize(index.size());
-            cm.resize(index.size());
             for (int j = 0; j < index.size(); j++) {
                 a[j].resize(index.size(), 0);
                 b[j].resize(index.size(), 0);
                 c[j].resize(index.size(), 0);
-                as[j].resize(index.size(), 0);
-                bs[j].resize(index.size(), 0);
-                cs[j].resize(index.size(), 0);
-                am[j].resize(index.size(), 0);
-                bm[j].resize(index.size(), 0);
-                cm[j].resize(index.size(), 0);
             }
             for (int j = 0; j < frankenstein_trajectory.size() - i - 1; j++) {
                 for (int f1 = 0; f1 < index.size(); f1++) {
@@ -525,47 +451,29 @@ Domains::finishAnalysis(int nframes)
                         rvec temp1, temp2;
                         rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j][f1], temp1);
                         rvec_dec(temp1, medx[f1]);
+
                         rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + i][f2], temp2);
                         rvec_dec(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]);
-
-                        rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j][f1], temp1);
-                        rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + i][f2], temp2);
-                        as[f1][f2] += (norm(temp1) - medx_scalar[f1]) * (norm(temp2) - medy_scalar[f2]);
-                        bs[f1][f2] += (norm(temp1) - medx_scalar[f1]) * (norm(temp1) - medx_scalar[f1]);
-                        cs[f1][f2] += (norm(temp2) - medy_scalar[f2]) * (norm(temp2) - medy_scalar[f2]);
-
-                        rvec_sub(frankenstein_trajectory[j][f1], frankenstein_trajectory[j + 1][f1], temp1);
-                        rvec_dec(temp1, medx_mini[f1]);
-                        rvec_sub(frankenstein_trajectory[j + i][f2], frankenstein_trajectory[j + i + 1][f2], temp2);
-                        rvec_dec(temp2, medy_mini[f2]);
-                        am[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
-                        bm[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
-                        cm[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
                     }
                 }
             }
             for (int j = 0; j < index.size(); j++) {
                 for (int f = 0; f < index.size(); f++) {
                     crltns[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
-                    crltns_scalars[i][j][f] = as[j][f] / (std::sqrt(bs[j][f] * cs[j][f]));
-                    crltns_mini[i][j][f] = am[j][f] / (std::sqrt(bm[j][f] * cm[j][f]));
                 }
             }
             sfree(medx);
             sfree(medy);
-            sfree(medx_mini);
-            sfree(medy_mini);
         }
     }
     std::cout << "\n" ;
 
     std::cout << "printing correlation's file\n" ;
     make_correlation_file(crltns, "corrs_vectors.txt", m);
-    make_correlation_file(crltns_scalars, "corrs_vectors_scalars.txt", m);
-    make_correlation_file(crltns_mini, "corrs_mini_vectors.txt", m);
     std::cout << "done\n" ;
 
     for (int i1 = 0; i1 < 100; i1++) {
@@ -579,17 +487,15 @@ Domains::finishAnalysis(int nframes)
                 }
             }
         }
-        std::cout << i1 << " - " << i5 << "   ";
+        std::cout << i1 << " - " << i5 << " | ";
         if (i1 % 10 == 0) {
             std::cout << "\n" ;
         }
     }
+    std::cout << "\n";
 
-    std::cout << "\n\n";
-
-    std::vector< std::vector< int > > graph, graph_tau;
+    std::vector< std::vector< int > > graph;
     graph.resize(index.size());
-    graph_tau.resize(index.size());
 
     for (int i = 1; i <= k; i++) {
         for (int j = 0; j < index.size(); j++) {
@@ -603,19 +509,12 @@ Domains::finishAnalysis(int nframes)
                     }
                     if (local_flag) {
                         graph[j].push_back(f);
-                        graph_tau[j].push_back(i);
                     }
                 }
             }
         }
     }
 
-    /*for (int i = 0; i < index.size(); i++) {
-        if (graph[i].size() > 0) {
-            std::cout << i << " - " << graph[i].size() << "\n";
-        }
-    }*/
-
     std::cout << "\n";
     std::vector< std::vector < int > > rout_old, rout_new, rout_out;
     bool flag = true;
@@ -677,26 +576,18 @@ Domains::finishAnalysis(int nframes)
         std::sort(rout_out[i].begin(), rout_out[i].end());
     }
 
-    std::cout << "\nfinished routs\n\n" ;
-
     rout_new.resize(0);
     for (int i = rout_old.size() - 1; i >= 0; i--) {
         bool lflag1 = true;
         for (int j = 0; j < i; j++) {
             if (rout_old[i].size() <= rout_old[j].size()) {
-                int la = 0/*, lb1 = 0, lb2 = 0*/;
+                int la = 0;
                 for (int f = 0; f < rout_old[j].size(); f++) {
                     if (rout_out[j][f] == rout_out[i][la]) {
                         la++;
-                        /*if (la < rout_old[i].size()) {
-                            lb1 += rout_old[i][la];
-                        }*/
                     }
-                    /*if (la > 0 && la < rout_old[i].size() && f + 1 < rout_old[j].size()) {
-                        lb2 += rout_old[j][f + 1];
-                    }*/
                 }
-                if (la == rout_out[i].size() /*&& lb1 >= lb2*/) {
+                if (la == rout_out[i].size()) {
                     lflag1 = false;
                 }
             }
@@ -705,6 +596,7 @@ Domains::finishAnalysis(int nframes)
             rout_new.push_back(rout_old[i]);
         }
     }
+    std::cout << "\nfinished routs\n\n" ;
 
     make_rout_file(crl_border, index, rout_new, "routs.txt");
 
@@ -717,7 +609,6 @@ Domains::writeOutput()
 
 }
 
-
 /*! \brief
  * The main function for the analysis template.
  */