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.
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
