1 /* 2 * Copyright (c) 2020 Robert Shaw 3 * This file is a part of Libecpint. 4 * 5 * Permission is hereby granted, free of charge, to any person obtaining 6 * a copy of this software and associated documentation files (the 7 * "Software"), to deal in the Software without restriction, including 8 * without limitation the rights to use, copy, modify, merge, publish, 9 * distribute, sublicense, and/or sell copies of the Software, and to 10 * permit persons to whom the Software is furnished to do so, subject to 11 * the following conditions: 12 * 13 * The above copyright notice and this permission notice shall be 14 * included in all copies or substantial portions of the Software. 15 * 16 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 17 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF 18 * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND 19 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 20 * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION 21 * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION 22 * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 23 */ 24 25 #ifndef GC_QUAD_HEAD 26 #define GC_QUAD_HEAD 27 28 #include <functional> 29 #include <vector> 30 31 namespace libecpint { 32 33 /// Different choices of integration algorithm, see references 34 enum GCTYPE { 35 ONEPOINT, ///< Described in Perez92 36 TWOPOINT ///< Described in Perez93 37 }; 38 39 /** 40 * \class GCQuadrature 41 * \brief Performs adaptive Gauss-Chebyshev quadrature of the second kind for any given function. 42 * 43 * Stores the weights and abscissae for the quadrature, and provides two different methods to integrate on [-1, 1] 44 * Also contains means to transform the region of integration to [0, infinity) and [rmin, rmax] 45 * 46 * REFERENCES: 47 * (Perez92) J.M. Perez-Jorda et al., Comput. Phys. Comm. 70 (1992), 271-284 48 * (Perez93) J.M. Perez-Jorda et al., Comput. Phys. Comm. 77 (1993), 46-56 49 * (Krack98) M. Krack, A.M. Koster, J. Chem. Phys. 108 (1998), 3226 - 3234 50 * (Flores06) R. Flores-Moreno et al., J. Comput. Chem. 27 (2006), 1009-1019 51 */ 52 class GCQuadrature { 53 private: 54 int maxN; ///< Maximum number of points to use in quadrature 55 int M; ///< Index of midpoint 56 57 std::vector<double> x; ///< Weights 58 std::vector<double> w; ///< Abscissae 59 60 double I; ///< Integral value 61 62 GCTYPE t; ///< Algorithm type to be used 63 64 /// Worker function for integration routines, should not be called directly. 65 double sumTerms(const std::function<double(double, const double*, int)> &f, 66 const double *p, int limit, int start, int end, int shift, int skip) const; 67 68 public: 69 70 /// Default constructor, creates empty object 71 GCQuadrature(); 72 73 /// Copy constructor, carbon copies all members 74 GCQuadrature(const GCQuadrature &other); 75 76 /** 77 * Intialises the integration grid to the given number of points, and integration type. 78 * ONEPOINT will choose N = 2^n - 1 closest to the given number of points, whilst 79 * TWOPOINT will choose N= 3*2^n - 1 in the same way. 80 * 81 * @param points - maximum number of quadrature points to be used 82 * @param t - the algorithm to be used (ONEPOINT / TWOPOINT) 83 */ 84 void initGrid(int points, GCTYPE t); 85 86 /** 87 * Integrates the given function (over [-1, 1] by default) to within the given tolerance. 88 * @param f - the function to be integrated 89 * @param params - array of parameters for the function to be integrated 90 * @param tolerance - change below which convergenced is considered to be achieved 91 * @param start - the index of the first point used in the integration 92 * @param end - the index of the last point used in the integration 93 * @returns the integral (first) and true if integration converged, false otherwise (second) 94 */ 95 std::pair<double, bool> integrate( 96 std::function<double(double, const double*, int)> &f, 97 const double *params, double tolerance, int start, int end) const; 98 99 /** 100 * Transforms the region of integration to [0, inf) using the logarithmic transformation of Krack98 101 */ 102 void transformZeroInf(); 103 104 /** 105 * Transforms region of integration to [rmin, rmax] using the linear transformation from Flores06, assuming 106 * a Gaussian envelope. rmin/rmax are the distances from the centre of the envelope such that the integrand is effectively zero. 107 * @param z - the exponent of the Gaussian envelope 108 * @param p - the centre of the Gaussian envelope 109 */ 110 void transformRMinMax(double z, double p); // Transfromation from [-1, 1] to [rmin, rmax] from Flores06 111 void untransformRMinMax(double z, double p); 112 113 /// @return the maximum number of quadrature points getN() const114 int getN() const { return maxN; } 115 116 /// @return a reference to the abscissae getX()117 std::vector<double>& getX() { return x; } getX() const118 const std::vector<double>& getX() const { return x; } 119 }; 120 } 121 122 #endif 123