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, 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/legacyheaders/copyrite.h"
54 #include "gromacs/analysisdata/analysisdata.h"
55 #include "gromacs/analysisdata/modules/average.h"
56 #include "gromacs/analysisdata/modules/plot.h"
57 #include "gromacs/fileio/confio.h"
58 #include "gromacs/fileio/pdbio.h"
59 #include "gromacs/fileio/trx.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/options/basicoptions.h"
63 #include "gromacs/options/filenameoption.h"
64 #include "gromacs/options/options.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/selection/selection.h"
67 #include "gromacs/selection/selectionoption.h"
68 #include "gromacs/topology/atomprop.h"
69 #include "gromacs/topology/symtab.h"
70 #include "gromacs/topology/topology.h"
71 #include "gromacs/trajectoryanalysis/analysismodule.h"
72 #include "gromacs/trajectoryanalysis/analysissettings.h"
73 #include "gromacs/utility/exceptions.h"
74 #include "gromacs/utility/futil.h"
75 #include "gromacs/utility/scoped_ptr_sfree.h"
76 #include "gromacs/utility/smalloc.h"
77 #include "gromacs/utility/stringutil.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[], atom_id i, atom_id j, real d2)
116 if (c[i].aa == NO_ATID)
121 else if (c[i].ab == NO_ATID)
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 = NO_ATID;
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 == NO_ATID) || (c[i].ab == NO_ATID))
189 fprintf(stderr, "Warning dot %d has no conections\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, int ndots, real dots[], rvec x[], t_atoms *atoms,
201 t_symtab *symtab, int ePBC, const matrix box, gmx_bool bIncludeSolute)
203 const char *const atomnm = "DOT";
204 const char *const resnm = "DOT";
205 const char *const title = "Connolly Dot Surface Generated by gmx sasa";
207 int i, i0, r0, ii0, k;
215 srenew(atoms->atom, atoms->nr+ndots);
216 memset(&atoms->atom[i0], 0, sizeof(*atoms->atom)*ndots);
217 srenew(atoms->atomname, atoms->nr+ndots);
218 srenew(atoms->resinfo, r0+1);
219 atoms->atom[i0].resind = r0;
220 t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0+1, ' ', 0, ' ');
221 if (atoms->pdbinfo != NULL)
223 srenew(atoms->pdbinfo, atoms->nr+ndots);
225 snew(xnew, atoms->nr+ndots);
226 for (i = 0; (i < atoms->nr); i++)
228 copy_rvec(x[i], xnew[i]);
230 for (i = k = 0; (i < ndots); i++)
233 atoms->atomname[ii0] = put_symtab(symtab, atomnm);
234 atoms->atom[ii0].resind = r0;
235 xnew[ii0][XX] = dots[k++];
236 xnew[ii0][YY] = dots[k++];
237 xnew[ii0][ZZ] = dots[k++];
238 if (atoms->pdbinfo != NULL)
240 atoms->pdbinfo[ii0].type = epdbATOM;
241 atoms->pdbinfo[ii0].atomnr = ii0;
242 atoms->pdbinfo[ii0].bfac = 0.0;
243 atoms->pdbinfo[ii0].occup = 0.0;
246 atoms->nr = i0+ndots;
248 write_sto_conf(fn, title, atoms, xnew, NULL, ePBC, const_cast<rvec *>(box));
254 init_t_atoms(&aaa, ndots, TRUE);
255 aaa.atom[0].resind = 0;
256 t_atoms_set_resinfo(&aaa, 0, symtab, resnm, 1, ' ', 0, ' ');
258 for (i = k = 0; (i < ndots); i++)
261 aaa.atomname[ii0] = put_symtab(symtab, atomnm);
262 aaa.pdbinfo[ii0].type = epdbATOM;
263 aaa.pdbinfo[ii0].atomnr = ii0;
264 aaa.atom[ii0].resind = 0;
265 xnew[ii0][XX] = dots[k++];
266 xnew[ii0][YY] = dots[k++];
267 xnew[ii0][ZZ] = dots[k++];
268 aaa.pdbinfo[ii0].bfac = 0.0;
269 aaa.pdbinfo[ii0].occup = 0.0;
272 write_sto_conf(fn, title, &aaa, xnew, NULL, ePBC, const_cast<rvec *>(box));
273 do_conect(fn, ndots, xnew);
274 free_t_atoms(&aaa, FALSE);
279 /********************************************************************
280 * Actual analysis module
284 * Implements `gmx sas` trajectory analysis module.
286 class Sasa : public TrajectoryAnalysisModule
291 virtual void initOptions(Options *options,
292 TrajectoryAnalysisSettings *settings);
293 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
294 const TopologyInformation &top);
296 virtual TrajectoryAnalysisModuleDataPointer startFrames(
297 const AnalysisDataParallelOptions &opt,
298 const SelectionCollection &selections);
299 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
300 TrajectoryAnalysisModuleData *pdata);
302 virtual void finishAnalysis(int nframes);
303 virtual void writeOutput();
307 * Surface areas as a function of time.
309 * First column is for the calculation group, and the rest for the
310 * output groups. This data is always produced.
314 * Per-atom surface areas as a function of time.
316 * Contains one data set for each column in `area_`.
317 * Each column corresponds to a selection position in `surfaceSel_`.
318 * This data is only produced if atom or residue areas have been
321 AnalysisData atomArea_;
323 * Per-residue surface areas as a function of time.
325 * Contains one data set for each column in `area_`.
326 * Each column corresponds to a distinct residue `surfaceSel_`.
327 * For example, if `surfaceSel_` selects residues 2, 5, and 7, there
328 * will be three columns here.
329 * This data is only produced if atom or residue areas have been
332 AnalysisData residueArea_;
334 * Free energy estimates as a function of time.
336 * Column layout is the same as for `area_`.
337 * This data is only produced if the output is requested.
339 AnalysisData dgSolv_;
341 * Total volume and density of the calculation group as a function of
344 * The first column is the volume and the second column is the density.
345 * This data is only produced if the output is requested.
347 AnalysisData volume_;
350 * The selection to calculate the surface for.
352 * Selection::originalId() and Selection::mappedId() store the mapping
353 * from the positions to the columns of `residueArea_`.
354 * The selection is computed with SelectionOption::dynamicMask(), i.e.,
355 * even in the presence of a dynamic selection, the number of returned
356 * positions is fixed, and SelectionPosition::selected() is used.
358 Selection surfaceSel_;
360 * List of optional additional output groups.
362 * Each of these must be a subset of the `surfaceSel_`.
363 * Selection::originalId() and Selection::mappedId() store the mapping
364 * from the positions to the corresponsing positions in `surfaceSel_`.
366 SelectionList outputSel_;
369 std::string fnAtomArea_;
370 std::string fnResidueArea_;
371 std::string fnDGSolv_;
372 std::string fnVolume_;
373 std::string fnConnolly_;
379 bool bIncludeSolute_;
382 //! Combined VdW and probe radii for each atom in the calculation group.
383 std::vector<real> radii_;
385 * Solvation free energy coefficients for each atom in the calculation
388 * Empty if the free energy output has not been requested.
390 std::vector<real> dgsFactor_;
392 // Copy and assign disallowed by base.
396 : TrajectoryAnalysisModule(SasaInfo::name, SasaInfo::shortDescription),
397 solsize_(0.14), ndots_(24), dgsDefault_(0), bIncludeSolute_(true), top_(NULL)
400 registerAnalysisDataset(&area_, "area");
401 registerAnalysisDataset(&atomArea_, "atomarea");
402 registerAnalysisDataset(&residueArea_, "resarea");
403 registerAnalysisDataset(&dgSolv_, "dgsolv");
404 registerAnalysisDataset(&volume_, "volume");
408 Sasa::initOptions(Options *options, TrajectoryAnalysisSettings *settings)
410 static const char *const desc[] = {
411 "[THISMODULE] computes solvent accessible surface areas.",
412 "See Eisenhaber F, Lijnzaad P, Argos P, Sander C, & Scharf M",
413 "(1995) J. Comput. Chem. 16, 273-284 for the algorithm used.",
414 "With [TT]-q[tt], the Connolly surface can be generated as well",
415 "in a [TT].pdb[tt] file where the nodes are represented as atoms",
416 "and the edges connecting the nearest nodes as CONECT records.",
417 "[TT]-odg[tt] allows for estimation of solvation free energies",
418 "from per-atom solvation energies per exposed surface area.[PAR]",
420 "The program requires a selection for the surface calculation to be",
421 "specified with [TT]-surface[tt]. This should always consist of all",
422 "non-solvent atoms in the system. The area of this group is always",
423 "calculated. Optionally, [TT]-output[tt] can specify additional",
424 "selections, which should be subsets of the calculation group.",
425 "The solvent-accessible areas for these groups are also extracted",
426 "from the full surface.[PAR]",
428 "The average and standard deviation of the area over the trajectory",
429 "can be calculated per residue and atom (options [TT]-or[tt] and",
430 "[TT]-oa[tt]).[PAR]",
431 //"In combination with the latter option an [TT].itp[tt] file can be",
432 //"generated (option [TT]-i[tt])",
433 //"which can be used to restrain surface atoms.[PAR]",
435 "With the [TT]-tv[tt] option the total volume and density of the",
436 "molecule can be computed.",
437 "Please consider whether the normal probe radius is appropriate",
438 "in this case or whether you would rather use, e.g., 0. It is good",
439 "to keep in mind that the results for volume and density are very",
440 "approximate. For example, in ice Ih, one can easily fit water molecules in the",
441 "pores which would yield a volume that is too low, and surface area and density",
442 "that are both too high."
445 options->setDescription(concatenateStrings(desc));
447 options->addOption(FileNameOption("o").filetype(eftPlot).outputFile().required()
448 .store(&fnArea_).defaultBasename("area")
449 .description("Total area as a function of time"));
450 options->addOption(FileNameOption("odg").filetype(eftPlot).outputFile()
451 .store(&fnDGSolv_).defaultBasename("dgsolv")
452 .description("Estimated solvation free energy as a function of time"));
453 options->addOption(FileNameOption("or").filetype(eftPlot).outputFile()
454 .store(&fnResidueArea_).defaultBasename("resarea")
455 .description("Average area per residue"));
456 options->addOption(FileNameOption("oa").filetype(eftPlot).outputFile()
457 .store(&fnAtomArea_).defaultBasename("atomarea")
458 .description("Average area per atom"));
459 options->addOption(FileNameOption("tv").filetype(eftPlot).outputFile()
460 .store(&fnVolume_).defaultBasename("volume")
461 .description("Total volume and density as a function of time"));
462 options->addOption(FileNameOption("q").filetype(eftPDB).outputFile()
463 .store(&fnConnolly_).defaultBasename("connolly")
464 .description("PDB file for Connolly surface"));
465 //options->addOption(FileNameOption("i").filetype(eftITP).outputFile()
466 // .store(&fnRestraints_).defaultBasename("surfat")
467 // .description("Topology file for position restraings on surface atoms"));
470 options->addOption(DoubleOption("probe").store(&solsize_)
471 .description("Radius of the solvent probe (nm)"));
472 options->addOption(IntegerOption("ndots").store(&ndots_)
473 .description("Number of dots per sphere, more dots means more accuracy"));
474 //options->addOption(DoubleOption("minarea").store(&minarea_)
475 // .description("The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file (see help)"));
476 options->addOption(BooleanOption("prot").store(&bIncludeSolute_)
477 .description("Output the protein to the Connolly [TT].pdb[tt] file too"));
478 options->addOption(DoubleOption("dgs").store(&dgsDefault_)
479 .description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
481 // Selections must select atoms for the VdW radii lookup to work.
482 // The calculation group uses dynamicMask() so that the coordinates
483 // match a static array of VdW radii.
484 options->addOption(SelectionOption("surface").store(&surfaceSel_)
485 .required().onlyAtoms().dynamicMask()
486 .description("Surface calculation selection"));
487 options->addOption(SelectionOption("output").storeVector(&outputSel_)
488 .onlyAtoms().multiValue()
489 .description("Output selection(s)"));
491 // Atom names etc. are required for the VdW radii lookup.
492 settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
496 Sasa::initAnalysis(const TrajectoryAnalysisSettings &settings,
497 const TopologyInformation &top)
499 const t_atoms &atoms = top.topology()->atoms;
500 top_ = top.topology();
502 //bITP = opt2bSet("-i", nfile, fnm);
504 !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
505 const bool bDGsol = !fnDGSolv_.empty();
510 fprintf(stderr, "Probe size too small, setting it to %g\n", solsize_);
515 fprintf(stderr, "Ndots too small, setting it to %d\n", ndots_);
518 please_cite(stderr, "Eisenhaber95");
519 //if ((top.ePBC() != epbcXYZ) || (TRICLINIC(fr.box)))
521 // fprintf(stderr, "\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
522 // "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
523 // "will certainly crash the analysis.\n\n");
528 if (!top.hasFullTopology())
530 GMX_THROW(InconsistentInputError("Cannot compute Delta G of solvation without a tpr file"));
534 if (strcmp(*(atoms.atomtype[0]), "?") == 0)
536 GMX_THROW(InconsistentInputError("Your input tpr file is too old (does not contain atom types). Cannot not compute Delta G of solvation"));
540 printf("Free energy of solvation predictions:\n");
541 please_cite(stdout, "Eisenberg86a");
546 // Now compute atomic radii including solvent probe size.
547 // Also, fetch solvation free energy coefficients and
548 // compute the residue indices that map the calculation atoms
549 // to the columns of residueArea_.
550 radii_.reserve(surfaceSel_.posCount());
553 dgsFactor_.reserve(surfaceSel_.posCount());
556 // TODO: Not exception-safe, but nice solution would be to have a C++
557 // atom properties class...
558 gmx_atomprop_t aps = gmx_atomprop_init();
560 ConstArrayRef<int> atomIndices = surfaceSel_.atomIndices();
561 int prevResind = atoms.atom[atomIndices[0]].resind;
564 for (int i = 0; i < surfaceSel_.posCount(); i++)
566 const int ii = atomIndices[i];
567 const int resind = atoms.atom[ii].resind;
568 if (resind != prevResind)
573 surfaceSel_.setOriginalId(i, resCount);
575 if (!gmx_atomprop_query(aps, epropVDW,
576 *(atoms.resinfo[resind].name),
577 *(atoms.atomname[ii]), &radius))
581 radii_.push_back(radius + solsize_);
585 if (!gmx_atomprop_query(aps, epropDGsol,
586 *(atoms.resinfo[resind].name),
587 *(atoms.atomtype[ii]), &dgsFactor))
589 dgsFactor = dgsDefault_;
591 dgsFactor_.push_back(dgsFactor);
596 fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
598 gmx_atomprop_destroy(aps);
600 // The loop above doesn't count the last residue.
603 // Pre-compute mapping from the output groups to the calculation group,
604 // and store it in the selection ID map for easy lookup.
605 for (size_t g = 0; g < outputSel_.size(); ++g)
607 ConstArrayRef<int> outputIndices = outputSel_[g].atomIndices();
608 for (int i = 0, j = 0; i < outputSel_[g].posCount(); ++i)
610 while (j < surfaceSel_.posCount() && outputIndices[i] > atomIndices[j])
614 if (j == surfaceSel_.posCount() || outputIndices[i] != atomIndices[j])
616 GMX_THROW(InconsistentInputError("Output selection is not a subset of the input selection"));
618 outputSel_[g].setOriginalId(i, j);
622 // Initialize all the output data objects and initialize the output plotters.
624 area_.setColumnCount(0, 1 + outputSel_.size());
626 AnalysisDataPlotModulePointer plotm(
627 new AnalysisDataPlotModule(settings.plotSettings()));
628 plotm->setFileName(fnArea_);
629 plotm->setTitle("Solvent Accessible Surface");
630 plotm->setXAxisIsTime();
631 plotm->setYLabel("Area (nm\\S2\\N)");
632 plotm->appendLegend("Total");
633 for (size_t i = 0; i < outputSel_.size(); ++i)
635 plotm->appendLegend(outputSel_[i].name());
637 area_.addModule(plotm);
642 atomArea_.setDataSetCount(1 + outputSel_.size());
643 residueArea_.setDataSetCount(1 + outputSel_.size());
644 for (size_t i = 0; i <= outputSel_.size(); ++i)
646 atomArea_.setColumnCount(i, surfaceSel_.posCount());
647 residueArea_.setColumnCount(i, resCount);
650 AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
651 for (int i = 0; i < surfaceSel_.posCount(); ++i)
653 avem->setXAxisValue(i, surfaceSel_.position(i).atomIndices()[0] + 1);
655 atomArea_.addModule(avem);
656 if (!fnAtomArea_.empty())
658 AnalysisDataPlotModulePointer plotm(
659 new AnalysisDataPlotModule(settings.plotSettings()));
660 plotm->setFileName(fnAtomArea_);
661 plotm->setTitle("Area per residue over the trajectory");
662 plotm->setXLabel("Atom");
663 plotm->setXFormat(8, 0);
664 plotm->setYLabel("Area (nm\\S2\\N)");
665 plotm->setErrorsAsSeparateColumn(true);
666 plotm->appendLegend("Average (nm\\S2\\N)");
667 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
668 avem->addModule(plotm);
672 AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
673 for (int i = 0; i < surfaceSel_.posCount(); ++i)
675 const int atomIndex = surfaceSel_.position(i).atomIndices()[0];
676 const int residueIndex = atoms.atom[atomIndex].resind;
677 avem->setXAxisValue(i, atoms.resinfo[residueIndex].nr);
679 residueArea_.addModule(avem);
680 if (!fnResidueArea_.empty())
682 AnalysisDataPlotModulePointer plotm(
683 new AnalysisDataPlotModule(settings.plotSettings()));
684 plotm->setFileName(fnResidueArea_);
685 plotm->setTitle("Area per atom over the trajectory");
686 plotm->setXLabel("Residue");
687 plotm->setXFormat(8, 0);
688 plotm->setYLabel("Area (nm\\S2\\N)");
689 plotm->setErrorsAsSeparateColumn(true);
690 plotm->appendLegend("Average (nm\\S2\\N)");
691 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
692 avem->addModule(plotm);
697 if (!fnDGSolv_.empty())
699 dgSolv_.setColumnCount(0, 1 + outputSel_.size());
700 AnalysisDataPlotModulePointer plotm(
701 new AnalysisDataPlotModule(settings.plotSettings()));
702 plotm->setFileName(fnDGSolv_);
703 plotm->setTitle("Free Energy of Solvation");
704 plotm->setXAxisIsTime();
705 plotm->setYLabel("D Gsolv");
706 plotm->appendLegend("Total");
707 for (size_t i = 0; i < outputSel_.size(); ++i)
709 plotm->appendLegend(outputSel_[i].name());
711 dgSolv_.addModule(plotm);
714 if (!fnVolume_.empty())
716 volume_.setColumnCount(0, 2);
717 AnalysisDataPlotModulePointer plotm(
718 new AnalysisDataPlotModule(settings.plotSettings()));
719 plotm->setFileName(fnVolume_);
720 plotm->setTitle("Volume and Density");
721 plotm->setXAxisIsTime();
722 plotm->appendLegend("Volume (nm\\S3\\N)");
723 plotm->appendLegend("Density (g/l)");
724 volume_.addModule(plotm);
729 * Temporary memory for use within a single-frame calculation.
731 class SasaModuleData : public TrajectoryAnalysisModuleData
735 * Reserves memory for the frame-local data.
737 * `residueCount` will be zero if per-residue data is not being
740 SasaModuleData(TrajectoryAnalysisModule *module,
741 const AnalysisDataParallelOptions &opt,
742 const SelectionCollection &selections,
743 int atomCount, int residueCount)
744 : TrajectoryAnalysisModuleData(module, opt, selections)
746 index_.reserve(atomCount);
747 // If the calculation group is not dynamic, pre-calculate
748 // the index, since it is not going to change.
749 for (int i = 0; i < atomCount; ++i)
753 atomAreas_.resize(atomCount);
754 res_a_.resize(residueCount);
757 virtual void finish() { finishDataHandles(); }
759 //! Indices of the calculation selection positions selected for the frame.
760 std::vector<int> index_;
762 * Atom areas for each calculation selection position for the frame.
764 * One entry for each position in the calculation group.
765 * Values for atoms not selected are set to zero.
767 std::vector<real> atomAreas_;
769 * Working array to accumulate areas for each residue.
771 * One entry for each distinct residue in the calculation group;
772 * indices are not directly residue numbers or residue indices.
774 * This vector is empty if residue area calculations are not being
777 std::vector<real> res_a_;
780 TrajectoryAnalysisModuleDataPointer Sasa::startFrames(
781 const AnalysisDataParallelOptions &opt,
782 const SelectionCollection &selections)
784 return TrajectoryAnalysisModuleDataPointer(
785 new SasaModuleData(this, opt, selections, surfaceSel_.posCount(),
786 residueArea_.columnCount(0)));
790 * Helper method to compute the areas for a single selection.
792 * \param[in] surfaceSel The calculation selection.
793 * \param[in] sel The selection to compute the areas for (can be
794 * `surfaceSel` or one of the output selections).
795 * \param[in] atomAreas Atom areas for each position in `surfaceSel`.
796 * \param[in] dgsFactor Free energy coefficients for each position in
797 * `surfaceSel`. If empty, free energies are not calculated.
798 * \param[out] totalAreaOut Total area of `sel` (sum of atom areas it selects).
799 * \param[out] dgsolvOut Solvation free energy.
800 * Will be zero of `dgsFactor` is empty.
801 * \param atomAreaHandle Data handle to use for storing atom areas for `sel`.
802 * \param resAreaHandle Data handle to use for storing residue areas for `sel`.
803 * \param resAreaWork Work array for accumulating the residue areas.
804 * If empty, atom and residue areas are not calculated.
806 * `atomAreaHandle` and `resAreaHandle` are not used if `resAreaWork` is empty.
808 void computeAreas(const Selection &surfaceSel, const Selection &sel,
809 const std::vector<real> &atomAreas,
810 const std::vector<real> &dgsFactor,
811 real *totalAreaOut, real *dgsolvOut,
812 AnalysisDataHandle atomAreaHandle,
813 AnalysisDataHandle resAreaHandle,
814 std::vector<real> *resAreaWork)
816 const bool bResAt = !resAreaWork->empty();
817 const bool bDGsolv = !dgsFactor.empty();
823 std::fill(resAreaWork->begin(), resAreaWork->end(),
824 static_cast<real>(0.0));
826 for (int i = 0; i < sel.posCount(); ++i)
828 // Get the index of the atom in the calculation group.
829 // For the output groups, the mapping has been precalculated in
831 const int ii = (sel != surfaceSel ? sel.position(i).mappedId() : i);
832 if (!surfaceSel.position(ii).selected())
834 // For the calculation group, skip unselected atoms.
835 if (sel == surfaceSel)
839 GMX_THROW(InconsistentInputError("Output selection is not a subset of the surface selection"));
841 // Get the internal index of the matching residue.
842 // These have been precalculated in initAnalysis().
843 const int ri = surfaceSel.position(ii).mappedId();
844 const real atomArea = atomAreas[ii];
845 totalArea += atomArea;
848 atomAreaHandle.setPoint(ii, atomArea);
849 (*resAreaWork)[ri] += atomArea;
853 dgsolv += atomArea * dgsFactor[ii];
858 for (size_t i = 0; i < (*resAreaWork).size(); ++i)
860 resAreaHandle.setPoint(i, (*resAreaWork)[i]);
863 *totalAreaOut = totalArea;
868 Sasa::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
869 TrajectoryAnalysisModuleData *pdata)
871 AnalysisDataHandle ah = pdata->dataHandle(area_);
872 AnalysisDataHandle dgh = pdata->dataHandle(dgSolv_);
873 AnalysisDataHandle aah = pdata->dataHandle(atomArea_);
874 AnalysisDataHandle rah = pdata->dataHandle(residueArea_);
875 AnalysisDataHandle vh = pdata->dataHandle(volume_);
876 const Selection &surfaceSel = pdata->parallelSelection(surfaceSel_);
877 const SelectionList &outputSel = pdata->parallelSelections(outputSel_);
878 SasaModuleData &frameData = *static_cast<SasaModuleData *>(pdata);
880 const bool bResAt = !frameData.res_a_.empty();
881 const bool bDGsol = !dgsFactor_.empty();
882 const bool bConnolly = (frnr == 0 && !fnConnolly_.empty());
884 // Update indices of selected atoms in the work array.
885 if (surfaceSel.isDynamic())
887 frameData.index_.clear();
888 for (int i = 0; i < surfaceSel.posCount(); ++i)
890 if (surfaceSel.position(i).selected())
892 frameData.index_.push_back(i);
897 // Determine what needs to be calculated.
899 if (bResAt || bDGsol || !outputSel.empty())
901 flag |= FLAG_ATOM_AREA;
907 if (volume_.columnCount() > 0)
912 // Do the low-level calculation.
913 // totarea and totvolume receive the values for the calculation group.
914 // area array contains the per-atom areas for the selected positions.
915 // surfacedots contains nsurfacedots entries, and contains the actual
917 real totarea, totvolume;
918 real *area = NULL, *surfacedots = NULL;
920 int retval = nsc_dclm_pbc(surfaceSel.coordinates().data(), &radii_[0],
921 frameData.index_.size(), ndots_, flag, &totarea,
922 &area, &totvolume, &surfacedots, &nsurfacedots,
923 &frameData.index_[0],
924 pbc != NULL ? pbc->ePBC : epbcNONE,
925 pbc != NULL ? pbc->box : NULL);
926 // Unpack the atomwise areas into the frameData.atomAreas_ array for easier
927 // indexing in the case of dynamic surfaceSel.
930 if (surfaceSel.isDynamic())
932 std::fill(frameData.atomAreas_.begin(), frameData.atomAreas_.end(),
933 static_cast<real>(0.0));
934 for (size_t i = 0; i < frameData.index_.size(); ++i)
936 frameData.atomAreas_[frameData.index_[i]] = area[i];
941 std::copy(area, area + surfaceSel.posCount(),
942 frameData.atomAreas_.begin());
946 scoped_ptr_sfree dotsGuard(surfacedots);
949 GMX_THROW(InternalError("nsc_dclm_pbc failed"));
954 // This is somewhat nasty, as it modifies the atoms and symtab
955 // structures. But since it is only used in the first frame, and no
956 // one else uses the topology after initialization, it may just work
957 // even with future parallelization.
958 connolly_plot(fnConnolly_.c_str(),
959 nsurfacedots, surfacedots, fr.x, &top_->atoms,
960 &top_->symtab, fr.ePBC, fr.box, bIncludeSolute_);
963 ah.startFrame(frnr, fr.time);
966 aah.startFrame(frnr, fr.time);
967 rah.startFrame(frnr, fr.time);
971 dgh.startFrame(frnr, fr.time);
974 ah.setPoint(0, totarea);
976 real totalArea, dgsolv;
977 if (bResAt || bDGsol)
979 computeAreas(surfaceSel, surfaceSel, frameData.atomAreas_, dgsFactor_,
980 &totalArea, &dgsolv, aah, rah, &frameData.res_a_);
983 dgh.setPoint(0, dgsolv);
986 for (size_t g = 0; g < outputSel.size(); ++g)
990 aah.selectDataSet(g + 1);
991 rah.selectDataSet(g + 1);
993 computeAreas(surfaceSel, outputSel[g], frameData.atomAreas_, dgsFactor_,
994 &totalArea, &dgsolv, aah, rah, &frameData.res_a_);
995 ah.setPoint(g + 1, totalArea);
998 dgh.setPoint(g + 1, dgsolv);
1016 for (int i = 0; i < surfaceSel.posCount(); ++i)
1018 totmass += surfaceSel.position(i).mass();
1020 const real density = totmass*AMU/(totvolume*NANO*NANO*NANO);
1021 vh.startFrame(frnr, fr.time);
1022 vh.setPoint(0, totvolume);
1023 vh.setPoint(1, density);
1029 Sasa::finishAnalysis(int /*nframes*/)
1033 // fp3 = ftp2FILE(efITP, nfile, fnm, "w");
1034 // fprintf(fp3, "[ position_restraints ]\n"
1035 // "#define FCX 1000\n"
1036 // "#define FCY 1000\n"
1037 // "#define FCZ 1000\n"
1038 // "; Atom Type fx fy fz\n");
1039 // for (i = 0; i < nx[0]; i++)
1041 // if (atom_area[i] > minarea)
1043 // fprintf(fp3, "%5d 1 FCX FCX FCZ\n", ii+1);
1059 const char SasaInfo::name[] = "sasa";
1060 const char SasaInfo::shortDescription[] =
1061 "Compute solvent accessible surface area";
1063 TrajectoryAnalysisModulePointer SasaInfo::create()
1065 return TrajectoryAnalysisModulePointer(new Sasa);
1068 } // namespace analysismodules