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              :  * Differentiation of contracted Gaussian functions.
      17              :  * Derivatives are then used to calculate analytic gradients of 1-el integrals.
      18              :  */
      19              : #include <math.h>
      20              : #include <stdlib.h>
      21              : 
      22              : #ifndef M_PI
      23              : #define M_PI 3.14159265358979323846
      24              : #endif
      25              : 
      26              : #include "grpp_diff_gaussian.h"
      27              : 
      28              : #include "libgrpp.h"
      29              : 
      30              : static double norm_factor(double alpha, int L);
      31              : 
      32              : /**
      33              :  * Performs differentiation of a contracted Gaussian.
      34              :  *
      35              :  * Note that the "2 alpha" factors are absorbed into coefficients, while the 'n'
      36              :  * factor is not. The latter must be accounted for explicitly at the stage of
      37              :  * gradient construction. For more details, see: T. Helgaker, P. Jorgensen, J.
      38              :  * Olsen, Molecular Electronic-Structure Theory, John Wiley & Sons Ltd, 2000.
      39              :  * Chapter 9.2.2, "Recurrence relations for Cartesian Gaussians"
      40              :  *
      41              :  */
      42            0 : void libgrpp_differentiate_shell(libgrpp_shell_t *shell,
      43              :                                  libgrpp_shell_t **shell_minus,
      44              :                                  libgrpp_shell_t **shell_plus) {
      45              :   // downwards
      46            0 :   if (shell->L > 0) {
      47            0 :     *shell_minus =
      48            0 :         libgrpp_new_shell(shell->origin, shell->L - 1, shell->num_primitives,
      49              :                           shell->coeffs, shell->alpha);
      50              : 
      51            0 :     for (int i = 0; i < shell->num_primitives; i++) {
      52              : 
      53            0 :       double alpha = shell->alpha[i];
      54            0 :       double L = shell->L;
      55              : 
      56            0 :       (*shell_minus)->coeffs[i] *=
      57            0 :           norm_factor(alpha, L) / norm_factor(alpha, L - 1);
      58              :     }
      59              :   } else {
      60            0 :     *shell_minus = NULL;
      61              :   }
      62              : 
      63              :   // upwards
      64            0 :   *shell_plus =
      65            0 :       libgrpp_new_shell(shell->origin, shell->L + 1, shell->num_primitives,
      66              :                         shell->coeffs, shell->alpha);
      67            0 :   for (int i = 0; i < shell->num_primitives; i++) {
      68              : 
      69            0 :     double alpha = shell->alpha[i];
      70            0 :     double L = shell->L;
      71              : 
      72            0 :     (*shell_plus)->coeffs[i] *=
      73            0 :         2.0 * alpha * norm_factor(alpha, L) / norm_factor(alpha, L + 1);
      74              :   }
      75            0 : }
      76              : 
      77              : /**
      78              :  * Calculates normalization factor for the primitive Gaussian
      79              :  * with the exponential parameter 'alpha' and angular momentum L.
      80              :  */
      81            0 : static double norm_factor(double alpha, int L) {
      82            0 :   return pow(2 * alpha / M_PI, 0.75) * pow(4 * alpha, 0.5 * L);
      83              : }