[gmx-users] tabulated potential energy=nan for r=0 nm and charged atoms

[chris.neale at utoronto.ca]

Dear users:

with assistance from Berk on the developers list, I am posting a  
work-around to this issue.

the nonbonded interactions calculate the nonbonded distance, r, using  
the gmx_invsqrt function. Thus, the optimized kernel from 4.5.4 and  
4.0.7 both give nan when two charges occupy the space, even when using  
tables that specify the force and potential energy to be zero.

I had expected that I would need to modify the generic kernel source,  
but that turned out to be unnecessary -- just using it was enough.

For some reason, using "export GMX_NOOPTIMIZEDKERNELS=1" prior to  
running mdrun solves the issue for 4.5.4 but not 4.0.7 (on my  
architecture). To solve the problem in v4.0.7 (and this also works for  
v4.5.4) I must "export GMX_NB_GENERIC=1; export GMX_NO_SOLV_OPT=1"  
prior to running mdrun. Note that the generic nonbonded kernel runs  
1.1x slower then the unoptimized kernel, which itself runs 2.4x slower  
than the optimized kernel (on my architecture for this system).

Thank you,
Chris.

-- original message --

Dear Users:

I am trying to work with tabulated potentials for the first time. I
would like to allow 2 charges to be at the same spot. Imagine taking a
sodium and chloride ion and removing the LJ entirely and I want the
system to evaluate energies correctly and be stable during the
simulation. It was my intention to add a sort-of soft-core to the
charge-charge interactions at very close distances so that the system
did not crash. However, my testing didn't even get that far. I find
that the Coulomb(SR) energy = nan even when using tables that suggest
it should be zero.

I suppose that there is some other place in the code where this
problem develops?

While looking at share/gromacs/top/table6-12.xvg it seems that this
should not lead to energies of nan, but it does.

$ head share/gromacs/top/table6-9.xvg
#
# Table AB1: ndisp=6 nrep=9
#
0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
2.0000000000e-03   0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
4.0000000000e-03   0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
6.0000000000e-03   0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
8.0000000000e-03   0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
1.0000000000e-02   0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00
1.2000000000e-02   0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00   0.0000000000e+00 0.0000000000e+00

and therefore, according to manual 4.5.4 page 151 equation 6.23, the
[(qi*qj)/(4*pi*epsilon)]*f(rij) coulomb SR value should equal zero. I
tried a few things, like modifying the values so that they are
constant from 0 to 0.04 nm in the table.xvg file, but this did not help.

#################### here is the .top file that I used along with
ffoplsaa to test out this idea:

#include "ffoplsaa.itp"
#include "na.itp"

[ system ]
; name
tables

[ molecules ]
; name  number
MOL_A          1
MOL_B          1

#################### and here is the .itp file that I used along with
ffoplsaa to test out this idea:

[ atomtypes ]
   aaa   AA  8     15.99940    0.0       A    0.00000      0.00000

[ moleculetype ]
; molname       nrexcl
MOL_A             1

[ atoms ]
; id    at type   res nr  residu name     at name         cg nr   charge
1       aaa       1       MLA             AA              1        1.0

[ moleculetype ]
; molname       nrexcl
MOL_B             1

[ atoms ]
; id    at type   res nr  residu name     at name         cg nr   charge
1       aaa       1       MLB             BB              1        -1.0

#################### And the .gro file:

title
      2
      1MLA     AA    1   0.000   0.000   0.000
      1MLB     BB    1   0.000   0.000   0.000
    5.00000  5.00000  5.00000

################### And my .mdp file is

nsteps = 0
ns_type = grid
rlist = 1
rcoulomb = 1
rvdw = 1
coulombtype = user
tc_grps             =  System
tau_t               =  1.0
ld_seed             =  -1
ref_t = 300
gen_temp = 300
gen_vel = yes
unconstrained_start = no
gen_seed = -1
Pcoupl = berendsen
pcoupltype = isotropic
tau_p = 4
compressibility = 4.5e-5
ref_p = 1.0

############################################################################

Still, when I run a 0-step mdrun (mdrun -deffnm md1 -table
mod2.table.xvg), I get:

     Energies (kJ/mol)
          LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
      0.00000e+00            nan            nan            nan            nan
      Temperature Pressure (bar)
              nan            nan

     # I get the same thing if the two charges are like or dislike...

But if I then modity my topologies to turn the +1 charge to zero, I
get sensible output:

          LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
      0.00000e+00    0.00000e+00    0.00000e+00    5.27174e+00    5.27174e+00
      Temperature Pressure (bar)
      4.22694e+02    4.66876e-01

Also, if I leave the charges on and separate the two ions by 0.001 nm,
I also don't get the nan values.

     Energies (kJ/mol)
          LJ (SR)   Coulomb (SR)      Potential    Kinetic En.   Total Energy
      0.00000e+00    0.00000e+00    0.00000e+00    6.64240e-01    6.64240e-01
      Temperature Pressure (bar)
      5.32595e+01    5.88265e-02


Thank you for your assistance,
Chris.