#+TITLE: Forces #+SETUPFILE: ../tools/theme.setup #+INCLUDE: ../tools/lib.org * Introduction To calculate interatomic forces, a derivative (with respect to nuclei) of the values, gradient and Laplacian of the Jastrow and the molecular orbitals (MO) is required. Furthermore, for the nonlocal pseudopotential, also the forces of the values of the single-electron move Jastrows are required. * Headers :noexport: #+begin_src elisp :noexport :results none (org-babel-lob-ingest "../tools/lib.org") #+end_src #+begin_src c :tangle (eval h_private_func) #ifndef QMCKL_FORCES_HPF #define QMCKL_FORCES_HPF #+end_src #+begin_src c :tangle (eval h_private_type) #ifndef QMCKL_FORCES_HPT #define QMCKL_FORCES_HPT #include #+end_src #+begin_src c :tangle (eval c_test) :noweb yes #include "qmckl.h" #include #include #ifdef HAVE_CONFIG_H #include "config.h" #endif #include #include #include "n2.h" #include "chbrclf.h" #include "qmckl_context_private_type.h" #include "qmckl_jastrow_champ_private_func.h" #include "qmckl_jastrow_champ_single_private_func.h" #include "qmckl_forces_private_func.h" int main() { qmckl_context context; context = qmckl_context_create(); #+end_src #+begin_src c :tangle (eval c) #ifdef HAVE_CONFIG_H #include "config.h" #endif #ifdef HAVE_STDINT_H #include #elif HAVE_INTTYPES_H #include #endif #include #include #include #include #include #include #include "qmckl.h" #include "qmckl_context_private_type.h" #include "qmckl_memory_private_type.h" #include "qmckl_memory_private_func.h" #include "qmckl_jastrow_champ_private_type.h" #include "qmckl_jastrow_champ_private_func.h" #include "qmckl_jastrow_champ_single_private_type.h" #include "qmckl_jastrow_champ_single_private_func.h" #include "qmckl_forces_private_type.h" #include "qmckl_forces_private_func.h" #+end_src * Context ** Data structure #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_forces_struct{ double * restrict forces_jastrow_en; uint64_t forces_jastrow_en_date; double * restrict forces_jastrow_en_g; uint64_t forces_jastrow_en_g_date; double * restrict forces_jastrow_en_l; uint64_t forces_jastrow_en_l_date; double * restrict forces_tmp_c; uint64_t forces_tmp_c_date; double * restrict forces_dtmp_c; uint64_t forces_dtmp_c_date; double * restrict forces_een_n; uint64_t forces_een_n_date; double * restrict forces_jastrow_een; uint64_t forces_jastrow_een_date; double * restrict forces_jastrow_een_g; uint64_t forces_jastrow_een_g_date; double * restrict forces_jastrow_een_l; uint64_t forces_jastrow_een_l_date; double * restrict forces_ao_value; uint64_t forces_ao_value_date; uint64_t forces_ao_value_maxsize; double * restrict forces_mo_value; uint64_t forces_mo_value_date; uint64_t forces_mo_value_maxsize; double * restrict forces_mo_g; uint64_t forces_mo_g_date; double * restrict forces_mo_l; uint64_t forces_mo_l_date; double * forces_jastrow_single_en; uint64_t forces_jastrow_single_en_date; double * forces_jastrow_single_een; uint64_t forces_jastrow_single_een_date; double * forces_delta_p; uint64_t forces_delta_p_date; } qmckl_forces_struct; #+end_src ** Test :noexport: #+begin_src c :tangle (eval c_test) /* Reference input data */ int64_t walk_num = n2_walk_num; int64_t elec_num = n2_elec_num; int64_t elec_up_num = n2_elec_up_num; int64_t elec_dn_num = n2_elec_dn_num; int64_t nucl_num = n2_nucl_num; double rescale_factor_ee = 0.6; double rescale_factor_en[2] = { 0.6, 0.6 }; double* elec_coord = &(n2_elec_coord[0][0][0]); const double* nucl_charge = n2_charge; double* nucl_coord = &(n2_nucl_coord[0][0]); /* Provide Electron data */ qmckl_exit_code rc; assert(!qmckl_electron_provided(context)); rc = qmckl_check(context, qmckl_set_electron_num (context, elec_up_num, elec_dn_num) ); assert(rc == QMCKL_SUCCESS); assert(qmckl_electron_provided(context)); rc = qmckl_check(context, qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num) ); assert(rc == QMCKL_SUCCESS); double elec_coord2[walk_num*3*elec_num]; rc = qmckl_check(context, qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num) ); assert(rc == QMCKL_SUCCESS); for (int64_t i=0 ; i<3*elec_num ; ++i) { assert( elec_coord[i] == elec_coord2[i] ); } /* Provide Nucleus data */ assert(!qmckl_nucleus_provided(context)); rc = qmckl_check(context, qmckl_set_nucleus_num (context, nucl_num) ); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); double nucl_coord2[3*nucl_num]; rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*nucl_num); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_check(context, qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num) ); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); rc = qmckl_check(context, qmckl_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3) ); assert(rc == QMCKL_SUCCESS); for (int64_t k=0 ; k<3 ; ++k) { for (int64_t i=0 ; inucleus.num * ctx->electron.walker.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_en", "Array too small. Expected 3*nucl_num*walk_num"); } memcpy(forces_jastrow_en, ctx->forces.forces_jastrow_en, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_en (context, & forces_jastrow_en, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: forces_jastrow_en(size_max) end function qmckl_get_forces_jastrow_en end interface #+end_src **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_en(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_en(qmckl_context context) { qmckl_exit_code rc; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_jastrow_en", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->jastrow_champ.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_forces_jastrow_en", NULL); } /* Check if en rescaled distance is provided */ rc = qmckl_provide_en_distance_rescaled(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_en_distance_rescaled_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->date > ctx->forces.forces_jastrow_en_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_jastrow_en != NULL) { rc = qmckl_free(context, ctx->forces.forces_jastrow_en); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_en", "Unable to free ctx->forces.forces_jastrow_en"); } ctx->forces.forces_jastrow_en = NULL; } } /* Allocate array */ if (ctx->forces.forces_jastrow_en == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * 3 * ctx->nucleus.num * sizeof(double); double* forces_jastrow_en = (double*) qmckl_malloc(context, mem_info); if (forces_jastrow_en == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_en", NULL); } ctx->forces.forces_jastrow_en = forces_jastrow_en; } rc = qmckl_compute_forces_jastrow_en(context, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.aord_num, ctx->jastrow_champ.a_vector, ctx->jastrow_champ.en_distance_rescaled, ctx->jastrow_champ.en_distance_rescaled_gl, ctx->forces.forces_jastrow_en); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_jastrow_en_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_en :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_jastrow_en_args | Variable | Type | In/Out | Description | |---------------------------+-------------------------------------------+--------+---------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~type_nucl_num~ | ~int64_t~ | in | Number of unique nuclei | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | IDs of unique nuclei | | ~aord_num~ | ~int64_t~ | in | Number of coefficients | | ~a_vector~ | ~double[type_nucl_num][aord_num+1]~ | in | List of coefficients | | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus distances | | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in | Electron-nucleus distance derivatives | | ~forces_jastrow_en~ | ~double[walk_num][nucl_num][3]~ | out | Electron-nucleus forces | |---------------------------+-------------------------------------------+--------+---------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_forces_jastrow_en_doc( & context, walk_num, elec_num, nucl_num, type_nucl_num, & type_nucl_vector, aord_num, a_vector, & en_distance_rescaled, en_distance_rescaled_gl, forces_jastrow_en) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) integer (c_int64_t) , intent(in) , value :: aord_num real (c_double ) , intent(in) :: a_vector(aord_num+1,type_nucl_num) real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num) real (c_double ) , intent(in) :: en_distance_rescaled_gl(4, elec_num,nucl_num,walk_num) real (c_double ) , intent(out) :: forces_jastrow_en(3,nucl_num,walk_num) integer(qmckl_exit_code) :: info integer*8 :: i, a, k, nw, ii double precision :: x, x1, kf double precision :: denom, invdenom, invdenom2, f double precision :: dx(3) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (aord_num < 0) then info = QMCKL_INVALID_ARG_7 return endif do nw =1, walk_num forces_jastrow_en(:,:,nw) = 0.0d0 do a = 1, nucl_num do i = 1, elec_num x = en_distance_rescaled(i,a,nw) if(abs(x) < 1.d-12) continue denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x invdenom = 1.0d0 / denom invdenom2 = invdenom*invdenom dx(1) = -en_distance_rescaled_gl(1,i,a,nw) dx(2) = -en_distance_rescaled_gl(2,i,a,nw) dx(3) = -en_distance_rescaled_gl(3,i,a,nw) f = a_vector(1, type_nucl_vector(a)+1) * invdenom2 forces_jastrow_en(1,a,nw) = forces_jastrow_en(1,a,nw) + f * dx(1) forces_jastrow_en(2,a,nw) = forces_jastrow_en(2,a,nw) + f * dx(2) forces_jastrow_en(3,a,nw) = forces_jastrow_en(3,a,nw) + f * dx(3) kf = 2.d0 x1 = x x = 1.d0 do k=2, aord_num f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x forces_jastrow_en(1,a,nw) = forces_jastrow_en(1,a,nw) + f * x1 * dx(1) forces_jastrow_en(2,a,nw) = forces_jastrow_en(2,a,nw) + f * x1 * dx(2) forces_jastrow_en(3,a,nw) = forces_jastrow_en(3,a,nw) + f * x1 * dx(3) x = x*x1 kf = kf + 1.d0 end do end do end do end do end function qmckl_compute_forces_jastrow_en_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_en_doc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en ); qmckl_exit_code qmckl_compute_forces_jastrow_en ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en ); qmckl_exit_code qmckl_compute_forces_jastrow_en ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_forces_jastrow_en (const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en ) { #ifdef HAVE_HPC return qmckl_compute_forces_jastrow_en_doc #else return qmckl_compute_forces_jastrow_en_doc #endif (context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, aord_num, a_vector, en_distance_rescaled, en_distance_rescaled_gl, forces_jastrow_en); } #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces Jastrow en\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double forces_jastrow_en[walk_num][nucl_num][3]; rc = qmckl_get_forces_jastrow_en(context, &forces_jastrow_en[0][0][0], 3*nucl_num*walk_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_en[walk_num][nucl_num][3]; rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_en, &(finite_difference_force_en[0][0][0]), 1); for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ //printf("%.10f\t", finite_difference_force_en[nw][a][k]); //printf("%.10f\n", forces_jastrow_en[nw][a][k]); } } } for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ assert(fabs(finite_difference_force_en[nw][a][k] - forces_jastrow_en[nw][a][k]) < 1.e-8); } } } printf("OK\n"); #+end_src *** Force of electron-nucleus Jastrow gradient To avoid having to create new arrays for the Hessian of the rescaled distances, we compute parts of it manually in this routine. \begin{eqnarray} \partial_{\alpha,m} \partial_{i,n} J_{en} & = \frac{a_1 \partial_{\alpha,m} \partial_{i,n} \widetilde{R}_{i\alpha} }{(1+a_2 \widetilde{R}_{i\alpha})^2} - 2 \frac{a_1 a_2 \partial_{\alpha,m} \widetilde{R}_{i\alpha} \partial_{i,n} \widetilde{R}_{i\alpha}}{(1+a_2 \widetilde{R}_{i\alpha})^3} + \sum_{k=2}^{N_\text{aord}} a_k k \left((k-1) \widetilde{R}_{i\alpha}^{k-2} \partial_{\alpha,m} \widetilde{R}_{i\alpha} \partial_{i,n} \widetilde{R}_{i\alpha} + \widetilde{R}_{i\alpha}^{k-1} \partial_{\alpha,m}\partial_{i,n} \widetilde{R}_{i\alpha} \right) \\ & = -\frac{a_1 e^{-\kappa R_{i\alpha}}}{R_{i\alpha }(1+a_2 \widetilde{R}_{i\alpha})^2} + \frac{a_1 \partial_{i,m} \widetilde{R}_{i\alpha} \partial_{i,n} \widetilde{R}_{i\alpha} }{(1+a_2 \widetilde{R}_{i\alpha})^2}\left (\frac{e^{\kappa R_{i\alpha}}}{R_{i\alpha}} + \kappa e^{\kappa R_{i\alpha}} \right) + 2 \frac{a_1 a_2 \partial_{i,m} \widetilde{R}_{i\alpha} \partial_{i,n} \widetilde{R}_{i\alpha}}{(1+a_2 \widetilde{R}_{i\alpha})^3} \\ &- \sum_{k=2}^{N_\text{aord}} a_k k \left((k-1) \widetilde{R}_{i\alpha}^{k-2} \partial_{\alpha,m} \widetilde{R}_{i\alpha} \partial_{i,n} \widetilde{R}_{i\alpha} - \widetilde{R}_{i\alpha}^{k-1} \partial_{i,m}\widetilde{R}_{i\alpha}\partial_{i,n} \widetilde{R}_{i\alpha} e^{\kappa R_{i\alpha}} (\kappa + \frac{1}{R_{i\alpha}}) + \delta_{n,m} e^{-\kappa R_{i\alpha}}\frac{1}{R_{i\alpha}} \right)\\ \end{eqnarray} **** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_jastrow_en_g(qmckl_context context, double* const forces_jastrow_en_g, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_jastrow_en_g(qmckl_context context, double* const forces_jastrow_en_g, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_forces_jastrow_en_g(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 3*3 * ctx->nucleus.num * ctx->electron.walker.num * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_en_g", "Array too small. Expected 3*3*nucl_num*walk_num_elec_num"); } memcpy(forces_jastrow_en_g, ctx->forces.forces_jastrow_en_g, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_en_g (context, & forces_jastrow_en_g, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: forces_jastrow_en_g(size_max) end function qmckl_get_forces_jastrow_en_g end interface #+end_src **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_en_g(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_en_g(qmckl_context context) { qmckl_exit_code rc; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_jastrow_en_g", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->jastrow_champ.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_forces_jastrow_en_g", NULL); } /* Check if en rescaled distance is provided */ rc = qmckl_provide_en_distance_rescaled(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_en_distance_rescaled_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->date > ctx->forces.forces_jastrow_en_g_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_jastrow_en_g != NULL) { rc = qmckl_free(context, ctx->forces.forces_jastrow_en_g); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_en", "Unable to free ctx->forces.forces_jastrow_en_g"); } ctx->forces.forces_jastrow_en_g = NULL; } } /* Allocate array */ if (ctx->forces.forces_jastrow_en_g == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * 3 * 3 * ctx->nucleus.num * ctx->electron.num * sizeof(double); double* forces_jastrow_en_g = (double*) qmckl_malloc(context, mem_info); if (forces_jastrow_en_g == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_en_g", NULL); } ctx->forces.forces_jastrow_en_g = forces_jastrow_en_g; } rc = qmckl_compute_forces_jastrow_en_g(context, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.aord_num, ctx->jastrow_champ.a_vector, ctx->jastrow_champ.rescale_factor_en, ctx->electron.en_distance, ctx->jastrow_champ.en_distance_rescaled, ctx->jastrow_champ.en_distance_rescaled_gl, ctx->forces.forces_jastrow_en_g); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_jastrow_en_g_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_en_g :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_jastrow_en_g_args | Variable | Type | In/Out | Description | |---------------------------+----------------------------------------------+--------+---------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~type_nucl_num~ | ~int64_t~ | in | Number of unique nuclei | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | IDs of unique nuclei | | ~aord_num~ | ~int64_t~ | in | Number of coefficients | | ~a_vector~ | ~double[type_nucl_num][aord_num+1]~ | in | List of coefficients | | ~rescale_factor_en~ | ~double[type_nucl_num]~ | in | Rescale factor for electron-nucleus | | ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | in | Electron-nucleus distances | | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus distances | | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in | Electron-nucleus distance derivatives | | ~forces_jastrow_en_g~ | ~double[walk_num][nucl_num][3][elec_num][3]~ | out | Force of electron-nucleus gradient | |---------------------------+----------------------------------------------+--------+---------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_forces_jastrow_en_g_doc( & context, walk_num, elec_num, nucl_num, type_nucl_num, & type_nucl_vector, aord_num, a_vector, rescale_factor_en, en_distance, & en_distance_rescaled, en_distance_rescaled_gl, forces_jastrow_en_g) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) integer (c_int64_t) , intent(in) , value :: aord_num real (c_double ) , intent(in) :: a_vector(aord_num+1,type_nucl_num) real (c_double ) , intent(in) :: rescale_factor_en(type_nucl_num) real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num) real (c_double ) , intent(in) :: en_distance_rescaled_gl(4, elec_num,nucl_num,walk_num) real (c_double ) , intent(out) :: forces_jastrow_en_g(3,elec_num,3,nucl_num,walk_num) integer(qmckl_exit_code) :: info integer*8 :: i, a, k, nw, ii, m,l double precision :: x, x1, kf double precision :: denom, invdenom, invdenom2, f, f2, expk, invdist double precision :: dx(4) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (aord_num < 0) then info = QMCKL_INVALID_ARG_7 return endif do nw =1, walk_num forces_jastrow_en_g(:,:,:,:,nw) = 0.0d0 do a = 1, nucl_num do i = 1, elec_num expk = dexp(rescale_factor_en(type_nucl_vector(a)+1) * en_distance(a,i,nw)) invdist = 1.d0 / en_distance(a,i,nw) x = en_distance_rescaled(i,a,nw) if(abs(x) < 1.d-12) continue denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x invdenom = 1.0d0 / denom invdenom2 = invdenom*invdenom f = a_vector(1, type_nucl_vector(a)+1) * invdenom2 do m = 1, 3 dx(m) = en_distance_rescaled_gl(m,i,a,nw) end do do m = 1, 3 do l = 1,3 if (m == l) then forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) - f * invdist / expk end if forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) + & f * dx(m) * dx(l) * expk * (invdist + rescale_factor_en(type_nucl_vector(a)+1)) forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) + 2.d0 * f * invdenom * & a_vector(2, type_nucl_vector(a)+1) * dx(m) * dx(l) end do end do kf = 2.d0 x1 = x x = 1.d0 do k=2, aord_num f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x f2 = a_vector(k+1,type_nucl_vector(a)+1) * kf * x * (kf-1.d0) do m = 1, 3 do l = 1, 3 if (m == l) then forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) - f * x1 * invdist / expk end if forces_jastrow_en_g(m,i,l,a,nw) = forces_jastrow_en_g(m,i,l,a,nw) - f2 * dx(m) * dx(l) & + f * x1 * dx(m) * dx(l) * rescale_factor_en(type_nucl_vector(a)+1) * expk & + f * x1 * dx(m) * dx(l) * invdist * expk end do end do x = x*x1 kf = kf + 1.d0 end do end do end do end do end function qmckl_compute_forces_jastrow_en_g_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_en_g_doc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* rescale_factor_en, const double* en_distance, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en_g ); qmckl_exit_code qmckl_compute_forces_jastrow_en_g ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* rescale_factor_en, const double* en_distance, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en_g ); qmckl_exit_code qmckl_compute_forces_jastrow_en_g ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* rescale_factor_en, const double* en_distance, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en_g ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_forces_jastrow_en_g (const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* rescale_factor_en, const double* en_distance, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en_g ) { #ifdef HAVE_HPC return qmckl_compute_forces_jastrow_en_g_doc #else return qmckl_compute_forces_jastrow_en_g_doc #endif (context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, aord_num, a_vector, rescale_factor_en, en_distance, en_distance_rescaled, en_distance_rescaled_gl, forces_jastrow_en_g); } #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces Jastrow en G\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double forces_jastrow_en_g[walk_num][nucl_num][3][elec_num][3]; rc = qmckl_get_forces_jastrow_en_g(context, &forces_jastrow_en_g[0][0][0][0][0], 3*3*nucl_num*walk_num*elec_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_en_g[walk_num][nucl_num][3][4][elec_num]; rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_en_gl, &finite_difference_force_en_g[0][0][0][0][0], 4*elec_num); for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ for (int i = 0; i < elec_num; i++){ for (int l = 1; l < 3; l++){ //printf("finite_difference_force_en_g[%i][%i][%i][%i][%i] %+3.10f \n", nw,a,k,l,i,finite_difference_force_en_g[nw][a][k][l][i]); //printf("forces_jastrow_en_g [%i][%i][%i][%i][%i] %+3.10f\n", nw,a,k,i,l,forces_jastrow_en_g[nw][a][k][i][l]); } } } } } for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ for (int i = 0; i < elec_num; i++){ for (int l = 1; l < 3; l++){ assert(fabs(finite_difference_force_en_g[nw][a][k][l][i] - forces_jastrow_en_g[nw][a][k][i][l]) < 1.e-8); } } } } } printf("OK\n"); #+end_src *** Force of electron-nucleus Jastrow Laplacian Calculates the force of the electron-nucleus Jastrow Laplacian. To avoid having to calculate higher order derivatives of the rescaled distances elsewhere, the neccessery terms are calcuculated on the fly. \begin{eqnarray} \partial_{\alpha,m} \nabla^2 J_{en} =& \sum_i \left(\partial_{\alpha,m}\frac{a_1 \nabla^2 \widetilde{R}_{i\alpha}}{(1+a_2 \widetilde{R}_{i\alpha})^2} - \partial_{\alpha,m} \frac{2a_1a_2 (\nabla\widetilde{R}_{i\alpha}\cdot \nabla \widetilde{R}_{i\alpha})}{(1+a_2\widetilde{R}_{i\alpha})^3} + \partial_{\alpha,m} \sum_{k=2}^{N_\text{aord}} a_k k \left( \widetilde{R}_{i\alpha}^{k-1} \nabla^2\widetilde{R}_{i\alpha} + (k-1)\widetilde{R}_{i\alpha}^{k-2} (\nabla \widetilde{R}_{i\alpha}\cdot \nabla \widetilde{R}_{i\alpha} ) \right) \right)\\ =&\sum_i \left( \frac{a_1 \partial_{i,m} \widetilde{R}_{i\alpha}}{(1+a_2 \widetilde{R}_{i\alpha})^2}\left( \frac{2\kappa}{R_{i\alpha}} - \kappa^2 + \frac{2}{R^2_{i\alpha}}\right) + \frac{2a_1a_2\partial_{i,m}\widetilde{R}_{i\alpha}\nabla^2\widetilde{R}_{i\alpha}}{(1+a_2 \widetilde{R}_{i\alpha})^3} - \frac{4a_1a_2\kappa e^{-\kappa R_{i\alpha}}\partial_{i,m}\widetilde{R}_{i\alpha}}{(1+a_2 \widetilde{R}_{i\alpha})^3} - \frac{6 a_1 a_2^2 (\nabla \widetilde{R}_{i\alpha} \cdot \nabla \widetilde{R}_{i\alpha})\partial_{i,m}\widetilde{R}_{i\alpha}}{(1+a_2 \widetilde{R}_{i\alpha})^4} \right.\\ &\left.+ \sum_{k=2}^{N_\text{aord}} a_k k \left( -(k-1) \widetilde{R}_{i\alpha}^{k-2}\partial_{i,m}\widetilde{R}_{i\alpha} \nabla^2\widetilde{R}_{i\alpha} + \widetilde{R}_{i\alpha}^{k-1}\partial_{i,m}\widetilde{R}_{i\alpha}\left(\frac{2\kappa}{R_{i\alpha}} - \kappa^2 + \frac{2}{R_{i\alpha}^2} \right) - (k-1)(k-2)\widetilde{R}_{i\alpha}^{k-3}\partial_{i,m}\widetilde{R}_{i\alpha}(\nabla\widetilde{R}_{i\alpha} \cdot \nabla \widetilde{R}_{i\alpha}) + 2\kappa(k-1)\widetilde{R}_{i\alpha}^{k-2} e^{-\kappa R_{i\alpha}} \partial_{i,m}\widetilde{R}_{i\alpha} \right) \right) \end{eqnarray} **** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_jastrow_en_l(qmckl_context context, double* const forces_jastrow_en_l, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_jastrow_en_l(qmckl_context context, double* const forces_jastrow_en_l, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_forces_jastrow_en_l(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 3 * ctx->nucleus.num * ctx->electron.walker.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_en_l", "Array too small. Expected 3*nucl_num*walk_num"); } memcpy(forces_jastrow_en_l, ctx->forces.forces_jastrow_en_l, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_en_l (context, & forces_jastrow_en_l, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: forces_jastrow_en_l(size_max) end function qmckl_get_forces_jastrow_en_l end interface #+end_src **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_en_l(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_en_l(qmckl_context context) { qmckl_exit_code rc; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_jastrow_en_l", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->jastrow_champ.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_forces_jastrow_en_l", NULL); } /* Check if en rescaled distance is provided */ rc = qmckl_provide_en_distance_rescaled(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_en_distance_rescaled_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->date > ctx->forces.forces_jastrow_en_l_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_jastrow_en_l != NULL) { rc = qmckl_free(context, ctx->forces.forces_jastrow_en_l); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_en_l", "Unable to free ctx->forces.forces_jastrow_en_l"); } ctx->forces.forces_jastrow_en_l = NULL; } } /* Allocate array */ if (ctx->forces.forces_jastrow_en_l == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * 3 * ctx->nucleus.num * sizeof(double); double* forces_jastrow_en_l = (double*) qmckl_malloc(context, mem_info); if (forces_jastrow_en_l == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_en_l", NULL); } ctx->forces.forces_jastrow_en_l = forces_jastrow_en_l; } rc = qmckl_compute_forces_jastrow_en_l(context, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.aord_num, ctx->jastrow_champ.a_vector, ctx->jastrow_champ.rescale_factor_en, ctx->electron.en_distance, ctx->jastrow_champ.en_distance_rescaled, ctx->jastrow_champ.en_distance_rescaled_gl, ctx->forces.forces_jastrow_en_l); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_jastrow_en_l_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_en_l :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_jastrow_en_l_args | Variable | Type | In/Out | Description | |---------------------------+-------------------------------------------+--------+---------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~type_nucl_num~ | ~int64_t~ | in | Number of unique nuclei | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | IDs of unique nuclei | | ~aord_num~ | ~int64_t~ | in | Number of coefficients | | ~a_vector~ | ~double[type_nucl_num][aord_num+1]~ | in | List of coefficients | | ~rescale_factor_en~ | ~double[type_nucl_num]~ | in | Rescale factor for electron-nucleus | | ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | in | Electron-nucleus distances | | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus distances | | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in | Electron-nucleus distance derivatives | | ~forces_jastrow_en_l~ | ~double[walk_num][nucl_num][3]~ | out | Forces of electron-nucleus Laplacian | |---------------------------+-------------------------------------------+--------+---------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_forces_jastrow_en_l_doc( & context, walk_num, elec_num, nucl_num, type_nucl_num, & type_nucl_vector, aord_num, a_vector, rescale_factor_en, en_distance, & en_distance_rescaled, en_distance_rescaled_gl, forces_jastrow_en_l) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) integer (c_int64_t) , intent(in) , value :: aord_num real (c_double ) , intent(in) :: a_vector(aord_num+1,type_nucl_num) real (c_double ) , intent(in) :: rescale_factor_en(type_nucl_num) real (c_double ) , intent(in) :: en_distance(nucl_num, elec_num, walk_num) real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num) real (c_double ) , intent(in) :: en_distance_rescaled_gl(4, elec_num,nucl_num,walk_num) real (c_double ) , intent(out) :: forces_jastrow_en_l(3,nucl_num,walk_num) integer(qmckl_exit_code) :: info integer*8 :: i, a, k, nw, ii, m,l double precision :: x, x1, kf, grad_c2 double precision :: denom, invdenom, invdenom2, f, f2, expk, invdist double precision :: dx(4) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (aord_num < 0) then info = QMCKL_INVALID_ARG_7 return endif do nw =1, walk_num forces_jastrow_en_l(:,:,nw) = 0.0d0 do a = 1, nucl_num do i = 1, elec_num expk = dexp(rescale_factor_en(type_nucl_vector(a)+1) * en_distance(a,i,nw)) invdist = 1.d0 / en_distance(a,i,nw) x = en_distance_rescaled(i,a,nw) if(abs(x) < 1.d-12) continue denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x invdenom = 1.0d0 / denom invdenom2 = invdenom*invdenom f = a_vector(1, type_nucl_vector(a)+1) * invdenom2 do m = 1, 4 dx(m) = en_distance_rescaled_gl(m,i,a,nw) end do grad_c2 = dx(1)*dx(1) + dx(2)*dx(2) + dx(3)*dx(3) do m = 1, 3 forces_jastrow_en_l(m,a,nw) = forces_jastrow_en_l(m,a,nw) + f * dx(m) * & (rescale_factor_en(type_nucl_vector(a)+1) & * (2 *invdist - rescale_factor_en(type_nucl_vector(a)+1)) + 2*invdist*invdist) forces_jastrow_en_l(m,a,nw) = forces_jastrow_en_l(m,a,nw) + 2.d0 * f & * invdenom * a_vector(2, type_nucl_vector(a)+1) * dx(m) * dx(4) forces_jastrow_en_l(m,a,nw) = forces_jastrow_en_l(m,a,nw) - 4.d0 * f & * invdenom * a_vector(2, type_nucl_vector(a)+1) * rescale_factor_en(type_nucl_vector(a)+1) & * dx(m)/expk forces_jastrow_en_l(m,a,nw) = forces_jastrow_en_l(m,a,nw) - 6.d0 * f & * invdenom * invdenom * a_vector(2, type_nucl_vector(a)+1) * a_vector(2, type_nucl_vector(a)+1) & * grad_c2 * dx(m) end do kf = 2.d0 x1 = x x = 1.d0 do k=2, aord_num f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x do m = 1, 3 forces_jastrow_en_l(m,a,nw) = forces_jastrow_en_l(m,a,nw) & - f * dx(m) * dx(4) * (kf-1.d0) & - f / x1 * (kf-1.d0) * (kf-2.d0) * dx(m) /expk /expk & + f * x1 * rescale_factor_en(type_nucl_vector(a)+1) * dx(m) * dx(4) * expk & + f * x1 * 2 * dx(m) * invdist * invdist & + 2 * f * (kf-1.d0) * dx(m) * rescale_factor_en(type_nucl_vector(a)+1) / expk end do x = x*x1 kf = kf + 1.d0 end do end do end do end do end function qmckl_compute_forces_jastrow_en_l_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_en_l_doc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* rescale_factor_en, const double* en_distance, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en_l ); qmckl_exit_code qmckl_compute_forces_jastrow_en_l ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* rescale_factor_en, const double* en_distance, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en_l ); qmckl_exit_code qmckl_compute_forces_jastrow_en_l ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* rescale_factor_en, const double* en_distance, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en_l ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_forces_jastrow_en_l (const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* rescale_factor_en, const double* en_distance, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, double* const forces_jastrow_en_l ) { #ifdef HAVE_HPC return qmckl_compute_forces_jastrow_en_l_doc #else return qmckl_compute_forces_jastrow_en_l_doc #endif (context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, aord_num, a_vector, rescale_factor_en, en_distance, en_distance_rescaled, en_distance_rescaled_gl, forces_jastrow_en_l); } #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces Jastrow en L\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double forces_jastrow_en_l[walk_num][nucl_num][3]; rc = qmckl_get_forces_jastrow_en_l(context, &forces_jastrow_en_l[0][0][0], 3*nucl_num*walk_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_en_l[walk_num][nucl_num][3][4][elec_num]; rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_en_gl, &finite_difference_force_en_l[0][0][0][0][0], 4*elec_num); double finite_difference_force_en_l_sum[walk_num][nucl_num][3]; for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ finite_difference_force_en_l_sum[nw][a][k] = 0; for (int i = 0; i < elec_num; i++){ finite_difference_force_en_l_sum[nw][a][k] = finite_difference_force_en_l_sum[nw][a][k] + finite_difference_force_en_l[nw][a][k][3][i]; } } } } for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ //printf("finite_difference_force_en_l_sum[%i][%i][%i] %+3.10f \n", nw,a,k,finite_difference_force_en_l_sum[nw][a][k]); //printf("forces_jastrow_en_l [%i][%i][%i] %+3.10f\n", nw,a,k,forces_jastrow_en_l[nw][a][k]); } } } for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ assert(fabs(finite_difference_force_en_l_sum[nw][a][k] - forces_jastrow_en_l[nw][a][k]) < 1.e-8); } } } printf("OK\n"); #+end_src *** Force of single electron-nucleus Jastrow value Calculate the force of the single-electron contribution ot the electron-nucleus Jastrow value. \[ \delta\partial_{\alpha,m} J_{en} = -\left( \frac{a_1 \partial_{i,m} \widetilde{R}_{i\alpha}^{\text{new}}}{(1+a_2 \widetilde{R}_{i\alpha}^{\text{new}} )^2} + \sum_{k=2}^{N_\text{aord}} a_k k {\widetilde{R}_{i\alpha}^{\text{new}}}^{k-1} \partial_{i,m} \widetilde{R}_{i\alpha}^{\text{new}} \right) + \left( \frac{a_1 \partial_{i,m} \widetilde{R}_{i\alpha}^{\text{old}}}{(1+a_2 \widetilde{R}_{i\alpha}^{\text{old}} )^2} + \sum_{k=2}^{N_\text{aord}} a_k k {\widetilde{R}_{i\alpha}^{\text{old}}}^{k-1} \partial_{i,m} \widetilde{R}_{i\alpha}^{\text{old}} \right) \] The function ~qmckl_set_single_point~ has to be called before this function to set the new electron position. **** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_jastrow_single_en(qmckl_context context, double* const forces_jastrow_single_en, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_jastrow_single_en(qmckl_context context, double* const forces_jastrow_single_en, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_forces_jastrow_single_en(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); const int64_t sze = ctx->electron.walker.num * 3 * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_single_en", "input array too small"); } memcpy(forces_jastrow_single_en, ctx->forces.forces_jastrow_single_en, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_single_en (context, & forces_jastrow_single_en, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: size_max real(c_double), intent(out) :: forces_jastrow_single_en(size_max) end function qmckl_get_forces_jastrow_single_en end interface #+end_src **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_single_en(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :export none qmckl_exit_code qmckl_provide_forces_jastrow_single_en(qmckl_context context) { qmckl_exit_code rc = QMCKL_SUCCESS; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_jastrow_single_en", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->jastrow_champ.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_forces_jastrow_single_en", NULL); } /* Check if en rescaled distance is provided */ rc = qmckl_provide_en_distance_rescaled(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_en_distance_rescaled_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance is provided */ rc = qmckl_provide_en_rescaled_single(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_en_rescaled_single_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->forces.forces_jastrow_single_en_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_jastrow_single_en != NULL) { rc = qmckl_free(context, ctx->forces.forces_jastrow_single_en); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_single_en", "Unable to free ctx->forces.forces_jastrow_single_en"); } ctx->forces.forces_jastrow_single_en = NULL; } } /* Allocate array */ if (ctx->forces.forces_jastrow_single_en == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * 3 * ctx->nucleus.num * sizeof(double); double* forces_jastrow_single_en = (double*) qmckl_malloc(context, mem_info); if (forces_jastrow_single_en == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_single_en", NULL); } ctx->forces.forces_jastrow_single_en = forces_jastrow_single_en; } rc = qmckl_compute_forces_jastrow_single_en(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.aord_num, ctx->jastrow_champ.a_vector, ctx->jastrow_champ.en_distance_rescaled, ctx->jastrow_champ.en_distance_rescaled_gl, ctx->single_point.en_rescaled_single, ctx->single_point.en_rescaled_single_gl, ctx->forces.forces_jastrow_single_en); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_jastrow_single_en_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_single_en :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_single_en_args | Variable | Type | In/Out | Description | |----------------------------+-------------------------------------------+--------+-------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of single electron | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~type_nucl_num~ | ~int64_t~ | in | Number of unique nuclei | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | IDs of unique nuclei | | ~aord_num~ | ~int64_t~ | in | Number of coefficients | | ~a_vector~ | ~double[type_nucl_num][aord_num+1]~ | in | List of coefficients | | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~en_distance_rescaled_gl~ | ~double[walk_num][nucl_num][elec_num][4]~ | in | Electron-nucleus rescaled distance derivatives | | ~en_rescaled_single~ | ~double[walk_num][nucl_num]~ | in | Electron-nucleus rescaled single distances | | ~en_rescaled_single_gl~ | ~double[walk_num][nucl_num][4]~ | in | Electron-nucleus rescaled single distance derivatives | | ~forces_jastrow_single_en~ | ~double[walk_num][nucl_num][3]~ | out | Single electron-nucleus forces | |----------------------------+-------------------------------------------+--------+-------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_forces_jastrow_single_en_doc( & context, num_in, walk_num, elec_num, nucl_num, type_nucl_num, & type_nucl_vector, aord_num, a_vector, & en_distance_rescaled, en_distance_rescaled_gl, en_rescaled_single, en_rescaled_single_gl, forces_jastrow_single_en) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num_in integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: type_nucl_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) integer (c_int64_t) , intent(in) , value :: aord_num real (c_double ) , intent(in) :: a_vector(aord_num+1,type_nucl_num) real (c_double ) , intent(in) :: en_distance_rescaled(elec_num,nucl_num,walk_num) real (c_double ) , intent(in) :: en_distance_rescaled_gl(4, elec_num,nucl_num,walk_num) real (c_double ) , intent(in) :: en_rescaled_single(nucl_num,walk_num) real (c_double ) , intent(in) :: en_rescaled_single_gl(4, nucl_num,walk_num) real (c_double ) , intent(out) :: forces_jastrow_single_en(3,nucl_num,walk_num) integer(qmckl_exit_code) :: info integer*8 :: i, a, k, nw, ii, num double precision :: x, x1, kf, x_old, x1_old double precision :: denom, invdenom, invdenom2, f double precision :: denom_old, invdenom_old, invdenom2_old, f_old double precision :: dx(3), dx_old(3) num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_5 return endif if (aord_num < 0) then info = QMCKL_INVALID_ARG_8 return endif do nw =1, walk_num forces_jastrow_single_en(:,:,nw) = 0.0d0 do a = 1, nucl_num x_old = en_distance_rescaled(num,a,nw) x = en_rescaled_single(a,nw) denom = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x invdenom = 1.0d0 / denom invdenom2 = invdenom*invdenom denom_old = 1.0d0 + a_vector(2, type_nucl_vector(a)+1) * x_old invdenom_old = 1.0d0 / denom_old invdenom2_old = invdenom_old*invdenom_old dx(1) = -en_rescaled_single_gl(1,a,nw) dx(2) = -en_rescaled_single_gl(2,a,nw) dx(3) = -en_rescaled_single_gl(3,a,nw) dx_old(1) = -en_distance_rescaled_gl(1,num,a,nw) dx_old(2) = -en_distance_rescaled_gl(2,num,a,nw) dx_old(3) = -en_distance_rescaled_gl(3,num,a,nw) f = a_vector(1, type_nucl_vector(a)+1) * invdenom2 f_old = a_vector(1, type_nucl_vector(a)+1) * invdenom2_old forces_jastrow_single_en(1,a,nw) = forces_jastrow_single_en(1,a,nw) + f * dx(1) - f_old * dx_old(1) forces_jastrow_single_en(2,a,nw) = forces_jastrow_single_en(2,a,nw) + f * dx(2) - f_old * dx_old(2) forces_jastrow_single_en(3,a,nw) = forces_jastrow_single_en(3,a,nw) + f * dx(3) - f_old * dx_old(3) kf = 2.d0 x1 = x x = 1.d0 x1_old = x_old x_old = 1.d0 do k=2, aord_num f = a_vector(k+1,type_nucl_vector(a)+1) * kf * x f_old = a_vector(k+1,type_nucl_vector(a)+1) * kf * x_old forces_jastrow_single_en(1,a,nw) = forces_jastrow_single_en(1,a,nw) + f * x1 * dx(1) - f_old * x1_old * dx_old(1) forces_jastrow_single_en(2,a,nw) = forces_jastrow_single_en(2,a,nw) + f * x1 * dx(2) - f_old * x1_old * dx_old(2) forces_jastrow_single_en(3,a,nw) = forces_jastrow_single_en(3,a,nw) + f * x1 * dx(3) - f_old * x1_old * dx_old(3) x = x*x1 x_old = x_old*x1_old kf = kf + 1.d0 end do end do end do end function qmckl_compute_forces_jastrow_single_en_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_single_en_doc ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, const double* en_rescaled_single, const double* en_rescaled_single_gl, double* const forces_jastrow_single_en ); qmckl_exit_code qmckl_compute_forces_jastrow_single_en ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, const double* en_rescaled_single, const double* en_rescaled_single_gl, double* const forces_jastrow_single_en ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_forces_jastrow_single_en (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, const int64_t* type_nucl_vector, const int64_t aord_num, const double* a_vector, const double* en_distance_rescaled, const double* en_distance_rescaled_gl, const double* en_rescaled_single, const double* en_rescaled_single_gl, double* const forces_jastrow_single_en ) { #ifdef HAVE_HPC return qmckl_compute_forces_jastrow_single_en_doc #else return qmckl_compute_forces_jastrow_single_en_doc #endif (context, num, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, aord_num, a_vector, en_distance_rescaled, en_distance_rescaled_gl, en_rescaled_single, en_rescaled_single_gl, forces_jastrow_single_en ); } #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces Jastrow single en\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double new_coords[6] = {1.0,2.0,3.0,4.0,5.0,6.0}; rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3*walk_num); assert (rc == QMCKL_SUCCESS); double forces_jastrow_single_en[walk_num][nucl_num][3]; rc = qmckl_get_forces_jastrow_single_en(context, &forces_jastrow_single_en[0][0][0], 3*nucl_num*walk_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_single_en[walk_num][nucl_num][3]; rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_single_en, &(finite_difference_force_single_en[0][0][0]), 1); for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ //printf("%.10f\t", finite_difference_force_single_en[nw][a][k]); //printf("%.10f\n", forces_jastrow_single_en[nw][a][k]); assert(fabs(finite_difference_force_single_en[nw][a][k] - forces_jastrow_single_en[nw][a][k]) < 1.e-8); } } } printf("OK\n"); #+end_src ** Electron-electron-nucleus Jastrow *** Force of P matrix Calculates the force of the P matrix. This is a required component to calculate the force on the 3-body Jastrow. \[ \partial_{\alpha,m} P_{i,\alpha,k,l} = \sum_j \widetilde{r}_{i,j,k} \partial_{\alpha,m}\widetilde{R}_{j,\alpha,l} = -\sum_j \widetilde{r}_{i,j,k} \partial_{j,m}\widetilde{R}_{j,\alpha,l} \] **** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_tmp_c(qmckl_context context, double* const forces_tmp_c, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_tmp_c(qmckl_context context, double* const forces_tmp_c, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_forces_tmp_c(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 4 * ctx->electron.walker.num * ctx->electron.num * ctx->nucleus.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num+1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_tmp_c", "Array too small. Expected 4*walk_num*elec_num*nucl_num*cord_num*(cord_num+1)"); } memcpy(forces_tmp_c, ctx->forces.forces_tmp_c, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_tmp_c(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_tmp_c(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->jastrow_champ.cord_num > 0) { /* Check if en rescaled distance is provided */ rc = qmckl_provide_een_rescaled_e(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance gl derivatives is provided */ rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; } /* Compute if necessary */ if (ctx->date > ctx->forces.forces_tmp_c_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_tmp_c != NULL) { rc = qmckl_free(context, ctx->forces.forces_tmp_c); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_tmp_c", "Unable to free ctx->forces.forces_tmp_c"); } ctx->forces.forces_tmp_c = NULL; } } /* Allocate array */ if (ctx->forces.forces_tmp_c == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 4 * ctx->electron.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num+1) * ctx->nucleus.num * ctx->electron.walker.num * sizeof(double); double* forces_tmp_c = (double*) qmckl_malloc(context, mem_info); if (forces_tmp_c == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_tmp_c", NULL); } ctx->forces.forces_tmp_c = forces_tmp_c; } rc = qmckl_compute_forces_tmp_c(context, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.een_rescaled_e, ctx->jastrow_champ.een_rescaled_n_gl, ctx->forces.forces_tmp_c); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_tmp_c_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_tmp_c :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_tmp_c_args | Variable | Type | In/Out | Description | |---------------------+---------------------------------------------------------------------+--------+---------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-nucleus rescaled | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled derivatives | | ~forces_tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | out | Force on the P matrix | |---------------------+---------------------------------------------------------------------+--------+---------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_tmp_c( & context, walk_num, elec_num, nucl_num, cord_num,& een_rescaled_e, een_rescaled_n_gl, forces_tmp_c) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer (c_int64_t) , intent(in), value :: walk_num, elec_num, cord_num, nucl_num real (c_double ) , intent(in) :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(out) :: forces_tmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) integer*8 :: nw, i integer*8 :: l, m, k, a,j info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return do nw=1, walk_num do i=0, cord_num-1 info = qmckl_dgemm(context,'N','N',elec_num*1_8,& nucl_num*(cord_num+1)*4_8, elec_num*1_8, -1.0d0, & een_rescaled_e(1,1,i,nw),1_8*elec_num, & een_rescaled_n_gl(1,1,1,0,nw),elec_num*1_8, & 0.0d0, & forces_tmp_c(1,1,1,0,i,nw),1_8*elec_num) end do end do end function qmckl_compute_forces_tmp_c #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_tmp_c ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_e, const double* een_rescaled_n_gl, double* const forces_tmp_c ); #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces tmp_c\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double forces_tmp_c[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num]; rc = qmckl_get_forces_tmp_c(context, &forces_tmp_c[0][0][0][0][0][0], 4*nucl_num*walk_num*elec_num*(cord_num+1)*cord_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_tmp_c[walk_num][cord_num][cord_num+1][nucl_num][3][elec_num]; double* nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (nucleus_coord == NULL) { return QMCKL_ALLOCATION_FAILED; } rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num); double* temp_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (temp_coord == NULL) { free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } double output[walk_num][cord_num][cord_num+1][nucl_num][elec_num]; // Copy original coordinates for (int i = 0; i < 3 * nucl_num; i++) { temp_coord[i] = nucleus_coord[i]; } for (int64_t a = 0; a < nucl_num; a++) { for (int64_t k = 0; k < 3; k++) { for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x; // Update coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); // Call the provided function rc = qmckl_get_jastrow_champ_tmp_c(context, &output[0][0][0][0][0], 4*nucl_num*walk_num*elec_num*(cord_num+1)*cord_num); assert(rc == QMCKL_SUCCESS); // Accumulate derivative using finite-difference coefficients for (int nw=0 ; nwelectron.walker.num * ctx->electron.num * ctx->nucleus.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num+1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_dtmp_c", "Array too small. Expected 4*4*walk_num*elec_num*nucl_num*cord_num*(cord_num+1)"); } memcpy(forces_dtmp_c, ctx->forces.forces_dtmp_c, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_dtmp_c(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_dtmp_c(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->jastrow_champ.cord_num > 0) { /* Check if en rescaled distance is provided */ rc = qmckl_provide_een_rescaled_e_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance gl derivatives is provided */ rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; } /* Compute if necessary */ if (ctx->date > ctx->forces.forces_dtmp_c_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_dtmp_c != NULL) { rc = qmckl_free(context, ctx->forces.forces_dtmp_c); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_dtmp_c", "Unable to free ctx->forces.forces_dtmp_c"); } ctx->forces.forces_dtmp_c = NULL; } } /* Allocate array */ if (ctx->forces.forces_dtmp_c == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 4* 4 * ctx->electron.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num+1) * ctx->nucleus.num * ctx->electron.walker.num * sizeof(double); double* forces_dtmp_c = (double*) qmckl_malloc(context, mem_info); if (forces_dtmp_c == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_dtmp_c", NULL); } ctx->forces.forces_dtmp_c = forces_dtmp_c; } rc = qmckl_compute_forces_dtmp_c(context, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.een_rescaled_e_gl, ctx->jastrow_champ.een_rescaled_n_gl, ctx->forces.forces_dtmp_c); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_dtmp_c_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_dtmp_c :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_dtmp_c_args | Variable | Type | In/Out | Description | |---------------------+------------------------------------------------------------------------+--------+---------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~een_rescaled_e_gl~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | in | Electron-nucleus rescaled | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled derivatives | | ~forces_dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][4][elec_num]~ | out | Force of derivative of the P matrix | |---------------------+------------------------------------------------------------------------+--------+---------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_dtmp_c_doc( & context, walk_num, elec_num, nucl_num, cord_num,& een_rescaled_e_gl, een_rescaled_n_gl, forces_dtmp_c) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer (c_int64_t) , intent(in), value :: walk_num, elec_num, cord_num, nucl_num real (c_double ) , intent(in) :: een_rescaled_e_gl(elec_num, 4, elec_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(out) :: forces_dtmp_c(elec_num, 4, nucl_num,4,0:cord_num, 0:cord_num-1, walk_num) integer*8 :: nw, i, k integer*8 :: l, m, a,j, kk info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return do nw = 1, walk_num do l = 0, cord_num-1 do m = 1, cord_num do k = 1, 4 do a = 1, nucl_num do kk = 1, 4 do i = 1, elec_num forces_dtmp_c(i,kk,a,k,m,l,nw) = 0.d0 do j = 1, elec_num forces_dtmp_c(i,kk,a,k,m,l,nw) = forces_dtmp_c(i,kk,a,k,m,l,nw) - & een_rescaled_e_gl(i,kk,j,l,nw) * een_rescaled_n_gl(j,k,a,m,nw) end do end do end do end do end do end do end do end do end function qmckl_compute_forces_dtmp_c_doc #+end_src #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_dtmp_c_hpc( & context, walk_num, elec_num, nucl_num, cord_num,& een_rescaled_e_gl, een_rescaled_n_gl, forces_dtmp_c) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer (c_int64_t) , intent(in), value :: walk_num, elec_num, cord_num, nucl_num real (c_double ) , intent(in) :: een_rescaled_e_gl(elec_num, 4, elec_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(out) :: forces_dtmp_c(elec_num, 4, nucl_num,4,0:cord_num, 0:cord_num-1, walk_num) integer*8 :: nw, i, k integer*8 :: l, m, a,j, kk double precision, allocatable :: tmp(:,:,:,:) info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return allocate(tmp(elec_num, 4, 4, nucl_num)) do nw = 1, walk_num do l = 0, cord_num-1 do m = 1, cord_num info = qmckl_dgemm(context,'N','N', elec_num*4, 4*nucl_num, elec_num, -1.d0, & een_rescaled_e_gl(1,1,1,l,nw), elec_num*4_8, & een_rescaled_n_gl(1,1,1,m,nw), elec_num, 0.d0, & tmp, elec_num*4_8) do k = 1, 4 do a = 1, nucl_num forces_dtmp_c(:,:,a,k,m,l,nw) = tmp(:,:,k,a) end do end do end do end do end do deallocate(tmp) end function qmckl_compute_forces_dtmp_c_hpc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_dtmp_c_doc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_e_gl, const double* een_rescaled_n_gl, double* const forces_dtmp_c ); qmckl_exit_code qmckl_compute_forces_dtmp_c_hpc ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_e_gl, const double* een_rescaled_n_gl, double* const forces_dtmp_c ); qmckl_exit_code qmckl_compute_forces_dtmp_c ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_e_gl, const double* een_rescaled_n_gl, double* const forces_dtmp_c ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_forces_dtmp_c ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_e_gl, const double* een_rescaled_n_gl, double* const forces_dtmp_c ) { #ifdef HAVE_HPC return qmckl_compute_forces_dtmp_c_hpc #else return qmckl_compute_forces_dtmp_c_doc #endif (context, walk_num, elec_num, nucl_num, cord_num, een_rescaled_e_gl, een_rescaled_n_gl, forces_dtmp_c ); } #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces dtmp_c\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double forces_dtmp_c[walk_num][cord_num][cord_num+1][4][nucl_num][4][elec_num]; rc = qmckl_get_forces_dtmp_c(context, &forces_dtmp_c[0][0][0][0][0][0][0], 4*4*nucl_num*walk_num*elec_num*(cord_num+1)*cord_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_dtmp_c[walk_num][cord_num][cord_num+1][3][nucl_num][4][elec_num]; nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (nucleus_coord == NULL) { return QMCKL_ALLOCATION_FAILED; } rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num); temp_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (temp_coord == NULL) { free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } double doutput[walk_num][cord_num][cord_num+1][nucl_num][4][elec_num]; // Copy original coordinates for (int i = 0; i < 3 * nucl_num; i++) { temp_coord[i] = nucleus_coord[i]; } for (int64_t a = 0; a < nucl_num; a++) { for (int64_t k = 0; k < 3; k++) { for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x; // Update coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); // Call the provided function rc = qmckl_get_jastrow_champ_dtmp_c(context, &doutput[0][0][0][0][0][0], 4*4*nucl_num*walk_num*elec_num*(cord_num+1)*cord_num); assert(rc == QMCKL_SUCCESS); // Accumulate derivative using finite-difference coefficients for (int nw=0 ; nwnucleus.num * ctx->electron.walker.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_een", "Array too small. Expected 3*nucl_num*walk_num"); } memcpy(forces_jastrow_een, ctx->forces.forces_jastrow_een, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_een (context, & forces_jastrow_een, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: forces_jastrow_een(size_max) end function qmckl_get_forces_jastrow_een end interface #+end_src **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_een(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_een(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->jastrow_champ.cord_num > 0) { /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance gl derivatives is provided */ rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if tmp_c is provided */ rc = qmckl_provide_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if forces_tmp_c is provided */ rc = qmckl_provide_forces_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_jastrow_champ_c_vector_full(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_lkpm_combined_index(context); if(rc != QMCKL_SUCCESS) return rc; } /* Compute if necessary */ if (ctx->date > ctx->forces.forces_jastrow_een_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_jastrow_een != NULL) { rc = qmckl_free(context, ctx->forces.forces_jastrow_een); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_een", "Unable to free ctx->forces.forces_jastrow_een"); } ctx->forces.forces_jastrow_een = NULL; } } /* Allocate array */ if (ctx->forces.forces_jastrow_een == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->nucleus.num * ctx->electron.walker.num * sizeof(double); double* forces_jastrow_een = (double*) qmckl_malloc(context, mem_info); if (forces_jastrow_een == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_een", NULL); } ctx->forces.forces_jastrow_een = forces_jastrow_een; } rc = qmckl_compute_forces_jastrow_een(context, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.dim_c_vector, ctx->jastrow_champ.c_vector_full, ctx->jastrow_champ.lkpm_combined_index, ctx->jastrow_champ.een_rescaled_n, ctx->jastrow_champ.een_rescaled_n_gl, ctx->jastrow_champ.tmp_c, ctx->forces.forces_tmp_c, ctx->forces.forces_jastrow_een); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_jastrow_een_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_een :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_jastrow_een_args | Variable | Type | In/Out | Description | |-----------------------+---------------------------------------------------------------------+--------+--------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~dim_c_vector~ | ~int64_t~ | in | dimension of full coefficient vector | | ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | in | full coefficient vector | | ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | in | combined indices | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled derivatives | | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | P matrix | | ~forces_tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | Force of P matrix | | ~forces_jastrow_een~ | ~double[walk_num][nucl_num][3]~ | out | Force of electron-electron-nucleus Jastrow value | |-----------------------+---------------------------------------------------------------------+--------+--------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_een( & context, walk_num, elec_num, nucl_num, cord_num,& dim_c_vector, c_vector_full, lkpm_combined_index, & een_rescaled_n, een_rescaled_n_gl, tmp_c, forces_tmp_c,forces_jastrow_een) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer (c_int64_t) , intent(in), value :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real (c_double ) , intent(in) :: c_vector_full(nucl_num, dim_c_vector) real (c_double ) , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_tmp_c(elec_num,4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(out) :: forces_jastrow_een(3, nucl_num, walk_num) integer*8 :: i, a, j, l, k, p, m, n, nw, ii double precision :: accu, accu2, cn info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return forces_jastrow_een = 0.d0 if (cord_num == 0) return do nw =1, walk_num do n = 1, dim_c_vector l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) p = lkpm_combined_index(n, 3) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = c_vector_full(a, n) if(cn == 0.d0) cycle do j = 1, elec_num do ii = 1, 3 accu = een_rescaled_n(j,a,m,nw) * forces_tmp_c(j,ii,a,m+l,k,nw) & - een_rescaled_n_gl(j,ii,a,m,nw) * tmp_c(j,a,m+l,k,nw) forces_jastrow_een(ii, a, nw) = forces_jastrow_een(ii, a, nw) + accu * cn end do end do end do end do end do end function qmckl_compute_forces_jastrow_een #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_een ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* een_rescaled_n, const double* een_rescaled_n_gl, const double* tmp_c, const double* forces_tmp_c, double* const forces_jastrow_een ); #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces Jastrow een\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double forces_jastrow_een[walk_num][nucl_num][3]; rc = qmckl_get_forces_jastrow_een(context, &forces_jastrow_een[0][0][0], 3*nucl_num*walk_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_een[walk_num][nucl_num][3]; rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_een, &(finite_difference_force_een[0][0][0]), 1); for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ //printf("%.10f\t", finite_difference_force_een[nw][a][k]); //printf("%.10f\n", forces_jastrow_een[nw][a][k]); } } } for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ assert(fabs(finite_difference_force_een[nw][a][k] - forces_jastrow_een[nw][a][k]) < 1.e-8); } } } printf("OK\n"); #+end_src *** Force of derivatives of electron-nucleus rescaled distance Calculates the force of the derivatives of the electron-nucleus rescaled distance. \begin{eqnarray} \partial_{\alpha,m}\partial_{i,n} \widetilde{R}_{i\alpha} =& \partial_{\alpha,m} \left( -\frac{r_{i,n} - R_{\alpha,n}}{R_{i\alpha}} \kappa l e^{-\kappa l R_{i\alpha}}\right)\\ =& - \frac{(r_{i,n} - R_{\alpha,n})(r_{i,m} - R_{\alpha,m})}{R_{i\alpha}^3}l\kappa e^{-\kappa l R_{i\alpha}} - \frac{(r_{i,n} - R_{\alpha,n})(r_{i,m} - R_{\alpha,m})}{R_{i\alpha}^2} \kappa^2 l^2 e^{-\kappa l R_{i\alpha}} + \delta_{n,m} \frac{1}{r}\kappa l e^{-\kappa l R_{i\alpha}}\\ \partial_{\alpha,m}\nabla^2 \widetilde{R}_{i\alpha} =& \partial_{\alpha,m} \left( \kappa l \left( \kappa l - \frac{2}{R_{i\alpha}}\right) e^{-\kappa l R_{i\alpha}}\right)\\ =&-2\kappa l \frac{r_{i,m} - R_{\alpha,m}}{R_{i\alpha}^3}e^{-\kappa l R_{i\alpha}} + \kappa^2 l^2 \left(\kappa l - \frac{2}{R_{i\alpha}} \right) \frac{r_{i,m}-R_{\alpha,m}}{R_{i\alpha}} e^{-\kappa l R_{i\alpha}} \end{eqnarray} **** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_een_rescaled_n_gl(qmckl_context context, double* const forces_een_n, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_een_rescaled_n_gl(qmckl_context context, double* const forces_een_n, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_forces_een_rescaled_n_gl(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 3* ctx->electron.num * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1); if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_een_rescaled_n_gl", "Array too small. Expected ctx->electron.num * 3 * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1)"); } memcpy(forces_een_n, ctx->forces.forces_een_n, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_een_rescaled_n_gl(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_een_rescaled_n_gl(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); /* Check if ee distance is provided */ qmckl_exit_code rc = qmckl_provide_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if ee distance is provided */ rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if ee gl distance is provided */ rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->date > ctx->forces.forces_een_n_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_een_n != NULL) { rc = qmckl_free(context, ctx->forces.forces_een_n); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_een_rescaled_n_gl", "Unable to free ctx->forces.forces_een_n"); } ctx->forces.forces_een_n = NULL; } } /* Allocate array */ if (ctx->forces.forces_een_n == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.num * 3 * 4 * ctx->nucleus.num * ctx->electron.walker.num * (ctx->jastrow_champ.cord_num + 1) * sizeof(double); double* forces_een_n = (double*) qmckl_malloc(context, mem_info); if (forces_een_n == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_een_rescaled_n_gl", NULL); } ctx->forces.forces_een_n = forces_een_n; } rc = qmckl_compute_forces_een_rescaled_n_gl(context, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.type_nucl_num, ctx->jastrow_champ.type_nucl_vector, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.rescale_factor_en, ctx->electron.en_distance, ctx->jastrow_champ.een_rescaled_n, ctx->jastrow_champ.een_rescaled_n_gl, ctx->forces.forces_een_n); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_een_n_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_een_rescaled_n_gl :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_compute_forces_een_rescaled_n_gl_args | Variable | Type | In/Out | Description | |---------------------+----------------------------------------------------------+--------+----------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of atoms | | ~type_nucl_num~ | ~int64_t~ | in | Number of atom types | | ~type_nucl_vector~ | ~int64_t[nucl_num]~ | in | Types of atoms | | ~cord_num~ | ~int64_t~ | in | Order of polynomials | | ~rescale_factor_en~ | ~double[nucl_num]~ | in | Factor to rescale ee distances | | ~en_distance~ | ~double[walk_num][elec_num][nucl_num]~ | in | Electron-nucleus distances | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances derivativies | | ~forces_een_n~ | ~double[walk_num][0:cord_num][3][nucl_num][4][elec_num]~ | out | Force of electron-nucleus rescaled distances derivatives | |---------------------+----------------------------------------------------------+--------+----------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_een_rescaled_n_gl( & context, walk_num, elec_num, nucl_num, type_nucl_num, type_nucl_vector, & cord_num, rescale_factor_en, & en_distance, een_rescaled_n, een_rescaled_n_gl,forces_een_n) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer (c_int64_t) , intent(in),value :: walk_num, elec_num, nucl_num, type_nucl_num, cord_num integer (c_int64_t) , intent(in) :: type_nucl_vector(nucl_num) real (c_double ) , intent(in) :: rescale_factor_en(type_nucl_num) real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num) real (c_double ) , intent(out) :: forces_een_n(elec_num,4,nucl_num,3,0:cord_num,walk_num) double precision :: x, ria_inv, kappa_l, temp integer*8 :: i, a, k, l, nw, m, n info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then info = QMCKL_INVALID_CONTEXT return endif if (walk_num <= 0) then info = QMCKL_INVALID_ARG_2 return endif if (elec_num <= 0) then info = QMCKL_INVALID_ARG_3 return endif if (nucl_num <= 0) then info = QMCKL_INVALID_ARG_4 return endif if (cord_num < 0) then info = QMCKL_INVALID_ARG_5 return endif !forces_een_n = 0.0d0 !do nw = 1, walk_num ! do l = 0, cord_num ! do a = 1, nucl_num ! kappa_l = dble(l)*rescale_factor_en(type_nucl_vector(a)+1) ! do i = 1, elec_num ! ria_inv = 1.0d0 / en_distance(a,i,nw) ! do n = 1, 3 ! do m = 1, 3 ! if (m == n) then ! forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) + & ! kappa_l*een_rescaled_n(i,a,l,nw) * ria_inv ! end if ! temp = een_rescaled_n_gl(i,m,a,l,nw) * een_rescaled_n_gl(i,n,a,l,nw) / een_rescaled_n(i,a,l,nw) ! forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) - & ! temp * ria_inv / kappa_l - temp ! end do ! forces_een_n(i,4,a,n,l,nw) = forces_een_n(i,4,a,n,l,nw) + & ! (2.0d0 * ria_inv * ria_inv & ! - een_rescaled_n_gl(i,4,a,l,nw)/een_rescaled_n(i,a,l,nw)) * & ! een_rescaled_n_gl(i,n,a,l,nw) ! end do ! end do ! end do ! end do !end do forces_een_n = 0.0d0 do nw = 1, walk_num do i = 1, elec_num do a = 1, nucl_num do l = 0, cord_num kappa_l = dble(l)*rescale_factor_en(type_nucl_vector(a)+1) do m = 1, 4 do n = 1, 3 if (l == 0) then forces_een_n(i,m,a,n,l,nw) = 0.0d0 else if (m == n) then forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) + & kappa_l*een_rescaled_n(i,a,l,nw)/en_distance(a,i,nw) end if if (m < 4) then forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) - & een_rescaled_n_gl(i,m,a,l,nw) * een_rescaled_n_gl(i,n,a,l,nw) / & (kappa_l*een_rescaled_n(i,a,l,nw)*en_distance(a,i,nw)) - & een_rescaled_n_gl(i,m,a,l,nw) * & een_rescaled_n_gl(i,n,a,l,nw)/een_rescaled_n(i,a,l,nw) else forces_een_n(i,m,a,n,l,nw) = forces_een_n(i,m,a,n,l,nw) + & 2.0d0 * een_rescaled_n_gl(i,n,a,l,nw) / & (en_distance(a,i,nw)*en_distance(a,i,nw)) - & een_rescaled_n_gl(i,m,a,l,nw) * & een_rescaled_n_gl(i,n,a,l,nw)/een_rescaled_n(i,a,l,nw) end if end if end do end do end do end do end do end do end function qmckl_compute_forces_een_rescaled_n_gl #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_een_rescaled_n_gl ( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t type_nucl_num, int64_t* const type_nucl_vector, const int64_t cord_num, const double* rescale_factor_en, const double* en_distance, const double* een_rescaled_n, const double* een_rescaled_n_gl, double* const forces_een_n ); #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces of een_rescaled_n_gl\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double forces_een_n[walk_num][cord_num+1][3][nucl_num][4][elec_num]; rc = qmckl_get_forces_een_rescaled_n_gl(context, &forces_een_n[0][0][0][0][0][0], 3*4*nucl_num*walk_num*elec_num*(cord_num+1)); assert(rc == QMCKL_SUCCESS); double finite_difference_force_een_n[walk_num][cord_num+1][3][nucl_num][4][elec_num]; nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (nucleus_coord == NULL) { return QMCKL_ALLOCATION_FAILED; } rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num); temp_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (temp_coord == NULL) { free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } double ddoutput[walk_num][cord_num+1][nucl_num][4][elec_num]; // Copy original coordinates for (int i = 0; i < 3 * nucl_num; i++) { temp_coord[i] = nucleus_coord[i]; } for (int64_t a = 0; a < nucl_num; a++) { for (int64_t k = 0; k < 3; k++) { for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x; // Update coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); // Call the provided function rc = qmckl_get_jastrow_champ_een_rescaled_n_gl(context,&ddoutput[0][0][0][0][0], 4*elec_num*nucl_num*(cord_num+1)*walk_num); assert(rc == QMCKL_SUCCESS); // Accumulate derivative using finite-difference coefficients for (int nw=0 ; nwelectron.walker.num * 3 * ctx->electron.num * 3 * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_een_g", "Array too small. Expected 3*3*walk_num*elec_num*nucl_num"); } memcpy(forces_jastrow_een_g, ctx->forces.forces_jastrow_een_g, sze*sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_een_g (context, & forces_jastrow_een_g, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: forces_jastrow_een_g(size_max) end function qmckl_get_forces_jastrow_een_g end interface #+end_src # **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_een_g(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_een_g(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->jastrow_champ.cord_num > 0) { /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_jastrow_champ_c_vector_full(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_lkpm_combined_index(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if tmp_c is provided */ rc = qmckl_provide_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if dtmp_c is provided */ rc = qmckl_provide_dtmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if forces_tmp_c is provided */ rc = qmckl_provide_forces_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if forces_dtmp_c is provided */ rc = qmckl_provide_forces_dtmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if ne distance is provided */ rc = qmckl_provide_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_forces_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; } /* Compute if necessary */ if (ctx->date > ctx->forces.forces_jastrow_een_g_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_jastrow_een_g != NULL) { rc = qmckl_free(context, ctx->forces.forces_jastrow_een_g); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_een_g", "Unable to free ctx->forces.forces_jastrow_een_g"); } ctx->forces.forces_jastrow_een_g = NULL; } } /* Allocate array */ if (ctx->forces.forces_jastrow_een_g == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * 3 * ctx->electron.num * ctx->electron.walker.num * ctx->nucleus.num * sizeof(double); double* forces_jastrow_een_g = (double*) qmckl_malloc(context, mem_info); if (forces_jastrow_een_g == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_een_g", NULL); } ctx->forces.forces_jastrow_een_g = forces_jastrow_een_g; } rc = qmckl_compute_forces_jastrow_een_g(context, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.dim_c_vector, ctx->jastrow_champ.c_vector_full, ctx->jastrow_champ.lkpm_combined_index, ctx->electron.en_distance, ctx->jastrow_champ.tmp_c, ctx->jastrow_champ.dtmp_c, ctx->forces.forces_tmp_c, ctx->forces.forces_dtmp_c, ctx->forces.forces_een_n, ctx->jastrow_champ.een_rescaled_n, ctx->jastrow_champ.een_rescaled_n_gl, ctx->forces.forces_jastrow_een_g); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_jastrow_een_g_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_een_g :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_jastrow_een_g_args | Variable | Type | In/Out | Description | |------------------------+------------------------------------------------------------------------+--------+----------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~dim_c_vector~ | ~int64_t~ | in | dimension of full coefficient vector | | ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | in | full coefficient vector | | ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | in | combined indices | | ~en_distance~ | ~double[elec_num][nucl_num]~ | in | Electron-nucleus distance | | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | P matrix | | ~dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | Derivative of P matrix | | ~forces_tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | Force of P matrix | | ~forces_dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][4][elec_num]~ | in | Force of derivative of P matrix | | ~forces_een_n~ | ~double[walk_num][0:cord_num][3][nucl_num][4][elec_num]~ | in | Force of derivatives of electron-nucleus rescaled factor | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled factor | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Derivative of electron-nucleus rescaled factor | | ~forces_jastrow_een_g~ | ~double[walk_num][3][nucl_num][3][elec_num]~ | out | Force of the electron-electron-nucleus Jastrow gradient | |------------------------+------------------------------------------------------------------------+--------+----------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_een_g( & context, walk_num, elec_num, nucl_num, & cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, & en_distance, tmp_c, dtmp_c, forces_tmp_c, forces_dtmp_c, forces_een_n, een_rescaled_n, & een_rescaled_n_gl, forces_jastrow_een_g)& result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer (c_int64_t) , intent(in),value :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real (c_double ) , intent(in) :: c_vector_full(nucl_num, dim_c_vector) real (c_double ) , intent(in) :: en_distance(nucl_num, elec_num) real (c_double ) , intent(in) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_tmp_c(elec_num, 4,nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_dtmp_c(elec_num, 4, nucl_num,4,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_een_n(elec_num, 4, nucl_num, 3, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(out) :: forces_jastrow_een_g(elec_num,3,nucl_num,3,walk_num) integer*8 :: i, a, j, l, k, m, n, nw, ii, jj double precision :: accu, accu2, cn info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return if (cord_num == 0) return forces_jastrow_een_g = 0.0d0 do nw =1, walk_num do n = 1, dim_c_vector l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = c_vector_full(a, n) if(cn == 0.d0) cycle do ii = 1, 3 do i = 1, elec_num do jj = 1, 3 forces_jastrow_een_g(i, ii, a, jj, nw) = forces_jastrow_een_g(i, ii, a, jj, nw) + (& tmp_c(i,a, m, k,nw) * forces_een_n(i,ii,a,jj ,m+l,nw) & + forces_tmp_c(i,jj,a, m, k,nw) * een_rescaled_n_gl(i,ii,a,m+l,nw) & + tmp_c(i,a, m+l,k,nw) * forces_een_n(i,ii,a,jj, m, nw) & + forces_tmp_c(i,jj,a, m+l,k,nw) * een_rescaled_n_gl(i,ii,a,m, nw) & - dtmp_c(i,ii,a, m, k,nw) * een_rescaled_n_gl(i,jj,a,m+l,nw) & + forces_dtmp_c(i,ii,a,jj,m, k,nw) * een_rescaled_n(i,a, m+l,nw) & - dtmp_c(i,ii,a, m+l,k,nw) * een_rescaled_n_gl(i,jj,a,m, nw) & + forces_dtmp_c(i,ii,a,jj,m+l,k,nw) * een_rescaled_n(i,a, m, nw) & ) * cn end do end do end do end do end do end do end function qmckl_compute_forces_jastrow_een_g #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_een_g( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* en_distance, const double* tmp_c, const double* dtmp_c, const double* forces_tmp_c, const double* forces_dtmp_c, const double* forces_een_n, const double* een_rescaled_n, const double* een_rescaled_n_gl, double* const forces_jastrow_een_g ); #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces Jastrow een G\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double forces_jastrow_een_g[walk_num][3][nucl_num][3][elec_num]; rc = qmckl_get_forces_jastrow_een_g(context, &forces_jastrow_een_g[0][0][0][0][0], 3*3*nucl_num*walk_num*elec_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_een_g[walk_num][nucl_num][3][4][elec_num]; rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_een_gl, &finite_difference_force_een_g[0][0][0][0][0], 4*elec_num); for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ for (int i = 0; i < elec_num; i++){ for (int l = 0; l < 3; l++){ //printf("finite_difference_force_een_g[%i][%i][%i][%i][%i] %+3.10f \n", nw,a,k,l,i,finite_difference_force_een_g[nw][a][k][l][i]); //printf("forces_jastrow_een_g [%i][%i][%i][%i][%i] %+3.10f\n", nw,k,a,l,i,forces_jastrow_een_g[nw][k][a][l][i]); //assert(fabs(finite_difference_force_een_g[nw][a][k][l][i] - forces_jastrow_een_g[nw][k][a][l][i]) < 1.e-8); } } } } } for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ for (int i = 0; i < elec_num; i++){ for (int l = 1; l < 3; l++){ assert(fabs(finite_difference_force_een_g[nw][a][k][l][i] - forces_jastrow_een_g[nw][k][a][l][i]) < 1.e-8); } } } } } printf("OK\n"); #+end_src *** Force of electron-electron-nucleus Jastrow Laplacian Calculates the force of the electron-electron-nucleus Jastrow Laplacian. \begin{eqnarray} \partial_{\alpha,m} \nabla^2 J_{een} &=& \partial_{\alpha,m} \sum_{i=1}^{N_\text{elec}}\sum_{p=2}^{N_\text{nord}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}} c_{l,k,p,\alpha} \left(\nabla^2\widetilde{R}_{i\alpha}^{(p-k+l)/2} P_{i\alpha,k,(p-k-l)/2} + \nabla^2\widetilde{R}_{i\alpha}^{(p-k-l)/2} P_{i\alpha,k,(p-k+l)/2} + \widetilde{R}_{i\alpha}^{(p-k+l)/2} \nabla^2 P_{i\alpha,k,(p-k-l)/2} + \widetilde{R}_{i\alpha}^{(p-k-l)/2} \nabla^2 P_{i\alpha,k,(p-k+l)/2} \right.\\ &&+ \left. 2 \nabla \widetilde{R}_{i\alpha}^{(p-k+l)/2} \cdot \nabla P_{i\alpha,k,(p-k-l)/2} + 2 \nabla \widetilde{R}_{i\alpha}^{(p-k-l)/2} \cdot \nabla P_{i\alpha,k,(p-k+l)/2} \right)\\ &=& \sum_{i=1}^{N_\text{elec}}\sum_{p=2}^{N_\text{nord}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}} c_{l,k,p,\alpha} \left( \partial_{\alpha,m} \nabla^2\widetilde{R}_{i\alpha}^{(p-k+l)/2} P_{i\alpha,k,(p-k-l)/2} + \partial_{\alpha,m} \nabla^2\widetilde{R}_{i\alpha}^{(p-k-l)/2} P_{i\alpha,k,(p-k+l)/2} + \partial_{\alpha,m} \widetilde{R}_{i\alpha}^{(p-k+l)/2} \nabla^2 P_{i\alpha,k,(p-k-l)/2} + \partial_{\alpha,m} \widetilde{R}_{i\alpha}^{(p-k-l)/2} \nabla^2 P_{i\alpha,k,(p-k+l)/2} \right.\\ &&+\nabla^2\widetilde{R}_{i\alpha}^{(p-k+l)/2} \partial_{\alpha,m} P_{i\alpha,k,(p-k-l)/2} + \nabla^2\widetilde{R}_{i\alpha}^{(p-k-l)/2} \partial_{\alpha,m} P_{i\alpha,k,(p-k+l)/2} + \widetilde{R}_{i\alpha}^{(p-k+l)/2} \partial_{\alpha,m} \nabla^2 P_{i\alpha,k,(p-k-l)/2} + \widetilde{R}_{i\alpha}^{(p-k-l)/2} \partial_{\alpha,m} \nabla^2 P_{i\alpha,k,(p-k+l)/2} \\ &&+ \left. 2 \partial_{\alpha,m} \nabla \widetilde{R}_{i\alpha}^{(p-k+l)/2} \cdot \nabla P_{i\alpha,k,(p-k-l)/2} + 2 \partial_{\alpha,m} \nabla \widetilde{R}_{i\alpha}^{(p-k-l)/2} \cdot \nabla P_{i\alpha,k,(p-k+l)/2} + 2 \nabla \widetilde{R}_{i\alpha}^{(p-k+l)/2} \cdot \partial_{\alpha,m} \nabla P_{i\alpha,k,(p-k-l)/2} + 2 \nabla \widetilde{R}_{i\alpha}^{(p-k-l)/2} \cdot \partial_{\alpha,m} \nabla P_{i\alpha,k,(p-k+l)/2} \right)\\ &=& \sum_{i=1}^{N_\text{elec}}\sum_{p=2}^{N_\text{nord}} \sum_{k=0}^{p-1} \sum_{l=0}^{p-k-2\delta_{k,0}} c_{l,k,p,\alpha} \left( \partial_{\alpha,m} \nabla^2\widetilde{R}_{i\alpha}^{(p-k+l)/2} P_{i\alpha,k,(p-k-l)/2} + \partial_{\alpha,m} \nabla^2\widetilde{R}_{i\alpha}^{(p-k-l)/2} P_{i\alpha,k,(p-k+l)/2} - \partial_{i,m} \widetilde{R}_{i\alpha}^{(p-k+l)/2} \nabla^2 P_{i\alpha,k,(p-k-l)/2} - \partial_{i,m} \widetilde{R}_{i\alpha}^{(p-k-l)/2} \nabla^2 P_{i\alpha,k,(p-k+l)/2} \right.\\ &&+\nabla^2\widetilde{R}_{i\alpha}^{(p-k+l)/2} \partial_{\alpha,m} P_{i\alpha,k,(p-k-l)/2} + \nabla^2\widetilde{R}_{i\alpha}^{(p-k-l)/2} \partial_{\alpha,m} P_{i\alpha,k,(p-k+l)/2} + \widetilde{R}_{i\alpha}^{(p-k+l)/2} \partial_{\alpha,m} \nabla^2 P_{i\alpha,k,(p-k-l)/2} + \widetilde{R}_{i\alpha}^{(p-k-l)/2} \partial_{\alpha,m} \nabla^2 P_{i\alpha,k,(p-k+l)/2} \\ &&+ \left. 2 \partial_{\alpha,m} \nabla \widetilde{R}_{i\alpha}^{(p-k+l)/2} \cdot \nabla P_{i\alpha,k,(p-k-l)/2} + 2 \partial_{\alpha,m} \nabla \widetilde{R}_{i\alpha}^{(p-k-l)/2} \cdot \nabla P_{i\alpha,k,(p-k+l)/2} + 2 \nabla \widetilde{R}_{i\alpha}^{(p-k+l)/2} \cdot \partial_{\alpha,m} \nabla P_{i\alpha,k,(p-k-l)/2} + 2 \nabla \widetilde{R}_{i\alpha}^{(p-k-l)/2} \cdot \partial_{\alpha,m} \nabla P_{i\alpha,k,(p-k+l)/2} \right) \end{eqnarray} **** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_jastrow_een_l(qmckl_context context, double* const forces_jastrow_een_l, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_jastrow_een_l(qmckl_context context, double* const forces_jastrow_een_l, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_forces_jastrow_een_l(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = ctx->electron.walker.num * 3 * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_een_l", "Array too small. Expected 3*walk_num*nucl_num"); } memcpy(forces_jastrow_een_l, ctx->forces.forces_jastrow_een_l, sze*sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_een_l (context, & forces_jastrow_een_l, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (qmckl_context) , intent(in), value :: context integer(c_int64_t), intent(in), value :: size_max real(c_double), intent(out) :: forces_jastrow_een_l(size_max) end function qmckl_get_forces_jastrow_een_l end interface #+end_src # **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_een_l(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_een_l(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->jastrow_champ.cord_num > 0) { /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_jastrow_champ_c_vector_full(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_lkpm_combined_index(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if tmp_c is provided */ rc = qmckl_provide_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if dtmp_c is provided */ rc = qmckl_provide_dtmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if forces_tmp_c is provided */ rc = qmckl_provide_forces_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if forces_dtmp_c is provided */ rc = qmckl_provide_forces_dtmp_c(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if ne distance is provided */ rc = qmckl_provide_en_distance(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_forces_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; } /* Compute if necessary */ if (ctx->date > ctx->forces.forces_jastrow_een_l_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_jastrow_een_l != NULL) { rc = qmckl_free(context, ctx->forces.forces_jastrow_een_l); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_een_l", "Unable to free ctx->forces.forces_jastrow_een_l"); } ctx->forces.forces_jastrow_een_l = NULL; } } /* Allocate array */ if (ctx->forces.forces_jastrow_een_l == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->electron.walker.num * ctx->nucleus.num * sizeof(double); double* forces_jastrow_een_l = (double*) qmckl_malloc(context, mem_info); if (forces_jastrow_een_l == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_een_l", NULL); } ctx->forces.forces_jastrow_een_l = forces_jastrow_een_l; } rc = qmckl_compute_forces_jastrow_een_l(context, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.dim_c_vector, ctx->jastrow_champ.c_vector_full, ctx->jastrow_champ.lkpm_combined_index, ctx->electron.en_distance, ctx->jastrow_champ.tmp_c, ctx->jastrow_champ.dtmp_c, ctx->forces.forces_tmp_c, ctx->forces.forces_dtmp_c, ctx->forces.forces_een_n, ctx->jastrow_champ.een_rescaled_n, ctx->jastrow_champ.een_rescaled_n_gl, ctx->forces.forces_jastrow_een_l); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_jastrow_een_l_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_een_l :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_jastrow_een_l_args | Variable | Type | In/Out | Description | |------------------------+------------------------------------------------------------------------+--------+---------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~dim_c_vector~ | ~int64_t~ | in | dimension of full coefficient vector | | ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | in | full coefficient vector | | ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | in | combined indices | | ~en_distance~ | ~double[elec_num][nucl_num]~ | in | Electron-nucleus distance | | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | P matrix | | ~dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | Derivative of P matrix | | ~forces_tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | Force of P matrix | | ~forces_dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][4][nucl_num][4][elec_num]~ | in | Force of derivative of P matrix | | ~forces_een_n~ | ~double[walk_num][0:cord_num][3][nucl_num][4][elec_num]~ | in | Force of derivative of electron-nucleus rescaled factor | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled factor | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Derivative of electron-nucleus rescaled factor | | ~forces_jastrow_een_l~ | ~double[walk_num][3][nucl_num]~ | out | Force of electron-electron-nucleus Jastrow Laplacian | |------------------------+------------------------------------------------------------------------+--------+---------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_een_l( & context, walk_num, elec_num, nucl_num, & cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, & en_distance, tmp_c, dtmp_c, forces_tmp_c, forces_dtmp_c, forces_een_n, een_rescaled_n, & een_rescaled_n_gl, forces_jastrow_een_l)& result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer (c_int64_t) , intent(in),value :: walk_num, elec_num, cord_num, nucl_num, dim_c_vector integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real (c_double ) , intent(in) :: c_vector_full(nucl_num, dim_c_vector) real (c_double ) , intent(in) :: en_distance(nucl_num, elec_num) real (c_double ) , intent(in) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: dtmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_tmp_c(elec_num, 4,nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_dtmp_c(elec_num, 4, nucl_num,4,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_een_n(elec_num, 4, nucl_num, 3, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(out) :: forces_jastrow_een_l(nucl_num,3,walk_num) integer*8 :: i, a, j, l, k, m, n, nw, ii, jj double precision :: accu, accu2, cn info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 if (cord_num < 0) info = QMCKL_INVALID_ARG_5 if (info /= QMCKL_SUCCESS) return if (cord_num == 0) return forces_jastrow_een_l = 0.0d0 do nw =1, walk_num do n = 1, dim_c_vector l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = c_vector_full(a, n) if(cn == 0.d0) cycle do ii = 1, 3 do i = 1, elec_num forces_jastrow_een_l(a, ii, nw) = forces_jastrow_een_l(a, ii, nw) + (& tmp_c(i,a, m, k,nw) * forces_een_n(i,4,a,ii ,m+l,nw) & + forces_tmp_c(i,ii,a, m, k,nw) * een_rescaled_n_gl(i,4,a,m+l,nw) & + tmp_c(i,a, m+l,k,nw) * forces_een_n(i,4,a,ii, m, nw) & + forces_tmp_c(i,ii,a, m+l,k,nw) * een_rescaled_n_gl(i,4,a,m, nw) & - dtmp_c(i,4,a, m, k,nw) * een_rescaled_n_gl(i,ii,a,m+l,nw) & + forces_dtmp_c(i,4,a,ii,m, k,nw) * een_rescaled_n(i,a, m+l,nw) & - dtmp_c(i,4,a, m+l,k,nw) * een_rescaled_n_gl(i,ii,a,m, nw) & + forces_dtmp_c(i,4,a,ii,m+l,k,nw) * een_rescaled_n(i,a, m, nw) & ) * cn end do end do cn = cn + cn do ii = 1, 3 do i = 1, elec_num do jj = 1, 3 forces_jastrow_een_l(a, ii, nw) = forces_jastrow_een_l(a, ii, nw) + (& dtmp_c(i,jj,a,m, k,nw) * forces_een_n(i,jj,a,ii ,m+l,nw) + & dtmp_c(i,jj,a,m+l, k,nw) * forces_een_n(i,jj,a,ii ,m,nw) + & forces_dtmp_c(i,jj,a,ii,m, k,nw) * een_rescaled_n_gl(i,jj,a,m+l,nw) + & forces_dtmp_c(i,jj,a,ii,m+l, k,nw) * een_rescaled_n_gl(i,jj,a,m,nw) & ) * cn end do end do end do end do end do end do end function qmckl_compute_forces_jastrow_een_l #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_een_l( const qmckl_context context, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* en_distance, const double* tmp_c, const double* dtmp_c, const double* forces_tmp_c, const double* forces_dtmp_c, const double* forces_een_n, const double* een_rescaled_n, const double* een_rescaled_n_gl, double* const forces_jastrow_een_l ); #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces Jastrow een l\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double forces_jastrow_een_l[walk_num][3][nucl_num]; rc = qmckl_get_forces_jastrow_een_l(context, &forces_jastrow_een_l[0][0][0], 3*nucl_num*walk_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_een_l[walk_num][nucl_num][3][4][elec_num]; rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_factor_een_gl, &finite_difference_force_een_l[0][0][0][0][0], 4*elec_num); double finite_difference_force_een_l_sum[walk_num][nucl_num][3]; for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ finite_difference_force_een_l_sum[nw][a][k] = 0; for (int i = 0; i < elec_num; i++){ finite_difference_force_een_l_sum[nw][a][k] = finite_difference_force_een_l_sum[nw][a][k] + finite_difference_force_een_l[nw][a][k][3][i]; } } } } for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ //printf("finite_difference_force_een_l_sum[%i][%i][%i] %+3.10f \n", nw,a,k,finite_difference_force_een_l_sum[nw][a][k]); //printf("forces_jastrow_een_l [%i][%i][%i] %+3.10f\n", nw,k,a,forces_jastrow_een_l[nw][k][a]); assert(fabs(finite_difference_force_een_l_sum[nw][a][k] - forces_jastrow_een_l[nw][k][a]) < 1.e-8); } } } for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ assert(fabs(finite_difference_force_een_l_sum[nw][a][k] - forces_jastrow_een_l[nw][k][a]) < 1.e-8); } } } printf("OK\n"); #+end_src *** Force of $\delta P$ matrix Calculates the force of the $\delta P$ matrix, required for force in a single-electron move. Here, the $r^\text{new}$ and $R^\text{new}$ are the new electron and nucleus positions of the electron that is being moved. \begin{eqnarray*} \partial_{\alpha,m}\delta P_{i,\alpha,k,l} &=& \partial_{\alpha,m} \left(\sum_{j=1}^{N_\text{elec}} \left( \delta \widetilde{r}_{ij}^k \widetilde{R}_{j\alpha}^l \delta_{i,\text{num}} \right) + \widetilde{r}_{i,\text{num}}^k \delta \widetilde{R}_{\text{num},\alpha}^l + \delta \widetilde{r}_{i,\text{num}}^k \left( \widetilde{R}_{\text{num},\alpha}^l + \delta \widetilde{R}_{\text{num},\alpha}^l \right)\right) \\ &=& \partial_{\alpha,m} \left(\sum_{j=1}^{N_\text{elec}} \left( \delta \widetilde{r}_{ij}^k \widetilde{R}_{j\alpha}^l \delta_{i,\text{num}} \right) + \widetilde{r}_{i,\text{num}}^k \widetilde{R}_{\text{num},\alpha}^l + {\widetilde{r}_{i,\text{num}}^{\text{new}}}^k {\widetilde{R}_{\text{num},\alpha}^{\text{new}}}^l \right) \\ &=& \left(\sum_{j=1}^{N_\text{elec}} \left( \delta \widetilde{r}_{ij}^k \partial_{\alpha,m} \widetilde{R}_{j\alpha}^l \delta_{i,\text{num}} \right) + \widetilde{r}_{i,\text{num}}^k \partial_{\alpha,m}\widetilde{R}_{\text{num},\alpha}^l + {\widetilde{r}_{i,\text{num}}^{\text{new}}}^k \partial_{\alpha,m}{\widetilde{R}_{\text{num},\alpha}^{\text{new}}}^l \right) \\ &=& \left(\sum_{j=1}^{N_\text{elec}} \left( -\delta \widetilde{r}_{ij}^k \partial_{i,m} \widetilde{R}_{j\alpha}^l \delta_{i,\text{num}} \right) - \widetilde{r}_{i,\text{num}}^k \partial_{i,m}\widetilde{R}_{\text{num},\alpha}^l - {\widetilde{r}_{i,\text{num}}^{\text{new}}}^k \partial_{i,m}{\widetilde{R}_{\text{num},\alpha}^{\text{new}}}^l \right) \\ \end{eqnarray*} The function ~qmckl_set_single_point~ has to be called before this function to set the new electron position. **** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_jastrow_delta_p(qmckl_context context, double* const forces_delta_p, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_jastrow_delta_p(qmckl_context context, double* const forces_delta_p, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_forces_jastrow_delta_p(context); if (rc != QMCKL_SUCCESS) return rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); int64_t sze = 3 * ctx->electron.walker.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_delta_p", "Array too small."); } memcpy(forces_delta_p, ctx->forces.forces_delta_p, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src **** provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_delta_p(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_delta_p(qmckl_context context) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); qmckl_exit_code rc = qmckl_provide_een_rescaled_single_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_single_n(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_single_e(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_een_rescaled_single_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Compute if necessary */ if (ctx->single_point.date > ctx->forces.forces_delta_p_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_delta_p != NULL) { rc = qmckl_free(context, ctx->forces.forces_delta_p); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_delta_p", "Unable to free ctx->forces.forces_delta_p"); } ctx->forces.forces_delta_p = NULL; } } /* Allocate array */ if (ctx->forces.forces_delta_p == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->electron.walker.num * ctx->jastrow_champ.cord_num * (ctx->jastrow_champ.cord_num + 1) * ctx->nucleus.num * ctx->electron.num * sizeof(double); double* forces_delta_p = (double*) qmckl_malloc(context, mem_info); if (forces_delta_p == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_delta_p", NULL); } ctx->forces.forces_delta_p = forces_delta_p; } rc = qmckl_compute_forces_jastrow_delta_p(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.een_rescaled_n, ctx->jastrow_champ.een_rescaled_e, ctx->single_point.een_rescaled_single_n, ctx->single_point.een_rescaled_single_e, ctx->jastrow_champ.een_rescaled_n_gl, ctx->single_point.een_rescaled_single_n_gl, ctx->forces.forces_delta_p); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_delta_p_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_delta_p_doc :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_delta_p_args | Variable | Type | In/Out | Description | |----------------------------+---------------------------------------------------------------------+--------+--------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of single electron | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-electron rescaled distances | | ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | in | Electron-nucleus single rescaled distances | | ~een_rescaled_single_e~ | ~double[walk_num][0:cord_num][elec_num]~ | in | Electron-electron single rescaled distances | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances derivatives | | ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~ | in | Electron-nucleus single rescaled distances derivatives | | ~forces_delta_p~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][3][elec_num]~ | out | Force of $\delta P$ matrix | |----------------------------+---------------------------------------------------------------------+--------+--------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_delta_p_doc( & context, num_in, walk_num, elec_num, nucl_num, cord_num, & een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, & een_rescaled_n_gl, een_rescaled_single_n_gl, forces_delta_p) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t) , intent(in), value :: num_in, walk_num, elec_num, cord_num, nucl_num real(c_double) , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) real(c_double) , intent(in) :: een_rescaled_e(elec_num, elec_num, 0:cord_num, walk_num) real(c_double) , intent(in) :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num) real(c_double) , intent(in) :: een_rescaled_single_e(elec_num, 0:cord_num, walk_num) real(c_double) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) real(c_double) , intent(in) :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num) real(c_double) , intent(out) :: forces_delta_p(elec_num,3,nucl_num,0:cord_num, 0:cord_num-1, walk_num) double precision :: tmp1, tmp2 double precision :: delta_e(elec_num) integer*8 :: i, a, j, l, k, p, m, n, nw, num double precision :: tmp, accu integer*8 :: LDA, LDB, LDC num = num_in + 1_8 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_3 if (elec_num <= 0) info = QMCKL_INVALID_ARG_4 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5 if (cord_num < 0) info = QMCKL_INVALID_ARG_6 if (info /= QMCKL_SUCCESS) return if (cord_num == 0) then forces_delta_p = 0.d0 return endif do nw=1, walk_num do m=0, cord_num-1 do j = 1, elec_num delta_e(j) = een_rescaled_e(j,num,m,nw) - een_rescaled_single_e(j,m,nw) end do do l=1, cord_num do a = 1, nucl_num do k = 1, 3 tmp1 = een_rescaled_n_gl(num, k, a, l, nw) tmp2 = een_rescaled_single_n_gl(k,a,l,nw) accu = 0.d0 do j = 1, elec_num forces_delta_p(j,k,a,l,m,nw) = een_rescaled_e(j,num,m,nw) * tmp1 - & een_rescaled_single_e(j,m,nw) * tmp2 accu = accu + delta_e(j) * een_rescaled_n_gl(j,k,a,l,nw) end do forces_delta_p(num,k,a,l,m,nw) = forces_delta_p(num,k,a,l,m,nw) + accu end do end do end do end do end do end function qmckl_compute_forces_jastrow_delta_p_doc #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes qmckl_exit_code qmckl_compute_forces_jastrow_delta_p_hpc (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const forces_delta_p ) { if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT; if (num < 0) return QMCKL_INVALID_ARG_2; if (walk_num <= 0) return QMCKL_INVALID_ARG_3; if (elec_num <= 0) return QMCKL_INVALID_ARG_4; if (nucl_num <= 0) return QMCKL_INVALID_ARG_5; if (cord_num < 0) return QMCKL_INVALID_ARG_6; if (een_rescaled_n == NULL) return QMCKL_INVALID_ARG_7; if (een_rescaled_e == NULL) return QMCKL_INVALID_ARG_8; if (een_rescaled_single_n == NULL) return QMCKL_INVALID_ARG_9; if (een_rescaled_single_e == NULL) return QMCKL_INVALID_ARG_10; if (een_rescaled_n_gl == NULL) return QMCKL_INVALID_ARG_11; if (een_rescaled_single_n_gl == NULL) return QMCKL_INVALID_ARG_12; if (cord_num == 0) { const size_t dim = elec_num*3*nucl_num*(cord_num+1)*cord_num; #pragma omp parallel for for (int64_t nw = 0; nw < walk_num; ++nw) { memset(&forces_delta_p[dim*nw],0,dim*sizeof(double)); } return QMCKL_SUCCESS; } #pragma omp parallel for for (int64_t nw = 0; nw < walk_num; ++nw) { for (int64_t m = 0; m <= cord_num-1; ++m) { double delta_e[elec_num]; const double* een_rescaled_e_ = &een_rescaled_e[elec_num*(num+elec_num*(m+(cord_num+1)*nw))]; const double* een_rescaled_single_e_ = &een_rescaled_single_e[elec_num*(m+(cord_num+1)*nw)]; for (int64_t j = 0; j < elec_num; ++j) { delta_e[j] = een_rescaled_e_[j] - een_rescaled_single_e_[j]; } for (int64_t l = 1; l <= cord_num; ++l) { for (int64_t a = 0; a < nucl_num; ++a) { for (int64_t k = 0; k < 3; ++k) { const double* een_rescaled_n_gl_ = &een_rescaled_n_gl[elec_num*(k+4*(a+nucl_num*(l+(cord_num+1)*nw)))]; const double tmp1 = een_rescaled_n_gl_[num]; const double tmp2 = een_rescaled_single_n_gl[k+4*(a+nucl_num*(l+(cord_num+1)*nw))]; double* forces_delta_p_ = &forces_delta_p[elec_num*(k+3*(a+nucl_num*(l+(cord_num+1)*(m+cord_num*nw))))]; double accu = 0.0; #pragma omp simd reduction(+:accu) for (int64_t j = 0; j < elec_num; ++j) { forces_delta_p_[j] = een_rescaled_e_[j] * tmp1 - een_rescaled_single_e_[j] * tmp2; accu += delta_e[j] * een_rescaled_n_gl_[j]; } forces_delta_p_[num] += accu; } } } } } return QMCKL_SUCCESS; } #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_delta_p_doc (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const forces_delta_p ); qmckl_exit_code qmckl_compute_forces_jastrow_delta_p_hpc (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const forces_delta_p ); qmckl_exit_code qmckl_compute_forces_jastrow_delta_p (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const forces_delta_p ); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_delta_p (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const double* een_rescaled_n, const double* een_rescaled_e, const double* een_rescaled_single_n, const double* een_rescaled_single_e, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const forces_delta_p ) { #ifdef HAVE_HPC return qmckl_compute_forces_jastrow_delta_p_hpc #else return qmckl_compute_forces_jastrow_delta_p_doc #endif (context, num, walk_num, elec_num, nucl_num, cord_num, een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, een_rescaled_n_gl, een_rescaled_single_n_gl, forces_delta_p); } #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces delta p\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3*walk_num); assert (rc == QMCKL_SUCCESS); double forces_delta_p[walk_num][cord_num][cord_num+1][nucl_num][3][elec_num]; rc = qmckl_get_forces_jastrow_delta_p(context, &forces_delta_p[0][0][0][0][0][0], 3*nucl_num*walk_num*elec_num*(cord_num+1)*cord_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_delta_p[walk_num][cord_num][cord_num+1][nucl_num][3][elec_num]; nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (nucleus_coord == NULL) { return QMCKL_ALLOCATION_FAILED; } rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num); temp_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (temp_coord == NULL) { free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } double output_delta_p[walk_num][cord_num][cord_num+1][nucl_num][elec_num]; // Copy original coordinates for (int i = 0; i < 3 * nucl_num; i++) { temp_coord[i] = nucleus_coord[i]; } for (int64_t a = 0; a < nucl_num; a++) { for (int64_t k = 0; k < 3; k++) { for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x; // Update coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); rc = qmckl_single_touch(context); assert(rc == QMCKL_SUCCESS); // Call the provided function rc = qmckl_get_jastrow_champ_delta_p(context, &output_delta_p[0][0][0][0][0], nucl_num*walk_num*elec_num*(cord_num+1)*cord_num); assert(rc == QMCKL_SUCCESS); // Accumulate derivative using finite-difference coefficients for (int nw=0 ; nwelectron.walker.num * 3 * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_jastrow_single_een", "input array too small"); } memcpy(forces_jastrow_single_een, ctx->forces.forces_jastrow_single_een, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_forces_jastrow_single_een (context, & forces_jastrow_single_een, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: size_max real(c_double), intent(out) :: forces_jastrow_single_een(size_max) end function qmckl_get_forces_jastrow_single_een end interface #+end_src **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_jastrow_single_een(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :export none qmckl_exit_code qmckl_provide_forces_jastrow_single_een(qmckl_context context) { qmckl_exit_code rc = QMCKL_SUCCESS; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_jastrow_single_een", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->jastrow_champ.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_forces_jastrow_single_een", NULL); } if (ctx->jastrow_champ.cord_num > 0) { /* Check if en rescaled distance is provided */ rc = qmckl_provide_een_rescaled_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_een_rescaled_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance is provided */ rc = qmckl_provide_een_rescaled_single_n(context); if(rc != QMCKL_SUCCESS) return rc; /* Check if en rescaled distance derivatives is provided */ rc = qmckl_provide_een_rescaled_single_n_gl(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_forces_tmp_c(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_jastrow_champ_delta_p(context); if(rc != QMCKL_SUCCESS) return rc; rc = qmckl_provide_forces_jastrow_delta_p(context); if(rc != QMCKL_SUCCESS) return rc; } /* Compute if necessary */ if (ctx->single_point.date > ctx->forces.forces_jastrow_single_een_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { if (ctx->forces.forces_jastrow_single_een != NULL) { rc = qmckl_free(context, ctx->forces.forces_jastrow_single_een); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_jastrow_single_een", "Unable to free ctx->forces.forces_jastrow_single_een"); } ctx->forces.forces_jastrow_single_een = NULL; } } /* Allocate array */ if (ctx->forces.forces_jastrow_single_een == NULL) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->electron.walker.num * 3 * ctx->nucleus.num * sizeof(double); double* forces_jastrow_single_een = (double*) qmckl_malloc(context, mem_info); if (forces_jastrow_single_een == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_provide_forces_jastrow_single_een", NULL); } ctx->forces.forces_jastrow_single_een = forces_jastrow_single_een; } rc = qmckl_compute_forces_jastrow_single_een(context, ctx->single_point.num, ctx->electron.walker.num, ctx->electron.num, ctx->nucleus.num, ctx->jastrow_champ.cord_num, ctx->jastrow_champ.dim_c_vector, ctx->jastrow_champ.c_vector_full, ctx->jastrow_champ.lkpm_combined_index, ctx->single_point.delta_p, ctx->forces.forces_delta_p, ctx->jastrow_champ.tmp_c, ctx->forces.forces_tmp_c, ctx->jastrow_champ.een_rescaled_n, ctx->single_point.een_rescaled_single_n, ctx->jastrow_champ.een_rescaled_n_gl, ctx->single_point.een_rescaled_single_n_gl, ctx->forces.forces_jastrow_single_een); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_jastrow_single_een_date = ctx->single_point.date; } return QMCKL_SUCCESS; } #+end_src **** Compute :PROPERTIES: :Name: qmckl_compute_forces_jastrow_single_een :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_single_een_args | Variable | Type | In/Out | Description | |-----------------------------+---------------------------------------------------------------------+--------+--------------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~num~ | ~int64_t~ | in | Index of single electron | | ~walk_num~ | ~int64_t~ | in | Number of walkers | | ~elec_num~ | ~int64_t~ | in | Number of electrons | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~cord_num~ | ~int64_t~ | in | order of polynomials | | ~dim_c_vector~ | ~int64_t~ | in | dimension of full coefficient vector | | ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | in | full coefficient vector | | ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | in | combined indices | | ~delta_p~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | $\delta P$ matrix | | ~forces_delta_p~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][3][elec_num]~ | in | Force of $\delta P$matrix | | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | P matrix | | ~forces_tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | Force of P matrix | | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled distances | | ~een_rescaled_single_n~ | ~double[walk_num][0:cord_num][nucl_num]~ | in | Electron-nucleus single rescaled distances | | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Electron-nucleus rescaled distances derivatives | | ~een_rescaled_single_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4]~ | in | Electron-nucleus single rescaled distances derivatives | | ~forces_jastrow_single_een~ | ~double[walk_num][nucl_num][3]~ | out | Single electron--electron-nucleus Jastrow forces | |-----------------------------+---------------------------------------------------------------------+--------+--------------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_forces_jastrow_single_een_doc( & context, num_in, walk_num, elec_num, nucl_num, cord_num, & dim_c_vector, c_vector_full, lkpm_combined_index, & delta_p, forces_delta_p, tmp_c, forces_tmp_c, & een_rescaled_n, een_rescaled_single_n, een_rescaled_n_gl, een_rescaled_single_n_gl, forces_jastrow_single_een) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num_in integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: dim_c_vector integer (c_int64_t) , intent(in) , value :: cord_num integer(c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real(c_double) , intent(in) :: c_vector_full(nucl_num, dim_c_vector) real (c_double ) , intent(in) :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_delta_p(elec_num, 3, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_tmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(out) :: forces_jastrow_single_een(3,nucl_num,walk_num) integer(qmckl_exit_code) :: info double precision :: een_rescaled_delta_n(nucl_num, 0:cord_num), een_rescaled_delta_n_gl(3,nucl_num, 0:cord_num) integer*8 :: i, a, j, l, k, p, m, n, nw, num, kk double precision :: accu, accu2, cn integer*8 :: LDA, LDB, LDC num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_3 if (elec_num <= 0) info = QMCKL_INVALID_ARG_4 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5 if (cord_num < 0) info = QMCKL_INVALID_ARG_6 if (info /= QMCKL_SUCCESS) return forces_jastrow_single_een = 0.0d0 if (cord_num == 0) return do nw =1, walk_num een_rescaled_delta_n(:,:) = een_rescaled_single_n(:,:,nw) - een_rescaled_n(num,:,:,nw) do kk = 1,3 een_rescaled_delta_n_gl(kk,:,:) = een_rescaled_single_n_gl(kk,:,:,nw) - een_rescaled_n_gl(num,kk,:,:,nw) end do do n = 1, dim_c_vector l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) p = lkpm_combined_index(n, 3) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = c_vector_full(a, n) if(cn == 0.d0) cycle do kk = 1, 3 accu = 0.0d0 do j = 1, elec_num accu = accu - een_rescaled_n_gl(j,kk,a,m,nw) * delta_p(j,a,m+l,k,nw) + & een_rescaled_n(j,a,m,nw) * forces_delta_p(j,kk,a,m+l,k,nw) end do accu = accu - een_rescaled_delta_n_gl(kk,a,m) * (tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw)) + & een_rescaled_delta_n(a,m) * (forces_tmp_c(num,kk,a,m+l,k,nw) + forces_delta_p(num,kk,a,m+l,k,nw)) forces_jastrow_single_een(kk,a,nw) = forces_jastrow_single_een(kk,a,nw) + accu * cn end do end do end do end do end function qmckl_compute_forces_jastrow_single_een_doc #+end_src #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_forces_jastrow_single_een_hpc( & context, num_in, walk_num, elec_num, nucl_num, cord_num, & dim_c_vector, c_vector_full, lkpm_combined_index, & delta_p, forces_delta_p, tmp_c, forces_tmp_c, & een_rescaled_n, een_rescaled_single_n, een_rescaled_n_gl, een_rescaled_single_n_gl, forces_jastrow_single_een) & bind(C) result(info) use qmckl implicit none integer (qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in) , value :: num_in integer (c_int64_t) , intent(in) , value :: walk_num integer (c_int64_t) , intent(in) , value :: elec_num integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) , value :: dim_c_vector integer (c_int64_t) , intent(in) , value :: cord_num integer(c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) real(c_double) , intent(in) :: c_vector_full(nucl_num, dim_c_vector) real (c_double ) , intent(in) :: delta_p(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_delta_p(elec_num, 3, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: tmp_c(elec_num, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: forces_tmp_c(elec_num, 4, nucl_num,0:cord_num, 0:cord_num-1, walk_num) real (c_double ) , intent(in) :: een_rescaled_n(elec_num, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num, 4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(in) :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num) real (c_double ) , intent(out) :: forces_jastrow_single_een(3,nucl_num,walk_num) integer(qmckl_exit_code) :: info double precision, allocatable :: een_rescaled_delta_n(:, :), een_rescaled_delta_n_gl(:,:,:) integer*8 :: i, a, j, l, k, p, m, n, nw, num, kk double precision :: accu2, cn double precision, allocatable :: accu(:,:), tmp(:) double precision, external :: ddot integer :: elec_num4 num = num_in + 1 info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT if (walk_num <= 0) info = QMCKL_INVALID_ARG_3 if (elec_num <= 0) info = QMCKL_INVALID_ARG_4 if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5 if (cord_num < 0) info = QMCKL_INVALID_ARG_6 if (info /= QMCKL_SUCCESS) return forces_jastrow_single_een = 0.0d0 if (cord_num == 0) return allocate(een_rescaled_delta_n(nucl_num, 0:cord_num), een_rescaled_delta_n_gl(3,nucl_num, 0:cord_num), & accu(3,nucl_num), tmp(nucl_num)) elec_num4 = int(elec_num,4) do nw =1, walk_num een_rescaled_delta_n(:,:) = een_rescaled_single_n(:,:,nw) - een_rescaled_n(num,:,:,nw) een_rescaled_delta_n_gl(1:3,:,:) = een_rescaled_single_n_gl(1:3,:,:,nw) - een_rescaled_n_gl(num,1:3,:,:,nw) do n = 1, dim_c_vector l = lkpm_combined_index(n, 1) k = lkpm_combined_index(n, 2) p = lkpm_combined_index(n, 3) m = lkpm_combined_index(n, 4) do a = 1, nucl_num cn = c_vector_full(a, n) accu(1:3,a) = 0.0d0 if(cn == 0.d0) cycle tmp(a) = tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw) call dgemv('T', elec_num4, 3, -1.d0, een_rescaled_n_gl(1,1,a,m,nw), elec_num4, & delta_p(1,a,m+l,k,nw), 1, 0.d0, accu(1,a), 1) call dgemv('T', elec_num4, 3, 1.d0, forces_delta_p(1,1,a,m+l,k,nw), elec_num4, & een_rescaled_n(1,a,m,nw), 1, 1.d0, accu(1,a), 1) enddo accu(1,:) = accu(1,:) - een_rescaled_delta_n_gl(1,:,m)*tmp(:) accu(2,:) = accu(2,:) - een_rescaled_delta_n_gl(2,:,m)*tmp(:) accu(3,:) = accu(3,:) - een_rescaled_delta_n_gl(3,:,m)*tmp(:) accu(1,:) = accu(1,:) + een_rescaled_delta_n(:,m)*forces_tmp_c(num,1,:,m+l,k,nw) accu(2,:) = accu(2,:) + een_rescaled_delta_n(:,m)*forces_tmp_c(num,2,:,m+l,k,nw) accu(3,:) = accu(3,:) + een_rescaled_delta_n(:,m)*forces_tmp_c(num,3,:,m+l,k,nw) accu(1,:) = accu(1,:) + een_rescaled_delta_n(:,m)*forces_delta_p(num,1,:,m+l,k,nw) accu(2,:) = accu(2,:) + een_rescaled_delta_n(:,m)*forces_delta_p(num,2,:,m+l,k,nw) accu(3,:) = accu(3,:) + een_rescaled_delta_n(:,m)*forces_delta_p(num,3,:,m+l,k,nw) forces_jastrow_single_een(1,:,nw) = forces_jastrow_single_een(1,:,nw) + accu(1,:) * c_vector_full(:,n) forces_jastrow_single_een(2,:,nw) = forces_jastrow_single_een(2,:,nw) + accu(2,:) * c_vector_full(:,n) forces_jastrow_single_een(3,:,nw) = forces_jastrow_single_een(3,:,nw) + accu(3,:) * c_vector_full(:,n) end do end do end function qmckl_compute_forces_jastrow_single_een_hpc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_single_een_doc ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* delta_p, const double* forces_delta_p, const double* tmp_c, const double* forces_tmp_c, const double* een_rescaled_n, const double* een_rescaled_single_n, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const forces_jastrow_single_een ); qmckl_exit_code qmckl_compute_forces_jastrow_single_een ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* delta_p, const double* forces_delta_p, const double* tmp_c, const double* forces_tmp_c, const double* een_rescaled_n, const double* een_rescaled_single_n, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const forces_jastrow_single_een ); qmckl_exit_code qmckl_compute_forces_jastrow_single_een_hpc ( const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* delta_p, const double* forces_delta_p, const double* tmp_c, const double* forces_tmp_c, const double* een_rescaled_n, const double* een_rescaled_single_n, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const forces_jastrow_single_een ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_forces_jastrow_single_een (const qmckl_context context, const int64_t num, const int64_t walk_num, const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, const int64_t dim_c_vector, const double* c_vector_full, const int64_t* lkpm_combined_index, const double* delta_p, const double* forces_delta_p, const double* tmp_c, const double* forces_tmp_c, const double* een_rescaled_n, const double* een_rescaled_single_n, const double* een_rescaled_n_gl, const double* een_rescaled_single_n_gl, double* const forces_jastrow_single_een ) { #ifdef HAVE_HPC return qmckl_compute_forces_jastrow_single_een_hpc #else return qmckl_compute_forces_jastrow_single_een_doc #endif (context, num, walk_num, elec_num, nucl_num, cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, delta_p, forces_delta_p, tmp_c, forces_tmp_c, een_rescaled_n, een_rescaled_single_n, een_rescaled_n_gl, een_rescaled_single_n_gl, forces_jastrow_single_een ); } #+end_src **** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces Jastrow single een\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_single_point(context, 'N', 2, new_coords, 3*walk_num); assert (rc == QMCKL_SUCCESS); double forces_jastrow_single_een[walk_num][nucl_num][3]; rc = qmckl_get_forces_jastrow_single_een(context, &forces_jastrow_single_een[0][0][0], 3*nucl_num*walk_num); assert(rc == QMCKL_SUCCESS); double finite_difference_force_single_een[walk_num][nucl_num][3]; rc = qmckl_finite_difference_deriv_n(context, delta_x, &qmckl_get_jastrow_champ_single_een, &(finite_difference_force_single_een[0][0][0]), 1); for (int nw = 0; nw < walk_num; nw++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ //printf("%.10f\t", finite_difference_force_single_een[nw][a][k]); //printf("%.10f\n", forces_jastrow_single_een[nw][a][k]); assert(fabs(finite_difference_force_single_een[nw][a][k] - forces_jastrow_single_een[nw][a][k]) < 1.e-8); } } } printf("OK\n"); #+end_src * Forces of the orbitals ** Reset test for orbitals :noexport: #+begin_src c :tangle (eval c_test) rc = qmckl_context_destroy(context); assert (rc == QMCKL_SUCCESS); context = qmckl_context_create(); nucl_num = chbrclf_nucl_num; nucl_charge = chbrclf_charge; nucl_coord = &(chbrclf_nucl_coord[0][0]); rc = qmckl_set_nucleus_num (context, nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_nucleus_charge(context, nucl_charge, nucl_num); assert(rc == QMCKL_SUCCESS); assert(qmckl_nucleus_provided(context)); const int64_t shell_num = chbrclf_shell_num; const int64_t prim_num = chbrclf_prim_num; const int64_t ao_num = chbrclf_ao_num; const int64_t * nucleus_index = &(chbrclf_basis_nucleus_index[0]); const int64_t * nucleus_shell_num = &(chbrclf_basis_nucleus_shell_num[0]); const int32_t * shell_ang_mom = &(chbrclf_basis_shell_ang_mom[0]); const int64_t * shell_prim_num = &(chbrclf_basis_shell_prim_num[0]); const int64_t * shell_prim_index = &(chbrclf_basis_shell_prim_index[0]); const double * shell_factor = &(chbrclf_basis_shell_factor[0]); const double * exponent = &(chbrclf_basis_exponent[0]); const double * coefficient = &(chbrclf_basis_coefficient[0]); const double * prim_factor = &(chbrclf_basis_prim_factor[0]); const double * ao_factor = &(chbrclf_basis_ao_factor[0]); const char typ = 'G'; assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_type (context, typ); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_num (context, shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_prim_num (context, prim_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, nucl_num); assert(rc == QMCKL_ALREADY_SET); rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, nucl_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, shell_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_exponent (context, exponent, prim_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_coefficient (context, coefficient, prim_num); assert(rc == QMCKL_SUCCESS); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, prim_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_ao_basis_ao_num(context, ao_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, ao_num); assert(rc == QMCKL_SUCCESS); assert(qmckl_ao_basis_provided(context)); int64_t shell_num_test ; int64_t prim_num_test ; int64_t ao_num_test ; int64_t * nucleus_index_test ; int64_t * nucleus_shell_num_test; int32_t * shell_ang_mom_test ; int64_t * shell_prim_num_test ; int64_t * shell_prim_index_test ; double * shell_factor_test ; double * exponent_test ; double * coefficient_test ; double * prim_factor_test ; double * ao_factor_test ; char typ_test ; rc = qmckl_get_ao_basis_type (context, &typ_test); assert (rc == QMCKL_SUCCESS); assert(typ == typ_test); rc = qmckl_get_ao_basis_shell_num (context, &shell_num_test); assert (rc == QMCKL_SUCCESS); assert(shell_num == shell_num_test); rc = qmckl_get_ao_basis_prim_num (context, &prim_num_test); assert (rc == QMCKL_SUCCESS); assert(prim_num == prim_num_test); nucleus_index_test = (int64_t*) malloc (nucl_num * sizeof(int64_t)); rc = qmckl_get_ao_basis_nucleus_index (context, nucleus_index_test, nucl_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < nucl_num ; ++i) { assert(nucleus_index_test[i] == nucleus_index[i]); } free(nucleus_index_test); nucleus_shell_num_test = (int64_t*) malloc ( nucl_num * sizeof(int64_t)); rc = qmckl_get_ao_basis_nucleus_shell_num (context, nucleus_shell_num_test, nucl_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < nucl_num ; ++i) { assert(nucleus_shell_num_test[i] == nucleus_shell_num[i]); } free(nucleus_shell_num_test); shell_ang_mom_test = (int32_t*) malloc ( shell_num * sizeof(int32_t)); rc = qmckl_get_ao_basis_shell_ang_mom (context, shell_ang_mom_test, shell_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < shell_num ; ++i) { assert(shell_ang_mom_test[i] == shell_ang_mom[i]); } free(shell_ang_mom_test); shell_factor_test = (double*) malloc ( shell_num * sizeof(double)); rc = qmckl_get_ao_basis_shell_factor (context, shell_factor_test, shell_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < shell_num ; ++i) { assert(shell_factor_test[i] == shell_factor[i]); } free(shell_factor_test); shell_prim_num_test = (int64_t*) malloc ( shell_num * sizeof(int64_t)); rc = qmckl_get_ao_basis_shell_prim_num (context, shell_prim_num_test, shell_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < shell_num ; ++i) { assert(shell_prim_num_test[i] == shell_prim_num[i]); } free(shell_prim_num_test); shell_prim_index_test = (int64_t*) malloc ( shell_num * sizeof(int64_t)); rc = qmckl_get_ao_basis_shell_prim_index (context, shell_prim_index_test, shell_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < shell_num ; ++i) { assert(shell_prim_index_test[i] == shell_prim_index[i]); } free(shell_prim_index_test); exponent_test = (double*) malloc ( prim_num * sizeof(double)); rc = qmckl_get_ao_basis_exponent(context, exponent_test, prim_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < prim_num ; ++i) { assert(exponent_test[i] == exponent[i]); } free(exponent_test); coefficient_test = (double*) malloc ( prim_num * sizeof(double)); rc = qmckl_get_ao_basis_coefficient(context, coefficient_test, prim_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < prim_num ; ++i) { assert(coefficient_test[i] == coefficient[i]); } free(coefficient_test); prim_factor_test = (double*) malloc ( prim_num * sizeof(double)); rc = qmckl_get_ao_basis_prim_factor (context, prim_factor_test, prim_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < prim_num ; ++i) { assert(prim_factor_test[i] == prim_factor[i]); } free(prim_factor_test); rc = qmckl_get_ao_basis_ao_num(context, &ao_num_test); assert(ao_num == ao_num_test); ao_factor_test = (double*) malloc ( ao_num * sizeof(double)); rc = qmckl_get_ao_basis_ao_factor (context, ao_factor_test, ao_num); assert (rc == QMCKL_SUCCESS); for (int64_t i=0 ; i < ao_num ; ++i) { assert(ao_factor_test[i] == ao_factor[i]); } free(ao_factor_test); #define walk_num 1 // chbrclf_walk_num #define elec_num chbrclf_elec_num #define prim_num chbrclf_prim_num elec_up_num = chbrclf_elec_up_num; elec_dn_num = chbrclf_elec_dn_num; elec_coord = &(chbrclf_elec_coord[0][0][0]); rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); assert (rc == QMCKL_SUCCESS); assert(qmckl_electron_provided(context)); int64_t point_num = elec_num; rc = qmckl_set_point(context, 'N', point_num, elec_coord, point_num*3); assert(rc == QMCKL_SUCCESS); int64_t mo_num = chbrclf_mo_num; rc = qmckl_set_mo_basis_mo_num(context, mo_num); assert (rc == QMCKL_SUCCESS); const double * mo_coefficient = &(chbrclf_mo_coef[0]); rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient, chbrclf_mo_num*chbrclf_ao_num); assert (rc == QMCKL_SUCCESS); #+end_src ** Force of AO value Computes the forces on the values of the atomic orbitals (AO). These are equal to the gradients of the AOs multiplied with a minus sign, and summed over different indices. *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_ao_value(qmckl_context context, double* const forces_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_forces_ao_value(qmckl_context context, double* const forces_ao_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_forces_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->nucleus.num * 3 * ctx->point.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_ao_value", "Array too small. Expected walk_num*nucl_num*point_num*3"); } memcpy(forces_ao_value, ctx->forces.forces_ao_value, sze*sizeof(double)); return QMCKL_SUCCESS; } #+end_src *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_ao_value(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_ao_value(qmckl_context context) { qmckl_exit_code rc = QMCKL_SUCCESS; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_ao_value", NULL); } 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_provide_forces_ao_value", NULL); } /* Compute if necessary */ if (ctx->point.date > ctx->forces.forces_ao_value_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->ao_basis.ao_num * 3 * ctx->nucleus.num * ctx->point.num * sizeof(double); if (mem_info.size > ctx->forces.forces_ao_value_maxsize) { if (ctx->forces.forces_ao_value != NULL) { rc = qmckl_free(context, ctx->forces.forces_ao_value); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_forces_ao_value", "Unable to free ctx->forces.forces_ao_value"); } ctx->forces.forces_ao_value = NULL; } } /* Allocate array */ if (ctx->forces.forces_ao_value == NULL) { double* forces_ao_value = (double*) qmckl_malloc(context, mem_info); ctx->forces.forces_ao_value_maxsize = mem_info.size; if (forces_ao_value == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_forces_ao_value", NULL); } ctx->forces.forces_ao_value = forces_ao_value; } rc = qmckl_provide_ao_basis_ao_vgl(context); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, rc, "qmckl_provide_ao_basis_ao_vgl", NULL); } //memset(ctx->forces.forces_ao_value, 0, mem_info.size); rc = qmckl_compute_forces_ao_value_doc(context, ctx->ao_basis.ao_num, ctx->ao_basis.shell_num, ctx->point.num, ctx->nucleus.num, ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.shell_ang_mom, ctx->ao_basis.ao_factor, ctx->ao_basis.ao_vgl, ctx->forces.forces_ao_value); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_ao_value_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_forces_ao_value :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_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 | | ~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 | | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | | ~ao_factor~ | ~double[ao_num]~ | in | Normalization factor of the AOs | | ~ao_vgl~ | ~double[point_num][5][shell_num]~ | in | Value, gradients and Laplacian of the shells | | ~forces_ao_value~ | ~double[nucl_num][3][point_num][ao_num]~ | out | Forces of the AOs | |---------------------+------------------------------------------+--------+----------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes function qmckl_compute_forces_ao_value_doc(context, & ao_num, shell_num, point_num, nucl_num, & nucleus_index, nucleus_shell_num, & shell_ang_mom, ao_factor, ao_vgl, forces_ao_value) & bind(C) result(info) use qmckl_constants use qmckl, only : qmckl_ao_polynomial_vgl, qmckl_get_numprec_precision implicit none integer (qmckl_context), intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: ao_num integer (c_int64_t) , intent(in) , value :: shell_num integer (c_int64_t) , intent(in) , value :: point_num integer (c_int64_t) , intent(in) , value :: nucl_num integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num) real (c_double ) , intent(in) :: ao_factor(ao_num) real (c_double ) , intent(in) :: ao_vgl(ao_num,5,point_num) real (c_double ) , intent(out) :: forces_ao_value(ao_num,point_num,3,nucl_num) integer(qmckl_exit_code) :: info integer :: l, il, k integer*8 :: ipoint, inucl integer*8 :: ishell_start, ishell_end, ishell integer :: lstart(0:20) double precision :: x, y, z, r2 double precision :: cutoff integer , allocatable :: ao_index(:) allocate(ao_index(ao_num)) !forces_ao_value = 0.d0 ! 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 do ipoint = 1, point_num do inucl=1,nucl_num ! 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 forces_ao_value(k,ipoint,1,inucl) = -ao_vgl(k,2,ipoint) forces_ao_value(k,ipoint,2,inucl) = -ao_vgl(k,3,ipoint) forces_ao_value(k,ipoint,3,inucl) = -ao_vgl(k,4,ipoint) k = k+1 end do end do end do end do deallocate(ao_index) end function qmckl_compute_forces_ao_value_doc #+end_src #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_ao_value_doc( const qmckl_context context, const int64_t ao_num, const int64_t shell_num, const int64_t point_num, const int64_t nucl_num, const int64_t* nucleus_index, const int64_t* nucleus_shell_num, const int32_t* shell_ang_mom, const double* ao_factor, const double* ao_vgl, double* const forces_ao_value ); #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces AO value\n"); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double * forces_ao_value = (double*) malloc(nucl_num * 3 *point_num * ao_num *sizeof(double)); rc = qmckl_get_forces_ao_value(context, &forces_ao_value[0], 3*nucl_num*ao_num*point_num); assert(rc == QMCKL_SUCCESS); double * finite_difference_force_ao_value = (double*) malloc(3 *nucl_num * point_num * ao_num *sizeof(double)); nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (nucleus_coord == NULL) { return QMCKL_ALLOCATION_FAILED; } rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num); temp_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (temp_coord == NULL) { free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } double * ao_output = (double*) malloc(point_num*ao_num*sizeof(double)); // Copy original coordinates for (int i = 0; i < 3 * nucl_num; i++) { temp_coord[i] = nucleus_coord[i]; } for (int64_t a = 0; a < nucl_num; a++) { for (int64_t k = 0; k < 3; k++) { for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x; // Update coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); // Call the provided function rc = qmckl_get_ao_basis_ao_value(context,&ao_output[0], point_num*ao_num); assert(rc == QMCKL_SUCCESS); // Accumulate derivative using finite-difference coefficients for (int i = 0; i < point_num; i++) { for (int j = 0; j < ao_num; j++) { if (m == -4) { finite_difference_force_ao_value[k*ao_num*point_num*nucl_num + a*ao_num*point_num + i*ao_num + j] = 0.0; } finite_difference_force_ao_value[k*ao_num*point_num*nucl_num + a*ao_num*point_num + i*ao_num + j] += coef[m + 4] * ao_output[i*ao_num + j]/delta_x; } } } temp_coord[k+a*3] = nucleus_coord[k+3*a]; } } // Reset coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); free(nucleus_coord); free(temp_coord); free(ao_output); for (int j = 0; j < ao_num; j++){ for (int i = 0; i < point_num; i++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ //printf("k=%i a=%i i=%i j=%i\n", k, a, i, j); //printf("%.10f\t", finite_difference_force_ao_value[k*ao_num*point_num*nucl_num + a*ao_num*point_num + i*ao_num + j]); //printf("%.10f\n", forces_ao_value[a*3*ao_num*point_num + k*ao_num*point_num + i*ao_num + j]); assert(fabs(finite_difference_force_ao_value[k*ao_num*point_num*nucl_num + a*ao_num*point_num + i*ao_num + j] - forces_ao_value[a*3*ao_num*point_num + k*ao_num*point_num + i*ao_num + j]) < 1.e-9); } } } } free(forces_ao_value); free(finite_difference_force_ao_value); printf("OK\n"); #+end_src ** Force of MO value Calculates the force on the molecular orbitals (MO) values. These are calculated by taking a linear combination of AOs. *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_mo_value(qmckl_context context, double* const forces_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_forces_mo_value(qmckl_context context, double* const forces_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_forces_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 * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_mo_value", "input array too small"); } memcpy(forces_mo_value, ctx->forces.forces_mo_value, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_forces_mo_value (context, & forces_mo_value, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real(c_double), intent(out) :: forces_mo_value(*) integer (c_int64_t) , intent(in) , value :: size_max end function qmckl_get_forces_mo_value end interface #+end_src #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_mo_value_inplace (qmckl_context context, double* const forces_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_forces_mo_value_inplace (qmckl_context context, double* const forces_mo_value, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_get_forces_mo_value", NULL); } qmckl_exit_code rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); const int64_t sze = ctx->point.num * 3 * ctx->mo_basis.mo_num * ctx->nucleus.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->forces.forces_mo_value; ctx->forces.forces_mo_value = forces_mo_value; rc = qmckl_provide_forces_mo_value(context); if (rc != QMCKL_SUCCESS) return rc; ctx->forces.forces_mo_value = old_array; return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_forces_mo_value_inplace (context, & forces_mo_value, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real(c_double), intent(out) :: forces_mo_value(*) integer (c_int64_t) , intent(in) , value :: size_max end function qmckl_get_forces_mo_value_inplace end interface #+end_src *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_mo_value(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :export none qmckl_exit_code qmckl_provide_forces_mo_value(qmckl_context context) { qmckl_exit_code rc = QMCKL_SUCCESS; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_mo_value", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->mo_basis.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_mo_basis_mo_vgl", NULL); } /* Compute if necessary */ if (ctx->point.date > ctx->forces.forces_mo_value_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->mo_basis.mo_num * ctx->point.num * ctx->nucleus.num * sizeof(double); if (mem_info.size > ctx->forces.forces_mo_value_maxsize) { if (ctx->forces.forces_mo_value != NULL) { rc = qmckl_free(context, ctx->forces.forces_mo_value); assert (rc == QMCKL_SUCCESS); ctx->forces.forces_mo_value = NULL; } } /* Allocate array */ if (ctx->forces.forces_mo_value == NULL) { double* forces_mo_value = (double*) qmckl_malloc(context, mem_info); ctx->forces.forces_mo_value_maxsize = mem_info.size; if (forces_mo_value == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_forces_mo_value", NULL); } ctx->forces.forces_mo_value = forces_mo_value; } rc = qmckl_provide_forces_ao_value(context); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_forces_ao_value", NULL); } rc = qmckl_compute_forces_mo_value_doc(context, ctx->nucleus.num, ctx->ao_basis.ao_num, ctx->mo_basis.mo_num, ctx->point.num, ctx->ao_basis.shell_num, ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.shell_ang_mom, ctx->mo_basis.coefficient_t, ctx->forces.forces_ao_value, ctx->forces.forces_mo_value); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_mo_value_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_forces_mo_value :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_forces_mo_value_args | Variable | Type | In/Out | Description | |---------------------+------------------------------------------+--------+-------------------------------------------------| | ~context~ | ~qmckl_context~ | in | Global state | | ~nucl_num~ | ~int64_t~ | in | Number of AOs | | ~ao_num~ | ~int64_t~ | in | Number of AOs | | ~mo_num~ | ~int64_t~ | in | Number of MOs | | ~point_num~ | ~int64_t~ | in | Number of points | | ~shell_num~ | ~int64_t~ | in | Number of shells | | ~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 | | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | | ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix | | ~forces_ao_value~ | ~double[nucl_num][3][point_num][ao_num]~ | in | Forces of AOs | | ~forces_mo_value~ | ~double[nucl_num][3][point_num][mo_num]~ | out | Forces of MOs | |---------------------+------------------------------------------+--------+-------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer(qmckl_exit_code) function qmckl_compute_forces_mo_value_doc(context, & nucl_num,ao_num, mo_num, point_num, shell_num, nucleus_index, nucleus_shell_num, & shell_ang_mom, coefficient_t, forces_ao_value, forces_mo_value) & result(info) bind(C) use, intrinsic :: iso_c_binding use qmckl implicit none integer(qmckl_context), intent(in) :: context integer(c_int64_t) , intent(in), value ::nucl_num, ao_num, mo_num, point_num, shell_num real(c_double) , intent(in) :: forces_ao_value(ao_num,point_num,3,nucl_num) real(c_double) , intent(in) :: coefficient_t(mo_num,ao_num) integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num) real(c_double) , intent(out) :: forces_mo_value(mo_num,point_num,3,nucl_num) integer*8 :: i,j,a, m, ishell_start, ishell_end, il, ishell integer :: l, k double precision :: c1, c2, c3 integer, allocatable :: ao_index(:) integer :: lstart(0:20) allocate(ao_index(ao_num)) do l=0,20 lstart(l) = l*(l+1)*(l+2)/6 +1 end do k=1 do a=1,nucl_num ishell_start = nucleus_index(a) + 1 ishell_end = nucleus_index(a) + nucleus_shell_num(a) 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 forces_mo_value = 0.0d0 do a=1,nucl_num ishell_start = nucleus_index(a) + 1 ishell_end = nucleus_index(a) + nucleus_shell_num(a) do j=1,point_num do ishell = ishell_start, ishell_end k = ao_index(ishell) l = shell_ang_mom(ishell) do il = lstart(l), lstart(l+1)-1 c1 = forces_ao_value(k,j,1,a) c2 = forces_ao_value(k,j,2,a) c3 = forces_ao_value(k,j,3,a) if (c1 == 0.d0 .and. c2 == 0.d0 .and. c3 == 0.d0) cycle do i=1,mo_num forces_mo_value(i,j,1,a) = forces_mo_value(i,j,1,a) + coefficient_t(i,k) * c1 forces_mo_value(i,j,2,a) = forces_mo_value(i,j,2,a) + coefficient_t(i,k) * c2 forces_mo_value(i,j,3,a) = forces_mo_value(i,j,3,a) + coefficient_t(i,k) * c3 end do k = k + 1 end do end do end do end do deallocate(ao_index) end function qmckl_compute_forces_mo_value_doc #+end_src #+begin_src c :tangle (eval h_private_func) :comments org qmckl_exit_code qmckl_compute_forces_mo_value_doc ( const qmckl_context context, const int64_t nucl_num, const int64_t ao_num, const int64_t mo_num, const int64_t point_num, const int64_t shell_num, const int64_t* nucleus_index, const int64_t* nucleus_shell_num, const int32_t* shell_ang_mom, const double* coefficient_t, const double* forces_ao_value, double* const forces_mo_value ); #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces MO value\n"); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double * forces_mo_value = (double*) malloc(nucl_num * point_num* 3 * mo_num *sizeof(double)); //double forces_mo_value[nucl_num][3][point_num][mo_num]; rc = qmckl_get_forces_mo_value(context, &forces_mo_value[0], 3*nucl_num*mo_num*point_num); assert(rc == QMCKL_SUCCESS); //double finite_difference_force_mo_value[3][nucl_num][point_num][mo_num]; double * finite_difference_force_mo_value = (double*) malloc(nucl_num* 3 * point_num * mo_num * sizeof(double)); nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (nucleus_coord == NULL) { return QMCKL_ALLOCATION_FAILED; } rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num); temp_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (temp_coord == NULL) { free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } double* mo_output = (double*) malloc(point_num * mo_num * sizeof(double)); if (mo_output == NULL) { free(temp_coord); free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } // Copy original coordinates for (int i = 0; i < 3 * nucl_num; i++) { temp_coord[i] = nucleus_coord[i]; } for (int64_t a = 0; a < nucl_num; a++) { for (int64_t k = 0; k < 3; k++) { for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x; // Update coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); // Call the provided function rc = qmckl_get_mo_basis_mo_value(context,&mo_output[0], point_num*mo_num); assert(rc == QMCKL_SUCCESS); // Accumulate derivative using finite-difference coefficients for (int i = 0; i < point_num; i++) { for (int j = 0; j < mo_num; j++) { if (m == -4) { finite_difference_force_mo_value[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j]= 0.0; } finite_difference_force_mo_value[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] += coef[m + 4] * mo_output[i*mo_num+ j]/delta_x; } } } temp_coord[k+a*3] = nucleus_coord[k+3*a]; } } // Reset coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); free(nucleus_coord); free(temp_coord); free(mo_output); for (int j = 0; j < mo_num; j++){ for (int i = 0; i < point_num; i++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ //printf("k=%i a=%i i=%i j=%i\n", k, a, i, j); //printf("%.10f\t", finite_difference_force_mo_value[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j]); //printf("%.10f\n", forces_mo_value[a*3*mo_num*point_num + k*mo_num*point_num + i*mo_num + j]); assert(fabs(finite_difference_force_mo_value[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] - forces_mo_value[a*3*mo_num*point_num + k*mo_num*point_num + i*mo_num + j]) < 1.e-9); } } } } free(forces_mo_value); free(finite_difference_force_mo_value); printf("OK\n"); #+end_src ** Force of MO gradient Calculates the forces of the gradients of the MOs. These are calculated by taking a linear combination of the required components of the Hessian of the AOs. *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_mo_g(qmckl_context context, double* const forces_mo_g, const int64_t size_max); qmckl_exit_code qmckl_get_forces_mo_g_inplace(qmckl_context context, double* const forces_mo_g, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_mo_g(qmckl_context context, double* const forces_mo_g, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_get_forces_mo_g", NULL); } qmckl_exit_code rc; rc = qmckl_provide_forces_mo_g(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 * 3 * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_mo_g", "input array too small"); } memcpy(forces_mo_g, ctx->forces.forces_mo_g, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_mo_g_inplace (qmckl_context context, double* const forces_mo_g, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_get_forces_mo_g", NULL); } qmckl_exit_code rc; qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); const int64_t sze = ctx->point.num * 3 * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_mo_g", "input array too small"); } rc = qmckl_context_touch(context); if (rc != QMCKL_SUCCESS) return rc; double* old_array = ctx->forces.forces_mo_g; ctx->forces.forces_mo_g = forces_mo_g; rc = qmckl_provide_forces_mo_g(context); if (rc != QMCKL_SUCCESS) return rc; ctx->forces.forces_mo_g = old_array; return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_forces_mo_g (context, & forces_mo_g, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real(c_double), intent(out) :: forces_mo_g(*) integer (c_int64_t) , intent(in) , value :: size_max end function qmckl_get_forces_mo_g end interface interface integer(qmckl_exit_code) function qmckl_get_forces_mo_g_inplace (context, & forces_mo_g, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real(c_double), intent(out) :: forces_mo_g(*) integer (c_int64_t) , intent(in) , value :: size_max end function qmckl_get_forces_mo_g_inplace end interface #+end_src *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_mo_g(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :export none qmckl_exit_code qmckl_provide_forces_mo_g(qmckl_context context) { qmckl_exit_code rc = QMCKL_SUCCESS; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_mo_g", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->mo_basis.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_forces_mo_g", NULL); } /* Compute if necessary */ if (ctx->point.date > ctx->forces.forces_mo_g_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * 3 * ctx->mo_basis.mo_num * ctx->point.num * ctx->nucleus.num * sizeof(double); if (ctx->forces.forces_mo_g != NULL) { qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero; rc = qmckl_get_malloc_info(context, ctx->forces.forces_mo_g, &mem_info_test); /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the memory was not allocated with qmckl_malloc */ if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) { rc = qmckl_free(context, ctx->forces.forces_mo_g); assert (rc == QMCKL_SUCCESS); ctx->forces.forces_mo_g = NULL; } } /* Allocate array */ if (ctx->forces.forces_mo_g == NULL) { double* forces_mo_g = (double*) qmckl_malloc(context, mem_info); if (forces_mo_g == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_forces_mo_g", NULL); } ctx->forces.forces_mo_g = forces_mo_g; } rc = qmckl_provide_ao_basis_ao_hessian(context); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_ao_basis_ao_hessian", NULL); } rc = qmckl_compute_forces_mo_g(context, ctx->ao_basis.ao_num, ctx->mo_basis.mo_num, ctx->point.num, ctx->nucleus.num, ctx->ao_basis.shell_num, ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.shell_ang_mom, ctx->mo_basis.coefficient_t, ctx->ao_basis.ao_hessian, ctx->forces.forces_mo_g); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_mo_g_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_forces_mo_g :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_compute_forces_mo_g_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 | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~shell_num~ | ~int64_t~ | in | Number of shells | | ~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 | | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | | ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix | | ~ao_hessian~ | ~double[nucl_num][3][point_num][4][ao_num]~ | in | Hessian of AOs | | ~forces_mo_g~ | ~double[nucl_num][3][point_num][3][mo_num]~ | out | Forces on gradients of MOs | |---------------------+---------------------------------------------+--------+-------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_forces_mo_g_doc(context, & ao_num, mo_num, point_num, nucl_num, & shell_num, nucleus_index, nucleus_shell_num, shell_ang_mom, & coefficient_t, ao_hessian, forces_mo_g) & bind(C) result(info) use qmckl implicit none integer(qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in), value :: nucl_num, ao_num, mo_num, point_num, shell_num integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num) real (c_double ) , intent(in) :: coefficient_t(mo_num,ao_num) real (c_double ) , intent(in) :: ao_hessian(ao_num,4,point_num,3,nucl_num) real (c_double ) , intent(out) :: forces_mo_g(mo_num,3,point_num,3,nucl_num) integer*8 :: i,j, m, n,a double precision :: c1 double precision, allocatable :: tmp(:,:,:,:) allocate(tmp(ao_num,3,point_num,3)) do a=1,nucl_num tmp(:,1:3,:,:) = ao_hessian(:,1:3,:,:,a) info = qmckl_dgemm(context, 'N', 'N', mo_num, 3*point_num*3, ao_num, & -1.d0, coefficient_t, mo_num, & tmp, ao_num, 0.d0, forces_mo_g(1,1,1,1,a), mo_num) if (info /= 0) return end do deallocate(tmp) end function qmckl_compute_forces_mo_g_doc #+end_src #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_forces_mo_g_hpc(context, & ao_num, mo_num, point_num, nucl_num, & shell_num, nucleus_index, nucleus_shell_num, shell_ang_mom, & coefficient_t, ao_hessian, forces_mo_g) & bind(C) result(info) use qmckl implicit none integer(qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in), value :: nucl_num, ao_num, mo_num, point_num, shell_num integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num) real (c_double ) , intent(in) :: coefficient_t(mo_num,ao_num) real (c_double ) , intent(in) :: ao_hessian(ao_num,4,point_num,3,nucl_num) real (c_double ) , intent(out) :: forces_mo_g(mo_num,3,point_num,3,nucl_num) integer*8 :: i,j, m, n,a double precision :: c1 integer :: l, il, k integer*8 :: ishell_start, ishell_end, ishell integer :: lstart(0:20) integer , allocatable :: ao_index(:) allocate(ao_index(ao_num)) do l=0,20 lstart(l) = l*(l+1)*(l+2)/6 +1 end do k=1 do a=1,nucl_num ishell_start = nucleus_index(a) + 1 ishell_end = nucleus_index(a) + nucleus_shell_num(a) 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 forces_mo_g = 0.d0 do a=1,nucl_num ishell_start = nucleus_index(a) + 1 ishell_end = nucleus_index(a) + nucleus_shell_num(a) do n = 1, 3 do j=1,point_num do m = 1, 3 do ishell = ishell_start, ishell_end l = shell_ang_mom(ishell) il = lstart(l+1)-lstart(l) do k = ao_index(ishell), ao_index(ishell) + il - 1 c1 = ao_hessian(k, m, j, n, a) if (c1 == 0.d0) cycle do i=1,mo_num forces_mo_g(i, m, j, n, a) = forces_mo_g(i, m, j, n, a) - coefficient_t(i,k) * c1 end do end do end do end do end do end do end do deallocate(ao_index) end function qmckl_compute_forces_mo_g_hpc #+end_src #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_compute_forces_mo_g ( const qmckl_context context, const int64_t ao_num, const int64_t mo_num, const int64_t point_num, const int64_t nucl_num, const int64_t shell_num, const int64_t* nucleus_index, const int64_t* nucleus_shell_num, const int32_t* shell_ang_mom, const double* coefficient_t, const double* ao_hessian, double* const forces_mo_g ); qmckl_exit_code qmckl_compute_forces_mo_g_doc ( const qmckl_context context, const int64_t ao_num, const int64_t mo_num, const int64_t point_num, const int64_t nucl_num, const int64_t shell_num, const int64_t* nucleus_index, const int64_t* nucleus_shell_num, const int32_t* shell_ang_mom, const double* coefficient_t, const double* ao_hessian, double* const forces_mo_g ); qmckl_exit_code qmckl_compute_forces_mo_g_hpc ( const qmckl_context context, const int64_t ao_num, const int64_t mo_num, const int64_t point_num, const int64_t nucl_num, const int64_t shell_num, const int64_t* nucleus_index, const int64_t* nucleus_shell_num, const int32_t* shell_ang_mom, const double* coefficient_t, const double* ao_hessian, double* const forces_mo_g ); #+end_src #+begin_src c :tangle (eval c) :comments org :exports none qmckl_exit_code qmckl_compute_forces_mo_g ( const qmckl_context context, const int64_t ao_num, const int64_t mo_num, const int64_t point_num, const int64_t nucl_num, const int64_t shell_num, const int64_t* nucleus_index, const int64_t* nucleus_shell_num, const int32_t* shell_ang_mom, const double* coefficient_t, const double* ao_hessian, double* const forces_mo_g ) { #ifdef HAVE_HPC return qmckl_compute_forces_mo_g_doc #else return qmckl_compute_forces_mo_g_doc #endif (context, ao_num, mo_num, point_num, nucl_num, shell_num, nucleus_index, nucleus_shell_num, shell_ang_mom, coefficient_t, ao_hessian, forces_mo_g ); } #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces MO gradient\n"); qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double * forces_mo_g_doc = (double*) malloc(nucl_num* 3 * point_num* 3 * mo_num *sizeof(double)); rc = qmckl_get_forces_mo_g(context, &forces_mo_g_doc[0], 3*3*nucl_num*mo_num*point_num); rc = qmckl_compute_forces_mo_g_doc(context, ctx->ao_basis.ao_num, ctx->mo_basis.mo_num, ctx->point.num, ctx->nucleus.num, ctx->ao_basis.shell_num, ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.shell_ang_mom, ctx->mo_basis.coefficient_t, ctx->ao_basis.ao_hessian, &forces_mo_g_doc[0]); assert(rc == QMCKL_SUCCESS); double * forces_mo_g = (double*) malloc(nucl_num* 3 * point_num* 3 * mo_num *sizeof(double)); rc = qmckl_compute_forces_mo_g_hpc(context, ctx->ao_basis.ao_num, ctx->mo_basis.mo_num, ctx->point.num, ctx->nucleus.num, ctx->ao_basis.shell_num, ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.shell_ang_mom, ctx->mo_basis.coefficient_t, ctx->ao_basis.ao_hessian, &forces_mo_g[0]); assert(rc == QMCKL_SUCCESS); double * finite_difference_force_mo_g = (double*) malloc(nucl_num* 3 * point_num* 3 * mo_num * sizeof(double)); nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (nucleus_coord == NULL) { return QMCKL_ALLOCATION_FAILED; } rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num); temp_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (temp_coord == NULL) { free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } mo_output = (double*) malloc(5 * point_num * mo_num * sizeof(double)); if (mo_output == NULL) { free(temp_coord); free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } // Copy original coordinates for (int i = 0; i < 3 * nucl_num; i++) { temp_coord[i] = nucleus_coord[i]; } for (int64_t a = 0; a < nucl_num; a++) { for (int64_t k = 0; k < 3; k++) { for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x; // Update coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); // Call the provided function rc = qmckl_get_mo_basis_mo_vgl(context,&mo_output[0], 5*point_num*mo_num); assert(rc == QMCKL_SUCCESS); // Accumulate derivative using finite-difference coefficients for (int i = 0; i < point_num; i++) { for (int n = 0; n < 3; n++){ for (int j = 0; j < mo_num; j++) { if (m == -4) { finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j] = 0.0; } finite_difference_force_mo_g[k*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + i*3*mo_num + n*mo_num + j] += coef[m + 4] * mo_output[i*mo_num*5 + (n+1)*mo_num + j]/delta_x; } } } } temp_coord[k+a*3] = nucleus_coord[k+3*a]; } } // Reset coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); free(nucleus_coord); free(temp_coord); free(mo_output); for (int a = 0; a < nucl_num; a++) { for (int n = 0; n < 3; n++){ for (int j = 0; j < point_num; j++){ for (int m = 0; m < 3; m++){ for (int i = 0; i < mo_num; i++){ printf("a=%i n=%i j=%i m=%i i=%i\n", a, n, j, m, i); printf("hpc %.10f\n", forces_mo_g[a*9*mo_num*point_num + n*3*mo_num*point_num + j*3*mo_num + m*mo_num + i]); printf("doc %.10f\n", forces_mo_g_doc[a*9*mo_num*point_num + n*3*mo_num*point_num + j*3*mo_num + m*mo_num + i]); printf("fd %.10f\n", finite_difference_force_mo_g[n*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + j*3*mo_num + m*mo_num + i]); assert(fabs(forces_mo_g[a*9*mo_num*point_num + n*3*mo_num*point_num + j*3*mo_num + m*mo_num + i] - forces_mo_g_doc[a*9*mo_num*point_num + n*3*mo_num*point_num + j*3*mo_num + m*mo_num + i] ) < 1.e-9); assert(fabs(finite_difference_force_mo_g[n*3*mo_num*point_num*nucl_num + a*3*mo_num*point_num + j*3*mo_num + m*mo_num + i] - forces_mo_g[a*9*mo_num*point_num + n*3*mo_num*point_num + j*3*mo_num + m*mo_num + i]) < 1.e-8); } } } } } free(forces_mo_g_doc); free(forces_mo_g); free(finite_difference_force_mo_g); printf("OK\n"); #+end_src ** Force of MO Laplacian Computes the forces on the Laplacian of the MOs. They are calculated by taking a linear combination of the ~ao_hessian[:,:,:,3,:]~ components of the AO Hessian. These store the derivatives of the Laplacian of the AO. *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes qmckl_exit_code qmckl_get_forces_mo_l(qmckl_context context, double* const forces_mo_l, const int64_t size_max); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_forces_mo_l(qmckl_context context, double* const forces_mo_l, const int64_t size_max) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_NULL_CONTEXT; } qmckl_exit_code rc; rc = qmckl_provide_forces_mo_l(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 * 3 * ctx->mo_basis.mo_num * ctx->nucleus.num; if (size_max < sze) { return qmckl_failwith( context, QMCKL_INVALID_ARG_3, "qmckl_get_forces_mo_l", "input array too small"); } memcpy(forces_mo_l, ctx->forces.forces_mo_l, sze * sizeof(double)); return QMCKL_SUCCESS; } #+end_src #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(qmckl_exit_code) function qmckl_get_forces_mo_l (context, & forces_mo_l, size_max) bind(C) use, intrinsic :: iso_c_binding import implicit none integer (c_int64_t) , intent(in) , value :: context real(c_double), intent(out) :: forces_mo_l(*) integer (c_int64_t) , intent(in) , value :: size_max end function qmckl_get_forces_mo_l end interface #+end_src *** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_forces_mo_l(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :export none qmckl_exit_code qmckl_provide_forces_mo_l(qmckl_context context) { qmckl_exit_code rc = QMCKL_SUCCESS; if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return qmckl_failwith( context, QMCKL_INVALID_CONTEXT, "qmckl_provide_forces_mo_l", NULL); } qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); if (!ctx->mo_basis.provided) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_provide_mo_basis_mo_vgl", NULL); } /* Compute if necessary */ if (ctx->point.date > ctx->forces.forces_mo_l_date) { qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = 3 * ctx->mo_basis.mo_num * ctx->point.num * ctx->nucleus.num * sizeof(double); if (ctx->forces.forces_mo_l != NULL) { qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero; rc = qmckl_get_malloc_info(context, ctx->forces.forces_mo_l, &mem_info_test); /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the memory was not allocated with qmckl_malloc */ if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) { rc = qmckl_free(context, ctx->forces.forces_mo_l); assert (rc == QMCKL_SUCCESS); ctx->forces.forces_mo_l= NULL; } } /* Allocate array */ if (ctx->forces.forces_mo_l == NULL) { double* forces_mo_l = (double*) qmckl_malloc(context, mem_info); if (forces_mo_l == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, "qmckl_forces_mo_l", NULL); } ctx->forces.forces_mo_l = forces_mo_l; } rc = qmckl_provide_ao_basis_ao_hessian(context); if (rc != QMCKL_SUCCESS) { return qmckl_failwith( context, QMCKL_NOT_PROVIDED, "qmckl_ao_basis_ao_hessian", NULL); } rc = qmckl_compute_forces_mo_l_doc(context, ctx->ao_basis.ao_num, ctx->mo_basis.mo_num, ctx->point.num, ctx->nucleus.num, ctx->ao_basis.shell_num, ctx->ao_basis.nucleus_index, ctx->ao_basis.nucleus_shell_num, ctx->ao_basis.shell_ang_mom, ctx->mo_basis.coefficient_t, ctx->ao_basis.ao_hessian, ctx->forces.forces_mo_l); if (rc != QMCKL_SUCCESS) { return rc; } ctx->forces.forces_mo_l_date = ctx->date; } return QMCKL_SUCCESS; } #+end_src *** Compute :PROPERTIES: :Name: qmckl_compute_forces_mo_l :CRetType: qmckl_exit_code :FRetType: qmckl_exit_code :END: #+NAME: qmckl_compute_forces_mo_l_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 | | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | | ~shell_num~ | ~int64_t~ | in | Number of shells | | ~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 | | ~shell_ang_mom~ | ~int32_t[shell_num]~ | in | Angular momentum of each shell | | ~coefficient_t~ | ~double[mo_num][ao_num]~ | in | Transpose of the AO to MO transformation matrix | | ~ao_hessian~ | ~double[nucl_num][3][point_num][4][ao_num]~ | in | Hessian of AOs | | ~forces_mo_l~ | ~double[nucl_num][3][point_num][mo_num]~ | out | Forces on MO Laplacian | |---------------------+---------------------------------------------+--------+-------------------------------------------------| #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_forces_mo_l_doc(context, & ao_num, mo_num, point_num, nucl_num, & shell_num, nucleus_index, nucleus_shell_num, shell_ang_mom, & coefficient_t, ao_hessian, forces_mo_l) & bind(C) result(info) use qmckl implicit none integer(qmckl_context), intent(in), value :: context integer (c_int64_t) , intent(in), value :: nucl_num, ao_num, mo_num, point_num, shell_num integer (c_int64_t) , intent(in) :: nucleus_index(nucl_num) integer (c_int64_t) , intent(in) :: nucleus_shell_num(nucl_num) integer (c_int32_t) , intent(in) :: shell_ang_mom(shell_num) real (c_double ) , intent(in) :: coefficient_t(mo_num,ao_num) real (c_double ) , intent(in) :: ao_hessian(ao_num,4,point_num,3,nucl_num) real (c_double ) , intent(out) :: forces_mo_l(mo_num,point_num,3,nucl_num) integer*8 :: i,j, m, n,a, il, ishell, ishell_start, ishell_end integer :: lstart(0:20) integer :: l, k double precision :: c1 integer , allocatable :: ao_index(:) allocate(ao_index(ao_num)) do l=0,20 lstart(l) = l*(l+1)*(l+2)/6 +1 end do k=1 do a=1,nucl_num ishell_start = nucleus_index(a) + 1 ishell_end = nucleus_index(a) + nucleus_shell_num(a) 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 forces_mo_l = 0.d0 info = QMCKL_SUCCESS do a=1, nucl_num ishell_start = nucleus_index(a) + 1 ishell_end = nucleus_index(a) + nucleus_shell_num(a) do j=1,point_num do n = 1, 3 do ishell = ishell_start, ishell_end k = ao_index(ishell) l = shell_ang_mom(ishell) do il = lstart(l), lstart(l+1)-1 c1 = ao_hessian(k, 4, j, n, a) if (c1 == 0.d0) then k = k + 1 cycle end if do i=1,mo_num forces_mo_l(i, j, n, a) = forces_mo_l(i, j, n, a) - coefficient_t(i,k) * c1 end do k = k+1 end do end do end do end do end do deallocate(ao_index) !do a=1,nucl_num ! do m = 1, 3 ! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, & ! -1.d0, coefficient_t, mo_num, & ! ao_hessian(:, 4, :, m, a), ao_num, & ! 1.d0, forces_mo_l(:, :, m, a), mo_num) ! end do !end do end function qmckl_compute_forces_mo_l_doc #+end_src #+RESULTS: #+begin_src c :tangle (eval h_func) :comments org qmckl_exit_code qmckl_compute_forces_mo_l_doc ( const qmckl_context context, const int64_t ao_num, const int64_t mo_num, const int64_t point_num, const int64_t nucl_num, const int64_t shell_num, const int64_t* nucleus_index, const int64_t* nucleus_shell_num, const int32_t* shell_ang_mom, const double* coefficient_t, const double* ao_hessian, double* const forces_mo_l ); #+end_src *** Test :noexport: #+begin_src c :tangle (eval c_test) printf("Forces MO laplacian\n"); rc = qmckl_set_nucleus_coord(context, 'T', &(nucl_coord[0]), 3*nucl_num); assert(rc == QMCKL_SUCCESS); double * forces_mo_l = (double*) malloc(nucl_num * point_num* 3 * mo_num *sizeof(double)); rc = qmckl_get_forces_mo_l(context, &forces_mo_l[0], 3*nucl_num*mo_num*point_num); assert(rc == QMCKL_SUCCESS); double * finite_difference_force_mo_l= (double*) malloc(nucl_num* 3 * point_num * mo_num * sizeof(double)); nucleus_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (nucleus_coord == NULL) { return QMCKL_ALLOCATION_FAILED; } rc = qmckl_get_nucleus_coord(context, 'N', nucleus_coord, 3*nucl_num); temp_coord = (double*) malloc(3 * nucl_num * sizeof(double)); if (temp_coord == NULL) { free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } mo_output = (double*) malloc(5 * point_num * mo_num * sizeof(double)); if (mo_output == NULL) { free(temp_coord); free(nucleus_coord); return QMCKL_ALLOCATION_FAILED; } // Copy original coordinates for (int i = 0; i < 3 * nucl_num; i++) { temp_coord[i] = nucleus_coord[i]; } for (int64_t a = 0; a < nucl_num; a++) { for (int64_t k = 0; k < 3; k++) { for (int64_t m = -4; m <= 4; m++) { // Apply finite difference displacement temp_coord[k+a*3] = nucleus_coord[k+3*a] + (double) m * delta_x; // Update coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); // Call the provided function rc = qmckl_get_mo_basis_mo_vgl(context,&mo_output[0], 5*point_num*mo_num); assert(rc == QMCKL_SUCCESS); // Accumulate derivative using finite-difference coefficients for (int i = 0; i < point_num; i++) { for (int j = 0; j < mo_num; j++) { if (m == -4) { finite_difference_force_mo_l[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] = 0.0; } finite_difference_force_mo_l[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] += coef[m + 4] * mo_output[i*mo_num*5 + 4*mo_num + j]/delta_x; } } } temp_coord[k+a*3] = nucleus_coord[k+3*a]; } } // Reset coordinates in the context rc = qmckl_set_nucleus_coord(context, 'N', temp_coord, 3*nucl_num); assert(rc == QMCKL_SUCCESS); rc = qmckl_context_touch(context); assert(rc == QMCKL_SUCCESS); free(nucleus_coord); free(temp_coord); free(mo_output); for (int j = 0; j < mo_num; j++){ for (int i = 0; i < point_num; i++){ for (int a = 0; a < nucl_num; a++) { for (int k = 0; k < 3; k++){ //printf("k=%i a=%i i=%i j=%i\n", k, a, i, j); //printf("%.10f\t", finite_difference_force_mo_l[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j]); //printf("%.10f\n", forces_mo_l[a*3*mo_num*point_num + k*mo_num*point_num + i*mo_num + j]); assert(fabs(finite_difference_force_mo_l[k*mo_num*point_num*nucl_num + a*mo_num*point_num + i*mo_num + j] - forces_mo_l[a*3*mo_num*point_num + k*mo_num*point_num + i*mo_num + j]) < 1.e-8); } } } } free(forces_mo_l); free(finite_difference_force_mo_l); printf("OK\n"); #+end_src * End of files :noexport: #+begin_src c :tangle (eval h_private_type) #endif #+end_src #+begin_src c :tangle (eval h_private_func) #endif #+end_src ** Test #+begin_src c :tangle (eval c_test) rc = qmckl_context_destroy(context); assert (rc == QMCKL_SUCCESS); return 0; } #+end_src