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