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