diff --git a/org/qmckl_jastrow_champ.org b/org/qmckl_jastrow_champ.org index cf4b37e..94bb814 100644 --- a/org/qmckl_jastrow_champ.org +++ b/org/qmckl_jastrow_champ.org @@ -3964,7 +3964,15 @@ qmckl_compute_jastrow_champ_factor_ee_gl_hpc(const qmckl_context context, double * restrict factor_ee_gl_3 = factor_ee_gl_2 + elec_num; for (int i = 0; i < elec_num; ++i) { - if (j == i) continue; + if (j == i) { + if (j == 0) { + factor_ee_gl_0[i] = 0.0; + factor_ee_gl_1[i] = 0.0; + factor_ee_gl_2[i] = 0.0; + factor_ee_gl_3[i] = 0.0; + } + continue; + } double x = xj[i]; @@ -10600,6 +10608,35 @@ qmckl_get_jastrow_champ_factor_een_gl(qmckl_context context, memcpy(factor_een_gl, ctx->jastrow_champ.factor_een_gl, sze*sizeof(double)); + return QMCKL_SUCCESS; +} + +qmckl_exit_code +qmckl_get_jastrow_champ_factor_een_grad(qmckl_context context, + double* const factor_een_grad, + const int64_t size_max) +{ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_exit_code rc; + + rc = qmckl_provide_jastrow_champ_factor_een_grad(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + int64_t sze = ctx->electron.walker.num * 3 * ctx->electron.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_champ_factor_een_grad", + "Array too small. Expected 3*walk_num*elec_num"); + } + memcpy(factor_een_grad, ctx->jastrow_champ.factor_een_grad, sze*sizeof(double)); + return QMCKL_SUCCESS; } #+end_src @@ -10618,11 +10655,24 @@ interface real(c_double), intent(out) :: factor_een_gl(size_max) end function qmckl_get_jastrow_champ_factor_een_gl end interface + +interface + integer(qmckl_exit_code) function qmckl_get_jastrow_champ_factor_een_grad (context, & + factor_een_grad, size_max) bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + integer (qmckl_context) , intent(in), value :: context + integer(c_int64_t), intent(in), value :: size_max + real(c_double), intent(out) :: factor_een_grad(size_max) + end function qmckl_get_jastrow_champ_factor_een_grad +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_jastrow_champ_factor_een_gl(qmckl_context context); +qmckl_exit_code qmckl_provide_jastrow_champ_factor_een_grad(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none @@ -10726,6 +10776,110 @@ qmckl_exit_code qmckl_provide_jastrow_champ_factor_een_gl(qmckl_context context) ctx->jastrow_champ.factor_een_gl_date = ctx->date; } + return QMCKL_SUCCESS; +} + + +qmckl_exit_code qmckl_provide_jastrow_champ_factor_een_grad(qmckl_context context) +{ + + qmckl_exit_code rc; + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return QMCKL_NULL_CONTEXT; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + if (ctx->jastrow_champ.cord_num > 0) { + + /* Check if en rescaled distance is provided */ + rc = qmckl_provide_een_rescaled_e(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_een_rescaled_n(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance is provided */ + rc = qmckl_provide_een_rescaled_e_gl(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_een_rescaled_n_gl(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_jastrow_champ_c_vector_full(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if en rescaled distance derivatives is provided */ + rc = qmckl_provide_lkpm_combined_index(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if tmp_c is provided */ + rc = qmckl_provide_tmp_c(context); + if(rc != QMCKL_SUCCESS) return rc; + + /* Check if dtmp_c is provided */ + rc = qmckl_provide_dtmp_c(context); + if(rc != QMCKL_SUCCESS) return rc; + + } + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow_champ.factor_een_grad_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow_champ.factor_een_grad != NULL) { + rc = qmckl_free(context, ctx->jastrow_champ.factor_een_grad); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_jastrow_champ_factor_een_grad", + "Unable to free ctx->jastrow_champ.factor_een_grad"); + } + ctx->jastrow_champ.factor_een_grad = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow_champ.factor_een_grad == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = 4 * ctx->electron.num * ctx->electron.walker.num * sizeof(double); + double* factor_een_grad = (double*) qmckl_malloc(context, mem_info); + + if (factor_een_grad == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_jastrow_champ_factor_een_grad", + NULL); + } + ctx->jastrow_champ.factor_een_grad = factor_een_grad; + } + + rc = qmckl_compute_jastrow_champ_factor_een_grad(context, + ctx->electron.walker.num, + ctx->electron.num, + ctx->nucleus.num, + ctx->jastrow_champ.cord_num, + ctx->jastrow_champ.dim_c_vector, + ctx->jastrow_champ.c_vector_full, + ctx->jastrow_champ.lkpm_combined_index, + ctx->jastrow_champ.tmp_c, + ctx->jastrow_champ.dtmp_c, + ctx->jastrow_champ.een_rescaled_n, + ctx->jastrow_champ.een_rescaled_n_gl, + ctx->jastrow_champ.factor_een_grad); + + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->jastrow_champ.factor_een_grad_date = ctx->date; + } + return QMCKL_SUCCESS; } #+end_src @@ -11280,25 +11434,27 @@ print("factor_een:",factor_een) /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); -double factor_een_gl[walk_num][4][elec_num]; -rc = qmckl_check(context, - rc = qmckl_get_jastrow_champ_factor_een_gl(context, &(factor_een_gl[0][0][0]),4*walk_num*elec_num) - ); -assert(rc == QMCKL_SUCCESS); - /* -printf("%20.15e\n", factor_een_gl[0][0][0]); -assert(fabs(8.967809309100624e-02 - factor_een_gl[0][0][0]) < 1e-12); +{ + double factor_een_gl[walk_num][4][elec_num]; + rc = qmckl_check(context, + rc = qmckl_get_jastrow_champ_factor_een_gl(context, &(factor_een_gl[0][0][0]),4*walk_num*elec_num) + ); + assert(rc == QMCKL_SUCCESS); -printf("%20.15e\n", factor_een_gl[0][1][1]); -assert(fabs(3.543090132452453e-02 - factor_een_gl[0][1][1]) < 1e-12); + printf("%20.15e\n", factor_een_gl[0][0][0]); + assert(fabs(8.967809309100624e-02 - factor_een_gl[0][0][0]) < 1e-12); -printf("%20.15e\n", factor_een_gl[0][2][2]); -assert(fabs(8.996044894431991e-04 - factor_een_gl[0][2][2]) < 1e-12); + printf("%20.15e\n", factor_een_gl[0][1][1]); + assert(fabs(3.543090132452453e-02 - factor_een_gl[0][1][1]) < 1e-12); -printf("%20.15e\n", factor_een_gl[0][3][3]); -assert(fabs(-1.175028308456619e+00 - factor_een_gl[0][3][3]) < 1e-12); -*/ + printf("%20.15e\n", factor_een_gl[0][2][2]); + assert(fabs(8.996044894431991e-04 - factor_een_gl[0][2][2]) < 1e-12); + + printf("%20.15e\n", factor_een_gl[0][3][3]); + assert(fabs(-1.175028308456619e+00 - factor_een_gl[0][3][3]) < 1e-12); +} +,*/ { printf("factor_een_gl_hpc\n"); @@ -11392,6 +11548,389 @@ assert(fabs(-1.175028308456619e+00 - factor_een_gl[0][3][3]) < 1e-12); } #+end_src +**** Compute Gradient only + :PROPERTIES: + :Name: qmckl_compute_jastrow_champ_factor_een_grad + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_factor_een_grad_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 nuclei | + | ~cord_num~ | ~int64_t~ | in | order of polynomials | + | ~dim_c_vector~ | ~int64_t~ | in | dimension of full coefficient vector | + | ~c_vector_full~ | ~double[dim_c_vector][nucl_num]~ | in | full coefficient vector | + | ~lkpm_combined_index~ | ~int64_t[4][dim_c_vector]~ | in | combined indices | + | ~tmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][elec_num]~ | in | Temporary intermediate tensor | + | ~dtmp_c~ | ~double[walk_num][0:cord_num-1][0:cord_num][nucl_num][4][elec_num]~ | in | vector of non-zero coefficients | + | ~een_rescaled_n~ | ~double[walk_num][0:cord_num][nucl_num][elec_num]~ | in | Electron-nucleus rescaled factor | + | ~een_rescaled_n_gl~ | ~double[walk_num][0:cord_num][nucl_num][4][elec_num]~ | in | Derivative of Electron-nucleus rescaled factor | + | ~factor_een_grad~ | ~double[walk_num][4][elec_num]~ | out | Derivative of Electron-nucleus jastrow | + + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +function qmckl_compute_jastrow_champ_factor_een_grad_doc( & + context, walk_num, elec_num, nucl_num, & + cord_num, dim_c_vector, c_vector_full, lkpm_combined_index, & + tmp_c, dtmp_c, een_rescaled_n, een_rescaled_n_gl, factor_een_grad) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer (c_int64_t) , intent(in) , value :: walk_num + integer (c_int64_t) , intent(in) , value :: elec_num + integer (c_int64_t) , intent(in) , value :: nucl_num + integer (c_int64_t) , intent(in) , value :: cord_num + integer (c_int64_t) , intent(in) , value :: dim_c_vector + real (c_double ) , intent(in) :: c_vector_full(nucl_num,dim_c_vector) + integer (c_int64_t) , intent(in) :: lkpm_combined_index(dim_c_vector,4) + real (c_double ) , intent(in) :: tmp_c(elec_num,nucl_num,0:cord_num,0:cord_num-1,walk_num) + real (c_double ) , intent(in) :: dtmp_c(elec_num,4,nucl_num,0:cord_num,0:cord_num-1,walk_num) + real (c_double ) , intent(in) :: een_rescaled_n(elec_num,nucl_num,0:cord_num,walk_num) + real (c_double ) , intent(in) :: een_rescaled_n_gl(elec_num,4,nucl_num,0:cord_num,walk_num) + real (c_double ) , intent(out) :: factor_een_grad(elec_num,3,walk_num) + integer(qmckl_exit_code) :: info + + integer*8 :: i, a, j, l, k, m, n, nw, ii + double precision :: accu, accu2, cn + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT + if (walk_num <= 0) info = QMCKL_INVALID_ARG_2 + if (elec_num <= 0) info = QMCKL_INVALID_ARG_3 + if (nucl_num <= 0) info = QMCKL_INVALID_ARG_4 + if (cord_num < 0) info = QMCKL_INVALID_ARG_5 + if (info /= QMCKL_SUCCESS) return + + + if (cord_num == 0) then + factor_een_grad = 0.0d0 + return + end if + + do nw =1, walk_num + factor_een_grad(:,:,nw) = 0.0d0 + do n = 1, dim_c_vector + l = lkpm_combined_index(n, 1) + k = lkpm_combined_index(n, 2) + m = lkpm_combined_index(n, 4) + + do a = 1, nucl_num + cn = c_vector_full(a, n) + if(cn == 0.d0) cycle + + do ii = 1, 3 + do j = 1, elec_num + factor_een_grad(j,ii,nw) = factor_een_grad(j,ii,nw) + ( & + tmp_c(j,a,m,k,nw) * een_rescaled_n_gl(j,ii,a,m+l,nw) + & + (dtmp_c(j,ii,a,m,k,nw)) * een_rescaled_n(j,a,m+l,nw) + & + (dtmp_c(j,ii,a,m+l,k,nw)) * een_rescaled_n(j,a,m ,nw) + & + tmp_c(j,a,m+l,k,nw) * een_rescaled_n_gl(j,ii,a,m,nw) & + ) * cn + end do + end do + + end do + end do + end do + +end function qmckl_compute_jastrow_champ_factor_een_grad_doc + #+end_src + + #+CALL: generate_private_c_header(table=qmckl_factor_een_grad_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_factor_een_grad_doc" ) + + #+RESULTS: + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_grad_doc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_c_vector, + const double* c_vector_full, + const int64_t* lkpm_combined_index, + const double* tmp_c, + const double* dtmp_c, + const double* een_rescaled_n, + const double* een_rescaled_n_gl, + double* const factor_een_grad ); + #+end_src + + #+CALL: generate_private_c_header(table=qmckl_factor_een_grad_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_factor_een_grad" ) + + #+RESULTS: + #+begin_src c :tangle (eval h_private_func) :comments org + qmckl_exit_code qmckl_compute_jastrow_champ_factor_een_grad ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_c_vector, + const double* c_vector_full, + const int64_t* lkpm_combined_index, + const double* tmp_c, + const double* dtmp_c, + const double* een_rescaled_n, + const double* een_rescaled_n_gl, + double* const factor_een_grad ); + #+end_src + + #+begin_src c :tangle (eval c) :comments org +qmckl_exit_code +qmckl_compute_jastrow_champ_factor_een_grad(const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_c_vector, + const double *c_vector_full, + const int64_t *lkpm_combined_index, + const double *tmp_c, + const double *dtmp_c, + const double *een_rescaled_n, + const double *een_rescaled_n_gl, + double* const factor_een_grad) +{ +#ifdef HAVE_HPC + return qmckl_compute_jastrow_champ_factor_een_grad_hpc +#else + return qmckl_compute_jastrow_champ_factor_een_grad_doc +#endif + (context, walk_num, elec_num, nucl_num, + cord_num, dim_c_vector, c_vector_full, + lkpm_combined_index, tmp_c, dtmp_c, + een_rescaled_n, een_rescaled_n_gl, + factor_een_grad); +} +#+end_src +***** HPC implementation :noexport: + #+CALL: generate_private_c_header(table=qmckl_factor_een_grad_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_factor_een_grad_hpc" ) + + #+RESULTS: + #+begin_src c :tangle (eval h_private_func) :comments org +qmckl_exit_code +qmckl_compute_jastrow_champ_factor_een_grad_hpc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_c_vector, + const double* c_vector_full, + const int64_t* lkpm_combined_index, + const double* tmp_c, + const double* dtmp_c, + const double* een_rescaled_n, + const double* een_rescaled_n_gl, + double* const factor_een_grad ); + #+end_src + + #+begin_src c :tangle (eval c) :comments org +qmckl_exit_code +qmckl_compute_jastrow_champ_factor_een_grad_hpc(const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const int64_t dim_c_vector, + const double *c_vector_full, + const int64_t *lkpm_combined_index, + const double *tmp_c, + const double *dtmp_c, + const double *een_rescaled_n, + const double *een_rescaled_n_gl, + double* const factor_een_grad) +{ + + int64_t info = QMCKL_SUCCESS; + + if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT; + if (walk_num <= 0) return QMCKL_INVALID_ARG_2; + if (elec_num <= 0) return QMCKL_INVALID_ARG_3; + if (nucl_num <= 0) return QMCKL_INVALID_ARG_4; + if (cord_num < 0) return QMCKL_INVALID_ARG_5; + + + if (cord_num == 0) { +#ifdef HAVE_OPENMP +#pragma omp parallel for +#endif + for (size_t nw = 0; nw < (size_t) walk_num; ++nw) { + memset(&factor_een_grad[elec_num*3*nw], 0, elec_num*3*sizeof(double)); + } + return QMCKL_SUCCESS; + } + + const size_t elec_num2 = elec_num << 1; + +#ifdef HAVE_OPENMP +#pragma omp parallel for +#endif + for (size_t nw = 0; nw < (size_t) walk_num; ++nw) { + bool touched = false; + double* const restrict factor_een_grad_0nw = &(factor_een_grad[elec_num*3*nw]); + for (size_t n = 0; n < (size_t) dim_c_vector; ++n) { + const size_t l = lkpm_combined_index[n]; + const size_t k = lkpm_combined_index[n+ dim_c_vector]; + const size_t m = lkpm_combined_index[n+3*dim_c_vector]; + + const size_t en = elec_num*nucl_num; + const size_t len = l*en; + const size_t len4 = len << 2; + const size_t cn = cord_num*nw; + const size_t c1 = cord_num+1; + const size_t addr0 = en*(m+c1*(k+cn)); + const size_t addr1 = en*(m+cn); + + const double* restrict tmp_c_mkn = &(tmp_c[addr0]); + const double* restrict tmp_c_mlkn = tmp_c_mkn + len; + const double* restrict een_rescaled_n_mnw = &(een_rescaled_n[addr1]); + const double* restrict een_rescaled_n_mlnw = een_rescaled_n_mnw + len; + const double* restrict dtmp_c_mknw = &(dtmp_c[addr0 << 2]); + const double* restrict dtmp_c_mlknw = dtmp_c_mknw + len4; + const double* restrict een_rescaled_n_gl_mnw = &(een_rescaled_n_gl[addr1 << 2]); + const double* restrict een_rescaled_n_gl_mlnw = een_rescaled_n_gl_mnw + len4; + for (size_t a = 0; a < (size_t) nucl_num; a++) { + double cn = c_vector_full[a+n*nucl_num]; + if (cn == 0.0) continue; + + const size_t ishift = elec_num*a; + const size_t ishift4 = ishift << 2; + + const double* restrict tmp_c_amlkn = tmp_c_mlkn + ishift; + const double* restrict tmp_c_amkn = tmp_c_mkn + ishift; + const double* restrict een_rescaled_n_amnw = een_rescaled_n_mnw + ishift; + const double* restrict een_rescaled_n_amlnw = een_rescaled_n_mlnw + ishift; + const double* restrict dtmp_c_0amknw = dtmp_c_mknw + ishift4; + const double* restrict dtmp_c_0amlknw = dtmp_c_mlknw + ishift4; + const double* restrict een_rescaled_n_gl_0amnw = een_rescaled_n_gl_mnw + ishift4; + const double* restrict een_rescaled_n_gl_0amlnw = een_rescaled_n_gl_mlnw + ishift4; + + const double* restrict dtmp_c_1amknw = dtmp_c_0amknw + elec_num; + const double* restrict dtmp_c_1amlknw = dtmp_c_0amlknw + elec_num; + const double* restrict dtmp_c_2amknw = dtmp_c_0amknw + elec_num2; + const double* restrict dtmp_c_2amlknw = dtmp_c_0amlknw + elec_num2; + const double* restrict een_rescaled_n_gl_1amnw = een_rescaled_n_gl_0amnw + elec_num; + const double* restrict een_rescaled_n_gl_1amlnw = een_rescaled_n_gl_0amlnw + elec_num; + const double* restrict een_rescaled_n_gl_2amnw = een_rescaled_n_gl_0amnw + elec_num2; + const double* restrict een_rescaled_n_gl_2amlnw = een_rescaled_n_gl_0amlnw + elec_num2; + double* const restrict factor_een_grad_1nw = factor_een_grad_0nw + elec_num; + double* const restrict factor_een_grad_2nw = factor_een_grad_0nw + elec_num2; + + if (touched) { +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (size_t j = 0; j < (size_t) elec_num; ++j) { + factor_een_grad_0nw[j] = factor_een_grad_0nw[j] + cn * + (tmp_c_amkn[j] * een_rescaled_n_gl_0amlnw[j] + + dtmp_c_0amknw[j] * een_rescaled_n_amlnw[j] + + dtmp_c_0amlknw[j] * een_rescaled_n_amnw[j] + + tmp_c_amlkn[j] * een_rescaled_n_gl_0amnw[j]); + } + +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (size_t j = 0; j < (size_t) elec_num; ++j) { + factor_een_grad_1nw[j] = factor_een_grad_1nw[j] + cn * + (tmp_c_amkn[j] * een_rescaled_n_gl_1amlnw[j] + + dtmp_c_1amknw[j] * een_rescaled_n_amlnw[j] + + dtmp_c_1amlknw[j] * een_rescaled_n_amnw[j] + + tmp_c_amlkn[j] * een_rescaled_n_gl_1amnw[j]); + } + +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (size_t j = 0; j < (size_t) elec_num; ++j) { + factor_een_grad_2nw[j] = factor_een_grad_2nw[j] + cn * + (tmp_c_amkn[j] * een_rescaled_n_gl_2amlnw[j] + + dtmp_c_2amknw[j] * een_rescaled_n_amlnw[j] + + dtmp_c_2amlknw[j] * een_rescaled_n_amnw[j] + + tmp_c_amlkn[j] * een_rescaled_n_gl_2amnw[j]); + } + + } else { + + touched = true; + +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (size_t j = 0; j < (size_t) elec_num; ++j) { + factor_een_grad_0nw[j] = cn * + (tmp_c_amkn[j] * een_rescaled_n_gl_0amlnw[j] + + dtmp_c_0amknw[j] * een_rescaled_n_amlnw[j] + + dtmp_c_0amlknw[j] * een_rescaled_n_amnw[j] + + tmp_c_amlkn[j] * een_rescaled_n_gl_0amnw[j]); + } + +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (size_t j = 0; j < (size_t) elec_num; ++j) { + factor_een_grad_1nw[j] = cn * + (tmp_c_amkn[j] * een_rescaled_n_gl_1amlnw[j] + + dtmp_c_1amknw[j] * een_rescaled_n_amlnw[j] + + dtmp_c_1amlknw[j] * een_rescaled_n_amnw[j] + + tmp_c_amlkn[j] * een_rescaled_n_gl_1amnw[j]); + } + +#ifdef HAVE_OPENMP +#pragma omp simd +#endif + for (size_t j = 0; j < (size_t) elec_num; ++j) { + factor_een_grad_2nw[j] = cn * + (tmp_c_amkn[j] * een_rescaled_n_gl_2amlnw[j] + + dtmp_c_2amknw[j] * een_rescaled_n_amlnw[j] + + dtmp_c_2amlknw[j] * een_rescaled_n_amnw[j] + + tmp_c_amlkn[j] * een_rescaled_n_gl_2amnw[j]); + } + + } + } + } + if (!touched) { + memset(factor_een_grad_0nw, 0, elec_num*3*sizeof(double)); + } + } + return info; +} + #+end_src + +***** Test + + #+begin_src c :tangle (eval c_test) +/* Check if Jastrow is properly initialized */ +assert(qmckl_jastrow_champ_provided(context)); + +{ + double factor_een_gl[walk_num][4][elec_num]; + rc = qmckl_get_jastrow_champ_factor_een_gl(context, &(factor_een_gl[0][0][0]),4*walk_num*elec_num); + + double factor_een_grad[walk_num][3][elec_num]; + rc = qmckl_get_jastrow_champ_factor_een_grad(context, &(factor_een_grad[0][0][0]),3*walk_num*elec_num); + + for (int nw=0 ; nwjastrow_champ.gl, sze*sizeof(double)); + return QMCKL_SUCCESS; +} + #+end_src + + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code +qmckl_get_jastrow_champ_grad(qmckl_context context, + double* const grad, + const int64_t size_max); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code +qmckl_get_jastrow_champ_grad(qmckl_context context, + double* const grad, + const int64_t size_max) +{ + qmckl_exit_code rc; + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_get_jastrow_champ_grad", + NULL); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + rc = qmckl_provide_jastrow_champ_grad(context); + if (rc != QMCKL_SUCCESS) return rc; + + int64_t sze = 3 * ctx->electron.walker.num * ctx->electron.num; + if (size_max < sze) { + return qmckl_failwith( context, + QMCKL_INVALID_ARG_3, + "qmckl_get_jastrow_champ_grad", + "Array too small. Expected walker.num * electron.num * 3"); + } + memcpy(grad, ctx->jastrow_champ.grad, sze*sizeof(double)); + return QMCKL_SUCCESS; } #+end_src @@ -11826,6 +12412,7 @@ end interface **** Provide :noexport: #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_provide_jastrow_champ_gl(qmckl_context context); +qmckl_exit_code qmckl_provide_jastrow_champ_grad(qmckl_context context); #+end_src #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none @@ -11907,11 +12494,92 @@ qmckl_exit_code qmckl_provide_jastrow_champ_gl(qmckl_context context) ctx->jastrow_champ.gl_date = ctx->date; } + return QMCKL_SUCCESS; +} + +qmckl_exit_code qmckl_provide_jastrow_champ_grad(qmckl_context context) +{ + + qmckl_exit_code rc; + + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return qmckl_failwith( context, + QMCKL_INVALID_CONTEXT, + "qmckl_provide_jastrow_champ_grad", + NULL); + } + + qmckl_context_struct* const ctx = (qmckl_context_struct*) context; + assert (ctx != NULL); + + if (!ctx->jastrow_champ.provided) { + return qmckl_failwith( context, + QMCKL_NOT_PROVIDED, + "qmckl_provide_jastrow_champ_grad", + NULL); + } + + + rc = qmckl_provide_jastrow_champ_value(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_jastrow_champ_factor_ee_gl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_jastrow_champ_factor_en_gl(context); + if (rc != QMCKL_SUCCESS) return rc; + + rc = qmckl_provide_jastrow_champ_factor_een_grad(context); + if (rc != QMCKL_SUCCESS) return rc; + + /* Compute if necessary */ + if (ctx->date > ctx->jastrow_champ.grad_date) { + + if (ctx->electron.walker.num > ctx->electron.walker_old.num) { + if (ctx->jastrow_champ.grad != NULL) { + rc = qmckl_free(context, ctx->jastrow_champ.grad); + if (rc != QMCKL_SUCCESS) { + return qmckl_failwith( context, rc, + "qmckl_provide_jastrow_champ_grad", + "Unable to free ctx->jastrow_champ.grad"); + } + ctx->jastrow_champ.grad = NULL; + } + } + + /* Allocate array */ + if (ctx->jastrow_champ.grad == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->electron.walker.num * ctx->electron.num * 3 * sizeof(double); + double* grad = (double*) qmckl_malloc(context, mem_info); + + if (grad == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_provide_jastrow_champ_grad", + NULL); + } + ctx->jastrow_champ.grad = grad; + } + + rc = qmckl_compute_jastrow_champ_grad_doc(context, + ctx->electron.walker.num, + ctx->electron.num, + ctx->jastrow_champ.value, + ctx->jastrow_champ.factor_ee_gl, + ctx->jastrow_champ.factor_en_gl, + ctx->jastrow_champ.factor_een_grad, + ctx->jastrow_champ.grad); + + ctx->jastrow_champ.grad_date = ctx->date; + } + return QMCKL_SUCCESS; } #+end_src -**** Compute +**** Compute GL :PROPERTIES: :Name: qmckl_compute_jastrow_champ_gl_doc :CRetType: qmckl_exit_code @@ -12091,14 +12759,14 @@ qmckl_exit_code qmckl_compute_jastrow_champ_gl ( } #+end_src -**** Test :noexport: +***** Test - #+begin_src c :tangle (eval c_test) + #+begin_src c :tangle (eval c_test) printf("Total Jastrow derivatives\n"); /* Check if Jastrow is properly initialized */ assert(qmckl_jastrow_champ_provided(context)); - { + double factor_ee_gl[walk_num][4][elec_num]; rc = qmckl_check(context, qmckl_get_jastrow_champ_factor_ee_gl(context, &(factor_ee_gl[0][0][0]), walk_num*elec_num*4) @@ -12138,19 +12806,222 @@ assert(qmckl_jastrow_champ_provided(context)); printf("k = %ld\n", k); printf("m = %ld\n", m); printf("e = %ld\n", e); - printf("total_j_deriv = %20.15e\n", total_j_deriv[k][m][e]/total_j[k]); + printf("total_j_deriv/total_j = %20.15e\n", total_j_deriv[k][m][e]/total_j[k]); printf("factor_ee_gl + factor_en_gl + factor_een_gl = %20.15e\n", factor_ee_gl[k][m][e] + factor_en_gl[k][m][e] + factor_een_gl[k][m][e]); fflush(stdout); } - assert (fabs(total_j_deriv[k][m][e]/total_j[k] - (factor_ee_gl[k][m][e] + factor_en_gl[k][m][e] + factor_een_gl[k][m][e])) < 1.e-12); + assert (fabs(total_j_deriv[k][m][e]/total_j[k] - (factor_ee_gl[k][m][e] + factor_en_gl[k][m][e] + factor_een_gl[k][m][e])) < 1.e-8); } } } } + + #+end_src + +**** Compute Gradient only + :PROPERTIES: + :Name: qmckl_compute_jastrow_champ_grad_doc + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + #+NAME: qmckl_jastrow_champ_grad_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 | + | ~value~ | ~double[walk_num]~ | in | Total Jastrow | + | ~gl_ee~ | ~double[walk_num][4][elec_num]~ | in | ee component | + | ~gl_en~ | ~double[walk_num][4][elec_num]~ | in | eN component | + | ~grad_een~ | ~double[walk_num][3][elec_num]~ | in | eeN component | + | ~grad~ | ~double[walk_num][3][elec_num]~ | out | Total Jastrow factor | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +function qmckl_compute_jastrow_champ_grad_doc(context, & + walk_num, elec_num, value, gl_ee, gl_en, grad_een, grad) & + result(info) bind(C) + use, intrinsic :: iso_c_binding + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer (c_int64_t) , intent(in) , value :: walk_num + integer (c_int64_t) , intent(in) , value :: elec_num + real (c_double ) , intent(in) :: value(walk_num) + real (c_double ) , intent(in) :: gl_ee(elec_num,4,walk_num) + real (c_double ) , intent(in) :: gl_en(elec_num,4,walk_num) + real (c_double ) , intent(in) :: grad_een(elec_num,3,walk_num) + real (c_double ) , intent(out) :: grad(elec_num,3,walk_num) + + integer(qmckl_exit_code) :: info + integer*8 :: i, j, k + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (walk_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + do k = 1, walk_num + do j=1,3 + do i = 1, elec_num + grad(i,j,k) = gl_ee(i,j,k) + gl_en(i,j,k) + grad_een(i,j,k) + end do + end do + grad(:,:,k) = grad(:,:,k) * value(k) + end do + + +end function qmckl_compute_jastrow_champ_grad_doc + #+end_src + +#+CALL: generate_private_c_header(table=qmckl_jastrow_champ_grad_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_grad") + +#+RESULTS: +#+begin_src c :tangle (eval h_private_func) :comments org +qmckl_exit_code qmckl_compute_jastrow_champ_grad ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const double* value, + const double* gl_ee, + const double* gl_en, + const double* grad_een, + double* const grad ); +#+end_src + +#+CALL: generate_private_c_header(table=qmckl_jastrow_champ_grad_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_grad_doc") + +#+RESULTS: +#+begin_src c :tangle (eval h_private_func) :comments org +qmckl_exit_code qmckl_compute_jastrow_champ_grad_doc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const double* value, + const double* gl_ee, + const double* gl_en, + const double* grad_een, + double* const grad ); +#+end_src + +#+CALL: generate_private_c_header(table=qmckl_jastrow_champ_grad_args,rettyp=get_value("CRetType"),fname="qmckl_compute_jastrow_champ_grad_hpc") + +#+RESULTS: +#+begin_src c :tangle (eval h_private_func) :comments org +qmckl_exit_code qmckl_compute_jastrow_champ_grad_hpc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const double* value, + const double* gl_ee, + const double* gl_en, + const double* grad_een, + double* const grad ); +#+end_src + + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +inline qmckl_exit_code +qmckl_compute_jastrow_champ_grad_hpc ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const double* value, + const double* gl_ee, + const double* gl_en, + const double* grad_een, + double* const grad) +{ + + if (context == QMCKL_NULL_CONTEXT) return QMCKL_INVALID_CONTEXT; + if (walk_num <= 0 ) return QMCKL_INVALID_ARG_2; + if (elec_num <= 0 ) return QMCKL_INVALID_ARG_3; + if (value == NULL ) return QMCKL_INVALID_ARG_4; + if (gl_ee == NULL ) return QMCKL_INVALID_ARG_5; + if (gl_en == NULL ) return QMCKL_INVALID_ARG_6; + if (grad_een == NULL ) return QMCKL_INVALID_ARG_7; + if (grad == NULL ) return QMCKL_INVALID_ARG_8; + + for (int64_t k = 0; k < walk_num; ++k) { + for (int64_t j = 0; j < 3; ++j) { + for (int64_t i = 0; i < elec_num; ++i) { + grad[i + elec_num*(j + k*3)] = ( gl_ee[i + elec_num*(j + k*4)] + + gl_en[i + elec_num*(j + k*4)] + grad_een[i + elec_num*(j + k*3)] )* value[k]; + } + } + } + + return QMCKL_SUCCESS; +} + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes +qmckl_exit_code qmckl_compute_jastrow_champ_grad ( + const qmckl_context context, + const int64_t walk_num, + const int64_t elec_num, + const double* value, + const double* gl_ee, + const double* gl_en, + const double* grad_een, + double* const grad) +{ + +#ifdef HAVE_HPC + return qmckl_compute_jastrow_champ_grad_hpc +#else + return qmckl_compute_jastrow_champ_grad_doc +#endif + (context, walk_num, elec_num, value, gl_ee, gl_en, grad_een, grad); +} + #+end_src + +***** Test + + #+begin_src c :tangle (eval c_test) +printf("Total Jastrow gradient only\n"); +/* Check if Jastrow is properly initialized */ +assert(qmckl_jastrow_champ_provided(context)); + +{ + double total_j_grad[walk_num][3][elec_num]; + rc = qmckl_check(context, + qmckl_get_jastrow_champ_grad(context, &(total_j_grad[0][0][0]), walk_num*elec_num*3) + ); + assert(rc == QMCKL_SUCCESS); + + + double total_j_deriv[walk_num][4][elec_num]; + rc = qmckl_check(context, + qmckl_get_jastrow_champ_gl(context, &(total_j_deriv[0][0][0]), walk_num*elec_num*4) + ); + assert(rc == QMCKL_SUCCESS); + + + + for (int64_t k=0 ; k< walk_num ; ++k) { + for (int64_t m=0 ; m<3; ++m) { + for (int64_t e=0 ; e 1e-12) { + printf("%ld %ld %ld\n", k, m, e); + printf("total_j_grad = %20.15e\n", total_j_grad[k][m][e]); + } + assert (fabs(total_j_deriv[k][m][e] - total_j_grad[k][m][e]) < 1.e-12); + } + } + } } - #+end_src +} + #+end_src * End of files :noexport: