src/gromacs/trajectoryanalysis/modules/sasa.cpp

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