[gmx-users] energy conservation with lookup table

[Zhe Wu]

Dear GROMACS developers and users,

I observed a problem with energy conservation in NVE simulation with
nonbond potential lookup table.
I use the same setup for the L-J particles, one in the "analytical"
way in the *.itp file as functional form 1
(with vdw_type= switch), the other using User define table provided in
/share/gromacs/top/table6-12.xvg (with vdw_type=User).

For unknown reason(s), the "analytical" one attains fairly good energy
conservation, while the lookup
table one doesn't (even after I changed the spacing interval from
0.002 nm to 0.001 nm).

Is it a bug or a silly mistake I made? I have been checking it for the
whole day, search through the
manual and mailing list and I just cannot find anything particularly
helpful. I am using GROMACS 4.0.5
and I use 3.3.3 (change the table and the mdp file of accordingly of
course) to make a test and the
energy still doesn't conserve with table. But I am use the single
precision for Gromacs executable and
haven't check the double precision yet.

Did anybody observe a similar problem before? While using potential
lookup table, will it still have shift
scheme or switch scheme? Any idea or suggestion will be really
appreciated! Thank you in advanced
for your help.


Below are the details for the simulation (in GMX 4.0.5):
4 core, red hat, for 1ns NVE simulation;

relevant lines in itp file from MARTINI (for both "analytical" and table case):
------------------------------------------------------------------------------------------
[ nonbond_params ]
 P1 	P1 	1 	0.19402E-00 	0.20914E-02 ; almost attractive

; test
;[ moleculetype ]
; molname       nrexcl
;SOL             2

;[ atoms ]
;   nr   type  resnr residue  atom   cgnr     charge       mass
;     1     P1      1    SOL     OW      1      0.0
--------------------------------------------------------------------------------------

mdp file for the "analytical" one, for the relevant lines:
----------------------------------------------------------------------
; VARIOUS PREPROCESSING OPTIONS =
title                    = Martini
cpp                      = /usr/bin/cpp

; RUN CONTROL PARAMETERS =
integrator               = md
; start time and timestep in ps =
tinit                    = 0.0
dt                       = 0.005   ; I am using a CG model
nsteps                   = 200000
; number of steps for center of mass motion removal =
nstcomm                  = 1

; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist                  = 10
; ns algorithm (simple or grid) =
ns_type                  = grid
; Periodic boundary conditions: xyz or none =
pbc                      = xyz
; nblist cut-off         =
rlist                    = 1.8   ; try to make the cut-off long enough

; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype              = shift   ; no electrostatic in the system
;fourierspacing           = 0.2
rcoulomb_switch          = 0.0
;pme_order                = 6
rcoulomb                 = 1.6
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r                = 1.4
; Method for doing Van der Waals =
vdw_type                 = switch   ; can still attain reasonable
energy conservation
; cut-off lengths        =
rvdw_switch              = 1.2
rvdw                     = 1.6
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr                 = No

; NVE simulation
; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling   =
tcoupl                   = no
; Groups to couple separately =
tc-grps                  = System
; Time constant (ps) and reference temperature (K) =
tau_t                    = 0.05
ref_t                    = 300
; Pressure coupling      =
Pcoupl                   = no
Pcoupltype               = semiisotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau_p                    = 0.2  0.2
compressibility          = 3e-5 3e-5
ref_p                    = 1.0  1.0

; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel                  = no
gen_temp                 = 105
gen_seed                 = 473529
------------------------------------------------------------------------------------------

results from g_energy for "analytical" one:
------------------------------------------------------------------------------------------
Energy                        Aver.              RMSD                      Drift
Potential                    -25618.5          150.676                0.00821594
Total Energy             -16529              0.976822             0.000626028
------------------------------------------------------------------------------------------

mdp file for the one using table. All lines are the same except the
following two:
------------------------------------------------------------------------------------------
table-extension                = 1.2  ; nm
vdw_type                         = User
--------------------------------------------------------------------------------------

The table is taken from the package 4.0.5 in /share/gromacs/top/table6-12.xvg.

results from g_energy for "analytical" one:
------------------------------------------------------------------------------------------
Energy                        Aver.              RMSD                      Drift
Potential                   -30257.3          1516.79                 -5.79857
Total Energy             -25641             3017.7                   -11.5046
------------------------------------------------------------------------------------------