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