2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2006, The GROMACS development team.
6 * Copyright (c) 2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 * Implements gmx::analysismodules::Sasa.
41 * \author Teemu Murtola <teemu.murtola@gmail.com> (C++ conversion)
42 * \ingroup module_trajectoryanalysis
52 #include "gromacs/analysisdata/analysisdata.h"
53 #include "gromacs/analysisdata/modules/average.h"
54 #include "gromacs/analysisdata/modules/plot.h"
55 #include "gromacs/fileio/confio.h"
56 #include "gromacs/fileio/pdbio.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/options/basicoptions.h"
60 #include "gromacs/options/filenameoption.h"
61 #include "gromacs/options/ioptionscontainer.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/selection/selection.h"
64 #include "gromacs/selection/selectionoption.h"
65 #include "gromacs/topology/atomprop.h"
66 #include "gromacs/topology/symtab.h"
67 #include "gromacs/topology/topology.h"
68 #include "gromacs/trajectory/trajectoryframe.h"
69 #include "gromacs/trajectoryanalysis/analysismodule.h"
70 #include "gromacs/trajectoryanalysis/analysissettings.h"
71 #include "gromacs/trajectoryanalysis/topologyinformation.h"
72 #include "gromacs/utility/exceptions.h"
73 #include "gromacs/utility/futil.h"
74 #include "gromacs/utility/pleasecite.h"
75 #include "gromacs/utility/smalloc.h"
76 #include "gromacs/utility/stringutil.h"
77 #include "gromacs/utility/unique_cptr.h"
79 #include "surfacearea.h"
84 namespace analysismodules
90 //! \addtogroup module_trajectoryanalysis
93 //! Tracks information on two nearest neighbors of a single surface dot.
96 //! Index of the second nearest neighbor dot.
98 //! Index of the nearest neighbor dot.
100 //! Squared distance to `aa`.
102 //! Squared distance to `ab`.
107 * Updates nearest neighbor information for a surface dot.
109 * \param[in,out] c Nearest neighbor information array to update.
110 * \param[in] i Index in `c` to update.
111 * \param[in] j Index of the other surface dot to add to the array.
112 * \param[in] d2 Squared distance between `i` and `j`.
114 void add_rec(t_conect c[], int i, int j, real d2)
121 else if (c[i].ab == -1)
126 else if (d2 < c[i].d2a)
131 else if (d2 < c[i].d2b)
136 /* Swap them if necessary: a must be larger than b */
137 if (c[i].d2a < c[i].d2b)
149 * Adds CONECT records for surface dots.
151 * \param[in] fn PDB file to append the CONECT records to.
152 * \param[in] n Number of dots in `x`.
153 * \param[in] x Array of surface dot positions.
155 * Adds a CONECT record that connects each surface dot to its two nearest
156 * neighbors. The function is copied verbatim from the old gmx_sas.c
159 void do_conect(const char* fn, int n, rvec x[])
167 fprintf(stderr, "Building CONECT records\n");
169 for (i = 0; (i < n); i++)
171 c[i].aa = c[i].ab = -1;
174 for (i = 0; (i < n); i++)
176 for (j = i + 1; (j < n); j++)
178 rvec_sub(x[i], x[j], dx);
180 add_rec(c, i, j, d2);
181 add_rec(c, j, i, d2);
184 fp = gmx_ffopen(fn, "a");
185 for (i = 0; (i < n); i++)
187 if ((c[i].aa == -1) || (c[i].ab == -1))
189 fprintf(stderr, "Warning dot %d has no connections\n", i + 1);
191 fprintf(fp, "CONECT%5d%5d%5d\n", i + 1, c[i].aa + 1, c[i].ab + 1);
198 * Plots the surface into a PDB file, optionally including the original atoms.
200 void connolly_plot(const char* fn,
208 gmx_bool bIncludeSolute)
210 const char* const atomnm = "DOT";
211 const char* const resnm = "DOT";
212 const char* const title = "Connolly Dot Surface Generated by gmx sasa";
214 int i, i0, r0, ii0, k;
222 srenew(atoms->atom, atoms->nr + ndots);
223 memset(&atoms->atom[i0], 0, sizeof(*atoms->atom) * ndots);
224 srenew(atoms->atomname, atoms->nr + ndots);
225 srenew(atoms->resinfo, r0 + 1);
226 atoms->atom[i0].resind = r0;
227 t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0 + 1, ' ', 0, ' ');
228 if (atoms->pdbinfo != nullptr)
230 srenew(atoms->pdbinfo, atoms->nr + ndots);
232 snew(xnew, atoms->nr + ndots);
233 for (i = 0; (i < atoms->nr); i++)
235 copy_rvec(x[i], xnew[i]);
237 for (i = k = 0; (i < ndots); i++)
240 atoms->atomname[ii0] = put_symtab(symtab, atomnm);
241 atoms->atom[ii0].resind = r0;
242 xnew[ii0][XX] = dots[k++];
243 xnew[ii0][YY] = dots[k++];
244 xnew[ii0][ZZ] = dots[k++];
245 if (atoms->pdbinfo != nullptr)
247 atoms->pdbinfo[ii0].type = epdbATOM;
248 atoms->pdbinfo[ii0].atomnr = ii0;
249 atoms->pdbinfo[ii0].bfac = 0.0;
250 atoms->pdbinfo[ii0].occup = 0.0;
253 atoms->nr = i0 + ndots;
254 atoms->nres = r0 + 1;
255 write_sto_conf(fn, title, atoms, xnew, nullptr, ePBC, const_cast<rvec*>(box));
261 init_t_atoms(&aaa, ndots, TRUE);
262 aaa.atom[0].resind = 0;
263 t_atoms_set_resinfo(&aaa, 0, symtab, resnm, 1, ' ', 0, ' ');
265 for (i = k = 0; (i < ndots); i++)
268 aaa.atomname[ii0] = put_symtab(symtab, atomnm);
269 aaa.pdbinfo[ii0].type = epdbATOM;
270 aaa.pdbinfo[ii0].atomnr = ii0;
271 aaa.atom[ii0].resind = 0;
272 xnew[ii0][XX] = dots[k++];
273 xnew[ii0][YY] = dots[k++];
274 xnew[ii0][ZZ] = dots[k++];
275 aaa.pdbinfo[ii0].bfac = 0.0;
276 aaa.pdbinfo[ii0].occup = 0.0;
279 write_sto_conf(fn, title, &aaa, xnew, nullptr, ePBC, const_cast<rvec*>(box));
280 do_conect(fn, ndots, xnew);
286 /********************************************************************
287 * Actual analysis module
291 * Implements `gmx sas` trajectory analysis module.
293 class Sasa : public TrajectoryAnalysisModule
298 void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
299 void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
301 TrajectoryAnalysisModuleDataPointer startFrames(const AnalysisDataParallelOptions& opt,
302 const SelectionCollection& selections) override;
303 void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
305 void finishAnalysis(int nframes) override;
306 void writeOutput() override;
310 * Surface areas as a function of time.
312 * First column is for the calculation group, and the rest for the
313 * output groups. This data is always produced.
317 * Per-atom surface areas as a function of time.
319 * Contains one data set for each column in `area_`.
320 * Each column corresponds to a selection position in `surfaceSel_`.
321 * This data is only produced if atom or residue areas have been
324 AnalysisData atomArea_;
326 * Per-residue surface areas as a function of time.
328 * Contains one data set for each column in `area_`.
329 * Each column corresponds to a distinct residue `surfaceSel_`.
330 * For example, if `surfaceSel_` selects residues 2, 5, and 7, there
331 * will be three columns here.
332 * This data is only produced if atom or residue areas have been
335 AnalysisData residueArea_;
337 * Free energy estimates as a function of time.
339 * Column layout is the same as for `area_`.
340 * This data is only produced if the output is requested.
342 AnalysisData dgSolv_;
344 * Total volume and density of the calculation group as a function of
347 * The first column is the volume and the second column is the density.
348 * This data is only produced if the output is requested.
350 AnalysisData volume_;
353 * The selection to calculate the surface for.
355 * Selection::originalId() and Selection::mappedId() store the mapping
356 * from the positions to the columns of `residueArea_`.
357 * The selection is computed with SelectionOption::dynamicMask(), i.e.,
358 * even in the presence of a dynamic selection, the number of returned
359 * positions is fixed, and SelectionPosition::selected() is used.
361 Selection surfaceSel_;
363 * List of optional additional output groups.
365 * Each of these must be a subset of the `surfaceSel_`.
366 * Selection::originalId() and Selection::mappedId() store the mapping
367 * from the positions to the corresponsing positions in `surfaceSel_`.
369 SelectionList outputSel_;
372 std::string fnAtomArea_;
373 std::string fnResidueArea_;
374 std::string fnDGSolv_;
375 std::string fnVolume_;
376 std::string fnConnolly_;
382 bool bIncludeSolute_;
384 //! Global topology corresponding to the input.
386 //! Per-atom data corresponding to the input.
388 //! Combined VdW and probe radii for each atom in the calculation group.
389 std::vector<real> radii_;
391 * Solvation free energy coefficients for each atom in the calculation
394 * Empty if the free energy output has not been requested.
396 std::vector<real> dgsFactor_;
397 //! Calculation algorithm.
398 SurfaceAreaCalculator calculator_;
400 // Copy and assign disallowed by base.
407 bIncludeSolute_(true),
412 registerAnalysisDataset(&area_, "area");
413 registerAnalysisDataset(&atomArea_, "atomarea");
414 registerAnalysisDataset(&residueArea_, "resarea");
415 registerAnalysisDataset(&dgSolv_, "dgsolv");
416 registerAnalysisDataset(&volume_, "volume");
419 void Sasa::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
421 static const char* const desc[] = {
422 "[THISMODULE] computes solvent accessible surface areas.",
423 "See Eisenhaber F, Lijnzaad P, Argos P, Sander C, & Scharf M",
424 "(1995) J. Comput. Chem. 16, 273-284 for the algorithm used.",
425 "With [TT]-q[tt], the Connolly surface can be generated as well",
426 "in a [REF].pdb[ref] file where the nodes are represented as atoms",
427 "and the edges connecting the nearest nodes as CONECT records.",
428 "[TT]-odg[tt] allows for estimation of solvation free energies",
429 "from per-atom solvation energies per exposed surface area.[PAR]",
431 "The program requires a selection for the surface calculation to be",
432 "specified with [TT]-surface[tt]. This should always consist of all",
433 "non-solvent atoms in the system. The area of this group is always",
434 "calculated. Optionally, [TT]-output[tt] can specify additional",
435 "selections, which should be subsets of the calculation group.",
436 "The solvent-accessible areas for these groups are also extracted",
437 "from the full surface.[PAR]",
439 "The average and standard deviation of the area over the trajectory",
440 "can be calculated per residue and atom (options [TT]-or[tt] and", "[TT]-oa[tt]).[PAR]",
441 //"In combination with the latter option an [REF].itp[ref] file can be",
442 //"generated (option [TT]-i[tt])",
443 //"which can be used to restrain surface atoms.[PAR]",
445 "With the [TT]-tv[tt] option the total volume and density of the",
446 "molecule can be computed. With [TT]-pbc[tt] (the default), you",
447 "must ensure that your molecule/surface group is not split across PBC.",
448 "Otherwise, you will get non-sensical results.",
449 "Please also consider whether the normal probe radius is appropriate",
450 "in this case or whether you would rather use, e.g., 0. It is good",
451 "to keep in mind that the results for volume and density are very",
452 "approximate. For example, in ice Ih, one can easily fit water molecules in the",
453 "pores which would yield a volume that is too low, and surface area and density",
454 "that are both too high."
457 settings->setHelpText(desc);
459 options->addOption(FileNameOption("o")
464 .defaultBasename("area")
465 .description("Total area as a function of time"));
467 FileNameOption("odg")
471 .defaultBasename("dgsolv")
472 .description("Estimated solvation free energy as a function of time"));
473 options->addOption(FileNameOption("or")
476 .store(&fnResidueArea_)
477 .defaultBasename("resarea")
478 .description("Average area per residue"));
479 options->addOption(FileNameOption("oa")
483 .defaultBasename("atomarea")
484 .description("Average area per atom"));
485 options->addOption(FileNameOption("tv")
489 .defaultBasename("volume")
490 .description("Total volume and density as a function of time"));
491 options->addOption(FileNameOption("q")
495 .defaultBasename("connolly")
496 .description("PDB file for Connolly surface"));
497 // options->addOption(FileNameOption("i").filetype(eftITP).outputFile()
498 // .store(&fnRestraints_).defaultBasename("surfat")
499 // .description("Topology file for position restraings on surface atoms"));
503 DoubleOption("probe").store(&solsize_).description("Radius of the solvent probe (nm)"));
504 options->addOption(IntegerOption("ndots").store(&ndots_).description(
505 "Number of dots per sphere, more dots means more accuracy"));
506 // options->addOption(DoubleOption("minarea").store(&minarea_)
507 // .description("The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file (see help)"));
509 BooleanOption("prot").store(&bIncludeSolute_).description("Output the protein to the Connolly [REF].pdb[ref] file too"));
511 DoubleOption("dgs").store(&dgsDefault_).description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
513 // Selections must select atoms for the VdW radii lookup to work.
514 // The calculation group uses dynamicMask() so that the coordinates
515 // match a static array of VdW radii.
516 options->addOption(SelectionOption("surface")
521 .description("Surface calculation selection"));
523 SelectionOption("output").storeVector(&outputSel_).onlySortedAtoms().multiValue().description("Output selection(s)"));
525 // Atom names etc. are required for the VdW radii lookup.
526 settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
529 void Sasa::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
532 atoms_ = top.copyAtoms();
534 // bITP = opt2bSet("-i", nfile, fnm);
535 const bool bResAt = !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
536 const bool bDGsol = !fnDGSolv_.empty();
541 fprintf(stderr, "Probe size too small, setting it to %g\n", solsize_);
546 fprintf(stderr, "Ndots too small, setting it to %d\n", ndots_);
549 please_cite(stderr, "Eisenhaber95");
550 // if ((top.ePBC() != epbcXYZ) || (TRICLINIC(fr.box)))
552 // fprintf(stderr, "\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
553 // "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
554 // "will certainly crash the analysis.\n\n");
559 if (!top.hasFullTopology())
561 GMX_THROW(InconsistentInputError(
562 "Cannot compute Delta G of solvation without a tpr file"));
566 if (strcmp(*(atoms_->atomtype[0]), "?") == 0)
568 GMX_THROW(InconsistentInputError(
569 "Your input tpr file is too old (does not contain atom types). Cannot not "
570 "compute Delta G of solvation"));
574 printf("Free energy of solvation predictions:\n");
575 please_cite(stdout, "Eisenberg86a");
580 // Now compute atomic radii including solvent probe size.
581 // Also, fetch solvation free energy coefficients and
582 // compute the residue indices that map the calculation atoms
583 // to the columns of residueArea_.
584 radii_.reserve(surfaceSel_.posCount());
587 dgsFactor_.reserve(surfaceSel_.posCount());
590 const int resCount = surfaceSel_.initOriginalIdsToGroup(top.mtop(), INDEX_RES);
592 // TODO: Not exception-safe, but nice solution would be to have a C++
593 // atom properties class...
596 ArrayRef<const int> atomIndices = surfaceSel_.atomIndices();
598 for (int i = 0; i < surfaceSel_.posCount(); i++)
600 const int ii = atomIndices[i];
601 const int resind = atoms_->atom[ii].resind;
603 if (!aps.setAtomProperty(epropVDW, *(atoms_->resinfo[resind].name), *(atoms_->atomname[ii]), &radius))
607 radii_.push_back(radius + solsize_);
611 if (!aps.setAtomProperty(epropDGsol, *(atoms_->resinfo[resind].name),
612 *(atoms_->atomtype[ii]), &dgsFactor))
614 dgsFactor = dgsDefault_;
616 dgsFactor_.push_back(dgsFactor);
621 fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
624 // Pre-compute mapping from the output groups to the calculation group,
625 // and store it in the selection ID map for easy lookup.
626 for (size_t g = 0; g < outputSel_.size(); ++g)
628 ArrayRef<const int> outputIndices = outputSel_[g].atomIndices();
629 for (int i = 0, j = 0; i < outputSel_[g].posCount(); ++i)
631 while (j < surfaceSel_.posCount() && outputIndices[i] > atomIndices[j])
635 if (j == surfaceSel_.posCount() || outputIndices[i] != atomIndices[j])
637 const std::string message = formatString(
638 "Output selection '%s' is not a subset of "
639 "the surface selection (atom %d is the first "
640 "atom not in the surface selection)",
641 outputSel_[g].name(), outputIndices[i] + 1);
642 GMX_THROW(InconsistentInputError(message));
644 outputSel_[g].setOriginalId(i, j);
648 calculator_.setDotCount(ndots_);
649 calculator_.setRadii(radii_);
651 // Initialize all the output data objects and initialize the output plotters.
653 area_.setColumnCount(0, 1 + outputSel_.size());
655 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
656 plotm->setFileName(fnArea_);
657 plotm->setTitle("Solvent Accessible Surface");
658 plotm->setXAxisIsTime();
659 plotm->setYLabel("Area (nm\\S2\\N)");
660 plotm->appendLegend("Total");
661 for (size_t i = 0; i < outputSel_.size(); ++i)
663 plotm->appendLegend(outputSel_[i].name());
665 area_.addModule(plotm);
670 atomArea_.setDataSetCount(1 + outputSel_.size());
671 residueArea_.setDataSetCount(1 + outputSel_.size());
672 for (size_t i = 0; i <= outputSel_.size(); ++i)
674 atomArea_.setColumnCount(i, surfaceSel_.posCount());
675 residueArea_.setColumnCount(i, resCount);
678 AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
679 for (int i = 0; i < surfaceSel_.posCount(); ++i)
681 avem->setXAxisValue(i, surfaceSel_.position(i).atomIndices()[0] + 1);
683 atomArea_.addModule(avem);
684 if (!fnAtomArea_.empty())
686 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
687 plotm->setFileName(fnAtomArea_);
688 plotm->setTitle("Area per atom over the trajectory");
689 plotm->setXLabel("Atom");
690 plotm->setXFormat(8, 0);
691 plotm->setYLabel("Area (nm\\S2\\N)");
692 plotm->setErrorsAsSeparateColumn(true);
693 plotm->appendLegend("Average (nm\\S2\\N)");
694 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
695 avem->addModule(plotm);
699 AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
701 for (int i = 0; i < surfaceSel_.posCount(); ++i)
703 const int residueGroup = surfaceSel_.position(i).mappedId();
704 if (residueGroup >= nextRow)
706 GMX_ASSERT(residueGroup == nextRow,
707 "Inconsistent (non-uniformly increasing) residue grouping");
708 const int atomIndex = surfaceSel_.position(i).atomIndices()[0];
709 const int residueIndex = atoms_->atom[atomIndex].resind;
710 avem->setXAxisValue(nextRow, atoms_->resinfo[residueIndex].nr);
714 residueArea_.addModule(avem);
715 if (!fnResidueArea_.empty())
717 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
718 plotm->setFileName(fnResidueArea_);
719 plotm->setTitle("Area per residue over the trajectory");
720 plotm->setXLabel("Residue");
721 plotm->setXFormat(8, 0);
722 plotm->setYLabel("Area (nm\\S2\\N)");
723 plotm->setErrorsAsSeparateColumn(true);
724 plotm->appendLegend("Average (nm\\S2\\N)");
725 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
726 avem->addModule(plotm);
731 if (!fnDGSolv_.empty())
733 dgSolv_.setColumnCount(0, 1 + outputSel_.size());
734 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
735 plotm->setFileName(fnDGSolv_);
736 plotm->setTitle("Free Energy of Solvation");
737 plotm->setXAxisIsTime();
738 plotm->setYLabel("D Gsolv");
739 plotm->appendLegend("Total");
740 for (size_t i = 0; i < outputSel_.size(); ++i)
742 plotm->appendLegend(outputSel_[i].name());
744 dgSolv_.addModule(plotm);
747 if (!fnVolume_.empty())
749 volume_.setColumnCount(0, 2);
750 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
751 plotm->setFileName(fnVolume_);
752 plotm->setTitle("Volume and Density");
753 plotm->setXAxisIsTime();
754 plotm->appendLegend("Volume (nm\\S3\\N)");
755 plotm->appendLegend("Density (g/l)");
756 volume_.addModule(plotm);
761 * Temporary memory for use within a single-frame calculation.
763 class SasaModuleData : public TrajectoryAnalysisModuleData
767 * Reserves memory for the frame-local data.
769 * `residueCount` will be zero if per-residue data is not being
772 SasaModuleData(TrajectoryAnalysisModule* module,
773 const AnalysisDataParallelOptions& opt,
774 const SelectionCollection& selections,
777 TrajectoryAnalysisModuleData(module, opt, selections)
779 index_.reserve(atomCount);
780 // If the calculation group is not dynamic, pre-calculate
781 // the index, since it is not going to change.
782 for (int i = 0; i < atomCount; ++i)
786 atomAreas_.resize(atomCount);
787 res_a_.resize(residueCount);
790 void finish() override { finishDataHandles(); }
792 //! Indices of the calculation selection positions selected for the frame.
793 std::vector<int> index_;
795 * Atom areas for each calculation selection position for the frame.
797 * One entry for each position in the calculation group.
798 * Values for atoms not selected are set to zero.
800 std::vector<real> atomAreas_;
802 * Working array to accumulate areas for each residue.
804 * One entry for each distinct residue in the calculation group;
805 * indices are not directly residue numbers or residue indices.
807 * This vector is empty if residue area calculations are not being
810 std::vector<real> res_a_;
813 TrajectoryAnalysisModuleDataPointer Sasa::startFrames(const AnalysisDataParallelOptions& opt,
814 const SelectionCollection& selections)
816 return TrajectoryAnalysisModuleDataPointer(new SasaModuleData(
817 this, opt, selections, surfaceSel_.posCount(), residueArea_.columnCount(0)));
821 * Helper method to compute the areas for a single selection.
823 * \param[in] surfaceSel The calculation selection.
824 * \param[in] sel The selection to compute the areas for (can be
825 * `surfaceSel` or one of the output selections).
826 * \param[in] atomAreas Atom areas for each position in `surfaceSel`.
827 * \param[in] dgsFactor Free energy coefficients for each position in
828 * `surfaceSel`. If empty, free energies are not calculated.
829 * \param[out] totalAreaOut Total area of `sel` (sum of atom areas it selects).
830 * \param[out] dgsolvOut Solvation free energy.
831 * Will be zero of `dgsFactor` is empty.
832 * \param atomAreaHandle Data handle to use for storing atom areas for `sel`.
833 * \param resAreaHandle Data handle to use for storing residue areas for `sel`.
834 * \param resAreaWork Work array for accumulating the residue areas.
835 * If empty, atom and residue areas are not calculated.
837 * `atomAreaHandle` and `resAreaHandle` are not used if `resAreaWork` is empty.
839 void computeAreas(const Selection& surfaceSel,
840 const Selection& sel,
841 const std::vector<real>& atomAreas,
842 const std::vector<real>& dgsFactor,
845 AnalysisDataHandle atomAreaHandle,
846 AnalysisDataHandle resAreaHandle,
847 std::vector<real>* resAreaWork)
849 const bool bResAt = !resAreaWork->empty();
850 const bool bDGsolv = !dgsFactor.empty();
856 std::fill(resAreaWork->begin(), resAreaWork->end(), 0.0_real);
858 for (int i = 0; i < sel.posCount(); ++i)
860 // Get the index of the atom in the calculation group.
861 // For the output groups, the mapping has been precalculated in
863 const int ii = (sel != surfaceSel ? sel.position(i).mappedId() : i);
864 if (!surfaceSel.position(ii).selected())
866 // For the calculation group, skip unselected atoms.
867 if (sel == surfaceSel)
871 GMX_THROW(InconsistentInputError(
872 "Output selection is not a subset of the surface selection"));
874 // Get the internal index of the matching residue.
875 // These have been precalculated in initAnalysis().
876 const int ri = surfaceSel.position(ii).mappedId();
877 const real atomArea = atomAreas[ii];
878 totalArea += atomArea;
881 atomAreaHandle.setPoint(ii, atomArea);
882 (*resAreaWork)[ri] += atomArea;
886 dgsolv += atomArea * dgsFactor[ii];
891 for (size_t i = 0; i < (*resAreaWork).size(); ++i)
893 resAreaHandle.setPoint(i, (*resAreaWork)[i]);
896 *totalAreaOut = totalArea;
900 void Sasa::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
902 AnalysisDataHandle ah = pdata->dataHandle(area_);
903 AnalysisDataHandle dgh = pdata->dataHandle(dgSolv_);
904 AnalysisDataHandle aah = pdata->dataHandle(atomArea_);
905 AnalysisDataHandle rah = pdata->dataHandle(residueArea_);
906 AnalysisDataHandle vh = pdata->dataHandle(volume_);
907 const Selection& surfaceSel = pdata->parallelSelection(surfaceSel_);
908 const SelectionList& outputSel = pdata->parallelSelections(outputSel_);
909 SasaModuleData& frameData = *static_cast<SasaModuleData*>(pdata);
911 const bool bResAt = !frameData.res_a_.empty();
912 const bool bDGsol = !dgsFactor_.empty();
913 const bool bConnolly = (frnr == 0 && !fnConnolly_.empty());
915 // Update indices of selected atoms in the work array.
916 if (surfaceSel.isDynamic())
918 frameData.index_.clear();
919 for (int i = 0; i < surfaceSel.posCount(); ++i)
921 if (surfaceSel.position(i).selected())
923 frameData.index_.push_back(i);
928 // Determine what needs to be calculated.
930 if (bResAt || bDGsol || !outputSel.empty())
932 flag |= FLAG_ATOM_AREA;
938 if (volume_.columnCount() > 0)
943 // Do the low-level calculation.
944 // totarea and totvolume receive the values for the calculation group.
945 // area array contains the per-atom areas for the selected positions.
946 // surfacedots contains nsurfacedots entries, and contains the actual
948 real totarea, totvolume;
949 real *area = nullptr, *surfacedots = nullptr;
951 calculator_.calculate(surfaceSel.coordinates().data(), pbc, frameData.index_.size(),
952 frameData.index_.data(), flag, &totarea, &totvolume, &area, &surfacedots,
954 // Unpack the atomwise areas into the frameData.atomAreas_ array for easier
955 // indexing in the case of dynamic surfaceSel.
958 if (surfaceSel.isDynamic())
960 std::fill(frameData.atomAreas_.begin(), frameData.atomAreas_.end(), 0.0_real);
961 for (size_t i = 0; i < frameData.index_.size(); ++i)
963 frameData.atomAreas_[frameData.index_[i]] = area[i];
968 std::copy(area, area + surfaceSel.posCount(), frameData.atomAreas_.begin());
972 const sfree_guard dotsGuard(surfacedots);
976 if (fr.natoms != mtop_->natoms)
979 InconsistentInputError("Connolly plot (-q) is only supported for trajectories "
980 "that contain all the atoms"));
982 // This is somewhat nasty, as it modifies the atoms and symtab
983 // structures. But since it is only used in the first frame, and no
984 // one else uses the topology after initialization, it may just work
985 // even with future parallelization.
986 connolly_plot(fnConnolly_.c_str(), nsurfacedots, surfacedots, fr.x, atoms_.get(),
987 &mtop_->symtab, fr.ePBC, fr.box, bIncludeSolute_);
990 ah.startFrame(frnr, fr.time);
993 aah.startFrame(frnr, fr.time);
994 rah.startFrame(frnr, fr.time);
998 dgh.startFrame(frnr, fr.time);
1001 ah.setPoint(0, totarea);
1003 real totalArea, dgsolv;
1004 if (bResAt || bDGsol)
1006 computeAreas(surfaceSel, surfaceSel, frameData.atomAreas_, dgsFactor_, &totalArea, &dgsolv,
1007 aah, rah, &frameData.res_a_);
1010 dgh.setPoint(0, dgsolv);
1013 for (size_t g = 0; g < outputSel.size(); ++g)
1017 aah.selectDataSet(g + 1);
1018 rah.selectDataSet(g + 1);
1020 computeAreas(surfaceSel, outputSel[g], frameData.atomAreas_, dgsFactor_, &totalArea,
1021 &dgsolv, aah, rah, &frameData.res_a_);
1022 ah.setPoint(g + 1, totalArea);
1025 dgh.setPoint(g + 1, dgsolv);
1043 for (int i = 0; i < surfaceSel.posCount(); ++i)
1045 totmass += surfaceSel.position(i).mass();
1047 const real density = totmass * AMU / (totvolume * NANO * NANO * NANO);
1048 vh.startFrame(frnr, fr.time);
1049 vh.setPoint(0, totvolume);
1050 vh.setPoint(1, density);
1055 void Sasa::finishAnalysis(int /*nframes*/)
1059 // fp3 = ftp2FILE(efITP, nfile, fnm, "w");
1060 // fprintf(fp3, "[ position_restraints ]\n"
1061 // "#define FCX 1000\n"
1062 // "#define FCY 1000\n"
1063 // "#define FCZ 1000\n"
1064 // "; Atom Type fx fy fz\n");
1065 // for (i = 0; i < nx[0]; i++)
1067 // if (atom_area[i] > minarea)
1069 // fprintf(fp3, "%5d 1 FCX FCX FCZ\n", ii+1);
1076 void Sasa::writeOutput() {}
1082 const char SasaInfo::name[] = "sasa";
1083 const char SasaInfo::shortDescription[] = "Compute solvent accessible surface area";
1085 TrajectoryAnalysisModulePointer SasaInfo::create()
1087 return TrajectoryAnalysisModulePointer(new Sasa);
1090 } // namespace analysismodules