# Gaussians¶

MRCPP provides some simple features for analytic Gaussian functions. These are meant to be used as a starting point for MW computations, they are not meant for heavy analytical computation, like GTO basis sets. The Gaussian features are available by including:

#include "MRCPP/Gaussians"

template<int D>
class mrcpp::GaussFunc : public mrcpp::Gaussian<D>

Gaussian function in D dimensions with a simple monomial in front.

• Monodimensional Gaussian (GaussFunc<1>):

$$g(x) = \alpha (x-x_0)^a e^{-\beta (x-x_0)^2}$$

• Multidimensional Gaussian (GaussFunc<D>):

$$G(x) = \prod_{d=1}^D g^d(x^d)$$

Public Functions

GaussFunc(double beta, double alpha, const Coord<D> &pos = {}, const std::array<int, D> &pow = {})

Return

New GaussFunc object

Parameters
• [in] beta: Exponent, $$e^{-\beta r^2}$$

• [in] alpha: Coefficient, $$\alpha e^{-r^2}$$

• [in] pos: Position $$(x - pos), (y - pos), ...$$

• [in] pow: Monomial power, $$x^{pow}, y^{pow}, ...$$

double calcCoulombEnergy(GaussFunc<D> &rhs)

Compute Coulomb repulsion energy between two GaussFuncs.

Return

Coulomb energy

Note

Both Gaussians must be normalized to unit charge $$\alpha = (\beta/\pi)^{D/2}$$ for this to be correct!

Parameters
• [in] this: Left hand GaussFunc

• [in] rhs: Right hand GaussFunc

double calcOverlap(GaussFunc<D> &rhs) override

Compute analytic overlap between two Gaussians.

Return

Overlap integral

Parameters
• [in] this: Left hand Gaussian

• [in] rhs: Right hand Gaussian

double calcOverlap(GaussPoly<D> &rhs) override

Compute analytic overlap between two Gaussians.

Return

Overlap integral

Parameters
• [in] this: Left hand Gaussian

• [in] rhs: Right hand Gaussian

GaussPoly<D> differentiate(int dir) override

Compute analytic derivative of Gaussian.

Return

New GaussPoly

Parameters
• [in] dir: Cartesian direction of derivative

GaussPoly<D> mult(const GaussFunc<D> &rhs)

Multiply two GaussFuncs.

Return

New GaussPoly

Parameters
• [in] this: Left hand side of multiply

• [in] rhs: Right hand side of multiply

GaussFunc<D> mult(double c)

Multiply GaussFunc by scalar.

Return

New GaussFunc

Parameters
• [in] c: Scalar to multiply

double evalf(const Coord<D> &r) const

Return

Function value in a point

Parameters
• [in] r: Cartesian coordinate

double getSquareNorm()

Return

Squared L2 norm of function $$||f||^2$$

void normalize()

Rescale function by its norm $$||f||^{-1}$$.

template<int D>
class mrcpp::GaussPoly : public mrcpp::Gaussian<D>

Gaussian function in D dimensions with a general polynomial in front.

• Monodimensional Gaussian (GaussPoly<1>):

$$g(x) = \alpha P(x-x_0) e^{-\beta (x-x_0)^2}$$

• Multidimensional Gaussian (GaussFunc<D>):

$$G(x) = \prod_{d=1}^D g^d(x^d)$$

Public Functions

GaussPoly(double alpha = 0.0, double coef = 1.0, const Coord<D> &pos = {}, const std::array<int, D> &power = {})

Return

New GaussPoly object

Parameters
• [in] beta: Exponent, $$e^{-\beta r^2}$$

• [in] alpha: Coefficient, $$\alpha e^{-r^2}$$

• [in] pos: Position $$(x - pos), (y - pos), ...$$

• [in] pow: Max polynomial degree, $$P_0(x), P_1(y), ...$$

double calcOverlap(GaussFunc<D> &b) override

Compute analytic overlap between two Gaussians.

Return

Overlap integral

Parameters
• [in] this: Left hand Gaussian

• [in] rhs: Right hand Gaussian

double calcOverlap(GaussPoly<D> &b) override

Compute analytic overlap between two Gaussians.

Return

Overlap integral

Parameters
• [in] this: Left hand Gaussian

• [in] rhs: Right hand Gaussian

GaussPoly differentiate(int dir) override

Compute analytic derivative of Gaussian.

Return

New GaussPoly

Parameters
• [in] dir: Cartesian direction of derivative

GaussPoly<D> mult(double c)

Multiply GaussPoly by scalar.

Return

New GaussPoly

Parameters
• [in] c: Scalar to multiply

void setPoly(int d, Polynomial &poly)

Set polynomial in given dimension.

Parameters
• [in] d: Cartesian direction

• [in] poly: Polynomial to set

double evalf(const Coord<D> &r) const

Return

Function value in a point

Parameters
• [in] r: Cartesian coordinate

double getSquareNorm()

Return

Squared L2 norm of function $$||f||^2$$

void normalize()

Rescale function by its norm $$||f||^{-1}$$.

template<int D>
class mrcpp::GaussExp : public mrcpp::RepresentableFunction<D>

Gaussian expansion in D dimensions.

• Monodimensional Gaussian expansion:

$$g(x) = \sum_{m=1}^M g_m(x) = \sum_{m=1}^M \alpha_m e^{-\beta (x-x^0)^2}$$

• Multidimensional Gaussian expansion:

$$G(x) = \sum_{m=1}^M G_m(x) = \sum_{m=1}^M \prod_{d=1}^D g_m^d(x^d)$$

Public Functions

double calcCoulombEnergy()

Return

Coulomb repulsion energy between all pairs in GaussExp, including self-interaction

Note

Each Gaussian must be normalized to unit charge $$c = (\alpha/\pi)^{D/2}$$ for this to be correct!

double evalf(const Coord<D> &r) const override

Return

Function value in a point

Parameters
• [in] r: Cartesian coordinate

void append(const Gaussian<D> &g)

Append Gaussian to expansion.

void append(const GaussExp<D> &g)

Append GaussExp to expansion.

## Examples¶

A GaussFunc is a simple D-dimensional Gaussian function with a Cartesian monomial in front, e.g. in 3D:

$f(r) = \alpha (x-x_0)^a (y-y_0)^b (z-z_0)^c e^{-\beta \|r-r_0\|^2}$
double alpha, beta;
std::array<int, 3> pow = {a, b, c};
mrcpp::Coord<3> pos = {x_0, y_0, z_0};
mrcpp::GaussFunc<3> gauss(beta, alpha, pos, pow);

double E = gauss.calcCoulombEnergy(gauss);              // Analytical energy


This Gaussian function can be used to build an empty grid based on the position and exponent. The grid will then be refined close to the center of the Gaussian, with deeper refinement for higher exponents (steeper function):

mrcpp::FunctionTree<3> g_tree(MRA);
mrcpp::build_grid(g_tree, gauss);                       // Build empty grid
mrcpp::project(prec, g_tree, gauss);                    // Project Gaussian


GaussPoly is a generalization of the GaussFunc, where there is an arbitrary polynomial in front of the exponential

$f(r) = \alpha P(r-r_0) e^{-\beta \|r-r_0\|^2}$

For instance, the following function can be constructed

$f(r) = \alpha (a_x + b_x x + c_x x^2) (a_y + b_y y + c_y y^2) (a_z + b_z z + c_z z^2)e^{-\beta \|r-r_0\|^2}$
auto gauss_poly = GaussPoly<D>(beta, alpha, pos, pow);

// Create polynomial in x, y and z direction
auto pol_x = Polynomial(2); // 2 is the degree of the polynomial
pol_x.getCoefs() << a_x, b_x, c_x;
auto pol_y = Polynomial(2);
pol_y.getCoefs() << a_y, b_y, c_y;
auto pol_z = Polynomial(2);
pol_z.getCoefs() << a_z, b_z, c_z;

guass_poly.setPoly(0, pol_x);
guass_poly.setPoly(1, pol_y);
guass_poly.setPoly(2, pol_z);


A GaussExp is a collection of Gaussians in the form

$G(r) = \sum_i c_i g_i(r)$

where $$g_i$$ can be either GaussFunc or GaussPoly

$g_i(r) = \alpha_i P_i(r-r_i)e^{-\beta_i\|r-r_i\|^2}$

Individual Gaussian functions can be appended to the GaussExp and treated as a single function:

mrcpp::GaussExp<3> g_exp;                               // Empty Gaussian expansion
for (int i = 0; i < N; i++) {
double alpha_i, beta_i;                             // Individual parameters
std::array<int, 3> pow_i;                           // Individual parameters
std::array<double, 3> pos_i;                        // Individual parameters
mrcpp::GaussFunc<3> gauss_i(beta_i, alpha_i, pos_i, pow_i);
g_exp.append(gauss_i);                              // Append Gaussian to expansion
}
mrcpp::project(prec, tree, g_exp);                      // Project full expansion