differences between pair_style table and non-tabulated potentials


I have been working on porting a dissociative water model (Mahadevan &
Garofalini, J Phys Chem B 111, 2007, 8919; in short, it's a modified
coul/wolf, LJ 1/r^6 attraction, short-range repulsion, and
Stillinger-Weber 3body, all non-bonded interactions with no bond
constraints) to LAMMPS (5Mar12) and have run into an odd problem.

I am finding that if I call my new pair_style (dwp/coul/wolf) directly
under NPT conditions, like

pair_style hybrid/overlay dwp/coul/wolf 0.224215247 5.0 10.0 sw
pair_coeff * * sw dwp.sw.kcal O H
pair_coeff 1 1 dwp/coul/wolf 6116.99331 608.245029 24.0 0.61000000
pair_coeff 1 2 dwp/coul/wolf 32859.0488 0.00000000 24.0 0.20009697
pair_coeff 2 2 dwp/coul/wolf 0.00000000 0.00000000 24.0 0.61500000

the system expands continuously and voids open up in my liquid water.
However, if I use this pair_style to generate tables using pair_write,
then read those same tables in using pair_style table, it works fine:

pair_style dwp/coul/wolf 0.224215247 5.0 10.0
# note that the sw component is omitted since pair_write doesn't like sw
pair_coeff 1 1 6116.99331 608.245029 24.0 0.61000000
pair_coeff 1 2 32859.0488 0.00000000 24.0 0.20009697
pair_coeff 2 2 0.00000000 0.00000000 24.0 0.61500000
pair_write 1 1 20000 r 0.4 10.0 tableOO.tab DWP_O_O -0.904 -0.904
pair_write 1 2 20000 r 0.4 10.0 tableHO.tab DWP_H_O +0.452 -0.904
pair_write 2 2 20000 r 0.4 10.0 tableHH.tab DWP_H_H +0.452 +0.452
run 0
pair_style hybrid/overlay table linear 20000 sw
pair_coeff * * sw dwp.sw.kcal O H
pair_coeff 1 1 table tableOO.tab DWP_O_O
pair_coeff 1 2 table tableHO.tab DWP_H_O
pair_coeff 2 2 table tableHH.tab DWP_H_H

In both cases, all other input parameters are the same, so shouldn't
the results be equivalent?

My system behaves similarly to how it would if external tension is
applied, so I am inclined to believe that the barostat I'm using here
(fix 1 all npt temp 298.0 298.0 100.0 aniso 1.0 1.0 1000) is doing
something funny. However, I'm using the same barostat in the
tabulated forces. Other than that, the only difference of which I can
think are the charges in my data file (-0.904 and +0.452, real units)
are somehow causing problems when the potential is calculated on the

Any insight would be greatly appreciated!

Glenn K. Lockwood

Does your new pair style compute the same pressure (virial)
as when you use pair_style table? If it doesn't then the box
will respond differently with NPT. If it does, then I don't see
how fix npt would care as to which pair style you are using.
Compare the thermo output from the 2 runs carefully. There
must be some difference or the 2 cases would evolve the same.


Sure enough, I was feeding the wrong values into ev_tally(). After
correcting that, our water model seems to produce reasonable values in
LAMMPS. Thank you very much.

Glenn K. Lockwood