From 9dc8c0653d4f74aa3165884d4996111c93519bbb Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 15 Feb 2024 20:37:56 +0100 Subject: [PATCH] added Boys & Handy's Jastrow --- plugins/local/jastrow/EZFIO.cfg | 46 +++- plugins/local/jastrow/bh_param.irp.f | 252 ++++++++++++++++++ .../non_h_ints_mu/jast_deriv_utils_vect.irp.f | 96 +++++++ 3 files changed, 392 insertions(+), 2 deletions(-) create mode 100644 plugins/local/jastrow/bh_param.irp.f diff --git a/plugins/local/jastrow/EZFIO.cfg b/plugins/local/jastrow/EZFIO.cfg index 23dde8ea..8fd2d05a 100644 --- a/plugins/local/jastrow/EZFIO.cfg +++ b/plugins/local/jastrow/EZFIO.cfg @@ -1,13 +1,13 @@ [j2e_type] type: character*(32) -doc: type of the 2e-Jastrow: [ None | Mu | Mur | Qmckl ] +doc: type of the 2e-Jastrow: [ None | Mu | Mu_Nu | Mur | Boys | Boys_Handy | Qmckl ] interface: ezfio,provider,ocaml default: Mu [j1e_type] type: character*(32) -doc: type of the 1e-Jastrow: [ None | Gauss | Charge_Harmonizer ] +doc: type of the 1e-Jastrow: [ None | Gauss | Charge_Harmonizer | Charge_Harmonizer_AO ] interface: ezfio,provider,ocaml default: None @@ -151,3 +151,45 @@ interface: ezfio,provider,ocaml default: 1.0 ezfio_name: nu_erf +[jBH_size] +type: integer +doc: number of terms per atom in Boys-Handy-Jastrow +interface: ezfio,provider,ocaml +default: 1 + +[jBH_c] +type: double precision +doc: coefficients of terms in Boys-Handy-Jastrow +interface: ezfio +size: (jastrow.jBH_size,nuclei.nucl_num) + +[jBH_m] +type: integer +doc: powers of terms in Boys-Handy-Jastrow +interface: ezfio +size: (jastrow.jBH_size,nuclei.nucl_num) + +[jBH_n] +type: integer +doc: powers of terms in Boys-Handy-Jastrow +interface: ezfio +size: (jastrow.jBH_size,nuclei.nucl_num) + +[jBH_o] +type: integer +doc: powers of terms in Boys-Handy-Jastrow +interface: ezfio +size: (jastrow.jBH_size,nuclei.nucl_num) + +[jBH_ee] +type: double precision +doc: parameters of e-e terms in Boys-Handy-Jastrow +interface: ezfio +size: (nuclei.nucl_num) + +[jBH_en] +type: double precision +doc: parameters of e-n terms in Boys-Handy-Jastrow +interface: ezfio +size: (nuclei.nucl_num) + diff --git a/plugins/local/jastrow/bh_param.irp.f b/plugins/local/jastrow/bh_param.irp.f new file mode 100644 index 00000000..790cf97c --- /dev/null +++ b/plugins/local/jastrow/bh_param.irp.f @@ -0,0 +1,252 @@ + + BEGIN_PROVIDER [double precision, jBH_ee, (nucl_num)] +&BEGIN_PROVIDER [double precision, jBH_en, (nucl_num)] +&BEGIN_PROVIDER [double precision, jBH_c , (jBH_size, nucl_num)] +&BEGIN_PROVIDER [integer , jBH_m , (jBH_size, nucl_num)] +&BEGIN_PROVIDER [integer , jBH_n , (jBH_size, nucl_num)] +&BEGIN_PROVIDER [integer , jBH_o , (jBH_size, nucl_num)] + + BEGIN_DOC + ! + ! parameters of Boys-Handy-Jastrow + ! + END_DOC + + implicit none + logical :: exists + integer :: i_nucl, p + integer :: ierr + + PROVIDE ezfio_filename + + ! --- + + if(mpi_master) then + call ezfio_has_jastrow_jBH_ee(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + include 'mpif.h' + call MPI_BCAST(jBH_ee, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read Boys-Handy e-e param with MPI' + endif + IRP_ENDIF + + if(exists) then + if(mpi_master) then + write(6,'(A)') '.. >>>>> [ IO READ: jBH_ee ] <<<<< ..' + call ezfio_get_jastrow_jBH_ee(jBH_ee) + IRP_IF MPI + call MPI_BCAST(jBH_ee, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read jBH_ee with MPI' + endif + IRP_ENDIF + endif + else + + jBH_ee = 1.d0 + call ezfio_set_jastrow_jBH_ee(jBH_ee) + endif + + ! --- + + if(mpi_master) then + call ezfio_has_jastrow_jBH_en(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + call MPI_BCAST(jBH_en, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read Boys-Handy e-n param with MPI' + endif + IRP_ENDIF + + if(exists) then + if(mpi_master) then + write(6,'(A)') '.. >>>>> [ IO READ: jBH_en ] <<<<< ..' + call ezfio_get_jastrow_jBH_en(jBH_en) + IRP_IF MPI + call MPI_BCAST(jBH_en, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read jBH_en with MPI' + endif + IRP_ENDIF + endif + else + + jBH_en = 1.d0 + call ezfio_set_jastrow_jBH_en(jBH_en) + endif + + ! --- + + if(mpi_master) then + call ezfio_has_jastrow_jBH_c(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + call MPI_BCAST(jBH_c, (jBH_size*nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read Boys-Handy coeff with MPI' + endif + IRP_ENDIF + + if(exists) then + if(mpi_master) then + write(6,'(A)') '.. >>>>> [ IO READ: jBH_c ] <<<<< ..' + call ezfio_get_jastrow_jBH_c(jBH_c) + IRP_IF MPI + call MPI_BCAST(jBH_c, (jBH_size*nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read jBH_c with MPI' + endif + IRP_ENDIF + endif + else + + jBH_c = 0.d0 + call ezfio_set_jastrow_jBH_c(jBH_c) + endif + + ! --- + + if(mpi_master) then + call ezfio_has_jastrow_jBH_m(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + call MPI_BCAST(jBH_m, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read Boys-Handy m powers with MPI' + endif + IRP_ENDIF + + if(exists) then + if(mpi_master) then + write(6,'(A)') '.. >>>>> [ IO READ: jBH_m ] <<<<< ..' + call ezfio_get_jastrow_jBH_m(jBH_m) + IRP_IF MPI + call MPI_BCAST(jBH_m, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read jBH_m with MPI' + endif + IRP_ENDIF + endif + else + + jBH_m = 0 + call ezfio_set_jastrow_jBH_m(jBH_m) + endif + + ! --- + + if(mpi_master) then + call ezfio_has_jastrow_jBH_n(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + call MPI_BCAST(jBH_n, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read Boys-Handy n powers with MPI' + endif + IRP_ENDIF + + if(exists) then + if(mpi_master) then + write(6,'(A)') '.. >>>>> [ IO READ: jBH_n ] <<<<< ..' + call ezfio_get_jastrow_jBH_n(jBH_n) + IRP_IF MPI + call MPI_BCAST(jBH_n, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read jBH_n with MPI' + endif + IRP_ENDIF + endif + else + + jBH_n = 0 + call ezfio_set_jastrow_jBH_n(jBH_n) + endif + + ! --- + + if(mpi_master) then + call ezfio_has_jastrow_jBH_o(exists) + endif + + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + + IRP_IF MPI + call MPI_BCAST(jBH_o, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read Boys-Handy o powers with MPI' + endif + IRP_ENDIF + + if(exists) then + if(mpi_master) then + write(6,'(A)') '.. >>>>> [ IO READ: jBH_o ] <<<<< ..' + call ezfio_get_jastrow_jBH_o(jBH_o) + IRP_IF MPI + call MPI_BCAST(jBH_o, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if(ierr /= MPI_SUCCESS) then + stop 'Unable to read jBH_o with MPI' + endif + IRP_ENDIF + endif + else + + jBH_o = 0 + call ezfio_set_jastrow_jBH_o(jBH_o) + endif + + ! --- + + print *, ' parameters for Boys-Handy Jastrow' + print *, ' nb of terms per nucleus = ', jBH_size + + do i_nucl = 1, nucl_num + print *, ' i_nucl = ', i_nucl + print *, ' ee-term = ', jBH_ee(i_nucl) + print *, ' en-term = ', jBH_en(i_nucl) + print *, 'm n o c' + do p = 1, jBH_size + write(*,'(3(I4,2x), E15.7)') jBH_m(p,i_nucl), jBH_n(p,i_nucl), jBH_o(p,i_nucl), jBH_c(p,i_nucl) + enddo + enddo + + +END_PROVIDER + +! --- + diff --git a/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f b/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f index 5777a44a..88778ee0 100644 --- a/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f +++ b/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f @@ -109,6 +109,16 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res) endif ! env_type + elseif(j2e_type .eq. "Boys_Handy") then + + PROVIDE jBH_size jBH_en jBH_ee jBH_m jBH_n jBH_o jBH_c + + if(env_type .ne. "None") then + + call grad1_j12_r1_seq(r1, n_grid2, resx, resy, resz) + + endif ! env_type + else print *, ' Error in get_grad1_u12_withsq_r1_seq: Unknown Jastrow' @@ -157,9 +167,13 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) double precision, intent(out) :: gradz(n_grid2) integer :: jpoint + integer :: i_nucl, p, mpA, npA, opA double precision :: r2(3) double precision :: dx, dy, dz, r12, tmp double precision :: mu_val, mu_tmp, mu_der(3) + double precision :: rn(3), f1A, gard1_f1A(3), f2A, gard2_f2A(3), g12, gard1_g12(3) + double precision :: tmp1, tmp2 + PROVIDE j2e_type @@ -267,6 +281,57 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz) gradz(jpoint) = tmp * dz enddo + elseif(j2e_type .eq. "Boys_Handy") then + + do jpoint = 1, n_points_extra_final_grid ! r2 + + r2(1) = final_grid_points_extra(1,jpoint) + r2(2) = final_grid_points_extra(2,jpoint) + r2(3) = final_grid_points_extra(3,jpoint) + + gradx(jpoint) = 0.d0 + grady(jpoint) = 0.d0 + gradz(jpoint) = 0.d0 + do i_nucl = 1, nucl_num + + rn(1) = nucl_coord(i_nucl,1) + rn(2) = nucl_coord(i_nucl,2) + rn(3) = nucl_coord(i_nucl,3) + + call jBH_elem_fct_grad(jBH_en(i_nucl), r1, rn, f1A, gard1_f1A) + call jBH_elem_fct_grad(jBH_en(i_nucl), r2, rn, f2A, gard2_f2A) + call jBH_elem_fct_grad(jBH_ee(i_nucl), r1, r2, g12, gard1_g12) + + do p = 1, jBH_size + mpA = jBH_m(p,i_nucl) + npA = jBH_n(p,i_nucl) + opA = jBH_o(p,i_nucl) + tmp = jBH_c(p,i_nucl) + if(mpA .eq. npA) then + tmp = tmp * 0.5d0 + endif + + tmp1 = 0.d0 + if(mpA .gt. 0) then + tmp1 = tmp1 + dble(mpA) * f1A**dble(mpA-1) * f2A**dble(npA) + endif + if(npA .gt. 0) then + tmp1 = tmp1 + dble(npA) * f1A**dble(npA-1) * f2A**dble(mpA) + endif + tmp1 = tmp1 * g12**dble(opA) + + tmp2 = 0.d0 + if(opA .gt. 0) then + tmp2 = tmp2 + dble(opA) * g12**dble(opA-1) * (f1A**dble(mpA) * f2A**dble(npA) + f1A**dble(npA) * f2A**dble(mpA)) + endif + + gradx(jpoint) = gradx(jpoint) + tmp * (tmp1 * gard1_f1A(1) + tmp2 * gard1_g12(1)) + grady(jpoint) = grady(jpoint) + tmp * (tmp1 * gard1_f1A(2) + tmp2 * gard1_g12(2)) + gradz(jpoint) = gradz(jpoint) + tmp * (tmp1 * gard1_f1A(3) + tmp2 * gard1_g12(3)) + enddo ! p + enddo ! i_nucl + enddo ! jpoint + else print *, ' Error in grad1_j12_r1_seq: Unknown j2e_type = ', j2e_type @@ -757,3 +822,34 @@ end ! --- +subroutine jBH_elem_fct_grad(alpha, r1, r2, fct, gard1_fct) + + implicit none + double precision, intent(in) :: alpha, r1(3), r2(3) + double precision, intent(out) :: fct, gard1_fct(3) + double precision :: dist, tmp1, tmp2 + + dist = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) & + + (r1(2) - r2(2)) * (r1(2) - r2(2)) & + + (r1(3) - r2(3)) * (r1(3) - r2(3)) ) + + tmp1 = 1.d0 / (1.d0 + alpha * dist) + + fct = alpha * dist * tmp1 + + if(dist .lt. 1d-10) then + gard1_fct(1) = 0.d0 + gard1_fct(2) = 0.d0 + gard1_fct(3) = 0.d0 + else + tmp2 = alpha * tmp1 * tmp1 / dist + gard1_fct(1) = tmp2 * (r1(1) - r2(1)) + gard1_fct(2) = tmp2 * (r1(2) - r2(2)) + gard1_fct(3) = tmp2 * (r1(3) - r2(3)) + endif + + return +end + +! --- +