src/tools/gmx_sas.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.3.2
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2007, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Groningen Machine for Chemical Simulation
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38 #include <math.h>
  39 #include <stdlib.h>
  40
  41 #include "sysstuff.h"
  42 #include "string.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "vec.h"
  47 #include "xvgr.h"
  48 #include "pbc.h"
  49 #include "copyrite.h"
  50 #include "futil.h"
  51 #include "statutil.h"
  52 #include "index.h"
  53 #include "nsc.h"
  54 #include "pdbio.h"
  55 #include "confio.h"
  56 #include "rmpbc.h"
  57 #include "names.h"
  58 #include "atomprop.h"
  59 #include "physics.h"
  60 #include "tpxio.h"
  61 #include "gmx_ana.h"
  62
  63
  64 typedef struct {
  65   atom_id  aa,ab;
  66   real     d2a,d2b;
  67 } t_conect;
  68
  69 void add_rec(t_conect c[],atom_id i,atom_id j,real d2)
  70 {
  71   if (c[i].aa == NO_ATID) {
  72     c[i].aa  = j;
  73     c[i].d2a = d2;
  74   }
  75   else if (c[i].ab == NO_ATID) {
  76     c[i].ab  = j;
  77     c[i].d2b = d2;
  78   }
  79   else if (d2 < c[i].d2a) {
  80     c[i].aa  = j;
  81     c[i].d2a = d2;
  82   }
  83   else if (d2 < c[i].d2b) {
  84     c[i].ab  = j;
  85     c[i].d2b = d2;
  86   }
  87   /* Swap them if necessary: a must be larger than b */
  88   if (c[i].d2a < c[i].d2b) {
  89     j        = c[i].ab;
  90     c[i].ab  = c[i].aa;
  91     c[i].aa  = j;
  92     d2       = c[i].d2b;
  93     c[i].d2b = c[i].d2a;
  94     c[i].d2a = d2;
  95   }
  96 }
  97
  98 void do_conect(const char *fn,int n,rvec x[])
  99 {
 100   FILE     *fp;
 101   int      i,j;
 102   t_conect *c;
 103   rvec     dx;
 104   real     d2;
 105
 106   fprintf(stderr,"Building CONECT records\n");
 107   snew(c,n);
 108   for(i=0; (i<n); i++)
 109     c[i].aa = c[i].ab = NO_ATID;
 110
 111   for(i=0; (i<n); i++) {
 112     for(j=i+1; (j<n); j++) {
 113       rvec_sub(x[i],x[j],dx);
 114       d2 = iprod(dx,dx);
 115       add_rec(c,i,j,d2);
 116       add_rec(c,j,i,d2);
 117     }
 118   }
 119   fp = ffopen(fn,"a");
 120   for(i=0; (i<n); i++) {
 121     if ((c[i].aa == NO_ATID) || (c[i].ab == NO_ATID))
 122       fprintf(stderr,"Warning dot %d has no conections\n",i+1);
 123     fprintf(fp,"CONECT%5d%5d%5d\n",i+1,c[i].aa+1,c[i].ab+1);
 124   }
 125   ffclose(fp);
 126   sfree(c);
 127 }
 128
 129 void connelly_plot(const char *fn,int ndots,real dots[],rvec x[],t_atoms *atoms,
 130                    t_symtab *symtab,int ePBC,matrix box,bool bSave)
 131 {
 132   static const char *atomnm="DOT";
 133   static const char *resnm ="DOT";
 134   static const char *title ="Connely Dot Surface Generated by g_sas";
 135
 136   int  i,i0,r0,ii0,k;
 137   rvec *xnew;
 138   t_atoms aaa;
 139
 140   if (bSave) {
 141     i0 = atoms->nr;
 142     r0 = atoms->nres;
 143     srenew(atoms->atom,atoms->nr+ndots);
 144     srenew(atoms->atomname,atoms->nr+ndots);
 145     srenew(atoms->resinfo,r0+1);
 146     atoms->atom[i0].resind = r0;
 147     t_atoms_set_resinfo(atoms,i0,symtab,resnm,r0+1,' ',' ');
 148     srenew(atoms->pdbinfo,atoms->nr+ndots);
 149     snew(xnew,atoms->nr+ndots);
 150     for(i=0; (i<atoms->nr); i++)
 151       copy_rvec(x[i],xnew[i]);
 152     for(i=k=0; (i<ndots); i++) {
 153       ii0 = i0+i;
 154       atoms->atomname[ii0] = put_symtab(symtab,atomnm);
 155       atoms->pdbinfo[ii0].type = epdbATOM;
 156       atoms->pdbinfo[ii0].atomnr= ii0;
 157       atoms->atom[ii0].resind = r0;
 158       xnew[ii0][XX] = dots[k++];
 159       xnew[ii0][YY] = dots[k++];
 160       xnew[ii0][ZZ] = dots[k++];
 161       atoms->pdbinfo[ii0].bfac  = 0.0;
 162       atoms->pdbinfo[ii0].occup = 0.0;
 163     }
 164     atoms->nr   = i0+ndots;
 165     atoms->nres = r0+1;
 166     write_sto_conf(fn,title,atoms,xnew,NULL,ePBC,box);
 167     atoms->nres = r0;
 168     atoms->nr   = i0;
 169   }
 170   else {
 171     init_t_atoms(&aaa,ndots,TRUE);
 172     aaa.atom[0].resind = 0;
 173     t_atoms_set_resinfo(&aaa,0,symtab,resnm,1,' ',' ');
 174     snew(xnew,ndots);
 175     for(i=k=0; (i<ndots); i++) {
 176       ii0 = i;
 177       aaa.atomname[ii0] = put_symtab(symtab,atomnm);
 178       aaa.pdbinfo[ii0].type = epdbATOM;
 179       aaa.pdbinfo[ii0].atomnr= ii0;
 180       aaa.atom[ii0].resind = 0;
 181       xnew[ii0][XX] = dots[k++];
 182       xnew[ii0][YY] = dots[k++];
 183       xnew[ii0][ZZ] = dots[k++];
 184       aaa.pdbinfo[ii0].bfac  = 0.0;
 185       aaa.pdbinfo[ii0].occup = 0.0;
 186     }
 187     aaa.nr = ndots;
 188     write_sto_conf(fn,title,&aaa,xnew,NULL,ePBC,box);
 189     do_conect(fn,ndots,xnew);
 190     free_t_atoms(&aaa,FALSE);
 191   }
 192   sfree(xnew);
 193 }
 194
 195 real calc_radius(char *atom)
 196 {
 197   real r;
 198
 199   switch (atom[0]) {
 200   case 'C':
 201     r = 0.16;
 202     break;
 203   case 'O':
 204     r = 0.13;
 205     break;
 206   case 'N':
 207     r = 0.14;
 208     break;
 209   case 'S':
 210     r = 0.2;
 211     break;
 212   case 'H':
 213     r = 0.1;
 214     break;
 215   default:
 216     r = 1e-3;
 217   }
 218   return r;
 219 }
 220
 221 void sas_plot(int nfile,t_filenm fnm[],real solsize,int ndots,
 222               real qcut,bool bSave,real minarea,bool bPBC,
 223               real dgs_default,bool bFindex, const output_env_t oenv)
 224 {
 225   FILE         *fp,*fp2,*fp3=NULL,*vp;
 226   char   *flegend[] = { "Hydrophobic", "Hydrophilic",
 227                               "Total", "D Gsolv" };
 228   char   *vlegend[] = { "Volume (nm\\S3\\N)", "Density (g/l)" };
 229   const char   *vfile;
 230   real         t;
 231   gmx_atomprop_t aps=NULL;
 232   int          status,ndefault;
 233   int          i,j,ii,nfr,natoms,flag,nsurfacedots,res;
 234   rvec         *xtop,*x;
 235   matrix       topbox,box;
 236   t_topology   top;
 237   char         title[STRLEN];
 238   int          ePBC;
 239   bool         bTop;
 240   t_atoms      *atoms;
 241   bool         *bOut,*bPhobic;
 242   bool         bConnelly;
 243   bool         bResAt,bITP,bDGsol;
 244   real         *radius,*dgs_factor=NULL,*area=NULL,*surfacedots=NULL;
 245   real         at_area,*atom_area=NULL,*atom_area2=NULL;
 246   real         *res_a=NULL,*res_area=NULL,*res_area2=NULL;
 247   real         totarea,totvolume,totmass=0,density,harea,tarea,fluc2;
 248   atom_id      **index,*findex;
 249   int          *nx,nphobic,npcheck,retval;
 250   char         **grpname,*fgrpname;
 251   real         dgsolv;
 252
 253   bITP   = opt2bSet("-i",nfile,fnm);
 254   bResAt = opt2bSet("-or",nfile,fnm) || opt2bSet("-oa",nfile,fnm) || bITP;
 255
 256   bTop = read_tps_conf(ftp2fn(efTPS,nfile,fnm),title,&top,&ePBC,
 257                        &xtop,NULL,topbox,FALSE);
 258   atoms = &(top.atoms);
 259
 260   if (!bTop) {
 261     fprintf(stderr,"No tpr file, will not compute Delta G of solvation\n");
 262     bDGsol = FALSE;
 263   } else {
 264     bDGsol = strcmp(*(atoms->atomtype[0]),"?") != 0;
 265     if (!bDGsol) {
 266       fprintf(stderr,"Warning: your tpr file is too old, will not compute "
 267               "Delta G of solvation\n");
 268     } else {
 269       printf("In case you use free energy of solvation predictions:\n");
 270       please_cite(stdout,"Eisenberg86a");
 271     }
 272   }
 273
 274   aps = gmx_atomprop_init();
 275
 276   if ((natoms=read_first_x(oenv,&status,ftp2fn(efTRX,nfile,fnm),
 277                            &t,&x,box))==0)
 278     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 279
 280   if ((ePBC != epbcXYZ) || (TRICLINIC(box))) {
 281     fprintf(stderr,"\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
 282             "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
 283             "will certainly crash the analysis.\n\n");
 284   }
 285   snew(nx,2);
 286   snew(index,2);
 287   snew(grpname,2);
 288   fprintf(stderr,"Select a group for calculation of surface and a group for output:\n");
 289   get_index(atoms,ftp2fn_null(efNDX,nfile,fnm),2,nx,index,grpname);
 290
 291   if (bFindex) {
 292     fprintf(stderr,"Select a group of hydrophobic atoms:\n");
 293     get_index(atoms,ftp2fn_null(efNDX,nfile,fnm),1,&nphobic,&findex,&fgrpname);
 294   }
 295   snew(bOut,natoms);
 296   for(i=0; i<nx[1]; i++)
 297     bOut[index[1][i]] = TRUE;
 298
 299   /* Now compute atomic readii including solvent probe size */
 300   snew(radius,natoms);
 301   snew(bPhobic,nx[0]);
 302   if (bResAt) {
 303     snew(atom_area,nx[0]);
 304     snew(atom_area2,nx[0]);
 305     snew(res_a,atoms->nres);
 306     snew(res_area,atoms->nres);
 307     snew(res_area2,atoms->nres);
 308   }
 309   if (bDGsol)
 310     snew(dgs_factor,nx[0]);
 311
 312   /* Get a Van der Waals radius for each atom */
 313   ndefault = 0;
 314   for(i=0; (i<natoms); i++) {
 315     if (!gmx_atomprop_query(aps,epropVDW,
 316                             *(atoms->resinfo[atoms->atom[i].resind].name),
 317                             *(atoms->atomname[i]),&radius[i]))
 318       ndefault++;
 319     /* radius[i] = calc_radius(*(top->atoms.atomname[i])); */
 320     radius[i] += solsize;
 321   }
 322   if (ndefault > 0)
 323     fprintf(stderr,"WARNING: could not find a Van der Waals radius for %d atoms\n",ndefault);
 324   /* Determine which atom is counted as hydrophobic */
 325   if (bFindex) {
 326     npcheck = 0;
 327     for(i=0; (i<nx[0]); i++) {
 328       ii = index[0][i];
 329       for(j=0; (j<nphobic); j++) {
 330         if (findex[j] == ii) {
 331           bPhobic[i] = TRUE;
 332           if (bOut[ii])
 333             npcheck++;
 334         }
 335       }
 336     }
 337     if (npcheck != nphobic)
 338       gmx_fatal(FARGS,"Consistency check failed: not all %d atoms in the hydrophobic index\n"
 339                   "found in the normal index selection (%d atoms)",nphobic,npcheck);
 340   }
 341   else
 342     nphobic = 0;
 343
 344   for(i=0; (i<nx[0]); i++) {
 345     ii = index[0][i];
 346     if (!bFindex) {
 347       bPhobic[i] = fabs(atoms->atom[ii].q) <= qcut;
 348       if (bPhobic[i] && bOut[ii])
 349         nphobic++;
 350     }
 351     if (bDGsol)
 352       if (!gmx_atomprop_query(aps,epropDGsol,
 353                               *(atoms->resinfo[atoms->atom[ii].resind].name),
 354                               *(atoms->atomtype[ii]),&(dgs_factor[i])))
 355         dgs_factor[i] = dgs_default;
 356     if (debug)
 357       fprintf(debug,"Atom %5d %5s-%5s: q= %6.3f, r= %6.3f, dgsol= %6.3f, hydrophobic= %s\n",
 358               ii+1,*(atoms->resinfo[atoms->atom[ii].resind].name),
 359               *(atoms->atomname[ii]),
 360               atoms->atom[ii].q,radius[ii]-solsize,dgs_factor[i],
 361               BOOL(bPhobic[i]));
 362   }
 363   fprintf(stderr,"%d out of %d atoms were classified as hydrophobic\n",
 364           nphobic,nx[1]);
 365
 366   fp=xvgropen(opt2fn("-o",nfile,fnm),"Solvent Accessible Surface","Time (ps)",
 367               "Area (nm\\S2\\N)",oenv);
 368   xvgr_legend(fp,asize(flegend) - (bDGsol ? 0 : 1),flegend,oenv);
 369   vfile = opt2fn_null("-tv",nfile,fnm);
 370   if (vfile) {
 371     if (!bTop) {
 372       gmx_fatal(FARGS,"Need a tpr file for option -tv");
 373     }
 374     vp=xvgropen(vfile,"Volume and Density","Time (ps)","",oenv);
 375     xvgr_legend(vp,asize(vlegend),vlegend,oenv);
 376     totmass  = 0;
 377     ndefault = 0;
 378     for(i=0; (i<nx[0]); i++) {
 379       real mm;
 380       ii = index[0][i];
 381       /*
 382       if (!query_atomprop(atomprop,epropMass,
 383                           *(top->atoms.resname[top->atoms.atom[ii].resnr]),
 384                           *(top->atoms.atomname[ii]),&mm))
 385         ndefault++;
 386       totmass += mm;
 387       */
 388       totmass += atoms->atom[ii].m;
 389     }
 390     if (ndefault)
 391       fprintf(stderr,"WARNING: Using %d default masses for density calculation, which most likely are inaccurate\n",ndefault);
 392   }
 393   else
 394     vp = NULL;
 395
 396   gmx_atomprop_destroy(aps);
 397
 398   nfr=0;
 399   do {
 400     if (bPBC)
 401       rm_pbc(&top.idef,ePBC,natoms,box,x,x);
 402
 403     bConnelly = (nfr==0 && opt2bSet("-q",nfile,fnm));
 404     if (bConnelly) {
 405       if (!bTop)
 406         gmx_fatal(FARGS,"Need a tpr file for Connelly plot");
 407       flag = FLAG_ATOM_AREA | FLAG_DOTS;
 408     } else {
 409       flag = FLAG_ATOM_AREA;
 410     }
 411     if (vp) {
 412       flag = flag | FLAG_VOLUME;
 413     }
 414
 415     if (debug)
 416       write_sto_conf("check.pdb","pbc check",atoms,x,NULL,ePBC,box);
 417
 418     retval = nsc_dclm_pbc(x,radius,nx[0],ndots,flag,&totarea,
 419                           &area,&totvolume,&surfacedots,&nsurfacedots,
 420                           index[0],ePBC,bPBC ? box : NULL);
 421     if (retval)
 422       gmx_fatal(FARGS,"Something wrong in nsc_dclm_pbc");
 423
 424     if (bConnelly)
 425       connelly_plot(ftp2fn(efPDB,nfile,fnm),
 426                     nsurfacedots,surfacedots,x,atoms,
 427                     &(top.symtab),ePBC,box,bSave);
 428     harea  = 0;
 429     tarea  = 0;
 430     dgsolv = 0;
 431     if (bResAt)
 432       for(i=0; i<atoms->nres; i++)
 433         res_a[i] = 0;
 434     for(i=0; (i<nx[0]); i++) {
 435       ii = index[0][i];
 436       if (bOut[ii]) {
 437         at_area = area[i];
 438         if (bResAt) {
 439           atom_area[i] += at_area;
 440           atom_area2[i] += sqr(at_area);
 441           res_a[atoms->atom[ii].resind] += at_area;
 442         }
 443         tarea += at_area;
 444         if (bDGsol)
 445           dgsolv += at_area*dgs_factor[i];
 446         if (bPhobic[i])
 447           harea += at_area;
 448       }
 449     }
 450     if (bResAt)
 451       for(i=0; i<atoms->nres; i++) {
 452         res_area[i] += res_a[i];
 453         res_area2[i] += sqr(res_a[i]);
 454       }
 455     fprintf(fp,"%10g  %10g  %10g  %10g",t,harea,tarea-harea,tarea);
 456     if (bDGsol)
 457       fprintf(fp,"  %10g\n",dgsolv);
 458     else
 459       fprintf(fp,"\n");
 460
 461     /* Print volume */
 462     if (vp) {
 463       density = totmass*AMU/(totvolume*NANO*NANO*NANO);
 464       fprintf(vp,"%12.5e  %12.5e  %12.5e\n",t,totvolume,density);
 465     }
 466     if (area) {
 467       sfree(area);
 468       area = NULL;
 469     }
 470     if (surfacedots) {
 471       sfree(surfacedots);
 472       surfacedots = NULL;
 473     }
 474     nfr++;
 475   } while (read_next_x(oenv,status,&t,natoms,x,box));
 476
 477   fprintf(stderr,"\n");
 478   close_trj(status);
 479   ffclose(fp);
 480   if (vp)
 481     ffclose(vp);
 482
 483   /* if necessary, print areas per atom to file too: */
 484   if (bResAt) {
 485     for(i=0; i<atoms->nres; i++) {
 486       res_area[i] /= nfr;
 487       res_area2[i] /= nfr;
 488     }
 489     for(i=0; i<nx[0]; i++) {
 490       atom_area[i] /= nfr;
 491       atom_area2[i] /= nfr;
 492     }
 493     fprintf(stderr,"Printing out areas per atom\n");
 494     fp  = xvgropen(opt2fn("-or",nfile,fnm),"Area per residue","Residue",
 495                    "Area (nm\\S2\\N)",oenv);
 496     fp2 = xvgropen(opt2fn("-oa",nfile,fnm),"Area per atom","Atom #",
 497                    "Area (nm\\S2\\N)",oenv);
 498     if (bITP) {
 499       fp3 = ftp2FILE(efITP,nfile,fnm,"w");
 500       fprintf(fp3,"[ position_restraints ]\n"
 501               "#define FCX 1000\n"
 502               "#define FCY 1000\n"
 503               "#define FCZ 1000\n"
 504               "; Atom  Type  fx   fy   fz\n");
 505     }
 506     for(i=0; i<nx[0]; i++) {
 507       ii = index[0][i];
 508       res = atoms->atom[ii].resind;
 509       if (i==nx[0]-1 || res!=atoms->atom[index[0][i+1]].resind) {
 510         fluc2 = res_area2[res]-sqr(res_area[res]);
 511         if (fluc2 < 0)
 512           fluc2 = 0;
 513         fprintf(fp,"%10d  %10g %10g\n",
 514                 atoms->resinfo[res].nr,res_area[res],sqrt(fluc2));
 515       }
 516       fluc2 = atom_area2[i]-sqr(atom_area[i]);
 517       if (fluc2 < 0)
 518         fluc2 = 0;
 519       fprintf(fp2,"%d %g %g\n",index[0][i]+1,atom_area[i],sqrt(fluc2));
 520       if (bITP && (atom_area[i] > minarea))
 521         fprintf(fp3,"%5d   1     FCX  FCX  FCZ\n",ii+1);
 522     }
 523     if (bITP)
 524       ffclose(fp3);
 525     ffclose(fp);
 526   }
 527
 528     /* Be a good citizen, keep our memory free! */
 529     sfree(x);
 530     sfree(nx);
 531     for(i=0;i<2;i++)
 532     {
 533         sfree(index[i]);
 534         sfree(grpname[i]);
 535     }
 536     sfree(bOut);
 537     sfree(radius);
 538     sfree(bPhobic);
 539
 540     if(bResAt)
 541     {
 542         sfree(atom_area);
 543         sfree(atom_area2);
 544         sfree(res_a);
 545         sfree(res_area);
 546         sfree(res_area2);
 547     }
 548     if(bDGsol)
 549     {
 550         sfree(dgs_factor);
 551     }
 552 }
 553
 554 int gmx_sas(int argc,char *argv[])
 555 {
 556   const char *desc[] = {
 557     "g_sas computes hydrophobic, hydrophilic and total solvent accessible surface area.",
 558     "As a side effect the Connolly surface can be generated as well in",
 559     "a pdb file where the nodes are represented as atoms and the vertices",
 560     "connecting the nearest nodes as CONECT records.",
 561     "The program will ask for a group for the surface calculation",
 562     "and a group for the output. The calculation group should always",
 563     "consists of all the non-solvent atoms in the system.",
 564     "The output group can be the whole or part of the calculation group.",
 565     "The area can be plotted",
 566     "per residue and atom as well (options [TT]-or[tt] and [TT]-oa[tt]).",
 567     "In combination with the latter option an [TT]itp[tt] file can be",
 568     "generated (option [TT]-i[tt])",
 569     "which can be used to restrain surface atoms.[PAR]",
 570     "By default, periodic boundary conditions are taken into account,",
 571     "this can be turned off using the [TT]-nopbc[tt] option.[PAR]",
 572     "With the [TT]-tv[tt] option the total volume and density of the molecule can be",
 573     "computed.",
 574     "Please consider whether the normal probe radius is appropriate",
 575     "in this case or whether you would rather use e.g. 0. It is good",
 576     "to keep in mind that the results for volume and density are very",
 577     "approximate, in e.g. ice Ih one can easily fit water molecules in the",
 578     "pores which would yield too low volume, too high surface area and too",
 579     "high density."
 580   };
 581
 582   output_env_t oenv;
 583   static real solsize = 0.14;
 584   static int  ndots   = 24;
 585   static real qcut    = 0.2;
 586   static real minarea = 0.5, dgs_default=0;
 587   static bool bSave   = TRUE,bPBC=TRUE,bFindex=FALSE;
 588   t_pargs pa[] = {
 589     { "-probe", FALSE, etREAL, {&solsize},
 590       "Radius of the solvent probe (nm)" },
 591     { "-ndots",   FALSE, etINT,  {&ndots},
 592       "Number of dots per sphere, more dots means more accuracy" },
 593     { "-qmax",    FALSE, etREAL, {&qcut},
 594       "The maximum charge (e, absolute value) of a hydrophobic atom" },
 595     { "-f_index", FALSE, etBOOL, {&bFindex},
 596       "Determine from a group in the index file what are the hydrophobic atoms rather than from the charge" },
 597     { "-minarea", FALSE, etREAL, {&minarea},
 598       "The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file  (see help)" },
 599     { "-pbc",     FALSE, etBOOL, {&bPBC},
 600       "Take periodicity into account" },
 601     { "-prot",    FALSE, etBOOL, {&bSave},
 602       "Output the protein to the connelly pdb file too" },
 603     { "-dgs",     FALSE, etREAL, {&dgs_default},
 604       "default value for solvation free energy per area (kJ/mol/nm^2)" }
 605   };
 606   t_filenm  fnm[] = {
 607     { efTRX, "-f",   NULL,       ffREAD },
 608     { efTPS, "-s",   NULL,       ffREAD },
 609     { efXVG, "-o",   "area",     ffWRITE },
 610     { efXVG, "-or",  "resarea",  ffOPTWR },
 611     { efXVG, "-oa",  "atomarea", ffOPTWR },
 612     { efXVG, "-tv",  "volume",   ffOPTWR },
 613     { efPDB, "-q",   "connelly", ffOPTWR },
 614     { efNDX, "-n",   "index",    ffOPTRD },
 615     { efITP, "-i",   "surfat",   ffOPTWR }
 616   };
 617 #define NFILE asize(fnm)
 618
 619   CopyRight(stderr,argv[0]);
 620   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 621                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 622   if (solsize < 0) {
 623     solsize=1e-3;
 624     fprintf(stderr,"Probe size too small, setting it to %g\n",solsize);
 625   }
 626   if (ndots < 20) {
 627     ndots = 20;
 628     fprintf(stderr,"Ndots too small, setting it to %d\n",ndots);
 629   }
 630
 631   please_cite(stderr,"Eisenhaber95");
 632
 633   sas_plot(NFILE,fnm,solsize,ndots,qcut,bSave,minarea,bPBC,dgs_default,bFindex,
 634           oenv);
 635
 636   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
 637   do_view(oenv,opt2fn_null("-or",NFILE,fnm),"-nxy");
 638   do_view(oenv,opt2fn_null("-oa",NFILE,fnm),"-nxy");
 639
 640   thanx(stderr);
 641
 642   return 0;
 643 }