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