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             :  * Calculation of momentum integrals.
      17             :  *
      18             :  * For details, see:
      19             :  * T. Helgaker, P. Jorgensen, J. Olsen, Molecular Electronic-Structure Theory,
      20             :  * John Wiley & Sons Ltd, 2000.
      21             :  * Chapter 9.3.4, "Momentum and kinetic-energy integrals"
      22             :  */
      23             : #include <math.h>
      24             : #include <stdlib.h>
      25             : #include <string.h>
      26             : #ifndef M_PI
      27             : #define M_PI 3.14159265358979323846
      28             : #endif
      29             : #include "grpp_momentum.h"
      30             : 
      31             : #include "grpp_norm_gaussian.h"
      32             : #include "grpp_utils.h"
      33             : #include "libgrpp.h"
      34             : 
      35             : static void momentum_integrals_shell_pair_obara_saika(
      36             :     libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double alpha_A,
      37             :     double alpha_B, double *momentum_x_matrix, double *momentum_y_matrix,
      38             :     double *momentum_z_matrix);
      39             : 
      40             : /**
      41             :  * returns imaginary(!) part of integrals over the momentum operator p = -i
      42             :  * \hbar \nabla. The "minus" sign is included.
      43             :  */
      44           0 : void libgrpp_momentum_integrals(libgrpp_shell_t *shell_A,
      45             :                                 libgrpp_shell_t *shell_B,
      46             :                                 double *momentum_x_matrix,
      47             :                                 double *momentum_y_matrix,
      48             :                                 double *momentum_z_matrix) {
      49           0 :   int size_A = libgrpp_get_shell_size(shell_A);
      50           0 :   int size_B = libgrpp_get_shell_size(shell_B);
      51             : 
      52           0 :   double *buf_x = calloc(size_A * size_B, sizeof(double));
      53           0 :   double *buf_y = calloc(size_A * size_B, sizeof(double));
      54           0 :   double *buf_z = calloc(size_A * size_B, sizeof(double));
      55             : 
      56           0 :   memset(momentum_x_matrix, 0, size_A * size_B * sizeof(double));
      57           0 :   memset(momentum_y_matrix, 0, size_A * size_B * sizeof(double));
      58           0 :   memset(momentum_z_matrix, 0, size_A * size_B * sizeof(double));
      59             : 
      60             :   // loop over primitives in contractions
      61           0 :   for (int i = 0; i < shell_A->num_primitives; i++) {
      62           0 :     for (int j = 0; j < shell_B->num_primitives; j++) {
      63           0 :       double alpha_i = shell_A->alpha[i];
      64           0 :       double alpha_j = shell_B->alpha[j];
      65           0 :       double coef_A_i = shell_A->coeffs[i];
      66           0 :       double coef_B_j = shell_B->coeffs[j];
      67             : 
      68           0 :       momentum_integrals_shell_pair_obara_saika(shell_A, shell_B, alpha_i,
      69             :                                                 alpha_j, buf_x, buf_y, buf_z);
      70             : 
      71           0 :       libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf_x,
      72             :                     momentum_x_matrix);
      73           0 :       libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf_y,
      74             :                     momentum_y_matrix);
      75           0 :       libgrpp_daxpy(size_A * size_B, coef_A_i * coef_B_j, buf_z,
      76             :                     momentum_z_matrix);
      77             :     }
      78             :   }
      79             : 
      80           0 :   free(buf_x);
      81           0 :   free(buf_y);
      82           0 :   free(buf_z);
      83           0 : }
      84             : 
      85           0 : static void momentum_integrals_shell_pair_obara_saika(
      86             :     libgrpp_shell_t *shell_A, libgrpp_shell_t *shell_B, double alpha_A,
      87             :     double alpha_B, double *momentum_x_matrix, double *momentum_y_matrix,
      88             :     double *momentum_z_matrix) {
      89           0 :   int size_A = libgrpp_get_shell_size(shell_A);
      90           0 :   int size_B = libgrpp_get_shell_size(shell_B);
      91           0 :   int L_A = shell_A->L;
      92           0 :   int L_B = shell_B->L;
      93           0 :   double N_A = libgrpp_gaussian_norm_factor(L_A, 0, 0, alpha_A);
      94           0 :   double N_B = libgrpp_gaussian_norm_factor(L_B, 0, 0, alpha_B);
      95             : 
      96           0 :   double p = alpha_A + alpha_B;
      97           0 :   double mu = alpha_A * alpha_B / (alpha_A + alpha_B);
      98           0 :   double *A = shell_A->origin;
      99           0 :   double *B = shell_B->origin;
     100             : 
     101             :   // calculate S_ij
     102           0 :   double S[3][LIBGRPP_MAX_BASIS_L + 1][LIBGRPP_MAX_BASIS_L + 1];
     103             : 
     104           0 :   for (int coord = 0; coord < 3; coord++) {
     105           0 :     double P = (alpha_A * A[coord] + alpha_B * B[coord]) / p;
     106             : 
     107           0 :     double X_AB = A[coord] - B[coord];
     108           0 :     double X_PA = P - A[coord];
     109           0 :     double X_PB = P - B[coord];
     110           0 :     double pfac = 1.0 / (2.0 * p);
     111             : 
     112           0 :     for (int i = 0; i <= L_A + 1; i++) {
     113           0 :       for (int j = 0; j <= L_B + 1; j++) {
     114           0 :         double S_ij = 0.0;
     115             : 
     116           0 :         if (i + j == 0) {
     117           0 :           S[coord][0][0] = sqrt(M_PI / p) * exp(-mu * X_AB * X_AB);
     118           0 :           continue;
     119             :         }
     120             : 
     121           0 :         if (i == 0) { // upward by j
     122           0 :           S_ij += X_PB * S[coord][i][j - 1];
     123           0 :           if (j - 1 > 0) {
     124           0 :             S_ij += (j - 1) * pfac * S[coord][i][j - 2];
     125             :           }
     126             :         } else { // upward by i
     127           0 :           S_ij += X_PA * S[coord][i - 1][j];
     128           0 :           if (i - 1 > 0) {
     129           0 :             S_ij += (i - 1) * pfac * S[coord][i - 2][j];
     130             :           }
     131           0 :           if (j > 0) {
     132           0 :             S_ij += j * pfac * S[coord][i - 1][j - 1];
     133             :           }
     134             :         }
     135             : 
     136           0 :         S[coord][i][j] = S_ij;
     137             :       }
     138             :     }
     139             :   }
     140             : 
     141             :   // calculate D^1_ij
     142             : 
     143             :   double D1[3][LIBGRPP_MAX_BASIS_L][LIBGRPP_MAX_BASIS_L];
     144             : 
     145           0 :   for (int coord = 0; coord < 3; coord++) {
     146           0 :     for (int i = 0; i <= L_A; i++) {
     147           0 :       for (int j = 0; j <= L_B; j++) {
     148             : 
     149           0 :         double D1_ij = 0.0;
     150           0 :         D1_ij += 2.0 * alpha_A * S[coord][i + 1][j];
     151           0 :         if (i >= 1) {
     152           0 :           D1_ij -= i * S[coord][i - 1][j];
     153             :         }
     154             : 
     155           0 :         D1[coord][i][j] = D1_ij;
     156             :       }
     157             :     }
     158             :   }
     159             : 
     160             :   // loop over cartesian functions inside the shells
     161           0 :   for (int m = 0; m < size_A; m++) {
     162           0 :     for (int n = 0; n < size_B; n++) {
     163           0 :       int n_A = shell_A->cart_list[3 * m + 0];
     164           0 :       int l_A = shell_A->cart_list[3 * m + 1];
     165           0 :       int m_A = shell_A->cart_list[3 * m + 2];
     166           0 :       int n_B = shell_B->cart_list[3 * n + 0];
     167           0 :       int l_B = shell_B->cart_list[3 * n + 1];
     168           0 :       int m_B = shell_B->cart_list[3 * n + 2];
     169             : 
     170           0 :       momentum_x_matrix[m * size_B + n] =
     171           0 :           -N_A * N_B * D1[0][n_A][n_B] * S[1][l_A][l_B] * S[2][m_A][m_B];
     172           0 :       momentum_y_matrix[m * size_B + n] =
     173           0 :           -N_A * N_B * S[0][n_A][n_B] * D1[1][l_A][l_B] * S[2][m_A][m_B];
     174           0 :       momentum_z_matrix[m * size_B + n] =
     175           0 :           -N_A * N_B * S[0][n_A][n_B] * S[1][l_A][l_B] * D1[2][m_A][m_B];
     176             :     }
     177             :   }
     178           0 : }