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             : #include "libgrpp.h"
      16             : 
      17             : #include <assert.h>
      18             : #include <stdio.h>
      19             : #include <stdlib.h>
      20             : #include <string.h>
      21             : 
      22             : #include "grpp_utils.h"
      23             : 
      24             : /**
      25             :  * Evaluates integrals over the full GRPP operator consisting of three parts:
      26             :  * - scalar relativistic (local)
      27             :  * - scalar relativistic (semi-local)
      28             :  * - effective spin-orbit (semi-local)
      29             :  * - outercore potentials (non-local)
      30             :  *
      31             :  * See libgrpp.h for the definition of the libgrpp_grpp_t structure.
      32             :  */
      33           0 : void libgrpp_full_grpp_integrals(libgrpp_shell_t *shell_A,
      34             :                                  libgrpp_shell_t *shell_B,
      35             :                                  libgrpp_grpp_t *grpp_operator,
      36             :                                  double *grpp_origin, double *arep_matrix,
      37             :                                  double *so_x_matrix, double *so_y_matrix,
      38             :                                  double *so_z_matrix) {
      39           0 :   assert(libgrpp_is_initialized());
      40             : 
      41           0 :   size_t size = shell_A->cart_size * shell_B->cart_size;
      42           0 :   double *buf_arep = (double *)calloc(size, sizeof(double));
      43           0 :   double *buf_so_x = (double *)calloc(size, sizeof(double));
      44           0 :   double *buf_so_y = (double *)calloc(size, sizeof(double));
      45           0 :   double *buf_so_z = (double *)calloc(size, sizeof(double));
      46             : 
      47           0 :   memset(arep_matrix, 0, sizeof(double) * size);
      48           0 :   memset(so_x_matrix, 0, sizeof(double) * size);
      49           0 :   memset(so_y_matrix, 0, sizeof(double) * size);
      50           0 :   memset(so_z_matrix, 0, sizeof(double) * size);
      51             : 
      52             :   /*
      53             :    * radially-local ("type-1") integrals
      54             :    */
      55           0 :   libgrpp_type1_integrals(shell_A, shell_B, grpp_origin, grpp_operator->U_L,
      56             :                           buf_arep);
      57           0 :   libgrpp_daxpy(size, 1.0, buf_arep, arep_matrix);
      58             : 
      59             :   /*
      60             :    * semilocal AREP ("type-2") integrals
      61             :    */
      62           0 :   for (int L = 0; L < grpp_operator->n_arep; L++) {
      63           0 :     libgrpp_type2_integrals(shell_A, shell_B, grpp_origin,
      64           0 :                             grpp_operator->U_arep[L], buf_arep);
      65           0 :     libgrpp_daxpy(size, 1.0, buf_arep, arep_matrix);
      66             :   }
      67             : 
      68             :   /*
      69             :    * semilocal SO ("type-3") integrals
      70             :    */
      71           0 :   for (int i_so = 0; i_so < grpp_operator->n_esop; i_so++) {
      72           0 :     libgrpp_potential_t *so_potential = grpp_operator->U_esop[i_so];
      73           0 :     libgrpp_spin_orbit_integrals(shell_A, shell_B, grpp_origin, so_potential,
      74             :                                  buf_so_x, buf_so_y, buf_so_z);
      75             : 
      76           0 :     int L = so_potential->L;
      77           0 :     libgrpp_daxpy(size, 2.0 / (2 * L + 1), buf_so_x, so_x_matrix);
      78           0 :     libgrpp_daxpy(size, 2.0 / (2 * L + 1), buf_so_y, so_y_matrix);
      79           0 :     libgrpp_daxpy(size, 2.0 / (2 * L + 1), buf_so_z, so_z_matrix);
      80             :   }
      81             : 
      82             :   /*
      83             :    * integrals over outercore non-local potentials,
      84             :    * the part specific for GRPP.
      85             :    *
      86             :    * note that proper pre-factors for the SO part are calculated inside
      87             :    * the libgrpp_outercore_potential_integrals() procedure.
      88             :    */
      89           0 :   libgrpp_outercore_potential_integrals(
      90             :       shell_A, shell_B, grpp_origin, grpp_operator->n_oc_shells,
      91             :       grpp_operator->U_oc, grpp_operator->oc_shells, buf_arep, buf_so_x,
      92             :       buf_so_y, buf_so_z);
      93             : 
      94           0 :   libgrpp_daxpy(size, 1.0, buf_arep, arep_matrix);
      95           0 :   libgrpp_daxpy(size, 1.0, buf_so_x, so_x_matrix);
      96           0 :   libgrpp_daxpy(size, 1.0, buf_so_y, so_y_matrix);
      97           0 :   libgrpp_daxpy(size, 1.0, buf_so_z, so_z_matrix);
      98             : 
      99             :   /*
     100             :    * cleanup
     101             :    */
     102           0 :   free(buf_arep);
     103           0 :   free(buf_so_x);
     104           0 :   free(buf_so_y);
     105           0 :   free(buf_so_z);
     106           0 : }