src/gromacs/trajectoryanalysis/modules/freevolume.cpp

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