[gmx-developers] likely an error in gromacs

[Dongsheng Zhang]

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





> 
> > 
> > 
> > All the best!
> > 
> > 
> > Dongsheng
> > 
> > 
> > _______________________________________________
> > 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.
> 
>