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