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