mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-22 03:23:29 +01:00
work in progress for get_d2_good
This commit is contained in:
eginer
parent
44ac34bc03
commit
b0d49f6351
@ -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
|
||||
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)
|
||||
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_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
|
||||
|
||||
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)
|
||||
|
||||
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 = 2
|
||||
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
|
||||
!DIR$ LOOP COUNT AVG(4)
|
||||
do k=1,N_states
|
||||
mat_p(k, putj, puti) = mat_p(k, putj, puti) + coefs(k,1) * hij
|
||||
mat_m(k, putj, puti) = mat_m(k, putj, puti) + coefs(k,2) * hji
|
||||
mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hij
|
||||
mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hji
|
||||
enddo
|
||||
else
|
||||
!DIR$ LOOP COUNT AVG(4)
|
||||
do k=1,N_states
|
||||
mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij
|
||||
mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji
|
||||
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hij
|
||||
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hji
|
||||
enddo
|
||||
end if
|
||||
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)
|
||||
!DIR$ LOOP COUNT AVG(4)
|
||||
do k=1,N_states
|
||||
mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij
|
||||
mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji
|
||||
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hij
|
||||
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hji
|
||||
enddo
|
||||
endif
|
||||
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)
|
||||
!DIR$ LOOP COUNT AVG(4)
|
||||
do k=1,N_states
|
||||
mat_p(k, puti, putj) = mat_p(k, puti, putj) +coefs(k,1) * hij
|
||||
mat_m(k, puti, putj) = mat_m(k, puti, putj) +coefs(k,2) * hji
|
||||
mat_l(k, puti, putj) = mat_l(k, puti, putj) +coefs(k,1) * hij
|
||||
mat_r(k, puti, putj) = mat_r(k, puti, putj) +coefs(k,2) * hji
|
||||
enddo
|
||||
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
|
||||
!DIR$ LOOP COUNT AVG(4)
|
||||
do k=1,N_states
|
||||
mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij
|
||||
mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji
|
||||
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hij
|
||||
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hji
|
||||
enddo
|
||||
else
|
||||
!DIR$ LOOP COUNT AVG(4)
|
||||
do k=1,N_states
|
||||
mat_p(k, putj, puti) = mat_p(k, putj, puti) + coefs(k,1) * hij
|
||||
mat_m(k, putj, puti) = mat_m(k, putj, puti) + coefs(k,2) * hji
|
||||
mat_l(k, putj, puti) = mat_l(k, putj, puti) + coefs(k,1) * hij
|
||||
mat_r(k, putj, puti) = mat_r(k, putj, puti) + coefs(k,2) * hji
|
||||
enddo
|
||||
endif
|
||||
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)
|
||||
!DIR$ LOOP COUNT AVG(4)
|
||||
do k=1,N_states
|
||||
mat_p(k, puti, putj) = mat_p(k, puti, putj) + coefs(k,1) * hij
|
||||
mat_m(k, puti, putj) = mat_m(k, puti, putj) + coefs(k,2) * hji
|
||||
mat_l(k, puti, putj) = mat_l(k, puti, putj) + coefs(k,1) * hij
|
||||
mat_r(k, puti, putj) = mat_r(k, puti, putj) + coefs(k,2) * hji
|
||||
enddo
|
||||
end if
|
||||
end if
|
||||
|
||||
@ -40,13 +40,13 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
|
||||
if(bannedOrb(p2,2)) cycle
|
||||
if(banned(p1, p2, bant)) cycle ! rentable?
|
||||
if(p1 == h1 .or. p2 == h2) then
|
||||
print*,'in hij 1'
|
||||
! print*,'in hij 1'
|
||||
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
|
||||
print*,'in chelou 1 !!!!!!!!!!!!!!!!!!!!!!!'
|
||||
! print*,'in chelou 1 !!!!!!!!!!!!!!!!!!!!!!!'
|
||||
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
|
||||
hij = hij_cache1(p2) * phase
|
||||
end if
|
||||
@ -61,13 +61,13 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
|
||||
if(bannedOrb(p2,2)) cycle
|
||||
if(banned(p1, p2, bant)) cycle ! rentable?
|
||||
if(p1 == h1 .or. p2 == h2) then
|
||||
print*,'in hji 1'
|
||||
! print*,'in hji 1'
|
||||
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
|
||||
print*,'in chelou 1 ji !!!!!!!!!!!!!!!!!!!!!!!'
|
||||
! print*,'in chelou 1 ji !!!!!!!!!!!!!!!!!!!!!!!'
|
||||
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
|
||||
hji = hji_cache1(p2) * phase
|
||||
end if
|
||||
@ -97,14 +97,14 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
|
||||
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
|
||||
print*,'in hij 2'
|
||||
! print*,'in hij 2'
|
||||
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,0.d0)) cycle
|
||||
else
|
||||
print*,'in chelou 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
! print*,'in chelou 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
! 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))
|
||||
@ -122,12 +122,12 @@ subroutine get_d0_new(gen, phasemask, bannedOrb, banned, mat_l, mat_r, mask, h,
|
||||
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
|
||||
print*,'in hji 2'
|
||||
call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
|
||||
! print*,'in hji 2'
|
||||
! 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,0.d0)) cycle
|
||||
else
|
||||
print*,'in chelou 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
! print*,'in chelou 2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
|
||||
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,0.d0)) cycle
|
||||
hji = (hji) * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
|
||||
|
||||
328
src/cipsi_tc_bi_ortho/get_d2_good.irp.f
Normal file
328
src/cipsi_tc_bi_ortho/get_d2_good.irp.f
Normal file
@ -0,0 +1,328 @@
|
||||
|
||||
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'
|
||||
|
||||
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)
|
||||
|
||||
! |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>
|
||||
! hji = <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
|
||||
! hji = 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
|
||||
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
|
||||
|
||||
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
|
||||
! hji = dconjg(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)
|
||||
|
||||
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,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(h1, h2,p1, p2) - mo_bi_ortho_tc_two_e(h1, h2, p2,p1))
|
||||
if (hji /= (0.d0,0.d0)) then
|
||||
! take conjugate to get contribution to <alpha|H|psi> instead of <psi|H|alpha>
|
||||
! hji = dconjg(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)
|
||||
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
|
||||
@ -579,11 +579,12 @@ 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)
|
||||
if(nt == 4) then
|
||||
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_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_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_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
|
||||
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_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_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
|
||||
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)))
|
||||
@ -786,8 +787,8 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
|
||||
do iii = 1, N_det
|
||||
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)
|
||||
!!! psi_h_alpha += i_h_alpha * leigvec_tc_bi_orth(iii,1)
|
||||
!!! alpha_h_psi += alpha_h_i * reigvec_tc_bi_orth(iii,1)
|
||||
!! psi_h_alpha += i_h_alpha * leigvec_tc_bi_orth(iii,1)
|
||||
!! alpha_h_psi += alpha_h_i * reigvec_tc_bi_orth(iii,1)
|
||||
psi_h_alpha += i_h_alpha * 1.d0
|
||||
alpha_h_psi += alpha_h_i * 1.d0
|
||||
enddo
|
||||
@ -795,20 +796,25 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
|
||||
!!! call debug_det(det,N_int)
|
||||
!!! print*,psi_h_alpha *alpha_h_psi, psi_h_alpha, alpha_h_psi
|
||||
!!! print*,psi_h_alpha_tmp*alpha_h_psi_tmp,psi_h_alpha_tmp,alpha_h_psi_tmp
|
||||
!!! 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
|
||||
!!! if(dabs(psi_h_alpha_tmp*alpha_h_psi_tmp).gt.1.d+10)then
|
||||
if(dabs(psi_h_alpha).gt.1.d-10.or.dabs(alpha_h_psi).gt.1.d-10)then
|
||||
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
|
||||
! if(dabs(psi_h_alpha_tmp).gt.1.d-10.or.dabs(alpha_h_psi_tmp).gt.1.d-10)then
|
||||
if(dabs(alpha_h_psi_tmp).gt.1.d-10)then
|
||||
if(degree==2)then
|
||||
! if(degree==2)then
|
||||
call debug_det(det,N_int)
|
||||
print*,'psi_h_alpha,alpha_h_psi'
|
||||
print*,psi_h_alpha,alpha_h_psi
|
||||
print*,psi_h_alpha_tmp, alpha_h_psi_tmp
|
||||
endif
|
||||
endif
|
||||
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)
|
||||
do iii = 1, N_det
|
||||
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(det, psi_det(1,1,iii), N_int, alpha_h_i)
|
||||
print*,iii,degree,i_h_alpha,alpha_h_i
|
||||
enddo
|
||||
stop
|
||||
endif
|
||||
endif
|
||||
! stop
|
||||
! endif
|
||||
|
||||
!if(alpha_h_psi*psi_h_alpha/delta_E.gt.1.d-10)then
|
||||
! print*, 'E0,Hii,E_shift'
|
||||
|
||||
@ -38,8 +38,8 @@ program fci
|
||||
|
||||
|
||||
my_grid_becke = .True.
|
||||
my_n_pt_r_grid = 30
|
||||
my_n_pt_a_grid = 50
|
||||
my_n_pt_r_grid = 10
|
||||
my_n_pt_a_grid = 14
|
||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||
pruning = -1.d0
|
||||
touch pruning
|
||||
|
||||
Loading…
Reference in New Issue
Block a user