mirror of
https://github.com/TREX-CoE/trexio.git
synced 2025-01-08 20:33:36 +01:00
Merge pull request #124 from joguenzl/master
Support for numerical atomic orbitals
This commit is contained in:
GitHub
commit
d011447f6f
@ -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
|
||||
|
||||
@ -42,6 +42,7 @@ typedef int32_t trexio_exit_code;
|
||||
#include "config.h"
|
||||
#endif
|
||||
|
||||
#include <math.h>
|
||||
#include <pthread.h>
|
||||
#include <assert.h>
|
||||
#include <stdlib.h>
|
||||
@ -5628,6 +5629,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
|
||||
@ -5646,6 +5653,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
|
||||
|
||||
|
||||
@ -5783,7 +5819,189 @@ 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;
|
||||
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]);
|
||||
|
||||
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;
|
||||
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]);
|
||||
|
||||
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;
|
||||
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]);
|
||||
|
||||
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.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;
|
||||
|
||||
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;
|
||||
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
|
||||
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;
|
||||
}
|
||||
|
||||
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.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) {
|
||||
*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];
|
||||
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;
|
||||
|
||||
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;
|
||||
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
|
||||
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;
|
||||
}
|
||||
#+end_src
|
||||
|
||||
#+begin_src c :tangle trexio_private.h
|
||||
/* Popcount and trailz */
|
||||
#if TREXIO_INT_SIZE == 64
|
||||
|
||||
@ -5910,6 +6128,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
|
||||
@ -6010,6 +6298,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 Error(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:
|
||||
- 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)
|
||||
|
||||
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:
|
||||
- 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:
|
||||
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 Error(rc)
|
||||
|
||||
return amplitudes
|
||||
#+end_src
|
||||
|
||||
* Fortran helper/wrapper functions
|
||||
|
||||
The function below adapts the original C-based ~trexio_open~ for Fortran.
|
||||
|
||||
135
trex.org
135
trex.org
@ -256,6 +256,60 @@
|
||||
|
||||
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 [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
|
||||
\[
|
||||
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 ~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:
|
||||
|
||||
\[
|
||||
\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
|
||||
~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
|
||||
\"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.
|
||||
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
|
||||
@ -272,20 +326,33 @@
|
||||
*** 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 | |
|
||||
| ~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$) | |
|
||||
| ~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 | |
|
||||
| ~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.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 | |
|
||||
|
||||
|
||||
|
||||
#+CALL: json(data=basis, title="basis")
|
||||
@ -293,20 +360,32 @@
|
||||
#+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" , [] ]
|
||||
, "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" ] ]
|
||||
, "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", [] ]
|
||||
, "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.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:
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user