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