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

[chris.neale at utoronto.ca]

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.