CFOUR has a utility program, xpolyreg, which can fit polynomials to multidimensional (up to six) data sets. This can be useful for fitting potentials or property surfaces, especially when the GRID option is used. The input is as follows:
First line: Information about weighting.
Second line: Order of polynomial to fit to the data.
Third line: Number of points to fit (n).
Fourth line: Dimensionality of the dataset.
Fifth line: Conversion factor applied to the polynomial.
Lines six through end: n lines, each containing the independent variables followed by the dependent variable and (perhaps) the weight.
If the first line of the input file is 0, uniform weighting is used.
If the first line of the input file is -1, weights are explicitly supplied in the input.
If the first line of the input file is -2, weights are given by 1 - (dk/R)x, where R is the range of ordinate values (z(max) - z(min)) and dk is (zk - min(z)). In this case, the value of x is read from the next line and the order of the polynomial therefore becomes the third line of the input.
Output from the program starts with the various data points, giving both the actual dependent variable value as well as that which is generated by the fitting polynomial. After this, the polynomial is printed according to the format given below.
Let us assume that a third-degree polynomial is fit to a two-dimensional data set (for example, these might be the bond length and bond angle of water) with 125 points. The fitting polynomial will be given in the following way:
0 0 value00
1 0 value10
2 0 value20
3 0 value30
0 1 value01
1 1 value11
2 1 value21
3 1 value31
0 2 value02
1 2 value12
2 2 value22
3 2 value32
0 3 value03
1 3 value13
2 3 value23
3 3 value33
where the values are the derivatives of the dependent variable with respect to the independent variables. That is, the polynomial would be
f= value00 + value10 * x1 + 1/2 value20 x12 + ...
where the derivatives are multiplied by whatever constant is given in the "fifth line" (see above) of the input. For example, if the dataset above (two dimensions, third-order unweighted fit) is for energies in Hartrees, but for a polynomial expressed in wavenumbers, the input would begin with
x1 value of first point x2 value of first point f value of first point ...
The executable is invoked with
xpolyreg < input_file
where input_file contains the data described above.
Sample input and output can be downloaded.