From ec7201783fbfd360d054891d7d093ec34577f98e Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 10 May 2022 19:18:19 +0200 Subject: [PATCH] Possibility to compute only values --- org/qmckl_ao.org | 1025 ++++++++++++++++++++++++++++++++++++++++++---- org/qmckl_mo.org | 508 ++++++++++++++++++++++- 2 files changed, 1448 insertions(+), 85 deletions(-) diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 03405bb..4120690 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -299,6 +299,8 @@ typedef struct qmckl_ao_basis_struct { uint64_t shell_vgl_date; double * restrict ao_vgl; uint64_t ao_vgl_date; + double * restrict ao_value; + uint64_t ao_value_date; int32_t uninitialized; bool provided; @@ -2490,8 +2492,10 @@ free(ao_factor_test); | ~primitive_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the primitives at current positions | | ~shell_vgl~ | ~double[point_num][5][shell_num]~ | Value, gradients, Laplacian of the primitives at current positions | | ~shell_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions | - | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | Value, gradients, Laplacian of the primitives at current positions | + | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | Value, gradients, Laplacian of the AOs at current positions | | ~ao_vgl_date~ | ~uint64_t~ | Last modification date of Value, gradients, Laplacian of the AOs at current positions | + | ~ao_value~ | ~double[point_num][ao_num]~ | Values of the the AOs at current positions | + | ~ao_value_date~ | ~uint64_t~ | Last modification date of the values of the AOs at current positions | |----------------------+-----------------------------------+----------------------------------------------------------------------------------------------| *** After initialization @@ -3022,7 +3026,7 @@ qmckl_get_ao_basis_ao_vgl (qmckl_context context, end interface #+end_src - Uses the give array to compute the VGL. + Uses the given array to compute the VGL. #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code @@ -3030,7 +3034,7 @@ qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context, double* const ao_vgl, const int64_t size_max); #+end_src - + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context, @@ -3088,6 +3092,133 @@ qmckl_get_ao_basis_ao_vgl_inplace (qmckl_context context, end interface #+end_src + + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_ao_basis_ao_value (qmckl_context context, + double* const ao_value, + const int64_t size_max); + #+end_src + + Returns the array of values of the atomic orbitals evaluated at + the current coordinates. See section [[Combining radial and polynomial parts]]. + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_ao_basis_ao_value (qmckl_context context, + double* const ao_value, + const int64_t size_max) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_ao_basis_ao_value", + NULL); + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ao_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + int64_t sze = ctx->ao_basis.ao_num * ctx->point.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_ao_basis_ao_value", + "input array too small"); + } + memcpy(ao_value, ctx->ao_basis.ao_value, (size_t) sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_get_ao_basis_ao_value (context, & + ao_value, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: ao_value(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function qmckl_get_ao_basis_ao_value + end interface + #+end_src + + Uses the given array to compute the value. + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_ao_basis_ao_value_inplace (qmckl_context context, + double* const ao_value, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_ao_basis_ao_value_inplace (qmckl_context context, + double* const ao_value, + const int64_t size_max) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_ao_basis_ao_value", + NULL); + } + + qmckl_exit_code rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + int64_t sze = ctx->ao_basis.ao_num * ctx->point.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_ao_basis_ao_value", + "input array too small"); + } + + rc = qmckl_context_touch(context); + if (rc != QMCKL_SUCCESS) return rc; + + double* old_array = ctx->ao_basis.ao_value; + + ctx->ao_basis.ao_value = ao_value; + + rc = qmckl_provide_ao_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + ctx->ao_basis.ao_value = old_array; + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_get_ao_basis_ao_value_inplace (context, & + ao_value, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: ao_value(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function qmckl_get_ao_basis_ao_value_inplace + end interface + #+end_src + * Radial part ** General functions for Gaussian basis functions @@ -5285,33 +5416,780 @@ for (int32_t ldl=3 ; ldl<=5 ; ++ldl) { #+end_src * Combining radial and polynomial parts + +** Values only + :PROPERTIES: + :Name: qmckl_compute_ao_value + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: +*** Unoptimized version + #+NAME: qmckl_ao_value_args_doc + | Variable | Type | In/Out | Description | + |-----------------------+-----------------------------------+--------+----------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~shell_num~ | ~int64_t~ | in | Number of shells | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~coord~ | ~double[3][point_num]~ | in | Coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | + | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | + | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | + | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | + | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | + | ~shell_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells | + | ~ao_value~ | ~double[point_num][ao_num]~ | out | Values of the AOs | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_ao_value_doc_f(context, & + ao_num, shell_num, point_num, nucl_num, & + coord, nucl_coord, nucleus_index, nucleus_shell_num, & + nucleus_range, nucleus_max_ang_mom, shell_ang_mom, & + ao_factor, shell_vgl, ao_value) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: ao_num + integer*8 , intent(in) :: shell_num + integer*8 , intent(in) :: point_num + integer*8 , intent(in) :: nucl_num + double precision , intent(in) :: coord(point_num,3) + double precision , intent(in) :: nucl_coord(nucl_num,3) + integer*8 , intent(in) :: nucleus_index(nucl_num) + integer*8 , intent(in) :: nucleus_shell_num(nucl_num) + double precision , intent(in) :: nucleus_range(nucl_num) + integer , intent(in) :: nucleus_max_ang_mom(nucl_num) + integer , intent(in) :: shell_ang_mom(shell_num) + double precision , intent(in) :: ao_factor(ao_num) + double precision , intent(in) :: shell_vgl(shell_num,5,point_num) + double precision , intent(out) :: ao_value(ao_num,point_num) + + double precision :: e_coord(3), n_coord(3) + integer*8 :: n_poly + integer :: l, il, k + integer*8 :: ipoint, inucl, ishell + integer*8 :: ishell_start, ishell_end + integer :: lstart(0:20) + double precision :: x, y, z, r2 + double precision :: cutoff + integer, external :: qmckl_ao_polynomial_vgl_doc_f + + double precision, allocatable :: poly_vgl(:,:) + integer , allocatable :: powers(:,:), ao_index(:) + + allocate(poly_vgl(5,ao_num), powers(3,ao_num), ao_index(ao_num)) + + ! Pre-computed data + do l=0,20 + lstart(l) = l*(l+1)*(l+2)/6 +1 + end do + + k=1 + do inucl=1,nucl_num + ishell_start = nucleus_index(inucl) + 1 + ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl) + do ishell = ishell_start, ishell_end + l = shell_ang_mom(ishell) + ao_index(ishell) = k + k = k + lstart(l+1) - lstart(l) + end do + end do + info = QMCKL_SUCCESS + + ! Don't compute polynomials when the radial part is zero. + cutoff = -dlog(1.d-12) + + do ipoint = 1, point_num + e_coord(1) = coord(ipoint,1) + e_coord(2) = coord(ipoint,2) + e_coord(3) = coord(ipoint,3) + do inucl=1,nucl_num + n_coord(1) = nucl_coord(inucl,1) + n_coord(2) = nucl_coord(inucl,2) + n_coord(3) = nucl_coord(inucl,3) + + ! Test if the point is in the range of the nucleus + x = e_coord(1) - n_coord(1) + y = e_coord(2) - n_coord(2) + z = e_coord(3) - n_coord(3) + + r2 = x*x + y*y + z*z + + if (r2 > cutoff*nucleus_range(inucl)) then + cycle + end if + + ! Compute polynomials + info = qmckl_ao_polynomial_vgl_doc_f(context, e_coord, n_coord, & + nucleus_max_ang_mom(inucl), n_poly, powers, 3_8, & + poly_vgl, 5_8) + + ! Loop over shells + ishell_start = nucleus_index(inucl) + 1 + ishell_end = nucleus_index(inucl) + nucleus_shell_num(inucl) + do ishell = ishell_start, ishell_end + k = ao_index(ishell) + l = shell_ang_mom(ishell) + do il = lstart(l), lstart(l+1)-1 + ! Value + ao_value(k,ipoint) = & + poly_vgl(1,il) * shell_vgl(ishell,1,ipoint) * ao_factor(k) + k = k+1 + end do + end do + end do + end do + + deallocate(poly_vgl, powers) +end function qmckl_compute_ao_value_doc_f + #+end_src + +*** HPC version + #+NAME: qmckl_ao_value_args_hpc_gaussian + | Variable | Type | In/Out | Description | + |-----------------------+-----------------------------+--------+----------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~shell_num~ | ~int64_t~ | in | Number of shells | + | ~prim_num~ | ~int64_t~ | in | Number of primitives | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~coord~ | ~double[3][point_num]~ | in | Coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | + | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | + | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | + | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | + | ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell | + | ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell | + | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | + | ~ao_expo~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | + | ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | + | ~ao_value~ | ~double[point_num][ao_num]~ | out | Values of the AOs | + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +#ifdef HAVE_HPC +qmckl_exit_code +qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context, + const int64_t ao_num, + const int64_t shell_num, + const int32_t* restrict prim_num_per_nucleus, + const int64_t point_num, + const int64_t nucl_num, + const double* restrict coord, + const double* restrict nucl_coord, + const int64_t* restrict nucleus_index, + const int64_t* restrict nucleus_shell_num, + const double* nucleus_range, + const int32_t* restrict nucleus_max_ang_mom, + const int32_t* restrict shell_ang_mom, + const double* restrict ao_factor, + const qmckl_matrix expo_per_nucleus, + const qmckl_tensor coef_per_nucleus, + double* restrict const ao_value ) +{ + int32_t lstart[32]; + for (int32_t l=0 ; l<32 ; ++l) { + lstart[l] = l*(l+1)*(l+2)/6; + } + + int64_t ao_index[shell_num+1]; + int64_t size_max = 0; + int64_t prim_max = 0; + int64_t shell_max = 0; + { + int64_t k=0; + for (int inucl=0 ; inucl < nucl_num ; ++inucl) { + prim_max = prim_num_per_nucleus[inucl] > prim_max ? + prim_num_per_nucleus[inucl] : prim_max; + shell_max = nucleus_shell_num[inucl] > shell_max ? + nucleus_shell_num[inucl] : shell_max; + const int64_t ishell_start = nucleus_index[inucl]; + const int64_t ishell_end = nucleus_index[inucl] + nucleus_shell_num[inucl]; + for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) { + const int l = shell_ang_mom[ishell]; + ao_index[ishell] = k; + k += lstart[l+1] - lstart[l]; + size_max = size_max < lstart[l+1] ? lstart[l+1] : size_max; + } + } + ao_index[shell_num] = ao_num+1; + } + + /* Don't compute polynomials when the radial part is zero. */ + double cutoff = -log(1.e-12); + +#ifdef HAVE_OPENMP +#pragma omp parallel +#endif + { + qmckl_exit_code rc; + double ar2[prim_max]; + int32_t powers[3*size_max]; + double poly_vgl[5*size_max]; + + double exp_mat[prim_max]; + double ce_mat[shell_max]; + + double coef_mat[nucl_num][shell_max][prim_max]; + for (int i=0 ; i cutoff * nucleus_range[inucl]) { + continue; + } + + int64_t n_poly; + switch (nucleus_max_ang_mom[inucl]) { + case 0: + break; + + case 1: + poly_vgl[0] = 0.; + poly_vgl[1] = x; + poly_vgl[2] = y; + poly_vgl[3] = z; + break; + + case 2: + poly_vgl[0] = 0.; + poly_vgl[1] = x; + poly_vgl[2] = y; + poly_vgl[3] = z; + poly_vgl[4] = x*x; + poly_vgl[5] = x*y; + poly_vgl[6] = x*z; + poly_vgl[7] = y*y; + poly_vgl[8] = y*z; + poly_vgl[9] = z*z; + break; + + default: + rc = qmckl_ao_polynomial_transp_vgl_hpc(context, e_coord, n_coord, + nucleus_max_ang_mom[inucl], + &n_poly, powers, (int64_t) 3, + poly_vgl, size_max); + assert (rc == QMCKL_SUCCESS); + break; + } + + /* Compute all exponents */ + + int64_t nidx = 0; + for (int64_t iprim = 0 ; iprim < prim_num_per_nucleus[inucl] ; ++iprim) { + const double v = qmckl_mat(expo_per_nucleus, iprim, inucl) * r2; + if (v <= cutoff) { + ar2[iprim] = v; + ++nidx; + } else { + break; + } + } + + for (int64_t iprim = 0 ; iprim < nidx ; ++iprim) { + exp_mat[iprim] = exp(-ar2[iprim]); + } + + for (int i=0 ; i 0) { + const double* restrict f = ao_factor + k; + const int64_t idx = lstart[l]; + + poly_vgl_1 = &(poly_vgl[idx]); + + switch (n) { + case(1): + ao_value_1[0] = s1 * f[0]; + break; + case (3): +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (int il=0 ; il<3 ; ++il) { + ao_value_1[il] = poly_vgl_1[il] * s1 * f[il]; + } + break; + case(6): +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (int il=0 ; il<6 ; ++il) { + ao_value_1[il] = poly_vgl_1[il] * s1 * f[il]; + } + break; + default: +#ifdef HAVE_OPENMP +#pragma omp simd simdlen(8) +#endif + for (int il=0 ; ilao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_value", + NULL); + } + + /* Compute if necessary */ + if (ctx->point.date > ctx->ao_basis.ao_value_date) { + + qmckl_exit_code rc; + + /* Provide required data */ +#ifndef HAVE_HPC + rc = qmckl_provide_ao_basis_shell_vgl(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_shell_vgl", NULL); + } +#endif + + /* Allocate array */ + if (ctx->ao_basis.ao_value == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->ao_basis.ao_num * 5 * ctx->point.num * sizeof(double); + double* ao_value = (double*) qmckl_malloc(context, mem_info); + + if (ao_value == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_ao_basis_ao_value", + NULL); + } + ctx->ao_basis.ao_value = ao_value; + } + + if (ctx->ao_basis.ao_vgl_date == ctx->point.date) { + + // ao_vgl has been computed at this step: Just copy the data. + + double * v = &(ctx->ao_basis.ao_value[0]); + double * vgl = &(ctx->ao_basis.ao_vgl[0]); + for (int i=0 ; ipoint.num ; ++i) { + for (int k=0 ; kao_basis.ao_num ; ++k) { + v[k] = vgl[k]; + } + v += ctx->ao_basis.ao_num; + vgl += ctx->ao_basis.ao_num * 5; + } + + } else { + +#ifdef HAVE_HPC + if (ctx->ao_basis.type == 'G') { + rc = qmckl_compute_ao_value_hpc_gaussian(context, + ctx->ao_basis.ao_num, + ctx->ao_basis.shell_num, + ctx->ao_basis.prim_num_per_nucleus, + ctx->point.num, + ctx->nucleus.num, + ctx->point.coord.data, + ctx->nucleus.coord.data, + ctx->ao_basis.nucleus_index, + ctx->ao_basis.nucleus_shell_num, + ctx->ao_basis.nucleus_range, + ctx->ao_basis.nucleus_max_ang_mom, + ctx->ao_basis.shell_ang_mom, + ctx->ao_basis.ao_factor, + ctx->ao_basis.expo_per_nucleus, + ctx->ao_basis.coef_per_nucleus, + ctx->ao_basis.ao_value); + /* + } else if (ctx->ao_basis.type == 'S') { + rc = qmck_compute_ao_value_hpc_slater(context, + ctx->ao_basis.ao_num, + ctx->ao_basis.shell_num, + ctx->ao_basis.prim_num, + ctx->point.num, + ctx->nucleus.num, + ctx->point.coord.data, + ctx->nucleus.coord.data, + ctx->ao_basis.nucleus_index, + ctx->ao_basis.nucleus_shell_num, + ctx->ao_basis.nucleus_range, + ctx->ao_basis.nucleus_max_ang_mom, + ctx->ao_basis.shell_ang_mom, + ctx->ao_basis.shell_prim_index, + ctx->ao_basis.shell_prim_num, + ctx->ao_basis.ao_factor, + ctx->ao_basis.exponent, + ctx->ao_basis.coefficient_normalized, + ctx->ao_basis.ao_value); + ,*/ + } else { + rc = qmckl_compute_ao_value_doc(context, + ctx->ao_basis.ao_num, + ctx->ao_basis.shell_num, + ctx->point.num, + ctx->nucleus.num, + ctx->point.coord.data, + ctx->nucleus.coord.data, + ctx->ao_basis.nucleus_index, + ctx->ao_basis.nucleus_shell_num, + ctx->ao_basis.nucleus_range, + ctx->ao_basis.nucleus_max_ang_mom, + ctx->ao_basis.shell_ang_mom, + ctx->ao_basis.ao_factor, + ctx->ao_basis.shell_vgl, + ctx->ao_basis.ao_value); + } +#else + rc = qmckl_compute_ao_value_doc(context, + ctx->ao_basis.ao_num, + ctx->ao_basis.shell_num, + ctx->point.num, + ctx->nucleus.num, + ctx->point.coord.data, + ctx->nucleus.coord.data, + ctx->ao_basis.nucleus_index, + ctx->ao_basis.nucleus_shell_num, + ctx->ao_basis.nucleus_range, + ctx->ao_basis.nucleus_max_ang_mom, + ctx->ao_basis.shell_ang_mom, + ctx->ao_basis.ao_factor, + ctx->ao_basis.shell_vgl, + ctx->ao_basis.ao_value); +#endif + if (rc != QMCKL_SUCCESS) { + return rc; + } + + } + + ctx->ao_basis.ao_value_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +**** Test :noexport: + + #+begin_src python :results output :exports none +import numpy as np +from math import sqrt + +h0 = 1.e-4 +def f(a,x,y): + return np.sum( [c * np.exp( -b*(np.linalg.norm(x-y))**2) for b,c in a] ) + +elec_26_w1 = np.array( [ 1.49050402641, 2.90106987953, -1.05920815468 ] ) +elec_15_w2 = np.array( [ -2.20180344582,-1.9113150239, 2.2193744778600002 ] ) +nucl_1 = np.array( [ -2.302574592081335e+00, -3.542027060505035e-01, -5.334129934317614e-02] ) + +#double ao_value[prim_num][5][elec_num]; +x = elec_26_w1 ; y = nucl_1 +a = [( 4.0382999999999998e+02, 1.4732000000000000e-03 * 5.9876577632594533e+04), + ( 1.2117000000000000e+02, 1.2672500000000000e-02 * 7.2836806319891484e+03), + ( 4.6344999999999999e+01, 5.8045100000000002e-02 * 1.3549226646722386e+03), + ( 1.9721000000000000e+01, 1.7051030000000000e-01 * 3.0376315094739988e+02), + ( 8.8623999999999992e+00, 3.1859579999999998e-01 * 7.4924579607137730e+01), + ( 3.9962000000000000e+00, 3.8450230000000002e-01 * 1.8590543353806009e+01), + ( 1.7636000000000001e+00, 2.7377370000000001e-01 * 4.4423176930919421e+00), + ( 7.0618999999999998e-01, 7.4396699999999996e-02 * 8.9541051939952665e-01)] + +norm = sqrt(3.) +# x^2 * g(r) +print ( "[26][0][219] : %25.15e"%(fx(a,x,y)) ) +print ( "[26][1][219] : %25.15e"%(df(a,x,y,1)) ) +print ( "[26][2][219] : %25.15e"%(df(a,x,y,2)) ) +print ( "[26][3][219] : %25.15e"%(df(a,x,y,3)) ) +print ( "[26][4][219] : %25.15e"%(lf(a,x,y)) ) + +print ( "[26][0][220] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[1] - y[1]) )) +print ( "[26][1][220] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[1] - y[1]) + norm*f(a,x,y) * (x[1] - y[1])) ) + +print ( "[26][0][221] : %25.15e"%(norm*f(a,x,y) * (x[0] - y[0]) * (x[2] - y[2])) ) +print ( "[26][1][221] : %25.15e"%(norm*df(a,x,y,1)* (x[0] - y[0]) * (x[2] - y[2]) + norm*f(a,x,y) * (x[2] - y[2])) ) + +print ( "[26][0][222] : %25.15e"%(f(a,x,y) * (x[1] - y[1]) * (x[1] - y[1])) ) +print ( "[26][1][222] : %25.15e"%(df(a,x,y,1)* (x[1] - y[1]) * (x[1] - y[1])) ) + +print ( "[26][0][223] : %25.15e"%(norm*f(a,x,y) * (x[1] - y[1]) * (x[2] - y[2])) ) +print ( "[26][1][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2])) ) + +print ( "[26][0][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) ) +print ( "[26][1][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) ) + + #+end_src + + #+RESULTS: + + #+begin_src c :tangle (eval c_test) :exports none +{ +#define walk_num 1 // chbrclf_walk_num +#define elec_num chbrclf_elec_num +#define shell_num chbrclf_shell_num +#define ao_num chbrclf_ao_num + +int64_t elec_up_num = chbrclf_elec_up_num; +int64_t elec_dn_num = chbrclf_elec_dn_num; +double* elec_coord = &(chbrclf_elec_coord[0][0][0]); + +rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_electron_walk_num (context, walk_num); +assert (rc == QMCKL_SUCCESS); + +assert(qmckl_electron_provided(context)); + +rc = qmckl_set_electron_coord (context, 'N', elec_coord, walk_num*elec_num*3); +assert(rc == QMCKL_SUCCESS); + + +double ao_value[elec_num][ao_num]; + +rc = qmckl_get_ao_basis_ao_value(context, &(ao_value[0][0]), + (int64_t) elec_num*ao_num); +assert (rc == QMCKL_SUCCESS); + +printf("\n"); +printf(" ao_value ao_value[26][219] %25.15e\n", ao_value[26][219]); +printf(" ao_value ao_value[26][220] %25.15e\n", ao_value[26][220]); +printf(" ao_value ao_value[26][221] %25.15e\n", ao_value[26][221]); +printf(" ao_value ao_value[26][222] %25.15e\n", ao_value[26][222]); +printf(" ao_value ao_value[26][223] %25.15e\n", ao_value[26][223]); +printf(" ao_value ao_value[26][224] %25.15e\n", ao_value[26][224]); +printf("\n"); + +assert( fabs(ao_value[26][219] - ( 1.020298798341620e-08)) < 1.e-14 ); +assert( fabs(ao_value[26][220] - ( 1.516643537739178e-08)) < 1.e-14 ); +assert( fabs(ao_value[26][221] - ( -4.686370882518819e-09)) < 1.e-14 ); +assert( fabs(ao_value[26][222] - ( 7.514816980753531e-09)) < 1.e-14 ); +assert( fabs(ao_value[26][223] - ( -4.021908374204471e-09)) < 1.e-14 ); +assert( fabs(ao_value[26][224] - ( 7.175045873560788e-10)) < 1.e-14 ); + +} + + #+end_src + +** Value, gradients, Laplacian :PROPERTIES: :Name: qmckl_compute_ao_vgl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: +*** Unoptimized version + #+NAME: qmckl_ao_vgl_args_doc + | Variable | Type | In/Out | Description | + |-----------------------+-----------------------------------+--------+----------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~shell_num~ | ~int64_t~ | in | Number of shells | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~coord~ | ~double[3][point_num]~ | in | Coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | + | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | + | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | + | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | + | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | + | ~shell_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells | + | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs | -** Unoptimized version - #+NAME: qmckl_ao_vgl_args_doc - | Variable | Type | In/Out | Description | - |-----------------------+-----------------------------------+--------+----------------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~ao_num~ | ~int64_t~ | in | Number of AOs | - | ~shell_num~ | ~int64_t~ | in | Number of shells | - | ~point_num~ | ~int64_t~ | in | Number of points | - | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | - | ~coord~ | ~double[3][point_num]~ | in | Coordinates | - | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | - | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | - | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | - | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | - | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | - | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | - | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | - | ~shell_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells | - | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs | - - #+begin_src f90 :comments org :tangle (eval f) :noweb yes + #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_ao_vgl_doc_f(context, & ao_num, shell_num, point_num, nucl_num, & coord, nucl_coord, nucleus_index, nucleus_shell_num, & @@ -5443,34 +6321,34 @@ integer function qmckl_compute_ao_vgl_doc_f(context, & deallocate(poly_vgl, powers) end function qmckl_compute_ao_vgl_doc_f - #+end_src + #+end_src -** HPC version - #+NAME: qmckl_ao_vgl_args_hpc_gaussian - | Variable | Type | In/Out | Description | - |-----------------------+--------------------------------+--------+----------------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~ao_num~ | ~int64_t~ | in | Number of AOs | - | ~shell_num~ | ~int64_t~ | in | Number of shells | - | ~prim_num~ | ~int64_t~ | in | Number of primitives | - | ~point_num~ | ~int64_t~ | in | Number of points | - | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | - | ~coord~ | ~double[3][point_num]~ | in | Coordinates | - | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | - | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | - | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | - | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | - | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | - | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | - | ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell | - | ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell | - | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | - | ~ao_expo~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | - | ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | - | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs | +*** HPC version + #+NAME: qmckl_ao_vgl_args_hpc_gaussian + | Variable | Type | In/Out | Description | + |-----------------------+--------------------------------+--------+----------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~shell_num~ | ~int64_t~ | in | Number of shells | + | ~prim_num~ | ~int64_t~ | in | Number of primitives | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~coord~ | ~double[3][point_num]~ | in | Coordinates | + | ~nucl_coord~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~nucleus_index~ | ~int64_t[nucl_num]~ | in | Index of the 1st shell of each nucleus | + | ~nucleus_shell_num~ | ~int64_t[nucl_num]~ | in | Number of shells per nucleus | + | ~nucleus_range~ | ~double[nucl_num]~ | in | Range beyond which all is zero | + | ~nucleus_max_ang_mom~ | ~int32_t[nucl_num]~ | in | Maximum angular momentum per nucleus | + | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | + | ~shell_prim_index~ | ~int64_t[shell_num]~ | in | Index of the 1st primitive of each shell | + | ~shell_prim_num~ | ~int64_t[shell_num]~ | in | Number of primitives per shell | + | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | + | ~ao_expo~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | + | ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | + | ~ao_vgl~ | ~double[point_num][5][ao_num]~ | out | Value, gradients and Laplacian of the AOs | - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none #ifdef HAVE_HPC qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian ( @@ -5784,14 +6662,13 @@ qmckl_compute_ao_vgl_hpc_gaussian ( return QMCKL_SUCCESS; } #endif - #+end_src + #+end_src -** Interfaces -# #+CALL: generate_c_header(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl")) -# (Commented because the header needs to go into h_private_func) +*** Interfaces + # #+CALL: generate_c_header(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl")) + # (Commented because the header needs to go into h_private_func) - #+RESULTS: - #+begin_src c :tangle (eval h_private_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_ao_vgl_doc ( const qmckl_context context, const int64_t ao_num, @@ -5808,8 +6685,9 @@ qmckl_compute_ao_vgl_hpc_gaussian ( const double* ao_factor, const double* shell_vgl, double* const ao_vgl ); - #+end_src - #+begin_src c :tangle (eval h_private_func) :comments org + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org #ifdef HAVE_HPC qmckl_exit_code qmckl_compute_ao_vgl_hpc_gaussian ( const qmckl_context context, @@ -5830,12 +6708,12 @@ qmckl_compute_ao_vgl_hpc_gaussian ( const qmckl_tensor coef_per_nucleus, double* const ao_vgl ); #endif - #+end_src + #+end_src - #+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc")) + #+CALL: generate_c_interface(table=qmckl_ao_vgl_args_doc,rettyp=get_value("CRetType"),fname="qmckl_compute_ao_vgl_doc")) - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_ao_vgl_doc & (context, & ao_num, & @@ -5892,16 +6770,15 @@ qmckl_compute_ao_vgl_hpc_gaussian ( ao_vgl) end function qmckl_compute_ao_vgl_doc - #+end_src + #+end_src +**** Provide :noexport: -*** Provide :noexport: - - #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context); - #+end_src + #+end_src - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) { @@ -6034,11 +6911,11 @@ qmckl_exit_code qmckl_provide_ao_vgl(qmckl_context context) return QMCKL_SUCCESS; } - #+end_src + #+end_src -*** Test :noexport: +**** Test :noexport: - #+begin_src python :results output :exports none + #+begin_src python :results output :exports none import numpy as np from math import sqrt @@ -6107,11 +6984,11 @@ print ( "[26][1][223] : %25.15e"%(norm*df(a,x,y,1)* (x[1] - y[1]) * (x[2] - y[2] print ( "[26][0][224] : %25.15e"%(f(a,x,y) * (x[2] - y[2]) * (x[2] - y[2])) ) print ( "[26][1][224] : %25.15e"%(df(a,x,y,1)* (x[2] - y[2]) * (x[2] - y[2])) ) - #+end_src + #+end_src - #+RESULTS: + #+RESULTS: - #+begin_src c :tangle (eval c_test) :exports none + #+begin_src c :tangle (eval c_test) :exports none { #define walk_num 1 // chbrclf_walk_num #define elec_num chbrclf_elec_num @@ -6206,7 +7083,7 @@ assert( fabs(ao_vgl[26][4][224] - ( 3.153244195820293e-08)) < 1.e-14 ); } - #+end_src + #+end_src * End of files :noexport: diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index ac7fdb9..4a148c6 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -92,10 +92,12 @@ int main() { Computed data: - |---------------+--------------------------+-------------------------------------------------------------------------------------| - | ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions | - | ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at point positions | - |---------------+--------------------------+-------------------------------------------------------------------------------------| + |-----------------+--------------------------+-------------------------------------------------------------------------------------| + | ~mo_value~ | ~[point_num][mo_num]~ | Value of the MOs at point positions | + | ~mo_value_date~ | ~uint64_t~ | Late modification date of the value of the MOs at point positions | + | ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions | + | ~mo_vgl_date~ | ~uint64_t~ | Late modification date of Value, gradients, Laplacian of the MOs at point positions | + |-----------------+--------------------------+-------------------------------------------------------------------------------------| ** Data structure @@ -106,7 +108,9 @@ typedef struct qmckl_mo_basis_struct { double * restrict coefficient_t; double * restrict mo_vgl; + double * restrict mo_value; uint64_t mo_vgl_date; + uint64_t mo_value_date; int32_t uninitialized; bool provided; @@ -418,7 +422,464 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { * Computation -** Computation of MOs +** Computation of MOs: values only + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_mo_basis_mo_value(qmckl_context context, + double* const mo_value, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_mo_basis_mo_value(qmckl_context context, + double* const mo_value, + const int64_t size_max) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ao_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_mo_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + const int64_t sze = ctx->point.num * ctx->mo_basis.mo_num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_mo_basis_mo_value", + "input array too small"); + } + memcpy(mo_value, ctx->mo_basis.mo_value, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_get_mo_basis_mo_value (context, & + mo_value, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: mo_value(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function qmckl_get_mo_basis_mo_value + end interface + #+end_src + + Uses the given array to compute the values. + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_mo_basis_mo_value_inplace (qmckl_context context, + double* const mo_value, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_mo_basis_mo_value_inplace (qmckl_context context, + double* const mo_value, + const int64_t size_max) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_mo_basis_mo_value", + NULL); + } + + qmckl_exit_code rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + const int64_t sze = ctx->mo_basis.mo_num * ctx->point.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_mo_basis_mo_value", + "input array too small"); + } + + rc = qmckl_context_touch(context); + if (rc != QMCKL_SUCCESS) return rc; + + double* old_array = ctx->mo_basis.mo_value; + + ctx->mo_basis.mo_value = mo_value; + + rc = qmckl_provide_mo_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + ctx->mo_basis.mo_value = old_array; + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_get_mo_basis_mo_value_inplace (context, & + mo_value, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: mo_value(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function qmckl_get_mo_basis_mo_value_inplace + end interface + #+end_src + +*** Provide + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_mo_value(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_mo_value(qmckl_context context) +{ + + qmckl_exit_code rc; + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + if (!ctx->ao_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_basis", + NULL); + } + + rc = qmckl_provide_ao_value(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_value", + NULL); + } + + if (!ctx->mo_basis.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_mo_basis", + NULL); + } + + /* Compute if necessary */ + if (ctx->point.date > ctx->mo_basis.mo_value_date) { + + /* Allocate array */ + if (ctx->mo_basis.mo_value == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->point.num * ctx->mo_basis.mo_num * sizeof(double); + double* mo_value = (double*) qmckl_malloc(context, mem_info); + + if (mo_value == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_mo_basis_mo_value", + NULL); + } + ctx->mo_basis.mo_value = mo_value; + } + + if (ctx->mo_basis.mo_vgl_date == ctx->point.date) { + + // mo_vgl has been computed at this step: Just copy the data. + + double * v = &(ctx->mo_basis.mo_value[0]); + double * vgl = &(ctx->mo_basis.mo_vgl[0]); + for (int i=0 ; ipoint.num ; ++i) { + for (int k=0 ; kmo_basis.mo_num ; ++k) { + v[k] = vgl[k]; + } + v += ctx->mo_basis.mo_num; + vgl += ctx->mo_basis.mo_num * 5; + } + + } else { + + rc = qmckl_compute_mo_basis_mo_value(context, + ctx->ao_basis.ao_num, + ctx->mo_basis.mo_num, + ctx->point.num, + ctx->mo_basis.coefficient_t, + ctx->ao_basis.ao_value, + ctx->mo_basis.mo_value); + + if (rc != QMCKL_SUCCESS) { + return rc; + } + + } + + ctx->mo_basis.mo_value_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_mo_basis_mo_value + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_mo_basis_mo_value_args + | Variable | Type | In/Out | Description | + |---------------------+-----------------------------+--------+-------------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~ao_num~ | ~int64_t~ | in | Number of AOs | + | ~mo_num~ | ~int64_t~ | in | Number of MOs | + | ~point_num~ | ~int64_t~ | in | Number of points | + | ~coef_normalized_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix | + | ~ao_value~ | ~double[point_num][ao_num]~ | in | Value of the AOs | + | ~mo_value~ | ~double[point_num][mo_num]~ | out | Value of the MOs | + + + The matrix of AO values is very sparse, so we use a sparse-dense + matrix multiplication instead of a dgemm, as exposed in + https://dx.doi.org/10.1007/978-3-642-38718-0_14. + + + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_mo_basis_mo_value_doc_f(context, & + ao_num, mo_num, point_num, & + coef_normalized_t, ao_value, mo_value) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: ao_num, mo_num + integer*8 , intent(in) :: point_num + double precision , intent(in) :: ao_value(ao_num,point_num) + double precision , intent(in) :: coef_normalized_t(mo_num,ao_num) + double precision , intent(out) :: mo_value(mo_num,point_num) + integer*8 :: i,j,k + double precision :: c1, c2, c3, c4, c5 + + integer*8 :: LDA, LDB, LDC + + info = QMCKL_SUCCESS + if (.True.) then ! fast algorithm + do j=1,point_num + mo_value(:,j) = 0.d0 + do k=1,ao_num + if (ao_value(k,j) /= 0.d0) then + c1 = ao_value(k,j) + do i=1,mo_num + mo_value(i,j) = mo_value(i,j) + coef_normalized_t(i,k) * c1 + end do + end if + end do + end do + + else ! dgemm + + LDA = size(coef_normalized_t,1) + LDB = size(ao_value,1) + LDC = size(mo_value,1) + + info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, & + coef_normalized_t, LDA, ao_value, LDB, & + 0.d0, mo_value, LDC) + + end if + +end function qmckl_compute_mo_basis_mo_value_doc_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_mo_basis_mo_value ( + const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* coef_normalized_t, + const double* ao_value, + double* const mo_value ); + #+end_src + + #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc ( + const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* coef_normalized_t, + const double* ao_value, + double* const mo_value ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_mo_basis_mo_value_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_doc")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_mo_basis_mo_value_doc & + (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: ao_num + integer (c_int64_t) , intent(in) , value :: mo_num + integer (c_int64_t) , intent(in) , value :: point_num + real (c_double ) , intent(in) :: coef_normalized_t(ao_num,mo_num) + real (c_double ) , intent(in) :: ao_value(ao_num,point_num) + real (c_double ) , intent(out) :: mo_value(mo_num,point_num) + + integer(c_int32_t), external :: qmckl_compute_mo_basis_mo_value_doc_f + info = qmckl_compute_mo_basis_mo_value_doc_f & + (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value) + + end function qmckl_compute_mo_basis_mo_value_doc + #+end_src + + #+begin_src c :tangle (eval c) :comments org +qmckl_exit_code +qmckl_compute_mo_basis_mo_value (const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* coef_normalized_t, + const double* ao_value, + double* const mo_value ) +{ +#ifdef HAVE_HPC + return qmckl_compute_mo_basis_mo_value_hpc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value); +#else + return qmckl_compute_mo_basis_mo_value_doc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_value, mo_value); +#endif +} + #+end_src + +*** HPC version + + + #+begin_src c :tangle (eval h_func) :comments org +#ifdef HAVE_HPC +qmckl_exit_code +qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* coef_normalized_t, + const double* ao_value, + double* const mo_value ); +#endif + #+end_src + + #+begin_src c :tangle (eval c) :comments org +#ifdef HAVE_HPC +qmckl_exit_code +qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context, + const int64_t ao_num, + const int64_t mo_num, + const int64_t point_num, + const double* restrict coef_normalized_t, + const double* restrict ao_value, + double* restrict const mo_value ) +{ + assert (context != QMCKL_NULL_CONTEXT); + +#ifdef HAVE_OPENMP + #pragma omp parallel for +#endif + for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) { + double* restrict const vgl1 = &(mo_value[ipoint*mo_num]); + const double* restrict avgl1 = &(ao_value[ipoint*ao_num]); + + for (int64_t i=0 ; i