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, 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     FILE*     fp;
 164     int       i, j;
 165     t_conect* c;
 166     rvec      dx;
 167     real      d2;
 168
 169     fprintf(stderr, "Building CONECT records\n");
 170     snew(c, n);
 171     for (i = 0; (i < n); i++)
 172     {
 173         c[i].aa = c[i].ab = -1;
 174     }
 175
 176     for (i = 0; (i < n); i++)
 177     {
 178         for (j = i + 1; (j < n); j++)
 179         {
 180             rvec_sub(x[i], x[j], dx);
 181             d2 = iprod(dx, dx);
 182             add_rec(c, i, j, d2);
 183             add_rec(c, j, i, d2);
 184         }
 185     }
 186     fp = gmx_ffopen(fn, "a");
 187     for (i = 0; (i < n); i++)
 188     {
 189         if ((c[i].aa == -1) || (c[i].ab == -1))
 190         {
 191             fprintf(stderr, "Warning dot %d has no connections\n", i + 1);
 192         }
 193         fprintf(fp, "CONECT%5d%5d%5d\n", i + 1, c[i].aa + 1, c[i].ab + 1);
 194     }
 195     gmx_ffclose(fp);
 196     sfree(c);
 197 }
 198
 199 /*! \brief
 200  * Plots the surface into a PDB file, optionally including the original atoms.
 201  */
 202 void connolly_plot(const char*  fn,
 203                    int          ndots,
 204                    const real   dots[],
 205                    rvec         x[],
 206                    t_atoms*     atoms,
 207                    t_symtab*    symtab,
 208                    PbcType      pbcType,
 209                    const matrix box,
 210                    gmx_bool     bIncludeSolute)
 211 {
 212     const char* const atomnm = "DOT";
 213     const char* const resnm  = "DOT";
 214     const char* const title  = "Connolly Dot Surface Generated by gmx sasa";
 215
 216     int     i, i0, r0, ii0, k;
 217     rvec*   xnew;
 218     t_atoms aaa;
 219
 220     if (bIncludeSolute)
 221     {
 222         i0 = atoms->nr;
 223         r0 = atoms->nres;
 224         srenew(atoms->atom, atoms->nr + ndots);
 225         memset(&atoms->atom[i0], 0, sizeof(*atoms->atom) * ndots);
 226         srenew(atoms->atomname, atoms->nr + ndots);
 227         srenew(atoms->resinfo, r0 + 1);
 228         atoms->atom[i0].resind = r0;
 229         t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0 + 1, ' ', 0, ' ');
 230         if (atoms->pdbinfo != nullptr)
 231         {
 232             srenew(atoms->pdbinfo, atoms->nr + ndots);
 233         }
 234         snew(xnew, atoms->nr + ndots);
 235         for (i = 0; (i < atoms->nr); i++)
 236         {
 237             copy_rvec(x[i], xnew[i]);
 238         }
 239         for (i = k = 0; (i < ndots); i++)
 240         {
 241             ii0                     = i0 + i;
 242             atoms->atomname[ii0]    = put_symtab(symtab, atomnm);
 243             atoms->atom[ii0].resind = r0;
 244             xnew[ii0][XX]           = dots[k++];
 245             xnew[ii0][YY]           = dots[k++];
 246             xnew[ii0][ZZ]           = dots[k++];
 247             if (atoms->pdbinfo != nullptr)
 248             {
 249                 atoms->pdbinfo[ii0].type   = epdbATOM;
 250                 atoms->pdbinfo[ii0].atomnr = ii0;
 251                 atoms->pdbinfo[ii0].bfac   = 0.0;
 252                 atoms->pdbinfo[ii0].occup  = 0.0;
 253             }
 254         }
 255         atoms->nr   = i0 + ndots;
 256         atoms->nres = r0 + 1;
 257         write_sto_conf(fn, title, atoms, xnew, nullptr, pbcType, const_cast<rvec*>(box));
 258         atoms->nres = r0;
 259         atoms->nr   = i0;
 260     }
 261     else
 262     {
 263         init_t_atoms(&aaa, ndots, TRUE);
 264         aaa.atom[0].resind = 0;
 265         t_atoms_set_resinfo(&aaa, 0, symtab, resnm, 1, ' ', 0, ' ');
 266         snew(xnew, ndots);
 267         for (i = k = 0; (i < ndots); i++)
 268         {
 269             ii0                     = i;
 270             aaa.atomname[ii0]       = put_symtab(symtab, atomnm);
 271             aaa.pdbinfo[ii0].type   = epdbATOM;
 272             aaa.pdbinfo[ii0].atomnr = ii0;
 273             aaa.atom[ii0].resind    = 0;
 274             xnew[ii0][XX]           = dots[k++];
 275             xnew[ii0][YY]           = dots[k++];
 276             xnew[ii0][ZZ]           = dots[k++];
 277             aaa.pdbinfo[ii0].bfac   = 0.0;
 278             aaa.pdbinfo[ii0].occup  = 0.0;
 279         }
 280         aaa.nr = ndots;
 281         write_sto_conf(fn, title, &aaa, xnew, nullptr, pbcType, const_cast<rvec*>(box));
 282         do_conect(fn, ndots, xnew);
 283         done_atom(&aaa);
 284     }
 285     sfree(xnew);
 286 }
 287
 288 /********************************************************************
 289  * Actual analysis module
 290  */
 291
 292 /*! \brief
 293  * Implements `gmx sas` trajectory analysis module.
 294  */
 295 class Sasa : public TrajectoryAnalysisModule
 296 {
 297 public:
 298     Sasa();
 299
 300     void initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings) override;
 301     void initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top) override;
 302
 303     TrajectoryAnalysisModuleDataPointer startFrames(const AnalysisDataParallelOptions& opt,
 304                                                     const SelectionCollection& selections) override;
 305     void                                analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, TrajectoryAnalysisModuleData* pdata) override;
 306
 307     void finishAnalysis(int nframes) override;
 308     void writeOutput() override;
 309
 310 private:
 311     /*! \brief
 312      * Surface areas as a function of time.
 313      *
 314      * First column is for the calculation group, and the rest for the
 315      * output groups.  This data is always produced.
 316      */
 317     AnalysisData area_;
 318     /*! \brief
 319      * Per-atom surface areas as a function of time.
 320      *
 321      * Contains one data set for each column in `area_`.
 322      * Each column corresponds to a selection position in `surfaceSel_`.
 323      * This data is only produced if atom or residue areas have been
 324      * requested.
 325      */
 326     AnalysisData atomArea_;
 327     /*! \brief
 328      * Per-residue surface areas as a function of time.
 329      *
 330      * Contains one data set for each column in `area_`.
 331      * Each column corresponds to a distinct residue `surfaceSel_`.
 332      * For example, if `surfaceSel_` selects residues 2, 5, and 7, there
 333      * will be three columns here.
 334      * This data is only produced if atom or residue areas have been
 335      * requested.
 336      */
 337     AnalysisData residueArea_;
 338     /*! \brief
 339      * Free energy estimates as a function of time.
 340      *
 341      * Column layout is the same as for `area_`.
 342      * This data is only produced if the output is requested.
 343      */
 344     AnalysisData dgSolv_;
 345     /*! \brief
 346      * Total volume and density of the calculation group as a function of
 347      * time.
 348      *
 349      * The first column is the volume and the second column is the density.
 350      * This data is only produced if the output is requested.
 351      */
 352     AnalysisData volume_;
 353
 354     /*! \brief
 355      * The selection to calculate the surface for.
 356      *
 357      * Selection::originalId() and Selection::mappedId() store the mapping
 358      * from the positions to the columns of `residueArea_`.
 359      * The selection is computed with SelectionOption::dynamicMask(), i.e.,
 360      * even in the presence of a dynamic selection, the number of returned
 361      * positions is fixed, and SelectionPosition::selected() is used.
 362      */
 363     Selection surfaceSel_;
 364     /*! \brief
 365      * List of optional additional output groups.
 366      *
 367      * Each of these must be a subset of the `surfaceSel_`.
 368      * Selection::originalId() and Selection::mappedId() store the mapping
 369      * from the positions to the corresponsing positions in `surfaceSel_`.
 370      */
 371     SelectionList outputSel_;
 372
 373     std::string fnArea_;
 374     std::string fnAtomArea_;
 375     std::string fnResidueArea_;
 376     std::string fnDGSolv_;
 377     std::string fnVolume_;
 378     std::string fnConnolly_;
 379
 380     double solsize_;
 381     int    ndots_;
 382     // double                  minarea_;
 383     double dgsDefault_;
 384     bool   bIncludeSolute_;
 385
 386     //! Global topology corresponding to the input.
 387     gmx_mtop_t* mtop_;
 388     //! Per-atom data corresponding to the input.
 389     AtomsDataPtr atoms_;
 390     //! Combined VdW and probe radii for each atom in the calculation group.
 391     std::vector<real> radii_;
 392     /*! \brief
 393      * Solvation free energy coefficients for each atom in the calculation
 394      * group.
 395      *
 396      * Empty if the free energy output has not been requested.
 397      */
 398     std::vector<real> dgsFactor_;
 399     //! Calculation algorithm.
 400     SurfaceAreaCalculator calculator_;
 401
 402     // Copy and assign disallowed by base.
 403 };
 404
 405 Sasa::Sasa() :
 406     solsize_(0.14),
 407     ndots_(24),
 408     dgsDefault_(0),
 409     bIncludeSolute_(true),
 410     mtop_(nullptr),
 411     atoms_(nullptr)
 412 {
 413     // minarea_ = 0.5;
 414     registerAnalysisDataset(&area_, "area");
 415     registerAnalysisDataset(&atomArea_, "atomarea");
 416     registerAnalysisDataset(&residueArea_, "resarea");
 417     registerAnalysisDataset(&dgSolv_, "dgsolv");
 418     registerAnalysisDataset(&volume_, "volume");
 419 }
 420
 421 void Sasa::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* settings)
 422 {
 423     static const char* const desc[] = {
 424         "[THISMODULE] computes solvent accessible surface areas.",
 425         "See Eisenhaber F, Lijnzaad P, Argos P, Sander C, & Scharf M",
 426         "(1995) J. Comput. Chem. 16, 273-284 for the algorithm used.",
 427         "With [TT]-q[tt], the Connolly surface can be generated as well",
 428         "in a [REF].pdb[ref] file where the nodes are represented as atoms",
 429         "and the edges connecting the nearest nodes as CONECT records.",
 430         "[TT]-odg[tt] allows for estimation of solvation free energies",
 431         "from per-atom solvation energies per exposed surface area.[PAR]",
 432
 433         "The program requires a selection for the surface calculation to be",
 434         "specified with [TT]-surface[tt]. This should always consist of all",
 435         "non-solvent atoms in the system. The area of this group is always",
 436         "calculated. Optionally, [TT]-output[tt] can specify additional",
 437         "selections, which should be subsets of the calculation group.",
 438         "The solvent-accessible areas for these groups are also extracted",
 439         "from the full surface.[PAR]",
 440
 441         "The average and standard deviation of the area over the trajectory",
 442         "can be calculated per residue and atom (options [TT]-or[tt] and", "[TT]-oa[tt]).[PAR]",
 443         //"In combination with the latter option an [REF].itp[ref] file can be",
 444         //"generated (option [TT]-i[tt])",
 445         //"which can be used to restrain surface atoms.[PAR]",
 446
 447         "With the [TT]-tv[tt] option the total volume and density of the",
 448         "molecule can be computed. With [TT]-pbc[tt] (the default), you",
 449         "must ensure that your molecule/surface group is not split across PBC.",
 450         "Otherwise, you will get non-sensical results.",
 451         "Please also consider whether the normal probe radius is appropriate",
 452         "in this case or whether you would rather use, e.g., 0. It is good",
 453         "to keep in mind that the results for volume and density are very",
 454         "approximate. For example, in ice Ih, one can easily fit water molecules in the",
 455         "pores which would yield a volume that is too low, and surface area and density",
 456         "that are both too high."
 457     };
 458
 459     settings->setHelpText(desc);
 460
 461     options->addOption(FileNameOption("o")
 462                                .filetype(eftPlot)
 463                                .outputFile()
 464                                .required()
 465                                .store(&fnArea_)
 466                                .defaultBasename("area")
 467                                .description("Total area as a function of time"));
 468     options->addOption(
 469             FileNameOption("odg")
 470                     .filetype(eftPlot)
 471                     .outputFile()
 472                     .store(&fnDGSolv_)
 473                     .defaultBasename("dgsolv")
 474                     .description("Estimated solvation free energy as a function of time"));
 475     options->addOption(FileNameOption("or")
 476                                .filetype(eftPlot)
 477                                .outputFile()
 478                                .store(&fnResidueArea_)
 479                                .defaultBasename("resarea")
 480                                .description("Average area per residue"));
 481     options->addOption(FileNameOption("oa")
 482                                .filetype(eftPlot)
 483                                .outputFile()
 484                                .store(&fnAtomArea_)
 485                                .defaultBasename("atomarea")
 486                                .description("Average area per atom"));
 487     options->addOption(FileNameOption("tv")
 488                                .filetype(eftPlot)
 489                                .outputFile()
 490                                .store(&fnVolume_)
 491                                .defaultBasename("volume")
 492                                .description("Total volume and density as a function of time"));
 493     options->addOption(FileNameOption("q")
 494                                .filetype(eftPDB)
 495                                .outputFile()
 496                                .store(&fnConnolly_)
 497                                .defaultBasename("connolly")
 498                                .description("PDB file for Connolly surface"));
 499     // options->addOption(FileNameOption("i").filetype(eftITP).outputFile()
 500     //                       .store(&fnRestraints_).defaultBasename("surfat")
 501     //                       .description("Topology file for position restraings on surface atoms"));
 502
 503
 504     options->addOption(
 505             DoubleOption("probe").store(&solsize_).description("Radius of the solvent probe (nm)"));
 506     options->addOption(IntegerOption("ndots").store(&ndots_).description(
 507             "Number of dots per sphere, more dots means more accuracy"));
 508     // options->addOption(DoubleOption("minarea").store(&minarea_)
 509     //                       .description("The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file (see help)"));
 510     options->addOption(
 511             BooleanOption("prot").store(&bIncludeSolute_).description("Output the protein to the Connolly [REF].pdb[ref] file too"));
 512     options->addOption(
 513             DoubleOption("dgs").store(&dgsDefault_).description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
 514
 515     // Selections must select atoms for the VdW radii lookup to work.
 516     // The calculation group uses dynamicMask() so that the coordinates
 517     // match a static array of VdW radii.
 518     options->addOption(SelectionOption("surface")
 519                                .store(&surfaceSel_)
 520                                .required()
 521                                .onlySortedAtoms()
 522                                .dynamicMask()
 523                                .description("Surface calculation selection"));
 524     options->addOption(
 525             SelectionOption("output").storeVector(&outputSel_).onlySortedAtoms().multiValue().description("Output selection(s)"));
 526
 527     // Atom names etc. are required for the VdW radii lookup.
 528     settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 529 }
 530
 531 void Sasa::initAnalysis(const TrajectoryAnalysisSettings& settings, const TopologyInformation& top)
 532 {
 533     mtop_  = top.mtop();
 534     atoms_ = top.copyAtoms();
 535
 536     // bITP   = opt2bSet("-i", nfile, fnm);
 537     const bool bResAt = !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
 538     const bool bDGsol = !fnDGSolv_.empty();
 539
 540     if (solsize_ < 0)
 541     {
 542         solsize_ = 1e-3;
 543         fprintf(stderr, "Probe size too small, setting it to %g\n", solsize_);
 544     }
 545     if (ndots_ < 20)
 546     {
 547         ndots_ = 20;
 548         fprintf(stderr, "Ndots too small, setting it to %d\n", ndots_);
 549     }
 550
 551     please_cite(stderr, "Eisenhaber95");
 552     // if ((top.pbcType() != PbcType::Xyz) || (TRICLINIC(fr.box)))
 553     //{
 554     //    fprintf(stderr, "\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
 555     //            "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
 556     //            "will certainly crash the analysis.\n\n");
 557     //}
 558
 559     if (bDGsol)
 560     {
 561         if (!top.hasFullTopology())
 562         {
 563             GMX_THROW(InconsistentInputError(
 564                     "Cannot compute Delta G of solvation without a tpr file"));
 565         }
 566         else
 567         {
 568             if (strcmp(*(atoms_->atomtype[0]), "?") == 0)
 569             {
 570                 GMX_THROW(InconsistentInputError(
 571                         "Your input tpr file is too old (does not contain atom types). Cannot not "
 572                         "compute Delta G of solvation"));
 573             }
 574             else
 575             {
 576                 printf("Free energy of solvation predictions:\n");
 577                 please_cite(stdout, "Eisenberg86a");
 578             }
 579         }
 580     }
 581
 582     // Now compute atomic radii including solvent probe size.
 583     // Also, fetch solvation free energy coefficients and
 584     // compute the residue indices that map the calculation atoms
 585     // to the columns of residueArea_.
 586     radii_.reserve(surfaceSel_.posCount());
 587     if (bDGsol)
 588     {
 589         dgsFactor_.reserve(surfaceSel_.posCount());
 590     }
 591
 592     const int resCount = surfaceSel_.initOriginalIdsToGroup(top.mtop(), INDEX_RES);
 593
 594     // TODO: Not exception-safe, but nice solution would be to have a C++
 595     // atom properties class...
 596     AtomProperties aps;
 597
 598     ArrayRef<const int> atomIndices = surfaceSel_.atomIndices();
 599     int                 ndefault    = 0;
 600     for (int i = 0; i < surfaceSel_.posCount(); i++)
 601     {
 602         const int ii     = atomIndices[i];
 603         const int resind = atoms_->atom[ii].resind;
 604         real      radius;
 605         if (!aps.setAtomProperty(epropVDW, *(atoms_->resinfo[resind].name), *(atoms_->atomname[ii]), &radius))
 606         {
 607             ndefault++;
 608         }
 609         radii_.push_back(radius + solsize_);
 610         if (bDGsol)
 611         {
 612             real dgsFactor;
 613             if (!aps.setAtomProperty(epropDGsol, *(atoms_->resinfo[resind].name),
 614                                      *(atoms_->atomtype[ii]), &dgsFactor))
 615             {
 616                 dgsFactor = dgsDefault_;
 617             }
 618             dgsFactor_.push_back(dgsFactor);
 619         }
 620     }
 621     if (ndefault > 0)
 622     {
 623         fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
 624     }
 625
 626     // Pre-compute mapping from the output groups to the calculation group,
 627     // and store it in the selection ID map for easy lookup.
 628     for (size_t g = 0; g < outputSel_.size(); ++g)
 629     {
 630         ArrayRef<const int> outputIndices = outputSel_[g].atomIndices();
 631         for (int i = 0, j = 0; i < outputSel_[g].posCount(); ++i)
 632         {
 633             while (j < surfaceSel_.posCount() && outputIndices[i] > atomIndices[j])
 634             {
 635                 ++j;
 636             }
 637             if (j == surfaceSel_.posCount() || outputIndices[i] != atomIndices[j])
 638             {
 639                 const std::string message = formatString(
 640                         "Output selection '%s' is not a subset of "
 641                         "the surface selection (atom %d is the first "
 642                         "atom not in the surface selection)",
 643                         outputSel_[g].name(), 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, totvolume;
 951     real *area = nullptr, *surfacedots = nullptr;
 952     int   nsurfacedots;
 953     calculator_.calculate(surfaceSel.coordinates().data(), pbc, frameData.index_.size(),
 954                           frameData.index_.data(), flag, &totarea, &totvolume, &area, &surfacedots,
 955                           &nsurfacedots);
 956     // Unpack the atomwise areas into the frameData.atomAreas_ array for easier
 957     // indexing in the case of dynamic surfaceSel.
 958     if (area != nullptr)
 959     {
 960         if (surfaceSel.isDynamic())
 961         {
 962             std::fill(frameData.atomAreas_.begin(), frameData.atomAreas_.end(), 0.0_real);
 963             for (size_t i = 0; i < frameData.index_.size(); ++i)
 964             {
 965                 frameData.atomAreas_[frameData.index_[i]] = area[i];
 966             }
 967         }
 968         else
 969         {
 970             std::copy(area, area + surfaceSel.posCount(), frameData.atomAreas_.begin());
 971         }
 972         sfree(area);
 973     }
 974     const sfree_guard dotsGuard(surfacedots);
 975
 976     if (bConnolly)
 977     {
 978         if (fr.natoms != mtop_->natoms)
 979         {
 980             GMX_THROW(
 981                     InconsistentInputError("Connolly plot (-q) is only supported for trajectories "
 982                                            "that contain all the atoms"));
 983         }
 984         // This is somewhat nasty, as it modifies the atoms and symtab
 985         // structures.  But since it is only used in the first frame, and no
 986         // one else uses the topology after initialization, it may just work
 987         // even with future parallelization.
 988         connolly_plot(fnConnolly_.c_str(), nsurfacedots, surfacedots, fr.x, atoms_.get(),
 989                       &mtop_->symtab, fr.pbcType, fr.box, bIncludeSolute_);
 990     }
 991
 992     ah.startFrame(frnr, fr.time);
 993     if (bResAt)
 994     {
 995         aah.startFrame(frnr, fr.time);
 996         rah.startFrame(frnr, fr.time);
 997     }
 998     if (bDGsol)
 999     {
1000         dgh.startFrame(frnr, fr.time);
1001     }
1002
1003     ah.setPoint(0, totarea);
1004
1005     real totalArea, dgsolv;
1006     if (bResAt || bDGsol)
1007     {
1008         computeAreas(surfaceSel, surfaceSel, frameData.atomAreas_, dgsFactor_, &totalArea, &dgsolv,
1009                      aah, rah, &frameData.res_a_);
1010         if (bDGsol)
1011         {
1012             dgh.setPoint(0, dgsolv);
1013         }
1014     }
1015     for (size_t g = 0; g < outputSel.size(); ++g)
1016     {
1017         if (bResAt)
1018         {
1019             aah.selectDataSet(g + 1);
1020             rah.selectDataSet(g + 1);
1021         }
1022         computeAreas(surfaceSel, outputSel[g], frameData.atomAreas_, dgsFactor_, &totalArea,
1023                      &dgsolv, aah, rah, &frameData.res_a_);
1024         ah.setPoint(g + 1, totalArea);
1025         if (bDGsol)
1026         {
1027             dgh.setPoint(g + 1, dgsolv);
1028         }
1029     }
1030
1031     ah.finishFrame();
1032     if (bResAt)
1033     {
1034         aah.finishFrame();
1035         rah.finishFrame();
1036     }
1037     if (bDGsol)
1038     {
1039         dgh.finishFrame();
1040     }
1041
1042     if (vh.isValid())
1043     {
1044         real totmass = 0;
1045         for (int i = 0; i < surfaceSel.posCount(); ++i)
1046         {
1047             totmass += surfaceSel.position(i).mass();
1048         }
1049         const real density = totmass * AMU / (totvolume * NANO * NANO * NANO);
1050         vh.startFrame(frnr, fr.time);
1051         vh.setPoint(0, totvolume);
1052         vh.setPoint(1, density);
1053         vh.finishFrame();
1054     }
1055 }
1056
1057 void Sasa::finishAnalysis(int /*nframes*/)
1058 {
1059     // if (bITP)
1060     //{
1061     //    fp3 = ftp2FILE(efITP, nfile, fnm, "w");
1062     //    fprintf(fp3, "[ position_restraints ]\n"
1063     //            "#define FCX 1000\n"
1064     //            "#define FCY 1000\n"
1065     //            "#define FCZ 1000\n"
1066     //            "; Atom  Type  fx   fy   fz\n");
1067     //    for (i = 0; i < nx[0]; i++)
1068     //    {
1069     //        if (atom_area[i] > minarea)
1070     //        {
1071     //            fprintf(fp3, "%5d   1     FCX  FCX  FCZ\n", ii+1);
1072     //        }
1073     //    }
1074     //    ffclose(fp3);
1075     //}
1076 }
1077
1078 void Sasa::writeOutput() {}
1079
1080 //! \}
1081
1082 } // namespace
1083
1084 const char SasaInfo::name[]             = "sasa";
1085 const char SasaInfo::shortDescription[] = "Compute solvent accessible surface area";
1086
1087 TrajectoryAnalysisModulePointer SasaInfo::create()
1088 {
1089     return TrajectoryAnalysisModulePointer(new Sasa);
1090 }
1091
1092 } // namespace analysismodules
1093
1094 } // namespace gmx