[gmx-users] unexpected large difference in potential energy calculation by different method: cut-off and switch

[Dongsheng Zhang]

Hi, 

I am comparing different methods to calculate the coulomb and VDW
potential energy. I think the results for cut-off and switch should be
close. However, the difference is very large. Could anyone tell me what
happen? Detailed information is as follows:

After I finish position restrained MD simulation, I use full2.mdp and
full3.mdp to make full2.tpr and full3.tpr, and run them. I get very
different results. Then, I think it maybe is a good idea to rerun a
trajectory file by using full2.tpr and full3.tpr, and compare the
results. I find the potential energies are very different, mainly from
Coulomb (SR).

The commands I used are:
mdrun -s full2.tpr -rerun full2.trr -e full2 -o full2 -x full2 -g full2
g_energy -f full2.edr -o energy2  (choose 5-10 )

mdrun -s full3.tpr -rerun full2.trr -e full3 -o full3 -x full3 -g full3
g_energy -f full3.edr -o energy3  (choose 5-10)

part of full2.mdp is:

; NEIGHBORSEARCHING PARAMETERS
nstlist                  = 5            ; nblist update frequency
ns_type                  = grid         ; ns algorithm (simple or grid)
pbc                      = xyz
rlist                    = 1.1          ; nblist cut-off        


; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics: Cut-off, 
coulombtype              = Switch
fourierspacing           = 0.1
pme_order                = 3
rcoulomb                 = 1.0
rcoulomb_switch          = 0.9
; Dielectric constant (DC) for cut-off or DC of reaction field
epsilon-r                = 1.0

; Method for doing Van der Waals
vdw-type                 = Switch
; cut-off lengths       
rvdw                     = 1.0
rvdw_switch              = 0.9

; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = no; EnerPres



; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling  
Tcoupl                   = nose-hoover ;berendsen
; Groups to couple separately
tc-grps                  = protein SOL
; Time constant (ps) and reference temperature (K)
tau_t                    = 0.1   0.1
ref_t                    = 300   300
; Pressure coupling     
Pcoupl                   = Parrinello-Rahman ;berendsen
Pcoupltype               = isotropic 
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p                    = 1.5
compressibility          = 4.5e-5
ref_p                    = 1.0   



part of full3.mdp is:

; NEIGHBORSEARCHING PARAMETERS
nstlist                  = 5            ; nblist update frequency
ns_type                  = grid         ; ns algorithm (simple or grid)
pbc                      = xyz
rlist                    = 1.0          ; nblist cut-off        


; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics: Cut-off, 
coulombtype              = Cut-off
fourierspacing           = 0.1
pme_order                = 3
rcoulomb                 = 1.0
rcoulomb_switch          = 0.9
; Dielectric constant (DC) for cut-off or DC of reaction field
epsilon-r                = 1.0

; Method for doing Van der Waals
vdw-type                 = cut-off
; cut-off lengths       
rvdw                     = 1.0
rvdw_switch              = 0.9

; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = no; EnerPres



; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling  
Tcoupl                   = nose-hoover ;berendsen
; Groups to couple separately
tc-grps                  = protein SOL
; Time constant (ps) and reference temperature (K)
tau_t                    = 0.1   0.1
ref_t                    = 300   300
; Pressure coupling     
Pcoupl                   = Parrinello-Rahman ;berendsen
Pcoupltype               = isotropic 
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau_p                    = 1.5
compressibility          = 4.5e-5
ref_p                    = 1.0   



I am using Gromacs 3.2.1 and amber94 FF. The water model is
ffamber_tip3p.itp


If anyone want more information, I will be more than happy to tell you. 


Thank you for your help!


-- 
Dongsheng Zhang <dong at pampas.chem.purdue.edu>