# How to carry out basis set fitting and extrapolation

I am playing around with basis set extrapolations, and i would like to perform some simple tests by doing some parameter fitting and checking the errors.

The functional form of the fitted curve is this

$$E(X) = E_\infty + Ae^{-\alpha X}$$

where $$X$$ is the cardinal number of the basis set used, and $$E$$ is the energy. This function makes sense because the rate of decay empirically fits the convergence (of correlation concistent basis sets) and the exponential term vanishes as $$x\rightarrow\infty$$.

We want to use energies from two low-cardinal numbers, and extrapolate to the $$E_\infty$$ as given by the function and the fitted parameters.

Assume that for the cardinal numbers $$X$$ and $$Y$$ we have the corresponding energies $$E_X$$ and $$E_Y$$. We then have the relationship

$$E_X = E_\infty + Ae^{-\alpha X} \\ E_Y = E_\infty + Ae^{-\alpha Y}$$

By using that $$A = \frac{E_X - E_\infty}{e^{-\alpha X}}$$, we can solve for $$E_\infty$$ to obtain

$$E_\infty = \frac{E_Ye^{-\alpha X} - E_Xe^{-\alpha Y}}{e^{-\alpha X} - e^{-\alpha Y}}$$

and solve for $$A$$ to obtain

$$A = \frac{E_Y - E_X}{e^{-\alpha Y} - e^{-\alpha X}}$$

From what I can see, if all I want is to obtain the extrapolated energy, then I do not need take $$A$$ into any consideration. I just need the cardinal numbers, the energies, and $$\alpha$$. In Python syntax this would be something like this

def extrapolate(X, Y, E_X, E_Y, alpha):
gamma = np.exp(-alpha*X)
omega = np.exp(-alpha*Y)
return (E_Y * gamma - E_X * omega) / (gamma - omega)


### Is the above approach the normal one for curve fitting and extrapolation?

In some tests I find that the extrapolated energy is virtually identical to the energy with the highest cardinal number used in the extrapolation. So either there is no benefit of my extrapolation, or I am doing something wrong.

(Note that I am not using Gaussian basis sets, but a MultiWavelet basis, and that the cardinal numbers are not really cardinal numbers, but a defined precision relative to the CBS. Still the approach should be the same regardless - even if the function I use for the fitting perhaps not is the perfect choice.)

• Hm, I would prefer 3 equations for 3 unknown parameters ( E_inf, A, Alfa) and therefore 3 degrees of freedom. Otherwise 1 parameter would be left for arbitrary choice and the other 2 would depend on its choice, still fitting both equations. Apr 11, 2021 at 12:37
• So given $X$, $Y$, and $Z$, and corresponding energies, I would derive expression for E and A in terms of these?
– Yoda
Apr 11, 2021 at 12:55
• How are you determining the alpha? Apr 11, 2021 at 12:55
• Alpha is fitted to multiple training cases. I guess I take the average value from these.
– Yoda
Apr 11, 2021 at 12:56
• Why not just fit the three points to the curve with scipy curve_fit ? Taking alpha from training means you are introducing an empirical parameters. Apr 11, 2021 at 12:59

If I well understand your question you have a set of experimental points $$(E_1,X_1), (E_2,X_2), ... , (E_n,X_n)$$. You want to fit to this set of points the function : $$E(X) = E_\infty + Ae^{-\alpha X}$$ So, you want to approximately determine the parameters $$E_\infty$$ , $$A$$ and $$\alpha$$.

This requires a non-linear regression calculus. For example see https://mathworld.wolfram.com/NonlinearLeastSquaresFitting.html

This calculus is iterative and requires to set some "guessed" values of the parameters to start the iterative process.

Another method (much simpler but for approximative result only) is not iterative and doesn't require "guessed" values is explained in : https://fr.scribd.com/doc/14674814/Regressions-et-equations-integrales

You could try it in folowing the numerical example below with $$\quad y(x)=E(X)$$ , $$a=E_\infty$$ , $$b=A$$ , $$c=-\alpha$$

Note : This simplified method is not reliable if the number of points is too low.

• +1 This is a good answer with maths, but I suspect OP was looking for a python based solution. Apr 29, 2021 at 8:21
• @Shoubhik R Maiti. Thank you for the judicious comment. In fact I was expecting that the OP could implement on pyton this very simple matrix calculus and gives a comparison of results as feeback. This method has been applied many times in other contexts but not in Chemistry as far as I know. Apr 29, 2021 at 8:33
• Another thing I would point out is that in general it would be difficult to get more than 4-5 points for basis set fitting, because as x (which is determined by basis set size) increases the computational cost increases rapidly, so its not practicaly possible. Does the non-iterative method work for that? I am a chemist mainly, so all of this math stuff is difficult for me :-) Apr 29, 2021 at 9:50
• I cannot answer without having an experimental data of 4 or 5 points. But a-priori I am pessimistic in this case. This depends also on the scatter. It is possible that the method be not convenient in this Chemistry context. In any case this would be interesting for me to see what happens and to look why. Apr 29, 2021 at 10:31
• I would add that even if the non-iterative method doesn't give good results in some cases (not enough points and too high scatter) generaly the values obtained for the parameters $a,b,c$ are at least of the correct magnitude. So they can be used as starting values for a classical iterative non-linear regression. This was used in a few cases instead of "guessed" initial values. Apr 29, 2021 at 10:49