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
50 #include "gromacs/legacyheaders/copyrite.h"
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/fileio/trx.h"
58 #include "gromacs/math/units.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/options/basicoptions.h"
61 #include "gromacs/options/filenameoption.h"
62 #include "gromacs/options/options.h"
63 #include "gromacs/pbcutil/pbc.h"
64 #include "gromacs/selection/selection.h"
65 #include "gromacs/selection/selectionoption.h"
66 #include "gromacs/topology/atomprop.h"
67 #include "gromacs/topology/symtab.h"
68 #include "gromacs/topology/topology.h"
69 #include "gromacs/trajectoryanalysis/analysismodule.h"
70 #include "gromacs/trajectoryanalysis/analysissettings.h"
71 #include "gromacs/utility/exceptions.h"
72 #include "gromacs/utility/futil.h"
73 #include "gromacs/utility/scoped_ptr_sfree.h"
74 #include "gromacs/utility/smalloc.h"
75 #include "gromacs/utility/stringutil.h"
82 namespace analysismodules
88 //! \addtogroup module_trajectoryanalysis
91 //! Tracks information on two nearest neighbors of a single surface dot.
94 //! Index of the second nearest neighbor dot.
96 //! Index of the nearest neighbor dot.
98 //! Squared distance to `aa`.
100 //! Squared distance to `ab`.
105 * Updates nearest neighbor information for a surface dot.
107 * \param[in,out] c Nearest neighbor information array to update.
108 * \param[in] i Index in `c` to update.
109 * \param[in] j Index of the other surface dot to add to the array.
110 * \param[in] d2 Squared distance between `i` and `j`.
112 void add_rec(t_conect c[], atom_id i, atom_id j, real d2)
114 if (c[i].aa == NO_ATID)
119 else if (c[i].ab == NO_ATID)
124 else if (d2 < c[i].d2a)
129 else if (d2 < c[i].d2b)
134 /* Swap them if necessary: a must be larger than b */
135 if (c[i].d2a < c[i].d2b)
147 * Adds CONECT records for surface dots.
149 * \param[in] fn PDB file to append the CONECT records to.
150 * \param[in] n Number of dots in `x`.
151 * \param[in] x Array of surface dot positions.
153 * Adds a CONECT record that connects each surface dot to its two nearest
154 * neighbors. The function is copied verbatim from the old gmx_sas.c
157 void do_conect(const char *fn, int n, rvec x[])
165 fprintf(stderr, "Building CONECT records\n");
167 for (i = 0; (i < n); i++)
169 c[i].aa = c[i].ab = NO_ATID;
172 for (i = 0; (i < n); i++)
174 for (j = i+1; (j < n); j++)
176 rvec_sub(x[i], x[j], dx);
178 add_rec(c, i, j, d2);
179 add_rec(c, j, i, d2);
182 fp = gmx_ffopen(fn, "a");
183 for (i = 0; (i < n); i++)
185 if ((c[i].aa == NO_ATID) || (c[i].ab == NO_ATID))
187 fprintf(stderr, "Warning dot %d has no conections\n", i+1);
189 fprintf(fp, "CONECT%5d%5d%5d\n", i+1, c[i].aa+1, c[i].ab+1);
196 * Plots the surface into a PDB file, optionally including the original atoms.
198 void connolly_plot(const char *fn, int ndots, real dots[], rvec x[], t_atoms *atoms,
199 t_symtab *symtab, int ePBC, const matrix box, gmx_bool bIncludeSolute)
201 const char *const atomnm = "DOT";
202 const char *const resnm = "DOT";
203 const char *const title = "Connolly Dot Surface Generated by gmx sasa";
205 int i, i0, r0, ii0, k;
213 srenew(atoms->atom, atoms->nr+ndots);
214 memset(&atoms->atom[i0], 0, sizeof(*atoms->atom)*ndots);
215 srenew(atoms->atomname, atoms->nr+ndots);
216 srenew(atoms->resinfo, r0+1);
217 atoms->atom[i0].resind = r0;
218 t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0+1, ' ', 0, ' ');
219 if (atoms->pdbinfo != NULL)
221 srenew(atoms->pdbinfo, atoms->nr+ndots);
223 snew(xnew, atoms->nr+ndots);
224 for (i = 0; (i < atoms->nr); i++)
226 copy_rvec(x[i], xnew[i]);
228 for (i = k = 0; (i < ndots); i++)
231 atoms->atomname[ii0] = put_symtab(symtab, atomnm);
232 atoms->atom[ii0].resind = r0;
233 xnew[ii0][XX] = dots[k++];
234 xnew[ii0][YY] = dots[k++];
235 xnew[ii0][ZZ] = dots[k++];
236 if (atoms->pdbinfo != NULL)
238 atoms->pdbinfo[ii0].type = epdbATOM;
239 atoms->pdbinfo[ii0].atomnr = ii0;
240 atoms->pdbinfo[ii0].bfac = 0.0;
241 atoms->pdbinfo[ii0].occup = 0.0;
244 atoms->nr = i0+ndots;
246 write_sto_conf(fn, title, atoms, xnew, NULL, ePBC, const_cast<rvec *>(box));
252 init_t_atoms(&aaa, ndots, TRUE);
253 aaa.atom[0].resind = 0;
254 t_atoms_set_resinfo(&aaa, 0, symtab, resnm, 1, ' ', 0, ' ');
256 for (i = k = 0; (i < ndots); i++)
259 aaa.atomname[ii0] = put_symtab(symtab, atomnm);
260 aaa.pdbinfo[ii0].type = epdbATOM;
261 aaa.pdbinfo[ii0].atomnr = ii0;
262 aaa.atom[ii0].resind = 0;
263 xnew[ii0][XX] = dots[k++];
264 xnew[ii0][YY] = dots[k++];
265 xnew[ii0][ZZ] = dots[k++];
266 aaa.pdbinfo[ii0].bfac = 0.0;
267 aaa.pdbinfo[ii0].occup = 0.0;
270 write_sto_conf(fn, title, &aaa, xnew, NULL, ePBC, const_cast<rvec *>(box));
271 do_conect(fn, ndots, xnew);
272 free_t_atoms(&aaa, FALSE);
277 /********************************************************************
278 * Actual analysis module
282 * Implements `gmx sas` trajectory analysis module.
284 class Sasa : public TrajectoryAnalysisModule
289 virtual void initOptions(Options *options,
290 TrajectoryAnalysisSettings *settings);
291 virtual void initAnalysis(const TrajectoryAnalysisSettings &settings,
292 const TopologyInformation &top);
294 virtual TrajectoryAnalysisModuleDataPointer startFrames(
295 const AnalysisDataParallelOptions &opt,
296 const SelectionCollection &selections);
297 virtual void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
298 TrajectoryAnalysisModuleData *pdata);
300 virtual void finishAnalysis(int nframes);
301 virtual void writeOutput();
305 * Surface areas as a function of time.
307 * First column is for the calculation group, and the rest for the
308 * output groups. This data is always produced.
312 * Per-atom surface areas as a function of time.
314 * Contains one data set for each column in `area_`.
315 * Each column corresponds to a selection position in `surfaceSel_`.
316 * This data is only produced if atom or residue areas have been
319 AnalysisData atomArea_;
321 * Per-residue surface areas as a function of time.
323 * Contains one data set for each column in `area_`.
324 * Each column corresponds to a distinct residue `surfaceSel_`.
325 * For example, if `surfaceSel_` selects residues 2, 5, and 7, there
326 * will be three columns here.
327 * This data is only produced if atom or residue areas have been
330 AnalysisData residueArea_;
332 * Free energy estimates as a function of time.
334 * Column layout is the same as for `area_`.
335 * This data is only produced if the output is requested.
337 AnalysisData dgSolv_;
339 * Total volume and density of the calculation group as a function of
342 * The first column is the volume and the second column is the density.
343 * This data is only produced if the output is requested.
345 AnalysisData volume_;
348 * The selection to calculate the surface for.
350 * Selection::originalId() and Selection::mappedId() store the mapping
351 * from the positions to the columns of `residueArea_`.
352 * The selection is computed with SelectionOption::dynamicMask(), i.e.,
353 * even in the presence of a dynamic selection, the number of returned
354 * positions is fixed, and SelectionPosition::selected() is used.
356 Selection surfaceSel_;
358 * List of optional additional output groups.
360 * Each of these must be a subset of the `surfaceSel_`.
361 * Selection::originalId() and Selection::mappedId() store the mapping
362 * from the positions to the corresponsing positions in `surfaceSel_`.
364 SelectionList outputSel_;
367 std::string fnAtomArea_;
368 std::string fnResidueArea_;
369 std::string fnDGSolv_;
370 std::string fnVolume_;
371 std::string fnConnolly_;
377 bool bIncludeSolute_;
380 //! Combined VdW and probe radii for each atom in the calculation group.
381 std::vector<real> radii_;
383 * Solvation free energy coefficients for each atom in the calculation
386 * Empty if the free energy output has not been requested.
388 std::vector<real> dgsFactor_;
390 // Copy and assign disallowed by base.
394 : TrajectoryAnalysisModule(SasaInfo::name, SasaInfo::shortDescription),
395 solsize_(0.14), ndots_(24), dgsDefault_(0), bIncludeSolute_(true), top_(NULL)
398 registerAnalysisDataset(&area_, "area");
399 registerAnalysisDataset(&atomArea_, "atomarea");
400 registerAnalysisDataset(&residueArea_, "resarea");
401 registerAnalysisDataset(&dgSolv_, "dgsolv");
402 registerAnalysisDataset(&volume_, "volume");
406 Sasa::initOptions(Options *options, TrajectoryAnalysisSettings *settings)
408 static const char *const desc[] = {
409 "[THISMODULE] computes solvent accessible surface areas.",
410 "See Eisenhaber F, Lijnzaad P, Argos P, Sander C, & Scharf M",
411 "(1995) J. Comput. Chem. 16, 273-284 for the algorithm used.",
412 "With [TT]-q[tt], the Connolly surface can be generated as well",
413 "in a [TT].pdb[tt] file where the nodes are represented as atoms",
414 "and the edges connecting the nearest nodes as CONECT records.",
415 "[TT]-odg[tt] allows for estimation of solvation free energies",
416 "from per-atom solvation energies per exposed surface area.[PAR]",
418 "The program requires a selection for the surface calculation to be",
419 "specified with [TT]-surface[tt]. This should always consist of all",
420 "non-solvent atoms in the system. The area of this group is always",
421 "calculated. Optionally, [TT]-output[tt] can specify additional",
422 "selections, which should be subsets of the calculation group.",
423 "The solvent-accessible areas for these groups are also extracted",
424 "from the full surface.[PAR]",
426 "The average and standard deviation of the area over the trajectory",
427 "can be calculated per residue and atom (options [TT]-or[tt] and",
428 "[TT]-oa[tt]).[PAR]",
429 //"In combination with the latter option an [TT].itp[tt] file can be",
430 //"generated (option [TT]-i[tt])",
431 //"which can be used to restrain surface atoms.[PAR]",
433 "With the [TT]-tv[tt] option the total volume and density of the",
434 "molecule can be computed.",
435 "Please consider whether the normal probe radius is appropriate",
436 "in this case or whether you would rather use, e.g., 0. It is good",
437 "to keep in mind that the results for volume and density are very",
438 "approximate. For example, in ice Ih, one can easily fit water molecules in the",
439 "pores which would yield a volume that is too low, and surface area and density",
440 "that are both too high."
443 options->setDescription(concatenateStrings(desc));
445 options->addOption(FileNameOption("o").filetype(eftPlot).outputFile().required()
446 .store(&fnArea_).defaultBasename("area")
447 .description("Total area as a function of time"));
448 options->addOption(FileNameOption("odg").filetype(eftPlot).outputFile()
449 .store(&fnDGSolv_).defaultBasename("dgsolv")
450 .description("Estimated solvation free energy as a function of time"));
451 options->addOption(FileNameOption("or").filetype(eftPlot).outputFile()
452 .store(&fnResidueArea_).defaultBasename("resarea")
453 .description("Average area per residue"));
454 options->addOption(FileNameOption("oa").filetype(eftPlot).outputFile()
455 .store(&fnAtomArea_).defaultBasename("atomarea")
456 .description("Average area per atom"));
457 options->addOption(FileNameOption("tv").filetype(eftPlot).outputFile()
458 .store(&fnVolume_).defaultBasename("volume")
459 .description("Total volume and density as a function of time"));
460 options->addOption(FileNameOption("q").filetype(eftPDB).outputFile()
461 .store(&fnConnolly_).defaultBasename("connolly")
462 .description("PDB file for Connolly surface"));
463 //options->addOption(FileNameOption("i").filetype(eftITP).outputFile()
464 // .store(&fnRestraints_).defaultBasename("surfat")
465 // .description("Topology file for position restraings on surface atoms"));
468 options->addOption(DoubleOption("probe").store(&solsize_)
469 .description("Radius of the solvent probe (nm)"));
470 options->addOption(IntegerOption("ndots").store(&ndots_)
471 .description("Number of dots per sphere, more dots means more accuracy"));
472 //options->addOption(DoubleOption("minarea").store(&minarea_)
473 // .description("The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file (see help)"));
474 options->addOption(BooleanOption("prot").store(&bIncludeSolute_)
475 .description("Output the protein to the Connolly [TT].pdb[tt] file too"));
476 options->addOption(DoubleOption("dgs").store(&dgsDefault_)
477 .description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
479 // Selections must select atoms for the VdW radii lookup to work.
480 // The calculation group uses dynamicMask() so that the coordinates
481 // match a static array of VdW radii.
482 options->addOption(SelectionOption("surface").store(&surfaceSel_)
483 .required().onlyAtoms().dynamicMask()
484 .description("Surface calculation selection"));
485 options->addOption(SelectionOption("output").storeVector(&outputSel_)
486 .onlyAtoms().multiValue()
487 .description("Output selection(s)"));
489 // Atom names etc. are required for the VdW radii lookup.
490 settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
494 Sasa::initAnalysis(const TrajectoryAnalysisSettings &settings,
495 const TopologyInformation &top)
497 const t_atoms &atoms = top.topology()->atoms;
498 top_ = top.topology();
500 //bITP = opt2bSet("-i", nfile, fnm);
502 !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
503 const bool bDGsol = !fnDGSolv_.empty();
508 fprintf(stderr, "Probe size too small, setting it to %g\n", solsize_);
513 fprintf(stderr, "Ndots too small, setting it to %d\n", ndots_);
516 please_cite(stderr, "Eisenhaber95");
517 //if ((top.ePBC() != epbcXYZ) || (TRICLINIC(fr.box)))
519 // fprintf(stderr, "\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
520 // "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
521 // "will certainly crash the analysis.\n\n");
526 if (!top.hasFullTopology())
528 GMX_THROW(InconsistentInputError("Cannot compute Delta G of solvation without a tpr file"));
532 if (strcmp(*(atoms.atomtype[0]), "?") == 0)
534 GMX_THROW(InconsistentInputError("Your input tpr file is too old (does not contain atom types). Cannot not compute Delta G of solvation"));
538 printf("Free energy of solvation predictions:\n");
539 please_cite(stdout, "Eisenberg86a");
544 // Now compute atomic radii including solvent probe size.
545 // Also, fetch solvation free energy coefficients and
546 // compute the residue indices that map the calculation atoms
547 // to the columns of residueArea_.
548 radii_.reserve(surfaceSel_.posCount());
551 dgsFactor_.reserve(surfaceSel_.posCount());
554 // TODO: Not exception-safe, but nice solution would be to have a C++
555 // atom properties class...
556 gmx_atomprop_t aps = gmx_atomprop_init();
558 ConstArrayRef<int> atomIndices = surfaceSel_.atomIndices();
559 int prevResind = atoms.atom[atomIndices[0]].resind;
562 for (int i = 0; i < surfaceSel_.posCount(); i++)
564 const int ii = atomIndices[i];
565 const int resind = atoms.atom[ii].resind;
566 if (resind != prevResind)
571 surfaceSel_.setOriginalId(i, resCount);
573 if (!gmx_atomprop_query(aps, epropVDW,
574 *(atoms.resinfo[resind].name),
575 *(atoms.atomname[ii]), &radius))
579 radii_.push_back(radius + solsize_);
583 if (!gmx_atomprop_query(aps, epropDGsol,
584 *(atoms.resinfo[resind].name),
585 *(atoms.atomtype[ii]), &dgsFactor))
587 dgsFactor = dgsDefault_;
589 dgsFactor_.push_back(dgsFactor);
594 fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
596 gmx_atomprop_destroy(aps);
598 // The loop above doesn't count the last residue.
601 // Pre-compute mapping from the output groups to the calculation group,
602 // and store it in the selection ID map for easy lookup.
603 for (size_t g = 0; g < outputSel_.size(); ++g)
605 ConstArrayRef<int> outputIndices = outputSel_[g].atomIndices();
606 for (int i = 0, j = 0; i < outputSel_[g].posCount(); ++i)
608 while (j < surfaceSel_.posCount() && outputIndices[i] > atomIndices[j])
612 if (j == surfaceSel_.posCount() || outputIndices[i] != atomIndices[j])
614 GMX_THROW(InconsistentInputError("Output selection is not a subset of the input selection"));
616 outputSel_[g].setOriginalId(i, j);
620 // Initialize all the output data objects and initialize the output plotters.
622 area_.setColumnCount(0, 1 + outputSel_.size());
624 AnalysisDataPlotModulePointer plotm(
625 new AnalysisDataPlotModule(settings.plotSettings()));
626 plotm->setFileName(fnArea_);
627 plotm->setTitle("Solvent Accessible Surface");
628 plotm->setXAxisIsTime();
629 plotm->setYLabel("Area (nm\\S2\\N)");
630 plotm->appendLegend("Total");
631 for (size_t i = 0; i < outputSel_.size(); ++i)
633 plotm->appendLegend(outputSel_[i].name());
635 area_.addModule(plotm);
640 atomArea_.setDataSetCount(1 + outputSel_.size());
641 residueArea_.setDataSetCount(1 + outputSel_.size());
642 for (size_t i = 0; i <= outputSel_.size(); ++i)
644 atomArea_.setColumnCount(i, surfaceSel_.posCount());
645 residueArea_.setColumnCount(i, resCount);
648 AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
649 for (int i = 0; i < surfaceSel_.posCount(); ++i)
651 avem->setXAxisValue(i, surfaceSel_.position(i).atomIndices()[0] + 1);
653 atomArea_.addModule(avem);
654 if (!fnAtomArea_.empty())
656 AnalysisDataPlotModulePointer plotm(
657 new AnalysisDataPlotModule(settings.plotSettings()));
658 plotm->setFileName(fnAtomArea_);
659 plotm->setTitle("Area per residue over the trajectory");
660 plotm->setXLabel("Atom");
661 plotm->setXFormat(8, 0);
662 plotm->setYLabel("Area (nm\\S2\\N)");
663 plotm->setErrorsAsSeparateColumn(true);
664 plotm->appendLegend("Average (nm\\S2\\N)");
665 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
666 avem->addModule(plotm);
670 AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
671 for (int i = 0; i < surfaceSel_.posCount(); ++i)
673 const int atomIndex = surfaceSel_.position(i).atomIndices()[0];
674 const int residueIndex = atoms.atom[atomIndex].resind;
675 avem->setXAxisValue(i, atoms.resinfo[residueIndex].nr);
677 residueArea_.addModule(avem);
678 if (!fnResidueArea_.empty())
680 AnalysisDataPlotModulePointer plotm(
681 new AnalysisDataPlotModule(settings.plotSettings()));
682 plotm->setFileName(fnResidueArea_);
683 plotm->setTitle("Area per atom over the trajectory");
684 plotm->setXLabel("Residue");
685 plotm->setXFormat(8, 0);
686 plotm->setYLabel("Area (nm\\S2\\N)");
687 plotm->setErrorsAsSeparateColumn(true);
688 plotm->appendLegend("Average (nm\\S2\\N)");
689 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
690 avem->addModule(plotm);
695 if (!fnDGSolv_.empty())
697 dgSolv_.setColumnCount(0, 1 + outputSel_.size());
698 AnalysisDataPlotModulePointer plotm(
699 new AnalysisDataPlotModule(settings.plotSettings()));
700 plotm->setFileName(fnDGSolv_);
701 plotm->setTitle("Free Energy of Solvation");
702 plotm->setXAxisIsTime();
703 plotm->setYLabel("D Gsolv");
704 plotm->appendLegend("Total");
705 for (size_t i = 0; i < outputSel_.size(); ++i)
707 plotm->appendLegend(outputSel_[i].name());
709 dgSolv_.addModule(plotm);
712 if (!fnVolume_.empty())
714 volume_.setColumnCount(0, 2);
715 AnalysisDataPlotModulePointer plotm(
716 new AnalysisDataPlotModule(settings.plotSettings()));
717 plotm->setFileName(fnVolume_);
718 plotm->setTitle("Volume and Density");
719 plotm->setXAxisIsTime();
720 plotm->appendLegend("Volume (nm\\S3\\N)");
721 plotm->appendLegend("Density (g/l)");
722 volume_.addModule(plotm);
727 * Temporary memory for use within a single-frame calculation.
729 class SasaModuleData : public TrajectoryAnalysisModuleData
733 * Reserves memory for the frame-local data.
735 * `residueCount` will be zero if per-residue data is not being
738 SasaModuleData(TrajectoryAnalysisModule *module,
739 const AnalysisDataParallelOptions &opt,
740 const SelectionCollection &selections,
741 int atomCount, int residueCount)
742 : TrajectoryAnalysisModuleData(module, opt, selections)
744 index_.reserve(atomCount);
745 // If the calculation group is not dynamic, pre-calculate
746 // the index, since it is not going to change.
747 for (int i = 0; i < atomCount; ++i)
751 atomAreas_.resize(atomCount);
752 res_a_.resize(residueCount);
755 virtual void finish() { finishDataHandles(); }
757 //! Indices of the calculation selection positions selected for the frame.
758 std::vector<int> index_;
760 * Atom areas for each calculation selection position for the frame.
762 * One entry for each position in the calculation group.
763 * Values for atoms not selected are set to zero.
765 std::vector<real> atomAreas_;
767 * Working array to accumulate areas for each residue.
769 * One entry for each distinct residue in the calculation group;
770 * indices are not directly residue numbers or residue indices.
772 * This vector is empty if residue area calculations are not being
775 std::vector<real> res_a_;
778 TrajectoryAnalysisModuleDataPointer Sasa::startFrames(
779 const AnalysisDataParallelOptions &opt,
780 const SelectionCollection &selections)
782 return TrajectoryAnalysisModuleDataPointer(
783 new SasaModuleData(this, opt, selections, surfaceSel_.posCount(),
784 residueArea_.columnCount(0)));
788 * Helper method to compute the areas for a single selection.
790 * \param[in] surfaceSel The calculation selection.
791 * \param[in] sel The selection to compute the areas for (can be
792 * `surfaceSel` or one of the output selections).
793 * \param[in] atomAreas Atom areas for each position in `surfaceSel`.
794 * \param[in] dgsFactor Free energy coefficients for each position in
795 * `surfaceSel`. If empty, free energies are not calculated.
796 * \param[out] totalAreaOut Total area of `sel` (sum of atom areas it selects).
797 * \param[out] dgsolvOut Solvation free energy.
798 * Will be zero of `dgsFactor` is empty.
799 * \param atomAreaHandle Data handle to use for storing atom areas for `sel`.
800 * \param resAreaHandle Data handle to use for storing residue areas for `sel`.
801 * \param resAreaWork Work array for accumulating the residue areas.
802 * If empty, atom and residue areas are not calculated.
804 * `atomAreaHandle` and `resAreaHandle` are not used if `resAreaWork` is empty.
806 void computeAreas(const Selection &surfaceSel, const Selection &sel,
807 const std::vector<real> &atomAreas,
808 const std::vector<real> &dgsFactor,
809 real *totalAreaOut, real *dgsolvOut,
810 AnalysisDataHandle atomAreaHandle,
811 AnalysisDataHandle resAreaHandle,
812 std::vector<real> *resAreaWork)
814 const bool bResAt = !resAreaWork->empty();
815 const bool bDGsolv = !dgsFactor.empty();
821 std::fill(resAreaWork->begin(), resAreaWork->end(),
822 static_cast<real>(0.0));
824 for (int i = 0; i < sel.posCount(); ++i)
826 // Get the index of the atom in the calculation group.
827 // For the output groups, the mapping has been precalculated in
829 const int ii = (sel != surfaceSel ? sel.position(i).mappedId() : i);
830 if (!surfaceSel.position(ii).selected())
832 // For the calculation group, skip unselected atoms.
833 if (sel == surfaceSel)
837 GMX_THROW(InconsistentInputError("Output selection is not a subset of the surface selection"));
839 // Get the internal index of the matching residue.
840 // These have been precalculated in initAnalysis().
841 const int ri = surfaceSel.position(ii).mappedId();
842 const real atomArea = atomAreas[ii];
843 totalArea += atomArea;
846 atomAreaHandle.setPoint(ii, atomArea);
847 (*resAreaWork)[ri] += atomArea;
851 dgsolv += atomArea * dgsFactor[ii];
856 for (size_t i = 0; i < (*resAreaWork).size(); ++i)
858 resAreaHandle.setPoint(i, (*resAreaWork)[i]);
861 *totalAreaOut = totalArea;
866 Sasa::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
867 TrajectoryAnalysisModuleData *pdata)
869 AnalysisDataHandle ah = pdata->dataHandle(area_);
870 AnalysisDataHandle dgh = pdata->dataHandle(dgSolv_);
871 AnalysisDataHandle aah = pdata->dataHandle(atomArea_);
872 AnalysisDataHandle rah = pdata->dataHandle(residueArea_);
873 AnalysisDataHandle vh = pdata->dataHandle(volume_);
874 const Selection &surfaceSel = pdata->parallelSelection(surfaceSel_);
875 const SelectionList &outputSel = pdata->parallelSelections(outputSel_);
876 SasaModuleData &frameData = *static_cast<SasaModuleData *>(pdata);
878 const bool bResAt = !frameData.res_a_.empty();
879 const bool bDGsol = !dgsFactor_.empty();
880 const bool bConnolly = (frnr == 0 && !fnConnolly_.empty());
882 // Update indices of selected atoms in the work array.
883 if (surfaceSel.isDynamic())
885 frameData.index_.clear();
886 for (int i = 0; i < surfaceSel.posCount(); ++i)
888 if (surfaceSel.position(i).selected())
890 frameData.index_.push_back(i);
895 // Determine what needs to be calculated.
897 if (bResAt || bDGsol || !outputSel.empty())
899 flag |= FLAG_ATOM_AREA;
905 if (volume_.columnCount() > 0)
910 // Do the low-level calculation.
911 // totarea and totvolume receive the values for the calculation group.
912 // area array contains the per-atom areas for the selected positions.
913 // surfacedots contains nsurfacedots entries, and contains the actual
915 real totarea, totvolume;
916 real *area = NULL, *surfacedots = NULL;
918 int retval = nsc_dclm_pbc(surfaceSel.coordinates().data(), &radii_[0],
919 frameData.index_.size(), ndots_, flag, &totarea,
920 &area, &totvolume, &surfacedots, &nsurfacedots,
921 &frameData.index_[0],
922 pbc != NULL ? pbc->ePBC : epbcNONE,
923 pbc != NULL ? pbc->box : NULL);
924 // Unpack the atomwise areas into the frameData.atomAreas_ array for easier
925 // indexing in the case of dynamic surfaceSel.
928 if (surfaceSel.isDynamic())
930 std::fill(frameData.atomAreas_.begin(), frameData.atomAreas_.end(),
931 static_cast<real>(0.0));
932 for (size_t i = 0; i < frameData.index_.size(); ++i)
934 frameData.atomAreas_[frameData.index_[i]] = area[i];
939 std::copy(area, area + surfaceSel.posCount(),
940 frameData.atomAreas_.begin());
944 scoped_ptr_sfree dotsGuard(surfacedots);
947 GMX_THROW(InternalError("nsc_dclm_pbc failed"));
952 // This is somewhat nasty, as it modifies the atoms and symtab
953 // structures. But since it is only used in the first frame, and no
954 // one else uses the topology after initialization, it may just work
955 // even with future parallelization.
956 connolly_plot(fnConnolly_.c_str(),
957 nsurfacedots, surfacedots, fr.x, &top_->atoms,
958 &top_->symtab, fr.ePBC, fr.box, bIncludeSolute_);
961 ah.startFrame(frnr, fr.time);
964 aah.startFrame(frnr, fr.time);
965 rah.startFrame(frnr, fr.time);
969 dgh.startFrame(frnr, fr.time);
972 ah.setPoint(0, totarea);
974 real totalArea, dgsolv;
975 if (bResAt || bDGsol)
977 computeAreas(surfaceSel, surfaceSel, frameData.atomAreas_, dgsFactor_,
978 &totalArea, &dgsolv, aah, rah, &frameData.res_a_);
981 dgh.setPoint(0, dgsolv);
984 for (size_t g = 0; g < outputSel.size(); ++g)
988 aah.selectDataSet(g + 1);
989 rah.selectDataSet(g + 1);
991 computeAreas(surfaceSel, outputSel[g], frameData.atomAreas_, dgsFactor_,
992 &totalArea, &dgsolv, aah, rah, &frameData.res_a_);
993 ah.setPoint(g + 1, totalArea);
996 dgh.setPoint(g + 1, dgsolv);
1014 for (int i = 0; i < surfaceSel.posCount(); ++i)
1016 totmass += surfaceSel.position(i).mass();
1018 const real density = totmass*AMU/(totvolume*NANO*NANO*NANO);
1019 vh.startFrame(frnr, fr.time);
1020 vh.setPoint(0, totvolume);
1021 vh.setPoint(1, density);
1027 Sasa::finishAnalysis(int /*nframes*/)
1031 // fp3 = ftp2FILE(efITP, nfile, fnm, "w");
1032 // fprintf(fp3, "[ position_restraints ]\n"
1033 // "#define FCX 1000\n"
1034 // "#define FCY 1000\n"
1035 // "#define FCZ 1000\n"
1036 // "; Atom Type fx fy fz\n");
1037 // for (i = 0; i < nx[0]; i++)
1039 // if (atom_area[i] > minarea)
1041 // fprintf(fp3, "%5d 1 FCX FCX FCZ\n", ii+1);
1057 const char SasaInfo::name[] = "sasa";
1058 const char SasaInfo::shortDescription[] =
1059 "Compute solvent accessible surface area";
1061 TrajectoryAnalysisModulePointer SasaInfo::create()
1063 return TrajectoryAnalysisModulePointer(new Sasa);
1066 } // namespace analysismodules