10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-23 12:56:14 +01:00
quantum_package/plugins/FOBOCI/corr_energy_2h2p.irp.f
2016-03-14 16:01:55 +01:00

426 lines
15 KiB
Fortran

BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_ab, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h2p_ab_2_orb, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h2p_bb_2_orb, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_a, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_b, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h2p_for_1h1p_double, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_aa, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h2p_per_orb_bb, (mo_tot_num)]
&BEGIN_PROVIDER [ double precision, total_corr_e_2h2p]
use bitmasks
print*,''
print*,'Providing the 2h2p correlation energy'
print*,''
implicit none
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,k,l
integer :: i_hole,j_hole,k_part,l_part
double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
double precision :: diag_H_mat_elem
integer :: i_ok,ispin
! Alpha - Beta correlation energy
total_corr_e_2h2p = 0.d0
corr_energy_2h2p_ab_2_orb = 0.d0
corr_energy_2h2p_bb_2_orb = 0.d0
corr_energy_2h2p_per_orb_ab = 0.d0
corr_energy_2h2p_per_orb_aa = 0.d0
corr_energy_2h2p_per_orb_bb = 0.d0
corr_energy_2h2p_for_1h1p_a = 0.d0
corr_energy_2h2p_for_1h1p_b = 0.d0
corr_energy_2h2p_for_1h1p_double = 0.d0
do i = 1, n_inact_orb ! beta
i_hole = list_inact(i)
do k = 1, n_virt_orb ! beta
k_part = list_virt(k)
do j = 1, n_inact_orb ! alpha
j_hole = list_inact(j)
do l = 1, n_virt_orb ! alpha
l_part = list_virt(l)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = (ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
contrib = hij*hij/delta_e
total_corr_e_2h2p += contrib
! Single orbital contribution
corr_energy_2h2p_per_orb_ab(i_hole) += contrib
corr_energy_2h2p_per_orb_ab(k_part) += contrib
! Couple of orbital contribution for the single 1h1p
corr_energy_2h2p_for_1h1p_a(j_hole,l_part) += contrib
corr_energy_2h2p_for_1h1p_a(l_part,j_hole) += contrib
corr_energy_2h2p_for_1h1p_b(j_hole,l_part) += contrib
corr_energy_2h2p_for_1h1p_b(l_part,j_hole) += contrib
! Couple of orbital contribution for the double 1h1p
corr_energy_2h2p_for_1h1p_double(i_hole,l_part) += contrib
corr_energy_2h2p_for_1h1p_double(l_part,i_hole) += contrib
corr_energy_2h2p_ab_2_orb(i_hole,j_hole) += contrib
corr_energy_2h2p_ab_2_orb(j_hole,i_hole) += contrib
corr_energy_2h2p_ab_2_orb(i_hole,k_part) += contrib
corr_energy_2h2p_ab_2_orb(k_part,i_hole) += contrib
corr_energy_2h2p_ab_2_orb(k_part,l_part) += contrib
corr_energy_2h2p_ab_2_orb(l_part,k_part) += contrib
enddo
enddo
enddo
enddo
! alpha alpha correlation energy
do i = 1, n_inact_orb
i_hole = list_inact(i)
do j = i+1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_virt_orb
k_part = list_virt(k)
do l = k+1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 1
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_2h2p += contrib
! Single orbital contribution
corr_energy_2h2p_per_orb_aa(i_hole) += contrib
corr_energy_2h2p_per_orb_aa(k_part) += contrib
! Couple of orbital contribution for the single 1h1p
corr_energy_2h2p_for_1h1p_a(i_hole,k_part) += contrib
corr_energy_2h2p_for_1h1p_a(k_part,i_hole) += contrib
enddo
enddo
enddo
enddo
! beta beta correlation energy
do i = 1, n_inact_orb
i_hole = list_inact(i)
do j = i+1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_virt_orb
k_part = list_virt(k)
do l = k+1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 2
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_2h2p += contrib
! Single orbital contribution
corr_energy_2h2p_per_orb_bb(i_hole) += contrib
corr_energy_2h2p_per_orb_bb(k_part) += contrib
corr_energy_2h2p_for_1h1p_b(i_hole,k_part) += contrib
corr_energy_2h2p_for_1h1p_b(k_part,i_hole) += contrib
! Two particle correlation energy
corr_energy_2h2p_bb_2_orb(i_hole,j_hole) += contrib
corr_energy_2h2p_bb_2_orb(j_hole,i_hole) += contrib
corr_energy_2h2p_bb_2_orb(i_hole,k_part) += contrib
corr_energy_2h2p_bb_2_orb(k_part,i_hole) += contrib
corr_energy_2h2p_bb_2_orb(k_part,l_part) += contrib
corr_energy_2h2p_bb_2_orb(l_part,k_part) += contrib
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, corr_energy_2h1p_ab_bb_per_2_orb, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_a, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_b, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h1p_for_1h1p_double, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_ab, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_aa, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_2h1p_per_orb_bb, (mo_tot_num)]
&BEGIN_PROVIDER [ double precision, total_corr_e_2h1p]
use bitmasks
implicit none
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,k,l
integer :: i_hole,j_hole,k_part,l_part
double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
double precision :: diag_H_mat_elem
integer :: i_ok,ispin
! Alpha - Beta correlation energy
total_corr_e_2h1p = 0.d0
corr_energy_2h1p_per_orb_ab = 0.d0
corr_energy_2h1p_per_orb_aa = 0.d0
corr_energy_2h1p_per_orb_bb = 0.d0
corr_energy_2h1p_ab_bb_per_2_orb = 0.d0
corr_energy_2h1p_for_1h1p_a = 0.d0
corr_energy_2h1p_for_1h1p_b = 0.d0
corr_energy_2h1p_for_1h1p_double = 0.d0
do i = 1, n_inact_orb
i_hole = list_inact(i)
do k = 1, n_act_orb
k_part = list_act(k)
do j = 1, n_inact_orb
j_hole = list_inact(j)
do l = 1, n_virt_orb
l_part = list_virt(l)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_2h1p += contrib
corr_energy_2h1p_ab_bb_per_2_orb(i_hole,j_hole) += contrib
corr_energy_2h1p_per_orb_ab(i_hole) += contrib
corr_energy_2h1p_per_orb_ab(l_part) += contrib
enddo
enddo
enddo
enddo
! Alpha Alpha spin correlation energy
do i = 1, n_inact_orb
i_hole = list_inact(i)
do j = i+1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_act_orb
k_part = list_act(k)
do l = 1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 1
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_2h1p += contrib
corr_energy_2h1p_per_orb_aa(i_hole) += contrib
corr_energy_2h1p_per_orb_aa(l_part) += contrib
enddo
enddo
enddo
enddo
! Beta Beta correlation energy
do i = 1, n_inact_orb
i_hole = list_inact(i)
do j = i+1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_act_orb
k_part = list_act(k)
do l = 1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 2
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
corr_energy_2h1p_ab_bb_per_2_orb(i_hole,j_hole) += contrib
total_corr_e_2h1p += contrib
corr_energy_2h1p_per_orb_bb(i_hole) += contrib
corr_energy_2h1p_per_orb_aa(l_part) += contrib
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_ab, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_1h2p_two_orb, (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_aa, (mo_tot_num)]
&BEGIN_PROVIDER [double precision, corr_energy_1h2p_per_orb_bb, (mo_tot_num)]
&BEGIN_PROVIDER [ double precision, total_corr_e_1h2p]
use bitmasks
implicit none
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,k,l
integer :: i_hole,j_hole,k_part,l_part
double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
double precision :: diag_H_mat_elem
integer :: i_ok,ispin
! Alpha - Beta correlation energy
total_corr_e_1h2p = 0.d0
corr_energy_1h2p_per_orb_ab = 0.d0
corr_energy_1h2p_per_orb_aa = 0.d0
corr_energy_1h2p_per_orb_bb = 0.d0
do i = 1, n_virt_orb
i_hole = list_virt(i)
do k = 1, n_act_orb
k_part = list_act(k)
do j = 1, n_inact_orb
j_hole = list_inact(j)
do l = 1, n_virt_orb
l_part = list_virt(l)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_1h2p += contrib
corr_energy_1h2p_per_orb_ab(i_hole) += contrib
corr_energy_1h2p_per_orb_ab(j_hole) += contrib
corr_energy_1h2p_two_orb(k_part,l_part) += contrib
corr_energy_1h2p_two_orb(l_part,k_part) += contrib
enddo
enddo
enddo
enddo
! Alpha Alpha correlation energy
do i = 1, n_virt_orb
i_hole = list_virt(i)
do j = 1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_act_orb
k_part = list_act(k)
do l = i+1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 1
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_1h2p += contrib
corr_energy_1h2p_per_orb_aa(i_hole) += contrib
corr_energy_1h2p_per_orb_ab(j_hole) += contrib
corr_energy_1h2p_two_orb(k_part,l_part) += contrib
corr_energy_1h2p_two_orb(l_part,k_part) += contrib
enddo
enddo
enddo
enddo
! Beta Beta correlation energy
do i = 1, n_virt_orb
i_hole = list_virt(i)
do j = 1, n_inact_orb
j_hole = list_inact(j)
do k = 1, n_act_orb
k_part = list_act(k)
do l = i+1,n_virt_orb
l_part = list_virt(l)
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
exc = get_mo_bielec_integral_schwartz(i_hole,j_hole,l_part,k_part,mo_integrals_map)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 2
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = hij - exc
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_1h2p += contrib
corr_energy_1h2p_per_orb_bb(i_hole) += contrib
corr_energy_1h2p_per_orb_ab(j_hole) += contrib
corr_energy_1h2p_two_orb(k_part,l_part) += contrib
corr_energy_1h2p_two_orb(l_part,k_part) += contrib
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, corr_energy_1h1p_spin_flip_per_orb, (mo_tot_num)]
&BEGIN_PROVIDER [ double precision, total_corr_e_1h1p_spin_flip]
use bitmasks
implicit none
integer(bit_kind) :: key_tmp(N_int,2)
integer :: i,j,k,l
integer :: i_hole,j_hole,k_part,l_part
double precision :: get_mo_bielec_integral_schwartz,hij,delta_e,exc,contrib
double precision :: diag_H_mat_elem
integer :: i_ok,ispin
! Alpha - Beta correlation energy
total_corr_e_1h1p_spin_flip = 0.d0
corr_energy_1h1p_spin_flip_per_orb = 0.d0
do i = 1, n_inact_orb
i_hole = list_inact(i)
do k = 1, n_act_orb
k_part = list_act(k)
do j = 1, n_act_orb
j_hole = list_act(j)
do l = 1, n_virt_orb
l_part = list_virt(l)
key_tmp = ref_bitmask
ispin = 2
call do_mono_excitation(key_tmp,i_hole,k_part,ispin,i_ok)
if(i_ok .ne.1)cycle
ispin = 1
call do_mono_excitation(key_tmp,j_hole,l_part,ispin,i_ok)
if(i_ok .ne.1)cycle
delta_e = -(ref_bitmask_energy - diag_H_mat_elem(key_tmp,N_int))
hij = get_mo_bielec_integral_schwartz(i_hole,j_hole,k_part,l_part,mo_integrals_map)
contrib = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij*hij))
total_corr_e_1h1p_spin_flip += contrib
corr_energy_1h1p_spin_flip_per_orb(i_hole) += contrib
enddo
enddo
enddo
enddo
END_PROVIDER