diff --git a/Makefile.am b/Makefile.am index f6d79a7..472ee62 100644 --- a/Makefile.am +++ b/Makefile.am @@ -118,7 +118,7 @@ python-install: $(qmckl_h) $(lib_LTLIBRARIES) $(dist_python_DATA) cp src/qmckl.py . ; \ export QMCKL_INCLUDEDIR="$(prefix)/include" ; \ export QMCKL_LIBDIR="$(prefix)/lib" ; \ - pip install . + python3 -m pip install . python-test: $(test_py) cd $(abs_srcdir)/python/test/ && \ diff --git a/README.md b/README.md index 0690294..e738b24 100644 --- a/README.md +++ b/README.md @@ -75,6 +75,27 @@ sudo make install sudo make installcheck ``` +## Python API + +- [SWIG](https://www.swig.org) (>= 4.0) is required to build the Python API for maintainers + +In order to install the `qmckl` Python package, first install the shared C library +`libqmckl` following the installation guide above and then run the following command: +``` +make python-install +``` + +To test the installation, run +``` +make python-test +``` + +Minimal example demonstrating the use of the `qmckl` Python API can be found in the +[test_api.py](https://github.com/TREX-CoE/qmckl/blob/master/python/test/test_api.py) file. + +We highly recommend to use +[virtual environments](https://docs.python.org/3/tutorial/venv.html) +to avoid compatibility issues and to improve reproducibility. ## Installation procedure for Guix users diff --git a/org/qmckl_ao.org b/org/qmckl_ao.org index 13e42d4..6187827 100644 --- a/org/qmckl_ao.org +++ b/org/qmckl_ao.org @@ -5648,33 +5648,14 @@ end function qmckl_compute_ao_value_doc_f | ~coef_normalized~ | ~double[prim_num]~ | in | Value, gradients and Laplacian of the shells | | ~ao_value~ | ~double[point_num][ao_num]~ | out | Values of the AOs | - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -#ifdef HAVE_HPC -qmckl_exit_code -qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context, - const int64_t ao_num, - const int64_t shell_num, - const int32_t* restrict prim_num_per_nucleus, - const int64_t point_num, - const int64_t nucl_num, - const double* restrict coord, - const double* restrict nucl_coord, - const int64_t* restrict nucleus_index, - const int64_t* restrict nucleus_shell_num, - const double* nucleus_range, - const int32_t* restrict nucleus_max_ang_mom, - const int32_t* restrict shell_ang_mom, - const double* restrict ao_factor, - const qmckl_matrix expo_per_nucleus, - const qmckl_tensor coef_per_nucleus, - double* restrict const ao_value ) -{ - int32_t lstart[32]; + #+NAME:ao_value_constants + #+begin_src c :exports none + int32_t lstart[32] __attribute__((aligned(64))); for (int32_t l=0 ; l<32 ; ++l) { lstart[l] = l*(l+1)*(l+2)/6; } - int64_t ao_index[shell_num+1]; + int64_t ao_index[shell_num+1] __attribute__((aligned(64))); int64_t size_max = 0; int64_t prim_max = 0; int64_t shell_max = 0; @@ -5700,17 +5681,42 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context, /* Don't compute polynomials when the radial part is zero. */ double cutoff = 27.631021115928547; // -log(1.e-12) + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +#ifdef HAVE_HPC +qmckl_exit_code +qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context, + const int64_t ao_num, + const int64_t shell_num, + const int32_t* restrict prim_num_per_nucleus, + const int64_t point_num, + const int64_t nucl_num, + const double* restrict coord, + const double* restrict nucl_coord, + const int64_t* restrict nucleus_index, + const int64_t* restrict nucleus_shell_num, + const double* nucleus_range, + const int32_t* restrict nucleus_max_ang_mom, + const int32_t* restrict shell_ang_mom, + const double* restrict ao_factor, + const qmckl_matrix expo_per_nucleus, + const qmckl_tensor coef_per_nucleus, + double* restrict const ao_value ) +{ + <> + #ifdef HAVE_OPENMP #pragma omp parallel #endif { qmckl_exit_code rc; - double ar2[prim_max]; - int32_t powers[3*size_max]; - double poly_vgl[5*size_max]; + double ar2[prim_max] __attribute__((aligned(64))); + int32_t powers[3*size_max] __attribute__((aligned(64))); + double poly_vgl[5*size_max] __attribute__((aligned(64))); - double exp_mat[prim_max]; - double ce_mat[shell_max]; + double exp_mat[prim_max] __attribute__((aligned(64))); + double ce_mat[shell_max] __attribute__((aligned(64))); double coef_mat[nucl_num][shell_max][prim_max]; for (int i=0 ; i 0) { const double* restrict f = ao_factor + k; @@ -5827,10 +5839,10 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context, poly_vgl_1 = &(poly_vgl[idx]); switch (n) { - case(1): + case 1: ao_value_1[0] = s1 * f[0]; break; - case (3): + case 3: #ifdef HAVE_OPENMP #pragma omp simd #endif @@ -5848,7 +5860,7 @@ qmckl_compute_ao_value_hpc_gaussian (const qmckl_context context, break; default: #ifdef HAVE_OPENMP -#pragma omp simd simdlen(8) +#pragma omp simd #endif for (int il=0 ; il prim_max ? - prim_num_per_nucleus[inucl] : prim_max; - shell_max = nucleus_shell_num[inucl] > shell_max ? - nucleus_shell_num[inucl] : shell_max; - const int64_t ishell_start = nucleus_index[inucl]; - const int64_t ishell_end = nucleus_index[inucl] + nucleus_shell_num[inucl]; - for (int64_t ishell = ishell_start ; ishell < ishell_end ; ++ishell) { - const int l = shell_ang_mom[ishell]; - ao_index[ishell] = k; - k += lstart[l+1] - lstart[l]; - size_max = size_max < lstart[l+1] ? lstart[l+1] : size_max; - } - } - ao_index[shell_num] = ao_num+1; - } - - /* Don't compute when the radial part is zero. */ - double cutoff = 27.631021115928547; // -log(1.e-12) + <> #ifdef HAVE_OPENMP #pragma omp parallel @@ -6508,18 +6491,6 @@ qmckl_compute_ao_vgl_hpc_gaussian ( qmckl_exit_code rc; double ar2[prim_max] __attribute__((aligned(64))); int32_t powers[3*size_max] __attribute__((aligned(64))); - double poly_vgl_l1[4][4] __attribute__((aligned(64))) = - {{1.0, 0.0, 0.0, 0.0}, - {0.0, 1.0, 0.0, 0.0}, - {0.0, 0.0, 1.0, 0.0}, - {0.0, 0.0, 0.0, 1.0}}; - double poly_vgl_l2[5][10]__attribute__((aligned(64))) = - {{1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, - {0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, - {0., 0., 0., 0., 2., 0., 0., 2., 0., 2.}}; - double poly_vgl[5][size_max] __attribute__((aligned(64))); double exp_mat[prim_max][8] __attribute__((aligned(64))) ; double ce_mat[shell_max][8] __attribute__((aligned(64))) ; @@ -6533,20 +6504,33 @@ qmckl_compute_ao_vgl_hpc_gaussian ( } } + double poly_vgl_l1[4][4] __attribute__((aligned(64))) = + {{1.0, 0.0, 0.0, 0.0}, + {0.0, 1.0, 0.0, 0.0}, + {0.0, 0.0, 1.0, 0.0}, + {0.0, 0.0, 0.0, 1.0}}; + double poly_vgl_l2[5][10]__attribute__((aligned(64))) = + {{1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, + {0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, + {0., 0., 0., 0., 2., 0., 0., 2., 0., 2.}}; + double poly_vgl[5][size_max] __attribute__((aligned(64))); + #ifdef HAVE_OPENMP #pragma omp for #endif for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) { const double e_coord[3] __attribute__((aligned(64))) = - { coord[ipoint], - coord[ipoint + point_num], - coord[ipoint + 2*point_num] }; - + { coord[ipoint], + coord[ipoint + point_num], + coord[ipoint + 2*point_num] }; + for (int64_t inucl=0 ; inucl < nucl_num ; ++inucl) { const double n_coord[3] __attribute__((aligned(64))) = - { nucl_coord[inucl], - nucl_coord[inucl + nucl_num], - nucl_coord[inucl + 2*nucl_num] }; + { nucl_coord[inucl], + nucl_coord[inucl + nucl_num], + nucl_coord[inucl + 2*nucl_num] }; /* Test if the point is in the range of the nucleus */ const double x = e_coord[0] - n_coord[0]; @@ -6616,6 +6600,7 @@ qmckl_compute_ao_vgl_hpc_gaussian ( for (int64_t iprim = 0 ; iprim < nidx ; ++iprim) { exp_mat[iprim][0] = exp(-ar2[iprim]); } + for (int64_t iprim = 0 ; iprim < nidx ; ++iprim) { double f = qmckl_mat(expo_per_nucleus, iprim, inucl) * exp_mat[iprim][0]; f = -f-f; @@ -6628,21 +6613,6 @@ qmckl_compute_ao_vgl_hpc_gaussian ( /* --- */ switch (8) { - case(5): - - for (int i=0 ; i 0~ + - ~n > 0~ + - ~k > 0~ + - ~lda >= m~ + - ~ldb >= n~ + - ~ldc >= n~ + - ~A~ is allocated with at least $m \times k \times 8$ bytes + - ~B~ is allocated with at least $k \times n \times 8$ bytes + - ~C~ is allocated with at least $m \times n \times 8$ bytes + - ~size_max_A >= m * k~ + - ~size_max_B >= k * n~ + - ~size_max_C >= m * n~ + + #+CALL: generate_c_header(table=qmckl_dgemm_safe_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm_safe") + + #+RESULTS: + #+BEGIN_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_dgemm_safe ( + const qmckl_context context, + const char TransA, + const char TransB, + const int64_t m, + const int64_t n, + const int64_t k, + const double alpha, + const double* A, + const int64_t size_max_A, + const int64_t lda, + const double* B, + const int64_t size_max_B, + const int64_t ldb, + const double beta, + double* const C, + const int64_t size_max_C, + const int64_t ldc ); + #+END_src + + #+begin_src f90 :tangle (eval f) :exports none +integer function qmckl_dgemm_safe_f(context, TransA, TransB, & + m, n, k, alpha, A, size_A, LDA, B, size_B, LDB, beta, C, size_C, LDC) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + character , intent(in) :: TransA, TransB + integer*8 , intent(in) :: m, n, k + double precision , intent(in) :: alpha, beta + integer*8 , intent(in) :: lda + integer*8 , intent(in) :: size_A + double precision , intent(in) :: A(lda,*) + integer*8 , intent(in) :: ldb + integer*8 , intent(in) :: size_B + double precision , intent(in) :: B(ldb,*) + integer*8 , intent(in) :: ldc + integer*8 , intent(in) :: size_C + double precision , intent(out) :: C(ldc,*) + + 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 (k <= 0_8) then + info = QMCKL_INVALID_ARG_6 + return + endif + + if (LDA <= 0) then + info = QMCKL_INVALID_ARG_10 + return + endif + + if (LDB <= 0) then + info = QMCKL_INVALID_ARG_13 + return + endif + + if (LDC <= 0) then + info = QMCKL_INVALID_ARG_17 + return + endif + + if (size_A <= 0) then + info = QMCKL_INVALID_ARG_9 + return + endif + + if (size_B <= 0) then + info = QMCKL_INVALID_ARG_12 + return + endif + + if (size_C <= 0) then + info = QMCKL_INVALID_ARG_16 + return + endif + + call dgemm(transA, transB, int(m,4), int(n,4), int(k,4), & + alpha, A, int(LDA,4), B, int(LDB,4), beta, C, int(LDC,4)) + +end function qmckl_dgemm_safe_f + #+end_src + +*** C interface :noexport: + + #+CALL: generate_c_interface(table=qmckl_dgemm_safe_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm_safe") + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_dgemm_safe & + (context, TransA, TransB, m, n, k, alpha, A, size_A, lda, B, size_B, ldb, beta, C, size_C, ldc) & + 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 + integer (c_int64_t) , intent(in) , value :: k + real (c_double ) , intent(in) , value :: alpha + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: size_A + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: size_B + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(in) , value :: beta + real (c_double ) , intent(out) :: C(ldc,*) + integer (c_int64_t) , intent(in) , value :: size_C + integer (c_int64_t) , intent(in) , value :: ldc + + integer(c_int32_t), external :: qmckl_dgemm_safe_f + info = qmckl_dgemm_safe_f & + (context, TransA, TransB, m, n, k, alpha, A, size_A, lda, B, size_B, ldb, beta, C, size_C, ldc) + + end function qmckl_dgemm_safe + #+end_src + + + #+CALL: generate_f_interface(table=qmckl_dgemm_safe_args,rettyp="qmckl_exit_code",fname="qmckl_dgemm_safe") + + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_dgemm_safe & + (context, TransA, TransB, m, n, k, alpha, A, size_A, lda, B, size_B, ldb, beta, C, size_C, ldc) & + 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 + integer (c_int64_t) , intent(in) , value :: k + real (c_double ) , intent(in) , value :: alpha + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: size_A + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: size_B + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(in) , value :: beta + real (c_double ) , intent(out) :: C(ldc,*) + integer (c_int64_t) , intent(in) , value :: size_C + integer (c_int64_t) , intent(in) , value :: ldc + + end function qmckl_dgemm_safe + end interface + #+end_src + ** ~qmckl_matmul~ Matrix multiplication using the =qmckl_matrix= data type: @@ -1452,8 +1671,6 @@ integer function qmckl_adjugate_f(context, na, A, LDA, B, ldb, det_l) & integer*8, intent(in) :: na double precision, intent(inout) :: det_l - integer :: i,j - info = QMCKL_SUCCESS if (context == QMCKL_NULL_CONTEXT) then @@ -2274,6 +2491,159 @@ qmckl_exit_code test_qmckl_adjugate(qmckl_context context); assert(QMCKL_SUCCESS == test_qmckl_adjugate(context)); #+end_src +** ~qmckl_adjugate_safe~ + + "Size-safe" proxy function with the same functionality as ~qmckl_adjugate~ + but with 2 additional arguments. These arguments ~size_max_M~ (where ~M~ is a matix) + are required primarily for the Python API, where compatibility with + NumPy arrays implies that sizes of the input and output arrays are provided. + + + #+NAME: qmckl_adjugate_safe_args + | Variable | Type | In/Out | Description | + |--------------+-----------------+--------+------------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~n~ | ~int64_t~ | in | Number of rows and columns of the input matrix | + | ~A~ | ~double[][lda]~ | in | Array containing the $n \times n$ matrix $A$ | + | ~size_max_A~ | ~int64_t~ | in | Size of the matrix $A$ | + | ~lda~ | ~int64_t~ | in | Leading dimension of array ~A~ | + | ~B~ | ~double[][ldb]~ | out | Adjugate of $A$ | + | ~size_max_B~ | ~int64_t~ | in | Size of the matrix $B$ | + | ~ldb~ | ~int64_t~ | in | Leading dimension of array ~B~ | + | ~det_l~ | ~double~ | inout | determinant of $A$ | + + Requirements: + + - ~context~ is not ~QMCKL_NULL_CONTEXT~ + - ~n > 0~ + - ~lda >= m~ + - ~A~ is allocated with at least $m \times m \times 8$ bytes + - ~ldb >= m~ + - ~B~ is allocated with at least $m \times m \times 8$ bytes + - ~size_max_A >= m * m~ + - ~size_max_B >= m * m~ + + #+CALL: generate_c_header(table=qmckl_adjugate_safe_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate_safe") + + #+RESULTS: + #+BEGIN_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_adjugate_safe ( + const qmckl_context context, + const int64_t n, + const double* A, + const int64_t size_max_A, + const int64_t lda, + double* const B, + const int64_t size_max_B, + const int64_t ldb, + double* det_l ); + #+END_src + + For small matrices (\le 5\times 5), we use specialized routines + for performance motivations. For larger sizes, we rely on the + LAPACK library. + + #+begin_src f90 :tangle (eval f) :exports none +integer function qmckl_adjugate_safe_f(context, & + na, A, size_A, LDA, B, size_B, LDB, det_l) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + double precision, intent(in) :: A (LDA,*) + integer*8, intent(in) :: size_A + integer*8, intent(in) :: LDA + double precision, intent(out) :: B (LDB,*) + integer*8, intent(in) :: size_B + integer*8, intent(in) :: LDB + integer*8, intent(in) :: na + double precision, intent(inout) :: det_l + + integer, external :: qmckl_adjugate_f + + info = QMCKL_SUCCESS + + if (size_A < na) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (size_B <= 0_8) then + info = QMCKL_INVALID_ARG_7 + return + endif + + info = qmckl_adjugate_f(context, na, A, LDA, B, LDB, det_l) + + if (info == QMCKL_INVALID_ARG_4) then + info = QMCKL_INVALID_ARG_5 + return + endif + + if (info == QMCKL_INVALID_ARG_6) then + info = QMCKL_INVALID_ARG_8 + return + endif + +end function qmckl_adjugate_safe_f + #+end_src + +*** C interface + + #+CALL: generate_c_interface(table=qmckl_adjugate_safe_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate_safe") + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_adjugate_safe & + (context, n, A, size_A, lda, B, size_B, ldb, det_l) & + 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 :: n + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + integer (c_int64_t) , intent(in) , value :: size_A + real (c_double ) , intent(out) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + integer (c_int64_t) , intent(in) , value :: size_B + real (c_double ) , intent(inout) :: det_l + + integer(c_int32_t), external :: qmckl_adjugate_safe_f + info = qmckl_adjugate_safe_f & + (context, n, A, size_A, lda, B, size_B, ldb, det_l) + + end function qmckl_adjugate_safe + #+end_src + + #+CALL: generate_f_interface(table=qmckl_adjugate_safe_args,rettyp="qmckl_exit_code",fname="qmckl_adjugate_safe") + + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_adjugate_safe & + (context, n, A, size_A, lda, B, size_B, ldb, det_l) & + 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 :: n + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + integer (c_int64_t) , intent(in) , value :: size_A + real (c_double ) , intent(out) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + integer (c_int64_t) , intent(in) , value :: size_B + real (c_double ) , intent(inout) :: det_l + + end function qmckl_adjugate_safe + end interface + #+end_src + ** ~qmckl_transpose~ Transposes a matrix: $A^\dagger_{ji} = A_{ij}$. @@ -2373,9 +2743,7 @@ qmckl_transpose (qmckl_context context, assert (qmckl_context_destroy(context) == QMCKL_SUCCESS); return 0; } - #+end_src - # -*- mode: org -*- # vim: syntax=c diff --git a/org/qmckl_context.org b/org/qmckl_context.org index dfcd906..cb00776 100644 --- a/org/qmckl_context.org +++ b/org/qmckl_context.org @@ -292,6 +292,9 @@ qmckl_context qmckl_context_create() { rc = qmckl_init_determinant(context); assert (rc == QMCKL_SUCCESS); + + rc = qmckl_init_jastrow(context); + assert (rc == QMCKL_SUCCESS); } /* Allocate qmckl_memory_struct */ diff --git a/org/qmckl_determinant.org b/org/qmckl_determinant.org index ba46707..78d54f1 100644 --- a/org/qmckl_determinant.org +++ b/org/qmckl_determinant.org @@ -1111,21 +1111,21 @@ const double* nucl_charge = chbrclf_charge; const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); rc = qmckl_set_electron_num (context, chbrclf_elec_up_num, chbrclf_elec_dn_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); + +rc = qmckl_set_electron_coord (context, 'N', chbrclf_walk_num, elec_coord, chbrclf_walk_num*chbrclf_elec_num*3); +qmckl_check(context, rc); assert(qmckl_electron_provided(context)); -rc = qmckl_set_electron_coord (context, 'N', chbrclf_walk_num, elec_coord, chbrclf_walk_num*chbrclf_elec_num*3); -assert(rc == QMCKL_SUCCESS); - rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), chbrclf_nucl_num*3); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_nucleus_charge(context, nucl_charge, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_nucleus_provided(context)); @@ -1145,57 +1145,57 @@ const char typ = 'G'; assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_type (context, typ); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_num (context, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_prim_num (context, chbrclf_prim_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_nucleus_index (context, nucleus_index, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_nucleus_shell_num (context, nucleus_shell_num, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_ang_mom (context, shell_ang_mom, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_factor (context, shell_factor, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_prim_num (context, shell_prim_num, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_shell_prim_index (context, shell_prim_index, chbrclf_shell_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_exponent (context, exponent, chbrclf_prim_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_coefficient (context, coefficient, chbrclf_prim_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_ao_basis_provided(context)); rc = qmckl_set_ao_basis_prim_factor (context, prim_factor, chbrclf_prim_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_ao_basis_ao_num(context, chbrclf_ao_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_ao_basis_ao_factor (context, ao_factor, chbrclf_ao_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_ao_basis_provided(context)); @@ -1203,22 +1203,22 @@ assert(qmckl_ao_basis_provided(context)); double ao_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_ao_num]; rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), (int64_t) 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_ao_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); /* Set up MO data */ rc = qmckl_set_mo_basis_mo_num(context, chbrclf_mo_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); const double * mo_coefficient = &(chbrclf_mo_coef[0]); rc = qmckl_set_mo_basis_coefficient(context, mo_coefficient); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_mo_basis_provided(context)); double mo_vgl[chbrclf_walk_num*chbrclf_elec_num][5][chbrclf_mo_num]; rc = qmckl_get_mo_basis_mo_vgl(context, &(mo_vgl[0][0][0]), 5*chbrclf_walk_num*chbrclf_elec_num*chbrclf_mo_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); /* Set up determinant data */ @@ -1238,19 +1238,19 @@ for(k = 0; k < det_num_beta; ++k) mo_index_beta[k][i][j] = j + 1; rc = qmckl_set_determinant_type (context, typ); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_determinant_det_num_alpha (context, det_num_alpha); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_determinant_det_num_beta (context, det_num_beta); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_determinant_mo_index_alpha (context, &(mo_index_alpha[0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_determinant_mo_index_beta (context, &(mo_index_beta[0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); // Get slater-determinant @@ -1258,10 +1258,10 @@ double det_vgl_alpha[det_num_alpha][chbrclf_walk_num][5][chbrclf_elec_up_num][ch double det_vgl_beta[det_num_beta][chbrclf_walk_num][5][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; rc = qmckl_get_det_vgl_alpha(context, &(det_vgl_alpha[0][0][0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_get_det_vgl_beta(context, &(det_vgl_beta[0][0][0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); #+end_src @@ -2018,10 +2018,10 @@ double det_inv_matrix_alpha[det_num_alpha][chbrclf_walk_num][chbrclf_elec_up_num double det_inv_matrix_beta[det_num_beta][chbrclf_walk_num][chbrclf_elec_dn_num][chbrclf_elec_dn_num]; rc = qmckl_get_det_inv_matrix_alpha(context, &(det_inv_matrix_alpha[0][0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_get_det_inv_matrix_beta(context, &(det_inv_matrix_beta[0][0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); #+end_src @@ -2034,7 +2034,7 @@ assert (rc == QMCKL_SUCCESS); *** Test #+begin_src c :tangle (eval c_test) rc = qmckl_context_destroy(context); - assert (rc == QMCKL_SUCCESS); + qmckl_check(context, rc); return 0; } diff --git a/org/qmckl_electron.org b/org/qmckl_electron.org index 3aeca6f..a7164d7 100644 --- a/org/qmckl_electron.org +++ b/org/qmckl_electron.org @@ -79,38 +79,28 @@ int main() { The following data stored in the context: - | Variable | Type | Description | - |---------------------------+---------------+---------------------------------------| - | ~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 | - | ~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 | - | ~walker~ | ~qmckl_point~ | Current set of walkers | - | ~walker_old~ | ~qmckl_point~ | Previous set of walkers | + | Variable | Type | Description | + |---------------------------+--------------------+--------------------------------------| + | ~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 | + | ~provided~ | ~bool~ | If true, ~electron~ is valid | + | ~walker~ | ~qmckl_point~ | Current set of walkers | + | ~walker_old~ | ~qmckl_point~ | Previous set of walkers | Computed data: - | Variable | Type | Description | - |-------------------------------------+--------------------------------------+----------------------------------------------------------------------| - | ~ee_distance~ | ~double[walker.num][num][num]~ | Electron-electron distances | - | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~en_distance~ | ~double[walker.num][nucl_num][num]~ | Electron-nucleus distances | - | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | - | ~ee_distance_rescaled~ | ~double[walker.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[walker.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 | - | ~ee_potential~ | ~double[walker.num]~ | Electron-electron potential energy | - | ~ee_potential_date~ | ~uint64_t~ | Last modification date of the electron-electron potential | - | ~en_potential~ | ~double[walker.num]~ | Electron-nucleus potential energy | - | ~en_potential_date~ | ~int64_t~ | Date when the electron-nucleus potential energy was computed | - | ~en_distance_rescaled~ | ~double[walker.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[walker.num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | - | ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + | Variable | Type | Description | + |-------------------------------------+----------------------------------------+----------------------------------------------------------------------| + | ~ee_distance~ | ~double[walker.num][num][num]~ | Electron-electron distances | + | ~ee_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~en_distance~ | ~double[walker.num][nucl_num][num]~ | Electron-nucleus distances | + | ~en_distance_date~ | ~uint64_t~ | Last modification date of the electron-electron distances | + | ~ee_potential~ | ~double[walker.num]~ | Electron-electron potential energy | + | ~ee_potential_date~ | ~uint64_t~ | Last modification date of the electron-electron potential | + | ~en_potential~ | ~double[walker.num]~ | Electron-nucleus potential energy | + | ~en_potential_date~ | ~int64_t~ | Date when the electron-nucleus potential energy was computed | ** Data structure @@ -127,24 +117,14 @@ typedef struct qmckl_electron_struct { int64_t down_num; qmckl_walker walker; qmckl_walker walker_old; - double rescale_factor_kappa_ee; - double rescale_factor_kappa_en; - uint64_t ee_distance_date; - uint64_t en_distance_date; - uint64_t ee_potential_date; - uint64_t en_potential_date; - uint64_t ee_distance_rescaled_date; - uint64_t ee_distance_rescaled_deriv_e_date; - uint64_t en_distance_rescaled_date; - uint64_t en_distance_rescaled_deriv_e_date; + uint64_t ee_distance_date; + uint64_t en_distance_date; + uint64_t ee_potential_date; + uint64_t en_potential_date; double* ee_distance; double* en_distance; double* ee_potential; double* en_potential; - 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; @@ -174,10 +154,6 @@ qmckl_exit_code qmckl_init_electron(qmckl_context context) { ctx->electron.uninitialized = (1 << 1) - 1; - /* Default values */ - ctx->electron.rescale_factor_kappa_ee = 1.0; - ctx->electron.rescale_factor_kappa_en = 1.0; - return QMCKL_SUCCESS; } #+end_src @@ -201,6 +177,182 @@ bool qmckl_electron_provided(const qmckl_context context) { } #+end_src +** Initialization functions + + To set the data relative to the electrons in the context, the + following functions need to be called. When the data structure is + initialized, the internal ~coord_new~ and ~coord_old~ arrays are + both not allocated. + + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num); +qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const int64_t walk_num, const double* coord, const int64_t size_max); + #+end_src + + #+NAME:pre2 + #+begin_src c :exports none +if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + +qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + +if (mask != 0 && !(ctx->electron.uninitialized & mask)) { + return qmckl_failwith( context, + QMCKL_ALREADY_SET, + "qmckl_set_electron_*", + NULL); +} + #+end_src + + #+NAME:post + #+begin_src c :exports none +ctx->electron.uninitialized &= ~mask; +ctx->electron.provided = (ctx->electron.uninitialized == 0); + +return QMCKL_SUCCESS; + #+end_src + + To set the number of electrons, we give the number of up-spin and + down-spin electrons to the context and we set the number of walkers. + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_set_electron_num(qmckl_context context, + const int64_t up_num, + const int64_t down_num) { + int32_t mask = 1 << 0; + + <> + + if (up_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_electron_num", + "up_num <= 0"); + } + + if (down_num < 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_electron_num", + "down_num < 0"); + } + + ctx->electron.up_num = up_num; + ctx->electron.down_num = down_num; + ctx->electron.num = up_num + down_num; + + <> +} + #+end_src + + + #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes +interface + integer(c_int32_t) function qmckl_set_electron_num(context, alpha, beta) 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 :: alpha + integer (c_int64_t) , intent(in) , value :: beta + end function +end interface + #+end_src + + + The following function sets the electron coordinates of all the + walkers. When this is done, the pointers to the old and new sets + of coordinates are swapped, and the new coordinates are + overwritten. This can be done only when the data relative to + electrons have been set. + + ~size_max~ should be equal equal or geater than ~elec_num * + walker.num * 3~, to be symmetric with ~qmckl_get_electron_coord~. + + Important: changing the electron coordinates increments the date + in the context. + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_set_electron_coord(qmckl_context context, + const char transp, + const int64_t walk_num, + const double* coord, + const int64_t size_max) +{ + + int32_t mask = 0; + + <> + + if (transp != 'N' && transp != 'T') { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_electron_coord", + "transp should be 'N' or 'T'"); + } + + if (walk_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_electron_coord", + "walk_num <= 0"); + } + + if (coord == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_electron_coord", + "coord is a null pointer"); + } + + const int64_t elec_num = ctx->electron.num; + + if (elec_num == 0L) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_electron_coord", + "elec_num is not set"); + } + + /* Swap pointers */ + qmckl_walker tmp = ctx->electron.walker_old; + ctx->electron.walker_old = ctx->electron.walker; + ctx->electron.walker = tmp; + + memcpy(&(ctx->point), &(ctx->electron.walker.point), sizeof(qmckl_point_struct)); + + qmckl_exit_code rc; + rc = qmckl_set_point(context, transp, walk_num*elec_num, coord, size_max); + if (rc != QMCKL_SUCCESS) return rc; + + ctx->electron.walker.num = walk_num; + memcpy(&(ctx->electron.walker.point), &(ctx->point), sizeof(qmckl_point_struct)); + + return QMCKL_SUCCESS; + +} + #+end_src + + #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes +interface + integer(c_int32_t) function qmckl_set_electron_coord(context, transp, walk_num, coord, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + character , intent(in) , value :: transp + integer (c_int64_t) , intent(in) , value :: walk_num + double precision , intent(in) :: coord(*) + integer (c_int64_t) , intent(in) , value :: size_max + end function +end interface + #+end_src + ** Access functions Access functions return ~QMCKL_SUCCESS~ when the data has been @@ -210,13 +362,6 @@ bool qmckl_electron_provided(const qmckl_context context) { successfully, the variable pointed by the pointer given in argument contains the requested data. Otherwise, this variable is untouched. - #+NAME:post - #+begin_src c :exports none -if ( (ctx->electron.uninitialized & mask) != 0) { - return NULL; -} - #+end_src - *** Number of electrons #+begin_src c :comments org :tangle (eval h_func) :exports none @@ -343,59 +488,6 @@ qmckl_get_electron_walk_num (const qmckl_context context, int64_t* const walk_nu } #+end_src -*** Scaling factors Kappa - - #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_kappa_ee); -qmckl_exit_code qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const rescale_factor_kappa_en); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_get_electron_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_kappa_ee) { - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_INVALID_CONTEXT; - } - - if (rescale_factor_kappa_ee == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_get_electron_rescale_factor_ee", - "rescale_factor_kappa_ee is a null pointer"); - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - assert (ctx->electron.rescale_factor_kappa_ee > 0.0); - - *rescale_factor_kappa_ee = ctx->electron.rescale_factor_kappa_ee; - return QMCKL_SUCCESS; -} - - -qmckl_exit_code -qmckl_get_electron_rescale_factor_en (const qmckl_context context, double* const rescale_factor_kappa_en) { - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_INVALID_CONTEXT; - } - - if (rescale_factor_kappa_en == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_get_electron_rescale_factor_en", - "rescale_factor_kappa_en is a null pointer"); - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - assert (ctx->electron.rescale_factor_kappa_en > 0.0); - *rescale_factor_kappa_en = ctx->electron.rescale_factor_kappa_en; - return QMCKL_SUCCESS; -} - #+end_src *** Electron coordinates @@ -470,230 +562,6 @@ qmckl_get_electron_coord (const qmckl_context context, } #+end_src - -** Initialization functions - - To set the data relative to the electrons in the context, the - following functions need to be called. When the data structure is - initialized, the internal ~coord_new~ and ~coord_old~ arrays are - both not allocated. - - #+begin_src c :comments org :tangle (eval h_func) -qmckl_exit_code qmckl_set_electron_num (qmckl_context context, const int64_t up_num, const int64_t down_num); -qmckl_exit_code qmckl_set_electron_coord (qmckl_context context, const char transp, const int64_t walk_num, const double* coord, const int64_t size_max); - -qmckl_exit_code qmckl_set_electron_rescale_factor_ee (qmckl_context context, const double kappa_ee); -qmckl_exit_code qmckl_set_electron_rescale_factor_en (qmckl_context context, const double kappa_en); - #+end_src - - #+NAME:pre2 - #+begin_src c :exports none -if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - -qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - -if (mask != 0 && !(ctx->electron.uninitialized & mask)) { - return qmckl_failwith( context, - QMCKL_ALREADY_SET, - "qmckl_set_electron_*", - NULL); -} - #+end_src - - #+NAME:post2 - #+begin_src c :exports none -ctx->electron.uninitialized &= ~mask; -ctx->electron.provided = (ctx->electron.uninitialized == 0); - -return QMCKL_SUCCESS; - #+end_src - - To set the number of electrons, we give the number of up-spin and - down-spin electrons to the context and we set the number of walkers. - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_set_electron_num(qmckl_context context, - const int64_t up_num, - const int64_t down_num) { - int32_t mask = 1 << 0; - - <> - - if (up_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_electron_num", - "up_num <= 0"); - } - - if (down_num < 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_electron_num", - "down_num < 0"); - } - - ctx->electron.up_num = up_num; - ctx->electron.down_num = down_num; - ctx->electron.num = up_num + down_num; - - <> -} - #+end_src - - - Next we set the rescale parameter for the rescaled distance metric. - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_set_electron_rescale_factor_ee(qmckl_context context, - const double rescale_factor_kappa_ee) { - - int32_t mask = 0; // can be changed - - <> - - if (rescale_factor_kappa_ee <= 0.0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_electron_rescale_factor_ee", - "rescale_factor_kappa_ee <= 0.0"); - } - - ctx->electron.rescale_factor_kappa_ee = rescale_factor_kappa_ee; - - return QMCKL_SUCCESS; -} - -qmckl_exit_code -qmckl_set_electron_rescale_factor_en(qmckl_context context, - const double rescale_factor_kappa_en) { - - int32_t mask = 0; // can be changed - - <> - - if (rescale_factor_kappa_en <= 0.0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_electron_rescale_factor_en", - "rescale_factor_kappa_en <= 0.0"); - } - - ctx->electron.rescale_factor_kappa_en = rescale_factor_kappa_en; - - return QMCKL_SUCCESS; -} - #+end_src - - #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes -interface - integer(c_int32_t) function qmckl_set_electron_num(context, alpha, beta) 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 :: alpha - integer (c_int64_t) , intent(in) , value :: beta - end function -end interface - #+end_src - - - The following function sets the electron coordinates of all the - walkers. When this is done, the pointers to the old and new sets - of coordinates are swapped, and the new coordinates are - overwritten. This can be done only when the data relative to - electrons have been set. - - ~size_max~ should be equal equal or geater than ~elec_num * - walker.num * 3~, to be symmetric with ~qmckl_get_electron_coord~. - - Important: changing the electron coordinates increments the date - in the context. - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_set_electron_coord(qmckl_context context, - const char transp, - const int64_t walk_num, - const double* coord, - const int64_t size_max) -{ - - int32_t mask = 0; // coord can be changed - - <> - - if (transp != 'N' && transp != 'T') { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_electron_coord", - "transp should be 'N' or 'T'"); - } - - if (walk_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_electron_coord", - "walk_num <= 0"); - } - - if (coord == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_electron_coord", - "coord is a null pointer"); - } - - const int64_t elec_num = ctx->electron.num; - - if (elec_num == 0L) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_electron_coord", - "elec_num is not set"); - } - - /* Swap pointers */ - qmckl_walker tmp = ctx->electron.walker_old; - ctx->electron.walker_old = ctx->electron.walker; - ctx->electron.walker = tmp; - - memcpy(&(ctx->point), &(ctx->electron.walker.point), sizeof(qmckl_point_struct)); - - qmckl_exit_code rc; - rc = qmckl_set_point(context, transp, walk_num*elec_num, coord, size_max); - if (rc != QMCKL_SUCCESS) return rc; - - ctx->electron.walker.num = walk_num; - memcpy(&(ctx->electron.walker.point), &(ctx->point), sizeof(qmckl_point_struct)); - - return QMCKL_SUCCESS; - -} - #+end_src - - #+begin_src f90 :comments org :tangle (eval fh_func) :noweb yes -interface - integer(c_int32_t) function qmckl_set_electron_coord(context, transp, walk_num, coord, size_max) bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - - integer (c_int64_t) , intent(in) , value :: context - character , intent(in) , value :: transp - integer (c_int64_t) , intent(in) , value :: walk_num - double precision , intent(in) :: coord(*) - integer (c_int64_t) , intent(in) , value :: size_max - end function -end interface - #+end_src - ** Test #+begin_src python :results output :exports none @@ -707,14 +575,12 @@ 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 = 1.0; -double rescale_factor_kappa_en = 1.0; -double nucl_rescale_factor_kappa = 1.0; +int64_t nucl_num = chbrclf_nucl_num; +double* charge = chbrclf_charge; +double* nucl_coord = &(chbrclf_nucl_coord[0][0]); + double* elec_coord = &(chbrclf_elec_coord[0][0][0]); -int64_t nucl_num = chbrclf_nucl_num; -double* charge = chbrclf_charge; -double* nucl_coord = &(chbrclf_nucl_coord[0][0]); /* --- */ @@ -734,65 +600,41 @@ assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_electron_provided(context)); rc = qmckl_get_electron_up_num (context, &n); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(n == elec_up_num); rc = qmckl_get_electron_down_num (context, &n); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(n == elec_dn_num); rc = qmckl_get_electron_num (context, &n); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(n == elec_num); -double k_ee = 0.; -double k_en = 0.; -rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); -assert(rc == QMCKL_SUCCESS); -assert(k_ee == 1.0); - -rc = qmckl_get_electron_rescale_factor_en (context, &k_en); -assert(rc == QMCKL_SUCCESS); -assert(k_en == 1.0); - -rc = qmckl_set_electron_rescale_factor_en(context, rescale_factor_kappa_en); -assert(rc == QMCKL_SUCCESS); - -rc = qmckl_set_electron_rescale_factor_ee(context, rescale_factor_kappa_ee); -assert(rc == QMCKL_SUCCESS); - -rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); -assert(rc == QMCKL_SUCCESS); -assert(k_ee == rescale_factor_kappa_ee); - -rc = qmckl_get_electron_rescale_factor_en (context, &k_en); -assert(rc == QMCKL_SUCCESS); -assert(k_en == rescale_factor_kappa_en); - int64_t w = 0; rc = qmckl_get_electron_walk_num (context, &w); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(w == 0); assert(qmckl_electron_provided(context)); rc = qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*elec_num*3); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_get_electron_walk_num (context, &w); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(w == walk_num); double elec_coord2[walk_num*3*elec_num]; rc = qmckl_get_electron_coord (context, 'N', elec_coord2, walk_num*3*elec_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); for (int64_t i=0 ; i<3*elec_num*walk_num ; ++i) { printf("%f %f\n", elec_coord[i], elec_coord2[i]); assert( elec_coord[i] == elec_coord2[i] ); @@ -1061,491 +903,6 @@ assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12); #+end_src -** 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 -qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_get_electron_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled) -{ - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - - qmckl_exit_code rc; - - rc = qmckl_provide_ee_distance_rescaled(context); - if (rc != QMCKL_SUCCESS) return rc; - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num; - memcpy(distance_rescaled, ctx->electron.ee_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_ee_distance_rescaled(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(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); - - - /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.ee_distance_rescaled_date) { - - if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->electron.ee_distance_rescaled); - ctx->electron.ee_distance_rescaled = NULL; - } - - /* Allocate array */ - if (ctx->electron.ee_distance_rescaled == NULL) { - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = ctx->electron.num * ctx->electron.num * - ctx->electron.walker.num * sizeof(double); - double* ee_distance_rescaled = (double*) qmckl_malloc(context, mem_info); - - if (ee_distance_rescaled == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_ee_distance_rescaled", - NULL); - } - ctx->electron.ee_distance_rescaled = ee_distance_rescaled; - } - - qmckl_exit_code rc = - qmckl_compute_ee_distance_rescaled(context, - ctx->electron.num, - ctx->electron.rescale_factor_kappa_en, - ctx->electron.walker.num, - ctx->electron.walker.point.coord.data, - ctx->electron.ee_distance_rescaled); - 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 - :CRetType: qmckl_exit_code - :FRetType: qmckl_exit_code - :END: - - #+NAME: qmckl_ee_distance_rescaled_args - | Variable | Type | In/Out | Description | - |---------------------------+----------------------------------------+--------+--------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances | - | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | - | ~ee_distance~ | ~double[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, & - coord, ee_distance_rescaled) & - 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,walk_num,3) - double precision , intent(out) :: ee_distance_rescaled(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(context, 'T', 'T', elec_num, elec_num, & - coord(1,k,1), elec_num * walk_num, & - coord(1,k,1), elec_num * walk_num, & - ee_distance_rescaled(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_f - #+end_src - - #+begin_src c :tangle (eval h_private_func) :comments org :exports none -qmckl_exit_code qmckl_compute_ee_distance_rescaled ( - 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 ); - #+end_src - - #+CALL: generate_c_interface(table=qmckl_ee_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_ee_distance_rescaled & - (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_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 - 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(elec_num,elec_num,walk_num) - - integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_f - info = qmckl_compute_ee_distance_rescaled_f & - (context, elec_num, rescale_factor_kappa_ee, walk_num, coord, ee_distance_rescaled) - - end function qmckl_compute_ee_distance_rescaled - #+end_src - -*** Test - - #+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] : ", (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] : 0.9992169566605263 - : [1][0][0] : 0.9992169566605263 - : [0][0][1] : 0.0 - : [0][1][1] : 0.9985724058042633 - : [1][0][1] : 0.9985724058042633 - - #+begin_src c :tangle (eval c_test) -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*) context; - assert (ctx != NULL); - - size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walker.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*) context; - assert (ctx != NULL); - - - /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.ee_distance_rescaled_deriv_e_date) { - - if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->electron.ee_distance_rescaled_deriv_e); - ctx->electron.ee_distance_rescaled_deriv_e = NULL; - } - - /* 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.walker.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.walker.num, - ctx->electron.walker.point.coord.data, - 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 - | Variable | Type | In/Out | Description | - |---------------------------+-------------------------------------------+--------+-------------------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances | - | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | - | ~ee_distance_deriv_e~ | ~double[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,walk_num,3) - 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,k,1), elec_num*walk_num, & - coord(1,k,1), elec_num*walk_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. -//assert(ee_distance[0] == 0.); -// -//// (1,0,0) == (0,1,0) -//assert(ee_distance[1] == ee_distance[elec_num]); -// -//// value of (1,0,0) -//assert(fabs(ee_distance[1]-7.152322512964209) < 1.e-12); -// -//// (0,0,1) == 0. -//assert(ee_distance[elec_num*elec_num] == 0.); -// -//// (1,0,1) == (0,1,1) -//assert(ee_distance[elec_num*elec_num+1] == ee_distance[elec_num*elec_num+elec_num]); -// -//// value of (1,0,1) -//assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12); - - #+end_src - ** Electron-electron potential ~ee_potential~ is given by @@ -1612,7 +969,6 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context) if (ctx->electron.walker.num > ctx->electron.walker_old.num) { free(ctx->electron.ee_potential); - ctx->electron.ee_distance_rescaled_deriv_e = NULL; } /* Allocate array */ @@ -1655,13 +1011,13 @@ qmckl_exit_code qmckl_provide_ee_potential(qmckl_context context) :END: #+NAME: qmckl_ee_potential_args - | Variable | Type | In/Out | Description | - |----------------+----------------------------------------+--------+--------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron rescaled distances | - | ~ee_potential~ | ~double[walk_num]~ | out | Electron-electron potential | + | Variable | Type | In/Out | Description | + |----------------+----------------------------------------+--------+-----------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~elec_num~ | ~int64_t~ | in | Number of electrons | + | ~walk_num~ | ~int64_t~ | in | Number of walkers | + | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances | + | ~ee_potential~ | ~double[walk_num]~ | out | Electron-electron potential | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_ee_potential_f(context, elec_num, walk_num, & @@ -1749,7 +1105,7 @@ end function qmckl_compute_ee_potential_f double ee_potential[walk_num]; rc = qmckl_get_electron_ee_potential(context, &(ee_potential[0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); #+end_src ** Electron-nucleus distances @@ -2002,20 +1358,20 @@ assert(!qmckl_nucleus_provided(context)); assert(qmckl_electron_provided(context)); rc = qmckl_set_nucleus_num (context, nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_nucleus_charge (context, charge, nucl_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); rc = qmckl_set_nucleus_coord (context, 'T', nucl_coord, 3*nucl_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(qmckl_nucleus_provided(context)); double en_distance[walk_num][nucl_num][elec_num]; rc = qmckl_get_electron_en_distance(context, &(en_distance[0][0][0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); // (e,n,w) in Fortran notation // (1,1,1) @@ -2038,559 +1394,6 @@ assert(fabs(en_distance[1][0][1] - 3.1804527583077356) < 1.e-12); #+end_src -** Electron-nucleus rescaled distances - - ~en_distance_rescaled~ stores the matrix of the rescaled distances between - electrons and nuclei. - - \[ - 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 -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*) context; - assert (ctx != NULL); - - size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.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*) context; - assert (ctx != NULL); - - if (!(ctx->nucleus.provided)) { - return QMCKL_NOT_PROVIDED; - } - - /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.en_distance_rescaled_date) { - - if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->electron.en_distance_rescaled); - ctx->electron.en_distance_rescaled = NULL; - } - - /* 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.walker.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.rescale_factor_kappa_en, - ctx->electron.walker.num, - ctx->electron.walker.point.coord.data, - ctx->nucleus.coord.data, - 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 - | Variable | Type | In/Out | Description | - |---------------------------+----------------------------------------+--------+-----------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | - | ~rescale_factor_kappa_en~ | ~double~ | in | The factor for rescaled distances | - | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | - | ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | - | ~en_distance_rescaled~ | ~double[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, rescale_factor_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) :: rescale_factor_kappa_en - integer*8 , intent(in) :: walk_num - double precision , intent(in) :: elec_coord(elec_num,walk_num,3) - 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 (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(context, 'T', 'T', elec_num, nucl_num, & - elec_coord(1,k,1), elec_num*walk_num, & - nucl_coord, nucl_num, & - en_distance_rescaled(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_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 rescale_factor_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, rescale_factor_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 :: rescale_factor_kappa_en - integer (c_int64_t) , intent(in) , value :: walk_num - real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) - 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, rescale_factor_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 - -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 ]) -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] : ", (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] : 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) - -assert(qmckl_electron_provided(context)); -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][nucl_num][elec_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*) context; - assert (ctx != NULL); - - size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.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*) context; - assert (ctx != NULL); - - if (!(ctx->nucleus.provided)) { - return QMCKL_NOT_PROVIDED; - } - - /* Compute if necessary */ - if (ctx->electron.walker.point.date > ctx->electron.en_distance_rescaled_deriv_e_date) { - - if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->electron.en_distance_rescaled_deriv_e); - ctx->electron.en_distance_rescaled_deriv_e = NULL; - } - - /* 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.walker.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.walker.num, - ctx->electron.walker.point.coord.data, - ctx->nucleus.coord.data, - 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 - | Variable | Type | In/Out | Description | - |--------------------------------+-------------------------------------------+--------+---------------------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~elec_num~ | ~int64_t~ | in | Number of electrons | - | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | - | ~rescale_factor_kappa_en~ | ~double~ | in | The factor for rescaled distances | - | ~walk_num~ | ~int64_t~ | in | Number of walkers | - | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | - | ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | - | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][nucl_num][elec_num][4]~ | 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,walk_num,3) - double precision , intent(in) :: nucl_coord(nucl_num,3) - double precision , intent(out) :: en_distance_rescaled_deriv_e(4,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,k,1), elec_num*walk_num, & - nucl_coord, nucl_num, & - en_distance_rescaled_deriv_e(1,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,walk_num,3) - 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_rescale_factor (context, nucl_rescale_factor_kappa); -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 -//// (e,n,w) in Fortran notation -//// (1,1,1) -//assert(fabs(en_distance_rescaled[0][0][0] - 7.546738741619978) < 1.e-12); -// -//// (1,2,1) -//assert(fabs(en_distance_rescaled[0][1][0] - 8.77102435246984) < 1.e-12); -// -//// (2,1,1) -//assert(fabs(en_distance_rescaled[0][0][1] - 3.698922010513608) < 1.e-12); -// -//// (1,1,2) -//assert(fabs(en_distance_rescaled[1][0][0] - 5.824059436060509) < 1.e-12); -// -//// (1,2,2) -//assert(fabs(en_distance_rescaled[1][1][0] - 7.080482110317645) < 1.e-12); -// -//// (2,1,2) -//assert(fabs(en_distance_rescaled[1][0][1] - 3.1804527583077356) < 1.e-12); - - #+end_src - ** Electron-nucleus potential ~en_potential~ stores the ~en~ potential energy @@ -2804,7 +1607,7 @@ end function qmckl_compute_en_potential_f double en_potential[walk_num]; rc = qmckl_get_electron_en_potential(context, &(en_potential[0])); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); #+end_src ** Generate initial coordinates diff --git a/org/qmckl_error.org b/org/qmckl_error.org index 33829e4..a80db17 100644 --- a/org/qmckl_error.org +++ b/org/qmckl_error.org @@ -608,6 +608,62 @@ qmckl_last_error(qmckl_context context, char* buffer) { end interface #+end_src +* Helper functions for debugging + + The following function prints to ~stderr~ an error message is the return code is + not ~QMCKL_SUCCESS~. + + # Header + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code +qmckl_check(qmckl_context context, qmckl_exit_code rc); + #+end_src + + # Source + #+begin_src c :tangle (eval c) :exports none +#include + +qmckl_exit_code +qmckl_check(qmckl_context context, qmckl_exit_code rc) +{ + + char fname[QMCKL_MAX_FUN_LEN]; + char message[QMCKL_MAX_MSG_LEN]; + + if (rc != QMCKL_SUCCESS) { + fprintf(stderr, "===========\nQMCKL ERROR\n%s\n", qmckl_string_of_error(rc)); + qmckl_get_error(context, &rc, fname, message); + fprintf(stderr, "Function: %s\nMessage: %s\n===========\n", fname, message); + } + + return rc; +} + #+end_src + + It should be used as: + + #+begin_src c +rc = qmckl_check(context, + qmckl_...(context, ...) + ); +assert (rc == QMCKL_SUCCESS); + #+end_src + +** Fortran inteface + + #+begin_src f90 :tangle (eval fh_func) :exports none :noweb yes + interface + function qmckl_check (context, rc) bind(C, name='qmckl_check') + use, intrinsic :: iso_c_binding + import + implicit none + integer(qmckl_exit_code) :: qmckl_check + integer (c_int64_t) , intent(in), value :: context + integer(qmckl_exit_code), intent(in) :: rc + end function qmckl_check + end interface + #+end_src + * End of files :noexport: #+begin_src c :comments link :tangle (eval h_private_type) @@ -623,12 +679,13 @@ qmckl_last_error(qmckl_context context, char* buffer) { qmckl_exit_code exit_code; exit_code = 1; - assert (qmckl_set_error(context, exit_code, "qmckl_transpose", "Success") == QMCKL_SUCCESS); + assert (qmckl_set_error(context, exit_code, "my_function", "Message") == QMCKL_SUCCESS); assert (qmckl_get_error(context, &exit_code, function_name, message) == QMCKL_SUCCESS); assert (exit_code == 1); - assert (strcmp(function_name,"qmckl_transpose") == 0); - assert (strcmp(message,"Success") == 0); + + assert (strcmp(function_name,"my_function") == 0); + assert (strcmp(message,"Message") == 0); exit_code = qmckl_context_destroy(context); assert(exit_code == QMCKL_SUCCESS); diff --git a/org/qmckl_jastrow.org b/org/qmckl_jastrow.org index edcb3e0..f078f66 100644 --- a/org/qmckl_jastrow.org +++ b/org/qmckl_jastrow.org @@ -12,15 +12,15 @@ J(\mathbf{r},\mathbf{R}) = J_{\text{eN}}(\mathbf{r},\mathbf{R}) + J_{\text{ee}}(\mathbf{r}) + J_{\text{eeN}}(\mathbf{r},\mathbf{R}) \] - In the following, we us the notations $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ and + In the following, we use the notations $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ and $R_{i\alpha} = |\mathbf{r}_i - \mathbf{R}_\alpha|$. $J_{\text{eN}}$ contains electron-nucleus terms: \[ J_{\text{eN}}(\mathbf{r},\mathbf{R}) = \sum_{i=1}^{N_\text{elec}} \sum_{\alpha=1}^{N_\text{nucl}} - \frac{a_1\, f(R_{i\alpha})}{1+a_2\, f(R_{i\alpha})} + - \sum_{p=2}^{N_\text{ord}^a} a_{p+1}\, [f(R_{i\alpha})]^p - J_{eN}^\infty + \frac{a_{1,\alpha}\, g_\alpha(R_{i\alpha})}{1+a_{2,\alpha}\, g_\alpha(R_{i\alpha})} + + \sum_{p=2}^{N_\text{ord}^a} a_{p+1,\alpha}\, [g_\alpha(R_{i\alpha})]^p - J_{eN}^\infty \] $J_{\text{ee}}$ contains electron-electron terms: @@ -33,18 +33,18 @@ and $J_{\text{eeN}}$ contains electron-electron-Nucleus terms: - \[ - J_{\text{eeN}}(\mathbf{r},\mathbf{R}) = - \sum_{\alpha=1}^{N_{\text{nucl}}} - \sum_{i=1}^{N_{\text{elec}}} - \sum_{j=1}^{i-1} - \sum_{p=2}^{N_{\text{ord}}} - \sum_{k=0}^{p-1} - \sum_{l=0}^{p-k-2\delta_{k,0}} - c_{lkp\alpha} \left[ g({r}_{ij}) \right]^k - \left[ \left[ g({R}_{i\alpha}) \right]^l + \left[ g({R}_{j\alpha}) \right]^l \right] - \left[ g({R}_{i\,\alpha}) \, g({R}_{j\alpha}) \right]^{(p-k-l)/2} - \] + \[ + J_{\text{eeN}}(\mathbf{r},\mathbf{R}) = + \sum_{\alpha=1}^{N_{\text{nucl}}} + \sum_{i=1}^{N_{\text{elec}}} + \sum_{j=1}^{i-1} + \sum_{p=2}^{N_{\text{ord}}} + \sum_{k=0}^{p-1} + \sum_{l=0}^{p-k-2\delta_{k,0}} + c_{lkp\alpha} \left[ f({r}_{ij}) \right]^k + \left[ \left[ g_\alpha({R}_{i\alpha}) \right]^l + \left[ g_\alpha({R}_{j\alpha}) \right]^l \right] + \left[ g_\alpha({R}_{i\,\alpha}) \, g_\alpha({R}_{j\alpha}) \right]^{(p-k-l)/2} + \] $c_{lkp\alpha}$ are non-zero only when $p-k-l$ is even. @@ -52,7 +52,7 @@ \[ f(r) = \frac{1-e^{-\kappa\, r}}{\kappa} \text{ and } - g(r) = e^{-\kappa\, r}. + g_\alpha(r) = e^{-\kappa_\alpha\, r}. \] The terms $J_{\text{ee}}^\infty$ and $J_{\text{eN}}^\infty$ are shifts to ensure that @@ -137,6 +137,8 @@ int main() { | Variable | Type | In/Out | Description | |---------------------------+---------------------------------------+--------+-------------------------------------------------------------------| | ~uninitialized~ | ~int32_t~ | in | Keeps bits set for uninitialized data | + | ~rescale_factor_ee~ | ~double~ | in | The distance scaling factor | + | ~rescale_factor_en~ | ~double[nucl_num]~ | in | The distance scaling factor | | ~aord_num~ | ~int64_t~ | in | The number of a coeffecients | | ~bord_num~ | ~int64_t~ | in | The number of b coeffecients | | ~cord_num~ | ~int64_t~ | in | The number of c coeffecients | @@ -145,39 +147,47 @@ int main() { | ~aord_vector~ | ~double[aord_num + 1][type_nucl_num]~ | in | Order of a polynomial coefficients | | ~bord_vector~ | ~double[bord_num + 1]~ | in | Order of b polynomial coefficients | | ~cord_vector~ | ~double[cord_num][type_nucl_num]~ | in | Order of c polynomial coefficients | - | ~factor_ee~ | ~double[walker.num]~ | out | Jastrow factor: electron-electron part | + | ~factor_ee~ | ~double[walker.num]~ | out | Jastrow factor: electron-electron part | | ~factor_ee_date~ | ~uint64_t~ | out | Jastrow factor: electron-electron part | - | ~factor_en~ | ~double[walker.num]~ | out | Jastrow factor: electron-nucleus part | + | ~factor_en~ | ~double[walker.num]~ | out | Jastrow factor: electron-nucleus part | | ~factor_en_date~ | ~uint64_t~ | out | Jastrow factor: electron-nucleus part | - | ~factor_een~ | ~double[walker.num]~ | out | Jastrow factor: electron-electron-nucleus part | + | ~factor_een~ | ~double[walker.num]~ | out | Jastrow factor: electron-electron-nucleus part | | ~factor_een_date~ | ~uint64_t~ | out | Jastrow factor: electron-electron-nucleus part | - | ~factor_ee_deriv_e~ | ~double[4][nelec][walker.num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | + | ~factor_ee_deriv_e~ | ~double[4][nelec][walker.num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | | ~factor_ee_deriv_e_date~ | ~uint64_t~ | out | Keep track of the date for the derivative | - | ~factor_en_deriv_e~ | ~double[4][nelec][walker.num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | + | ~factor_en_deriv_e~ | ~double[4][nelec][walker.num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | | ~factor_en_deriv_e_date~ | ~uint64_t~ | out | Keep track of the date for the en derivative | - | ~factor_een_deriv_e~ | ~double[4][nelec][walker.num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | + | ~factor_een_deriv_e~ | ~double[4][nelec][walker.num]~ | out | Derivative of the Jastrow factor: electron-electron-nucleus part | | ~factor_een_deriv_e_date~ | ~uint64_t~ | out | Keep track of the date for the een derivative | computed data: - | Variable | Type | In/Out | - |----------------------------+-----------------------------------------------------------------+-------------------------------------------------| - | ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | - | ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | - | ~asymp_jasb~ | ~double[2]~ | Asymptotic component | - | ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | - | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | - | ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | - | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | - | ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | - | ~tmp_c~ | ~double[walker.num][cord_num][cord_num+1][nucl_num][elec_num]~ | vector of non-zero coefficients | - | ~dtmp_c~ | ~double[walker.num][elec_num][4][nucl_num][cord_num+1][cord_num]~ | vector of non-zero coefficients | - | ~een_rescaled_n~ | ~double[walker.num][cord_num+1][nucl_num][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | - | ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | - | ~een_rescaled_e_deriv_e~ | ~double[walker.num][cord_num+1][elec_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | - | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | - | ~een_rescaled_n_deriv_e~ | ~double[walker.num][cord_num+1][nucl_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | - | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | + | Variable | Type | In/Out | + |-------------------------------------+-------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------| + | ~dim_cord_vect~ | ~int64_t~ | Number of unique C coefficients | + | ~dim_cord_vect_date~ | ~uint64_t~ | Number of unique C coefficients | + | ~asymp_jasb~ | ~double[2]~ | Asymptotic component | + | ~asymp_jasb_date~ | ~uint64_t~ | Asymptotic component | + | ~cord_vect_full~ | ~double[dim_cord_vect][nucl_num]~ | vector of non-zero coefficients | + | ~cord_vect_full_date~ | ~uint64_t~ | Keep track of changes here | + | ~lkpm_combined_index~ | ~int64_t[4][dim_cord_vect]~ | Transform l,k,p, and m into consecutive indices | + | ~lkpm_combined_index_date~ | ~uint64_t~ | Transform l,k,p, and m into consecutive indices | + | ~tmp_c~ | ~double[walker.num][cord_num][cord_num+1][nucl_num][elec_num]~ | vector of non-zero coefficients | + | ~dtmp_c~ | ~double[walker.num][elec_num][4][nucl_num][cord_num+1][cord_num]~ | vector of non-zero coefficients | + | ~ee_distance_rescaled~ | ~double[walker.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[walker.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[walker.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[walker.num][4][nucl_num][num]~ | Electron-electron rescaled distances derivatives | + | ~en_distance_rescaled_deriv_e_date~ | ~uint64_t~ | Last modification date of the electron-electron distance derivatives | + | ~een_rescaled_n~ | ~double[walker.num][cord_num+1][nucl_num][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord | + | ~een_rescaled_n_date~ | ~uint64_t~ | Keep track of the date of creation | + | ~een_rescaled_e_deriv_e~ | ~double[walker.num][cord_num+1][elec_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | + | ~een_rescaled_e_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | + | ~een_rescaled_n_deriv_e~ | ~double[walker.num][cord_num+1][nucl_num][4][elec_num]~ | The electron-electron rescaled distances raised to the powers defined by cord derivatives wrt electrons | + | ~een_rescaled_n_deriv_e_date~ | ~uint64_t~ | Keep track of the date of creation | #+NAME: jastrow_data #+BEGIN_SRC python :results none :exports none @@ -334,11 +344,11 @@ kappa_inv = 1.0/kappa #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_jastrow_struct{ - int32_t uninitialized; - int64_t aord_num; - int64_t bord_num; - int64_t cord_num; - int64_t type_nucl_num; + int32_t uninitialized; + int64_t aord_num; + int64_t bord_num; + int64_t cord_num; + int64_t type_nucl_num; uint64_t asymp_jasb_date; uint64_t tmp_c_date; uint64_t dtmp_c_date; @@ -348,35 +358,45 @@ typedef struct qmckl_jastrow_struct{ uint64_t factor_ee_deriv_e_date; uint64_t factor_en_deriv_e_date; uint64_t factor_een_deriv_e_date; - int64_t* type_nucl_vector; - double * aord_vector; - double * bord_vector; - double * cord_vector; - double * asymp_jasb; - double * factor_ee; - double * factor_en; - double * factor_een; - double * factor_ee_deriv_e; - double * factor_en_deriv_e; - double * factor_een_deriv_e; - int64_t dim_cord_vect; + double rescale_factor_ee; + double* rescale_factor_en; + int64_t* type_nucl_vector; + double * aord_vector; + double * bord_vector; + double * cord_vector; + double * asymp_jasb; + double * factor_ee; + double * factor_en; + double * factor_een; + double * factor_ee_deriv_e; + double * factor_en_deriv_e; + double * factor_een_deriv_e; + int64_t dim_cord_vect; uint64_t dim_cord_vect_date; - double * cord_vect_full; + double * cord_vect_full; uint64_t cord_vect_full_date; - int64_t* lkpm_combined_index; + int64_t* lkpm_combined_index; uint64_t lkpm_combined_index_date; - double * tmp_c; - double * dtmp_c; - double * een_rescaled_e; - double * een_rescaled_n; + double * tmp_c; + double * dtmp_c; + uint64_t ee_distance_rescaled_date; + uint64_t ee_distance_rescaled_deriv_e_date; + uint64_t en_distance_rescaled_date; + uint64_t en_distance_rescaled_deriv_e_date; + double* ee_distance_rescaled; + double* ee_distance_rescaled_deriv_e; + double* en_distance_rescaled; + double* en_distance_rescaled_deriv_e; + double * een_rescaled_e; + double * een_rescaled_n; uint64_t een_rescaled_e_date; uint64_t een_rescaled_n_date; - double * een_rescaled_e_deriv_e; - double * een_rescaled_n_deriv_e; + double * een_rescaled_e_deriv_e; + double * een_rescaled_n_deriv_e; uint64_t een_rescaled_e_deriv_e_date; uint64_t een_rescaled_n_deriv_e_date; - bool provided; - char * type; + bool provided; + char * type; #ifdef HAVE_HPC bool gpu_offload; @@ -406,7 +426,7 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); - ctx->jastrow.uninitialized = (1 << 6) - 1; + ctx->jastrow.uninitialized = (1 << 8) - 1; /* Default values */ @@ -414,6 +434,494 @@ qmckl_exit_code qmckl_init_jastrow(qmckl_context context) { } #+end_src +** Initialization functions + + To prepare for the Jastrow and its derivative, all the following functions need to be + called. + + #+begin_src c :comments org :tangle (eval h_func) +qmckl_exit_code qmckl_set_jastrow_rescale_factor_ee (qmckl_context context, const double kappa_ee); +qmckl_exit_code qmckl_set_jastrow_rescale_factor_en (qmckl_context context, const double* kappa_en, const int64_t size_max); +qmckl_exit_code qmckl_set_jastrow_ord_num (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num); +qmckl_exit_code qmckl_set_jastrow_type_nucl_num (qmckl_context context, const int64_t type_nucl_num); +qmckl_exit_code qmckl_set_jastrow_type_nucl_vector (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num); +qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector, const int64_t size_max); +qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector, const int64_t size_max); +qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector, const int64_t size_max); + #+end_src + + #+NAME:pre2 + #+begin_src c :exports none +if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + +qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + +if (mask != 0 && !(ctx->jastrow.uninitialized & mask)) { + printf("%d %d\n", mask, ctx->jastrow.uninitialized ); + return qmckl_failwith( context, + QMCKL_ALREADY_SET, + "qmckl_set_jastrow_*", + NULL); + } + #+end_src + + #+NAME:post2 + #+begin_src c :exports none +ctx->jastrow.uninitialized &= ~mask; +ctx->jastrow.provided = (ctx->jastrow.uninitialized == 0); +if (ctx->jastrow.provided) { + qmckl_exit_code rc_ = qmckl_finalize_jastrow(context); + if (rc_ != QMCKL_SUCCESS) return rc_; + } + +return QMCKL_SUCCESS; + #+end_src + + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_set_jastrow_ord_num(qmckl_context context, + const int64_t aord_num, + const int64_t bord_num, + const int64_t cord_num) +{ + + int32_t mask = 1 << 0; + +<> + + if (aord_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_ord_num", + "aord_num <= 0"); + } + + if (bord_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_ord_num", + "bord_num <= 0"); + } + + if (cord_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_ord_num", + "cord_num <= 0"); + } + + ctx->jastrow.aord_num = aord_num; + ctx->jastrow.bord_num = bord_num; + ctx->jastrow.cord_num = cord_num; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) +{ + int32_t mask = 1 << 1; + +<> + + if (type_nucl_num <= 0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_type_nucl_num", + "type_nucl_num < 0"); + } + + ctx->jastrow.type_nucl_num = type_nucl_num; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_type_nucl_vector(qmckl_context context, + int64_t const * type_nucl_vector, + const int64_t nucl_num) +{ + + int32_t mask = 1 << 2; + +<> + + int64_t type_nucl_num; + qmckl_exit_code rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (type_nucl_num == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_type_nucl_vector", + "type_nucl_num is not set"); + } + + if (type_nucl_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_type_nucl_vector", + "type_nucl_vector = NULL"); + } + + if (ctx->jastrow.type_nucl_vector != NULL) { + rc = qmckl_free(context, ctx->jastrow.type_nucl_vector); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_type_nucl_vector", + "Unable to free ctx->jastrow.type_nucl_vector"); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = nucl_num * sizeof(int64_t); + int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_type_nucl_vector", + NULL); + } + + memcpy(new_array, type_nucl_vector, mem_info.size); + + ctx->jastrow.type_nucl_vector = new_array; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_aord_vector(qmckl_context context, + double const * aord_vector, + const int64_t size_max) +{ + int32_t mask = 1 << 3; + +<> + + int64_t aord_num; + qmckl_exit_code rc = qmckl_get_jastrow_aord_num(context, &aord_num); + if (rc != QMCKL_SUCCESS) return rc; + + int64_t type_nucl_num; + rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (aord_num == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_coefficient", + "aord_num is not set"); + } + + if (aord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_aord_vector", + "aord_vector = NULL"); + } + + if (ctx->jastrow.aord_vector != NULL) { + rc = qmckl_free(context, ctx->jastrow.aord_vector); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_jastrow_aord_vector", + "Unable to free ctx->jastrow.aord_vector"); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double); + + if ((size_t) size_max < mem_info.size/sizeof(double)) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_aord_vector", + "Array too small. Expected (aord_num+1)*type_nucl_num"); + } + + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_coefficient", + NULL); + } + + memcpy(new_array, aord_vector, mem_info.size); + + ctx->jastrow.aord_vector = new_array; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_bord_vector(qmckl_context context, + double const * bord_vector, + const int64_t size_max) +{ + int32_t mask = 1 << 4; + +<> + + int64_t bord_num; + qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (bord_num == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_coefficient", + "bord_num is not set"); + } + + if (bord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_bord_vector", + "bord_vector = NULL"); + } + + if (ctx->jastrow.bord_vector != NULL) { + rc = qmckl_free(context, ctx->jastrow.bord_vector); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_jastrow_bord_vector", + "Unable to free ctx->jastrow.bord_vector"); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = (bord_num + 1) * sizeof(double); + + if ((size_t) size_max < mem_info.size/sizeof(double)) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_bord_vector", + "Array too small. Expected (bord_num+1)"); + } + + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_coefficient", + NULL); + } + + memcpy(new_array, bord_vector, mem_info.size); + + ctx->jastrow.bord_vector = new_array; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_cord_vector(qmckl_context context, + double const * cord_vector, + const int64_t size_max) +{ + int32_t mask = 1 << 5; + +<> + + qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context); + if (rc != QMCKL_SUCCESS) return rc; + + int64_t dim_cord_vect; + rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect); + if (rc != QMCKL_SUCCESS) return rc; + + int64_t type_nucl_num; + rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); + if (rc != QMCKL_SUCCESS) return rc; + + if (dim_cord_vect == 0) { + return qmckl_failwith( context, + QMCKL_FAILURE, + "qmckl_set_jastrow_coefficient", + "dim_cord_vect is not set"); + } + + if (cord_vector == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_cord_vector", + "cord_vector = NULL"); + } + + if (ctx->jastrow.cord_vector != NULL) { + rc = qmckl_free(context, ctx->jastrow.cord_vector); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_set_jastrow_cord_vector", + "Unable to free ctx->jastrow.cord_vector"); + } + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = dim_cord_vect * type_nucl_num * sizeof(double); + + if ((size_t) size_max < mem_info.size/sizeof(double)) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_cord_vector", + "Array too small. Expected dim_cord_vect * type_nucl_num"); + } + + double* new_array = (double*) qmckl_malloc(context, mem_info); + + if(new_array == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_set_jastrow_coefficient", + NULL); + } + + memcpy(new_array, cord_vector, mem_info.size); + + ctx->jastrow.cord_vector = new_array; + + <> +} + +qmckl_exit_code +qmckl_set_jastrow_rescale_factor_ee(qmckl_context context, + const double rescale_factor_ee) { + + int32_t mask = 1 << 6; + + <> + + if (rescale_factor_ee <= 0.0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_rescale_factor_ee", + "rescale_factor_ee <= 0.0"); + } + + ctx->jastrow.rescale_factor_ee = rescale_factor_ee; + + <> +} + + +qmckl_exit_code +qmckl_set_jastrow_rescale_factor_en(qmckl_context context, + const double* rescale_factor_en, + const int64_t size_max) { + + int32_t mask = 1 << 7; + + <> + + if (rescale_factor_en == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_rescale_factor_en", + "Null pointer"); + } + + if (size_max < ctx->nucleus.num) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_rescale_factor_en", + "Array too small"); + } + + + if (ctx->jastrow.rescale_factor_en != NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_set_jastrow_rescale_factor_en", + "Already set"); + } + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->nucleus.num * sizeof(double); + ctx->jastrow.rescale_factor_en = (double*) qmckl_malloc(context, mem_info); + + for (int64_t i=0 ; inucleus.num ; ++i) { + if (rescale_factor_en[i] <= 0.0) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_set_jastrow_rescale_factor_en", + "rescale_factor_en <= 0.0"); + } + ctx->jastrow.rescale_factor_en[i] = rescale_factor_en[i]; + } + + <> +} + #+end_src + + When the required information is completely entered, other data structures are + computed to accelerate the calculations. The intermediates factors + are precontracted using BLAS LEVEL 3 operations for an optimal flop count. + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + /* ----------------------------------- */ + /* Check for the necessary information */ + /* ----------------------------------- */ + + /* Check for the electron data + 1. elec_num + 2. ee_distances_rescaled + ,*/ + if (!(ctx->electron.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_electron", + NULL); + } + + /* Check for the nucleus data + 1. nucl_num + 2. en_distances_rescaled + ,*/ + if (!(ctx->nucleus.provided)) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_nucleus", + NULL); + } + + /* Decide if the Jastrow if offloaded on GPU or not */ +#if defined(HAVE_HPC) && (defined(HAVE_CUBLAS_OFFLOAD) || defined(HAVE_OPENACC_OFFLOAD) || defined(HAVE_OPENMP_OFFLOAD)) + ctx->jastrow.gpu_offload = true; // ctx->electron.num > 100; +#endif + + qmckl_exit_code rc = QMCKL_SUCCESS; + return rc; + + +} + #+end_src + ** Access functions #+begin_src c :comments org :tangle (eval h_func) :exports none @@ -425,7 +933,9 @@ qmckl_exit_code qmckl_get_jastrow_type_nucl_vector (qmckl_context context, int qmckl_exit_code qmckl_get_jastrow_aord_vector (qmckl_context context, double * const aord_vector, const int64_t size_max); qmckl_exit_code qmckl_get_jastrow_bord_vector (qmckl_context context, double * const bord_vector, const int64_t size_max); qmckl_exit_code qmckl_get_jastrow_cord_vector (qmckl_context context, double * const cord_vector, const int64_t size_max); - #+end_src +qmckl_exit_code qmckl_get_jastrow_rescale_factor_ee (const qmckl_context context, double* const rescale_factor_ee); +qmckl_exit_code qmckl_get_jastrow_rescale_factor_en (const qmckl_context context, double* const rescale_factor_en, const int64_t size_max); + #+end_src Along with these core functions, calculation of the jastrow factor requires the following additional information to be set: @@ -452,13 +962,6 @@ bool qmckl_jastrow_provided(const qmckl_context context) { } #+end_src - #+NAME:post - #+begin_src c :exports none -if ( (ctx->jastrow.uninitialized & mask) != 0) { - return NULL; -} - #+end_src - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code qmckl_get_jastrow_aord_num (const qmckl_context context, int64_t* const aord_num) { @@ -722,425 +1225,75 @@ qmckl_get_jastrow_cord_vector (const qmckl_context context, return QMCKL_SUCCESS; } - #+end_src - -** Initialization functions - - To prepare for the Jastrow and its derivative, all the following functions need to be - called. - - #+begin_src c :comments org :tangle (eval h_func) -qmckl_exit_code qmckl_set_jastrow_ord_num (qmckl_context context, const int64_t aord_num, const int64_t bord_num, const int64_t cord_num); -qmckl_exit_code qmckl_set_jastrow_type_nucl_num (qmckl_context context, const int64_t type_nucl_num); -qmckl_exit_code qmckl_set_jastrow_type_nucl_vector (qmckl_context context, const int64_t* type_nucl_vector, const int64_t nucl_num); -qmckl_exit_code qmckl_set_jastrow_aord_vector (qmckl_context context, const double * aord_vector, const int64_t size_max); -qmckl_exit_code qmckl_set_jastrow_bord_vector (qmckl_context context, const double * bord_vector, const int64_t size_max); -qmckl_exit_code qmckl_set_jastrow_cord_vector (qmckl_context context, const double * cord_vector, const int64_t size_max); - #+end_src - - #+NAME:pre2 - #+begin_src c :exports none -if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_NULL_CONTEXT; - } - -qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - -if (mask != 0 && !(ctx->jastrow.uninitialized & mask)) { - printf("%d %d\n", mask, ctx->jastrow.uninitialized ); - return qmckl_failwith( context, - QMCKL_ALREADY_SET, - "qmckl_set_jastrow_*", - NULL); - } - #+end_src - - #+NAME:post2 - #+begin_src c :exports none -ctx->jastrow.uninitialized &= ~mask; -ctx->jastrow.provided = (ctx->jastrow.uninitialized == 0); -if (ctx->jastrow.provided) { - qmckl_exit_code rc_ = qmckl_finalize_jastrow(context); - if (rc_ != QMCKL_SUCCESS) return rc_; - } - -return QMCKL_SUCCESS; - #+end_src - - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none qmckl_exit_code -qmckl_set_jastrow_ord_num(qmckl_context context, - const int64_t aord_num, - const int64_t bord_num, - const int64_t cord_num) -{ - - int32_t mask = 1 << 0; - -<> - - if (aord_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_ord_num", - "aord_num <= 0"); - } - - if (bord_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_ord_num", - "bord_num <= 0"); - } - - if (cord_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_ord_num", - "cord_num <= 0"); - } - - ctx->jastrow.aord_num = aord_num; - ctx->jastrow.bord_num = bord_num; - ctx->jastrow.cord_num = cord_num; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_type_nucl_num(qmckl_context context, const int64_t type_nucl_num) -{ - int32_t mask = 1 << 1; - -<> - - if (type_nucl_num <= 0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_type_nucl_num", - "type_nucl_num < 0"); - } - - ctx->jastrow.type_nucl_num = type_nucl_num; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_type_nucl_vector(qmckl_context context, - int64_t const * type_nucl_vector, - const int64_t nucl_num) -{ - - int32_t mask = 1 << 2; - -<> - - int64_t type_nucl_num; - qmckl_exit_code rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); - if (rc != QMCKL_SUCCESS) return rc; - - if (type_nucl_num == 0) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_jastrow_type_nucl_vector", - "type_nucl_num is not set"); - } - - if (type_nucl_vector == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_type_nucl_vector", - "type_nucl_vector = NULL"); - } - - if (ctx->jastrow.type_nucl_vector != NULL) { - rc = qmckl_free(context, ctx->jastrow.type_nucl_vector); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_set_type_nucl_vector", - NULL); - } - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = nucl_num * sizeof(int64_t); - int64_t* new_array = (int64_t*) qmckl_malloc(context, mem_info); - - if(new_array == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_type_nucl_vector", - NULL); - } - - memcpy(new_array, type_nucl_vector, mem_info.size); - - ctx->jastrow.type_nucl_vector = new_array; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_aord_vector(qmckl_context context, - double const * aord_vector, - const int64_t size_max) -{ - int32_t mask = 1 << 3; - -<> - - int64_t aord_num; - qmckl_exit_code rc = qmckl_get_jastrow_aord_num(context, &aord_num); - if (rc != QMCKL_SUCCESS) return rc; - - int64_t type_nucl_num; - rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); - if (rc != QMCKL_SUCCESS) return rc; - - if (aord_num == 0) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_jastrow_coefficient", - "aord_num is not set"); - } - - if (aord_vector == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_aord_vector", - "aord_vector = NULL"); - } - - if (ctx->jastrow.aord_vector != NULL) { - rc = qmckl_free(context, ctx->jastrow.aord_vector); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_set_ord_vector", - NULL); - } - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = (aord_num + 1) * type_nucl_num * sizeof(double); - - if ((size_t) size_max < mem_info.size/sizeof(double)) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_jastrow_aord_vector", - "Array too small. Expected (aord_num+1)*type_nucl_num"); - } - - double* new_array = (double*) qmckl_malloc(context, mem_info); - - if(new_array == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_coefficient", - NULL); - } - - memcpy(new_array, aord_vector, mem_info.size); - - ctx->jastrow.aord_vector = new_array; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_bord_vector(qmckl_context context, - double const * bord_vector, - const int64_t size_max) -{ - int32_t mask = 1 << 4; - -<> - - int64_t bord_num; - qmckl_exit_code rc = qmckl_get_jastrow_bord_num(context, &bord_num); - if (rc != QMCKL_SUCCESS) return rc; - - if (bord_num == 0) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_jastrow_coefficient", - "bord_num is not set"); - } - - if (bord_vector == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_bord_vector", - "bord_vector = NULL"); - } - - if (ctx->jastrow.bord_vector != NULL) { - rc = qmckl_free(context, ctx->jastrow.bord_vector); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_set_ord_vector", - NULL); - } - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = (bord_num + 1) * sizeof(double); - - if ((size_t) size_max < mem_info.size/sizeof(double)) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_jastrow_bord_vector", - "Array too small. Expected (bord_num+1)"); - } - - double* new_array = (double*) qmckl_malloc(context, mem_info); - - if(new_array == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_coefficient", - NULL); - } - - memcpy(new_array, bord_vector, mem_info.size); - - ctx->jastrow.bord_vector = new_array; - - <> -} - - -qmckl_exit_code -qmckl_set_jastrow_cord_vector(qmckl_context context, - double const * cord_vector, - const int64_t size_max) -{ - int32_t mask = 1 << 5; - -<> - - qmckl_exit_code rc = qmckl_provide_dim_cord_vect(context); - if (rc != QMCKL_SUCCESS) return rc; - - int64_t dim_cord_vect; - rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect); - if (rc != QMCKL_SUCCESS) return rc; - - int64_t type_nucl_num; - rc = qmckl_get_jastrow_type_nucl_num(context, &type_nucl_num); - if (rc != QMCKL_SUCCESS) return rc; - - if (dim_cord_vect == 0) { - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_set_jastrow_coefficient", - "dim_cord_vect is not set"); - } - - if (cord_vector == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_jastrow_cord_vector", - "cord_vector = NULL"); - } - - if (ctx->jastrow.cord_vector != NULL) { - rc = qmckl_free(context, ctx->jastrow.cord_vector); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_set_cord_vector", - NULL); - } - } - - qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; - mem_info.size = dim_cord_vect * type_nucl_num * sizeof(double); - - if ((size_t) size_max < mem_info.size/sizeof(double)) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_set_jastrow_cord_vector", - "Array too small. Expected dim_cord_vect * type_nucl_num"); - } - - double* new_array = (double*) qmckl_malloc(context, mem_info); - - if(new_array == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_set_jastrow_coefficient", - NULL); - } - - memcpy(new_array, cord_vector, mem_info.size); - - ctx->jastrow.cord_vector = new_array; - - <> -} - - #+end_src - - When the required information is completely entered, other data structures are - computed to accelerate the calculations. The intermediates factors - are precontracted using BLAS LEVEL 3 operations for an optimal flop count. - - #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none -qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_finalize_jastrow(qmckl_context context) { +qmckl_get_jastrow_rescale_factor_ee (const qmckl_context context, + double* const rescale_factor_ee) { if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { return QMCKL_INVALID_CONTEXT; } + if (rescale_factor_ee == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_rescale_factor_ee", + "rescale_factor_ee is a null pointer"); + } + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); - /* ----------------------------------- */ - /* Check for the necessary information */ - /* ----------------------------------- */ + int32_t mask = 1 << 6; - /* Check for the electron data - 1. elec_num - 2. ee_distances_rescaled - ,*/ - if (!(ctx->electron.provided)) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_electron", - NULL); + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; } + assert (ctx->jastrow.rescale_factor_ee > 0.0); - /* Check for the nucleus data - 1. nucl_num - 2. en_distances_rescaled - ,*/ - if (!(ctx->nucleus.provided)) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_nucleus", - NULL); - } - - /* Decide if the Jastrow if offloaded on GPU or not */ -#if defined(HAVE_HPC) && (defined(HAVE_CUBLAS_OFFLOAD) || defined(HAVE_OPENACC_OFFLOAD) || defined(HAVE_OPENMP_OFFLOAD)) - ctx->jastrow.gpu_offload = true; // ctx->electron.num > 100; -#endif - - qmckl_exit_code rc = QMCKL_SUCCESS; - return rc; - - + ,*rescale_factor_ee = ctx->jastrow.rescale_factor_ee; + return QMCKL_SUCCESS; } - #+end_src + + +qmckl_exit_code +qmckl_get_jastrow_rescale_factor_en (const qmckl_context context, + double* const rescale_factor_en, + const int64_t size_max) { + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_INVALID_CONTEXT; + } + + if (rescale_factor_en == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_jastrow_rescale_factor_en", + "rescale_factor_en is a null pointer"); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + int32_t mask = 1 << 7; + + if ( (ctx->jastrow.uninitialized & mask) != 0) { + return QMCKL_NOT_PROVIDED; + } + + if (size_max < ctx->nucleus.num) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_rescale_factor_en", + "Array to small"); + } + + assert(ctx->jastrow.rescale_factor_en != NULL); + for (int64_t i=0 ; inucleus.num ; ++i) { + rescale_factor_en[i] = ctx->jastrow.rescale_factor_en[i]; + } + + return QMCKL_SUCCESS; +} + #+end_src ** Test #+begin_src c :tangle (eval c_test) @@ -1149,13 +1302,12 @@ 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; -double rescale_factor_kappa_ee = 1.0; -double rescale_factor_kappa_en = 1.0; -double nucl_rescale_factor_kappa = 1.0; +int64_t nucl_num = n2_nucl_num; +double rescale_factor_ee = 1.0; +double rescale_factor_en[2] = { 1.0, 1.0 }; double* elec_coord = &(n2_elec_coord[0][0][0]); const double* nucl_charge = n2_charge; -int64_t nucl_num = n2_nucl_num; double* nucl_coord = &(n2_nucl_coord[0][0]); int64_t size_max; @@ -1165,42 +1317,23 @@ qmckl_exit_code rc; assert(!qmckl_electron_provided(context)); -rc = qmckl_set_electron_num (context, elec_up_num, elec_dn_num); +rc = qmckl_check(context, + qmckl_set_electron_num (context, elec_up_num, elec_dn_num) + ); assert(rc == QMCKL_SUCCESS); assert(qmckl_electron_provided(context)); -double k_ee = 0.; -double k_en = 0.; -rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); -assert(rc == QMCKL_SUCCESS); -assert(k_ee == 1.0); - -rc = qmckl_get_electron_rescale_factor_en (context, &k_en); -assert(rc == QMCKL_SUCCESS); -assert(k_en == 1.0); - -rc = qmckl_set_electron_rescale_factor_en(context, rescale_factor_kappa_en); -assert(rc == QMCKL_SUCCESS); - -rc = qmckl_set_electron_rescale_factor_ee(context, rescale_factor_kappa_ee); -assert(rc == QMCKL_SUCCESS); - -rc = qmckl_get_electron_rescale_factor_ee (context, &k_ee); -assert(rc == QMCKL_SUCCESS); -assert(k_ee == rescale_factor_kappa_ee); - -rc = qmckl_get_electron_rescale_factor_en (context, &k_en); -assert(rc == QMCKL_SUCCESS); -assert(k_en == rescale_factor_kappa_en); - - -rc = qmckl_set_electron_coord (context, 'N', walk_num, elec_coord, walk_num*3*elec_num); +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_get_electron_coord (context, 'N', 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] ); @@ -1211,34 +1344,27 @@ for (int64_t i=0 ; i<3*elec_num ; ++i) { assert(!qmckl_nucleus_provided(context)); -rc = qmckl_set_nucleus_num (context, nucl_num); +rc = qmckl_check(context, + qmckl_set_nucleus_num (context, nucl_num) + ); assert(rc == QMCKL_SUCCESS); assert(!qmckl_nucleus_provided(context)); -double k; -rc = qmckl_get_nucleus_rescale_factor (context, &k); -assert(rc == QMCKL_SUCCESS); -assert(k == 1.0); - - -rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa); -assert(rc == QMCKL_SUCCESS); - -rc = qmckl_get_nucleus_rescale_factor (context, &k); -assert(rc == QMCKL_SUCCESS); -assert(k == nucl_rescale_factor_kappa); - 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_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*nucl_num); +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_get_nucleus_coord (context, 'N', nucl_coord2, nucl_num*3); +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 ; idate > ctx->jastrow.asymp_jasb_date) { @@ -1373,7 +1500,7 @@ qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context) rc = qmckl_compute_asymp_jasb(context, ctx->jastrow.bord_num, ctx->jastrow.bord_vector, - rescale_factor_kappa_ee, + ctx->jastrow.rescale_factor_ee, ctx->jastrow.asymp_jasb); if (rc != QMCKL_SUCCESS) { return rc; @@ -1394,28 +1521,28 @@ qmckl_exit_code qmckl_provide_asymp_jasb(qmckl_context context) :END: #+NAME: qmckl_asymp_jasb_args - | Variable | Type | In/Out | Description | - |---------------------------+----------------------+--------+-------------------------| - | ~context~ | ~qmckl_context~ | in | Global state | - | ~bord_num~ | ~int64_t~ | in | Order of the polynomial | - | ~bord_vector~ | ~double[bord_num+1]~ | in | Values of b | - | ~rescale_factor_kappa_ee~ | ~double~ | in | Electron coordinates | - | ~asymp_jasb~ | ~double[2]~ | out | Asymptotic value | + | Variable | Type | In/Out | Description | + |---------------------+----------------------+--------+-------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~bord_num~ | ~int64_t~ | in | Order of the polynomial | + | ~bord_vector~ | ~double[bord_num+1]~ | in | Values of b | + | ~rescale_factor_ee~ | ~double~ | in | Electron coordinates | + | ~asymp_jasb~ | ~double[2]~ | out | Asymptotic value | #+begin_src f90 :comments org :tangle (eval f) :noweb yes -integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_kappa_ee, asymp_jasb) & +integer function qmckl_compute_asymp_jasb_f(context, bord_num, bord_vector, rescale_factor_ee, asymp_jasb) & result(info) use qmckl implicit none integer(qmckl_context), intent(in) :: context integer*8 , intent(in) :: bord_num double precision , intent(in) :: bord_vector(bord_num + 1) - double precision , intent(in) :: rescale_factor_kappa_ee + double precision , intent(in) :: rescale_factor_ee double precision , intent(out) :: asymp_jasb(2) integer*8 :: i, p double precision :: kappa_inv, x, asym_one - kappa_inv = 1.0d0 / rescale_factor_kappa_ee + kappa_inv = 1.0d0 / rescale_factor_ee info = QMCKL_SUCCESS @@ -1448,7 +1575,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( const qmckl_context context, const int64_t bord_num, const double* bord_vector, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, double* const asymp_jasb ) { if (context == QMCKL_NULL_CONTEXT){ @@ -1459,7 +1586,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( return QMCKL_INVALID_ARG_2; } - const double kappa_inv = 1.0 / rescale_factor_kappa_ee; + const double kappa_inv = 1.0 / rescale_factor_ee; const double asym_one = bord_vector[0] * kappa_inv / (1.0 + bord_vector[1] * kappa_inv); asymp_jasb[0] = asym_one; asymp_jasb[1] = 0.5 * asym_one; @@ -1483,7 +1610,7 @@ qmckl_exit_code qmckl_compute_asymp_jasb ( const qmckl_context context, const int64_t bord_num, const double* bord_vector, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, double* const asymp_jasb ); #+end_src @@ -1527,26 +1654,64 @@ double* bord_vector = &(n2_bord_vector[0]); double* cord_vector = &(n2_cord_vector[0][0]); int64_t dim_cord_vect=0; -/* Initialize the Jastrow data */ -rc = qmckl_init_jastrow(context); assert(!qmckl_jastrow_provided(context)); /* Set the data */ -rc = qmckl_set_jastrow_ord_num(context, aord_num, bord_num, cord_num); +rc = qmckl_check(context, + qmckl_set_jastrow_ord_num(context, aord_num, bord_num, cord_num) + ); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_jastrow_type_nucl_num(context, type_nucl_num); +rc = qmckl_check(context, + qmckl_set_jastrow_type_nucl_num(context, type_nucl_num) + ); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_jastrow_type_nucl_vector(context, type_nucl_vector, nucl_num); +rc = qmckl_check(context, + qmckl_set_jastrow_type_nucl_vector(context, type_nucl_vector, nucl_num) + ); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_jastrow_aord_vector(context, aord_vector,(aord_num+1)*type_nucl_num); +rc = qmckl_check(context, + qmckl_set_jastrow_aord_vector(context, aord_vector,(aord_num+1)*type_nucl_num) + ); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_jastrow_bord_vector(context, bord_vector,(bord_num+1)); +rc = qmckl_check(context, + qmckl_set_jastrow_bord_vector(context, bord_vector,(bord_num+1)) + ); assert(rc == QMCKL_SUCCESS); -rc = qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect); +rc = qmckl_check(context, + qmckl_get_jastrow_dim_cord_vect(context, &dim_cord_vect) + ); assert(rc == QMCKL_SUCCESS); -rc = qmckl_set_jastrow_cord_vector(context, cord_vector,dim_cord_vect*type_nucl_num); +rc = qmckl_check(context, + qmckl_set_jastrow_cord_vector(context, cord_vector,dim_cord_vect*type_nucl_num) + ); assert(rc == QMCKL_SUCCESS); +double k_ee = 0.; +double k_en[2] = { 0., 0. }; +rc = qmckl_check(context, + qmckl_set_jastrow_rescale_factor_en(context, rescale_factor_en, nucl_num) + ); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_check(context, + qmckl_set_jastrow_rescale_factor_ee(context, rescale_factor_ee) + ); +assert(rc == QMCKL_SUCCESS); + +rc = qmckl_check(context, + qmckl_get_jastrow_rescale_factor_ee (context, &k_ee) + ); +assert(rc == QMCKL_SUCCESS); +assert(k_ee == rescale_factor_ee); + +rc = qmckl_check(context, + qmckl_get_jastrow_rescale_factor_en (context, &(k_en[0]), nucl_num) + ); +assert(rc == QMCKL_SUCCESS); +for (int i=0 ; idate > ctx->jastrow.factor_ee_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_ee); - ctx->jastrow.factor_ee = NULL; + if (ctx->jastrow.factor_ee != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_ee); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_ee", + "Unable to free ctx->jastrow.factor_ee"); + } + ctx->jastrow.factor_ee = NULL; + } } /* Allocate array */ @@ -1660,7 +1832,7 @@ qmckl_exit_code qmckl_provide_factor_ee(qmckl_context context) ctx->electron.up_num, ctx->jastrow.bord_num, ctx->jastrow.bord_vector, - ctx->electron.ee_distance_rescaled, + ctx->jastrow.ee_distance_rescaled, ctx->jastrow.asymp_jasb, ctx->jastrow.factor_ee); if (rc != QMCKL_SUCCESS) { @@ -1892,9 +2064,12 @@ print("factor_ee :",factor_ee) assert(qmckl_jastrow_provided(context)); double factor_ee[walk_num]; -rc = qmckl_get_jastrow_factor_ee(context, factor_ee, walk_num); +rc = qmckl_check(context, + qmckl_get_jastrow_factor_ee(context, factor_ee, walk_num) + ); // calculate factor_ee +printf("%e\n%e\n\n",factor_ee[0],-4.282760865958113 ); assert(fabs(factor_ee[0]+4.282760865958113) < 1.e-12); #+end_src @@ -1979,8 +2154,15 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.factor_ee_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_ee_deriv_e); - ctx->jastrow.factor_ee_deriv_e = NULL; + if (ctx->jastrow.factor_ee_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_ee_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_ee_deriv_e", + "Unable to free ctx->jastrow.factor_ee_deriv_e"); + } + ctx->jastrow.factor_ee_deriv_e = NULL; + } } /* Allocate array */ @@ -2005,8 +2187,8 @@ qmckl_exit_code qmckl_provide_factor_ee_deriv_e(qmckl_context context) 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.ee_distance_rescaled, + ctx->jastrow.ee_distance_rescaled_deriv_e, ctx->jastrow.factor_ee_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; @@ -2562,8 +2744,15 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context) if (ctx->date > ctx->jastrow.factor_en_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_en); - ctx->jastrow.factor_en = NULL; + if (ctx->jastrow.factor_en != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_en); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_en", + "Unable to free ctx->jastrow.factor_en"); + } + ctx->jastrow.factor_en = NULL; + } } /* Allocate array */ @@ -2590,7 +2779,7 @@ qmckl_exit_code qmckl_provide_factor_en(qmckl_context context) ctx->jastrow.type_nucl_vector, ctx->jastrow.aord_num, ctx->jastrow.aord_vector, - ctx->electron.en_distance_rescaled, + ctx->jastrow.en_distance_rescaled, ctx->jastrow.factor_en); if (rc != QMCKL_SUCCESS) { return rc; @@ -2913,8 +3102,15 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.factor_en_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_en_deriv_e); - ctx->jastrow.factor_en_deriv_e = NULL; + if (ctx->jastrow.factor_en_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_en_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_en_deriv_e", + "Unable to free ctx->jastrow.factor_en_deriv_e"); + } + ctx->jastrow.factor_en_deriv_e = NULL; + } } /* Allocate array */ @@ -2941,8 +3137,8 @@ qmckl_exit_code qmckl_provide_factor_en_deriv_e(qmckl_context context) ctx->jastrow.type_nucl_vector, ctx->jastrow.aord_num, ctx->jastrow.aord_vector, - ctx->electron.en_distance_rescaled, - ctx->electron.en_distance_rescaled_deriv_e, + ctx->jastrow.en_distance_rescaled, + ctx->jastrow.en_distance_rescaled_deriv_e, ctx->jastrow.factor_en_deriv_e); if (rc != QMCKL_SUCCESS) { return rc; @@ -3258,7 +3454,506 @@ assert(fabs(factor_en_deriv_e[0][3][0]+0.9667363412285741 ) < 1.e-12); #+end_src -** Electron-electron rescaled distances for each order +** 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 +qmckl_exit_code qmckl_get_jastrow_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_jastrow_ee_distance_rescaled(qmckl_context context, double* const distance_rescaled) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_ee_distance_rescaled(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.num * ctx->electron.num * ctx->electron.walker.num; + memcpy(distance_rescaled, ctx->jastrow.ee_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_ee_distance_rescaled(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(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); + + + /* Compute if necessary */ + if (ctx->electron.walker.point.date > ctx->jastrow.ee_distance_rescaled_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow.ee_distance_rescaled != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.ee_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_ee_distance_rescaled", + "Unable to free ctx->jastrow.ee_distance_rescaled"); + } + ctx->jastrow.ee_distance_rescaled = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow.ee_distance_rescaled == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.num * ctx->electron.num * + ctx->electron.walker.num * sizeof(double); + double* ee_distance_rescaled = (double*) qmckl_malloc(context, mem_info); + + if (ee_distance_rescaled == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_ee_distance_rescaled", + NULL); + } + ctx->jastrow.ee_distance_rescaled = ee_distance_rescaled; + } + + qmckl_exit_code rc = + qmckl_compute_ee_distance_rescaled(context, + ctx->electron.num, + ctx->jastrow.rescale_factor_ee, + ctx->electron.walker.num, + ctx->electron.walker.point.coord.data, + ctx->jastrow.ee_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.ee_distance_rescaled_date = ctx->date; + } + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + :PROPERTIES: + :Name: qmckl_compute_ee_distance_rescaled + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_ee_distance_rescaled_args + | Variable | Type | In/Out | Description | + |---------------------+----------------------------------------+--------+--------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~elec_num~ | ~int64_t~ | in | Number of electrons | + | ~rescale_factor_ee~ | ~double~ | in | Factor to rescale ee distances | + | ~walk_num~ | ~int64_t~ | in | Number of walkers | + | ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | + | ~ee_distance~ | ~double[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_ee, walk_num, & + coord, ee_distance_rescaled) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: elec_num + double precision , intent(in) :: rescale_factor_ee + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: coord(elec_num,walk_num,3) + double precision , intent(out) :: ee_distance_rescaled(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(context, 'T', 'T', elec_num, elec_num, & + coord(1,k,1), elec_num * walk_num, & + coord(1,k,1), elec_num * walk_num, & + ee_distance_rescaled(1,1,k), elec_num, rescale_factor_ee) + if (info /= QMCKL_SUCCESS) then + exit + endif + end do + +end function qmckl_compute_ee_distance_rescaled_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_ee_distance_rescaled ( + const qmckl_context context, + const int64_t elec_num, + const double rescale_factor_ee, + const int64_t walk_num, + const double* coord, + double* const ee_distance_rescaled ); + #+end_src + + #+CALL: generate_c_interface(table=qmckl_ee_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_ee_distance_rescaled & + (context, elec_num, rescale_factor_ee, walk_num, coord, ee_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 + real (c_double ) , intent(in) , value :: rescale_factor_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(elec_num,elec_num,walk_num) + + integer(c_int32_t), external :: qmckl_compute_ee_distance_rescaled_f + info = qmckl_compute_ee_distance_rescaled_f & + (context, elec_num, rescale_factor_ee, walk_num, coord, ee_distance_rescaled) + + end function qmckl_compute_ee_distance_rescaled + #+end_src + +*** Test + + #+begin_src python :results output :exports none +import numpy as np + +kappa = 1.0 + +elec_1_w1 = np.array( [-0.250655104764153, 0.503070975550133 , -0.166554344502303]) +elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274 , 0.187773606533075]) +elec_5_w1 = np.array( [-0.127732483187947, -0.138975497694196 , -8.669850480215846E-002]) +elec_6_w1 = np.array( [-0.232271834949124, -1.059321673434182E-002 , -0.504862241464867]) + +print ( "[0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_1_w1)) )/kappa ) +print ( "[0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-elec_2_w1)) )/kappa ) +print ( "[1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-elec_1_w1)) )/kappa ) +print ( "[5][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-elec_5_w1)) )/kappa ) +print ( "[5][6] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-elec_6_w1)) )/kappa ) +print ( "[6][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-elec_5_w1)) )/kappa ) + #+end_src + + #+RESULTS: + : [0][0] : 0.0 + : [0][1] : 0.5502278003524018 + : [1][0] : 0.5502278003524018 + : [5][5] : 0.0 + : [5][6] : 0.3622098222364193 + : [6][5] : 0.3622098222364193 + + #+begin_src c :tangle (eval c_test) +assert(qmckl_electron_provided(context)); + + +double ee_distance_rescaled[walk_num * elec_num * elec_num]; +rc = qmckl_get_jastrow_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.5502278003524018) < 1.e-12); + +// (0,0,1) == 0. +assert(ee_distance_rescaled[5*elec_num + 5] == 0.); + +// (1,0,1) == (0,1,1) +assert(ee_distance_rescaled[5*elec_num+6] == ee_distance_rescaled[6*elec_num+5]); + +// value of (1,0,1) +assert(fabs(ee_distance_rescaled[5*elec_num+6]-0.3622098222364193) < 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_jastrow_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_jastrow_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*) context; + assert (ctx != NULL); + + size_t sze = 4 * ctx->electron.num * ctx->electron.num * ctx->electron.walker.num; + memcpy(distance_rescaled_deriv_e, ctx->jastrow.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*) context; + assert (ctx != NULL); + + + /* Compute if necessary */ + if (ctx->electron.walker.point.date > ctx->jastrow.ee_distance_rescaled_deriv_e_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow.ee_distance_rescaled_deriv_e != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.ee_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_ee_distance_rescaled_deriv_e", + "Unable to free ctx->jastrow.ee_distance_rescaled_deriv_e"); + } + ctx->jastrow.ee_distance_rescaled_deriv_e = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow.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.walker.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_provide_ee_distance_rescaled_deriv_e", + NULL); + } + ctx->jastrow.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->jastrow.rescale_factor_ee, + ctx->electron.walker.num, + ctx->electron.walker.point.coord.data, + ctx->jastrow.ee_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.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 + | Variable | Type | In/Out | Description | + |-----------------------+-------------------------------------------+--------+-------------------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~elec_num~ | ~int64_t~ | in | Number of electrons | + | ~rescale_factor_ee~ | ~double~ | in | Factor to rescale ee distances | + | ~walk_num~ | ~int64_t~ | in | Number of walkers | + | ~coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | + | ~ee_distance_deriv_e~ | ~double[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_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_ee + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: coord(elec_num,walk_num,3) + 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,k,1), elec_num*walk_num, & + coord(1,k,1), elec_num*walk_num, & + ee_distance_rescaled_deriv_e(1,1,1,k), elec_num, rescale_factor_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_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_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_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_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_jastrow_ee_distance_rescaled_deriv_e(context, ee_distance_rescaled_deriv_e); + +// TODO: Get exact values +//// (e1,e2,w) +//// (0,0,0) == 0. +//assert(ee_distance[0] == 0.); +// +//// (1,0,0) == (0,1,0) +//assert(ee_distance[1] == ee_distance[elec_num]); +// +//// value of (1,0,0) +//assert(fabs(ee_distance[1]-7.152322512964209) < 1.e-12); +// +//// (0,0,1) == 0. +//assert(ee_distance[elec_num*elec_num] == 0.); +// +//// (1,0,1) == (0,1,1) +//assert(ee_distance[elec_num*elec_num+1] == ee_distance[elec_num*elec_num+elec_num]); +// +//// value of (1,0,1) +//assert(fabs(ee_distance[elec_num*elec_num+1]-6.5517646321055665) < 1.e-12); + + #+end_src + +** Electron-electron-nucleus rescaled distances for each order ~een_rescaled_e~ stores the table of the rescaled distances between all pairs of electrons and raised to the power \(p\) defined by ~cord_num~: @@ -3334,8 +4029,15 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) if (ctx->date > ctx->jastrow.een_rescaled_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.een_rescaled_e); - ctx->jastrow.een_rescaled_e = NULL; + if (ctx->jastrow.een_rescaled_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.een_rescaled_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_een_rescaled_e", + "Unable to free ctx->jastrow.een_rescaled_e"); + } + ctx->jastrow.een_rescaled_e = NULL; + } } /* Allocate array */ @@ -3349,7 +4051,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) if (een_rescaled_e == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_een_rescaled_e", + "qmckl_provide_een_rescaled_e", NULL); } ctx->jastrow.een_rescaled_e = een_rescaled_e; @@ -3359,7 +4061,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) ctx->electron.walker.num, ctx->electron.num, ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_ee, + ctx->jastrow.rescale_factor_ee, ctx->electron.ee_distance, ctx->jastrow.een_rescaled_e); if (rc != QMCKL_SUCCESS) { @@ -3381,19 +4083,19 @@ qmckl_exit_code qmckl_provide_een_rescaled_e(qmckl_context context) :END: #+NAME: qmckl_factor_een_rescaled_e_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 | - | ~cord_num~ | ~int64_t~ | in | Order of polynomials | - | ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances | - | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances | - | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances | + | 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 | + | ~cord_num~ | ~int64_t~ | in | Order of polynomials | + | ~rescale_factor_ee~ | ~double~ | in | Factor to rescale ee distances | + | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances | + | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_e_doc_f( & - context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & + context, walk_num, elec_num, cord_num, rescale_factor_ee, & ee_distance, een_rescaled_e) & result(info) use qmckl @@ -3402,7 +4104,7 @@ integer function qmckl_compute_een_rescaled_e_doc_f( & integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: cord_num - double precision , intent(in) :: rescale_factor_kappa_ee + double precision , intent(in) :: rescale_factor_ee double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) double precision , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) double precision,dimension(:,:),allocatable :: een_rescaled_e_ij @@ -3444,7 +4146,7 @@ integer function qmckl_compute_een_rescaled_e_doc_f( & do j = 1, elec_num do i = 1, j - 1 k = k + 1 - een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_kappa_ee * ee_distance(i, j, nw)) + een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_ee * ee_distance(i, j, nw)) end do end do @@ -3489,7 +4191,7 @@ end function qmckl_compute_een_rescaled_e_doc_f const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ); #+end_src @@ -3499,7 +4201,7 @@ end function qmckl_compute_een_rescaled_e_doc_f #+RESULTS: #+begin_src f90 :tangle (eval f) :comments org :exports none integer(c_int32_t) function qmckl_compute_een_rescaled_e_doc & - (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & + (context, walk_num, elec_num, cord_num, rescale_factor_ee, & ee_distance, een_rescaled_e) & bind(C) result(info) @@ -3510,13 +4212,13 @@ end function qmckl_compute_een_rescaled_e_doc_f 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 :: cord_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee + real (c_double ) , intent(in) , value :: rescale_factor_ee real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) real (c_double ) , intent(out) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) integer(c_int32_t), external :: qmckl_compute_een_rescaled_e_doc_f info = qmckl_compute_een_rescaled_e_doc_f & - (context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e) + (context, walk_num, elec_num, cord_num, rescale_factor_ee, ee_distance, een_rescaled_e) end function qmckl_compute_een_rescaled_e_doc #+end_src @@ -3527,7 +4229,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ) { @@ -3589,8 +4291,8 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( k = 0; for (int i = 0; i < elec_num; ++i) { for (int j = 0; j < i; ++j) { - // een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_kappa_ee * ee_distance(i, j, nw)); - een_rescaled_e_ij[k + elec_pairs] = exp(-rescale_factor_kappa_ee * \ + // een_rescaled_e_ij(k, 2) = dexp(-rescale_factor_ee * ee_distance(i, j, nw)); + een_rescaled_e_ij[k + elec_pairs] = exp(-rescale_factor_ee * \ ee_distance[j + i*elec_num + nw*(elec_num*elec_num)]); k = k + 1; } @@ -3648,7 +4350,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ); #+end_src @@ -3659,7 +4361,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ); #+end_src @@ -3670,7 +4372,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ); #+end_src @@ -3681,14 +4383,14 @@ qmckl_exit_code qmckl_compute_een_rescaled_e_hpc ( const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* ee_distance, double* const een_rescaled_e ) { #ifdef HAVE_HPC - return qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e); + return qmckl_compute_een_rescaled_e_hpc(context, walk_num, elec_num, cord_num, rescale_factor_ee, ee_distance, een_rescaled_e); #else - return qmckl_compute_een_rescaled_e_doc(context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, ee_distance, een_rescaled_e); + return qmckl_compute_een_rescaled_e_doc(context, walk_num, elec_num, cord_num, rescale_factor_ee, ee_distance, een_rescaled_e); #endif } #+end_src @@ -3770,7 +4472,7 @@ assert(fabs(een_rescaled_e[0][2][1][4]-0.0005700154999202759) < 1.e-12); assert(fabs(een_rescaled_e[0][2][1][5]-0.3424402276009091) < 1.e-12); #+end_src -** Electron-electron rescaled distances for each order and derivatives +** Electron-electron-nucleus rescaled distances for each order and derivatives ~een_rescaled_e_deriv_e~ stores the table of the derivatives of the rescaled distances between all pairs of electrons and raised to the @@ -3845,8 +4547,15 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.een_rescaled_e_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.een_rescaled_e_deriv_e); - ctx->jastrow.een_rescaled_e_deriv_e = NULL; + if (ctx->jastrow.een_rescaled_e_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.een_rescaled_e_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_een_rescaled_e_deriv_e", + "Unable to free ctx->jastrow.een_rescaled_e_deriv_e"); + } + ctx->jastrow.een_rescaled_e_deriv_e = NULL; + } } /* Allocate array */ @@ -3860,7 +4569,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) if (een_rescaled_e_deriv_e == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_een_rescaled_e_deriv_e", + "qmckl_provide_een_rescaled_e_deriv_e", NULL); } ctx->jastrow.een_rescaled_e_deriv_e = een_rescaled_e_deriv_e; @@ -3870,7 +4579,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) ctx->electron.walker.num, ctx->electron.num, ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_ee, + ctx->jastrow.rescale_factor_ee, ctx->electron.walker.point.coord.data, ctx->electron.ee_distance, ctx->jastrow.een_rescaled_e, @@ -3894,21 +4603,21 @@ qmckl_exit_code qmckl_provide_een_rescaled_e_deriv_e(qmckl_context context) :END: #+NAME: qmckl_factor_een_rescaled_e_deriv_e_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 | - | ~cord_num~ | ~int64_t~ | in | Order of polynomials | - | ~rescale_factor_kappa_ee~ | ~double~ | in | Factor to rescale ee distances | - | ~coord_ee~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | - | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances | - | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-electron distances | - | ~een_rescaled_e_deriv_e~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | out | Electron-electron rescaled distances | + | 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 | + | ~cord_num~ | ~int64_t~ | in | Order of polynomials | + | ~rescale_factor_ee~ | ~double~ | in | Factor to rescale ee distances | + | ~coord_ee~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~ee_distance~ | ~double[walk_num][elec_num][elec_num]~ | in | Electron-electron distances | + | ~een_rescaled_e~ | ~double[walk_num][0:cord_num][elec_num][elec_num]~ | in | Electron-electron distances | + | ~een_rescaled_e_deriv_e~ | ~double[walk_num][0:cord_num][elec_num][4][elec_num]~ | out | Electron-electron rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( & - context, walk_num, elec_num, cord_num, rescale_factor_kappa_ee, & + context, walk_num, elec_num, cord_num, rescale_factor_ee, & coord_ee, ee_distance, een_rescaled_e, een_rescaled_e_deriv_e) & result(info) use qmckl @@ -3917,7 +4626,7 @@ integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( & integer*8 , intent(in) :: walk_num integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: cord_num - double precision , intent(in) :: rescale_factor_kappa_ee + double precision , intent(in) :: rescale_factor_ee double precision , intent(in) :: coord_ee(elec_num,3,walk_num) double precision , intent(in) :: ee_distance(elec_num,elec_num,walk_num) double precision , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) @@ -3966,7 +4675,7 @@ integer function qmckl_compute_factor_een_rescaled_e_deriv_e_f( & ! prepare the actual een table do l = 1, cord_num - kappa_l = - dble(l) * rescale_factor_kappa_ee + kappa_l = - dble(l) * rescale_factor_ee do j = 1, elec_num do i = 1, elec_num een_rescaled_e_deriv_e(i, 1, j, l, nw) = kappa_l * elec_dist_deriv_e(1, i, j) @@ -4003,7 +4712,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f const int64_t walk_num, const int64_t elec_num, const int64_t cord_num, - const double rescale_factor_kappa_ee, + const double rescale_factor_ee, const double* coord_ee, const double* ee_distance, const double* een_rescaled_e, @@ -4020,7 +4729,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f walk_num, & elec_num, & cord_num, & - rescale_factor_kappa_ee, & + rescale_factor_ee, & coord_ee, & ee_distance, & een_rescaled_e, & @@ -4034,7 +4743,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f 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 :: cord_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_ee + real (c_double ) , intent(in) , value :: rescale_factor_ee real (c_double ) , intent(in) :: coord_ee(elec_num,3,walk_num) real (c_double ) , intent(in) :: ee_distance(elec_num,elec_num,walk_num) real (c_double ) , intent(in) :: een_rescaled_e(elec_num,elec_num,0:cord_num,walk_num) @@ -4046,7 +4755,7 @@ end function qmckl_compute_factor_een_rescaled_e_deriv_e_f walk_num, & elec_num, & cord_num, & - rescale_factor_kappa_ee, & + rescale_factor_ee, & coord_ee, & ee_distance, & een_rescaled_e, & @@ -4137,7 +4846,6 @@ for l in range(0,cord_num+1): : een_rescaled_e_deriv_e[2, 1, 6, 2] = 0.5880599146214673 #+begin_src c :tangle (eval c_test) -//assert(qmckl_electron_provided(context)); double een_rescaled_e_deriv_e[walk_num][(cord_num + 1)][elec_num][4][elec_num]; size_max=walk_num*(cord_num + 1)*elec_num*4*elec_num; rc = qmckl_get_jastrow_een_rescaled_e_deriv_e(context, @@ -4152,7 +4860,596 @@ assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][4] + 0.00041342909359870457) < 1. assert(fabs(een_rescaled_e_deriv_e[0][2][1][0][5] + 0.5880599146214673 ) < 1.e-12); #+end_src -** Electron-nucleus rescaled distances for each order +** Electron-nucleus rescaled distances + + ~en_distance_rescaled~ stores the matrix of the rescaled distances between + electrons and nuclei. + + \[ + 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 +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*) context; + assert (ctx != NULL); + + size_t sze = ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num; + memcpy(distance_rescaled, ctx->jastrow.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*) context; + assert (ctx != NULL); + + if (!(ctx->nucleus.provided)) { + return QMCKL_NOT_PROVIDED; + } + + /* Compute if necessary */ + if (ctx->electron.walker.point.date > ctx->jastrow.en_distance_rescaled_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow.en_distance_rescaled != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.en_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_en_distance_rescaled", + "Unable to free ctx->jastrow.en_distance_rescaled"); + } + ctx->jastrow.en_distance_rescaled = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow.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.walker.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_provide_en_distance_rescaled", + NULL); + } + ctx->jastrow.en_distance_rescaled = en_distance_rescaled; + } + + qmckl_exit_code rc = + qmckl_compute_en_distance_rescaled(context, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow.rescale_factor_en, + ctx->electron.walker.num, + ctx->electron.walker.point.coord.data, + ctx->nucleus.coord.data, + ctx->jastrow.en_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.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 + | Variable | Type | In/Out | Description | + |------------------------+----------------------------------------+--------+-----------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~elec_num~ | ~int64_t~ | in | Number of electrons | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~rescale_factor_en~ | ~double[nucl_num]~ | in | The factor for rescaled distances | + | ~walk_num~ | ~int64_t~ | in | Number of walkers | + | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | + | ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | + | ~en_distance_rescaled~ | ~double[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, rescale_factor_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) :: rescale_factor_en(nucl_num) + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: elec_coord(elec_num,walk_num,3) + double precision , intent(in) :: nucl_coord(nucl_num,3) + double precision , intent(out) :: en_distance_rescaled(elec_num,nucl_num,walk_num) + + integer*8 :: i, k + double precision :: coord(3) + + 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 + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + do i=1, nucl_num + coord(1:3) = nucl_coord(i,1:3) + do k=1,walk_num + info = qmckl_distance_rescaled(context, 'T', 'T', elec_num, 1_8, & + elec_coord(1,k,1), elec_num*walk_num, coord, 1_8, & + en_distance_rescaled(1,i,k), elec_num, rescale_factor_en(i)) + if (info /= QMCKL_SUCCESS) then + return + endif + end do + 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* rescale_factor_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, & + rescale_factor_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) :: rescale_factor_en(nucl_num) + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) + 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, & + rescale_factor_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 + +kappa = 1.0 + +elec_1_w1 = np.array( [-0.250655104764153, 0.503070975550133 , -0.166554344502303]) +elec_2_w1 = np.array( [-0.587812193472177, -0.128751981129274 , 0.187773606533075]) +elec_5_w1 = np.array( [-0.127732483187947, -0.138975497694196 , -8.669850480215846E-002]) +elec_6_w1 = np.array( [-0.232271834949124, -1.059321673434182E-002 , -0.504862241464867]) +nucl_1 = np.array( [ 0., 0., 0. ]) +nucl_2 = np.array( [ 0., 0., 2.059801 ]) + +print ( "[0][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_1)) )/kappa ) +print ( "[1][0] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_1_w1-nucl_2)) )/kappa ) +print ( "[0][1] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_2_w1-nucl_1)) )/kappa ) +print ( "[0][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-nucl_1)) )/kappa ) +print ( "[1][5] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_5_w1-nucl_2)) )/kappa ) +print ( "[0][6] : ", (1.0 - np.exp(-kappa * np.linalg.norm(elec_6_w1-nucl_1)) )/kappa ) + + #+end_src + + #+RESULTS: + : [0][0] : 0.4435709484118112 + : [1][0] : 0.8993601506374442 + : [0][1] : 0.46760219699910477 + : [0][5] : 0.1875631834682101 + : [1][5] : 0.8840716589810682 + : [0][6] : 0.42640469987268914 + + #+begin_src c :tangle (eval c_test) + +assert(qmckl_electron_provided(context)); +assert(qmckl_nucleus_provided(context)); + +double en_distance_rescaled[walk_num][nucl_num][elec_num]; + +rc = qmckl_check(context, + 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.4435709484118112) < 1.e-12); + +// (1,2,1) +assert(fabs(en_distance_rescaled[0][1][0] - 0.8993601506374442) < 1.e-12); + +// (2,1,1) +assert(fabs(en_distance_rescaled[0][0][1] - 0.46760219699910477) < 1.e-12); + +// (1,1,2) +assert(fabs(en_distance_rescaled[0][0][5] - 0.1875631834682101) < 1.e-12); + +// (1,2,2) +assert(fabs(en_distance_rescaled[0][1][5] - 0.8840716589810682) < 1.e-12); + +// (2,1,2) +assert(fabs(en_distance_rescaled[0][0][6] - 0.42640469987268914) < 1.e-12); + + #+end_src + +** Electron-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][nucl_num][elec_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*) context; + assert (ctx != NULL); + + size_t sze = 4 * ctx->electron.num * ctx->nucleus.num * ctx->electron.walker.num; + memcpy(distance_rescaled_deriv_e, ctx->jastrow.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*) context; + assert (ctx != NULL); + + if (!(ctx->nucleus.provided)) { + return QMCKL_NOT_PROVIDED; + } + + /* Compute if necessary */ + if (ctx->electron.walker.point.date > ctx->jastrow.en_distance_rescaled_deriv_e_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow.en_distance_rescaled_deriv_e != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->jastrow.en_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_en_distance_rescaled_deriv_e", + "Unable to free ctx->jastrow.en_distance_rescaled_deriv_e"); + } + ctx->jastrow.en_distance_rescaled_deriv_e = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow.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.walker.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_provide_en_distance_rescaled_deriv_e", + NULL); + } + ctx->jastrow.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->jastrow.rescale_factor_en, + ctx->electron.walker.num, + ctx->electron.walker.point.coord.data, + ctx->nucleus.coord.data, + ctx->jastrow.en_distance_rescaled_deriv_e); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow.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 + | Variable | Type | In/Out | Description | + |--------------------------------+-------------------------------------------+--------+---------------------------------------| + | ~context~ | ~qmckl_context~ | in | Global state | + | ~elec_num~ | ~int64_t~ | in | Number of electrons | + | ~nucl_num~ | ~int64_t~ | in | Number of nuclei | + | ~rescale_factor_en~ | ~double[nucl_num]~ | in | The factors for rescaled distances | + | ~walk_num~ | ~int64_t~ | in | Number of walkers | + | ~elec_coord~ | ~double[3][walk_num][elec_num]~ | in | Electron coordinates | + | ~nucl_coord~ | ~double[3][elec_num]~ | in | Nuclear coordinates | + | ~en_distance_rescaled_deriv_e~ | ~double[walk_num][nucl_num][elec_num][4]~ | 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_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_en(nucl_num) + integer*8 , intent(in) :: walk_num + double precision , intent(in) :: elec_coord(elec_num,walk_num,3) + double precision , intent(in) :: nucl_coord(nucl_num,3) + double precision , intent(out) :: en_distance_rescaled_deriv_e(4,elec_num,nucl_num,walk_num) + + integer*8 :: i, k + double precision :: coord(3) + + 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 + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_5 + return + endif + + do i=1, nucl_num + coord(1:3) = nucl_coord(i,1:3) + do k=1,walk_num + info = qmckl_distance_rescaled_deriv_e(context, 'T', 'T', elec_num, 1_8, & + elec_coord(1,k,1), elec_num*walk_num, coord, 1_8, & + en_distance_rescaled_deriv_e(1,1,i,k), elec_num, rescale_factor_en(i)) + if (info /= QMCKL_SUCCESS) then + return + endif + end do + 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_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_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) :: rescale_factor_en(nucl_num) + integer (c_int64_t) , intent(in) , value :: walk_num + real (c_double ) , intent(in) :: elec_coord(elec_num,walk_num,3) + real (c_double ) , intent(in) :: nucl_coord(elec_num,3) + real (c_double ) , intent(out) :: en_distance_rescaled_deriv_e(4,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_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)); + +assert(qmckl_nucleus_provided(context)); + +double en_distance_rescaled_deriv_e[walk_num][4][nucl_num][elec_num]; + +rc = qmckl_check(context, + 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 +//// (e,n,w) in Fortran notation +//// (1,1,1) +//assert(fabs(en_distance_rescaled[0][0][0] - 7.546738741619978) < 1.e-12); +// +//// (1,2,1) +//assert(fabs(en_distance_rescaled[0][1][0] - 8.77102435246984) < 1.e-12); +// +//// (2,1,1) +//assert(fabs(en_distance_rescaled[0][0][1] - 3.698922010513608) < 1.e-12); +// +//// (1,1,2) +//assert(fabs(en_distance_rescaled[1][0][0] - 5.824059436060509) < 1.e-12); +// +//// (1,2,2) +//assert(fabs(en_distance_rescaled[1][1][0] - 7.080482110317645) < 1.e-12); +// +//// (2,1,2) +//assert(fabs(en_distance_rescaled[1][0][1] - 3.1804527583077356) < 1.e-12); + + #+end_src + +** Electron-electron-nucleus rescaled distances for each order ~een_rescaled_n~ stores the table of the rescaled distances between electrons and nucleii raised to the power \(p\) defined by ~cord_num~: @@ -4228,8 +5525,15 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) if (ctx->date > ctx->jastrow.een_rescaled_n_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.een_rescaled_n); - ctx->jastrow.een_rescaled_n = NULL; + if (ctx->jastrow.een_rescaled_n != NULL) { + rc = qmckl_free(context, ctx->jastrow.een_rescaled_n); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_een_rescaled_n", + "Unable to free ctx->jastrow.een_rescaled_n"); + } + ctx->jastrow.een_rescaled_n = NULL; + } } /* Allocate array */ @@ -4243,7 +5547,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) if (een_rescaled_n == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_een_rescaled_n", + "qmckl_provide_een_rescaled_n", NULL); } ctx->jastrow.een_rescaled_n = een_rescaled_n; @@ -4254,7 +5558,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) ctx->electron.num, ctx->nucleus.num, ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_en, + ctx->jastrow.rescale_factor_en, ctx->electron.en_distance, ctx->jastrow.een_rescaled_n); if (rc != QMCKL_SUCCESS) { @@ -4276,20 +5580,20 @@ qmckl_exit_code qmckl_provide_een_rescaled_n(qmckl_context context) :END: #+NAME: qmckl_factor_een_rescaled_n_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 | - | ~cord_num~ | ~int64_t~ | in | Order of polynomials | - | ~rescale_factor_kappa_en~ | ~double~ | 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]~ | out | Electron-nucleus rescaled distances | + | 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 | + | ~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]~ | out | Electron-nucleus rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_een_rescaled_n_f( & - context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_kappa_en, & + context, walk_num, elec_num, nucl_num, cord_num, rescale_factor_en, & en_distance, een_rescaled_n) & result(info) use qmckl @@ -4299,7 +5603,7 @@ integer function qmckl_compute_een_rescaled_n_f( & integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: cord_num - double precision , intent(in) :: rescale_factor_kappa_en + double precision , intent(in) :: rescale_factor_en(nucl_num) double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) double precision , intent(out) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) double precision :: x @@ -4341,7 +5645,7 @@ integer function qmckl_compute_een_rescaled_n_f( & do a = 1, nucl_num do i = 1, elec_num - een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw)) + een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_en(a) * en_distance(i, a, nw)) end do end do @@ -4364,7 +5668,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_n ( const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, - const double rescale_factor_kappa_en, + const double* rescale_factor_en, const double* en_distance, double* const een_rescaled_n ) { @@ -4400,8 +5704,8 @@ qmckl_exit_code qmckl_compute_een_rescaled_n ( // prepare the actual een table //een_rescaled_n(:, :, 0, nw) = 1.0d0 een_rescaled_n[i + a * elec_num + 0 + nw * elec_num*nucl_num*(cord_num+1)] = 1.0; - //een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_kappa_en * en_distance(i, a, nw)) - een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = exp(-rescale_factor_kappa_en * \ + //een_rescaled_n(i, a, 1, nw) = dexp(-rescale_factor_en(a) * en_distance(i, a, nw)) + een_rescaled_n[i + a*elec_num + elec_num*nucl_num + nw*elec_num*nucl_num*(cord_num+1)] = exp(-rescale_factor_en[a] * \ en_distance[i + a*elec_num + nw*elec_num*nucl_num]); } } @@ -4430,7 +5734,7 @@ qmckl_exit_code qmckl_compute_een_rescaled_n ( const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, - const double rescale_factor_kappa_en, + const double* rescale_factor_en, const double* en_distance, double* const een_rescaled_n ); #+end_src @@ -4570,8 +5874,15 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.een_rescaled_n_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.een_rescaled_n_deriv_e); - ctx->jastrow.een_rescaled_n_deriv_e = NULL; + if (ctx->jastrow.een_rescaled_n_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.een_rescaled_n_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_een_rescaled_n_deriv_e", + "Unable to free ctx->jastrow.een_rescaled_n_deriv_e"); + } + ctx->jastrow.een_rescaled_n_deriv_e = NULL; + } } /* Allocate array */ @@ -4585,7 +5896,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) if (een_rescaled_n_deriv_e == NULL) { return qmckl_failwith( context, QMCKL_ALLOCATION_FAILED, - "qmckl_een_rescaled_n_deriv_e", + "qmckl_provide_een_rescaled_n_deriv_e", NULL); } ctx->jastrow.een_rescaled_n_deriv_e = een_rescaled_n_deriv_e; @@ -4596,7 +5907,7 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) ctx->electron.num, ctx->nucleus.num, ctx->jastrow.cord_num, - ctx->electron.rescale_factor_kappa_en, + ctx->jastrow.rescale_factor_en, ctx->electron.walker.point.coord.data, ctx->nucleus.coord.data, ctx->electron.en_distance, @@ -4621,24 +5932,24 @@ qmckl_exit_code qmckl_provide_een_rescaled_n_deriv_e(qmckl_context context) :END: #+NAME: qmckl_compute_factor_een_rescaled_n_deriv_e_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 | - | ~cord_num~ | ~int64_t~ | in | Order of polynomials | - | ~rescale_factor_kappa_en~ | ~double~ | in | Factor to rescale ee distances | - | ~coord_ee~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | - | ~coord_en~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | - | ~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 distances | - | ~een_rescaled_n_deriv_e~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | out | Electron-nucleus rescaled distances | + | 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 | + | ~cord_num~ | ~int64_t~ | in | Order of polynomials | + | ~rescale_factor_en~ | ~double[nucl_num]~ | in | Factor to rescale ee distances | + | ~coord_ee~ | ~double[walk_num][3][elec_num]~ | in | Electron coordinates | + | ~coord_en~ | ~double[3][nucl_num]~ | in | Nuclear coordinates | + | ~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 distances | + | ~een_rescaled_n_deriv_e~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | out | Electron-nucleus rescaled distances | #+begin_src f90 :comments org :tangle (eval f) :noweb yes integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( & context, walk_num, elec_num, nucl_num, & - cord_num, rescale_factor_kappa_en, & + cord_num, rescale_factor_en, & coord_ee, coord_en, en_distance, een_rescaled_n, een_rescaled_n_deriv_e) & result(info) use qmckl @@ -4648,7 +5959,7 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( & integer*8 , intent(in) :: elec_num integer*8 , intent(in) :: nucl_num integer*8 , intent(in) :: cord_num - double precision , intent(in) :: rescale_factor_kappa_en + double precision , intent(in) :: rescale_factor_en(nucl_num) double precision , intent(in) :: coord_ee(elec_num,3,walk_num) double precision , intent(in) :: coord_en(nucl_num,3) double precision , intent(in) :: en_distance(elec_num,nucl_num,walk_num) @@ -4703,8 +6014,8 @@ integer function qmckl_compute_factor_een_rescaled_n_deriv_e_f( & end do do l = 0, cord_num - kappa_l = - dble(l) * rescale_factor_kappa_en do a = 1, nucl_num + kappa_l = - dble(l) * rescale_factor_en(a) do i = 1, elec_num een_rescaled_n_deriv_e(i, 1, a, l, nw) = kappa_l * elnuc_dist_deriv_e(1, i, a) een_rescaled_n_deriv_e(i, 2, a, l, nw) = kappa_l * elnuc_dist_deriv_e(2, i, a) @@ -4741,7 +6052,7 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f const int64_t elec_num, const int64_t nucl_num, const int64_t cord_num, - const double rescale_factor_kappa_en, + const double* rescale_factor_en, const double* coord_ee, const double* coord_en, const double* en_distance, @@ -4759,7 +6070,7 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f elec_num, & nucl_num, & cord_num, & - rescale_factor_kappa_en, & + rescale_factor_en, & coord_ee, & coord_en, & en_distance, & @@ -4775,7 +6086,7 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f 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 :: cord_num - real (c_double ) , intent(in) , value :: rescale_factor_kappa_en + real (c_double ) , intent(in) :: rescale_factor_en(nucl_num) real (c_double ) , intent(in) :: coord_ee(elec_num,3,walk_num) real (c_double ) , intent(in) :: coord_en(nucl_num,3) real (c_double ) , intent(in) :: en_distance(nucl_num,elec_num,walk_num) @@ -4789,7 +6100,7 @@ end function qmckl_compute_factor_een_rescaled_n_deriv_e_f elec_num, & nucl_num, & cord_num, & - rescale_factor_kappa_en, & + rescale_factor_en, & coord_ee, & coord_en, & en_distance, & @@ -5172,8 +6483,15 @@ qmckl_exit_code qmckl_provide_tmp_c(qmckl_context context) if (ctx->date > ctx->jastrow.tmp_c_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.tmp_c); - ctx->jastrow.tmp_c = NULL; + if (ctx->jastrow.tmp_c != NULL) { + rc = qmckl_free(context, ctx->jastrow.tmp_c); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_tmp_c", + "Unable to free ctx->jastrow.tmp_c"); + } + ctx->jastrow.tmp_c = NULL; + } } /* Allocate array */ @@ -5267,8 +6585,15 @@ qmckl_exit_code qmckl_provide_dtmp_c(qmckl_context context) if (ctx->date > ctx->jastrow.dtmp_c_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.dtmp_c); - ctx->jastrow.dtmp_c = NULL; + if (ctx->jastrow.dtmp_c != NULL) { + rc = qmckl_free(context, ctx->jastrow.dtmp_c); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_dtmp_c", + "Unable to free ctx->jastrow.dtmp_c"); + } + ctx->jastrow.dtmp_c = NULL; + } } /* Allocate array */ @@ -7037,8 +8362,15 @@ qmckl_exit_code qmckl_provide_factor_een(qmckl_context context) if (ctx->date > ctx->jastrow.factor_een_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_een); - ctx->jastrow.factor_een = NULL; + if (ctx->jastrow.factor_een != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_een); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_een", + "Unable to free ctx->jastrow.factor_een"); + } + ctx->jastrow.factor_een = NULL; + } } /* Allocate array */ @@ -7552,8 +8884,15 @@ qmckl_exit_code qmckl_provide_factor_een_deriv_e(qmckl_context context) if (ctx->date > ctx->jastrow.factor_een_deriv_e_date) { if (ctx->electron.walker.num > ctx->electron.walker_old.num) { - free(ctx->jastrow.factor_een_deriv_e); - ctx->jastrow.factor_een_deriv_e = NULL; + if (ctx->jastrow.factor_een_deriv_e != NULL) { + rc = qmckl_free(context, ctx->jastrow.factor_een_deriv_e); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_factor_een_deriv_e", + "Unable to free ctx->jastrow.factor_een_deriv_e"); + } + ctx->jastrow.factor_een_deriv_e = NULL; + } } /* Allocate array */ diff --git a/org/qmckl_mo.org b/org/qmckl_mo.org index 972cef0..dd9eca4 100644 --- a/org/qmckl_mo.org +++ b/org/qmckl_mo.org @@ -3,16 +3,16 @@ #+INCLUDE: ../tools/lib.org The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO -coefficient matrix \[C\]. Using these coefficients (e.g. from Hartree Fock SCF method) +coefficient matrix $C$. Using these coefficients (e.g. from Hartree Fock method) the MOs are defined as follows: \[ \phi_i(\mathbf{r}) = C_i * \chi_i (\mathbf{r}) \] - -In this section we demonstrate how to use the QMCkl specific DGEMM -function to calculate the MOs. +By default, all the MOs are computed. A subset of MOs can be selected +for computation, for example to remove computation of the useless +virtual MOs present in a MO coefficient matrix. * Headers :noexport: @@ -84,11 +84,11 @@ int main() { The following arrays are stored in the context: - |-----------------+--------------------+----------------------------------------| - | ~mo_num~ | | Number of MOs | - | ~coefficient~ | ~[mo_num][ao_num]~ | Orbital coefficients | - | ~coefficient_t~ | ~[ao_num][mo_num]~ | Transposed of the Orbital coefficients | - |-----------------+--------------------+----------------------------------------| + |-------------------+--------------------+----------------------------------------| + | ~mo_num~ | | Number of MOs | + | ~coefficient~ | ~[mo_num][ao_num]~ | MO coefficients | + | ~coefficient_t~ | ~[ao_num][mo_num]~ | Transposed of the Orbital coefficients | + |-------------------+--------------------+----------------------------------------| Computed data: @@ -103,7 +103,7 @@ int main() { #+begin_src c :comments org :tangle (eval h_private_type) typedef struct qmckl_mo_basis_struct { - int64_t mo_num; + int64_t mo_num; double * restrict coefficient; double * restrict coefficient_t; @@ -144,150 +144,6 @@ qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) { } #+end_src -** Access functions - - #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code -qmckl_get_mo_basis_mo_num (const qmckl_context context, - int64_t* mo_num); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_get_mo_basis_mo_num (const qmckl_context context, - int64_t* mo_num) -{ - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return qmckl_failwith( context, - QMCKL_INVALID_CONTEXT, - "qmckl_get_mo_basis_mo_num", - NULL); - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - int32_t mask = 1; - - if ( (ctx->mo_basis.uninitialized & mask) != 0) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_get_mo_basis_mo_num", - NULL); - } - - assert (ctx->mo_basis.mo_num > (int64_t) 0); - ,*mo_num = ctx->mo_basis.mo_num; - return QMCKL_SUCCESS; -} - - #+end_src - - #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code -qmckl_get_mo_basis_coefficient (const qmckl_context context, - double* const coefficient, - const int64_t size_max); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :exports none -qmckl_exit_code -qmckl_get_mo_basis_coefficient (const qmckl_context context, - double* const coefficient, - const int64_t size_max) -{ - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return qmckl_failwith( context, - QMCKL_INVALID_CONTEXT, - "qmckl_get_mo_basis_coefficient", - NULL); - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - int32_t mask = 1 << 1; - - if ( (ctx->ao_basis.uninitialized & mask) != 0) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_get_mo_basis_coefficient", - NULL); - } - - if (coefficient == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_get_mo_basis_coefficient", - "NULL pointer"); - } - - if (size_max < ctx->ao_basis.ao_num * ctx->mo_basis.mo_num) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_3, - "qmckl_get_mo_basis_coefficient", - "Array too small. Expected mo_num * ao_num"); - } - - assert (ctx->mo_basis.coefficient != NULL); - memcpy(coefficient, ctx->mo_basis.coefficient, - ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double)); - - return QMCKL_SUCCESS; -} - #+end_src - - When all the data for the AOs have been provided, the following - function returns ~true~. - - #+begin_src c :comments org :tangle (eval h_func) -bool qmckl_mo_basis_provided (const qmckl_context context); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :exports none -bool qmckl_mo_basis_provided(const qmckl_context context) { - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return false; - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - return ctx->mo_basis.provided; -} - - - #+end_src - -*** Fortran interfaces - - #+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_get_mo_basis_mo_num (context, & - mo_num) bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - integer (c_int64_t) , intent(out) :: mo_num - end function qmckl_get_mo_basis_mo_num -end interface - -interface - integer(c_int32_t) function qmckl_get_mo_basis_coefficient(context, & - coefficient, size_max) bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - double precision, intent(out) :: coefficient(*) - integer (c_int64_t) , intent(in), value :: size_max - end function qmckl_get_mo_basis_coefficient -end interface - - #+end_src - ** Initialization functions To set the basis set, all the following functions need to be @@ -398,6 +254,15 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { qmckl_context_struct* const ctx = (qmckl_context_struct*) context; assert (ctx != NULL); + if (ctx->mo_basis.coefficient_t != NULL) { + qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient_t); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_finalize_mo_basis", + NULL); + } + } + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; mem_info.size = ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double); double* new_array = (double*) qmckl_malloc(context, mem_info); @@ -410,15 +275,6 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { assert (ctx->mo_basis.coefficient != NULL); - if (ctx->mo_basis.coefficient_t != NULL) { - qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, rc, - "qmckl_finalize_mo_basis", - NULL); - } - } - for (int64_t i=0 ; iao_basis.ao_num ; ++i) { for (int64_t j=0 ; jmo_basis.mo_num ; ++j) { new_array[i*ctx->mo_basis.mo_num + j] = ctx->mo_basis.coefficient[j*ctx->ao_basis.ao_num + i]; @@ -426,11 +282,269 @@ qmckl_exit_code qmckl_finalize_mo_basis(qmckl_context context) { } ctx->mo_basis.coefficient_t = new_array; - qmckl_exit_code rc = QMCKL_SUCCESS; - return rc; + + qmckl_exit_code rc; + if (ctx->mo_basis.mo_vgl != NULL) { + rc = qmckl_free(context, ctx->mo_basis.mo_vgl); + if (rc != QMCKL_SUCCESS) return rc; + ctx->mo_basis.mo_vgl = NULL; + ctx->mo_basis.mo_vgl_date = 0; + } + + if (ctx->mo_basis.mo_value != NULL) { + rc = qmckl_free(context, ctx->mo_basis.mo_value); + if (rc != QMCKL_SUCCESS) return rc; + ctx->mo_basis.mo_value = NULL; + ctx->mo_basis.mo_value_date = 0; + } + + return QMCKL_SUCCESS; } #+end_src +** Access functions + + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code +qmckl_get_mo_basis_mo_num (const qmckl_context context, + int64_t* mo_num); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_mo_basis_mo_num (const qmckl_context context, + int64_t* mo_num) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_mo_basis_mo_num", + NULL); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + int32_t mask = 1; + + if ( (ctx->mo_basis.uninitialized & mask) != 0) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_get_mo_basis_mo_num", + NULL); + } + + assert (ctx->mo_basis.mo_num > (int64_t) 0); + ,*mo_num = ctx->mo_basis.mo_num; + return QMCKL_SUCCESS; +} + + #+end_src + + #+begin_src c :comments org :tangle (eval h_func) :exports none +qmckl_exit_code +qmckl_get_mo_basis_coefficient (const qmckl_context context, + double* const coefficient, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :exports none +qmckl_exit_code +qmckl_get_mo_basis_coefficient (const qmckl_context context, + double* const coefficient, + const int64_t size_max) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_mo_basis_coefficient", + NULL); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + int32_t mask = 1 << 1; + + if ( (ctx->mo_basis.uninitialized & mask) != 0) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_get_mo_basis_coefficient", + NULL); + } + + if (coefficient == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_mo_basis_coefficient", + "NULL pointer"); + } + + if (size_max < ctx->ao_basis.ao_num * ctx->mo_basis.mo_num) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_mo_basis_coefficient", + "Array too small: expected mo_num * ao_num."); + } + + assert (ctx->mo_basis.coefficient != NULL); + memcpy(coefficient, ctx->mo_basis.coefficient, + ctx->ao_basis.ao_num * ctx->mo_basis.mo_num * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + + When all the data for the AOs have been provided, the following + function returns ~true~. + + #+begin_src c :comments org :tangle (eval h_func) +bool qmckl_mo_basis_provided (const qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :exports none +bool qmckl_mo_basis_provided(const qmckl_context context) { + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return false; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + return ctx->mo_basis.provided; +} + + + #+end_src + +*** Fortran interfaces + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_get_mo_basis_mo_num (context, & + mo_num) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(out) :: mo_num + end function qmckl_get_mo_basis_mo_num +end interface + +interface + integer(c_int32_t) function qmckl_get_mo_basis_coefficient(context, & + coefficient, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in) , value :: context + double precision, intent(out) :: coefficient(*) + integer (c_int64_t) , intent(in), value :: size_max + end function qmckl_get_mo_basis_coefficient +end interface + + #+end_src + +** Update + + Useless MOs can be removed, for instance virtual MOs in a single + determinant calculation. + + To select a subset of MOs that will be kept, create an array of + integers of size =mo_num=. If the integer is zero, the MO is dropped, + otherwise it is kept. + + #+begin_src c :comments org :tangle (eval h_func) +bool qmckl_mo_basis_select_mo (const qmckl_context context, + const int32_t* keep, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :exports none +bool qmckl_mo_basis_select_mo (const qmckl_context context, + const int32_t* keep, + const int64_t size_max) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_mo_basis_select_mo", + NULL); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + if ( !(qmckl_mo_basis_provided(context)) ) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_get_mo_basis_select_mo", + NULL); + } + + if (keep == NULL) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_2, + "qmckl_get_mo_basis_select_mo", + "NULL pointer"); + } + + const int64_t mo_num = ctx->mo_basis.mo_num; + const int64_t ao_num = ctx->ao_basis.ao_num; + + if (size_max < mo_num) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_mo_basis_select_mo", + "Array too small: expected mo_num."); + } + + int64_t mo_num_new = 0; + for (int64_t i=0 ; imo_basis.coefficient[i*ao_num]), ao_num*sizeof(double)); + ++k; + } + } + + qmckl_exit_code rc = qmckl_free(context, ctx->mo_basis.coefficient); + if (rc != QMCKL_SUCCESS) return rc; + + ctx->mo_basis.coefficient = coefficient; + ctx->mo_basis.mo_num = mo_num_new; + + rc = qmckl_finalize_mo_basis(context); + return rc; +} + + #+end_src + +*** Fortran interface + + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none +interface + integer(c_int32_t) function qmckl_mo_basis_select_mo (context, & + keep, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (c_int64_t) , intent(in), value :: context + integer (c_int32_t) , intent(in) :: keep(*) + integer (c_int64_t) , intent(in), value :: size_max + end function qmckl_mo_basis_select_mo +end interface + #+end_src + * Computation ** Computation of MOs: values only @@ -628,15 +742,8 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context) #+end_src + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none - rc = qmckl_provide_ao_basis_ao_value(context); - if (rc != QMCKL_SUCCESS) { - return qmckl_failwith( context, - QMCKL_NOT_PROVIDED, - "qmckl_ao_value", - NULL); - } - if (ctx->mo_basis.mo_vgl_date == ctx->point.date) { // mo_vgl has been computed at this step: Just copy the data. @@ -653,6 +760,14 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context) } else { + rc = qmckl_provide_ao_basis_ao_value(context); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_ao_value", + NULL); + } + rc = qmckl_compute_mo_basis_mo_value(context, ctx->ao_basis.ao_num, ctx->mo_basis.mo_num, @@ -664,6 +779,8 @@ qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context) } #+end_src + + #+CALL: write_provider_post( group="mo_basis", data="mo_value" ) #+RESULTS: @@ -727,8 +844,8 @@ integer function qmckl_compute_mo_basis_mo_value_doc_f(context, & do j=1,point_num mo_value(:,j) = 0.d0 do k=1,ao_num - if (ao_value(k,j) /= 0.d0) then - c1 = ao_value(k,j) + c1 = ao_value(k,j) + if (c1 /= 0.d0) then do i=1,mo_num mo_value(i,j) = mo_value(i,j) + coefficient_t(i,k) * c1 end do @@ -736,7 +853,7 @@ integer function qmckl_compute_mo_basis_mo_value_doc_f(context, & end do end do - else ! dgemm + else ! dgemm for checking LDA = size(coefficient_t,1) LDB = size(ao_value,1) @@ -874,7 +991,7 @@ qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context, } } - int64_t n; + int64_t n=0; for (n=0 ; n < nidx-4 ; n+=4) { const double* restrict ck1 = coefficient_t + idx[n ]*mo_num; @@ -895,8 +1012,7 @@ qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context, } } - const int64_t n0 = n < 0 ? 0 : n; - for (int64_t m=n0 ; m < nidx ; m+=1) { + for (int64_t m=n ; m < nidx ; m+=1) { const double* restrict ck = coefficient_t + idx[m]*mo_num; const double a1 = av1[m]; @@ -1355,7 +1471,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context, } } - int64_t n; + int64_t n=0; for (n=0 ; n < nidx-4 ; n+=4) { const double* restrict ck1 = coefficient_t + idx[n ]*mo_num; @@ -1400,8 +1516,7 @@ qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context, } } - const int64_t n0 = n < 0 ? 0 : n; - for (int64_t m=n0 ; m < nidx ; m+=1) { + for (int64_t m=n ; m < nidx ; m+=1) { const double* restrict ck = coefficient_t + idx[m]*mo_num; const double a1 = av1[m]; const double a2 = av2[m]; @@ -1603,7 +1718,7 @@ rc = qmckl_get_ao_basis_ao_vgl(context, &(ao_vgl[0][0][0]), assert (rc == QMCKL_SUCCESS); /* Set up MO data */ -const int64_t mo_num = chbrclf_mo_num; +int64_t mo_num = chbrclf_mo_num; rc = qmckl_set_mo_basis_mo_num(context, mo_num); assert (rc == QMCKL_SUCCESS); @@ -1699,6 +1814,32 @@ printf(" mo_vgl mo_vgl[1][26][223] %25.15e\n", mo_vgl[2][1][3]); printf(" mo_vgl mo_vgl[0][26][224] %25.15e\n", mo_vgl[2][0][3]); printf(" mo_vgl mo_vgl[1][26][224] %25.15e\n", mo_vgl[2][1][3]); printf("\n"); + + +/* Check selection of MOs */ + + int32_t keep[mo_num]; + for (int i=0 ; inucleus.uninitialized = (1 << 3) - 1; /* Default values */ - ctx->nucleus.rescale_factor_kappa = 1.0; return QMCKL_SUCCESS; } #+end_src ** Access functions - #+NAME:post - #+begin_src c :exports none -if ( (ctx->nucleus.uninitialized & mask) != 0) { - return NULL; -} - #+end_src - #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_exit_code @@ -272,53 +258,6 @@ end interface #+begin_src c :comments org :tangle (eval h_func) :exports none qmckl_exit_code -qmckl_get_nucleus_rescale_factor(const qmckl_context context, - double* const rescale_factor_kappa); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_get_nucleus_rescale_factor (const qmckl_context context, - double* const rescale_factor_kappa) -{ - - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return QMCKL_INVALID_CONTEXT; - } - - if (rescale_factor_kappa == NULL) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_get_nucleus_rescale_factor", - "rescale_factor_kappa is a null pointer"); - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - assert (ctx->nucleus.rescale_factor_kappa > 0.0); - - (*rescale_factor_kappa) = ctx->nucleus.rescale_factor_kappa; - - return QMCKL_SUCCESS; -} - #+end_src - - #+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_get_nucleus_rescale_factor(context, kappa) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - real (c_double) , intent(out) :: kappa - end function qmckl_get_nucleus_rescale_factor -end interface - #+end_src - - #+begin_src c :comments org :tangle (eval h_func) :exports none -qmckl_exit_code qmckl_get_nucleus_coord(const qmckl_context context, const char transp, double* const coord, @@ -438,7 +377,7 @@ if (mask != 0 && !(ctx->nucleus.uninitialized & mask)) { } #+end_src - #+NAME:post2 + #+NAME:post #+begin_src c :exports none ctx->nucleus.uninitialized &= ~mask; ctx->nucleus.provided = (ctx->nucleus.uninitialized == 0); @@ -475,7 +414,7 @@ qmckl_set_nucleus_num(qmckl_context context, ctx->nucleus.num = num; - <> + <> } #+end_src @@ -541,7 +480,7 @@ qmckl_set_nucleus_charge(qmckl_context context, "Error in vector->double* conversion"); } - <> + <> } #+end_src @@ -619,7 +558,7 @@ qmckl_set_nucleus_coord(qmckl_context context, } if (rc != QMCKL_SUCCESS) return rc; - <> + <> } #+end_src @@ -638,55 +577,12 @@ interface end interface #+end_src - #+begin_src c :comments org :tangle (eval h_func) -qmckl_exit_code -qmckl_set_nucleus_rescale_factor(qmckl_context context, - const double kappa); - #+end_src - - Sets the rescale parameter for the nuclear distances. - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code -qmckl_set_nucleus_rescale_factor(qmckl_context context, - const double rescale_factor_kappa) -{ - int32_t mask = 0; // Can be updated - - <> - - if (rescale_factor_kappa <= 0.0) { - return qmckl_failwith( context, - QMCKL_INVALID_ARG_2, - "qmckl_set_nucleus_rescale_factor", - "rescale_factor_kappa cannot be <= 0."); - } - - ctx->nucleus.rescale_factor_kappa = rescale_factor_kappa; - - return QMCKL_SUCCESS; -} - #+end_src - - #+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_set_nucleus_rescale_factor(context, kappa) & - bind(C) - use, intrinsic :: iso_c_binding - import - implicit none - integer (c_int64_t) , intent(in) , value :: context - real (c_double) , intent(in) , value :: kappa - end function qmckl_set_nucleus_rescale_factor -end interface - #+end_src ** Test #+begin_src c :tangle (eval c_test) const double* nucl_charge = chbrclf_charge; const double* nucl_coord = &(chbrclf_nucl_coord[0][0]); -const double nucl_rescale_factor_kappa = 2.0; /* --- */ @@ -700,38 +596,25 @@ assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_nucleus_num (context, chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_nucleus_provided(context)); rc = qmckl_get_nucleus_num (context, &n); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(n == chbrclf_nucl_num); -double k; -rc = qmckl_get_nucleus_rescale_factor (context, &k); -assert(rc == QMCKL_SUCCESS); -assert(k == 1.0); - - -rc = qmckl_set_nucleus_rescale_factor (context, nucl_rescale_factor_kappa); -assert(rc == QMCKL_SUCCESS); - -rc = qmckl_get_nucleus_rescale_factor (context, &k); -assert(rc == QMCKL_SUCCESS); -assert(k == nucl_rescale_factor_kappa); - double nucl_coord2[3*chbrclf_nucl_num]; rc = qmckl_get_nucleus_coord (context, 'T', nucl_coord2, 3*chbrclf_nucl_num); assert(rc == QMCKL_NOT_PROVIDED); rc = qmckl_set_nucleus_coord (context, 'T', &(nucl_coord[0]), 3*chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert(!qmckl_nucleus_provided(context)); rc = qmckl_get_nucleus_coord (context, 'N', nucl_coord2, 3*chbrclf_nucl_num); -assert(rc == QMCKL_SUCCESS); +qmckl_check(context, rc); for (size_t k=0 ; k<3 ; ++k) { for (int64_t i=0 ; inucleus.num * ctx->nucleus.num; - if (sze > size_max) { - return qmckl_failwith(context, - QMCKL_INVALID_ARG_3, - "qmckl_get_nucleus_nn_distance_rescaled", - "Array too small"); - } - rc = qmckl_double_of_matrix(context, - ctx->nucleus.nn_distance_rescaled, - distance_rescaled, - size_max); - return rc; -} - #+end_src - - #+begin_src f90 :tangle (eval fh_func) :comments org :exports none -interface - integer(c_int32_t) function qmckl_get_nucleus_nn_distance_rescaled(context, distance_rescaled, 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) :: distance_rescaled(*) - integer (c_int64_t) , intent(in) , value :: size_max - end function -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_nn_distance_rescaled(qmckl_context context); - #+end_src - - #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none -qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context) -{ - /* Check input parameters */ - if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { - return (char) 0; - } - - qmckl_context_struct* const ctx = (qmckl_context_struct*) context; - assert (ctx != NULL); - - if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED; - - /* Allocate array */ - if (ctx->nucleus.nn_distance_rescaled.data == NULL) { - ctx->nucleus.nn_distance_rescaled = - qmckl_matrix_alloc(context, ctx->nucleus.num, ctx->nucleus.num); - - if (ctx->nucleus.nn_distance_rescaled.data == NULL) { - return qmckl_failwith( context, - QMCKL_ALLOCATION_FAILED, - "qmckl_nn_distance_rescaled", - NULL); - } - } - - qmckl_exit_code rc = - qmckl_compute_nn_distance_rescaled(context, - ctx->nucleus.num, - ctx->nucleus.rescale_factor_kappa, - ctx->nucleus.coord.data, - ctx->nucleus.nn_distance_rescaled.data); - if (rc != QMCKL_SUCCESS) { - return rc; - } - - ctx->nucleus.nn_distance_rescaled_date = ctx->date; - - return QMCKL_SUCCESS; -} - #+end_src - -*** Compute - - #+NAME: qmckl_nn_distance_rescaled_args - | qmckl_context | context | in | Global state | - | int64_t | nucl_num | in | Number of nuclei | - | double | coord[3][nucl_num] | in | Nuclear coordinates (au) | - | 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, rescale_factor_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) :: rescale_factor_kappa - double precision , intent(in) :: coord(nucl_num,3) - double precision , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num) - - integer*8 :: k - - info = QMCKL_SUCCESS - - if (context == QMCKL_NULL_CONTEXT) then - info = QMCKL_INVALID_CONTEXT - return - endif - - if (nucl_num <= 0) then - info = QMCKL_INVALID_ARG_2 - return - endif - - info = qmckl_distance_rescaled(context, 'T', 'T', nucl_num, nucl_num, & - coord, nucl_num, & - coord, nucl_num, & - nn_distance_rescaled, nucl_num, rescale_factor_kappa) - -end function qmckl_compute_nn_distance_rescaled_f - #+end_src - - #+begin_src c :tangle (eval h_private_func) :comments org :exports none -qmckl_exit_code qmckl_compute_nn_distance_rescaled ( - const qmckl_context context, - const int64_t nucl_num, - const double rescale_factor_kappa, - const double* coord, - double* const nn_distance_rescaled ); - #+end_src - - - #+CALL: generate_c_interface(table=qmckl_nn_distance_rescaled_args,rettyp="qmckl_exit_code",fname="qmckl_compute_nn_distance") - - #+RESULTS: - #+begin_src f90 :tangle (eval f) :comments org :exports none - integer(c_int32_t) function qmckl_compute_nn_distance_rescaled & - (context, nucl_num, rescale_factor_kappa, coord, nn_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 :: nucl_num - real (c_double ) , intent(in) , value :: rescale_factor_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, rescale_factor_kappa, coord, nn_distance_rescaled) - - end function qmckl_compute_nn_distance_rescaled - #+end_src - -*** Test - - #+begin_src c :tangle (eval c_test) -/* Reference input data */ -/* TODO */ - -//assert(qmckl_nucleus_provided(context)); -// -//double distance[nucl_num*nucl_num]; -//rc = qmckl_get_nucleus_nn_distance(context, distance, nucl_num*nucl_num); -//assert(distance[0] == 0.); -//assert(distance[1] == distance[nucl_num]); -//assert(fabs(distance[1]-2.070304721365169) < 1.e-12); - - #+end_src - ** Nuclear repulsion energy \[ diff --git a/org/qmckl_trexio.org b/org/qmckl_trexio.org index a94bd85..4490957 100644 --- a/org/qmckl_trexio.org +++ b/org/qmckl_trexio.org @@ -1097,37 +1097,25 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name, const int6 rc = qmckl_trexio_read_electron_X(context, file); if (rc != QMCKL_SUCCESS) { trexio_close(file); - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_trexio_read", - "Error reading electron"); + return rc; } rc = qmckl_trexio_read_nucleus_X(context, file); if (rc != QMCKL_SUCCESS) { trexio_close(file); - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_trexio_read", - "Error reading nucleus"); + return rc; } rc = qmckl_trexio_read_ao_X(context, file); if (rc != QMCKL_SUCCESS) { trexio_close(file); - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_trexio_read", - "Error reading AOs"); + return rc; } rc = qmckl_trexio_read_mo_X(context, file); if (rc != QMCKL_SUCCESS) { trexio_close(file); - return qmckl_failwith( context, - QMCKL_FAILURE, - "qmckl_trexio_read", - "Error reading MOs"); + return rc; } trexio_close(file); @@ -1149,27 +1137,19 @@ qmckl_trexio_read(const qmckl_context context, const char* file_name, const int6 #ifdef HAVE_TREXIO qmckl_exit_code rc; -char fname[256]; -char message[256]; +char filename[256]; #ifndef QMCKL_TEST_DIR #error "QMCKL_TEST_DIR is not defined" #endif -strncpy(fname, QMCKL_TEST_DIR,255); -strncat(fname, "/chbrclf", 255); -printf("Test file: %s\n", fname); +strncpy(filename, QMCKL_TEST_DIR,255); +strncat(filename, "/chbrclf", 255); +printf("Test file: %s\n", filename); -rc = qmckl_trexio_read(context, fname, 255); +rc = qmckl_trexio_read(context, filename, 255); -if (rc != QMCKL_SUCCESS) { - printf("%s\n", qmckl_string_of_error(rc)); - qmckl_get_error(context, &rc, fname, message); - printf("%s\n", fname); - printf("%s\n", message); - } - -assert ( rc == QMCKL_SUCCESS ); +qmckl_check(context, rc); #+end_src @@ -1179,11 +1159,11 @@ assert ( rc == QMCKL_SUCCESS ); printf("Electrons\n"); int64_t up_num, dn_num; rc = qmckl_get_electron_up_num(context, &up_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert (up_num == chbrclf_elec_up_num); rc = qmckl_get_electron_down_num(context, &dn_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert (dn_num == chbrclf_elec_dn_num); #+end_src @@ -1195,13 +1175,13 @@ printf("Nuclei\n"); int64_t nucl_num; rc = qmckl_get_nucleus_num(context, &nucl_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); assert (nucl_num == chbrclf_nucl_num); printf("Nuclear charges\n"); double * charge = (double*) malloc (nucl_num * sizeof(double)); rc = qmckl_get_nucleus_charge(context, charge, nucl_num); -assert (rc == QMCKL_SUCCESS); +qmckl_check(context, rc); for (int i=0 ; i=42", "wheel", - "numpy>=1.17.3" + "oldest-supported-numpy" ] build-backend = "setuptools.build_meta" diff --git a/python/setup.py b/python/setup.py index d2abecf..1a58a9e 100644 --- a/python/setup.py +++ b/python/setup.py @@ -71,5 +71,7 @@ setup(name = MODULE_NAME, "Operating System :: MacOS" ], python_requires = ">=3.0", - install_requires = ["numpy>=1.17.3"] + # The ABI incompatibility of numpy is a pain, for now set the + # min numpy version such that we have wheels for CPython 3.5 & 3.6 + install_requires = ["numpy>=1.13.3"] ) diff --git a/python/src/qmckl.i b/python/src/qmckl.i index e8b2d89..5573bdb 100644 --- a/python/src/qmckl.i +++ b/python/src/qmckl.i @@ -20,6 +20,7 @@ %include typemaps.i %apply int *OUTPUT { qmckl_exit_code *exit_code}; +%apply double *OUTPUT { double* det_l}; /* Avoid passing file_name length as an additiona argument */ %apply (char *STRING, int LENGTH) { (const char* file_name, const int64_t size_max) }; @@ -29,6 +30,9 @@ %cstring_bounded_output(char* function_name, 1024); %cstring_bounded_output(char* message, 1024); +%cstring_bounded_output(char* const basis_type, 2); +/* Required for qmckl_last_error function to work */ +%cstring_bounded_output(char* buffer, 1024); /* This block is needed make SWIG convert NumPy arrays to/from from the C pointer and size_max argument. NOTE: `numpy.i` interface file is not part of SWIG but it is included in the numpy distribution (under numpy/tools/swig/numpy.i) @@ -48,6 +52,12 @@ import_array(); /* Include typemaps generated by the process_header.py script */ %include "qmckl_include.i" +/* Some custom array typemaps which are not generated by process_header.py */ +%apply ( double* IN_ARRAY1 , int64_t DIM1 ) { ( const double * A, const int64_t size_max_A) }; +%apply ( double* IN_ARRAY1 , int64_t DIM1 ) { ( const double * B, const int64_t size_max_B) }; +%apply ( double* ARGOUT_ARRAY1 , int64_t DIM1 ) { ( double* const C, const int64_t size_max_C) }; +%apply ( double* ARGOUT_ARRAY1 , int64_t DIM1 ) { ( double* const B, const int64_t size_max_B) }; + /* Handle properly get_point */ diff --git a/python/test/test_api.py b/python/test/test_api.py index ccffd45..58b8476 100644 --- a/python/test/test_api.py +++ b/python/test/test_api.py @@ -33,6 +33,9 @@ assert mo_num == 404 pq.set_electron_coord(ctx, 'T', walk_num, coord) +ao_type = pq.get_ao_basis_type(ctx) +assert 'G' in ao_type + size_max = 5*walk_num*elec_num*mo_num mo_vgl = pq.get_mo_basis_mo_vgl(ctx, size_max)