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, 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, int ndots, const real dots[], rvec x[], t_atoms *atoms,
 201                    t_symtab *symtab, int ePBC, const matrix box, gmx_bool bIncludeSolute)
 202 {
 203     const char *const  atomnm = "DOT";
 204     const char *const  resnm  = "DOT";
 205     const char *const  title  = "Connolly Dot Surface Generated by gmx sasa";
 206
 207     int                i, i0, r0, ii0, k;
 208     rvec              *xnew;
 209     t_atoms            aaa;
 210
 211     if (bIncludeSolute)
 212     {
 213         i0 = atoms->nr;
 214         r0 = atoms->nres;
 215         srenew(atoms->atom, atoms->nr+ndots);
 216         memset(&atoms->atom[i0], 0, sizeof(*atoms->atom)*ndots);
 217         srenew(atoms->atomname, atoms->nr+ndots);
 218         srenew(atoms->resinfo, r0+1);
 219         atoms->atom[i0].resind = r0;
 220         t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0+1, ' ', 0, ' ');
 221         if (atoms->pdbinfo != nullptr)
 222         {
 223             srenew(atoms->pdbinfo, atoms->nr+ndots);
 224         }
 225         snew(xnew, atoms->nr+ndots);
 226         for (i = 0; (i < atoms->nr); i++)
 227         {
 228             copy_rvec(x[i], xnew[i]);
 229         }
 230         for (i = k = 0; (i < ndots); i++)
 231         {
 232             ii0                        = i0+i;
 233             atoms->atomname[ii0]       = put_symtab(symtab, atomnm);
 234             atoms->atom[ii0].resind    = r0;
 235             xnew[ii0][XX]              = dots[k++];
 236             xnew[ii0][YY]              = dots[k++];
 237             xnew[ii0][ZZ]              = dots[k++];
 238             if (atoms->pdbinfo != nullptr)
 239             {
 240                 atoms->pdbinfo[ii0].type   = epdbATOM;
 241                 atoms->pdbinfo[ii0].atomnr = ii0;
 242                 atoms->pdbinfo[ii0].bfac   = 0.0;
 243                 atoms->pdbinfo[ii0].occup  = 0.0;
 244             }
 245         }
 246         atoms->nr   = i0+ndots;
 247         atoms->nres = r0+1;
 248         write_sto_conf(fn, title, atoms, xnew, nullptr, ePBC, const_cast<rvec *>(box));
 249         atoms->nres = r0;
 250         atoms->nr   = i0;
 251     }
 252     else
 253     {
 254         init_t_atoms(&aaa, ndots, TRUE);
 255         aaa.atom[0].resind = 0;
 256         t_atoms_set_resinfo(&aaa, 0, symtab, resnm, 1, ' ', 0, ' ');
 257         snew(xnew, ndots);
 258         for (i = k = 0; (i < ndots); i++)
 259         {
 260             ii0                     = i;
 261             aaa.atomname[ii0]       = put_symtab(symtab, atomnm);
 262             aaa.pdbinfo[ii0].type   = epdbATOM;
 263             aaa.pdbinfo[ii0].atomnr = ii0;
 264             aaa.atom[ii0].resind    = 0;
 265             xnew[ii0][XX]           = dots[k++];
 266             xnew[ii0][YY]           = dots[k++];
 267             xnew[ii0][ZZ]           = dots[k++];
 268             aaa.pdbinfo[ii0].bfac   = 0.0;
 269             aaa.pdbinfo[ii0].occup  = 0.0;
 270         }
 271         aaa.nr = ndots;
 272         write_sto_conf(fn, title, &aaa, xnew, nullptr, ePBC, const_cast<rvec *>(box));
 273         do_conect(fn, ndots, xnew);
 274         done_atom(&aaa);
 275     }
 276     sfree(xnew);
 277 }
 278
 279 /********************************************************************
 280  * Actual analysis module
 281  */
 282
 283 /*! \brief
 284  * Implements `gmx sas` trajectory analysis module.
 285  */
 286 class Sasa : public TrajectoryAnalysisModule
 287 {
 288     public:
 289         Sasa();
 290
 291         void initOptions(IOptionsContainer          *options,
 292                          TrajectoryAnalysisSettings *settings) override;
 293         void initAnalysis(const TrajectoryAnalysisSettings &settings,
 294                           const TopologyInformation        &top) override;
 295
 296         TrajectoryAnalysisModuleDataPointer startFrames(
 297             const AnalysisDataParallelOptions &opt,
 298             const SelectionCollection         &selections) override;
 299         void analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 300                           TrajectoryAnalysisModuleData *pdata) override;
 301
 302         void finishAnalysis(int nframes) override;
 303         void writeOutput() override;
 304
 305     private:
 306         /*! \brief
 307          * Surface areas as a function of time.
 308          *
 309          * First column is for the calculation group, and the rest for the
 310          * output groups.  This data is always produced.
 311          */
 312         AnalysisData            area_;
 313         /*! \brief
 314          * Per-atom surface areas as a function of time.
 315          *
 316          * Contains one data set for each column in `area_`.
 317          * Each column corresponds to a selection position in `surfaceSel_`.
 318          * This data is only produced if atom or residue areas have been
 319          * requested.
 320          */
 321         AnalysisData            atomArea_;
 322         /*! \brief
 323          * Per-residue surface areas as a function of time.
 324          *
 325          * Contains one data set for each column in `area_`.
 326          * Each column corresponds to a distinct residue `surfaceSel_`.
 327          * For example, if `surfaceSel_` selects residues 2, 5, and 7, there
 328          * will be three columns here.
 329          * This data is only produced if atom or residue areas have been
 330          * requested.
 331          */
 332         AnalysisData            residueArea_;
 333         /*! \brief
 334          * Free energy estimates as a function of time.
 335          *
 336          * Column layout is the same as for `area_`.
 337          * This data is only produced if the output is requested.
 338          */
 339         AnalysisData            dgSolv_;
 340         /*! \brief
 341          * Total volume and density of the calculation group as a function of
 342          * time.
 343          *
 344          * The first column is the volume and the second column is the density.
 345          * This data is only produced if the output is requested.
 346          */
 347         AnalysisData            volume_;
 348
 349         /*! \brief
 350          * The selection to calculate the surface for.
 351          *
 352          * Selection::originalId() and Selection::mappedId() store the mapping
 353          * from the positions to the columns of `residueArea_`.
 354          * The selection is computed with SelectionOption::dynamicMask(), i.e.,
 355          * even in the presence of a dynamic selection, the number of returned
 356          * positions is fixed, and SelectionPosition::selected() is used.
 357          */
 358         Selection               surfaceSel_;
 359         /*! \brief
 360          * List of optional additional output groups.
 361          *
 362          * Each of these must be a subset of the `surfaceSel_`.
 363          * Selection::originalId() and Selection::mappedId() store the mapping
 364          * from the positions to the corresponsing positions in `surfaceSel_`.
 365          */
 366         SelectionList           outputSel_;
 367
 368         std::string             fnArea_;
 369         std::string             fnAtomArea_;
 370         std::string             fnResidueArea_;
 371         std::string             fnDGSolv_;
 372         std::string             fnVolume_;
 373         std::string             fnConnolly_;
 374
 375         double                  solsize_;
 376         int                     ndots_;
 377         //double                  minarea_;
 378         double                  dgsDefault_;
 379         bool                    bIncludeSolute_;
 380
 381         t_topology             *top_;
 382         //! Combined VdW and probe radii for each atom in the calculation group.
 383         std::vector<real>       radii_;
 384         /*! \brief
 385          * Solvation free energy coefficients for each atom in the calculation
 386          * group.
 387          *
 388          * Empty if the free energy output has not been requested.
 389          */
 390         std::vector<real>       dgsFactor_;
 391         //! Calculation algorithm.
 392         SurfaceAreaCalculator   calculator_;
 393
 394         // Copy and assign disallowed by base.
 395 };
 396
 397 Sasa::Sasa()
 398     : solsize_(0.14), ndots_(24), dgsDefault_(0), bIncludeSolute_(true), top_(nullptr)
 399 {
 400     //minarea_ = 0.5;
 401     registerAnalysisDataset(&area_, "area");
 402     registerAnalysisDataset(&atomArea_, "atomarea");
 403     registerAnalysisDataset(&residueArea_, "resarea");
 404     registerAnalysisDataset(&dgSolv_, "dgsolv");
 405     registerAnalysisDataset(&volume_, "volume");
 406 }
 407
 408 void
 409 Sasa::initOptions(IOptionsContainer *options, TrajectoryAnalysisSettings *settings)
 410 {
 411     static const char *const desc[] = {
 412         "[THISMODULE] computes solvent accessible surface areas.",
 413         "See Eisenhaber F, Lijnzaad P, Argos P, Sander C, & Scharf M",
 414         "(1995) J. Comput. Chem. 16, 273-284 for the algorithm used.",
 415         "With [TT]-q[tt], the Connolly surface can be generated as well",
 416         "in a [REF].pdb[ref] file where the nodes are represented as atoms",
 417         "and the edges connecting the nearest nodes as CONECT records.",
 418         "[TT]-odg[tt] allows for estimation of solvation free energies",
 419         "from per-atom solvation energies per exposed surface area.[PAR]",
 420
 421         "The program requires a selection for the surface calculation to be",
 422         "specified with [TT]-surface[tt]. This should always consist of all",
 423         "non-solvent atoms in the system. The area of this group is always",
 424         "calculated. Optionally, [TT]-output[tt] can specify additional",
 425         "selections, which should be subsets of the calculation group.",
 426         "The solvent-accessible areas for these groups are also extracted",
 427         "from the full surface.[PAR]",
 428
 429         "The average and standard deviation of the area over the trajectory",
 430         "can be calculated per residue and atom (options [TT]-or[tt] and",
 431         "[TT]-oa[tt]).[PAR]",
 432         //"In combination with the latter option an [REF].itp[ref] file can be",
 433         //"generated (option [TT]-i[tt])",
 434         //"which can be used to restrain surface atoms.[PAR]",
 435
 436         "With the [TT]-tv[tt] option the total volume and density of the",
 437         "molecule can be computed. With [TT]-pbc[tt] (the default), you",
 438         "must ensure that your molecule/surface group is not split across PBC.",
 439         "Otherwise, you will get non-sensical results.",
 440         "Please also consider whether the normal probe radius is appropriate",
 441         "in this case or whether you would rather use, e.g., 0. It is good",
 442         "to keep in mind that the results for volume and density are very",
 443         "approximate. For example, in ice Ih, one can easily fit water molecules in the",
 444         "pores which would yield a volume that is too low, and surface area and density",
 445         "that are both too high."
 446     };
 447
 448     settings->setHelpText(desc);
 449
 450     options->addOption(FileNameOption("o").filetype(eftPlot).outputFile().required()
 451                            .store(&fnArea_).defaultBasename("area")
 452                            .description("Total area as a function of time"));
 453     options->addOption(FileNameOption("odg").filetype(eftPlot).outputFile()
 454                            .store(&fnDGSolv_).defaultBasename("dgsolv")
 455                            .description("Estimated solvation free energy as a function of time"));
 456     options->addOption(FileNameOption("or").filetype(eftPlot).outputFile()
 457                            .store(&fnResidueArea_).defaultBasename("resarea")
 458                            .description("Average area per residue"));
 459     options->addOption(FileNameOption("oa").filetype(eftPlot).outputFile()
 460                            .store(&fnAtomArea_).defaultBasename("atomarea")
 461                            .description("Average area per atom"));
 462     options->addOption(FileNameOption("tv").filetype(eftPlot).outputFile()
 463                            .store(&fnVolume_).defaultBasename("volume")
 464                            .description("Total volume and density as a function of time"));
 465     options->addOption(FileNameOption("q").filetype(eftPDB).outputFile()
 466                            .store(&fnConnolly_).defaultBasename("connolly")
 467                            .description("PDB file for Connolly surface"));
 468     //options->addOption(FileNameOption("i").filetype(eftITP).outputFile()
 469     //                       .store(&fnRestraints_).defaultBasename("surfat")
 470     //                       .description("Topology file for position restraings on surface atoms"));
 471
 472
 473     options->addOption(DoubleOption("probe").store(&solsize_)
 474                            .description("Radius of the solvent probe (nm)"));
 475     options->addOption(IntegerOption("ndots").store(&ndots_)
 476                            .description("Number of dots per sphere, more dots means more accuracy"));
 477     //options->addOption(DoubleOption("minarea").store(&minarea_)
 478     //                       .description("The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file (see help)"));
 479     options->addOption(BooleanOption("prot").store(&bIncludeSolute_)
 480                            .description("Output the protein to the Connolly [REF].pdb[ref] file too"));
 481     options->addOption(DoubleOption("dgs").store(&dgsDefault_)
 482                            .description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
 483
 484     // Selections must select atoms for the VdW radii lookup to work.
 485     // The calculation group uses dynamicMask() so that the coordinates
 486     // match a static array of VdW radii.
 487     options->addOption(SelectionOption("surface").store(&surfaceSel_)
 488                            .required().onlySortedAtoms().dynamicMask()
 489                            .description("Surface calculation selection"));
 490     options->addOption(SelectionOption("output").storeVector(&outputSel_)
 491                            .onlySortedAtoms().multiValue()
 492                            .description("Output selection(s)"));
 493
 494     // Atom names etc. are required for the VdW radii lookup.
 495     settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 496 }
 497
 498 void
 499 Sasa::initAnalysis(const TrajectoryAnalysisSettings &settings,
 500                    const TopologyInformation        &top)
 501 {
 502     const t_atoms &atoms = top.topology()->atoms;
 503     top_ = top.topology();
 504
 505     //bITP   = opt2bSet("-i", nfile, fnm);
 506     const bool bResAt =
 507         !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
 508     const bool bDGsol = !fnDGSolv_.empty();
 509
 510     if (solsize_ < 0)
 511     {
 512         solsize_ = 1e-3;
 513         fprintf(stderr, "Probe size too small, setting it to %g\n", solsize_);
 514     }
 515     if (ndots_ < 20)
 516     {
 517         ndots_ = 20;
 518         fprintf(stderr, "Ndots too small, setting it to %d\n", ndots_);
 519     }
 520
 521     please_cite(stderr, "Eisenhaber95");
 522     //if ((top.ePBC() != epbcXYZ) || (TRICLINIC(fr.box)))
 523     //{
 524     //    fprintf(stderr, "\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
 525     //            "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
 526     //            "will certainly crash the analysis.\n\n");
 527     //}
 528
 529     if (bDGsol)
 530     {
 531         if (!top.hasFullTopology())
 532         {
 533             GMX_THROW(InconsistentInputError("Cannot compute Delta G of solvation without a tpr file"));
 534         }
 535         else
 536         {
 537             if (strcmp(*(atoms.atomtype[0]), "?") == 0)
 538             {
 539                 GMX_THROW(InconsistentInputError("Your input tpr file is too old (does not contain atom types). Cannot not compute Delta G of solvation"));
 540             }
 541             else
 542             {
 543                 printf("Free energy of solvation predictions:\n");
 544                 please_cite(stdout, "Eisenberg86a");
 545             }
 546         }
 547     }
 548
 549     // Now compute atomic radii including solvent probe size.
 550     // Also, fetch solvation free energy coefficients and
 551     // compute the residue indices that map the calculation atoms
 552     // to the columns of residueArea_.
 553     radii_.reserve(surfaceSel_.posCount());
 554     if (bDGsol)
 555     {
 556         dgsFactor_.reserve(surfaceSel_.posCount());
 557     }
 558
 559     const int resCount = surfaceSel_.initOriginalIdsToGroup(top.mtop(), INDEX_RES);
 560
 561     // TODO: Not exception-safe, but nice solution would be to have a C++
 562     // atom properties class...
 563     gmx_atomprop_t      aps = gmx_atomprop_init();
 564
 565     ArrayRef<const int> atomIndices = surfaceSel_.atomIndices();
 566     int                 ndefault    = 0;
 567     for (int i = 0; i < surfaceSel_.posCount(); i++)
 568     {
 569         const int ii     = atomIndices[i];
 570         const int resind = atoms.atom[ii].resind;
 571         real      radius;
 572         if (!gmx_atomprop_query(aps, epropVDW,
 573                                 *(atoms.resinfo[resind].name),
 574                                 *(atoms.atomname[ii]), &radius))
 575         {
 576             ndefault++;
 577         }
 578         radii_.push_back(radius + solsize_);
 579         if (bDGsol)
 580         {
 581             real dgsFactor;
 582             if (!gmx_atomprop_query(aps, epropDGsol,
 583                                     *(atoms.resinfo[resind].name),
 584                                     *(atoms.atomtype[ii]), &dgsFactor))
 585             {
 586                 dgsFactor = dgsDefault_;
 587             }
 588             dgsFactor_.push_back(dgsFactor);
 589         }
 590     }
 591     if (ndefault > 0)
 592     {
 593         fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
 594     }
 595     gmx_atomprop_destroy(aps);
 596
 597     // Pre-compute mapping from the output groups to the calculation group,
 598     // and store it in the selection ID map for easy lookup.
 599     for (size_t g = 0; g < outputSel_.size(); ++g)
 600     {
 601         ArrayRef<const int> outputIndices = outputSel_[g].atomIndices();
 602         for (int i = 0, j = 0; i < outputSel_[g].posCount(); ++i)
 603         {
 604             while (j < surfaceSel_.posCount() && outputIndices[i] > atomIndices[j])
 605             {
 606                 ++j;
 607             }
 608             if (j == surfaceSel_.posCount() || outputIndices[i] != atomIndices[j])
 609             {
 610                 const std::string message
 611                     = formatString("Output selection '%s' is not a subset of "
 612                                    "the surface selection (atom %d is the first "
 613                                    "atom not in the surface selection)",
 614                                    outputSel_[g].name(), outputIndices[i] + 1);
 615                 GMX_THROW(InconsistentInputError(message));
 616             }
 617             outputSel_[g].setOriginalId(i, j);
 618         }
 619     }
 620
 621     calculator_.setDotCount(ndots_);
 622     calculator_.setRadii(radii_);
 623
 624     // Initialize all the output data objects and initialize the output plotters.
 625
 626     area_.setColumnCount(0, 1 + outputSel_.size());
 627     {
 628         AnalysisDataPlotModulePointer plotm(
 629                 new AnalysisDataPlotModule(settings.plotSettings()));
 630         plotm->setFileName(fnArea_);
 631         plotm->setTitle("Solvent Accessible Surface");
 632         plotm->setXAxisIsTime();
 633         plotm->setYLabel("Area (nm\\S2\\N)");
 634         plotm->appendLegend("Total");
 635         for (size_t i = 0; i < outputSel_.size(); ++i)
 636         {
 637             plotm->appendLegend(outputSel_[i].name());
 638         }
 639         area_.addModule(plotm);
 640     }
 641
 642     if (bResAt)
 643     {
 644         atomArea_.setDataSetCount(1 + outputSel_.size());
 645         residueArea_.setDataSetCount(1 + outputSel_.size());
 646         for (size_t i = 0; i <= outputSel_.size(); ++i)
 647         {
 648             atomArea_.setColumnCount(i, surfaceSel_.posCount());
 649             residueArea_.setColumnCount(i, resCount);
 650         }
 651         {
 652             AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
 653             for (int i = 0; i < surfaceSel_.posCount(); ++i)
 654             {
 655                 avem->setXAxisValue(i, surfaceSel_.position(i).atomIndices()[0] + 1);
 656             }
 657             atomArea_.addModule(avem);
 658             if (!fnAtomArea_.empty())
 659             {
 660                 AnalysisDataPlotModulePointer plotm(
 661                         new AnalysisDataPlotModule(settings.plotSettings()));
 662                 plotm->setFileName(fnAtomArea_);
 663                 plotm->setTitle("Area per atom over the trajectory");
 664                 plotm->setXLabel("Atom");
 665                 plotm->setXFormat(8, 0);
 666                 plotm->setYLabel("Area (nm\\S2\\N)");
 667                 plotm->setErrorsAsSeparateColumn(true);
 668                 plotm->appendLegend("Average (nm\\S2\\N)");
 669                 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
 670                 avem->addModule(plotm);
 671             }
 672         }
 673         {
 674             AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
 675             int nextRow = 0;
 676             for (int i = 0; i < surfaceSel_.posCount(); ++i)
 677             {
 678                 const int residueGroup = surfaceSel_.position(i).mappedId();
 679                 if (residueGroup >= nextRow)
 680                 {
 681                     GMX_ASSERT(residueGroup == nextRow,
 682                                "Inconsistent (non-uniformly increasing) residue grouping");
 683                     const int atomIndex    = surfaceSel_.position(i).atomIndices()[0];
 684                     const int residueIndex = atoms.atom[atomIndex].resind;
 685                     avem->setXAxisValue(nextRow, atoms.resinfo[residueIndex].nr);
 686                     ++nextRow;
 687                 }
 688             }
 689             residueArea_.addModule(avem);
 690             if (!fnResidueArea_.empty())
 691             {
 692                 AnalysisDataPlotModulePointer plotm(
 693                         new AnalysisDataPlotModule(settings.plotSettings()));
 694                 plotm->setFileName(fnResidueArea_);
 695                 plotm->setTitle("Area per residue over the trajectory");
 696                 plotm->setXLabel("Residue");
 697                 plotm->setXFormat(8, 0);
 698                 plotm->setYLabel("Area (nm\\S2\\N)");
 699                 plotm->setErrorsAsSeparateColumn(true);
 700                 plotm->appendLegend("Average (nm\\S2\\N)");
 701                 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
 702                 avem->addModule(plotm);
 703             }
 704         }
 705     }
 706
 707     if (!fnDGSolv_.empty())
 708     {
 709         dgSolv_.setColumnCount(0, 1 + outputSel_.size());
 710         AnalysisDataPlotModulePointer plotm(
 711                 new AnalysisDataPlotModule(settings.plotSettings()));
 712         plotm->setFileName(fnDGSolv_);
 713         plotm->setTitle("Free Energy of Solvation");
 714         plotm->setXAxisIsTime();
 715         plotm->setYLabel("D Gsolv");
 716         plotm->appendLegend("Total");
 717         for (size_t i = 0; i < outputSel_.size(); ++i)
 718         {
 719             plotm->appendLegend(outputSel_[i].name());
 720         }
 721         dgSolv_.addModule(plotm);
 722     }
 723
 724     if (!fnVolume_.empty())
 725     {
 726         volume_.setColumnCount(0, 2);
 727         AnalysisDataPlotModulePointer plotm(
 728                 new AnalysisDataPlotModule(settings.plotSettings()));
 729         plotm->setFileName(fnVolume_);
 730         plotm->setTitle("Volume and Density");
 731         plotm->setXAxisIsTime();
 732         plotm->appendLegend("Volume (nm\\S3\\N)");
 733         plotm->appendLegend("Density (g/l)");
 734         volume_.addModule(plotm);
 735     }
 736 }
 737
 738 /*! \brief
 739  * Temporary memory for use within a single-frame calculation.
 740  */
 741 class SasaModuleData : public TrajectoryAnalysisModuleData
 742 {
 743     public:
 744         /*! \brief
 745          * Reserves memory for the frame-local data.
 746          *
 747          * `residueCount` will be zero if per-residue data is not being
 748          * calculated.
 749          */
 750         SasaModuleData(TrajectoryAnalysisModule          *module,
 751                        const AnalysisDataParallelOptions &opt,
 752                        const SelectionCollection         &selections,
 753                        int atomCount, int residueCount)
 754             : TrajectoryAnalysisModuleData(module, opt, selections)
 755         {
 756             index_.reserve(atomCount);
 757             // If the calculation group is not dynamic, pre-calculate
 758             // the index, since it is not going to change.
 759             for (int i = 0; i < atomCount; ++i)
 760             {
 761                 index_.push_back(i);
 762             }
 763             atomAreas_.resize(atomCount);
 764             res_a_.resize(residueCount);
 765         }
 766
 767         void finish() override { finishDataHandles(); }
 768
 769         //! Indices of the calculation selection positions selected for the frame.
 770         std::vector<int>        index_;
 771         /*! \brief
 772          * Atom areas for each calculation selection position for the frame.
 773          *
 774          * One entry for each position in the calculation group.
 775          * Values for atoms not selected are set to zero.
 776          */
 777         std::vector<real>       atomAreas_;
 778         /*! \brief
 779          * Working array to accumulate areas for each residue.
 780          *
 781          * One entry for each distinct residue in the calculation group;
 782          * indices are not directly residue numbers or residue indices.
 783          *
 784          * This vector is empty if residue area calculations are not being
 785          * performed.
 786          */
 787         std::vector<real>       res_a_;
 788 };
 789
 790 TrajectoryAnalysisModuleDataPointer Sasa::startFrames(
 791         const AnalysisDataParallelOptions &opt,
 792         const SelectionCollection         &selections)
 793 {
 794     return TrajectoryAnalysisModuleDataPointer(
 795             new SasaModuleData(this, opt, selections, surfaceSel_.posCount(),
 796                                residueArea_.columnCount(0)));
 797 }
 798
 799 /*! \brief
 800  * Helper method to compute the areas for a single selection.
 801  *
 802  * \param[in]  surfaceSel     The calculation selection.
 803  * \param[in]  sel            The selection to compute the areas for (can be
 804  *     `surfaceSel` or one of the output selections).
 805  * \param[in]  atomAreas      Atom areas for each position in `surfaceSel`.
 806  * \param[in]  dgsFactor      Free energy coefficients for each position in
 807  *     `surfaceSel`. If empty, free energies are not calculated.
 808  * \param[out] totalAreaOut Total area of `sel` (sum of atom areas it selects).
 809  * \param[out] dgsolvOut      Solvation free energy.
 810  *     Will be zero of `dgsFactor` is empty.
 811  * \param      atomAreaHandle Data handle to use for storing atom areas for `sel`.
 812  * \param      resAreaHandle  Data handle to use for storing residue areas for `sel`.
 813  * \param      resAreaWork    Work array for accumulating the residue areas.
 814  *     If empty, atom and residue areas are not calculated.
 815  *
 816  * `atomAreaHandle` and `resAreaHandle` are not used if `resAreaWork` is empty.
 817  */
 818 void computeAreas(const Selection &surfaceSel, const Selection &sel,
 819                   const std::vector<real> &atomAreas,
 820                   const std::vector<real> &dgsFactor,
 821                   real *totalAreaOut, real *dgsolvOut,
 822                   AnalysisDataHandle atomAreaHandle,
 823                   AnalysisDataHandle resAreaHandle,
 824                   std::vector<real> *resAreaWork)
 825 {
 826     const bool bResAt    = !resAreaWork->empty();
 827     const bool bDGsolv   = !dgsFactor.empty();
 828     real       totalArea = 0;
 829     real       dgsolv    = 0;
 830
 831     if (bResAt)
 832     {
 833         std::fill(resAreaWork->begin(), resAreaWork->end(), 0.0_real);
 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                       0.0_real);
 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