src/gromacs/trajectoryanalysis/modules/sasa.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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 by the GROMACS development team.
   7  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   8  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   9  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  10  * and including many others, as listed in the AUTHORS file in the
  11  * top-level source directory and at http://www.gromacs.org.
  12  *
  13  * GROMACS is free software; you can redistribute it and/or
  14  * modify it under the terms of the GNU Lesser General Public License
  15  * as published by the Free Software Foundation; either version 2.1
  16  * of the License, or (at your option) any later version.
  17  *
  18  * GROMACS is distributed in the hope that it will be useful,
  19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  21  * Lesser General Public License for more details.
  22  *
  23  * You should have received a copy of the GNU Lesser General Public
  24  * License along with GROMACS; if not, see
  25  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  26  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  27  *
  28  * If you want to redistribute modifications to GROMACS, please
  29  * consider that scientific software is very special. Version
  30  * control is crucial - bugs must be traceable. We will be happy to
  31  * consider code for inclusion in the official distribution, but
  32  * derived work must not be called official GROMACS. Details are found
  33  * in the README & COPYING files - if they are missing, get the
  34  * official version at http://www.gromacs.org.
  35  *
  36  * To help us fund GROMACS development, we humbly ask that you cite
  37  * the research papers on the package. Check out http://www.gromacs.org.
  38  */
  39 /*! \internal \file
  40  * \brief
  41  * Implements gmx::analysismodules::Sasa.
  42  *
  43  * \author Teemu Murtola <teemu.murtola@gmail.com> (C++ conversion)
  44  * \ingroup module_trajectoryanalysis
  45  */
  46 #include "gmxpre.h"
  47
  48 #include "sasa.h"
  49
  50 #include <algorithm>
  51 #include <string>
  52 #include <vector>
  53
  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/math/units.h"
  60 #include "gromacs/math/vec.h"
  61 #include "gromacs/options/basicoptions.h"
  62 #include "gromacs/options/filenameoption.h"
  63 #include "gromacs/options/ioptionscontainer.h"
  64 #include "gromacs/pbcutil/pbc.h"
  65 #include "gromacs/selection/selection.h"
  66 #include "gromacs/selection/selectionoption.h"
  67 #include "gromacs/topology/atomprop.h"
  68 #include "gromacs/topology/symtab.h"
  69 #include "gromacs/topology/topology.h"
  70 #include "gromacs/trajectory/trajectoryframe.h"
  71 #include "gromacs/trajectoryanalysis/analysismodule.h"
  72 #include "gromacs/trajectoryanalysis/analysissettings.h"
  73 #include "gromacs/trajectoryanalysis/topologyinformation.h"
  74 #include "gromacs/utility/exceptions.h"
  75 #include "gromacs/utility/futil.h"
  76 #include "gromacs/utility/pleasecite.h"
  77 #include "gromacs/utility/smalloc.h"
  78 #include "gromacs/utility/stringutil.h"
  79 #include "gromacs/utility/unique_cptr.h"
  80
  81 #include "surfacearea.h"
  82
  83 namespace gmx
  84 {
  85
  86 namespace analysismodules
  87 {
  88
  89 namespace
  90 {
  91
  92 //! \addtogroup module_trajectoryanalysis
  93 //! \{
  94
  95 //! Tracks information on two nearest neighbors of a single surface dot.
  96 struct t_conect
  97 {
  98     //! Index of the second nearest neighbor dot.
  99     int aa;
 100     //! Index of the nearest neighbor dot.
 101     int ab;
 102     //! Squared distance to `aa`.
 103     real d2a;
 104     //! Squared distance to `ab`.
 105     real d2b;
 106 };
 107
 108 /*! \brief
 109  * Updates nearest neighbor information for a surface dot.
 110  *
 111  * \param[in,out] c  Nearest neighbor information array to update.
 112  * \param[in]     i  Index in `c` to update.
 113  * \param[in]     j  Index of the other surface dot to add to the array.
 114  * \param[in]     d2 Squared distance between `i` and `j`.
 115  */
 116 void add_rec(t_conect c[], int i, int j, real d2)
 117 {
 118     if (c[i].aa == -1)
 119     { // NOLINT bugprone-branch-clone
 120         c[i].aa  = j;
 121         c[i].d2a = d2;
 122     }
 123     else if (c[i].ab == -1)
 124     { // NOLINT bugprone-branch-clone
 125         c[i].ab  = j;
 126         c[i].d2b = d2;
 127     }
 128     else if (d2 < c[i].d2a)
 129     {
 130         c[i].aa  = j;
 131         c[i].d2a = d2;
 132     }
 133     else if (d2 < c[i].d2b)
 134     {
 135         c[i].ab  = j;
 136         c[i].d2b = d2;
 137     }
 138     /* Swap them if necessary: a must be larger than b */
 139     if (c[i].d2a < c[i].d2b)
 140     {
 141         j        = c[i].ab;
 142         c[i].ab  = c[i].aa;
 143         c[i].aa  = j;
 144         d2       = c[i].d2b;
 145         c[i].d2b = c[i].d2a;
 146         c[i].d2a = d2;
 147     }
 148 }
 149
 150 /*! \brief
 151  * Adds CONECT records for surface dots.
 152  *
 153  * \param[in] fn  PDB file to append the CONECT records to.
 154  * \param[in] n   Number of dots in `x`.
 155  * \param[in] x   Array of surface dot positions.
 156  *
 157  * Adds a CONECT record that connects each surface dot to its two nearest
 158  * neighbors.  The function is copied verbatim from the old gmx_sas.c
 159  * implementation.
 160  */
 161 void do_conect(const char* fn, int n, rvec x[])
 162 {
 163     t_conect* c = nullptr;
 164
 165     fprintf(stderr, "Building CONECT records\n");
 166     snew(c, n);
 167     for (int i = 0; (i < n); i++)
 168     {
 169         c[i].aa = c[i].ab = -1;
 170     }
 171
 172     for (int i = 0; (i < n); i++)
 173     {
 174         for (int j = i + 1; (j < n); j++)
 175         {
 176             rvec dx;
 177             rvec_sub(x[i], x[j], dx);
 178             const real d2 = iprod(dx, dx);
 179             add_rec(c, i, j, d2);
 180             add_rec(c, j, i, d2);
 181         }
 182     }
 183     FILE* fp = gmx_ffopen(fn, "a");
 184     for (int i = 0; (i < n); i++)
 185     {
 186         if ((c[i].aa == -1) || (c[i].ab == -1))
 187         {
 188             fprintf(stderr, "Warning dot %d has no connections\n", i + 1);
 189         }
 190         fprintf(fp, "CONECT%5d%5d%5d\n", i + 1, c[i].aa + 1, c[i].ab + 1);
 191     }
 192     gmx_ffclose(fp);
 193     sfree(c);
 194 }
 195
 196 /*! \brief
 197  * Plots the surface into a PDB file, optionally including the original atoms.
 198  */
 199 void connolly_plot(const char*  fn,
 200                    int          ndots,
 201                    const real   dots[],
 202                    rvec         x[],
 203                    t_atoms*     atoms,
 204                    t_symtab*    symtab,
 205                    PbcType      pbcType,
 206                    const matrix box,
 207                    gmx_bool     bIncludeSolute)
 208 {
 209     const char* const atomnm = "DOT";
 210     const char* const resnm  = "DOT";
 211     const char* const title  = "Connolly Dot Surface Generated by gmx sasa";
 212
 213     rvec* xnew = nullptr;
 214
 215     if (bIncludeSolute)
 216     {
 217         int i0 = atoms->nr;
 218         int r0 = atoms->nres;
 219         srenew(atoms->atom, atoms->nr + ndots);
 220         memset(&atoms->atom[i0], 0, sizeof(*atoms->atom) * ndots);
 221         srenew(atoms->atomname, atoms->nr + ndots);
 222         srenew(atoms->resinfo, r0 + 1);
 223         atoms->atom[i0].resind = r0;
 224         t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0 + 1, ' ', 0, ' ');
 225         if (atoms->pdbinfo != nullptr)
 226         {
 227             srenew(atoms->pdbinfo, atoms->nr + ndots);
 228         }
 229         snew(xnew, atoms->nr + ndots);
 230         for (int i = 0; (i < atoms->nr); i++)
 231         {
 232             copy_rvec(x[i], xnew[i]);
 233         }
 234         int k = 0;
 235         for (int i = 0; (i < ndots); i++)
 236         {
 237             int ii0                 = i0 + i;
 238             atoms->atomname[ii0]    = put_symtab(symtab, atomnm);
 239             atoms->atom[ii0].resind = r0;
 240             xnew[ii0][XX]           = dots[k++];
 241             xnew[ii0][YY]           = dots[k++];
 242             xnew[ii0][ZZ]           = dots[k++];
 243             if (atoms->pdbinfo != nullptr)
 244             {
 245                 atoms->pdbinfo[ii0].type   = PdbRecordType::Atom;
 246                 atoms->pdbinfo[ii0].atomnr = ii0;
 247                 atoms->pdbinfo[ii0].bfac   = 0.0;
 248                 atoms->pdbinfo[ii0].occup  = 0.0;
 249             }
 250         }
 251         atoms->nr   = i0 + ndots;
 252         atoms->nres = r0 + 1;
 253         write_sto_conf(fn, title, atoms, xnew, nullptr, pbcType, const_cast<rvec*>(box));
 254         atoms->nres = r0;
 255         atoms->nr   = i0;
 256     }
 257     else
 258     {
 259         t_atoms aaa;
 260         init_t_atoms(&aaa, ndots, TRUE);
 261         aaa.atom[0].resind = 0;
 262         t_atoms_set_resinfo(&aaa, 0, symtab, resnm, 1, ' ', 0, ' ');
 263         snew(xnew, ndots);
 264         int k = 0;
 265         for (int i = 0; (i < ndots); i++)
 266         {
 267             int ii0                 = i;
 268             aaa.atomname[ii0]       = put_symtab(symtab, atomnm);
 269             aaa.pdbinfo[ii0].type   = PdbRecordType::Atom;
 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;
 277         }
 278         aaa.nr = ndots;
 279         write_sto_conf(fn, title, &aaa, xnew, nullptr, pbcType, const_cast<rvec*>(box));
 280         do_conect(fn, ndots, xnew);
 281         done_atom(&aaa);
 282     }
 283     sfree(xnew);
 284 }
 285
 286 /********************************************************************
 287  * Actual analysis module
 288  */
 289
 290 /*! \brief
 291  * Implements `gmx sas` trajectory analysis module.
 292  */
 293 class Sasa : public TrajectoryAnalysisModule
 294 {
 295 public:
 296     Sasa();
 297
 298     void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
 299     void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
 300
 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;
 304
 305     void finishAnalysis(int nframes) override;
 306     void writeOutput() override;
 307
 308 private:
 309     /*! \brief
 310      * Surface areas as a function of time.
 311      *
 312      * First column is for the calculation group, and the rest for the
 313      * output groups.  This data is always produced.
 314      */
 315     AnalysisData area_;
 316     /*! \brief
 317      * Per-atom surface areas as a function of time.
 318      *
 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
 322      * requested.
 323      */
 324     AnalysisData atomArea_;
 325     /*! \brief
 326      * Per-residue surface areas as a function of time.
 327      *
 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
 333      * requested.
 334      */
 335     AnalysisData residueArea_;
 336     /*! \brief
 337      * Free energy estimates as a function of time.
 338      *
 339      * Column layout is the same as for `area_`.
 340      * This data is only produced if the output is requested.
 341      */
 342     AnalysisData dgSolv_;
 343     /*! \brief
 344      * Total volume and density of the calculation group as a function of
 345      * time.
 346      *
 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.
 349      */
 350     AnalysisData volume_;
 351
 352     /*! \brief
 353      * The selection to calculate the surface for.
 354      *
 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.
 360      */
 361     Selection surfaceSel_;
 362     /*! \brief
 363      * List of optional additional output groups.
 364      *
 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_`.
 368      */
 369     SelectionList outputSel_;
 370
 371     std::string fnArea_;
 372     std::string fnAtomArea_;
 373     std::string fnResidueArea_;
 374     std::string fnDGSolv_;
 375     std::string fnVolume_;
 376     std::string fnConnolly_;
 377
 378     double solsize_;
 379     int    ndots_;
 380     // double                  minarea_;
 381     double dgsDefault_;
 382     bool   bIncludeSolute_;
 383
 384     //! Global topology corresponding to the input.
 385     gmx_mtop_t* mtop_;
 386     //! Per-atom data corresponding to the input.
 387     AtomsDataPtr atoms_;
 388     //! Combined VdW and probe radii for each atom in the calculation group.
 389     std::vector<real> radii_;
 390     /*! \brief
 391      * Solvation free energy coefficients for each atom in the calculation
 392      * group.
 393      *
 394      * Empty if the free energy output has not been requested.
 395      */
 396     std::vector<real> dgsFactor_;
 397     //! Calculation algorithm.
 398     SurfaceAreaCalculator calculator_;
 399
 400     // Copy and assign disallowed by base.
 401 };
 402
 403 Sasa::Sasa() :
 404     solsize_(0.14),
 405     ndots_(24),
 406     dgsDefault_(0),
 407     bIncludeSolute_(true),
 408     mtop_(nullptr),
 409     atoms_(nullptr)
 410 {
 411     // minarea_ = 0.5;
 412     registerAnalysisDataset(&area_, "area");
 413     registerAnalysisDataset(&atomArea_, "atomarea");
 414     registerAnalysisDataset(&residueArea_, "resarea");
 415     registerAnalysisDataset(&dgSolv_, "dgsolv");
 416     registerAnalysisDataset(&volume_, "volume");
 417 }
 418
 419 void Sasa::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 420 {
 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]",
 430
 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]",
 438
 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",
 441         "[TT]-oa[tt]).[PAR]",
 442         //"In combination with the latter option an [REF].itp[ref] file can be",
 443         //"generated (option [TT]-i[tt])",
 444         //"which can be used to restrain surface atoms.[PAR]",
 445
 446         "With the [TT]-tv[tt] option the total volume and density of the",
 447         "molecule can be computed. With [TT]-pbc[tt] (the default), you",
 448         "must ensure that your molecule/surface group is not split across PBC.",
 449         "Otherwise, you will get non-sensical results.",
 450         "Please also consider whether the normal probe radius is appropriate",
 451         "in this case or whether you would rather use, e.g., 0. It is good",
 452         "to keep in mind that the results for volume and density are very",
 453         "approximate. For example, in ice Ih, one can easily fit water molecules in the",
 454         "pores which would yield a volume that is too low, and surface area and density",
 455         "that are both too high."
 456     };
 457
 458     settings->setHelpText(desc);
 459
 460     options->addOption(FileNameOption("o")
 461                                .filetype(OptionFileType::Plot)
 462                                .outputFile()
 463                                .required()
 464                                .store(&fnArea_)
 465                                .defaultBasename("area")
 466                                .description("Total area as a function of time"));
 467     options->addOption(
 468             FileNameOption("odg")
 469                     .filetype(OptionFileType::Plot)
 470                     .outputFile()
 471                     .store(&fnDGSolv_)
 472                     .defaultBasename("dgsolv")
 473                     .description("Estimated solvation free energy as a function of time"));
 474     options->addOption(FileNameOption("or")
 475                                .filetype(OptionFileType::Plot)
 476                                .outputFile()
 477                                .store(&fnResidueArea_)
 478                                .defaultBasename("resarea")
 479                                .description("Average area per residue"));
 480     options->addOption(FileNameOption("oa")
 481                                .filetype(OptionFileType::Plot)
 482                                .outputFile()
 483                                .store(&fnAtomArea_)
 484                                .defaultBasename("atomarea")
 485                                .description("Average area per atom"));
 486     options->addOption(FileNameOption("tv")
 487                                .filetype(OptionFileType::Plot)
 488                                .outputFile()
 489                                .store(&fnVolume_)
 490                                .defaultBasename("volume")
 491                                .description("Total volume and density as a function of time"));
 492     options->addOption(FileNameOption("q")
 493                                .filetype(OptionFileType::PDB)
 494                                .outputFile()
 495                                .store(&fnConnolly_)
 496                                .defaultBasename("connolly")
 497                                .description("PDB file for Connolly surface"));
 498
 499     options->addOption(
 500             DoubleOption("probe").store(&solsize_).description("Radius of the solvent probe (nm)"));
 501     options->addOption(IntegerOption("ndots").store(&ndots_).description(
 502             "Number of dots per sphere, more dots means more accuracy"));
 503     options->addOption(
 504             BooleanOption("prot").store(&bIncludeSolute_).description("Output the protein to the Connolly [REF].pdb[ref] file too"));
 505     options->addOption(
 506             DoubleOption("dgs").store(&dgsDefault_).description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
 507
 508     // Selections must select atoms for the VdW radii lookup to work.
 509     // The calculation group uses dynamicMask() so that the coordinates
 510     // match a static array of VdW radii.
 511     options->addOption(SelectionOption("surface")
 512                                .store(&surfaceSel_)
 513                                .required()
 514                                .onlySortedAtoms()
 515                                .dynamicMask()
 516                                .description("Surface calculation selection"));
 517     options->addOption(
 518             SelectionOption("output").storeVector(&outputSel_).onlySortedAtoms().multiValue().description("Output selection(s)"));
 519
 520     // Atom names etc. are required for the VdW radii lookup.
 521     settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 522 }
 523
 524 void Sasa::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
 525 {
 526     mtop_  = top.mtop();
 527     atoms_ = top.copyAtoms();
 528
 529     // bITP   = opt2bSet("-i", nfile, fnm);
 530     const bool bResAt = !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
 531     const bool bDGsol = !fnDGSolv_.empty();
 532
 533     if (solsize_ < 0)
 534     {
 535         solsize_ = 1e-3;
 536         fprintf(stderr, "Probe size too small, setting it to %g\n", solsize_);
 537     }
 538     if (ndots_ < 20)
 539     {
 540         ndots_ = 20;
 541         fprintf(stderr, "Ndots too small, setting it to %d\n", ndots_);
 542     }
 543
 544     please_cite(stderr, "Eisenhaber95");
 545
 546     if (bDGsol)
 547     {
 548         if (!top.hasFullTopology())
 549         {
 550             GMX_THROW(InconsistentInputError(
 551                     "Cannot compute Delta G of solvation without a tpr file"));
 552         }
 553         else
 554         {
 555             if (strcmp(*(atoms_->atomtype[0]), "?") == 0)
 556             {
 557                 GMX_THROW(InconsistentInputError(
 558                         "Your input tpr file is too old (does not contain atom types). Cannot not "
 559                         "compute Delta G of solvation"));
 560             }
 561             else
 562             {
 563                 printf("Free energy of solvation predictions:\n");
 564                 please_cite(stdout, "Eisenberg86a");
 565             }
 566         }
 567     }
 568
 569     // Now compute atomic radii including solvent probe size.
 570     // Also, fetch solvation free energy coefficients and
 571     // compute the residue indices that map the calculation atoms
 572     // to the columns of residueArea_.
 573     radii_.reserve(surfaceSel_.posCount());
 574     if (bDGsol)
 575     {
 576         dgsFactor_.reserve(surfaceSel_.posCount());
 577     }
 578
 579     const int resCount = surfaceSel_.initOriginalIdsToGroup(top.mtop(), INDEX_RES);
 580
 581     // TODO: Not exception-safe, but nice solution would be to have a C++
 582     // atom properties class...
 583     AtomProperties aps;
 584
 585     ArrayRef<const int> atomIndices = surfaceSel_.atomIndices();
 586     int                 ndefault    = 0;
 587     for (int i = 0; i < surfaceSel_.posCount(); i++)
 588     {
 589         const int ii     = atomIndices[i];
 590         const int resind = atoms_->atom[ii].resind;
 591         real      radius = 0;
 592         if (!aps.setAtomProperty(epropVDW, *(atoms_->resinfo[resind].name), *(atoms_->atomname[ii]), &radius))
 593         {
 594             ndefault++;
 595         }
 596         radii_.push_back(radius + solsize_);
 597         if (bDGsol)
 598         {
 599             real dgsFactor = 0;
 600             if (!aps.setAtomProperty(
 601                         epropDGsol, *(atoms_->resinfo[resind].name), *(atoms_->atomtype[ii]), &dgsFactor))
 602             {
 603                 dgsFactor = dgsDefault_;
 604             }
 605             dgsFactor_.push_back(dgsFactor);
 606         }
 607     }
 608     if (ndefault > 0)
 609     {
 610         fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
 611     }
 612
 613     // Pre-compute mapping from the output groups to the calculation group,
 614     // and store it in the selection ID map for easy lookup.
 615     for (size_t g = 0; g < outputSel_.size(); ++g)
 616     {
 617         ArrayRef<const int> outputIndices = outputSel_[g].atomIndices();
 618         for (int i = 0, j = 0; i < outputSel_[g].posCount(); ++i)
 619         {
 620             while (j < surfaceSel_.posCount() && outputIndices[i] > atomIndices[j])
 621             {
 622                 ++j;
 623             }
 624             if (j == surfaceSel_.posCount() || outputIndices[i] != atomIndices[j])
 625             {
 626                 const std::string message = formatString(
 627                         "Output selection '%s' is not a subset of "
 628                         "the surface selection (atom %d is the first "
 629                         "atom not in the surface selection)",
 630                         outputSel_[g].name(),
 631                         outputIndices[i] + 1);
 632                 GMX_THROW(InconsistentInputError(message));
 633             }
 634             outputSel_[g].setOriginalId(i, j);
 635         }
 636     }
 637
 638     calculator_.setDotCount(ndots_);
 639     calculator_.setRadii(radii_);
 640
 641     // Initialize all the output data objects and initialize the output plotters.
 642
 643     area_.setColumnCount(0, 1 + outputSel_.size());
 644     {
 645         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 646         plotm->setFileName(fnArea_);
 647         plotm->setTitle("Solvent Accessible Surface");
 648         plotm->setXAxisIsTime();
 649         plotm->setYLabel("Area (nm\\S2\\N)");
 650         plotm->appendLegend("Total");
 651         for (size_t i = 0; i < outputSel_.size(); ++i)
 652         {
 653             plotm->appendLegend(outputSel_[i].name());
 654         }
 655         area_.addModule(plotm);
 656     }
 657
 658     if (bResAt)
 659     {
 660         atomArea_.setDataSetCount(1 + outputSel_.size());
 661         residueArea_.setDataSetCount(1 + outputSel_.size());
 662         for (size_t i = 0; i <= outputSel_.size(); ++i)
 663         {
 664             atomArea_.setColumnCount(i, surfaceSel_.posCount());
 665             residueArea_.setColumnCount(i, resCount);
 666         }
 667         {
 668             AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
 669             for (int i = 0; i < surfaceSel_.posCount(); ++i)
 670             {
 671                 avem->setXAxisValue(i, surfaceSel_.position(i).atomIndices()[0] + 1);
 672             }
 673             atomArea_.addModule(avem);
 674             if (!fnAtomArea_.empty())
 675             {
 676                 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 677                 plotm->setFileName(fnAtomArea_);
 678                 plotm->setTitle("Area per atom over the trajectory");
 679                 plotm->setXLabel("Atom");
 680                 plotm->setXFormat(8, 0);
 681                 plotm->setYLabel("Area (nm\\S2\\N)");
 682                 plotm->setErrorsAsSeparateColumn(true);
 683                 plotm->appendLegend("Average (nm\\S2\\N)");
 684                 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
 685                 avem->addModule(plotm);
 686             }
 687         }
 688         {
 689             AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
 690             int                              nextRow = 0;
 691             for (int i = 0; i < surfaceSel_.posCount(); ++i)
 692             {
 693                 const int residueGroup = surfaceSel_.position(i).mappedId();
 694                 if (residueGroup >= nextRow)
 695                 {
 696                     GMX_ASSERT(residueGroup == nextRow,
 697                                "Inconsistent (non-uniformly increasing) residue grouping");
 698                     const int atomIndex    = surfaceSel_.position(i).atomIndices()[0];
 699                     const int residueIndex = atoms_->atom[atomIndex].resind;
 700                     avem->setXAxisValue(nextRow, atoms_->resinfo[residueIndex].nr);
 701                     ++nextRow;
 702                 }
 703             }
 704             residueArea_.addModule(avem);
 705             if (!fnResidueArea_.empty())
 706             {
 707                 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 708                 plotm->setFileName(fnResidueArea_);
 709                 plotm->setTitle("Area per residue over the trajectory");
 710                 plotm->setXLabel("Residue");
 711                 plotm->setXFormat(8, 0);
 712                 plotm->setYLabel("Area (nm\\S2\\N)");
 713                 plotm->setErrorsAsSeparateColumn(true);
 714                 plotm->appendLegend("Average (nm\\S2\\N)");
 715                 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
 716                 avem->addModule(plotm);
 717             }
 718         }
 719     }
 720
 721     if (!fnDGSolv_.empty())
 722     {
 723         dgSolv_.setColumnCount(0, 1 + outputSel_.size());
 724         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 725         plotm->setFileName(fnDGSolv_);
 726         plotm->setTitle("Free Energy of Solvation");
 727         plotm->setXAxisIsTime();
 728         plotm->setYLabel("D Gsolv");
 729         plotm->appendLegend("Total");
 730         for (size_t i = 0; i < outputSel_.size(); ++i)
 731         {
 732             plotm->appendLegend(outputSel_[i].name());
 733         }
 734         dgSolv_.addModule(plotm);
 735     }
 736
 737     if (!fnVolume_.empty())
 738     {
 739         volume_.setColumnCount(0, 2);
 740         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 741         plotm->setFileName(fnVolume_);
 742         plotm->setTitle("Volume and Density");
 743         plotm->setXAxisIsTime();
 744         plotm->appendLegend("Volume (nm\\S3\\N)");
 745         plotm->appendLegend("Density (g/l)");
 746         volume_.addModule(plotm);
 747     }
 748 }
 749
 750 /*! \brief
 751  * Temporary memory for use within a single-frame calculation.
 752  */
 753 class SasaModuleData : public TrajectoryAnalysisModuleData
 754 {
 755 public:
 756     /*! \brief
 757      * Reserves memory for the frame-local data.
 758      *
 759      * `residueCount` will be zero if per-residue data is not being
 760      * calculated.
 761      */
 762     SasaModuleData(TrajectoryAnalysisModule*          module,
 763                    const AnalysisDataParallelOptions& opt,
 764                    const SelectionCollection&         selections,
 765                    int                                atomCount,
 766                    int                                residueCount) :
 767         TrajectoryAnalysisModuleData(module, opt, selections)
 768     {
 769         index_.reserve(atomCount);
 770         // If the calculation group is not dynamic, pre-calculate
 771         // the index, since it is not going to change.
 772         for (int i = 0; i < atomCount; ++i)
 773         {
 774             index_.push_back(i);
 775         }
 776         atomAreas_.resize(atomCount);
 777         res_a_.resize(residueCount);
 778     }
 779
 780     void finish() override { finishDataHandles(); }
 781
 782     //! Indices of the calculation selection positions selected for the frame.
 783     std::vector<int> index_;
 784     /*! \brief
 785      * Atom areas for each calculation selection position for the frame.
 786      *
 787      * One entry for each position in the calculation group.
 788      * Values for atoms not selected are set to zero.
 789      */
 790     std::vector<real> atomAreas_;
 791     /*! \brief
 792      * Working array to accumulate areas for each residue.
 793      *
 794      * One entry for each distinct residue in the calculation group;
 795      * indices are not directly residue numbers or residue indices.
 796      *
 797      * This vector is empty if residue area calculations are not being
 798      * performed.
 799      */
 800     std::vector<real> res_a_;
 801 };
 802
 803 TrajectoryAnalysisModuleDataPointer Sasa::startFrames(const AnalysisDataParallelOptions& opt,
 804                                                       const SelectionCollection&         selections)
 805 {
 806     return TrajectoryAnalysisModuleDataPointer(new SasaModuleData(
 807             this, opt, selections, surfaceSel_.posCount(), residueArea_.columnCount(0)));
 808 }
 809
 810 /*! \brief
 811  * Helper method to compute the areas for a single selection.
 812  *
 813  * \param[in]  surfaceSel     The calculation selection.
 814  * \param[in]  sel            The selection to compute the areas for (can be
 815  *     `surfaceSel` or one of the output selections).
 816  * \param[in]  atomAreas      Atom areas for each position in `surfaceSel`.
 817  * \param[in]  dgsFactor      Free energy coefficients for each position in
 818  *     `surfaceSel`. If empty, free energies are not calculated.
 819  * \param[out] totalAreaOut Total area of `sel` (sum of atom areas it selects).
 820  * \param[out] dgsolvOut      Solvation free energy.
 821  *     Will be zero of `dgsFactor` is empty.
 822  * \param      atomAreaHandle Data handle to use for storing atom areas for `sel`.
 823  * \param      resAreaHandle  Data handle to use for storing residue areas for `sel`.
 824  * \param      resAreaWork    Work array for accumulating the residue areas.
 825  *     If empty, atom and residue areas are not calculated.
 826  *
 827  * `atomAreaHandle` and `resAreaHandle` are not used if `resAreaWork` is empty.
 828  */
 829 void computeAreas(const Selection&         surfaceSel,
 830                   const Selection&         sel,
 831                   const std::vector<real>& atomAreas,
 832                   const std::vector<real>& dgsFactor,
 833                   real*                    totalAreaOut,
 834                   real*                    dgsolvOut,
 835                   AnalysisDataHandle       atomAreaHandle,
 836                   AnalysisDataHandle       resAreaHandle,
 837                   std::vector<real>*       resAreaWork)
 838 {
 839     const bool bResAt    = !resAreaWork->empty();
 840     const bool bDGsolv   = !dgsFactor.empty();
 841     real       totalArea = 0;
 842     real       dgsolv    = 0;
 843
 844     if (bResAt)
 845     {
 846         std::fill(resAreaWork->begin(), resAreaWork->end(), 0.0_real);
 847     }
 848     for (int i = 0; i < sel.posCount(); ++i)
 849     {
 850         // Get the index of the atom in the calculation group.
 851         // For the output groups, the mapping has been precalculated in
 852         // initAnalysis().
 853         const int ii = (sel != surfaceSel ? sel.position(i).mappedId() : i);
 854         if (!surfaceSel.position(ii).selected())
 855         {
 856             // For the calculation group, skip unselected atoms.
 857             if (sel == surfaceSel)
 858             {
 859                 continue;
 860             }
 861             GMX_THROW(InconsistentInputError(
 862                     "Output selection is not a subset of the surface selection"));
 863         }
 864         // Get the internal index of the matching residue.
 865         // These have been precalculated in initAnalysis().
 866         const int  ri       = surfaceSel.position(ii).mappedId();
 867         const real atomArea = atomAreas[ii];
 868         totalArea += atomArea;
 869         if (bResAt)
 870         {
 871             atomAreaHandle.setPoint(ii, atomArea);
 872             (*resAreaWork)[ri] += atomArea;
 873         }
 874         if (bDGsolv)
 875         {
 876             dgsolv += atomArea * dgsFactor[ii];
 877         }
 878     }
 879     if (bResAt)
 880     {
 881         for (size_t i = 0; i < (*resAreaWork).size(); ++i)
 882         {
 883             resAreaHandle.setPoint(i, (*resAreaWork)[i]);
 884         }
 885     }
 886     *totalAreaOut = totalArea;
 887     *dgsolvOut    = dgsolv;
 888 }
 889
 890 void Sasa::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
 891 {
 892     AnalysisDataHandle   ah         = pdata->dataHandle(area_);
 893     AnalysisDataHandle   dgh        = pdata->dataHandle(dgSolv_);
 894     AnalysisDataHandle   aah        = pdata->dataHandle(atomArea_);
 895     AnalysisDataHandle   rah        = pdata->dataHandle(residueArea_);
 896     AnalysisDataHandle   vh         = pdata->dataHandle(volume_);
 897     const Selection&     surfaceSel = pdata->parallelSelection(surfaceSel_);
 898     const SelectionList& outputSel  = pdata->parallelSelections(outputSel_);
 899     SasaModuleData&      frameData  = *static_cast<SasaModuleData*>(pdata);
 900
 901     const bool bResAt    = !frameData.res_a_.empty();
 902     const bool bDGsol    = !dgsFactor_.empty();
 903     const bool bConnolly = (frnr == 0 && !fnConnolly_.empty());
 904
 905     // Update indices of selected atoms in the work array.
 906     if (surfaceSel.isDynamic())
 907     {
 908         frameData.index_.clear();
 909         for (int i = 0; i < surfaceSel.posCount(); ++i)
 910         {
 911             if (surfaceSel.position(i).selected())
 912             {
 913                 frameData.index_.push_back(i);
 914             }
 915         }
 916     }
 917
 918     // Determine what needs to be calculated.
 919     int flag = 0;
 920     if (bResAt || bDGsol || !outputSel.empty())
 921     {
 922         flag |= FLAG_ATOM_AREA;
 923     }
 924     if (bConnolly)
 925     {
 926         flag |= FLAG_DOTS;
 927     }
 928     if (volume_.columnCount() > 0)
 929     {
 930         flag |= FLAG_VOLUME;
 931     }
 932
 933     // Do the low-level calculation.
 934     // totarea and totvolume receive the values for the calculation group.
 935     // area array contains the per-atom areas for the selected positions.
 936     // surfacedots contains nsurfacedots entries, and contains the actual
 937     // points.
 938     real  totarea   = 0;
 939     real  totvolume = 0;
 940     real *area = nullptr, *surfacedots = nullptr;
 941     int   nsurfacedots = 0;
 942     calculator_.calculate(surfaceSel.coordinates().data(),
 943                           pbc,
 944                           frameData.index_.size(),
 945                           frameData.index_.data(),
 946                           flag,
 947                           &totarea,
 948                           &totvolume,
 949                           &area,
 950                           &surfacedots,
 951                           &nsurfacedots);
 952     // Unpack the atomwise areas into the frameData.atomAreas_ array for easier
 953     // indexing in the case of dynamic surfaceSel.
 954     if (area != nullptr)
 955     {
 956         if (surfaceSel.isDynamic())
 957         {
 958             std::fill(frameData.atomAreas_.begin(), frameData.atomAreas_.end(), 0.0_real);
 959             for (size_t i = 0; i < frameData.index_.size(); ++i)
 960             {
 961                 frameData.atomAreas_[frameData.index_[i]] = area[i];
 962             }
 963         }
 964         else
 965         {
 966             std::copy(area, area + surfaceSel.posCount(), frameData.atomAreas_.begin());
 967         }
 968         sfree(area);
 969     }
 970     const sfree_guard dotsGuard(surfacedots);
 971
 972     if (bConnolly)
 973     {
 974         if (fr.natoms != mtop_->natoms)
 975         {
 976             GMX_THROW(
 977                     InconsistentInputError("Connolly plot (-q) is only supported for trajectories "
 978                                            "that contain all the atoms"));
 979         }
 980         // This is somewhat nasty, as it modifies the atoms and symtab
 981         // structures.  But since it is only used in the first frame, and no
 982         // one else uses the topology after initialization, it may just work
 983         // even with future parallelization.
 984         connolly_plot(fnConnolly_.c_str(),
 985                       nsurfacedots,
 986                       surfacedots,
 987                       fr.x,
 988                       atoms_.get(),
 989                       &mtop_->symtab,
 990                       fr.pbcType,
 991                       fr.box,
 992                       bIncludeSolute_);
 993     }
 994
 995     ah.startFrame(frnr, fr.time);
 996     if (bResAt)
 997     {
 998         aah.startFrame(frnr, fr.time);
 999         rah.startFrame(frnr, fr.time);
1000     }
1001     if (bDGsol)
1002     {
1003         dgh.startFrame(frnr, fr.time);
1004     }
1005
1006     ah.setPoint(0, totarea);
1007
1008     real totalArea = 0;
1009     real dgsolv    = 0;
1010     if (bResAt || bDGsol)
1011     {
1012         computeAreas(surfaceSel,
1013                      surfaceSel,
1014                      frameData.atomAreas_,
1015                      dgsFactor_,
1016                      &totalArea,
1017                      &dgsolv,
1018                      aah,
1019                      rah,
1020                      &frameData.res_a_);
1021         if (bDGsol)
1022         {
1023             dgh.setPoint(0, dgsolv);
1024         }
1025     }
1026     for (size_t g = 0; g < outputSel.size(); ++g)
1027     {
1028         if (bResAt)
1029         {
1030             aah.selectDataSet(g + 1);
1031             rah.selectDataSet(g + 1);
1032         }
1033         computeAreas(surfaceSel,
1034                      outputSel[g],
1035                      frameData.atomAreas_,
1036                      dgsFactor_,
1037                      &totalArea,
1038                      &dgsolv,
1039                      aah,
1040                      rah,
1041                      &frameData.res_a_);
1042         ah.setPoint(g + 1, totalArea);
1043         if (bDGsol)
1044         {
1045             dgh.setPoint(g + 1, dgsolv);
1046         }
1047     }
1048
1049     ah.finishFrame();
1050     if (bResAt)
1051     {
1052         aah.finishFrame();
1053         rah.finishFrame();
1054     }
1055     if (bDGsol)
1056     {
1057         dgh.finishFrame();
1058     }
1059
1060     if (vh.isValid())
1061     {
1062         real totmass = 0;
1063         for (int i = 0; i < surfaceSel.posCount(); ++i)
1064         {
1065             totmass += surfaceSel.position(i).mass();
1066         }
1067         const real density = totmass * gmx::c_amu / (totvolume * gmx::c_nano * gmx::c_nano * gmx::c_nano);
1068         vh.startFrame(frnr, fr.time);
1069         vh.setPoint(0, totvolume);
1070         vh.setPoint(1, density);
1071         vh.finishFrame();
1072     }
1073 }
1074
1075 void Sasa::finishAnalysis(int /*nframes*/) {}
1076
1077 void Sasa::writeOutput() {}
1078
1079 //! \}
1080
1081 } // namespace
1082
1083 const char SasaInfo::name[]             = "sasa";
1084 const char SasaInfo::shortDescription[] = "Compute solvent accessible surface area";
1085
1086 TrajectoryAnalysisModulePointer SasaInfo::create()
1087 {
1088     return TrajectoryAnalysisModulePointer(new Sasa);
1089 }
1090
1091 } // namespace analysismodules
1092
1093 } // namespace gmx