From 376e4940db043ea19ed2c040552f387f08a06a52 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 19 Sep 2016 18:06:34 +0200 Subject: [PATCH] second order works for 2p --- plugins/MRPT_Utils/second_order_new.irp.f | 13 +- plugins/MRPT_Utils/second_order_new_2p.irp.f | 283 +++++++++++++++++++ src/Determinants/slater_rules.irp.f | 2 + 3 files changed, 292 insertions(+), 6 deletions(-) create mode 100644 plugins/MRPT_Utils/second_order_new_2p.irp.f diff --git a/plugins/MRPT_Utils/second_order_new.irp.f b/plugins/MRPT_Utils/second_order_new.irp.f index 54b284f6..bcd08bf5 100644 --- a/plugins/MRPT_Utils/second_order_new.irp.f +++ b/plugins/MRPT_Utils/second_order_new.irp.f @@ -198,7 +198,8 @@ subroutine give_1h2p_new(matrix_1h2p) if(ispin == kspin .and. vorb.le.rorb)then cycle_same_spin_first_order = .True. endif - if(ispin .ne. kspin .and. cycle_same_spin_first_order == .False. )then ! condition not to double count +! if(ispin .ne. kspin .and. cycle_same_spin_first_order == .False. )then ! condition not to double count + if(cycle_same_spin_first_order == .False. )then ! condition not to double count ! FIRST ORDER CONTRIBUTION @@ -234,7 +235,7 @@ subroutine give_1h2p_new(matrix_1h2p) 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 + if(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 @@ -291,7 +292,7 @@ subroutine give_1h2p_new(matrix_1h2p) 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(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) @@ -325,7 +326,7 @@ subroutine give_1h2p_new(matrix_1h2p) 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(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) @@ -364,7 +365,7 @@ subroutine give_1h2p_new(matrix_1h2p) 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(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) @@ -399,7 +400,7 @@ subroutine give_1h2p_new(matrix_1h2p) 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(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) diff --git a/plugins/MRPT_Utils/second_order_new_2p.irp.f b/plugins/MRPT_Utils/second_order_new_2p.irp.f new file mode 100644 index 00000000..2e94527c --- /dev/null +++ b/plugins/MRPT_Utils/second_order_new_2p.irp.f @@ -0,0 +1,283 @@ + +subroutine give_2p_new(matrix_2p) + use bitmasks + implicit none + double precision , intent(inout) :: matrix_2p(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,n_act_orb,2,2) + double precision :: perturb_dets_phase(n_act_orb,n_act_orb,2,2) + double precision :: perturb_dets_hij(n_act_orb,n_act_orb,2,2) + double precision :: perturb_dets_hpsi0(n_act_orb,n_act_orb,2,2,N_states) + 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,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,n_act_orb,2,2,N_states) + double precision :: delta_e_inv(n_act_orb,n_act_orb,2,2,N_states) + double precision :: delta_e_inactive_virt(N_states) + integer :: istate + integer :: index_orb_act_mono(N_det,6) + 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_2p = 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 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 r,v fixed + do a = 1, n_act_orb + aorb = list_act(a) + do b = 1, n_act_orb + borb = list_act(b) + active_int(a,b,1) = get_mo_bielec_integral_schwartz(aorb,borb,rorb,vorb,mo_integrals_map) ! direct ( a--> r | b--> v ) + active_int(a,b,2) = get_mo_bielec_integral_schwartz(aorb,borb,vorb,rorb,mo_integrals_map) ! exchange ( b--> r | a--> v ) + perturb_dets_phase(a,b,1,1) = -1000.d0 + perturb_dets_phase(a,b,1,2) = -1000.d0 + perturb_dets_phase(a,b,2,2) = -1000.d0 + perturb_dets_phase(a,b,2,1) = -1000.d0 + perturb_dets_phase(b,a,1,1) = -1000.d0 + perturb_dets_phase(b,a,1,2) = -1000.d0 + perturb_dets_phase(b,a,2,2) = -1000.d0 + perturb_dets_phase(b,a,2,1) = -1000.d0 + enddo + 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) + enddo + 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(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 (aorb,rorb) + do jspin = 1, 2 ! spin of the couple a-a^dagger (borb,vorb) + do b = 1, n_act_orb ! First active + borb = list_act(b) + do a = 1, n_act_orb ! First active + aorb = list_act(a) +! if(ispin == 2.and. jspin ==1)then +! perturb_dets_phase(a,b,ispin,jspin) = -1000.0d0 +! perturb_dets_hij(a,b,ispin,jspin) = 0.d0 +! cycle ! condition not to double count +! endif + + if(ispin == jspin .and. vorb.le.rorb)then + perturb_dets_phase(a,b,ispin,jspin) = -1000.0d0 + perturb_dets_hij(a,b,ispin,jspin) = 0.d0 + cycle ! condition not to double count + endif + if(ispin == jspin .and. aorb.le.borb) then + perturb_dets_phase(a,b,ispin,jspin) = -1000.0d0 + perturb_dets_hij(a,b,ispin,jspin) = 0.d0 + cycle ! condition not to double count + endif + 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 (aorb,ispin) --> (rorb,ispin) + call clear_bit_to_integer(aorb,det_tmp(1,ispin),N_int) ! hole in "aorb" of spin Ispin + call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin + + ! Do the excitation (borb,jspin) --> (vorb,jspin) + call clear_bit_to_integer(borb,det_tmp(1,jspin),N_int) ! hole in "borb" 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,1)) + popcnt(det_tmp(inint,2)) + enddo + if(accu_elec .ne. elec_num_tab_local(2)+elec_num_tab_local(1))then + perturb_dets_phase(a,b,ispin,jspin) = -1000.0d0 + perturb_dets_hij(a,b,ispin,jspin) = 0.d0 + cycle + endif + do inint = 1, N_int + perturb_dets(inint,1,a,b,ispin,jspin) = det_tmp(inint,1) + perturb_dets(inint,2,a,b,ispin,jspin) = det_tmp(inint,2) + enddo + + call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) + perturb_dets_phase(a,b,ispin,jspin) = phase + + do istate = 1, N_states + delta_e(a,b,ispin,jspin,istate) = two_anhil(a,b,ispin,jspin,istate) + delta_e_inactive_virt(istate) + delta_e_inv(a,b,ispin,jspin,istate) = 1.d0 / delta_e(a,b,ispin,jspin,istate) + enddo + if(ispin == jspin)then + perturb_dets_hij(a,b,ispin,jspin) = phase * (active_int(a,b,2) - active_int(a,b,1) ) + else + perturb_dets_hij(a,b,ispin,jspin) = phase * active_int(a,b,1) + endif + call i_H_j(psi_det(1,1,idet),det_tmp,N_int,hij) + if(hij.ne.perturb_dets_hij(a,b,ispin,jspin))then + print*, active_int(a,b,1) , active_int(b,a,1) + double precision :: hmono,hdouble + call i_H_j_verbose(psi_det(1,1,idet),det_tmp,N_int,hij,hmono,hdouble) + print*, 'pb !! hij.ne.perturb_dets_hij(a,b,ispin,jspin)' + print*, ispin,jspin + print*, aorb,rorb,borb,vorb + print*, hij,perturb_dets_hij(a,b,ispin,jspin) + call debug_det(psi_det(1,1,idet),N_int) + call debug_det(det_tmp,N_int) + stop + endif + enddo ! b + enddo ! a + enddo ! jspin + enddo ! ispin + +!!!!!!!!!!!!!!!!!!!!!!!!!!! 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(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) + if (exc(0,1,1) == 1) then + ! 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 + else if (exc(0,1,1) == 2) then + 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 + index_orb_act_mono(idx(jdet),4) = list_act_reverse(exc(2,1,1)) !!! a_c ALPHA + index_orb_act_mono(idx(jdet),5) = list_act_reverse(exc(2,2,1)) !!! a^{\dagger}_{d} ALPHA + index_orb_act_mono(idx(jdet),6) = 1 + else if (exc(0,1,2) == 2) then + index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,2)) !!! a_a BETA + index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(2,1,2)) !!! a^{\dagger}_{b} BETA + index_orb_act_mono(idx(jdet),3) = 2 + index_orb_act_mono(idx(jdet),4) = list_act_reverse(exc(1,2,2)) !!! a_c BETA + index_orb_act_mono(idx(jdet),5) = list_act_reverse(exc(2,2,2)) !!! a^{\dagger}_{d} BETA + index_orb_act_mono(idx(jdet),6) = 2 + endif + 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_{a} (ispin) + !!!! SECOND ORDER CONTRIBTIONS + ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,jspin} a_{corb,jspin} a_{iorb,ispin} | Idet > + do jspin = 1, 2 + if(ispin == 2 .and. jspin ==1)cycle + do b = 1, n_act_orb + do a = 1, n_act_orb + logical :: cycle_same_spin_second_order(2) + cycle_same_spin_second_order(1) = .False. + cycle_same_spin_second_order(2) = .False. + if(perturb_dets_phase(a,b,ispin,jspin).le.-10d0)cycle + if(ispin == jspin .and. vorb.le.rorb)then + cycle_same_spin_second_order(1) = .True. + endif + if(ispin == jspin .and. aorb.le.borb)then + cycle_same_spin_second_order(2) = .True. + endif + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,a,b,ispin,jspin) + det_tmp(inint,2) = perturb_dets(inint,2,a,b,ispin,jspin) + enddo + do jdet = 1, idx(0) +! if(idx(jdet).gt.idet)cycle + do istate = 1, N_states + call i_H_j(psi_det(1,1,idx(jdet)),det_tmp,N_int,hij) + matrix_2p(idx(jdet),idet,istate) += hij * perturb_dets_hij(a,b,ispin,jspin) * delta_e_inv(a,b,ispin,jspin,istate) + enddo + enddo ! jdet + enddo ! b + enddo ! a + enddo ! jspin + enddo ! ispin + +! else if (degree(jdet) == 0)then +! +! endif +! enddo !! jdet + + + enddo + enddo + enddo +end diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index b1666864..cb96957a 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -747,6 +747,8 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) exc(1,1,2), & exc(1,2,1), & exc(1,2,2) ,mo_integrals_map) + print*, 'hij verbose ',hij * phase + print*, 'phase verbose',phase else if (exc(0,1,1) == 2) then ! Double alpha print*,'phase hij = ',phase