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