src/gromacs/gmxana/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 "gromacs/fileio/futil.h"
  51 #include "statutil.h"
  52 #include "index.h"
  53 #include "nsc.h"
  54 #include "gromacs/fileio/pdbio.h"
  55 #include "gromacs/fileio/confio.h"
  56 #include "rmpbc.h"
  57 #include "names.h"
  58 #include "atomprop.h"
  59 #include "physics.h"
  60 #include "gromacs/fileio/tpxio.h"
  61 #include "gromacs/fileio/trxio.h"
  62 #include "gmx_ana.h"
  63
  64
  65 typedef struct {
  66     atom_id  aa, ab;
  67     real     d2a, d2b;
  68 } t_conect;
  69
  70 void add_rec(t_conect c[], atom_id i, atom_id j, real d2)
  71 {
  72     if (c[i].aa == NO_ATID)
  73     {
  74         c[i].aa  = j;
  75         c[i].d2a = d2;
  76     }
  77     else if (c[i].ab == NO_ATID)
  78     {
  79         c[i].ab  = j;
  80         c[i].d2b = d2;
  81     }
  82     else if (d2 < c[i].d2a)
  83     {
  84         c[i].aa  = j;
  85         c[i].d2a = d2;
  86     }
  87     else if (d2 < c[i].d2b)
  88     {
  89         c[i].ab  = j;
  90         c[i].d2b = d2;
  91     }
  92     /* Swap them if necessary: a must be larger than b */
  93     if (c[i].d2a < c[i].d2b)
  94     {
  95         j        = c[i].ab;
  96         c[i].ab  = c[i].aa;
  97         c[i].aa  = j;
  98         d2       = c[i].d2b;
  99         c[i].d2b = c[i].d2a;
 100         c[i].d2a = d2;
 101     }
 102 }
 103
 104 void do_conect(const char *fn, int n, rvec x[])
 105 {
 106     FILE     *fp;
 107     int       i, j;
 108     t_conect *c;
 109     rvec      dx;
 110     real      d2;
 111
 112     fprintf(stderr, "Building CONECT records\n");
 113     snew(c, n);
 114     for (i = 0; (i < n); i++)
 115     {
 116         c[i].aa = c[i].ab = NO_ATID;
 117     }
 118
 119     for (i = 0; (i < n); i++)
 120     {
 121         for (j = i+1; (j < n); j++)
 122         {
 123             rvec_sub(x[i], x[j], dx);
 124             d2 = iprod(dx, dx);
 125             add_rec(c, i, j, d2);
 126             add_rec(c, j, i, d2);
 127         }
 128     }
 129     fp = ffopen(fn, "a");
 130     for (i = 0; (i < n); i++)
 131     {
 132         if ((c[i].aa == NO_ATID) || (c[i].ab == NO_ATID))
 133         {
 134             fprintf(stderr, "Warning dot %d has no conections\n", i+1);
 135         }
 136         fprintf(fp, "CONECT%5d%5d%5d\n", i+1, c[i].aa+1, c[i].ab+1);
 137     }
 138     ffclose(fp);
 139     sfree(c);
 140 }
 141
 142 void connelly_plot(const char *fn, int ndots, real dots[], rvec x[], t_atoms *atoms,
 143                    t_symtab *symtab, int ePBC, matrix box, gmx_bool bSave)
 144 {
 145     static const char *atomnm = "DOT";
 146     static const char *resnm  = "DOT";
 147     static const char *title  = "Connely Dot Surface Generated by g_sas";
 148
 149     int                i, i0, r0, ii0, k;
 150     rvec              *xnew;
 151     t_atoms            aaa;
 152
 153     if (bSave)
 154     {
 155         i0 = atoms->nr;
 156         r0 = atoms->nres;
 157         srenew(atoms->atom, atoms->nr+ndots);
 158         srenew(atoms->atomname, atoms->nr+ndots);
 159         srenew(atoms->resinfo, r0+1);
 160         atoms->atom[i0].resind = r0;
 161         t_atoms_set_resinfo(atoms, i0, symtab, resnm, r0+1, ' ', 0, ' ');
 162         srenew(atoms->pdbinfo, atoms->nr+ndots);
 163         snew(xnew, atoms->nr+ndots);
 164         for (i = 0; (i < atoms->nr); i++)
 165         {
 166             copy_rvec(x[i], xnew[i]);
 167         }
 168         for (i = k = 0; (i < ndots); i++)
 169         {
 170             ii0                        = i0+i;
 171             atoms->atomname[ii0]       = put_symtab(symtab, atomnm);
 172             atoms->pdbinfo[ii0].type   = epdbATOM;
 173             atoms->pdbinfo[ii0].atomnr = ii0;
 174             atoms->atom[ii0].resind    = r0;
 175             xnew[ii0][XX]              = dots[k++];
 176             xnew[ii0][YY]              = dots[k++];
 177             xnew[ii0][ZZ]              = dots[k++];
 178             atoms->pdbinfo[ii0].bfac   = 0.0;
 179             atoms->pdbinfo[ii0].occup  = 0.0;
 180         }
 181         atoms->nr   = i0+ndots;
 182         atoms->nres = r0+1;
 183         write_sto_conf(fn, title, atoms, xnew, NULL, ePBC, box);
 184         atoms->nres = r0;
 185         atoms->nr   = i0;
 186     }
 187     else
 188     {
 189         init_t_atoms(&aaa, ndots, TRUE);
 190         aaa.atom[0].resind = 0;
 191         t_atoms_set_resinfo(&aaa, 0, symtab, resnm, 1, ' ', 0, ' ');
 192         snew(xnew, ndots);
 193         for (i = k = 0; (i < ndots); i++)
 194         {
 195             ii0                     = i;
 196             aaa.atomname[ii0]       = put_symtab(symtab, atomnm);
 197             aaa.pdbinfo[ii0].type   = epdbATOM;
 198             aaa.pdbinfo[ii0].atomnr = ii0;
 199             aaa.atom[ii0].resind    = 0;
 200             xnew[ii0][XX]           = dots[k++];
 201             xnew[ii0][YY]           = dots[k++];
 202             xnew[ii0][ZZ]           = dots[k++];
 203             aaa.pdbinfo[ii0].bfac   = 0.0;
 204             aaa.pdbinfo[ii0].occup  = 0.0;
 205         }
 206         aaa.nr = ndots;
 207         write_sto_conf(fn, title, &aaa, xnew, NULL, ePBC, box);
 208         do_conect(fn, ndots, xnew);
 209         free_t_atoms(&aaa, FALSE);
 210     }
 211     sfree(xnew);
 212 }
 213
 214 real calc_radius(char *atom)
 215 {
 216     real r;
 217
 218     switch (atom[0])
 219     {
 220         case 'C':
 221             r = 0.16;
 222             break;
 223         case 'O':
 224             r = 0.13;
 225             break;
 226         case 'N':
 227             r = 0.14;
 228             break;
 229         case 'S':
 230             r = 0.2;
 231             break;
 232         case 'H':
 233             r = 0.1;
 234             break;
 235         default:
 236             r = 1e-3;
 237     }
 238     return r;
 239 }
 240
 241 void sas_plot(int nfile, t_filenm fnm[], real solsize, int ndots,
 242               real qcut, gmx_bool bSave, real minarea, gmx_bool bPBC,
 243               real dgs_default, gmx_bool bFindex, const output_env_t oenv)
 244 {
 245     FILE             *fp, *fp2, *fp3 = NULL, *vp;
 246     const char       *flegend[] = {
 247         "Hydrophobic", "Hydrophilic",
 248         "Total", "D Gsolv"
 249     };
 250     const char       *vlegend[]          = { "Volume (nm\\S3\\N)", "Density (g/l)" };
 251     const char       *or_and_oa_legend[] = { "Average (nm\\S2\\N)", "Standard deviation (nm\\S2\\N)" };
 252     const char       *vfile;
 253     real              t;
 254     gmx_atomprop_t    aps  = NULL;
 255     gmx_rmpbc_t       gpbc = NULL;
 256     t_trxstatus      *status;
 257     int               ndefault;
 258     int               i, j, ii, nfr, natoms, flag, nsurfacedots, res;
 259     rvec             *xtop, *x;
 260     matrix            topbox, box;
 261     t_topology        top;
 262     char              title[STRLEN];
 263     int               ePBC;
 264     gmx_bool          bTop;
 265     t_atoms          *atoms;
 266     gmx_bool         *bOut, *bPhobic;
 267     gmx_bool          bConnelly;
 268     gmx_bool          bResAt, bITP, bDGsol;
 269     real             *radius, *dgs_factor = NULL, *area = NULL, *surfacedots = NULL;
 270     real              at_area, *atom_area = NULL, *atom_area2 = NULL;
 271     real             *res_a = NULL, *res_area = NULL, *res_area2 = NULL;
 272     real              totarea, totvolume, totmass = 0, density, harea, tarea, fluc2;
 273     atom_id         **index, *findex;
 274     int              *nx, nphobic, npcheck, retval;
 275     char            **grpname, *fgrpname;
 276     real              dgsolv;
 277
 278     bITP   = opt2bSet("-i", nfile, fnm);
 279     bResAt = opt2bSet("-or", nfile, fnm) || opt2bSet("-oa", nfile, fnm) || bITP;
 280
 281     bTop = read_tps_conf(ftp2fn(efTPS, nfile, fnm), title, &top, &ePBC,
 282                          &xtop, NULL, topbox, FALSE);
 283     atoms = &(top.atoms);
 284
 285     if (!bTop)
 286     {
 287         fprintf(stderr, "No tpr file, will not compute Delta G of solvation\n");
 288         bDGsol = FALSE;
 289     }
 290     else
 291     {
 292         bDGsol = strcmp(*(atoms->atomtype[0]), "?") != 0;
 293         if (!bDGsol)
 294         {
 295             fprintf(stderr, "Warning: your tpr file is too old, will not compute "
 296                     "Delta G of solvation\n");
 297         }
 298         else
 299         {
 300             printf("In case you use free energy of solvation predictions:\n");
 301             please_cite(stdout, "Eisenberg86a");
 302         }
 303     }
 304
 305     aps = gmx_atomprop_init();
 306
 307     if ((natoms = read_first_x(oenv, &status, ftp2fn(efTRX, nfile, fnm),
 308                                &t, &x, box)) == 0)
 309     {
 310         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 311     }
 312
 313     if ((ePBC != epbcXYZ) || (TRICLINIC(box)))
 314     {
 315         fprintf(stderr, "\n\nWARNING: non-rectangular boxes may give erroneous results or crashes.\n"
 316                 "Analysis based on vacuum simulations (with the possibility of evaporation)\n"
 317                 "will certainly crash the analysis.\n\n");
 318     }
 319     snew(nx, 2);
 320     snew(index, 2);
 321     snew(grpname, 2);
 322     fprintf(stderr, "Select a group for calculation of surface and a group for output:\n");
 323     get_index(atoms, ftp2fn_null(efNDX, nfile, fnm), 2, nx, index, grpname);
 324
 325     if (bFindex)
 326     {
 327         fprintf(stderr, "Select a group of hydrophobic atoms:\n");
 328         get_index(atoms, ftp2fn_null(efNDX, nfile, fnm), 1, &nphobic, &findex, &fgrpname);
 329     }
 330     snew(bOut, natoms);
 331     for (i = 0; i < nx[1]; i++)
 332     {
 333         bOut[index[1][i]] = TRUE;
 334     }
 335
 336     /* Now compute atomic readii including solvent probe size */
 337     snew(radius, natoms);
 338     snew(bPhobic, nx[0]);
 339     if (bResAt)
 340     {
 341         snew(atom_area, nx[0]);
 342         snew(atom_area2, nx[0]);
 343         snew(res_a, atoms->nres);
 344         snew(res_area, atoms->nres);
 345         snew(res_area2, atoms->nres);
 346     }
 347     if (bDGsol)
 348     {
 349         snew(dgs_factor, nx[0]);
 350     }
 351
 352     /* Get a Van der Waals radius for each atom */
 353     ndefault = 0;
 354     for (i = 0; (i < natoms); i++)
 355     {
 356         if (!gmx_atomprop_query(aps, epropVDW,
 357                                 *(atoms->resinfo[atoms->atom[i].resind].name),
 358                                 *(atoms->atomname[i]), &radius[i]))
 359         {
 360             ndefault++;
 361         }
 362         /* radius[i] = calc_radius(*(top->atoms.atomname[i])); */
 363         radius[i] += solsize;
 364     }
 365     if (ndefault > 0)
 366     {
 367         fprintf(stderr, "WARNING: could not find a Van der Waals radius for %d atoms\n", ndefault);
 368     }
 369     /* Determine which atom is counted as hydrophobic */
 370     if (bFindex)
 371     {
 372         npcheck = 0;
 373         for (i = 0; (i < nx[0]); i++)
 374         {
 375             ii = index[0][i];
 376             for (j = 0; (j < nphobic); j++)
 377             {
 378                 if (findex[j] == ii)
 379                 {
 380                     bPhobic[i] = TRUE;
 381                     if (bOut[ii])
 382                     {
 383                         npcheck++;
 384                     }
 385                 }
 386             }
 387         }
 388         if (npcheck != nphobic)
 389         {
 390             gmx_fatal(FARGS, "Consistency check failed: not all %d atoms in the hydrophobic index\n"
 391                       "found in the normal index selection (%d atoms)", nphobic, npcheck);
 392         }
 393     }
 394     else
 395     {
 396         nphobic = 0;
 397     }
 398
 399     for (i = 0; (i < nx[0]); i++)
 400     {
 401         ii = index[0][i];
 402         if (!bFindex)
 403         {
 404             bPhobic[i] = fabs(atoms->atom[ii].q) <= qcut;
 405             if (bPhobic[i] && bOut[ii])
 406             {
 407                 nphobic++;
 408             }
 409         }
 410         if (bDGsol)
 411         {
 412             if (!gmx_atomprop_query(aps, epropDGsol,
 413                                     *(atoms->resinfo[atoms->atom[ii].resind].name),
 414                                     *(atoms->atomtype[ii]), &(dgs_factor[i])))
 415             {
 416                 dgs_factor[i] = dgs_default;
 417             }
 418         }
 419         if (debug)
 420         {
 421             fprintf(debug, "Atom %5d %5s-%5s: q= %6.3f, r= %6.3f, dgsol= %6.3f, hydrophobic= %s\n",
 422                     ii+1, *(atoms->resinfo[atoms->atom[ii].resind].name),
 423                     *(atoms->atomname[ii]),
 424                     atoms->atom[ii].q, radius[ii]-solsize, dgs_factor[i],
 425                     EBOOL(bPhobic[i]));
 426         }
 427     }
 428     fprintf(stderr, "%d out of %d atoms were classified as hydrophobic\n",
 429             nphobic, nx[1]);
 430
 431     fp = xvgropen(opt2fn("-o", nfile, fnm), "Solvent Accessible Surface", "Time (ps)",
 432                   "Area (nm\\S2\\N)", oenv);
 433     xvgr_legend(fp, asize(flegend) - (bDGsol ? 0 : 1), flegend, oenv);
 434     vfile = opt2fn_null("-tv", nfile, fnm);
 435     if (vfile)
 436     {
 437         if (!bTop)
 438         {
 439             gmx_fatal(FARGS, "Need a tpr file for option -tv");
 440         }
 441         vp = xvgropen(vfile, "Volume and Density", "Time (ps)", "", oenv);
 442         xvgr_legend(vp, asize(vlegend), vlegend, oenv);
 443         totmass  = 0;
 444         ndefault = 0;
 445         for (i = 0; (i < nx[0]); i++)
 446         {
 447             real mm;
 448             ii = index[0][i];
 449             /*
 450                if (!query_atomprop(atomprop,epropMass,
 451              *(top->atoms.resname[top->atoms.atom[ii].resnr]),
 452              *(top->atoms.atomname[ii]),&mm))
 453                ndefault++;
 454                totmass += mm;
 455              */
 456             totmass += atoms->atom[ii].m;
 457         }
 458         if (ndefault)
 459         {
 460             fprintf(stderr, "WARNING: Using %d default masses for density calculation, which most likely are inaccurate\n", ndefault);
 461         }
 462     }
 463     else
 464     {
 465         vp = NULL;
 466     }
 467
 468     gmx_atomprop_destroy(aps);
 469
 470     if (bPBC)
 471     {
 472         gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
 473     }
 474
 475     nfr = 0;
 476     do
 477     {
 478         if (bPBC)
 479         {
 480             gmx_rmpbc(gpbc, natoms, box, x);
 481         }
 482
 483         bConnelly = (nfr == 0 && opt2bSet("-q", nfile, fnm));
 484         if (bConnelly)
 485         {
 486             if (!bTop)
 487             {
 488                 gmx_fatal(FARGS, "Need a tpr file for Connelly plot");
 489             }
 490             flag = FLAG_ATOM_AREA | FLAG_DOTS;
 491         }
 492         else
 493         {
 494             flag = FLAG_ATOM_AREA;
 495         }
 496         if (vp)
 497         {
 498             flag = flag | FLAG_VOLUME;
 499         }
 500
 501         if (debug)
 502         {
 503             write_sto_conf("check.pdb", "pbc check", atoms, x, NULL, ePBC, box);
 504         }
 505
 506         retval = nsc_dclm_pbc(x, radius, nx[0], ndots, flag, &totarea,
 507                               &area, &totvolume, &surfacedots, &nsurfacedots,
 508                               index[0], ePBC, bPBC ? box : NULL);
 509         if (retval)
 510         {
 511             gmx_fatal(FARGS, "Something wrong in nsc_dclm_pbc");
 512         }
 513
 514         if (bConnelly)
 515         {
 516             connelly_plot(ftp2fn(efPDB, nfile, fnm),
 517                           nsurfacedots, surfacedots, x, atoms,
 518                           &(top.symtab), ePBC, box, bSave);
 519         }
 520         harea  = 0;
 521         tarea  = 0;
 522         dgsolv = 0;
 523         if (bResAt)
 524         {
 525             for (i = 0; i < atoms->nres; i++)
 526             {
 527                 res_a[i] = 0;
 528             }
 529         }
 530         for (i = 0; (i < nx[0]); i++)
 531         {
 532             ii = index[0][i];
 533             if (bOut[ii])
 534             {
 535                 at_area = area[i];
 536                 if (bResAt)
 537                 {
 538                     atom_area[i]                  += at_area;
 539                     atom_area2[i]                 += sqr(at_area);
 540                     res_a[atoms->atom[ii].resind] += at_area;
 541                 }
 542                 tarea += at_area;
 543                 if (bDGsol)
 544                 {
 545                     dgsolv += at_area*dgs_factor[i];
 546                 }
 547                 if (bPhobic[i])
 548                 {
 549                     harea += at_area;
 550                 }
 551             }
 552         }
 553         if (bResAt)
 554         {
 555             for (i = 0; i < atoms->nres; i++)
 556             {
 557                 res_area[i]  += res_a[i];
 558                 res_area2[i] += sqr(res_a[i]);
 559             }
 560         }
 561         fprintf(fp, "%10g  %10g  %10g  %10g", t, harea, tarea-harea, tarea);
 562         if (bDGsol)
 563         {
 564             fprintf(fp, "  %10g\n", dgsolv);
 565         }
 566         else
 567         {
 568             fprintf(fp, "\n");
 569         }
 570
 571         /* Print volume */
 572         if (vp)
 573         {
 574             density = totmass*AMU/(totvolume*NANO*NANO*NANO);
 575             fprintf(vp, "%12.5e  %12.5e  %12.5e\n", t, totvolume, density);
 576         }
 577         if (area)
 578         {
 579             sfree(area);
 580             area = NULL;
 581         }
 582         if (surfacedots)
 583         {
 584             sfree(surfacedots);
 585             surfacedots = NULL;
 586         }
 587         nfr++;
 588     }
 589     while (read_next_x(oenv, status, &t, x, box));
 590
 591     if (bPBC)
 592     {
 593         gmx_rmpbc_done(gpbc);
 594     }
 595
 596     fprintf(stderr, "\n");
 597     close_trj(status);
 598     ffclose(fp);
 599     if (vp)
 600     {
 601         ffclose(vp);
 602     }
 603
 604     /* if necessary, print areas per atom to file too: */
 605     if (bResAt)
 606     {
 607         for (i = 0; i < atoms->nres; i++)
 608         {
 609             res_area[i]  /= nfr;
 610             res_area2[i] /= nfr;
 611         }
 612         for (i = 0; i < nx[0]; i++)
 613         {
 614             atom_area[i]  /= nfr;
 615             atom_area2[i] /= nfr;
 616         }
 617         fprintf(stderr, "Printing out areas per atom\n");
 618         fp  = xvgropen(opt2fn("-or", nfile, fnm), "Area per residue over the trajectory", "Residue",
 619                        "Area (nm\\S2\\N)", oenv);
 620         xvgr_legend(fp, asize(or_and_oa_legend), or_and_oa_legend, oenv);
 621         fp2 = xvgropen(opt2fn("-oa", nfile, fnm), "Area per atom over the trajectory", "Atom #",
 622                        "Area (nm\\S2\\N)", oenv);
 623         xvgr_legend(fp2, asize(or_and_oa_legend), or_and_oa_legend, oenv);
 624         if (bITP)
 625         {
 626             fp3 = ftp2FILE(efITP, nfile, fnm, "w");
 627             fprintf(fp3, "[ position_restraints ]\n"
 628                     "#define FCX 1000\n"
 629                     "#define FCY 1000\n"
 630                     "#define FCZ 1000\n"
 631                     "; Atom  Type  fx   fy   fz\n");
 632         }
 633         for (i = 0; i < nx[0]; i++)
 634         {
 635             ii  = index[0][i];
 636             res = atoms->atom[ii].resind;
 637             if (i == nx[0]-1 || res != atoms->atom[index[0][i+1]].resind)
 638             {
 639                 fluc2 = res_area2[res]-sqr(res_area[res]);
 640                 if (fluc2 < 0)
 641                 {
 642                     fluc2 = 0;
 643                 }
 644                 fprintf(fp, "%10d  %10g %10g\n",
 645                         atoms->resinfo[res].nr, res_area[res], sqrt(fluc2));
 646             }
 647             fluc2 = atom_area2[i]-sqr(atom_area[i]);
 648             if (fluc2 < 0)
 649             {
 650                 fluc2 = 0;
 651             }
 652             fprintf(fp2, "%d %g %g\n", index[0][i]+1, atom_area[i], sqrt(fluc2));
 653             if (bITP && (atom_area[i] > minarea))
 654             {
 655                 fprintf(fp3, "%5d   1     FCX  FCX  FCZ\n", ii+1);
 656             }
 657         }
 658         if (bITP)
 659         {
 660             ffclose(fp3);
 661         }
 662         ffclose(fp);
 663     }
 664
 665     /* Be a good citizen, keep our memory free! */
 666     sfree(x);
 667     sfree(nx);
 668     for (i = 0; i < 2; i++)
 669     {
 670         sfree(index[i]);
 671         sfree(grpname[i]);
 672     }
 673     sfree(bOut);
 674     sfree(radius);
 675     sfree(bPhobic);
 676
 677     if (bResAt)
 678     {
 679         sfree(atom_area);
 680         sfree(atom_area2);
 681         sfree(res_a);
 682         sfree(res_area);
 683         sfree(res_area2);
 684     }
 685     if (bDGsol)
 686     {
 687         sfree(dgs_factor);
 688     }
 689 }
 690
 691 int gmx_sas(int argc, char *argv[])
 692 {
 693     const char     *desc[] = {
 694         "[TT]g_sas[tt] computes hydrophobic, hydrophilic and total solvent",
 695         "accessible surface area. See Eisenhaber F, Lijnzaad P, Argos P,",
 696         "Sander C, & Scharf M (1995) J. Comput. Chem. 16, 273-284.",
 697         "As a side effect, the Connolly surface can be generated as well in",
 698         "a [TT].pdb[tt] file where the nodes are represented as atoms and the",
 699         "vertice connecting the nearest nodes as CONECT records.",
 700         "The program will ask for a group for the surface calculation",
 701         "and a group for the output. The calculation group should always",
 702         "consists of all the non-solvent atoms in the system.",
 703         "The output group can be the whole or part of the calculation group.",
 704         "The average and standard deviation of the area over the trajectory can be plotted",
 705         "per residue and atom as well (options [TT]-or[tt] and [TT]-oa[tt]).",
 706         "In combination with the latter option an [TT].itp[tt] file can be",
 707         "generated (option [TT]-i[tt])",
 708         "which can be used to restrain surface atoms.[PAR]",
 709
 710         "By default, periodic boundary conditions are taken into account,",
 711         "this can be turned off using the [TT]-nopbc[tt] option.[PAR]",
 712
 713         "With the [TT]-tv[tt] option the total volume and density of the",
 714         "molecule can be computed.",
 715         "Please consider whether the normal probe radius is appropriate",
 716         "in this case or whether you would rather use e.g. 0. It is good",
 717         "to keep in mind that the results for volume and density are very",
 718         "approximate. For example, in ice Ih, one can easily fit water molecules in the",
 719         "pores which would yield a volume that is too low, and surface area and density",
 720         "that are both too high."
 721     };
 722
 723     output_env_t    oenv;
 724     static real     solsize = 0.14;
 725     static int      ndots   = 24;
 726     static real     qcut    = 0.2;
 727     static real     minarea = 0.5, dgs_default = 0;
 728     static gmx_bool bSave   = TRUE, bPBC = TRUE, bFindex = FALSE;
 729     t_pargs         pa[]    = {
 730         { "-probe", FALSE, etREAL, {&solsize},
 731           "Radius of the solvent probe (nm)" },
 732         { "-ndots",   FALSE, etINT,  {&ndots},
 733           "Number of dots per sphere, more dots means more accuracy" },
 734         { "-qmax",    FALSE, etREAL, {&qcut},
 735           "The maximum charge (e, absolute value) of a hydrophobic atom" },
 736         { "-f_index", FALSE, etBOOL, {&bFindex},
 737           "Determine from a group in the index file what are the hydrophobic atoms rather than from the charge" },
 738         { "-minarea", FALSE, etREAL, {&minarea},
 739           "The minimum area (nm^2) to count an atom as a surface atom when writing a position restraint file  (see help)" },
 740         { "-pbc",     FALSE, etBOOL, {&bPBC},
 741           "Take periodicity into account" },
 742         { "-prot",    FALSE, etBOOL, {&bSave},
 743           "Output the protein to the Connelly [TT].pdb[tt] file too" },
 744         { "-dgs",     FALSE, etREAL, {&dgs_default},
 745           "Default value for solvation free energy per area (kJ/mol/nm^2)" }
 746     };
 747     t_filenm        fnm[] = {
 748         { efTRX, "-f",   NULL,       ffREAD },
 749         { efTPS, "-s",   NULL,       ffREAD },
 750         { efXVG, "-o",   "area",     ffWRITE },
 751         { efXVG, "-or",  "resarea",  ffOPTWR },
 752         { efXVG, "-oa",  "atomarea", ffOPTWR },
 753         { efXVG, "-tv",  "volume",   ffOPTWR },
 754         { efPDB, "-q",   "connelly", ffOPTWR },
 755         { efNDX, "-n",   "index",    ffOPTRD },
 756         { efITP, "-i",   "surfat",   ffOPTWR }
 757     };
 758 #define NFILE asize(fnm)
 759
 760     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 761                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 762     {
 763         return 0;
 764     }
 765     if (solsize < 0)
 766     {
 767         solsize = 1e-3;
 768         fprintf(stderr, "Probe size too small, setting it to %g\n", solsize);
 769     }
 770     if (ndots < 20)
 771     {
 772         ndots = 20;
 773         fprintf(stderr, "Ndots too small, setting it to %d\n", ndots);
 774     }
 775
 776     please_cite(stderr, "Eisenhaber95");
 777
 778     sas_plot(NFILE, fnm, solsize, ndots, qcut, bSave, minarea, bPBC, dgs_default, bFindex,
 779              oenv);
 780
 781     do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
 782     do_view(oenv, opt2fn_null("-or", NFILE, fnm), "-nxy");
 783     do_view(oenv, opt2fn_null("-oa", NFILE, fnm), "-nxy");
 784
 785     return 0;
 786 }