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