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