[gmx-users] dipole correction in 3.1-beta

[David L. Bostick]

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	
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-