[gmx-users] Re: gmx-users Digest, Vol 20, Issue 9

[Dongsheng Zhang]

Dear David,

I have tried it. The difference in potential energy is very large. 

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

Could you please give me a clue where the difference in potential energy is from?

Thanks! 

> Date: Thu, 01 Dec 2005 22:44:16 +0100
> From: David van der Spoel <spoel at xray.bmc.uu.se>
> Subject: Re: [gmx-users] Re: Re:unexpected large difference in
> 	potential energy	calculation by different method: cut-off and switch
> To: Discussion list for GROMACS users <gmx-users at gromacs.org>
> Message-ID: <438F6EB0.9020003 at xray.bmc.uu.se>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
> 
> Dongsheng Zhang wrote:
> > Dear David,
> > 
> > Thank you for your reply. When I read your email at the first time. I
> > think you have answered my question. But when I think it again today, I
> > find I still have a question:
> > 
> > I can understand that there will be some artifacts to calculate forces
> > by using switch function in MD simulation because switch function only
> > smooths the potential. My question is: If I calculate the total
> > potential energy of a system by rerunning a trajectory file with two
> > different methods (Cut-off and switch), should I get close results
> > because I don't need to calculate the force to determine the next
> > configuration of the system? 
> Interesting question. Please try it.
> It depends on the options etc.
> > 
> > 
> > 
> > Dongsheng Zhang wrote:
> > 
> >>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:
> > 
> > 
> > Switched potentials should be avoided. In GROMACS the switch is atom 
> > based which gives huges artifacts, in other codes you can choose group 
> > based switches which has less artifacts but introduces a horrible 
> > arbitrariness in the definition of the groups. For more info see:
> > 
> > The Origin of Layer Structure Artifacts in Simulations of Liquid Water
> > David van der Spoel and Paul J. van Maaren, in press
> > J. Chem. Theor. Comput. http://dx.doi.org/10.1021/ct0502256
> > 
> > 
> 
> 
> -- 
> David.
> ________________________________________________________________________
> David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
> Dept. of Cell and Molecular Biology, Uppsala University.
> Husargatan 3, Box 596,  	75124 Uppsala, Sweden
> phone:	46 18 471 4205		fax: 46 18 511 755
> spoel at xray.bmc.uu.se	spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> 

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