diff --git a/src/bi_ort_ints/semi_num_ints_mo.irp.f b/src/bi_ort_ints/semi_num_ints_mo.irp.f index e32d4707..33f512cf 100644 --- a/src/bi_ort_ints/semi_num_ints_mo.irp.f +++ b/src/bi_ort_ints/semi_num_ints_mo.irp.f @@ -127,6 +127,9 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, implicit none integer :: ipoint + print*,'providing int2_grad1_u12_bimo_transp' + double precision :: wall0, wall1 + call wall_time(wall0) !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint) & @@ -142,6 +145,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, enddo !$OMP END DO !$OMP END PARALLEL + call wall_time(wall1) + print*,'Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0 END_PROVIDER diff --git a/src/tc_bi_ortho/slater_tc_3e.irp.f b/src/tc_bi_ortho/slater_tc_3e.irp.f index 0d5f8542..4bfb2da3 100644 --- a/src/tc_bi_ortho/slater_tc_3e.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e.irp.f @@ -49,6 +49,8 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) if(Ne(1)+Ne(2).ge.3)then !! ! alpha/alpha/beta three-body + double precision :: accu + accu = 0.d0 do i = 1, Ne(1) ii = occ(i,1) do j = i+1, Ne(1) @@ -60,11 +62,14 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR hthree += direct_int - exchange_int + accu += direct_int - exchange_int enddo enddo enddo + print*,'aab = ',accu ! beta/beta/alpha three-body + accu = 0.d0 do i = 1, Ne(2) ii = occ(i,2) do j = i+1, Ne(2) @@ -74,11 +79,14 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) hthree += direct_int - exchange_int + accu += direct_int - exchange_int enddo enddo enddo + print*,'abb = ',accu ! alpha/alpha/alpha three-body + accu = 0.d0 do i = 1, Ne(1) ii = occ(i,1) ! 1 do j = i+1, Ne(1) @@ -87,11 +95,14 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) mm = occ(m,1) ! 3 ! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS + accu += three_e_diag_parrallel_spin(mm,jj,ii) enddo enddo enddo + print*,'aaa = ',accu ! beta/beta/beta three-body + accu = 0.d0 do i = 1, Ne(2) ii = occ(i,2) ! 1 do j = i+1, Ne(2) @@ -100,9 +111,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) mm = occ(m,2) ! 3 ! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS + accu += three_e_diag_parrallel_spin(mm,jj,ii) enddo enddo enddo + print*,'bbb = ',accu endif end diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f index d585ce6c..a49a5958 100644 --- a/src/tc_bi_ortho/test_tc_fock.irp.f +++ b/src/tc_bi_ortho/test_tc_fock.irp.f @@ -13,105 +13,44 @@ program test_tc_fock !call routine_1 !call routine_2 - call routine_3() +! call routine_3() + call test_3e end ! --- -subroutine routine_0 +subroutine test_3e implicit none - use bitmasks ! you need to include the bitmasks_module.f90 features - integer :: i,a,j,m,i_ok - integer :: exc(0:2,2,2),h1,p1,s1,h2,p2,s2,degree - - integer(bit_kind), allocatable :: det_i(:,:) - double precision :: hmono,htwoe,hthree,htilde_ij,phase - double precision :: same, op, tot, accu - allocate(det_i(N_int,2)) - s1 = 1 - accu = 0.d0 - do i = 1, elec_alpha_num ! occupied - do a = elec_alpha_num+1, mo_num ! virtual - det_i = ref_bitmask - call do_single_excitation(det_i,i,a,s1,i_ok) - if(i_ok == -1)then - print*,'PB !!' - print*,i,a - stop - endif -! call debug_det(det_i,N_int) - call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) - call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) - op = fock_3_mat_a_op_sh_bi_orth(a,i) - same = fock_3_mat_a_sa_sh_bi_orth(a,i) -! same = 0.d0 - tot = same + op - if(dabs(tot - phase*hthree).gt.1.d-10)then - print*,'------' - print*,i,a,phase - print*,'hthree = ',phase*hthree - print*,'fock = ',tot - print*,'same,op= ',same,op - print*,dabs(tot - phase*hthree) - stop - endif - accu += dabs(tot - phase*hthree) - enddo - enddo + double precision :: integral_aaa,integral_aab,integral_abb,integral_bbb,accu + double precision :: hmono, htwoe, hthree, htot + call htilde_mu_mat_bi_ortho(ref_bitmask, ref_bitmask, N_int, hmono, htwoe, hthree, htot) +! call diag_htilde_three_body_ints_bi_ort(N_int, ref_bitmask, hthree) + print*,'hmono = ',hmono + print*,'htwoe = ',htwoe + print*,'hthree= ',hthree + print*,'htot = ',htot + print*,'' + print*,'' + print*,'TC_one= ',TC_HF_one_electron_energy + print*,'TC_two= ',TC_HF_two_e_energy + print*,'TC_3e = ',diag_three_elem_hf + print*,'TC_tot= ',TC_HF_energy + print*,'' + print*,'' + call give_aaa_contrib(integral_aaa) + print*,'integral_aaa = ',integral_aaa + call give_aab_contrib(integral_aab) + print*,'integral_aab = ',integral_aab + call give_abb_contrib(integral_abb) + print*,'integral_abb = ',integral_abb + call give_bbb_contrib(integral_bbb) + print*,'integral_bbb = ',integral_bbb + accu = integral_aaa + integral_aab + integral_abb + integral_bbb print*,'accu = ',accu + print*,'delta = ',hthree - accu -end subroutine routine_0 - -! --- - -subroutine routine_1 - - implicit none - integer :: i, a - double precision :: accu - - accu = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - accu += dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) - !if(dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) .gt. 1.d-10)then - print*, i, a - print*, dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) & - , fock_3_mat_a_op_sh_bi_orth_old(a,i), fock_3_mat_a_op_sh_bi_orth(a,i) - !endif - enddo - enddo - - print *, 'accu = ', accu - -end subroutine routine_1 - -! --- - -subroutine routine_2 - - implicit none - integer :: i, a - double precision :: accu - - accu = 0.d0 - do i = 1, mo_num - do a = 1, mo_num - accu += dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) - !if(dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) .gt. 1.d-10)then - print*, i, a - print*, dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) & - , fock_3_mat_a_sa_sh_bi_orth_old(a,i), fock_3_mat_a_sa_sh_bi_orth(a,i) - !endif - enddo - enddo - - print *, 'accu = ', accu - -end subroutine routine_2 - -! --- +end subroutine routine_3() diff --git a/src/tc_scf/fock_three.irp.f b/src/tc_scf/fock_three.irp.f index e9ad4ce5..f73a5049 100644 --- a/src/tc_scf/fock_three.irp.f +++ b/src/tc_scf/fock_three.irp.f @@ -74,39 +74,46 @@ BEGIN_PROVIDER [double precision, diag_three_elem_hf] implicit none integer :: i,j,k,ipoint,mm double precision :: contrib,weight,four_third,one_third,two_third,exchange_int_231 - if(.not.bi_ortho)then - if(three_body_h_tc)then - one_third = 1.d0/3.d0 - two_third = 2.d0/3.d0 - four_third = 4.d0/3.d0 + print*,'providing diag_three_elem_hf' + if(.not.three_body_h_tc)then diag_three_elem_hf = 0.d0 - do i = 1, elec_beta_num - do j = 1, elec_beta_num - do k = 1, elec_beta_num - call give_integrals_3_body(k,j,i,j,i,k,exchange_int_231) - diag_three_elem_hf += two_third * exchange_int_231 + else + if(.not.bi_ortho)then + one_third = 1.d0/3.d0 + two_third = 2.d0/3.d0 + four_third = 4.d0/3.d0 + diag_three_elem_hf = 0.d0 + do i = 1, elec_beta_num + do j = 1, elec_beta_num + do k = 1, elec_beta_num + call give_integrals_3_body(k,j,i,j,i,k,exchange_int_231) + diag_three_elem_hf += two_third * exchange_int_231 + enddo enddo enddo - enddo - do mm = 1, 3 - do ipoint = 1, n_points_final_grid - weight = final_weight_at_r_vector(ipoint) - contrib = 3.d0 * fock_3_w_kk_sum(ipoint,mm) * fock_3_rho_beta(ipoint) * fock_3_w_kk_sum(ipoint,mm) & - -2.d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * fock_3_w_kk_sum(ipoint,mm) & - -1.d0 * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) - contrib *= four_third - contrib += -two_third * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) & - - four_third * fock_3_w_kk_sum(ipoint,mm) * fock_3_w_kl_mo_k_mo_l(ipoint,mm) - diag_three_elem_hf += weight * contrib + do mm = 1, 3 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + contrib = 3.d0 * fock_3_w_kk_sum(ipoint,mm) * fock_3_rho_beta(ipoint) * fock_3_w_kk_sum(ipoint,mm) & + -2.d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * fock_3_w_kk_sum(ipoint,mm) & + -1.d0 * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) + contrib *= four_third + contrib += -two_third * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) & + - four_third * fock_3_w_kk_sum(ipoint,mm) * fock_3_w_kl_mo_k_mo_l(ipoint,mm) + diag_three_elem_hf += weight * contrib + enddo enddo - enddo - diag_three_elem_hf = - diag_three_elem_hf - else - diag_three_elem_hf = 0.D0 + diag_three_elem_hf = - diag_three_elem_hf + else + double precision :: integral_aaa,hthree, integral_aab,integral_abb,integral_bbb + provide mo_l_coef mo_r_coef + call give_aaa_contrib(integral_aaa) + call give_aab_contrib(integral_aab) + call give_abb_contrib(integral_abb) + call give_bbb_contrib(integral_bbb) + diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb endif - else - diag_three_elem_hf = 0.D0 - endif + endif END_PROVIDER diff --git a/src/tc_scf/fock_three_bi_ortho_new_new.irp.f b/src/tc_scf/fock_three_bi_ortho_new_new.irp.f index b56fc1ff..b0345957 100644 --- a/src/tc_scf/fock_three_bi_ortho_new_new.irp.f +++ b/src/tc_scf/fock_three_bi_ortho_new_new.irp.f @@ -154,6 +154,7 @@ BEGIN_PROVIDER [double precision, fock_b_tmp2_bi_ortho, (mo_num, mo_num)] END_PROVIDER subroutine contrib_3e_sss(a,i,j,k,integral) + implicit none integer, intent(in) :: a,i,j,k BEGIN_DOC ! returns the pure same spin contribution to F(a,i) from two orbitals j,k @@ -173,6 +174,7 @@ subroutine contrib_3e_sss(a,i,j,k,integral) end subroutine contrib_3e_soo(a,i,j,k,integral) + implicit none integer, intent(in) :: a,i,j,k BEGIN_DOC ! returns the same spin / opposite spin / opposite spin contribution to F(a,i) from two orbitals j,k diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index bb382c7c..4a875b59 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -81,8 +81,8 @@ subroutine routine_scf() print*,'TC HF total energy = ', TC_HF_energy print*,'TC HF 1 e energy = ', TC_HF_one_electron_energy print*,'TC HF 2 e energy = ', TC_HF_two_e_energy - if(.not. bi_ortho)then - print*,'TC HF 3 body = ', diag_three_elem_hf + if(three_body_h_tc)then + print*,'TC HF 3 body = ', diag_three_elem_hf endif print*,'***' e_delta = 10.d0 @@ -124,6 +124,9 @@ subroutine routine_scf() print*,'TC HF total energy = ', TC_HF_energy print*,'TC HF 1 e energy = ', TC_HF_one_electron_energy print*,'TC HF 2 non hermit = ', TC_HF_two_e_energy + if(three_body_h_tc)then + print*,'TC HF 3 body = ', diag_three_elem_hf + endif print*,'***' e_delta = dabs( TC_HF_energy - e_save ) print*, 'it, delta E = ', it, e_delta diff --git a/src/tc_scf/three_e_energy_bi_ortho.irp.f b/src/tc_scf/three_e_energy_bi_ortho.irp.f new file mode 100644 index 00000000..64212da8 --- /dev/null +++ b/src/tc_scf/three_e_energy_bi_ortho.irp.f @@ -0,0 +1,174 @@ + +subroutine contrib_3e_diag_sss(i,j,k,integral) + implicit none + integer, intent(in) :: i,j,k + BEGIN_DOC + ! returns the pure same spin contribution to diagonal matrix element of 3e term + END_DOC + double precision, intent(out) :: integral + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + call give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int )!!! < i k j | i k j > + call give_integrals_3_body_bi_ort(i, k, j, j, i, k, c_3_int) ! < i k j | j i k > + call give_integrals_3_body_bi_ort(i, k, j, k, j, i, c_minus_3_int)! < i k j | k j i > + integral = direct_int + c_3_int + c_minus_3_int + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(i, k, j, j, k, i, exch_13_int)!!! < i k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)!!! < i k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(i, k, j, k, i, j, exch_12_int)!!! < i k j | k i j > : E_12 + integral += - exch_13_int - exch_23_int - exch_12_int + integral = -integral +end + +subroutine contrib_3e_diag_soo(i,j,k,integral) + implicit none + integer, intent(in) :: i,j,k + BEGIN_DOC + ! returns the pure same spin contribution to diagonal matrix element of 3e term + END_DOC + double precision, intent(out) :: integral + double precision :: direct_int, exch_23_int + call give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int) ! < i k j | i k j > + call give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)! < i k j | i j k > : E_23 + integral = direct_int - exch_23_int + integral = -integral +end + + +subroutine give_aaa_contrib_bis(integral_aaa) + implicit none + double precision, intent(out) :: integral_aaa + double precision :: integral + integer :: i,j,k + integral_aaa = 0.d0 + do i = 1, elec_alpha_num + do j = i+1, elec_alpha_num + do k = j+1, elec_alpha_num + call contrib_3e_diag_sss(i,j,k,integral) + integral_aaa += integral + enddo + enddo + enddo + +end + +subroutine give_aaa_contrib(integral_aaa) + implicit none + double precision, intent(out) :: integral_aaa + double precision :: integral + integer :: i,j,k + integral_aaa = 0.d0 + do i = 1, elec_alpha_num + do j = 1, elec_alpha_num + do k = 1, elec_alpha_num + call contrib_3e_diag_sss(i,j,k,integral) + integral_aaa += integral + enddo + enddo + enddo + integral_aaa *= 1.d0/6.d0 +end + + +subroutine give_aab_contrib(integral_aab) + implicit none + double precision, intent(out) :: integral_aab + double precision :: integral + integer :: i,j,k + integral_aab = 0.d0 + do i = 1, elec_beta_num + do j = 1, elec_alpha_num + do k = 1, elec_alpha_num + call contrib_3e_diag_soo(i,j,k,integral) + integral_aab += integral + enddo + enddo + enddo + integral_aab *= 0.5d0 +end + + +subroutine give_aab_contrib_bis(integral_aab) + implicit none + double precision, intent(out) :: integral_aab + double precision :: integral + integer :: i,j,k + integral_aab = 0.d0 + do i = 1, elec_beta_num + do j = 1, elec_alpha_num + do k = j+1, elec_alpha_num + call contrib_3e_diag_soo(i,j,k,integral) + integral_aab += integral + enddo + enddo + enddo +end + + +subroutine give_abb_contrib(integral_abb) + implicit none + double precision, intent(out) :: integral_abb + double precision :: integral + integer :: i,j,k + integral_abb = 0.d0 + do i = 1, elec_alpha_num + do j = 1, elec_beta_num + do k = 1, elec_beta_num + call contrib_3e_diag_soo(i,j,k,integral) + integral_abb += integral + enddo + enddo + enddo + integral_abb *= 0.5d0 +end + +subroutine give_abb_contrib_bis(integral_abb) + implicit none + double precision, intent(out) :: integral_abb + double precision :: integral + integer :: i,j,k + integral_abb = 0.d0 + do i = 1, elec_alpha_num + do j = 1, elec_beta_num + do k = j+1, elec_beta_num + call contrib_3e_diag_soo(i,j,k,integral) + integral_abb += integral + enddo + enddo + enddo +end + +subroutine give_bbb_contrib_bis(integral_bbb) + implicit none + double precision, intent(out) :: integral_bbb + double precision :: integral + integer :: i,j,k + integral_bbb = 0.d0 + do i = 1, elec_beta_num + do j = i+1, elec_beta_num + do k = j+1, elec_beta_num + call contrib_3e_diag_sss(i,j,k,integral) + integral_bbb += integral + enddo + enddo + enddo + +end + +subroutine give_bbb_contrib(integral_bbb) + implicit none + double precision, intent(out) :: integral_bbb + double precision :: integral + integer :: i,j,k + integral_bbb = 0.d0 + do i = 1, elec_beta_num + do j = 1, elec_beta_num + do k = 1, elec_beta_num + call contrib_3e_diag_sss(i,j,k,integral) + integral_bbb += integral + enddo + enddo + enddo + integral_bbb *= 1.d0/6.d0 +end + +