src/gromacs/gmxana/gmx_do_dssp.c

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