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