diff --git a/plugins/MRPT_Utils/new_way.irp.f b/plugins/MRPT_Utils/new_way.irp.f index 83831e41..8cf6ec5f 100644 --- a/plugins/MRPT_Utils/new_way.irp.f +++ b/plugins/MRPT_Utils/new_way.irp.f @@ -112,8 +112,8 @@ subroutine give_2h1p_contrib(matrix_2h1p) 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),2) = list_act_reverse(exc(1,1,1)) !!! a^{\dagger}_a - index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,1)) !!! a_{b} + 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 @@ -189,4 +189,205 @@ subroutine give_2h1p_contrib(matrix_2h1p) +end + + +subroutine give_1h2p_contrib(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 :: 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 :: coef_perturb_from_idet(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 :: 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 +!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 + enddo + + integer :: degree(N_det) + 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) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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) + 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) = 0.0 + perturb_dets_hij(a,jspin,ispin) = 0.d0 + do istate = 1, N_states + coef_perturb_from_idet(a,jspin,ispin,istate) = 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 + 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 + + 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) & + - fock_virt_total_spin_trace(rorb,istate) & + - fock_virt_total_spin_trace(vorb,istate) & + + fock_core_inactive_total_spin_trace(iorb,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 +!!!!!!!!!!!!!!!!!!!!!1 Computation of the coefficient at first order coming from idet +!!!!!!!!!!!!!!!!!!!!! for the excitation (i,j)(ispin,jspin) ---> (r,a)(ispin,jspin) + do istate = 1, N_states + coef_perturb_from_idet(a,jspin,ispin,istate) = perturb_dets_hij(a,jspin,ispin) / delta_e(a,jspin,istate) + enddo + + 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 + 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,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 + else + ! Mono beta + index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,1,2)) !!! a_a + index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_{b} + index_orb_act_mono(idx(jdet),3) = 2 + endif + else + index_orb_act_mono(idx(jdet),1) = -1 + endif + enddo + + integer :: kspin + 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_{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 > + + 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 + + ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{vorb,kspin} a_{aorb,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}_{vorb,kspin} a_{borb,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 + + do istate = 1, N_states + matrix_1h2p(idx(jdet),idet,istate) += hja * coef_perturb_from_idet(aorb,kspin,ispin,istate) + 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}_{vorb,kspin} a_{aorb,kspin} a_{iorb,ispin} | Idet > + do ispin = 1, 2 + do kspin = 1, 2 + do a = 1, n_act_orb ! First active + aorb = list_act(a) + if(ispin == kspin .and. vorb.le.rorb)cycle ! condition not to double count + do istate = 1, N_states + matrix_1h2p(idet,idet,istate) += coef_perturb_from_idet(a,kspin,ispin,istate) * perturb_dets_hij(a,kspin,ispin) + enddo + enddo + enddo + enddo + + endif + + enddo + enddo + enddo + enddo + enddo + + + + + end diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index da13ec44..5a60d093 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -8,14 +8,14 @@ BEGIN_PROVIDER [integer(bit_kind), psi_active, (N_int,2,psi_det_size)] use bitmasks integer :: i,j,k,l provide cas_bitmask - print*, 'psi_active ' +!print*, 'psi_active ' do i = 1, N_det do j = 1, N_int psi_active(j,1,i) = iand(psi_det(j,1,i),cas_bitmask(j,1,1)) psi_active(j,2,i) = iand(psi_det(j,2,i),cas_bitmask(j,1,1)) enddo - call debug_det(psi_active(1,1,i),N_int) +! call debug_det(psi_active(1,1,i),N_int) enddo END_PROVIDER diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 0357fb88..cb4e12b4 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -773,6 +773,18 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) else if (exc(0,1,2) == 2) then ! Double beta + print*,'phase hij = ',phase + print*, get_mo_bielec_integral_schwartz( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map ) + print*, get_mo_bielec_integral_schwartz( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) + hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,2), & exc(2,1,2), &