LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 11:33:29 +01:00
Code Releases Activity

Merge branch 'dev' of https://github.com/eginer/qp2 into dev

Browse Source
This commit is contained in:
Anthony Scemama 2023-02-02 11:06:15 +01:00
commit e58c42ebdb
9 changed files with 1003 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

35
src/cipsi_tc_bi_ortho/get_d.irp.f
View File

@ -194,7 +194,7 @@ end subroutine get_d3_h
! --- ! ---
subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, coefs) subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
use bitmasks use bitmasks
implicit none implicit none
@ -203,7 +203,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s
integer(bit_kind), intent(in) :: phasemask(N_int,2) integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2) logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
double precision, intent(in) :: coefs(N_states,2) double precision, intent(in) :: coefs(N_states,2)
double precision, intent(inout) :: mat_p(N_states, mo_num, mo_num), mat_m(N_states, mo_num, mo_num) double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, external :: get_phase_bi double precision, external :: get_phase_bi
@ -222,7 +222,8 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s
tip = p(0,1) * p(0,2) tip = p(0,1) * p(0,2)
ma = sp ma = sp
print*,'in get d2' print*,'in get_d2'
stop
if(p(0,1) > p(0,2)) ma = 1 if(p(0,1) > p(0,2)) ma = 1
if(p(0,1) < p(0,2)) ma = 2 if(p(0,1) < p(0,2)) ma = 2
mi = mod(ma, 2) + 1 mi = mod(ma, 2) + 1
@ -259,14 +260,14 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s
if(ma == 1) then if(ma == 1) then
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
do k=1,N_states do k=1,N_states
mat_p(k, putj, puti) = mat_p(k, putj, puti) + coefs(k,1) * hij mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hij
mat_m(k, putj, puti) = mat_m(k, putj, puti) + coefs(k,2) * hji mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hji
enddo enddo
else else
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
do k=1,N_states do k=1,N_states
mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hij
mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hji
enddo enddo
end if end if
end do end do
@ -290,8 +291,8 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s
hji = hji * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) hji = hji * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
do k=1,N_states do k=1,N_states
mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hij
mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hji
enddo enddo
endif endif
end do end do
@ -323,8 +324,8 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s
hji = hji * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) hji = hji * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
do k=1,N_states do k=1,N_states
mat_p(k, puti, putj) = mat_p(k, puti, putj) +coefs(k,1) * hij mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,1) * hij
mat_m(k, puti, putj) = mat_m(k, puti, putj) +coefs(k,2) * hji mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,2) * hji
enddo enddo
end do end do
end do end do
@ -349,14 +350,14 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s
if (puti < putj) then if (puti < putj) then
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
do k=1,N_states do k=1,N_states
mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hij
mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hji
enddo enddo
else else
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
do k=1,N_states do k=1,N_states
mat_p(k, putj, puti) = mat_p(k, putj, puti) + coefs(k,1) * hij mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hij
mat_m(k, putj, puti) = mat_m(k, putj, puti) + coefs(k,2) * hji mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hji
enddo enddo
endif endif
end do end do
@ -375,8 +376,8 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s
hji = hji * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int) hji = hji * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4) !DIR$ LOOP COUNT AVG(4)
do k=1,N_states do k=1,N_states
mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hij
mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hji
enddo enddo
end if end if
end if end if

139
src/cipsi_tc_bi_ortho/get_d0_good.irp.f Normal file
View File

@ -0,0 +1,139 @@
subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
!todo: indices/conjg should be okay for complex
use bitmasks
implicit none
integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
integer(bit_kind) :: det(N_int, 2)
double precision, intent(in) :: coefs(N_states,2)
double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num)
double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
integer :: i, j, k, s, h1, h2, p1, p2, puti, putj, mm
double precision :: phase
double precision :: hij,hji
double precision, external :: get_phase_bi
logical :: ok
integer, parameter :: bant=1
double precision, allocatable :: hij_cache1(:), hij_cache2(:)
allocate (hij_cache1(mo_num),hij_cache2(mo_num))
double precision, allocatable :: hji_cache1(:), hji_cache2(:)
allocate (hji_cache1(mo_num),hji_cache2(mo_num))
! print*,'in get_d0_new'
! call debug_det(gen,N_int)
! print*,'coefs',coefs(1,:)
if(sp == 3) then ! AB
h1 = p(1,1)
h2 = p(1,2)
do p1=1, mo_num
if(bannedOrb(p1, 1)) cycle
! call get_mo_two_e_integrals_complex(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map)
do mm = 1, mo_num
hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,p1,h2,h1)
hji_cache1(mm) = mo_bi_ortho_tc_two_e(h2,h1,mm,p1)
enddo
!!!!!!!!!! <alpha|H|psi>
do p2=1, mo_num
if(bannedOrb(p2,2)) cycle
if(banned(p1, p2, bant)) cycle ! rentable?
if(p1 == h1 .or. p2 == h2) then
call apply_particles(mask, 1,p1,2,p2, det, ok, N_int)
! call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this
! call i_h_j_complex(det, gen, N_int, hij)
call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hij)
else
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = hij_cache1(p2) * phase
end if
if (hij == (0.d0,0.d0)) cycle
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * hij ! HOTSPOT
enddo
end do
!!!!!!!!!! <phi|H|alpha>
do p2=1, mo_num
if(bannedOrb(p2,2)) cycle
if(banned(p1, p2, bant)) cycle ! rentable?
if(p1 == h1 .or. p2 == h2) then
call apply_particles(mask, 1,p1,2,p2, det, ok, N_int)
! call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this
! call i_h_j_complex(det, gen, N_int, hij)
call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hji)
else
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hji = hji_cache1(p2) * phase
end if
if (hji == (0.d0,0.d0)) cycle
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,2) * hji ! HOTSPOT
enddo
end do
end do
else ! AA BB
p1 = p(1,sp)
p2 = p(2,sp)
do puti=1, mo_num
if(bannedOrb(puti, sp)) cycle
! call get_mo_two_e_integrals_complex(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,puti,p2,p1)
hij_cache2(mm) = mo_bi_ortho_tc_two_e(mm,puti,p1,p2)
hji_cache1(mm) = mo_bi_ortho_tc_two_e(p2,p1,mm,puti)
hji_cache2(mm) = mo_bi_ortho_tc_two_e(p1,p2,mm,puti)
enddo
!!!!!!!!!! <alpha|H|psi>
do putj=puti+1, mo_num
if(bannedOrb(putj, sp)) cycle
if(banned(puti, putj, bant)) cycle ! rentable?
if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
!call i_h_j_complex(gen, det, N_int, hij) ! need to take conjugate of this
! call i_h_j_complex(det, gen, N_int, hij)
call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hij)
if (hij == 0.d0) cycle
else
! hij = (mo_two_e_integral_complex(p1, p2, puti, putj) - mo_two_e_integral_complex(p2, p1, puti, putj))
! hij = (mo_bi_ortho_tc_two_e(p1, p2, puti, putj) - mo_bi_ortho_tc_two_e(p2, p1, puti, putj))
hij = (mo_bi_ortho_tc_two_e(puti, putj, p1, p2) - mo_bi_ortho_tc_two_e(puti, putj, p2, p1))
if (hij == 0.d0) cycle
hij = (hij) * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
end if
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
enddo
end do
!!!!!!!!!! <phi|H|alpha>
do putj=puti+1, mo_num
if(bannedOrb(putj, sp)) cycle
if(banned(puti, putj, bant)) cycle ! rentable?
if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hji)
if (hji == 0.d0) cycle
else
hji = (mo_bi_ortho_tc_two_e( p1, p2, puti, putj) - mo_bi_ortho_tc_two_e( p2, p1, puti, putj))
if (hji == 0.d0) cycle
hji = (hji) * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
end if
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
enddo
end do
end do
end if
deallocate(hij_cache1,hij_cache2)
end

454
src/cipsi_tc_bi_ortho/get_d1_good.irp.f Normal file
View File

@ -0,0 +1,454 @@
subroutine get_d1_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
!todo: indices should be okay for complex?
use bitmasks
implicit none
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
integer(bit_kind) :: det(N_int, 2)
double precision, intent(in) :: coefs(N_states,2)
double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num)
double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, external :: get_phase_bi
double precision, external :: mo_two_e_integral_complex
logical :: ok
logical, allocatable :: lbanned(:,:)
integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j
integer :: hfix, pfix, h1, h2, p1, p2, ib, k, l, mm
integer, parameter :: turn2(2) = (/2,1/)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant
double precision, allocatable :: hij_cache(:,:)
double precision :: hij, tmp_rowij(N_states, mo_num), tmp_rowij2(N_states, mo_num)
double precision, allocatable :: hji_cache(:,:)
double precision :: hji, tmp_rowji(N_states, mo_num), tmp_rowji2(N_states, mo_num)
! PROVIDE mo_integrals_map N_int
! print*,'in get_d1_new'
! call debug_det(gen,N_int)
! print*,'coefs',coefs(1,:)
allocate (lbanned(mo_num, 2))
allocate (hij_cache(mo_num,2))
allocate (hji_cache(mo_num,2))
lbanned = bannedOrb
do i=1, p(0,1)
lbanned(p(i,1), 1) = .true.
end do
do i=1, p(0,2)
lbanned(p(i,2), 2) = .true.
end do
ma = 1
if(p(0,2) >= 2) ma = 2
mi = turn2(ma)
bant = 1
if(sp == 3) then
!move MA
if(ma == 2) bant = 2
puti = p(1,mi)
hfix = h(1,ma)
p1 = p(1,ma)
p2 = p(2,ma)
if(.not. bannedOrb(puti, mi)) then
! call get_mo_two_e_integrals_complex(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,p2)
hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,p1)
hji_cache(mm,1) = mo_bi_ortho_tc_two_e(p1,p2,mm,hfix)
hji_cache(mm,2) = mo_bi_ortho_tc_two_e(p2,p1,mm,hfix)
enddo
!! <alpha|H|psi>
tmp_rowij = 0.d0
do putj=1, hfix-1
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hij = hij_cache(putj,1) - hij_cache(putj,2)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,1)
enddo
endif
end do
do putj=hfix+1, mo_num
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hij = hij_cache(putj,2) - hij_cache(putj,1)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowij(k,putj) = tmp_rowij(k,putj) + hij * coefs(k,1)
enddo
endif
end do
if(ma == 1) then
mat_r(1:N_states,1:mo_num,puti) = mat_r(1:N_states,1:mo_num,puti) + tmp_rowij(1:N_states,1:mo_num)
else
do l=1,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k,puti,l) = mat_r(k,puti,l) + tmp_rowij(k,l)
enddo
enddo
end if
!! <phi|H|alpha>
tmp_rowji = 0.d0
do putj=1, hfix-1
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hji = hji_cache(putj,1) - hji_cache(putj,2)
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,2)
enddo
endif
end do
do putj=hfix+1, mo_num
if(lbanned(putj, ma)) cycle
if(banned(putj, puti,bant)) cycle
hji = hji_cache(putj,2) - hji_cache(putj,1)
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowji(k,putj) = tmp_rowji(k,putj) + hji * coefs(k,2)
enddo
endif
end do
if(ma == 1) then
mat_l(1:N_states,1:mo_num,puti) = mat_l(1:N_states,1:mo_num,puti) + tmp_rowji(1:N_states,1:mo_num)
else
do l=1,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k,puti,l) = mat_l(k,puti,l) + tmp_rowji(k,l)
enddo
enddo
end if
end if
!MOVE MI
pfix = p(1,mi)
tmp_rowij = 0.d0
tmp_rowij2 = 0.d0
tmp_rowji = 0.d0
tmp_rowji2 = 0.d0
! call get_mo_two_e_integrals_complex(hfix,pfix,p1,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(hfix,pfix,p2,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p1)
hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p2)
hji_cache(mm,1) = mo_bi_ortho_tc_two_e(pfix,p1,mm,hfix)
hji_cache(mm,2) = mo_bi_ortho_tc_two_e(pfix,p2,mm,hfix)
enddo
putj = p1
!! <alpha|H|psi>
do puti=1,mo_num !HOT
if(lbanned(puti,mi)) cycle
!p1 fixed
putj = p1
if(.not. banned(putj,puti,bant)) then
hij = hij_cache(puti,2)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,1)
enddo
endif
end if
!
putj = p2
if(.not. banned(putj,puti,bant)) then
hij = hij_cache(puti,1)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
do k=1,N_states
tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,1)
enddo
endif
end if
end do
if(mi == 1) then
mat_r(:,:,p1) = mat_r(:,:,p1) + tmp_rowij(:,:)
mat_r(:,:,p2) = mat_r(:,:,p2) + tmp_rowij2(:,:)
else
do l=1,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k,p1,l) = mat_r(k,p1,l) + tmp_rowij(k,l)
mat_r(k,p2,l) = mat_r(k,p2,l) + tmp_rowij2(k,l)
enddo
enddo
end if
putj = p1
!! <phi|H|alpha>
do puti=1,mo_num !HOT
if(lbanned(puti,mi)) cycle
!p1 fixed
putj = p1
if(.not. banned(putj,puti,bant)) then
hji = hji_cache(puti,2)
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,2)
enddo
endif
end if
!
putj = p2
if(.not. banned(putj,puti,bant)) then
hji = hji_cache(puti,1)
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
do k=1,N_states
tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,2)
enddo
endif
end if
end do
if(mi == 1) then
mat_l(:,:,p1) = mat_l(:,:,p1) + tmp_rowji(:,:)
mat_l(:,:,p2) = mat_l(:,:,p2) + tmp_rowji2(:,:)
else
do l=1,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k,p1,l) = mat_l(k,p1,l) + tmp_rowji(k,l)
mat_l(k,p2,l) = mat_l(k,p2,l) + tmp_rowji2(k,l)
enddo
enddo
end if
else ! sp /= 3
if(p(0,ma) == 3) then
do i=1,3
hfix = h(1,ma)
puti = p(i, ma)
p1 = p(turn3(1,i), ma)
p2 = p(turn3(2,i), ma)
! call get_mo_two_e_integrals_complex(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,p2)
hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,p1)
hji_cache(mm,1) = mo_bi_ortho_tc_two_e(p1,p2,mm,hfix)
hji_cache(mm,2) = mo_bi_ortho_tc_two_e(p2,p1,mm,hfix)
enddo
!! <alpha|H|psi>
tmp_rowij = 0.d0
do putj=1,hfix-1
if(banned(putj,puti,1)) cycle
if(lbanned(putj,ma)) cycle
hij = hij_cache(putj,1) - hij_cache(putj,2)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,1)
endif
end do
do putj=hfix+1,mo_num
if(banned(putj,puti,1)) cycle
if(lbanned(putj,ma)) cycle
hij = hij_cache(putj,2) - hij_cache(putj,1)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
tmp_rowij(:,putj) = tmp_rowij(:,putj) + hij * coefs(:,1)
endif
end do
mat_r(:, :puti-1, puti) = mat_r(:, :puti-1, puti) + tmp_rowij(:,:puti-1)
do l=puti,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, l) = mat_r(k, puti,l) + tmp_rowij(k,l)
enddo
enddo
!! <phi|H|alpha>
tmp_rowji = 0.d0
do putj=1,hfix-1
if(banned(putj,puti,1)) cycle
if(lbanned(putj,ma)) cycle
hji = hji_cache(putj,1) - hji_cache(putj,2)
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,2)
endif
end do
do putj=hfix+1,mo_num
if(banned(putj,puti,1)) cycle
if(lbanned(putj,ma)) cycle
hji = hji_cache(putj,2) - hji_cache(putj,1)
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
tmp_rowji(:,putj) = tmp_rowji(:,putj) + hji * coefs(:,2)
endif
end do
mat_l(:, :puti-1, puti) = mat_l(:, :puti-1, puti) + tmp_rowji(:,:puti-1)
do l=puti,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, puti, l) = mat_l(k, puti,l) + tmp_rowji(k,l)
enddo
enddo
end do
else
hfix = h(1,mi)
pfix = p(1,mi)
p1 = p(1,ma)
p2 = p(2,ma)
tmp_rowij = 0.d0
tmp_rowij2 = 0.d0
tmp_rowji = 0.d0
tmp_rowji2 = 0.d0
! call get_mo_two_e_integrals_complex(hfix,p1,pfix,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
! call get_mo_two_e_integrals_complex(hfix,p2,pfix,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
do mm = 1, mo_num
hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,p1,pfix)
hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,p2,pfix)
hji_cache(mm,1) = mo_bi_ortho_tc_two_e(p1,pfix,mm,hfix)
hji_cache(mm,2) = mo_bi_ortho_tc_two_e(p2,pfix,mm,hfix)
enddo
putj = p2
!! <alpha|H|psi>
do puti=1,mo_num
if(lbanned(puti,ma)) cycle
putj = p2
if(.not. banned(puti,putj,1)) then
hij = hij_cache(puti,1)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowij(k,puti) = tmp_rowij(k,puti) + hij * coefs(k,1)
enddo
endif
end if
putj = p1
if(.not. banned(puti,putj,1)) then
hij = hij_cache(puti,2)
if (hij /= 0.d0) then
hij = hij * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
do k=1,N_states
tmp_rowij2(k,puti) = tmp_rowij2(k,puti) + hij * coefs(k,1)
enddo
endif
end if
end do
mat_r(:,:p2-1,p2) = mat_r(:,:p2-1,p2) + tmp_rowij(:,:p2-1)
do l=p2,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k,p2,l) = mat_r(k,p2,l) + tmp_rowij(k,l)
enddo
enddo
mat_r(:,:p1-1,p1) = mat_r(:,:p1-1,p1) + tmp_rowij2(:,:p1-1)
do l=p1,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k,p1,l) = mat_r(k,p1,l) + tmp_rowij2(k,l)
enddo
enddo
!! <phi|H|alpha>
putj = p2
do puti=1,mo_num
if(lbanned(puti,ma)) cycle
putj = p2
if(.not. banned(puti,putj,1)) then
hji = hji_cache(puti,1)
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
tmp_rowji(k,puti) = tmp_rowji(k,puti) + hji * coefs(k,2)
enddo
endif
end if
putj = p1
if(.not. banned(puti,putj,1)) then
hji = hji_cache(puti,2)
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
do k=1,N_states
tmp_rowji2(k,puti) = tmp_rowji2(k,puti) + hji * coefs(k,2)
enddo
endif
end if
end do
mat_l(:,:p2-1,p2) = mat_l(:,:p2-1,p2) + tmp_rowji(:,:p2-1)
do l=p2,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k,p2,l) = mat_l(k,p2,l) + tmp_rowji(k,l)
enddo
enddo
mat_l(:,:p1-1,p1) = mat_l(:,:p1-1,p1) + tmp_rowji2(:,:p1-1)
do l=p1,mo_num
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k,p1,l) = mat_l(k,p1,l) + tmp_rowji2(k,l)
enddo
enddo
end if
end if
deallocate(lbanned,hij_cache, hji_cache)
!! MONO
if(sp == 3) then
s1 = 1
s2 = 2
else
s1 = sp
s2 = sp
end if
do i1=1,p(0,s1)
ib = 1
if(s1 == s2) ib = i1+1
do i2=ib,p(0,s2)
p1 = p(i1,s1)
p2 = p(i2,s2)
if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
! gen is a selector; mask is ionized generator; det is alpha
! hij is contribution to <psi|H|alpha>
! call i_h_j_complex(gen, det, N_int, hij)
call htilde_mu_mat_opt_bi_ortho_no_3e(det, gen, N_int, hij)
call htilde_mu_mat_opt_bi_ortho_no_3e(gen, det, N_int, hji)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * dconjg(hij)
mat_r(k, p1, p2) = mat_r(k, p1, p2) + coefs(k,1) * hij
mat_l(k, p1, p2) = mat_l(k, p1, p2) + coefs(k,2) * hji
enddo
end do
end do
end

308
src/cipsi_tc_bi_ortho/get_d2_good.irp.f Normal file
View File

@ -0,0 +1,308 @@
subroutine get_d2_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, coefs)
!todo: indices/conjg should be correct for complex
use bitmasks
implicit none
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
double precision, intent(in) :: coefs(N_states,2)
double precision, intent(inout) :: mat_r(N_states, mo_num, mo_num)
double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, external :: get_phase_bi
integer :: i, j, k, tip, ma, mi, puti, putj
integer :: h1, h2, p1, p2, i1, i2
double precision :: phase
double precision :: hij,hji
integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0, 3,4, 0,0, 0,0, 2,4, 1,4, 0,0, 2,3, 1,3, 1,2 /), (/2,3,4/))
integer, parameter :: turn2(2) = (/2, 1/)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant
bant = 1
! print*, 'in get_d2_new'
! call debug_det(gen,N_int)
! print*,'coefs',coefs(1,:)
tip = p(0,1) * p(0,2) ! number of alpha particles times number of beta particles
ma = sp !1:(alpha,alpha); 2:(b,b); 3:(a,b)
if(p(0,1) > p(0,2)) ma = 1 ! more alpha particles than beta particles
if(p(0,1) < p(0,2)) ma = 2 ! fewer alpha particles than beta particles
mi = mod(ma, 2) + 1
if(sp == 3) then ! if one alpha and one beta xhole
!(where xholes refer to the ionizations from the generator, not the holes occupied in the ionized generator)
if(ma == 2) bant = 2 ! if more beta particles than alpha particles
if(tip == 3) then ! if 3 of one particle spin and 1 of the other particle spin
puti = p(1, mi)
if(bannedOrb(puti, mi)) return
h1 = h(1, ma)
h2 = h(2, ma)
!! <alpha|H|psi>
do i = 1, 3 ! loop over all 3 combinations of 2 particles with spin ma
putj = p(i, ma)
if(banned(putj,puti,bant)) cycle
i1 = turn3(1,i)
i2 = turn3(2,i)
p1 = p(i1, ma)
p2 = p(i2, ma)
! |G> = |psi_{gen,i}>
! |G'> = a_{x1} a_{x2} |G>
! |alpha> = a_{puti}^{\dagger} a_{putj}^{\dagger} |G'>
! |alpha> = t_{x1,x2}^{puti,putj} |G>
! hij = <psi_{selectors,i}|H|alpha>
! |alpha> = t_{p1,p2}^{h1,h2}|psi_{selectors,i}>
!todo: <i|H|j> = (<h1,h2|p1,p2> - <h1,h2|p2,p1>) * phase
! <psi|H|j> += dconjg(c_i) * <i|H|j>
! <j|H|i> = (<p1,p2|h1,h2> - <p2,p1|h1,h2>) * phase
! <j|H|psi> += <j|H|i> * c_i
! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2, p1, h1, h2)
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e( h1, h2, p2, p1)
if (hij == 0.d0) cycle
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
if(ma == 1) then ! if particle spins are (alpha,alpha,alpha,beta), then puti is beta and putj is alpha
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,1) * hij
enddo
else ! if particle spins are (beta,beta,beta,alpha), then puti is alpha and putj is beta
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
enddo
end if
end do
!! <phi|H|alpha>
do i = 1, 3 ! loop over all 3 combinations of 2 particles with spin ma
putj = p(i, ma)
if(banned(putj,puti,bant)) cycle
i1 = turn3(1,i)
i2 = turn3(2,i)
p1 = p(i1, ma)
p2 = p(i2, ma)
hji = mo_bi_ortho_tc_two_e(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e( p2, p1, h1, h2)
if (hji == 0.d0) cycle
hji = hji * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
if(ma == 1) then ! if particle spins are (alpha,alpha,alpha,beta), then puti is beta and putj is alpha
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,2) * hji
enddo
else ! if particle spins are (beta,beta,beta,alpha), then puti is alpha and putj is beta
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
enddo
end if
end do
else ! if 2 alpha and 2 beta particles
h1 = h(1,1)
h2 = h(1,2)
!! <alpha|H|psi>
do j = 1,2 ! loop over all 4 combinations of one alpha and one beta particle
putj = p(j, 2)
if(bannedOrb(putj, 2)) cycle
p2 = p(turn2(j), 2)
do i = 1,2
puti = p(i, 1)
if(banned(puti,putj,bant) .or. bannedOrb(puti,1)) cycle
p1 = p(turn2(i), 1)
! hij = <psi_{selectors,i}|H|alpha>
! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2)
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2 )
if (hij /= 0.d0) then
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = hij * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
enddo
endif
end do
end do
!! <phi|H|alpha>
do j = 1,2 ! loop over all 4 combinations of one alpha and one beta particle
putj = p(j, 2)
if(bannedOrb(putj, 2)) cycle
p2 = p(turn2(j), 2)
do i = 1,2
puti = p(i, 1)
if(banned(puti,putj,bant) .or. bannedOrb(puti,1)) cycle
p1 = p(turn2(i), 1)
hji = mo_bi_ortho_tc_two_e( p1, p2, h1, h2)
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
enddo
endif
end do
end do
end if
else ! if holes are (a,a) or (b,b)
if(tip == 0) then ! if particles are (a,a,a,a) or (b,b,b,b)
h1 = h(1, ma)
h2 = h(2, ma)
!! <alpha|H|psi>
do i=1,3
puti = p(i, ma)
if(bannedOrb(puti,ma)) cycle
do j=i+1,4
putj = p(j, ma)
if(bannedOrb(putj,ma)) cycle
if(banned(puti,putj,1)) cycle
i1 = turn2d(1, i, j)
i2 = turn2d(2, i, j)
p1 = p(i1, ma)
p2 = p(i2, ma)
! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2)
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
hij = mo_bi_ortho_tc_two_e(h1, h2, p1, p2) - mo_bi_ortho_tc_two_e(h1, h2, p2,p1 )
if (hij == 0.d0) cycle
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
hij = hij * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,1) * hij
enddo
end do
end do
!! <phi|H|alpha>
do i=1,3
puti = p(i, ma)
if(bannedOrb(puti,ma)) cycle
do j=i+1,4
putj = p(j, ma)
if(bannedOrb(putj,ma)) cycle
if(banned(puti,putj,1)) cycle
i1 = turn2d(1, i, j)
i2 = turn2d(2, i, j)
p1 = p(i1, ma)
p2 = p(i2, ma)
hji = mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1,h1, h2 )
if (hji == 0.d0) cycle
hji = hji * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,2) * hji
enddo
end do
end do
else if(tip == 3) then ! if particles are (a,a,a,b) (ma=1,mi=2) or (a,b,b,b) (ma=2,mi=1)
h1 = h(1, mi)
h2 = h(1, ma)
p1 = p(1, mi)
!! <alpha|H|psi>
do i=1,3
puti = p(turn3(1,i), ma)
if(bannedOrb(puti,ma)) cycle
putj = p(turn3(2,i), ma)
if(bannedOrb(putj,ma)) cycle
if(banned(puti,putj,1)) cycle
p2 = p(i, ma)
! hij = mo_bi_ortho_tc_two_e(p1, p2, h1, h2)
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
hij = mo_bi_ortho_tc_two_e(h1, h2,p1, p2 )
if (hij == 0.d0) cycle
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int)
hij = hij * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int)
if (puti < putj) then
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
enddo
else
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,1) * hij
enddo
endif
end do
!! <phi|H|alpha>
do i=1,3
puti = p(turn3(1,i), ma)
if(bannedOrb(puti,ma)) cycle
putj = p(turn3(2,i), ma)
if(bannedOrb(putj,ma)) cycle
if(banned(puti,putj,1)) cycle
p2 = p(i, ma)
hji = mo_bi_ortho_tc_two_e(p1, p2,h1, h2)
if (hji == 0.d0) cycle
hji = hji * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2, N_int)
if (puti < putj) then
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
enddo
else
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,2) * hji
enddo
endif
end do
else ! tip == 4 (a,a,b,b)
puti = p(1, sp)
putj = p(2, sp)
if(.not. banned(puti,putj,1)) then
p1 = p(1, mi)
p2 = p(2, mi)
h1 = h(1, mi)
h2 = h(2, mi)
!! <alpha|H|psi>
! hij = (mo_bi_ortho_tc_two_e(p1, p2, h1, h2) - mo_bi_ortho_tc_two_e(p2,p1, h1, h2))
!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!
! take the transpose of what's written above because later use the complex conjugate
hij = (mo_bi_ortho_tc_two_e(h1, h2,p1, p2) - mo_bi_ortho_tc_two_e(h1, h2, p2,p1))
if (hij /= 0.d0) then
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
! hij = dconjg(hij) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
hij = hij * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,1) * hij
enddo
end if
!! <phi|H|alpha>
hji = (mo_bi_ortho_tc_two_e(p1, p2,h1, h2) - mo_bi_ortho_tc_two_e( p2,p1, h1, h2))
if (hji /= 0.d0) then
hji = hji * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2, N_int)
!DIR$ LOOP COUNT AVG(4)
do k=1,N_states
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,2) * hji
enddo
end if
end if
end if
end if
end

153
src/cipsi_tc_bi_ortho/selection.irp.f
View File

@ -75,7 +75,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :) integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :)
logical, allocatable :: banned(:,:,:), bannedOrb(:,:) logical, allocatable :: banned(:,:,:), bannedOrb(:,:)
double precision, allocatable :: coef_fullminilist_rev(:,:) double precision, allocatable :: coef_fullminilist_rev(:,:)
double precision, allocatable :: mat(:,:,:), mat_p(:,:,:), mat_m(:,:,:) double precision, allocatable :: mat(:,:,:), mat_l(:,:,:), mat_r(:,:,:)
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
@ -208,7 +208,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
allocate( banned(mo_num, mo_num,2), bannedOrb(mo_num, 2) ) allocate( banned(mo_num, mo_num,2), bannedOrb(mo_num, 2) )
allocate( mat(N_states, mo_num, mo_num) ) allocate( mat(N_states, mo_num, mo_num) )
allocate( mat_p(N_states, mo_num, mo_num), mat_m(N_states, mo_num, mo_num) ) allocate( mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num) )
maskInd = -1 maskInd = -1
do s1 = 1, 2 do s1 = 1, 2
@ -411,9 +411,9 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting) call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
if(fullMatch) cycle if(fullMatch) cycle
call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting, mat_p, mat_m) call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting, mat_l, mat_r)
call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, mat_p, mat_m) call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, mat_l, mat_r)
endif endif
enddo enddo
@ -428,7 +428,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask,particle_mask, fock
deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
deallocate(banned, bannedOrb,mat) deallocate(banned, bannedOrb,mat)
deallocate(mat_p, mat_m) deallocate(mat_l, mat_r)
end subroutine select_singles_and_doubles end subroutine select_singles_and_doubles
@ -488,7 +488,7 @@ end subroutine spot_isinwf
! --- ! ---
subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, mat_p, mat_m) subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, mat_l, mat_r)
BEGIN_DOC BEGIN_DOC
! Computes the contributions A(r,s) by ! Computes the contributions A(r,s) by
@ -504,7 +504,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel) integer(bit_kind),intent(in) :: mask(N_int, 2), det(N_int, 2, N_sel)
logical, intent(inout) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2) logical, intent(inout) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num) double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
double precision, intent(inout) :: mat_p(N_states, mo_num, mo_num), mat_m(N_states, mo_num, mo_num) double precision, intent(inout) :: mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num)
integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt integer :: i, ii, j, k, l, h(0:2,2), p(0:4,2), nt
integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2) integer(bit_kind) :: perMask(N_int, 2), mobMask(N_int, 2), negMask(N_int, 2)
@ -516,8 +516,8 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
mat = 0d0 mat = 0d0
mat_p = 0d0 mat_l = 0d0
mat_m = 0d0 mat_r = 0d0
do i = 1, N_int do i = 1, N_int
negMask(i,1) = not(mask(i,1)) negMask(i,1) = not(mask(i,1))
@ -571,7 +571,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
perMask(j,2) = iand(mask(j,2), not(det(j,2,i))) perMask(j,2) = iand(mask(j,2), not(det(j,2,i)))
end do end do
! call get_d3_h ( det(1,1,i), bannedOrb, banned, mat , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) ) ! call get_d3_h ( det(1,1,i), bannedOrb, banned, mat , mask, p, sp, psi_selectors_coef_transp_tc (1, interesting(i)) )
! call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_m, mat_p, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) & ! call get_d3_htc( det(1,1,i), bannedOrb, banned, mat_r, mat_l, mask, p, sp, psi_selectors_rcoef_bi_orth_transp(1, interesting(i)) &
! , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) ) ! , psi_selectors_lcoef_bi_orth_transp(1, interesting(i)) )
call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int) call bitstring_to_list_in_selection(perMask(1,1), h(1,1), h(0,1), N_int)
@ -579,14 +579,16 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int) call get_mask_phase(psi_det_sorted_tc(1,1,interesting(i)), phasemask,N_int)
if(nt == 4) then if(nt == 4) then
call get_d2 (det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) ! call get_d2 (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
! call get_pm2(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i))) call get_d2_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
! call get_pm2(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i)))
elseif(nt == 3) then elseif(nt == 3) then
call get_d1 (det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) ! call get_d1 (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
! call get_pm1(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i))) call get_d1_new(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
! call get_pm1(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i)))
else else
call get_d0 (det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i))) call get_d0_new (det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, 1, interesting(i)))
! call get_pm0(det(1,1,i), phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i))) ! call get_pm0(det(1,1,i), phasemask, bannedOrb, banned, mat_l, mat_r, mask, h, p, sp, psi_selectors_coef_transp_tc(1, interesting(i)))
endif endif
elseif(nt == 4) then elseif(nt == 4) then
call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int) call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
@ -603,7 +605,7 @@ end subroutine splash_pq
! --- ! ---
subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, mat_p, mat_m) subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, mat_l, mat_r)
use bitmasks use bitmasks
use selection_types use selection_types
@ -611,7 +613,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
integer, intent(in) :: i_generator, sp, h1, h2 integer, intent(in) :: i_generator, sp, h1, h2
double precision, intent(in) :: mat(N_states, mo_num, mo_num) double precision, intent(in) :: mat(N_states, mo_num, mo_num)
double precision, intent(in) :: mat_p(N_states, mo_num, mo_num), mat_m(N_states, mo_num, mo_num) double precision, intent(in) :: mat_l(N_states, mo_num, mo_num), mat_r(N_states, mo_num, mo_num)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num) logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num)
double precision, intent(in) :: fock_diag_tmp(mo_num) double precision, intent(in) :: fock_diag_tmp(mo_num)
double precision, intent(in) :: E0(N_states) double precision, intent(in) :: E0(N_states)
@ -770,69 +772,64 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
! ------------------------------------------- ! -------------------------------------------
istate = 1 istate = 1
call htilde_mu_mat_bi_ortho_tot(det, det, N_int, Hii) ! call htilde_mu_mat_bi_ortho_tot(det, det, N_int, Hii)
double precision :: hmono, htwoe, hthree
call diag_htilde_mu_mat_fock_bi_ortho(N_int, det, hmono, htwoe, hthree, hii)
delta_E = E0(istate) - Hii + E_shift delta_E = E0(istate) - Hii + E_shift
!delta_E = 1.d0
! call get_excitation_degree( HF_bitmask, det, degree, N_int)
double precision :: alpha_h_psi_tmp, psi_h_alpha_tmp double precision :: alpha_h_psi_tmp, psi_h_alpha_tmp
psi_h_alpha_tmp = mat_m(istate, p1, p2) if(debug_tc_pt2 == 1)then !! Using the old version
alpha_h_psi_tmp = mat_p(istate, p1, p2) psi_h_alpha = 0.d0
! alpha_h_psi = 0.d0
psi_h_alpha = 0.d0 do iii = 1, N_det
alpha_h_psi = 0.d0 call htilde_mu_mat_bi_ortho_tot(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
do iii = 1, N_det call htilde_mu_mat_bi_ortho_tot(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha) psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function
call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i) alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function
! psi_h_alpha += i_h_alpha * leigvec_tc_bi_orth(iii,1) enddo
! alpha_h_psi += alpha_h_i * reigvec_tc_bi_orth(iii,1) else if(debug_tc_pt2 == 2)then !! debugging the new version
psi_h_alpha += i_h_alpha * 1.d0 psi_h_alpha_tmp = mat_l(istate, p1, p2) ! new version
alpha_h_psi += alpha_h_i * 1.d0 alpha_h_psi_tmp = mat_r(istate, p1, p2) ! new version
enddo psi_h_alpha = 0.d0
! print*,'---',p1,p2 alpha_h_psi = 0.d0
! call debug_det(det,N_int) do iii = 1, N_det ! old version
! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi call htilde_mu_mat_opt_bi_ortho_no_3e(psi_selectors(1,1,iii), det, N_int, i_h_alpha)
! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp call htilde_mu_mat_opt_bi_ortho_no_3e(det, psi_selectors(1,1,iii), N_int, alpha_h_i)
! if(dabs(psi_h_alpha - psi_h_alpha_tmp).gt.1.d-10 .or. dabs(alpha_h_psi - alpha_h_psi_tmp).gt.1.d-10)then ! psi_h_alpha += i_h_alpha * psi_selectors_coef_tc(iii,2,1) ! left function
! if(dabs(psi_h_alpha_tmp*alpha_h_psi_tmp).gt.1.d+10)then ! alpha_h_psi += alpha_h_i * psi_selectors_coef_tc(iii,1,1) ! right function
if(dabs(psi_h_alpha*alpha_h_psi - psi_h_alpha_tmp*alpha_h_psi_tmp).gt.1.d-10)then psi_h_alpha += i_h_alpha * 1.d0 ! left function
! print*,'---' alpha_h_psi += alpha_h_i * 1.d0 ! right function
! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi enddo
! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp if(dabs(psi_h_alpha).gt.1.d-10.or.dabs(alpha_h_psi).gt.1.d-10)then
call debug_det(det,N_int) if(dabs(psi_h_alpha - psi_h_alpha_tmp).gt.1.d-08 .or. dabs(alpha_h_psi - alpha_h_psi_tmp).gt.1.d-08)then
print*,dabs(psi_h_alpha*alpha_h_psi - psi_h_alpha_tmp*alpha_h_psi_tmp),psi_h_alpha *alpha_h_psi,psi_h_alpha_tmp*alpha_h_psi_tmp call debug_det(det,N_int)
print*,'-- Good ' print*,'psi_h_alpha,alpha_h_psi'
print*, psi_h_alpha, alpha_h_psi print*,psi_h_alpha,alpha_h_psi
print*,'-- bad ' print*,psi_h_alpha_tmp,alpha_h_psi_tmp
print*,psi_h_alpha_tmp,alpha_h_psi_tmp print*,dabs(psi_h_alpha - psi_h_alpha_tmp),dabs(alpha_h_psi - alpha_h_psi_tmp)
print*,'-- details good' alpha_h_psi = 0.d0
double precision :: accu_1, accu_2 psi_h_alpha = 0.d0
accu_1 = 0.d0 do iii = 1, N_det
accu_2 = 0.d0
do iii = 1, N_det call get_excitation_degree( psi_det(1,1,iii), det, degree, N_int)
call get_excitation_degree( psi_det(1,1,iii), det, degree, N_int) call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha)
call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha) call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i)
call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i) alpha_h_psi += alpha_h_i
print*,iii,degree,i_h_alpha,alpha_h_i psi_h_alpha += i_h_alpha
accu_1 += i_h_alpha if(dabs(i_h_alpha).gt.1.d-10.or.dabs(alpha_h_i).gt.1.d-10)then
accu_2 += alpha_h_i call debug_det(psi_det(1,1,iii),N_int)
print*,accu_1,accu_2 print*,iii,degree,i_h_alpha,alpha_h_i
print*,psi_h_alpha,alpha_h_psi
enddo ! print*,leigvec_tc_bi_orth(iii,1),reigvec_tc_bi_orth(iii,1)
! if(dabs(psi_h_alpha*alpha_h_psi).gt.1.d-10)then endif
! print*,p1,p2 enddo
! print*,det(1,1), det(1,2) stop
! call debug_det(det,N_int) endif
! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi endif
! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp else
! print*, dabs(psi_h_alpha*alpha_h_psi - psi_h_alpha_tmp*alpha_h_psi_tmp),& psi_h_alpha = mat_l(istate, p1, p2)
! psi_h_alpha *alpha_h_psi,psi_h_alpha_tmp*alpha_h_psi_tmp alpha_h_psi = mat_r(istate, p1, p2)
stop endif
endif
! endif
! stop
! endif
!if(alpha_h_psi*psi_h_alpha/delta_E.gt.1.d-10)then !if(alpha_h_psi*psi_h_alpha/delta_E.gt.1.d-10)then
! print*, 'E0,Hii,E_shift' ! print*, 'E0,Hii,E_shift'

2
src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f
View File

@ -90,7 +90,7 @@ subroutine run_stochastic_cipsi
call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states) call pt2_alloc(pt2_data_err, N_states)
call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection
stop ! stop
N_iter += 1 N_iter += 1

4
src/fci_tc_bi/selectors.irp.f
View File

@ -48,8 +48,8 @@ END_PROVIDER
do k=1,N_states do k=1,N_states
do i=1,N_det_selectors do i=1,N_det_selectors
psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k) psi_selectors_coef(i,k) = psi_coef_sorted_tc_gen(i,k)
! psi_selectors_coef_tc(i,1,k) = psi_r_coef_sorted_bi_ortho(i,k) ! psi_selectors_coef_tc(i,1,k) = psi_l_coef_sorted_bi_ortho(i,k)
! psi_selectors_coef_tc(i,2,k) = psi_l_coef_sorted_bi_ortho(i,k) ! psi_selectors_coef_tc(i,2,k) = psi_r_coef_sorted_bi_ortho(i,k)
psi_selectors_coef_tc(i,1,k) = 1.d0 psi_selectors_coef_tc(i,1,k) = 1.d0
psi_selectors_coef_tc(i,2,k) = 1.d0 psi_selectors_coef_tc(i,2,k) = 1.d0
enddo enddo

1
src/mo_two_e_ints/map_integrals.irp.f
View File

@ -309,6 +309,7 @@ end
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call map_get(mo_integrals_map,idx,tmp) call map_get(mo_integrals_map,idx,tmp)
banned_excitation(i,j) = dabs(tmp) < 1.d-14 banned_excitation(i,j) = dabs(tmp) < 1.d-14
! banned_excitation(i,j) = .False.
banned_excitation(j,i) = banned_excitation(i,j) banned_excitation(j,i) = banned_excitation(i,j)
if (banned_excitation(i,j)) icount = icount+2 if (banned_excitation(i,j)) icount = icount+2
enddo enddo

5
src/tc_keywords/EZFIO.cfg
View File

@ -184,3 +184,8 @@ doc: Thresholds on the Imag part of energy
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: 1.e-7 default: 1.e-7
[debug_tc_pt2]
type: integer
doc: If :: 1 then you compute the TC-PT2 the old way, :: 2 then you check with the new version but without three-body
interface: ezfio,provider,ocaml
default: -1