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. [TT]do_dssp[tt] assumes that the dssp executable is",
 392     "/usr/local/bin/dssp. If this is not the case, then you should",
 393     "set an environment variable [TT]DSSP[tt] pointing to the dssp",
 394     "executable, e.g.: [PAR]",
 395     "[TT]setenv DSSP /opt/dssp/bin/dssp[tt][PAR]",
 396     "The structure assignment for each residue and time is written to an",
 397     "[TT].xpm[tt] matrix file. This file can be visualized with for instance",
 398     "[TT]xv[tt] and can be converted to postscript with [TT]xpm2ps[tt].",
 399     "Individual chains are separated by light grey lines in the [TT].xpm[tt] and",
 400     "postscript files.",
 401     "The number of residues with each secondary structure type and the",
 402     "total secondary structure ([TT]-sss[tt]) count as a function of",
 403     "time are also written to file ([TT]-sc[tt]).[PAR]",
 404     "Solvent accessible surface (SAS) per residue can be calculated, both in",
 405     "absolute values (A^2) and in fractions of the maximal accessible",
 406     "surface of a residue. The maximal accessible surface is defined as",
 407     "the accessible surface of a residue in a chain of glycines.",
 408     "[BB]Note[bb] that the program [TT]g_sas[tt] can also compute SAS",
 409     "and that is more efficient.[PAR]",
 410     "Finally, this program can dump the secondary structure in a special file",
 411     "[TT]ssdump.dat[tt] for usage in the program [TT]g_chi[tt]. Together",
 412     "these two programs can be used to analyze dihedral properties as a",
 413     "function of secondary structure type."
 414   };
 415   static gmx_bool bVerbose;
 416   static const char *ss_string="HEBT";
 417   t_pargs pa[] = {
 418     { "-v",  FALSE, etBOOL, {&bVerbose},
 419       "HIDDENGenerate miles of useless information" },
 420     { "-sss", FALSE, etSTR, {&ss_string},
 421       "Secondary structures for structure count"}
 422   };
 423
 424   t_trxstatus *status;
 425   FILE       *tapein;
 426   FILE       *ss,*acc,*fTArea,*tmpf;
 427   const char *fnSCount,*fnArea,*fnTArea,*fnAArea;
 428   const char *leg[] = { "Phobic", "Phylic" };
 429   t_topology top;
 430   int        ePBC;
 431   t_atoms    *atoms;
 432   t_matrix   mat;
 433   int        nres,nr0,naccr,nres_plus_separators;
 434   gmx_bool       *bPhbres,bDoAccSurf;
 435   real       t;
 436   int        i,j,natoms,nframe=0;
 437   matrix     box;
 438   int        gnx;
 439   char       *grpnm,*ss_str;
 440   atom_id    *index;
 441   rvec       *xp,*x;
 442   int        *average_area;
 443   real       **accr,*accr_ptr=NULL,*av_area, *norm_av_area;
 444   char       pdbfile[32],tmpfile[32],title[256];
 445   char       dssp[256];
 446   const char *dptr;
 447   output_env_t oenv;
 448   gmx_rmpbc_t  gpbc=NULL;
 449
 450   t_filenm   fnm[] = {
 451     { efTRX, "-f",   NULL,      ffREAD },
 452     { efTPS, NULL,   NULL,      ffREAD },
 453     { efNDX, NULL,   NULL,      ffOPTRD },
 454     { efDAT, "-ssdump", "ssdump", ffOPTWR },
 455     { efMAP, "-map", "ss",      ffLIBRD },
 456     { efXPM, "-o",   "ss",      ffWRITE },
 457     { efXVG, "-sc",  "scount",  ffWRITE },
 458     { efXPM, "-a",   "area",    ffOPTWR },
 459     { efXVG, "-ta",  "totarea", ffOPTWR },
 460     { efXVG, "-aa",  "averarea",ffOPTWR }
 461   };
 462 #define NFILE asize(fnm)
 463
 464   CopyRight(stderr,argv[0]);
 465   parse_common_args(&argc,argv,
 466                     PCA_CAN_TIME | PCA_CAN_VIEW | PCA_TIME_UNIT | PCA_BE_NICE ,
 467                     NFILE,fnm, asize(pa),pa, asize(desc),desc,0,NULL,&oenv);
 468   fnSCount= opt2fn("-sc",NFILE,fnm);
 469   fnArea  = opt2fn_null("-a", NFILE,fnm);
 470   fnTArea = opt2fn_null("-ta",NFILE,fnm);
 471   fnAArea = opt2fn_null("-aa",NFILE,fnm);
 472   bDoAccSurf=(fnArea || fnTArea || fnAArea);
 473
 474   read_tps_conf(ftp2fn(efTPS,NFILE,fnm),title,&top,&ePBC,&xp,NULL,box,FALSE);
 475   atoms=&(top.atoms);
 476   check_oo(atoms);
 477   bPhbres=bPhobics(atoms);
 478
 479   get_index(atoms,ftp2fn_null(efNDX,NFILE,fnm),1,&gnx,&index,&grpnm);
 480   nres=0;
 481   nr0=-1;
 482   for(i=0; (i<gnx); i++) {
 483     if (atoms->atom[index[i]].resind != nr0) {
 484       nr0=atoms->atom[index[i]].resind;
 485       nres++;
 486     }
 487   }
 488   fprintf(stderr,"There are %d residues in your selected group\n",nres);
 489
 490   strcpy(pdbfile,"ddXXXXXX");
 491   gmx_tmpnam(pdbfile);
 492   if ((tmpf = fopen(pdbfile,"w")) == NULL) {
 493     sprintf(pdbfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);
 494     gmx_tmpnam(pdbfile);
 495     if ((tmpf = fopen(pdbfile,"w")) == NULL)
 496       gmx_fatal(FARGS,"Can not open tmp file %s",pdbfile);
 497   }
 498   else
 499     fclose(tmpf);
 500
 501   strcpy(tmpfile,"ddXXXXXX");
 502   gmx_tmpnam(tmpfile);
 503   if ((tmpf = fopen(tmpfile,"w")) == NULL) {
 504     sprintf(tmpfile,"%ctmp%cfilterXXXXXX",DIR_SEPARATOR,DIR_SEPARATOR);
 505     gmx_tmpnam(tmpfile);
 506     if ((tmpf = fopen(tmpfile,"w")) == NULL)
 507       gmx_fatal(FARGS,"Can not open tmp file %s",tmpfile);
 508   }
 509   else
 510     fclose(tmpf);
 511
 512   if ((dptr=getenv("DSSP")) == NULL)
 513     dptr="/usr/local/bin/dssp";
 514   if (!gmx_fexist(dptr))
 515     gmx_fatal(FARGS,"DSSP executable (%s) does not exist (use setenv DSSP)",
 516                 dptr);
 517   sprintf(dssp,"%s %s %s %s > /dev/null %s",
 518           dptr,bDoAccSurf?"":"-na",pdbfile,tmpfile,bVerbose?"":"2> /dev/null");
 519   if (bVerbose)
 520     fprintf(stderr,"dssp cmd='%s'\n",dssp);
 521
 522   if (fnTArea) {
 523     fTArea=xvgropen(fnTArea,"Solvent Accessible Surface Area",
 524                     output_env_get_xvgr_tlabel(oenv),"Area (nm\\S2\\N)",oenv);
 525     xvgr_legend(fTArea,2,leg,oenv);
 526   } else
 527     fTArea=NULL;
 528
 529   mat.map=NULL;
 530   mat.nmap=getcmap(libopen(opt2fn("-map",NFILE,fnm)),
 531                    opt2fn("-map",NFILE,fnm),&(mat.map));
 532
 533   natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
 534   if (natoms > atoms->nr)
 535     gmx_fatal(FARGS,"\nTrajectory does not match topology!");
 536   if (gnx > natoms)
 537     gmx_fatal(FARGS,"\nTrajectory does not match selected group!");
 538
 539   snew(average_area, atoms->nres);
 540   snew(av_area     , atoms->nres);
 541   snew(norm_av_area, atoms->nres);
 542   accr=NULL;
 543   naccr=0;
 544
 545   gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
 546   do {
 547     t = output_env_conv_time(oenv,t);
 548     if (bDoAccSurf && nframe>=naccr) {
 549       naccr+=10;
 550       srenew(accr,naccr);
 551       for(i=naccr-10; i<naccr; i++)
 552         snew(accr[i],2*atoms->nres-1);
 553     }
 554     gmx_rmpbc(gpbc,natoms,box,x);
 555     tapein=ffopen(pdbfile,"w");
 556     write_pdbfile_indexed(tapein,NULL,atoms,x,ePBC,box,' ',-1,gnx,index,NULL,TRUE);
 557     ffclose(tapein);
 558
 559 #ifdef GMX_NO_SYSTEM
 560     printf("Warning-- No calls to system(3) supported on this platform.");
 561     printf("Warning-- Skipping execution of 'system(\"%s\")'.", dssp);
 562     exit(1);
 563 #else
 564     if(0 != system(dssp))
 565     {
 566         gmx_fatal(FARGS,"Failed to execute command: %s",dssp);
 567     }
 568 #endif
 569
 570     /* strip_dssp returns the number of lines found in the dssp file, i.e.
 571      * the number of residues plus the separator lines */
 572
 573     if (bDoAccSurf)
 574         accr_ptr = accr[nframe];
 575
 576     nres_plus_separators = strip_dssp(tmpfile,nres,bPhbres,t,
 577                                       accr_ptr,fTArea,&mat,average_area,oenv);
 578     remove(tmpfile);
 579     remove(pdbfile);
 580     nframe++;
 581   } while(read_next_x(oenv,status,&t,natoms,x,box));
 582   fprintf(stderr,"\n");
 583   close_trj(status);
 584   if (fTArea)
 585     ffclose(fTArea);
 586   gmx_rmpbc_done(gpbc);
 587
 588   prune_ss_legend(&mat);
 589
 590   ss=opt2FILE("-o",NFILE,fnm,"w");
 591   mat.flags = 0;
 592   write_xpm_m(ss,mat);
 593   ffclose(ss);
 594
 595   if (opt2bSet("-ssdump",NFILE,fnm))
 596   {
 597       ss = opt2FILE("-ssdump",NFILE,fnm,"w");
 598       snew(ss_str,nres+1);
 599       fprintf(ss,"%d\n",nres);
 600       for(j=0; j<mat.nx; j++)
 601       {
 602           for(i=0; (i<mat.ny); i++)
 603           {
 604               ss_str[i] = mat.map[mat.matrix[j][i]].code.c1;
 605           }
 606           ss_str[i] = '\0';
 607           fprintf(ss,"%s\n",ss_str);
 608       }
 609       ffclose(ss);
 610       sfree(ss_str);
 611   }
 612   analyse_ss(fnSCount,&mat,ss_string,oenv);
 613
 614   if (bDoAccSurf) {
 615     write_sas_mat(fnArea,accr,nframe,nres_plus_separators,&mat);
 616
 617     for(i=0; i<atoms->nres; i++)
 618       av_area[i] = (average_area[i] / (real)nframe);
 619
 620     norm_acc(atoms, nres, av_area, norm_av_area);
 621
 622     if (fnAArea) {
 623       acc=xvgropen(fnAArea,"Average Accessible Area",
 624                    "Residue","A\\S2",oenv);
 625       for(i=0; (i<nres); i++)
 626         fprintf(acc,"%5d  %10g %10g\n",i+1,av_area[i], norm_av_area[i]);
 627       ffclose(acc);
 628     }
 629   }
 630
 631   view_all(oenv, NFILE, fnm);
 632
 633   thanx(stderr);
 634
 635   return 0;
 636 }