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, 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/fileio/trx.h"
  58 #include "gromacs/legacyheaders/copyrite.h"
  59 #include "gromacs/math/units.h"
  60 #include "gromacs/math/vec.h"
  61 #include "gromacs/options/basicoptions.h"
  62 #include "gromacs/options/filenameoption.h"
  63 #include "gromacs/options/options.h"
  64 #include "gromacs/pbcutil/pbc.h"
  65 #include "gromacs/selection/selection.h"
  66 #include "gromacs/selection/selectionoption.h"
  67 #include "gromacs/topology/atomprop.h"
  68 #include "gromacs/topology/symtab.h"
  69 #include "gromacs/topology/topology.h"
  70 #include "gromacs/trajectoryanalysis/analysismodule.h"
  71 #include "gromacs/trajectoryanalysis/analysissettings.h"
  72 #include "gromacs/utility/exceptions.h"
  73 #include "gromacs/utility/futil.h"
  74 #include "gromacs/utility/scoped_ptr_sfree.h"
  75 #include "gromacs/utility/smalloc.h"
  76 #include "gromacs/utility/stringutil.h"
  77
  78 #include "nsc.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     atom_id  aa;
  97     //! Index of the nearest neighbor dot.
  98     atom_id  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[], atom_id i, atom_id j, real d2)
 114 {
 115     if (c[i].aa == NO_ATID)
 116     {
 117         c[i].aa  = j;
 118         c[i].d2a = d2;
 119     }
 120     else if (c[i].ab == NO_ATID)
 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 = NO_ATID;
 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 == NO_ATID) || (c[i].ab == NO_ATID))
 187         {
 188             fprintf(stderr, "Warning dot %d has no conections\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 != NULL)
 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 != NULL)
 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, NULL, 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, NULL, ePBC, const_cast<rvec *>(box));
 272         do_conect(fn, ndots, xnew);
 273         free_t_atoms(&aaa, FALSE);
 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(Options                    *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
 391         // Copy and assign disallowed by base.
 392 };
 393
 394 Sasa::Sasa()
 395     : TrajectoryAnalysisModule(SasaInfo::name, SasaInfo::shortDescription),
 396       solsize_(0.14), ndots_(24), dgsDefault_(0), bIncludeSolute_(true), top_(NULL)
 397 {
 398     //minarea_ = 0.5;
 399     registerAnalysisDataset(&area_, "area");
 400     registerAnalysisDataset(&atomArea_, "atomarea");
 401     registerAnalysisDataset(&residueArea_, "resarea");
 402     registerAnalysisDataset(&dgSolv_, "dgsolv");
 403     registerAnalysisDataset(&volume_, "volume");
 404 }
 405
 406 void
 407 Sasa::initOptions(Options *options, TrajectoryAnalysisSettings *settings)
 408 {
 409     static const char *const desc[] = {
 410         "[THISMODULE] computes solvent accessible surface areas.",
 411         "See Eisenhaber F, Lijnzaad P, Argos P, Sander C, & Scharf M",
 412         "(1995) J. Comput. Chem. 16, 273-284 for the algorithm used.",
 413         "With [TT]-q[tt], the Connolly surface can be generated as well",
 414         "in a [TT].pdb[tt] file where the nodes are represented as atoms",
 415         "and the edges connecting the nearest nodes as CONECT records.",
 416         "[TT]-odg[tt] allows for estimation of solvation free energies",
 417         "from per-atom solvation energies per exposed surface area.[PAR]",
 418
 419         "The program requires a selection for the surface calculation to be",
 420         "specified with [TT]-surface[tt]. This should always consist of all",
 421         "non-solvent atoms in the system. The area of this group is always",
 422         "calculated. Optionally, [TT]-output[tt] can specify additional",
 423         "selections, which should be subsets of the calculation group.",
 424         "The solvent-accessible areas for these groups are also extracted",
 425         "from the full surface.[PAR]",
 426
 427         "The average and standard deviation of the area over the trajectory",
 428         "can be calculated per residue and atom (options [TT]-or[tt] and",
 429         "[TT]-oa[tt]).[PAR]",
 430         //"In combination with the latter option an [TT].itp[tt] file can be",
 431         //"generated (option [TT]-i[tt])",
 432         //"which can be used to restrain surface atoms.[PAR]",
 433
 434         "With the [TT]-tv[tt] option the total volume and density of the",
 435         "molecule can be computed.",
 436         "Please consider whether the normal probe radius is appropriate",
 437         "in this case or whether you would rather use, e.g., 0. It is good",
 438         "to keep in mind that the results for volume and density are very",
 439         "approximate. For example, in ice Ih, one can easily fit water molecules in the",
 440         "pores which would yield a volume that is too low, and surface area and density",
 441         "that are both too high."
 442     };
 443
 444     options->setDescription(concatenateStrings(desc));
 445
 446     options->addOption(FileNameOption("o").filetype(eftPlot).outputFile().required()
 447                            .store(&fnArea_).defaultBasename("area")
 448                            .description("Total area as a function of time"));
 449     options->addOption(FileNameOption("odg").filetype(eftPlot).outputFile()
 450                            .store(&fnDGSolv_).defaultBasename("dgsolv")
 451                            .description("Estimated solvation free energy as a function of time"));
 452     options->addOption(FileNameOption("or").filetype(eftPlot).outputFile()
 453                            .store(&fnResidueArea_).defaultBasename("resarea")
 454                            .description("Average area per residue"));
 455     options->addOption(FileNameOption("oa").filetype(eftPlot).outputFile()
 456                            .store(&fnAtomArea_).defaultBasename("atomarea")
 457                            .description("Average area per atom"));
 458     options->addOption(FileNameOption("tv").filetype(eftPlot).outputFile()
 459                            .store(&fnVolume_).defaultBasename("volume")
 460                            .description("Total volume and density as a function of time"));
 461     options->addOption(FileNameOption("q").filetype(eftPDB).outputFile()
 462                            .store(&fnConnolly_).defaultBasename("connolly")
 463                            .description("PDB file for Connolly surface"));
 464     //options->addOption(FileNameOption("i").filetype(eftITP).outputFile()
 465     //                       .store(&fnRestraints_).defaultBasename("surfat")
 466     //                       .description("Topology file for position restraings on surface atoms"));
 467
 468
 469     options->addOption(DoubleOption("probe").store(&solsize_)
 470                            .description("Radius of the solvent probe (nm)"));
 471     options->addOption(IntegerOption("ndots").store(&ndots_)
 472                            .description("Number of dots per sphere, more dots means more accuracy"));
 473     //options->addOption(DoubleOption("minarea").store(&minarea_)
 474     //                       .description("The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file (see help)"));
 475     options->addOption(BooleanOption("prot").store(&bIncludeSolute_)
 476                            .description("Output the protein to the Connolly [TT].pdb[tt] file too"));
 477     options->addOption(DoubleOption("dgs").store(&dgsDefault_)
 478                            .description("Default value for solvation free energy per area (kJ/mol/nm^2)"));
 479
 480     // Selections must select atoms for the VdW radii lookup to work.
 481     // The calculation group uses dynamicMask() so that the coordinates
 482     // match a static array of VdW radii.
 483     options->addOption(SelectionOption("surface").store(&surfaceSel_)
 484                            .required().onlyAtoms().dynamicMask()
 485                            .description("Surface calculation selection"));
 486     options->addOption(SelectionOption("output").storeVector(&outputSel_)
 487                            .onlyAtoms().multiValue()
 488                            .description("Output selection(s)"));
 489
 490     // Atom names etc. are required for the VdW radii lookup.
 491     settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);
 492 }
 493
 494 void
 495 Sasa::initAnalysis(const TrajectoryAnalysisSettings &settings,
 496                    const TopologyInformation        &top)
 497 {
 498     const t_atoms &atoms = top.topology()->atoms;
 499     top_ = top.topology();
 500
 501     //bITP   = opt2bSet("-i", nfile, fnm);
 502     const bool bResAt =
 503         !fnResidueArea_.empty() || !fnAtomArea_.empty(); // || bITP;
 504     const bool bDGsol = !fnDGSolv_.empty();
 505
 506     if (solsize_ < 0)
 507     {
 508         solsize_ = 1e-3;
 509         fprintf(stderr, "Probe size too small, setting it to %g\n", solsize_);
 510     }
 511     if (ndots_ < 20)
 512     {
 513         ndots_ = 20;
 514         fprintf(stderr, "Ndots too small, setting it to %d\n", ndots_);
 515     }
 516
 517     please_cite(stderr, "Eisenhaber95");
 518     //if ((top.ePBC() != epbcXYZ) || (TRICLINIC(fr.box)))
 519     //{
 520     //    fprintf(stderr, "\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
 521     //            "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
 522     //            "will certainly crash the analysis.\n\n");
 523     //}
 524
 525     if (bDGsol)
 526     {
 527         if (!top.hasFullTopology())
 528         {
 529             GMX_THROW(InconsistentInputError("Cannot compute Delta G of solvation without a tpr file"));
 530         }
 531         else
 532         {
 533             if (strcmp(*(atoms.atomtype[0]), "?") == 0)
 534             {
 535                 GMX_THROW(InconsistentInputError("Your input tpr file is too old (does not contain atom types). Cannot not compute Delta G of solvation"));
 536             }
 537             else
 538             {
 539                 printf("Free energy of solvation predictions:\n");
 540                 please_cite(stdout, "Eisenberg86a");
 541             }
 542         }
 543     }
 544
 545     // Now compute atomic radii including solvent probe size.
 546     // Also, fetch solvation free energy coefficients and
 547     // compute the residue indices that map the calculation atoms
 548     // to the columns of residueArea_.
 549     radii_.reserve(surfaceSel_.posCount());
 550     if (bDGsol)
 551     {
 552         dgsFactor_.reserve(surfaceSel_.posCount());
 553     }
 554
 555     // TODO: Not exception-safe, but nice solution would be to have a C++
 556     // atom properties class...
 557     gmx_atomprop_t     aps = gmx_atomprop_init();
 558
 559     ConstArrayRef<int> atomIndices = surfaceSel_.atomIndices();
 560     int                prevResind  = atoms.atom[atomIndices[0]].resind;
 561     int                resCount    = 0;
 562     int                ndefault    = 0;
 563     for (int i = 0; i < surfaceSel_.posCount(); i++)
 564     {
 565         const int ii     = atomIndices[i];
 566         const int resind = atoms.atom[ii].resind;
 567         if (resind != prevResind)
 568         {
 569             ++resCount;
 570             prevResind = resind;
 571         }
 572         surfaceSel_.setOriginalId(i, resCount);
 573         real      radius;
 574         if (!gmx_atomprop_query(aps, epropVDW,
 575                                 *(atoms.resinfo[resind].name),
 576                                 *(atoms.atomname[ii]), &radius))
 577         {
 578             ndefault++;
 579         }
 580         radii_.push_back(radius + solsize_);
 581         if (bDGsol)
 582         {
 583             real dgsFactor;
 584             if (!gmx_atomprop_query(aps, epropDGsol,
 585                                     *(atoms.resinfo[resind].name),
 586                                     *(atoms.atomtype[ii]), &dgsFactor))
 587             {
 588                 dgsFactor = dgsDefault_;
 589             }
 590             dgsFactor_.push_back(dgsFactor);
 591         }
 592     }
 593     if (ndefault > 0)
 594     {
 595         fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
 596     }
 597     gmx_atomprop_destroy(aps);
 598
 599     // The loop above doesn't count the last residue.
 600     ++resCount;
 601
 602     // Pre-compute mapping from the output groups to the calculation group,
 603     // and store it in the selection ID map for easy lookup.
 604     for (size_t g = 0; g < outputSel_.size(); ++g)
 605     {
 606         ConstArrayRef<int> outputIndices = outputSel_[g].atomIndices();
 607         for (int i = 0, j = 0; i < outputSel_[g].posCount(); ++i)
 608         {
 609             while (j < surfaceSel_.posCount() && outputIndices[i] > atomIndices[j])
 610             {
 611                 ++j;
 612             }
 613             if (j == surfaceSel_.posCount() || outputIndices[i] != atomIndices[j])
 614             {
 615                 GMX_THROW(InconsistentInputError("Output selection is not a subset of the input selection"));
 616             }
 617             outputSel_[g].setOriginalId(i, j);
 618         }
 619     }
 620
 621     // Initialize all the output data objects and initialize the output plotters.
 622
 623     area_.setColumnCount(0, 1 + outputSel_.size());
 624     {
 625         AnalysisDataPlotModulePointer plotm(
 626                 new AnalysisDataPlotModule(settings.plotSettings()));
 627         plotm->setFileName(fnArea_);
 628         plotm->setTitle("Solvent Accessible Surface");
 629         plotm->setXAxisIsTime();
 630         plotm->setYLabel("Area (nm\\S2\\N)");
 631         plotm->appendLegend("Total");
 632         for (size_t i = 0; i < outputSel_.size(); ++i)
 633         {
 634             plotm->appendLegend(outputSel_[i].name());
 635         }
 636         area_.addModule(plotm);
 637     }
 638
 639     if (bResAt)
 640     {
 641         atomArea_.setDataSetCount(1 + outputSel_.size());
 642         residueArea_.setDataSetCount(1 + outputSel_.size());
 643         for (size_t i = 0; i <= outputSel_.size(); ++i)
 644         {
 645             atomArea_.setColumnCount(i, surfaceSel_.posCount());
 646             residueArea_.setColumnCount(i, resCount);
 647         }
 648         {
 649             AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
 650             for (int i = 0; i < surfaceSel_.posCount(); ++i)
 651             {
 652                 avem->setXAxisValue(i, surfaceSel_.position(i).atomIndices()[0] + 1);
 653             }
 654             atomArea_.addModule(avem);
 655             if (!fnAtomArea_.empty())
 656             {
 657                 AnalysisDataPlotModulePointer plotm(
 658                         new AnalysisDataPlotModule(settings.plotSettings()));
 659                 plotm->setFileName(fnAtomArea_);
 660                 plotm->setTitle("Area per residue over the trajectory");
 661                 plotm->setXLabel("Atom");
 662                 plotm->setXFormat(8, 0);
 663                 plotm->setYLabel("Area (nm\\S2\\N)");
 664                 plotm->setErrorsAsSeparateColumn(true);
 665                 plotm->appendLegend("Average (nm\\S2\\N)");
 666                 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
 667                 avem->addModule(plotm);
 668             }
 669         }
 670         {
 671             AnalysisDataAverageModulePointer avem(new AnalysisDataAverageModule);
 672             for (int i = 0; i < surfaceSel_.posCount(); ++i)
 673             {
 674                 const int atomIndex     = surfaceSel_.position(i).atomIndices()[0];
 675                 const int residueIndex  = atoms.atom[atomIndex].resind;
 676                 avem->setXAxisValue(i, atoms.resinfo[residueIndex].nr);
 677             }
 678             residueArea_.addModule(avem);
 679             if (!fnResidueArea_.empty())
 680             {
 681                 AnalysisDataPlotModulePointer plotm(
 682                         new AnalysisDataPlotModule(settings.plotSettings()));
 683                 plotm->setFileName(fnResidueArea_);
 684                 plotm->setTitle("Area per atom over the trajectory");
 685                 plotm->setXLabel("Residue");
 686                 plotm->setXFormat(8, 0);
 687                 plotm->setYLabel("Area (nm\\S2\\N)");
 688                 plotm->setErrorsAsSeparateColumn(true);
 689                 plotm->appendLegend("Average (nm\\S2\\N)");
 690                 plotm->appendLegend("Standard deviation (nm\\S2\\N)");
 691                 avem->addModule(plotm);
 692             }
 693         }
 694     }
 695
 696     if (!fnDGSolv_.empty())
 697     {
 698         dgSolv_.setColumnCount(0, 1 + outputSel_.size());
 699         AnalysisDataPlotModulePointer plotm(
 700                 new AnalysisDataPlotModule(settings.plotSettings()));
 701         plotm->setFileName(fnDGSolv_);
 702         plotm->setTitle("Free Energy of Solvation");
 703         plotm->setXAxisIsTime();
 704         plotm->setYLabel("D Gsolv");
 705         plotm->appendLegend("Total");
 706         for (size_t i = 0; i < outputSel_.size(); ++i)
 707         {
 708             plotm->appendLegend(outputSel_[i].name());
 709         }
 710         dgSolv_.addModule(plotm);
 711     }
 712
 713     if (!fnVolume_.empty())
 714     {
 715         volume_.setColumnCount(0, 2);
 716         AnalysisDataPlotModulePointer plotm(
 717                 new AnalysisDataPlotModule(settings.plotSettings()));
 718         plotm->setFileName(fnVolume_);
 719         plotm->setTitle("Volume and Density");
 720         plotm->setXAxisIsTime();
 721         plotm->appendLegend("Volume (nm\\S3\\N)");
 722         plotm->appendLegend("Density (g/l)");
 723         volume_.addModule(plotm);
 724     }
 725 }
 726
 727 /*! \brief
 728  * Temporary memory for use within a single-frame calculation.
 729  */
 730 class SasaModuleData : public TrajectoryAnalysisModuleData
 731 {
 732     public:
 733         /*! \brief
 734          * Reserves memory for the frame-local data.
 735          *
 736          * `residueCount` will be zero if per-residue data is not being
 737          * calculated.
 738          */
 739         SasaModuleData(TrajectoryAnalysisModule          *module,
 740                        const AnalysisDataParallelOptions &opt,
 741                        const SelectionCollection         &selections,
 742                        int atomCount, int residueCount)
 743             : TrajectoryAnalysisModuleData(module, opt, selections)
 744         {
 745             index_.reserve(atomCount);
 746             // If the calculation group is not dynamic, pre-calculate
 747             // the index, since it is not going to change.
 748             for (int i = 0; i < atomCount; ++i)
 749             {
 750                 index_.push_back(i);
 751             }
 752             atomAreas_.resize(atomCount);
 753             res_a_.resize(residueCount);
 754         }
 755
 756         virtual void finish() { finishDataHandles(); }
 757
 758         //! Indices of the calculation selection positions selected for the frame.
 759         std::vector<int>        index_;
 760         /*! \brief
 761          * Atom areas for each calculation selection position for the frame.
 762          *
 763          * One entry for each position in the calculation group.
 764          * Values for atoms not selected are set to zero.
 765          */
 766         std::vector<real>       atomAreas_;
 767         /*! \brief
 768          * Working array to accumulate areas for each residue.
 769          *
 770          * One entry for each distinct residue in the calculation group;
 771          * indices are not directly residue numbers or residue indices.
 772          *
 773          * This vector is empty if residue area calculations are not being
 774          * performed.
 775          */
 776         std::vector<real>       res_a_;
 777 };
 778
 779 TrajectoryAnalysisModuleDataPointer Sasa::startFrames(
 780         const AnalysisDataParallelOptions &opt,
 781         const SelectionCollection         &selections)
 782 {
 783     return TrajectoryAnalysisModuleDataPointer(
 784             new SasaModuleData(this, opt, selections, surfaceSel_.posCount(),
 785                                residueArea_.columnCount(0)));
 786 }
 787
 788 /*! \brief
 789  * Helper method to compute the areas for a single selection.
 790  *
 791  * \param[in]  surfaceSel     The calculation selection.
 792  * \param[in]  sel            The selection to compute the areas for (can be
 793  *     `surfaceSel` or one of the output selections).
 794  * \param[in]  atomAreas      Atom areas for each position in `surfaceSel`.
 795  * \param[in]  dgsFactor      Free energy coefficients for each position in
 796  *     `surfaceSel`. If empty, free energies are not calculated.
 797  * \param[out] totalAreaOut Total area of `sel` (sum of atom areas it selects).
 798  * \param[out] dgsolvOut      Solvation free energy.
 799  *     Will be zero of `dgsFactor` is empty.
 800  * \param      atomAreaHandle Data handle to use for storing atom areas for `sel`.
 801  * \param      resAreaHandle  Data handle to use for storing residue areas for `sel`.
 802  * \param      resAreaWork    Work array for accumulating the residue areas.
 803  *     If empty, atom and residue areas are not calculated.
 804  *
 805  * `atomAreaHandle` and `resAreaHandle` are not used if `resAreaWork` is empty.
 806  */
 807 void computeAreas(const Selection &surfaceSel, const Selection &sel,
 808                   const std::vector<real> &atomAreas,
 809                   const std::vector<real> &dgsFactor,
 810                   real *totalAreaOut, real *dgsolvOut,
 811                   AnalysisDataHandle atomAreaHandle,
 812                   AnalysisDataHandle resAreaHandle,
 813                   std::vector<real> *resAreaWork)
 814 {
 815     const bool bResAt    = !resAreaWork->empty();
 816     const bool bDGsolv   = !dgsFactor.empty();
 817     real       totalArea = 0;
 818     real       dgsolv    = 0;
 819
 820     if (bResAt)
 821     {
 822         std::fill(resAreaWork->begin(), resAreaWork->end(),
 823                   static_cast<real>(0.0));
 824     }
 825     for (int i = 0; i < sel.posCount(); ++i)
 826     {
 827         // Get the index of the atom in the calculation group.
 828         // For the output groups, the mapping has been precalculated in
 829         // initAnalysis().
 830         const int  ii = (sel != surfaceSel ? sel.position(i).mappedId() : i);
 831         if (!surfaceSel.position(ii).selected())
 832         {
 833             // For the calculation group, skip unselected atoms.
 834             if (sel == surfaceSel)
 835             {
 836                 continue;
 837             }
 838             GMX_THROW(InconsistentInputError("Output selection is not a subset of the surface selection"));
 839         }
 840         // Get the internal index of the matching residue.
 841         // These have been precalculated in initAnalysis().
 842         const int  ri       = surfaceSel.position(ii).mappedId();
 843         const real atomArea = atomAreas[ii];
 844         totalArea += atomArea;
 845         if (bResAt)
 846         {
 847             atomAreaHandle.setPoint(ii, atomArea);
 848             (*resAreaWork)[ri] += atomArea;
 849         }
 850         if (bDGsolv)
 851         {
 852             dgsolv += atomArea * dgsFactor[ii];
 853         }
 854     }
 855     if (bResAt)
 856     {
 857         for (size_t i = 0; i < (*resAreaWork).size(); ++i)
 858         {
 859             resAreaHandle.setPoint(i, (*resAreaWork)[i]);
 860         }
 861     }
 862     *totalAreaOut = totalArea;
 863     *dgsolvOut    = dgsolv;
 864 }
 865
 866 void
 867 Sasa::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
 868                    TrajectoryAnalysisModuleData *pdata)
 869 {
 870     AnalysisDataHandle   ah         = pdata->dataHandle(area_);
 871     AnalysisDataHandle   dgh        = pdata->dataHandle(dgSolv_);
 872     AnalysisDataHandle   aah        = pdata->dataHandle(atomArea_);
 873     AnalysisDataHandle   rah        = pdata->dataHandle(residueArea_);
 874     AnalysisDataHandle   vh         = pdata->dataHandle(volume_);
 875     const Selection     &surfaceSel = pdata->parallelSelection(surfaceSel_);
 876     const SelectionList &outputSel  = pdata->parallelSelections(outputSel_);
 877     SasaModuleData      &frameData  = *static_cast<SasaModuleData *>(pdata);
 878
 879     const bool           bResAt    = !frameData.res_a_.empty();
 880     const bool           bDGsol    = !dgsFactor_.empty();
 881     const bool           bConnolly = (frnr == 0 && !fnConnolly_.empty());
 882
 883     // Update indices of selected atoms in the work array.
 884     if (surfaceSel.isDynamic())
 885     {
 886         frameData.index_.clear();
 887         for (int i = 0; i < surfaceSel.posCount(); ++i)
 888         {
 889             if (surfaceSel.position(i).selected())
 890             {
 891                 frameData.index_.push_back(i);
 892             }
 893         }
 894     }
 895
 896     // Determine what needs to be calculated.
 897     int                  flag      = 0;
 898     if (bResAt || bDGsol || !outputSel.empty())
 899     {
 900         flag |= FLAG_ATOM_AREA;
 901     }
 902     if (bConnolly)
 903     {
 904         flag |= FLAG_DOTS;
 905     }
 906     if (volume_.columnCount() > 0)
 907     {
 908         flag |= FLAG_VOLUME;
 909     }
 910
 911     // Do the low-level calculation.
 912     // totarea and totvolume receive the values for the calculation group.
 913     // area array contains the per-atom areas for the selected positions.
 914     // surfacedots contains nsurfacedots entries, and contains the actual
 915     // points.
 916     real  totarea, totvolume;
 917     real *area = NULL, *surfacedots = NULL;
 918     int   nsurfacedots;
 919     int   retval = nsc_dclm_pbc(surfaceSel.coordinates().data(), &radii_[0],
 920                                 frameData.index_.size(), ndots_, flag, &totarea,
 921                                 &area, &totvolume, &surfacedots, &nsurfacedots,
 922                                 &frameData.index_[0],
 923                                 pbc != NULL ? pbc->ePBC : epbcNONE,
 924                                 pbc != NULL ? pbc->box : NULL);
 925     // Unpack the atomwise areas into the frameData.atomAreas_ array for easier
 926     // indexing in the case of dynamic surfaceSel.
 927     if (area != NULL)
 928     {
 929         if (surfaceSel.isDynamic())
 930         {
 931             std::fill(frameData.atomAreas_.begin(), frameData.atomAreas_.end(),
 932                       static_cast<real>(0.0));
 933             for (size_t i = 0; i < frameData.index_.size(); ++i)
 934             {
 935                 frameData.atomAreas_[frameData.index_[i]] = area[i];
 936             }
 937         }
 938         else
 939         {
 940             std::copy(area, area + surfaceSel.posCount(),
 941                       frameData.atomAreas_.begin());
 942         }
 943         sfree(area);
 944     }
 945     scoped_ptr_sfree dotsGuard(surfacedots);
 946     if (retval != 0)
 947     {
 948         GMX_THROW(InternalError("nsc_dclm_pbc failed"));
 949     }
 950
 951     if (bConnolly)
 952     {
 953         // This is somewhat nasty, as it modifies the atoms and symtab
 954         // structures.  But since it is only used in the first frame, and no
 955         // one else uses the topology after initialization, it may just work
 956         // even with future parallelization.
 957         connolly_plot(fnConnolly_.c_str(),
 958                       nsurfacedots, surfacedots, fr.x, &top_->atoms,
 959                       &top_->symtab, fr.ePBC, fr.box, bIncludeSolute_);
 960     }
 961
 962     ah.startFrame(frnr, fr.time);
 963     if (bResAt)
 964     {
 965         aah.startFrame(frnr, fr.time);
 966         rah.startFrame(frnr, fr.time);
 967     }
 968     if (bDGsol)
 969     {
 970         dgh.startFrame(frnr, fr.time);
 971     }
 972
 973     ah.setPoint(0, totarea);
 974
 975     real totalArea, dgsolv;
 976     if (bResAt || bDGsol)
 977     {
 978         computeAreas(surfaceSel, surfaceSel, frameData.atomAreas_, dgsFactor_,
 979                      &totalArea, &dgsolv, aah, rah, &frameData.res_a_);
 980         if (bDGsol)
 981         {
 982             dgh.setPoint(0, dgsolv);
 983         }
 984     }
 985     for (size_t g = 0; g < outputSel.size(); ++g)
 986     {
 987         if (bResAt)
 988         {
 989             aah.selectDataSet(g + 1);
 990             rah.selectDataSet(g + 1);
 991         }
 992         computeAreas(surfaceSel, outputSel[g], frameData.atomAreas_, dgsFactor_,
 993                      &totalArea, &dgsolv, aah, rah, &frameData.res_a_);
 994         ah.setPoint(g + 1, totalArea);
 995         if (bDGsol)
 996         {
 997             dgh.setPoint(g + 1, dgsolv);
 998         }
 999     }
1000
1001     ah.finishFrame();
1002     if (bResAt)
1003     {
1004         aah.finishFrame();
1005         rah.finishFrame();
1006     }
1007     if (bDGsol)
1008     {
1009         dgh.finishFrame();
1010     }
1011
1012     if (vh.isValid())
1013     {
1014         real totmass = 0;
1015         for (int i = 0; i < surfaceSel.posCount(); ++i)
1016         {
1017             totmass += surfaceSel.position(i).mass();
1018         }
1019         const real density = totmass*AMU/(totvolume*NANO*NANO*NANO);
1020         vh.startFrame(frnr, fr.time);
1021         vh.setPoint(0, totvolume);
1022         vh.setPoint(1, density);
1023         vh.finishFrame();
1024     }
1025 }
1026
1027 void
1028 Sasa::finishAnalysis(int /*nframes*/)
1029 {
1030     //if (bITP)
1031     //{
1032     //    fp3 = ftp2FILE(efITP, nfile, fnm, "w");
1033     //    fprintf(fp3, "[ position_restraints ]\n"
1034     //            "#define FCX 1000\n"
1035     //            "#define FCY 1000\n"
1036     //            "#define FCZ 1000\n"
1037     //            "; Atom  Type  fx   fy   fz\n");
1038     //    for (i = 0; i < nx[0]; i++)
1039     //    {
1040     //        if (atom_area[i] > minarea)
1041     //        {
1042     //            fprintf(fp3, "%5d   1     FCX  FCX  FCZ\n", ii+1);
1043     //        }
1044     //    }
1045     //    ffclose(fp3);
1046     //}
1047 }
1048
1049 void
1050 Sasa::writeOutput()
1051 {
1052 }
1053
1054 //! \}
1055
1056 }       // namespace
1057
1058 const char SasaInfo::name[]             = "sasa";
1059 const char SasaInfo::shortDescription[] =
1060     "Compute solvent accessible surface area";
1061
1062 TrajectoryAnalysisModulePointer SasaInfo::create()
1063 {
1064     return TrajectoryAnalysisModulePointer(new Sasa);
1065 }
1066
1067 } // namespace analysismodules
1068
1069 } // namespace gmx