From c6b7acbc4e68dc22acaa5e5198c96314392d2311 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 19 Sep 2016 13:38:37 +0200 Subject: [PATCH] definitive version for the second order of 1h2p and 2h1p --- plugins/MRPT_Utils/second_order_new.irp.f | 419 ++++++++++++++++------ src/Determinants/slater_rules.irp.f | 12 +- 2 files changed, 321 insertions(+), 110 deletions(-) diff --git a/plugins/MRPT_Utils/second_order_new.irp.f b/plugins/MRPT_Utils/second_order_new.irp.f index 57439580..54b284f6 100644 --- a/plugins/MRPT_Utils/second_order_new.irp.f +++ b/plugins/MRPT_Utils/second_order_new.irp.f @@ -137,7 +137,6 @@ subroutine give_1h2p_new(matrix_1h2p) !!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 @@ -163,10 +162,17 @@ subroutine give_1h2p_new(matrix_1h2p) 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 ALPHA + index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,2,1)) !!! a^{\dagger}_{b} ALPHA + 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 BETA + index_orb_act_mono(idx(jdet),5) = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_{b} BETA + index_orb_act_mono(idx(jdet),6) = 2 endif - else - index_orb_act_mono(idx(jdet),1) = -1 - endif enddo @@ -185,65 +191,83 @@ subroutine give_1h2p_new(matrix_1h2p) 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) + +! if(ispin == kspin .and. vorb.le.rorb)cycle ! condition not to double count + logical :: cycle_same_spin_first_order + cycle_same_spin_first_order = .False. + if(ispin == kspin .and. vorb.le.rorb)then + cycle_same_spin_first_order = .True. 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) + if(ispin .ne. kspin .and. cycle_same_spin_first_order == .False. )then ! 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 + + contrib_hij = hja * hia + do istate = 1, N_states + accu_contrib(istate) += contrib_hij * delta_e_inv(aorb,kspin,istate) + enddo 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 + logical :: cycle_same_spin_second_order + cycle_same_spin_second_order = .False. + if(ispin == jspin .and. vorb.le.rorb)then + cycle_same_spin_second_order = .True. + endif + if(ispin .ne. jspin .or. cycle_same_spin_second_order == .False.)then + 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 + if(dabs(hia).le.1.d-12)cycle + if(dabs(hab).le.1.d-12)cycle + + 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 + if(dabs(hjb).le.1.d-12)cycle + 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 + endif enddo @@ -253,6 +277,176 @@ subroutine give_1h2p_new(matrix_1h2p) matrix_1h2p(idet,idx(jdet),istate) += accu_contrib(istate) enddo + else if (degree(jdet) == 2)then + ! CASE OF THE DOUBLE EXCITATIONS, ONLY THIRD ORDER EFFECTS + accu_contrib = 0.d0 + do ispin = 1, 2 ! you loop on all possible spin for the excitation + ! a^{\dagger}_r a_{i} (ispin) + ! if it is standard exchange case, the hole ALPHA == the part. BETA + if (index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),5))then + aorb = index_orb_act_mono(idx(jdet),1) !! the HOLE of the ALPHA electron + borb = index_orb_act_mono(idx(jdet),4) !! the HOLE of the BETA electron + ! first case :: | det_tmp > == a_{borb,\beta} | Idet > + cycle_same_spin_second_order = .False. + if(ispin == 2 .and. vorb.le.rorb)then + cycle_same_spin_second_order = .True. + endif + if(ispin .ne. 2 .or. cycle_same_spin_second_order == .False.)then ! condition not to double count + if(perturb_dets_phase(borb,2,ispin) .le. -10.d0)cycle + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,borb,2,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,borb,2,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,borb,2,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,borb,2,ispin) + enddo + hia = perturb_dets_hij(borb,2,ispin) + if(dabs(hia).le.1.d-12)cycle + call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),1,i_ok) + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) + hab = fock_operator_local(aorb,borb,1) * phase + + if(dabs(hab).le.1.d-12)cycle + call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int) + if(ispin == 2)then + hjb = phase * (active_int(aorb,1) - active_int(aorb,2) ) + else if (ispin == 1)then + hjb = phase * active_int(aorb,1) + endif + if(dabs(hjb).le.1.d-12)cycle + do istate = 1, N_states + accu_contrib(istate) += hia * delta_e_inv(borb,2,istate) & ! | Idet > --> | det_tmp > + ! | det_tmp > --> | det_tmp_bis > + * hab / (delta_e(borb,2,istate) + one_anhil_one_creat(aorb,borb,1,1,istate)) & + * hjb + enddo + endif + ! second case :: | det_tmp > == a_{aorb,\alpha} | Idet > + cycle_same_spin_second_order = .False. + if(ispin == 1 .and. vorb.le.rorb)then + cycle_same_spin_second_order = .True. + endif + if(ispin .ne. 1 .or. cycle_same_spin_second_order == .False.)then ! condition not to double count + if(perturb_dets_phase(aorb,1,ispin) .le. -10.d0)cycle + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,aorb,1,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,aorb,1,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,1,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,1,ispin) + enddo + hia = perturb_dets_hij(aorb,1,ispin) + if(dabs(hia).le.1.d-12)cycle + call do_mono_excitation(det_tmp_bis,list_act(borb),list_act(aorb),2,i_ok) + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) + hab = fock_operator_local(aorb,borb,2) * phase + + if(dabs(hab).le.1.d-12)cycle + call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int) + if(ispin == 1)then + hjb = phase * (active_int(borb,1) - active_int(borb,2) ) + else if (ispin == 2)then + hjb = phase * active_int(borb,1) + endif + if(dabs(hjb).le.1.d-12)cycle + do istate = 1, N_states + accu_contrib(istate) += hia * delta_e_inv(aorb,1,istate) & ! | Idet > --> | det_tmp > + ! | det_tmp > --> | det_tmp_bis > + * hab / (delta_e(aorb,1,istate) + one_anhil_one_creat(borb,aorb,2,2,istate)) & + * hjb + enddo + endif + + ! if it is a closed shell double excitation, the hole ALPHA == the hole BETA + else if (index_orb_act_mono(idx(jdet),1) == index_orb_act_mono(idx(jdet),4))then + aorb = index_orb_act_mono(idx(jdet),1) !! the HOLE of the ALPHA electron + borb = index_orb_act_mono(idx(jdet),2) !! the PART of the ALPHA electron + ! first case :: | det_tmp > == a_{aorb,\beta} | Idet > + cycle_same_spin_second_order = .False. + if(ispin == 2 .and. vorb.le.rorb)then + cycle_same_spin_second_order = .True. + endif + if(ispin .ne. 2 .or. cycle_same_spin_second_order == .False.)then ! condition not to double count + if(perturb_dets_phase(aorb,2,ispin) .le. -10.d0)cycle + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,aorb,2,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,aorb,2,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,2,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,2,ispin) + enddo + hia = perturb_dets_hij(aorb,2,ispin) + if(dabs(hia).le.1.d-12)cycle + call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),1,i_ok) + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) + hab = fock_operator_local(aorb,borb,1) * phase + + if(dabs(hab).le.1.d-12)cycle + call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int) + if(ispin == 2)then + hjb = phase * (active_int(borb,1) - active_int(borb,2) ) + else if (ispin == 1)then + hjb = phase * active_int(borb,1) + endif + if(dabs(hjb).le.1.d-12)cycle + do istate = 1, N_states + accu_contrib(istate) += hia * delta_e_inv(aorb,2,istate) & ! | Idet > --> | det_tmp > + ! | det_tmp > --> | det_tmp_bis > + * hab / (delta_e(aorb,2,istate) + one_anhil_one_creat(aorb,borb,1,1,istate)) & + * hjb + enddo + endif + + ! second case :: | det_tmp > == a_{aorb,\alpha} | Idet > + cycle_same_spin_second_order = .False. + if(ispin == 1 .and. vorb.le.rorb)then + cycle_same_spin_second_order = .True. + endif + if(ispin .ne. 1 .or. cycle_same_spin_second_order == .False.)then ! condition not to double count + if(perturb_dets_phase(aorb,1,ispin) .le. -10.d0)cycle + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,aorb,1,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,aorb,1,ispin) + det_tmp_bis(inint,1) = perturb_dets(inint,1,aorb,1,ispin) + det_tmp_bis(inint,2) = perturb_dets(inint,2,aorb,1,ispin) + enddo + hia = perturb_dets_hij(aorb,1,ispin) + if(dabs(hia).le.1.d-12)cycle + call do_mono_excitation(det_tmp_bis,list_act(aorb),list_act(borb),2,i_ok) + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) + hab = fock_operator_local(aorb,borb,2) * phase + + if(dabs(hab).le.1.d-12)cycle + call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp_bis,exc,phase,N_int) + if(ispin == 1)then + hjb = phase * (active_int(borb,1) - active_int(borb,2) ) + else if (ispin == 2)then + hjb = phase * active_int(borb,1) + endif + if(dabs(hjb).le.1.d-12)cycle + do istate = 1, N_states + accu_contrib(istate) += hia * delta_e_inv(aorb,1,istate) & ! | Idet > --> | det_tmp > + ! | det_tmp > --> | det_tmp_bis > + * hab / (delta_e(aorb,1,istate) + one_anhil_one_creat(aorb,borb,2,2,istate)) & + * hjb + enddo + endif + + + else + ! one should not fall in this case ... + call debug_det(psi_det(1,1,i),N_int) + call debug_det(psi_det(1,1,idx(jdet)),N_int) + call get_double_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + call decode_exc(exc,2,h1,p1,h2,p2,s1,s2) + integer :: h1, p1, h2, p2, s1, s2 + print*, h1, p1, h2, p2, s1, s2 + + print*, 'pb !!! it is a double but not an exchange case ....' + stop + endif + 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 ! @@ -447,59 +641,77 @@ subroutine give_2h1p_new(matrix_2h1p) 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. iorb.le.jorb)cycle ! condition not to double count - - ! | 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 - ! 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,1) - active_int(borb,2) ) - else - hja = phase * active_int(borb,1) +! if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count + logical :: cycle_same_spin_first_order + cycle_same_spin_first_order = .False. + if(ispin == kspin .and. iorb.le.jorb)then + cycle_same_spin_first_order = .True. endif - + if(ispin .ne. kspin .or. cycle_same_spin_first_order == .False. )then! condition not to double count - do istate = 1, N_states - accu_contrib(istate) += hja * perturb_dets_hij(aorb,kspin,ispin) * delta_e_inv(aorb,kspin,istate) - enddo + ! | 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 + ! 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,1) - active_int(borb,2) ) + else + hja = phase * active_int(borb,1) + endif +!! if(dabs(hja).le.1.d-10)cycle + + + do istate = 1, N_states + accu_contrib(istate) += hja * perturb_dets_hij(aorb,kspin,ispin) * delta_e_inv(aorb,kspin,istate) + enddo + endif + logical :: cycle_same_spin_second_order !!!! SECOND ORDER CONTRIBUTIONS !!!! SECOND ORDER CONTRIBTIONS ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{corb,jspin} a_{jorb,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}_{aorb,kspin} a_{borb,kspin} a_{iorb,kspin} | Idet > - call do_mono_excitation(det_tmp_bis,list_act(borb),list_act(aorb),kspin,i_ok) - if(i_ok .ne. 1)cycle - hia = perturb_dets_hij(corb,jspin,ispin) - call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) - hab = fock_operator_local(borb,aorb,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(borb,aorb,kspin,kspin,istate)) & - *hjb - enddo - enddo + cycle_same_spin_second_order = .False. + if(ispin == jspin .and. iorb.le.jorb)then + cycle_same_spin_second_order = .True. + endif + if(ispin .ne. jspin .or. cycle_same_spin_second_order == .False. )then! condition not to double count + 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}_{aorb,kspin} a_{borb,kspin} a_{iorb,kspin} | Idet > + call do_mono_excitation(det_tmp_bis,list_act(borb),list_act(aorb),kspin,i_ok) + if(i_ok .ne. 1)cycle + hia = perturb_dets_hij(corb,jspin,ispin) + if(dabs(hia).le.1.d-10)cycle + call get_mono_excitation(det_tmp,det_tmp_bis,exc,phase,N_int) + hab = fock_operator_local(borb,aorb,kspin) * phase + if(dabs(hab).le.1.d-10)cycle + + 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 + if(dabs(hjb).le.1.d-10)cycle + 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(borb,aorb,kspin,kspin,istate)) & + *hjb + enddo + enddo ! jspin + endif enddo enddo ! ispin do istate = 1, N_states @@ -516,6 +728,7 @@ subroutine give_2h1p_new(matrix_2h1p) if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count do a = 1, n_act_orb ! First active contrib_hij = perturb_dets_hij(a,kspin,ispin) * perturb_dets_hij(a,kspin,ispin) + if(dabs(contrib_hij).le.1.d-10)cycle do istate = 1, N_states accu_contrib(istate) += contrib_hij * delta_e_inv(a,kspin,istate) enddo diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index d153d008..b1666864 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1353,6 +1353,7 @@ subroutine get_excitation_degree_vector_mono_or_exchange(key1,key2,degree,Nint,s integer(bit_kind), intent(in) :: key2(Nint,2) integer, intent(out) :: degree(sze) integer, intent(out) :: idx(0:sze) + integer(bit_kind) :: key_tmp(Nint,2) integer :: i,l,d,m integer :: exchange_1,exchange_2 @@ -1367,15 +1368,12 @@ 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(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)))) + key_tmp(1,1) = xor(key1(1,1,i),key2(1,1)) + key_tmp(1,2) = xor(key1(1,2,i),key2(1,2)) + if(popcnt(key_tmp(1,1)) .gt.3 .or. popcnt(key_tmp(1,2)) .gt.3 )cycle !! no double excitations of same spin 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 + if(popcnt(xor(key_tmp(1,1),key_tmp(1,2))) == 0)then degree(l) = ishft(d,-1) idx(l) = i l = l+1