src/tools/g_sas.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Great Red Oystrich Makes All Chemists Sane
  28  */
  29 static char *SRCID_g_sas_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <stdlib.h>
  33 #include "sysstuff.h"
  34 #include "string.h"
  35 #include "typedefs.h"
  36 #include "smalloc.h"
  37 #include "macros.h"
  38 #include "vec.h"
  39 #include "xvgr.h"
  40 #include "pbc.h"
  41 #include "copyrite.h"
  42 #include "futil.h"
  43 #include "statutil.h"
  44 #include "rdgroup.h"
  45 #include "nsc.h"
  46 #include "pdbio.h"
  47 #include "confio.h"
  48 #include "rmpbc.h"
  49
  50 void connelly_plot(char *fn,int ndots,real dots[],rvec x[],t_atoms *atoms,
  51                    t_symtab *symtab,matrix box)
  52 {
  53   static char *atomnm="DOT";
  54   static char *resnm ="DOT";
  55   static char *title ="Connely Dot Surface Generated by g_sas";
  56
  57   int  i,i0,ii0,k;
  58   rvec *xnew;
  59
  60   i0 = atoms->nr;
  61   srenew(atoms->atom,atoms->nr+ndots);
  62   srenew(atoms->atomname,atoms->nr+ndots);
  63   srenew(atoms->resname,atoms->nr+ndots);
  64   srenew(atoms->pdbinfo,atoms->nr+ndots);
  65   snew(xnew,atoms->nr+ndots);
  66   for(i=0; (i<atoms->nr); i++)
  67     copy_rvec(x[i],xnew[i]);
  68   for(i=k=0; (i<ndots); i++) {
  69     ii0 = i0+i;
  70     atoms->resname[ii0]  = put_symtab(symtab,resnm);
  71     atoms->atomname[ii0] = put_symtab(symtab,atomnm);
  72     strcpy(atoms->pdbinfo[ii0].pdbresnr,"1");
  73     atoms->pdbinfo[ii0].type = epdbATOM;
  74     atoms->atom[ii0].chain = ' ';
  75     atoms->pdbinfo[ii0].atomnr= ii0;
  76     atoms->atom[ii0].resnr = 1;
  77     xnew[ii0][XX] = dots[k++];
  78     xnew[ii0][YY] = dots[k++];
  79     xnew[ii0][ZZ] = dots[k++];
  80     atoms->pdbinfo[ii0].bfac  = 0.0;
  81     atoms->pdbinfo[ii0].occup = 0.0;
  82   }
  83   atoms->nr = i0+ndots;
  84   write_sto_conf(fn,title,atoms,xnew,NULL,box);
  85   atoms->nr = i0;
  86
  87   sfree(xnew);
  88 }
  89
  90 real calc_radius(char *atom)
  91 {
  92   real r;
  93
  94   switch (atom[0]) {
  95   case 'C':
  96     r = 0.16;
  97     break;
  98   case 'O':
  99     r = 0.13;
 100     break;
 101   case 'N':
 102     r = 0.14;
 103     break;
 104   case 'S':
 105     r = 0.2;
 106     break;
 107   case 'H':
 108     r = 0.1;
 109     break;
 110   default:
 111     r = 1e-3;
 112   }
 113   return r;
 114 }
 115
 116 void sas_plot(int nfile,t_filenm fnm[],real solsize,int ndots,
 117               real qcut,int nskip)
 118 {
 119   FILE         *fp;
 120   real         t;
 121   int          status;
 122   int          i,j,natoms,flag,nsurfacedots;
 123   rvec         *x;
 124   matrix       box;
 125   t_topology   *top;
 126   bool         *bPhobic;
 127   bool         bConnelly;
 128   real         *radius,*area=NULL,*surfacedots=NULL;
 129   real         totarea,totvolume,harea;
 130
 131   if ((natoms=read_first_x(&status,ftp2fn(efTRX,nfile,fnm),
 132                            &t,&x,box))==0)
 133     fatal_error(0,"Could not read coordinates from statusfile\n");
 134   top=read_top(ftp2fn(efTPX,nfile,fnm));
 135
 136   /* Now comput atomic readii including solvent probe size */
 137   snew(radius,natoms);
 138   snew(bPhobic,natoms);
 139   for(i=0; (i<natoms); i++) {
 140     radius[i]  = calc_radius(*(top->atoms.atomname[i])) + solsize;
 141     bPhobic[i] = fabs(top->atoms.atom[i].q) <= qcut;
 142   }
 143   fp=xvgropen(ftp2fn(efXVG,nfile,fnm),"Solvent Accessible Surface","Time (ps)",
 144               "Area (nm\\S2\\N)");
 145
 146   j=0;
 147   do {
 148     if ((j++ % 10) == 0)
 149       fprintf(stderr,"\rframe: %5d",j-1);
 150
 151     if ((nskip > 0) && (((j-1) % nskip) == 0)) {
 152       rm_pbc(&top->idef,natoms,box,x,x);
 153
 154       bConnelly = ((j == 1) && (opt2bSet("-q",nfile,fnm)));
 155       if (bConnelly)
 156         flag = FLAG_ATOM_AREA | FLAG_DOTS;
 157       else
 158         flag = FLAG_ATOM_AREA;
 159
 160       if (NSC(x[0],radius,natoms,ndots,flag,&totarea,
 161               &area,&totvolume,&surfacedots,&nsurfacedots))
 162         fatal_error(0,"Something wrong in NSC");
 163
 164       if (bConnelly)
 165         connelly_plot(ftp2fn(efPDB,nfile,fnm),
 166                       nsurfacedots,surfacedots,x,&(top->atoms),
 167                       &(top->symtab),box);
 168
 169       harea = 0;
 170       for(i=0; (i<natoms); i++) {
 171         if (bPhobic[i])
 172           harea += area[i];
 173       }
 174       fprintf(fp,"%10g  %10g  %10g\n",t,harea,totarea);
 175
 176       if (area)
 177         sfree(area);
 178       if (surfacedots)
 179         sfree(surfacedots);
 180     }
 181   } while (read_next_x(status,&t,natoms,x,box));
 182
 183   fprintf(stderr,"\n");
 184   close_trj(status);
 185
 186   sfree(x);
 187 }
 188
 189 int main(int argc,char *argv[])
 190 {
 191   static char *desc[] = {
 192     "g_sas computes hydrophobic and total solvent accessible surface area."
 193   };
 194
 195   static real solsize = 0.14;
 196   static int  ndots   = 24,nskip=1;
 197   static real qcut    = 0.2;
 198   t_pargs pa[] = {
 199     { "-solsize", FALSE, etREAL, {&solsize},
 200         "Radius of the solvent probe (nm)" },
 201     { "-ndots",   FALSE, etINT,  {&ndots},
 202         "Number of dots per sphere, more dots means more accuracy" },
 203     { "-qmax",    FALSE, etREAL, {&qcut},
 204         "The maximum charge (e, absolute value) of a hydrophobic atom" },
 205     { "-skip",    FALSE, etINT,  {&nskip},
 206       "Do only every nth frame" }
 207   };
 208   t_filenm  fnm[] = {
 209     { efTRX, "-f",   NULL,       ffREAD },
 210     { efTPX, "-s",   NULL,       ffREAD },
 211     { efXVG, "-o",   "area",     ffWRITE },
 212     { efPDB, "-q",   "connelly", ffOPTWR }
 213   };
 214 #define NFILE asize(fnm)
 215
 216   CopyRight(stderr,argv[0]);
 217   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 218                     NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
 219   if (solsize <= 0) {
 220     solsize=1e-3;
 221     fprintf(stderr,"Solsize too small, setting it to %g\n",solsize);
 222   }
 223   if (ndots < 20) {
 224     ndots = 20;
 225     fprintf(stderr,"Ndots too small, setting it to %d\n",ndots);
 226   }
 227
 228   sas_plot(NFILE,fnm,solsize,ndots,qcut,nskip);
 229
 230   xvgr_file(opt2fn("-o",NFILE,fnm),"-nxy");
 231
 232   thanx(stdout);
 233
 234   return 0;
 235 }