2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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"
63 #define EPS0 8.85419E-12
64 #define ELC 1.60219E-19
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*/
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 /****************************************************************************/
76 static int ce = 0, cb = 0;
78 /* this routine integrates the array data and returns the resulting array */
79 /* routine uses simple trapezoid rule */
80 static void p_integrate(double* result, const double data[], int ndata, double slWidth)
88 "Warning: nr of slices very small. This will result"
92 fprintf(stderr, "Integrating from slice %d to slice %d\n", cb, ndata - ce);
94 for (slice = cb; slice < (ndata - ce); slice++)
97 for (i = cb; i < slice; i++)
99 sum += slWidth * (data[i] + 0.5 * (data[i + 1] - data[i]));
105 static void calc_potential(const char* fn,
108 double*** slPotential,
112 const t_topology* top,
120 const gmx_output_env_t* oenv)
122 rvec* x0; /* coordinates without pbc */
123 matrix box; /* box (3x3) */
124 int natoms; /* nr. atoms in trj */
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 */
134 gmx_rmpbc_t gpbc = nullptr;
150 default: gmx_fatal(FARGS, "Invalid axes. Terminating\n");
153 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
155 gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
160 *nslices = static_cast<int>(box[axis][axis] * 10.0); /* default value */
162 fprintf(stderr, "\nDividing the box in %d slices\n", *nslices);
164 snew(*slField, nr_grps);
165 snew(*slCharge, nr_grps);
166 snew(*slPotential, nr_grps);
168 for (i = 0; i < nr_grps; i++)
170 snew((*slField)[i], *nslices);
171 snew((*slCharge)[i], *nslices);
172 snew((*slPotential)[i], *nslices);
176 gpbc = gmx_rmpbc_init(&top->idef, pbcType, natoms);
178 /*********** Start processing trajectory ***********/
181 *slWidth = box[axis][axis] / static_cast<real>((*nslices));
183 gmx_rmpbc(gpbc, natoms, box, x0);
185 /* calculate position of center of mass based on group 1 */
186 calc_xcm(x0, gnx[0], index[0], top->atoms.atom, xcm, FALSE);
189 for (n = 0; n < nr_grps; n++)
191 /* Check whether we actually have all positions of the requested index
192 * group in the trajectory file */
196 "You selected a group with %d atoms, but only %d atoms\n"
197 "were found in the trajectory.\n",
201 for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */
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));
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
212 if (slice > (*nslices))
214 fprintf(stderr,"Warning: slice = %d\n",slice);
217 (*slCharge)[n][slice] += top->atoms.atom[index[n][i]].q;
221 z = x0[index[n][i]][axis];
225 z += box[axis][axis];
227 if (z > box[axis][axis])
229 z -= box[axis][axis];
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;
238 } while (read_next_x(oenv, status, &t, x0, box));
240 gmx_rmpbc_done(gpbc);
242 /*********** done with status file **********/
245 /* slCharge now contains the total charge per slice, summed over all
246 frames. Now divide by nr_frames and integrate twice
253 "\n\nRead %d frames from trajectory. Calculating potential"
254 "in spherical coordinates\n",
259 fprintf(stderr, "\n\nRead %d frames from trajectory. Calculating potential\n", nr_frames);
262 for (n = 0; n < nr_grps; n++)
264 for (i = 0; i < *nslices; i++)
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
271 slVolume = 4 * M_PI * gmx::square(i) * gmx::square(*slWidth) * *slWidth;
274 (*slCharge)[n][i] = 0;
278 (*slCharge)[n][i] = (*slCharge)[n][i] / (nr_frames * slVolume);
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]);
289 /* Now we have charge densities */
292 if (bCorrect && !bSpherical)
294 for (n = 0; n < nr_grps; n++)
298 for (i = 0; i < *nslices; i++)
300 if (std::abs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN)
303 qsum += (*slCharge)[n][i];
307 for (i = 0; i < *nslices; i++)
309 if (std::abs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN)
311 (*slCharge)[n][i] -= qsum;
317 for (n = 0; n < nr_grps; n++)
319 /* integrate twice to get field and potential */
320 p_integrate((*slField)[n], (*slCharge)[n], *nslices, *slWidth);
324 if (bCorrect && !bSpherical)
326 for (n = 0; n < nr_grps; n++)
330 for (i = 0; i < *nslices; i++)
332 if (std::abs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN)
335 qsum += (*slField)[n][i];
339 for (i = 0; i < *nslices; i++)
341 if (std::abs((*slCharge)[n][i]) >= GMX_DOUBLE_MIN)
343 (*slField)[n][i] -= qsum;
349 for (n = 0; n < nr_grps; n++)
351 p_integrate((*slPotential)[n], (*slField)[n], *nslices, *slWidth);
354 /* Now correct for eps0 and in spherical case for r*/
355 for (n = 0; n < nr_grps; n++)
357 for (i = 0; i < *nslices; i++)
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));
366 (*slPotential)[n][i] = ELC * (*slPotential)[n][i] * -1.0E9 / EPS0;
367 (*slField)[n][i] = ELC * (*slField)[n][i] * 1E18 / EPS0;
372 sfree(x0); /* free memory used by coordinate array */
375 static void plot_potential(double* potential[],
383 const char* const grpname[],
385 const gmx_output_env_t* oenv)
387 FILE *pot, /* xvgr file with potential */
388 *cha, /* xvgr file with charges */
389 *fie; /* xvgr files with fields */
390 char buf[256]; /* for xvgr title */
393 sprintf(buf, "Electrostatic Potential");
394 pot = xvgropen(afile, buf, "Box (nm)", "Potential (V)", oenv);
395 xvgr_legend(pot, nr_grps, grpname, oenv);
397 sprintf(buf, "Charge Distribution");
398 cha = xvgropen(bfile, buf, "Box (nm)", "Charge density (q/nm\\S3\\N)", oenv);
399 xvgr_legend(cha, nr_grps, grpname, oenv);
401 sprintf(buf, "Electric Field");
402 fie = xvgropen(cfile, buf, "Box (nm)", "Field (V/nm)", oenv);
403 xvgr_legend(fie, nr_grps, grpname, oenv);
405 for (slice = cb; slice < (nslices - ce); slice++)
407 fprintf(pot, "%20.16g ", slice * slWidth);
408 fprintf(cha, "%20.16g ", slice * slWidth);
409 fprintf(fie, "%20.16g ", slice * slWidth);
410 for (n = 0; n < nr_grps; n++)
412 fprintf(pot, " %20.16g", potential[n][slice]);
413 fprintf(fie, " %20.16g", field[n][slice] / 1e9); /* convert to V/nm */
414 fprintf(cha, " %20.16g", charge[n][slice]);
426 int gmx_potential(int argc, char* argv[])
428 const char* desc[] = {
429 "[THISMODULE] computes the electrostatical potential across the box. The potential is",
430 "calculated by first summing the charges per slice and then integrating",
431 "twice of this charge distribution. Periodic boundaries are not taken",
432 "into account. Reference of potential is taken to be the left side of",
433 "the box. It is also possible to calculate the potential in spherical",
434 "coordinates as function of r by calculating a charge distribution in",
435 "spherical slices and twice integrating them. epsilon_r is taken as 1,",
436 "but 2 is more appropriate in many cases."
438 gmx_output_env_t* oenv;
439 static int axis = 2; /* normal to memb. default z */
440 static const char* axtitle = "Z";
441 static int nslices = 10; /* nr of slices defined */
442 static int ngrps = 1;
443 static gmx_bool bSpherical = FALSE; /* default is bilayer types */
444 static real fudge_z = 0; /* translate coordinates */
445 static gmx_bool bCorrect = false;
451 "Take the normal on the membrane in direction X, Y or Z." },
456 "Calculate potential as function of boxlength, dividing the box"
457 " in this number of slices." },
462 "Discard this number of first slices of box for integration" },
467 "Discard this number of last slices of box for integration" },
472 "Translate all coordinates by this distance in the direction of the box" },
473 { "-spherical", FALSE, etBOOL, { &bSpherical }, "Calculate in spherical coordinates" },
474 { "-ng", FALSE, etINT, { &ngrps }, "Number of groups to consider" },
479 "Assume net zero charge of groups to improve accuracy" }
481 const char* bugs[] = { "Discarding slices for integration should not be necessary." };
483 double **potential, /* potential per slice */
484 **charge, /* total charge per slice */
485 **field, /* field per slice */
486 slWidth; /* width of one slice */
487 char** grpname; /* groupnames */
488 int* ngx; /* sizes of groups */
489 t_topology* top; /* topology */
491 int** index; /* indices for all groups */
493 /* files for g_order */
494 { efTRX, "-f", nullptr, ffREAD }, /* trajectory file */
495 { efNDX, nullptr, nullptr, ffREAD }, /* index file */
496 { efTPR, nullptr, nullptr, ffREAD }, /* topology file */
497 { efXVG, "-o", "potential", ffWRITE }, /* xvgr output file */
498 { efXVG, "-oc", "charge", ffWRITE }, /* xvgr output file */
499 { efXVG, "-of", "field", ffWRITE }, /* xvgr output file */
502 #define NFILE asize(fnm)
504 if (!parse_common_args(
505 &argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
511 axis = toupper(axtitle[0]) - 'X';
513 top = read_top(ftp2fn(efTPR, NFILE, fnm), &pbcType); /* read topology file */
515 snew(grpname, ngrps);
519 rd_index(ftp2fn(efNDX, NFILE, fnm), ngrps, ngx, index, grpname);
522 calc_potential(ftp2fn(efTRX, NFILE, fnm),
539 plot_potential(potential,
542 opt2fn("-o", NFILE, fnm),
543 opt2fn("-oc", NFILE, fnm),
544 opt2fn("-of", NFILE, fnm),
551 do_view(oenv, opt2fn("-o", NFILE, fnm), nullptr); /* view xvgr file */
552 do_view(oenv, opt2fn("-oc", NFILE, fnm), nullptr); /* view xvgr file */
553 do_view(oenv, opt2fn("-of", NFILE, fnm), nullptr); /* view xvgr file */