[gmx-users] dipole correction in 3.1-beta
David L. Bostick
dbostick at physics.unc.edu
Wed Feb 20 23:26:04 CET 2002
Hello,
This question is directed more at the developers than the users, however..
here goes. I am interested in employing a simulation with slab geometry
and therefore went delving into the code which computes the dipole
correction for the force and energy. Another user Rui Quiao had expressed
interest to the user list for the same thing. As I was looking into this I
noticed the following:
In the file ewald_util.c for version 3.0 the factor placed in front of the
square dipole moment in the electrostatic energy , called dipole_coeff, is
assigned like this:
if(epsilon_surface==0)
dipole_coeff=0;
else
dipole_coeff=M_PI/(1.5*epsilon_surface*ONE_4PI_EPS0*vol);
The obvious error here of ONE_4PI_EPS0 being in the denominator has been
corrected in 3.1-beta like this:
if(epsilon_surface==0)
dipole_coeff=0;
else
dipole_coeff=M_PI*ONE_4PI_EPS0/(1.5*epsilon_surface*vol);
However, this looks a lot like the prefactor for the dipole correction for
the force if epsilon_surface is taken to be the permittivity of free
space..i.e. "epsilon_zero". It seems that still there is something wrong,
however, I am not sure. The dipole term or "shape term" for 3D ewald in
terms of gromacs variables looks like.
Vdipole = dipole_coeff*(dipole)^2
.. leaving out the conversions .. DEBYE2ENM*DEBYE2ENM
where,
dipole_coeff = ONE_4PI_EPS0*2*M_PI / (2*epsilon_surface + 1) / vol;
this leaves the correction to the force to be,
f_pme[i][j]-=2.0*dipole_coeff*mu_tot[j]*charge[i]
again, leaving out the conversion factor DEBYE2ENM
If the albebra is done for the negative gradient of Vdipole with
epsilon_surface = 1 for a vacuum, you obtain a dipole_coeff, that looks
exactly like,
dipole_coeff = ONE_4PI_EPS0*2*M_PI / (3*vol) or,
dipole_coeff = ONE_4PI_EPS0*M_PI / (1.5 * vol)
If you take the gradient of this, you get
f_pme[i][j]-=2.0*charge[i]*dipole_coeff*mu_tot[j]
neglecting conversion,
or,
f_pme[i][j]-=prefactor*mu_tot[j]
where prefactor = 2*charge[i]/(3*vol*epsilon_zero)
= charge[i] / (1.5*vol*epsilon_zero)
The question is why is there no factor of (2*epsilon_surface + 1) which
would account for epsilon_surface = something other than 1.. Is there a
mistake here, or are there things I am neglecting.
Thanks in advance,
David
-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
David Bostick Office: 262 Venable Hall
Dept. of Physics and Astronomy Phone: (919)962-0165
Program in Molecular and Cellular Biophysics
UNC-Chapel Hill
CB #3255 Phillips Hall dbostick at physics.unc.edu
Chapel Hill, NC 27599 http://www.unc.edu/~dbostick
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
More information about the gromacs.org_gmx-users
mailing list