LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-28 21:34:01 +02:00
Code Releases Activity

PROVIDING bug

Browse Source
This commit is contained in:
Abdallah Ammar 2024-03-01 13:37:46 +01:00
parent 5226bb2ef4
commit b0fdf35983
5 changed files with 281 additions and 204 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
plugins/local/tc_bi_ortho/slater_tc_opt.irp.f
View File

@ -8,8 +8,13 @@ subroutine provide_all_three_ints_bi_ortho()
END_DOC END_DOC
implicit none implicit none
double precision :: t1, t2
PROVIDE ao_two_e_integrals_in_map 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_body_h_tc) then
if(three_e_3_idx_term) then if(three_e_3_idx_term) then
@ -32,6 +37,9 @@ subroutine provide_all_three_ints_bi_ortho()
endif endif
call wall_time(t2)
print *, ' end provide_all_three_ints_bi_ortho after (min) = ', (t2-t1)/60.d0
return return
end 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, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot double precision, intent(out) :: hmono, htwoe, hthree, htot
integer :: degree integer :: degree
PROVIDE pure_three_body_h_tc
hmono = 0.d0 hmono = 0.d0
htwoe = 0.d0 htwoe = 0.d0
htot = 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 if(degree == 0) then
call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot) call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot)
else if (degree == 1) then 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 else if(degree == 2) then
call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
endif endif
@ -111,7 +122,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
if(degree == 0) then if(degree == 0) then
call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot) call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot)
else if (degree == 1) then 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 else if(degree == 2) then
call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot) call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
else 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 double precision, intent(out) :: htot
integer :: degree integer :: degree
htot = 0.d0 htot = 0.d0
call get_excitation_degree(key_i, key_j, degree, Nint) call get_excitation_degree(key_i, key_j, degree, Nint)
if(degree.gt.2) return 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) 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) 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) call double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
endif endif

2
plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f
View File

@ -15,6 +15,8 @@
implicit none implicit none
double precision :: hmono, htwoe, htot, hthree double precision :: hmono, htwoe, htot, hthree
PROVIDE N_int
PROVIDE HF_bitmask
PROVIDE mo_l_coef mo_r_coef PROVIDE mo_l_coef mo_r_coef
call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot) call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot)

419
plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f
View File

@ -19,6 +19,7 @@ subroutine single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
integer, intent(in) :: Nint integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot double precision, intent(out) :: hmono, htwoe, hthree, htot
integer :: occ(Nint*bit_kind_size,2) integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer :: degree,exc(0:2,2,2) 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 bitstring_to_list_ab(key_i, occ, Ne, Nint)
call get_single_excitation(key_i, key_j, exc, phase, Nint) call get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc, 1, h1, p1, h2, p2, s1, s2) 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 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 use bitmasks
implicit none implicit none
integer, intent(in) :: Nint
integer, intent(in) :: h, p, spin integer, intent(in) :: h, p, spin
double precision, intent(in) :: phase 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 double precision, intent(out) :: hmono, htwoe, hthree, htot
integer(bit_kind) :: differences(N_int,2) integer(bit_kind) :: differences(Nint,2)
integer(bit_kind) :: hole(N_int,2) integer(bit_kind) :: hole(Nint,2)
integer(bit_kind) :: partcl(N_int,2) integer(bit_kind) :: partcl(Nint,2)
integer :: occ_hole(N_int*bit_kind_size,2) integer :: occ_hole(Nint*bit_kind_size,2)
integer :: occ_partcl(N_int*bit_kind_size,2) integer :: occ_partcl(Nint*bit_kind_size,2)
integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
integer :: i0,i integer :: i0,i
double precision :: buffer_c(mo_num),buffer_x(mo_num) 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) buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
enddo 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,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)) 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,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)) partcl (i,2) = iand(differences(i,2), key_i(i,2))
enddo enddo
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, 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, N_int) call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint)
hmono = mo_bi_ortho_tc_one_e(p,h) hmono = mo_bi_ortho_tc_one_e(p,h)
htwoe = fock_op_2_e_tc_closed_shell(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 hthree = 0.d0
if (three_body_h_tc .and. elec_num.gt.2 .and. three_e_4_idx_term) then 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 endif
htwoe = htwoe * phase htwoe = htwoe * phase
@ -134,24 +136,27 @@ end
! --- ! ---
subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree) subroutine three_comp_fock_elem(Nint, 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)
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(1) = 0
nexc(2) = 0 nexc(2) = 0
!! Get all the holes and particles of key_i with respect to the ROHF determinant !! 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,1) = xor(key_i(i,1),ref_bitmask(i,1))
hole(i,2) = xor(key_i(i,2),ref_bitmask(i,2)) hole(i,2) = xor(key_i(i,2),ref_bitmask(i,2))
particle(i,1) = iand(hole(i,1),key_i(i,1)) 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(1) = nexc(1) + popcnt(hole(i,1))
nexc(2) = nexc(2) + popcnt(hole(i,2)) nexc(2) = nexc(2) + popcnt(hole(i,2))
enddo enddo
integer :: tmp(2) integer :: tmp(2)
!DIR$ FORCEINLINE !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(1) == nexc(1)) ! Number of particles alpha
ASSERT (tmp(2) == nexc(2)) ! Number of particle beta ASSERT (tmp(2) == nexc(2)) ! Number of particle beta
!DIR$ FORCEINLINE !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(1) == nexc(1)) ! Number of holes alpha
ASSERT (tmp(2) == nexc(2)) ! Number of holes beta 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 do ispin=1,2
na = elec_num_tab(ispin) na = elec_num_tab(ispin)
nb = elec_num_tab(iand(ispin,1)+1) nb = elec_num_tab(iand(ispin,1)+1)
do i=1,nexc(ispin) do i = 1, nexc(ispin)
!DIR$ FORCEINLINE !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 !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
enddo enddo
end end
! ---
subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb) subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb)
use bitmasks use bitmasks
implicit none implicit none
@ -365,111 +374,118 @@ subroutine fock_a_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,N
end end
! ---
BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) ] 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
fock_op_2_e_tc_closed_shell = -1000.d0 BEGIN_DOC
call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) ! Closed-shell part of the Fock operator for the TC operator
do i = 1, N_int END_DOC
key_virt(i,1) = full_ijkl_bitmask(i)
key_virt(i,2) = full_ijkl_bitmask(i) implicit none
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)) PROVIDE N_int
enddo
call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) integer :: h0,p0,h,p,k0,k,i
! docc ---> virt single excitations integer :: n_occ_ab(2)
do h0 = 1, n_occ_ab(1) integer :: occ(N_int*bit_kind_size,2)
h=occ(h0,1) integer :: n_occ_ab_virt(2)
do p0 = 1, n_occ_ab_virt(1) integer :: occ_virt(N_int*bit_kind_size,2)
p = occ_virt(p0,1) integer(bit_kind) :: key_test(N_int)
accu = 0.d0 integer(bit_kind) :: key_virt(N_int,2)
do k0 = 1, n_occ_ab(1) double precision :: accu
k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) fock_op_2_e_tc_closed_shell = -1000.d0
enddo call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int)
fock_op_2_e_tc_closed_shell(p,h) = accu
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
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_virt(1) do h0 = 1, n_occ_ab(1)
h = occ_virt(h0,1) h = occ(h0,1)
do p0 = 1, n_occ_ab(1) do p0 = 1, n_occ_ab_virt(1)
p=occ(p0,1) p = occ_virt(p0,1)
accu = 0.d0 accu = 0.d0
do k0 = 1, n_occ_ab(1) do k0 = 1, n_occ_ab(1)
k = occ(k0,1) k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo enddo
fock_op_2_e_tc_closed_shell(p,h) = accu fock_op_2_e_tc_closed_shell(p,h) = accu
enddo
enddo enddo
enddo
do h0 = 1, n_occ_ab_virt(1)
! virt ---> virt single excitations h = occ_virt(h0,1)
do h0 = 1, n_occ_ab_virt(1) do p0 = 1, n_occ_ab(1)
h=occ_virt(h0,1) p = occ(p0,1)
do p0 = 1, n_occ_ab_virt(1) accu = 0.d0
p = occ_virt(p0,1) do k0 = 1, n_occ_ab(1)
accu = 0.d0 k = occ(k0,1)
do k0 = 1, n_occ_ab(1) accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
k = occ(k0,1) enddo
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) fock_op_2_e_tc_closed_shell(p,h) = accu
enddo 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) do h0 = 1, n_occ_ab_virt(1)
h = occ_virt(h0,1) h=occ_virt(h0,1)
do p0 = 1, n_occ_ab_virt(1) do p0 = 1, n_occ_ab_virt(1)
p=occ_virt(p0,1) p = occ_virt(p0,1)
accu = 0.d0 accu = 0.d0
do k0 = 1, n_occ_ab(1) do k0 = 1, n_occ_ab(1)
k = occ(k0,1) k = occ(k0,1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) 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
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)
! docc ---> docc single excitations do p0 = 1, n_occ_ab_virt(1)
do h0 = 1, n_occ_ab(1) p=occ_virt(p0,1)
h=occ(h0,1) accu = 0.d0
do p0 = 1, n_occ_ab(1) do k0 = 1, n_occ_ab(1)
p = occ(p0,1) k = occ(k0,1)
accu = 0.d0 accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
do k0 = 1, n_occ_ab(1) enddo
k = occ(k0,1) fock_op_2_e_tc_closed_shell(p,h) = accu
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
enddo enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo enddo
enddo
do h0 = 1, n_occ_ab(1) ! docc ---> docc single excitations
h = occ(h0,1) do h0 = 1, n_occ_ab(1)
do p0 = 1, n_occ_ab(1) h=occ(h0,1)
p=occ(p0,1) do p0 = 1, n_occ_ab(1)
accu = 0.d0 p = occ(p0,1)
do k0 = 1, n_occ_ab(1) accu = 0.d0
k = occ(k0,1) do k0 = 1, n_occ_ab(1)
accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h) 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 enddo
fock_op_2_e_tc_closed_shell(p,h) = accu
enddo enddo
enddo
! do i = 1, mo_num ! do i = 1, mo_num
! write(*,'(100(F10.5,X))')fock_op_2_e_tc_closed_shell(:,i) ! 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 END_PROVIDER
! ---
subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot) subroutine single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
BEGIN_DOC BEGIN_DOC
! <key_j |H_tilde | key_i> for single excitation ONLY FOR ONE- AND TWO-BODY TERMS ! <key_j |H_tilde | key_i> 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 implicit none
integer, intent(in) :: Nint integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
double precision, intent(out) :: htot double precision, intent(out) :: htot
double precision :: hmono, htwoe
double precision :: hmono, htwoe
integer :: occ(Nint*bit_kind_size,2) integer :: occ(Nint*bit_kind_size,2)
integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk
integer :: degree,exc(0:2,2,2) 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 get_single_excitation(key_i, key_j, exc, phase, Nint)
call decode_exc(exc,1,h1,p1,h2,p2,s1,s2) 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) call get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h1, p1, s1, phase, hmono, htwoe, htot)
end
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 subroutine get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h, p, spin, phase, hmono, htwoe, htot)
integer,intent(in) :: h,p,spin
double precision, intent(in) :: phase use bitmasks
integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2)
double precision, intent(out) :: hmono,htwoe,htot implicit none
integer(bit_kind) :: differences(N_int,2) integer, intent(in) :: Nint
integer(bit_kind) :: hole(N_int,2) integer, intent(in) :: h, p, spin
integer(bit_kind) :: partcl(N_int,2) double precision, intent(in) :: phase
integer :: occ_hole(N_int*bit_kind_size,2) integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
integer :: occ_partcl(N_int*bit_kind_size,2) double precision, intent(out) :: hmono,htwoe,htot
integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
integer :: i0,i integer(bit_kind) :: differences(Nint,2)
double precision :: buffer_c(mo_num),buffer_x(mo_num) integer(bit_kind) :: hole(Nint,2)
do i=1, mo_num integer(bit_kind) :: partcl(Nint,2)
buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h) integer :: occ_hole(Nint*bit_kind_size,2)
buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h) integer :: occ_partcl(Nint*bit_kind_size,2)
enddo integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
do i = 1, N_int integer :: i0,i
differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1)) double precision :: buffer_c(mo_num), buffer_x(mo_num)
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)) do i = 1, mo_num
hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h)
partcl(i,1) = iand(differences(i,1),key_i(i,1)) buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
partcl(i,2) = iand(differences(i,2),key_i(i,2)) enddo
enddo
call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) do i = 1, Nint
call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1))
hmono = mo_bi_ortho_tc_one_e(p,h) differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2))
htwoe = fock_op_2_e_tc_closed_shell(p,h) hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1))
! holes :: direct terms hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
do i0 = 1, n_occ_ab_hole(1) partcl(i,1) = iand(differences(i,1),key_i(i,1))
i = occ_hole(i0,1) partcl(i,2) = iand(differences(i,2),key_i(i,2))
htwoe -= buffer_c(i) enddo
enddo
do i0 = 1, n_occ_ab_hole(2) call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint)
i = occ_hole(i0,2) call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint)
htwoe -= buffer_c(i) hmono = mo_bi_ortho_tc_one_e(p,h)
enddo htwoe = fock_op_2_e_tc_closed_shell(p,h)
! holes :: exchange terms ! holes :: direct terms
do i0 = 1, n_occ_ab_hole(spin) do i0 = 1, n_occ_ab_hole(1)
i = occ_hole(i0,spin) i = occ_hole(i0,1)
htwoe += buffer_x(i) htwoe -= buffer_c(i)
enddo enddo
do i0 = 1, n_occ_ab_hole(2)
! particles :: direct terms i = occ_hole(i0,2)
do i0 = 1, n_occ_ab_partcl(1) htwoe -= buffer_c(i)
i = occ_partcl(i0,1) enddo
htwoe += buffer_c(i)
enddo ! holes :: exchange terms
do i0 = 1, n_occ_ab_partcl(2) do i0 = 1, n_occ_ab_hole(spin)
i = occ_partcl(i0,2) i = occ_hole(i0,spin)
htwoe += buffer_c(i) htwoe += buffer_x(i)
enddo enddo
! particles :: exchange terms ! particles :: direct terms
do i0 = 1, n_occ_ab_partcl(spin) do i0 = 1, n_occ_ab_partcl(1)
i = occ_partcl(i0,spin) i = occ_partcl(i0,1)
htwoe -= buffer_x(i) htwoe += buffer_c(i)
enddo enddo
htwoe = htwoe * phase do i0 = 1, n_occ_ab_partcl(2)
hmono = hmono * phase i = occ_partcl(i0,2)
htot = htwoe + hmono 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 end

13
plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f
View File

@ -7,6 +7,10 @@ program tc_bi_ortho
! !
END_DOC END_DOC
implicit none
PROVIDE N_int
my_grid_becke = .True. my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r my_n_pt_r_grid = tc_grid1_r
@ -66,6 +70,15 @@ subroutine routine_diag()
! provide overlap_bi_ortho ! provide overlap_bi_ortho
! provide htilde_matrix_elmt_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 if(N_states .eq. 1) then
print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1) print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1)

28
plugins/local/tc_bi_ortho/tc_hmat.irp.f
View File

@ -13,16 +13,34 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)]
implicit none implicit none
integer :: i, j integer :: i, j
double precision :: t1, t2
double precision :: htot 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 i = 1
j = 1 j = 1
call htilde_mu_mat_opt_bi_ortho_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot) 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 SHARED (N_det, psi_det, N_int, htilde_matrix_elmt_bi_ortho)
!$OMP DO
do i = 1, N_det do i = 1, N_det
do j = 1, N_det do j = 1, N_det
! < J |Htilde | I > ! < 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 htilde_matrix_elmt_bi_ortho(j,i) = htot
enddo enddo
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 END_PROVIDER