src/gromacs/gmxana/gmx_do_dssp.c

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