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              :  * 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       350780 : libgrpp_potential_t *libgrpp_new_potential(int L, int J, int num_primitives,
      33              :                                            int *powers, double *coeffs,
      34              :                                            double *alpha) {
      35       350780 :   libgrpp_potential_t *pot =
      36       350780 :       (libgrpp_potential_t *)calloc(1, sizeof(libgrpp_potential_t));
      37              : 
      38       350780 :   pot->L = L;
      39       350780 :   pot->J = J;
      40       350780 :   pot->powers = (int *)calloc(num_primitives, sizeof(int));
      41       350780 :   pot->coeffs = (double *)calloc(num_primitives, sizeof(double));
      42       350780 :   pot->alpha = (double *)calloc(num_primitives, sizeof(double));
      43              : 
      44       350780 :   pot->num_primitives = 0;
      45      1235765 :   for (int i = 0; i < num_primitives; i++) {
      46       884985 :     if (fabs(coeffs[i]) < LIBGRPP_ZERO_THRESH) {
      47          941 :       continue;
      48              :     }
      49              : 
      50       884044 :     pot->coeffs[pot->num_primitives] = coeffs[i];
      51       884044 :     pot->powers[pot->num_primitives] = powers[i];
      52       884044 :     pot->alpha[pot->num_primitives] = alpha[i];
      53              : 
      54       884044 :     pot->num_primitives++;
      55              :   }
      56              : 
      57       350780 :   return pot;
      58              : }
      59              : 
      60              : /*
      61              :  * destructor for the pseudopotential
      62              :  */
      63       350780 : void libgrpp_delete_potential(libgrpp_potential_t *potential) {
      64       350780 :   if (potential == NULL) {
      65              :     return;
      66              :   }
      67              : 
      68       350780 :   free(potential->powers);
      69       350780 :   free(potential->coeffs);
      70       350780 :   free(potential->alpha);
      71       350780 :   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     38916548 : double libgrpp_potential_value(libgrpp_potential_t *potential, double r) {
      80     38916548 :   double val = 0.0;
      81     38916548 :   double r_2 = r * r;
      82              : 
      83    206574668 :   for (int i = 0; i < potential->num_primitives; i++) {
      84    167658120 :     int n = potential->powers[i];
      85    167658120 :     val +=
      86    167658120 :         potential->coeffs[i] * pow(r, n - 2) * exp(-potential->alpha[i] * r_2);
      87              :   }
      88              : 
      89     38916548 :   return val;
      90              : }
      91              : 
      92              : /*
      93              :  * removes redundant (zero) primitives from the RPP.
      94              :  * argument remains constant.
      95              :  */
      96              : libgrpp_potential_t *
      97       156604 : libgrpp_shrink_potential(libgrpp_potential_t *src_potential) {
      98       156604 :   int n = src_potential->num_primitives;
      99       156604 :   int *new_powers = calloc(n, sizeof(int));
     100       156604 :   double *new_coeffs = calloc(n, sizeof(double));
     101       156604 :   double *new_alpha = calloc(n, sizeof(double));
     102              : 
     103       156604 :   int n_nonzero_primitives = 0;
     104       516500 :   for (int i = 0; i < n; i++) {
     105       359896 :     if (fabs(src_potential->coeffs[i]) > LIBGRPP_ZERO_THRESH) {
     106       359896 :       new_powers[n_nonzero_primitives] = src_potential->powers[i];
     107       359896 :       new_coeffs[n_nonzero_primitives] = src_potential->coeffs[i];
     108       359896 :       new_alpha[n_nonzero_primitives] = src_potential->alpha[i];
     109       359896 :       n_nonzero_primitives++;
     110              :     }
     111              :   }
     112              : 
     113       156604 :   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       156604 :   free(new_powers);
     118       156604 :   free(new_coeffs);
     119       156604 :   free(new_alpha);
     120              : 
     121       156604 :   return new_pot;
     122              : }