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