From c9decf482fc6b68a9ba68dea91b6875ed8782e04 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Mon, 5 Jul 2021 22:58:04 +0530 Subject: [PATCH] Working on factor_ee_deriv_e. #22 --- org/qmckl_jastrow.org | 408 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 404 insertions(+), 4 deletions(-) diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index df09b9c..52e8d90 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -1297,12 +1297,10 @@ assert(fabs(asymp_jasb[1]-0.31567342786262853) < 1.e-12); #+end_src - ** Electron-electron component \(f_{ee}\) - Calculate the asymptotic component ~asymp_jasb~ to be substracted from the final - electron-electron jastrow factor \(f_{ee}\). The asymptotic componenet is calculated - via the ~bord_vector~ and the electron-electron rescale factor ~rescale_factor_kappa~. + Calculate the electron-electron jastrow component ~factor_ee~ using the ~asymp_jasb~ + componenet and the electron-electron rescaled distances ~ee_distance_rescaled~. \[ f_{ee} = \sum_{i,jjastrow.factor_ee_deriv_e, ctx->electron.walk_num*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_factor_ee_deriv_e(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_factor_ee_deriv_e(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* const) context; + assert (ctx != NULL); + + /* Check if ee rescaled distance is provided */ + rc = qmckl_provide_ee_distance_rescaled(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow.factor_ee_date) { + + /* Allocate array */ + if (ctx->jastrow.factor_ee_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walk_num * sizeof(double); + double* factor_ee_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (factor_ee_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_factor_ee_deriv_e", + NULL); + } + ctx->jastrow.factor_ee_deriv_e = factor_ee_deriv_e; + } + + qmckl_exit_code rc = + qmckl_compute_factor_ee_deriv_e(context, + ctx->electron.walk_num, + ctx->electron.num, + ctx->electron.up_num, + ctx->jastrow.bord_num, + ctx->jastrow.bord_vector, + ctx->electron.ee_distance_rescaled, + ctx->electron.ee_distance_rescaled_deriv_e, + ctx->jastrow.asymp_jasb, + ctx->jastrow.factor_ee_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.factor_ee_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_factor_ee_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_ee_deriv_e_args + | qmckl_context | context | in | Global state | + | int64_t | walk_num | in | Number of walkers | + | int64_t | elec_num | in | Number of electrons | + | int64_t | up_num | in | Number of alpha electrons | + | int64_t | bord_num | in | Number of coefficients | + | double | bord_vector[bord_num + 1] | in | List of coefficients | + | double | ee_distance_rescaled[walk_num][elec_num][elec_num] | in | Electron-electron distances | + | double | ee_distance_rescaled_deriv_e[4][walk_num][elec_num][elec_num] | in | Electron-electron distances | + | double | asymp_jasb[2] | in | Electron-electron distances | + | double | factor_ee_deriv_e[walk_num] | out | Electron-electron distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_factor_ee_deriv_e_f(context, walk_num, elec_num, up_num, bord_num, & + bord_vector, ee_distance_rescaled, ee_distance_rescaled_deriv_e, & + asymp_jasb, factor_ee_deriv_e) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: walk_num, elec_num, bord_num, up_num + double precision , intent(in) :: bord_vector(bord_num) + double precision , intent(in) :: ee_distance_rescaled(walk_num, elec_num, elec_num) + double precision , intent(in) :: ee_distance_rescaled_deriv_e(walk_num, 4, elec_num, elec_num) + double precision , intent(in) :: asymp_jasb(2) + double precision , intent(out) :: factor_ee_deriv_e(walk_num, 4, elec_num) + + integer*8 :: i, j, p, ipar, nw, ii + double precision :: x, spin_fact, y + double precision :: den, invden, invden2, invden3, xinv + double precision :: lap1, lap2, lap3, third + double precision, dimension(3) :: pow_ser_g + double precision, dimension(4) :: dx + + 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 (bord_num <= 0) then + info = QMCKL_INVALID_ARG_4 + return + endif + + factor_ee_deriv_e = 0.0d0 + third = 1.0d0 / 3.0d0 + + do nw =1, walk_num + do j = 1, elec_num + do i = 1, elec_num + x = ee_distance_rescaled(nw,i,j) + pow_ser_g = 0.0d0 + spin_fact = 1.0d0 + den = 1.0d0 + bord_vector(2) * ee_distance_rescaled(nw, i, j) + invden = 1.0d0 / den + invden2 = invden * invden + invden3 = invden2 * invden + xinv = 1.0d0 / (ee_distance_rescaled(nw, i, j) + 1.0d0-18) + ipar = 1 + + do ii = 1, 4 + dx(ii) = ee_distance_rescaled_deriv_e(nw, ii, j, i) + end do + + if((i .LE. up_num .AND. j .LE. up_num ) .OR. & + (i .GT. up_num .AND. j .GT. up_num)) then + spin_fact = 0.5d0 + endif + + lap1 = 0.0d0 + lap2 = 0.0d0 + lap3 = 0.0d0 + do ii = 1, 3 + x = ee_distance_rescaled(nw, i, j) + do p = 2, bord_num + y = p * bord_vector(p + 1) * x + pow_ser_g(ii) = pow_ser_g(ii) + y * dx(ii) + lap1 = lap1 + (p - 1) * y * xinv * dx(ii) * dx(ii) + lap2 = lap2 + y + x = x * ee_distance_rescaled(nw, i, j) + end do + + lap3 = lap3 - 2.0d0 * bord_vector(2) * dx(ii) * dx(ii) + + factor_ee_deriv_e(nw, ii, j) = factor_ee_deriv_e(nw, ii, j) + spin_fact * bord_vector(1) * & + dx(ii) * invden + pow_ser_g(ii) + end do + + ii = 4 + lap2 = lap2 * dx(ii) * third + lap3 = lap3 + den * dx(ii) + lap3 = lap3 + spin_fact * bord_vector(1) * invden3 + factor_ee_deriv_e(nw, ii, j) = factor_ee_deriv_e(nw, ii, j) + lap1 + lap2 + lap3 + + end do + end do + end do + +end function qmckl_compute_factor_ee_deriv_e_f + #+end_src + + #+CALL: generate_c_header(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_compute_factor_ee_deriv_e ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t up_num, + const int64_t bord_num, + const double* bord_vector, + const double* ee_distance_rescaled, + const double* ee_distance_rescaled_deriv_e, + const double* asymp_jasb, + double* const factor_ee_deriv_e ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_factor_ee_deriv_e_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_factor_ee_deriv_e & + (context, & + walk_num, & + elec_num, & + up_num, & + bord_num, & + bord_vector, & + ee_distance_rescaled, & + ee_distance_rescaled_deriv_e, & + asymp_jasb, & + factor_ee_deriv_e) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: walk_num + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: up_num + integer (c_int64_t) , intent(in) , value :: bord_num + real (c_double ) , intent(in) :: bord_vector(bord_num + 1) + real (c_double ) , intent(in) :: ee_distance_rescaled(elec_num,elec_num,walk_num) + real (c_double ) , intent(in) :: ee_distance_rescaled_deriv_e(elec_num,elec_num,walk_num,4) + real (c_double ) , intent(in) :: asymp_jasb(2) + real (c_double ) , intent(out) :: factor_ee_deriv_e(walk_num) + + integer(c_int32_t), external :: qmckl_compute_factor_ee_deriv_e_f + info = qmckl_compute_factor_ee_deriv_e_f & + (context, & + walk_num, & + elec_num, & + up_num, & + bord_num, & + bord_vector, & + ee_distance_rescaled, & + ee_distance_rescaled_deriv_e, & + asymp_jasb, & + factor_ee_deriv_e) + + end function qmckl_compute_factor_ee_deriv_e + #+end_src + +*** Test + #+begin_src python :results output :exports none :noweb yes +import numpy as np + +<> + +<> + +kappa = 1.0 + +elec_coord = np.array(elec_coord)[0] +elec_dist = np.zeros(shape=(elec_num, elec_num),dtype=float) +for i in range(elec_num): + for j in range(elec_num): + elec_dist[i, j] = np.linalg.norm(np.abs(elec_coord[i] - elec_coord[j])) + +elec_dist_deriv_e = np.zeros(shape=(4,elec_num, elec_num),dtype=float) +for j in range(elec_num): + for i in range(elec_num): + rij_inv = 1.0 / elec_dist[i, j] + for ii in range(3): + elec_dist_deriv_e[ii, i, j] = (elec_coord[i][ii] - elec_coord[j][ii]) * rij_inv + elec_dist_deriv_e[3, i, j] = 2.0 * rij_inv + elec_dist_deriv_e[:, j, j] = 0.0 + +ee_distance_rescaled_deriv_e = np.zeros(shape=(4,elec_num,elec_num),dtype=float) +for i in range(elec_num): + for j in range(elec_num): + f = 1.0 - kappa * ee_distance_rescaled[i][j] + for ii in range(4): + ee_distance_rescaled_deriv_e[ii][i][j] = elec_dist_deriv_e[ii][i][j] + ee_distance_rescaled_deriv_e[3][i][j] = ee_distance_rescaled_deriv_e[3][i][j] + \ + (-kappa * ee_distance_rescaled_deriv_e[0][i][j]**2) + \ + (-kappa * ee_distance_rescaled_deriv_e[1][i][j]**2) + \ + (-kappa * ee_distance_rescaled_deriv_e[2][i][j]**2) + for ii in range(4): + ee_distance_rescaled_deriv_e[ii][i][j] = ee_distance_rescaled_deriv_e[ii][i][j] * f + + +third = 1.0 / 3.0 +factor_ee_deriv_e = np.zeros(shape=(4,elec_num),dtype=float) +dx = np.zeros(shape=(4),dtype=float) +pow_ser_g = np.zeros(shape=(4),dtype=float) +for i in range(elec_num): + for j in range(elec_num): + x = ee_distance_rescaled[i][j] + pow_ser_g = np.zeros(shape=(4),dtype=float) + spin_fact = 1.0 + den = 1.0 + bord_vector[1] * ee_distance_rescaled[i][j] + invden = 1.0 / den + invden2 = invden * invden + invden3 = invden2 * invden + xinv = 1.0 / (ee_distance_rescaled[i][j] + 1.0-18) + ipar = 1 + + for ii in range(4): + dx[ii] = ee_distance_rescaled_deriv_e[ii][j][i] + + if((i <= up_num and j <= up_num ) or \ + (i > up_num and j > up_num)): + spin_fact = 0.5 + + lap1 = 0.0 + lap2 = 0.0 + lap3 = 0.0 + for ii in range(3): + x = ee_distance_rescaled[i][j] + for p in range(1,bord_num): + y = p * bord_vector[p + 1] * x + pow_ser_g[ii] = pow_ser_g[ii] + y * dx[ii] + lap1 = lap1 + (p - 1) * y * xinv * dx[ii] * dx[ii] + lap2 = lap2 + y + x = x * ee_distance_rescaled[i][j] + + lap3 = lap3 - 2.0 * bord_vector[1] * dx[ii] * dx[ii] + + factor_ee_deriv_e[ii][j] = factor_ee_deriv_e[ii][j] + spin_fact * bord_vector[0] * \ + dx[ii] * invden + pow_ser_g[ii] + + ii = 3 + lap2 = lap2 * dx[ii] * third + lap3 = lap3 + den * dx[ii] + lap3 = lap3 + spin_fact * bord_vector[1] * invden3 + factor_ee_deriv_e[ii][j] = factor_ee_deriv_e[ii][j] + lap1 + lap2 + lap3 + +print("factor_ee_deriv_e:",factor_ee_deriv_e) + + + #+end_src + + #+RESULTS: + #+begin_example + asym_one : 0.43340325572525706 + asymp_jasb[0] : 0.5323750557252571 + asymp_jasb[1] : 0.31567342786262853 + factor_ee_deriv_e: [[-0.98933561 -1.90335162 1.55808557 -0.4175529 0.61492377 -0.39831866 + -0.39927944 0.65028576 -0.3615353 0.18529135] + [ 1.63379825 -0.20386723 -0.68681515 -0.12852739 -0.26752663 -0.38482884 + -0.08399636 0.0271138 0.47700861 -0.15444407] + [-1.50784905 -0.77857073 -2.65602338 0.47954785 2.47647536 -0.73979109 + -1.45880891 -0.57274462 0.20895041 0.22564652] + [ 3.33386195 2.70852681 -4.1135207 7.24473889 -1.68083113 6.11612348 + 3.74311012 0.3229593 7.21394604 2.7578531 ]] + #+end_example + + + #+begin_src c :tangle (eval c_test) +/* Check if Jastrow is properly initialized */ +assert(qmckl_jastrow_provided(context)); + +double factor_ee_deriv_e[walk_num]; +rc = qmckl_get_jastrow_factor_ee_deriv_e(context, factor_ee_deriv_e); + +// calculate factor_ee_deriv_e +assert(fabs(factor_ee_deriv_e[0]+4.282760865958113) < 1.e-12); + + #+end_src + * End of files :noexport: