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