PrimitiveOverlapIntegralEngine.hpp

Go to the documentation of this file.

// This file is part of GQCG-GQCP.
//
// Copyright (C) 2017-2020  the GQCG developers
//
// GQCG-GQCP is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// GQCG-GQCP is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with GQCG-GQCP.  If not, see <http://www.gnu.org/licenses/>.
 
#pragma once
 
#include "Basis/Integrals/Primitive/BaseScalarPrimitiveIntegralEngine.hpp"
#include "Basis/Integrals/Primitive/McMurchieDavidsonCoefficient.hpp"
#include "Basis/ScalarBasis/GTOShell.hpp"
#include "Basis/ScalarBasis/LondonGTOShell.hpp"
#include "Mathematical/Functions/CartesianGTO.hpp"
#include "Mathematical/Functions/LondonCartesianGTO.hpp"
#include "Operator/FirstQuantized/OverlapOperator.hpp"
#include "Utilities/complex.hpp"
#include "Utilities/type_traits.hpp"
 
#include <boost/math/constants/constants.hpp>
 
 
namespace GQCP {
 
 
template <typename _Shell>
class PrimitiveOverlapIntegralEngine:
    public BaseScalarPrimitiveIntegralEngine {
public:
    // The type of shell that this integral engine is related to.
    using Shell = _Shell;
 
    // The type of primitive that underlies the type of shell.
    using Primitive = typename Shell::Primitive;
 
    // The scalar representation of an overlap integral.
    using IntegralScalar = product_t<OverlapOperator::Scalar, typename Primitive::OutputType>;
 
 
public:
    /*
     *  MARK: CartesianGTO integrals
     */
 
    template <typename Z = Shell>
    enable_if_t<std::is_same<Z, GTOShell>::value, IntegralScalar> calculate(const CartesianGTO& left, const CartesianGTO& right) {
 
        // The 3D integral is separable in three 1D integrals.
        IntegralScalar primitive_integral {1.0};
        for (const auto& direction : {GQCP::CartesianDirection::x, GQCP::CartesianDirection::y, GQCP::CartesianDirection::z}) {
            const auto i = left.cartesianExponents().value(direction);
            const auto j = right.cartesianExponents().value(direction);
 
            primitive_integral *= this->calculate1D(left.gaussianExponent(), left.center()(direction), i, right.gaussianExponent(), right.center()(direction), j);
        }
 
        return primitive_integral;
    }
 
 
    template <typename Z = Shell>
    enable_if_t<std::is_same<Z, GTOShell>::value, IntegralScalar> calculate1D(const double a, const double K, const int i, const double b, const double L, const int j) {
 
        // Negative Cartesian exponents should be ignored: the correct value for the corresponding integral is 0.
        if ((i < 0) || (j < 0)) {
            return 0.0;
        }
 
        // Use the McMurchie-Davidson recursion to calculate the overlap integral.
        const auto p = a + b;
        const McMurchieDavidsonCoefficient E {K, a, L, b};
 
        return std::pow(boost::math::constants::pi<double>() / p, 0.5) * E(i, j, 0);
    }
 
 
    /*
     *  MARK: LondonCartesianGTO integrals
     */
 
    template <typename Z = Shell>
    enable_if_t<std::is_same<Z, LondonGTOShell>::value, IntegralScalar> calculate(const LondonCartesianGTO& left, const LondonCartesianGTO& right) {
 
        const Vector<double, 3> k1 = right.kVector() - left.kVector();  // The k-vector of the London overlap distribution.
 
        // The 3D integral is separable in three 1D integrals.
        IntegralScalar primitive_integral {1.0};
        for (const auto& direction : {GQCP::CartesianDirection::x, GQCP::CartesianDirection::y, GQCP::CartesianDirection::z}) {
            const auto a = left.cartesianGTO().gaussianExponent();
            const auto K = left.cartesianGTO().center()(direction);
            const auto i = left.cartesianGTO().cartesianExponents().value(direction);
 
            const auto b = right.cartesianGTO().gaussianExponent();
            const auto L = right.cartesianGTO().center()(direction);
            const auto j = right.cartesianGTO().cartesianExponents().value(direction);
 
            const auto k1_component = k1(direction);
            primitive_integral *= this->calculate1D(k1_component, a, K, i, b, L, j);
        }
 
        return primitive_integral;
    }
 
 
    template <typename Z = Shell>
    enable_if_t<std::is_same<Z, LondonGTOShell>::value, IntegralScalar> calculate1D(const complex k1, const double a, const double K, const int i, const double b, const double L, const int j) {
 
        // Negative Cartesian exponents should be ignored: the correct value for the corresponding integral is 0.
        if ((i < 0) || (j < 0)) {
            return 0.0;
        }
 
        using namespace GQCP::literals;
 
 
        // Use the McMurchie-Davidson recursion to calculate the overlap integral.
        const auto p = a + b;
        const McMurchieDavidsonCoefficient E {K, a, L, b};
        const auto P = E.centerOfMass();
 
        IntegralScalar integral {0.0};
        for (int t = 0; t <= i + j; t++) {
            integral += E(i, j, t) *
                        std::pow(-1.0_ii * k1, t) *
                        std::pow(boost::math::constants::pi<double>() / p, 0.5) *
                        std::exp(-1.0_ii * k1 * P) *
                        std::exp(-std::pow(k1, 2) / (4 * p));
        }
 
        return integral;
    }
};
 
 
}  // namespace GQCP