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