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