When dealing with charge density data or wave function data stored in a cube file, I sometimes would like to find a charge density in a specific point in the space. However, only density on grid is stored in the file. Interpolation is necessary for this purpose. In scipy, scipy.interpolate.RegularGridInterpolator can be used for 3d interpolation, but only "nearest" and "linear" method is available. In order to get high presicion, Tricubic interpolation could help. In two packages, this function is available. pytricubic and eqtools, I tried both and both works fine. I would recommend the eqtools since it is easy to install and doesnt have a dependence on Eigen or boost.
Here I give a small test script for comparison of nearst, linear, and tricubic interpolation accuracy.
from scipy.interpolate import RegularGridInterpolator
from eqtools.trispline import Spline
n = 30
xgrid = np.linspace(-1.0, 1.0, n)
ygrid = np.linspace(-1.0, 1.0, n)
zgrid = np.linspace(-1.0, 1.0, n)
f = np.zeros((n,n,n), dtype = 'float')
for i in range(n):
for j in range(n):
for k in range(n):
f[i][j][k] = np.sin(5*xgrid[i])*\
np.sin(5*ygrid[j])*\
np.sin(5*zgrid[k])
interpolating_function_tricubic = Spline(xgrid,ygrid,zgrid, f)
interpolating_function = RegularGridInterpolator((xgrid,ygrid,zgrid), f, method="linear")
interpolating_function_nearest = RegularGridInterpolator((xgrid,ygrid,zgrid), f, method="nearest")
m = 20
error_nearest = 0.0
error = 0.0
error_tricubic = 0.0
xgrid_d = np.linspace(-0.5, 0.5, m)
ygrid_d = np.linspace(-0.5, 0.5, m)
zgrid_d = np.linspace(-0.5, 0.5, m)
for i in range(m):
for j in range(m):
for k in range(m):
error += np.abs(interpolating_function([xgrid_d[i], ygrid_d[j], zgrid_d[k]]) - \
np.sin(5*xgrid_d[i])*\
np.sin(5*ygrid_d[j])*\
np.sin(5*zgrid_d[k]))[0]
error_nearest += np.abs(interpolating_function_nearest([xgrid_d[i], ygrid_d[j], zgrid_d[k]]) - \
np.sin(5*xgrid_d[i])*\
np.sin(5*ygrid_d[j])*\
np.sin(5*zgrid_d[k]))[0]
error_tricubic += np.abs(interpolating_function_tricubic.ev(\
xgrid_d[i], ygrid_d[j], zgrid_d[k]) -\
np.sin(5*xgrid_d[i])*\
np.sin(5*ygrid_d[j])*\
np.sin(5*zgrid_d[k]))[0]
print "average error for nearest, linear, and tricubic interpolation"
print error_nearest / m**3, error / m**3, error_tricubic / m **3
It gave such result, showing better accuracy for tricubic interpolation.
average error for nearest, linear, and tricubic interpolation
0.0479790688641 0.0103316526679 0.000350221662438