src/gromacs/gmxana/gmx_potential.cpp

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