diff --git a/plugins/MRPT_Utils/new_way_second_order_coef.irp.f b/plugins/MRPT_Utils/new_way_second_order_coef.irp.f index ae356e11..5c4b562f 100644 --- a/plugins/MRPT_Utils/new_way_second_order_coef.irp.f +++ b/plugins/MRPT_Utils/new_way_second_order_coef.irp.f @@ -19,6 +19,7 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p) double precision :: get_mo_bielec_integral_schwartz double precision :: active_int(n_act_orb,2) double precision :: hij,phase + integer :: index_orb_act_mono(N_det,6) !matrix_2h1p = 0.d0 elec_num_tab_local = 0 @@ -47,10 +48,12 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p) integer :: idx(0:N_det) double precision :: delta_e(n_act_orb,2,N_states) integer :: istate - integer :: index_orb_act_mono(N_det,3) do idet = 1, N_det - call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + call get_excitation_degree_vector_mono_or_exchange(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) +! if(idet == 81)then +! call get_excitation_degree_vector_mono_or_exchange_verbose(psi_det(1,1,1),psi_det(1,1,idet),degree,N_int,N_det,idx) +! endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) do jspin = 1, 2 ! spin of the couple z-a^dagger (j,a) @@ -109,82 +112,52 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p) enddo enddo enddo -!!!!!!!!!!!!!!!!!!!!!!!!!!!! Second order coefficient : interactions between the perturbers throw the active space - do a = 1, n_act_orb - do jspin = 1, 2 - do ispin = 1, 2 - if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle - ! determinant perturber | det_tmp > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{a,jspin} a_{i,ispin} | Idet > - do inint = 1, N_int - det_tmp(inint,1) = iand(perturb_dets(inint,1,a,jspin,ispin),cas_bitmask(inint,1,1)) - det_tmp(inint,2) = iand(perturb_dets(inint,2,a,jspin,ispin),cas_bitmask(inint,1,1)) - enddo - do istate = 1, N_states - coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0 - enddo - do b = 1, n_act_orb - do kspin = jspin , jspin - integer :: degree_scalar - if( perturb_dets_phase(b,kspin,ispin) .le. -10.d0)cycle - do inint = 1, N_int - det_tmp_j(inint,1) = iand(perturb_dets(inint,1,b,kspin,ispin),cas_bitmask(inint,1,1)) - det_tmp_j(inint,2) = iand(perturb_dets(inint,2,b,kspin,ispin),cas_bitmask(inint,1,1)) - enddo - call get_excitation_degree(det_tmp,det_tmp_j,degree_scalar,N_int) - if (degree_scalar > 2 .or. degree_scalar == 0)cycle - ! determinant perturber | det_tmp_j > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{b,jspin} a_{i,ispin} | Idet > -! print*, '**********************' -! integer(bit_kind) :: det_bis(N_int,2) -! call debug_det(det_tmp,N_int) -! call debug_det(det_tmp_j,N_int) -! do inint = 1, N_int -! det_bis(inint,1) = perturb_dets(inint,1,b,kspin,ispin) -! det_bis(inint,2) = perturb_dets(inint,2,b,kspin,ispin) -! enddo -! call debug_det(det_bis,N_int) - call i_H_j_dyall(det_tmp,det_tmp_j,N_int,hij) - do istate = 1, N_states - coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(b,kspin,ispin,istate,1) & - * hij / delta_e(a,jspin,istate) - if(dabs(hij).gt.0.01d0)then - print*,degree_scalar, hij - print*, coef_perturb_from_idet(b,kspin,ispin,istate,1)* hij / delta_e(a,jspin,istate),coef_perturb_from_idet(a,jspin,ispin,istate,1) - - endif - enddo - enddo - enddo - enddo - enddo - enddo - do a = 1, n_act_orb - do jspin = 1, 2 - do ispin = 1, 2 - if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle - do istate = 1, N_states -! print*, coef_perturb_from_idet(a,jspin,ispin,istate,1),coef_perturb_from_idet(a,jspin,ispin,istate,2) - coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(a,jspin,ispin,istate,1) - enddo - enddo - enddo - enddo -! stop + + !!!!!!!!!!!!!!!!!!!!!!!!!!! determination of the connections between I and the other J determinants mono excited in the CAS !!!!!!!!!!!!!!!!!!!!!!!!!!!! the determinants I and J must be connected by the following operator -!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!! + integer :: i_hole,i_part + double precision :: hij_test + double precision :: fock_operator_local(n_act_orb,n_act_orb,2) do jdet = 1, idx(0) if(idx(jdet).ne.idet)then - call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) - if (exc(0,1,1) == 1) then - ! Mono alpha - index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_a - index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,1,1)) !!! a_{b} - index_orb_act_mono(idx(jdet),3) = 1 - else - ! Mono beta - index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_a - index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,1,2)) !!! a_{b} - index_orb_act_mono(idx(jdet),3) = 2 + if(degree(jdet)==1)then + call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + if (exc(0,1,1) == 1) then + ! Mono alpha + i_hole = list_act_reverse(exc(1,1,1)) !!! a_a + i_part = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b} + kspin = 1 !!! kspin + index_orb_act_mono(idx(jdet),1) = i_hole + index_orb_act_mono(idx(jdet),2) = i_part + index_orb_act_mono(idx(jdet),3) = kspin + call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij) + fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator + fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator + else + ! Mono beta + i_hole = list_act_reverse(exc(1,1,2)) !!! a_a + i_part = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_{b} + kspin = 2 !!! kspin + index_orb_act_mono(idx(jdet),1) = i_hole + index_orb_act_mono(idx(jdet),2) = i_part + index_orb_act_mono(idx(jdet),3) = kspin + call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij) + fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator + fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator + endif + else if(degree(jdet)==2)then + call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + ! Mono alpha + index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,1)) !!! a_a + index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b} + index_orb_act_mono(idx(jdet),3) = 1 + ! Mono beta + index_orb_act_mono(idx(jdet),4) = list_act_reverse(exc(1,1,2)) !!! a_a + index_orb_act_mono(idx(jdet),5) = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_{b} + index_orb_act_mono(idx(jdet),6) = 2 endif else index_orb_act_mono(idx(jdet),1) = -1 @@ -192,56 +165,181 @@ subroutine give_2h1p_contrib_sec_order(matrix_2h1p) enddo integer :: kspin + integer :: corb,i_ok + integer(bit_kind) :: det_tmp_bis(N_int,2) + double precision :: hib , hab , hja + double precision :: delta_e_ab(N_states) + double precision :: hib_test,hja_test,hab_test do jdet = 1, idx(0) if(idx(jdet).ne.idet)then - ! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator - ! are connected by the presence of the perturbers determinants |det_tmp> - aorb = index_orb_act_mono(idx(jdet),1) ! a^{\dagger}_{aorb} - borb = index_orb_act_mono(idx(jdet),2) ! a_{borb} - kspin = index_orb_act_mono(idx(jdet),3) ! spin of the excitation - ! the determinants Idet and Jdet interact throw the following operator - ! | Jdet > = a_{borb,kspin} a^{\dagger}_{aorb, kspin} | Idet > +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE OF THE MONO EXCITATIONS + if(degree(jdet) == 1)then + ! ! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator + ! ! are connected by the presence of the perturbers determinants |det_tmp> + aorb = index_orb_act_mono(idx(jdet),1) ! a^{\dagger}_{aorb} + borb = index_orb_act_mono(idx(jdet),2) ! a_{borb} + kspin = index_orb_act_mono(idx(jdet),3) ! spin of the excitation + do ispin = 1, 2 ! you loop on all possible spin for the excitation + ! a^{\dagger}_r a_{i} (ispin) + ! ! the determinants Idet and Jdet interact throw the following operator + ! ! | Jdet > = a_{borb,kspin} a^{\dagger}_{aorb, kspin} | Idet > + do jspin = 1, 2 + if (jspin .ne. kspin)then - do ispin = 1, 2 ! you loop on all possible spin for the excitation - ! a^{\dagger}_r a_{i} (ispin) - if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count + do corb = 1, n_act_orb + if(perturb_dets_phase(corb,jspin,ispin).le.-100d0)cycle + ! ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{corb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet > + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,corb,jspin,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,corb,jspin,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,corb,jspin,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin) + enddo + ! ! < idet | H | det_tmp > = phase * (ir|cv) + call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) + if(ispin == jspin)then + hib= phase * (active_int(corb,1) - active_int(corb,2)) + else + hib= phase * active_int(corb,1) + endif + + ! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} | det_tmp > + call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok) + if(i_ok .ne. 1)cycle + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) - ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet > - do inint = 1, N_int - det_tmp(inint,1) = perturb_dets(inint,1,aorb,kspin,ispin) - det_tmp(inint,2) = perturb_dets(inint,2,aorb,kspin,ispin) - enddo - double precision :: hja - ! you determine the interaction between the excited determinant and the other parent | Jdet > - ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{borb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Jdet > - ! hja = < det_tmp | H | Jdet > - call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int) - if(kspin == ispin)then - hja = phase * (active_int(borb,2) - active_int(borb,1) ) - else - hja = phase * active_int(borb,1) + ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb} + hab = (fock_operator_local(aorb,borb,kspin) ) * phase + if(isnan(hab))then + print*, '1' + stop + endif + ! < jdet | H | det_tmp_bis > = phase * (ir|cv) + call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int) + if(ispin == jspin)then + hja= phase * (active_int(corb,1) - active_int(corb,2)) + else + hja= phase * (active_int(corb,1)) + endif + do istate = 1, N_states + delta_e_ab(istate) = delta_e(corb,jspin,istate) + one_anhil_one_creat(borb,aorb,kspin,kspin,istate) + matrix_2h1p(idx(jdet),idet,istate) = matrix_2h1p(idx(jdet),idet,istate) + & + hib / delta_e(corb,jspin,istate) * hab / delta_e_ab(istate) * hja + ! ! < det_tmp | H | Idet > / delta_E (Idet --> det_tmp ) + ! ! < det_tmp | H | det_tmp_bis > / delta_E (Idet --> det_tmp --> det_tmp_bis) + ! ! < det_tmp_bis | H | Jdet > + enddo + enddo ! corb + else + if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count + do corb = 1, n_act_orb + if(corb == aorb .or. corb == borb) cycle + if(perturb_dets_phase(corb,jspin,ispin).le.-100d0)cycle + ! ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{corb,jspin} a_{iorb,ispin} | Idet > + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,corb,jspin,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,corb,jspin,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,corb,jspin,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin) + enddo + ! < idet | H | det_tmp > = phase * ( (ir|cv) - (iv|cr) ) + call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) + if(ispin == jspin)then + hib= phase * (active_int(corb,1) - active_int(corb,2)) + else + hib= phase * active_int(corb,1) + endif + ! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} | det_tmp > + call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok) + if(i_ok .ne. 1)cycle + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) +! ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb} + hab = fock_operator_local(aorb,borb,kspin) * phase + if(isnan(hab))then + print*, '2' + stop + endif + ! < jdet | H | det_tmp_bis > = phase * ( (ir|cv) - (iv|cr) ) + call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int) + if(ispin == jspin)then + hja= phase * (active_int(corb,1) - active_int(corb,2)) + else + hja= phase * (active_int(corb,1)) + endif + do istate = 1, N_states + delta_e_ab(istate) = delta_e(corb,jspin,istate) + one_anhil_one_creat(borb,aorb,kspin,kspin,istate) + matrix_2h1p(idx(jdet),idet,istate) = matrix_2h1p(idx(jdet),idet,istate) + & + hib / delta_e(corb,jspin,istate) * hab / delta_e_ab(istate) * hja + ! ! < det_tmp | H | Idet > / delta_E (Idet --> det_tmp ) + ! ! < det_tmp | H | det_tmp_bis > / delta_E (Idet --> det_tmp --> det_tmp_bis) + ! ! < det_tmp_bis | H | Jdet > + enddo + enddo ! corb + endif + enddo + enddo + ! + else !! Double excitation operators + ! + if (index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),5))then !! spin exchange + do ispin = 1, 2 ! you loop on all possible spin for the excitation + ! a^{\dagger}_r a_{i} (ispin) + !!! ! first combination of spin :: | det_tmp > = a^{\dagger}_{aorb,beta} | Idet > + jspin = 2 + aorb = index_orb_act_mono(idx(jdet),1) ! hole of the alpha electron + borb = index_orb_act_mono(idx(jdet),2) ! particle of the alpha electron + if(perturb_dets_phase(aorb,jspin,ispin).le.-100d0)cycle + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin) + enddo + ! | det_tmp > = a^{\dagger}_{aorb,beta} | Idet > + call get_double_excitation(det_tmp,psi_det(1,1,idet),exc,phase,N_int) + if(ispin == jspin)then + hib= phase * (active_int(aorb,1) - active_int(aorb,2)) + else + hib= phase * (active_int(aorb,1)) + endif + if(hib .ne. perturb_dets_hij(aorb,jspin,ispin))then + print*, 'pb !!' + print*, 'hib .ne. perturb_dets_hij(aorb,jspin,ispin)' + stop + endif + enddo !! ispin + + else if(index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),4))then !! closed shell double excitation + + else + call get_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,degree_scalar,phase,N_int) + integer :: h1,h2,p1,p2,s1,s2 , degree_scalar + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + print*, h1,p1,h2,p2,s1,s2 + call debug_det(psi_det(1,1,idet),N_int) + call debug_det(psi_det(1,1,idx(jdet)),N_int) + print*, idet,idx(jdet) + print*, 'pb !!!!!!!!!!!!!' + call get_excitation_degree_vector_mono_or_exchange_verbose(psi_det(1,1,1),psi_det(1,1,idet),degree,N_int,N_det,idx) + stop + endif endif - do istate = 1, N_states - matrix_2h1p(idx(jdet),idet,istate) += hja * coef_perturb_from_idet(aorb,kspin,ispin,istate,2) - enddo - enddo ! ispin - else - ! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations - ! - ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet > - do ispin = 1, 2 - do kspin = 1, 2 - if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count - do a = 1, n_act_orb ! First active - do istate = 1, N_states - matrix_2h1p(idet,idet,istate) += coef_perturb_from_idet(a,kspin,ispin,istate,2) * perturb_dets_hij(a,kspin,ispin) - enddo - enddo - enddo - enddo - + !! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations + !! + !! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet > + !!do ispin = 1, 2 + !! do kspin = 1, 2 + !! if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count + !! do a = 1, n_act_orb ! First active + !! do istate = 1, N_states + !! matrix_2h1p(idet,idet,istate) += coef_perturb_from_idet(a,kspin,ispin,istate,2) * perturb_dets_hij(a,kspin,ispin) + !! enddo + !! enddo + !! enddo + !!enddo + ! endif enddo @@ -261,8 +359,8 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) use bitmasks implicit none double precision , intent(inout) :: matrix_1h2p(N_det,N_det,*) - integer :: i,v,r,a,b - integer :: iorb, vorb, rorb, aorb, borb + integer :: i,v,r,a,b,c + integer :: iorb, vorb, rorb, aorb, borb,corb integer :: ispin,jspin integer :: idet,jdet integer(bit_kind) :: perturb_dets(N_int,2,n_act_orb,2,2) @@ -292,6 +390,9 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) double precision :: hja double precision :: contrib_hij double precision :: fock_operator_local(n_act_orb,n_act_orb,2) + double precision :: fock_operator_from_core(n_act_orb,n_act_orb) + double precision :: fock_operator_from_virt(n_act_orb,n_act_orb) + double precision :: fock_operator_from_act(n_act_orb,n_act_orb,n_act_orb,2) accu_contrib = 0.d0 !matrix_1h2p = 0.d0 @@ -394,10 +495,11 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) !!!!!!!!!!!!!!!!!!!!!!!!!!! determination of the connections between I and the other J determinants mono excited in the CAS !!!!!!!!!!!!!!!!!!!!!!!!!!!! the determinants I and J must be connected by the following operator !!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!! integer :: i_hole,i_part + double precision :: hij_test do jdet = 1, idx(0) if(idx(jdet).ne.idet)then -! print*, degree(jdet) if(degree(jdet)==1)then call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) if (exc(0,1,1) == 1) then @@ -408,7 +510,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) index_orb_act_mono(idx(jdet),1) = i_hole index_orb_act_mono(idx(jdet),2) = i_part index_orb_act_mono(idx(jdet),3) = kspin - call i_H_j_dyall(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),N_int,hij) + call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij) fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator else @@ -419,7 +521,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) index_orb_act_mono(idx(jdet),1) = i_hole index_orb_act_mono(idx(jdet),2) = i_part index_orb_act_mono(idx(jdet),3) = kspin - call i_H_j_dyall(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),N_int,hij) + call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij) fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator endif @@ -439,11 +541,17 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) endif enddo + integer ::dorb,i_ok + integer(bit_kind) :: det_tmp_bis(N_int,2) + double precision :: hib , hab + double precision :: delta_e_ab(N_states) + double precision :: hib_test,hja_test,hab_test do jdet = 1, idx(0) if(idx(jdet).ne.idet)then - if(degree(jdet) == 1)then +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE OF THE MONO EXCITATIONS + if(degree(jdet) == 1)then ! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator ! are connected by the presence of the perturbers determinants |det_tmp> aorb = index_orb_act_mono(idx(jdet),1) ! a_{aorb} @@ -454,11 +562,6 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) do ispin = 1, 2 ! you loop on all possible spin for the excitation ! a^{\dagger}_r a_{i} (ispin) - integer ::corb,dorb,i_ok - integer(bit_kind) :: det_tmp_bis(N_int,2) - double precision :: hib , hab - double precision :: delta_e_ab(N_states) - double precision :: hib_test,hja_test,hab_test if(ispin == kspin .and. vorb.le.rorb)cycle ! condition not to double count do jspin = 1, 2 if (jspin .ne. kspin)then @@ -472,18 +575,12 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin) enddo ! < idet | H | det_tmp > = phase * (ir|cv) -! call i_H_j(det_tmp,psi_det(1,1,idet),N_int,hib) call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) if(ispin == jspin)then hib= phase * (active_int(corb,1) - active_int(corb,2)) else hib= phase * active_int(corb,1) endif -! if(hib_test .ne. hib)then -! print*, 'hib_test .ne. hib' -! print*, hib, hib_test -! stop -! endif ! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} | det_tmp > call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok) @@ -491,21 +588,14 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb} -! call i_H_j(det_tmp_bis,det_tmp,N_int,hab) - hab = fock_operator_local(aorb,borb,kspin) * phase + hab = (fock_operator_local(aorb,borb,kspin) ) * phase ! < jdet | H | det_tmp_bis > = phase * (ir|cv) -! call i_H_j(det_tmp_bis,psi_det(1,1,idx(jdet)),N_int,hja) call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int) if(ispin == jspin)then hja= phase * (active_int(corb,1) - active_int(corb,2)) else hja= phase * (active_int(corb,1)) endif -! if(hja_test .ne. hja)then -! print*, 'hja_test .ne. hja' -! print*, hja, hja_test -! stop -! endif do istate = 1, N_states delta_e_ab(istate) = delta_e(corb,jspin,istate) + one_anhil_one_creat(borb,aorb,kspin,kspin,istate) matrix_1h2p(idx(jdet),idet,istate) = matrix_1h2p(idx(jdet),idet,istate) + & @@ -519,7 +609,7 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) do corb = 1, n_act_orb if(corb == aorb .or. corb == borb) cycle if(perturb_dets_phase(corb,jspin,ispin).le.-100d0)cycle - ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{corb,kspin} a_{iorb,ispin} | Idet > + ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{corb,jspin} a_{iorb,ispin} | Idet > do inint = 1, N_int det_tmp(inint,1) = perturb_dets(inint,1,corb,jspin,ispin) det_tmp(inint,2) = perturb_dets(inint,2,corb,jspin,ispin) @@ -527,38 +617,25 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin) enddo ! < idet | H | det_tmp > = phase * ( (ir|cv) - (iv|cr) ) -! call i_H_j(det_tmp,psi_det(1,1,idet),N_int,hib) call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) if(ispin == jspin)then hib= phase * (active_int(corb,1) - active_int(corb,2)) else hib= phase * active_int(corb,1) endif -! if(hib_test .ne. hib)then -! print*, 'hib_test .ne. hib jspin == kspin' -! print*, hib, hib_test -! stop -! endif ! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} | det_tmp > call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok) if(i_ok .ne. 1)cycle call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) ! ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb} hab = fock_operator_local(aorb,borb,kspin) * phase -! call i_H_j(det_tmp_bis,det_tmp,N_int,hab) ! < jdet | H | det_tmp_bis > = phase * ( (ir|cv) - (iv|cr) ) -! call i_H_j(det_tmp_bis,psi_det(1,1,idx(jdet)),N_int,hja) call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int) if(ispin == jspin)then hja= phase * (active_int(corb,1) - active_int(corb,2)) else hja= phase * (active_int(corb,1)) endif -! if(hja_test .ne. hja)then -! print*, 'hja_test .ne. hja' -! print*, hja, hja_test -! stop -! endif do istate = 1, N_states delta_e_ab(istate) = delta_e(corb,jspin,istate) + one_anhil_one_creat(borb,aorb,kspin,kspin,istate) matrix_1h2p(idx(jdet),idet,istate) = matrix_1h2p(idx(jdet),idet,istate) + & @@ -574,14 +651,114 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) enddo ! ispin else !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Case of double excitations !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! call debug_det(psi_det(1,1,idet),N_int) -! call debug_det(psi_det(1,1,idx(jdet)),N_int) -! pause + ! a^{\dagger}_r a_{i} (ispin) + aorb = index_orb_act_mono(idx(jdet),4) ! hole of a beta electron + borb = index_orb_act_mono(idx(jdet),5) ! propagation of the hole :: mono excitation of alpha spin + do ispin = 1, 2 ! you loop on all possible spin for the excitation + ! a^{\dagger}_r a_{i} (ispin) + ! ! first combination of spin :: | det_tmp > = a_{aorb,beta} | Idet > + jspin = 2 + if(perturb_dets_phase(aorb,jspin,ispin).le.-100d0)cycle + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin) + enddo + call get_double_excitation(det_tmp,psi_det(1,1,idet),exc,phase,N_int) + if(ispin == jspin)then + hib= phase * (active_int(borb,1) - active_int(borb,2)) + else + hib= phase * (active_int(borb,1)) + endif + if( index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),5))then + call do_mono_excitation(det_tmp_bis,list_act(borb),list_act(aorb),1,i_ok) + if(i_ok .ne. 1)then + call debug_det(psi_det(1,1,idet),N_int) + call debug_det(psi_det(1,1,idx(jdet)),N_int) + print*, aorb, borb + call debug_det(det_tmp,N_int) + stop + endif + else + call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),1,i_ok) + endif + + if(i_ok .ne. 1)cycle + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) + ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb} + if (aorb == borb)then + print*, 'iahaha' + stop + endif + hab = fock_operator_local(aorb,borb,1) * phase + call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int) + if(ispin == jspin)then + hja= phase * (active_int(borb,1) - active_int(borb,2)) + else + hja= phase * (active_int(borb,1)) + endif + do istate = 1, N_states + delta_e_ab(istate) = delta_e(aorb,jspin,istate) + one_anhil_one_creat(borb,aorb,1,1,istate) + matrix_1h2p(idx(jdet),idet,istate) = matrix_1h2p(idx(jdet),idet,istate) + & + hib / delta_e(aorb,jspin,istate) * hab / delta_e_ab(istate) * hja + ! < det_tmp | H | Idet > / delta_E (Idet --> det_tmp ) + ! < det_tmp | H | det_tmp_bis > / delta_E (Idet --> det_tmp --> det_tmp_bis) + ! < det_tmp_bis | H | Jdet > + enddo !! istate + + ! ! second combination of spin :: | det_tmp > = a_{aorb,alpha} | Idet > + jspin = 1 + if(perturb_dets_phase(aorb,jspin,ispin).le.-100d0)cycle + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,jspin,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,jspin,ispin) + enddo + call get_double_excitation(det_tmp,psi_det(1,1,idet),exc,phase,N_int) + if(ispin == jspin)then + hib= phase * (active_int(borb,1) - active_int(borb,2)) + else + hib= phase * (active_int(borb,1)) + endif + if( index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),5))then + call do_mono_excitation(det_tmp_bis,list_act(borb),list_act(aorb),2,i_ok) + if(i_ok .ne. 1)then + call debug_det(psi_det(1,1,idet),N_int) + call debug_det(psi_det(1,1,idx(jdet)),N_int) + print*, aorb, borb + call debug_det(det_tmp,N_int) + stop + endif + else + call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),2,i_ok) + endif + + if(i_ok .ne. 1)cycle + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) + ! < det_tmp | H | det_tmp_bis > = F_{aorb,borb} + hab = fock_operator_local(aorb,borb,2) * phase + call get_double_excitation(det_tmp_bis,psi_det(1,1,idx(jdet)),exc,phase,N_int) + if(ispin == jspin)then + hja= phase * (active_int(borb,1) - active_int(borb,2)) + else + hja= phase * (active_int(borb,1)) + endif + do istate = 1, N_states + delta_e_ab(istate) = delta_e(aorb,jspin,istate) + one_anhil_one_creat(borb,aorb,1,1,istate) + matrix_1h2p(idx(jdet),idet,istate) = matrix_1h2p(idx(jdet),idet,istate) + & + hib / delta_e(aorb,jspin,istate) * hab / delta_e_ab(istate) * hja + ! < det_tmp | H | Idet > / delta_E (Idet --> det_tmp ) + ! < det_tmp | H | det_tmp_bis > / delta_E (Idet --> det_tmp --> det_tmp_bis) + ! < det_tmp_bis | H | Jdet > + enddo !! istate + enddo !! ispin - endif + endif !! en of test if jdet is a single or a double excitation of type K_ab - else + else !! jdet is idet ! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations ! ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{aorb,kspin} a_{iorb,ispin} | Idet > @@ -603,14 +780,13 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) endif - enddo + enddo !! jdet enddo enddo enddo enddo - print* , 'accu_contrib = ',accu_contrib @@ -618,104 +794,3 @@ subroutine give_1h2p_contrib_sec_order(matrix_1h2p) end - - - - - - -! do a = 1, n_act_orb -! do jspin = 1, 2 -! do ispin = 1, 2 -! if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle -! ! determinant perturber | det_tmp > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{a,jspin} a_{i,ispin} | Idet > -! do inint = 1, N_int -! det_tmp(inint,1) = perturb_dets(inint,1,a,jspin,ispin) -! det_tmp(inint,2) = perturb_dets(inint,2,a,jspin,ispin) -! enddo -! do istate = 1, N_states -! coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0 -! enddo -! do b = 1, n_act_orb -! do kspin = jspin , jspin -! integer :: degree_scalar -! if( perturb_dets_phase(b,kspin,ispin) .le. -10.d0)cycle -! do inint = 1, N_int -! det_tmp_j(inint,1) = perturb_dets(inint,1,b,kspin,ispin) -! det_tmp_j(inint,2) = perturb_dets(inint,2,b,kspin,ispin) -! enddo -! call get_excitation_degree(det_tmp,det_tmp_j,degree_scalar,N_int) -! if (degree_scalar > 2 .or. degree_scalar == 0)cycle -! ! determinant perturber | det_tmp_j > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{b,jspin} a_{i,ispin} | Idet > -! call i_H_j(det_tmp,det_tmp_j,N_int,hij) -! do istate = 1, N_states -! coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(b,kspin,ispin,istate,1) & -! * hij / delta_e(a,jspin,istate) -! endif -! enddo -! enddo -! enddo -! enddo -! enddo -! enddo - - - -! do a = 1, n_act_orb -! do jspin = 1, 2 -! do ispin = 1, 2 -! if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle -! ! determinant perturber | det_tmp > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{a,jspin} a_{i,ispin} | Idet > -! do inint = 1, N_int -! det_tmp(inint,1) = iand(perturb_dets(inint,1,a,jspin,ispin),cas_bitmask(inint,1,1)) -! det_tmp(inint,2) = iand(perturb_dets(inint,2,a,jspin,ispin),cas_bitmask(inint,1,1)) -! enddo -! do istate = 1, N_states -! coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0 -! enddo -! do b = 1, n_act_orb -! do kspin = jspin , jspin -! integer :: degree_scalar -! if( perturb_dets_phase(b,kspin,ispin) .le. -10.d0)cycle -! do inint = 1, N_int -! det_tmp_j(inint,1) = iand(perturb_dets(inint,1,b,kspin,ispin),cas_bitmask(inint,1,1)) -! det_tmp_j(inint,2) = iand(perturb_dets(inint,2,b,kspin,ispin),cas_bitmask(inint,1,1)) -! enddo -! call get_excitation_degree(det_tmp,det_tmp_j,degree_scalar,N_int) -! if (degree_scalar > 2 .or. degree_scalar == 0)cycle -! ! determinant perturber | det_tmp_j > = a^{\dagger}_{r,ispin} a^{\dagger}_{v,jspin} a_{b,jspin} a_{i,ispin} | Idet > -!! print*, '**********************' -!! integer(bit_kind) :: det_bis(N_int,2) -!! call debug_det(det_tmp,N_int) -!! call debug_det(det_tmp_j,N_int) -!! do inint = 1, N_int -!! det_bis(inint,1) = perturb_dets(inint,1,b,kspin,ispin) -!! det_bis(inint,2) = perturb_dets(inint,2,b,kspin,ispin) -!! enddo -!! call debug_det(det_bis,N_int) -! call i_H_j_dyall(det_tmp,det_tmp_j,N_int,hij) -! do istate = 1, N_states -! coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(b,kspin,ispin,istate,1) & -! * hij / delta_e(a,jspin,istate) -! if(dabs(hij).gt.0.01d0)then -! print*,degree_scalar, hij -! print*, coef_perturb_from_idet(b,kspin,ispin,istate,1)* hij / delta_e(a,jspin,istate),coef_perturb_from_idet(a,jspin,ispin,istate,1) -! -! endif -! enddo -! enddo -! enddo -! enddo -! enddo -! enddo - -! do a = 1, n_act_orb -! do jspin = 1, 2 -! do ispin = 1, 2 -! if( perturb_dets_phase(a,jspin,ispin) .le. -10.d0)cycle -! do istate = 1, N_states -! coef_perturb_from_idet(a,jspin,ispin,istate,2) += coef_perturb_from_idet(a,jspin,ispin,istate,1) -! enddo -! enddo -! enddo -! enddo diff --git a/plugins/MRPT_Utils/print_1h2p.irp.f b/plugins/MRPT_Utils/print_1h2p.irp.f index 7d2d6c23..a03f2659 100644 --- a/plugins/MRPT_Utils/print_1h2p.irp.f +++ b/plugins/MRPT_Utils/print_1h2p.irp.f @@ -17,7 +17,8 @@ subroutine routine enddo enddo enddo - call give_1h2p_contrib_sec_order(matrix_1h2p) + if(.False.)then + call give_1h2p_contrib(matrix_1h2p) double precision :: accu accu = 0.d0 do i = 1, N_det @@ -25,7 +26,25 @@ subroutine routine accu += matrix_1h2p(i,j,1) * psi_coef(i,1) * psi_coef(j,1) enddo enddo - print*, 'accu', accu + print*, 'second order ', accu + endif + + if(.True.)then + do i = 1, N_det + do j = 1, N_det + do istate = 1, N_states + matrix_1h2p(i,j,istate) = 0.d0 + enddo + enddo + enddo + call give_1h2p_new(matrix_1h2p) + accu = 0.d0 + do i = 1, N_det + do j = 1, N_det + accu += matrix_1h2p(i,j,1) * psi_coef(i,1) * psi_coef(j,1) + enddo + enddo + print*, 'third order ', accu deallocate (matrix_1h2p) end diff --git a/plugins/MRPT_Utils/second_order_new.irp.f b/plugins/MRPT_Utils/second_order_new.irp.f new file mode 100644 index 00000000..7bfeeb9c --- /dev/null +++ b/plugins/MRPT_Utils/second_order_new.irp.f @@ -0,0 +1,303 @@ + +subroutine give_1h2p_new(matrix_1h2p) + use bitmasks + implicit none + double precision , intent(inout) :: matrix_1h2p(N_det,N_det,*) + integer :: i,v,r,a,b,c + integer :: iorb, vorb, rorb, aorb, borb,corb + integer :: ispin,jspin + integer :: idet,jdet + integer(bit_kind) :: perturb_dets(N_int,2,n_act_orb,2,2) + double precision :: perturb_dets_phase(n_act_orb,2,2) + double precision :: perturb_dets_hij(n_act_orb,2,2) + double precision :: perturb_dets_hpsi0(n_act_orb,2,2,N_states) + double precision :: coef_perturb_from_idet(n_act_orb,2,2,N_states,2) + logical :: already_generated(n_act_orb,2,2) + integer :: inint + integer :: elec_num_tab_local(2),acu_elec + integer(bit_kind) :: det_tmp(N_int,2) + integer(bit_kind) :: det_tmp_j(N_int,2) + integer :: exc(0:2,2,2) + integer :: accu_elec + double precision :: get_mo_bielec_integral_schwartz + double precision :: active_int(n_act_orb,2) + double precision :: hij,phase + double precision :: accu_contrib(N_states) + integer :: degree(N_det) + integer :: idx(0:N_det) + double precision :: delta_e(n_act_orb,2,N_states) + double precision :: delta_e_inv(n_act_orb,2,N_states) + integer :: istate + integer :: index_orb_act_mono(N_det,6) + double precision :: delta_e_inactive_virt(N_states) + integer :: kspin + double precision :: delta_e_ja(N_states) + double precision :: hja + double precision :: contrib_hij + double precision :: fock_operator_local(n_act_orb,n_act_orb,2) + double precision :: hij_test + integer ::i_ok + integer(bit_kind) :: det_tmp_bis(N_int,2) + double precision :: hib , hab + double precision :: delta_e_ab(N_states) + double precision :: hib_test,hja_test,hab_test + integer :: i_hole,i_part + double precision :: hia,hjb + integer :: other_spin(2) + other_spin(1) = 2 + other_spin(2) = 1 + + accu_contrib = 0.d0 +!matrix_1h2p = 0.d0 + + elec_num_tab_local = 0 + do inint = 1, N_int + elec_num_tab_local(1) += popcnt(psi_det(inint,1,1)) + elec_num_tab_local(2) += popcnt(psi_det(inint,2,1)) + enddo + do i = 1, n_inact_orb ! First inactive + iorb = list_inact(i) + do v = 1, n_virt_orb ! First virtual + vorb = list_virt(v) + do r = 1, n_virt_orb ! Second virtual + rorb = list_virt(r) + ! take all the integral you will need for i,j,r fixed + do a = 1, n_act_orb + aorb = list_act(a) + active_int(a,1) = get_mo_bielec_integral_schwartz(iorb,aorb,rorb,vorb,mo_integrals_map) ! direct + active_int(a,2) = get_mo_bielec_integral_schwartz(iorb,aorb,vorb,rorb,mo_integrals_map) ! exchange + perturb_dets_phase(a,1,1) = -1000.d0 + perturb_dets_phase(a,1,2) = -1000.d0 + perturb_dets_phase(a,2,2) = -1000.d0 + perturb_dets_phase(a,2,1) = -1000.d0 + enddo + + + do istate = 1, N_states + delta_e_inactive_virt(istate) = & + - fock_virt_total_spin_trace(rorb,istate) & + - fock_virt_total_spin_trace(vorb,istate) & + + fock_core_inactive_total_spin_trace(iorb,istate) + enddo + do idet = 1, N_det + call get_excitation_degree_vector_mono_or_exchange(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements + do ispin = 1, 2 ! spin of the couple a-a^dagger (iorb,rorb) + do jspin = 1, 2 ! spin of the couple a-a^dagger (aorb,vorb) + do a = 1, n_act_orb ! First active + aorb = list_act(a) + do istate = 1, N_states + perturb_dets_hpsi0(a,jspin,ispin,istate) = 0.d0 + coef_perturb_from_idet(a,jspin,ispin,istate,1) = 0.d0 + coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0 + enddo + if(ispin == jspin .and. vorb.le.rorb)cycle ! condition not to double count + do inint = 1, N_int + det_tmp(inint,1) = psi_det(inint,1,idet) + det_tmp(inint,2) = psi_det(inint,2,idet) + enddo + ! Do the excitation inactive -- > virtual + call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin + call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin + + ! Do the excitation active -- > virtual + call clear_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! hole in "aorb" of spin Jspin + call set_bit_to_integer(vorb,det_tmp(1,jspin),N_int) ! particle in "vorb" of spin Jspin + + ! Check if the excitation is possible or not on psi_det(idet) + accu_elec= 0 + do inint = 1, N_int + accu_elec+= popcnt(det_tmp(inint,jspin)) + enddo + if(accu_elec .ne. elec_num_tab_local(jspin))then + perturb_dets_phase(a,jspin,ispin) = -1000.0d0 + perturb_dets_hij(a,jspin,ispin) = 0.d0 + do istate = 1, N_states + coef_perturb_from_idet(a,jspin,ispin,istate,1) = 0.d0 + coef_perturb_from_idet(a,jspin,ispin,istate,2) = 0.d0 + enddo + cycle + endif + do inint = 1, N_int + perturb_dets(inint,1,a,jspin,ispin) = det_tmp(inint,1) + perturb_dets(inint,2,a,jspin,ispin) = det_tmp(inint,2) + enddo + + call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) + perturb_dets_phase(a,jspin,ispin) = phase + + do istate = 1, N_states + delta_e(a,jspin,istate) = one_anhil(a,jspin,istate) + delta_e_inactive_virt(istate) + delta_e_inv(a,jspin,istate) = 1.d0 / delta_e(a,jspin,istate) + enddo + if(ispin == jspin)then + perturb_dets_hij(a,jspin,ispin) = phase * (active_int(a,1) - active_int(a,2) ) + else + perturb_dets_hij(a,jspin,ispin) = phase * active_int(a,1) + endif + enddo + enddo + enddo + +!!!!!!!!!!!!!!!!!!!!!!!!!!! determination of the connections between I and the other J determinants mono excited in the CAS +!!!!!!!!!!!!!!!!!!!!!!!!!!!! the determinants I and J must be connected by the following operator +!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!! + do jdet = 1, idx(0) + if(idx(jdet).ne.idet)then + if(degree(jdet)==1)then + call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + if (exc(0,1,1) == 1) then + ! Mono alpha + i_hole = list_act_reverse(exc(1,1,1)) !!! a_a + i_part = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b} + kspin = 1 !!! kspin + index_orb_act_mono(idx(jdet),1) = i_hole + index_orb_act_mono(idx(jdet),2) = i_part + index_orb_act_mono(idx(jdet),3) = kspin + call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij) + fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator + fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator + else + ! Mono beta + i_hole = list_act_reverse(exc(1,1,2)) !!! a_a + i_part = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_{b} + kspin = 2 !!! kspin + index_orb_act_mono(idx(jdet),1) = i_hole + index_orb_act_mono(idx(jdet),2) = i_part + index_orb_act_mono(idx(jdet),3) = kspin + call i_H_j_dyall(psi_active(1,1,idet),psi_active(1,1,idx(jdet)),N_int,hij) + fock_operator_local(i_hole,i_part,kspin) = hij * phase ! phase less fock operator + fock_operator_local(i_part,i_hole,kspin) = hij * phase ! phase less fock operator + endif + endif + else + index_orb_act_mono(idx(jdet),1) = -1 + endif + enddo + + + + do jdet = 1, idx(0) +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CASE OF THE MONO EXCITATIONS + if(degree(jdet) == 1)then + ! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator + ! are connected by the presence of the perturbers determinants |det_tmp> + aorb = index_orb_act_mono(idx(jdet),1) ! a_{aorb} + borb = index_orb_act_mono(idx(jdet),2) ! a^{\dagger}_{borb} + kspin = index_orb_act_mono(idx(jdet),3) ! spin of the excitation + ! the determinants Idet and Jdet interact throw the following operator + ! | Jdet > = a^{\dagger}_{borb,kspin} a_{aorb, kspin} | Idet > + + accu_contrib = 0.d0 + do ispin = 1, 2 ! you loop on all possible spin for the excitation + ! a^{\dagger}_r a_{i} (ispin) + if(ispin == kspin .and. vorb.le.rorb)cycle ! condition not to double count + + ! FIRST ORDER CONTRIBUTION + + ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{aorb,kspin} a_{iorb,ispin} | Idet > + if(perturb_dets_phase(aorb,kspin,ispin) .le. -10.d0)cycle + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,aorb,kspin,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,aorb,kspin,ispin) + enddo + call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) + if(kspin == ispin)then + hia = phase * (active_int(aorb,1) - active_int(aorb,2) ) + else + hia = phase * active_int(aorb,1) + endif + call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int) + if(kspin == ispin)then + hja = phase * (active_int(borb,1) - active_int(borb,2) ) + else + hja = phase * active_int(borb,1) + endif + + call i_H_j(det_tmp,psi_det(1,1,idx(jdet)),N_int,hij_test) + if(hij_test .ne. hja)then + print*, 'hij_test .ne. hja' + stop + endif + contrib_hij = hja * hia + do istate = 1, N_states + accu_contrib(istate) += contrib_hij * delta_e_inv(aorb,kspin,istate) + enddo + + !!!! SECOND ORDER CONTRIBTIONS + ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,jspin} a_{corb,jspin} a_{iorb,ispin} | Idet > + do jspin = 1, 2 + do corb = 1, n_act_orb + if(perturb_dets_phase(corb,jspin,ispin) .le. -10.d0)cycle + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,corb,jspin,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,corb,jspin,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,corb,jspin,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,corb,jspin,ispin) + enddo + ! | det_tmp_bis > = a^{\dagger}_{borb,kspin} a_{aorb,kspin} a_{iorb,ispin} | Idet > + call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),kspin,i_ok) + if(i_ok .ne. 1)cycle + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) + hia = perturb_dets_hij(corb,jspin,ispin) + hab = fock_operator_local(aorb,borb,kspin) * phase + + call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int) + if(jspin == ispin)then + hjb = phase * (active_int(corb,1) - active_int(corb,2) ) + else + hjb = phase * active_int(corb,1) + endif + do istate = 1, N_states + accu_contrib(istate)+=hia * delta_e_inv(corb,jspin,istate) & ! | Idet > --> | det_tmp > + ! | det_tmp > --> | det_tmp_bis > + *hab / (delta_e(corb,jspin,istate) + one_anhil_one_creat(aorb,borb,kspin,kspin,istate)) & + *hjb + enddo + enddo + enddo + + + + enddo ! ispin + do istate = 1, N_states + matrix_1h2p(idet,idx(jdet),istate) += accu_contrib(istate) + enddo + + else if (degree(jdet) == 0)then + ! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations + ! + ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{aorb,kspin} a_{iorb,ispin} | Idet > + accu_contrib = 0.d0 + do ispin = 1, 2 + do kspin = 1, 2 + do a = 1, n_act_orb ! First active + if( perturb_dets_phase(a,kspin,ispin) .le. -10.d0)cycle + if(ispin == kspin .and. vorb.le.rorb)cycle ! condition not to double count + contrib_hij = perturb_dets_hij(a,kspin,ispin) * perturb_dets_hij(a,kspin,ispin) + do istate = 1, N_states + accu_contrib(istate) += contrib_hij * delta_e_inv(a,kspin,istate) + enddo + enddo + enddo + enddo + do istate = 1, N_states + matrix_1h2p(idet,idet,istate) += accu_contrib(istate) + enddo + + endif + enddo !! jdet + + + enddo + enddo + enddo + enddo + + + + + +end + diff --git a/src/Determinants/create_excitations.irp.f b/src/Determinants/create_excitations.irp.f index b2a78216..6af49681 100644 --- a/src/Determinants/create_excitations.irp.f +++ b/src/Determinants/create_excitations.irp.f @@ -31,8 +31,8 @@ subroutine do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2)) enddo if(n_elec_tmp .ne. elec_num)then - print*, n_elec_tmp,elec_num - call debug_det(key_in,N_int) + !print*, n_elec_tmp,elec_num + !call debug_det(key_in,N_int) i_ok = -1 endif end diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 01ac8e8d..d153d008 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1367,8 +1367,178 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s do i=1,sze d = popcnt(xor( key1(1,1,i), key2(1,1))) + & popcnt(xor( key1(1,2,i), key2(1,2))) - exchange_1 = popcnt(xor(iand(key1(1,1,i),key1(1,2,i)),iand(key2(1,1),key2(1,2)))) - exchange_2 = popcnt(iand(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2)))) + exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + if (d > 4)cycle + if (d ==4)then + if(exchange_1 .eq. 0 ) then + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + else + cycle + endif +! pause + else + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + endif + enddo + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + d = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,2),key2(1,2)))) + & + popcnt(xor(ior(key1(2,1,i),key1(2,2,i)),ior(key2(2,2),key2(2,2)))) + exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + & + popcnt(ior(xor(key1(2,1,i),key2(2,1)),xor(key1(2,2,i),key2(2,2)))) + if (d > 4)cycle + if (d ==4)then + if(exchange_1 .eq. 0 ) then + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + else + cycle + endif +! pause + else + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + d = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2)))) + & + popcnt(xor(ior(key1(2,1,i),key1(2,2,i)),ior(key2(2,1),key2(2,2)))) + & + popcnt(xor(ior(key1(3,1,i),key1(3,2,i)),ior(key2(3,1),key2(3,2)))) + exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + & + popcnt(ior(xor(key1(2,1,i),key2(2,1)),xor(key1(2,2,i),key2(2,2)))) + & + popcnt(ior(xor(key1(3,1,i),key2(3,1)),xor(key1(3,2,i),key2(3,2)))) + if (d > 4)cycle + if (d ==4)then + if(exchange_1 .eq. 0 ) then + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + else + cycle + endif +! pause + else + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + d = 0 + exchange_1 = 0 + !DIR$ LOOP COUNT MIN(4) + do m=1,Nint + d = d + popcnt(xor( key1(m,1,i), key2(m,1))) & + + popcnt(xor( key1(m,2,i), key2(m,2))) + exchange_1 = popcnt(xor(ior(key1(m,1,i),key1(m,2,i)),ior(key2(m,1),key2(m,2)))) + exchange_2 = popcnt(ior(xor(key1(m,1,i),key2(m,1)),xor(key1(m,2,i),key2(m,2)))) + enddo + if (d > 4)cycle + if (d ==4)then + if(exchange_1 .eq. 0 ) then + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + else if (exchange_1 .eq. 2 .and. exchange_2.eq.2)then + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + else + cycle + endif +! pause + else + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + endif + enddo + + endif + idx(0) = l-1 +end + +subroutine get_excitation_degree_vector_mono_or_exchange_verbose(key1,key2,degree,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Applies get_excitation_degree to an array of determinants and return only the mono excitations + ! and the connections through exchange integrals + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: degree(sze) + integer, intent(out) :: idx(0:sze) + + integer :: i,l,d,m + integer :: exchange_1,exchange_2 + + ASSERT (Nint > 0) + ASSERT (sze > 0) + + l=1 + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + d = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + exchange_1 = popcnt(xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2)))) + exchange_2 = popcnt(ior(xor(key1(1,1,i),key2(1,1)),xor(key1(1,2,i),key2(1,2)))) + if(i==99)then + integer(bit_kind) :: key_test(N_int,2) + key_test(1,2) = 0_bit_kind + call debug_det(key2,N_int) + key_test(1,1) = ior(key2(1,1),key2(1,2)) + call debug_det(key_test,N_int) + key_test(1,1) = ior(key1(1,1,i),key1(1,2,i)) + call debug_det(key1(1,1,i),N_int) + call debug_det(key_test,N_int) + key_test(1,1) = xor(ior(key1(1,1,i),key1(1,2,i)),ior(key2(1,1),key2(1,2))) + call debug_det(key_test,N_int) + print*, exchange_1 , exchange_2 + stop + endif if (d > 4)cycle if (d ==4)then if(exchange_1 .eq. 0 ) then