Fix: Implement numerically safe switching function

[jhenin]

replacing function and derivative with order 1 Taylor expansions around indeterminate form 0/0 for original expression

fixes #903

[jhenin]

With epsilon at 1e-7:

and at 1e-8:

[HanatoK]

How about using the second-order derivatives (the same as the PLUMED code)?

Let $f(x) = \frac{1-x^n}{1-x^m}$, then we have

$$ \begin{dcases} \lim_{x \rightarrow1} f(x) &= \frac{n}{m}\\ \lim_{x \rightarrow1} \frac{\rm d}{{\rm d}x} f(x) &= \frac{n^2-nm}{2m}\\ \lim_{x \rightarrow1} \frac{\rm d^2}{{\rm d}x^2} f(x) &= \frac{2 {{n}^{3}}\mathop{+}\left( \mathop{-}\left( 3 m\right) \mathop{-}3\right) {{n}^{2}}\mathop{+}\left( {{m}^{2}}\mathop{+}3 m\right) n}{6 m} \end{dcases} $$

The limitation of the second-order derivative is calculated by the following maxima code (thanks for Gemini!):

limit(diff((1 - x^n) / (1 - x^m), x, 2), x, 1);

$$ \begin{split} f(x\pm \epsilon) &= \lim_{a \rightarrow1} f(a) + ((x-1)\pm \epsilon)\lim_{a \rightarrow1} f'(a)+((x-1)\pm \epsilon)^2\lim_{a \rightarrow1} \frac{f''(a)}{2}\\ &=\frac{n}{m} + ((x-1)\pm \epsilon)\frac{n^2-nm}{2m}+((x-1)\pm \epsilon)^2 \frac{2 {{n}^{3}}\mathop{+}\left( \mathop{-}\left( 3 m\right) \mathop{-}3\right) {{n}^{2}}\mathop{+}\left( {{m}^{2}}\mathop{+}3 m\right) n}{12 m} \end{split} $$

The calculation of the derivative should be also straight forward:

$$ f'(x\pm \epsilon)=\frac{n^2-nm}{2m}+2((x-1)\pm \epsilon)\frac{2 {{n}^{3}}\mathop{+}\left( \mathop{-}\left( 3 m\right) \mathop{-}3\right) {{n}^{2}}\mathop{+}\left( {{m}^{2}}\mathop{+}3 m\right) n}{12 m} $$

[giacomofiorin]

Looks good to me already with the order 1 (the interval is tiny).

[giacomofiorin]

Years ago, we had to wrap all the STL math functions because of some exotic platforms where VMD was being built.

It doesn't seem much of an issue any more, and few instances of unwrapped functions got undetected. I'll fix them in this branch, but we could perhaps consider removing those wrappers.

[jhenin]

Agreed on removing them.

[HanatoK]

I attach my experiment results with cutoff = 6.0 and second-order Taylor expansion here:

My test code can be found in https://godbolt.org/z/fnjnKoTqd.