src/gromacs/gmxana/gmx_do_dssp.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-2004, The GROMACS development team.
   6  * Copyright (c) 2012,2013,2014,2015,2017,2018,2019, 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 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45 #include <numeric>
  46
  47 #include "gromacs/commandline/pargs.h"
  48 #include "gromacs/commandline/viewit.h"
  49 #include "gromacs/fileio/confio.h"
  50 #include "gromacs/fileio/matio.h"
  51 #include "gromacs/fileio/pdbio.h"
  52 #include "gromacs/fileio/trxio.h"
  53 #include "gromacs/fileio/xvgr.h"
  54 #include "gromacs/gmxana/gmx_ana.h"
  55 #include "gromacs/gmxana/gstat.h"
  56 #include "gromacs/pbcutil/rmpbc.h"
  57 #include "gromacs/topology/index.h"
  58 #include "gromacs/topology/topology.h"
  59 #include "gromacs/utility/arraysize.h"
  60 #include "gromacs/utility/cstringutil.h"
  61 #include "gromacs/utility/dir_separator.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/futil.h"
  64 #include "gromacs/utility/smalloc.h"
  65 #include "gromacs/utility/strdb.h"
  66
  67 #if GMX_NATIVE_WINDOWS
  68     #define NULL_DEVICE  "nul"
  69     #define popen  _popen
  70     #define pclose _pclose
  71 #else
  72     #define NULL_DEVICE  "/dev/null"
  73 #endif
  74
  75 struct DsspInputStrings
  76 {
  77     std::string dptr;
  78     std::string pdbfile;
  79     std::string tmpfile;
  80 };
  81
  82 static const char* prepareToRedirectStdout(bool bVerbose)
  83 {
  84     return bVerbose ? "" : "2>" NULL_DEVICE;
  85 }
  86
  87 static void printDsspResult(char *dssp, const DsspInputStrings &strings,
  88                             const std::string &redirectionString)
  89 {
  90 #if HAVE_PIPES || GMX_NATIVE_WINDOWS
  91     sprintf(dssp, "%s -i %s %s",
  92             strings.dptr.c_str(), strings.pdbfile.c_str(), redirectionString.c_str());
  93 #else
  94     sprintf(dssp, "%s -i %s -o %s > %s %s",
  95             strings.dptr.c_str(), strings.pdbfile.c_str(), strings.tmpfile.c_str(), NULL_DEVICE, redirectionString.c_str());
  96 #endif
  97 }
  98
  99
 100 static int strip_dssp(FILE *tapeout, int nres,
 101                       const gmx_bool bPhobres[], real t,
 102                       real *acc, FILE *fTArea,
 103                       t_matrix *mat, int average_area[],
 104                       const gmx_output_env_t *oenv)
 105 {
 106     static gmx_bool bFirst = TRUE;
 107     static char    *ssbuf;
 108     char            buf[STRLEN+1];
 109     char            SSTP;
 110     int             nr, iacc, nresidues;
 111     int             naccf, naccb; /* Count hydrophobic and hydrophilic residues */
 112     real            iaccf, iaccb;
 113
 114
 115     /* Skip header */
 116     do
 117     {
 118         fgets2(buf, STRLEN, tapeout);
 119     }
 120     while (std::strstr(buf, "KAPPA") == nullptr);
 121     if (bFirst)
 122     {
 123         /* Since we also have empty lines in the dssp output (temp) file,
 124          * and some line content is saved to the ssbuf variable,
 125          * we need more memory than just nres elements. To be shure,
 126          * we allocate 2*nres-1, since for each chain there is a
 127          * separating line in the temp file. (At most each residue
 128          * could have been defined as a separate chain.) */
 129         snew(ssbuf, 2*nres-1);
 130     }
 131
 132     iaccb     = iaccf = 0;
 133     nresidues = 0;
 134     naccf     = 0;
 135     naccb     = 0;
 136     for (nr = 0; (fgets2(buf, STRLEN, tapeout) != nullptr); nr++)
 137     {
 138         if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
 139         {
 140             SSTP = '=';     /* Chain separator sign '=' */
 141         }
 142         else
 143         {
 144             SSTP = buf[16] == ' ' ? '~' : buf[16];
 145         }
 146         ssbuf[nr] = SSTP;
 147
 148         buf[39] = '\0';
 149
 150         /* Only calculate solvent accessible area if needed */
 151         if ((nullptr != acc) && (buf[13] != '!'))
 152         {
 153             sscanf(&(buf[34]), "%d", &iacc);
 154             acc[nr] = iacc;
 155             /* average_area and bPhobres are counted from 0...nres-1 */
 156             average_area[nresidues] += iacc;
 157             if (bPhobres[nresidues])
 158             {
 159                 naccb++;
 160                 iaccb += iacc;
 161             }
 162             else
 163             {
 164                 naccf++;
 165                 iaccf += iacc;
 166             }
 167             /* Keep track of the residue number (do not count chain separator lines '!') */
 168             nresidues++;
 169         }
 170     }
 171     ssbuf[nr] = '\0';
 172
 173     if (bFirst)
 174     {
 175         if (nullptr != acc)
 176         {
 177             fprintf(stderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n", naccb, naccf);
 178         }
 179
 180         mat->title     = "Secondary structure";
 181         mat->legend    = "";
 182         mat->label_x   = output_env_get_time_label(oenv);
 183         mat->label_y   = "Residue";
 184         mat->bDiscrete = true;
 185         mat->ny        = nr;
 186         mat->axis_y.resize(nr);
 187         std::iota(mat->axis_y.begin(), mat->axis_y.end(), 1);
 188         mat->axis_x.resize(0);
 189         mat->matrix.resize(0, 0);
 190         bFirst      = false;
 191     }
 192     mat->axis_x.push_back(t);
 193     mat->matrix.resize(mat->matrix.extent(0), nr);
 194     mat->nx = mat->matrix.extent(0);
 195     auto columnIndex = mat->nx - 1;
 196     for (int i = 0; i < nr; i++)
 197     {
 198         t_xpmelmt c = {ssbuf[i], 0};
 199         mat->matrix(columnIndex, i) = std::max(static_cast<t_matelmt>(0), searchcmap(mat->map, c));
 200     }
 201
 202     if (fTArea)
 203     {
 204         fprintf(fTArea, "%10g  %10g  %10g\n", t, 0.01*iaccb, 0.01*iaccf);
 205     }
 206
 207     /* Return the number of lines found in the dssp file (i.e. number
 208      * of redidues plus chain separator lines).
 209      * This is the number of y elements needed for the area xpm file */
 210     return nr;
 211 }
 212
 213 static gmx_bool *bPhobics(t_atoms *atoms)
 214 {
 215     int         j, i, nb;
 216     char      **cb;
 217     gmx_bool   *bb;
 218     int         n_surf;
 219     char        surffn[] = "surface.dat";
 220     char      **surf_res, **surf_lines;
 221
 222
 223     nb = get_lines("phbres.dat", &cb);
 224     snew(bb, atoms->nres);
 225
 226     n_surf = get_lines(surffn, &surf_lines);
 227     snew(surf_res, n_surf);
 228     for (i = 0; (i < n_surf); i++)
 229     {
 230         snew(surf_res[i], 5);
 231         sscanf(surf_lines[i], "%s", surf_res[i]);
 232     }
 233
 234
 235     for (i = 0, j = 0; (i < atoms->nres); i++)
 236     {
 237         if (-1 != search_str(n_surf, surf_res, *atoms->resinfo[i].name) )
 238         {
 239             bb[j++] = (-1 != search_str(nb, cb, *atoms->resinfo[i].name));
 240         }
 241     }
 242
 243     if (i != j)
 244     {
 245         fprintf(stderr, "Not all residues were recognized (%d from %d), the result may be inaccurate!\n", j, i);
 246     }
 247
 248     for (i = 0; (i < n_surf); i++)
 249     {
 250         sfree(surf_res[i]);
 251     }
 252     sfree(surf_res);
 253
 254     return bb;
 255 }
 256
 257 static void check_oo(t_atoms *atoms)
 258 {
 259     char *OOO;
 260
 261     int   i;
 262
 263     OOO = gmx_strdup("O");
 264
 265     for (i = 0; (i < atoms->nr); i++)
 266     {
 267         if (std::strcmp(*(atoms->atomname[i]), "OXT") == 0)
 268         {
 269             *atoms->atomname[i] = OOO;
 270         }
 271         else if (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
 272         {
 273             *atoms->atomname[i] = OOO;
 274         }
 275         else if (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
 276         {
 277             *atoms->atomname[i] = OOO;
 278         }
 279     }
 280 }
 281
 282 static void norm_acc(t_atoms *atoms, int nres,
 283                      const real av_area[], real norm_av_area[])
 284 {
 285     int     i, n, n_surf;
 286
 287     char    surffn[] = "surface.dat";
 288     char  **surf_res, **surf_lines;
 289     double *surf;
 290
 291     n_surf = get_lines(surffn, &surf_lines);
 292     snew(surf, n_surf);
 293     snew(surf_res, n_surf);
 294     for (i = 0; (i < n_surf); i++)
 295     {
 296         snew(surf_res[i], 5);
 297         sscanf(surf_lines[i], "%s %lf", surf_res[i], &surf[i]);
 298     }
 299
 300     for (i = 0; (i < nres); i++)
 301     {
 302         n = search_str(n_surf, surf_res, *atoms->resinfo[i].name);
 303         if (n != -1)
 304         {
 305             norm_av_area[i] = av_area[i] / surf[n];
 306         }
 307         else
 308         {
 309             fprintf(stderr, "Residue %s not found in surface database (%s)\n",
 310                     *atoms->resinfo[i].name, surffn);
 311         }
 312     }
 313 }
 314
 315 static void prune_ss_legend(t_matrix *mat)
 316 {
 317     std::vector<bool>      isPresent(mat->map.size());
 318     std::vector<int>       newnum(mat->map.size());
 319     std::vector<t_mapping> newmap;
 320
 321     for (int f = 0; f < mat->nx; f++)
 322     {
 323         for (int r = 0; r < mat->ny; r++)
 324         {
 325             isPresent[mat->matrix(f, r)] = true;
 326         }
 327     }
 328
 329     for (size_t i = 0; i < mat->map.size(); i++)
 330     {
 331         newnum[i] = -1;
 332         if (isPresent[i])
 333         {
 334             newnum[i] = newmap.size();
 335             newmap.emplace_back(mat->map[i]);
 336         }
 337     }
 338     if (newmap.size() != mat->map.size())
 339     {
 340         std::swap(mat->map, newmap);
 341         for (int f = 0; f < mat->nx; f++)
 342         {
 343             for (int r = 0; r < mat->ny; r++)
 344             {
 345                 mat->matrix(f, r) = newnum[mat->matrix(f, r)];
 346             }
 347         }
 348     }
 349 }
 350
 351 static void write_sas_mat(const char *fn, real **accr, int nframe, int nres, t_matrix *mat)
 352 {
 353     real  lo, hi;
 354     int   i, j, nlev;
 355     t_rgb rlo = {1, 1, 1}, rhi = {0, 0, 0};
 356     FILE *fp;
 357
 358     if (fn)
 359     {
 360         hi = lo = accr[0][0];
 361         for (i = 0; i < nframe; i++)
 362         {
 363             for (j = 0; j < nres; j++)
 364             {
 365                 lo = std::min(lo, accr[i][j]);
 366                 hi = std::max(hi, accr[i][j]);
 367             }
 368         }
 369         fp   = gmx_ffopen(fn, "w");
 370         nlev = static_cast<int>(hi-lo+1);
 371         write_xpm(fp, 0, "Solvent Accessible Surface", "Surface (A^2)",
 372                   "Time", "Residue Index", nframe, nres,
 373                   mat->axis_x.data(), mat->axis_y.data(), accr, lo, hi, rlo, rhi, &nlev);
 374         gmx_ffclose(fp);
 375     }
 376 }
 377
 378 static void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,
 379                        const gmx_output_env_t *oenv)
 380 {
 381     FILE                    *fp;
 382     int                      ss_count, total_count;
 383
 384     gmx::ArrayRef<t_mapping> map = mat->map;
 385     std::vector<int>         count(map.size());
 386     std::vector<int>         total(map.size(), 0);
 387     // This copying would not be necessary if xvgr_legend could take a
 388     // view of string views
 389     std::vector<std::string> leg;
 390     leg.reserve(map.size()+1);
 391     leg.emplace_back("Structure");
 392     for (const auto &m : map)
 393     {
 394         leg.emplace_back(m.desc);
 395     }
 396
 397     fp = xvgropen(outfile, "Secondary Structure",
 398                   output_env_get_xvgr_tlabel(oenv), "Number of Residues", oenv);
 399     if (output_env_get_print_xvgr_codes(oenv))
 400     {
 401         fprintf(fp, "@ subtitle \"Structure = ");
 402     }
 403     for (size_t s = 0; s < std::strlen(ss_string); s++)
 404     {
 405         if (s > 0)
 406         {
 407             fprintf(fp, " + ");
 408         }
 409         for (const auto &m : map)
 410         {
 411             if (ss_string[s] == m.code.c1)
 412             {
 413                 fprintf(fp, "%s", m.desc);
 414             }
 415         }
 416     }
 417     fprintf(fp, "\"\n");
 418     xvgrLegend(fp, leg, oenv);
 419
 420     total_count = 0;
 421     for (int f = 0; f < mat->nx; f++)
 422     {
 423         ss_count = 0;
 424         for (auto &c : count)
 425         {
 426             c = 0;
 427         }
 428         for (int r = 0; r < mat->ny; r++)
 429         {
 430             count[mat->matrix(f, r)]++;
 431             total[mat->matrix(f, r)]++;
 432         }
 433         for (gmx::index s = 0; s != gmx::ssize(map); ++s)
 434         {
 435             if (std::strchr(ss_string, map[s].code.c1))
 436             {
 437                 ss_count    += count[s];
 438                 total_count += count[s];
 439             }
 440         }
 441         fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);
 442         for (const auto &c : count)
 443         {
 444             fprintf(fp, " %5d", c);
 445         }
 446         fprintf(fp, "\n");
 447     }
 448     /* now print column totals */
 449     fprintf(fp, "%-8s %5d", "# Totals", total_count);
 450     for (const auto &t : total)
 451     {
 452         fprintf(fp, " %5d", t);
 453     }
 454     fprintf(fp, "\n");
 455
 456     /* now print probabilities */
 457     fprintf(fp, "%-8s %5.2f", "# SS pr.", total_count / static_cast<real>(mat->nx * mat->ny));
 458     for (const auto &t : total)
 459     {
 460         fprintf(fp, " %5.2f", t / static_cast<real>(mat->nx * mat->ny));
 461     }
 462     fprintf(fp, "\n");
 463
 464     xvgrclose(fp);
 465 }
 466
 467 int gmx_do_dssp(int argc, char *argv[])
 468 {
 469     const char        *desc[] = {
 470         "[THISMODULE] ",
 471         "reads a trajectory file and computes the secondary structure for",
 472         "each time frame ",
 473         "calling the dssp program. If you do not have the dssp program,",
 474         "get it from http://swift.cmbi.ru.nl/gv/dssp. [THISMODULE] assumes ",
 475         "that the dssp executable is located in ",
 476         "[TT]/usr/local/bin/dssp[tt]. If this is not the case, then you should",
 477         "set an environment variable [TT]DSSP[tt] pointing to the dssp",
 478         "executable, e.g.: [PAR]",
 479         "[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
 480         "Since version 2.0.0, dssp is invoked with a syntax that differs",
 481         "from earlier versions. If you have an older version of dssp,",
 482         "use the [TT]-ver[tt] option to direct do_dssp to use the older syntax.",
 483         "By default, do_dssp uses the syntax introduced with version 2.0.0.",
 484         "Even newer versions (which at the time of writing are not yet released)",
 485         "are assumed to have the same syntax as 2.0.0.[PAR]",
 486         "The structure assignment for each residue and time is written to an",
 487         "[REF].xpm[ref] matrix file. This file can be visualized with for instance",
 488         "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].",
 489         "Individual chains are separated by light grey lines in the [REF].xpm[ref] and",
 490         "postscript files.",
 491         "The number of residues with each secondary structure type and the",
 492         "total secondary structure ([TT]-sss[tt]) count as a function of",
 493         "time are also written to file ([TT]-sc[tt]).[PAR]",
 494         "Solvent accessible surface (SAS) per residue can be calculated, both in",
 495         "absolute values (A^2) and in fractions of the maximal accessible",
 496         "surface of a residue. The maximal accessible surface is defined as",
 497         "the accessible surface of a residue in a chain of glycines.",
 498         "[BB]Note[bb] that the program [gmx-sas] can also compute SAS",
 499         "and that is more efficient.[PAR]",
 500         "Finally, this program can dump the secondary structure in a special file",
 501         "[TT]ssdump.dat[tt] for usage in the program [gmx-chi]. Together",
 502         "these two programs can be used to analyze dihedral properties as a",
 503         "function of secondary structure type."
 504     };
 505     static gmx_bool    bVerbose;
 506     static const char *ss_string   = "HEBT";
 507     static int         dsspVersion = 2;
 508     t_pargs            pa[]        = {
 509         { "-v",  FALSE, etBOOL, {&bVerbose},
 510           "HIDDENGenerate miles of useless information" },
 511         { "-sss", FALSE, etSTR, {&ss_string},
 512           "Secondary structures for structure count"},
 513         { "-ver", FALSE, etINT, {&dsspVersion},
 514           "DSSP major version. Syntax changed with version 2"}
 515     };
 516
 517     t_trxstatus       *status;
 518     FILE              *tapein,  *tapeout;
 519     FILE              *ss, *acc, *fTArea, *tmpf;
 520     const char        *fnSCount, *fnArea, *fnTArea, *fnAArea;
 521     const char        *leg[] = { "Phobic", "Phylic" };
 522     t_topology         top;
 523     int                ePBC;
 524     t_atoms           *atoms;
 525     t_matrix           mat;
 526     int                nres, nr0, naccr, nres_plus_separators;
 527     gmx_bool          *bPhbres, bDoAccSurf;
 528     real               t;
 529     int                natoms, nframe = 0;
 530     matrix             box = {{0}};
 531     int                gnx;
 532     char              *grpnm;
 533     int               *index;
 534     rvec              *xp, *x;
 535     int               *average_area;
 536     real             **accr, *accr_ptr = nullptr, *av_area, *norm_av_area;
 537     char               pdbfile[32], tmpfile[32];
 538     char               dssp[256];
 539     const char        *dptr;
 540     gmx_output_env_t  *oenv;
 541     gmx_rmpbc_t        gpbc = nullptr;
 542
 543     t_filenm           fnm[] = {
 544         { efTRX, "-f",   nullptr,      ffREAD },
 545         { efTPS, nullptr,   nullptr,      ffREAD },
 546         { efNDX, nullptr,   nullptr,      ffOPTRD },
 547         { efDAT, "-ssdump", "ssdump", ffOPTWR },
 548         { efMAP, "-map", "ss",      ffLIBRD },
 549         { efXPM, "-o",   "ss",      ffWRITE },
 550         { efXVG, "-sc",  "scount",  ffWRITE },
 551         { efXPM, "-a",   "area",    ffOPTWR },
 552         { efXVG, "-ta",  "totarea", ffOPTWR },
 553         { efXVG, "-aa",  "averarea", ffOPTWR }
 554     };
 555 #define NFILE asize(fnm)
 556
 557     if (!parse_common_args(&argc, argv,
 558                            PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT,
 559                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
 560     {
 561         return 0;
 562     }
 563     fnSCount   = opt2fn("-sc", NFILE, fnm);
 564     fnArea     = opt2fn_null("-a", NFILE, fnm);
 565     fnTArea    = opt2fn_null("-ta", NFILE, fnm);
 566     fnAArea    = opt2fn_null("-aa", NFILE, fnm);
 567     bDoAccSurf = ((fnArea != nullptr) || (fnTArea != nullptr) || (fnAArea != nullptr));
 568
 569     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xp, nullptr, box, FALSE);
 570     atoms = &(top.atoms);
 571     check_oo(atoms);
 572     bPhbres = bPhobics(atoms);
 573
 574     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpnm);
 575     nres = 0;
 576     nr0  = -1;
 577     for (int i = 0; (i < gnx); i++)
 578     {
 579         if (atoms->atom[index[i]].resind != nr0)
 580         {
 581             nr0 = atoms->atom[index[i]].resind;
 582             nres++;
 583         }
 584     }
 585     fprintf(stderr, "There are %d residues in your selected group\n", nres);
 586
 587     std::strcpy(pdbfile, "ddXXXXXX");
 588     gmx_tmpnam(pdbfile);
 589     if ((tmpf = fopen(pdbfile, "w")) == nullptr)
 590     {
 591         sprintf(pdbfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
 592         gmx_tmpnam(pdbfile);
 593         if ((tmpf = fopen(pdbfile, "w")) == nullptr)
 594         {
 595             gmx_fatal(FARGS, "Can not open tmp file %s", pdbfile);
 596         }
 597     }
 598     fclose(tmpf);
 599
 600     std::strcpy(tmpfile, "ddXXXXXX");
 601     gmx_tmpnam(tmpfile);
 602     if ((tmpf = fopen(tmpfile, "w")) == nullptr)
 603     {
 604         sprintf(tmpfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
 605         gmx_tmpnam(tmpfile);
 606         if ((tmpf = fopen(tmpfile, "w")) == nullptr)
 607         {
 608             gmx_fatal(FARGS, "Can not open tmp file %s", tmpfile);
 609         }
 610     }
 611     fclose(tmpf);
 612
 613     if ((dptr = getenv("DSSP")) == nullptr)
 614     {
 615         dptr = "/usr/local/bin/dssp";
 616     }
 617     if (!gmx_fexist(dptr))
 618     {
 619         gmx_fatal(FARGS, "DSSP executable (%s) does not exist (use setenv DSSP)",
 620                   dptr);
 621     }
 622     std::string      redirectionString;
 623     redirectionString = prepareToRedirectStdout(bVerbose);
 624     DsspInputStrings dsspStrings;
 625     dsspStrings.dptr    = dptr;
 626     dsspStrings.pdbfile = pdbfile;
 627     dsspStrings.tmpfile = tmpfile;
 628     if (dsspVersion >= 2)
 629     {
 630         if (dsspVersion > 2)
 631         {
 632             printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by do_dssp. Assuming version 2 syntax.\n\n", dsspVersion);
 633         }
 634
 635         printDsspResult(dssp, dsspStrings, redirectionString);
 636     }
 637     else
 638     {
 639         if (bDoAccSurf)
 640         {
 641             dsspStrings.dptr.clear();
 642         }
 643         else
 644         {
 645             dsspStrings.dptr = "-na";
 646         }
 647         printDsspResult(dssp, dsspStrings, redirectionString);
 648     }
 649     fprintf(stderr, "dssp cmd='%s'\n", dssp);
 650
 651     if (fnTArea)
 652     {
 653         fTArea = xvgropen(fnTArea, "Solvent Accessible Surface Area",
 654                           output_env_get_xvgr_tlabel(oenv), "Area (nm\\S2\\N)", oenv);
 655         xvgr_legend(fTArea, 2, leg, oenv);
 656     }
 657     else
 658     {
 659         fTArea = nullptr;
 660     }
 661
 662     mat.map = readcmap(opt2fn("-map", NFILE, fnm));
 663
 664     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 665     if (natoms > atoms->nr)
 666     {
 667         gmx_fatal(FARGS, "\nTrajectory does not match topology!");
 668     }
 669     if (gnx > natoms)
 670     {
 671         gmx_fatal(FARGS, "\nTrajectory does not match selected group!");
 672     }
 673
 674     snew(average_area, atoms->nres);
 675     snew(av_area, atoms->nres);
 676     snew(norm_av_area, atoms->nres);
 677     accr  = nullptr;
 678     naccr = 0;
 679
 680     gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
 681     do
 682     {
 683         t = output_env_conv_time(oenv, t);
 684         if (bDoAccSurf && nframe >= naccr)
 685         {
 686             naccr += 10;
 687             srenew(accr, naccr);
 688             for (int i = naccr-10; i < naccr; i++)
 689             {
 690                 snew(accr[i], 2*atoms->nres-1);
 691             }
 692         }
 693         gmx_rmpbc(gpbc, natoms, box, x);
 694         tapein = gmx_ffopen(pdbfile, "w");
 695         write_pdbfile_indexed(tapein, nullptr, atoms, x, ePBC, box, ' ', -1, gnx, index, nullptr, FALSE);
 696         gmx_ffclose(tapein);
 697         /* strip_dssp returns the number of lines found in the dssp file, i.e.
 698          * the number of residues plus the separator lines */
 699
 700 #if HAVE_PIPES || GMX_NATIVE_WINDOWS
 701         if (nullptr == (tapeout = popen(dssp, "r")))
 702 #else
 703         if (0 != system(dssp) || nullptr == (tapeout = gmx_ffopen(tmpfile, "r")))
 704 #endif
 705         {
 706             remove(pdbfile);
 707             remove(tmpfile);
 708             gmx_fatal(FARGS, "Failed to execute command: %s\n"
 709                       "Try specifying your dssp version with the -ver option.", dssp);
 710         }
 711         if (bDoAccSurf)
 712         {
 713             accr_ptr = accr[nframe];
 714         }
 715         /* strip_dssp returns the number of lines found in the dssp file, i.e.
 716          * the number of residues plus the separator lines */
 717         nres_plus_separators = strip_dssp(tapeout, nres, bPhbres, t,
 718                                           accr_ptr, fTArea, &mat, average_area, oenv);
 719 #if HAVE_PIPES || GMX_NATIVE_WINDOWS
 720         pclose(tapeout);
 721 #else
 722         gmx_ffclose(tapeout);
 723 #endif
 724         remove(tmpfile);
 725         remove(pdbfile);
 726         nframe++;
 727     }
 728     while (read_next_x(oenv, status, &t, x, box));
 729     fprintf(stderr, "\n");
 730     close_trx(status);
 731     if (fTArea)
 732     {
 733         xvgrclose(fTArea);
 734     }
 735     gmx_rmpbc_done(gpbc);
 736
 737     prune_ss_legend(&mat);
 738
 739     ss        = opt2FILE("-o", NFILE, fnm, "w");
 740     mat.flags = 0;
 741     write_xpm_m(ss, mat);
 742     gmx_ffclose(ss);
 743
 744     if (opt2bSet("-ssdump", NFILE, fnm))
 745     {
 746         ss = opt2FILE("-ssdump", NFILE, fnm, "w");
 747         fprintf(ss, "%d\n", nres);
 748         for (gmx::index j = 0; j != mat.matrix.extent(0); ++j)
 749         {
 750             auto row = mat.matrix.asView()[j];
 751             for (gmx::index i = 0; i != row.extent(0); ++i)
 752             {
 753                 fputc(mat.map[row[i]].code.c1, ss);
 754             }
 755             fputc('\n', ss);
 756         }
 757         gmx_ffclose(ss);
 758     }
 759     analyse_ss(fnSCount, &mat, ss_string, oenv);
 760
 761     if (bDoAccSurf)
 762     {
 763         write_sas_mat(fnArea, accr, nframe, nres_plus_separators, &mat);
 764
 765         for (int i = 0; i < atoms->nres; i++)
 766         {
 767             av_area[i] = (average_area[i] / static_cast<real>(nframe));
 768         }
 769
 770         norm_acc(atoms, nres, av_area, norm_av_area);
 771
 772         if (fnAArea)
 773         {
 774             acc = xvgropen(fnAArea, "Average Accessible Area",
 775                            "Residue", "A\\S2", oenv);
 776             for (int i = 0; (i < nres); i++)
 777             {
 778                 fprintf(acc, "%5d  %10g %10g\n", i+1, av_area[i], norm_av_area[i]);
 779             }
 780             xvgrclose(acc);
 781         }
 782     }
 783
 784     view_all(oenv, NFILE, fnm);
 785
 786     return 0;
 787 }