From 2438df280f66bf7198664b60ba58e45f39425a92 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Wed, 19 Feb 2025 17:20:42 +0100 Subject: [PATCH] Introduced qmckl_compute_forces_jastrow_delta_p_hpc --- org/qmckl_forces.org | 111 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 107 insertions(+), 4 deletions(-) diff --git a/org/qmckl_forces.org b/org/qmckl_forces.org index 47128c5..cab567d 100644 --- a/org/qmckl_forces.org +++ b/org/qmckl_forces.org @@ -4482,6 +4482,93 @@ integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_delta_p_doc( & end function qmckl_compute_forces_jastrow_delta_p_doc #+end_src + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer(qmckl_exit_code) function qmckl_compute_forces_jastrow_delta_p_hpc( & + context, num_in, walk_num, elec_num, nucl_num, cord_num, & + een_rescaled_n, een_rescaled_e, een_rescaled_single_n, een_rescaled_single_e, & + een_rescaled_n_gl, een_rescaled_single_n_gl, forces_delta_p) & + 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 :: num_in, walk_num, elec_num, cord_num, nucl_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_e(elec_num, elec_num, 0:cord_num, walk_num) + real(c_double) , intent(in) :: een_rescaled_single_n(nucl_num, 0:cord_num, walk_num) + real(c_double) , intent(in) :: een_rescaled_single_e(elec_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(in) :: een_rescaled_single_n_gl(4, nucl_num, 0:cord_num, walk_num) + real(c_double) , intent(out) :: forces_delta_p(elec_num,3,nucl_num,0:cord_num, 0:cord_num-1, walk_num) + + double precision :: tmp1, tmp2 + double precision :: delta_e(elec_num) + + integer*8 :: i, a, j, l, k, p, m, n, nw, num + double precision :: tmp, accu + integer :: elec_num4, size4 + + double precision, external :: ddot + double precision, allocatable :: een_rescaled_n_g(:,:,:,:) + + num = num_in + 1 + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) info = QMCKL_INVALID_CONTEXT + if (walk_num <= 0) info = QMCKL_INVALID_ARG_3 + if (elec_num <= 0) info = QMCKL_INVALID_ARG_4 + if (nucl_num <= 0) info = QMCKL_INVALID_ARG_5 + if (cord_num < 0) info = QMCKL_INVALID_ARG_6 + if (info /= QMCKL_SUCCESS) return + + + if (cord_num == 0) then + forces_delta_p = 0.d0 + return + endif + + + elec_num4 = elec_num + size4 = 3*nucl_num*cord_num + + allocate(een_rescaled_n_g(elec_num, 3, nucl_num, cord_num)) + + do nw=1, walk_num + do l=1,cord_num + do a=1,nucl_num + een_rescaled_n_g(:,1:3,a,l) = een_rescaled_n_gl(:,1:3,a,l,nw) + enddo + enddo + + do m=0, cord_num-1 + + do l=1, cord_num + do a = 1, nucl_num + do k = 1, 3 + tmp1 = een_rescaled_n_gl(num, k, a, l, nw) + tmp2 = een_rescaled_single_n_gl(k,a,l,nw) + do j = 1, elec_num + forces_delta_p(j,k,a,l,m,nw) = een_rescaled_e(j,num,m,nw) * tmp1 - & + een_rescaled_single_e(j,m,nw) * tmp2 + + end do + end do + end do + end do + + delta_e(:) = een_rescaled_e(:,num,m,nw) - een_rescaled_single_e(:,m,nw) + + call dgemv('T', elec_num4, size4, 1.d0, & + een_rescaled_n_g, elec_num4, delta_e, 1, 1.d0, & + forces_delta_p(num,1,1,1,m,nw), elec_num4) + + end do + end do + +end function qmckl_compute_forces_jastrow_delta_p_hpc + #+end_src + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none qmckl_exit_code qmckl_compute_forces_jastrow_delta_p_doc (const qmckl_context context, @@ -4499,7 +4586,21 @@ qmckl_compute_forces_jastrow_delta_p_doc (const qmckl_context context, double* const forces_delta_p ); qmckl_exit_code +qmckl_compute_forces_jastrow_delta_p_hpc (const qmckl_context context, + const int64_t num, + const int64_t walk_num, + const int64_t elec_num, + const int64_t nucl_num, + const int64_t cord_num, + const double* een_rescaled_n, + const double* een_rescaled_e, + const double* een_rescaled_single_n, + const double* een_rescaled_single_e, + const double* een_rescaled_n_gl, + const double* een_rescaled_single_n_gl, + double* const forces_delta_p ); +qmckl_exit_code qmckl_compute_forces_jastrow_delta_p (const qmckl_context context, const int64_t num, const int64_t walk_num, @@ -4532,7 +4633,7 @@ qmckl_compute_forces_jastrow_delta_p (const qmckl_context context, double* const forces_delta_p ) { #ifdef HAVE_HPC - return qmckl_compute_forces_jastrow_delta_p_doc + return qmckl_compute_forces_jastrow_delta_p_hpc #else return qmckl_compute_forces_jastrow_delta_p_doc #endif @@ -5001,9 +5102,10 @@ function qmckl_compute_forces_jastrow_single_een_hpc( & integer*8 :: i, a, j, l, k, p, m, n, nw, num, kk double precision :: accu2, cn double precision, allocatable :: accu(:,:), tmp(:) - integer*8 :: LDA, LDB, LDC double precision, external :: ddot + integer :: elec_num4 + num = num_in + 1 info = QMCKL_SUCCESS @@ -5022,6 +5124,7 @@ function qmckl_compute_forces_jastrow_single_een_hpc( & allocate(een_rescaled_delta_n(nucl_num, 0:cord_num), een_rescaled_delta_n_gl(3,nucl_num, 0:cord_num), & accu(3,nucl_num), tmp(nucl_num)) + elec_num4 = elec_num do nw =1, walk_num een_rescaled_delta_n(:,:) = een_rescaled_single_n(:,:,nw) - een_rescaled_n(num,:,:,nw) een_rescaled_delta_n_gl(1:3,:,:) = een_rescaled_single_n_gl(1:3,:,:,nw) - een_rescaled_n_gl(num,1:3,:,:,nw) @@ -5036,9 +5139,9 @@ function qmckl_compute_forces_jastrow_single_een_hpc( & accu(1:3,a) = 0.0d0 if(cn == 0.d0) cycle tmp(a) = tmp_c(num,a,m+l,k,nw) + delta_p(num,a,m+l,k,nw) - call dgemv('T', elec_num, 3, -1.d0, een_rescaled_n_gl(1,1,a,m,nw), elec_num, & + call dgemv('T', elec_num4, 3, -1.d0, een_rescaled_n_gl(1,1,a,m,nw), elec_num4, & delta_p(1,a,m+l,k,nw), 1, 0.d0, accu(1,a), 1) - call dgemv('T', elec_num, 3, 1.d0, forces_delta_p(1,1,a,m+l,k,nw), elec_num, & + call dgemv('T', elec_num4, 3, 1.d0, forces_delta_p(1,1,a,m+l,k,nw), elec_num4, & een_rescaled_n(1,a,m,nw), 1, 1.d0, accu(1,a), 1) enddo