The source is here.

Points of interest

use of a binary operation to apply a complex operation to two functions
use of asymptotic analysis to ensure good behavior in presence of numerical noise

Background

In nuclear physics density functional theory it is necessary to compute functions of the form

$U(r) = \frac{\Delta^2 (r)}{\rho^{2/3} (r)}$

The functions $\Delta$ and $\rho$ are both expected to go to zero at large as is the ratio (i.e., $\Delta$ goes to zero faster than $\rho$ ). Moreover, $\rho$ should everywhere be positive and is not expected to be zero in the interior region (for ground states only?).

Implementation

The first problem is how to compose this operation inside MADNESS. One could square , use unaryop() to compute the negative fractional power of $\rho()$ , and then multiply the two. With care (see below) this should work. Easier, faster, and more accurate is to do all of the above at once. This is accomplished with a binary operation that acts upon the input function values and returns the result.

The second and most significant problem is numerical noise in the functions that can lead to $\rho$ being zero and even negative while is non-zero. However, the expected asymptotics tell us that if either or $\rho$ are so small that noise is dominating, that the result of the binary operation should be zero. [Aside. To accomplish the same using a unary operation the operation that computes $rho^{-2/3}$ should return zero if is small. But this precludes us from simultaneously using information about the size of $\Delta$ and does not ensure that both are computed at the same level of refinement.]

Analysis is necessary. The threshold at which to screen values to their asymptotic form depends on the problem, the accuracy of computation, and possibly the box size. In this problem we choose

$\Delta(r) = exp(- | r | )$

and

$\rho(r) = exp(- 2 | r | ) = \Delta^2(r)$

Thus, the exact result is

$U(r) = \frac{\Delta^2 (r)}{\rho^{2/3} (r)} = exp( - 2 | r | / 3)$

Note that the result has a smaller exponent than the two input functions and is therefore significant at a longer range. Since we cannot generate information we do not have, once the input functions degenerate into numerical noise we must expect that the ratio is also just noise. In the physical application, the potential is applied to another function that is also decaying expoentially, which makes small noise at long range not significant. By screening to the physically expected value of zero we therefore ensure correct physics.