src/rcore.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2011,2012,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
  36 #include <fstream>
  37 #include <omp.h>
  38 #include <iostream>
  39 #include <iomanip>
  40
  41 #include <newfit.h>
  42
  43 #include <gromacs/trajectoryanalysis/topologyinformation.h>
  44
  45 /*! \brief
  46  * Template class to serve as a basis for user analysis tools.
  47  */
  48 class AnalysisTemplate : public TrajectoryAnalysisModule
  49 {
  50     public:
  51
  52         virtual void initOptions(IOptionsContainer          *options,
  53                                  TrajectoryAnalysisSettings *settings);
  54         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
  55                                   const TopologyInformation        &top);
  56
  57         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
  58                                   TrajectoryAnalysisModuleData *pdata);
  59
  60         virtual void finishAnalysis(int nframes);
  61         virtual void writeOutput();
  62
  63     private:
  64         //class ModuleData;
  65
  66         std::string                         outputFName;        // selectable
  67         Selection                           sel_;
  68         std::vector< RVec >                 reference;
  69         std::vector< std::vector < RVec > > trajectory;
  70         std::vector< int >                  index;
  71         double                              epsi        {0.25};  // selectable
  72         double                              prt         {0.66}; // selectable
  73 };
  74
  75 void
  76 AnalysisTemplate::initOptions(IOptionsContainer          *options,
  77                               TrajectoryAnalysisSettings *settings)
  78 {
  79     static const char *const desc[] = {
  80         "Analysis tool for finding molecular core."
  81     };
  82
  83     // Add the descriptive text (program help text) to the options
  84     settings->setHelpText(desc);
  85     // Add option for selecting a subset of atoms
  86     options->addOption(SelectionOption("select")
  87                         .store(&sel_).required()
  88                         .description("Atoms that are considered as part of the excluded volume"));
  89     // Add option for output file name
  90     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
  91                         .store(&outputFName).defaultBasename("rcore")
  92                         .description("Index file from the rcore"));
  93     // Add option for epsi constant
  94     options->addOption(DoubleOption("eps")
  95                         .store(&epsi)
  96                         .description("thermal vibrations' constant"));
  97     // Add option for particion constant
  98     options->addOption(DoubleOption("-prt")
  99                         .store(&prt)
 100                         .description("blob particion"));
 101     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
 102     settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
 103     settings->setPBC(true);
 104 }
 105
 106 void
 107 AnalysisTemplate::initAnalysis(const TrajectoryAnalysisSettings &settings,
 108                                const TopologyInformation        &top)
 109 {
 110     index.resize(0);
 111     for (ArrayRef< const int >::iterator ai {sel_.atomIndices().begin()}; ai < sel_.atomIndices().end(); ++ai) {
 112         index.push_back(*ai);
 113     }
 114
 115     reference.resize(0);
 116     if (top.hasFullTopology()) {
 117         for (size_t i {0}; i < index.size(); ++i) {
 118             reference.push_back(top.x().at(index[i]));
 119         }
 120     }
 121
 122     trajectory.resize(0);
 123 }
 124
 125 void
 126 AnalysisTemplate::analyzeFrame(int                          frnr,
 127                                const t_trxframe             &fr,
 128                                t_pbc                        *pbc,
 129                                TrajectoryAnalysisModuleData *pdata)
 130 {
 131     trajectory.resize(trajectory.size() + 1);
 132     trajectory.back().resize(0);
 133     for (size_t i {0}; i < index.size(); ++i)
 134     {
 135         trajectory.back().push_back(fr.x[index[i]]);
 136     }
 137 }
 138
 139 void
 140 AnalysisTemplate::finishAnalysis(int nframes)
 141 {
 142     std::vector< std::pair< size_t, size_t > > fitPairs;
 143     fitPairs.resize(0);
 144     for (size_t i {0}; i < index.size(); ++i) {
 145         fitPairs.push_back(std::make_pair(i, i));
 146     }
 147
 148     std::vector< std::pair< size_t, double > > noise;
 149
 150     size_t old = index.size(), nw = 0;
 151     while (old != nw) {
 152         noise.resize(0);
 153         for (size_t i {0}; i < index.size(); ++i) {
 154             noise.push_back(std::make_pair(i, 0.));
 155         }
 156         for (size_t i {0}; i < trajectory.size(); ++i) {
 157             MyFitNew(reference, trajectory[i], fitPairs, 0.000'001);
 158         }
 159         for (size_t i {0}; i < trajectory.size(); ++i) {
 160             for (size_t j {0}; j < trajectory[i].size(); ++j) {
 161                 noise[j].second += (reference[j] - trajectory[i][j]).norm();
 162             }
 163         }
 164         for (size_t i {0}; i < noise.size(); ++i) {
 165             noise[i].second /= trajectory.size();
 166         }
 167         std::sort(noise.begin(), noise.end(), [](const std::pair< size_t, double > a, const std::pair< size_t, double > b){ return a.second < b.second;});
 168         if (size_t temp {fitPairs.size() * prt}; noise[temp].second > epsi) {
 169             fitPairs.resize(0);
 170             for (size_t i {0}; i <= temp; ++i) {
 171                 fitPairs.push_back(std::make_pair(noise[i].first, noise[i].first));
 172             }
 173             old = nw;
 174             nw = temp;
 175         } else {
 176             old = nw;
 177             std::ofstream file(outputFName, std::ofstream::app);
 178             file << "[ rcore ]\n";
 179             for (size_t i {0}; i < temp; ++i) {
 180                 file << std::setw(7) << noise[i].first;
 181                 if (temp % 7 == 0) {
 182                     file << "\n";
 183                 }
 184             }
 185         }
 186     }
 187 }
 188
 189 void
 190 AnalysisTemplate::writeOutput()
 191 {
 192
 193 }
 194
 195 /*! \brief
 196  * The main function for the analysis template.
 197  */
 198 int
 199 main(int argc, char *argv[])
 200 {
 201     return gmx::TrajectoryAnalysisCommandLineRunner::runAsMain<AnalysisTemplate>(argc, argv);
 202 }