diff --git a/org/qmckl_distance.org b/org/qmckl_distance.org index bccf454..082209b 100644 --- a/org/qmckl_distance.org +++ b/org/qmckl_distance.org @@ -1105,6 +1105,360 @@ end function qmckl_distance_rescaled_f #+end_src *** Test :noexport: +* Rescaled Distance Derivatives + +** ~qmckl_distance_rescaled_deriv_e~ + :PROPERTIES: + :Name: qmckl_distance_rescaled_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + ~qmckl_distance_rescaled_deriv_e~ computes the matrix of the gradient and laplacian of the + rescaled distance with respect to the electron coordinates. The derivative is a rank 3 tensor. + The first dimension has a dimension of 4 of which the first three coordinates + contains the gradient vector and the last index is the laplacian. + + + \[ + C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa + \] + + Here the gradient is defined as follows: + + \[ + \nabla (C_{ij}(\mathbf{r}_{ee})) = \left(\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta x},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta y},\frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta z} \right) + \] + and the laplacian is defined as follows: + + \[ + \triangle (C_{ij}(r_{ee})) = \frac{\delta^2}{\delta x^2} + \frac{\delta^2}{\delta y^2} + \frac{\delta^2}{\delta z^2} + \] + + Using the above three formulae, the expression for the gradient and laplacian is + as follows: + + \[ + \frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta x} = \frac{|(x_i - x_j)|}{r_{ij}} (1 - \kappa R_{ij}) + \] + + \[ + \frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta y} = \frac{|(y_i - y_j)|}{r_{ij}} (1 - \kappa R_{ij}) + \] + + \[ + \frac{\delta C_{ij}(\mathbf{r}_{ee})}{\delta z} = \frac{|(z_i - z_j)|}{r_{ij}} (1 - \kappa R_{ij}) + \] + + \[ + \Delta(C_{ij}(r_{ee}) = \left[ \frac{2}{r_{ij}} - \kappa \right] (1-\kappa R_{ij}) + \] + + If the input array is normal (~'N'~), the xyz coordinates are in + the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran. + + #+NAME: qmckl_distance_rescaled_deriv_e_args + | qmckl_context | context | in | Global state | + | char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed | + | char | transb | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed | + | int64_t | m | in | Number of points in the first set | + | int64_t | n | in | Number of points in the second set | + | double | A[][lda] | in | Array containing the $m \times 3$ matrix $A$ | + | int64_t | lda | in | Leading dimension of array ~A~ | + | double | B[][ldb] | in | Array containing the $n \times 3$ matrix $B$ | + | int64_t | ldb | in | Leading dimension of array ~B~ | + | double | C[4][n][ldc] | out | Array containing the $4 \times m \times n$ matrix $C$ | + | int64_t | ldc | in | Leading dimension of array ~C~ | + | double | rescale_factor_kappa | in | Factor for calculating rescaled distances derivatives | + +*** Requirements + + - ~context~ is not ~QMCKL_NULL_CONTEXT~ + - ~m > 0~ + - ~n > 0~ + - ~lda >= 3~ if ~transa == 'N'~ + - ~lda >= m~ if ~transa == 'T'~ + - ~ldb >= 3~ if ~transb == 'N'~ + - ~ldb >= n~ if ~transb == 'T'~ + - ~ldc >= m~ + - ~A~ is allocated with at least $3 \times m \times 8$ bytes + - ~B~ is allocated with at least $3 \times n \times 8$ bytes + - ~C~ is allocated with at least $4 \times m \times n \times 8$ bytes + +*** C header + + #+CALL: generate_c_header(table=qmckl_distance_rescaled_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_distance_rescaled_deriv_e ( + const qmckl_context context, + const char transa, + const char transb, + const int64_t m, + const int64_t n, + const double* A, + const int64_t lda, + const double* B, + const int64_t ldb, + double* const C, + const int64_t ldc, + const double rescale_factor_kappa); + #+end_src + +*** Source + #+begin_src f90 :tangle (eval f) +integer function qmckl_distance_rescaled_deriv_e_f(context, transa, transb, m, n, & + A, LDA, B, LDB, C, LDC, rescale_factor_kappa) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + character , intent(in) :: transa, transb + integer*8 , intent(in) :: m, n + integer*8 , intent(in) :: lda + real*8 , intent(in) :: A(lda,*) + integer*8 , intent(in) :: ldb + real*8 , intent(in) :: B(ldb,*) + integer*8 , intent(in) :: ldc + real*8 , intent(out) :: C(4,ldc,*) + real*8 , intent(in) :: rescale_factor_kappa + + integer*8 :: i,j + real*8 :: x, y, z, dist, dist_inv + real*8 :: rescale_factor_kappa_inv, rij + integer :: transab + + rescale_factor_kappa_inv = 1.0d0/rescale_factor_kappa; + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (m <= 0_8) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (n <= 0_8) then + info = QMCKL_INVALID_ARG_5 + return + endif + + if (transa == 'N' .or. transa == 'n') then + transab = 0 + else if (transa == 'T' .or. transa == 't') then + transab = 1 + else + transab = -100 + endif + + if (transb == 'N' .or. transb == 'n') then + continue + else if (transa == 'T' .or. transa == 't') then + transab = transab + 2 + else + transab = -100 + endif + + ! check for LDA + if (transab < 0) then + info = QMCKL_INVALID_ARG_1 + return + endif + + if (iand(transab,1) == 0 .and. LDA < 3) then + info = QMCKL_INVALID_ARG_7 + return + endif + + if (iand(transab,1) == 1 .and. LDA < m) then + info = QMCKL_INVALID_ARG_7 + return + endif + + if (iand(transab,2) == 0 .and. LDA < 3) then + info = QMCKL_INVALID_ARG_7 + return + endif + + if (iand(transab,2) == 2 .and. LDA < m) then + info = QMCKL_INVALID_ARG_7 + return + endif + + ! check for LDB + if (iand(transab,1) == 0 .and. LDB < 3) then + info = QMCKL_INVALID_ARG_9 + return + endif + + if (iand(transab,1) == 1 .and. LDB < n) then + info = QMCKL_INVALID_ARG_9 + return + endif + + if (iand(transab,2) == 0 .and. LDB < 3) then + info = QMCKL_INVALID_ARG_9 + return + endif + + if (iand(transab,2) == 2 .and. LDB < n) then + info = QMCKL_INVALID_ARG_9 + return + endif + + ! check for LDC + if (LDC < m) then + info = QMCKL_INVALID_ARG_11 + return + endif + + + select case (transab) + + case(0) + + do j=1,n + do i=1,m + x = A(1,i) - B(1,j) + y = A(2,i) - B(2,j) + z = A(3,i) - B(3,j) + dist = dsqrt(x*x + y*y + z*z) + dist_inv = 1.0d0/dist + rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv + C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij) + end do + end do + + case(1) + + do j=1,n + do i=1,m + x = A(i,1) - B(1,j) + y = A(i,2) - B(2,j) + z = A(i,3) - B(3,j) + dist = dsqrt(x*x + y*y + z*z) + dist_inv = 1.0d0/dist + rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv + C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij) + end do + end do + + case(2) + + do j=1,n + do i=1,m + x = A(1,i) - B(j,1) + y = A(2,i) - B(j,2) + z = A(3,i) - B(j,3) + dist = dsqrt(x*x + y*y + z*z) + dist_inv = 1.0d0/dist + rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv + C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij) + end do + end do + + case(3) + + do j=1,n + do i=1,m + x = A(i,1) - B(j,1) + y = A(i,2) - B(j,2) + z = A(i,3) - B(j,3) + dist = dsqrt(x*x + y*y + z*z) + dist_inv = 1.0d0/dist + rij = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv + C(1,i,j) = x * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(2,i,j) = y * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(3,i,j) = z * dist_inv * ( 1.0d0 - rescale_factor_kappa_inv * rij) + C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa_inv) * ( 1.0d0 - rescale_factor_kappa_inv * rij) + end do + end do + + end select + +end function qmckl_distance_rescaled_deriv_e_f + #+end_src + +*** Performance + + This function is more efficient when ~A~ and ~B~ are transposed. + +** C interface :noexport: + + #+CALL: generate_c_interface(table=qmckl_distance_rescaled_deriv_e_args,rettyp=get_value("FRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_distance_rescaled_deriv_e & + (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + character , intent(in) , value :: transa + character , intent(in) , value :: transb + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(out) :: C(4,ldc,n) + integer (c_int64_t) , intent(in) , value :: ldc + real (c_double ) , intent(in) , value :: rescale_factor_kappa + + integer(c_int32_t), external :: qmckl_distance_rescaled_deriv_e_f + info = qmckl_distance_rescaled_deriv_e_f & + (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) + + end function qmckl_distance_rescaled_deriv_e + #+end_src + + #+CALL: generate_f_interface(table=qmckl_distance_rescaled_deriv_e_args,rettyp=get_value("FRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_distance_rescaled_deriv_e & + (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, rescale_factor_kappa) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + character , intent(in) , value :: transa + character , intent(in) , value :: transb + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(out) :: C(4,ldc,n) + integer (c_int64_t) , intent(in) , value :: ldc + real (c_double ) , intent(in) , value :: rescale_factor_kappa + + end function qmckl_distance_rescaled_deriv_e + end interface + #+end_src + * End of files :noexport: #+begin_src c :comments link :tangle (eval c_test) diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 9cdd100..742ba56 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -63,25 +63,30 @@ int main() { The following data stored in the context: - | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | - | ~num~ | ~int64_t~ | Total number of electrons | - | ~up_num~ | ~int64_t~ | Number of up-spin electrons | - | ~down_num~ | ~int64_t~ | Number of down-spin electrons | - | ~walk_num~ | ~int64_t~ | Number of walkers | - | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | - | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | - | ~provided~ | ~bool~ | If true, ~electron~ is valid | - | ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates | - | ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates | - | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | - | ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances | - | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | - | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron distances | - | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | - | ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~uninitialized~ | ~int32_t~ | Keeps bit set for uninitialized data | + | ~num~ | ~int64_t~ | Total number of electrons | + | ~up_num~ | ~int64_t~ | Number of up-spin electrons | + | ~down_num~ | ~int64_t~ | Number of down-spin electrons | + | ~walk_num~ | ~int64_t~ | Number of walkers | + | ~rescale_factor_kappa_ee~ | ~double~ | The distance scaling factor | + | ~rescale_factor_kappa_en~ | ~double~ | The distance scaling factor | + | ~provided~ | ~bool~ | If true, ~electron~ is valid | + | ~coord_new~ | ~double[walk_num][3][num]~ | New set of electron coordinates | + | ~coord_old~ | ~double[walk_num][3][num]~ | Old set of electron coordinates | + | ~coord_new_date~ | ~uint64_t~ | Last modification date of the coordinates | + | ~ee_distance~ | ~double[walk_num][num][num]~ | Electron-electron distances | + | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~en_distance~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | + | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~ee_distance_rescaled~ | ~double[walk_num][num][num]~ | Electron-electron rescaled distances | + | ~ee_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~ee_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives | + | ~ee_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + | ~en_distance_rescaled~ | ~double[walk_num][nucl_num][num]~ | Electron-nucleus distances | + | ~en_distance_rescaled_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][4][num][num]~ | Electron-electron rescaled distances derivatives | + | ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + ** Data structure @@ -97,13 +102,17 @@ typedef struct qmckl_electron_struct { int64_t ee_distance_date; int64_t en_distance_date; int64_t ee_distance_rescaled_date; + int64_t ee_distance_rescaled_deriv_e_date; int64_t en_distance_rescaled_date; + int64_t en_distance_rescaled_deriv_e_date; double* coord_new; double* coord_old; double* ee_distance; double* en_distance; double* ee_distance_rescaled; + double* ee_distance_rescaled_deriv_e; double* en_distance_rescaled; + double* en_distance_rescaled_deriv_e; int32_t uninitialized; bool provided; } qmckl_electron_struct; @@ -701,8 +710,9 @@ int64_t walk_num = chbrclf_walk_num; int64_t elec_num = chbrclf_elec_num; int64_t elec_up_num = chbrclf_elec_up_num; int64_t elec_dn_num = chbrclf_elec_dn_num; -double rescale_factor_kappa_ee = 2.0; -double rescale_factor_kappa_en = 3.0; +double rescale_factor_kappa_ee = 1.0; +double rescale_factor_kappa_en = 1.0; +double nucl_rescale_factor_kappa = 1.0; double* elec_coord = &(chbrclf_elec_coord[0][0][0]); int64_t nucl_num = chbrclf_nucl_num; @@ -803,7 +813,7 @@ for (int64_t i=0 ; i<3*elec_num ; ++i) { the dependencies are more recent than the date of the data to compute. If it is the case, then the data is recomputed and the current date is stored. - + ** Electron-electron distances *** Get @@ -980,7 +990,7 @@ qmckl_exit_code qmckl_compute_ee_distance ( #+end_src *** Test - + #+begin_src python :results output :exports none import numpy as np @@ -1035,6 +1045,15 @@ assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12); ** Electron-electron rescaled distances + ~ee_distance_rescaled~ stores the matrix of the rescaled distances between all + pairs of electrons: + + \[ + C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa + \] + + where \(C_{ij}\) is the matrix of electron-electron distances. + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes @@ -1127,12 +1146,12 @@ qmckl_exit_code qmckl_provide_ee_distance_rescaled(qmckl_context context) :END: #+NAME: qmckl_ee_distance_rescaled_args - | qmckl_context | context | in | Global state | - | int64_t | elec_num | in | Number of electrons | - | double | rescale_factor_kappa_ee | in | Factor to rescale ee distances | - | int64_t | walk_num | in | Number of walkers | - | double | coord[walk_num][3][elec_num] | in | Electron coordinates | - | double | ee_distance[walk_num][elec_num][elec_num] | out | Electron-electron distances | + | qmckl_context | context | in | Global state | + | int64_t | elec_num | in | Number of electrons | + | double | rescale_factor_kappa_ee | in | Factor to rescale ee distances | + | int64_t | walk_num | in | Number of walkers | + | double | coord[walk_num][3][elec_num] | in | Electron coordinates | + | double | ee_distance[walk_num][elec_num][elec_num] | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_ee_distance_rescaled_f(context, elec_num, rescale_factor_kappa_ee, walk_num, & @@ -1219,26 +1238,28 @@ qmckl_exit_code qmckl_compute_ee_distance_rescaled ( #+begin_src python :results output :exports none import numpy as np +kappa = 1.0 + elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ]) elec_2_w1 = np.array( [ 3.51983380318, -1.08717381954, -1.19617708027 ]) elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ]) elec_2_w2 = np.array( [ 3.17996025085, -1.40260577202, 1.49473607540 ]) -print ( "[0][0][0] : ", np.linalg.norm(elec_1_w1-elec_1_w1) ) -print ( "[0][1][0] : ", np.linalg.norm(elec_1_w1-elec_2_w1) ) -print ( "[1][0][0] : ", np.linalg.norm(elec_2_w1-elec_1_w1) ) -print ( "[0][0][1] : ", np.linalg.norm(elec_1_w2-elec_1_w2) ) -print ( "[0][1][1] : ", np.linalg.norm(elec_1_w2-elec_2_w2) ) -print ( "[1][0][1] : ", np.linalg.norm(elec_2_w2-elec_1_w2) ) +print ( "[0][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_1_w1)) )/kappa ) +print ( "[0][1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_2_w1)) )/kappa ) +print ( "[1][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-elec_1_w1)) )/kappa ) +print ( "[0][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-elec_1_w2)) )/kappa ) +print ( "[0][1][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-elec_2_w2)) )/kappa ) +print ( "[1][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w2-elec_1_w2)) )/kappa ) #+end_src #+RESULTS: : [0][0][0] : 0.0 - : [0][1][0] : 7.152322512964209 - : [1][0][0] : 7.152322512964209 + : [0][1][0] : 0.9992169566605263 + : [1][0][0] : 0.9992169566605263 : [0][0][1] : 0.0 - : [0][1][1] : 6.5517646321055665 - : [1][0][1] : 6.5517646321055665 + : [0][1][1] : 0.9985724058042633 + : [1][0][1] : 0.9985724058042633 #+begin_src c :tangle (eval c_test) assert(qmckl_electron_provided(context)); @@ -1247,6 +1268,230 @@ assert(qmckl_electron_provided(context)); double ee_distance_rescaled[walk_num * elec_num * elec_num]; rc = qmckl_get_electron_ee_distance_rescaled(context, ee_distance_rescaled); +// (e1,e2,w) +// (0,0,0) == 0. +assert(ee_distance_rescaled[0] == 0.); + +// (1,0,0) == (0,1,0) +assert(ee_distance_rescaled[1] == ee_distance_rescaled[elec_num]); + +// value of (1,0,0) +assert(fabs(ee_distance_rescaled[1]-0.9992169566605263) < 1.e-12); + +// (0,0,1) == 0. +assert(ee_distance_rescaled[elec_num*elec_num] == 0.); + +// (1,0,1) == (0,1,1) +assert(ee_distance_rescaled[elec_num*elec_num+1] == ee_distance_rescaled[elec_num*elec_num+elec_num]); + +// value of (1,0,1) +assert(fabs(ee_distance_rescaled[elec_num*elec_num+1]-0.9985724058042633) < 1.e-12); + + #+end_src + +** Electron-electron rescaled distance gradients and laplacian with respect to electron coords + + The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ + needs to be perturbed with respect to the electorn coordinates. + This data is stored in the ~ee_distance_rescaled_deriv_e~ tensor. The + The first three elements of this three index tensor ~[4][num][num]~ gives the + derivatives in the x, y, and z directions $dx, dy, dz$ and the last index + gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context context, double* const distance_rescaled_deriv_e); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_electron_ee_distance_rescaled_deriv_e(qmckl_context context, double* const distance_rescaled_deriv_e) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ee_distance_rescaled_deriv_e(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walk_num; + memcpy(distance_rescaled_deriv_e, ctx->electron.ee_distance_rescaled_deriv_e, 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_ee_distance_rescaled_deriv_e(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_ee_distance_rescaled_deriv_e(qmckl_context context) +{ + + 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); + + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->electron.ee_distance_rescaled_deriv_e_date) { + + /* Allocate array */ + if (ctx->electron.ee_distance_rescaled_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 4 * ctx->electron.num * ctx->electron.num * + ctx->electron.walk_num * sizeof(double); + double* ee_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (ee_distance_rescaled_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_ee_distance_rescaled_deriv_e", + NULL); + } + ctx->electron.ee_distance_rescaled_deriv_e = ee_distance_rescaled_deriv_e; + } + + qmckl_exit_code rc = + qmckl_compute_ee_distance_rescaled_deriv_e(context, + ctx->electron.num, + ctx->electron.rescale_factor_kappa_en, + ctx->electron.walk_num, + ctx->electron.coord_new, + ctx->electron.ee_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->electron.ee_distance_rescaled_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_ee_distance_rescaled_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_ee_distance_rescaled_deriv_e_args + | qmckl_context | context | in | Global state | + | int64_t | elec_num | in | Number of electrons | + | double | rescale_factor_kappa_ee | in | Factor to rescale ee distances | + | int64_t | walk_num | in | Number of walkers | + | double | coord[walk_num][3][elec_num] | in | Electron coordinates | + | double | ee_distance_deriv_e[walk_num][4][elec_num][elec_num] | out | Electron-electron rescaled distance derivatives | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_ee_distance_rescaled_deriv_e_f(context, elec_num, rescale_factor_kappa_ee, walk_num, & + coord, ee_distance_rescaled_deriv_e) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + double precision , intent(in) :: rescale_factor_kappa_ee + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: coord(elec_num,3,walk_num) + double precision , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num) + + integer*8 :: k + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + do k=1,walk_num + info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, elec_num, & + coord(1,1,k), elec_num, & + coord(1,1,k), elec_num, & + ee_distance_rescaled_deriv_e(1,1,1,k), elec_num, rescale_factor_kappa_ee) + if (info /= QMCKL_SUCCESS) then + exit + endif + end do + +end function qmckl_compute_ee_distance_rescaled_deriv_e_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_ee_distance_rescaled_deriv_e ( + const qmckl_context context, + const int64_t elec_num, + const double rescale_factor_kappa_ee, + const int64_t walk_num, + const double* coord, + double* const ee_distance_rescaled_deriv_e ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_ee_distance_rescaled_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_ee_distance_rescaled_deriv_e & + (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled_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 :: elec_num + real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: coord(elec_num,3,walk_num) + real (c_double ) , intent(out) :: ee_distance_rescaled_deriv_e(4,elec_num,elec_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_deriv_e_f + info = qmckl_compute_ee_distance_rescaled_deriv_e_f & + (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled_deriv_e) + + end function qmckl_compute_ee_distance_rescaled_deriv_e + #+end_src + +*** Test + + #+begin_src python :results output :exports none +import numpy as np + +# TODO + #+end_src + + #+begin_src c :tangle (eval c_test) +assert(qmckl_electron_provided(context)); + + +double ee_distance_rescaled_deriv_e[4 * walk_num * elec_num * elec_num]; +rc = qmckl_get_electron_ee_distance_rescaled_deriv_e(context, ee_distance_rescaled_deriv_e); + // TODO: Get exact values //// (e1,e2,w) //// (0,0,0) == 0. @@ -1269,6 +1514,7 @@ rc = qmckl_get_electron_ee_distance_rescaled(context, ee_distance_rescaled); #+end_src + ** Electron-nucleus distances *** Get @@ -1538,6 +1784,15 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); ** Electron-nucleus rescaled distances + ~en_distance_rescaled~ stores the matrix of the rescaled distances between + electrons and nucleii. + + \[ + C_{ij} = \left( 1 - \exp{-\kappa C_{ij}}\right)/\kappa + \] + + where \(C_{ij}\) is the matrix of electron-nucleus distances. + *** Get #+begin_src c :comments org :tangle (eval h_func) :noweb yes @@ -1747,6 +2002,8 @@ qmckl_exit_code qmckl_compute_en_distance_rescaled ( #+begin_src python :results output :exports none import numpy as np +kappa = 1.0 + elec_1_w1 = np.array( [ -2.26995253563, -5.15737533569, -2.22940072417 ]) elec_2_w1 = np.array( [ 3.51983380318, -1.08717381954, -1.19617708027 ]) elec_1_w2 = np.array( [ -2.34410619736, -3.20016115904, -1.53496759012 ]) @@ -1754,21 +2011,22 @@ elec_2_w2 = np.array( [ 3.17996025085, -1.40260577202, 1.49473607540 ]) nucl_1 = np.array( [ 1.096243353458458e+00, 8.907054016973815e-01, 7.777092280258892e-01 ] ) nucl_2 = np.array( [ 1.168459237342663e+00, 1.125660720053393e+00, 2.833370314829343e+00 ] ) -print ( "[0][0][0] : ", np.linalg.norm(elec_1_w1-nucl_1) ) -print ( "[0][1][0] : ", np.linalg.norm(elec_1_w1-nucl_2) ) -print ( "[0][0][1] : ", np.linalg.norm(elec_2_w1-nucl_1) ) -print ( "[1][0][0] : ", np.linalg.norm(elec_1_w2-nucl_1) ) -print ( "[1][1][0] : ", np.linalg.norm(elec_1_w2-nucl_2) ) -print ( "[1][0][1] : ", np.linalg.norm(elec_2_w2-nucl_1) ) +print ( "[0][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_1)) )/kappa ) +print ( "[0][1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_2)) )/kappa ) +print ( "[0][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-nucl_1)) )/kappa ) +print ( "[1][0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-nucl_1)) )/kappa ) +print ( "[1][1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w2-nucl_2)) )/kappa ) +print ( "[1][0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w2-nucl_1)) )/kappa ) + #+end_src #+RESULTS: - : [0][0][0] : 7.546738741619978 - : [0][1][0] : 8.77102435246984 - : [0][0][1] : 3.698922010513608 - : [1][0][0] : 5.824059436060509 - : [1][1][0] : 7.080482110317645 - : [1][0][1] : 3.1804527583077356 + : [0][0][0] : 0.9994721712909764 + : [0][1][0] : 0.9998448354439821 + : [0][0][1] : 0.9752498074577688 + : [1][0][0] : 0.9970444172399963 + : [1][1][0] : 0.9991586325813303 + : [1][0][1] : 0.9584331688679852 #+begin_src c :tangle (eval c_test) @@ -1788,6 +2046,275 @@ assert(qmckl_nucleus_provided(context)); double en_distance_rescaled[walk_num][nucl_num][elec_num]; rc = qmckl_get_electron_en_distance_rescaled(context, &(en_distance_rescaled[0][0][0])); + +assert (rc == QMCKL_SUCCESS); + +// (e,n,w) in Fortran notation +// (1,1,1) +assert(fabs(en_distance_rescaled[0][0][0] - 0.9994721712909764) < 1.e-12); + +// (1,2,1) +assert(fabs(en_distance_rescaled[0][1][0] - 0.9998448354439821) < 1.e-12); + +// (2,1,1) +assert(fabs(en_distance_rescaled[0][0][1] - 0.9752498074577688) < 1.e-12); + +// (1,1,2) +assert(fabs(en_distance_rescaled[1][0][0] - 0.9970444172399963) < 1.e-12); + +// (1,2,2) +assert(fabs(en_distance_rescaled[1][1][0] - 0.9991586325813303) < 1.e-12); + +// (2,1,2) +assert(fabs(en_distance_rescaled[1][0][1] - 0.9584331688679852) < 1.e-12); + + #+end_src + +** Electron-nucleus rescaled distance gradients and laplacian with respect to electron coords + + The rescaled distances which is given as $R = (1 - \exp{-\kappa r})/\kappa$ + needs to be perturbed with respect to the nuclear coordinates. + This data is stored in the ~en_distance_rescaled_deriv_e~ tensor. The + The first three elements of this three index tensor ~[4][num][num]~ gives the + derivatives in the x, y, and z directions $dx, dy, dz$ and the last index + gives the Laplacian $\partial x^2 + \partial y^2 + \partial z^2$. + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_electron_en_distance_rescaled_deriv_e(qmckl_context context, double* distance_rescaled_deriv_e) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_en_distance_rescaled_deriv_e(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num; + memcpy(distance_rescaled_deriv_e, ctx->electron.en_distance_rescaled_deriv_e, 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_en_distance_rescaled_deriv_e(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_en_distance_rescaled_deriv_e(qmckl_context context) +{ + + 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); + + if (!(ctx->nucleus.provided)) { + return QMCKL_NOT_PROVIDED; + } + + /* Compute if necessary */ + if (ctx->electron.coord_new_date > ctx->electron.en_distance_rescaled_deriv_e_date) { + + /* Allocate array */ + if (ctx->electron.en_distance_rescaled_deriv_e == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 4 * ctx->electron.num * ctx->nucleus.num * + ctx->electron.walk_num * sizeof(double); + double* en_distance_rescaled_deriv_e = (double*) qmckl_malloc(context, mem_info); + + if (en_distance_rescaled_deriv_e == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_en_distance_rescaled_deriv_e", + NULL); + } + ctx->electron.en_distance_rescaled_deriv_e = en_distance_rescaled_deriv_e; + } + + qmckl_exit_code rc = + qmckl_compute_en_distance_rescaled_deriv_e(context, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.rescale_factor_kappa_en, + ctx->electron.walk_num, + ctx->electron.coord_new, + ctx->nucleus.coord, + ctx->electron.en_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->electron.en_distance_rescaled_deriv_e_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_en_distance_rescaled_deriv_e + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_en_distance_rescaled_deriv_e_args + | qmckl_context | context | in | Global state | + | int64_t | elec_num | in | Number of electrons | + | int64_t | nucl_num | in | Number of nuclei | + | double | rescale_factor_kappa_en | in | The factor for rescaled distances | + | int64_t | walk_num | in | Number of walkers | + | double | elec_coord[walk_num][3][elec_num] | in | Electron coordinates | + | double | nucl_coord[3][elec_num] | in | Nuclear coordinates | + | double | en_distance_rescaled_deriv_e_date[walk_num][4][nucl_num][elec_num] | out | Electron-nucleus distance derivatives | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_en_distance_rescaled_deriv_e_f(context, elec_num, nucl_num, & + rescale_factor_kappa_en, walk_num, elec_coord, & + nucl_coord, en_distance_rescaled_deriv_e) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + integer*8 , intent(in) :: nucl_num + double precision , intent(in) :: rescale_factor_kappa_en + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: elec_coord(elec_num,3,walk_num) + double precision , intent(in) :: nucl_coord(nucl_num,3) + double precision , intent(out) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,walk_num) + + integer*8 :: k + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (elec_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + if (nucl_num <= 0) then + info = QMCKL_INVALID_ARG_3 + return + endif + + ! TODO: comparison with 0 + !if (rescale_factor_kappa_en <= 0) then + ! info = QMCKL_INVALID_ARG_4 + ! return + !endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + do k=1,walk_num + info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, nucl_num, & + elec_coord(1,1,k), elec_num, & + nucl_coord, nucl_num, & + en_distance_rescaled_deriv_e(1,1,k), elec_num, rescale_factor_kappa_en) + if (info /= QMCKL_SUCCESS) then + exit + endif + end do + +end function qmckl_compute_en_distance_rescaled_deriv_e_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_en_distance_rescaled_deriv_e ( + const qmckl_context context, + const int64_t elec_num, + const int64_t nucl_num, + const double rescale_factor_kappa_en, + const int64_t walk_num, + const double* elec_coord, + const double* nucl_coord, + double* const en_distance_rescaled_deriv_e ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_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_en_distance_rescaled_deriv_e & + (context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled_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 :: elec_num + integer (c_int64_t) , intent(in) , value :: nucl_num + real (c_double ) , intent(in) , value :: rescale_factor_kappa_en + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: elec_coord(elec_num,3,walk_num) + real (c_double ) , intent(in) :: nucl_coord(elec_num,3) + real (c_double ) , intent(out) :: en_distance_rescaled_deriv_e(elec_num,nucl_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_deriv_e_f + info = qmckl_compute_en_distance_rescaled_deriv_e_f & + (context, elec_num, nucl_num, rescale_factor_kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled_deriv_e) + + end function qmckl_compute_en_distance_rescaled_deriv_e + #+end_src + +*** Test + + #+begin_src python :results output :exports none +import numpy as np + +# TODO + #+end_src + + + #+begin_src c :tangle (eval c_test) + +assert(qmckl_electron_provided(context)); + +rc = qmckl_set_nucleus_num (context, nucl_num); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_set_nucleus_charge (context, charge); +assert (rc == QMCKL_SUCCESS); + +rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord); +assert (rc == QMCKL_SUCCESS); + +assert(qmckl_nucleus_provided(context)); + +double en_distance_rescaled_deriv_e[walk_num][4][nucl_num][elec_num]; + +rc = qmckl_get_electron_en_distance_rescaled_deriv_e(context, &(en_distance_rescaled_deriv_e[0][0][0][0])); + assert (rc == QMCKL_SUCCESS); // TODO: check exact values diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index 99237af..cc7dce6 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -426,7 +426,7 @@ qmckl_set_nucleus_rescale_factor(qmckl_context context, const double rescale_fac if (rescale_factor_kappa <= 0.0) { return qmckl_failwith( context, QMCKL_INVALID_ARG_2, - "qmckl_set_nucleus_kappa", + "qmckl_set_nucleus_rescale_factor", "rescale_factor_kappa cannot be <= 0."); }