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