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