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