From fb97b44a5ff63bfc2ae9fb7465ce6ac833b49d08 Mon Sep 17 00:00:00 2001 From: joguenzl Date: Fri, 31 Mar 2023 13:01:51 +0200 Subject: [PATCH 1/6] Initial support for numerical atom-centered orbitals --- src/pytrexio.i | 8 + src/templates_front/templator_front.org | 369 +++++++++++++++++++++++- trex.org | 113 ++++++-- 3 files changed, 461 insertions(+), 29 deletions(-) diff --git a/src/pytrexio.i b/src/pytrexio.i index b6e5485..13c9e10 100644 --- a/src/pytrexio.i +++ b/src/pytrexio.i @@ -108,6 +108,14 @@ import_array(); /* For some reasons SWIG does not apply the proper bitfield_t typemap, so one has to manually specify int64_t* ARGOUT_ARRAY1 below */ %apply (int64_t* ARGOUT_ARRAY1, int32_t DIM1) {(bitfield_t* const bit_list, const int32_t N_int)}; +/* NAO functions */ +%apply double *OUTPUT { double* const log_r_out, double* const amplitude}; +%apply (double *ARGOUT_ARRAY1, int DIM1) {(double* const amplitudes, int amplitude_cnt)}; +%apply (double* IN_ARRAY1, int DIM1) {(double* grid_r, int n_grid_r), + (double* interpolator, int n_interp), (double* nucleus_coords, int n_nuc_co), (double* normalization, int n_norm)}; +%apply (int64_t* IN_ARRAY1, int DIM1) {(int64_t* grid_start, int n_grid_st), + (int64_t* grid_size, int n_grid_si), (int64_t* nucleus_index, int n_nuc_id)}; + /* This tells SWIG to treat char ** dset_in pattern as a special case Enables access to trexio_[...]_write_dset_str set of functions directly, i.e. by converting input list of strings from Python into char ** of C diff --git a/src/templates_front/templator_front.org b/src/templates_front/templator_front.org index 91bfdf8..7a56e80 100644 --- a/src/templates_front/templator_front.org +++ b/src/templates_front/templator_front.org @@ -42,6 +42,8 @@ typedef int32_t trexio_exit_code; #include "config.h" #endif +#define _USE_MATH_DEFINES +#include #include #include #include @@ -5585,6 +5587,12 @@ def has_determinant_list(trexio_file) -> bool: bitfield representation of the determinant. If the creation of the bitfield requires a change of sign, the return code is ~TREXIO_PHASE_CHANGE~. + ~trexio_convert_nao_radius~ functions convert a radius from the center of an NAO to its logarithmic interpolation grid. + + ~trexio_evaluate_nao_radial~ function evaluates the radial function of an NAO for a given distance from the center. + + ~trexio_evaluate_nao_radial_all~ evaluates all radial functions at a given position in space. + ** C #+begin_src c :tangle prefix_front.h @@ -5603,6 +5611,35 @@ trexio_exit_code trexio_to_orbital_list_up_dn (const int32_t N_int, const bitfie trexio_exit_code trexio_safe_to_orbital_list (const int32_t N_int, const bitfield_t* dset_in, const int64_t dim_in, int32_t* const dset_out, const int64_t dim_out, int32_t* const num); trexio_exit_code trexio_safe_to_orbital_list_up_dn (const int32_t N_int, const bitfield_t* dset_in, const int64_t dim_in, int32_t* const dset_up_out, const int64_t dim_up_out, int32_t* const dset_dn_out, const int64_t dim_dn_out, int32_t* const num_up, int32_t* const num_dn); trexio_exit_code trexio_to_bitfield_list (const int32_t* orb_list, const int32_t occupied_num, bitfield_t* const bit_list, const int32_t N_int); + +trexio_exit_code trexio_convert_nao_radius_32 (const float r, + const float* const grid_r, float* const log_r_out); +trexio_exit_code trexio_convert_nao_radius_64 (const double r, + const double* const grid_r, double* const log_r_out); +trexio_exit_code trexio_convert_nao_radius_py (const double r, + double* grid_r, int32_t n_grid_r, double* const log_r_out); +trexio_exit_code trexio_evaluate_nao_radial (const int32_t shell_index, + const double r, const int32_t* const grid_start, const int32_t* const grid_size, + const double* const grid_r, const double* const interpolator, + const double* const normalization, double* const amplitude); + +trexio_exit_code trexio_evaluate_nao_radial_all (const int32_t shell_num, + const int32_t* const nucleus_index, const double* const nucleus_coords, + const int32_t* const grid_start, const int32_t* const grid_size, + const double* const grid_r, const double* const interpolator, const double* const normalization, + const double rx, const double ry, const double rz, double* const amplitude); + +trexio_exit_code trexio_evaluate_nao_radial_py (const int shell_index, + const double r, int64_t* grid_start, int n_grid_st, int64_t* grid_size, + int n_grid_si, double* grid_r, int n_grid_r, double* interpolator, + int n_interp, double* normalization, int n_norm, double* const amplitude); + +trexio_exit_code trexio_evaluate_nao_radial_all_py (const int32_t shell_num, + int64_t* nucleus_index, int n_nuc_id, double* nucleus_coords, int n_nuc_co, + int64_t* grid_start, int n_grid_st, int64_t* grid_size, int n_grid_si, + double* grid_r, int n_grid_r, double* interpolator, int n_interp, double* normalization, int n_norm, + const double rx, const double ry, const double rz, double* const amplitudes, int amplitude_cnt); + #+end_src @@ -5740,7 +5777,174 @@ trexio_safe_to_orbital_list_up_dn (const int32_t N_int, } #+end_src - #+begin_src c :tangle trexio_private.h + #+begin_src c :tangle prefix_front.c +trexio_exit_code +trexio_convert_nao_radius_32 (const float r, const float* const grid_r, float* const log_r_out) +{ + if (r < 0) return TREXIO_INVALID_ARG_1; + if (grid_r == NULL) return TREXIO_INVALID_ARG_2; + + *log_r_out = log(r / grid_r[0]) / log(grid_r[1] / grid_r[0]); + + return TREXIO_SUCCESS; +} + +trexio_exit_code +trexio_convert_nao_radius_64 (const double r, const double* const grid_r, double* const log_r_out) +{ + if (r < 0) return TREXIO_INVALID_ARG_1; + if (grid_r == NULL) return TREXIO_INVALID_ARG_2; + + *log_r_out = log(r / grid_r[0]) / log(grid_r[1] / grid_r[0]); + + return TREXIO_SUCCESS; +} + +trexio_exit_code +trexio_convert_nao_radius_py (const double r, double* grid_r, int32_t n_grid, double* const log_r_out) +{ + if (r < 0) return TREXIO_INVALID_ARG_1; + if (grid_r == NULL) return TREXIO_INVALID_ARG_2; + + *log_r_out = log(r / grid_r[0]) / log(grid_r[1] / grid_r[0]); + + return TREXIO_SUCCESS; +} + + #+end_src + + #+begin_src c :tangle prefix_front.c +trexio_exit_code +trexio_evaluate_nao_radial (const int32_t shell_index, const double r, const int32_t* const grid_start, + const int32_t* const grid_size, const double* const grid_r, const double* const interpolator, const double* const normalization, double* const amplitude) +{ + if (shell_index < 0) return TREXIO_INVALID_ARG_1; + if (r < 0) return TREXIO_INVALID_ARG_2; + if (grid_start == 0) return TREXIO_INVALID_ARG_3; + if (grid_size == 0) return TREXIO_INVALID_ARG_4; + if (grid_r == NULL) return TREXIO_INVALID_ARG_5; + if (interpolator == 0) return TREXIO_INVALID_ARG_6; + if (normalization == 0) return TREXIO_INVALID_ARG_7; + + const int32_t i0 = 4*grid_start[shell_index]; + + // Convert radius to logarithmic units + double r_log = 0; + trexio_convert_nao_radius_64 (r, grid_r + grid_start[shell_index], &r_log); + int32_t i_log = (int32_t) r_log; + if (i_log < 0) + i_log = 0; + else if (i_log >= grid_size[shell_index]) + return 0; // NAOs vanish at the boundary by definition + + double t = r_log - (double) i_log; + double val_spline = interpolator[i0 + 4*i_log + 0]; + val_spline += t * interpolator[i0 + 4*i_log + 1]; + val_spline += t * t * interpolator[i0 + 4*i_log + 2]; + val_spline += t * t * t * interpolator[i0 + 4*i_log + 3]; + + *amplitude = val_spline * normalization[shell_index] / r; + return TREXIO_SUCCESS; +} + +trexio_exit_code +trexio_evaluate_nao_radial_all (const int32_t shell_num, const int32_t* const nucleus_index, const double* const nucleus_coords, const int32_t* const grid_start, + const int32_t* const grid_size, const double* const grid_r, const double* const interpolator, + const double* const normalization, const double rx, const double ry, const double rz, double* const amplitude) +{ + if (shell_num < 0) return TREXIO_INVALID_ARG_1; + if (nucleus_index == 0) return TREXIO_INVALID_ARG_2; + if (nucleus_coords == 0) return TREXIO_INVALID_ARG_3; + if (grid_start == 0) return TREXIO_INVALID_ARG_4; + if (grid_size == 0) return TREXIO_INVALID_ARG_5; + if (grid_r == NULL) return TREXIO_INVALID_ARG_6; + if (interpolator == 0) return TREXIO_INVALID_ARG_7; + if (normalization == 0) return TREXIO_INVALID_ARG_8; + + for (int shell_index = 0; shell_index < shell_num; shell_index++) { + const int32_t nuc_index = nucleus_index[shell_index]; + const double dx = nucleus_coords[3*nuc_index + 0] - rx; + const double dy = nucleus_coords[3*nuc_index + 1] - ry; + const double dz = nucleus_coords[3*nuc_index + 2] - rz; + const double r = sqrt(dx*dx + dy*dy + dz*dz); + + // All possibly reported errors have been caught above + trexio_evaluate_nao_radial(shell_index, r, grid_start, + grid_size, grid_r, interpolator, normalization, &litude[shell_index]); + } + + return TREXIO_SUCCESS; +} + +trexio_exit_code trexio_evaluate_nao_radial_py (const int shell_index, + const double r, int64_t* grid_start, int n_grid_st, + int64_t* grid_size, int n_grid_si, double* grid_r, int n_grid_r, + double* interpolator, int n_interp, double* normalization, int n_norm, double* const amplitude) +{ + // Code needs to be copied because of the use of int64_t mandated by Python + // If a 64-bit version is implemented, this can be avoided + if (shell_index < 0) return TREXIO_INVALID_ARG_1; + if (r < 0) return TREXIO_INVALID_ARG_2; + if (grid_start == 0) return TREXIO_INVALID_ARG_3; + if (grid_size == 0) return TREXIO_INVALID_ARG_4; + if (grid_r == NULL) return TREXIO_INVALID_ARG_5; + if (interpolator == 0) return TREXIO_INVALID_ARG_6; + if (normalization == 0) return TREXIO_INVALID_ARG_7; + + const int32_t i0 = 4*grid_start[shell_index]; + + // Convert radius to logarithmic units + double r_log = 0; + trexio_convert_nao_radius_64 (r, grid_r + grid_start[shell_index], &r_log); + int32_t i_log = (int32_t) r_log; + if (i_log < 0) + i_log = 0; + else if (i_log >= grid_size[shell_index]) + return 0; // NAOs vanish at the boundary by definition + + double t = r_log - (double) i_log; + double val_spline = interpolator[i0 + 4*i_log + 0]; + val_spline += t * interpolator[i0 + 4*i_log + 1]; + val_spline += t * t * interpolator[i0 + 4*i_log + 2]; + val_spline += t * t * t * interpolator[i0 + 4*i_log + 3]; + + *amplitude = val_spline * normalization[shell_index] / r; + return TREXIO_SUCCESS; +} + +trexio_exit_code trexio_evaluate_nao_radial_all_py (const int32_t shell_num, + int64_t* nucleus_index, int n_nuc_id, double* nucleus_coords, int n_nuc_co, + int64_t* grid_start, int n_grid_st, int64_t* grid_size, int n_grid_si, + double* grid_r, int n_grid_r, double* interpolator, int n_interp, + double* normalization, int n_norm, + const double rx, const double ry, const double rz, double* const amplitudes, int amplitude_cnt) +{ + if (shell_num < 0) return TREXIO_INVALID_ARG_1; + if (nucleus_index == 0) return TREXIO_INVALID_ARG_2; + if (nucleus_coords == 0) return TREXIO_INVALID_ARG_3; + if (grid_start == 0) return TREXIO_INVALID_ARG_4; + if (grid_size == 0) return TREXIO_INVALID_ARG_5; + if (grid_r == NULL) return TREXIO_INVALID_ARG_6; + if (interpolator == 0) return TREXIO_INVALID_ARG_7; + if (normalization == 0) return TREXIO_INVALID_ARG_8; + + for (int shell_index = 0; shell_index < shell_num; shell_index++) { + const int32_t nuc_index = nucleus_index[shell_index]; + const double dx = nucleus_coords[3*nuc_index + 0] - rx; + const double dy = nucleus_coords[3*nuc_index + 1] - ry; + const double dz = nucleus_coords[3*nuc_index + 2] - rz; + const double r = sqrt(dx*dx + dy*dy + dz*dz); + + // All possibly reported errors have been caught above + trexio_evaluate_nao_radial_py(shell_index, r, grid_start, n_grid_st, + grid_size, n_grid_si, grid_r, n_grid_r, interpolator, n_interp, normalization, n_norm, &litudes[shell_index]); + } + + return TREXIO_SUCCESS; +} + #+end_src + + #+begin_src c :tangle trexio_private.h /* Popcount and trailz */ #if TREXIO_INT_SIZE == 64 @@ -5867,6 +6071,76 @@ interface end interface #+end_src + #+begin_src f90 :tangle prefix_fortran.f90 +interface + integer(trexio_exit_code) function trexio_convert_nao_radius_32(r, grid_r, log_r_out) & + bind(C, name="trexio_convert_nao_radius_32") + use, intrinsic :: iso_c_binding + import + real(c_float), intent(in), value :: r + real(c_float), intent(in) :: grid_r(*) + real(c_float), intent(out) :: log_r_out + end function trexio_convert_nao_radius_32 +end interface + #+end_src + + #+begin_src f90 :tangle prefix_fortran.f90 +interface + integer(trexio_exit_code) function trexio_convert_nao_radius_64(r, grid_r, log_r_out) & + bind(C, name="trexio_convert_nao_radius_64") + use, intrinsic :: iso_c_binding + import + real(c_double), intent(in), value :: r + real(c_double), intent(in) :: grid_r(*) + real(c_double), intent(out) :: log_r_out + end function trexio_convert_nao_radius_64 +end interface + #+end_src + + #+begin_src f90 :tangle prefix_fortran.f90 +interface + integer(trexio_exit_code) function trexio_evaluate_nao_radial (shell_index, r, & + grid_start, grid_size, grid_r, interpolator, normalization, amplitude) & + bind(C, name="trexio_evaluate_nao_radial") + use, intrinsic :: iso_c_binding + import + integer(c_int32_t), intent(in), value :: shell_index + real(c_double), intent(in), value :: r + integer(c_int32_t), intent(in) :: grid_start(*) + integer(c_int32_t), intent(in) :: grid_size(*) + real(c_double), intent(in) :: grid_r(*) + real(c_double), intent(in) :: interpolator(*) + real(c_double), intent(in) :: normalization(*) + real(c_double), intent(out) :: amplitude + end function trexio_evaluate_nao_radial +end interface + #+end_src + + #+begin_src f90 :tangle prefix_fortran.f90 +interface + integer(trexio_exit_code) function trexio_evaluate_nao_radial_all (shell_num, & + nucleus_index, nucleus_coords, & + grid_start, grid_size, grid_r, interpolator, & + normalization, rx, ry, rz, amplitudes) & + bind(C, name="trexio_evaluate_nao_radial_all") + use, intrinsic :: iso_c_binding + import + integer(c_int32_t), intent(in), value :: shell_num + integer(c_int32_t), intent(in) :: nucleus_index(*) + real(c_double), intent(in) :: nucleus_coords(*) + integer(c_int32_t), intent(in) :: grid_start(*) + integer(c_int32_t), intent(in) :: grid_size(*) + real(c_double), intent(in) :: grid_r(*) + real(c_double), intent(in) :: interpolator(*) + real(c_double), intent(in) :: normalization(*) + real(c_double), intent(in), value :: rx + real(c_double), intent(in), value :: ry + real(c_double), intent(in), value :: rz + real(c_double), intent(out) :: amplitudes(*) + end function trexio_evaluate_nao_radial_all +end interface + #+end_src + ** Python #+begin_src python :tangle basic_python.py @@ -5967,6 +6241,99 @@ def to_orbital_list_up_dn(n_int: int, determinant: list) -> tuple: return (orbital_list_up[0:occ_num_up], orbital_list_dn[0:occ_num_dn]) #+end_src + #+begin_src python :tangle basic_python.py +def convert_nao_radius(r: float, grid_r) -> float: + """Converts the radius r as a distance from a nucleus to the shell + s logarithmic grid. + + Input: + ~r~ - the radius to be converted + ~grid_r~ - The radial grid of the shell. Note that this is only the + grid of the shell of interest, not the array of all shells. + + Returns: + Float that gives the radius in the shell's logarithmic units + + Raises: + - Exception from AssertionError if TREXIO return code ~rc~ is different from TREXIO_SUCCESS and prints the error message using trexio_string_of_error. + - Exception from some other error (e.g. RuntimeError). + + """ + + rc, r_log = pytr.trexio_convert_nao_radius_py(r, grid_r) + + if rc != TREXIO_SUCCESS: + raise Exception(rc) + + return r_log + +def evaluate_nao_radial(shell_index, r, grid_start, grid_size, grid_r, interpolator, normalization) -> float: + """Evaluates the radial function of a given NAO shell at a distance from its center. + + Input: + ~shell_index~ - index of the shell of interest + ~r~ - distance from the shell center + ~grid_start~ - array of starting points of the individual splines + ~grid_size~ - array of number of data points per spline + ~grid_r~ - The radial grid of the shell. Note that this is only the + grid of the shell of interest, not the array of all shells. + ~interpolator~ - Interpolator of the spline. It is assumed to contain a cubic spline. + An interpolator for a kinetic energy spline can also be passed. + ~normalization~ - array of radial function normalization constants. + + Returns: + Value of the spline at the given radius + + Raises: + - Exception from AssertionError if TREXIO return code ~rc~ is different from TREXIO_SUCCESS and prints the error message using trexio_string_of_error. + - Exception from some other error (e.g. RuntimeError). + + """ + rc, amplitude = pytr.trexio_evaluate_nao_radial_py(shell_index, r, grid_start, grid_size, grid_r, interpolator.flatten(), normalization) + + if rc != TREXIO_SUCCESS: + raise Error(rc) + + return amplitude + +def evaluate_nao_radial_all(nucleus_index, nucleus_coords, grid_start, + grid_size, grid_r, interpolator, normalization, r) -> float: + """Evaluates the radial functions of all NAO shells at a given position in space. + + Input: + ~nucleus_index~ - array giving the centers of the NAO + ~nucleus_coords~ - array giving the coordinates of the NAO centers + ~grid_start~ - array of starting points of the individual splines + ~grid_size~ - array of number of data points per spline + ~grid_r~ - The radial grid of the shell. Note that this is only the + grid of the shell of interest, not the array of all shells. + ~interpolator~ - Interpolator of the spline. It is assumed to contain a cubic spline. + An interpolator for a kinetic energy spline can also be passed. + ~normalization~ - array of radial function normalization constants. + ~r~ - the position in space at which the functions are evaluated + + Returns: + Array of spline values at ~r~ + + Raises: + - Exception if ~r~ is not three dimensional + - Exception from AssertionError if TREXIO return code ~rc~ is different from TREXIO_SUCCESS and prints the error message using trexio_string_of_error. + - Exception from some other error (e.g. RuntimeError). + + """ + if len(r) != 3: + raise Error(TREXIO_INVALID_ARG7) + shell_cnt = len(nucleus_index) + rc, amplitudes = pytr.trexio_evaluate_nao_radial_all_py(shell_cnt, \ + nucleus_index, nucleus_coords.flatten(), grid_start, grid_size, grid_r, \ + interpolator.flatten(), normalization, r[0], r[1], r[2], shell_cnt) + + if rc != TREXIO_SUCCESS: + raise Exception(rc) + + return amplitudes + #+end_src + * Fortran helper/wrapper functions The function below adapts the original C-based ~trexio_open~ for Fortran. diff --git a/trex.org b/trex.org index ba901a0..857fdb0 100644 --- a/trex.org +++ b/trex.org @@ -254,6 +254,44 @@ All the basis set parameters are stored in one-dimensional arrays. +*** Numerical orbitals + + Trexio supports numerical atom centered orbitals. The implementation is + based on the approach of FHI-aims [Reference to paper]. These orbitals are + defined by the atom they are centered on, their angular momentum and a + radial function $R_s$, which is of the form + \[ + R_s(\mathbf{r}) = \mathcal{N}_s \frac{u_i(\mathbf{r})}{r^{-n_s}}. + \] + Where $u_i(\mathbf{r})$ is numerically tabulated on a dense logarithmic + grid. It is constructed to vanish for any $\mathbf{r}$ + outside of the grid. The reference points are stored in ~numgrid_r~ + and ~numgrid_phi~. Additionaly, a separate spline for the kinetic energy + can be stored in ~numgrid_kin~. The index of the first data point for + each shell is stored in ~numgrid_start~. + + FHI-aims uses cubic interpolation between the reference points. + The interpolation coefficients are given in the ~interpolator~ array. + The number of coefficients per interpolation point is stored as + ~interp_coeff_cnt~ so that other interpolation functions can be + implemented if needed. The ~interpolator_kin~ array provides a spline for + the kinetic energy. + The argument passed to the interpolants is on the logarithmic scale of + the reference points: If the argument is an integer $i$, the interpolant + will return the value of $u(\mathbf{r})$ at the $i$th reference point. + A radius is converted to this scale by (note the zero-indexing) + \[ + i_{log} = \frac{1}{c} \cdot \log{\frac{r}{r_0}} + \] + where + \[ + c = \log\left \frac{r_1}{r_0}\right. + \] + For convenience, this conversion and functions to evaluate the splines + are provided with trexio. Since these implementations are not adapted to + a specific software architecture, a programm using these orbitals should + reimplement them with consideration for its specific needs. + *** Plane waves A plane wave is defined as @@ -270,20 +308,30 @@ *** Data definitions #+NAME: basis - | Variable | Type | Dimensions | Description | - |-----------------+---------+---------------------+-----------------------------------------------------------------| - | ~type~ | ~str~ | | Type of basis set: "Gaussian", "Slater" or "PW" for plane waves | - | ~prim_num~ | ~dim~ | | Total number of primitives | - | ~shell_num~ | ~dim~ | | Total number of shells | - | ~nucleus_index~ | ~index~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and atomic indices | - | ~shell_ang_mom~ | ~int~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and angular momenta | - | ~shell_factor~ | ~float~ | ~(basis.shell_num)~ | Normalization factor of each shell ($\mathcal{N}_s$) | - | ~r_power~ | ~int~ | ~(basis.shell_num)~ | Power to which $r$ is raised ($n_s$) | - | ~shell_index~ | ~index~ | ~(basis.prim_num)~ | One-to-one correspondence between primitives and shell index | - | ~exponent~ | ~float~ | ~(basis.prim_num)~ | Exponents of the primitives ($\gamma_{ks}$) | - | ~coefficient~ | ~float~ | ~(basis.prim_num)~ | Coefficients of the primitives ($a_{ks}$) | - | ~prim_factor~ | ~float~ | ~(basis.prim_num)~ | Normalization coefficients for the primitives ($f_{ks}$) | - | ~e_cut~ | ~float~ | | Energy cut-off for plane-wave calculations | + | Variable | Type | Dimensions | Description | | + |--------------------+---------+----------------------------------------------+------------------------------------------------------------------------------+---| + | ~type~ | ~str~ | | Type of basis set: "Gaussian", "Slater", "Numerical" or "PW" for plane waves | | + | ~prim_num~ | ~dim~ | | Total number of primitives | | + | ~shell_num~ | ~dim~ | | Total number of shells | | + | ~numgrid_num~ | ~dim~ | | Total number of grid points for numerical orbitals | | + | ~interp_coeff_cnt~ | ~dim~ | | Number of coefficients for the numerical orbital interpolator | | + | ~nucleus_index~ | ~index~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and atomic indices | | + | ~shell_ang_mom~ | ~int~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and angular momenta | | + | ~shell_factor~ | ~float~ | ~(basis.shell_num)~ | Normalization factor of each shell ($\mathcal{N}_s$) | | + | ~r_power~ | ~int~ | ~(basis.shell_num)~ | Power to which $r$ is raised ($n_s$) | | + | ~numgrid_start~ | ~index~ | ~(basis.shell_num)~ | Index of the first data point for a given numerical orbital | | + | ~numgrid_size~ | ~dim~ | ~(basis.shell_num)~ | Number of data points per numerical orbital | | + | ~shell_index~ | ~index~ | ~(basis.prim_num)~ | One-to-one correspondence between primitives and shell index | | + | ~exponent~ | ~float~ | ~(basis.prim_num)~ | Exponents of the primitives ($\gamma_{ks}$) | | + | ~coefficient~ | ~float~ | ~(basis.prim_num)~ | Coefficients of the primitives ($a_{ks}$) | | + | ~prim_factor~ | ~float~ | ~(basis.prim_num)~ | Normalization coefficients for the primitives ($f_{ks}$) | | + | ~e_cut~ | ~float~ | | Energy cut-off for plane-wave calculations | | + | ~numgrid_radius~ | ~float~ | ~(basis.numgrid_num)~ | Radii of grid points for numerical orbitals | | + | ~numgrid_phi~ | ~float~ | ~(basis.numgrid_num)~ | Wave function values for numerical orbitals | | + | ~numgrid_kin~ | ~float~ | ~(basis.numgrid_num)~ | Kinetic energy of numerical orbitals | | + | ~interpolator~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital interpolation function | | + | ~interpolator_kin~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital kinetic energy interpolation function | | + #+CALL: json(data=basis, title="basis") @@ -291,20 +339,29 @@ #+RESULTS: :results: #+begin_src python :tangle trex.json - "basis": { - "type" : [ "str" , [] ] - , "prim_num" : [ "dim" , [] ] - , "shell_num" : [ "dim" , [] ] - , "nucleus_index" : [ "index", [ "basis.shell_num" ] ] - , "shell_ang_mom" : [ "int" , [ "basis.shell_num" ] ] - , "shell_factor" : [ "float", [ "basis.shell_num" ] ] - , "r_power" : [ "int" , [ "basis.shell_num" ] ] - , "shell_index" : [ "index", [ "basis.prim_num" ] ] - , "exponent" : [ "float", [ "basis.prim_num" ] ] - , "coefficient" : [ "float", [ "basis.prim_num" ] ] - , "prim_factor" : [ "float", [ "basis.prim_num" ] ] - , "e_cut" : [ "float", [] ] - } , + "basis": { + "type" : [ "str" , [] ] + , "prim_num" : [ "dim" , [] ] + , "shell_num" : [ "dim" , [] ] + , "numgrid_num" : [ "dim" , [] ] + , "interp_coeff_cnt" : [ "dim" , [] ] + , "nucleus_index" : [ "index", [ "basis.shell_num" ] ] + , "shell_ang_mom" : [ "int" , [ "basis.shell_num" ] ] + , "shell_factor" : [ "float", [ "basis.shell_num" ] ] + , "r_power" : [ "int" , [ "basis.shell_num" ] ] + , "numgrid_start" : [ "index", [ "basis.shell_num" ] ] + , "numgrid_size" : [ "dim" , [ "basis.shell_num" ] ] + , "shell_index" : [ "index", [ "basis.prim_num" ] ] + , "exponent" : [ "float", [ "basis.prim_num" ] ] + , "coefficient" : [ "float", [ "basis.prim_num" ] ] + , "prim_factor" : [ "float", [ "basis.prim_num" ] ] + , "e_cut" : [ "float", [] ] + , "numgrid_radius" : [ "float", [ "basis.numgrid_num" ] ] + , "numgrid_phi" : [ "float", [ "basis.numgrid_num" ] ] + , "numgrid_kin" : [ "float", [ "basis.numgrid_num" ] ] + , "interpolator" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] + , "interpolator_kin" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] + } , #+end_src :end: From b7d9f32bfb12e656fceaf83d7a7329f22dc6c2f7 Mon Sep 17 00:00:00 2001 From: joguenzl Date: Tue, 16 May 2023 14:13:45 +0200 Subject: [PATCH 2/6] Change NAO 'kinetic energy' fields to 'Laplacian' --- trex.org | 115 ++++++++++++++++++++++++++++++------------------------- 1 file changed, 62 insertions(+), 53 deletions(-) diff --git a/trex.org b/trex.org index 857fdb0..c7014f3 100644 --- a/trex.org +++ b/trex.org @@ -266,16 +266,19 @@ Where $u_i(\mathbf{r})$ is numerically tabulated on a dense logarithmic grid. It is constructed to vanish for any $\mathbf{r}$ outside of the grid. The reference points are stored in ~numgrid_r~ - and ~numgrid_phi~. Additionaly, a separate spline for the kinetic energy - can be stored in ~numgrid_kin~. The index of the first data point for - each shell is stored in ~numgrid_start~. + and ~numgrid_phi~. Additionaly, a separate spline for the gradient and Laplacian + can be stored in ~numgrid_grad~ and ~numgrid_lap~. The index of the first data point for + each shell is stored in ~numgrid_start~, the number of data points per spline + is stored in ~numgrid_size~ for convenience. - FHI-aims uses cubic interpolation between the reference points. - The interpolation coefficients are given in the ~interpolator~ array. - The number of coefficients per interpolation point is stored as - ~interp_coeff_cnt~ so that other interpolation functions can be - implemented if needed. The ~interpolator_kin~ array provides a spline for - the kinetic energy. + What kind of spline is used can be provided in the ~interpolator_kind~ field. + For example, FHI-aims uses a cubic spline, so the ~interpolator_kind~ is + \"polynomial\" and the ~interp_coeff_cnt~ is $4$. In this case, the first + interpolation coefficient per data point is the absolute term, the second is for + the linear term etc. + The interpolation coefficients for the wave function are given in the + ~interpolator_phi~ array. The ~interpolator_grad~ and ~interpolator_lap~ + arrays provide a spline for the gradient and Laplacian, respectively. The argument passed to the interpolants is on the logarithmic scale of the reference points: If the argument is an integer $i$, the interpolant will return the value of $u(\mathbf{r})$ at the $i$th reference point. @@ -308,29 +311,32 @@ *** Data definitions #+NAME: basis - | Variable | Type | Dimensions | Description | | - |--------------------+---------+----------------------------------------------+------------------------------------------------------------------------------+---| - | ~type~ | ~str~ | | Type of basis set: "Gaussian", "Slater", "Numerical" or "PW" for plane waves | | - | ~prim_num~ | ~dim~ | | Total number of primitives | | - | ~shell_num~ | ~dim~ | | Total number of shells | | - | ~numgrid_num~ | ~dim~ | | Total number of grid points for numerical orbitals | | - | ~interp_coeff_cnt~ | ~dim~ | | Number of coefficients for the numerical orbital interpolator | | - | ~nucleus_index~ | ~index~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and atomic indices | | - | ~shell_ang_mom~ | ~int~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and angular momenta | | - | ~shell_factor~ | ~float~ | ~(basis.shell_num)~ | Normalization factor of each shell ($\mathcal{N}_s$) | | - | ~r_power~ | ~int~ | ~(basis.shell_num)~ | Power to which $r$ is raised ($n_s$) | | - | ~numgrid_start~ | ~index~ | ~(basis.shell_num)~ | Index of the first data point for a given numerical orbital | | - | ~numgrid_size~ | ~dim~ | ~(basis.shell_num)~ | Number of data points per numerical orbital | | - | ~shell_index~ | ~index~ | ~(basis.prim_num)~ | One-to-one correspondence between primitives and shell index | | - | ~exponent~ | ~float~ | ~(basis.prim_num)~ | Exponents of the primitives ($\gamma_{ks}$) | | - | ~coefficient~ | ~float~ | ~(basis.prim_num)~ | Coefficients of the primitives ($a_{ks}$) | | - | ~prim_factor~ | ~float~ | ~(basis.prim_num)~ | Normalization coefficients for the primitives ($f_{ks}$) | | - | ~e_cut~ | ~float~ | | Energy cut-off for plane-wave calculations | | - | ~numgrid_radius~ | ~float~ | ~(basis.numgrid_num)~ | Radii of grid points for numerical orbitals | | - | ~numgrid_phi~ | ~float~ | ~(basis.numgrid_num)~ | Wave function values for numerical orbitals | | - | ~numgrid_kin~ | ~float~ | ~(basis.numgrid_num)~ | Kinetic energy of numerical orbitals | | - | ~interpolator~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital interpolation function | | - | ~interpolator_kin~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital kinetic energy interpolation function | | + | Variable | Type | Dimensions | Description | | + |---------------------+---------+----------------------------------------------+------------------------------------------------------------------------------+---| + | ~type~ | ~str~ | | Type of basis set: "Gaussian", "Slater", "Numerical" or "PW" for plane waves | | + | ~prim_num~ | ~dim~ | | Total number of primitives | | + | ~shell_num~ | ~dim~ | | Total number of shells | | + | ~numgrid_num~ | ~dim~ | | Total number of grid points for numerical orbitals | | + | ~interp_coeff_cnt~ | ~dim~ | | Number of coefficients for the numerical orbital interpolator | | + | ~nucleus_index~ | ~index~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and atomic indices | | + | ~shell_ang_mom~ | ~int~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and angular momenta | | + | ~shell_factor~ | ~float~ | ~(basis.shell_num)~ | Normalization factor of each shell ($\mathcal{N}_s$) | | + | ~r_power~ | ~int~ | ~(basis.shell_num)~ | Power to which $r$ is raised ($n_s$) | | + | ~numgrid_start~ | ~index~ | ~(basis.shell_num)~ | Index of the first data point for a given numerical orbital | | + | ~numgrid_size~ | ~dim~ | ~(basis.shell_num)~ | Number of data points per numerical orbital | | + | ~shell_index~ | ~index~ | ~(basis.prim_num)~ | One-to-one correspondence between primitives and shell index | | + | ~exponent~ | ~float~ | ~(basis.prim_num)~ | Exponents of the primitives ($\gamma_{ks}$) | | + | ~coefficient~ | ~float~ | ~(basis.prim_num)~ | Coefficients of the primitives ($a_{ks}$) | | + | ~prim_factor~ | ~float~ | ~(basis.prim_num)~ | Normalization coefficients for the primitives ($f_{ks}$) | | + | ~e_cut~ | ~float~ | | Energy cut-off for plane-wave calculations | | + | ~numgrid_radius~ | ~float~ | ~(basis.numgrid_num)~ | Radii of grid points for numerical orbitals | | + | ~numgrid_phi~ | ~float~ | ~(basis.numgrid_num)~ | Wave function values for numerical orbitals | | + | ~numgrid_grad~ | ~float~ | ~(basis.numgrid_num)~ | Radial gradient of numerical orbitals | | + | ~numgrid_lap~ | ~float~ | ~(basis.numgrid_num)~ | Laplacian of numerical orbitals | | + | ~interpolator_kind~ | ~str~ | | Kind of spline, e.g. "polynomial" | | + | ~interpolator_phi~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital interpolation function | | + | ~interpolator_grad~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital gradient interpolation function | | + | ~interpolator_lap~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital laplacian interpolation function | | @@ -340,27 +346,30 @@ :results: #+begin_src python :tangle trex.json "basis": { - "type" : [ "str" , [] ] - , "prim_num" : [ "dim" , [] ] - , "shell_num" : [ "dim" , [] ] - , "numgrid_num" : [ "dim" , [] ] - , "interp_coeff_cnt" : [ "dim" , [] ] - , "nucleus_index" : [ "index", [ "basis.shell_num" ] ] - , "shell_ang_mom" : [ "int" , [ "basis.shell_num" ] ] - , "shell_factor" : [ "float", [ "basis.shell_num" ] ] - , "r_power" : [ "int" , [ "basis.shell_num" ] ] - , "numgrid_start" : [ "index", [ "basis.shell_num" ] ] - , "numgrid_size" : [ "dim" , [ "basis.shell_num" ] ] - , "shell_index" : [ "index", [ "basis.prim_num" ] ] - , "exponent" : [ "float", [ "basis.prim_num" ] ] - , "coefficient" : [ "float", [ "basis.prim_num" ] ] - , "prim_factor" : [ "float", [ "basis.prim_num" ] ] - , "e_cut" : [ "float", [] ] - , "numgrid_radius" : [ "float", [ "basis.numgrid_num" ] ] - , "numgrid_phi" : [ "float", [ "basis.numgrid_num" ] ] - , "numgrid_kin" : [ "float", [ "basis.numgrid_num" ] ] - , "interpolator" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] - , "interpolator_kin" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] + "type" : [ "str" , [] ] + , "prim_num" : [ "dim" , [] ] + , "shell_num" : [ "dim" , [] ] + , "numgrid_num" : [ "dim" , [] ] + , "interp_coeff_cnt" : [ "dim" , [] ] + , "nucleus_index" : [ "index", [ "basis.shell_num" ] ] + , "shell_ang_mom" : [ "int" , [ "basis.shell_num" ] ] + , "shell_factor" : [ "float", [ "basis.shell_num" ] ] + , "r_power" : [ "int" , [ "basis.shell_num" ] ] + , "numgrid_start" : [ "index", [ "basis.shell_num" ] ] + , "numgrid_size" : [ "dim" , [ "basis.shell_num" ] ] + , "shell_index" : [ "index", [ "basis.prim_num" ] ] + , "exponent" : [ "float", [ "basis.prim_num" ] ] + , "coefficient" : [ "float", [ "basis.prim_num" ] ] + , "prim_factor" : [ "float", [ "basis.prim_num" ] ] + , "e_cut" : [ "float", [] ] + , "numgrid_radius" : [ "float", [ "basis.numgrid_num" ] ] + , "numgrid_phi" : [ "float", [ "basis.numgrid_num" ] ] + , "numgrid_grad" : [ "float", [ "basis.numgrid_num" ] ] + , "numgrid_lap" : [ "float", [ "basis.numgrid_num" ] ] + , "interpolator_kind" : [ "str" , [] ] + , "interpolator_phi" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] + , "interpolator_grad" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] + , "interpolator_lap" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] } , #+end_src :end: From c5516b0afb3d1923d263f58c76c6398b46e4a11d Mon Sep 17 00:00:00 2001 From: joguenzl Date: Thu, 25 May 2023 13:01:24 +0200 Subject: [PATCH 3/6] Fix upper vs. lower case --- trex.org | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/trex.org b/trex.org index c7014f3..84a7b69 100644 --- a/trex.org +++ b/trex.org @@ -273,7 +273,7 @@ What kind of spline is used can be provided in the ~interpolator_kind~ field. For example, FHI-aims uses a cubic spline, so the ~interpolator_kind~ is - \"polynomial\" and the ~interp_coeff_cnt~ is $4$. In this case, the first + \"Polynomial\" and the ~interp_coeff_cnt~ is $4$. In this case, the first interpolation coefficient per data point is the absolute term, the second is for the linear term etc. The interpolation coefficients for the wave function are given in the @@ -333,7 +333,7 @@ | ~numgrid_phi~ | ~float~ | ~(basis.numgrid_num)~ | Wave function values for numerical orbitals | | | ~numgrid_grad~ | ~float~ | ~(basis.numgrid_num)~ | Radial gradient of numerical orbitals | | | ~numgrid_lap~ | ~float~ | ~(basis.numgrid_num)~ | Laplacian of numerical orbitals | | - | ~interpolator_kind~ | ~str~ | | Kind of spline, e.g. "polynomial" | | + | ~interpolator_kind~ | ~str~ | | Kind of spline, e.g. "Polynomial" | | | ~interpolator_phi~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital interpolation function | | | ~interpolator_grad~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital gradient interpolation function | | | ~interpolator_lap~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital laplacian interpolation function | | From e62dde53be8183bd1a692ccef5ac9dccd98ebe68 Mon Sep 17 00:00:00 2001 From: joguenzl Date: Tue, 30 May 2023 10:24:58 +0200 Subject: [PATCH 4/6] Expand NAO documentation --- trex.org | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/trex.org b/trex.org index 84a7b69..f3d8e65 100644 --- a/trex.org +++ b/trex.org @@ -257,7 +257,9 @@ *** Numerical orbitals Trexio supports numerical atom centered orbitals. The implementation is - based on the approach of FHI-aims [Reference to paper]. These orbitals are + based on the approach of FHI-aims [Blum, V. et al; Ab initio molecular + simulations with numeric atom-centered orbitals; Computer Physics + Communications 2009]. These orbitals are defined by the atom they are centered on, their angular momentum and a radial function $R_s$, which is of the form \[ @@ -266,9 +268,20 @@ Where $u_i(\mathbf{r})$ is numerically tabulated on a dense logarithmic grid. It is constructed to vanish for any $\mathbf{r}$ outside of the grid. The reference points are stored in ~numgrid_r~ - and ~numgrid_phi~. Additionaly, a separate spline for the gradient and Laplacian - can be stored in ~numgrid_grad~ and ~numgrid_lap~. The index of the first data point for - each shell is stored in ~numgrid_start~, the number of data points per spline + and ~numgrid_phi~. Additionaly, a separate spline for the first and second + derivative of $u(\mathbf{r})$ can be stored in ~numgrid_grad~ and ~numgrid_lap~. + Storing them in this form allows to calculate the actual gradients and + Laplacian easily as follows: + + \[ + \grad_{x_i} \phi = \frac{x_i}{r^2}\left u\prime(r) - \frac{u(r)}{r}\right + \] + \[ + \Delta \phi = \frac{1}{r^3}\left x^2 u''(r) + \left3x^2-r^2\right \left \frac{u(r)}{r^2} - \frac{u'(r)}{r}\right \right + \] + + The index of the first data point for each shell is stored in + ~numgrid_start~, the number of data points per spline is stored in ~numgrid_size~ for convenience. What kind of spline is used can be provided in the ~interpolator_kind~ field. From a8fec4d6c75a9c024cf84f758ca4df9e72be3c83 Mon Sep 17 00:00:00 2001 From: joguenzl Date: Mon, 5 Jun 2023 16:21:09 +0200 Subject: [PATCH 5/6] Incorporate PR feedback --- src/templates_front/templator_front.org | 46 ++++++++++++++++--------- 1 file changed, 30 insertions(+), 16 deletions(-) diff --git a/src/templates_front/templator_front.org b/src/templates_front/templator_front.org index bc2eae4..af771ae 100644 --- a/src/templates_front/templator_front.org +++ b/src/templates_front/templator_front.org @@ -42,7 +42,6 @@ typedef int32_t trexio_exit_code; #include "config.h" #endif -#define _USE_MATH_DEFINES #include #include #include @@ -5818,6 +5817,7 @@ trexio_convert_nao_radius_32 (const float r, const float* const grid_r, float* c { if (r < 0) return TREXIO_INVALID_ARG_1; if (grid_r == NULL) return TREXIO_INVALID_ARG_2; + if (log_r_out == NULL) return TREXIO_INVALID_ARG_3; *log_r_out = log(r / grid_r[0]) / log(grid_r[1] / grid_r[0]); @@ -5829,6 +5829,7 @@ trexio_convert_nao_radius_64 (const double r, const double* const grid_r, double { if (r < 0) return TREXIO_INVALID_ARG_1; if (grid_r == NULL) return TREXIO_INVALID_ARG_2; + if (log_r_out == NULL) return TREXIO_INVALID_ARG_3; *log_r_out = log(r / grid_r[0]) / log(grid_r[1] / grid_r[0]); @@ -5840,6 +5841,7 @@ trexio_convert_nao_radius_py (const double r, double* grid_r, int32_t n_grid, do { if (r < 0) return TREXIO_INVALID_ARG_1; if (grid_r == NULL) return TREXIO_INVALID_ARG_2; + if (log_r_out == NULL) return TREXIO_INVALID_ARG_3; *log_r_out = log(r / grid_r[0]) / log(grid_r[1] / grid_r[0]); @@ -5864,7 +5866,7 @@ trexio_evaluate_nao_radial (const int32_t shell_index, const double r, const int const int32_t i0 = 4*grid_start[shell_index]; // Convert radius to logarithmic units - double r_log = 0; + double r_log = 0.0; trexio_convert_nao_radius_64 (r, grid_r + grid_start[shell_index], &r_log); int32_t i_log = (int32_t) r_log; if (i_log < 0) @@ -5896,6 +5898,8 @@ trexio_evaluate_nao_radial_all (const int32_t shell_num, const int32_t* const nu if (interpolator == 0) return TREXIO_INVALID_ARG_7; if (normalization == 0) return TREXIO_INVALID_ARG_8; + trexio_exit_code rc; + for (int shell_index = 0; shell_index < shell_num; shell_index++) { const int32_t nuc_index = nucleus_index[shell_index]; const double dx = nucleus_coords[3*nuc_index + 0] - rx; @@ -5904,8 +5908,11 @@ trexio_evaluate_nao_radial_all (const int32_t shell_num, const int32_t* const nu const double r = sqrt(dx*dx + dy*dy + dz*dz); // All possibly reported errors have been caught above - trexio_evaluate_nao_radial(shell_index, r, grid_start, + rc = trexio_evaluate_nao_radial(shell_index, r, grid_start, grid_size, grid_r, interpolator, normalization, &litude[shell_index]); + + if (rc != TREXIO_SUCCESS) + return rc; } return TREXIO_SUCCESS; @@ -5929,13 +5936,16 @@ trexio_exit_code trexio_evaluate_nao_radial_py (const int shell_index, const int32_t i0 = 4*grid_start[shell_index]; // Convert radius to logarithmic units - double r_log = 0; + double r_log = 0.0; trexio_convert_nao_radius_64 (r, grid_r + grid_start[shell_index], &r_log); int32_t i_log = (int32_t) r_log; - if (i_log < 0) - i_log = 0; - else if (i_log >= grid_size[shell_index]) - return 0; // NAOs vanish at the boundary by definition + if (i_log < 0) { + *amplitude = interpolator[i0] * normalization[shell_index] / r; + return TREXIO_SUCCESS; + } else if (i_log >= grid_size[shell_index]) { + *amplitude = 0.0; + return TREXIO_SUCCESS; // NAOs vanish at the boundary by definition + } double t = r_log - (double) i_log; double val_spline = interpolator[i0 + 4*i_log + 0]; @@ -5963,6 +5973,8 @@ trexio_exit_code trexio_evaluate_nao_radial_all_py (const int32_t shell_num, if (interpolator == 0) return TREXIO_INVALID_ARG_7; if (normalization == 0) return TREXIO_INVALID_ARG_8; + trexio_exit_code rc; + for (int shell_index = 0; shell_index < shell_num; shell_index++) { const int32_t nuc_index = nucleus_index[shell_index]; const double dx = nucleus_coords[3*nuc_index + 0] - rx; @@ -5971,8 +5983,10 @@ trexio_exit_code trexio_evaluate_nao_radial_all_py (const int32_t shell_num, const double r = sqrt(dx*dx + dy*dy + dz*dz); // All possibly reported errors have been caught above - trexio_evaluate_nao_radial_py(shell_index, r, grid_start, n_grid_st, + rc = trexio_evaluate_nao_radial_py(shell_index, r, grid_start, n_grid_st, grid_size, n_grid_si, grid_r, n_grid_r, interpolator, n_interp, normalization, n_norm, &litudes[shell_index]); + if (rc != TREXIO_SUCCESS) + return rc; } return TREXIO_SUCCESS; @@ -6298,7 +6312,7 @@ def convert_nao_radius(r: float, grid_r) -> float: rc, r_log = pytr.trexio_convert_nao_radius_py(r, grid_r) if rc != TREXIO_SUCCESS: - raise Exception(rc) + raise Error(rc) return r_log @@ -6320,8 +6334,8 @@ def evaluate_nao_radial(shell_index, r, grid_start, grid_size, grid_r, interpola Value of the spline at the given radius Raises: - - Exception from AssertionError if TREXIO return code ~rc~ is different from TREXIO_SUCCESS and prints the error message using trexio_string_of_error. - - Exception from some other error (e.g. RuntimeError). + - Error from AssertionError if TREXIO return code ~rc~ is different from TREXIO_SUCCESS and prints the error message using trexio_string_of_error. + - Error from some other error (e.g. RuntimeError). """ rc, amplitude = pytr.trexio_evaluate_nao_radial_py(shell_index, r, grid_start, grid_size, grid_r, interpolator.flatten(), normalization) @@ -6351,9 +6365,9 @@ def evaluate_nao_radial_all(nucleus_index, nucleus_coords, grid_start, Array of spline values at ~r~ Raises: - - Exception if ~r~ is not three dimensional - - Exception from AssertionError if TREXIO return code ~rc~ is different from TREXIO_SUCCESS and prints the error message using trexio_string_of_error. - - Exception from some other error (e.g. RuntimeError). + - Error if ~r~ is not three dimensional + - Error from AssertionError if TREXIO return code ~rc~ is different from TREXIO_SUCCESS and prints the error message using trexio_string_of_error. + - Error from some other error (e.g. RuntimeError). """ if len(r) != 3: @@ -6364,7 +6378,7 @@ def evaluate_nao_radial_all(nucleus_index, nucleus_coords, grid_start, interpolator.flatten(), normalization, r[0], r[1], r[2], shell_cnt) if rc != TREXIO_SUCCESS: - raise Exception(rc) + raise Error(rc) return amplitudes #+end_src From 1faa478025df4a8e6cf5f99e487dca8b454c38c4 Mon Sep 17 00:00:00 2001 From: joguenzl Date: Mon, 5 Jun 2023 16:59:20 +0200 Subject: [PATCH 6/6] Rename numgrid -> nao_grid --- trex.org | 50 +++++++++++++++++++++++++------------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/trex.org b/trex.org index 2431903..769a172 100644 --- a/trex.org +++ b/trex.org @@ -269,9 +269,9 @@ \] Where $u_i(\mathbf{r})$ is numerically tabulated on a dense logarithmic grid. It is constructed to vanish for any $\mathbf{r}$ - outside of the grid. The reference points are stored in ~numgrid_r~ - and ~numgrid_phi~. Additionaly, a separate spline for the first and second - derivative of $u(\mathbf{r})$ can be stored in ~numgrid_grad~ and ~numgrid_lap~. + outside of the grid. The reference points are stored in ~nao_grid_r~ + and ~nao_grid_phi~. Additionaly, a separate spline for the first and second + derivative of $u(\mathbf{r})$ can be stored in ~nao_grid_grad~ and ~nao_grid_lap~. Storing them in this form allows to calculate the actual gradients and Laplacian easily as follows: @@ -283,8 +283,8 @@ \] The index of the first data point for each shell is stored in - ~numgrid_start~, the number of data points per spline - is stored in ~numgrid_size~ for convenience. + ~nao_grid_start~, the number of data points per spline + is stored in ~nao_grid_size~ for convenience. What kind of spline is used can be provided in the ~interpolator_kind~ field. For example, FHI-aims uses a cubic spline, so the ~interpolator_kind~ is @@ -331,27 +331,27 @@ | ~type~ | ~str~ | | Type of basis set: "Gaussian", "Slater", "Numerical" or "PW" for plane waves | | | ~prim_num~ | ~dim~ | | Total number of primitives | | | ~shell_num~ | ~dim~ | | Total number of shells | | - | ~numgrid_num~ | ~dim~ | | Total number of grid points for numerical orbitals | | + | ~nao_grid_num~ | ~dim~ | | Total number of grid points for numerical orbitals | | | ~interp_coeff_cnt~ | ~dim~ | | Number of coefficients for the numerical orbital interpolator | | | ~nucleus_index~ | ~index~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and atomic indices | | | ~shell_ang_mom~ | ~int~ | ~(basis.shell_num)~ | One-to-one correspondence between shells and angular momenta | | | ~shell_factor~ | ~float~ | ~(basis.shell_num)~ | Normalization factor of each shell ($\mathcal{N}_s$) | | | ~r_power~ | ~int~ | ~(basis.shell_num)~ | Power to which $r$ is raised ($n_s$) | | - | ~numgrid_start~ | ~index~ | ~(basis.shell_num)~ | Index of the first data point for a given numerical orbital | | - | ~numgrid_size~ | ~dim~ | ~(basis.shell_num)~ | Number of data points per numerical orbital | | + | ~nao_grid_start~ | ~index~ | ~(basis.shell_num)~ | Index of the first data point for a given numerical orbital | | + | ~nao_grid_size~ | ~dim~ | ~(basis.shell_num)~ | Number of data points per numerical orbital | | | ~shell_index~ | ~index~ | ~(basis.prim_num)~ | One-to-one correspondence between primitives and shell index | | | ~exponent~ | ~float~ | ~(basis.prim_num)~ | Exponents of the primitives ($\gamma_{ks}$) | | | ~coefficient~ | ~float~ | ~(basis.prim_num)~ | Coefficients of the primitives ($a_{ks}$) | | | ~prim_factor~ | ~float~ | ~(basis.prim_num)~ | Normalization coefficients for the primitives ($f_{ks}$) | | | ~e_cut~ | ~float~ | | Energy cut-off for plane-wave calculations | | - | ~numgrid_radius~ | ~float~ | ~(basis.numgrid_num)~ | Radii of grid points for numerical orbitals | | - | ~numgrid_phi~ | ~float~ | ~(basis.numgrid_num)~ | Wave function values for numerical orbitals | | - | ~numgrid_grad~ | ~float~ | ~(basis.numgrid_num)~ | Radial gradient of numerical orbitals | | - | ~numgrid_lap~ | ~float~ | ~(basis.numgrid_num)~ | Laplacian of numerical orbitals | | + | ~nao_grid_radius~ | ~float~ | ~(basis.nao_grid_num)~ | Radii of grid points for numerical orbitals | | + | ~nao_grid_phi~ | ~float~ | ~(basis.nao_grid_num)~ | Wave function values for numerical orbitals | | + | ~nao_grid_grad~ | ~float~ | ~(basis.nao_grid_num)~ | Radial gradient of numerical orbitals | | + | ~nao_grid_lap~ | ~float~ | ~(basis.nao_grid_num)~ | Laplacian of numerical orbitals | | | ~interpolator_kind~ | ~str~ | | Kind of spline, e.g. "Polynomial" | | - | ~interpolator_phi~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital interpolation function | | - | ~interpolator_grad~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital gradient interpolation function | | - | ~interpolator_lap~ | ~float~ | ~(basis.interp_coeff_cnt,basis.numgrid_num)~ | Coefficients for numerical orbital laplacian interpolation function | | + | ~interpolator_phi~ | ~float~ | ~(basis.interp_coeff_cnt,basis.nao_grid_num)~ | Coefficients for numerical orbital interpolation function | | + | ~interpolator_grad~ | ~float~ | ~(basis.interp_coeff_cnt,basis.nao_grid_num)~ | Coefficients for numerical orbital gradient interpolation function | | + | ~interpolator_lap~ | ~float~ | ~(basis.interp_coeff_cnt,basis.nao_grid_num)~ | Coefficients for numerical orbital laplacian interpolation function | | @@ -364,27 +364,27 @@ "type" : [ "str" , [] ] , "prim_num" : [ "dim" , [] ] , "shell_num" : [ "dim" , [] ] - , "numgrid_num" : [ "dim" , [] ] + , "nao_grid_num" : [ "dim" , [] ] , "interp_coeff_cnt" : [ "dim" , [] ] , "nucleus_index" : [ "index", [ "basis.shell_num" ] ] , "shell_ang_mom" : [ "int" , [ "basis.shell_num" ] ] , "shell_factor" : [ "float", [ "basis.shell_num" ] ] , "r_power" : [ "int" , [ "basis.shell_num" ] ] - , "numgrid_start" : [ "index", [ "basis.shell_num" ] ] - , "numgrid_size" : [ "dim" , [ "basis.shell_num" ] ] + , "nao_grid_start" : [ "index", [ "basis.shell_num" ] ] + , "nao_grid_size" : [ "dim" , [ "basis.shell_num" ] ] , "shell_index" : [ "index", [ "basis.prim_num" ] ] , "exponent" : [ "float", [ "basis.prim_num" ] ] , "coefficient" : [ "float", [ "basis.prim_num" ] ] , "prim_factor" : [ "float", [ "basis.prim_num" ] ] , "e_cut" : [ "float", [] ] - , "numgrid_radius" : [ "float", [ "basis.numgrid_num" ] ] - , "numgrid_phi" : [ "float", [ "basis.numgrid_num" ] ] - , "numgrid_grad" : [ "float", [ "basis.numgrid_num" ] ] - , "numgrid_lap" : [ "float", [ "basis.numgrid_num" ] ] + , "nao_grid_radius" : [ "float", [ "basis.nao_grid_num" ] ] + , "nao_grid_phi" : [ "float", [ "basis.nao_grid_num" ] ] + , "nao_grid_grad" : [ "float", [ "basis.nao_grid_num" ] ] + , "nao_grid_lap" : [ "float", [ "basis.nao_grid_num" ] ] , "interpolator_kind" : [ "str" , [] ] - , "interpolator_phi" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] - , "interpolator_grad" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] - , "interpolator_lap" : [ "float", [ "basis.numgrid_num", "basis.interp_coeff_cnt" ] ] + , "interpolator_phi" : [ "float", [ "basis.nao_grid_num", "basis.interp_coeff_cnt" ] ] + , "interpolator_grad" : [ "float", [ "basis.nao_grid_num", "basis.interp_coeff_cnt" ] ] + , "interpolator_lap" : [ "float", [ "basis.nao_grid_num", "basis.interp_coeff_cnt" ] ] } , #+end_src :end: