From cb70c1f5684761a781be885eddd5f82f0dc52b96 Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Tue, 1 Jun 2021 13:01:14 +0530 Subject: [PATCH] Created a function to provide the derivative functions along with Doc. #17 --- org/qmckl_distance.org | 354 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 354 insertions(+) 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)