src/spacetimecorr.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements gmx::analysismodules::Freevolume.
  38  *
  39  * \author Titov Anatoly <Wapuk-cobaka@yandex.ru>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42
  43
  44 #include <iostream>
  45 #include <algorithm>
  46 #include <vector>
  47 #include <math.h>
  48 #include <omp.h>
  49 #include <string>
  50
  51 #include <gromacs/trajectoryanalysis.h>
  52 #include <gromacs/utility/smalloc.h>
  53 #include <gromacs/math/do_fit.h>
  54
  55 using namespace gmx;
  56
  57 using gmx::RVec;
  58
  59 void make_correlation_matrix_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start)
  60 {
  61     FILE *file;
  62     file = std::fopen(file_name, "w+");
  63     for (int i = start; i < correlations.size(); i++) {
  64         if (correlations.size() - start > 1) {
  65             std::fprintf(file, "correlation between 'n' and 'n + %d' frames\n", i);
  66         }
  67         if (i % 50 == 0) {
  68             std::cout << "correlation between 'n' and 'n + " << i << "' frames\n";
  69         }
  70         for (int j = 0; j < correlations[i].size(); j++) {
  71             for (int f = 0; f < correlations[i][j].size(); f++) {
  72                 std::fprintf(file, "%3.4f ", correlations[i][j][f]);
  73             }
  74             std::fprintf(file, "\n");
  75         }
  76     }
  77     std::fclose(file);
  78 }
  79
  80 void make_correlation_pairs_file(std::vector< std::vector< std::vector< float > > > correlations, const char* file_name, int start)
  81 {
  82     FILE *file;
  83     file = std::fopen(file_name, "w+");
  84     for (int i = 0; i < correlations.front().size(); i++) {
  85         for (int j = 0; j < correlations.front().front().size(); j++) {
  86             //std::fprintf(file, "correlation between point '%d' and point '%d'\n", i, j);
  87             std::fprintf(file, "%d %d\n", i, j);
  88             for (int k = 0; k < correlations.size(); k++) {
  89                 std::fprintf(file, "%3.4f ", correlations[k][i][j]);
  90             }
  91             std::fprintf(file, "\n");
  92         }
  93         if (i % 10 == 0) {
  94             std::cout << "correlations in row " << i << " out of " << correlations.front().size() << " completed\n";
  95         }
  96     }
  97     std::fclose(file);
  98 }
  99
 100 void correlation_evaluation(std::vector< std::vector< RVec > > traj, int window_size, std::vector< std::vector< std::vector< float > > > &crl) {
 101     crl.resize(traj.size() - window_size + 1);
 102     for (int i = 0; i < crl.size(); i++) {
 103         crl[i].resize(traj.front().size());
 104         for (int j = 0; j < crl[i].size(); j++) {
 105             crl[i][j].resize(traj.front().size(), 0);
 106         }
 107     }
 108
 109     RVec temp_zero;
 110     temp_zero[0] = 0;
 111     temp_zero[1] = 0;
 112     temp_zero[2] = 0;
 113
 114     std::vector< float > d;
 115     d.resize(traj.front().size(), 0);
 116     //int nth = omp_get_thread_num();
 117
 118
 119     //changed parrallel
 120     //#pragma omp parallel shared(crl, temp_zero, d)
 121     //{
 122     #pragma omp parallel for schedule(dynamic)
 123     for (int i = 0; i < crl.size(); i++) {
 124         std::vector< RVec > medx, medy;
 125         RVec temp1, temp2;
 126         float tmp;
 127         medx.resize(traj.front().size(), temp_zero);
 128         medy.resize(traj.front().size(), temp_zero);
 129             //#pragma omp for schedule(dynamic)
 130         for (int j = i; j < i + window_size - 1; j++) {
 131             for (int f = 0; f < traj.front().size(); f++) {
 132                 medx[f] += (traj[i][f] - traj[j][f]);
 133                 medy[f] += (traj[i][f] - traj[j][f]);
 134             }
 135         }
 136             //#pragma omp barrier
 137             //#pragma omp for schedule(dynamic)
 138         for (int j = 0; j < traj.front().size(); j++) {
 139             tmp = traj.size() - 1;
 140             tmp -= i;
 141             tmp = 1 / tmp;
 142                 //tmp = 1 / (traj.size() - 1 - i);
 143             medx[j] *= tmp;
 144             medy[j] *= tmp;
 145         }
 146             //#pragma omp barrier
 147         std::vector< std::vector< float > > a, b, c;
 148         a.resize(traj.front().size(), d);
 149         b.resize(traj.front().size(), d);
 150         c.resize(traj.front().size(), d);
 151         for (int j = i; j < i + window_size - 1; j++) {
 152             for (int f1 = 0; f1 < traj.front().size(); f1++) {
 153                     //#pragma omp for schedule(dynamic)
 154                 for (int f2 = 0; f2 < traj.front().size(); f2++) {
 155                     temp1 = traj[i][f1] - traj[j][f1] - medx[f1];
 156                     temp2 = traj[i][f2] - traj[j][f2] - medy[f2];
 157                     a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
 158                     b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
 159                     c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
 160                 }
 161                     //#pragma omp barrier
 162             }
 163         }
 164
 165             //#pragma omp for schedule(dynamic)
 166         for (int j = 0; j < traj.front().size(); j++) {
 167             for (int f = 0; f < traj.front().size(); f++) {
 168                 crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
 169             }
 170         }
 171             //#pragma omp barrier
 172         medx.resize(0);
 173         medy.resize(0);
 174         std::cout << i << " corr done\n";
 175     }
 176     #pragma omp barrier
 177     //}
 178 }
 179
 180 /*! \brief
 181  * \ingroup module_trajectoryanalysis
 182  */
 183
 184 class SimpleCorr : public TrajectoryAnalysisModule
 185 {
 186     public:
 187
 188         SimpleCorr();
 189         virtual ~SimpleCorr();
 190
 191         //! Set the options and setting
 192         virtual void initOptions(IOptionsContainer          *options,
 193                                  TrajectoryAnalysisSettings *settings);
 194
 195         //! First routine called by the analysis framework
 196         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 197         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 198                                   const TopologyInformation        &top);
 199
 200         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 201                                          const t_trxframe                   &fr);
 202
 203         //! Call for each frame of the trajectory
 204         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 205         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 206                                   TrajectoryAnalysisModuleData *pdata);
 207
 208         //! Last routine called by the analysis framework
 209         // virtual void finishAnalysis(t_pbc *pbc);
 210         virtual void finishAnalysis(int nframes);
 211
 212         //! Routine to write output, that is additional over the built-in
 213         virtual void writeOutput();
 214
 215     private:
 216
 217         SelectionList                                               sel_;
 218
 219         std::vector< std::vector< RVec > >                          trajectory;
 220
 221         std::vector< int >                                          index;
 222         int                                                         frames              = 0;
 223         int                                                         W_size              = 0;
 224         std::string                                                 OutPutName; // selectable
 225         // Copy and assign disallowed by base.
 226 };
 227
 228 SimpleCorr::SimpleCorr(): TrajectoryAnalysisModule()
 229 {
 230 }
 231
 232 SimpleCorr::~SimpleCorr()
 233 {
 234 }
 235
 236 void
 237 SimpleCorr::initOptions(IOptionsContainer          *options,
 238                      TrajectoryAnalysisSettings *settings)
 239 {
 240     static const char *const desc[] = {
 241         "[THISMODULE] to be done"
 242     };
 243     // Add the descriptive text (program help text) to the options
 244     settings->setHelpText(desc);
 245     // Add option for output file name
 246     //options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 247     //                        .store(&fnNdx_).defaultBasename("domains")
 248     //                        .description("Index file from the domains"));
 249     // Add option for time window size constant
 250     options->addOption(gmx::IntegerOption("Window")
 251                             .store(&W_size)
 252                             .description("number of frames for time travel"));
 253     // Add option for selection list
 254     options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
 255                            .required().dynamicMask().multiValue()
 256                            .description("Domains to form rigid skeleton"));
 257     options->addOption(StringOption("out_put")
 258                             .store(&OutPutName)
 259                             .description("<your name here> + <local file tag>.txt"));
 260     // Control input settings
 261     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 262     settings->setPBC(true);
 263 }
 264
 265 void
 266 SimpleCorr::initAnalysis(const TrajectoryAnalysisSettings &settings,
 267                       const TopologyInformation        &top)
 268 {
 269     ArrayRef< const int > atomind = sel_[0].atomIndices();
 270     index.resize(0);
 271     for (ArrayRef< const int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 272         index.push_back(*ai);
 273     }
 274     trajectory.resize(0);
 275 }
 276
 277 void
 278 SimpleCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 279                              const t_trxframe                       &fr)
 280 {
 281 }
 282
 283 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.md_npt.no-sol.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/CA.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/SLca' -Window 1000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/window'
 284 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/trp-cage.md_npt.no-sol.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/CA.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/SLca' -Window 1000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/window'
 285
 286 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.md_npt.no-sol.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/CA.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/SLca' -Window 10000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/window10'
 287 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/trp-cage.md_npt.no-sol.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/CA.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/SLca' -Window 10000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/window10'
 288
 289 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/trp-cage.md_npt.no-sol.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/CA.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/SLca' -Window 100000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv1/window100'
 290 // -s '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/trp-cage.pdb' -f '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/trp-cage.md_npt.no-sol.xtc' -n '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/CA.ndx' -sf '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/SLca' -Window 100000 -out_put '/home/titov_ai/BioStore/titov_ai/trp_cage/fullTrj/kv2/window100'
 291
 292 void
 293 SimpleCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 294                       TrajectoryAnalysisModuleData *pdata)
 295 {
 296     trajectory.resize(trajectory.size() + 1);
 297     trajectory.back().resize(index.size());
 298     for (int i = 0; i < index.size(); i++) {
 299         trajectory.back()[i] = fr.x[index[i]];
 300     }
 301     frames++;
 302 }
 303
 304 void
 305 SimpleCorr::finishAnalysis(int nframes)
 306 {
 307     std::vector< std::vector< std::vector< float > > > crltns;
 308
 309     std::cout << "\nCorrelation's evaluation - start\n";
 310     correlation_evaluation(trajectory, W_size, crltns);
 311     std::cout << "\nCorrelation's evaluation - end\n";
 312
 313     std::cout << "corelation matrix - start printing\n";
 314     make_correlation_matrix_file(crltns, (OutPutName + "_matrix.txt").c_str(), 0);
 315     std::cout << "corelation matrix - printed\n";
 316
 317     std::cout << "corelation pairs - start printing\n";
 318     make_correlation_pairs_file(crltns, (OutPutName + "_pairs.txt").c_str(), 0);
 319     std::cout << "corelation pairs - printed\n";
 320
 321 }
 322
 323 void
 324 SimpleCorr::writeOutput()
 325 {
 326     std::cout << "GAME OVER\n";
 327 }
 328
 329 /*! \brief
 330  * The main function for the analysis template.
 331  */
 332 int
 333 main(int argc, char *argv[])
 334 {
 335     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<SimpleCorr>(argc, argv);
 336 }