Line data    Source code

       1             : /*----------------------------------------------------------------------------*/
       2             : /*  CP2K: A general program to perform molecular dynamics simulations         */
       3             : /*  Copyright 2000-2025 CP2K developers group <https://cp2k.org>              */
       4             : /*                                                                            */
       5             : /*  SPDX-License-Identifier: MIT                                              */
       6             : /*----------------------------------------------------------------------------*/
       7             : 
       8             : /*
       9             :  *  libgrpp - a library for the evaluation of integrals over
      10             :  *            generalized relativistic pseudopotentials.
      11             :  *
      12             :  *  Copyright (C) 2021-2023 Alexander Oleynichenko
      13             :  */
      14             : 
      15             : /*
      16             :  * Implementation of the Gn(x) auxiliary function.
      17             :  * This function is required to calculate integrals over the 1/r^2 operator.
      18             :  *
      19             :  * More on the Gn(x) function evaluation:
      20             :  * (1) J. 0. Jensen, A. H. Cameri, C. P. Vlahacos, D. Zeroka, H. F. Hameka, C.
      21             :  * N. Merrow, Evaluation of one-electron integrals for arbitrary operators V(r)
      22             :  * over cartesian Gaussians: Application to inverse-square distance and Yukawa
      23             :  * operators. J. Comput. Chem. 14(8), 986 (1993). doi: 10.1002/jcc.540140814 (2)
      24             :  * B. Gao, A. J. Thorvaldsen, K. Ruud, GEN1INT: A unified procedure for the
      25             :  * evaluation of one-electron integrals over Gaussian basis functions and their
      26             :  * geometric derivatives. Int. J. Quantum Chem. 111(4), 858 (2011).
      27             :  *     doi: 10.1002/qua.22886
      28             :  */
      29             : #include <math.h>
      30             : #include <string.h>
      31             : 
      32             : #ifndef M_PI
      33             : #define M_PI 3.1415926535897932384626433
      34             : #endif
      35             : 
      36             : #include "grpp_factorial.h"
      37             : #include "grpp_specfunc.h"
      38             : 
      39             : static double gfun_taylor(int n, double x);
      40             : 
      41             : /**
      42             :  * Calculates values of the Gn(x) auxiliary function for n = 0, ..., nmax
      43             :  * and stores them into the g[] array.
      44             :  */
      45           0 : void libgrpp_gfun_values(double x, int nmax, double *g) {
      46           0 :   memset(g, 0, (nmax + 1) * sizeof(double));
      47             : 
      48           0 :   if (x <= 12.0) {
      49             :     /*
      50             :      * downward recursion
      51             :      */
      52           0 :     g[nmax] = gfun_taylor(nmax, x);
      53             : 
      54           0 :     for (int n = nmax; n > 0; n--) {
      55           0 :       g[n - 1] = (1.0 - 2.0 * x * g[n]) / (2.0 * n - 1.0);
      56             :     }
      57             :   } else {
      58             :     /*
      59             :      * upward recursion
      60             :      */
      61           0 :     double sqrt_x = sqrt(x);
      62           0 :     g[0] = libgrpp_Dawsons_Integral(sqrt_x) / sqrt_x;
      63             : 
      64           0 :     for (int n = 0; n < nmax; n++) {
      65           0 :       g[n + 1] = (1.0 - (2 * n + 1) * g[n]) / (2.0 * x);
      66             :     }
      67             :   }
      68           0 : }
      69             : 
      70             : /**
      71             :  * Calculates value of the Gn(x) auxiliary function using the Taylor expansion.
      72             :  * The Taylor series converges for x <= 30.
      73             :  */
      74           0 : static double gfun_taylor(int n, double x) {
      75           0 :   const double thresh = 1e-15;
      76           0 :   double sum = 0.0;
      77             : 
      78           0 :   for (int k = 0; k < 100; k++) {
      79             : 
      80           0 :     double y_exp = exp(-x);
      81           0 :     double y_pow = pow(x, k);
      82           0 :     double y_fac = libgrpp_factorial(k);
      83           0 :     double y_nk1 = 2.0 * n + 2.0 * k + 1.0;
      84             : 
      85           0 :     double contrib = y_exp * y_pow / y_fac / y_nk1;
      86           0 :     sum += contrib;
      87             : 
      88           0 :     if (fabs(contrib) < thresh) {
      89             :       break;
      90             :     }
      91             :   }
      92             : 
      93           0 :   return sum;
      94             : }