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