src/tools/gmx_sas.c

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