From 2ecfc55dbce7f60af365683f7debcb3e02d20954 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 17 Mar 2023 14:54:58 +0100 Subject: [PATCH] Added eN cusp fitting in MOs --- org/qmckl_mo.org | 360 +++++++++++++++++++++++++++++++++-------------- 1 file changed, 253 insertions(+), 107 deletions(-) diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index bf5b1ab..f8b6b20 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -71,6 +71,8 @@ cusp is exact. #+begin_src c :tangle (eval h_private_func) #ifndef QMCKL_MO_HPF #define QMCKL_MO_HPF + +#include "qmckl_blas_private_type.h" #+end_src #+begin_src c :tangle (eval h_private_type) @@ -78,6 +80,7 @@ cusp is exact. #define QMCKL_MO_HPT #include +#include "qmckl_blas_private_type.h" #+end_src #+begin_src c :tangle (eval c_test) :noweb yes @@ -140,11 +143,11 @@ int main() { Computed data: - |-----------------+--------------------------+----------------------------------------------------------------------------------| - | ~r_cusp_param~ | ~[nucl_num][mo_num][4]~ | Parameters of the functions for Cusp adjustments | - | ~mo_value~ | ~[point_num][mo_num]~ | Value of the MOs at point positions | - | ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions | - |-----------------+--------------------------+----------------------------------------------------------------------------------| + |--------------+--------------------------+-----------------------------------------------------------| + | ~cusp_param~ | ~[nucl_num][4][mo_num]~ | Parameters of the functions for Cusp adjustments | + | ~mo_value~ | ~[point_num][mo_num]~ | Value of the MOs at point positions | + | ~mo_vgl~ | ~[point_num][5][mo_num]~ | Value, gradients, Laplacian of the MOs at point positions | + |--------------+--------------------------+-----------------------------------------------------------| ** Data structure @@ -157,7 +160,7 @@ typedef struct qmckl_mo_basis_struct { double * restrict mo_vgl; double * restrict mo_value; - qmckl_tensor r_cusp_param; + qmckl_tensor cusp_param; uint64_t mo_vgl_date; uint64_t mo_value_date; @@ -431,36 +434,6 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context, double* coord = (double*) malloc(ctx->nucleus.num * 3 * sizeof(double)); - // Evaluate MO vgl at r_cusp - qmckl_tensor mo_vgl_at_r_cusp; - { - int64_t sze[3] = { ctx->mo_basis.mo_num, 5, ctx->nucleus.num }; - mo_vgl_at_r_cusp = qmckl_tensor_alloc(context, 3, &(sze[0])); - - rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3); - if (rc != QMCKL_SUCCESS) return rc; - - int64_t i=0; - for (int64_t inucl=0 ; inuclnucleus.num ; ++inucl) { - coord[i] += r_cusp[inucl]; - i+=3; - } - - rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3)); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_mo_basis_r_cusp", - "Unable to set coordinates at r_cusp"); - } - - rc = qmckl_get_mo_basis_mo_vgl(context, - &(qmckl_ten3(mo_vgl_at_r_cusp,0,0,0)), - ctx->mo_basis.mo_num * 5 * ctx->nucleus.num); - if (rc != QMCKL_SUCCESS) return rc; - - } - // Set r_cusp { assert (ctx->mo_basis.r_cusp == NULL); @@ -480,14 +453,14 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context, // Allocate cusp parameters and set them to zero { - if (ctx->mo_basis.r_cusp_param.size[0] != 0) { - rc = qmckl_tensor_free(context, &(ctx->mo_basis.r_cusp_param)); + if (ctx->mo_basis.cusp_param.size[0] != 0) { + rc = qmckl_tensor_free(context, &(ctx->mo_basis.cusp_param)); if (rc != QMCKL_SUCCESS) return rc; } - int64_t sze[3] = { 4, ctx->mo_basis.mo_num, ctx->nucleus.num }; - ctx->mo_basis.r_cusp_param = qmckl_tensor_alloc(context, 3, &(sze[0])); - ctx->mo_basis.r_cusp_param = qmckl_tensor_set(ctx->mo_basis.r_cusp_param, 0.); + int64_t sze[3] = { ctx->mo_basis.mo_num, 4, ctx->nucleus.num }; + ctx->mo_basis.cusp_param = qmckl_tensor_alloc(context, 3, &(sze[0])); + ctx->mo_basis.cusp_param = qmckl_tensor_set(ctx->mo_basis.cusp_param, 0.); } @@ -515,20 +488,24 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context, // Evaluate MO vgl at r_cusp without s components - qmckl_tensor mo_vgl_at_r_cusp_no_s; + qmckl_tensor mo_vgl_at_r_cusp_s; { - int64_t sze[3] = { ctx->mo_basis.mo_num, 5, ctx->nucleus.num }; - mo_vgl_at_r_cusp_no_s = qmckl_tensor_alloc(context, 3, &(sze[0])); + int64_t sze[3] = { ctx->mo_basis.mo_num, 3, ctx->nucleus.num }; + mo_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0])); + } + + { + qmckl_tensor ao_vgl_at_r_cusp_s; + int64_t sze[3] = { ctx->ao_basis.ao_num, 5, ctx->nucleus.num }; + ao_vgl_at_r_cusp_s = qmckl_tensor_alloc(context, 3, &(sze[0])); rc = qmckl_double_of_matrix(context, ctx->nucleus.coord, coord, ctx->nucleus.num * 3); if (rc != QMCKL_SUCCESS) return rc; - int64_t i=0; - for (int64_t inucl=0 ; inuclnucleus.num ; ++inucl) { - coord[i] += r_cusp[inucl]; - i+=3; + for (int64_t i=0 ; inucleus.num ; ++i) { + coord[2*ctx->nucleus.num + i] += r_cusp[i]; } - + rc = qmckl_set_point(context, 'T', ctx->nucleus.num, coord, (ctx->nucleus.num * 3)); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, @@ -537,11 +514,28 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context, "Unable to set coordinates at r_cusp"); } - rc = qmckl_get_mo_basis_mo_vgl(context, - &(qmckl_ten3(mo_vgl_at_r_cusp_no_s,0,0,0)), - ctx->mo_basis.mo_num * 5 * ctx->nucleus.num); - if (rc != QMCKL_SUCCESS) return rc; + rc = qmckl_get_ao_basis_ao_vgl(context, + &(qmckl_ten3(ao_vgl_at_r_cusp_s,0,0,0)), + ctx->ao_basis.ao_num * 5 * ctx->point.num); + for (int64_t inucl=0 ; inuclnucleus.num ; ++inucl) { + for (int64_t i=0 ; imo_basis.mo_num ; ++i) { + qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl) = 0.; + qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl) = 0.; + qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl) = 0.; + for (int64_t k=0 ; kao_basis.ao_num ; ++k) { + if ( ctx->ao_basis.ao_nucl[k] == inucl && ctx->ao_basis.ao_ang_mom[k] == 0) { + const double ck = ctx->mo_basis.coefficient[k + i*ctx->ao_basis.ao_num]; + qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl) += ck * qmckl_ten3(ao_vgl_at_r_cusp_s,k,0,inucl); + qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl) += ck * qmckl_ten3(ao_vgl_at_r_cusp_s,k,3,inucl); + qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl) += ck * qmckl_ten3(ao_vgl_at_r_cusp_s,k,4,inucl); + } + } + } + } + rc = qmckl_tensor_free(context,&ao_vgl_at_r_cusp_s); + + if (rc != QMCKL_SUCCESS) return rc; } // Compute parameters @@ -551,36 +545,30 @@ qmckl_set_mo_basis_r_cusp(qmckl_context context, if (Z < 0.1) continue; // Avoid dummy atoms const double R = r_cusp[inucl]; for (int64_t i=0 ; imo_basis.mo_num ; ++i) { - const double phi = qmckl_ten3(mo_vgl_at_r_cusp,i,0,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,0,inucl); - const double grad_phi = qmckl_ten3(mo_vgl_at_r_cusp,i,1,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,1,inucl); - const double lap_phi = qmckl_ten3(mo_vgl_at_r_cusp,i,4,inucl) - qmckl_ten3(mo_vgl_at_r_cusp_no_s,i,4,inucl); + const double phi = qmckl_ten3(mo_vgl_at_r_cusp_s,i,0,inucl); + const double grad_phi = qmckl_ten3(mo_vgl_at_r_cusp_s,i,1,inucl); + const double lap_phi = qmckl_ten3(mo_vgl_at_r_cusp_s,i,2,inucl); const double eta = qmckl_mat(mo_value_at_nucl_no_s,i,inucl); - qmckl_ten3(ctx->mo_basis.r_cusp_param,0,i,inucl) = + qmckl_ten3(ctx->mo_basis.cusp_param,i,0,inucl) = -(R*(2.*eta*Z-6.*grad_phi)+lap_phi*R*R+6.*phi)/(2.*R*Z-6.); - qmckl_ten3(ctx->mo_basis.r_cusp_param,1,i,inucl) = + qmckl_ten3(ctx->mo_basis.cusp_param,i,1,inucl) = (lap_phi*R*R*Z-6.*grad_phi*R*Z+6.*phi*Z+6.*eta*Z)/(2.*R*Z-6.); - qmckl_ten3(ctx->mo_basis.r_cusp_param,2,i,inucl) = + qmckl_ten3(ctx->mo_basis.cusp_param,i,2,inucl) = -(R*(-5.*grad_phi*Z-1.5*lap_phi)+lap_phi*R*R*Z+3.*phi*Z+3.*eta*Z+6.*grad_phi)/(R*R*Z-3.*R); - qmckl_ten3(ctx->mo_basis.r_cusp_param,3,i,inucl) = + qmckl_ten3(ctx->mo_basis.cusp_param,i,3,inucl) = (R*(-2.*grad_phi*Z-lap_phi)+0.5*lap_phi*R*R*Z+phi*Z+eta*Z+3.*grad_phi)/(R*R*R*Z-3.*R*R); -printf("%ld %ld %f %f %f %f\n",i, inucl, - qmckl_ten3(ctx->mo_basis.r_cusp_param,0,i,inucl), - qmckl_ten3(ctx->mo_basis.r_cusp_param,1,i,inucl), - qmckl_ten3(ctx->mo_basis.r_cusp_param,2,i,inucl), - qmckl_ten3(ctx->mo_basis.r_cusp_param,3,i,inucl)); } } } free(coord); qmckl_matrix_free(context, &mo_value_at_nucl_no_s); - qmckl_tensor_free(context, &mo_vgl_at_r_cusp_no_s); - qmckl_tensor_free(context, &mo_vgl_at_r_cusp); + qmckl_tensor_free(context, &mo_vgl_at_r_cusp_s); // Restore old points if (old_point_num > 0) { @@ -1115,6 +1103,7 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context) ctx->ao_basis.ao_ang_mom, ctx->electron.en_distance, ctx->mo_basis.r_cusp, + ctx->mo_basis.cusp_param, ctx->mo_basis.coefficient_t, ctx->ao_basis.ao_value, ctx->mo_basis.mo_value); @@ -1195,7 +1184,7 @@ end function qmckl_compute_mo_basis_mo_value_doc_f #+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 + #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_mo_basis_mo_value ( const qmckl_context context, const int64_t ao_num, @@ -1209,7 +1198,7 @@ end function qmckl_compute_mo_basis_mo_value_doc_f #+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 + #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc ( const qmckl_context context, const int64_t ao_num, @@ -1267,7 +1256,7 @@ qmckl_compute_mo_basis_mo_value (const qmckl_context context, *** HPC version - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org #ifdef HAVE_HPC qmckl_exit_code qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context, @@ -1582,7 +1571,10 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context) ctx->ao_basis.ao_nucl, ctx->ao_basis.ao_ang_mom, ctx->electron.en_distance, + ctx->nucleus.coord, + ctx->point.coord, ctx->mo_basis.r_cusp, + ctx->mo_basis.cusp_param, ctx->mo_basis.coefficient_t, ctx->ao_basis.ao_vgl, ctx->mo_basis.mo_vgl); @@ -1678,7 +1670,7 @@ end function qmckl_compute_mo_basis_mo_vgl_doc_f #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl")) #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_mo_basis_mo_vgl ( const qmckl_context context, const int64_t ao_num, @@ -1692,7 +1684,7 @@ end function qmckl_compute_mo_basis_mo_vgl_doc_f #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_doc")) #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_doc ( const qmckl_context context, const int64_t ao_num, @@ -1750,7 +1742,7 @@ qmckl_compute_mo_basis_mo_vgl (const qmckl_context context, *** HPC version - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org #ifdef HAVE_HPC qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context, @@ -1910,6 +1902,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context, | ~ao_ang_mom~ | ~int32_t[ao_num]~ | in | Angular momentum of the shell | | ~en_distance~ | ~double[point_num][nucl_num]~ | in | Electron-nucleus distances | | ~r_cusp~ | ~double[nucl_num]~ | in | Cusp-adjustment radius | + | ~cusp_param~ | ~double[nucl_num][4][mo_num]~ | in | Cusp-adjustment parameters | | ~coefficient_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 | Cusp correction for the values of the MOs | @@ -1920,7 +1913,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context, #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_mo_basis_mo_value_cusp_doc_f(context, & nucl_num, ao_num, mo_num, point_num, ao_nucl, ao_ang_mom, en_distance, & - r_cusp, coefficient_t, ao_value, mo_value) & + r_cusp, cusp_param, coefficient_t, ao_value, mo_value) & result(info) use qmckl implicit none @@ -1930,31 +1923,44 @@ integer function qmckl_compute_mo_basis_mo_value_cusp_doc_f(context, & integer*4 , intent(in) :: ao_ang_mom(ao_num) double precision , intent(in) :: en_distance(nucl_num, point_num) double precision , intent(in) :: r_cusp(nucl_num) - double precision , intent(in) :: coefficient_t(mo_num,ao_num) - double precision , intent(in) :: ao_value(ao_num,point_num) - double precision , intent(out) :: mo_value(mo_num,point_num) + double precision , intent(in) :: cusp_param(mo_num, 4, nucl_num) + double precision , intent(in) :: coefficient_t(mo_num, ao_num) + double precision , intent(in) :: ao_value(ao_num, point_num) + double precision , intent(out) :: mo_value(mo_num, point_num) - integer*8 :: i,j,k, inucl + integer*8 :: i, j, k, inucl + double precision :: r info = QMCKL_SUCCESS - do j=1,point_num - mo_value(:,j) = 0.d0 + do i=1,point_num + mo_value(:,i) = 0.d0 do k=1,ao_num - if (ao_value(k,j) == 0.d0) cycle + if (ao_value(k,i) == 0.d0) cycle inucl = ao_nucl(k)+1 - if ( (en_distance(inucl,j) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle - mo_value(:,j) = mo_value(:,j) + coefficient_t(:,k) * ao_value(k,j) - end do - end do + if ( (en_distance(inucl,i) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle + mo_value(:,i) = mo_value(:,i) + coefficient_t(:,k) * ao_value(k,i) + end do ! k + do inucl=1,nucl_num + r = en_distance(inucl,i) + if (r > r_cusp(inucl)) cycle + + do j=1,mo_num + mo_value(j,i) = mo_value(j,i) + & + cusp_param(j,1,inucl) + r*(cusp_param(j,2,inucl) + r*( & + cusp_param(j,3,inucl) + r* cusp_param(j,4,inucl) )) + enddo + enddo ! inucl + enddo ! i + end function qmckl_compute_mo_basis_mo_value_cusp_doc_f #+end_src - #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp")) + #+CALL: generate_c_header(table=qmckl_mo_basis_mo_value_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_value_cusp")) #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_mo_basis_mo_value_cusp ( const qmckl_context context, const int64_t nucl_num, @@ -1965,6 +1971,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f const int32_t* ao_ang_mom, const double* en_distance, const double* r_cusp, + const qmckl_tensor cusp_param, const double* coefficient_t, const double* ao_value, double* const mo_value ); @@ -1984,6 +1991,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f const int32_t* ao_ang_mom, const double* en_distance, const double* r_cusp, + const double* cusp_param, const double* coefficient_t, const double* ao_value, double* const mo_value ); @@ -2003,6 +2011,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f ao_ang_mom, & en_distance, & r_cusp, & + cusp_param, & coefficient_t, & ao_value, & mo_value) & @@ -2020,6 +2029,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f integer (c_int32_t) , intent(in) :: ao_ang_mom(ao_num) real (c_double ) , intent(in) :: en_distance(nucl_num,point_num) real (c_double ) , intent(in) :: r_cusp(nucl_num) + real (c_double ) , intent(in) :: cusp_param(mo_num,4,nucl_num) real (c_double ) , intent(in) :: coefficient_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) @@ -2035,6 +2045,7 @@ end function qmckl_compute_mo_basis_mo_value_cusp_doc_f ao_ang_mom, & en_distance, & r_cusp, & + cusp_param, & coefficient_t, & ao_value, & mo_value) @@ -2053,26 +2064,34 @@ qmckl_compute_mo_basis_mo_value_cusp (const qmckl_context context, const int32_t* ao_ang_mom, const double* en_distance, const double* r_cusp, + const qmckl_tensor cusp_param_tensor, const double* coefficient_t, const double* ao_value, double* const mo_value ) { + qmckl_exit_code rc; + #ifdef HAVE_HPC - return qmckl_compute_mo_basis_mo_value_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num, - ao_nucl, ao_ang_mom, en_distance, r_cusp, - coefficient_t, ao_value, mo_value ); + rc = qmckl_compute_mo_basis_mo_value_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num, + ao_nucl, ao_ang_mom, en_distance, r_cusp, + cusp_param_tensor, coefficient_t, ao_value, mo_value ); #else - return qmckl_compute_mo_basis_mo_value_cusp_doc (context, nucl_num, ao_num, mo_num, point_num, - ao_nucl, ao_ang_mom, en_distance, r_cusp, - coefficient_t, ao_value, mo_value ); + double * cusp_param = qmckl_alloc_double_of_tensor(context, cusp_param_tensor); + + rc = qmckl_compute_mo_basis_mo_value_cusp_doc (context, nucl_num, ao_num, mo_num, point_num, + ao_nucl, ao_ang_mom, en_distance, r_cusp, + cusp_param, coefficient_t, ao_value, mo_value ); + + qmckl_free(context, cusp_param); #endif + return rc; } #+end_src *** HPC version - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org #ifdef HAVE_HPC qmckl_exit_code qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context, @@ -2084,6 +2103,7 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context, const int32_t* ao_ang_mom, const double* en_distance, const double* r_cusp, + const qmckl_tensor cusp_param, const double* coefficient_t, const double* ao_value, double* const mo_value ); @@ -2102,6 +2122,7 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context, const int32_t* ao_ang_mom, const double* en_distance, const double* r_cusp, + const qmckl_tensor cusp_param, const double* coefficient_t, const double* ao_value, double* const mo_value) @@ -2158,14 +2179,30 @@ qmckl_compute_mo_basis_mo_value_cusp_hpc (const qmckl_context context, for (int64_t m=n ; m < nidx ; m+=1) { const double* restrict ck = coefficient_t + idx[m]*mo_num; const double a1 = av1[m]; - + #ifdef HAVE_OPENMP - #pragma omp simd +#pragma omp simd #endif for (int64_t i=0 ; i r_cusp(inucl)) cycle + + r_vec(1:3) = point_coord(j,1:3) - nucl_coord(inucl,1:3) + r_inv = 1.d0/r + + + do i=1,mo_num + mo_vgl(i,1,j) = mo_vgl(i,1,j) + & + cusp_param(i,1,inucl) + r*(cusp_param(i,2,inucl) + r*( & + cusp_param(i,3,inucl) + r* cusp_param(i,4,inucl) )) + + c1 = r_inv * cusp_param(i,2,inucl) + 2.d0*cusp_param(i,3,inucl) + & + r * 3.d0 * cusp_param(i,4,inucl) + + mo_vgl(i,2,j) = mo_vgl(i,2,j) + r_vec(1) * c1 + mo_vgl(i,3,j) = mo_vgl(i,3,j) + r_vec(2) * c1 + mo_vgl(i,4,j) = mo_vgl(i,4,j) + r_vec(3) * c1 + + mo_vgl(i,5,j) = mo_vgl(i,5,j) + & + 2.d0*cusp_param(i,2,inucl)*r_inv + & + 6.d0*cusp_param(i,3,inucl) + & + 12.d0*cusp_param(i,4,inucl)*r + + enddo + enddo ! inucl end do info = QMCKL_SUCCESS end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f #+end_src - #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp")) - +# #+CALL: generate_c_header(table=qmckl_mo_basis_mo_vgl_cusp_args,rettyp=get_value("CRetType"),fname="qmckl_compute_mo_basis_mo_vgl_cusp")) + #+RESULTS: - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp ( const qmckl_context context, const int64_t nucl_num, @@ -2256,7 +2331,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f const int64_t* ao_nucl, const int32_t* ao_ang_mom, const double* en_distance, + const qmckl_matrix nucl_coord, + const qmckl_matrix point_coord, const double* r_cusp, + const qmckl_tensor cusp_param, const double* coefficient_t, const double* ao_vgl, double* const mo_vgl ); @@ -2275,7 +2353,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f const int64_t* ao_nucl, const int32_t* ao_ang_mom, const double* en_distance, + const double* nucl_coord, + const double* point_coord, const double* r_cusp, + const double* cusp_param, const double* coefficient_t, const double* ao_vgl, double* const mo_vgl ); @@ -2294,7 +2375,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f ao_nucl, & ao_ang_mom, & en_distance, & + nucl_coord, & + point_coord, & r_cusp, & + cusp_param, & coefficient_t, & ao_vgl, & mo_vgl) & @@ -2311,7 +2395,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f integer (c_int64_t) , intent(in) :: ao_nucl(ao_num) integer (c_int32_t) , intent(in) :: ao_ang_mom(ao_num) real (c_double ) , intent(in) :: en_distance(nucl_num,point_num) + real (c_double ) , intent(in) :: nucl_coord(nucl_num,3) + real (c_double ) , intent(in) :: point_coord(point_num,3) real (c_double ) , intent(in) :: r_cusp(nucl_num) + real (c_double ) , intent(in) :: cusp_param(mo_num,4,nucl_num) real (c_double ) , intent(in) :: coefficient_t(ao_num,mo_num) real (c_double ) , intent(in) :: ao_vgl(ao_num,5,point_num) real (c_double ) , intent(out) :: mo_vgl(mo_num,5,point_num) @@ -2326,7 +2413,10 @@ end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f ao_nucl, & ao_ang_mom, & en_distance, & + nucl_coord, & + point_coord, & r_cusp, & + cusp_param, & coefficient_t, & ao_vgl, & mo_vgl) @@ -2344,27 +2434,43 @@ qmckl_compute_mo_basis_mo_vgl_cusp (const qmckl_context context, const int64_t* ao_nucl, const int32_t* ao_ang_mom, const double* en_distance, + const qmckl_matrix nucl_coord_matrix, + const qmckl_matrix point_coord_matrix, const double* r_cusp, + const qmckl_tensor cusp_param_tensor, const double* coefficient_t, const double* ao_vgl, double* const mo_vgl ) { + qmckl_exit_code rc; + #ifdef HAVE_HPC - return qmckl_compute_mo_basis_mo_vgl_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num, - ao_nucl, ao_ang_mom, en_distance, r_cusp, - coefficient_t, ao_vgl, mo_vgl ); + rc = qmckl_compute_mo_basis_mo_vgl_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num, + ao_nucl, ao_ang_mom, en_distance, nucl_coord_matrix, + point_coord_matrix, r_cusp, cusp_param_tensor, + coefficient_t, ao_vgl, mo_vgl ); #else - return qmckl_compute_mo_basis_mo_vgl_cusp_doc (context, nucl_num, ao_num, mo_num, point_num, - ao_nucl, ao_ang_mom, en_distance, r_cusp, - coefficient_t, ao_vgl, mo_vgl ); + double * nucl_coord = qmckl_alloc_double_of_matrix(context, nucl_coord_matrix); + double * point_coord = qmckl_alloc_double_of_matrix(context, point_coord_matrix); + double * cusp_param = qmckl_alloc_double_of_tensor(context, cusp_param_tensor); + + rc = qmckl_compute_mo_basis_mo_vgl_cusp_doc (context, nucl_num, ao_num, mo_num, point_num, + ao_nucl, ao_ang_mom, en_distance, nucl_coord, + point_coord, r_cusp, cusp_param, coefficient_t, + ao_vgl, mo_vgl ); + + qmckl_free(context, nucl_coord); + qmckl_free(context, point_coord); + qmckl_free(context, cusp_param); #endif + return rc; } #+end_src *** HPC version - #+begin_src c :tangle (eval h_func) :comments org + #+begin_src c :tangle (eval h_private_func) :comments org #ifdef HAVE_HPC qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context, @@ -2375,7 +2481,10 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context, const int64_t* ao_nucl, const int32_t* ao_ang_mom, const double* en_distance, + const qmckl_matrix nucl_coord, + const qmckl_matrix point_coord, const double* r_cusp, + const qmckl_tensor cusp_param, const double* coefficient_t, const double* ao_vgl, double* const mo_vgl ); @@ -2393,7 +2502,10 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context, const int64_t* ao_nucl, const int32_t* ao_ang_mom, const double* en_distance, + const qmckl_matrix nucl_coord, + const qmckl_matrix point_coord, const double* r_cusp, + const qmckl_tensor cusp_param, const double* coefficient_t, const double* ao_vgl, double* const mo_vgl ) @@ -2512,6 +2624,40 @@ qmckl_compute_mo_basis_mo_vgl_cusp_hpc (const qmckl_context context, vgl5[i] += ck[i] * a5; } } + + // TODO + for (int64_t inucl=0 ; inucl