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