src/domains.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 #include <iostream>
  44 #include <chrono>
  45 #include <omp.h>
  46 #include <thread>
  47 #include <string>
  48
  49 #include <gromacs/trajectoryanalysis.h>
  50 #include <gromacs/pbcutil/pbc.h>
  51 #include <gromacs/pbcutil/rmpbc.h>
  52 #include <gromacs/utility/smalloc.h>
  53 #include <gromacs/math/vectypes.h>
  54 #include <gromacs/math/vec.h>
  55 #include <gromacs/math/do_fit.h>
  56
  57 #define MAX_NTRICVEC 12
  58
  59 using namespace gmx;
  60
  61 using gmx::RVec;
  62
  63 struct node {
  64     int  n;
  65     RVec r;
  66     bool yep;
  67 };
  68
  69 void make_graph(int mgwi_natoms, rvec *mgwi_x, std::vector< std::vector< node > > &mgwi_graph)
  70 {
  71     mgwi_graph.resize(mgwi_natoms);
  72     for (int i = 0; i < mgwi_natoms; i++) {
  73         mgwi_graph[i].resize(mgwi_natoms);
  74     }
  75     for (int i = 0; i < mgwi_natoms; i++) {
  76         for (int j = 0; j < mgwi_natoms; j++) {
  77             rvec_sub(mgwi_x[i], mgwi_x[j], mgwi_graph[i][j].r);
  78             mgwi_graph[i][j].n = 0;
  79         }
  80     }
  81 }
  82
  83 void update_graph(std::vector< std::vector< node > > &ugwi_graph, rvec *ugwi_x, long double ugwi_epsi) {
  84     rvec ugwi_temp;
  85     int ugwi_for = ugwi_graph.size();
  86     for (int i = 0; i < ugwi_for; i++) {
  87         for (int j = i; j < ugwi_for; j++) {
  88             rvec_sub(ugwi_x[i], ugwi_x[j], ugwi_temp);
  89             rvec_dec(ugwi_temp, ugwi_graph[i][j].r.as_vec());
  90             if (norm(ugwi_temp) <= ugwi_epsi) {
  91                 if (i == j) {
  92                     ugwi_graph[i][j].n++;
  93                 }
  94                 else {
  95                     ugwi_graph[i][j].n++;
  96                     ugwi_graph[j][i].n++;
  97                 }
  98             }
  99         }
 100     }
 101 }
 102
 103 void check_domains(long double cd_delta, int cd_frames, std::vector< std::vector< std::vector< node > > > &cd_graph) {
 104     int cd_for1 = cd_graph.size(), cd_for2 = cd_graph[1].size();
 105     for (int k = 0; k < cd_for1; k++) {
 106         for (int i = 0; i < cd_for2; i++) {
 107             for (int j = 0; j < cd_for2; j++) {
 108                 if (cd_graph[k][i][j].n >= cd_frames * cd_delta) {
 109                     cd_graph[k][i][j].yep = true;
 110                 }
 111                 else {
 112                     cd_graph[k][i][j].yep = false;
 113                 }
 114             }
 115         }
 116     }
 117 }
 118
 119 void find_domain_sizes(std::vector< std::vector< std::vector< node > > > fds_graph, std::vector< std::vector< int > > &fds_domsizes) {
 120     fds_domsizes.resize(fds_graph.size());
 121     int fds_for1 = fds_graph.size(), fds_for2 = fds_graph[1].size();
 122     for (int i = 0; i < fds_for1; i++) {
 123         fds_domsizes[i].resize(fds_for2, 0);
 124         for (int j = 0; j < fds_for2; j++) {
 125             for (int k = 0; k < fds_for2; k++) {
 126                 if (fds_graph[i][j][k].yep) {
 127                     fds_domsizes[i][j]++;
 128                 }
 129             }
 130         }
 131     }
 132 }
 133
 134 void get_maxsized_domain(std::vector< int > &gmd_max_d, std::vector< std::vector< std::vector< node > > > gmd_graph, std::vector< std::vector< int > > gmd_domsizes) {
 135     int gmd_number1 = 0, gmd_number2 = 0;
 136     int gmd_for1 = gmd_domsizes.size(), gmd_for2 = gmd_domsizes[0].size(), gmd_for3 = gmd_graph[1][1].size();
 137     for (int i = 0; i < gmd_for1; i++) {
 138         for (int j = 0; j < gmd_for2; j++) {
 139             if (gmd_domsizes[i][j] >= gmd_domsizes[gmd_number1][gmd_number2]) {
 140                 gmd_number1 = i;
 141                 gmd_number2 = j;
 142             }
 143         }
 144     }
 145     gmd_max_d.resize(0);
 146     for (int i = 0; i < gmd_for3; i++) {
 147         if (gmd_graph[gmd_number1][gmd_number2][i].yep) {
 148             gmd_max_d.push_back(i);
 149         }
 150     }
 151 }
 152
 153 void delete_domain_from_graph(std::vector< std::vector< std::vector< node > > > &ddf_graph, std::vector< int > ddf_domain) {
 154     int ddfg_for1 = ddf_domain.size(), ddfg_for2 = ddf_graph.size(), ddfg_for3 = ddf_graph[1].size();
 155     for (int i = 0; i < ddfg_for1; i++) {
 156         for (int k = 0; k < ddfg_for2; k++) {
 157             for (int j = 0; j < ddfg_for3; j++) {
 158                 if (ddf_graph[k][ddf_domain[i]][j].yep) {
 159                     ddf_graph[k][ddf_domain[i]][j].yep = false;
 160                 }
 161                 if (ddf_graph[k][j][ddf_domain[i]].yep) {
 162                     ddf_graph[k][j][ddf_domain[i]].yep = false;
 163                 }
 164             }
 165         }
 166     }
 167 }
 168
 169 bool check_domsizes(std::vector< std::vector< int > > cd_domsizes, int cd_domain_min_size) {
 170     int cd_for1 = cd_domsizes.size(), cd_for2 = cd_domsizes[0].size();
 171     for (int i = 0; i < cd_for1; i++) {
 172         for (int j = 0; j < cd_for2; j++) {
 173             if (cd_domsizes[i][j] >= cd_domain_min_size) {
 174                 return true;
 175             }
 176         }
 177     }
 178     return false;
 179 }
 180
 181 void print_domains(std::vector< std::vector< int > > pd_domains, std::vector< int > index, std::string fnNdx_) {
 182     FILE                *fpNdx_;
 183     fpNdx_ = std::fopen(fnNdx_.c_str(), "w+");
 184     int write_count;
 185     for (int i = 0; i < pd_domains.size(); i++) {
 186         std::fprintf(fpNdx_, "[domain_%d]\n", i + 1);
 187         write_count = 0;
 188         for (int j = 0; j < pd_domains[i].size(); j++) {
 189             write_count++;
 190             if (write_count > 20) {
 191                 write_count -= 20;
 192                 std::fprintf(fpNdx_, "\n");
 193             }
 194             std::fprintf(fpNdx_, "%5d ", index[pd_domains[i][j]] + 1);
 195         }
 196         std::fprintf(fpNdx_,"\n\n");
 197     }
 198     std::fprintf(fpNdx_,"\n");
 199     std::fclose(fpNdx_);
 200 }
 201
 202
 203 /*! \brief
 204  * Class used to compute free volume in a simulations box.
 205  *
 206  * Inherits TrajectoryAnalysisModule and all functions from there.
 207  * Does not implement any new functionality.
 208  *
 209  * \ingroup module_trajectoryanalysis
 210  */
 211 class Domains : public TrajectoryAnalysisModule
 212 {
 213     public:
 214
 215         Domains();
 216         virtual ~Domains();
 217
 218         //! Set the options and setting
 219         virtual void initOptions(IOptionsContainer          *options,
 220                                  TrajectoryAnalysisSettings *settings);
 221
 222         //! First routine called by the analysis framework
 223         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
 224         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
 225                                   const TopologyInformation        &top);
 226
 227         virtual void initAfterFirstFrame(const TrajectoryAnalysisSettings   &settings,
 228                                          const t_trxframe                   &fr);
 229
 230         //! Call for each frame of the trajectory
 231         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
 232         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 233                                   TrajectoryAnalysisModuleData *pdata);
 234
 235         //! Last routine called by the analysis framework
 236         // virtual void finishAnalysis(t_pbc *pbc);
 237         virtual void finishAnalysis(int nframes);
 238
 239         //! Routine to write output, that is additional over the built-in
 240         virtual void writeOutput();
 241
 242     private:
 243
 244         std::string                                                 fnNdx_;
 245         std::vector< std::vector< std::vector< node > > >           graph;
 246         std::vector< std::vector< int > >                           domains;
 247         std::vector< std::vector< int > >                           domsizes;
 248         std::vector< int >                                          index;
 249         std::vector< int >                                          numbers;
 250         std::vector< std::vector < RVec > >                         trajectory;
 251         Selection                                                   selec;
 252         int                                                         frames              = 0;
 253         int                                                         etalon_frame        = 0; // should be selectable
 254         int                                                         domain_min_size     = 5; // should be selectable
 255         double                                                      delta               = 0.90; //0.95 // should be selectable
 256         double                                                      epsi                = 0.15; //0.3 колебания внутри домена // should be selectable
 257         int                                                         domains_ePBC;
 258         // Copy and assign disallowed by base.
 259 };
 260
 261 Domains::Domains(): TrajectoryAnalysisModule()
 262 {
 263 }
 264
 265 Domains::~Domains()
 266 {
 267 }
 268
 269 void
 270 Domains::initOptions(IOptionsContainer          *options,
 271                      TrajectoryAnalysisSettings *settings)
 272 {
 273     static const char *const desc[] = {
 274         "[THISMODULE] to be done"
 275     };
 276     // Add the descriptive text (program help text) to the options
 277     settings->setHelpText(desc);
 278     // Add option for selecting a subset of atoms
 279     options->addOption(SelectionOption("select")
 280                            .store(&selec).required()
 281                            .description("Atoms that are considered as part of the excluded volume"));
 282     // Add option for output file name
 283     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 284                             .store(&fnNdx_).defaultBasename("domains")
 285                             .description("Index file from the domains"));
 286     // Add option for etalon_frame constant
 287     options->addOption(gmx::IntegerOption("etalon_frame")
 288                             .store(&etalon_frame)
 289                             .description("basic frame to base evaluation on"));
 290     // Add option for etalon_frame constant
 291     options->addOption(gmx::IntegerOption("dms")
 292                             .store(&domain_min_size)
 293                             .description("minimum domain size"));
 294     // Add option for epsi constant
 295     options->addOption(DoubleOption("epsilon")
 296                             .store(&epsi)
 297                             .description("thermal vibrations' constant"));
 298     // Add option for delta constant
 299     options->addOption(DoubleOption("delta")
 300                             .store(&delta)
 301                             .description("domain membership probability"));
 302     // Control input settings
 303     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 304     settings->setPBC(true);
 305 }
 306
 307 void
 308 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 309                       const TopologyInformation        &top)
 310 {
 311     domains_ePBC = top.ePBC();
 312 }
 313
 314 void
 315 Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
 316                              const t_trxframe                       &fr)
 317 {
 318     t_pbc  pbc;
 319     t_pbc *ppbc = settings.hasPBC() ? &pbc : NULL;
 320     matrix boxx;
 321     copy_mat(fr.box, boxx);
 322     if (ppbc != NULL) {
 323         set_pbc(ppbc, domains_ePBC, boxx);
 324     }
 325     ConstArrayRef< int > atomind  = selec.atomIndices();
 326     index.resize(0);
 327     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
 328         index.push_back(*ai);
 329     }
 330     trajectory.resize(frames + 1);
 331     trajectory[frames].resize(selec.atomCount());
 332
 333     for (int i = 0; i < selec.atomCount(); i++) {
 334         trajectory[frames][i] = fr.x[index[i]];
 335     }
 336     frames++;
 337     graph.resize(0);
 338 }
 339
 340 void
 341 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 342                       TrajectoryAnalysisModuleData *pdata)
 343 {
 344     trajectory.resize(frames + 1);
 345     trajectory[frames].resize(index.size());
 346     for (int i = 0; i < index.size(); i++) {
 347         trajectory[frames][i] = fr.x[index[i]];
 348     }
 349     frames++;
 350 }
 351
 352 //domains -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA'
 353 //domains -s '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.tpr' -f '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.10th.xtc' -select 'name CA'
 354 //domains -s '/home/toluk/Data/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Data/reca_rd/reca_rd.mono.xtc' - select 'name CA'
 355
 356 void
 357 Domains::finishAnalysis(int nframes)
 358 {
 359     frames--;
 360     int bone = index.size() - domain_min_size + 1;
 361     graph.resize(bone);
 362     real **w_rls;
 363
 364     snew(w_rls, bone);
 365     for (int i = 0; i < bone; i++) {
 366         snew(w_rls[i], index.size());
 367         for (int j = 0; j < index.size(); j++) {
 368             if (j >= i && j <= i + domain_min_size - 1) {
 369                 w_rls[i][j] = 1;
 370             } else {
 371                 w_rls[i][j] = 0;
 372             }
 373         }
 374     }
 375
 376     /*
 377     rvec **etalon;
 378     snew(etalon, bone);
 379     for (int i = 0; i < bone; i++) {
 380         snew(etalon[i], index.size());
 381         for (int j = 0; j < index.size(); j++) {
 382             copy_rvec(trajectory[etalon_frame][j].as_vec(), etalon[i][j]);
 383         }
 384     }
 385     */
 386
 387     for (int i = 0; i < bone; i++) {
 388         rvec *etalon;
 389         snew(etalon, index.size());
 390         for (int j = 0; j < index.size(); j++) {
 391             copy_rvec(trajectory[etalon_frame][j].as_vec(), etalon[j]);
 392         }
 393         reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[i]);
 394         make_graph(index.size(), etalon, graph[i]);
 395         sfree(etalon);
 396     }
 397
 398     std::chrono::time_point<std::chrono::system_clock> start1, end1;
 399     start1 = std::chrono::system_clock::now();
 400
 401     /*
 402     rvec **traj;
 403     snew(traj, bone);
 404
 405     for (int i = 1; i < frames; i++) {
 406         for (int j = 0; j < bone; j++) {
 407             sfree(traj[j]);
 408             snew(traj[j], index.size());
 409             for (int k = 0; k < index.size(); k++) {
 410                 copy_rvec(trajectory[i][k].as_vec(), traj[j][k]);
 411             }
 412         }
 413         #pragma omp parallel
 414         {
 415             #pragma omp for schedule(dynamic)
 416             for (int j = 0; j < bone; j++) {
 417                 //#pragma omp ordered
 418                 //{
 419                 reset_x(index.size(), NULL, index.size(), NULL, etalon[j], w_rls[j]);
 420                 reset_x(index.size(), NULL, index.size(), NULL, traj[j], w_rls[j]);
 421                 do_fit(index.size(), w_rls[j], etalon[j], traj[j]);
 422                 update_graph(graph[j], traj[j], epsi);
 423                 //}
 424             }
 425         }
 426         std::cout << i << " / " << frames << " processed\n";
 427         if (i % 25 == 0) {
 428             end1 = std::chrono::system_clock::now();
 429             int seconds = 0, minutes = 0, hours = 0;
 430             seconds = std::chrono::duration_cast<std::chrono::seconds>(end1-start1).count();
 431             seconds = (float)seconds * (float)( (float)(frames - i) / (float)i );
 432             hours = seconds / 3600;
 433             seconds -= hours * 3600;
 434             minutes = seconds / 60;
 435             seconds %= 60;
 436             std::cout << "foretold time: " << hours << " hour(s) " << minutes << " minute(s) " << seconds << "second(s)\n";
 437         }
 438     }
 439     */
 440
 441     for (int i = 1; i < frames; i++) {
 442         #pragma omp parallel
 443         {
 444             #pragma omp for schedule(dynamic)
 445             for (int j = 0; j < bone; j++) {
 446                 rvec *etalon, *traj;
 447                 snew(etalon, index.size());
 448                 for (int k = 0; k < index.size(); k++) {
 449                     copy_rvec(trajectory[etalon_frame][k].as_vec(), etalon[k]);
 450                 }
 451                 snew(traj, index.size());
 452                 for (int k = 0; k < index.size(); k++) {
 453                     copy_rvec(trajectory[i][k].as_vec(), traj[k]);
 454                 }
 455                 reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[j]);
 456                 reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[j]);
 457                 do_fit(index.size(), w_rls[j], etalon, traj);
 458                 update_graph(graph[j], traj, epsi);
 459                 sfree(etalon);
 460                 sfree(traj);
 461             }
 462         }
 463         std::cout << i << " / " << frames << " processed\n";
 464         if (i % 25 == 0) {
 465             end1 = std::chrono::system_clock::now();
 466             int seconds = 0, minutes = 0, hours = 0;
 467             seconds = std::chrono::duration_cast<std::chrono::seconds>(end1-start1).count();
 468             seconds = (float)seconds * (float)( (float)(frames - i) / (float)i );
 469             hours = seconds / 3600;
 470             seconds -= hours * 3600;
 471             minutes = seconds / 60;
 472             seconds %= 60;
 473             std::cout << "foretold time: " << hours << " hour(s) " << minutes << " minute(s) " << seconds << "second(s)\n";
 474         }
 475     }
 476     std::cout << "finale cheking\n";
 477     check_domains(delta, frames, graph);
 478     std::cout << "finding domains' sizes\n";
 479     find_domain_sizes(graph, domsizes);
 480     std::vector< int > a;
 481     a.resize(0);
 482     std::cout << "finding domains\n";
 483     while (check_domsizes(domsizes, domain_min_size)) {
 484         domains.push_back(a);
 485         get_maxsized_domain(domains.back(), graph, domsizes);
 486         delete_domain_from_graph(graph, domains.back());
 487         domsizes.resize(0);
 488         find_domain_sizes(graph, domsizes);
 489     }
 490     for (int i = 0; i < bone; i++) {
 491         sfree(w_rls[i]);
 492         //sfree(etalon[i]);
 493         //sfree(traj[i]);
 494     }
 495     sfree(w_rls);
 496     //sfree(etalon);
 497     //sfree(traj);
 498 }
 499
 500 void
 501 Domains::writeOutput()
 502 {
 503     std::cout << "making output file\n";
 504     print_domains(domains, index, fnNdx_); // see function for details | numbers from index
 505     std::cout << "\n END \n";
 506 }
 507
 508
 509 /*! \brief
 510  * The main function for the analysis template.
 511  */
 512 int
 513 main(int argc, char *argv[])
 514 {
 515     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
 516 }