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,2021, 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(OptionFileType::Plot)
 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 = 0;
 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,
 269                         "Could not determine VDW radius for %s-%s. Set to zero.\n",
 270                         *(atoms->resinfo[resnr].name),
 271                         *(atoms->atomname[i]));
 272             }
 273             vdw_radius_.push_back(0.0);
 274         }
 275     }
 276
 277     // Increase cutoff by proberadius to make sure we do not miss
 278     // anything
 279     cutoff_ += probeRadius_;
 280
 281     if (nnovdw >= maxnovdw)
 282     {
 283         fprintf(stderr, "Could not determine VDW radius for %d particles. These were set to zero.\n", nnovdw);
 284     }
 285
 286     if (seed_ == 0)
 287     {
 288         seed_ = static_cast<int>(gmx::makeRandomSeed());
 289     }
 290
 291     // Print parameters to output. Maybe should make dependent on
 292     // verbosity flag?
 293     printf("cutoff       = %g nm\n", cutoff_);
 294     printf("probe_radius = %g nm\n", probeRadius_);
 295     printf("seed         = %d\n", seed_);
 296     printf("ninsert      = %d probes per nm^3\n", ninsert_);
 297
 298     // Initiate the random number generator
 299     rng_.seed(seed_);
 300
 301     // Initiate the neighborsearching code
 302     nb_.setCutoff(cutoff_);
 303 }
 304
 305 void FreeVolume::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
 306 {
 307     AnalysisDataHandle                 dh  = pdata->dataHandle(data_);
 308     const Selection&                   sel = pdata->parallelSelection(sel_);
 309     gmx::UniformRealDistribution<real> dist;
 310
 311     GMX_RELEASE_ASSERT(nullptr != pbc, "You have no periodic boundary conditions");
 312
 313     // Analysis framework magic
 314     dh.startFrame(frnr, fr.time);
 315
 316     // Compute volume and number of insertions to perform
 317     real V       = det(fr.box);
 318     int  Ninsert = static_cast<int>(ninsert_ * V);
 319
 320     // Use neighborsearching tools!
 321     AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, sel);
 322
 323     // Then loop over insertions
 324     int NinsTot = 0;
 325     for (int i = 0; (i < Ninsert); i++)
 326     {
 327         rvec rand, ins, dx;
 328
 329         for (int m = 0; (m < DIM); m++)
 330         {
 331             // Generate random number between 0 and 1
 332             rand[m] = dist(rng_);
 333         }
 334         // Generate random 3D position within the box
 335         mvmul(fr.box, rand, ins);
 336
 337         // Find the first reference position within the cutoff.
 338         bool                           bOverlap = false;
 339         AnalysisNeighborhoodPair       pair;
 340         AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(ins);
 341         while (!bOverlap && pairSearch.findNextPair(&pair))
 342         {
 343             int jp = pair.refIndex();
 344             // Compute distance vector to first atom in the neighborlist
 345             pbc_dx(pbc, ins, sel.position(jp).x(), dx);
 346
 347             // See whether the distance is smaller than allowed
 348             bOverlap = (norm(dx) < probeRadius_ + vdw_radius_[sel.position(jp).refId()]);
 349         }
 350
 351         if (!bOverlap)
 352         {
 353             // We found some free volume!
 354             NinsTot++;
 355         }
 356     }
 357     // Compute total free volume for this frame
 358     double frac = 0;
 359     if (Ninsert > 0)
 360     {
 361         frac = (100.0 * NinsTot) / Ninsert;
 362     }
 363     // Add the free volume fraction to the data set in column 0
 364     dh.setPoint(0, frac);
 365     // Add the total volume to the data set in column 1
 366     dh.setPoint(1, V);
 367
 368     // Magic
 369     dh.finishFrame();
 370 }
 371
 372
 373 void FreeVolume::finishAnalysis(int /* nframes */)
 374 {
 375     please_cite(stdout, "Bondi1964a");
 376     please_cite(stdout, "Lourenco2013a");
 377 }
 378
 379 void FreeVolume::writeOutput()
 380 {
 381     // Final results come from statistics module in analysis framework
 382     double FVaver  = adata_->average(0, 0);
 383     double FVerror = adata_->standardDeviation(0, 0);
 384     printf("Free volume %.2f +/- %.2f %%\n", FVaver, FVerror);
 385
 386     double Vaver  = adata_->average(0, 1);
 387     double Verror = adata_->standardDeviation(0, 1);
 388     printf("Total volume %.2f +/- %.2f nm^3\n", Vaver, Verror);
 389
 390     printf("Number of molecules %d total mass %.2f Dalton\n", nmol_, mtot_);
 391     double RhoAver  = mtot_ / (Vaver * 1e-24 * gmx::c_avogadro);
 392     double RhoError = gmx::square(RhoAver / Vaver) * Verror;
 393     printf("Average molar mass: %.2f Dalton\n", mtot_ / nmol_);
 394
 395     double VmAver  = Vaver / nmol_;
 396     double VmError = Verror / nmol_;
 397     printf("Density rho: %.2f +/- %.2f nm^3\n", RhoAver, RhoError);
 398     printf("Molecular volume Vm assuming homogeneity: %.4f +/- %.4f nm^3\n", VmAver, VmError);
 399
 400     double VvdWaver  = (1 - FVaver / 100) * VmAver;
 401     double VvdWerror = 0;
 402     printf("Molecular van der Waals volume assuming homogeneity:  %.4f +/- %.4f nm^3\n", VvdWaver, VvdWerror);
 403
 404     double FFVaver  = 1 - 1.3 * ((100 - FVaver) / 100);
 405     double FFVerror = (FVerror / FVaver) * FFVaver;
 406     printf("Fractional free volume %.3f +/- %.3f\n", FFVaver, FFVerror);
 407 }
 408
 409 } // namespace
 410
 411 const char FreeVolumeInfo::name[]             = "freevolume";
 412 const char FreeVolumeInfo::shortDescription[] = "Calculate free volume";
 413
 414 TrajectoryAnalysisModulePointer FreeVolumeInfo::create()
 415 {
 416     return TrajectoryAnalysisModulePointer(new FreeVolume);
 417 }
 418
 419 } // namespace analysismodules
 420
 421 } // namespace gmx