/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2013, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
*
* GROMACS is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public License
* \author David van der Spoel <david.vanderspoel@icm.uu.se>
* \ingroup module_trajectoryanalysis
*/
+#include "gmxpre.h"
+
#include "freevolume.h"
-#include "gromacs/legacyheaders/pbc.h"
-#include "gromacs/legacyheaders/vec.h"
-#include "gromacs/legacyheaders/atomprop.h"
-#include "gromacs/legacyheaders/copyrite.h"
+#include <string>
#include "gromacs/analysisdata/analysisdata.h"
+#include "gromacs/analysisdata/modules/average.h"
#include "gromacs/analysisdata/modules/plot.h"
+#include "gromacs/fileio/trx.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/math/vec.h"
#include "gromacs/options/basicoptions.h"
#include "gromacs/options/filenameoption.h"
#include "gromacs/options/options.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/random/random.h"
+#include "gromacs/selection/nbsearch.h"
#include "gromacs/selection/selection.h"
#include "gromacs/selection/selectionoption.h"
+#include "gromacs/topology/atomprop.h"
+#include "gromacs/topology/topology.h"
#include "gromacs/trajectoryanalysis/analysissettings.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/exceptions.h"
-#include "gromacs/utility/stringutil.h"
namespace gmx
{
namespace analysismodules
{
-const char FreeVolume::name[] = "freevolume";
-const char FreeVolume::shortDescription[] =
- "Calculate free volume";
+namespace
+{
+
+/*! \brief
+ * Class used to compute free volume in a simulations box.
+ *
+ * Inherits TrajectoryAnalysisModule and all functions from there.
+ * Does not implement any new functionality.
+ *
+ * \ingroup module_trajectoryanalysis
+ */
+class FreeVolume : public TrajectoryAnalysisModule
+{
+ public:
+ //! Constructor
+ FreeVolume();
+
+ //! Destructor
+ virtual ~FreeVolume();
+
+ //! Set the options and setting
+ virtual void initOptions(Options *options,
+ TrajectoryAnalysisSettings *settings);
+
+ //! First routine called by the analysis framework
+ virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
+ const TopologyInformation &top);
+
+ //! Call for each frame of the trajectory
+ virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+ TrajectoryAnalysisModuleData *pdata);
+
+ //! Last routine called by the analysis framework
+ virtual void finishAnalysis(int nframes);
+
+ //! Routine to write output, that is additional over the built-in
+ virtual void writeOutput();
+
+ private:
+ std::string fnFreevol_;
+ Selection sel_;
+ AnalysisData data_;
+ AnalysisDataAverageModulePointer adata_;
+
+ int nmol_;
+ double mtot_;
+ double cutoff_;
+ double probeRadius_;
+ gmx_rng_t rng_;
+ int seed_, ninsert_;
+ AnalysisNeighborhood nb_;
+ //! The van der Waals radius per atom
+ std::vector<double> vdw_radius_;
+
+ // Copy and assign disallowed by base.
+};
// Constructor. Here it is important to initialize the pointer to
// subclasses that are elements of the main class. Here we have only
// one. The type of this depends on what kind of tool you need.
// Here we only have simple value/time kind of data.
FreeVolume::FreeVolume()
- : TrajectoryAnalysisModule(name, shortDescription),
+ : TrajectoryAnalysisModule(FreeVolumeInfo::name, FreeVolumeInfo::shortDescription),
adata_(new AnalysisDataAverageModule())
{
// We only compute two numbers per frame
- data_.setColumnCount(2);
+ data_.setColumnCount(0, 2);
// Tell the analysis framework that this component exists
registerAnalysisDataset(&data_, "freevolume");
rng_ = NULL;
- nbsearch_ = NULL;
nmol_ = 0;
mtot_ = 0;
cutoff_ = 0;
{
gmx_rng_destroy(rng_);
}
- if (NULL != nbsearch_)
- {
- gmx_ana_nbsearch_free(nbsearch_);
- }
}
TrajectoryAnalysisSettings *settings)
{
static const char *const desc[] = {
- "g_freevol can calculate the free volume in a box as",
+ "[THISMODULE] calculates the free volume in a box as",
"a function of time. The free volume is",
"plotted as a fraction of the total volume.",
"The program tries to insert a probe with a given radius,",
};
// Add the descriptive text (program help text) to the options
- options->setDescription(concatenateStrings(desc));
+ options->setDescription(desc);
// Add option for optional output file
options->addOption(FileNameOption("o").filetype(eftPlot).outputFile()
.description("Computed free volume"));
// Add option for selecting a subset of atoms
- options->addOption(SelectionOption("select").required().valueCount(1)
- .store(&sel_)
+ options->addOption(SelectionOption("select")
+ .store(&sel_).defaultSelectionText("all")
.onlyAtoms());
// Add option for the probe radius and initialize it
rng_ = gmx_rng_init(seed_);
// Initiate the neighborsearching code
- nbsearch_ = gmx_ana_nbsearch_create(cutoff_, atoms->nr);
+ nb_.setCutoff(cutoff_);
}
void
int Ninsert = static_cast<int>(ninsert_*V);
// Use neighborsearching tools!
- gmx_ana_nbsearch_pos_init(nbsearch_, pbc, sel.positions());
+ AnalysisNeighborhoodSearch nbsearch = nb_.initSearch(pbc, sel);
// Then loop over insertions
int NinsTot = 0;
mvmul(fr.box, rand, ins);
// Find the first reference position within the cutoff.
- bool bOverlap = false;
- int jp = NOTSET;
- if (gmx_ana_nbsearch_first_within(nbsearch_, ins, &jp))
+ bool bOverlap = false;
+ AnalysisNeighborhoodPair pair;
+ AnalysisNeighborhoodPairSearch pairSearch = nbsearch.startPairSearch(ins);
+ while (!bOverlap && pairSearch.findNextPair(&pair))
{
- do
- {
- // Compute distance vector to first atom in the neighborlist
- pbc_dx(pbc, ins, sel.position(jp).x(), dx);
+ int jp = pair.refIndex();
+ // Compute distance vector to first atom in the neighborlist
+ pbc_dx(pbc, ins, sel.position(jp).x(), dx);
- // See whether the distance is smaller than allowed
- bOverlap = (norm(dx) <
- probeRadius_+vdw_radius_[sel.position(jp).refId()]);
+ // See whether the distance is smaller than allowed
+ bOverlap = (norm(dx) <
+ probeRadius_+vdw_radius_[sel.position(jp).refId()]);
- // Finds the next reference position within the cutoff.
- }
- while (!bOverlap &&
- (gmx_ana_nbsearch_next_within(nbsearch_, &jp)));
}
if (!bOverlap)
void
-FreeVolume::finishAnalysis(int nframes)
+FreeVolume::finishAnalysis(int /* nframes */)
{
please_cite(stdout, "Bondi1964a");
+ please_cite(stdout, "Lourenco2013a");
}
void
FreeVolume::writeOutput()
{
// Final results come from statistics module in analysis framework
- double FVaver = adata_->average(0);
- double FVerror = adata_->stddev(0);
+ double FVaver = adata_->average(0, 0);
+ double FVerror = adata_->standardDeviation(0, 0);
printf("Free volume %.2f +/- %.2f %%\n", FVaver, FVerror);
- double Vaver = adata_->average(1);
- double Verror = adata_->stddev(1);
+ double Vaver = adata_->average(0, 1);
+ double Verror = adata_->standardDeviation(0, 1);
printf("Total volume %.2f +/- %.2f nm^3\n", Vaver, Verror);
printf("Number of molecules %d total mass %.2f Dalton\n", nmol_, mtot_);
- double RhoAver = mtot_ / (adata_->average(1) * 1e-24 * AVOGADRO);
+ double RhoAver = mtot_ / (Vaver * 1e-24 * AVOGADRO);
double RhoError = sqr(RhoAver / Vaver)*Verror;
printf("Average molar mass: %.2f Dalton\n", mtot_/nmol_);
printf("Fractional free volume %.3f +/- %.3f\n", FFVaver, FFVerror);
}
+} // namespace
+
+const char FreeVolumeInfo::name[] = "freevolume";
+const char FreeVolumeInfo::shortDescription[] =
+ "Calculate free volume";
+
+TrajectoryAnalysisModulePointer FreeVolumeInfo::create()
+{
+ return TrajectoryAnalysisModulePointer(new FreeVolume);
+}
+
} // namespace analysismodules
} // namespace gmx
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff