PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-23 04:43:45 +01:00
Code Releases Activity

changed the factor 2 in basis_correction and mu_of_r in order to adapt to new normalization factor

Browse Source
This commit is contained in:
eginer 2023-02-27 18:35:51 +01:00
parent 8515fcf93f
commit b51d72d5b9
4 changed files with 18 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -41,13 +41,15 @@
if(mu_of_r_potential == "cas_ful")then if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r) ! 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 else
! You take the on-top of the CAS wave function computed separately ! 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) on_top = total_cas_on_top_density(ipoint,istate)
endif endif
! We take the extrapolated on-top pair density * 2 because of normalization ! We take the extrapolated on-top pair density
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top) 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) 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 if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r) ! 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 else
! You take the on-top of the CAS wave function computed separately ! 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) on_top = total_cas_on_top_density(ipoint,istate)
endif endif
! We take the extrapolated on-top pair density * 2 because of normalization ! We take the extrapolated on-top pair density
on_top_extrap = 2.d0 * mu_correction_of_on_top(mu,on_top) 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) 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 if(mu_of_r_potential == "cas_ful")then
! You take the on-top of the CAS wave function which is computed with mu(r) ! 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 else
! You take the on-top of the CAS wave function computed separately ! 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) on_top = total_cas_on_top_density(ipoint,istate)
endif endif
! We DO NOT take the extrapolated on-top pair density, but there is * 2 because of normalization ! We DO NOT take the extrapolated on-top pair density
on_top_extrap = 2.d0 * on_top 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) 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)

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 b = 1, n_act_orb ! 2
do c = 1, n_act_orb ! 1 do c = 1, n_act_orb ! 1
do d = 1, n_act_orb ! 2 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 rho_tilde(b,a) += rho
two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b) two_bod_dens += rho * mos_array_act_r1(a) * mos_array_act_r2(b)
enddo 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 j = 1, n_core_inact_act_orb
do i = 1, n_core_inact_act_orb do i = 1, n_core_inact_act_orb
! 1 2 1 2 ! 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 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 istate = 1, N_states
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
f_psi = f_psi_cas_ab(ipoint,istate) 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 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 w_psi = 1.d+10
else else