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