src/gromacs/trajectoryanalysis/modules/freevolume.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014, 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 David van der Spoel <david.vanderspoel@icm.uu.se>
  40  * \ingroup module_trajectoryanalysis
  41  */
  42 #include "freevolume.h"
  43
  44 #include <string>
  45
  46 #include "gromacs/legacyheaders/atomprop.h"
  47 #include "gromacs/legacyheaders/copyrite.h"
  48 #include "gromacs/random/random.h"
  49 #include "gromacs/legacyheaders/pbc.h"
  50 #include "gromacs/legacyheaders/vec.h"
  51
  52 #include "gromacs/analysisdata/analysisdata.h"
  53 #include "gromacs/analysisdata/modules/average.h"
  54 #include "gromacs/analysisdata/modules/plot.h"
  55 #include "gromacs/options/basicoptions.h"
  56 #include "gromacs/options/filenameoption.h"
  57 #include "gromacs/options/options.h"
  58 #include "gromacs/selection/nbsearch.h"
  59 #include "gromacs/selection/selection.h"
  60 #include "gromacs/selection/selectionoption.h"
  61 #include "gromacs/trajectoryanalysis/analysissettings.h"
  62 #include "gromacs/utility/arrayref.h"
  63 #include "gromacs/utility/exceptions.h"
  64
  65 namespace gmx
  66 {
  67
  68 namespace analysismodules
  69 {
  70
  71 namespace
  72 {
  73
  74 /*! \brief
  75  * Class used to compute free volume in a simulations box.
  76  *
  77  * Inherits TrajectoryAnalysisModule and all functions from there.
  78  * Does not implement any new functionality.
  79  *
  80  * \ingroup module_trajectoryanalysis
  81  */
  82 class FreeVolume : public TrajectoryAnalysisModule
  83 {
  84     public:
  85         //! Constructor
  86         FreeVolume();
  87
  88         //! Destructor
  89         virtual ~FreeVolume();
  90
  91         //! Set the options and setting
  92         virtual void initOptions(Options                    *options,
  93                                  TrajectoryAnalysisSettings *settings);
  94
  95         //! First routine called by the analysis framework
  96         virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
  97                                   const TopologyInformation        &top);
  98
  99         //! Call for each frame of the trajectory
 100         virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 101                                   TrajectoryAnalysisModuleData *pdata);
 102
 103         //! Last routine called by the analysis framework
 104         virtual void finishAnalysis(int nframes);
 105
 106         //! Routine to write output, that is additional over the built-in
 107         virtual void writeOutput();
 108
 109     private:
 110         std::string                       fnFreevol_;
 111         Selection                         sel_;
 112         AnalysisData                      data_;
 113         AnalysisDataAverageModulePointer  adata_;
 114
 115         int                               nmol_;
 116         double                            mtot_;
 117         double                            cutoff_;
 118         double                            probeRadius_;
 119         gmx_rng_t                         rng_;
 120         int                               seed_, ninsert_;
 121         AnalysisNeighborhood              nb_;
 122         //! The van der Waals radius per atom
 123         std::vector<double>               vdw_radius_;
 124
 125         // Copy and assign disallowed by base.
 126 };
 127
 128 // Constructor. Here it is important to initialize the pointer to
 129 // subclasses that are elements of the main class. Here we have only
 130 // one. The type of this depends on what kind of tool you need.
 131 // Here we only have simple value/time kind of data.
 132 FreeVolume::FreeVolume()
 133     : TrajectoryAnalysisModule(FreeVolumeInfo::name, FreeVolumeInfo::shortDescription),
 134       adata_(new AnalysisDataAverageModule())
 135 {
 136     // We only compute two numbers per frame
 137     data_.setColumnCount(0, 2);
 138     // Tell the analysis framework that this component exists
 139     registerAnalysisDataset(&data_, "freevolume");
 140     rng_         = NULL;
 141     nmol_        = 0;
 142     mtot_        = 0;
 143     cutoff_      = 0;
 144     probeRadius_ = 0;
 145     seed_        = -1;
 146     ninsert_     = 1000;
 147 }
 148
 149
 150 FreeVolume::~FreeVolume()
 151 {
 152     // Destroy C structures where there is no automatic memory release
 153     // C++ takes care of memory in classes (hopefully)
 154     if (NULL != rng_)
 155     {
 156         gmx_rng_destroy(rng_);
 157     }
 158 }
 159
 160
 161 void
 162 FreeVolume::initOptions(Options                    *options,
 163                         TrajectoryAnalysisSettings *settings)
 164 {
 165     static const char *const desc[] = {
 166         "[THISMODULE] calculates the free volume in a box as",
 167         "a function of time. The free volume is",
 168         "plotted as a fraction of the total volume.",
 169         "The program tries to insert a probe with a given radius,",
 170         "into the simulations box and if the distance between the",
 171         "probe and any atom is less than the sums of the",
 172         "van der Waals radii of both atoms, the position is",
 173         "considered to be occupied, i.e. non-free. By using a",
 174         "probe radius of 0, the true free volume is computed.",
 175         "By using a larger radius, e.g. 0.14 nm, roughly corresponding",
 176         "to a water molecule, the free volume for a hypothetical",
 177         "particle with that size will be produced.",
 178         "Note however, that since atoms are treated as hard-spheres",
 179         "these number are very approximate, and typically only",
 180         "relative changes are meaningful, for instance by doing a",
 181         "series of simulations at different temperature.[PAR]",
 182         "The group specified by the selection is considered to",
 183         "delineate non-free volume.",
 184         "The number of insertions per unit of volume is important",
 185         "to get a converged result. About 1000/nm^3 yields an overall",
 186         "standard deviation that is determined by the fluctuations in",
 187         "the trajectory rather than by the fluctuations due to the",
 188         "random numbers.[PAR]",
 189         "The results are critically dependent on the van der Waals radii;",
 190         "we recommend to use the values due to Bondi (1964).[PAR]",
 191         "The Fractional Free Volume (FFV) that some authors like to use",
 192         "is given by 1 - 1.3*(1-Free Volume). This value is printed on",
 193         "the terminal."
 194     };
 195
 196     // Add the descriptive text (program help text) to the options
 197     options->setDescription(desc);
 198
 199     // Add option for optional output file
 200     options->addOption(FileNameOption("o").filetype(eftPlot).outputFile()
 201                            .store(&fnFreevol_).defaultBasename("freevolume")
 202                            .description("Computed free volume"));
 203
 204     // Add option for selecting a subset of atoms
 205     options->addOption(SelectionOption("select")
 206                            .store(&sel_).defaultSelectionText("all")
 207                            .onlyAtoms());
 208
 209     // Add option for the probe radius and initialize it
 210     options->addOption(DoubleOption("radius").store(&probeRadius_)
 211                            .description("Radius of the probe to be inserted (nm, 0 yields the true free volume)"));
 212
 213     // Add option for the random number seed and initialize it to
 214     // generate a value automatically
 215     options->addOption(IntegerOption("seed").store(&seed_)
 216                            .description("Seed for random number generator."));
 217
 218     // Add option to determine number of insertion trials per frame
 219     options->addOption(IntegerOption("ninsert").store(&ninsert_)
 220                            .description("Number of probe insertions per cubic nm to try for each frame in the trajectory."));
 221
 222     // Control input settings
 223     settings->setFlags(TrajectoryAnalysisSettings::efRequireTop |
 224                        TrajectoryAnalysisSettings::efNoUserPBC);
 225     settings->setPBC(true);
 226 }
 227
 228
 229 void
 230 FreeVolume::initAnalysis(const TrajectoryAnalysisSettings &settings,
 231                          const TopologyInformation        &top)
 232 {
 233     // Add the module that will contain the averaging and the time series
 234     // for our calculation
 235     data_.addModule(adata_);
 236
 237     // Add a module for plotting the data automatically at the end of
 238     // the calculation. With this in place you only have to add data
 239     // points to the data et.
 240     AnalysisDataPlotModulePointer plotm_(new AnalysisDataPlotModule());
 241     plotm_->setSettings(settings.plotSettings());
 242     plotm_->setFileName(fnFreevol_);
 243     plotm_->setTitle("Free Volume");
 244     plotm_->setXAxisIsTime();
 245     plotm_->setYLabel("Free Volume (%)");
 246     plotm_->appendLegend("Free Volume");
 247     plotm_->appendLegend("Volume");
 248
 249     data_.addModule(plotm_);
 250
 251     // Initiate variable
 252     cutoff_               = 0;
 253     int            nnovdw = 0;
 254     gmx_atomprop_t aps    = gmx_atomprop_init();
 255     t_atoms       *atoms  = &(top.topology()->atoms);
 256
 257     // Compute total mass
 258     mtot_ = 0;
 259     for (int i = 0; (i < atoms->nr); i++)
 260     {
 261         mtot_ += atoms->atom[i].m;
 262     }
 263
 264     // Extracts number of molecules
 265     nmol_ = top.topology()->mols.nr;
 266
 267     // Loop over atoms in the selection using an iterator
 268     const int          maxnovdw = 10;
 269     ConstArrayRef<int> atomind  = sel_.atomIndices();
 270     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ++ai)
 271     {
 272         // Dereference the iterator to obtain an atom number
 273         int  i = *ai;
 274         real value;
 275
 276         // Lookup the Van der Waals radius of this atom
 277         int resnr = atoms->atom[i].resind;
 278         if (TRUE == gmx_atomprop_query(aps, epropVDW,
 279                                        *(atoms->resinfo[resnr].name),
 280                                        *(atoms->atomname[i]),
 281                                        &value))
 282         {
 283             vdw_radius_.push_back(value);
 284             if (value > cutoff_)
 285             {
 286                 cutoff_ = value;
 287             }
 288         }
 289         else
 290         {
 291             nnovdw++;
 292             if (nnovdw < maxnovdw)
 293             {
 294                 fprintf(stderr, "Could not determine VDW radius for %s-%s. Set to zero.\n",
 295                         *(atoms->resinfo[resnr].name),
 296                         *(atoms->atomname[i]));
 297             }
 298             vdw_radius_.push_back(0.0);
 299         }
 300     }
 301     gmx_atomprop_destroy(aps);
 302
 303     // Increase cutoff by proberadius to make sure we do not miss
 304     // anything
 305     cutoff_ += probeRadius_;
 306
 307     if (nnovdw >= maxnovdw)
 308     {
 309         fprintf(stderr, "Could not determine VDW radius for %d particles. These were set to zero.\n", nnovdw);
 310     }
 311
 312     // Print parameters to output. Maybe should make dependent on
 313     // verbosity flag?
 314     printf("cutoff       = %g nm\n", cutoff_);
 315     printf("probe_radius = %g nm\n", probeRadius_);
 316     printf("seed         = %d\n", seed_);
 317     printf("ninsert      = %d probes per nm^3\n", ninsert_);
 318
 319     // Initiate the random number generator
 320     rng_ = gmx_rng_init(seed_);
 321
 322     // Initiate the neighborsearching code
 323     nb_.setCutoff(cutoff_);
 324 }
 325
 326 void
 327 FreeVolume::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 328                          TrajectoryAnalysisModuleData *pdata)
 329 {
 330     AnalysisDataHandle       dh   = pdata->dataHandle(data_);
 331     const Selection         &sel  = pdata->parallelSelection(sel_);
 332
 333     GMX_RELEASE_ASSERT(NULL != pbc, "You have no periodic boundary conditions");
 334
 335     // Analysis framework magic
 336     dh.startFrame(frnr, fr.time);
 337
 338     // Compute volume and number of insertions to perform
 339     real V       = det(fr.box);
 340     int  Ninsert = static_cast<int>(ninsert_*V);
 341
 342     // Use neighborsearching tools!
 343     AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, sel);
 344
 345     // Then loop over insertions
 346     int NinsTot = 0;
 347     for (int i = 0; (i < Ninsert); i++)
 348     {
 349         rvec rand, ins, dx;
 350
 351         for (int m = 0; (m < DIM); m++)
 352         {
 353             // Generate random number between 0 and 1
 354             rand[m] = gmx_rng_uniform_real(rng_);
 355         }
 356         // Generate random 3D position within the box
 357         mvmul(fr.box, rand, ins);
 358
 359         // Find the first reference position within the cutoff.
 360         bool                           bOverlap = false;
 361         AnalysisNeighborhoodPair       pair;
 362         AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(ins);
 363         while (!bOverlap && pairSearch.findNextPair(&pair))
 364         {
 365             int jp = pair.refIndex();
 366             // Compute distance vector to first atom in the neighborlist
 367             pbc_dx(pbc, ins, sel.position(jp).x(), dx);
 368
 369             // See whether the distance is smaller than allowed
 370             bOverlap = (norm(dx) <
 371                         probeRadius_+vdw_radius_[sel.position(jp).refId()]);
 372
 373         }
 374
 375         if (!bOverlap)
 376         {
 377             // We found some free volume!
 378             NinsTot++;
 379         }
 380     }
 381     // Compute total free volume for this frame
 382     double frac = 0;
 383     if (Ninsert > 0)
 384     {
 385         frac = (100.0*NinsTot)/Ninsert;
 386     }
 387     // Add the free volume fraction to the data set in column 0
 388     dh.setPoint(0, frac);
 389     // Add the total volume to the data set in column 1
 390     dh.setPoint(1, V);
 391
 392     // Magic
 393     dh.finishFrame();
 394 }
 395
 396
 397 void
 398 FreeVolume::finishAnalysis(int /* nframes */)
 399 {
 400     please_cite(stdout, "Bondi1964a");
 401     please_cite(stdout, "Lourenco2013a");
 402 }
 403
 404 void
 405 FreeVolume::writeOutput()
 406 {
 407     // Final results come from statistics module in analysis framework
 408     double FVaver  = adata_->average(0, 0);
 409     double FVerror = adata_->standardDeviation(0, 0);
 410     printf("Free volume %.2f +/- %.2f %%\n", FVaver, FVerror);
 411
 412     double Vaver  = adata_->average(0, 1);
 413     double Verror = adata_->standardDeviation(0, 1);
 414     printf("Total volume %.2f +/- %.2f nm^3\n", Vaver, Verror);
 415
 416     printf("Number of molecules %d total mass %.2f Dalton\n", nmol_, mtot_);
 417     double RhoAver  = mtot_ / (Vaver * 1e-24 * AVOGADRO);
 418     double RhoError = sqr(RhoAver / Vaver)*Verror;
 419     printf("Average molar mass: %.2f Dalton\n", mtot_/nmol_);
 420
 421     double VmAver  = Vaver/nmol_;
 422     double VmError = Verror/nmol_;
 423     printf("Density rho: %.2f +/- %.2f nm^3\n", RhoAver, RhoError);
 424     printf("Molecular volume Vm assuming homogeneity: %.4f +/- %.4f nm^3\n",
 425            VmAver, VmError);
 426
 427     double VvdWaver  = (1-FVaver/100)*VmAver;
 428     double VvdWerror = 0;
 429     printf("Molecular van der Waals volume assuming homogeneity:  %.4f +/- %.4f nm^3\n",
 430            VvdWaver, VvdWerror);
 431
 432     double FFVaver  = 1-1.3*((100-FVaver)/100);
 433     double FFVerror = (FVerror/FVaver)*FFVaver;
 434     printf("Fractional free volume %.3f +/- %.3f\n", FFVaver, FFVerror);
 435 }
 436
 437 }       // namespace
 438
 439 const char FreeVolumeInfo::name[]             = "freevolume";
 440 const char FreeVolumeInfo::shortDescription[] =
 441     "Calculate free volume";
 442
 443 TrajectoryAnalysisModulePointer FreeVolumeInfo::create()
 444 {
 445     return TrajectoryAnalysisModulePointer(new FreeVolume);
 446 }
 447
 448 } // namespace analysismodules
 449
 450 } // namespace gmx