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