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, 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.
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.
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.
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.
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.
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.
46 #include "gromacs/utility/smalloc.h"
49 #include "gromacs/pbcutil/rmpbc.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/fileio/xvgr.h"
53 #include "gromacs/utility/futil.h"
54 #include "gromacs/commandline/pargs.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/fileio/trxio.h"
60 #include "gromacs/utility/fatalerror.h"
62 #define EPS0 8.85419E-12
63 #define ELC 1.60219E-19
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*/
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 /****************************************************************************/
75 static int ce
= 0, cb
= 0;
77 /* this routine integrates the array data and returns the resulting array */
78 /* routine uses simple trapezoid rule */
79 void p_integrate(double *result
, double data
[], int ndata
, double slWidth
)
86 fprintf(stderr
, "Warning: nr of slices very small. This will result"
90 fprintf(stderr
, "Integrating from slice %d to slice %d\n", cb
, ndata
-ce
);
92 for (slice
= cb
; slice
< (ndata
-ce
); slice
++)
95 for (i
= cb
; i
< slice
; i
++)
97 sum
+= slWidth
* (data
[i
] + 0.5 * (data
[i
+1] - data
[i
]));
104 void calc_potential(const char *fn
, atom_id
**index
, int gnx
[],
105 double ***slPotential
, double ***slCharge
,
106 double ***slField
, int *nslices
,
107 t_topology
*top
, int ePBC
,
108 int axis
, int nr_grps
, double *slWidth
,
109 double fudge_z
, gmx_bool bSpherical
, gmx_bool bCorrect
,
110 const output_env_t oenv
)
112 rvec
*x0
; /* coordinates without pbc */
113 matrix box
; /* box (3x3) */
114 int natoms
; /* nr. atoms in trj */
116 int **slCount
, /* nr. of atoms in one slice for a group */
117 i
, j
, n
, /* loop indices */
120 nr_frames
= 0, /* number of frames */
121 slice
; /* current slice */
122 double slVolume
; /* volume of slice for spherical averaging */
127 gmx_rmpbc_t gpbc
= NULL
;
141 gmx_fatal(FARGS
, "Invalid axes. Terminating\n");
144 if ((natoms
= read_first_x(oenv
, &status
, fn
, &t
, &x0
, box
)) == 0)
146 gmx_fatal(FARGS
, "Could not read coordinates from statusfile\n");
151 *nslices
= (int)(box
[axis
][axis
] * 10); /* default value */
154 fprintf(stderr
, "\nDividing the box in %d slices\n", *nslices
);
156 snew(*slField
, nr_grps
);
157 snew(*slCharge
, nr_grps
);
158 snew(*slPotential
, nr_grps
);
160 for (i
= 0; i
< nr_grps
; i
++)
162 snew((*slField
)[i
], *nslices
);
163 snew((*slCharge
)[i
], *nslices
);
164 snew((*slPotential
)[i
], *nslices
);
168 gpbc
= gmx_rmpbc_init(&top
->idef
, ePBC
, natoms
);
170 /*********** Start processing trajectory ***********/
173 *slWidth
= box
[axis
][axis
]/(*nslices
);
175 gmx_rmpbc(gpbc
, natoms
, box
, x0
);
177 /* calculate position of center of mass based on group 1 */
178 calc_xcm(x0
, gnx
[0], index
[0], top
->atoms
.atom
, xcm
, FALSE
);
181 for (n
= 0; n
< nr_grps
; n
++)
183 /* Check whether we actually have all positions of the requested index
184 * group in the trajectory file */
187 gmx_fatal(FARGS
, "You selected a group with %d atoms, but only %d atoms\n"
188 "were found in the trajectory.\n", gnx
[n
], natoms
);
190 for (i
= 0; i
< gnx
[n
]; i
++) /* loop over all atoms in index file */
194 rvec_add(x0
[index
[n
][i
]], xcm
, x0
[index
[n
][i
]]);
195 /* only distance from origin counts, not sign */
196 slice
= norm(x0
[index
[n
][i
]])/(*slWidth
);
198 /* this is a nice check for spherical groups but not for
199 all water in a cubic box since a lot will fall outside
201 if (slice > (*nslices))
203 fprintf(stderr,"Warning: slice = %d\n",slice);
206 (*slCharge
)[n
][slice
] += top
->atoms
.atom
[index
[n
][i
]].q
;
210 z
= x0
[index
[n
][i
]][axis
];
214 z
+= box
[axis
][axis
];
216 if (z
> box
[axis
][axis
])
218 z
-= box
[axis
][axis
];
220 /* determine which slice atom is in */
221 slice
= (z
/ (*slWidth
));
222 (*slCharge
)[n
][slice
] += top
->atoms
.atom
[index
[n
][i
]].q
;
228 while (read_next_x(oenv
, status
, &t
, x0
, box
));
230 gmx_rmpbc_done(gpbc
);
232 /*********** done with status file **********/
235 /* slCharge now contains the total charge per slice, summed over all
236 frames. Now divide by nr_frames and integrate twice
242 fprintf(stderr
, "\n\nRead %d frames from trajectory. Calculating potential"
243 "in spherical coordinates\n", nr_frames
);
247 fprintf(stderr
, "\n\nRead %d frames from trajectory. Calculating potential\n",
251 for (n
= 0; n
< nr_grps
; n
++)
253 for (i
= 0; i
< *nslices
; i
++)
257 /* charge per volume is now the summed charge, divided by the nr
258 of frames and by the volume of the slice it's in, 4pi r^2 dr
260 slVolume
= 4*M_PI
* sqr(i
) * sqr(*slWidth
) * *slWidth
;
263 (*slCharge
)[n
][i
] = 0;
267 (*slCharge
)[n
][i
] = (*slCharge
)[n
][i
] / (nr_frames
* slVolume
);
272 /* get charge per volume */
273 (*slCharge
)[n
][i
] = (*slCharge
)[n
][i
] * (*nslices
) /
274 (nr_frames
* box
[axis
][axis
] * box
[ax1
][ax1
] * box
[ax2
][ax2
]);
277 /* Now we have charge densities */
280 if (bCorrect
&& !bSpherical
)
282 for (n
= 0; n
< nr_grps
; n
++)
286 for (i
= 0; i
< *nslices
; i
++)
288 if (fabs((*slCharge
)[n
][i
]) >= GMX_DOUBLE_MIN
)
291 qsum
+= (*slCharge
)[n
][i
];
295 for (i
= 0; i
< *nslices
; i
++)
297 if (fabs((*slCharge
)[n
][i
]) >= GMX_DOUBLE_MIN
)
299 (*slCharge
)[n
][i
] -= qsum
;
305 for (n
= 0; n
< nr_grps
; n
++)
307 /* integrate twice to get field and potential */
308 p_integrate((*slField
)[n
], (*slCharge
)[n
], *nslices
, *slWidth
);
312 if (bCorrect
&& !bSpherical
)
314 for (n
= 0; n
< nr_grps
; n
++)
318 for (i
= 0; i
< *nslices
; i
++)
320 if (fabs((*slCharge
)[n
][i
]) >= GMX_DOUBLE_MIN
)
323 qsum
+= (*slField
)[n
][i
];
327 for (i
= 0; i
< *nslices
; i
++)
329 if (fabs((*slCharge
)[n
][i
]) >= GMX_DOUBLE_MIN
)
331 (*slField
)[n
][i
] -= qsum
;
337 for (n
= 0; n
< nr_grps
; n
++)
339 p_integrate((*slPotential
)[n
], (*slField
)[n
], *nslices
, *slWidth
);
342 /* Now correct for eps0 and in spherical case for r*/
343 for (n
= 0; n
< nr_grps
; n
++)
345 for (i
= 0; i
< *nslices
; i
++)
349 (*slPotential
)[n
][i
] = ELC
* (*slPotential
)[n
][i
] * -1.0E9
/
350 (EPS0
* i
* (*slWidth
));
351 (*slField
)[n
][i
] = ELC
* (*slField
)[n
][i
] * 1E18
/
352 (EPS0
* i
* (*slWidth
));
356 (*slPotential
)[n
][i
] = ELC
* (*slPotential
)[n
][i
] * -1.0E9
/ EPS0
;
357 (*slField
)[n
][i
] = ELC
* (*slField
)[n
][i
] * 1E18
/ EPS0
;
362 sfree(x0
); /* free memory used by coordinate array */
365 void plot_potential(double *potential
[], double *charge
[], double *field
[],
366 const char *afile
, const char *bfile
, const char *cfile
,
367 int nslices
, int nr_grps
, const char *grpname
[], double slWidth
,
368 const output_env_t oenv
)
370 FILE *pot
, /* xvgr file with potential */
371 *cha
, /* xvgr file with charges */
372 *fie
; /* xvgr files with fields */
373 char buf
[256]; /* for xvgr title */
376 sprintf(buf
, "Electrostatic Potential");
377 pot
= xvgropen(afile
, buf
, "Box (nm)", "Potential (V)", oenv
);
378 xvgr_legend(pot
, nr_grps
, grpname
, oenv
);
380 sprintf(buf
, "Charge Distribution");
381 cha
= xvgropen(bfile
, buf
, "Box (nm)", "Charge density (q/nm\\S3\\N)", oenv
);
382 xvgr_legend(cha
, nr_grps
, grpname
, oenv
);
384 sprintf(buf
, "Electric Field");
385 fie
= xvgropen(cfile
, buf
, "Box (nm)", "Field (V/nm)", oenv
);
386 xvgr_legend(fie
, nr_grps
, grpname
, oenv
);
388 for (slice
= cb
; slice
< (nslices
- ce
); slice
++)
390 fprintf(pot
, "%20.16g ", slice
* slWidth
);
391 fprintf(cha
, "%20.16g ", slice
* slWidth
);
392 fprintf(fie
, "%20.16g ", slice
* slWidth
);
393 for (n
= 0; n
< nr_grps
; n
++)
395 fprintf(pot
, " %20.16g", potential
[n
][slice
]);
396 fprintf(fie
, " %20.16g", field
[n
][slice
]/1e9
); /* convert to V/nm */
397 fprintf(cha
, " %20.16g", charge
[n
][slice
]);
409 int gmx_potential(int argc
, char *argv
[])
411 const char *desc
[] = {
412 "[THISMODULE] computes the electrostatical potential across the box. The potential is",
413 "calculated by first summing the charges per slice and then integrating",
414 "twice of this charge distribution. Periodic boundaries are not taken",
415 "into account. Reference of potential is taken to be the left side of",
416 "the box. It is also possible to calculate the potential in spherical",
417 "coordinates as function of r by calculating a charge distribution in",
418 "spherical slices and twice integrating them. epsilon_r is taken as 1,",
419 "but 2 is more appropriate in many cases."
422 static int axis
= 2; /* normal to memb. default z */
423 static const char *axtitle
= "Z";
424 static int nslices
= 10; /* nr of slices defined */
425 static int ngrps
= 1;
426 static gmx_bool bSpherical
= FALSE
; /* default is bilayer types */
427 static real fudge_z
= 0; /* translate coordinates */
428 static gmx_bool bCorrect
= 0;
430 { "-d", FALSE
, etSTR
, {&axtitle
},
431 "Take the normal on the membrane in direction X, Y or Z." },
432 { "-sl", FALSE
, etINT
, {&nslices
},
433 "Calculate potential as function of boxlength, dividing the box"
434 " in this number of slices." },
435 { "-cb", FALSE
, etINT
, {&cb
},
436 "Discard this number of first slices of box for integration" },
437 { "-ce", FALSE
, etINT
, {&ce
},
438 "Discard this number of last slices of box for integration" },
439 { "-tz", FALSE
, etREAL
, {&fudge_z
},
440 "Translate all coordinates by this distance in the direction of the box" },
441 { "-spherical", FALSE
, etBOOL
, {&bSpherical
},
442 "Calculate spherical thingie" },
443 { "-ng", FALSE
, etINT
, {&ngrps
},
444 "Number of groups to consider" },
445 { "-correct", FALSE
, etBOOL
, {&bCorrect
},
446 "Assume net zero charge of groups to improve accuracy" }
448 const char *bugs
[] = {
449 "Discarding slices for integration should not be necessary."
452 double **potential
, /* potential per slice */
453 **charge
, /* total charge per slice */
454 **field
, /* field per slice */
455 slWidth
; /* width of one slice */
456 char **grpname
; /* groupnames */
457 int *ngx
; /* sizes of groups */
458 t_topology
*top
; /* topology */
460 atom_id
**index
; /* indices for all groups */
461 t_filenm fnm
[] = { /* files for g_order */
462 { efTRX
, "-f", NULL
, ffREAD
}, /* trajectory file */
463 { efNDX
, NULL
, NULL
, ffREAD
}, /* index file */
464 { efTPX
, NULL
, NULL
, ffREAD
}, /* topology file */
465 { efXVG
, "-o", "potential", ffWRITE
}, /* xvgr output file */
466 { efXVG
, "-oc", "charge", ffWRITE
}, /* xvgr output file */
467 { efXVG
, "-of", "field", ffWRITE
}, /* xvgr output file */
470 #define NFILE asize(fnm)
472 if (!parse_common_args(&argc
, argv
, PCA_CAN_VIEW
| PCA_CAN_TIME
| PCA_BE_NICE
,
473 NFILE
, fnm
, asize(pa
), pa
, asize(desc
), desc
, asize(bugs
), bugs
,
480 axis
= toupper(axtitle
[0]) - 'X';
482 top
= read_top(ftp2fn(efTPX
, NFILE
, fnm
), &ePBC
); /* read topology file */
484 snew(grpname
, ngrps
);
488 rd_index(ftp2fn(efNDX
, NFILE
, fnm
), ngrps
, ngx
, index
, grpname
);
491 calc_potential(ftp2fn(efTRX
, NFILE
, fnm
), index
, ngx
,
492 &potential
, &charge
, &field
,
493 &nslices
, top
, ePBC
, axis
, ngrps
, &slWidth
, fudge_z
,
494 bSpherical
, bCorrect
, oenv
);
496 plot_potential(potential
, charge
, field
, opt2fn("-o", NFILE
, fnm
),
497 opt2fn("-oc", NFILE
, fnm
), opt2fn("-of", NFILE
, fnm
),
498 nslices
, ngrps
, (const char**)grpname
, slWidth
, oenv
);
500 do_view(oenv
, opt2fn("-o", NFILE
, fnm
), NULL
); /* view xvgr file */
501 do_view(oenv
, opt2fn("-oc", NFILE
, fnm
), NULL
); /* view xvgr file */
502 do_view(oenv
, opt2fn("-of", NFILE
, fnm
), NULL
); /* view xvgr file */