src/tools/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 #include "gmx_ana.h"
  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     CopyRight(stderr, argv[0]);
 525     parse_common_args(&argc, argv,
 526                       PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT | PCA_BE_NICE,
 527                       NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
 528     fnSCount   = opt2fn("-sc", NFILE, fnm);
 529     fnArea     = opt2fn_null("-a", NFILE, fnm);
 530     fnTArea    = opt2fn_null("-ta", NFILE, fnm);
 531     fnAArea    = opt2fn_null("-aa", NFILE, fnm);
 532     bDoAccSurf = (fnArea || fnTArea || fnAArea);
 533
 534     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xp, NULL, box, FALSE);
 535     atoms = &(top.atoms);
 536     check_oo(atoms);
 537     bPhbres = bPhobics(atoms);
 538
 539     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &gnx, &index, &grpnm);
 540     nres = 0;
 541     nr0  = -1;
 542     for (i = 0; (i < gnx); i++)
 543     {
 544         if (atoms->atom[index[i]].resind != nr0)
 545         {
 546             nr0 = atoms->atom[index[i]].resind;
 547             nres++;
 548         }
 549     }
 550     fprintf(stderr, "There are %d residues in your selected group\n", nres);
 551
 552     strcpy(pdbfile, "ddXXXXXX");
 553     gmx_tmpnam(pdbfile);
 554     if ((tmpf = fopen(pdbfile, "w")) == NULL)
 555     {
 556         sprintf(pdbfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
 557         gmx_tmpnam(pdbfile);
 558         if ((tmpf = fopen(pdbfile, "w")) == NULL)
 559         {
 560             gmx_fatal(FARGS, "Can not open tmp file %s", pdbfile);
 561         }
 562     }
 563     else
 564     {
 565         fclose(tmpf);
 566     }
 567
 568     strcpy(tmpfile, "ddXXXXXX");
 569     gmx_tmpnam(tmpfile);
 570     if ((tmpf = fopen(tmpfile, "w")) == NULL)
 571     {
 572         sprintf(tmpfile, "%ctmp%cfilterXXXXXX", DIR_SEPARATOR, DIR_SEPARATOR);
 573         gmx_tmpnam(tmpfile);
 574         if ((tmpf = fopen(tmpfile, "w")) == NULL)
 575         {
 576             gmx_fatal(FARGS, "Can not open tmp file %s", tmpfile);
 577         }
 578     }
 579     else
 580     {
 581         fclose(tmpf);
 582     }
 583
 584     if ((dptr = getenv("DSSP")) == NULL)
 585     {
 586         dptr = "/usr/local/bin/dssp";
 587     }
 588     if (!gmx_fexist(dptr))
 589     {
 590         gmx_fatal(FARGS, "DSSP executable (%s) does not exist (use setenv DSSP)",
 591                   dptr);
 592     }
 593     if (dsspVersion >= 2)
 594     {
 595         if (dsspVersion > 2)
 596         {
 597             printf("\nWARNING: You use DSSP version %d, which is not explicitly\nsupported by do_dssp. Assuming version 2 syntax.\n\n", dsspVersion);
 598         }
 599
 600         sprintf(dssp, "%s -i %s -o %s > /dev/null %s",
 601                 dptr, pdbfile, tmpfile, bVerbose ? "" : "2> /dev/null");
 602     }
 603     else
 604     {
 605         sprintf(dssp, "%s %s %s %s > /dev/null %s",
 606                 dptr, bDoAccSurf ? "" : "-na", pdbfile, tmpfile, bVerbose ? "" : "2> /dev/null");
 607
 608     }
 609     fprintf(stderr, "dssp cmd='%s'\n", dssp);
 610
 611     if (fnTArea)
 612     {
 613         fTArea = xvgropen(fnTArea, "Solvent Accessible Surface Area",
 614                           output_env_get_xvgr_tlabel(oenv), "Area (nm\\S2\\N)", oenv);
 615         xvgr_legend(fTArea, 2, leg, oenv);
 616     }
 617     else
 618     {
 619         fTArea = NULL;
 620     }
 621
 622     mat.map  = NULL;
 623     mat.nmap = getcmap(libopen(opt2fn("-map", NFILE, fnm)),
 624                        opt2fn("-map", NFILE, fnm), &(mat.map));
 625
 626     natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
 627     if (natoms > atoms->nr)
 628     {
 629         gmx_fatal(FARGS, "\nTrajectory does not match topology!");
 630     }
 631     if (gnx > natoms)
 632     {
 633         gmx_fatal(FARGS, "\nTrajectory does not match selected group!");
 634     }
 635
 636     snew(average_area, atoms->nres);
 637     snew(av_area, atoms->nres);
 638     snew(norm_av_area, atoms->nres);
 639     accr  = NULL;
 640     naccr = 0;
 641
 642     gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms, box);
 643     do
 644     {
 645         t = output_env_conv_time(oenv, t);
 646         if (bDoAccSurf && nframe >= naccr)
 647         {
 648             naccr += 10;
 649             srenew(accr, naccr);
 650             for (i = naccr-10; i < naccr; i++)
 651             {
 652                 snew(accr[i], 2*atoms->nres-1);
 653             }
 654         }
 655         gmx_rmpbc(gpbc, natoms, box, x);
 656         tapein = ffopen(pdbfile, "w");
 657         write_pdbfile_indexed(tapein, NULL, atoms, x, ePBC, box, ' ', -1, gnx, index, NULL, TRUE);
 658         ffclose(tapein);
 659
 660 #ifdef GMX_NO_SYSTEM
 661         printf("Warning-- No calls to system(3) supported on this platform.");
 662         printf("Warning-- Skipping execution of 'system(\"%s\")'.", dssp);
 663         exit(1);
 664 #else
 665         if (0 != system(dssp))
 666         {
 667             gmx_fatal(FARGS, "Failed to execute command: %s\n",
 668                       "Try specifying your dssp version with the -ver option.", dssp);
 669         }
 670 #endif
 671
 672         /* strip_dssp returns the number of lines found in the dssp file, i.e.
 673          * the number of residues plus the separator lines */
 674
 675         if (bDoAccSurf)
 676         {
 677             accr_ptr = accr[nframe];
 678         }
 679
 680         nres_plus_separators = strip_dssp(tmpfile, nres, bPhbres, t,
 681                                           accr_ptr, fTArea, &mat, average_area, oenv);
 682         remove(tmpfile);
 683         remove(pdbfile);
 684         nframe++;
 685     }
 686     while (read_next_x(oenv, status, &t, natoms, x, box));
 687     fprintf(stderr, "\n");
 688     close_trj(status);
 689     if (fTArea)
 690     {
 691         ffclose(fTArea);
 692     }
 693     gmx_rmpbc_done(gpbc);
 694
 695     prune_ss_legend(&mat);
 696
 697     ss        = opt2FILE("-o", NFILE, fnm, "w");
 698     mat.flags = 0;
 699     write_xpm_m(ss, mat);
 700     ffclose(ss);
 701
 702     if (opt2bSet("-ssdump", NFILE, fnm))
 703     {
 704         ss = opt2FILE("-ssdump", NFILE, fnm, "w");
 705         snew(ss_str, nres+1);
 706         fprintf(ss, "%d\n", nres);
 707         for (j = 0; j < mat.nx; j++)
 708         {
 709             for (i = 0; (i < mat.ny); i++)
 710             {
 711                 ss_str[i] = mat.map[mat.matrix[j][i]].code.c1;
 712             }
 713             ss_str[i] = '\0';
 714             fprintf(ss, "%s\n", ss_str);
 715         }
 716         ffclose(ss);
 717         sfree(ss_str);
 718     }
 719     analyse_ss(fnSCount, &mat, ss_string, oenv);
 720
 721     if (bDoAccSurf)
 722     {
 723         write_sas_mat(fnArea, accr, nframe, nres_plus_separators, &mat);
 724
 725         for (i = 0; i < atoms->nres; i++)
 726         {
 727             av_area[i] = (average_area[i] / (real)nframe);
 728         }
 729
 730         norm_acc(atoms, nres, av_area, norm_av_area);
 731
 732         if (fnAArea)
 733         {
 734             acc = xvgropen(fnAArea, "Average Accessible Area",
 735                            "Residue", "A\\S2", oenv);
 736             for (i = 0; (i < nres); i++)
 737             {
 738                 fprintf(acc, "%5d  %10g %10g\n", i+1, av_area[i], norm_av_area[i]);
 739             }
 740             ffclose(acc);
 741         }
 742     }
 743
 744     view_all(oenv, NFILE, fnm);
 745
 746     thanx(stderr);
 747
 748     return 0;
 749 }