LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-26 06:22:04 +02:00
Code Releases Activity

Merge pull request #252 from QuantumPackage/dev-stable-rdm
Some checks failed
continuous-integration/drone/push Build is failing
Details

Browse Source 
Only one convention for TWO-rDM and all factors updated for basis_correction
This commit is contained in:
Anthony Scemama 2023-02-27 22:55:44 +01:00 committed by GitHub
commit 67e1477af8
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
15 changed files with 223 additions and 225 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
external/qp2-dependencies vendored

@ -1 +1 @@
Subproject commit ce14f57b50511825a9fedb096749200779d3f4d4
Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8

3
src/basis_correction/TODO
View File

@ -1 +1,4 @@
change all correlation functionals with the pbe_on_top general
factor 2 in two-rdm involved in:
on-top, mu(r), pbe-on-top, sc_basis_corr and so on
casscf : state_av_act_2_rdm_spin_trace_mo

24
src/basis_correction/pbe_on_top.irp.f
View File

@ -41,13 +41,15 @@
if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
! factor 2 because convention N(N-1)/ 2 in provider on_top_cas_mu_r
on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate)
else
! You take the on-top of the CAS wave function computed separately
! No factor 2 because convention N(N-1) in provider total_cas_on_top_density
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We take the extrapolated on-top pair density * 2 because of normalization
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
! We take the extrapolated on-top pair density
on_top_extrap = mu_correction_of_on_top(mu,on_top)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
@ -103,13 +105,15 @@
if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
! factor 2 because convention N(N-1)/ 2 in provider on_top_cas_mu_r
on_top = 2.d0 * on_top_cas_mu_r(ipoint,istate)
else
! You take the on-top of the CAS wave function computed separately
! No factor 2 because convention N(N-1) in provider total_cas_on_top_density
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We take the extrapolated on-top pair density * 2 because of normalization
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top)
! We take the extrapolated on-top pair density
on_top_extrap = mu_correction_of_on_top(mu,on_top)
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)
@ -165,13 +169,15 @@
if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
! factor 2 because convention N(N-1)/ 2 in provider on_top_cas_mu_r
on_top = 1.d0 * on_top_cas_mu_r(ipoint,istate)
else
! You take the on-top of the CAS wave function computed separately
! No factor 2 because convention N(N-1) in provider total_cas_on_top_density
on_top = total_cas_on_top_density(ipoint,istate)
endif
! We DO NOT take the extrapolated on-top pair density, but there is * 2 because of normalization
on_top_extrap = 2.d0 * on_top
! We DO NOT take the extrapolated on-top pair density
on_top_extrap = on_top
call ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top_extrap,eps_c_md_on_top_PBE)

6
src/basis_correction/print_routine.irp.f
View File

@ -18,7 +18,7 @@ subroutine print_basis_correction
print*, ''
print*, 'For more details look at Journal of Chemical Physics 149, 194301 1-15 (2018) '
print*, ' Journal of Physical Chemistry Letters 10, 2931-2937 (2019) '
print*, ' ???REF SC?'
print*, ' Journal of Chemical Physics 152, 174104 (2020) '
print*, '****************************************'
print*, '****************************************'
print*, 'mu_of_r_potential = ',mu_of_r_potential
@ -56,14 +56,14 @@ subroutine print_basis_correction
print*,''
print*,'********************************************'
print*,'********************************************'
print*,'+) PBE-on-top Ecmd functional : (??????? REF-SCF ??????????)'
print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) '
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate)
enddo
print*,''
print*,'********************************************'
print*,'+) PBE-on-top no spin polarization Ecmd functional : (??????? REF-SCF ??????????)'
print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)'
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION'
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)

2
src/mu_of_r/f_psi_i_a_v_utils.irp.f
View File

@ -194,7 +194,7 @@ subroutine give_f_aa_val_ab(r1,r2,f_aa_val_ab,two_bod_dens,istate)
do b = 1, n_act_orb ! 2
do c = 1, n_act_orb ! 1
do d = 1, n_act_orb ! 2
rho = mos_array_act_r1(c) * mos_array_act_r2(d) * act_2_rdm_ab_mo(d,c,b,a,istate)
rho = mos_array_act_r1(c) * mos_array_act_r2(d) * 0.5d0 * act_2_rdm_ab_mo(d,c,b,a,istate)
rho_tilde(b,a) += rho
two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b)
enddo

2
src/mu_of_r/f_psi_utils.irp.f
View File

@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, full_occ_2_rdm_cntrctd_trans, (n_points_final_
do j = 1, n_core_inact_act_orb
do i = 1, n_core_inact_act_orb
! 1 2 1 2
full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint)
full_occ_2_rdm_cntrctd_trans(ipoint,k,l,istate) += 0.5d0 * full_occ_2_rdm_ab_mo(i,j,k,l,istate) * core_inact_act_mos_in_r_array(j,ipoint) * core_inact_act_mos_in_r_array(i,ipoint)
enddo
enddo
enddo

2
src/mu_of_r/mu_of_r_conditions.irp.f
View File

@ -110,7 +110,7 @@
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
f_psi = f_psi_cas_ab(ipoint,istate)
on_top = on_top_cas_mu_r(ipoint,istate)
on_top = on_top_cas_mu_r(ipoint,istate)
if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then
w_psi = 1.d+10
else

23
src/two_body_rdm/act_2_rdm.irp.f
View File

@ -4,21 +4,18 @@
BEGIN_DOC
! 12 12
! 1 2 1 2 == <ij|kl>
! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons
! act_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta+beta/alpha electrons
!
! <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>
! <Psi_{istate}| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi_{istate}>
!
! + <Psi_{istate}| a^{\dagger}_{i \beta} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \beta} |Psi_{istate}>
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act}
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is alpha, electron 2 is beta
!
! act_2_rdm_ab_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta
!
! Therefore you don't necessary have symmetry between electron 1 and 2
END_DOC
integer :: ispin
double precision :: wall_1, wall_2
@ -44,6 +41,7 @@
endif
call wall_time(wall_2)
print*,'Wall time to provide act_2_rdm_ab_mo',wall_2 - wall_1
act_2_rdm_ab_mo *= 2.d0
END_PROVIDER
@ -56,7 +54,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
@ -84,6 +82,7 @@
call wall_time(wall_2)
print*,'Wall time to provide act_2_rdm_aa_mo',wall_2 - wall_1
act_2_rdm_aa_mo *= 2.d0
END_PROVIDER
@ -96,7 +95,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
@ -124,6 +123,7 @@
call wall_time(wall_2)
print*,'Wall time to provide act_2_rdm_bb_mo',wall_2 - wall_1
act_2_rdm_bb_mo *= 2.d0
END_PROVIDER
BEGIN_PROVIDER [double precision, act_2_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
@ -135,7 +135,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec}^{act} * (N_{elec}^{act} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec}^{act} * (N_{elec}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
@ -161,6 +161,7 @@
call ezfio_set_two_body_rdm_io_two_body_rdm_spin_trace("Read")
endif
act_2_rdm_spin_trace_mo *= 2.d0
call wall_time(wall_2)
print*,'Wall time to provide act_2_rdm_spin_trace_mo',wall_2 - wall_1
END_PROVIDER

52
src/two_body_rdm/example.irp.f
View File

@ -117,10 +117,10 @@ subroutine routine_active_only
endif
wee_ab(istate) += vijkl * rdmab
wee_aa(istate) += vijkl * rdmaa
wee_bb(istate) += vijkl * rdmbb
wee_tot(istate) += vijkl * rdmtot
wee_ab(istate) += 0.5d0 * vijkl * rdmab
wee_aa(istate) += 0.5d0 * vijkl * rdmaa
wee_bb(istate) += 0.5d0 * vijkl * rdmbb
wee_tot(istate) += 0.5d0 * vijkl * rdmtot
enddo
enddo
@ -144,13 +144,13 @@ subroutine routine_active_only
print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate)
print*,'--------------------------'
print*,'accu_aa = ',accu_aa
print*,'N_a (N_a-1)/2 = ', elec_alpha_num*(elec_alpha_num-1)*0.5
print*,'N_a (N_a-1) = ', elec_alpha_num*(elec_alpha_num-1)
print*,'accu_bb = ',accu_bb
print*,'N_b (N_b-1)/2 = ', elec_beta_num*(elec_beta_num-1)*0.5
print*,'2 N_b (N_b-1) = ', elec_beta_num*(elec_beta_num-1)*2
print*,'accu_ab = ',accu_ab
print*,'N_a N_b = ', elec_beta_num*elec_alpha_num
print*,'accu_tot = ',accu_tot
print*,'Ne(Ne-1)/2 = ',(elec_num-1)*elec_num * 0.5
print*,'Ne(Ne-1) = ',(elec_num-1)*elec_num
enddo
wee_aa_st_av = 0.d0
wee_bb_st_av = 0.d0
@ -192,10 +192,10 @@ subroutine routine_active_only
print*,spin_trace,rdm_tot_st_av,dabs(spin_trace - rdm_tot_st_av)
endif
wee_aa_st_av += vijkl * rdm_aa_st_av
wee_bb_st_av += vijkl * rdm_bb_st_av
wee_ab_st_av += vijkl * rdm_ab_st_av
wee_tot_st_av += vijkl * rdm_tot_st_av
wee_aa_st_av += 0.5d0 * vijkl * rdm_aa_st_av
wee_bb_st_av += 0.5d0 * vijkl * rdm_bb_st_av
wee_ab_st_av += 0.5d0 * vijkl * rdm_ab_st_av
wee_tot_st_av += 0.5d0 * vijkl * rdm_tot_st_av
enddo
enddo
enddo
@ -279,10 +279,10 @@ subroutine routine_full_mos
rdmbb = full_occ_2_rdm_bb_mo(l,k,j,i,istate)
rdmtot = full_occ_2_rdm_spin_trace_mo(l,k,j,i,istate)
wee_ab(istate) += vijkl * rdmab
wee_aa(istate) += vijkl * rdmaa
wee_bb(istate) += vijkl * rdmbb
wee_tot(istate)+= vijkl * rdmtot
wee_ab(istate) += 0.5d0 * vijkl * rdmab
wee_aa(istate) += 0.5d0 * vijkl * rdmaa
wee_bb(istate) += 0.5d0 * vijkl * rdmbb
wee_tot(istate)+= 0.5d0 * vijkl * rdmtot
enddo
enddo
aa_norm(istate) += full_occ_2_rdm_aa_mo(j,i,j,i,istate)
@ -310,18 +310,19 @@ subroutine routine_full_mos
print*,'Normalization of two-rdms '
print*,''
print*,'aa_norm(istate) = ',aa_norm(istate)
print*,'N_alpha(N_alpha-1)/2 = ',elec_num_tab(1) * (elec_num_tab(1) - 1)/2
print*,'N_alpha(N_alpha-1) = ',elec_num_tab(1) * (elec_num_tab(1) - 1)
print*,''
print*,'bb_norm(istate) = ',bb_norm(istate)
print*,'N_alpha(N_alpha-1)/2 = ',elec_num_tab(2) * (elec_num_tab(2) - 1)/2
print*,'N_alpha(N_alpha-1) = ',elec_num_tab(2) * (elec_num_tab(2) - 1)
print*,''
print*,'ab_norm(istate) = ',ab_norm(istate)
print*,'N_alpha * N_beta = ',elec_num_tab(1) * elec_num_tab(2)
print*,'N_alpha * N_beta *2 = ',elec_num_tab(1) * elec_num_tab(2) * 2
print*,''
print*,'tot_norm(istate) = ',tot_norm(istate)
print*,'N(N-1)/2 = ',elec_num*(elec_num - 1)/2
print*,'N(N-1) = ',elec_num*(elec_num - 1)
enddo
! return
wee_aa_st_av = 0.d0
wee_bb_st_av = 0.d0
wee_ab_st_av = 0.d0
@ -341,10 +342,10 @@ subroutine routine_full_mos
rdm_ab_st_av = state_av_full_occ_2_rdm_ab_mo(l,k,j,i)
rdm_tot_st_av = state_av_full_occ_2_rdm_spin_trace_mo(l,k,j,i)
wee_aa_st_av += vijkl * rdm_aa_st_av
wee_bb_st_av += vijkl * rdm_bb_st_av
wee_ab_st_av += vijkl * rdm_ab_st_av
wee_tot_st_av += vijkl * rdm_tot_st_av
wee_aa_st_av += 0.5d0 * vijkl * rdm_aa_st_av
wee_bb_st_av += 0.5d0 * vijkl * rdm_bb_st_av
wee_ab_st_av += 0.5d0 * vijkl * rdm_ab_st_av
wee_tot_st_av += 0.5d0 * vijkl * rdm_tot_st_av
enddo
enddo
enddo
@ -354,17 +355,12 @@ subroutine routine_full_mos
print*,''
print*,'STATE AVERAGE ENERGY '
print*,'wee_aa_st_av = ',wee_aa_st_av
print*,'wee_aa_st_av_2 = ',wee_aa_st_av_2
print*,'wee_bb_st_av = ',wee_bb_st_av
print*,'wee_bb_st_av_2 = ',wee_bb_st_av_2
print*,'wee_ab_st_av = ',wee_ab_st_av
print*,'wee_ab_st_av_2 = ',wee_ab_st_av_2
print*,'Sum of components = ',wee_aa_st_av + wee_bb_st_av + wee_ab_st_av
print*,'Sum of components_2 = ',wee_aa_st_av_2 + wee_bb_st_av_2 + wee_ab_st_av_2
print*,''
print*,'Full energy '
print*,'wee_tot_st_av = ',wee_tot_st_av
print*,'wee_tot_st_av_2 = ',wee_tot_st_av_2
print*,'wee_tot_st_av_3 = ',wee_tot_st_av_3
end

154
src/two_body_rdm/full_orb_2_rdm.irp.f
View File

@ -4,22 +4,18 @@
full_occ_2_rdm_ab_mo = 0.d0
integer :: i,j,k,l,iorb,jorb,korb,lorb,istate
BEGIN_DOC
! full_occ_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons
! full_occ_2_rdm_ab_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta + beta/alpha electrons
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
!
! + <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \beta} |Psi>
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta}
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} * 2
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA
!
! full_occ_2_rdm_ab_mo(i,j,k,l,istate) = i:alpha, j:beta, j:alpha, l:beta
!
! Therefore you don't necessary have symmetry between electron 1 and 2
!
! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO ARE SET TO ZERO
END_DOC
full_occ_2_rdm_ab_mo = 0.d0
@ -50,7 +46,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = 2.d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -64,7 +60,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = 2.d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -76,7 +72,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0
full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 2.D0
enddo
enddo
@ -93,7 +89,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb,istate) = 2.d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -107,7 +103,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb,istate) = 2.d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -119,7 +115,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 1.D0
full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb,istate) = 2.D0
enddo
enddo
endif
@ -139,7 +135,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -172,11 +168,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -186,8 +182,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
@ -203,11 +199,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -217,8 +213,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
endif
@ -237,7 +233,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -270,11 +266,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -284,8 +280,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
@ -301,11 +297,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -315,8 +311,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
endif
@ -335,14 +331,14 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero
! The two-electron energy of each state can be computed as:
!
! \sum_{i,j,k,l = 1, n_core_inact_act_orb} full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) * < ii jj | kk ll >
! \sum_{i,j,k,l = 1, n_core_inact_act_orb} full_occ_2_rdm_spin_trace_mo(i,j,k,l,istate) * 1/2 * < ii jj | kk ll >
!
! with ii = list_core_inact_act(i), jj = list_core_inact_act(j), kk = list_core_inact_act(k), ll = list_core_inact_act(l)
END_DOC
@ -377,11 +373,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -390,8 +386,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
if (.not.no_core_density)then
@ -403,11 +399,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
enddo
enddo
enddo
@ -416,8 +412,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
endif
@ -433,11 +429,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -446,8 +442,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
if (.not.no_core_density)then
@ -459,11 +455,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! 1 2 1 2 : EXCHANGE TERM
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb,istate) += -1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -472,8 +468,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 0.5d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0d0
full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb,istate) -= 1.0d0
enddo
enddo
endif
@ -489,14 +485,14 @@
korb = list_inact(k)
! ALPHA INACTIVE - BETA ACTIVE
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! beta alph beta alph
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_beta(jorb,iorb,istate)
! BETA INACTIVE - ALPHA ACTIVE
! beta alph beta alpha
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0d0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -506,8 +502,8 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0
full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0D0
full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 1.0D0
enddo
enddo
@ -523,14 +519,14 @@
korb = list_core(k)
!! BETA ACTIVE - ALPHA CORE
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0D0 * one_e_dm_mo_beta(jorb,iorb,istate)
! beta alph beta alph
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_beta(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0D0 * one_e_dm_mo_beta(jorb,iorb,istate)
!! ALPHA ACTIVE - BETA CORE
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb,istate) += 1.0D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
! beta alph beta alph
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 0.5D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb,istate) += 1.0D0 * one_e_dm_mo_alpha(jorb,iorb,istate)
enddo
enddo
enddo
@ -540,8 +536,8 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 0.5D0
full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 0.5D0
full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb,istate) += 1.0D0
full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb,istate) += 1.0D0
enddo
enddo

21
src/two_body_rdm/state_av_act_2rdm.irp.f
View File

@ -2,21 +2,16 @@
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! state_av_act_2_rdm_ab_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs
! state_av_act_2_rdm_ab_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha/beta+beta/alpha electron pairs
!
! = \sum_{istate} w(istate) * <Psi_{istate}| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi_{istate}>
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act}
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * N_{\beta}^{act} * 2
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is alpha, electron 2 is beta
!
! state_av_act_2_rdm_ab_mo(i,j,k,l) = i:alpha, j:beta, j:alpha, l:beta
!
! Therefore you don't necessary have symmetry between electron 1 and 2
END_DOC
allocate(state_weights(N_states))
state_weights = state_average_weight
@ -34,6 +29,8 @@
call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(wall_2)
print*,'Wall time to provide state_av_act_2_rdm_ab_mo',wall_2 - wall_1
! factor 2 to have the correct normalization factor
state_av_act_2_rdm_ab_mo *= 2.d0
END_PROVIDER
@ -48,7 +45,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha}^{act} * (N_{\alpha}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
@ -63,6 +60,8 @@
call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(wall_2)
print*,'Wall time to provide state_av_act_2_rdm_aa_mo',wall_2 - wall_1
! factor 2 to have the correct normalization factor
state_av_act_2_rdm_aa_mo *= 2.d0
END_PROVIDER
@ -76,7 +75,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta}^{act} * (N_{\beta}^{act} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
@ -91,6 +90,8 @@
call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(wall_2)
print*,'Wall time to provide state_av_act_2_rdm_bb_mo',wall_2 - wall_1
! factor 2 to have the correct normalization factor
state_av_act_2_rdm_bb_mo *= 2.d0
END_PROVIDER
@ -104,7 +105,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC

148
src/two_body_rdm/state_av_full_orb_2_rdm.irp.f
View File

@ -4,22 +4,16 @@
state_av_full_occ_2_rdm_ab_mo = 0.d0
integer :: i,j,k,l,iorb,jorb,korb,lorb
BEGIN_DOC
! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha/beta electrons
! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha/beta + beta/alpha electrons
!
! = \sum_{istate} w(istate) * <Psi_{istate}| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi_{istate}>
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta}
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * N_{\beta} * 2
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! !!!!! WARNING !!!!! For efficiency reasons, electron 1 is ALPHA, electron 2 is BETA
!
! state_av_full_occ_2_rdm_ab_mo(i,j,k,l) = i:alpha, j:beta, j:alpha, l:beta
!
! Therefore you don't necessary have symmetry between electron 1 and 2
!
! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO is set to zero
END_DOC
state_av_full_occ_2_rdm_ab_mo = 0.d0
@ -47,7 +41,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = 2.d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -61,7 +55,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = 2.d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -73,7 +67,7 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 1.D0
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0
enddo
enddo
@ -90,7 +84,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,iorb) = 2.d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -104,7 +98,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_ab_mo(jorb,korb,iorb,korb) = 2.d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -116,7 +110,7 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 1.D0
state_av_full_occ_2_rdm_ab_mo(korb,jorb,korb,jorb) = 2.D0
enddo
enddo
endif
@ -135,7 +129,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\alpha} * (N_{\alpha} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -167,11 +161,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -181,8 +175,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
@ -198,11 +192,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -212,8 +206,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_aa_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
endif
@ -231,7 +225,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{\beta} * (N_{\beta} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -263,11 +257,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -277,8 +271,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
@ -294,11 +288,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -308,8 +302,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_bb_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
endif
@ -328,7 +322,7 @@
!
! WHERE ALL ORBITALS (i,j,k,l) BELONGS TO ALL OCCUPIED ORBITALS : core, inactive and active
!
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)/2
! THE NORMALIZATION (i.e. sum of diagonal elements) IS SET TO N_{elec} * (N_{elec} - 1)
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
@ -364,11 +358,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -377,8 +371,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
if (.not.no_core_density)then
@ -390,11 +384,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
enddo
enddo
enddo
@ -403,8 +397,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
endif
@ -420,11 +414,11 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -433,8 +427,8 @@
jorb = list_inact(j)
do k = 1, n_inact_orb
korb = list_inact(k)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
if (.not.no_core_density)then
@ -446,11 +440,11 @@
do k = 1, n_core_orb
korb = list_core(k)
! 1 2 1 2 : DIRECT TERM
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! 1 2 1 2 : EXCHANGE TERM
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,korb,iorb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,iorb,korb) += -1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -459,8 +453,8 @@
jorb = list_core(j)
do k = 1, n_core_orb
korb = list_core(k)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 0.5d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,jorb,korb) -= 1.0d0
enddo
enddo
endif
@ -476,14 +470,14 @@
korb = list_inact(k)
! ALPHA INACTIVE - BETA ACTIVE
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! beta alph beta alph
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_beta_average(jorb,iorb)
! BETA INACTIVE - ALPHA ACTIVE
! beta alph beta alpha
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5d0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0d0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -493,8 +487,8 @@
do k = 1, n_inact_orb
korb = list_inact(k)
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 0.5D0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0d0
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0d0
enddo
enddo
@ -510,14 +504,14 @@
korb = list_core(k)
!! BETA ACTIVE - ALPHA CORE
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0D0 * one_e_dm_mo_beta_average(jorb,iorb)
! beta alph beta alph
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_beta_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0D0 * one_e_dm_mo_beta_average(jorb,iorb)
!! ALPHA ACTIVE - BETA CORE
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,iorb,korb) += 1.0D0 * one_e_dm_mo_alpha_average(jorb,iorb)
! beta alph beta alph
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 0.5D0 * one_e_dm_mo_alpha_average(jorb,iorb)
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,iorb) += 1.0D0 * one_e_dm_mo_alpha_average(jorb,iorb)
enddo
enddo
enddo
@ -527,8 +521,8 @@
do k = 1, n_core_orb
korb = list_core(k)
! alph beta alph beta
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 0.5D0
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 0.5D0
state_av_full_occ_2_rdm_spin_trace_mo(korb,jorb,korb,jorb) += 1.0D0
state_av_full_occ_2_rdm_spin_trace_mo(jorb,korb,jorb,korb) += 1.0D0
enddo
enddo

2
src/two_body_rdm/test_2_rdm.irp.f
View File

@ -2,7 +2,7 @@ program test_2_rdm
implicit none
read_wf = .True.
touch read_wf
call routine_active_only
! call routine_active_only
call routine_full_mos
end

3
src/two_body_rdm/two_e_dm_mo.irp.f
View File

@ -29,7 +29,8 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)]
enddo
enddo
enddo
two_e_dm_mo(:,:,:,:) = two_e_dm_mo(:,:,:,:) * 2.d0
two_e_dm_mo(:,:,:,:) = two_e_dm_mo(:,:,:,:)
! * 2.d0
END_PROVIDER

4
src/two_rdm_routines/davidson_like_2rdm.irp.f
View File

@ -4,8 +4,8 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
BEGIN_DOC
! if ispin == 1 :: alpha/alpha 2rdm
! == 2 :: beta /beta 2rdm
! == 3 :: alpha/beta 2rdm
! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba))
! == 3 :: alpha/beta + beta/alpha 2rdm
! == 4 :: spin traced 2rdm :: aa + bb + ab + ba
!
! Assumes that the determinants are in psi_det
!