From fe27080069d6685ff1ff94fa9ce1d78d6cc853c2 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 15:27:39 +0100 Subject: [PATCH 1/6] beginning to introduce a factor 2 in two-rdm --- external/qp2-dependencies | 2 +- src/basis_correction/TODO | 2 ++ src/basis_correction/print_routine.irp.f | 6 +++--- src/two_body_rdm/example.irp.f | 2 +- src/two_body_rdm/test_2_rdm.irp.f | 2 +- src/two_rdm_routines/davidson_like_2rdm.irp.f | 4 ++-- 6 files changed, 10 insertions(+), 8 deletions(-) diff --git a/external/qp2-dependencies b/external/qp2-dependencies index ce14f57b..f40bde09 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit ce14f57b50511825a9fedb096749200779d3f4d4 +Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8 diff --git a/src/basis_correction/TODO b/src/basis_correction/TODO index e28d593a..36c438c8 100644 --- a/src/basis_correction/TODO +++ b/src/basis_correction/TODO @@ -1 +1,3 @@ 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 diff --git a/src/basis_correction/print_routine.irp.f b/src/basis_correction/print_routine.irp.f index 67c5c6c2..c2558d22 100644 --- a/src/basis_correction/print_routine.irp.f +++ b/src/basis_correction/print_routine.irp.f @@ -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) diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index de3d97b9..67de9df4 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -120,7 +120,7 @@ subroutine routine_active_only wee_ab(istate) += vijkl * rdmab wee_aa(istate) += vijkl * rdmaa wee_bb(istate) += vijkl * rdmbb - wee_tot(istate) += vijkl * rdmtot + wee_tot(istate) += vijkl * rdmtot enddo enddo diff --git a/src/two_body_rdm/test_2_rdm.irp.f b/src/two_body_rdm/test_2_rdm.irp.f index 123261d8..4eb8f9f0 100644 --- a/src/two_body_rdm/test_2_rdm.irp.f +++ b/src/two_body_rdm/test_2_rdm.irp.f @@ -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 diff --git a/src/two_rdm_routines/davidson_like_2rdm.irp.f b/src/two_rdm_routines/davidson_like_2rdm.irp.f index 2e5aa4d1..ad7a3b21 100644 --- a/src/two_rdm_routines/davidson_like_2rdm.irp.f +++ b/src/two_rdm_routines/davidson_like_2rdm.irp.f @@ -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 ! From 3ba5d3b540424df9cacfe578d24b209ed93e018e Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 15:45:09 +0100 Subject: [PATCH 2/6] factor two introduced in non active only non state average two-rdm, it works with example.irp.f --- src/basis_correction/TODO | 1 + src/two_body_rdm/act_2_rdm.irp.f | 4 ++ src/two_body_rdm/example.irp.f | 47 ++++++++++++------------ src/two_body_rdm/state_av_act_2rdm.irp.f | 9 +++-- src/two_body_rdm/two_e_dm_mo.irp.f | 3 +- 5 files changed, 37 insertions(+), 27 deletions(-) diff --git a/src/basis_correction/TODO b/src/basis_correction/TODO index 36c438c8..64a6ddeb 100644 --- a/src/basis_correction/TODO +++ b/src/basis_correction/TODO @@ -1,3 +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 diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f index 41b28aea..61ae4e47 100644 --- a/src/two_body_rdm/act_2_rdm.irp.f +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -44,6 +44,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 @@ -84,6 +85,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 @@ -124,6 +126,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)] @@ -161,6 +164,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 diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 67de9df4..985f6a5d 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -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)/2 = ',(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)/2 = ',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 diff --git a/src/two_body_rdm/state_av_act_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f index db793047..9e3d1df0 100644 --- a/src/two_body_rdm/state_av_act_2rdm.irp.f +++ b/src/two_body_rdm/state_av_act_2rdm.irp.f @@ -34,6 +34,7 @@ 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 +! state_av_act_2_rdm_ab_mo *= 2.d0 END_PROVIDER @@ -48,7 +49,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 +64,7 @@ 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 +! state_av_act_2_rdm_aa_mo *= 2.d0 END_PROVIDER @@ -76,7 +78,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 +93,7 @@ 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 +! state_av_act_2_rdm_bb_mo *= 2.d0 END_PROVIDER @@ -104,7 +107,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 diff --git a/src/two_body_rdm/two_e_dm_mo.irp.f b/src/two_body_rdm/two_e_dm_mo.irp.f index a4dea15f..7e35fc7b 100644 --- a/src/two_body_rdm/two_e_dm_mo.irp.f +++ b/src/two_body_rdm/two_e_dm_mo.irp.f @@ -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 From fd63ab1355d481e90a3cef7401e934da63f04a19 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 17:33:43 +0100 Subject: [PATCH 3/6] state_av_full_occ_2_rdm_aa_mo work --- src/two_body_rdm/example.irp.f | 7 +- src/two_body_rdm/state_av_act_2rdm.irp.f | 6 +- .../state_av_full_orb_2_rdm.irp.f | 132 +++++++++--------- 3 files changed, 70 insertions(+), 75 deletions(-) diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 985f6a5d..01e971ba 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -322,7 +322,7 @@ subroutine routine_full_mos print*,'N(N-1)/2 = ',elec_num*(elec_num - 1) enddo - return +! return wee_aa_st_av = 0.d0 wee_bb_st_av = 0.d0 wee_ab_st_av = 0.d0 @@ -355,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 diff --git a/src/two_body_rdm/state_av_act_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f index 9e3d1df0..a97ec3f9 100644 --- a/src/two_body_rdm/state_av_act_2rdm.irp.f +++ b/src/two_body_rdm/state_av_act_2rdm.irp.f @@ -34,7 +34,7 @@ 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 -! state_av_act_2_rdm_ab_mo *= 2.d0 + state_av_act_2_rdm_ab_mo *= 2.d0 END_PROVIDER @@ -64,7 +64,7 @@ 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 -! state_av_act_2_rdm_aa_mo *= 2.d0 + state_av_act_2_rdm_aa_mo *= 2.d0 END_PROVIDER @@ -93,7 +93,7 @@ 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 -! state_av_act_2_rdm_bb_mo *= 2.d0 + state_av_act_2_rdm_bb_mo *= 2.d0 END_PROVIDER diff --git a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f index b3a5fe65..251b4d96 100644 --- a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f +++ b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f @@ -47,7 +47,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 +61,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 +73,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 +90,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 +104,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 +116,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 @@ -167,11 +167,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 +181,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 +198,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 +212,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 @@ -263,11 +263,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 +277,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 +294,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 +308,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 @@ -364,11 +364,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 +377,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 +390,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 +403,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 +420,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 +433,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 +446,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 +459,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 +476,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 +493,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 +510,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 +527,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 From c95f1ee0ac22e79aa0f7b785e7503b1853079c80 Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 17:42:44 +0100 Subject: [PATCH 4/6] changed all two-rdm with the normalization convtion to N(N-1) and not N(N-1)/2 --- src/two_body_rdm/example.irp.f | 4 +- src/two_body_rdm/full_orb_2_rdm.irp.f | 132 +++++++++++++------------- 2 files changed, 68 insertions(+), 68 deletions(-) diff --git a/src/two_body_rdm/example.irp.f b/src/two_body_rdm/example.irp.f index 01e971ba..30e2685a 100644 --- a/src/two_body_rdm/example.irp.f +++ b/src/two_body_rdm/example.irp.f @@ -150,7 +150,7 @@ subroutine routine_active_only 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 + print*,'Ne(Ne-1) = ',(elec_num-1)*elec_num enddo wee_aa_st_av = 0.d0 wee_bb_st_av = 0.d0 @@ -319,7 +319,7 @@ subroutine routine_full_mos 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) + print*,'N(N-1) = ',elec_num*(elec_num - 1) enddo ! return diff --git a/src/two_body_rdm/full_orb_2_rdm.irp.f b/src/two_body_rdm/full_orb_2_rdm.irp.f index fba88172..a3de8ea9 100644 --- a/src/two_body_rdm/full_orb_2_rdm.irp.f +++ b/src/two_body_rdm/full_orb_2_rdm.irp.f @@ -50,7 +50,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 +64,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 +76,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 +93,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 +107,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 +119,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 @@ -172,11 +172,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 +186,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 +203,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 +217,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 @@ -270,11 +270,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 +284,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 +301,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 +315,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 @@ -377,11 +377,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 +390,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 +403,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 +416,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 +433,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 +446,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 +459,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 +472,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 +489,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 +506,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 +523,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 +540,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 From 8515fcf93fae51d4658f2a9911f39d6bf1a2fced Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 17:48:42 +0100 Subject: [PATCH 5/6] updated DOC for correct normalization of two-rdms --- src/two_body_rdm/act_2_rdm.irp.f | 19 +++++++--------- src/two_body_rdm/full_orb_2_rdm.irp.f | 22 ++++++++----------- src/two_body_rdm/state_av_act_2rdm.irp.f | 12 +++++----- .../state_av_full_orb_2_rdm.irp.f | 16 +++++--------- 4 files changed, 27 insertions(+), 42 deletions(-) diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f index 61ae4e47..c550e991 100644 --- a/src/two_body_rdm/act_2_rdm.irp.f +++ b/src/two_body_rdm/act_2_rdm.irp.f @@ -4,21 +4,18 @@ BEGIN_DOC ! 12 12 ! 1 2 1 2 == -! 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 ! -! +! +! +! + ! ! 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 @@ -57,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 @@ -98,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 @@ -138,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 diff --git a/src/two_body_rdm/full_orb_2_rdm.irp.f b/src/two_body_rdm/full_orb_2_rdm.irp.f index a3de8ea9..78a28acf 100644 --- a/src/two_body_rdm/full_orb_2_rdm.irp.f +++ b/src/two_body_rdm/full_orb_2_rdm.irp.f @@ -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 ! -! +! +! +! + ! ! 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 @@ -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" ! @@ -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" ! @@ -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 diff --git a/src/two_body_rdm/state_av_act_2rdm.irp.f b/src/two_body_rdm/state_av_act_2rdm.irp.f index a97ec3f9..cd417a9d 100644 --- a/src/two_body_rdm/state_av_act_2rdm.irp.f +++ b/src/two_body_rdm/state_av_act_2rdm.irp.f @@ -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) * ! ! 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,7 @@ 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 @@ -64,6 +60,7 @@ 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 @@ -93,6 +90,7 @@ 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 diff --git a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f index 251b4d96..2e44665d 100644 --- a/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f +++ b/src/two_body_rdm/state_av_full_orb_2_rdm.irp.f @@ -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) * ! ! 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 @@ -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" ! @@ -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" ! @@ -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" ! From b51d72d5b905ef5ac198d42cc4b53ca90192d90d Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 18:35:51 +0100 Subject: [PATCH 6/6] changed the factor 2 in basis_correction and mu_of_r in order to adapt to new normalization factor --- src/basis_correction/pbe_on_top.irp.f | 24 +++++++++++++++--------- src/mu_of_r/f_psi_i_a_v_utils.irp.f | 2 +- src/mu_of_r/f_psi_utils.irp.f | 2 +- src/mu_of_r/mu_of_r_conditions.irp.f | 2 +- 4 files changed, 18 insertions(+), 12 deletions(-) diff --git a/src/basis_correction/pbe_on_top.irp.f b/src/basis_correction/pbe_on_top.irp.f index cc9cdec9..a25fd61b 100644 --- a/src/basis_correction/pbe_on_top.irp.f +++ b/src/basis_correction/pbe_on_top.irp.f @@ -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) diff --git a/src/mu_of_r/f_psi_i_a_v_utils.irp.f b/src/mu_of_r/f_psi_i_a_v_utils.irp.f index 427da199..0d08e193 100644 --- a/src/mu_of_r/f_psi_i_a_v_utils.irp.f +++ b/src/mu_of_r/f_psi_i_a_v_utils.irp.f @@ -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 diff --git a/src/mu_of_r/f_psi_utils.irp.f b/src/mu_of_r/f_psi_utils.irp.f index bdd76f18..95f10f36 100644 --- a/src/mu_of_r/f_psi_utils.irp.f +++ b/src/mu_of_r/f_psi_utils.irp.f @@ -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 diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f index 5c41acdc..f9c3b3b3 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -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