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