src/gromacs/gmxana/gmx_spatial.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2007,2008,2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39
  40
  41 #include "gromacs/commandline/pargs.h"
  42 #include "typedefs.h"
  43 #include "gromacs/utility/smalloc.h"
  44 #include "vec.h"
  45 #include "gromacs/fileio/tpxio.h"
  46 #include "gromacs/fileio/trxio.h"
  47 #include <math.h>
  48 #include "index.h"
  49 #include "pbc.h"
  50 #include "rmpbc.h"
  51 #include "gmx_ana.h"
  52 #include "macros.h"
  53
  54
  55 static const double bohr = 0.529177249;  /* conversion factor to compensate for VMD plugin conversion... */
  56
  57 static void mequit(void)
  58 {
  59     printf("Memory allocation error\n");
  60     exit(1);
  61 }
  62
  63 int gmx_spatial(int argc, char *argv[])
  64 {
  65     const char     *desc[] = {
  66         "[THISMODULE] calculates the spatial distribution function and ",
  67         "outputs it in a form that can be read by VMD as Gaussian98 cube format. ",
  68         "For a system of 32,000 atoms and a 50 ns trajectory, the SDF can be generated ",
  69         "in about 30 minutes, with most of the time dedicated to the two runs through ",
  70         "[TT]trjconv[tt] that are required to center everything properly. ",
  71         "This also takes a whole bunch of space (3 copies of the trajectory file). ",
  72         "Still, the pictures are pretty and very informative when the fitted selection is properly made. ",
  73         "3-4 atoms in a widely mobile group (like a free amino acid in solution) works ",
  74         "well, or select the protein backbone in a stable folded structure to get the SDF ",
  75         "of solvent and look at the time-averaged solvation shell. ",
  76         "It is also possible using this program to generate the SDF based on some arbitrary ",
  77         "Cartesian coordinate. To do that, simply omit the preliminary [gmx-trjconv] steps. \n",
  78         "USAGE: \n",
  79         "1. Use [gmx-make_ndx] to create a group containing the atoms around which you want the SDF \n",
  80         "2. [TT]gmx trjconv -s a.tpr -f a.tng -o b.tng -boxcenter tric -ur compact -pbc none[tt] \n",
  81         "3. [TT]gmx trjconv -s a.tpr -f b.tng -o c.tng -fit rot+trans[tt] \n",
  82         "4. run [THISMODULE] on the [TT]c.tng[tt] output of step #3. \n",
  83         "5. Load [TT]grid.cube[tt] into VMD and view as an isosurface. \n",
  84         "[BB]Note[bb] that systems such as micelles will require [TT]gmx trjconv -pbc cluster[tt] between steps 1 and 2\n",
  85         "WARNINGS:[BR]",
  86         "The SDF will be generated for a cube that contains all bins that have some non-zero occupancy. ",
  87         "However, the preparatory [TT]-fit rot+trans[tt] option to [gmx-trjconv] implies that your system will be rotating ",
  88         "and translating in space (in order that the selected group does not). Therefore the values that are ",
  89         "returned will only be valid for some region around your central group/coordinate that has full overlap ",
  90         "with system volume throughout the entire translated/rotated system over the course of the trajectory. ",
  91         "It is up to the user to ensure that this is the case. \n",
  92         "BUGS:[BR]",
  93         "When the allocated memory is not large enough, a segmentation fault may occur. This is usually detected ",
  94         "and the program is halted prior to the fault while displaying a warning message suggesting the use of the [TT]-nab[tt] (Number of Additional Bins)",
  95         "option. However, the program does not detect all such events. If you encounter a segmentation fault, run it again ",
  96         "with an increased [TT]-nab[tt] value. \n",
  97         "RISKY OPTIONS:[BR]",
  98         "To reduce the amount of space and time required, you can output only the coords ",
  99         "that are going to be used in the first and subsequent run through [gmx-trjconv]. ",
 100         "However, be sure to set the [TT]-nab[tt] option to a sufficiently high value since ",
 101         "memory is allocated for cube bins based on the initial coordinates and the [TT]-nab[tt] ",
 102         "option value. \n"
 103     };
 104
 105     static gmx_bool bPBC         = FALSE;
 106     static gmx_bool bSHIFT       = FALSE;
 107     static int      iIGNOREOUTER = -1;   /*Positive values may help if the surface is spikey */
 108     static gmx_bool bCUTDOWN     = TRUE;
 109     static real     rBINWIDTH    = 0.05; /* nm */
 110     static gmx_bool bCALCDIV     = TRUE;
 111     static int      iNAB         = 4;
 112
 113     t_pargs         pa[] = {
 114         { "-pbc",      FALSE, etBOOL, {&bPBC},
 115           "Use periodic boundary conditions for computing distances" },
 116         { "-div",      FALSE, etBOOL, {&bCALCDIV},
 117           "Calculate and apply the divisor for bin occupancies based on atoms/minimal cube size. Set as TRUE for visualization and as FALSE ([TT]-nodiv[tt]) to get accurate counts per frame" },
 118         { "-ign",      FALSE, etINT, {&iIGNOREOUTER},
 119           "Do not display this number of outer cubes (positive values may reduce boundary speckles; -1 ensures outer surface is visible)" },
 120         /*    { "-cut",      bCUTDOWN, etBOOL, {&bCUTDOWN},*/
 121         /*      "Display a total cube that is of minimal size" }, */
 122         { "-bin",      FALSE, etREAL, {&rBINWIDTH},
 123           "Width of the bins (nm)" },
 124         { "-nab",      FALSE, etINT, {&iNAB},
 125           "Number of additional bins to ensure proper memory allocation" }
 126     };
 127
 128     double          MINBIN[3];
 129     double          MAXBIN[3];
 130     t_topology      top;
 131     int             ePBC;
 132     char            title[STRLEN];
 133     t_trxframe      fr;
 134     rvec           *xtop, *shx[26];
 135     matrix          box, box_pbc;
 136     t_trxstatus    *status;
 137     int             flags = TRX_READ_X;
 138     t_pbc           pbc;
 139     t_atoms        *atoms;
 140     int             natoms;
 141     char           *grpnm, *grpnmp;
 142     atom_id        *index, *indexp;
 143     int             i, nidx, nidxp;
 144     int             v;
 145     int             j, k;
 146     long         ***bin = (long ***)NULL;
 147     long            nbin[3];
 148     FILE           *flp;
 149     long            x, y, z, minx, miny, minz, maxx, maxy, maxz;
 150     long            numfr, numcu;
 151     long            tot, max, min;
 152     double          norm;
 153     output_env_t    oenv;
 154     gmx_rmpbc_t     gpbc = NULL;
 155
 156     t_filenm        fnm[] = {
 157         { efTPS,  NULL,  NULL, ffREAD }, /* this is for the topology */
 158         { efTRX, "-f", NULL, ffREAD },   /* and this for the trajectory */
 159         { efNDX, NULL, NULL, ffOPTRD }
 160     };
 161
 162 #define NFILE asize(fnm)
 163
 164     /* This is the routine responsible for adding default options,
 165      * calling the X/motif interface, etc. */
 166     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_CAN_VIEW,
 167                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
 168     {
 169         return 0;
 170     }
 171
 172     read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &xtop, NULL, box, TRUE);
 173     sfree(xtop);
 174
 175     atoms = &(top.atoms);
 176     printf("Select group to generate SDF:\n");
 177     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nidx, &index, &grpnm);
 178     printf("Select group to output coords (e.g. solute):\n");
 179     get_index(atoms, ftp2fn_null(efNDX, NFILE, fnm), 1, &nidxp, &indexp, &grpnmp);
 180
 181     /* The first time we read data is a little special */
 182     natoms = read_first_frame(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &fr, flags);
 183
 184     /* Memory Allocation */
 185     MINBIN[XX] = MAXBIN[XX] = fr.x[0][XX];
 186     MINBIN[YY] = MAXBIN[YY] = fr.x[0][YY];
 187     MINBIN[ZZ] = MAXBIN[ZZ] = fr.x[0][ZZ];
 188     for (i = 1; i < top.atoms.nr; ++i)
 189     {
 190         if (fr.x[i][XX] < MINBIN[XX])
 191         {
 192             MINBIN[XX] = fr.x[i][XX];
 193         }
 194         if (fr.x[i][XX] > MAXBIN[XX])
 195         {
 196             MAXBIN[XX] = fr.x[i][XX];
 197         }
 198         if (fr.x[i][YY] < MINBIN[YY])
 199         {
 200             MINBIN[YY] = fr.x[i][YY];
 201         }
 202         if (fr.x[i][YY] > MAXBIN[YY])
 203         {
 204             MAXBIN[YY] = fr.x[i][YY];
 205         }
 206         if (fr.x[i][ZZ] < MINBIN[ZZ])
 207         {
 208             MINBIN[ZZ] = fr.x[i][ZZ];
 209         }
 210         if (fr.x[i][ZZ] > MAXBIN[ZZ])
 211         {
 212             MAXBIN[ZZ] = fr.x[i][ZZ];
 213         }
 214     }
 215     for (i = ZZ; i >= XX; --i)
 216     {
 217         MAXBIN[i]  = (ceil((MAXBIN[i]-MINBIN[i])/rBINWIDTH)+(double)iNAB)*rBINWIDTH+MINBIN[i];
 218         MINBIN[i] -= (double)iNAB*rBINWIDTH;
 219         nbin[i]    = (long)ceil((MAXBIN[i]-MINBIN[i])/rBINWIDTH);
 220     }
 221     bin = (long ***)malloc(nbin[XX]*sizeof(long **));
 222     if (!bin)
 223     {
 224         mequit();
 225     }
 226     for (i = 0; i < nbin[XX]; ++i)
 227     {
 228         bin[i] = (long **)malloc(nbin[YY]*sizeof(long *));
 229         if (!bin[i])
 230         {
 231             mequit();
 232         }
 233         for (j = 0; j < nbin[YY]; ++j)
 234         {
 235             bin[i][j] = (long *)calloc(nbin[ZZ], sizeof(long));
 236             if (!bin[i][j])
 237             {
 238                 mequit();
 239             }
 240         }
 241     }
 242     copy_mat(box, box_pbc);
 243     numfr = 0;
 244     minx  = miny = minz = 999;
 245     maxx  = maxy = maxz = 0;
 246
 247     if (bPBC)
 248     {
 249         gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
 250     }
 251     /* This is the main loop over frames */
 252     do
 253     {
 254         /* Must init pbc every step because of pressure coupling */
 255
 256         copy_mat(box, box_pbc);
 257         if (bPBC)
 258         {
 259             gmx_rmpbc_trxfr(gpbc, &fr);
 260             set_pbc(&pbc, ePBC, box_pbc);
 261         }
 262
 263         for (i = 0; i < nidx; i++)
 264         {
 265             if (fr.x[index[i]][XX] < MINBIN[XX] || fr.x[index[i]][XX] > MAXBIN[XX] ||
 266                 fr.x[index[i]][YY] < MINBIN[YY] || fr.x[index[i]][YY] > MAXBIN[YY] ||
 267                 fr.x[index[i]][ZZ] < MINBIN[ZZ] || fr.x[index[i]][ZZ] > MAXBIN[ZZ])
 268             {
 269                 printf("There was an item outside of the allocated memory. Increase the value given with the -nab option.\n");
 270                 printf("Memory was allocated for [%f,%f,%f]\tto\t[%f,%f,%f]\n", MINBIN[XX], MINBIN[YY], MINBIN[ZZ], MAXBIN[XX], MAXBIN[YY], MAXBIN[ZZ]);
 271                 printf("Memory was required for [%f,%f,%f]\n", fr.x[index[i]][XX], fr.x[index[i]][YY], fr.x[index[i]][ZZ]);
 272                 exit(1);
 273             }
 274             x = (long)ceil((fr.x[index[i]][XX]-MINBIN[XX])/rBINWIDTH);
 275             y = (long)ceil((fr.x[index[i]][YY]-MINBIN[YY])/rBINWIDTH);
 276             z = (long)ceil((fr.x[index[i]][ZZ]-MINBIN[ZZ])/rBINWIDTH);
 277             ++bin[x][y][z];
 278             if (x < minx)
 279             {
 280                 minx = x;
 281             }
 282             if (x > maxx)
 283             {
 284                 maxx = x;
 285             }
 286             if (y < miny)
 287             {
 288                 miny = y;
 289             }
 290             if (y > maxy)
 291             {
 292                 maxy = y;
 293             }
 294             if (z < minz)
 295             {
 296                 minz = z;
 297             }
 298             if (z > maxz)
 299             {
 300                 maxz = z;
 301             }
 302         }
 303         numfr++;
 304         /* printf("%f\t%f\t%f\n",box[XX][XX],box[YY][YY],box[ZZ][ZZ]); */
 305
 306     }
 307     while (read_next_frame(oenv, status, &fr));
 308
 309     if (bPBC)
 310     {
 311         gmx_rmpbc_done(gpbc);
 312     }
 313
 314     if (!bCUTDOWN)
 315     {
 316         minx = miny = minz = 0;
 317         maxx = nbin[XX];
 318         maxy = nbin[YY];
 319         maxz = nbin[ZZ];
 320     }
 321
 322     /* OUTPUT */
 323     flp = gmx_ffopen("grid.cube", "w");
 324     fprintf(flp, "Spatial Distribution Function\n");
 325     fprintf(flp, "test\n");
 326     fprintf(flp, "%5d%12.6f%12.6f%12.6f\n", nidxp, (MINBIN[XX]+(minx+iIGNOREOUTER)*rBINWIDTH)*10./bohr, (MINBIN[YY]+(miny+iIGNOREOUTER)*rBINWIDTH)*10./bohr, (MINBIN[ZZ]+(minz+iIGNOREOUTER)*rBINWIDTH)*10./bohr);
 327     fprintf(flp, "%5ld%12.6f%12.6f%12.6f\n", maxx-minx+1-(2*iIGNOREOUTER), rBINWIDTH*10./bohr, 0., 0.);
 328     fprintf(flp, "%5ld%12.6f%12.6f%12.6f\n", maxy-miny+1-(2*iIGNOREOUTER), 0., rBINWIDTH*10./bohr, 0.);
 329     fprintf(flp, "%5ld%12.6f%12.6f%12.6f\n", maxz-minz+1-(2*iIGNOREOUTER), 0., 0., rBINWIDTH*10./bohr);
 330     for (i = 0; i < nidxp; i++)
 331     {
 332         v = 2;
 333         if (*(top.atoms.atomname[indexp[i]][0]) == 'C')
 334         {
 335             v = 6;
 336         }
 337         if (*(top.atoms.atomname[indexp[i]][0]) == 'N')
 338         {
 339             v = 7;
 340         }
 341         if (*(top.atoms.atomname[indexp[i]][0]) == 'O')
 342         {
 343             v = 8;
 344         }
 345         if (*(top.atoms.atomname[indexp[i]][0]) == 'H')
 346         {
 347             v = 1;
 348         }
 349         if (*(top.atoms.atomname[indexp[i]][0]) == 'S')
 350         {
 351             v = 16;
 352         }
 353         fprintf(flp, "%5d%12.6f%12.6f%12.6f%12.6f\n", v, 0., (double)fr.x[indexp[i]][XX]*10./bohr, (double)fr.x[indexp[i]][YY]*10./bohr, (double)fr.x[indexp[i]][ZZ]*10./bohr);
 354     }
 355
 356     tot = 0;
 357     for (k = 0; k < nbin[XX]; k++)
 358     {
 359         if (!(k < minx || k > maxx))
 360         {
 361             continue;
 362         }
 363         for (j = 0; j < nbin[YY]; j++)
 364         {
 365             if (!(j < miny || j > maxy))
 366             {
 367                 continue;
 368             }
 369             for (i = 0; i < nbin[ZZ]; i++)
 370             {
 371                 if (!(i < minz || i > maxz))
 372                 {
 373                     continue;
 374                 }
 375                 if (bin[k][j][i] != 0)
 376                 {
 377                     printf("A bin was not empty when it should have been empty. Programming error.\n");
 378                     printf("bin[%d][%d][%d] was = %ld\n", k, j, i, bin[k][j][i]);
 379                     exit(1);
 380                 }
 381             }
 382         }
 383     }
 384
 385     min = 999;
 386     max = 0;
 387     for (k = 0; k < nbin[XX]; k++)
 388     {
 389         if (k < minx+iIGNOREOUTER || k > maxx-iIGNOREOUTER)
 390         {
 391             continue;
 392         }
 393         for (j = 0; j < nbin[YY]; j++)
 394         {
 395             if (j < miny+iIGNOREOUTER || j > maxy-iIGNOREOUTER)
 396             {
 397                 continue;
 398             }
 399             for (i = 0; i < nbin[ZZ]; i++)
 400             {
 401                 if (i < minz+iIGNOREOUTER || i > maxz-iIGNOREOUTER)
 402                 {
 403                     continue;
 404                 }
 405                 tot += bin[k][j][i];
 406                 if (bin[k][j][i] > max)
 407                 {
 408                     max = bin[k][j][i];
 409                 }
 410                 if (bin[k][j][i] < min)
 411                 {
 412                     min = bin[k][j][i];
 413                 }
 414             }
 415         }
 416     }
 417
 418     numcu = (maxx-minx+1-(2*iIGNOREOUTER))*(maxy-miny+1-(2*iIGNOREOUTER))*(maxz-minz+1-(2*iIGNOREOUTER));
 419     if (bCALCDIV)
 420     {
 421         norm = ((double)numcu*(double)numfr) / (double)tot;
 422     }
 423     else
 424     {
 425         norm = 1.0;
 426     }
 427
 428     for (k = 0; k < nbin[XX]; k++)
 429     {
 430         if (k < minx+iIGNOREOUTER || k > maxx-iIGNOREOUTER)
 431         {
 432             continue;
 433         }
 434         for (j = 0; j < nbin[YY]; j++)
 435         {
 436             if (j < miny+iIGNOREOUTER || j > maxy-iIGNOREOUTER)
 437             {
 438                 continue;
 439             }
 440             for (i = 0; i < nbin[ZZ]; i++)
 441             {
 442                 if (i < minz+iIGNOREOUTER || i > maxz-iIGNOREOUTER)
 443                 {
 444                     continue;
 445                 }
 446                 fprintf(flp, "%12.6f ", norm*(double)bin[k][j][i]/(double)numfr);
 447             }
 448             fprintf(flp, "\n");
 449         }
 450         fprintf(flp, "\n");
 451     }
 452     gmx_ffclose(flp);
 453
 454     /* printf("x=%d to %d\n",minx,maxx); */
 455     /* printf("y=%d to %d\n",miny,maxy); */
 456     /* printf("z=%d to %d\n",minz,maxz); */
 457
 458     if (bCALCDIV)
 459     {
 460         printf("Counts per frame in all %ld cubes divided by %le\n", numcu, 1.0/norm);
 461         printf("Normalized data: average %le, min %le, max %le\n", 1.0, norm*(double)min/(double)numfr, norm*(double)max/(double)numfr);
 462     }
 463     else
 464     {
 465         printf("grid.cube contains counts per frame in all %ld cubes\n", numcu);
 466         printf("Raw data: average %le, min %le, max %le\n", 1.0/norm, (double)min/(double)numfr, (double)max/(double)numfr);
 467     }
 468
 469     return 0;
 470 }