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 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

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

Parameters:

beta – [in] Exponent, \( e^{-\beta r^2} \)
alpha – [in] Coefficient, \( \alpha e^{-r^2} \)
pos – [in] Position \( (x - pos[0]), (y - pos[1]), ... \)
pow – [in] Monomial power, \( x^{pow[0]}, y^{pow[1]}, ... \)

Returns:

New GaussFunc object

double calcCoulombEnergy(const GaussFunc<D> &rhs) const¶

Compute Coulomb repulsion energy between two GaussFuncs.

Note

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

Parameters:

this – [in] Left hand GaussFunc
rhs – [in] Right hand GaussFunc

Returns:

Coulomb energy

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

Parameters:: r – [in] Cartesian coordinate
Returns:: Function value in a point

virtual GaussPoly<D> differentiate(int dir) const override¶

Compute analytic derivative of Gaussian.

Parameters:: dir – [in] Cartesian direction of derivative
Returns:: New GaussPoly

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

Multiply two GaussFuncs.

Parameters:

this – [in] Left hand side of multiply
rhs – [in] Right hand side of multiply

Returns:

New GaussPoly

GaussFunc<D> mult(double c)¶

Multiply GaussFunc by scalar.

Parameters:: c – [in] Scalar to multiply
Returns:: New GaussFunc

GaussExp<D> periodify(const std::array<double, D> &period, double nStdDev = 4.0) const¶

Generates a GaussExp that is semi-periodic around a unit-cell.

nStdDev = 1, 2, 3 and 4 ensures atleast 68.27%, 95.45%, 99.73% and 99.99% of the integral is conserved with respect to the integration limits.

Parameters:

period – [in] The period of the unit cell
nStdDev – [in] Number of standard diviations covered in each direction. Default 4.0

Returns:

Semi-periodic version of a Gaussian around a unit-cell

inline void normalize()¶: Rescale function by its norm \( ||f||^{-1} \).

template<int D> class 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 = {})¶

Parameters:

beta – [in] Exponent, \( e^{-\beta r^2} \)
alpha – [in] Coefficient, \( \alpha e^{-r^2} \)
pos – [in] Position \( (x - pos[0]), (y - pos[1]), ... \)
pow – [in] Max polynomial degree, \( P_0(x), P_1(y), ... \)

Returns:

New GaussPoly object

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

Parameters:: r – [in] Cartesian coordinate
Returns:: Function value in a point

virtual GaussPoly differentiate(int dir) const override¶

Compute analytic derivative of Gaussian.

Parameters:: dir – [in] Cartesian direction of derivative
Returns:: New GaussPoly

GaussPoly<D> mult(double c)¶

Multiply GaussPoly by scalar.

Parameters:: c – [in] Scalar to multiply
Returns:: New GaussPoly

void setPoly(int d, Polynomial &poly)¶

Set polynomial in given dimension.

Parameters:

d – [in] Cartesian direction
poly – [in] Polynomial to set

GaussExp<D> periodify(const std::array<double, D> &period, double nStdDev = 4.0) const¶

Generates a GaussExp that is semi-periodic around a unit-cell.

nStdDev = 1, 2, 3 and 4 ensures atleast 68.27%, 95.45%, 99.73% and 99.99% of the integral is conserved with respect to the integration limits.

Parameters:

period – [in] The period of the unit cell
nStdDev – [in] Number of standard diviations covered in each direction. Default 4.0

Returns:

Semi-periodic version of a Gaussian around a unit-cell

inline void normalize()¶: Rescale function by its norm \( ||f||^{-1} \).

template<int D> class 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() const¶

Note

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

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

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

Parameters:: r – [in] Cartesian coordinate
Returns:: Function value in a point

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;

// Add polynomials to gauss_poly
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