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