From 3ba5d3b540424df9cacfe578d24b209ed93e018e Mon Sep 17 00:00:00 2001 From: eginer Date: Mon, 27 Feb 2023 15:45:09 +0100 Subject: [PATCH] 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