Fix link for dssp download
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_do_dssp.cpp
index cc8c87c6925a13e911f7761c991c4620aa3b89c3..36fb601e93df2e5dde55686f2c0ae6abf661f9f4 100644 (file)
@@ -3,7 +3,8 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2013,2014,2015,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
  */
 #include "gmxpre.h"
 
+#include "config.h"
+
 #include <cstdlib>
 #include <cstring>
 
 #include <algorithm>
+#include <numeric>
 
 #include "gromacs/commandline/pargs.h"
 #include "gromacs/commandline/viewit.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/strdb.h"
 
-static int strip_dssp(char *dsspfile, int nres,
-                      const gmx_bool bPhobres[], real t,
-                      real *acc, FILE *fTArea,
-                      t_matrix *mat, int average_area[],
-                      const gmx_output_env_t *oenv)
+#if GMX_NATIVE_WINDOWS
+#    define NULL_DEVICE "nul"
+#    define popen _popen
+#    define pclose _pclose
+#else
+#    define NULL_DEVICE "/dev/null"
+#endif
+
+struct DsspInputStrings
 {
-    static gmx_bool bFirst = TRUE;
-    static char    *ssbuf;
-    FILE           *tapeout;
-    static int      xsize, frame;
-    char            buf[STRLEN+1];
-    char            SSTP;
-    int             i, nr, iacc, nresidues;
-    int             naccf, naccb; /* Count hydrophobic and hydrophilic residues */
-    real            iaccf, iaccb;
-    t_xpmelmt       c;
-
-    tapeout = gmx_ffopen(dsspfile, "r");
+    std::string dptr;
+    std::string pdbfile;
+    std::string tmpfile;
+};
+
+static const char* prepareToRedirectStdout(bool bVerbose)
+{
+    return bVerbose ? "" : "2>" NULL_DEVICE;
+}
+
+static void printDsspResult(char* dssp, const DsspInputStrings& strings, const std::string& redirectionString)
+{
+#if HAVE_PIPES || GMX_NATIVE_WINDOWS
+    sprintf(dssp, "%s -i %s %s", strings.dptr.c_str(), strings.pdbfile.c_str(), redirectionString.c_str());
+#else
+    sprintf(dssp,
+            "%s -i %s -o %s > %s %s",
+            strings.dptr.c_str(),
+            strings.pdbfile.c_str(),
+            strings.tmpfile.c_str(),
+            NULL_DEVICE,
+            redirectionString.c_str());
+#endif
+}
+
+
+static int strip_dssp(FILE*                   tapeout,
+                      int                     nres,
+                      const gmx_bool          bPhobres[],
+                      real                    t,
+                      real*                   acc,
+                      FILE*                   fTArea,
+                      t_matrix*               mat,
+                      int                     average_area[],
+                      const gmx_output_env_t* oenv)
+{
+    static gmx_bool    bFirst = TRUE;
+    static std::string ssbuf;
+    char               buf[STRLEN + 1];
+    char               SSTP;
+    int                nr, iacc, nresidues;
+    int                naccf, naccb; /* Count hydrophobic and hydrophilic residues */
+    real               iaccf, iaccb;
+
 
     /* Skip header */
     do
     {
         fgets2(buf, STRLEN, tapeout);
-    }
-    while (std::strstr(buf, "KAPPA") == nullptr);
+    } while (std::strstr(buf, "KAPPA") == nullptr);
     if (bFirst)
     {
         /* Since we also have empty lines in the dssp output (temp) file,
@@ -94,18 +133,18 @@ static int strip_dssp(char *dsspfile, int nres,
          * we allocate 2*nres-1, since for each chain there is a
          * separating line in the temp file. (At most each residue
          * could have been defined as a separate chain.) */
-        snew(ssbuf, 2*nres-1);
+        ssbuf.resize(2 * nres - 1);
     }
 
-    iaccb     = iaccf = 0;
-    nresidues = 0;
-    naccf     = 0;
-    naccb     = 0;
+    iaccb = iaccf = 0;
+    nresidues     = 0;
+    naccf         = 0;
+    naccb         = 0;
     for (nr = 0; (fgets2(buf, STRLEN, tapeout) != nullptr); nr++)
     {
         if (buf[13] == '!') /* Chain separator line has '!' at pos. 13 */
         {
-            SSTP = '=';     /* Chain separator sign '=' */
+            SSTP = '='; /* Chain separator sign '=' */
         }
         else
         {
@@ -145,45 +184,31 @@ static int strip_dssp(char *dsspfile, int nres,
             fprintf(stderr, "%d residues were classified as hydrophobic and %d as hydrophilic.\n", naccb, naccf);
         }
 
-        sprintf(mat->title, "Secondary structure");
-        mat->legend[0] = 0;
-        sprintf(mat->label_x, "%s", output_env_get_time_label(oenv).c_str());
-        sprintf(mat->label_y, "Residue");
-        mat->bDiscrete = TRUE;
+        mat->title     = "Secondary structure";
+        mat->legend    = "";
+        mat->label_x   = output_env_get_time_label(oenv);
+        mat->label_y   = "Residue";
+        mat->bDiscrete = true;
         mat->ny        = nr;
-        snew(mat->axis_y, nr);
-        for (i = 0; i < nr; i++)
-        {
-            mat->axis_y[i] = i+1;
-        }
-        mat->axis_x = nullptr;
-        mat->matrix = nullptr;
-        xsize       = 0;
-        frame       = 0;
-        bFirst      = FALSE;
-    }
-    if (frame >= xsize)
-    {
-        xsize += 10;
-        srenew(mat->axis_x, xsize);
-        srenew(mat->matrix, xsize);
+        mat->axis_y.resize(nr);
+        std::iota(mat->axis_y.begin(), mat->axis_y.end(), 1);
+        mat->axis_x.resize(0);
+        mat->matrix.resize(1, 1);
+        bFirst = false;
     }
-    mat->axis_x[frame] = t;
-    snew(mat->matrix[frame], nr);
-    c.c2 = 0;
-    for (i = 0; i < nr; i++)
+    mat->axis_x.push_back(t);
+    mat->matrix.resize(++(mat->nx), nr);
+    auto columnIndex = mat->nx - 1;
+    for (int i = 0; i < nr; i++)
     {
-        c.c1                  = ssbuf[i];
-        mat->matrix[frame][i] = std::max(static_cast<t_matelmt>(0), searchcmap(mat->nmap, mat->map, c));
+        t_xpmelmt c                 = { ssbuf[i], 0 };
+        mat->matrix(columnIndex, i) = std::max(static_cast<t_matelmt>(0), searchcmap(mat->map, c));
     }
-    frame++;
-    mat->nx = frame;
 
     if (fTArea)
     {
-        fprintf(fTArea, "%10g  %10g  %10g\n", t, 0.01*iaccb, 0.01*iaccf);
+        fprintf(fTArea, "%10g  %10g  %10g\n", t, 0.01 * iaccb, 0.01 * iaccf);
     }
-    gmx_ffclose(tapeout);
 
     /* Return the number of lines found in the dssp file (i.e. number
      * of redidues plus chain separator lines).
@@ -191,59 +216,76 @@ static int strip_dssp(char *dsspfile, int nres,
     return nr;
 }
 
-static gmx_bool *bPhobics(t_atoms *atoms)
+static gmx_bool* bPhobics(t_atoms* atoms)
 {
-    int         i, nb;
-    char      **cb;
-    gmx_bool   *bb;
+    int       j, i, nb;
+    char**    cb;
+    gmx_bool* bb;
+    int       n_surf;
+    char      surffn[] = "surface.dat";
+    char **   surf_res, **surf_lines;
 
 
     nb = get_lines("phbres.dat", &cb);
     snew(bb, atoms->nres);
 
-    for (i = 0; (i < atoms->nres); i++)
+    n_surf = get_lines(surffn, &surf_lines);
+    snew(surf_res, n_surf);
+    for (i = 0; (i < n_surf); i++)
     {
-        if (-1 != search_str(nb, cb, *atoms->resinfo[i].name) )
+        snew(surf_res[i], 5);
+        sscanf(surf_lines[i], "%s", surf_res[i]);
+    }
+
+
+    for (i = 0, j = 0; (i < atoms->nres); i++)
+    {
+        if (-1 != search_str(n_surf, surf_res, *atoms->resinfo[i].name))
         {
-            bb[i] = TRUE;
+            bb[j++] = (-1 != search_str(nb, cb, *atoms->resinfo[i].name));
         }
     }
+
+    if (i != j)
+    {
+        fprintf(stderr,
+                "Not all residues were recognized (%d from %d), the result may be inaccurate!\n",
+                j,
+                i);
+    }
+
+    for (i = 0; (i < n_surf); i++)
+    {
+        sfree(surf_res[i]);
+    }
+    sfree(surf_res);
+
     return bb;
 }
 
-static void check_oo(t_atoms *atoms)
+static void check_oo(t_atomsatoms)
 {
-    char *OOO;
-
-    int   i;
+    char* OOO = gmx_strdup("O");
 
-    OOO = gmx_strdup("O");
-
-    for (i = 0; (i < atoms->nr); i++)
+    for (int i = 0; (i < atoms->nr); i++)
     {
-        if (std::strcmp(*(atoms->atomname[i]), "OXT") == 0)
-        {
-            *atoms->atomname[i] = OOO;
-        }
-        else if (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
-        {
-            *atoms->atomname[i] = OOO;
-        }
-        else if (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
+        if ((std::strcmp(*(atoms->atomname[i]), "OXT") == 0)
+            || (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
+            || (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
+            || (std::strcmp(*(atoms->atomname[i]), "OT1") == 0))
         {
             *atoms->atomname[i] = OOO;
         }
     }
 }
 
-static void norm_acc(t_atoms *atoms, int nres,
-                     const real av_area[], real norm_av_area[])
+static void norm_acc(t_atoms* atoms, int nres, const real av_area[], real norm_av_area[])
 {
-    int     i, n, n_surf;
+    int i, n, n_surf;
 
     char    surffn[] = "surface.dat";
-    char  **surf_res, **surf_lines;
-    double *surf;
+    char ** surf_res, **surf_lines;
+    doublesurf;
 
     n_surf = get_lines(surffn, &surf_lines);
     snew(surf, n_surf);
@@ -263,63 +305,53 @@ static void norm_acc(t_atoms *atoms, int nres,
         }
         else
         {
-            fprintf(stderr, "Residue %s not found in surface database (%s)\n",
-                    *atoms->resinfo[i].name, surffn);
+            fprintf(stderr, "Residue %s not found in surface database (%s)\n", *atoms->resinfo[i].name, surffn);
         }
     }
 }
 
-static void prune_ss_legend(t_matrix *mat)
+static void prune_ss_legend(t_matrixmat)
 {
-    gmx_bool  *present;
-    int       *newnum;
-    int        i, r, f, newnmap;
-    t_mapping *newmap;
+    std::vector<bool>      isPresent(mat->map.size());
+    std::vector<int>       newnum(mat->map.size());
+    std::vector<t_mapping> newmap;
 
-    snew(present, mat->nmap);
-    snew(newnum, mat->nmap);
-
-    for (f = 0; f < mat->nx; f++)
+    for (int f = 0; f < mat->nx; f++)
     {
-        for (r = 0; r < mat->ny; r++)
+        for (int r = 0; r < mat->ny; r++)
         {
-            present[mat->matrix[f][r]] = TRUE;
+            isPresent[mat->matrix(f, r)] = true;
         }
     }
 
-    newnmap = 0;
-    newmap  = nullptr;
-    for (i = 0; i < mat->nmap; i++)
+    for (size_t i = 0; i < mat->map.size(); i++)
     {
         newnum[i] = -1;
-        if (present[i])
+        if (isPresent[i])
         {
-            newnum[i] = newnmap;
-            newnmap++;
-            srenew(newmap, newnmap);
-            newmap[newnmap-1] = mat->map[i];
+            newnum[i] = newmap.size();
+            newmap.emplace_back(mat->map[i]);
         }
     }
-    if (newnmap != mat->nmap)
+    if (newmap.size() != mat->map.size())
     {
-        mat->nmap = newnmap;
-        mat->map  = newmap;
-        for (f = 0; f < mat->nx; f++)
+        std::swap(mat->map, newmap);
+        for (int f = 0; f < mat->nx; f++)
         {
-            for (r = 0; r < mat->ny; r++)
+            for (int r = 0; r < mat->ny; r++)
             {
-                mat->matrix[f][r] = newnum[mat->matrix[f][r]];
+                mat->matrix(f, r) = newnum[mat->matrix(f, r)];
             }
         }
     }
 }
 
-static void write_sas_mat(const char *fn, real **accr, int nframe, int nres, t_matrix *mat)
+static void write_sas_mat(const char* fn, real** accr, int nframe, int nres, t_matrix* mat)
 {
     real  lo, hi;
     int   i, j, nlev;
-    t_rgb rlo = {1, 1, 1}, rhi = {0, 0, 0};
-    FILE *fp;
+    t_rgb rlo = { 1, 1, 1 }, rhi = { 0, 0, 0 };
+    FILEfp;
 
     if (fn)
     {
@@ -333,119 +365,126 @@ static void write_sas_mat(const char *fn, real **accr, int nframe, int nres, t_m
             }
         }
         fp   = gmx_ffopen(fn, "w");
-        nlev = static_cast<int>(hi-lo+1);
-        write_xpm(fp, 0, "Solvent Accessible Surface", "Surface (A^2)",
-                  "Time", "Residue Index", nframe, nres,
-                  mat->axis_x, mat->axis_y, accr, lo, hi, rlo, rhi, &nlev);
+        nlev = static_cast<int>(hi - lo + 1);
+        write_xpm(fp,
+                  0,
+                  "Solvent Accessible Surface",
+                  "Surface (A^2)",
+                  "Time",
+                  "Residue Index",
+                  nframe,
+                  nres,
+                  mat->axis_x.data(),
+                  mat->axis_y.data(),
+                  accr,
+                  lo,
+                  hi,
+                  rlo,
+                  rhi,
+                  &nlev);
         gmx_ffclose(fp);
     }
 }
 
-static void analyse_ss(const char *outfile, t_matrix *mat, const char *ss_string,
-                       const gmx_output_env_t *oenv)
+static void analyse_ss(const char* outfile, t_matrix* mat, const char* ss_string, const gmx_output_env_t* oenv)
 {
-    FILE        *fp;
-    t_mapping   *map;
-    int          f, r, *count, *total, ss_count, total_count;
-    size_t       s;
-    const char** leg;
+    FILE* fp;
+    int   ss_count, total_count;
 
-    map = mat->map;
-    snew(count, mat->nmap);
-    snew(total, mat->nmap);
-    snew(leg, mat->nmap+1);
-    leg[0] = "Structure";
-    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+    gmx::ArrayRef<t_mapping> map = mat->map;
+    std::vector<int>         count(map.size());
+    std::vector<int>         total(map.size(), 0);
+    // This copying would not be necessary if xvgr_legend could take a
+    // view of string views
+    std::vector<std::string> leg;
+    leg.reserve(map.size() + 1);
+    leg.emplace_back("Structure");
+    for (const auto& m : map)
     {
-        leg[s+1] = gmx_strdup(map[s].desc);
+        leg.emplace_back(m.desc);
     }
 
-    fp = xvgropen(outfile, "Secondary Structure",
-                  output_env_get_xvgr_tlabel(oenv), "Number of Residues", oenv);
+    fp = xvgropen(
+            outfile, "Secondary Structure", output_env_get_xvgr_tlabel(oenv), "Number of Residues", oenv);
     if (output_env_get_print_xvgr_codes(oenv))
     {
         fprintf(fp, "@ subtitle \"Structure = ");
     }
-    for (s = 0; s < std::strlen(ss_string); s++)
+    for (size_t s = 0; s < std::strlen(ss_string); s++)
     {
         if (s > 0)
         {
             fprintf(fp, " + ");
         }
-        for (f = 0; f < mat->nmap; f++)
+        for (const auto& m : map)
         {
-            if (ss_string[s] == map[f].code.c1)
+            if (ss_string[s] == m.code.c1)
             {
-                fprintf(fp, "%s", map[f].desc);
+                fprintf(fp, "%s", m.desc);
             }
         }
     }
     fprintf(fp, "\"\n");
-    xvgr_legend(fp, mat->nmap+1, leg, oenv);
+    xvgrLegend(fp, leg, oenv);
 
     total_count = 0;
-    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
-    {
-        total[s] = 0;
-    }
-    for (f = 0; f < mat->nx; f++)
+    for (int f = 0; f < mat->nx; f++)
     {
         ss_count = 0;
-        for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+        for (auto& c : count)
         {
-            count[s] = 0;
+            c = 0;
         }
-        for (r = 0; r < mat->ny; r++)
+        for (int r = 0; r < mat->ny; r++)
         {
-            count[mat->matrix[f][r]]++;
-            total[mat->matrix[f][r]]++;
+            count[mat->matrix(f, r)]++;
+            total[mat->matrix(f, r)]++;
         }
-        for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+        for (gmx::index s = 0; s != gmx::ssize(map); ++s)
         {
             if (std::strchr(ss_string, map[s].code.c1))
             {
-                ss_count    += count[s];
+                ss_count += count[s];
                 total_count += count[s];
             }
         }
         fprintf(fp, "%8g %5d", mat->axis_x[f], ss_count);
-        for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+        for (const auto& c : count)
         {
-            fprintf(fp, " %5d", count[s]);
+            fprintf(fp, " %5d", c);
         }
         fprintf(fp, "\n");
     }
     /* now print column totals */
     fprintf(fp, "%-8s %5d", "# Totals", total_count);
-    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+    for (const auto& t : total)
     {
-        fprintf(fp, " %5d", total[s]);
+        fprintf(fp, " %5d", t);
     }
     fprintf(fp, "\n");
 
     /* now print probabilities */
     fprintf(fp, "%-8s %5.2f", "# SS pr.", total_count / static_cast<real>(mat->nx * mat->ny));
-    for (s = 0; s < static_cast<size_t>(mat->nmap); s++)
+    for (const auto& t : total)
     {
-        fprintf(fp, " %5.2f", total[s] / static_cast<real>(mat->nx * mat->ny));
+        fprintf(fp, " %5.2f", t / static_cast<real>(mat->nx * mat->ny));
     }
     fprintf(fp, "\n");
 
     xvgrclose(fp);
-    sfree(leg);
-    sfree(count);
 }
 
-int gmx_do_dssp(int argc, char *argv[])
+int gmx_do_dssp(int argc, charargv[])
 {
-    const char        *desc[] = {
+    const chardesc[] = {
         "[THISMODULE] ",
         "reads a trajectory file and computes the secondary structure for",
         "each time frame ",
         "calling the dssp program. If you do not have the dssp program,",
-        "get it from http://swift.cmbi.ru.nl/gv/dssp. [THISMODULE] assumes ",
+        "get it from https://swift.cmbi.umcn.nl/gv/dssp. [THISMODULE] assumes ",
         "that the dssp executable is located in ",
-        "[TT]/usr/local/bin/dssp[tt]. If this is not the case, then you should",
+        // NOLINTNEXTLINE(bugprone-suspicious-missing-comma)
+        "[TT]" GMX_DSSP_PROGRAM_PATH "[tt]. If this is not the case, then you should",
         "set an environment variable [TT]DSSP[tt] pointing to the dssp",
         "executable, e.g.: [PAR]",
         "[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
@@ -475,60 +514,65 @@ int gmx_do_dssp(int argc, char *argv[])
         "function of secondary structure type."
     };
     static gmx_bool    bVerbose;
-    static const char *ss_string   = "HEBT";
+    static const charss_string   = "HEBT";
     static int         dsspVersion = 2;
     t_pargs            pa[]        = {
-        { "-v",  FALSE, etBOOL, {&bVerbose},
-          "HIDDENGenerate miles of useless information" },
-        { "-sss", FALSE, etSTR, {&ss_string},
-          "Secondary structures for structure count"},
-        { "-ver", FALSE, etINT, {&dsspVersion},
-          "DSSP major version. Syntax changed with version 2"}
+        { "-v", FALSE, etBOOL, { &bVerbose }, "HIDDENGenerate miles of useless information" },
+        { "-sss", FALSE, etSTR, { &ss_string }, "Secondary structures for structure count" },
+        { "-ver",
+          FALSE,
+          etINT,
+          { &dsspVersion },
+          "DSSP major version. Syntax changed with version 2" }
     };
 
-    t_trxstatus       *status;
-    FILE              *tapein;
-    FILE              *ss, *acc, *fTArea, *tmpf;
-    const char        *fnSCount, *fnArea, *fnTArea, *fnAArea;
-    const char        *leg[] = { "Phobic", "Phylic" };
-    t_topology         top;
-    int                ePBC;
-    t_atoms           *atoms;
-    t_matrix           mat;
-    int                nres, nr0, naccr, nres_plus_separators;
-    gmx_bool          *bPhbres, bDoAccSurf;
-    real               t;
-    int                i, j, natoms, nframe = 0;
-    matrix             box = {{0}};
-    int                gnx;
-    char              *grpnm, *ss_str;
-    int               *index;
-    rvec              *xp, *x;
-    int               *average_area;
-    real             **accr, *accr_ptr = nullptr, *av_area, *norm_av_area;
-    char               pdbfile[32], tmpfile[32];
-    char               dssp[256];
-    const char        *dptr;
-    gmx_output_env_t  *oenv;
-    gmx_rmpbc_t        gpbc = nullptr;
-
-    t_filenm           fnm[] = {
-        { efTRX, "-f",   nullptr,      ffREAD },
-        { efTPS, nullptr,   nullptr,      ffREAD },
-        { efNDX, nullptr,   nullptr,      ffOPTRD },
-        { efDAT, "-ssdump", "ssdump", ffOPTWR },
-        { efMAP, "-map", "ss",      ffLIBRD },
-        { efXPM, "-o",   "ss",      ffWRITE },
-        { efXVG, "-sc",  "scount",  ffWRITE },
-        { efXPM, "-a",   "area",    ffOPTWR },
-        { efXVG, "-ta",  "totarea", ffOPTWR },
-        { efXVG, "-aa",  "averarea", ffOPTWR }
+    t_trxstatus*      status;
+    FILE *            tapein, *tapeout;
+    FILE *            ss, *acc, *fTArea, *tmpf;
+    const char *      fnSCount, *fnArea, *fnTArea, *fnAArea;
+    const char*       leg[] = { "Phobic", "Phylic" };
+    t_topology        top;
+    PbcType           pbcType;
+    t_atoms*          atoms;
+    t_matrix          mat;
+    int               nres, nr0, naccr, nres_plus_separators;
+    gmx_bool *        bPhbres, bDoAccSurf;
+    real              t;
+    int               natoms, nframe = 0;
+    matrix            box = { { 0 } };
+    int               gnx;
+    char*             grpnm;
+    int*              index;
+    rvec *            xp, *x;
+    int*              average_area;
+    real **           accr, *accr_ptr = nullptr, *av_area, *norm_av_area;
+    char              pdbfile[32], tmpfile[32];
+    char              dssp[256];
+    const char*       dptr;
+    gmx_output_env_t* oenv;
+    gmx_rmpbc_t       gpbc = nullptr;
+
+    t_filenm fnm[] = {
+        { efTRX, "-f", nullptr, ffREAD },     { efTPS, nullptr, nullptr, ffREAD },
+        { efNDX, nullptr, nullptr, ffOPTRD }, { efDAT, "-ssdump", "ssdump", ffOPTWR },
+        { efMAP, "-map", "ss", ffLIBRD },     { efXPM, "-o", "ss", ffWRITE },
+        { efXVG, "-sc", "scount", ffWRITE },  { efXPM, "-a", "area", ffOPTWR },
+        { efXVG, "-ta", "totarea", ffOPTWR }, { efXVG, "-aa", "averarea", ffOPTWR }
     };
 #define NFILE asize(fnm)
 
-    if (!parse_common_args(&argc, argv,
+    if (!parse_common_args(&argc,
+                           argv,
                            PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT,
-                           NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
+                           NFILE,
+                           fnm,
+                           asize(pa),
+                           pa,
+                           asize(desc),
+                           desc,
+                           0,
+                           nullptr,
+                           &oenv))
     {
         return 0;
     }
@@ -538,7 +582,7 @@ int gmx_do_dssp(int argc, char *argv[])
     fnAArea    = opt2fn_null("-aa", NFILE, fnm);
     bDoAccSurf = ((fnArea != nullptr) || (fnTArea != nullptr) || (fnAArea != nullptr));
 
-    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xp, nullptr, box, FALSE);
+    read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xp, nullptr, box, FALSE);
     atoms = &(top.atoms);
     check_oo(atoms);
     bPhbres = bPhobics(atoms);
@@ -546,7 +590,7 @@ int gmx_do_dssp(int argc, char *argv[])
     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpnm);
     nres = 0;
     nr0  = -1;
-    for (i = 0; (i < gnx); i++)
+    for (int i = 0; (i < gnx); i++)
     {
         if (atoms->atom[index[i]].resind != nr0)
         {
@@ -567,10 +611,7 @@ int gmx_do_dssp(int argc, char *argv[])
             gmx_fatal(FARGS, "Can not open tmp file %s", pdbfile);
         }
     }
-    else
-    {
-        fclose(tmpf);
-    }
+    fclose(tmpf);
 
     std::strcpy(tmpfile, "ddXXXXXX");
     gmx_tmpnam(tmpfile);
@@ -583,42 +624,55 @@ int gmx_do_dssp(int argc, char *argv[])
             gmx_fatal(FARGS, "Can not open tmp file %s", tmpfile);
         }
     }
-    else
-    {
-        fclose(tmpf);
-    }
+    fclose(tmpf);
 
+    const std::string defpathenv = GMX_DSSP_PROGRAM_PATH;
     if ((dptr = getenv("DSSP")) == nullptr)
     {
-        dptr = "/usr/local/bin/dssp";
+        dptr = defpathenv.c_str();
     }
     if (!gmx_fexist(dptr))
     {
-        gmx_fatal(FARGS, "DSSP executable (%s) does not exist (use setenv DSSP)",
-                  dptr);
+        gmx_fatal(FARGS, "DSSP executable (%s) does not exist (use setenv DSSP)", dptr);
     }
+    std::string redirectionString;
+    redirectionString = prepareToRedirectStdout(bVerbose);
+    DsspInputStrings dsspStrings;
+    dsspStrings.dptr    = dptr;
+    dsspStrings.pdbfile = pdbfile;
+    dsspStrings.tmpfile = tmpfile;
     if (dsspVersion >= 2)
     {
         if (dsspVersion > 2)
         {
-            printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by do_dssp. Assuming version 2 syntax.\n\n", dsspVersion);
+            printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by "
+                   "do_dssp. Assuming version 2 syntax.\n\n",
+                   dsspVersion);
         }
 
-        sprintf(dssp, "%s -i %s -o %s > /dev/null %s",
-                dptr, pdbfile, tmpfile, bVerbose ? "" : "2> /dev/null");
+        printDsspResult(dssp, dsspStrings, redirectionString);
     }
     else
     {
-        sprintf(dssp, "%s %s %s %s > /dev/null %s",
-                dptr, bDoAccSurf ? "" : "-na", pdbfile, tmpfile, bVerbose ? "" : "2> /dev/null");
-
+        if (bDoAccSurf)
+        {
+            dsspStrings.dptr.clear();
+        }
+        else
+        {
+            dsspStrings.dptr = "-na";
+        }
+        printDsspResult(dssp, dsspStrings, redirectionString);
     }
     fprintf(stderr, "dssp cmd='%s'\n", dssp);
 
     if (fnTArea)
     {
-        fTArea = xvgropen(fnTArea, "Solvent Accessible Surface Area",
-                          output_env_get_xvgr_tlabel(oenv), "Area (nm\\S2\\N)", oenv);
+        fTArea = xvgropen(fnTArea,
+                          "Solvent Accessible Surface Area",
+                          output_env_get_xvgr_tlabel(oenv),
+                          "Area (nm\\S2\\N)",
+                          oenv);
         xvgr_legend(fTArea, 2, leg, oenv);
     }
     else
@@ -626,8 +680,7 @@ int gmx_do_dssp(int argc, char *argv[])
         fTArea = nullptr;
     }
 
-    mat.map  = nullptr;
-    mat.nmap = readcmap(opt2fn("-map", NFILE, fnm), &(mat.map));
+    mat.map = readcmap(opt2fn("-map", NFILE, fnm));
 
     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
     if (natoms > atoms->nr)
@@ -645,7 +698,7 @@ int gmx_do_dssp(int argc, char *argv[])
     accr  = nullptr;
     naccr = 0;
 
-    gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
+    gpbc = gmx_rmpbc_init(&top.idef, pbcType, natoms);
     do
     {
         t = output_env_conv_time(oenv, t);
@@ -653,37 +706,48 @@ int gmx_do_dssp(int argc, char *argv[])
         {
             naccr += 10;
             srenew(accr, naccr);
-            for (i = naccr-10; i < naccr; i++)
+            for (int i = naccr - 10; i < naccr; i++)
             {
-                snew(accr[i], 2*atoms->nres-1);
+                snew(accr[i], 2 * atoms->nres - 1);
             }
         }
         gmx_rmpbc(gpbc, natoms, box, x);
         tapein = gmx_ffopen(pdbfile, "w");
-        write_pdbfile_indexed(tapein, nullptr, atoms, x, ePBC, box, ' ', -1, gnx, index, nullptr, TRUE, FALSE);
+        write_pdbfile_indexed(tapein, nullptr, atoms, x, pbcType, box, ' ', -1, gnx, index, nullptr, FALSE);
         gmx_ffclose(tapein);
-
-        if (0 != system(dssp))
-        {
-            gmx_fatal(FARGS, "Failed to execute command: %s\n"
-                      "Try specifying your dssp version with the -ver option.", dssp);
-        }
-
         /* strip_dssp returns the number of lines found in the dssp file, i.e.
          * the number of residues plus the separator lines */
 
+#if HAVE_PIPES || GMX_NATIVE_WINDOWS
+        if (nullptr == (tapeout = popen(dssp, "r")))
+#else
+        if (0 != system(dssp) || nullptr == (tapeout = gmx_ffopen(tmpfile, "r")))
+#endif
+        {
+            remove(pdbfile);
+            remove(tmpfile);
+            gmx_fatal(FARGS,
+                      "Failed to execute command: %s\n"
+                      "Try specifying your dssp version with the -ver option.",
+                      dssp);
+        }
         if (bDoAccSurf)
         {
             accr_ptr = accr[nframe];
         }
-
-        nres_plus_separators = strip_dssp(tmpfile, nres, bPhbres, t,
-                                          accr_ptr, fTArea, &mat, average_area, oenv);
+        /* strip_dssp returns the number of lines found in the dssp file, i.e.
+         * the number of residues plus the separator lines */
+        nres_plus_separators =
+                strip_dssp(tapeout, nres, bPhbres, t, accr_ptr, fTArea, &mat, average_area, oenv);
+#if HAVE_PIPES || GMX_NATIVE_WINDOWS
+        pclose(tapeout);
+#else
+        gmx_ffclose(tapeout);
+#endif
         remove(tmpfile);
         remove(pdbfile);
         nframe++;
-    }
-    while (read_next_x(oenv, status, &t, x, box));
+    } while (read_next_x(oenv, status, &t, x, box));
     fprintf(stderr, "\n");
     close_trx(status);
     if (fTArea)
@@ -702,19 +766,17 @@ int gmx_do_dssp(int argc, char *argv[])
     if (opt2bSet("-ssdump", NFILE, fnm))
     {
         ss = opt2FILE("-ssdump", NFILE, fnm, "w");
-        snew(ss_str, nres+1);
         fprintf(ss, "%d\n", nres);
-        for (j = 0; j < mat.nx; j++)
+        for (gmx::index j = 0; j != mat.matrix.extent(0); ++j)
         {
-            for (i = 0; (i < mat.ny); i++)
+            auto row = mat.matrix.asView()[j];
+            for (gmx::index i = 0; i != row.extent(0); ++i)
             {
-                ss_str[i] = mat.map[mat.matrix[j][i]].code.c1;
+                fputc(mat.map[row[i]].code.c1, ss);
             }
-            ss_str[i] = '\0';
-            fprintf(ss, "%s\n", ss_str);
+            fputc('\n', ss);
         }
         gmx_ffclose(ss);
-        sfree(ss_str);
     }
     analyse_ss(fnSCount, &mat, ss_string, oenv);
 
@@ -722,7 +784,7 @@ int gmx_do_dssp(int argc, char *argv[])
     {
         write_sas_mat(fnArea, accr, nframe, nres_plus_separators, &mat);
 
-        for (i = 0; i < atoms->nres; i++)
+        for (int i = 0; i < atoms->nres; i++)
         {
             av_area[i] = (average_area[i] / static_cast<real>(nframe));
         }
@@ -731,11 +793,10 @@ int gmx_do_dssp(int argc, char *argv[])
 
         if (fnAArea)
         {
-            acc = xvgropen(fnAArea, "Average Accessible Area",
-                           "Residue", "A\\S2", oenv);
-            for (i = 0; (i < nres); i++)
+            acc = xvgropen(fnAArea, "Average Accessible Area", "Residue", "A\\S2", oenv);
+            for (int i = 0; (i < nres); i++)
             {
-                fprintf(acc, "%5d  %10g %10g\n", i+1, av_area[i], norm_av_area[i]);
+                fprintf(acc, "%5d  %10g %10g\n", i + 1, av_area[i], norm_av_area[i]);
             }
             xvgrclose(acc);
         }