From bb155c0dfd0b12dadf34b321f0aaa6eb25ac3055 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Fri, 15 Sep 2023 11:30:10 +0200 Subject: [PATCH] J_qmckl en & ee --- src/non_h_ints_mu/NEED | 1 + src/non_h_ints_mu/jast_deriv.irp.f | 95 ++++++++++++++++++++++++++ src/non_h_ints_mu/qmckl.irp.f | 102 ++++++++++++++++++++++++++++ src/non_h_ints_mu/tc_integ_an.irp.f | 14 ++-- src/qmckl/README.md | 4 ++ src/qmckl/qmckl.F90 | 1 + 6 files changed, 213 insertions(+), 4 deletions(-) create mode 100644 src/non_h_ints_mu/qmckl.irp.f create mode 100644 src/qmckl/README.md create mode 100644 src/qmckl/qmckl.F90 diff --git a/src/non_h_ints_mu/NEED b/src/non_h_ints_mu/NEED index d09ab4a5..ecde6390 100644 --- a/src/non_h_ints_mu/NEED +++ b/src/non_h_ints_mu/NEED @@ -1,2 +1,3 @@ +qmckl ao_tc_eff_map bi_ortho_mos diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f index ee01886c..6b8445b1 100644 --- a/src/non_h_ints_mu/jast_deriv.irp.f +++ b/src/non_h_ints_mu/jast_deriv.irp.f @@ -24,11 +24,15 @@ double precision :: v1b_r1, v1b_r2, u2b_r12 double precision :: grad1_v1b(3), grad1_u2b(3) double precision :: dx, dy, dz + double precision :: time0, time1 double precision, external :: j12_mu, j1b_nucl PROVIDE j1b_type PROVIDE final_grid_points_extra + print*, ' providing grad1_u12_num & grad1_u12_squared_num ...' + call wall_time(time0) + grad1_u12_num = 0.d0 grad1_u12_squared_num = 0.d0 @@ -112,6 +116,94 @@ !$OMP END DO !$OMP END PARALLEL + elseif (j1b_type .eq. 1000) then + + double precision :: f + f = 1.d0 / dble(elec_num - 1) + + double precision, allocatable :: rij(:,:,:) + allocate( rij(3, 2, n_points_extra_final_grid) ) + + use qmckl + integer(qmckl_exit_code) :: rc + + integer*8 :: npoints + npoints = n_points_extra_final_grid + + double precision, allocatable :: gl(:,:,:) + allocate( gl(2,4,n_points_extra_final_grid) ) + + do ipoint = 1, n_points_final_grid ! r1 + + do jpoint = 1, n_points_extra_final_grid ! r2 + rij(1:3, 1, jpoint) = final_grid_points (1:3, ipoint) + rij(1:3, 2, jpoint) = final_grid_points_extra(1:3, jpoint) + enddo + + + rc = qmckl_set_electron_coord(qmckl_ctx_jastrow, 'N', npoints, rij, npoints*6_8) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in set_electron_coord' + stop -1 + endif + + + ! --- + ! e-e term + + rc = qmckl_get_jastrow_champ_factor_ee_gl(qmckl_ctx_jastrow, gl, 8_8*npoints) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in fact_ee_gl' + stop -1 + endif + + do jpoint = 1, n_points_extra_final_grid ! r2 + grad1_u12_num(jpoint,ipoint,1) = gl(1,1,jpoint) + grad1_u12_num(jpoint,ipoint,2) = gl(1,2,jpoint) + grad1_u12_num(jpoint,ipoint,3) = gl(1,3,jpoint) + enddo + + ! --- + ! e-e-n term + +! rc = qmckl_get_jastrow_champ_factor_een_gl(qmckl_ctx_jastrow, gl, 8_8*npoints) +! if (rc /= QMCKL_SUCCESS) then +! print *, irp_here, 'qmckl error in fact_een_gl' +! stop -1 +! endif +! +! do jpoint = 1, n_points_extra_final_grid ! r2 +! grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + gl(1,1,jpoint) +! grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + gl(1,2,jpoint) +! grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + gl(1,3,jpoint) +! enddo + + ! --- + ! e-n term + + rc = qmckl_get_jastrow_champ_factor_en_gl(qmckl_ctx_jastrow, gl, 8_8*npoints) + if (rc /= QMCKL_SUCCESS) then + print *, irp_here, 'qmckl error in fact_en_gl' + stop -1 + endif + + do jpoint = 1, n_points_extra_final_grid ! r2 + grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + f * gl(1,1,jpoint) + grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + f * gl(1,2,jpoint) + grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + f * gl(1,3,jpoint) + enddo + + do jpoint = 1, n_points_extra_final_grid ! r2 + dx = grad1_u12_num(jpoint,ipoint,1) + dy = grad1_u12_num(jpoint,ipoint,2) + dz = grad1_u12_num(jpoint,ipoint,3) + grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz + enddo + + enddo + + deallocate(gl, rij) + else print *, ' j1b_type = ', j1b_type, 'not implemented yet' @@ -119,6 +211,9 @@ endif + call wall_time(time1) + print*, ' Wall time for grad1_u12_num & grad1_u12_squared_num (min) =', (time1-time0)/60.d0 + END_PROVIDER ! --- diff --git a/src/non_h_ints_mu/qmckl.irp.f b/src/non_h_ints_mu/qmckl.irp.f new file mode 100644 index 00000000..d83de4dc --- /dev/null +++ b/src/non_h_ints_mu/qmckl.irp.f @@ -0,0 +1,102 @@ +BEGIN_PROVIDER [ integer*8, qmckl_ctx_jastrow ] + use qmckl + implicit none + BEGIN_DOC + ! Context for the QMCKL library + END_DOC + integer(qmckl_exit_code) :: rc + + qmckl_ctx_jastrow = qmckl_context_create() + + rc = qmckl_set_nucleus_num(qmckl_ctx_jastrow, nucl_num*1_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_nucleus_charge(qmckl_ctx_jastrow, nucl_charge, nucl_num*1_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_nucleus_coord(qmckl_ctx_jastrow, 'T', nucl_coord, nucl_num*3_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_electron_num(qmckl_ctx_jastrow, 1_8, 1_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + + ! Jastrow parameters + rc = qmckl_set_jastrow_champ_type_nucl_num (qmckl_ctx_jastrow, 2_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_jastrow_champ_type_nucl_vector (qmckl_ctx_jastrow, (/0_8,1_8,1_8/), 1_8*nucl_num) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_jastrow_champ_rescale_factor_ee (qmckl_ctx_jastrow, 0.6d0) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_jastrow_champ_rescale_factor_en (qmckl_ctx_jastrow, (/0.6d0, 0.6d0 /), 2_8 ) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_jastrow_champ_aord_num (qmckl_ctx_jastrow, 5_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_jastrow_champ_bord_num (qmckl_ctx_jastrow, 5_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_jastrow_champ_cord_num (qmckl_ctx_jastrow, 0_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + +! double precision :: a_vector(12) = dble(& +! (/ 0.00000000, 0.00000000, -0.71168405, -0.44415699, -0.13865109, 0.07002267 , & +! 0.00000000, 0.00000000, -0.11379992, 0.04542846, 0.01696997, -0.01809299 /) ) + +! double precision :: b_vector(6) = dble(& +! (/ 0.00000000, 0.65603311, 0.14581988, 0.03138163, 0.00153156, -0.00447302 /) ) + +! double precision :: c_vector(46) = & +! (/ 1.06384279d0, -1.44303973d0, -0.92409833d0, 0.11845356d0, -0.02980776d0, & +! 1.07048863d0, 0.06009623d0, -0.01854872d0, -0.00915398d0, 0.01324198d0, & +! -0.00504959d0, -0.01202497d0, -0.00531644d0, 0.15101629d0, -0.00723831d0, & +! -0.00384182d0, -0.00295036d0, -0.00114583d0, 0.00158107d0, -0.00078107d0, & +! -0.00080000d0, -0.14140576d0, -0.00237271d0, -0.03006706d0, 0.01537009d0, & +! -0.02327226d0, 0.16502789d0, -0.01458259d0, -0.09946065d0, 0.00850029d0, & +! -0.02969361d0, -0.01159547d0, 0.00516313d0, 0.00405247d0, -0.02200886d0, & +! 0.03376709d0, 0.01277767d0, -0.01523013d0, -0.00739224d0, -0.00463953d0, & +! 0.00003174d0, -0.01421128d0, 0.00808140d0, 0.00612988d0, -0.00610632d0, & +! 0.01926215d0 /) + +! a_vector = 0.d0 +! b_vector = 0.d0 +! c_vector = 0.d0 + + double precision :: a_vector(12) = dble(& + (/ 0.00000000 , 0.00000000, -0.45105821, -0.23519218, -0.03825391, 0.10072866, & + 0.00000000 , 0.00000000, -0.06930592, -0.02909224, -0.00134650, 0.01477242 /) ) + + double precision :: b_vector(6) = dble(& + (/ 0.00000000, 0.00000000, 0.29217862, -0.00450671, -0.02925982, -0.01381532 /) ) + + double precision :: c_vector(46) + c_vector = 0.d0 + + rc = qmckl_set_jastrow_champ_a_vector(qmckl_ctx_jastrow, a_vector, 12_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + + rc = qmckl_set_jastrow_champ_b_vector(qmckl_ctx_jastrow, b_vector, 6_8) + rc = qmckl_check(qmckl_ctx_jastrow, rc) + if (rc /= QMCKL_SUCCESS) stop -1 + +! rc = qmckl_set_jastrow_champ_c_vector(qmckl_ctx_jastrow, c_vector, 46_8) +! rc = qmckl_check(qmckl_ctx_jastrow, rc) +! if (rc /= QMCKL_SUCCESS) stop -1 + +END_PROVIDER diff --git a/src/non_h_ints_mu/tc_integ_an.irp.f b/src/non_h_ints_mu/tc_integ_an.irp.f index ae7af987..a6459761 100644 --- a/src/non_h_ints_mu/tc_integ_an.irp.f +++ b/src/non_h_ints_mu/tc_integ_an.irp.f @@ -106,8 +106,11 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f elseif(j1b_type .ge. 100) then - PROVIDE int2_grad1_u12_ao_num - int2_grad1_u12_ao = int2_grad1_u12_ao_num +! PROVIDE int2_grad1_u12_ao_num +! int2_grad1_u12_ao = int2_grad1_u12_ao_num + + PROVIDE int2_grad1_u12_ao_num_1shot + int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot else @@ -222,8 +225,11 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p elseif(j1b_type .ge. 100) then - PROVIDE int2_grad1_u12_square_ao_num - int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num + ! PROVIDE int2_grad1_u12_square_ao_num + ! int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num + + PROVIDE int2_grad1_u12_square_ao_num_1shot + int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num_1shot else diff --git a/src/qmckl/README.md b/src/qmckl/README.md new file mode 100644 index 00000000..ebc4b089 --- /dev/null +++ b/src/qmckl/README.md @@ -0,0 +1,4 @@ +#QMCkl + +Info related to the QMCkl library. + diff --git a/src/qmckl/qmckl.F90 b/src/qmckl/qmckl.F90 new file mode 100644 index 00000000..94ac962f --- /dev/null +++ b/src/qmckl/qmckl.F90 @@ -0,0 +1 @@ +#include