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