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