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 "domaintype.h"
  44 #include "newfit.h"
  45
  46 #include <gromacs/trajectoryanalysis.h>
  47 #include <gromacs/trajectoryanalysis/topologyinformation.h>
  48 //#include "/home/titov_ai/Develop/gromacs-original/src/gromacs/trajectoryanalysis/topologyinformation.h"
  49
  50 #include <iostream>
  51 //#include <chrono>
  52 #include <cfloat>
  53 #include <omp.h>
  54 #include <vector>
  55 #include <cstdio>
  56
  57 using namespace gmx;
  58
  59 using gmx::RVec;
  60
  61 /*! \brief
  62  * \ingroup module_trajectoryanalysis
  63  */
  64 class Domains : public TrajectoryAnalysisModule
  65 {
  66     public:
  67
  68         Domains();
  69         virtual ~Domains();
  70
  71         //! Set the options and setting
  72         virtual void initOptions(IOptionsContainer          *options,
  73                                  TrajectoryAnalysisSettings *settings);
  74
  75         //! First routine called by the analysis framework
  76         // virtual void initAnalysis(const t_trxframe &fr, t_pbc *pbc);
  77         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
  78                                   const TopologyInformation        &top);
  79
  80         //! Call for each frame of the trajectory
  81         // virtual void analyzeFrame(const t_trxframe &fr, t_pbc *pbc);
  82         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
  83                                   TrajectoryAnalysisModuleData *pdata);
  84
  85         //! Last routine called by the analysis framework
  86         // virtual void finishAnalysis(t_pbc *pbc);
  87         virtual void finishAnalysis(int nframes);
  88
  89         //! Routine to write output, that is additional over the built-in
  90         virtual void writeOutput();
  91
  92     private:
  93
  94         std::string                                                 fnNdx_;
  95         domainType                                                  testSubject;
  96         std::vector< size_t >                                       index;
  97         std::vector< RVec >                                         trajectory;
  98         std::vector< RVec >                                         reference;
  99         Selection                                                   selec;
 100         std::vector< std::vector< std::pair< size_t, size_t > > >   fitPairs;
 101         unsigned int                                                bone;
 102         int                                                         window                      {-1}; // selectable
 103         int                                                         domain_min_size             {-1}; // selectable
 104         double                                                      delta                       {-1}; // selectable
 105         double                                                      epsi                        {-1}; // selectable
 106         int                                                         twStep                      {-1}; // selectable
 107         int                                                         DomainSearchingAlgorythm    {-1}; // selectable
 108         // Copy and assign disallowed by base.
 109 };
 110
 111 Domains::Domains(): TrajectoryAnalysisModule()
 112 {
 113 }
 114
 115 Domains::~Domains()
 116 {
 117 }
 118
 119 void
 120 Domains::initOptions(IOptionsContainer          *options,
 121                      TrajectoryAnalysisSettings *settings)
 122 {
 123     static const char *const desc[] = {
 124         "[THISMODULE] to be done"
 125         "-f trajectory.xtc -s structure.tpr -n index.ndx -sf selectionListFile -on outPutName -dms 4+ -DSA 0/1/2 -eps 0.2 -dlt 0.95 -wSize 1000 -twStep 100"
 126     };
 127     // Add the descriptive text (program help text) to the options
 128     settings->setHelpText(desc);
 129     // Add option for selecting a subset of atoms
 130     options->addOption(SelectionOption("select")
 131                            .store(&selec).required()
 132                            .description("Atoms that are considered as part of the excluded volume."));
 133     // Add option for output file name
 134     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
 135                             .store(&fnNdx_).defaultBasename("domains")
 136                             .description("Index file for the algorithm output."));
 137     // Add option for domain min size constant
 138     options->addOption(gmx::IntegerOption("dms")
 139                             .store(&domain_min_size)
 140                             .description("minimum domain size, should be >= 4."));
 141     // Add option for Domain's Searching Algorythm
 142     options->addOption(gmx::IntegerOption("DSA")
 143                             .store(&DomainSearchingAlgorythm)
 144                             .description("Domain's Searching Algorythm: 0 == (default) from bigger to smaller | 1 == just the first biggest domain | 2 == from smaller to bigger."));
 145     // Add option for epsi constant
 146     options->addOption(DoubleOption("eps")
 147                             .store(&epsi)
 148                             .description("thermal vibrations' constant."));
 149     // Add option for delta constant
 150     options->addOption(DoubleOption("dlt")
 151                             .store(&delta)
 152                             .description("domain membership probability."));
 153     // Add option for domain min size constant
 154     options->addOption(gmx::IntegerOption("wSize")
 155                             .store(&window)
 156                             .description("flowing window to evaluate domains from."));
 157     // Add option for time step between windows' starts
 158     options->addOption(gmx::IntegerOption("twStep")
 159                             .store(&twStep)
 160                             .description("time step between windows' starting positions."));
 161     // Control input settings
 162     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 163     settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 164     //settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 165     settings->setPBC(true);
 166 }
 167
 168 // -f trajectory.xtc -s structure.tpr -n index.ndx -sf selectionListFile -on outPutName -dms 4+ -DSA 0/1/2 -eps 0.2 -dlt 0.95 -wSize 1000 -twStep 100
 169
 170 void
 171 Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 172                       const TopologyInformation        &top)
 173 {
 174     // считывание индекса
 175     index.resize(0);
 176     for (ArrayRef< const int >::iterator ai {selec.atomIndices().begin()}; (ai < selec.atomIndices().end()); ++ai) {
 177         index.push_back(static_cast< size_t >(*ai));
 178     }
 179
 180     // считывание референсной структуры
 181     reference.resize(0);
 182
 183     if (top.hasFullTopology()) {
 184         for (const auto &i : index) {
 185             reference.push_back(top.x().at(i));
 186         }
 187     }
 188
 189     // вычисление числа основных (от основание) последовательностей
 190     // в идеале нужно брать все сочетания, но ограничиваемся последовательными наборами
 191     bone = static_cast< size_t >(index.size() - static_cast< size_t >(domain_min_size) + 1);
 192
 193     // создание пар для фитирования структур
 194     fitPairs.resize(0);
 195     fitPairs.resize(bone);
 196     for (size_t i {0}; i < bone; ++i) {
 197         fitPairs[i].resize(0);
 198         for (size_t j {0}; j < domain_min_size; ++j) {
 199             fitPairs[i].push_back(std::make_pair(j + i, j + i));
 200         }
 201     }
 202
 203     trajectory.resize(0);
 204     trajectory.resize(index.size());
 205
 206     // задание необходимых параметров для вычисления доменов
 207     testSubject.setDefaults(index, reference, window, domain_min_size, DomainSearchingAlgorythm, twStep, bone, epsi, delta, fnNdx_);
 208 }
 209
 210 void
 211 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
 212 {
 213     // считывания текущего фрейма траектории
 214     for (size_t i {0}; i < index.size(); ++i) {
 215         trajectory[i] = fr.x[index[i]];
 216     }
 217
 218     // создание копий фрейма для каждой основной последовательности
 219     std::vector< std::vector< RVec > > tsTemp;
 220     tsTemp.resize(0);
 221     tsTemp.resize(bone, trajectory);
 222
 223     // фитирование каждой копии
 224     #pragma omp parallel for ordered schedule(dynamic) firstprivate(reference)
 225     for (size_t i = 0; i < bone; ++i) {
 226         MyFitNew(reference, tsTemp[i], fitPairs[i], 0.000'001);
 227     }
 228     #pragma omp barrier
 229
 230     // (обновление/потенциальное выделение) доменов
 231     testSubject.update(tsTemp, frnr);
 232 }
 233
 234 void
 235 Domains::finishAnalysis(int nframes)
 236 {
 237
 238 }
 239
 240 void
 241 Domains::writeOutput()
 242 {
 243     std::cout << "\n END \n";
 244 }
 245
 246 /*! \brief
 247  * The main function for the analysis template.
 248  */
 249 int
 250 main(int argc, char *argv[])
 251 {
 252     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<Domains>(argc, argv);
 253 }