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,2016,2017 by the GROMACS development team.
   5  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Implements gmx::analysismodules::Freevolume.
  39  *
  40  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  41  * \ingroup module_trajectoryanalysis
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "freevolume.h"
  46
  47 #include <string>
  48
  49 #include "gromacs/analysisdata/analysisdata.h"
  50 #include "gromacs/analysisdata/modules/average.h"
  51 #include "gromacs/analysisdata/modules/plot.h"
  52 #include "gromacs/math/functions.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/math/vec.h"
  55 #include "gromacs/options/basicoptions.h"
  56 #include "gromacs/options/filenameoption.h"
  57 #include "gromacs/options/ioptionscontainer.h"
  58 #include "gromacs/pbcutil/pbc.h"
  59 #include "gromacs/random/threefry.h"
  60 #include "gromacs/random/uniformrealdistribution.h"
  61 #include "gromacs/selection/nbsearch.h"
  62 #include "gromacs/selection/selection.h"
  63 #include "gromacs/selection/selectionoption.h"
  64 #include "gromacs/topology/atomprop.h"
  65 #include "gromacs/topology/mtop_util.h"
  66 #include "gromacs/topology/topology.h"
  67 #include "gromacs/trajectory/trajectoryframe.h"
  68 #include "gromacs/trajectoryanalysis/analysissettings.h"
  69 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  70 #include "gromacs/utility/arrayref.h"
  71 #include "gromacs/utility/exceptions.h"
  72 #include "gromacs/utility/pleasecite.h"
  73
  74 namespace gmx
  75 {
  76
  77 namespace analysismodules
  78 {
  79
  80 namespace
  81 {
  82
  83 /*! \brief
  84  * Class used to compute free volume in a simulations box.
  85  *
  86  * Inherits TrajectoryAnalysisModule and all functions from there.
  87  * Does not implement any new functionality.
  88  *
  89  * \ingroup module_trajectoryanalysis
  90  */
  91 class FreeVolume : public TrajectoryAnalysisModule
  92 {
  93 public:
  94     FreeVolume();
  95     ~FreeVolume() override {}
  96
  97     void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
  98     void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
  99     void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
 100     void finishAnalysis(int nframes) override;
 101     void writeOutput() override;
 102
 103 private:
 104     std::string                      fnFreevol_;
 105     Selection                        sel_;
 106     AnalysisData                     data_;
 107     AnalysisDataAverageModulePointer adata_;
 108
 109     int                      nmol_;
 110     double                   mtot_;
 111     double                   cutoff_;
 112     double                   probeRadius_;
 113     gmx::DefaultRandomEngine rng_;
 114     int                      seed_, ninsert_;
 115     AnalysisNeighborhood     nb_;
 116     //! The van der Waals radius per atom
 117     std::vector<double> vdw_radius_;
 118
 119     // Copy and assign disallowed by base.
 120 };
 121
 122 // Constructor. Here it is important to initialize the pointer to
 123 // subclasses that are elements of the main class. Here we have only
 124 // one. The type of this depends on what kind of tool you need.
 125 // Here we only have simple value/time kind of data.
 126 FreeVolume::FreeVolume() : adata_(new AnalysisDataAverageModule())
 127 {
 128     // We only compute two numbers per frame
 129     data_.setColumnCount(0, 2);
 130     // Tell the analysis framework that this component exists
 131     registerAnalysisDataset(&data_, "freevolume");
 132     nmol_        = 0;
 133     mtot_        = 0;
 134     cutoff_      = 0;
 135     probeRadius_ = 0;
 136     seed_        = 0;
 137     ninsert_     = 1000;
 138 }
 139
 140
 141 void FreeVolume::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 142 {
 143     static const char* const desc[] = {
 144         "[THISMODULE] calculates the free volume in a box as",
 145         "a function of time. The free volume is",
 146         "plotted as a fraction of the total volume.",
 147         "The program tries to insert a probe with a given radius,",
 148         "into the simulations box and if the distance between the",
 149         "probe and any atom is less than the sums of the",
 150         "van der Waals radii of both atoms, the position is",
 151         "considered to be occupied, i.e. non-free. By using a",
 152         "probe radius of 0, the true free volume is computed.",
 153         "By using a larger radius, e.g. 0.14 nm, roughly corresponding",
 154         "to a water molecule, the free volume for a hypothetical",
 155         "particle with that size will be produced.",
 156         "Note however, that since atoms are treated as hard-spheres",
 157         "these number are very approximate, and typically only",
 158         "relative changes are meaningful, for instance by doing a",
 159         "series of simulations at different temperature.[PAR]",
 160         "The group specified by the selection is considered to",
 161         "delineate non-free volume.",
 162         "The number of insertions per unit of volume is important",
 163         "to get a converged result. About 1000/nm^3 yields an overall",
 164         "standard deviation that is determined by the fluctuations in",
 165         "the trajectory rather than by the fluctuations due to the",
 166         "random numbers.[PAR]",
 167         "The results are critically dependent on the van der Waals radii;",
 168         "we recommend to use the values due to Bondi (1964).[PAR]",
 169         "The Fractional Free Volume (FFV) that some authors like to use",
 170         "is given by 1 - 1.3*(1-Free Volume). This value is printed on",
 171         "the terminal."
 172     };
 173
 174     settings->setHelpText(desc);
 175
 176     // Add option for optional output file
 177     options->addOption(FileNameOption("o")
 178                                .filetype(eftPlot)
 179                                .outputFile()
 180                                .store(&fnFreevol_)
 181                                .defaultBasename("freevolume")
 182                                .description("Computed free volume"));
 183
 184     // Add option for selecting a subset of atoms
 185     options->addOption(
 186             SelectionOption("select").store(&sel_).defaultSelectionText("all").onlyAtoms().description(
 187                     "Atoms that are considered as part of the excluded volume"));
 188
 189     // Add option for the probe radius and initialize it
 190     options->addOption(
 191             DoubleOption("radius").store(&probeRadius_).description("Radius of the probe to be inserted (nm, 0 yields the true free volume)"));
 192
 193     // Add option for the random number seed and initialize it to
 194     // generate a value automatically
 195     options->addOption(IntegerOption("seed").store(&seed_).description(
 196             "Seed for random number generator (0 means generate)."));
 197
 198     // Add option to determine number of insertion trials per frame
 199     options->addOption(IntegerOption("ninsert").store(&ninsert_).description(
 200             "Number of probe insertions per cubic nm to try for each frame in the trajectory."));
 201
 202     // Control input settings
 203     settings->setFlags(TrajectoryAnalysisSettings::efRequireTop | TrajectoryAnalysisSettings::efNoUserPBC);
 204     settings->setPBC(true);
 205 }
 206
 207
 208 void FreeVolume::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
 209 {
 210     // Add the module that will contain the averaging and the time series
 211     // for our calculation
 212     data_.addModule(adata_);
 213
 214     // Add a module for plotting the data automatically at the end of
 215     // the calculation. With this in place you only have to add data
 216     // points to the data et.
 217     AnalysisDataPlotModulePointer plotm_(new AnalysisDataPlotModule());
 218     plotm_->setSettings(settings.plotSettings());
 219     plotm_->setFileName(fnFreevol_);
 220     plotm_->setTitle("Free Volume");
 221     plotm_->setXAxisIsTime();
 222     plotm_->setYLabel("Free Volume (%)");
 223     plotm_->appendLegend("Free Volume");
 224     plotm_->appendLegend("Volume");
 225
 226     data_.addModule(plotm_);
 227
 228     // Initiate variable
 229     cutoff_               = 0;
 230     int            nnovdw = 0;
 231     AtomProperties aps;
 232     auto           atoms = top.copyAtoms();
 233
 234     // Compute total mass
 235     mtot_ = 0;
 236     for (int i = 0; (i < atoms->nr); i++)
 237     {
 238         mtot_ += atoms->atom[i].m;
 239     }
 240
 241     // Extracts number of molecules
 242     nmol_ = gmx_mtop_num_molecules(*top.mtop());
 243
 244     // Loop over atoms in the selection using an iterator
 245     const int           maxnovdw = 10;
 246     ArrayRef<const int> atomind  = sel_.atomIndices();
 247     for (ArrayRef<const int>::iterator ai = atomind.begin(); (ai < atomind.end()); ++ai)
 248     {
 249         // Dereference the iterator to obtain an atom number
 250         int  i = *ai;
 251         real value;
 252
 253         // Lookup the Van der Waals radius of this atom
 254         int resnr = atoms->atom[i].resind;
 255         if (aps.setAtomProperty(epropVDW, *(atoms->resinfo[resnr].name), *(atoms->atomname[i]), &value))
 256         {
 257             vdw_radius_.push_back(value);
 258             if (value > cutoff_)
 259             {
 260                 cutoff_ = value;
 261             }
 262         }
 263         else
 264         {
 265             nnovdw++;
 266             if (nnovdw < maxnovdw)
 267             {
 268                 fprintf(stderr, "Could not determine VDW radius for %s-%s. Set to zero.\n",
 269                         *(atoms->resinfo[resnr].name), *(atoms->atomname[i]));
 270             }
 271             vdw_radius_.push_back(0.0);
 272         }
 273     }
 274
 275     // Increase cutoff by proberadius to make sure we do not miss
 276     // anything
 277     cutoff_ += probeRadius_;
 278
 279     if (nnovdw >= maxnovdw)
 280     {
 281         fprintf(stderr, "Could not determine VDW radius for %d particles. These were set to zero.\n", nnovdw);
 282     }
 283
 284     if (seed_ == 0)
 285     {
 286         seed_ = static_cast<int>(gmx::makeRandomSeed());
 287     }
 288
 289     // Print parameters to output. Maybe should make dependent on
 290     // verbosity flag?
 291     printf("cutoff       = %g nm\n", cutoff_);
 292     printf("probe_radius = %g nm\n", probeRadius_);
 293     printf("seed         = %d\n", seed_);
 294     printf("ninsert      = %d probes per nm^3\n", ninsert_);
 295
 296     // Initiate the random number generator
 297     rng_.seed(seed_);
 298
 299     // Initiate the neighborsearching code
 300     nb_.setCutoff(cutoff_);
 301 }
 302
 303 void FreeVolume::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
 304 {
 305     AnalysisDataHandle                 dh  = pdata->dataHandle(data_);
 306     const Selection&                   sel = pdata->parallelSelection(sel_);
 307     gmx::UniformRealDistribution<real> dist;
 308
 309     GMX_RELEASE_ASSERT(nullptr != pbc, "You have no periodic boundary conditions");
 310
 311     // Analysis framework magic
 312     dh.startFrame(frnr, fr.time);
 313
 314     // Compute volume and number of insertions to perform
 315     real V       = det(fr.box);
 316     int  Ninsert = static_cast<int>(ninsert_ * V);
 317
 318     // Use neighborsearching tools!
 319     AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, sel);
 320
 321     // Then loop over insertions
 322     int NinsTot = 0;
 323     for (int i = 0; (i < Ninsert); i++)
 324     {
 325         rvec rand, ins, dx;
 326
 327         for (int m = 0; (m < DIM); m++)
 328         {
 329             // Generate random number between 0 and 1
 330             rand[m] = dist(rng_);
 331         }
 332         // Generate random 3D position within the box
 333         mvmul(fr.box, rand, ins);
 334
 335         // Find the first reference position within the cutoff.
 336         bool                           bOverlap = false;
 337         AnalysisNeighborhoodPair       pair;
 338         AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(ins);
 339         while (!bOverlap && pairSearch.findNextPair(&pair))
 340         {
 341             int jp = pair.refIndex();
 342             // Compute distance vector to first atom in the neighborlist
 343             pbc_dx(pbc, ins, sel.position(jp).x(), dx);
 344
 345             // See whether the distance is smaller than allowed
 346             bOverlap = (norm(dx) < probeRadius_ + vdw_radius_[sel.position(jp).refId()]);
 347         }
 348
 349         if (!bOverlap)
 350         {
 351             // We found some free volume!
 352             NinsTot++;
 353         }
 354     }
 355     // Compute total free volume for this frame
 356     double frac = 0;
 357     if (Ninsert > 0)
 358     {
 359         frac = (100.0 * NinsTot) / Ninsert;
 360     }
 361     // Add the free volume fraction to the data set in column 0
 362     dh.setPoint(0, frac);
 363     // Add the total volume to the data set in column 1
 364     dh.setPoint(1, V);
 365
 366     // Magic
 367     dh.finishFrame();
 368 }
 369
 370
 371 void FreeVolume::finishAnalysis(int /* nframes */)
 372 {
 373     please_cite(stdout, "Bondi1964a");
 374     please_cite(stdout, "Lourenco2013a");
 375 }
 376
 377 void FreeVolume::writeOutput()
 378 {
 379     // Final results come from statistics module in analysis framework
 380     double FVaver  = adata_->average(0, 0);
 381     double FVerror = adata_->standardDeviation(0, 0);
 382     printf("Free volume %.2f +/- %.2f %%\n", FVaver, FVerror);
 383
 384     double Vaver  = adata_->average(0, 1);
 385     double Verror = adata_->standardDeviation(0, 1);
 386     printf("Total volume %.2f +/- %.2f nm^3\n", Vaver, Verror);
 387
 388     printf("Number of molecules %d total mass %.2f Dalton\n", nmol_, mtot_);
 389     double RhoAver  = mtot_ / (Vaver * 1e-24 * AVOGADRO);
 390     double RhoError = gmx::square(RhoAver / Vaver) * Verror;
 391     printf("Average molar mass: %.2f Dalton\n", mtot_ / nmol_);
 392
 393     double VmAver  = Vaver / nmol_;
 394     double VmError = Verror / nmol_;
 395     printf("Density rho: %.2f +/- %.2f nm^3\n", RhoAver, RhoError);
 396     printf("Molecular volume Vm assuming homogeneity: %.4f +/- %.4f nm^3\n", VmAver, VmError);
 397
 398     double VvdWaver  = (1 - FVaver / 100) * VmAver;
 399     double VvdWerror = 0;
 400     printf("Molecular van der Waals volume assuming homogeneity:  %.4f +/- %.4f nm^3\n", VvdWaver, VvdWerror);
 401
 402     double FFVaver  = 1 - 1.3 * ((100 - FVaver) / 100);
 403     double FFVerror = (FVerror / FVaver) * FFVaver;
 404     printf("Fractional free volume %.3f +/- %.3f\n", FFVaver, FFVerror);
 405 }
 406
 407 } // namespace
 408
 409 const char FreeVolumeInfo::name[]             = "freevolume";
 410 const char FreeVolumeInfo::shortDescription[] = "Calculate free volume";
 411
 412 TrajectoryAnalysisModulePointer FreeVolumeInfo::create()
 413 {
 414     return TrajectoryAnalysisModulePointer(new FreeVolume);
 415 }
 416
 417 } // namespace analysismodules
 418
 419 } // namespace gmx