struct t_conect
{
//! Index of the second nearest neighbor dot.
- int aa;
+ int aa;
//! Index of the nearest neighbor dot.
- int ab;
+ int ab;
//! Squared distance to `aa`.
- real d2a;
+ real d2a;
//! Squared distance to `ab`.
- real d2b;
+ real d2b;
};
/*! \brief
* neighbors. The function is copied verbatim from the old gmx_sas.c
* implementation.
*/
-void do_conect(const char *fn, int n, rvec x[])
+void do_conect(const char* fn, int n, rvec x[])
{
- FILE *fp;
+ FILE* fp;
int i, j;
- t_conect *c;
+ t_conect* c;
rvec dx;
real d2;
for (i = 0; (i < n); i++)
{
- for (j = i+1; (j < n); j++)
+ for (j = i + 1; (j < n); j++)
{
rvec_sub(x[i], x[j], dx);
d2 = iprod(dx, dx);
{
if ((c[i].aa == -1) || (c[i].ab == -1))
{
- fprintf(stderr, "Warning dot %d has no connections\n", i+1);
+ fprintf(stderr, "Warning dot %d has no connections\n", i + 1);
}
- fprintf(fp, "CONECT%5d%5d%5d\n", i+1, c[i].aa+1, c[i].ab+1);
+ fprintf(fp, "CONECT%5d%5d%5d\n", i + 1, c[i].aa + 1, c[i].ab + 1);
}
gmx_ffclose(fp);
sfree(c);
/*! \brief
* Plots the surface into a PDB file, optionally including the original atoms.
*/
-void connolly_plot(const char *fn, int ndots, const real dots[], rvec x[], t_atoms *atoms,
- t_symtab *symtab, int ePBC, const matrix box, gmx_bool bIncludeSolute)
+void connolly_plot(const char* fn,
+ int ndots,
+ const real dots[],
+ rvec x[],
+ t_atoms* atoms,
+ t_symtab* symtab,
+ int ePBC,
+ const matrix box,
+ gmx_bool bIncludeSolute)
{
- const char *const atomnm = "DOT";
- const char *const resnm = "DOT";
- const char *const title = "Connolly Dot Surface Generated by gmx sasa";
+ const char* const atomnm = "DOT";
+ const char* const resnm = "DOT";
+ const char* const title = "Connolly Dot Surface Generated by gmx sasa";
- int i, i0, r0, ii0, k;
- rvec *xnew;
- t_atoms aaa;
+ int i, i0, r0, ii0, k;
+ rvec* xnew;
+ t_atoms aaa;
if (bIncludeSolute)
{
i0 = atoms->nr;
r0 = atoms->nres;
- srenew(atoms->atom, atoms->nr+ndots);
- memset(&atoms->atom[i0], 0, sizeof(*atoms->atom)*ndots);
- srenew(atoms->atomname, atoms->nr+ndots);
- srenew(atoms->resinfo, r0+1);
+ srenew(atoms->atom, atoms->nr + ndots);
+ memset(&atoms->atom[i0], 0, sizeof(*atoms->atom) * ndots);
+ srenew(atoms->atomname, atoms->nr + ndots);
+ srenew(atoms->resinfo, r0 + 1);
atoms->atom[i0].resind = r0;
- t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0+1, ' ', 0, ' ');
+ t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0 + 1, ' ', 0, ' ');
if (atoms->pdbinfo != nullptr)
{
- srenew(atoms->pdbinfo, atoms->nr+ndots);
+ srenew(atoms->pdbinfo, atoms->nr + ndots);
}
- snew(xnew, atoms->nr+ndots);
+ snew(xnew, atoms->nr + ndots);
for (i = 0; (i < atoms->nr); i++)
{
copy_rvec(x[i], xnew[i]);
}
for (i = k = 0; (i < ndots); i++)
{
- ii0 = i0+i;
- atoms->atomname[ii0] = put_symtab(symtab, atomnm);
- atoms->atom[ii0].resind = r0;
- xnew[ii0][XX] = dots[k++];
- xnew[ii0][YY] = dots[k++];
- xnew[ii0][ZZ] = dots[k++];
+ ii0 = i0 + i;
+ atoms->atomname[ii0] = put_symtab(symtab, atomnm);
+ atoms->atom[ii0].resind = r0;
+ xnew[ii0][XX] = dots[k++];
+ xnew[ii0][YY] = dots[k++];
+ xnew[ii0][ZZ] = dots[k++];
if (atoms->pdbinfo != nullptr)
{
atoms->pdbinfo[ii0].type = epdbATOM;
atoms->pdbinfo[ii0].occup = 0.0;
}
}
- atoms->nr = i0+ndots;
- atoms->nres = r0+1;
- write_sto_conf(fn, title, atoms, xnew, nullptr, ePBC, const_cast<rvec *>(box));
+ atoms->nr = i0 + ndots;
+ atoms->nres = r0 + 1;
+ write_sto_conf(fn, title, atoms, xnew, nullptr, ePBC, const_cast<rvec*>(box));
atoms->nres = r0;
atoms->nr = i0;
}
aaa.pdbinfo[ii0].occup = 0.0;
}
aaa.nr = ndots;
- write_sto_conf(fn, title, &aaa, xnew, nullptr, ePBC, const_cast<rvec *>(box));
+ write_sto_conf(fn, title, &aaa, xnew, nullptr, ePBC, const_cast<rvec*>(box));
do_conect(fn, ndots, xnew);
done_atom(&aaa);
}
*/
class Sasa : public TrajectoryAnalysisModule
{
- public:
- Sasa();
-
- void initOptions(IOptionsContainer *options,
- TrajectoryAnalysisSettings *settings) override;
- void initAnalysis(const TrajectoryAnalysisSettings &settings,
- const TopologyInformation &top) override;
-
- TrajectoryAnalysisModuleDataPointer startFrames(
- const AnalysisDataParallelOptions &opt,
- const SelectionCollection &selections) override;
- void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
- TrajectoryAnalysisModuleData *pdata) override;
-
- void finishAnalysis(int nframes) override;
- void writeOutput() override;
-
- private:
- /*! \brief
- * Surface areas as a function of time.
- *
- * First column is for the calculation group, and the rest for the
- * output groups. This data is always produced.
- */
- AnalysisData area_;
- /*! \brief
- * Per-atom surface areas as a function of time.
- *
- * Contains one data set for each column in `area_`.
- * Each column corresponds to a selection position in `surfaceSel_`.
- * This data is only produced if atom or residue areas have been
- * requested.
- */
- AnalysisData atomArea_;
- /*! \brief
- * Per-residue surface areas as a function of time.
- *
- * Contains one data set for each column in `area_`.
- * Each column corresponds to a distinct residue `surfaceSel_`.
- * For example, if `surfaceSel_` selects residues 2, 5, and 7, there
- * will be three columns here.
- * This data is only produced if atom or residue areas have been
- * requested.
- */
- AnalysisData residueArea_;
- /*! \brief
- * Free energy estimates as a function of time.
- *
- * Column layout is the same as for `area_`.
- * This data is only produced if the output is requested.
- */
- AnalysisData dgSolv_;
- /*! \brief
- * Total volume and density of the calculation group as a function of
- * time.
- *
- * The first column is the volume and the second column is the density.
- * This data is only produced if the output is requested.
- */
- AnalysisData volume_;
-
- /*! \brief
- * The selection to calculate the surface for.
- *
- * Selection::originalId() and Selection::mappedId() store the mapping
- * from the positions to the columns of `residueArea_`.
- * The selection is computed with SelectionOption::dynamicMask(), i.e.,
- * even in the presence of a dynamic selection, the number of returned
- * positions is fixed, and SelectionPosition::selected() is used.
- */
- Selection surfaceSel_;
- /*! \brief
- * List of optional additional output groups.
- *
- * Each of these must be a subset of the `surfaceSel_`.
- * Selection::originalId() and Selection::mappedId() store the mapping
- * from the positions to the corresponsing positions in `surfaceSel_`.
- */
- SelectionList outputSel_;
-
- std::string fnArea_;
- std::string fnAtomArea_;
- std::string fnResidueArea_;
- std::string fnDGSolv_;
- std::string fnVolume_;
- std::string fnConnolly_;
-
- double solsize_;
- int ndots_;
- //double minarea_;
- double dgsDefault_;
- bool bIncludeSolute_;
-
- //! Global topology corresponding to the input.
- gmx_mtop_t *mtop_;
- //! Per-atom data corresponding to the input.
- AtomsDataPtr atoms_;
- //! Combined VdW and probe radii for each atom in the calculation group.
- std::vector<real> radii_;
- /*! \brief
- * Solvation free energy coefficients for each atom in the calculation
- * group.
- *
- * Empty if the free energy output has not been requested.
- */
- std::vector<real> dgsFactor_;
- //! Calculation algorithm.
- SurfaceAreaCalculator calculator_;
-
- // Copy and assign disallowed by base.
+public:
+ Sasa();
+
+ void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
+ void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
+
+ TrajectoryAnalysisModuleDataPointer startFrames(const AnalysisDataParallelOptions& opt,
+ const SelectionCollection& selections) override;
+ void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
+
+ void finishAnalysis(int nframes) override;
+ void writeOutput() override;
+
+private:
+ /*! \brief
+ * Surface areas as a function of time.
+ *
+ * First column is for the calculation group, and the rest for the
+ * output groups. This data is always produced.
+ */
+ AnalysisData area_;
+ /*! \brief
+ * Per-atom surface areas as a function of time.
+ *
+ * Contains one data set for each column in `area_`.
+ * Each column corresponds to a selection position in `surfaceSel_`.
+ * This data is only produced if atom or residue areas have been
+ * requested.
+ */
+ AnalysisData atomArea_;
+ /*! \brief
+ * Per-residue surface areas as a function of time.
+ *
+ * Contains one data set for each column in `area_`.
+ * Each column corresponds to a distinct residue `surfaceSel_`.
+ * For example, if `surfaceSel_` selects residues 2, 5, and 7, there
+ * will be three columns here.
+ * This data is only produced if atom or residue areas have been
+ * requested.
+ */
+ AnalysisData residueArea_;
+ /*! \brief
+ * Free energy estimates as a function of time.
+ *
+ * Column layout is the same as for `area_`.
+ * This data is only produced if the output is requested.
+ */
+ AnalysisData dgSolv_;
+ /*! \brief
+ * Total volume and density of the calculation group as a function of
+ * time.
+ *
+ * The first column is the volume and the second column is the density.
+ * This data is only produced if the output is requested.
+ */
+ AnalysisData volume_;
+
+ /*! \brief
+ * The selection to calculate the surface for.
+ *
+ * Selection::originalId() and Selection::mappedId() store the mapping
+ * from the positions to the columns of `residueArea_`.
+ * The selection is computed with SelectionOption::dynamicMask(), i.e.,
+ * even in the presence of a dynamic selection, the number of returned
+ * positions is fixed, and SelectionPosition::selected() is used.
+ */
+ Selection surfaceSel_;
+ /*! \brief
+ * List of optional additional output groups.
+ *
+ * Each of these must be a subset of the `surfaceSel_`.
+ * Selection::originalId() and Selection::mappedId() store the mapping
+ * from the positions to the corresponsing positions in `surfaceSel_`.
+ */
+ SelectionList outputSel_;
+
+ std::string fnArea_;
+ std::string fnAtomArea_;
+ std::string fnResidueArea_;
+ std::string fnDGSolv_;
+ std::string fnVolume_;
+ std::string fnConnolly_;
+
+ double solsize_;
+ int ndots_;
+ // double minarea_;
+ double dgsDefault_;
+ bool bIncludeSolute_;
+
+ //! Global topology corresponding to the input.
+ gmx_mtop_t* mtop_;
+ //! Per-atom data corresponding to the input.
+ AtomsDataPtr atoms_;
+ //! Combined VdW and probe radii for each atom in the calculation group.
+ std::vector<real> radii_;
+ /*! \brief
+ * Solvation free energy coefficients for each atom in the calculation
+ * group.
+ *
+ * Empty if the free energy output has not been requested.
+ */
+ std::vector<real> dgsFactor_;
+ //! Calculation algorithm.
+ SurfaceAreaCalculator calculator_;
+
+ // Copy and assign disallowed by base.
};
-Sasa::Sasa()
- : solsize_(0.14), ndots_(24), dgsDefault_(0), bIncludeSolute_(true),
- mtop_(nullptr), atoms_(nullptr)
+Sasa::Sasa() :
+ solsize_(0.14),
+ ndots_(24),
+ dgsDefault_(0),
+ bIncludeSolute_(true),
+ mtop_(nullptr),
+ atoms_(nullptr)
{
- //minarea_ = 0.5;
+ // minarea_ = 0.5;
registerAnalysisDataset(&area_, "area");
registerAnalysisDataset(&atomArea_, "atomarea");
registerAnalysisDataset(&residueArea_, "resarea");
registerAnalysisDataset(&volume_, "volume");
}
-void
-Sasa::initOptions(IOptionsContainer *options, TrajectoryAnalysisSettings *settings)
+void Sasa::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
{
- static const char *const desc[] = {
+ static const char* const desc[] = {
"[THISMODULE] computes solvent accessible surface areas.",
"See Eisenhaber F, Lijnzaad P, Argos P, Sander C, & Scharf M",
"(1995) J. Comput. Chem. 16, 273-284 for the algorithm used.",
"from the full surface.[PAR]",
"The average and standard deviation of the area over the trajectory",
- "can be calculated per residue and atom (options [TT]-or[tt] and",
- "[TT]-oa[tt]).[PAR]",
+ "can be calculated per residue and atom (options [TT]-or[tt] and", "[TT]-oa[tt]).[PAR]",
//"In combination with the latter option an [REF].itp[ref] file can be",
//"generated (option [TT]-i[tt])",
//"which can be used to restrain surface atoms.[PAR]",
settings->setHelpText(desc);
- options->addOption(FileNameOption("o").filetype(eftPlot).outputFile().required()
- .store(&fnArea_).defaultBasename("area")
- .description("Total area as a function of time"));
- options->addOption(FileNameOption("odg").filetype(eftPlot).outputFile()
- .store(&fnDGSolv_).defaultBasename("dgsolv")
- .description("Estimated solvation free energy as a function of time"));
- options->addOption(FileNameOption("or").filetype(eftPlot).outputFile()
- .store(&fnResidueArea_).defaultBasename("resarea")
- .description("Average area per residue"));
- options->addOption(FileNameOption("oa").filetype(eftPlot).outputFile()
- .store(&fnAtomArea_).defaultBasename("atomarea")
- .description("Average area per atom"));
- options->addOption(FileNameOption("tv").filetype(eftPlot).outputFile()
- .store(&fnVolume_).defaultBasename("volume")
- .description("Total volume and density as a function of time"));
- options->addOption(FileNameOption("q").filetype(eftPDB).outputFile()
- .store(&fnConnolly_).defaultBasename("connolly")
- .description("PDB file for Connolly surface"));
- //options->addOption(FileNameOption("i").filetype(eftITP).outputFile()
+ options->addOption(FileNameOption("o")
+ .filetype(eftPlot)
+ .outputFile()
+ .required()
+ .store(&fnArea_)
+ .defaultBasename("area")
+ .description("Total area as a function of time"));
+ options->addOption(
+ FileNameOption("odg")
+ .filetype(eftPlot)
+ .outputFile()
+ .store(&fnDGSolv_)
+ .defaultBasename("dgsolv")
+ .description("Estimated solvation free energy as a function of time"));
+ options->addOption(FileNameOption("or")
+ .filetype(eftPlot)
+ .outputFile()
+ .store(&fnResidueArea_)
+ .defaultBasename("resarea")
+ .description("Average area per residue"));
+ options->addOption(FileNameOption("oa")
+ .filetype(eftPlot)
+ .outputFile()
+ .store(&fnAtomArea_)
+ .defaultBasename("atomarea")
+ .description("Average area per atom"));
+ options->addOption(FileNameOption("tv")
+ .filetype(eftPlot)
+ .outputFile()
+ .store(&fnVolume_)
+ .defaultBasename("volume")
+ .description("Total volume and density as a function of time"));
+ options->addOption(FileNameOption("q")
+ .filetype(eftPDB)
+ .outputFile()
+ .store(&fnConnolly_)
+ .defaultBasename("connolly")
+ .description("PDB file for Connolly surface"));
+ // options->addOption(FileNameOption("i").filetype(eftITP).outputFile()
// .store(&fnRestraints_).defaultBasename("surfat")
// .description("Topology file for position restraings on surface atoms"));
- options->addOption(DoubleOption("probe").store(&solsize_)
- .description("Radius of the solvent probe (nm)"));
- options->addOption(IntegerOption("ndots").store(&ndots_)
- .description("Number of dots per sphere, more dots means more accuracy"));
- //options->addOption(DoubleOption("minarea").store(&minarea_)
+ options->addOption(
+ DoubleOption("probe").store(&solsize_).description("Radius of the solvent probe (nm)"));
+ options->addOption(IntegerOption("ndots").store(&ndots_).description(
+ "Number of dots per sphere, more dots means more accuracy"));
+ // options->addOption(DoubleOption("minarea").store(&minarea_)
// .description("The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file (see help)"));
- options->addOption(BooleanOption("prot").store(&bIncludeSolute_)
- .description("Output the protein to the Connolly [REF].pdb[ref] file too"));
- options->addOption(DoubleOption("dgs").store(&dgsDefault_)
- .description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
+ options->addOption(
+ BooleanOption("prot").store(&bIncludeSolute_).description("Output the protein to the Connolly [REF].pdb[ref] file too"));
+ options->addOption(
+ DoubleOption("dgs").store(&dgsDefault_).description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
// Selections must select atoms for the VdW radii lookup to work.
// The calculation group uses dynamicMask() so that the coordinates
// match a static array of VdW radii.
- options->addOption(SelectionOption("surface").store(&surfaceSel_)
- .required().onlySortedAtoms().dynamicMask()
- .description("Surface calculation selection"));
- options->addOption(SelectionOption("output").storeVector(&outputSel_)
- .onlySortedAtoms().multiValue()
- .description("Output selection(s)"));
+ options->addOption(SelectionOption("surface")
+ .store(&surfaceSel_)
+ .required()
+ .onlySortedAtoms()
+ .dynamicMask()
+ .description("Surface calculation selection"));
+ options->addOption(
+ SelectionOption("output").storeVector(&outputSel_).onlySortedAtoms().multiValue().description("Output selection(s)"));
// Atom names etc. are required for the VdW radii lookup.
settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
}
-void
-Sasa::initAnalysis(const TrajectoryAnalysisSettings &settings,
- const TopologyInformation &top)
+void Sasa::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
{
mtop_ = top.mtop();
atoms_ = top.copyAtoms();
- //bITP = opt2bSet("-i", nfile, fnm);
- const bool bResAt =
- !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
+ // bITP = opt2bSet("-i", nfile, fnm);
+ const bool bResAt = !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
const bool bDGsol = !fnDGSolv_.empty();
if (solsize_ < 0)
}
please_cite(stderr, "Eisenhaber95");
- //if ((top.ePBC() != epbcXYZ) || (TRICLINIC(fr.box)))
+ // if ((top.ePBC() != epbcXYZ) || (TRICLINIC(fr.box)))
//{
// fprintf(stderr, "\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
// "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
{
if (!top.hasFullTopology())
{
- GMX_THROW(InconsistentInputError("Cannot compute Delta G of solvation without a tpr file"));
+ GMX_THROW(InconsistentInputError(
+ "Cannot compute Delta G of solvation without a tpr file"));
}
else
{
if (strcmp(*(atoms_->atomtype[0]), "?") == 0)
{
- GMX_THROW(InconsistentInputError("Your input tpr file is too old (does not contain atom types). Cannot not compute Delta G of solvation"));
+ GMX_THROW(InconsistentInputError(
+ "Your input tpr file is too old (does not contain atom types). Cannot not "
+ "compute Delta G of solvation"));
}
else
{
// TODO: Not exception-safe, but nice solution would be to have a C++
// atom properties class...
- AtomProperties aps;
+ AtomProperties aps;
ArrayRef<const int> atomIndices = surfaceSel_.atomIndices();
int ndefault = 0;
const int ii = atomIndices[i];
const int resind = atoms_->atom[ii].resind;
real radius;
- if (!aps.setAtomProperty(epropVDW,
- *(atoms_->resinfo[resind].name),
- *(atoms_->atomname[ii]), &radius))
+ if (!aps.setAtomProperty(epropVDW, *(atoms_->resinfo[resind].name), *(atoms_->atomname[ii]), &radius))
{
ndefault++;
}
if (bDGsol)
{
real dgsFactor;
- if (!aps.setAtomProperty(epropDGsol,
- *(atoms_->resinfo[resind].name),
+ if (!aps.setAtomProperty(epropDGsol, *(atoms_->resinfo[resind].name),
*(atoms_->atomtype[ii]), &dgsFactor))
{
dgsFactor = dgsDefault_;
}
if (j == surfaceSel_.posCount() || outputIndices[i] != atomIndices[j])
{
- const std::string message
- = formatString("Output selection '%s' is not a subset of "
- "the surface selection (atom %d is the first "
- "atom not in the surface selection)",
- outputSel_[g].name(), outputIndices[i] + 1);
+ const std::string message = formatString(
+ "Output selection '%s' is not a subset of "
+ "the surface selection (atom %d is the first "
+ "atom not in the surface selection)",
+ outputSel_[g].name(), outputIndices[i] + 1);
GMX_THROW(InconsistentInputError(message));
}
outputSel_[g].setOriginalId(i, j);
area_.setColumnCount(0, 1 + outputSel_.size());
{
- AnalysisDataPlotModulePointer plotm(
- new AnalysisDataPlotModule(settings.plotSettings()));
+ AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
plotm->setFileName(fnArea_);
plotm->setTitle("Solvent Accessible Surface");
plotm->setXAxisIsTime();
atomArea_.addModule(avem);
if (!fnAtomArea_.empty())
{
- AnalysisDataPlotModulePointer plotm(
- new AnalysisDataPlotModule(settings.plotSettings()));
+ AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
plotm->setFileName(fnAtomArea_);
plotm->setTitle("Area per atom over the trajectory");
plotm->setXLabel("Atom");
}
{
AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
- int nextRow = 0;
+ int nextRow = 0;
for (int i = 0; i < surfaceSel_.posCount(); ++i)
{
const int residueGroup = surfaceSel_.position(i).mappedId();
residueArea_.addModule(avem);
if (!fnResidueArea_.empty())
{
- AnalysisDataPlotModulePointer plotm(
- new AnalysisDataPlotModule(settings.plotSettings()));
+ AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
plotm->setFileName(fnResidueArea_);
plotm->setTitle("Area per residue over the trajectory");
plotm->setXLabel("Residue");
if (!fnDGSolv_.empty())
{
dgSolv_.setColumnCount(0, 1 + outputSel_.size());
- AnalysisDataPlotModulePointer plotm(
- new AnalysisDataPlotModule(settings.plotSettings()));
+ AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
plotm->setFileName(fnDGSolv_);
plotm->setTitle("Free Energy of Solvation");
plotm->setXAxisIsTime();
if (!fnVolume_.empty())
{
volume_.setColumnCount(0, 2);
- AnalysisDataPlotModulePointer plotm(
- new AnalysisDataPlotModule(settings.plotSettings()));
+ AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
plotm->setFileName(fnVolume_);
plotm->setTitle("Volume and Density");
plotm->setXAxisIsTime();
*/
class SasaModuleData : public TrajectoryAnalysisModuleData
{
- public:
- /*! \brief
- * Reserves memory for the frame-local data.
- *
- * `residueCount` will be zero if per-residue data is not being
- * calculated.
- */
- SasaModuleData(TrajectoryAnalysisModule *module,
- const AnalysisDataParallelOptions &opt,
- const SelectionCollection &selections,
- int atomCount, int residueCount)
- : TrajectoryAnalysisModuleData(module, opt, selections)
+public:
+ /*! \brief
+ * Reserves memory for the frame-local data.
+ *
+ * `residueCount` will be zero if per-residue data is not being
+ * calculated.
+ */
+ SasaModuleData(TrajectoryAnalysisModule* module,
+ const AnalysisDataParallelOptions& opt,
+ const SelectionCollection& selections,
+ int atomCount,
+ int residueCount) :
+ TrajectoryAnalysisModuleData(module, opt, selections)
+ {
+ index_.reserve(atomCount);
+ // If the calculation group is not dynamic, pre-calculate
+ // the index, since it is not going to change.
+ for (int i = 0; i < atomCount; ++i)
{
- index_.reserve(atomCount);
- // If the calculation group is not dynamic, pre-calculate
- // the index, since it is not going to change.
- for (int i = 0; i < atomCount; ++i)
- {
- index_.push_back(i);
- }
- atomAreas_.resize(atomCount);
- res_a_.resize(residueCount);
+ index_.push_back(i);
}
+ atomAreas_.resize(atomCount);
+ res_a_.resize(residueCount);
+ }
- void finish() override { finishDataHandles(); }
-
- //! Indices of the calculation selection positions selected for the frame.
- std::vector<int> index_;
- /*! \brief
- * Atom areas for each calculation selection position for the frame.
- *
- * One entry for each position in the calculation group.
- * Values for atoms not selected are set to zero.
- */
- std::vector<real> atomAreas_;
- /*! \brief
- * Working array to accumulate areas for each residue.
- *
- * One entry for each distinct residue in the calculation group;
- * indices are not directly residue numbers or residue indices.
- *
- * This vector is empty if residue area calculations are not being
- * performed.
- */
- std::vector<real> res_a_;
+ void finish() override { finishDataHandles(); }
+
+ //! Indices of the calculation selection positions selected for the frame.
+ std::vector<int> index_;
+ /*! \brief
+ * Atom areas for each calculation selection position for the frame.
+ *
+ * One entry for each position in the calculation group.
+ * Values for atoms not selected are set to zero.
+ */
+ std::vector<real> atomAreas_;
+ /*! \brief
+ * Working array to accumulate areas for each residue.
+ *
+ * One entry for each distinct residue in the calculation group;
+ * indices are not directly residue numbers or residue indices.
+ *
+ * This vector is empty if residue area calculations are not being
+ * performed.
+ */
+ std::vector<real> res_a_;
};
-TrajectoryAnalysisModuleDataPointer Sasa::startFrames(
- const AnalysisDataParallelOptions &opt,
- const SelectionCollection &selections)
+TrajectoryAnalysisModuleDataPointer Sasa::startFrames(const AnalysisDataParallelOptions& opt,
+ const SelectionCollection& selections)
{
- return TrajectoryAnalysisModuleDataPointer(
- new SasaModuleData(this, opt, selections, surfaceSel_.posCount(),
- residueArea_.columnCount(0)));
+ return TrajectoryAnalysisModuleDataPointer(new SasaModuleData(
+ this, opt, selections, surfaceSel_.posCount(), residueArea_.columnCount(0)));
}
/*! \brief
*
* `atomAreaHandle` and `resAreaHandle` are not used if `resAreaWork` is empty.
*/
-void computeAreas(const Selection &surfaceSel, const Selection &sel,
- const std::vector<real> &atomAreas,
- const std::vector<real> &dgsFactor,
- real *totalAreaOut, real *dgsolvOut,
- AnalysisDataHandle atomAreaHandle,
- AnalysisDataHandle resAreaHandle,
- std::vector<real> *resAreaWork)
+void computeAreas(const Selection& surfaceSel,
+ const Selection& sel,
+ const std::vector<real>& atomAreas,
+ const std::vector<real>& dgsFactor,
+ real* totalAreaOut,
+ real* dgsolvOut,
+ AnalysisDataHandle atomAreaHandle,
+ AnalysisDataHandle resAreaHandle,
+ std::vector<real>* resAreaWork)
{
const bool bResAt = !resAreaWork->empty();
const bool bDGsolv = !dgsFactor.empty();
// Get the index of the atom in the calculation group.
// For the output groups, the mapping has been precalculated in
// initAnalysis().
- const int ii = (sel != surfaceSel ? sel.position(i).mappedId() : i);
+ const int ii = (sel != surfaceSel ? sel.position(i).mappedId() : i);
if (!surfaceSel.position(ii).selected())
{
// For the calculation group, skip unselected atoms.
{
continue;
}
- GMX_THROW(InconsistentInputError("Output selection is not a subset of the surface selection"));
+ GMX_THROW(InconsistentInputError(
+ "Output selection is not a subset of the surface selection"));
}
// Get the internal index of the matching residue.
// These have been precalculated in initAnalysis().
*dgsolvOut = dgsolv;
}
-void
-Sasa::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
- TrajectoryAnalysisModuleData *pdata)
+void Sasa::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
{
AnalysisDataHandle ah = pdata->dataHandle(area_);
AnalysisDataHandle dgh = pdata->dataHandle(dgSolv_);
AnalysisDataHandle aah = pdata->dataHandle(atomArea_);
AnalysisDataHandle rah = pdata->dataHandle(residueArea_);
AnalysisDataHandle vh = pdata->dataHandle(volume_);
- const Selection &surfaceSel = pdata->parallelSelection(surfaceSel_);
- const SelectionList &outputSel = pdata->parallelSelections(outputSel_);
- SasaModuleData &frameData = *static_cast<SasaModuleData *>(pdata);
+ const Selection& surfaceSel = pdata->parallelSelection(surfaceSel_);
+ const SelectionList& outputSel = pdata->parallelSelections(outputSel_);
+ SasaModuleData& frameData = *static_cast<SasaModuleData*>(pdata);
- const bool bResAt = !frameData.res_a_.empty();
- const bool bDGsol = !dgsFactor_.empty();
- const bool bConnolly = (frnr == 0 && !fnConnolly_.empty());
+ const bool bResAt = !frameData.res_a_.empty();
+ const bool bDGsol = !dgsFactor_.empty();
+ const bool bConnolly = (frnr == 0 && !fnConnolly_.empty());
// Update indices of selected atoms in the work array.
if (surfaceSel.isDynamic())
}
// Determine what needs to be calculated.
- int flag = 0;
+ int flag = 0;
if (bResAt || bDGsol || !outputSel.empty())
{
flag |= FLAG_ATOM_AREA;
real totarea, totvolume;
real *area = nullptr, *surfacedots = nullptr;
int nsurfacedots;
- calculator_.calculate(surfaceSel.coordinates().data(), pbc,
- frameData.index_.size(), frameData.index_.data(), flag,
- &totarea, &totvolume, &area,
- &surfacedots, &nsurfacedots);
+ calculator_.calculate(surfaceSel.coordinates().data(), pbc, frameData.index_.size(),
+ frameData.index_.data(), flag, &totarea, &totvolume, &area, &surfacedots,
+ &nsurfacedots);
// Unpack the atomwise areas into the frameData.atomAreas_ array for easier
// indexing in the case of dynamic surfaceSel.
if (area != nullptr)
{
if (surfaceSel.isDynamic())
{
- std::fill(frameData.atomAreas_.begin(), frameData.atomAreas_.end(),
- 0.0_real);
+ std::fill(frameData.atomAreas_.begin(), frameData.atomAreas_.end(), 0.0_real);
for (size_t i = 0; i < frameData.index_.size(); ++i)
{
frameData.atomAreas_[frameData.index_[i]] = area[i];
}
else
{
- std::copy(area, area + surfaceSel.posCount(),
- frameData.atomAreas_.begin());
+ std::copy(area, area + surfaceSel.posCount(), frameData.atomAreas_.begin());
}
sfree(area);
}
{
if (fr.natoms != mtop_->natoms)
{
- GMX_THROW(InconsistentInputError("Connolly plot (-q) is only supported for trajectories that contain all the atoms"));
+ GMX_THROW(
+ InconsistentInputError("Connolly plot (-q) is only supported for trajectories "
+ "that contain all the atoms"));
}
// This is somewhat nasty, as it modifies the atoms and symtab
// structures. But since it is only used in the first frame, and no
// one else uses the topology after initialization, it may just work
// even with future parallelization.
- connolly_plot(fnConnolly_.c_str(),
- nsurfacedots, surfacedots, fr.x, atoms_.get(),
+ connolly_plot(fnConnolly_.c_str(), nsurfacedots, surfacedots, fr.x, atoms_.get(),
&mtop_->symtab, fr.ePBC, fr.box, bIncludeSolute_);
}
real totalArea, dgsolv;
if (bResAt || bDGsol)
{
- computeAreas(surfaceSel, surfaceSel, frameData.atomAreas_, dgsFactor_,
- &totalArea, &dgsolv, aah, rah, &frameData.res_a_);
+ computeAreas(surfaceSel, surfaceSel, frameData.atomAreas_, dgsFactor_, &totalArea, &dgsolv,
+ aah, rah, &frameData.res_a_);
if (bDGsol)
{
dgh.setPoint(0, dgsolv);
aah.selectDataSet(g + 1);
rah.selectDataSet(g + 1);
}
- computeAreas(surfaceSel, outputSel[g], frameData.atomAreas_, dgsFactor_,
- &totalArea, &dgsolv, aah, rah, &frameData.res_a_);
+ computeAreas(surfaceSel, outputSel[g], frameData.atomAreas_, dgsFactor_, &totalArea,
+ &dgsolv, aah, rah, &frameData.res_a_);
ah.setPoint(g + 1, totalArea);
if (bDGsol)
{
{
totmass += surfaceSel.position(i).mass();
}
- const real density = totmass*AMU/(totvolume*NANO*NANO*NANO);
+ const real density = totmass * AMU / (totvolume * NANO * NANO * NANO);
vh.startFrame(frnr, fr.time);
vh.setPoint(0, totvolume);
vh.setPoint(1, density);
}
}
-void
-Sasa::finishAnalysis(int /*nframes*/)
+void Sasa::finishAnalysis(int /*nframes*/)
{
- //if (bITP)
+ // if (bITP)
//{
// fp3 = ftp2FILE(efITP, nfile, fnm, "w");
// fprintf(fp3, "[ position_restraints ]\n"
//}
}
-void
-Sasa::writeOutput()
-{
-}
+void Sasa::writeOutput() {}
//! \}
-} // namespace
+} // namespace
const char SasaInfo::name[] = "sasa";
-const char SasaInfo::shortDescription[] =
- "Compute solvent accessible surface area";
+const char SasaInfo::shortDescription[] = "Compute solvent accessible surface area";
TrajectoryAnalysisModulePointer SasaInfo::create()
{
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff