From c02d8871d853c31cc78a4b18dfe8f9a02924d4ec Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Wed, 26 May 2021 12:52:28 +0530 Subject: [PATCH] Function to calculate rescaled distances for elec-nucl. #15 --- org/qmckl_distance.org | 15 ++- org/qmckl_electron.org | 279 +++++++++++++++++++++++++++++++++++++++++ org/qmckl_nucleus.org | 12 +- 3 files changed, 297 insertions(+), 9 deletions(-) diff --git a/org/qmckl_distance.org b/org/qmckl_distance.org index c9ef84c..5b9737f 100644 --- a/org/qmckl_distance.org +++ b/org/qmckl_distance.org @@ -807,6 +807,7 @@ assert(test_qmckl_dist(context) == QMCKL_SUCCESS); | int64_t | ldb | in | Leading dimension of array ~B~ | | double | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$ | | int64_t | ldc | in | Leading dimension of array ~C~ | + | double | kappa | in | Factor for calculating rescaled distances | *** Requirements @@ -839,13 +840,14 @@ assert(test_qmckl_dist(context) == QMCKL_SUCCESS); const double* B, const int64_t ldb, double* const C, - const int64_t ldc ); + const int64_t ldc, + const double kappa); #+end_src *** Source #+begin_src f90 :tangle (eval f) integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, & - A, LDA, B, LDB, C, LDC) & + A, LDA, B, LDB, C, LDC, kappa) & result(info) use qmckl implicit none @@ -858,6 +860,7 @@ integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, & real*8 , intent(in) :: B(ldb,*) integer*8 , intent(in) :: ldc real*8 , intent(out) :: C(ldc,*) + real*8 , intent(in) :: kappa integer*8 :: i,j real*8 :: x, y, z @@ -988,7 +991,7 @@ end function qmckl_distance_rescaled_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_distance_rescaled & - (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) & + (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, kappa) & bind(C) result(info) use, intrinsic :: iso_c_binding @@ -1005,10 +1008,11 @@ end function qmckl_distance_rescaled_f integer (c_int64_t) , intent(in) , value :: ldb real (c_double ) , intent(out) :: C(ldc,n) integer (c_int64_t) , intent(in) , value :: ldc + real (c_double ) , intent(in) , value :: kappa integer(c_int32_t), external :: qmckl_distance_rescaled_f info = qmckl_distance_rescaled_f & - (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) + (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, kappa) end function qmckl_distance_rescaled #+end_src @@ -1019,7 +1023,7 @@ end function qmckl_distance_rescaled_f #+begin_src f90 :tangle (eval fh_func) :comments org :exports none interface integer(c_int32_t) function qmckl_distance_rescaled & - (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) & + (context, transa, transb, m, n, A, lda, B, ldb, C, ldc, kappa) & bind(C) use, intrinsic :: iso_c_binding import @@ -1036,6 +1040,7 @@ end function qmckl_distance_rescaled_f integer (c_int64_t) , intent(in) , value :: ldb real (c_double ) , intent(out) :: C(ldc,n) integer (c_int64_t) , intent(in) , value :: ldc + real (c_double ) , intent(in) , value :: kappa end function qmckl_distance_rescaled end interface diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index aa6b600..153f48b 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -1276,6 +1276,285 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); #+end_src +** Electron-nucleus rescaled distances + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_electron_en_distance_rescaled(qmckl_context context, double* distance_rescaled) +{ + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_en_distance_rescaled(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walk_num; + memcpy(distance_rescaled, ctx->electron.en_distance_rescaled, 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(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(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_date) { + + /* Allocate array */ + if (ctx->electron.en_distance_rescaled == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.num * ctx->nucleus.num * + ctx->electron.walk_num * sizeof(double); + double* en_distance_rescaled = (double*) qmckl_malloc(context, mem_info); + + if (en_distance_rescaled == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_en_distance_rescaled", + NULL); + } + ctx->electron.en_distance_rescaled = en_distance_rescaled; + } + + qmckl_exit_code rc = + qmckl_compute_en_distance_rescaled(context, + ctx->electron.num, + ctx->nucleus.num, + ctx->electron.kappa_en, + ctx->electron.walk_num, + ctx->electron.coord_new, + ctx->nucleus.coord, + ctx->electron.en_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->electron.en_distance_rescaled_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_en_distance_rescaled + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_en_distance_rescaled_args + | qmckl_context | context | in | Global state | + | int64_t | elec_num | in | Number of electrons | + | int64_t | nucl_num | in | Number of nuclei | + | double | 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_date[walk_num][nucl_num][elec_num] | out | Electron-nucleus distances | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_en_distance_rescaled_f(context, elec_num, nucl_num, kappa_en, walk_num, elec_coord, & + nucl_coord, en_distance_rescaled) & + 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) :: 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(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 (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(context, 'T', 'T', elec_num, nucl_num, & + elec_coord(1,1,k), elec_num, & + nucl_coord, nucl_num, & + en_distance_rescaled(1,1,k), elec_num, kappa_en) + if (info /= QMCKL_SUCCESS) then + exit + endif + end do + +end function qmckl_compute_en_distance_rescaled_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_en_distance_rescaled ( + const qmckl_context context, + const int64_t elec_num, + const int64_t nucl_num, + const double kappa_en, + const int64_t walk_num, + const double* elec_coord, + const double* nucl_coord, + double* const en_distance_rescaled ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_en_distance_rescaled_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 & + (context, elec_num, nucl_num, kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled) & + 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 :: 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(elec_num,nucl_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_en_distance_rescaled_f + info = qmckl_compute_en_distance_rescaled_f & + (context, elec_num, nucl_num, kappa_en, walk_num, elec_coord, nucl_coord, en_distance_rescaled) + + end function qmckl_compute_en_distance_rescaled + #+end_src + +*** Test + + #+begin_src python :results output :exports none +import numpy as np + +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 ]) +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) ) + #+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 + + #+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_kappa (context, nucl_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[walk_num][nucl_num][elec_num]; + +rc = qmckl_get_electron_en_distance_rescaled(context, &(en_distance[0][0][0])); +assert (rc == QMCKL_SUCCESS); + +// TODO: check exact values +//// (e,n,w) in Fortran notation +//// (1,1,1) +//assert(fabs(en_distance[0][0][0] - 7.546738741619978) < 1.e-12); +// +//// (1,2,1) +//assert(fabs(en_distance[0][1][0] - 8.77102435246984) < 1.e-12); +// +//// (2,1,1) +//assert(fabs(en_distance[0][0][1] - 3.698922010513608) < 1.e-12); +// +//// (1,1,2) +//assert(fabs(en_distance[1][0][0] - 5.824059436060509) < 1.e-12); +// +//// (1,2,2) +//assert(fabs(en_distance[1][1][0] - 7.080482110317645) < 1.e-12); +// +//// (2,1,2) +//assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); + + #+end_src + * End of files :noexport: #+begin_src c :tangle (eval h_private_type) diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index 418c368..1c890c9 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -780,6 +780,7 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context) qmckl_exit_code rc = qmckl_compute_nn_distance_rescaled(context, ctx->nucleus.num, + ctx->nucleus.kappa, ctx->nucleus.coord, ctx->nucleus.nn_distance_rescaled); if (rc != QMCKL_SUCCESS) { @@ -801,12 +802,13 @@ qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context) | double | nn_distance_rescaled[nucl_num][nucl_num] | out | Nucleus-nucleus rescaled distances (au) | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, coord, nn_distance_rescaled) & +integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, kappa, coord, nn_distance_rescaled) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: nucl_num + double precision , intent(in) :: kappa double precision , intent(in) :: coord(nucl_num,3) double precision , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num) @@ -827,7 +829,7 @@ integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, coord, info = qmckl_distance_rescaled(context, 'T', 'T', nucl_num, nucl_num, & coord, nucl_num, & coord, nucl_num, & - nn_distance_rescaled, nucl_num) + nn_distance_rescaled, nucl_num, kappa) end function qmckl_compute_nn_distance_rescaled_f #+end_src @@ -836,6 +838,7 @@ end function qmckl_compute_nn_distance_rescaled_f qmckl_exit_code qmckl_compute_nn_distance_rescaled ( const qmckl_context context, const int64_t nucl_num, + const double kappa, const double* coord, double* const nn_distance_rescaled ); #+end_src @@ -846,7 +849,7 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled ( #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_nn_distance_rescaled & - (context, nucl_num, coord, nn_distance_rescaled) & + (context, nucl_num, kappa, coord, nn_distance_rescaled) & bind(C) result(info) use, intrinsic :: iso_c_binding @@ -854,12 +857,13 @@ qmckl_exit_code qmckl_compute_nn_distance_rescaled ( integer (c_int64_t) , intent(in) , value :: context integer (c_int64_t) , intent(in) , value :: nucl_num + real (c_double ) , intent(in) , value :: kappa real (c_double ) , intent(in) :: coord(nucl_num,3) real (c_double ) , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num) integer(c_int32_t), external :: qmckl_compute_nn_distance_rescaled_f info = qmckl_compute_nn_distance_rescaled_f & - (context, nucl_num, coord, nn_distance_rescaled) + (context, nucl_num, kappa, coord, nn_distance_rescaled) end function qmckl_compute_nn_distance_rescaled #+end_src