[gmx-developers] Re: Cubic splines - again

[Mathias PUETZ]

Hi Berk,

I don't really understand what you are getting at with your "as well" 
statement below.

With cubic splines you simply can't have 
  - the function values
  - the first order derivates
  - and the second order derivatives
right with the same set of coefficients. if you want all at the same time 
with the same set of coefficients,
you need quintic splines, which are more expensive (it wouldn't be a 
catastrophy for the
computation time, but would mean a tremendous amount of code hacking).

True, getting the 2nd order derivates right minimizes the error in the 
potential calculation to o(h**3), h=x_i-x_i+1,
but the point is, it does *not* minimize the error in the forces at the 
same time.

So it essentially boils down two these four choices:

 1. get better force accuracy and keep the current accuracy for the 
potential, but you have to pay for it with quintic splines
 2. have either accurate potentials and inaccurate forces (the way it is 
right now)
 3. improve force accuracy at the expense of some potential accuracy (what 
I proposed)
 4. approximate potentials and forces using separate cubic spline 
polynomials

I don't consider the 4th choice a good one, but I listed it for 
completeness anyway.
It means you would have two different sets of coefficient tables, one for 
the potentials and
one for the forces, which, in turn, would lead to a change in the inner 
loops of the force kernels.
Such an approach would also introduce a small systematical bias, since the 
simulated forces
are no longer exact derivatives of the measured potential function, which 
I, personally, would
stay away from.

So overall I would go for the 3rd choice, because

 - I get better energy conservation during the integration
 - I keep force and potential consistent both in absolute values and in 
accuracy (both to o(h**2))

the only thing I loose is some accuracy in the potential, which
 - is approximate anyway
 - can always be improved by increasing the number of table sites

I consider the first two arguments more important, but then again, I am 
not certain, that this really is the best way to go,
since I don't really have strong experience using tabulated potentials and 
have not performed any studies on this.

Mathias PUETZ wrote:

>
> Hi,
>
> the cubic spline function
>
> spline_i (x) = a+b(x - x_i) + c (x - x_i)**2 + d (x - x_i)**3
>
>  has four variable coefficients: a,b,c,d.
>
>  For each table segment i you have four linear equations which 
> determine the coefficients for the spline function spline_i(x)
>
> v(x_i) = spline_i(x_i)
> v(x_i+1) = spline_i(x_i+1)
> v''(x_i) = spline_i''(x_i)
> v''(x_i+1) = spline_i''(x_i+1)
>
> solving these equations for a,b,c,d gives exactly the spline formula 
> you are using.
>
> Since you force v and v'' to be continous, v' (the forces) cannot be 
> continuous at x_i
> unless v(x) itself is a polynomial function with a degree of less than 
> 4, which it typically isn't.
> The error will be small in most cases, but it is noticable if you 
> print out the forces
> calculated by your spline formula and the exact force function in the 
> vicinity of the x_i
>
> I believe it would be more appropriate to enforce
>
> v(x_i) = spline_i(x_i)
> v(x_i+1) = spline_i(x_i+1)
> v'(x_i) = spline_i'(x_i)
> v'(x_i+1) = spline_i'(x_i+1)
>
> which will give you slightly different spline coefficients, but 
> guarantees that both v(x) and v'(x)
> will be continuous at x_i.
>
You are right.

But I think we would still prefer to have continuous second derivatives 
as well.
So the problem is then not in the interloops, but in the generation of 
the tables.
There we should not use the second derivative of the potential,
but rather the second derivative of the cubic spline fit of the potential.

Berk.
Viele Grüsse / Best regards,
Dr. Mathias Pütz

IT Specialist for Application Perfomance

Deep Computing - Strategic Growth BusinessDeep 
IBM Systems & Technology Group

e-mail:  mpuetz at de.ibm.com
mobile: + 49-(0)160-7120602
fax:         + 49-(0)6131-84-6660

snailmail:
  IBM Deutschland GmbH
  Department B458
  Hechtsheimer Str. 2 / Building 12
  55131 Mainz
  Germany
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20060822/f624d925/attachment.html>