diff --git a/src/cipsi_tc_bi_ortho/get_d.irp.f b/src/cipsi_tc_bi_ortho/get_d.irp.f index 7fdc5e12..9421787e 100644 --- a/src/cipsi_tc_bi_ortho/get_d.irp.f +++ b/src/cipsi_tc_bi_ortho/get_d.irp.f @@ -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 diff --git a/src/cipsi_tc_bi_ortho/get_d0_good.irp.f b/src/cipsi_tc_bi_ortho/get_d0_good.irp.f index b390ac5c..59068e46 100644 --- a/src/cipsi_tc_bi_ortho/get_d0_good.irp.f +++ b/src/cipsi_tc_bi_ortho/get_d0_good.irp.f @@ -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) diff --git a/src/cipsi_tc_bi_ortho/get_d2_good.irp.f b/src/cipsi_tc_bi_ortho/get_d2_good.irp.f new file mode 100644 index 00000000..0e17e306 --- /dev/null +++ b/src/cipsi_tc_bi_ortho/get_d2_good.irp.f @@ -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) + + !! + 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 = + ! |alpha> = t_{p1,p2}^{h1,h2}|psi_{selectors,i}> + !todo: = ( - ) * phase + ! += dconjg(c_i) * + ! = ( - ) * phase + ! += * 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 instead of +! 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 + !! + 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 = + ! |alpha> = t_{p1,p2}^{h1,h2}|psi_{selectors,i}> + !todo: = ( - ) * phase + ! += dconjg(c_i) * + ! = ( - ) * phase + ! += * 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 instead of +! 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) + !! + 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 = +! 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 instead of +! 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 + !! + 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) + !! + 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 instead of +! 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 + !! + 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) + !! + 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 instead of +! 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 + !! + 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) + !! +! 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 instead of +! 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 + !! + 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 instead of +! 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 diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f index 33043e74..51853fcd 100644 --- a/src/cipsi_tc_bi_ortho/selection.irp.f +++ b/src/cipsi_tc_bi_ortho/selection.irp.f @@ -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' diff --git a/src/fci_tc_bi/fci_tc_bi_ortho.irp.f b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f index 84ac8166..f15c14da 100644 --- a/src/fci_tc_bi/fci_tc_bi_ortho.irp.f +++ b/src/fci_tc_bi/fci_tc_bi_ortho.irp.f @@ -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