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