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