[gmx-developers] Re: Cubic splines - again

[Berk Hess]

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.