[gmx-developers] Re: Cubic splines - again

[Berk Hess]

Mathias PUETZ wrote:

>
> 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 think the spline interpolation of the potential I suggested is 
impractical.
The interpolation will always depend on the endpoints of the table.
We would not want the forces at a certain distance to depend on how
long our table is (even if the effect is very small).

I think Mathias' suggestion is the best solution.
I don't agree with the phrasing of point 3 though.
The force should always be the exact derivative of the potential.
Currently this is not the case, which is bad.
So with cubic splines this leaves only one solution.

If we prescribe V and V' at table points the relative error in the second
derivative at a table point is:
V'''' h^2/(3 V'')
where h is the table spacing (if I did not make a mistake in my math).
The jump in the second derivative is double this.
Even for the r^-12 LJ repulsion this error is very small with h=1/500 nm.

I guess we only need to supply different table values to the innerloops.

But the user tables would need to be changed, as they currently
contain V and V''.

Berk.