src/gromacs/gmxana/gmx_potential.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-2004, 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
  41 #include <math.h>
  42 #include <ctype.h>
  43
  44 #include "sysstuff.h"
  45 #include <string.h>
  46 #include "typedefs.h"
  47 #include "gromacs/utility/smalloc.h"
  48 #include "macros.h"
  49 #include "princ.h"
  50 #include "rmpbc.h"
  51 #include "vec.h"
  52 #include "xvgr.h"
  53 #include "pbc.h"
  54 #include "gromacs/fileio/futil.h"
  55 #include "gromacs/commandline/pargs.h"
  56 #include "index.h"
  57 #include "gmx_ana.h"
  58 #include "gromacs/utility/cstringutil.h"
  59 #include "gromacs/fileio/trxio.h"
  60
  61 #include "gromacs/legacyheaders/gmx_fatal.h"
  62
  63 #define EPS0 8.85419E-12
  64 #define ELC 1.60219E-19
  65
  66 /****************************************************************************/
  67 /* This program calculates the electrostatic potential across the box by    */
  68 /* determining the charge density in slices of the box and integrating these*/
  69 /* twice.                                                                   */
  70 /* Peter Tieleman, April 1995                                               */
  71 /* It now also calculates electrostatic potential in spherical micelles,    */
  72 /* using \frac{1}{r}\frac{d^2r\Psi}{r^2} = - \frac{\rho}{\epsilon_0}        */
  73 /* This probably sucks but it seems to work.                                */
  74 /****************************************************************************/
  75
  76 static int ce = 0, cb = 0;
  77
  78 /* this routine integrates the array data and returns the resulting array */
  79 /* routine uses simple trapezoid rule                                     */
  80 void p_integrate(double *result, double data[], int ndata, double slWidth)
  81 {
  82     int    i, slice;
  83     double sum;
  84
  85     if (ndata <= 2)
  86     {
  87         fprintf(stderr, "Warning: nr of slices very small. This will result"
  88                 "in nonsense.\n");
  89     }
  90
  91     fprintf(stderr, "Integrating from slice %d to slice %d\n", cb, ndata-ce);
  92
  93     for (slice = cb; slice < (ndata-ce); slice++)
  94     {
  95         sum = 0;
  96         for (i = cb; i < slice; i++)
  97         {
  98             sum += slWidth * (data[i] + 0.5 * (data[i+1] - data[i]));
  99         }
 100         result[slice] = sum;
 101     }
 102     return;
 103 }
 104
 105 void calc_potential(const char *fn, atom_id **index, int gnx[],
 106                     double ***slPotential, double ***slCharge,
 107                     double ***slField, int *nslices,
 108                     t_topology *top, int ePBC,
 109                     int axis, int nr_grps, double *slWidth,
 110                     double fudge_z, gmx_bool bSpherical, gmx_bool bCorrect,
 111                     const output_env_t oenv)
 112 {
 113     rvec        *x0;      /* coordinates without pbc */
 114     matrix       box;     /* box (3x3) */
 115     int          natoms;  /* nr. atoms in trj */
 116     t_trxstatus *status;
 117     int        **slCount, /* nr. of atoms in one slice for a group */
 118                  i, j, n, /* loop indices */
 119                  teller    = 0,
 120                  ax1       = 0, ax2 = 0,
 121                  nr_frames = 0, /* number of frames */
 122                  slice;         /* current slice */
 123     double       slVolume;      /* volume of slice for spherical averaging */
 124     double       qsum, nn;
 125     real         t;
 126     double       z;
 127     rvec         xcm;
 128     gmx_rmpbc_t  gpbc = NULL;
 129
 130     switch (axis)
 131     {
 132         case 0:
 133             ax1 = 1; ax2 = 2;
 134             break;
 135         case 1:
 136             ax1 = 0; ax2 = 2;
 137             break;
 138         case 2:
 139             ax1 = 0; ax2 = 1;
 140             break;
 141         default:
 142             gmx_fatal(FARGS, "Invalid axes. Terminating\n");
 143     }
 144
 145     if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
 146     {
 147         gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
 148     }
 149
 150     if (!*nslices)
 151     {
 152         *nslices = (int)(box[axis][axis] * 10); /* default value */
 153
 154     }
 155     fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
 156
 157     snew(*slField, nr_grps);
 158     snew(*slCharge, nr_grps);
 159     snew(*slPotential, nr_grps);
 160
 161     for (i = 0; i < nr_grps; i++)
 162     {
 163         snew((*slField)[i], *nslices);
 164         snew((*slCharge)[i], *nslices);
 165         snew((*slPotential)[i], *nslices);
 166     }
 167
 168
 169     gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
 170
 171     /*********** Start processing trajectory ***********/
 172     do
 173     {
 174         *slWidth = box[axis][axis]/(*nslices);
 175         teller++;
 176         gmx_rmpbc(gpbc, natoms, box, x0);
 177
 178         /* calculate position of center of mass based on group 1 */
 179         calc_xcm(x0, gnx[0], index[0], top->atoms.atom, xcm, FALSE);
 180         svmul(-1, xcm, xcm);
 181
 182         for (n = 0; n < nr_grps; n++)
 183         {
 184             /* Check whether we actually have all positions of the requested index
 185              * group in the trajectory file */
 186             if (gnx[n] > natoms)
 187             {
 188                 gmx_fatal(FARGS, "You selected a group with %d atoms, but only %d atoms\n"
 189                           "were found in the trajectory.\n", gnx[n], natoms);
 190             }
 191             for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
 192             {
 193                 if (bSpherical)
 194                 {
 195                     rvec_add(x0[index[n][i]], xcm, x0[index[n][i]]);
 196                     /* only distance from origin counts, not sign */
 197                     slice = norm(x0[index[n][i]])/(*slWidth);
 198
 199                     /* this is a nice check for spherical groups but not for
 200                        all water in a cubic box since a lot will fall outside
 201                        the sphere
 202                        if (slice > (*nslices))
 203                        {
 204                        fprintf(stderr,"Warning: slice = %d\n",slice);
 205                        }
 206                      */
 207                     (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
 208                 }
 209                 else
 210                 {
 211                     z = x0[index[n][i]][axis];
 212                     z = z + fudge_z;
 213                     if (z < 0)
 214                     {
 215                         z += box[axis][axis];
 216                     }
 217                     if (z > box[axis][axis])
 218                     {
 219                         z -= box[axis][axis];
 220                     }
 221                     /* determine which slice atom is in */
 222                     slice                  = (z / (*slWidth));
 223                     (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
 224                 }
 225             }
 226         }
 227         nr_frames++;
 228     }
 229     while (read_next_x(oenv, status, &t, x0, box));
 230
 231     gmx_rmpbc_done(gpbc);
 232
 233     /*********** done with status file **********/
 234     close_trj(status);
 235
 236     /* slCharge now contains the total charge per slice, summed over all
 237        frames. Now divide by nr_frames and integrate twice
 238      */
 239
 240
 241     if (bSpherical)
 242     {
 243         fprintf(stderr, "\n\nRead %d frames from trajectory. Calculating potential"
 244                 "in spherical coordinates\n", nr_frames);
 245     }
 246     else
 247     {
 248         fprintf(stderr, "\n\nRead %d frames from trajectory. Calculating potential\n",
 249                 nr_frames);
 250     }
 251
 252     for (n = 0; n < nr_grps; n++)
 253     {
 254         for (i = 0; i < *nslices; i++)
 255         {
 256             if (bSpherical)
 257             {
 258                 /* charge per volume is now the summed charge, divided by the nr
 259                    of frames and by the volume of the slice it's in, 4pi r^2 dr
 260                  */
 261                 slVolume = 4*M_PI * sqr(i) * sqr(*slWidth) * *slWidth;
 262                 if (slVolume == 0)
 263                 {
 264                     (*slCharge)[n][i] = 0;
 265                 }
 266                 else
 267                 {
 268                     (*slCharge)[n][i] = (*slCharge)[n][i] / (nr_frames * slVolume);
 269                 }
 270             }
 271             else
 272             {
 273                 /* get charge per volume */
 274                 (*slCharge)[n][i] = (*slCharge)[n][i] * (*nslices) /
 275                     (nr_frames * box[axis][axis] * box[ax1][ax1] * box[ax2][ax2]);
 276             }
 277         }
 278         /* Now we have charge densities */
 279     }
 280
 281     if (bCorrect && !bSpherical)
 282     {
 283         for (n = 0; n < nr_grps; n++)
 284         {
 285             nn   = 0;
 286             qsum = 0;
 287             for (i = 0; i < *nslices; i++)
 288             {
 289                 if (fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN)
 290                 {
 291                     nn++;
 292                     qsum += (*slCharge)[n][i];
 293                 }
 294             }
 295             qsum /= nn;
 296             for (i = 0; i < *nslices; i++)
 297             {
 298                 if (fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN)
 299                 {
 300                     (*slCharge)[n][i] -= qsum;
 301                 }
 302             }
 303         }
 304     }
 305
 306     for (n = 0; n < nr_grps; n++)
 307     {
 308         /* integrate twice to get field and potential */
 309         p_integrate((*slField)[n], (*slCharge)[n], *nslices, *slWidth);
 310     }
 311
 312
 313     if (bCorrect && !bSpherical)
 314     {
 315         for (n = 0; n < nr_grps; n++)
 316         {
 317             nn   = 0;
 318             qsum = 0;
 319             for (i = 0; i < *nslices; i++)
 320             {
 321                 if (fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN)
 322                 {
 323                     nn++;
 324                     qsum += (*slField)[n][i];
 325                 }
 326             }
 327             qsum /= nn;
 328             for (i = 0; i < *nslices; i++)
 329             {
 330                 if (fabs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN)
 331                 {
 332                     (*slField)[n][i] -= qsum;
 333                 }
 334             }
 335         }
 336     }
 337
 338     for (n = 0; n < nr_grps; n++)
 339     {
 340         p_integrate((*slPotential)[n], (*slField)[n], *nslices, *slWidth);
 341     }
 342
 343     /* Now correct for eps0 and in spherical case for r*/
 344     for (n = 0; n < nr_grps; n++)
 345     {
 346         for (i = 0; i < *nslices; i++)
 347         {
 348             if (bSpherical)
 349             {
 350                 (*slPotential)[n][i] = ELC * (*slPotential)[n][i] * -1.0E9 /
 351                     (EPS0 * i * (*slWidth));
 352                 (*slField)[n][i] = ELC * (*slField)[n][i] * 1E18 /
 353                     (EPS0 * i * (*slWidth));
 354             }
 355             else
 356             {
 357                 (*slPotential)[n][i] = ELC * (*slPotential)[n][i] * -1.0E9 / EPS0;
 358                 (*slField)[n][i]     = ELC * (*slField)[n][i] * 1E18 / EPS0;
 359             }
 360         }
 361     }
 362
 363     sfree(x0); /* free memory used by coordinate array */
 364 }
 365
 366 void plot_potential(double *potential[], double *charge[], double *field[],
 367                     const char *afile, const char *bfile, const char *cfile,
 368                     int nslices, int nr_grps, const char *grpname[], double slWidth,
 369                     const output_env_t oenv)
 370 {
 371     FILE       *pot,     /* xvgr file with potential */
 372     *cha,                /* xvgr file with charges   */
 373     *fie;                /* xvgr files with fields   */
 374     char       buf[256]; /* for xvgr title */
 375     int        slice, n;
 376
 377     sprintf(buf, "Electrostatic Potential");
 378     pot = xvgropen(afile, buf, "Box (nm)", "Potential (V)", oenv);
 379     xvgr_legend(pot, nr_grps, grpname, oenv);
 380
 381     sprintf(buf, "Charge Distribution");
 382     cha = xvgropen(bfile, buf, "Box (nm)", "Charge density (q/nm\\S3\\N)", oenv);
 383     xvgr_legend(cha, nr_grps, grpname, oenv);
 384
 385     sprintf(buf, "Electric Field");
 386     fie = xvgropen(cfile, buf, "Box (nm)", "Field (V/nm)", oenv);
 387     xvgr_legend(fie, nr_grps, grpname, oenv);
 388
 389     for (slice = cb; slice < (nslices - ce); slice++)
 390     {
 391         fprintf(pot, "%20.16g  ", slice * slWidth);
 392         fprintf(cha, "%20.16g  ", slice * slWidth);
 393         fprintf(fie, "%20.16g  ", slice * slWidth);
 394         for (n = 0; n < nr_grps; n++)
 395         {
 396             fprintf(pot, "   %20.16g", potential[n][slice]);
 397             fprintf(fie, "   %20.16g", field[n][slice]/1e9); /* convert to V/nm */
 398             fprintf(cha, "   %20.16g", charge[n][slice]);
 399         }
 400         fprintf(pot, "\n");
 401         fprintf(cha, "\n");
 402         fprintf(fie, "\n");
 403     }
 404
 405     gmx_ffclose(pot);
 406     gmx_ffclose(cha);
 407     gmx_ffclose(fie);
 408 }
 409
 410 int gmx_potential(int argc, char *argv[])
 411 {
 412     const char        *desc[] = {
 413         "[THISMODULE] computes the electrostatical potential across the box. The potential is",
 414         "calculated by first summing the charges per slice and then integrating",
 415         "twice of this charge distribution. Periodic boundaries are not taken",
 416         "into account. Reference of potential is taken to be the left side of",
 417         "the box. It is also possible to calculate the potential in spherical",
 418         "coordinates as function of r by calculating a charge distribution in",
 419         "spherical slices and twice integrating them. epsilon_r is taken as 1,",
 420         "but 2 is more appropriate in many cases."
 421     };
 422     output_env_t       oenv;
 423     static int         axis       = 2;       /* normal to memb. default z  */
 424     static const char *axtitle    = "Z";
 425     static int         nslices    = 10;      /* nr of slices defined       */
 426     static int         ngrps      = 1;
 427     static gmx_bool    bSpherical = FALSE;   /* default is bilayer types   */
 428     static real        fudge_z    = 0;       /* translate coordinates      */
 429     static gmx_bool    bCorrect   = 0;
 430     t_pargs            pa []      = {
 431         { "-d",   FALSE, etSTR, {&axtitle},
 432           "Take the normal on the membrane in direction X, Y or Z." },
 433         { "-sl",  FALSE, etINT, {&nslices},
 434           "Calculate potential as function of boxlength, dividing the box"
 435           " in this number of slices." },
 436         { "-cb",  FALSE, etINT, {&cb},
 437           "Discard this number of  first slices of box for integration" },
 438         { "-ce",  FALSE, etINT, {&ce},
 439           "Discard this number of last slices of box for integration" },
 440         { "-tz",  FALSE, etREAL, {&fudge_z},
 441           "Translate all coordinates by this distance in the direction of the box" },
 442         { "-spherical", FALSE, etBOOL, {&bSpherical},
 443           "Calculate spherical thingie" },
 444         { "-ng",       FALSE, etINT, {&ngrps},
 445           "Number of groups to consider" },
 446         { "-correct",  FALSE, etBOOL, {&bCorrect},
 447           "Assume net zero charge of groups to improve accuracy" }
 448     };
 449     const char        *bugs[] = {
 450         "Discarding slices for integration should not be necessary."
 451     };
 452
 453     double           **potential,              /* potential per slice        */
 454     **charge,                                  /* total charge per slice     */
 455     **field,                                   /* field per slice            */
 456                        slWidth;                /* width of one slice         */
 457     char      **grpname;                       /* groupnames                 */
 458     int        *ngx;                           /* sizes of groups            */
 459     t_topology *top;                           /* topology        */
 460     int         ePBC;
 461     atom_id   **index;                         /* indices for all groups     */
 462     t_filenm    fnm[] = {                      /* files for g_order       */
 463         { efTRX, "-f", NULL,  ffREAD },        /* trajectory file             */
 464         { efNDX, NULL, NULL,  ffREAD },        /* index file          */
 465         { efTPX, NULL, NULL,  ffREAD },        /* topology file               */
 466         { efXVG, "-o", "potential", ffWRITE }, /* xvgr output file    */
 467         { efXVG, "-oc", "charge", ffWRITE },   /* xvgr output file    */
 468         { efXVG, "-of", "field", ffWRITE },    /* xvgr output file    */
 469     };
 470
 471 #define NFILE asize(fnm)
 472
 473     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 474                            NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
 475                            &oenv))
 476     {
 477         return 0;
 478     }
 479
 480     /* Calculate axis */
 481     axis = toupper(axtitle[0]) - 'X';
 482
 483     top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
 484
 485     snew(grpname, ngrps);
 486     snew(index, ngrps);
 487     snew(ngx, ngrps);
 488
 489     rd_index(ftp2fn(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
 490
 491
 492     calc_potential(ftp2fn(efTRX, NFILE, fnm), index, ngx,
 493                    &potential, &charge, &field,
 494                    &nslices, top, ePBC, axis, ngrps, &slWidth, fudge_z,
 495                    bSpherical, bCorrect, oenv);
 496
 497     plot_potential(potential, charge, field, opt2fn("-o", NFILE, fnm),
 498                    opt2fn("-oc", NFILE, fnm), opt2fn("-of", NFILE, fnm),
 499                    nslices, ngrps, (const char**)grpname, slWidth, oenv);
 500
 501     do_view(oenv, opt2fn("-o", NFILE, fnm), NULL);  /* view xvgr file */
 502     do_view(oenv, opt2fn("-oc", NFILE, fnm), NULL); /* view xvgr file */
 503     do_view(oenv, opt2fn("-of", NFILE, fnm), NULL); /* view xvgr file */
 504
 505     return 0;
 506 }