src/gromacs/gmxana/gmx_do_dssp.cpp

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