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