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), ndots_(24), dgsDefault_(0), bIncludeSolute_(true), mtop_(nullptr), atoms_(nullptr)
 405 {
 406     // minarea_ = 0.5;
 407     registerAnalysisDataset(&area_, "area");
 408     registerAnalysisDataset(&atomArea_, "atomarea");
 409     registerAnalysisDataset(&residueArea_, "resarea");
 410     registerAnalysisDataset(&dgSolv_, "dgsolv");
 411     registerAnalysisDataset(&volume_, "volume");
 412 }
 413
 414 void Sasa::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 415 {
 416     static const char* const desc[] = {
 417         "[THISMODULE] computes solvent accessible surface areas.",
 418         "See Eisenhaber F, Lijnzaad P, Argos P, Sander C, & Scharf M",
 419         "(1995) J. Comput. Chem. 16, 273-284 for the algorithm used.",
 420         "With [TT]-q[tt], the Connolly surface can be generated as well",
 421         "in a [REF].pdb[ref] file where the nodes are represented as atoms",
 422         "and the edges connecting the nearest nodes as CONECT records.",
 423         "[TT]-odg[tt] allows for estimation of solvation free energies",
 424         "from per-atom solvation energies per exposed surface area.[PAR]",
 425
 426         "The program requires a selection for the surface calculation to be",
 427         "specified with [TT]-surface[tt]. This should always consist of all",
 428         "non-solvent atoms in the system. The area of this group is always",
 429         "calculated. Optionally, [TT]-output[tt] can specify additional",
 430         "selections, which should be subsets of the calculation group.",
 431         "The solvent-accessible areas for these groups are also extracted",
 432         "from the full surface.[PAR]",
 433
 434         "The average and standard deviation of the area over the trajectory",
 435         "can be calculated per residue and atom (options [TT]-or[tt] and",
 436         "[TT]-oa[tt]).[PAR]",
 437         //"In combination with the latter option an [REF].itp[ref] file can be",
 438         //"generated (option [TT]-i[tt])",
 439         //"which can be used to restrain surface atoms.[PAR]",
 440
 441         "With the [TT]-tv[tt] option the total volume and density of the",
 442         "molecule can be computed. With [TT]-pbc[tt] (the default), you",
 443         "must ensure that your molecule/surface group is not split across PBC.",
 444         "Otherwise, you will get non-sensical results.",
 445         "Please also consider whether the normal probe radius is appropriate",
 446         "in this case or whether you would rather use, e.g., 0. It is good",
 447         "to keep in mind that the results for volume and density are very",
 448         "approximate. For example, in ice Ih, one can easily fit water molecules in the",
 449         "pores which would yield a volume that is too low, and surface area and density",
 450         "that are both too high."
 451     };
 452
 453     settings->setHelpText(desc);
 454
 455     options->addOption(FileNameOption("o")
 456                                .filetype(OptionFileType::Plot)
 457                                .outputFile()
 458                                .required()
 459                                .store(&fnArea_)
 460                                .defaultBasename("area")
 461                                .description("Total area as a function of time"));
 462     options->addOption(
 463             FileNameOption("odg")
 464                     .filetype(OptionFileType::Plot)
 465                     .outputFile()
 466                     .store(&fnDGSolv_)
 467                     .defaultBasename("dgsolv")
 468                     .description("Estimated solvation free energy as a function of time"));
 469     options->addOption(FileNameOption("or")
 470                                .filetype(OptionFileType::Plot)
 471                                .outputFile()
 472                                .store(&fnResidueArea_)
 473                                .defaultBasename("resarea")
 474                                .description("Average area per residue"));
 475     options->addOption(FileNameOption("oa")
 476                                .filetype(OptionFileType::Plot)
 477                                .outputFile()
 478                                .store(&fnAtomArea_)
 479                                .defaultBasename("atomarea")
 480                                .description("Average area per atom"));
 481     options->addOption(FileNameOption("tv")
 482                                .filetype(OptionFileType::Plot)
 483                                .outputFile()
 484                                .store(&fnVolume_)
 485                                .defaultBasename("volume")
 486                                .description("Total volume and density as a function of time"));
 487     options->addOption(FileNameOption("q")
 488                                .filetype(OptionFileType::PDB)
 489                                .outputFile()
 490                                .store(&fnConnolly_)
 491                                .defaultBasename("connolly")
 492                                .description("PDB file for Connolly surface"));
 493
 494     options->addOption(
 495             DoubleOption("probe").store(&solsize_).description("Radius of the solvent probe (nm)"));
 496     options->addOption(IntegerOption("ndots").store(&ndots_).description(
 497             "Number of dots per sphere, more dots means more accuracy"));
 498     options->addOption(
 499             BooleanOption("prot").store(&bIncludeSolute_).description("Output the protein to the Connolly [REF].pdb[ref] file too"));
 500     options->addOption(
 501             DoubleOption("dgs").store(&dgsDefault_).description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
 502
 503     // Selections must select atoms for the VdW radii lookup to work.
 504     // The calculation group uses dynamicMask() so that the coordinates
 505     // match a static array of VdW radii.
 506     options->addOption(SelectionOption("surface")
 507                                .store(&surfaceSel_)
 508                                .required()
 509                                .onlySortedAtoms()
 510                                .dynamicMask()
 511                                .description("Surface calculation selection"));
 512     options->addOption(
 513             SelectionOption("output").storeVector(&outputSel_).onlySortedAtoms().multiValue().description("Output selection(s)"));
 514
 515     // Atom names etc. are required for the VdW radii lookup.
 516     settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 517 }
 518
 519 void Sasa::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
 520 {
 521     mtop_  = top.mtop();
 522     atoms_ = top.copyAtoms();
 523
 524     // bITP   = opt2bSet("-i", nfile, fnm);
 525     const bool bResAt = !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
 526     const bool bDGsol = !fnDGSolv_.empty();
 527
 528     if (solsize_ < 0)
 529     {
 530         solsize_ = 1e-3;
 531         fprintf(stderr, "Probe size too small, setting it to %g\n", solsize_);
 532     }
 533     if (ndots_ < 20)
 534     {
 535         ndots_ = 20;
 536         fprintf(stderr, "Ndots too small, setting it to %d\n", ndots_);
 537     }
 538
 539     please_cite(stderr, "Eisenhaber95");
 540
 541     if (bDGsol)
 542     {
 543         if (!top.hasFullTopology())
 544         {
 545             GMX_THROW(InconsistentInputError(
 546                     "Cannot compute Delta G of solvation without a tpr file"));
 547         }
 548         else
 549         {
 550             if (strcmp(*(atoms_->atomtype[0]), "?") == 0)
 551             {
 552                 GMX_THROW(InconsistentInputError(
 553                         "Your input tpr file is too old (does not contain atom types). Cannot not "
 554                         "compute Delta G of solvation"));
 555             }
 556             else
 557             {
 558                 printf("Free energy of solvation predictions:\n");
 559                 please_cite(stdout, "Eisenberg86a");
 560             }
 561         }
 562     }
 563
 564     // Now compute atomic radii including solvent probe size.
 565     // Also, fetch solvation free energy coefficients and
 566     // compute the residue indices that map the calculation atoms
 567     // to the columns of residueArea_.
 568     radii_.reserve(surfaceSel_.posCount());
 569     if (bDGsol)
 570     {
 571         dgsFactor_.reserve(surfaceSel_.posCount());
 572     }
 573
 574     const int resCount = surfaceSel_.initOriginalIdsToGroup(top.mtop(), INDEX_RES);
 575
 576     // TODO: Not exception-safe, but nice solution would be to have a C++
 577     // atom properties class...
 578     AtomProperties aps;
 579
 580     ArrayRef<const int> atomIndices = surfaceSel_.atomIndices();
 581     int                 ndefault    = 0;
 582     for (int i = 0; i < surfaceSel_.posCount(); i++)
 583     {
 584         const int ii     = atomIndices[i];
 585         const int resind = atoms_->atom[ii].resind;
 586         real      radius = 0;
 587         if (!aps.setAtomProperty(epropVDW, *(atoms_->resinfo[resind].name), *(atoms_->atomname[ii]), &radius))
 588         {
 589             ndefault++;
 590         }
 591         radii_.push_back(radius + solsize_);
 592         if (bDGsol)
 593         {
 594             real dgsFactor = 0;
 595             if (!aps.setAtomProperty(
 596                         epropDGsol, *(atoms_->resinfo[resind].name), *(atoms_->atomtype[ii]), &dgsFactor))
 597             {
 598                 dgsFactor = dgsDefault_;
 599             }
 600             dgsFactor_.push_back(dgsFactor);
 601         }
 602     }
 603     if (ndefault > 0)
 604     {
 605         fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
 606     }
 607
 608     // Pre-compute mapping from the output groups to the calculation group,
 609     // and store it in the selection ID map for easy lookup.
 610     for (size_t g = 0; g < outputSel_.size(); ++g)
 611     {
 612         ArrayRef<const int> outputIndices = outputSel_[g].atomIndices();
 613         for (int i = 0, j = 0; i < outputSel_[g].posCount(); ++i)
 614         {
 615             while (j < surfaceSel_.posCount() && outputIndices[i] > atomIndices[j])
 616             {
 617                 ++j;
 618             }
 619             if (j == surfaceSel_.posCount() || outputIndices[i] != atomIndices[j])
 620             {
 621                 const std::string message = formatString(
 622                         "Output selection '%s' is not a subset of "
 623                         "the surface selection (atom %d is the first "
 624                         "atom not in the surface selection)",
 625                         outputSel_[g].name(),
 626                         outputIndices[i] + 1);
 627                 GMX_THROW(InconsistentInputError(message));
 628             }
 629             outputSel_[g].setOriginalId(i, j);
 630         }
 631     }
 632
 633     calculator_.setDotCount(ndots_);
 634     calculator_.setRadii(radii_);
 635
 636     // Initialize all the output data objects and initialize the output plotters.
 637
 638     area_.setColumnCount(0, 1 + outputSel_.size());
 639     {
 640         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 641         plotm->setFileName(fnArea_);
 642         plotm->setTitle("Solvent Accessible Surface");
 643         plotm->setXAxisIsTime();
 644         plotm->setYLabel("Area (nm\\S2\\N)");
 645         plotm->appendLegend("Total");
 646         for (size_t i = 0; i < outputSel_.size(); ++i)
 647         {
 648             plotm->appendLegend(outputSel_[i].name());
 649         }
 650         area_.addModule(plotm);
 651     }
 652
 653     if (bResAt)
 654     {
 655         atomArea_.setDataSetCount(1 + outputSel_.size());
 656         residueArea_.setDataSetCount(1 + outputSel_.size());
 657         for (size_t i = 0; i <= outputSel_.size(); ++i)
 658         {
 659             atomArea_.setColumnCount(i, surfaceSel_.posCount());
 660             residueArea_.setColumnCount(i, resCount);
 661         }
 662         {
 663             AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
 664             for (int i = 0; i < surfaceSel_.posCount(); ++i)
 665             {
 666                 avem->setXAxisValue(i, surfaceSel_.position(i).atomIndices()[0] + 1);
 667             }
 668             atomArea_.addModule(avem);
 669             if (!fnAtomArea_.empty())
 670             {
 671                 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 672                 plotm->setFileName(fnAtomArea_);
 673                 plotm->setTitle("Area per atom over the trajectory");
 674                 plotm->setXLabel("Atom");
 675                 plotm->setXFormat(8, 0);
 676                 plotm->setYLabel("Area (nm\\S2\\N)");
 677                 plotm->setErrorsAsSeparateColumn(true);
 678                 plotm->appendLegend("Average (nm\\S2\\N)");
 679                 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
 680                 avem->addModule(plotm);
 681             }
 682         }
 683         {
 684             AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
 685             int                              nextRow = 0;
 686             for (int i = 0; i < surfaceSel_.posCount(); ++i)
 687             {
 688                 const int residueGroup = surfaceSel_.position(i).mappedId();
 689                 if (residueGroup >= nextRow)
 690                 {
 691                     GMX_ASSERT(residueGroup == nextRow,
 692                                "Inconsistent (non-uniformly increasing) residue grouping");
 693                     const int atomIndex    = surfaceSel_.position(i).atomIndices()[0];
 694                     const int residueIndex = atoms_->atom[atomIndex].resind;
 695                     avem->setXAxisValue(nextRow, atoms_->resinfo[residueIndex].nr);
 696                     ++nextRow;
 697                 }
 698             }
 699             residueArea_.addModule(avem);
 700             if (!fnResidueArea_.empty())
 701             {
 702                 AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 703                 plotm->setFileName(fnResidueArea_);
 704                 plotm->setTitle("Area per residue over the trajectory");
 705                 plotm->setXLabel("Residue");
 706                 plotm->setXFormat(8, 0);
 707                 plotm->setYLabel("Area (nm\\S2\\N)");
 708                 plotm->setErrorsAsSeparateColumn(true);
 709                 plotm->appendLegend("Average (nm\\S2\\N)");
 710                 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
 711                 avem->addModule(plotm);
 712             }
 713         }
 714     }
 715
 716     if (!fnDGSolv_.empty())
 717     {
 718         dgSolv_.setColumnCount(0, 1 + outputSel_.size());
 719         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 720         plotm->setFileName(fnDGSolv_);
 721         plotm->setTitle("Free Energy of Solvation");
 722         plotm->setXAxisIsTime();
 723         plotm->setYLabel("D Gsolv");
 724         plotm->appendLegend("Total");
 725         for (size_t i = 0; i < outputSel_.size(); ++i)
 726         {
 727             plotm->appendLegend(outputSel_[i].name());
 728         }
 729         dgSolv_.addModule(plotm);
 730     }
 731
 732     if (!fnVolume_.empty())
 733     {
 734         volume_.setColumnCount(0, 2);
 735         AnalysisDataPlotModulePointer plotm(new AnalysisDataPlotModule(settings.plotSettings()));
 736         plotm->setFileName(fnVolume_);
 737         plotm->setTitle("Volume and Density");
 738         plotm->setXAxisIsTime();
 739         plotm->appendLegend("Volume (nm\\S3\\N)");
 740         plotm->appendLegend("Density (g/l)");
 741         volume_.addModule(plotm);
 742     }
 743 }
 744
 745 /*! \brief
 746  * Temporary memory for use within a single-frame calculation.
 747  */
 748 class SasaModuleData : public TrajectoryAnalysisModuleData
 749 {
 750 public:
 751     /*! \brief
 752      * Reserves memory for the frame-local data.
 753      *
 754      * `residueCount` will be zero if per-residue data is not being
 755      * calculated.
 756      */
 757     SasaModuleData(TrajectoryAnalysisModule*          module,
 758                    const AnalysisDataParallelOptions& opt,
 759                    const SelectionCollection&         selections,
 760                    int                                atomCount,
 761                    int                                residueCount) :
 762         TrajectoryAnalysisModuleData(module, opt, selections)
 763     {
 764         index_.reserve(atomCount);
 765         // If the calculation group is not dynamic, pre-calculate
 766         // the index, since it is not going to change.
 767         for (int i = 0; i < atomCount; ++i)
 768         {
 769             index_.push_back(i);
 770         }
 771         atomAreas_.resize(atomCount);
 772         res_a_.resize(residueCount);
 773     }
 774
 775     void finish() override { finishDataHandles(); }
 776
 777     //! Indices of the calculation selection positions selected for the frame.
 778     std::vector<int> index_;
 779     /*! \brief
 780      * Atom areas for each calculation selection position for the frame.
 781      *
 782      * One entry for each position in the calculation group.
 783      * Values for atoms not selected are set to zero.
 784      */
 785     std::vector<real> atomAreas_;
 786     /*! \brief
 787      * Working array to accumulate areas for each residue.
 788      *
 789      * One entry for each distinct residue in the calculation group;
 790      * indices are not directly residue numbers or residue indices.
 791      *
 792      * This vector is empty if residue area calculations are not being
 793      * performed.
 794      */
 795     std::vector<real> res_a_;
 796 };
 797
 798 TrajectoryAnalysisModuleDataPointer Sasa::startFrames(const AnalysisDataParallelOptions& opt,
 799                                                       const SelectionCollection&         selections)
 800 {
 801     return TrajectoryAnalysisModuleDataPointer(new SasaModuleData(
 802             this, opt, selections, surfaceSel_.posCount(), residueArea_.columnCount(0)));
 803 }
 804
 805 /*! \brief
 806  * Helper method to compute the areas for a single selection.
 807  *
 808  * \param[in]  surfaceSel     The calculation selection.
 809  * \param[in]  sel            The selection to compute the areas for (can be
 810  *     `surfaceSel` or one of the output selections).
 811  * \param[in]  atomAreas      Atom areas for each position in `surfaceSel`.
 812  * \param[in]  dgsFactor      Free energy coefficients for each position in
 813  *     `surfaceSel`. If empty, free energies are not calculated.
 814  * \param[out] totalAreaOut Total area of `sel` (sum of atom areas it selects).
 815  * \param[out] dgsolvOut      Solvation free energy.
 816  *     Will be zero of `dgsFactor` is empty.
 817  * \param      atomAreaHandle Data handle to use for storing atom areas for `sel`.
 818  * \param      resAreaHandle  Data handle to use for storing residue areas for `sel`.
 819  * \param      resAreaWork    Work array for accumulating the residue areas.
 820  *     If empty, atom and residue areas are not calculated.
 821  *
 822  * `atomAreaHandle` and `resAreaHandle` are not used if `resAreaWork` is empty.
 823  */
 824 void computeAreas(const Selection&         surfaceSel,
 825                   const Selection&         sel,
 826                   const std::vector<real>& atomAreas,
 827                   const std::vector<real>& dgsFactor,
 828                   real*                    totalAreaOut,
 829                   real*                    dgsolvOut,
 830                   AnalysisDataHandle       atomAreaHandle,
 831                   AnalysisDataHandle       resAreaHandle,
 832                   std::vector<real>*       resAreaWork)
 833 {
 834     const bool bResAt    = !resAreaWork->empty();
 835     const bool bDGsolv   = !dgsFactor.empty();
 836     real       totalArea = 0;
 837     real       dgsolv    = 0;
 838
 839     if (bResAt)
 840     {
 841         std::fill(resAreaWork->begin(), resAreaWork->end(), 0.0_real);
 842     }
 843     for (int i = 0; i < sel.posCount(); ++i)
 844     {
 845         // Get the index of the atom in the calculation group.
 846         // For the output groups, the mapping has been precalculated in
 847         // initAnalysis().
 848         const int ii = (sel != surfaceSel ? sel.position(i).mappedId() : i);
 849         if (!surfaceSel.position(ii).selected())
 850         {
 851             // For the calculation group, skip unselected atoms.
 852             if (sel == surfaceSel)
 853             {
 854                 continue;
 855             }
 856             GMX_THROW(InconsistentInputError(
 857                     "Output selection is not a subset of the surface selection"));
 858         }
 859         // Get the internal index of the matching residue.
 860         // These have been precalculated in initAnalysis().
 861         const int  ri       = surfaceSel.position(ii).mappedId();
 862         const real atomArea = atomAreas[ii];
 863         totalArea += atomArea;
 864         if (bResAt)
 865         {
 866             atomAreaHandle.setPoint(ii, atomArea);
 867             (*resAreaWork)[ri] += atomArea;
 868         }
 869         if (bDGsolv)
 870         {
 871             dgsolv += atomArea * dgsFactor[ii];
 872         }
 873     }
 874     if (bResAt)
 875     {
 876         for (size_t i = 0; i < (*resAreaWork).size(); ++i)
 877         {
 878             resAreaHandle.setPoint(i, (*resAreaWork)[i]);
 879         }
 880     }
 881     *totalAreaOut = totalArea;
 882     *dgsolvOut    = dgsolv;
 883 }
 884
 885 void Sasa::analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata)
 886 {
 887     AnalysisDataHandle   ah         = pdata->dataHandle(area_);
 888     AnalysisDataHandle   dgh        = pdata->dataHandle(dgSolv_);
 889     AnalysisDataHandle   aah        = pdata->dataHandle(atomArea_);
 890     AnalysisDataHandle   rah        = pdata->dataHandle(residueArea_);
 891     AnalysisDataHandle   vh         = pdata->dataHandle(volume_);
 892     const Selection&     surfaceSel = pdata->parallelSelection(surfaceSel_);
 893     const SelectionList& outputSel  = pdata->parallelSelections(outputSel_);
 894     SasaModuleData&      frameData  = *static_cast<SasaModuleData*>(pdata);
 895
 896     const bool bResAt    = !frameData.res_a_.empty();
 897     const bool bDGsol    = !dgsFactor_.empty();
 898     const bool bConnolly = (frnr == 0 && !fnConnolly_.empty());
 899
 900     // Update indices of selected atoms in the work array.
 901     if (surfaceSel.isDynamic())
 902     {
 903         frameData.index_.clear();
 904         for (int i = 0; i < surfaceSel.posCount(); ++i)
 905         {
 906             if (surfaceSel.position(i).selected())
 907             {
 908                 frameData.index_.push_back(i);
 909             }
 910         }
 911     }
 912
 913     // Determine what needs to be calculated.
 914     int flag = 0;
 915     if (bResAt || bDGsol || !outputSel.empty())
 916     {
 917         flag |= FLAG_ATOM_AREA;
 918     }
 919     if (bConnolly)
 920     {
 921         flag |= FLAG_DOTS;
 922     }
 923     if (volume_.columnCount() > 0)
 924     {
 925         flag |= FLAG_VOLUME;
 926     }
 927
 928     // Do the low-level calculation.
 929     // totarea and totvolume receive the values for the calculation group.
 930     // area array contains the per-atom areas for the selected positions.
 931     // surfacedots contains nsurfacedots entries, and contains the actual
 932     // points.
 933     real  totarea   = 0;
 934     real  totvolume = 0;
 935     real *area = nullptr, *surfacedots = nullptr;
 936     int   nsurfacedots = 0;
 937     calculator_.calculate(surfaceSel.coordinates().data(),
 938                           pbc,
 939                           frameData.index_.size(),
 940                           frameData.index_.data(),
 941                           flag,
 942                           &totarea,
 943                           &totvolume,
 944                           &area,
 945                           &surfacedots,
 946                           &nsurfacedots);
 947     // Unpack the atomwise areas into the frameData.atomAreas_ array for easier
 948     // indexing in the case of dynamic surfaceSel.
 949     if (area != nullptr)
 950     {
 951         if (surfaceSel.isDynamic())
 952         {
 953             std::fill(frameData.atomAreas_.begin(), frameData.atomAreas_.end(), 0.0_real);
 954             for (size_t i = 0; i < frameData.index_.size(); ++i)
 955             {
 956                 frameData.atomAreas_[frameData.index_[i]] = area[i];
 957             }
 958         }
 959         else
 960         {
 961             std::copy(area, area + surfaceSel.posCount(), frameData.atomAreas_.begin());
 962         }
 963         sfree(area);
 964     }
 965     const sfree_guard dotsGuard(surfacedots);
 966
 967     if (bConnolly)
 968     {
 969         if (fr.natoms != mtop_->natoms)
 970         {
 971             GMX_THROW(
 972                     InconsistentInputError("Connolly plot (-q) is only supported for trajectories "
 973                                            "that contain all the atoms"));
 974         }
 975         // This is somewhat nasty, as it modifies the atoms and symtab
 976         // structures.  But since it is only used in the first frame, and no
 977         // one else uses the topology after initialization, it may just work
 978         // even with future parallelization.
 979         connolly_plot(fnConnolly_.c_str(),
 980                       nsurfacedots,
 981                       surfacedots,
 982                       fr.x,
 983                       atoms_.get(),
 984                       &mtop_->symtab,
 985                       fr.pbcType,
 986                       fr.box,
 987                       bIncludeSolute_);
 988     }
 989
 990     ah.startFrame(frnr, fr.time);
 991     if (bResAt)
 992     {
 993         aah.startFrame(frnr, fr.time);
 994         rah.startFrame(frnr, fr.time);
 995     }
 996     if (bDGsol)
 997     {
 998         dgh.startFrame(frnr, fr.time);
 999     }
1000
1001     ah.setPoint(0, totarea);
1002
1003     real totalArea = 0;
1004     real dgsolv    = 0;
1005     if (bResAt || bDGsol)
1006     {
1007         computeAreas(surfaceSel,
1008                      surfaceSel,
1009                      frameData.atomAreas_,
1010                      dgsFactor_,
1011                      &totalArea,
1012                      &dgsolv,
1013                      aah,
1014                      rah,
1015                      &frameData.res_a_);
1016         if (bDGsol)
1017         {
1018             dgh.setPoint(0, dgsolv);
1019         }
1020     }
1021     for (size_t g = 0; g < outputSel.size(); ++g)
1022     {
1023         if (bResAt)
1024         {
1025             aah.selectDataSet(g + 1);
1026             rah.selectDataSet(g + 1);
1027         }
1028         computeAreas(surfaceSel,
1029                      outputSel[g],
1030                      frameData.atomAreas_,
1031                      dgsFactor_,
1032                      &totalArea,
1033                      &dgsolv,
1034                      aah,
1035                      rah,
1036                      &frameData.res_a_);
1037         ah.setPoint(g + 1, totalArea);
1038         if (bDGsol)
1039         {
1040             dgh.setPoint(g + 1, dgsolv);
1041         }
1042     }
1043
1044     ah.finishFrame();
1045     if (bResAt)
1046     {
1047         aah.finishFrame();
1048         rah.finishFrame();
1049     }
1050     if (bDGsol)
1051     {
1052         dgh.finishFrame();
1053     }
1054
1055     if (vh.isValid())
1056     {
1057         real totmass = 0;
1058         for (int i = 0; i < surfaceSel.posCount(); ++i)
1059         {
1060             totmass += surfaceSel.position(i).mass();
1061         }
1062         const real density = totmass * gmx::c_amu / (totvolume * gmx::c_nano * gmx::c_nano * gmx::c_nano);
1063         vh.startFrame(frnr, fr.time);
1064         vh.setPoint(0, totvolume);
1065         vh.setPoint(1, density);
1066         vh.finishFrame();
1067     }
1068 }
1069
1070 void Sasa::finishAnalysis(int /*nframes*/) {}
1071
1072 void Sasa::writeOutput() {}
1073
1074 //! \}
1075
1076 } // namespace
1077
1078 const char SasaInfo::name[]             = "sasa";
1079 const char SasaInfo::shortDescription[] = "Compute solvent accessible surface area";
1080
1081 TrajectoryAnalysisModulePointer SasaInfo::create()
1082 {
1083     return TrajectoryAnalysisModulePointer(new Sasa);
1084 }
1085
1086 } // namespace analysismodules
1087
1088 } // namespace gmx