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