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,gmx_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,' ',0,' ');
 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,' ',0,' ');
 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,gmx_bool bSave,real minarea,gmx_bool bPBC,
 223               real dgs_default,gmx_bool bFindex, const output_env_t oenv)
 224 {
 225   FILE         *fp,*fp2,*fp3=NULL,*vp;
 226   const char   *flegend[] = { "Hydrophobic", "Hydrophilic",
 227                               "Total", "D Gsolv" };
 228   const char   *vlegend[] = { "Volume (nm\\S3\\N)", "Density (g/l)" };
 229   const char   *vfile;
 230   real         t;
 231   gmx_atomprop_t aps=NULL;
 232   gmx_rmpbc_t  gpbc=NULL;
 233   t_trxstatus  *status;
 234   int          ndefault;
 235   int          i,j,ii,nfr,natoms,flag,nsurfacedots,res;
 236   rvec         *xtop,*x;
 237   matrix       topbox,box;
 238   t_topology   top;
 239   char         title[STRLEN];
 240   int          ePBC;
 241   gmx_bool         bTop;
 242   t_atoms      *atoms;
 243   gmx_bool         *bOut,*bPhobic;
 244   gmx_bool         bConnelly;
 245   gmx_bool         bResAt,bITP,bDGsol;
 246   real         *radius,*dgs_factor=NULL,*area=NULL,*surfacedots=NULL;
 247   real         at_area,*atom_area=NULL,*atom_area2=NULL;
 248   real         *res_a=NULL,*res_area=NULL,*res_area2=NULL;
 249   real         totarea,totvolume,totmass=0,density,harea,tarea,fluc2;
 250   atom_id      **index,*findex;
 251   int          *nx,nphobic,npcheck,retval;
 252   char         **grpname,*fgrpname;
 253   real         dgsolv;
 254
 255   bITP   = opt2bSet("-i",nfile,fnm);
 256   bResAt = opt2bSet("-or",nfile,fnm) || opt2bSet("-oa",nfile,fnm) || bITP;
 257
 258   bTop = read_tps_conf(ftp2fn(efTPS,nfile,fnm),title,&top,&ePBC,
 259                        &xtop,NULL,topbox,FALSE);
 260   atoms = &(top.atoms);
 261
 262   if (!bTop) {
 263     fprintf(stderr,"No tpr file, will not compute Delta G of solvation\n");
 264     bDGsol = FALSE;
 265   } else {
 266     bDGsol = strcmp(*(atoms->atomtype[0]),"?") != 0;
 267     if (!bDGsol) {
 268       fprintf(stderr,"Warning: your tpr file is too old, will not compute "
 269               "Delta G of solvation\n");
 270     } else {
 271       printf("In case you use free energy of solvation predictions:\n");
 272       please_cite(stdout,"Eisenberg86a");
 273     }
 274   }
 275
 276   aps = gmx_atomprop_init();
 277
 278   if ((natoms=read_first_x(oenv,&status,ftp2fn(efTRX,nfile,fnm),
 279                            &t,&x,box))==0)
 280     gmx_fatal(FARGS,"Could not read coordinates from statusfile\n");
 281
 282   if ((ePBC != epbcXYZ) || (TRICLINIC(box))) {
 283     fprintf(stderr,"\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
 284             "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
 285             "will certainly crash the analysis.\n\n");
 286   }
 287   snew(nx,2);
 288   snew(index,2);
 289   snew(grpname,2);
 290   fprintf(stderr,"Select a group for calculation of surface and a group for output:\n");
 291   get_index(atoms,ftp2fn_null(efNDX,nfile,fnm),2,nx,index,grpname);
 292
 293   if (bFindex) {
 294     fprintf(stderr,"Select a group of hydrophobic atoms:\n");
 295     get_index(atoms,ftp2fn_null(efNDX,nfile,fnm),1,&nphobic,&findex,&fgrpname);
 296   }
 297   snew(bOut,natoms);
 298   for(i=0; i<nx[1]; i++)
 299     bOut[index[1][i]] = TRUE;
 300
 301   /* Now compute atomic readii including solvent probe size */
 302   snew(radius,natoms);
 303   snew(bPhobic,nx[0]);
 304   if (bResAt) {
 305     snew(atom_area,nx[0]);
 306     snew(atom_area2,nx[0]);
 307     snew(res_a,atoms->nres);
 308     snew(res_area,atoms->nres);
 309     snew(res_area2,atoms->nres);
 310   }
 311   if (bDGsol)
 312     snew(dgs_factor,nx[0]);
 313
 314   /* Get a Van der Waals radius for each atom */
 315   ndefault = 0;
 316   for(i=0; (i<natoms); i++) {
 317     if (!gmx_atomprop_query(aps,epropVDW,
 318                             *(atoms->resinfo[atoms->atom[i].resind].name),
 319                             *(atoms->atomname[i]),&radius[i]))
 320       ndefault++;
 321     /* radius[i] = calc_radius(*(top->atoms.atomname[i])); */
 322     radius[i] += solsize;
 323   }
 324   if (ndefault > 0)
 325     fprintf(stderr,"WARNING: could not find a Van der Waals radius for %d atoms\n",ndefault);
 326   /* Determine which atom is counted as hydrophobic */
 327   if (bFindex) {
 328     npcheck = 0;
 329     for(i=0; (i<nx[0]); i++) {
 330       ii = index[0][i];
 331       for(j=0; (j<nphobic); j++) {
 332         if (findex[j] == ii) {
 333           bPhobic[i] = TRUE;
 334           if (bOut[ii])
 335             npcheck++;
 336         }
 337       }
 338     }
 339     if (npcheck != nphobic)
 340       gmx_fatal(FARGS,"Consistency check failed: not all %d atoms in the hydrophobic index\n"
 341                   "found in the normal index selection (%d atoms)",nphobic,npcheck);
 342   }
 343   else
 344     nphobic = 0;
 345
 346   for(i=0; (i<nx[0]); i++) {
 347     ii = index[0][i];
 348     if (!bFindex) {
 349       bPhobic[i] = fabs(atoms->atom[ii].q) <= qcut;
 350       if (bPhobic[i] && bOut[ii])
 351         nphobic++;
 352     }
 353     if (bDGsol)
 354       if (!gmx_atomprop_query(aps,epropDGsol,
 355                               *(atoms->resinfo[atoms->atom[ii].resind].name),
 356                               *(atoms->atomtype[ii]),&(dgs_factor[i])))
 357         dgs_factor[i] = dgs_default;
 358     if (debug)
 359       fprintf(debug,"Atom %5d %5s-%5s: q= %6.3f, r= %6.3f, dgsol= %6.3f, hydrophobic= %s\n",
 360               ii+1,*(atoms->resinfo[atoms->atom[ii].resind].name),
 361               *(atoms->atomname[ii]),
 362               atoms->atom[ii].q,radius[ii]-solsize,dgs_factor[i],
 363               BOOL(bPhobic[i]));
 364   }
 365   fprintf(stderr,"%d out of %d atoms were classified as hydrophobic\n",
 366           nphobic,nx[1]);
 367
 368   fp=xvgropen(opt2fn("-o",nfile,fnm),"Solvent Accessible Surface","Time (ps)",
 369               "Area (nm\\S2\\N)",oenv);
 370   xvgr_legend(fp,asize(flegend) - (bDGsol ? 0 : 1),flegend,oenv);
 371   vfile = opt2fn_null("-tv",nfile,fnm);
 372   if (vfile) {
 373     if (!bTop) {
 374       gmx_fatal(FARGS,"Need a tpr file for option -tv");
 375     }
 376     vp=xvgropen(vfile,"Volume and Density","Time (ps)","",oenv);
 377     xvgr_legend(vp,asize(vlegend),vlegend,oenv);
 378     totmass  = 0;
 379     ndefault = 0;
 380     for(i=0; (i<nx[0]); i++) {
 381       real mm;
 382       ii = index[0][i];
 383       /*
 384       if (!query_atomprop(atomprop,epropMass,
 385                           *(top->atoms.resname[top->atoms.atom[ii].resnr]),
 386                           *(top->atoms.atomname[ii]),&mm))
 387         ndefault++;
 388       totmass += mm;
 389       */
 390       totmass += atoms->atom[ii].m;
 391     }
 392     if (ndefault)
 393       fprintf(stderr,"WARNING: Using %d default masses for density calculation, which most likely are inaccurate\n",ndefault);
 394   }
 395   else
 396     vp = NULL;
 397
 398   gmx_atomprop_destroy(aps);
 399
 400   if (bPBC)
 401     gpbc = gmx_rmpbc_init(&top.idef,ePBC,natoms,box);
 402
 403   nfr=0;
 404   do {
 405     if (bPBC)
 406       gmx_rmpbc(gpbc,natoms,box,x);
 407
 408     bConnelly = (nfr==0 && opt2bSet("-q",nfile,fnm));
 409     if (bConnelly) {
 410       if (!bTop)
 411         gmx_fatal(FARGS,"Need a tpr file for Connelly plot");
 412       flag = FLAG_ATOM_AREA | FLAG_DOTS;
 413     } else {
 414       flag = FLAG_ATOM_AREA;
 415     }
 416     if (vp) {
 417       flag = flag | FLAG_VOLUME;
 418     }
 419
 420     if (debug)
 421       write_sto_conf("check.pdb","pbc check",atoms,x,NULL,ePBC,box);
 422
 423     retval = nsc_dclm_pbc(x,radius,nx[0],ndots,flag,&totarea,
 424                           &area,&totvolume,&surfacedots,&nsurfacedots,
 425                           index[0],ePBC,bPBC ? box : NULL);
 426     if (retval)
 427       gmx_fatal(FARGS,"Something wrong in nsc_dclm_pbc");
 428
 429     if (bConnelly)
 430       connelly_plot(ftp2fn(efPDB,nfile,fnm),
 431                     nsurfacedots,surfacedots,x,atoms,
 432                     &(top.symtab),ePBC,box,bSave);
 433     harea  = 0;
 434     tarea  = 0;
 435     dgsolv = 0;
 436     if (bResAt)
 437       for(i=0; i<atoms->nres; i++)
 438         res_a[i] = 0;
 439     for(i=0; (i<nx[0]); i++) {
 440       ii = index[0][i];
 441       if (bOut[ii]) {
 442         at_area = area[i];
 443         if (bResAt) {
 444           atom_area[i] += at_area;
 445           atom_area2[i] += sqr(at_area);
 446           res_a[atoms->atom[ii].resind] += at_area;
 447         }
 448         tarea += at_area;
 449         if (bDGsol)
 450           dgsolv += at_area*dgs_factor[i];
 451         if (bPhobic[i])
 452           harea += at_area;
 453       }
 454     }
 455     if (bResAt)
 456       for(i=0; i<atoms->nres; i++) {
 457         res_area[i] += res_a[i];
 458         res_area2[i] += sqr(res_a[i]);
 459       }
 460     fprintf(fp,"%10g  %10g  %10g  %10g",t,harea,tarea-harea,tarea);
 461     if (bDGsol)
 462       fprintf(fp,"  %10g\n",dgsolv);
 463     else
 464       fprintf(fp,"\n");
 465
 466     /* Print volume */
 467     if (vp) {
 468       density = totmass*AMU/(totvolume*NANO*NANO*NANO);
 469       fprintf(vp,"%12.5e  %12.5e  %12.5e\n",t,totvolume,density);
 470     }
 471     if (area) {
 472       sfree(area);
 473       area = NULL;
 474     }
 475     if (surfacedots) {
 476       sfree(surfacedots);
 477       surfacedots = NULL;
 478     }
 479     nfr++;
 480   } while (read_next_x(oenv,status,&t,natoms,x,box));
 481
 482   if (bPBC)
 483     gmx_rmpbc_done(gpbc);
 484
 485   fprintf(stderr,"\n");
 486   close_trj(status);
 487   ffclose(fp);
 488   if (vp)
 489     ffclose(vp);
 490
 491   /* if necessary, print areas per atom to file too: */
 492   if (bResAt) {
 493     for(i=0; i<atoms->nres; i++) {
 494       res_area[i] /= nfr;
 495       res_area2[i] /= nfr;
 496     }
 497     for(i=0; i<nx[0]; i++) {
 498       atom_area[i] /= nfr;
 499       atom_area2[i] /= nfr;
 500     }
 501     fprintf(stderr,"Printing out areas per atom\n");
 502     fp  = xvgropen(opt2fn("-or",nfile,fnm),"Area per residue","Residue",
 503                    "Area (nm\\S2\\N)",oenv);
 504     fp2 = xvgropen(opt2fn("-oa",nfile,fnm),"Area per atom","Atom #",
 505                    "Area (nm\\S2\\N)",oenv);
 506     if (bITP) {
 507       fp3 = ftp2FILE(efITP,nfile,fnm,"w");
 508       fprintf(fp3,"[ position_restraints ]\n"
 509               "#define FCX 1000\n"
 510               "#define FCY 1000\n"
 511               "#define FCZ 1000\n"
 512               "; Atom  Type  fx   fy   fz\n");
 513     }
 514     for(i=0; i<nx[0]; i++) {
 515       ii = index[0][i];
 516       res = atoms->atom[ii].resind;
 517       if (i==nx[0]-1 || res!=atoms->atom[index[0][i+1]].resind) {
 518         fluc2 = res_area2[res]-sqr(res_area[res]);
 519         if (fluc2 < 0)
 520           fluc2 = 0;
 521         fprintf(fp,"%10d  %10g %10g\n",
 522                 atoms->resinfo[res].nr,res_area[res],sqrt(fluc2));
 523       }
 524       fluc2 = atom_area2[i]-sqr(atom_area[i]);
 525       if (fluc2 < 0)
 526         fluc2 = 0;
 527       fprintf(fp2,"%d %g %g\n",index[0][i]+1,atom_area[i],sqrt(fluc2));
 528       if (bITP && (atom_area[i] > minarea))
 529         fprintf(fp3,"%5d   1     FCX  FCX  FCZ\n",ii+1);
 530     }
 531     if (bITP)
 532       ffclose(fp3);
 533     ffclose(fp);
 534   }
 535
 536     /* Be a good citizen, keep our memory free! */
 537     sfree(x);
 538     sfree(nx);
 539     for(i=0;i<2;i++)
 540     {
 541         sfree(index[i]);
 542         sfree(grpname[i]);
 543     }
 544     sfree(bOut);
 545     sfree(radius);
 546     sfree(bPhobic);
 547
 548     if(bResAt)
 549     {
 550         sfree(atom_area);
 551         sfree(atom_area2);
 552         sfree(res_a);
 553         sfree(res_area);
 554         sfree(res_area2);
 555     }
 556     if(bDGsol)
 557     {
 558         sfree(dgs_factor);
 559     }
 560 }
 561
 562 int gmx_sas(int argc,char *argv[])
 563 {
 564   const char *desc[] = {
 565     "g_sas computes hydrophobic, hydrophilic and total solvent accessible surface area.",
 566     "As a side effect the Connolly surface can be generated as well in",
 567     "a pdb file where the nodes are represented as atoms and the vertices",
 568     "connecting the nearest nodes as CONECT records.",
 569     "The program will ask for a group for the surface calculation",
 570     "and a group for the output. The calculation group should always",
 571     "consists of all the non-solvent atoms in the system.",
 572     "The output group can be the whole or part of the calculation group.",
 573     "The area can be plotted",
 574     "per residue and atom as well (options [TT]-or[tt] and [TT]-oa[tt]).",
 575     "In combination with the latter option an [TT]itp[tt] file can be",
 576     "generated (option [TT]-i[tt])",
 577     "which can be used to restrain surface atoms.[PAR]",
 578     "By default, periodic boundary conditions are taken into account,",
 579     "this can be turned off using the [TT]-nopbc[tt] option.[PAR]",
 580     "With the [TT]-tv[tt] option the total volume and density of the molecule can be",
 581     "computed.",
 582     "Please consider whether the normal probe radius is appropriate",
 583     "in this case or whether you would rather use e.g. 0. It is good",
 584     "to keep in mind that the results for volume and density are very",
 585     "approximate, in e.g. ice Ih one can easily fit water molecules in the",
 586     "pores which would yield too low volume, too high surface area and too",
 587     "high density."
 588   };
 589
 590   output_env_t oenv;
 591   static real solsize = 0.14;
 592   static int  ndots   = 24;
 593   static real qcut    = 0.2;
 594   static real minarea = 0.5, dgs_default=0;
 595   static gmx_bool bSave   = TRUE,bPBC=TRUE,bFindex=FALSE;
 596   t_pargs pa[] = {
 597     { "-probe", FALSE, etREAL, {&solsize},
 598       "Radius of the solvent probe (nm)" },
 599     { "-ndots",   FALSE, etINT,  {&ndots},
 600       "Number of dots per sphere, more dots means more accuracy" },
 601     { "-qmax",    FALSE, etREAL, {&qcut},
 602       "The maximum charge (e, absolute value) of a hydrophobic atom" },
 603     { "-f_index", FALSE, etBOOL, {&bFindex},
 604       "Determine from a group in the index file what are the hydrophobic atoms rather than from the charge" },
 605     { "-minarea", FALSE, etREAL, {&minarea},
 606       "The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file  (see help)" },
 607     { "-pbc",     FALSE, etBOOL, {&bPBC},
 608       "Take periodicity into account" },
 609     { "-prot",    FALSE, etBOOL, {&bSave},
 610       "Output the protein to the connelly pdb file too" },
 611     { "-dgs",     FALSE, etREAL, {&dgs_default},
 612       "default value for solvation free energy per area (kJ/mol/nm^2)" }
 613   };
 614   t_filenm  fnm[] = {
 615     { efTRX, "-f",   NULL,       ffREAD },
 616     { efTPS, "-s",   NULL,       ffREAD },
 617     { efXVG, "-o",   "area",     ffWRITE },
 618     { efXVG, "-or",  "resarea",  ffOPTWR },
 619     { efXVG, "-oa",  "atomarea", ffOPTWR },
 620     { efXVG, "-tv",  "volume",   ffOPTWR },
 621     { efPDB, "-q",   "connelly", ffOPTWR },
 622     { efNDX, "-n",   "index",    ffOPTRD },
 623     { efITP, "-i",   "surfat",   ffOPTWR }
 624   };
 625 #define NFILE asize(fnm)
 626
 627   CopyRight(stderr,argv[0]);
 628   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 629                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
 630   if (solsize < 0) {
 631     solsize=1e-3;
 632     fprintf(stderr,"Probe size too small, setting it to %g\n",solsize);
 633   }
 634   if (ndots < 20) {
 635     ndots = 20;
 636     fprintf(stderr,"Ndots too small, setting it to %d\n",ndots);
 637   }
 638
 639   please_cite(stderr,"Eisenhaber95");
 640
 641   sas_plot(NFILE,fnm,solsize,ndots,qcut,bSave,minarea,bPBC,dgs_default,bFindex,
 642           oenv);
 643
 644   do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
 645   do_view(oenv,opt2fn_null("-or",NFILE,fnm),"-nxy");
 646   do_view(oenv,opt2fn_null("-oa",NFILE,fnm),"-nxy");
 647
 648   thanx(stderr);
 649
 650   return 0;
 651 }