radial.h

Go to the documentation of this file.

//  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  //
//                                                                                   //
//                       / /   / /    __    \ \  / /                                 //
//                      / /__ / /   / _ \    \ \/ /                                  //
//                     /  ___  /   | |/_/    / /\ \                                  //
//                    / /   / /    \_\      / /  \ \                                 //
//                                                                                   //
//                                                                                   //
//  Copyright (c) 2017, Jakub Benda, Charles University in Prague                    //
//                                                                                   //
// MIT License:                                                                      //
//                                                                                   //
//  Permission is hereby granted, free of charge, to any person obtaining a          //
// copy of this software and associated documentation files (the "Software"),        //
// to deal in the Software without restriction, including without limitation         //
// the rights to use, copy, modify, merge, publish, distribute, sublicense,          //
// and/or sell copies of the Software, and to permit persons to whom the             //
// Software is furnished to do so, subject to the following conditions:              //
//                                                                                   //
//  The above copyright notice and this permission notice shall be included          //
// in all copies or substantial portions of the Software.                            //
//                                                                                   //
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS          //
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE       //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, //
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF         //
// OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.  //
//                                                                                   //
//  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  //

#ifndef HEX_ECS_RADIAL_H
#define HEX_ECS_RADIAL_H

// --------------------------------------------------------------------------------- //

#include <cmath>
#include <functional>
#include <vector>

// --------------------------------------------------------------------------------- //

#include "hex-arrays.h"
#include "hex-symbandmatrix.h"
#include "hex-special.h"

// --------------------------------------------------------------------------------- //

#include "bspline.h"
#include "gauss.h"
#include "io.h"
#include "parallel.h"

// --------------------------------------------------------------------------------- //

#define EXPANSION_QUADRATURE_POINTS 20

// --------------------------------------------------------------------------------- //

// Complex -> Complex function.
typedef std::function<Complex(Complex)> CCFunction;

// Complex x Complex -> Complex function.
typedef std::function<Complex(Complex,Complex)> C2CFunction;

// --------------------------------------------------------------------------------- //

class RadialIntegrals
{
    public:
        
        RadialIntegrals
        (
            Bspline const & bspline_x,
            Bspline const & bspline_y,
            int Nlambdas
        );
        
        void setupOneElectronIntegrals (bool sharedscratch = false, bool iammaster = true);
        
        void setupOneElectronIntegrals (Parallel const & par, CommandLine const & cmd);
        
        void setupTwoElectronIntegrals (Parallel const & par, CommandLine const & cmd);
        
        void verbose (bool v) { verbose_ = v; }
        
        int maxlambda () const { return Nlambdas_ - 1; }
        
        // compute overlap matrix of radial function
        Complex computeOverlapMatrixElement_iknot (Bspline const & bspline, GaussLegendre const & g, int i, int j, CCFunction func, int iknot) const;
        Complex computeOverlapMatrixElement (Bspline const & bspline, GaussLegendre const & g, int i, int j, CCFunction func) const;
        SymBandMatrix<Complex> computeOverlapMatrix (Bspline const & bspline, CCFunction func) const;
        
        static CsrMatrix<LU_int_t,Complex> computeOverlapMatrix
        (
            Bspline const & bspline1,
            Bspline const & bspline2,
            Real R_left,
            Real R_right
        );
        
        Complex computeD_iknot
        (
            Bspline const & bspline,
            GaussLegendre const & g,
            int i, int j, int iknot
        ) const;
        
        Complex computeD
        (
            Bspline const & bspline,
            GaussLegendre const & g,
            int i, int j
        ) const;
        
        Complex computeM_iknot
        (
            Bspline const & bspline,
            GaussLegendre const & g,
            int a, int i, int j,
            int iknot,
            Real rmin, Real rmax,
            Real scale
        ) const;
        
        Complex computeM
        (
            Bspline const & bspline,
            GaussLegendre const & g,
            int a, int i, int j,
            Real rmin, Real rmax,
            bool scale = false
        ) const;
        
        Complex computeM
        (
            Bspline const & bspline,
            GaussLegendre const & g,
            int a, int i, int j,
            int begin_knot, int end_knot,
            Real rmin, Real rmax,
            bool scale = false
        ) const;
        
        cArray computeMi
        (
            Bspline const & bspline,
            GaussLegendre const & g,
            int a,
            Real rmin, Real rmax
        ) const;
        
        Complex computeR (int lambda, int i, int j, int k, int l) const;
        
        void diagonalR (int lambda);
        void coupledR (int lambda);
        
        Complex computeRtri
        (
            int lambda,
            Bspline const & xspline, GaussLegendre const & xg,
            int k, int l, int iknotx, Real tx1, Real tx2,
            Bspline const & yspline, GaussLegendre const & yg,
            int m, int n, int iknoty, Real ty1, Real ty2
        ) const;
        
        SymBandMatrix<Complex> calc_R_tr_dia_block
        (
            unsigned lambda,
            int i, int k
        ) const;
        
        void apply_R_matrix
        (
            unsigned int lambda,
            Complex a, const cArrayView src,
            Complex b,       cArrayView dst,
            int x_row_limit = -1,
            int y_row_limit = -1
        ) const;
        
        cArray overlap
        (
            Bspline const & bspline,
            GaussLegendre const & g,
            CCFunction funct
        ) const;
        
        cArray overlap
        (
            Bspline const & bsplinex,
            GaussLegendre const & gx,
            Bspline const & bspliney,
            GaussLegendre const & gy,
            C2CFunction funct
        ) const;
        
        static cArray overlapP
        (
            Bspline const & bspline,
            GaussLegendre const & g,
            Real Z, int n, int l
        );
        
        cArray overlapj
        (
            Bspline const & bspline,
            GaussLegendre const & g,
            int maxell,
            const rArrayView vk,
            bool fast_bessel = false
        ) const;
        
        // Return reference to the B-spline object.
        Bspline const & bspline () const { return bspline_x_; }
        Bspline const & bspline_x () const { return bspline_x_; }
        Bspline const & bspline_y () const { return bspline_y_; }
        
        // Return the Gauss-Legendre integrator object.
        GaussLegendre const & gaussleg () const { return g_x_; }
        GaussLegendre const & gaussleg_x () const { return g_x_; }
        GaussLegendre const & gaussleg_y () const { return g_y_; }
        
        #define OneElectronMatrixAccessors(M) \
            SymBandMatrix<Complex> const & M () const { return M##_x_; } \
            SymBandMatrix<Complex> const & M##_x () const { return M##_x_; } \
            SymBandMatrix<Complex> const & M##_y () const { return M##_y_; } \
            Complex M##_x (std::size_t i, std::size_t j) const { return M##_x_(i,j); } \
            Complex M##_y (std::size_t i, std::size_t j) const { return M##_y_(i,j); } \
            
        // Access the precomputed one-electron overlap matrices.
        OneElectronMatrixAccessors(D)
        OneElectronMatrixAccessors(S)
        OneElectronMatrixAccessors(Mm1)
        OneElectronMatrixAccessors(Mm2)
        
        #define OneElectronMatrixArrayAccessors(M) \
            SymBandMatrix<Complex> const & M (int L) const { return M##_x_[L]; } \
            SymBandMatrix<Complex> const & M##_x (int L) const { return M##_x_[L]; } \
            SymBandMatrix<Complex> const & M##_y (int L) const { return M##_y_[L]; } \
            Complex M##_x (int L, std::size_t i, std::size_t j) const { return M##_x_[L](i,j); } \
            Complex M##_y (int L, std::size_t i, std::size_t j) const { return M##_y_[L](i,j); } \
        
        // Access the precomputed scaled full integral moments of order L / -L-1.
        OneElectronMatrixArrayAccessors(Mtr_L)
        OneElectronMatrixArrayAccessors(Mtr_mLm1)
        
        #define OneElectronPartialMatrixAccessors(M) \
            cArrayView M##_x (int L = -1) const \
            { \
                if (L < 0) return M##_x_; \
                std::size_t mi_size = bspline_x_.Nspline() * (2 * bspline_x_.order() + 1) * (bspline_x_.order() + 1); \
                return cArrayView (M##_x_, L * mi_size, mi_size); \
            } \
            cArrayView M##_y (int L = -1) const \
            { \
                if (L < 0) return M##_y_; \
                std::size_t mi_size = bspline_y_.Nspline() * (2 * bspline_y_.order() + 1) * (bspline_y_.order() + 1); \
                return cArrayView (M##_y_, L * mi_size, mi_size); \
            }
        
        // Access the precomputed scaled partial overlap matrices of order L / -L-1.
        OneElectronPartialMatrixAccessors(Mitr_L)
        OneElectronPartialMatrixAccessors(Mitr_mLm1)
        
        // Return reference to the precomputed matrix of two-electron integrals for given multipole.
        BlockSymBandMatrix<Complex> const & R_tr_dia (unsigned i) const { return R_tr_dia_[i]; }
        CsrMatrix<LU_int_t,Complex> const & R_coupled (unsigned i) const { return R_coupled_[i]; }
        
    private:
        
        // B-spline environment
        Bspline bspline_x_;
        Bspline bspline_y_;
        
        // Asymptotic radiuses
        // - These are somewhat arbitrary radiuses that restrict position dependence of the Coulombic
        //   interaction. When either of the coordinates exceeds the limit, it is assumed to be exactly equal
        //   to the bound instead, making the potential constant with respect to further change of that
        //   coordinate. These bounds are sufficiently far from the real grid to avoid corruption of the
        //   physical solution. In single-domain problem it just stops position dependence of the potential
        //   for large distances, which is irrelevant. The importance lies in the multi-domain solution,
        //   where the potential may diverge in prepended complex grid.
        Real rxmin_, rymin_, rxmax_, rymax_;
        
        // Gauss-Legendre integrator
        GaussLegendre g_x_;
        GaussLegendre g_y_;
        
        // one-electron moment and overlap matrices
        SymBandMatrix<Complex> D_x_, S_x_, Mm1_x_, Mm2_x_;
        SymBandMatrix<Complex> D_y_, S_y_, Mm1_y_, Mm2_y_;
        
        // one-electron full integral moments for various orders (used to calculate R-integrals)
        std::vector<SymBandMatrix<Complex>> Mtr_L_x_, Mtr_mLm1_x_;
        std::vector<SymBandMatrix<Complex>> Mtr_L_y_, Mtr_mLm1_y_;
        
        // partial one-electron integral moments for various orders (used to calculate R-integrals)
        cArray Mitr_L_x_, Mitr_mLm1_x_;
        cArray Mitr_L_y_, Mitr_mLm1_y_;
        
        // two-electron integral matrices
        std::vector<BlockSymBandMatrix<Complex>> R_tr_dia_;
        std::vector<CsrMatrix<LU_int_t,Complex>> R_coupled_;
        
        // verbose output
        bool verbose_;
        
        // number of multipole matrices
        int Nlambdas_;
};

#endif // HEX_ECS_RADIAL_H