Line data    Source code

       1              : /*----------------------------------------------------------------------------*/
       2              : /*  CP2K: A general program to perform molecular dynamics simulations         */
       3              : /*  Copyright 2000-2026 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              :  * representation of (generalized) effective core potentials
      17              :  */
      18              : 
      19              : #include <math.h>
      20              : #include <stdlib.h>
      21              : #include <string.h>
      22              : 
      23              : #ifndef M_PI
      24              : #define M_PI 3.1415926535897932384626433
      25              : #endif
      26              : 
      27              : #include "libgrpp.h"
      28              : 
      29              : /**
      30              :  * constructor for the pseudopotential
      31              :  */
      32       661520 : libgrpp_potential_t *libgrpp_new_potential(int L, int J, int num_primitives,
      33              :                                            int *powers, double *coeffs,
      34              :                                            double *alpha) {
      35       661520 :   libgrpp_potential_t *pot =
      36       661520 :       (libgrpp_potential_t *)calloc(1, sizeof(libgrpp_potential_t));
      37              : 
      38       661520 :   pot->L = L;
      39       661520 :   pot->J = J;
      40       661520 :   pot->powers = (int *)calloc(num_primitives, sizeof(int));
      41       661520 :   pot->coeffs = (double *)calloc(num_primitives, sizeof(double));
      42       661520 :   pot->alpha = (double *)calloc(num_primitives, sizeof(double));
      43              : 
      44       661520 :   pot->num_primitives = 0;
      45      2070482 :   for (int i = 0; i < num_primitives; i++) {
      46      1408962 :     if (fabs(coeffs[i]) < LIBGRPP_ZERO_THRESH) {
      47        34050 :       continue;
      48              :     }
      49              : 
      50      1374912 :     pot->coeffs[pot->num_primitives] = coeffs[i];
      51      1374912 :     pot->powers[pot->num_primitives] = powers[i];
      52      1374912 :     pot->alpha[pot->num_primitives] = alpha[i];
      53              : 
      54      1374912 :     pot->num_primitives++;
      55              :   }
      56              : 
      57       661520 :   return pot;
      58              : }
      59              : 
      60              : /*
      61              :  * destructor for the pseudopotential
      62              :  */
      63       661520 : void libgrpp_delete_potential(libgrpp_potential_t *potential) {
      64       661520 :   if (potential == NULL) {
      65              :     return;
      66              :   }
      67              : 
      68       661520 :   free(potential->powers);
      69       661520 :   free(potential->coeffs);
      70       661520 :   free(potential->alpha);
      71       661520 :   free(potential);
      72              : }
      73              : 
      74              : /*
      75              :  * calculates value of the pseudopotential at the point 'r'
      76              :  *
      77              :  * TODO: remove the invocation of 'pow()'
      78              :  */
      79     54819711 : double libgrpp_potential_value(libgrpp_potential_t *potential, double r) {
      80     54819711 :   double val = 0.0;
      81     54819711 :   double r_2 = r * r;
      82              : 
      83    253365703 :   for (int i = 0; i < potential->num_primitives; i++) {
      84    198545992 :     int n = potential->powers[i];
      85    198545992 :     val +=
      86    198545992 :         potential->coeffs[i] * pow(r, n - 2) * exp(-potential->alpha[i] * r_2);
      87              :   }
      88              : 
      89     54819711 :   return val;
      90              : }
      91              : 
      92              : /*
      93              :  * removes redundant (zero) primitives from the RPP.
      94              :  * argument remains constant.
      95              :  */
      96              : libgrpp_potential_t *
      97       309601 : libgrpp_shrink_potential(libgrpp_potential_t *src_potential) {
      98       309601 :   int n = src_potential->num_primitives;
      99       309601 :   int *new_powers = calloc(n, sizeof(int));
     100       309601 :   double *new_coeffs = calloc(n, sizeof(double));
     101       309601 :   double *new_alpha = calloc(n, sizeof(double));
     102              : 
     103       309601 :   int n_nonzero_primitives = 0;
     104       911097 :   for (int i = 0; i < n; i++) {
     105       601496 :     if (fabs(src_potential->coeffs[i]) > LIBGRPP_ZERO_THRESH) {
     106       601496 :       new_powers[n_nonzero_primitives] = src_potential->powers[i];
     107       601496 :       new_coeffs[n_nonzero_primitives] = src_potential->coeffs[i];
     108       601496 :       new_alpha[n_nonzero_primitives] = src_potential->alpha[i];
     109       601496 :       n_nonzero_primitives++;
     110              :     }
     111              :   }
     112              : 
     113       309601 :   libgrpp_potential_t *new_pot = libgrpp_new_potential(
     114              :       src_potential->L, src_potential->J, n_nonzero_primitives, new_powers,
     115              :       new_coeffs, new_alpha);
     116              : 
     117       309601 :   free(new_powers);
     118       309601 :   free(new_coeffs);
     119       309601 :   free(new_alpha);
     120              : 
     121       309601 :   return new_pot;
     122              : }