[gmx-developers] function gmx_pme_do will cause segmentation fault for mdp option: epsilon-r=0

[Mark Tianwu Zang]

Dear all,
I found that in force.c line 1324(for gromacs-4.5.7) or line 1966(for
gromacs-4.6.2), a variable called "eterm"(for gromacs -4.5) or elfac(for
gromacs-4.6) is defined as

eterm = ONE_4PI_EPS0/pme->epsilon_r*tmp1[kx]*denom[kx];

or

elfac = ONE_4PI_EPS0/pme->epsilon_r;

where pme->epsilon_r = ir->epsilon_r, which corresponds to the mdp option:
epsilon-r.
I have not tested on gromacs-4.6 yet but I have already tested on
gromacs-4.5, and as expected, if epsilon-r is made zero in *.mdp, a
segmentation fault will appear at function gmx_pme_do. (Since the two
pieces of code above are almost the same, I think I probably get the same
result for groamcs-4.6).
The physical meaning of epsilon-r=0 is that I am only interested in an
infinite dieletric system(no Coulomb interactions despite non-zero
charges).

On the other hand, in forcerec.c:2623(version 4.5) or 1609(verion 4.6),
there is another factor: fr->epsfac defined as below:

if(fr->epsilon_r != 0)
 fr->epsfac = ONE_4PI_EPS0/fr->epsilon_r;
else
 fr->epsfac = 0;

which has already considered this possible situation.

However, although the meaning of fr->epsfac is the same as the ones of
variable eterm(the first half) and elfac, this factor (fr->epsfac) will
only be used for calculating non-bonded (short-range) Coulomb interaction,
bonded 1-4 Coulomb interaction and long-range interaction using reaction
field method. For long-range interaction using PME/Ewald
method, fr->epsfac will never be used, and the introduction(actually
recalculation..) of eterm/elfac causes the failure of computing PME.

Thanks,
-Mark
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20130701/da00758c/attachment.html>