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