From 9dc8c0653d4f74aa3165884d4996111c93519bbb Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Thu, 15 Feb 2024 20:37:56 +0100 Subject: [PATCH 1/5] 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 + +! --- + From 6fa207a9fbebf79d3f0e36f2196bc3260878d45d Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Fri, 16 Feb 2024 16:45:54 +0100 Subject: [PATCH 2/5] Boys & Handy Jastrow: OK --- .../non_h_ints_mu/jast_deriv_utils_vect.irp.f | 4 +--- plugins/local/tc_scf/rh_tcscf_diis.irp.f | 23 ++++++------------- 2 files changed, 8 insertions(+), 19 deletions(-) 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 88778ee0..31ad5756 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 @@ -113,10 +113,8 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res) PROVIDE jBH_size jBH_en jBH_ee jBH_m jBH_n jBH_o jBH_c - if(env_type .ne. "None") then - + if(env_type .eq. "None") then call grad1_j12_r1_seq(r1, n_grid2, resx, resy, resz) - endif ! env_type else diff --git a/plugins/local/tc_scf/rh_tcscf_diis.irp.f b/plugins/local/tc_scf/rh_tcscf_diis.irp.f index 12678500..431b6e08 100644 --- a/plugins/local/tc_scf/rh_tcscf_diis.irp.f +++ b/plugins/local/tc_scf/rh_tcscf_diis.irp.f @@ -22,6 +22,9 @@ subroutine rh_tcscf_diis() logical, external :: qp_stop + PROVIDE level_shift_TCSCF + PROVIDE mo_l_coef mo_r_coef + it = 0 e_save = 0.d0 dim_DIIS = 0 @@ -41,19 +44,6 @@ subroutine rh_tcscf_diis() ! --- - PROVIDE level_shift_TCSCF - PROVIDE mo_l_coef mo_r_coef - - !write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & - ! '====', '================', '================', '================', '================', '================' & - ! , '================', '================', '================', '====', '========' - !write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & - ! ' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' & - ! , ' gradient ', ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)' - !write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & - ! '====', '================', '================', '================', '================', '================' & - ! , '================', '================', '================', '====', '========' - write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & '====', '================', '================', '================', '================', '================' & , '================', '================', '====', '========' @@ -81,8 +71,6 @@ subroutine rh_tcscf_diis() er_save = er_DIIS call wall_time(t1) - !write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & - ! it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 @@ -91,6 +79,8 @@ subroutine rh_tcscf_diis() PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot converged = .false. + call ezfio_set_tc_scf_converged_tcscf(converged) + !do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. dsqrt(thresh_tcscf))) do while(.not. converged) @@ -253,8 +243,9 @@ subroutine rh_tcscf_diis() endif call lock_io - if (converged) then + if(converged) then write(json_unit, json_true_fmtx) 'converged' + call ezfio_set_tc_scf_converged_tcscf(converged) else write(json_unit, json_false_fmtx) 'converged' endif From fcde51ea002e01e6b2d109c17320b2a877673845 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Thu, 22 Feb 2024 08:15:01 +0100 Subject: [PATCH 3/5] small modif --- plugins/local/jastrow/bh_param.irp.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/plugins/local/jastrow/bh_param.irp.f b/plugins/local/jastrow/bh_param.irp.f index 790cf97c..167d8814 100644 --- a/plugins/local/jastrow/bh_param.irp.f +++ b/plugins/local/jastrow/bh_param.irp.f @@ -236,7 +236,7 @@ print *, ' nb of terms per nucleus = ', jBH_size do i_nucl = 1, nucl_num - print *, ' i_nucl = ', i_nucl + print *, ' nucl = ', nucl_label(i_nucl) print *, ' ee-term = ', jBH_ee(i_nucl) print *, ' en-term = ', jBH_en(i_nucl) print *, 'm n o c' From 7fc4c8148552105b08088b2e3c79c59f48c86f50 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Fri, 23 Feb 2024 23:51:03 +0100 Subject: [PATCH 4/5] few modif --- plugins/local/jastrow/bh_param.irp.f | 2 +- plugins/local/tc_scf/EZFIO.cfg | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/plugins/local/jastrow/bh_param.irp.f b/plugins/local/jastrow/bh_param.irp.f index 167d8814..1ed871bc 100644 --- a/plugins/local/jastrow/bh_param.irp.f +++ b/plugins/local/jastrow/bh_param.irp.f @@ -239,7 +239,7 @@ print *, ' nucl = ', nucl_label(i_nucl) print *, ' ee-term = ', jBH_ee(i_nucl) print *, ' en-term = ', jBH_en(i_nucl) - print *, 'm n o c' + 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 diff --git a/plugins/local/tc_scf/EZFIO.cfg b/plugins/local/tc_scf/EZFIO.cfg index 313d6f2b..3dfa9a71 100644 --- a/plugins/local/tc_scf/EZFIO.cfg +++ b/plugins/local/tc_scf/EZFIO.cfg @@ -2,3 +2,10 @@ type: Threshold doc: Energy bi-tc HF interface: ezfio + +[converged_tcscf] +type: logical +doc: If |true|, tc-scf has converged +interface: ezfio,provider,ocaml +default: False + From b0fdf35983c6458805923c1a440cff72ad7fc150 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Fri, 1 Mar 2024 13:37:46 +0100 Subject: [PATCH 5/5] PROVIDING bug --- plugins/local/tc_bi_ortho/slater_tc_opt.irp.f | 23 +- .../tc_bi_ortho/slater_tc_opt_diag.irp.f | 2 + .../tc_bi_ortho/slater_tc_opt_single.irp.f | 419 ++++++++++-------- plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f | 13 + plugins/local/tc_bi_ortho/tc_hmat.irp.f | 28 +- 5 files changed, 281 insertions(+), 204 deletions(-) diff --git a/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f b/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f index 9901a853..59efc943 100644 --- a/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f +++ b/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f @@ -8,8 +8,13 @@ subroutine provide_all_three_ints_bi_ortho() END_DOC implicit none + double precision :: t1, t2 + PROVIDE ao_two_e_integrals_in_map + print *, ' start provide_all_three_ints_bi_ortho' + call wall_time(t1) + if(three_body_h_tc) then if(three_e_3_idx_term) then @@ -32,6 +37,9 @@ subroutine provide_all_three_ints_bi_ortho() endif + call wall_time(t2) + print *, ' end provide_all_three_ints_bi_ortho after (min) = ', (t2-t1)/60.d0 + return end @@ -83,8 +91,11 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) double precision, intent(out) :: hmono, htwoe, hthree, htot + integer :: degree + PROVIDE pure_three_body_h_tc + hmono = 0.d0 htwoe = 0.d0 htot = 0.d0 @@ -99,7 +110,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, if(degree == 0) then call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot) else if (degree == 1) then - call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot) + call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) else if(degree == 2) then call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) endif @@ -111,7 +122,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, if(degree == 0) then call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot) else if (degree == 1) then - call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot) + call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) else if(degree == 2) then call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) else @@ -149,16 +160,16 @@ subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot) double precision, intent(out) :: htot integer :: degree - htot = 0.d0 + htot = 0.d0 call get_excitation_degree(key_i, key_j, degree, Nint) if(degree.gt.2) return - if(degree == 0)then + if(degree == 0) then call diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_i,htot) - else if (degree == 1)then + else if (degree == 1) then call single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint,key_j, key_i , htot) - else if(degree == 2)then + else if(degree == 2) then call double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) endif diff --git a/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f b/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f index cc1a0603..78f9dc66 100644 --- a/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f +++ b/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f @@ -15,6 +15,8 @@ implicit none double precision :: hmono, htwoe, htot, hthree + PROVIDE N_int + PROVIDE HF_bitmask PROVIDE mo_l_coef mo_r_coef call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot) diff --git a/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f b/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f index 81bf69f4..e57cb05c 100644 --- a/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f +++ b/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f @@ -19,6 +19,7 @@ subroutine single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) double precision, intent(out) :: hmono, htwoe, hthree, htot + integer :: occ(Nint*bit_kind_size,2) integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk integer :: degree,exc(0:2,2,2) @@ -44,27 +45,28 @@ subroutine single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, call bitstring_to_list_ab(key_i, occ, Ne, Nint) call get_single_excitation(key_i, key_j, exc, phase, Nint) call decode_exc(exc, 1, h1, p1, h2, p2, s1, s2) - call get_single_excitation_from_fock_tc(key_i, key_j, h1, p1, s1, phase, hmono, htwoe, hthree, htot) + call get_single_excitation_from_fock_tc(Nint, key_i, key_j, h1, p1, s1, phase, hmono, htwoe, hthree, htot) end ! --- -subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, hmono, htwoe, hthree, htot) +subroutine get_single_excitation_from_fock_tc(Nint, key_i, key_j, h, p, spin, phase, hmono, htwoe, hthree, htot) use bitmasks implicit none + integer, intent(in) :: Nint integer, intent(in) :: h, p, spin double precision, intent(in) :: phase - integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2) + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) double precision, intent(out) :: hmono, htwoe, hthree, htot - integer(bit_kind) :: differences(N_int,2) - integer(bit_kind) :: hole(N_int,2) - integer(bit_kind) :: partcl(N_int,2) - integer :: occ_hole(N_int*bit_kind_size,2) - integer :: occ_partcl(N_int*bit_kind_size,2) + integer(bit_kind) :: differences(Nint,2) + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: partcl(Nint,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer :: occ_partcl(Nint*bit_kind_size,2) integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) integer :: i0,i double precision :: buffer_c(mo_num),buffer_x(mo_num) @@ -74,7 +76,7 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h) enddo - do i = 1, N_int + do i = 1, Nint differences(i,1) = xor(key_i(i,1), ref_closed_shell_bitmask(i,1)) differences(i,2) = xor(key_i(i,2), ref_closed_shell_bitmask(i,2)) hole (i,1) = iand(differences(i,1), ref_closed_shell_bitmask(i,1)) @@ -83,8 +85,8 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h partcl (i,2) = iand(differences(i,2), key_i(i,2)) enddo - call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) - call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) + call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint) + call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint) hmono = mo_bi_ortho_tc_one_e(p,h) htwoe = fock_op_2_e_tc_closed_shell(p,h) @@ -122,7 +124,7 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h hthree = 0.d0 if (three_body_h_tc .and. elec_num.gt.2 .and. three_e_4_idx_term) then - call three_comp_fock_elem(key_i, h, p, spin, hthree) + call three_comp_fock_elem(Nint, key_i, h, p, spin, hthree) endif htwoe = htwoe * phase @@ -134,24 +136,27 @@ end ! --- -subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree) - implicit none - integer,intent(in) :: h_fock,p_fock,ispin_fock - integer(bit_kind), intent(in) :: key_i(N_int,2) - double precision, intent(out) :: hthree - integer :: nexc(2),i,ispin,na,nb - integer(bit_kind) :: hole(N_int,2) - integer(bit_kind) :: particle(N_int,2) - integer :: occ_hole(N_int*bit_kind_size,2) - integer :: occ_particle(N_int*bit_kind_size,2) - integer :: n_occ_ab_hole(2),n_occ_ab_particle(2) - integer(bit_kind) :: det_tmp(N_int,2) +subroutine three_comp_fock_elem(Nint, key_i, h_fock, p_fock, ispin_fock, hthree) + implicit none + integer, intent(in) :: Nint + integer, intent(in) :: h_fock, p_fock, ispin_fock + integer(bit_kind), intent(in) :: key_i(Nint,2) + double precision, intent(out) :: hthree + + integer :: nexc(2),i,ispin,na,nb + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: particle(Nint,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer :: occ_particle(Nint*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_particle(2) + integer(bit_kind) :: det_tmp(Nint,2) nexc(1) = 0 nexc(2) = 0 + !! Get all the holes and particles of key_i with respect to the ROHF determinant - do i=1,N_int + do i = 1, Nint hole(i,1) = xor(key_i(i,1),ref_bitmask(i,1)) hole(i,2) = xor(key_i(i,2),ref_bitmask(i,2)) particle(i,1) = iand(hole(i,1),key_i(i,1)) @@ -161,13 +166,14 @@ subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree) nexc(1) = nexc(1) + popcnt(hole(i,1)) nexc(2) = nexc(2) + popcnt(hole(i,2)) enddo + integer :: tmp(2) !DIR$ FORCEINLINE - call bitstring_to_list_ab(particle, occ_particle, tmp, N_int) + call bitstring_to_list_ab(particle, occ_particle, tmp, Nint) ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha ASSERT (tmp(2) == nexc(2)) ! Number of particle beta !DIR$ FORCEINLINE - call bitstring_to_list_ab(hole, occ_hole, tmp, N_int) + call bitstring_to_list_ab(hole, occ_hole, tmp, Nint) ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha ASSERT (tmp(2) == nexc(2)) ! Number of holes beta @@ -181,15 +187,18 @@ subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree) do ispin=1,2 na = elec_num_tab(ispin) nb = elec_num_tab(iand(ispin,1)+1) - do i=1,nexc(ispin) + do i = 1, nexc(ispin) !DIR$ FORCEINLINE - call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb) + call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, Nint, na, nb) !DIR$ FORCEINLINE - call fock_a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb) + call fock_a_tc_operator ( occ_hole (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, Nint, na, nb) enddo enddo + end +! --- + subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb) use bitmasks implicit none @@ -365,111 +374,118 @@ subroutine fock_a_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,N end +! --- -BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) ] - implicit none - BEGIN_DOC -! Closed-shell part of the Fock operator for the TC operator - END_DOC - integer :: h0,p0,h,p,k0,k,i - integer :: n_occ_ab(2) - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab_virt(2) - integer :: occ_virt(N_int*bit_kind_size,2) - integer(bit_kind) :: key_test(N_int) - integer(bit_kind) :: key_virt(N_int,2) - double precision :: accu +BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num)] - fock_op_2_e_tc_closed_shell = -1000.d0 - call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) - do i = 1, N_int - key_virt(i,1) = full_ijkl_bitmask(i) - key_virt(i,2) = full_ijkl_bitmask(i) - key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1)) - key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2)) - enddo - call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) - ! docc ---> virt single excitations - do h0 = 1, n_occ_ab(1) - h=occ(h0,1) - do p0 = 1, n_occ_ab_virt(1) - p = occ_virt(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) - enddo - fock_op_2_e_tc_closed_shell(p,h) = accu + BEGIN_DOC + ! Closed-shell part of the Fock operator for the TC operator + END_DOC + + implicit none + + PROVIDE N_int + + integer :: h0,p0,h,p,k0,k,i + integer :: n_occ_ab(2) + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab_virt(2) + integer :: occ_virt(N_int*bit_kind_size,2) + integer(bit_kind) :: key_test(N_int) + integer(bit_kind) :: key_virt(N_int,2) + double precision :: accu + + fock_op_2_e_tc_closed_shell = -1000.d0 + call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) + + do i = 1, N_int + key_virt(i,1) = full_ijkl_bitmask(i) + key_virt(i,2) = full_ijkl_bitmask(i) + key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1)) + key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2)) enddo - enddo - - do h0 = 1, n_occ_ab_virt(1) - h = occ_virt(h0,1) - do p0 = 1, n_occ_ab(1) - p=occ(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) - enddo - fock_op_2_e_tc_closed_shell(p,h) = accu + call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) + ! docc ---> virt single excitations + do h0 = 1, n_occ_ab(1) + h = occ(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p = occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo enddo - enddo - - ! virt ---> virt single excitations - do h0 = 1, n_occ_ab_virt(1) - h=occ_virt(h0,1) - do p0 = 1, n_occ_ab_virt(1) - p = occ_virt(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) - enddo - fock_op_2_e_tc_closed_shell(p,h) = accu + + do h0 = 1, n_occ_ab_virt(1) + h = occ_virt(h0,1) + do p0 = 1, n_occ_ab(1) + p = occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo enddo - enddo - - do h0 = 1, n_occ_ab_virt(1) - h = occ_virt(h0,1) - do p0 = 1, n_occ_ab_virt(1) - p=occ_virt(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + + ! virt ---> virt single excitations + do h0 = 1, n_occ_ab_virt(1) + h=occ_virt(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p = occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu enddo - fock_op_2_e_tc_closed_shell(p,h) = accu enddo - enddo - - - ! docc ---> docc single excitations - do h0 = 1, n_occ_ab(1) - h=occ(h0,1) - do p0 = 1, n_occ_ab(1) - p = occ(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + + do h0 = 1, n_occ_ab_virt(1) + h = occ_virt(h0,1) + do p0 = 1, n_occ_ab_virt(1) + p=occ_virt(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu enddo - fock_op_2_e_tc_closed_shell(p,h) = accu enddo - enddo - - do h0 = 1, n_occ_ab(1) - h = occ(h0,1) - do p0 = 1, n_occ_ab(1) - p=occ(p0,1) - accu = 0.d0 - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + + + ! docc ---> docc single excitations + do h0 = 1, n_occ_ab(1) + h=occ(h0,1) + do p0 = 1, n_occ_ab(1) + p = occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu + enddo + enddo + + do h0 = 1, n_occ_ab(1) + h = occ(h0,1) + do p0 = 1, n_occ_ab(1) + p=occ(p0,1) + accu = 0.d0 + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) + enddo + fock_op_2_e_tc_closed_shell(p,h) = accu enddo - fock_op_2_e_tc_closed_shell(p,h) = accu enddo - enddo ! do i = 1, mo_num ! write(*,'(100(F10.5,X))')fock_op_2_e_tc_closed_shell(:,i) @@ -477,8 +493,10 @@ BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) END_PROVIDER +! --- subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) + BEGIN_DOC ! for single excitation ONLY FOR ONE- AND TWO-BODY TERMS !! @@ -492,8 +510,9 @@ subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) implicit none integer, intent(in) :: Nint integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) - double precision, intent(out) :: htot - double precision :: hmono, htwoe + double precision, intent(out) :: htot + + double precision :: hmono, htwoe integer :: occ(Nint*bit_kind_size,2) integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk integer :: degree,exc(0:2,2,2) @@ -517,75 +536,85 @@ subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) call get_single_excitation(key_i, key_j, exc, phase, Nint) call decode_exc(exc,1,h1,p1,h2,p2,s1,s2) - call get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h1,p1,s1,phase,hmono,htwoe,htot) -end - - -subroutine get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h,p,spin,phase,hmono,htwoe,htot) - use bitmasks - implicit none - integer,intent(in) :: h,p,spin - double precision, intent(in) :: phase - integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2) - double precision, intent(out) :: hmono,htwoe,htot - integer(bit_kind) :: differences(N_int,2) - integer(bit_kind) :: hole(N_int,2) - integer(bit_kind) :: partcl(N_int,2) - integer :: occ_hole(N_int*bit_kind_size,2) - integer :: occ_partcl(N_int*bit_kind_size,2) - integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) - integer :: i0,i - double precision :: buffer_c(mo_num),buffer_x(mo_num) - do i=1, mo_num - buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h) - buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h) - enddo - do i = 1, N_int - differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1)) - differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2)) - hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1)) - hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) - partcl(i,1) = iand(differences(i,1),key_i(i,1)) - partcl(i,2) = iand(differences(i,2),key_i(i,2)) - enddo - call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) - call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) - hmono = mo_bi_ortho_tc_one_e(p,h) - htwoe = fock_op_2_e_tc_closed_shell(p,h) - ! holes :: direct terms - do i0 = 1, n_occ_ab_hole(1) - i = occ_hole(i0,1) - htwoe -= buffer_c(i) - enddo - do i0 = 1, n_occ_ab_hole(2) - i = occ_hole(i0,2) - htwoe -= buffer_c(i) - enddo - - ! holes :: exchange terms - do i0 = 1, n_occ_ab_hole(spin) - i = occ_hole(i0,spin) - htwoe += buffer_x(i) - enddo - - ! particles :: direct terms - do i0 = 1, n_occ_ab_partcl(1) - i = occ_partcl(i0,1) - htwoe += buffer_c(i) - enddo - do i0 = 1, n_occ_ab_partcl(2) - i = occ_partcl(i0,2) - htwoe += buffer_c(i) - enddo - - ! particles :: exchange terms - do i0 = 1, n_occ_ab_partcl(spin) - i = occ_partcl(i0,spin) - htwoe -= buffer_x(i) - enddo - htwoe = htwoe * phase - hmono = hmono * phase - htot = htwoe + hmono + call get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h1, p1, s1, phase, hmono, htwoe, htot) + +end + +! --- + +subroutine get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h, p, spin, phase, hmono, htwoe, htot) + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer, intent(in) :: h, p, spin + double precision, intent(in) :: phase + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hmono,htwoe,htot + + integer(bit_kind) :: differences(Nint,2) + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: partcl(Nint,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer :: occ_partcl(Nint*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) + integer :: i0,i + double precision :: buffer_c(mo_num), buffer_x(mo_num) + + do i = 1, mo_num + buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h) + buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h) + enddo + + do i = 1, Nint + differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1)) + differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2)) + hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1)) + hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) + partcl(i,1) = iand(differences(i,1),key_i(i,1)) + partcl(i,2) = iand(differences(i,2),key_i(i,2)) + enddo + + call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint) + call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint) + hmono = mo_bi_ortho_tc_one_e(p,h) + htwoe = fock_op_2_e_tc_closed_shell(p,h) + + ! holes :: direct terms + do i0 = 1, n_occ_ab_hole(1) + i = occ_hole(i0,1) + htwoe -= buffer_c(i) + enddo + do i0 = 1, n_occ_ab_hole(2) + i = occ_hole(i0,2) + htwoe -= buffer_c(i) + enddo + + ! holes :: exchange terms + do i0 = 1, n_occ_ab_hole(spin) + i = occ_hole(i0,spin) + htwoe += buffer_x(i) + enddo + + ! particles :: direct terms + do i0 = 1, n_occ_ab_partcl(1) + i = occ_partcl(i0,1) + htwoe += buffer_c(i) + enddo + do i0 = 1, n_occ_ab_partcl(2) + i = occ_partcl(i0,2) + htwoe += buffer_c(i) + enddo + + ! particles :: exchange terms + do i0 = 1, n_occ_ab_partcl(spin) + i = occ_partcl(i0,spin) + htwoe -= buffer_x(i) + enddo + htwoe = htwoe * phase + hmono = hmono * phase + htot = htwoe + hmono end diff --git a/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f b/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f index 64982ab6..398e96db 100644 --- a/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f +++ b/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f @@ -7,6 +7,10 @@ program tc_bi_ortho ! END_DOC + implicit none + + PROVIDE N_int + my_grid_becke = .True. PROVIDE tc_grid1_a tc_grid1_r my_n_pt_r_grid = tc_grid1_r @@ -66,6 +70,15 @@ subroutine routine_diag() ! provide overlap_bi_ortho ! provide htilde_matrix_elmt_bi_ortho + if(noL_standard) then + PROVIDE noL_0e + PROVIDE noL_1e + PROVIDE noL_2e + endif + + PROVIDE htilde_matrix_elmt_bi_ortho + return + if(N_states .eq. 1) then print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1) diff --git a/plugins/local/tc_bi_ortho/tc_hmat.irp.f b/plugins/local/tc_bi_ortho/tc_hmat.irp.f index 88652caa..abec410d 100644 --- a/plugins/local/tc_bi_ortho/tc_hmat.irp.f +++ b/plugins/local/tc_bi_ortho/tc_hmat.irp.f @@ -13,16 +13,34 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)] implicit none integer :: i, j + double precision :: t1, t2 double precision :: htot + + PROVIDE N_int + PROVIDE psi_det + PROVIDE three_e_3_idx_term - call provide_all_three_ints_bi_ortho + if(noL_standard) then + PROVIDE noL_0e + PROVIDE noL_1e + PROVIDE noL_2e + endif + + print *, ' PROVIDING htilde_matrix_elmt_bi_ortho ...' + call wall_time(t1) + + call provide_all_three_ints_bi_ortho() i = 1 j = 1 call htilde_mu_mat_opt_bi_ortho_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot) - !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j, htot) & + + !$OMP PARALLEL & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, j, htot) & !$OMP SHARED (N_det, psi_det, N_int, htilde_matrix_elmt_bi_ortho) + !$OMP DO do i = 1, N_det do j = 1, N_det ! < J |Htilde | I > @@ -31,7 +49,11 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)] htilde_matrix_elmt_bi_ortho(j,i) = htot enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(t2) + print *, ' wall time for htilde_matrix_elmt_bi_ortho (min) =', (t2-t1)/60.d0 END_PROVIDER