Tricubic Interpolation - ProkopHapala/FireCore GitHub Wiki

Motivation

Motivation 1: accuracy, consistency & convergnece

  • FireCore use grid accelerated forcefield (GridFF) for acceleration of interaction of molecules with rigid substrate and for simulation of AFM images.
  • GridFF aims to inteprpolate on only nergy (potential) but more importaintly the forces (i.e. derivative of energy with respect to position of particle) which are used for optimization of molecular geometry or running molecular dynamics.
  • Original version of GridFF uses trilinear interpolation of both Energy and Forces, where each component of force is interpolated independnetly.
    • The advantage of this algorithm is that it is numerically simple to implement, fast and highly paralel (vectorized).
    • A disadvantage is that energy and force are not exactly consistent (i.e. force is not exactly derivative of the energy), because both energy and force are piecewise linar functions while, derivatiwe of piecewise energy would be piecewise constant (i.e. step-function like staircase; such function would be not sooth making convergence of relaxation almost impossible)
    • The problem with this inconsistency is that the interpolated forcefield is not necessarily conservative (resp. that it has zero rotation in terms of vector fied analysis). Unlike, forcefield which is exact derivative of a scalar potential function, which is mathematically guarantied to be conservative with zero rotation.
      • This means that particle in such non-conservative forcefield may forever move around (e.g. in some closed orbit) and never converge to minimim (e.g. like some dust particle in vortex of a liquid).
    • I was suspiciios there may be such problem for long time, but I did not encounter it in practice up to recently.
    • Recencly I found I'm not able to converge some organic molecules on NaCl substrate with GridFF to lower residual forces than $F_{conv} < 10^{-3}$ [eV/A].

Motivation 2 - memory requirement & speed

  • Trilinear interpolator of force needs to store both energy and force independently, although a the physical information is already in the potential (if we assume the forces are exactly derivative of the potential). This means that the forcefield require 4x more memory. Especially on GPU this may be relevant (assuming each $1 nm^3$ of grid with spacing $0.1$ A require 1 million gridpoints, with 8 bytes per double, and 4 different forces (Pauli, Londo, Coulomb, Hbond), this is $100^3 * 16 * 8 * 4 * 4 = 2$ GB of VRAM per $1 nm^3$).
    • If the limiting factro is reading data from memory, 4x higher memory size means also 4x higher memory bandwidth (theoretically).
    • Larger memory size can also decrease cashe efficinelcy because size of cache can be limiting factor.

Tricubic interpolation

  • While linear interpolation require 2 parametr per each interval in 1D (typically values at the endpoints of the interval),interpolation by cubic polynominal require 4 parameters.
  • In 3D the interpolation splines are constructed as tensor product of 1D splines, therefore require $n^3$ parameters ( 2^3=8 for trilinear and 4^3=64 for triubic)
    • Form this analysis it is immedietely clear that tricubic interpolation is computationally much more demanding cosidering both number of arithmetic operations but also number of data (parameters) which has to be read from grid for each interpolation.
      • Cashing of the data can be useful if many particles are collocated in the same interval (voxel of the grid), or if the particle moves slowly and stay long time in the same voxel. Nevertheless we did no yet implemented any explicit cashing scheme.

Hermite cubic spline

  • In 1D the Hermite cubic spline seems ideal for interpolation of forcefield which we know exactly at the sampling points including its derivative, because this is exactly the scenario for which it was developed. I.e. 1D Hermite cubic spline on an interval $(x_1 ... x_2)$ is linar combindation of cubic polynomical where the weight coefficients of each basis polynominal are values of the inrerpolated function and the boundary of interpolated interval $(w_1 = f(x_1)$ and $w_2 = f(x_2) )$, and value of the derivative at the same point ( $w_3 = \partial f(x_2)/\partial x$ and $w_2 = \partial f(x_2)/\partial x$ ).

    • This require to explicitly store not only the potneital energy but also the forces (i.e. we will not save any memory).

      • This becomes the problem in 2D and 3D since where we need to store not only the derivatives along one dimension $\partial f(x,y,z)/\partial x$ but also the mixed derivatives like $\partial \partial f(x,y,z)/(\partial x \partial y)$ and $\partial \partial \partial f(x,y,z)/(\partial x \partial y \partial z) $

    • There are however variants of Hermite spline, (such as Natural Hermite cubic spline) where we do not store the derivatives, but we approximate the derivatives by finite differences of the function value in neighboring grid points.

      • In case of 3D (tricubic) Natural Hermite spline interpolation this means we read potential energy values from 4x4x4 stencil neighborhood of the interpolated point from the global grid (i.e. 64 floating point values). We do not need to store any derivatives.
      • However the forces obtained by hermite spline with derivatives approximate by finite differences are less smooth than canonical Hermite spline with explicit analytical deivatives in control points. In particular, the piecewise polynominals for forces are parabolas (as can be expected a parabola is derivative of cubic polynominal), but with sharp kinks at the grid-points which does not match the analytical forces.
        • Therefore based on our numerical experiments we do not consider this stratedy (using appriximate derivatives by finite diferences) very atractive.
        • In fact our numerical experiments shows these approximate Natural Hermite cubic splines have higher average error in the force than the liner interpolation of forces. Nevertheless, cubic splines are ate least conservative (as they are true derivative of the potential), which gives them significant advantage over triliner interpolation, even if the accuracy is lower.

Cubic B-spline

  • B-splines are typically used for approximation rather then interpolation. This suggest they are sub-optimal choice for our case where we know exact value of our function at the grid points.
  • Nevertheless B-splines are generaly more smooth than Hermite splines because even the derivative of cubic B-spline is C2 smooth (see image )
  • Another advantage of B-splines (over Hermite spline) is that we need to store just one grid weight coeficients (each grid point store weight coefficient atributed to bell-shaped cubic spline centered at that grid point). Therefore we do not need to store derivatives explicitly.
  • The disadvantage of using B-splines is that we cannot use the known function values (i.e. potential energy) at the grid points directly but instead we need to fit the values of the weight coeficients attributed to bell-splines to reproduce them.
    • This means we need to implement of fitting algorithm which optimize values of ~1 million parameters (grid weights) depending on ~1 million sampled function values. Fortunately linear system corespoding to this linear fit is very sparse matrix as each bell-spline affects function values only in nearest neightbor grid points (3x3x3 stencil = 27 points in 3D).
      • Neverthless solving even this sparse system may be rather demainding and problematic especially when it is stiff, which is the case when fit e.g. Lenard-Jones potential which is very flat far from atoms, but has shar spike (almost singularity) near the atoms.
        • Currently I implemented rather only rather inefficient dynamicall dampled fitting algorithm, which requies several thousands interations to convege and even than it is sometimes prone to instability.

Code

References