[gmx-developers] likely an error in gromacs

[Dongsheng Zhang]

On Sat, 2006-08-19 at 19:54 +0200, David van der Spoel wrote:
> Dongsheng Zhang wrote:
> > On Sat, 2006-08-19 at 08:22 +0200, David van der Spoel wrote: 
> >> Dongsheng Zhang wrote:
> >>> Dear GMX developers,
> >>>
> >>> I used a user-defined potential for non-bonded interactions. By chance I
> >>> found there was likely an error in potential energy calculation. To make
> >>> the story short, I will just tell you what different tests I have done.
> >>>
> >>> In the test, I calculated the potential energy between one pair by using
> >>> user-defined potential (tabulated potential). The coordinates of these
> >>> two particles are:
> >>> 42.04700       37.95200       42.47800
> >>> 42.99200       37.88300       42.77700
> >>>
> >>> The distance between them is 0.99357250826 (got from my own code), it is
> >>> between 0.992 and 0.994
> >>>
> >>> I varied the tabulated potential to see how gromacs and my own code to
> >>> correspond:
> >>>
> >>> Table one (only relevant part): (To make it simple, I set the second
> >>> derivative term zero)
> >>>
> >>> 0.992	0	0	0		0		0	0
> >>> 0.994	0	0	10		0		0	0
> >>>
> >>> The result from my code is 7.86254128, The result from gromacs (by
> >>> defining energy groups to get the interaction of one pair) is 7.802608
> >>>
> >>> Table two:
> >>>
> >>> 0.992	0	0	0		0		0	0
> >>> 0.994	0	0	100		0		0	0
> >>>
> >>> The result from my code is 7.86254128e+01. The result from gromacs is
> >>> 78.562798, which is not 10 times of 7.802608.
> >>>
> >>>
> >>>
> >>> Could anyone give me a hint what happens in gromacs? Thank you a lot!
> >> have you implemented the same cubic spline interpolation scheme as in 
> >> gromacs?
> >> see chapter 6 in the manual.
> > 
> > Yes. I did. I have read Chapter 6 in the manual to build the potential
> > table and write my own code. The only reason to set the second
> > derivative term zero is to make the calculation simple and easy
> > checking.
> > 
> > 
> > If I vary the second derivative term, both gromacs and my own code can
> > responds correctly.
> > 
> > Continue the test from my previous email,
> > 
> > Table 3:
> > 
> > 0.992	0	0	0		0	0	0
> > 0.994	0	0	10		100000	0	0
> > 
> > 
> > The result from my code is 7.84252826 (the result from table 1 is
> > 7.86254128). The result from gromacs is 7.782593 (the result from table
> > 1 is 7.802608).
> > 
> > 
> > Table 4:
> > 0.992	0	0	0		100000	0	0
> > 0.994	0	0	10		100000	0	0
> > 
> > 
> > The result from my code is 7.82892957, the result from gromacs is
> > 7.768992.
> > 
> > Table 5:
> > 
> > 
> > 
> >      0.978     1.02249     2.13803    -1.14279    -50.1809     1.30597
> > 213.001
> >       0.98     1.02041     2.12496    -1.12887    -49.3674     1.27435
> > 206.995
> >      0.982     1.01833     2.11201    -1.11514    -48.5688     1.24355
> > 201.17
> >      0.984     1.01626     2.09916    -1.10161    -47.7847     1.21355
> > 195.521
> >      0.986      1.0142     2.08641    -1.08827    -47.0147     1.18434
> > 190.041
> >      0.988     1.01215     2.07376    -1.07512    -46.2587     1.15589
> > 184.726
> >       0.99      1.0101     2.06122    -1.06216    -45.5164     1.12818
> > 179.569
> >      0.992     0		  0    0	   0	       0		  0
> >      0.994     0		  0    10	   0	       0		  0
> >      0.996     1.00402     2.02419    -1.02434    -43.3685     1.04927
> > 165.004
> >      0.998       1.002     2.01205    -1.01208    -42.6781     1.02431
> > 160.434
> > 
> > 
> > The result from my code is 7.86254128 (the same as that from table 1).
> > The result from gromacs is 7.212971, which is different from the result
> > from table 1 (7.802608). Table 5 looks a little messy. The key point is
> > that the entries for r = 0.992 and 0.994 is the same as table 1.
> > Therefore I expect that I can get the same result as table 1, but
> > gromacs gives me different results. I hope my test results can give you
> > some clues where gromacs made a mistake or I had made a mistake. If you
> > need more information, please let me know.
> > 
> > 
> > All the best!
> > 
> > 
> > 
> > Dongsheng
> > 
> Please submit a bugzilla and upload the different tables and also a tpr 
> file.
> 
I just report it. Thank you for help!



> -- 
> 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://folding.bmc.uu.se
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the 
> www interface or send it to gmx-developers-request at gromacs.org.