MWOperators¶

The MW operators discussed in this chapter is available to the application program by including:

#include "MRCPP/MWOperators"

template<int D> class MWOperator¶

Fixme.

Fixme

Subclassed by mrcpp::ConvolutionOperator< 3 >, mrcpp::ConvolutionOperator< D >, mrcpp::DerivativeOperator< D >

ConvolutionOperator¶

Note

The convolution operators have separate precision parameters for their construction and application. The build_prec argument to the operator constructors will affect e.g. the number of terms in the separated representations of the Poisson/Helmholtz approximations, as well as the operator bandwidth. The apply_prec argument to the apply function relates only to the adaptive construction of the output function, based on a wavelet norm error estimate.

template<int D> class IdentityConvolution : public mrcpp::ConvolutionOperator<D>¶

Convolution with an identity kernel.

The identity kernel (Dirac’s delta function) is approximated by a narrow Gaussian function: \( I(r-r') = \delta(r-r') \approx \alpha e^{-\beta (r-r')^2} \)

Public Functions

IdentityConvolution(const MultiResolutionAnalysis<D> &mra, double prec)¶

Constructor of the IdentityConvolution object.

This will project a kernel of a single gaussian with exponent sqrt(10/build_prec).

Parameters:

mra – [in] Which MRA the operator is defined
prec – [in] Build precision, closeness to delta function

Returns:

New IdentityConvolution object

IdentityConvolution(const MultiResolutionAnalysis<D> &mra, double prec, int root, int reach = 1)¶

Constructor of the IdentityConvolution object in case of Periodic Boundary Conditions (PBC)

This will project a kernel of a single gaussian with exponent sqrt(10/build_prec). This version of the constructor is used for calculations within periodic boundary conditions (PBC). The root parameter is the coarsest negative scale at wich the operator is applied. The reach parameter is the bandwidth of the operator at the root scale. For details see MWOperator

Parameters:

mra – [in] Which MRA the operator is defined
prec – [in] Build precision, closeness to delta function
root – [in] root scale of operator.
reach – [in] width at root scale (applies to periodic boundary conditions)

Returns:

New IdentityConvolution object

template<int D> class DerivativeConvolution : public mrcpp::ConvolutionOperator<D>¶

Convolution with a derivative kernel.

Derivative operator written as a convolution. The derivative kernel (derivative of Dirac’s delta function) is approximated by the derivative of a narrow Gaussian function: \( D^x(r-r') = \frac{d}{dx}\delta(r-r') \approx \frac{d}{dx} \alpha e^{-\beta (r-r')^2} \)

NOTE: This is not the recommended derivative operator for practial calculations, it’s a proof-of-concept operator. Use the ABGVOperator for “cuspy” functions and the BSOperator for smooth functions.

Public Functions

DerivativeConvolution(const MultiResolutionAnalysis<D> &mra, double prec)¶

This will project a kernel of a single differentiated gaussian with exponent sqrt(10/build_prec).

Parameters:

mra – [in] Which MRA the operator is defined
pr – [in] Build precision, closeness to delta function

Returns:

New DerivativeConvolution object

class PoissonOperator : public mrcpp::ConvolutionOperator<3>¶

Convolution with the Poisson Green’s function kernel.

The Poisson kernel is approximated as a sum of Gaussian functions in order to allow for separated application of the operator in the Cartesian directions: \( P(r-r') = \frac{1}{|r-r'|} \approx \sum_m^M \alpha_m e^{-\beta_m (r-r')^2} \)

Public Functions

PoissonOperator(const MultiResolutionAnalysis<3> &mra, double prec)¶

This will construct a gaussian expansion to approximate 1/r, and project each term into a one-dimensional MW operator. Subsequent application of this operator will apply each of the terms to the input function in all Cartesian directions.

Parameters:

mra – [in] Which MRA the operator is defined
pr – [in] Build precision, closeness to 1/r

Returns:

New PoissonOperator object

class HelmholtzOperator : public mrcpp::ConvolutionOperator<3>¶

Convolution with the Helmholtz Green’s function kernel.

The Helmholtz kernel is approximated as a sum of gaussian functions in order to allow for separated application of the operator in the Cartesian directions: \( H(r-r') = \frac{e^{-\mu|r-r'|}}{|r-r'|} \approx \sum_m^M \alpha_m e^{-\beta_m (r-r')^2} \)

Public Functions

HelmholtzOperator(const MultiResolutionAnalysis<3> &mra, double m, double prec)¶

This will construct a gaussian expansion to approximate exp(-mu*r)/r, and project each term into a one-dimensional MW operator. Subsequent application of this operator will apply each of the terms to the input function in all Cartesian directions.

Parameters:

mra – [in] Which MRA the operator is defined
m – [in] Exponential parameter of the operator
pr – [in] Build precision, closeness to exp(-mu*r)/r

Returns:

New HelmholtzOperator object

Warning

doxygenfunction: Unable to resolve function “mrcpp::apply” with arguments “(double, FunctionTree<D>&, ConvolutionOperator<D>&, FunctionTree<D>&, int, bool)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template void apply (CompFunction< 3 > &out, DerivativeOperator< 3 > &oper, CompFunction< 3 > &inp, int dir=-1, const ComplexDouble(*metric)[4]) ———^

Warning

doxygenfunction: Unable to resolve function “mrcpp::apply” with arguments “(double, FunctionTree<D>&, ConvolutionOperator<D>&, FunctionTree<D>&, FunctionTreeVector<D>&, int, bool)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template void apply (CompFunction< 3 > &out, DerivativeOperator< 3 > &oper, CompFunction< 3 > &inp, int dir=-1, const ComplexDouble(*metric)[4]) ———^

DerivativeOperators¶

Note

The derivative operators have clearly defined requirements on the output grid structure, based on the grid of the input function. This means that there is no real grid adaptivity, and thus no precision parameter is needed for the application of such an operator.

template<int D> class ABGVOperator : public mrcpp::DerivativeOperator<D>¶

Derivative operator as defined by Alpert, Beylkin, Ginez and Vozovoi, J Comp Phys 182, 149-190 (2002).

NOTE: This is the recommended derivative operator for “cuspy” or discontinuous functions. The BSOperator is recommended for smooth functions.

Public Functions

ABGVOperator(const MultiResolutionAnalysis<D> &mra, double a, double b)¶

Boundary parameters correspond to:

a=0.0 b=0.0: Strictly local “center” difference
a=0.5 b=0.5: Semi-local central difference
a=1.0 b=0.0: Semi-local forward difference
a=0.0 b=1.0: Semi-local backward difference

Parameters:

mra – [in] Which MRA the operator is defined
a – [in] Left boundary condition
b – [in] Right boundary condition

Returns:

New ABGVOperator object

template<int D> class PHOperator : public mrcpp::DerivativeOperator<D>¶

Derivative operator based on the smoothing derivative of Pavel Holoborodko .

NOTE: This is not the recommended derivative operator for practial calculations, it’s a proof-of-concept operator. Use the ABGVOperator for “cuspy” functions and the BSOperator for smooth functions.

Public Functions

PHOperator(const MultiResolutionAnalysis<D> &mra, int order)¶

Parameters:

mra – [in] Which MRA the operator is defined
order – [in] Derivative order, defined for 1 and 2

Returns:

New PHOperator object

template<int D> class BSOperator : public mrcpp::DerivativeOperator<D>¶

B-spline derivative operator as defined by Anderson etal, J Comp Phys X 4, 100033 (2019).

NOTE: This is the recommended derivative operator only for smooth functions. Use the ABGVOperator if the function has known cusps or discontinuities.

Public Functions

explicit BSOperator(const MultiResolutionAnalysis<D> &mra, int order)¶

Parameters:

mra – [in] Which MRA the operator is defined
order – [in] Derivative order, defined for 1, 2 and 3

Returns:

New BSOperator object

Warning

doxygenfunction: Unable to resolve function “mrcpp::apply” with arguments “(FunctionTree<D>&, DerivativeOperator<D>&, FunctionTree<D>&, int)”. Candidate function could not be parsed. Parsing error is Invalid C++ declaration: Expected ‘<’ after ‘template’ [error at 9] template void apply (CompFunction< 3 > &out, DerivativeOperator< 3 > &oper, CompFunction< 3 > &inp, int dir=-1, const ComplexDouble(*metric)[4]) ———^

Warning

doxygenfunction: Unable to resolve function “mrcpp::divergence” with arguments (FunctionTree<D>&, DerivativeOperator<D>&, FunctionTreeVector<D>&) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches:

- template<int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D, T> *inp, const ComplexDouble (*metric)[4])
- template<int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T>*> *inp, const ComplexDouble (*metric)[4])
- template<int D, typename T> void divergence(FunctionTree<D, T> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D, T> &inp)
- template<int D, typename T> void divergence(FunctionTree<D, T> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T>*> &inp)

Warning

doxygenfunction: Unable to resolve function “mrcpp::gradient” with arguments (DerivativeOperator<D>&, FunctionTree<D>&) in doxygen xml output for project “MRCPP” from directory: _build/xml. Potential matches:

- std::vector<CompFunction<3>*> gradient(DerivativeOperator<3> &oper, CompFunction<3> &inp, const ComplexDouble (*metric)[4])
- template<int D, typename T> FunctionTreeVector<D, T> gradient(DerivativeOperator<D> &oper, FunctionTree<D, T> &inp)

Examples¶

PoissonOperator¶

The electrostatic potential \(g\) arising from a charge distribution \(f\) are related through the Poisson equation

\[-\nabla^2 g(r) = f(r)\]

This equation can be solved with respect to the potential by inverting the differential operator into the Green’s function integral convolution operator

\[g(r) = \int \frac{1}{4\pi\|r-r'\|} f(r') dr'\]

This operator is available in the MW representation, and can be solved with arbitrary (finite) precision in linear complexity with respect to system size. Given an arbitrary charge dirtribution f_tree in the MW representation, the potential is computed in the following way:

double apply_prec;                              // Precision for operator application
double build_prec;                              // Precision for operator construction

mrcpp::PoissonOperator P(MRA, build_prec);      // MW representation of Poisson operator
mrcpp::FunctionTree<3> f_tree(MRA);             // Input function
mrcpp::FunctionTree<3> g_tree(MRA);             // Output function

mrcpp::apply(apply_prec, g_tree, P, f_tree);    // Apply operator adaptively

The Coulomb self-interaction energy can now be computed as the dot product:

double E = mrcpp::dot(g_tree, f_tree);

HelmholtzOperator¶

The Helmholtz operator is a generalization of the Poisson operator and is given as the integral convolution

\[g(r) = \int \frac{e^{-\mu\|r-r'\|}}{4\pi\|r-r'\|} f(r') dr'\]

The operator is the inverse of the shifted Laplacian

\[\big[-\nabla^2 + \mu^2 \big] g(r) = f(r)\]

and appears e.g. when solving the SCF equations. The construction and application is similar to the Poisson operator, with an extra argument for the \(\mu\) parameter

double apply_prec;                              // Precision for operator application
double build_prec;                              // Precision for operator construction
double mu;                                      // Must be a positive real number

mrcpp::HelmholtzOperator H(MRA, mu, build_prec);// MW representation of Helmholtz operator
mrcpp::FunctionTree<3> f_tree(MRA);             // Input function
mrcpp::FunctionTree<3> g_tree(MRA);             // Output function

mrcpp::apply(apply_prec, g_tree, H, f_tree);    // Apply operator adaptively

ABGVOperator¶

The ABGV (Alpert, Beylkin, Gines, Vozovoi) derivative operator is initialized with two parameters \(a\) and \(b\) accounting for the boundary conditions between adjacent nodes, see Alpert et al.

double a = 0.0, b = 0.0;                        // Boundary conditions for operator
mrcpp::ABGVOperator<3> D(MRA, a, b);            // MW derivative operator
mrcpp::FunctionTree<3> f(MRA);                  // Input function
mrcpp::FunctionTree<3> f_x(MRA);                // Output function
mrcpp::FunctionTree<3> f_y(MRA);                // Output function
mrcpp::FunctionTree<3> f_z(MRA);                // Output function

mrcpp::apply(f_x, D, f, 0);                     // Operator application in x direction
mrcpp::apply(f_y, D, f, 1);                     // Operator application in y direction
mrcpp::apply(f_z, D, f, 2);                     // Operator application in z direction

The tree structure of the output function will depend on the choice of parameters \(a\) and \(b\): if both are zero, the output grid will be identical to the input grid; otherwise the grid will be widened by one node (on each side) in the direction of application.

PHOperator¶

The PH derivative operator is based on the noise reducing derivative of Pavel Holoborodko. This operator is also available as a direct second derivative.

mrcpp::PHOperator<3> D1(MRA, 1);                // MW 1st derivative operator
mrcpp::PHOperator<3> D2(MRA, 2);                // MW 2nd derivative operator
mrcpp::FunctionTree<3> f(MRA);                  // Input function
mrcpp::FunctionTree<3> f_x(MRA);                // Output function
mrcpp::FunctionTree<3> f_xx(MRA);               // Output function

mrcpp::apply(f_x, D1, f, 0);                    // Operator application in x direction
mrcpp::apply(f_xx, D2, f, 0);                   // Operator application in x direction

Special thanks to Prof. Robert J. Harrison (Stony Brook University) for sharing the operator coefficients.

BSOperator¶

The BS derivative operator is based on a pre-projection onto B-splines in order to remove the discontinuities in the MW basis, see Anderson et al. This operator is also available as a direct second and third derivative.

mrcpp::BSOperator<3> D1(MRA, 1);                // MW 1st derivative operator
mrcpp::BSOperator<3> D2(MRA, 2);                // MW 2nd derivative operator
mrcpp::BSOperator<3> D3(MRA, 3);                // MW 3nd derivative operator
mrcpp::FunctionTree<3> f(MRA);                  // Input function
mrcpp::FunctionTree<3> f_x(MRA);                // Output function
mrcpp::FunctionTree<3> f_yy(MRA);               // Output function
mrcpp::FunctionTree<3> f_zzz(MRA);              // Output function

mrcpp::apply(f_x, D1, f, 0);                    // Operator application in x direction
mrcpp::apply(f_yy, D2, f, 1);                   // Operator application in x direction
mrcpp::apply(f_zzz, D3, f, 2);                  // Operator application in x direction