diff --git a/src/bi_ort_ints/three_body_ijmk.irp.f b/src/bi_ort_ints/three_body_ijmk.irp.f index 0c02e4c5..669861b7 100644 --- a/src/bi_ort_ints/three_body_ijmk.irp.f +++ b/src/bi_ort_ints/three_body_ijmk.irp.f @@ -86,6 +86,13 @@ tmp_loc_1 = mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) tmp_loc_2 = tmp_aux_2(ipoint,n) + tmp1(ipoint,1,n) = int2_grad1_u12_bimo_t(ipoint,1,n,n) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,1,k,i) * tmp_loc_2 + tmp1(ipoint,2,n) = int2_grad1_u12_bimo_t(ipoint,2,n,n) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,2,k,i) * tmp_loc_2 + tmp1(ipoint,3,n) = int2_grad1_u12_bimo_t(ipoint,3,n,n) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,3,k,i) * tmp_loc_2 + tmp1(ipoint,4,n) = int2_grad1_u12_bimo_t(ipoint,1,n,n) * int2_grad1_u12_bimo_t(ipoint,1,k,i) & + + int2_grad1_u12_bimo_t(ipoint,2,n,n) * int2_grad1_u12_bimo_t(ipoint,2,k,i) & + + int2_grad1_u12_bimo_t(ipoint,3,n,n) * int2_grad1_u12_bimo_t(ipoint,3,k,i) + enddo enddo !$OMP END DO diff --git a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f index 1d1b26cc..8524253a 100644 --- a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f +++ b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f @@ -90,3 +90,96 @@ subroutine htcdag_bi_ortho_calc_tdav_slow(v, u, N_st, sze) end +subroutine i_H_tc_psi_phi(key,keys,coef_l,coef_r,Nint,Ndet,Ndet_max,Nstate,i_H_chi_array,i_H_phi_array) + use bitmasks + implicit none + BEGIN_DOC +! Computes $\langle i|H|Phi \rangle = \sum_J c^R_J \langle i | H | J \rangle$. +! +! AND $\langle Chi|H| i \rangle = \sum_J c^L_J \langle J | H | i \rangle$. +! +! CONVENTION: i_H_phi_array(0) = total matrix element, +! +! i_H_phi_array(1) = one-electron matrix element, +! +! i_H_phi_array(2) = two-electron matrix element, +! +! i_H_phi_array(3) = three-electron matrix element, +! +! Uses filter_connected_i_H_psi0 to get all the $|J \rangle$ to which $|i \rangle$ +! is connected. +! +! The i_H_psi_minilist is much faster but requires to build the +! minilists. + END_DOC + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef_l(Ndet_max,Nstate),coef_r(Ndet_max,Nstate) + double precision, intent(out) :: i_H_chi_array(0:3,Nstate),i_H_phi_array(0:3,Nstate) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hmono, htwoe, hthree, htot + integer, allocatable :: idx(:) + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + allocate(idx(0:Ndet)) + + i_H_chi_array = 0.d0 + i_H_phi_array = 0.d0 + + call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) + if (Nstate == 1) then + + do ii=1,idx(0) + i = idx(ii) + ! computes + !DIR$ FORCEINLINE + call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot) + i_H_chi_array(0,1) = i_H_chi_array(0,1) + coef_l(i,1)*htot + i_H_chi_array(1,1) = i_H_chi_array(1,1) + coef_l(i,1)*hmono + i_H_chi_array(2,1) = i_H_chi_array(2,1) + coef_l(i,1)*htwoe + i_H_chi_array(3,1) = i_H_chi_array(3,1) + coef_l(i,1)*hthree + ! computes + !DIR$ FORCEINLINE + call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot) + i_H_phi_array(0,1) = i_H_phi_array(0,1) + coef_r(i,1)*htot + i_H_phi_array(1,1) = i_H_phi_array(1,1) + coef_r(i,1)*hmono + i_H_phi_array(2,1) = i_H_phi_array(2,1) + coef_r(i,1)*htwoe + i_H_phi_array(3,1) = i_H_phi_array(3,1) + coef_r(i,1)*hthree + enddo + + else + + do ii=1,idx(0) + i = idx(ii) + ! computes + !DIR$ FORCEINLINE + call htilde_mu_mat_opt_bi_ortho(keys(1,1,i), key, Nint, hmono, htwoe, hthree, htot) + do j = 1, Nstate + i_H_chi_array(0,j) = i_H_chi_array(0,j) + coef_l(i,j)*htot + i_H_chi_array(1,j) = i_H_chi_array(1,j) + coef_l(i,j)*hmono + i_H_chi_array(2,j) = i_H_chi_array(2,j) + coef_l(i,j)*htwoe + i_H_chi_array(3,j) = i_H_chi_array(3,j) + coef_l(i,j)*hthree + enddo + ! computes + !DIR$ FORCEINLINE + call htilde_mu_mat_opt_bi_ortho(key,keys(1,1,i), Nint, hmono, htwoe, hthree, htot) + do j = 1, Nstate + i_H_phi_array(0,j) = i_H_phi_array(0,j) + coef_r(i,j)*htot + i_H_phi_array(1,j) = i_H_phi_array(1,j) + coef_r(i,j)*hmono + i_H_phi_array(2,j) = i_H_phi_array(2,j) + coef_r(i,j)*htwoe + i_H_phi_array(3,j) = i_H_phi_array(3,j) + coef_r(i,j)*hthree + enddo + enddo + + endif + +end + diff --git a/src/tc_bi_ortho/slater_tc_3e_slow.irp.f b/src/tc_bi_ortho/slater_tc_3e_slow.irp.f index 35abbbc4..cb33d343 100644 --- a/src/tc_bi_ortho/slater_tc_3e_slow.irp.f +++ b/src/tc_bi_ortho/slater_tc_3e_slow.irp.f @@ -184,7 +184,7 @@ subroutine single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree) ii = occ(i,s1) do j = i+1, Ne(s1) jj = occ(j,s1) -! ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1) +! !ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1) hthree += three_e_single_parrallel_spin(jj,ii,p1,h1) ! USES THE 4-IDX TENSOR enddo enddo diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f index 182c03d7..f1a7cc0a 100644 --- a/src/tc_bi_ortho/test_tc_fock.irp.f +++ b/src/tc_bi_ortho/test_tc_fock.irp.f @@ -152,9 +152,7 @@ subroutine routine_tot() ! do i = 1, elec_num_tab(s1) ! do a = elec_num_tab(s1)+1, mo_num ! virtual do i = 1, elec_beta_num - do a = elec_beta_num+1, elec_alpha_num! virtual -! do i = elec_beta_num+1, elec_alpha_num -! do a = elec_alpha_num+1, mo_num! virtual + do a = elec_beta_num+1, mo_num! virtual print*,i,a det_i = ref_bitmask @@ -167,7 +165,7 @@ subroutine routine_tot() call htilde_mu_mat_bi_ortho_slow(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij) print*,htilde_ij - if(dabs(htilde_ij).lt.1.d-10)cycle +! if(dabs(htilde_ij).lt.1.d-10)cycle print*, ' excited det' call debug_det(det_i, N_int) @@ -184,9 +182,12 @@ subroutine routine_tot() ! endif err_ai = dabs(dabs(ref) - dabs(new)) if(err_ai .gt. 1d-7) then + print*,'---------' print*,'s1 = ',s1 print*, ' warning on', i, a print*, ref,new,err_ai + print*,hmono, htwoe, hthree + print*,'---------' endif print*, ref,new,err_ai err_tot += err_ai diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index 0ae515bb..f4553f3e 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -208,10 +208,10 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ] if(three_body_h_tc) then !call wall_time(tt0) - !PROVIDE fock_a_tot_3e_bi_orth - !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth - PROVIDE fock_3e_uhf_mo_a - Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a + PROVIDE fock_a_tot_3e_bi_orth + Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth +! PROVIDE fock_3e_uhf_mo_a +! Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a !call wall_time(tt1) !print*, ' 3-e term:', tt1-tt0 endif @@ -241,21 +241,13 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ] if(bi_ortho) then - !allocate(tmp(ao_num,ao_num)) - !tmp = Fock_matrix_tc_ao_beta - !if(three_body_h_tc) then - ! tmp += fock_3e_uhf_ao_b - !endif - !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1)) - !deallocate(tmp) - call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) & , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) ) if(three_body_h_tc) then - !PROVIDE fock_b_tot_3e_bi_orth - !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth - PROVIDE fock_3e_uhf_mo_b - Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b + PROVIDE fock_b_tot_3e_bi_orth + Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth +! PROVIDE fock_3e_uhf_mo_b +! Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b endif else