From 61ddf14ab9665ac8af4f2bd0aac07b4d071a5c69 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 28 Feb 2025 18:04:16 +0100 Subject: [PATCH] Added average mu(r) with rho^2 --- .../basis_correction/print_routine.irp.f | 3 +- src/mu_of_r/mu_of_r_conditions.irp.f | 30 +++++++++++++++++-- 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/plugins/local/basis_correction/print_routine.irp.f b/plugins/local/basis_correction/print_routine.irp.f index 8879fd5d..cc9744f2 100644 --- a/plugins/local/basis_correction/print_routine.irp.f +++ b/plugins/local/basis_correction/print_routine.irp.f @@ -79,7 +79,8 @@ subroutine print_basis_correction print*,'' print*,'**************' do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) , state ',istate,' = ',mu_average_prov(istate) + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) [rho ], state ',istate,' = ',mu_average_prov(istate) + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' Average mu(r) [rho^2], state ',istate,' = ',mu_average_prov2(istate) enddo end 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 ba3675d4..d13605f6 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -211,7 +211,7 @@ END_PROVIDER - BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)] +BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)] implicit none BEGIN_DOC ! average value of mu(r) weighted with the total one-e density and divided by the number of electrons @@ -233,7 +233,33 @@ enddo mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate) enddo - END_PROVIDER +END_PROVIDER + +BEGIN_PROVIDER [double precision, mu_average_prov2, (N_states)] + implicit none + BEGIN_DOC + ! average value of mu(r) weighted with square of the total one-e density + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the one- and two-body density matrix are excluded + END_DOC + integer :: ipoint,istate + double precision :: weight,density,norm + mu_average_prov2 = 0.d0 + do istate = 1, N_states + norm = 0.d0 + do ipoint = 1, n_points_final_grid + weight =final_weight_at_r_vector(ipoint) + density = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate) & + + one_e_dm_and_grad_beta_in_r(4,ipoint,istate) + if(mu_of_r_prov(ipoint,istate).gt.1.d+09)cycle + mu_average_prov2(istate) += mu_of_r_prov(ipoint,istate) * weight * density*density + norm = norm + density*density*weight + enddo + mu_average_prov2(istate) = mu_average_prov2(istate) / norm + enddo +END_PROVIDER BEGIN_PROVIDER [double precision, mu_of_r_projector_mo, (n_points_final_grid) ]