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