From ee267e27e9b6027ea76a05e2b43271cc8906fd55 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Mon, 21 Sep 2020 15:38:26 +0200 Subject: [PATCH] minor modifs --- config/gfortran_debug.cfg | 3 +- src/basis_correction/print_routine.irp.f | 1 + src/dft_utils_in_r/mo_in_r.irp.f | 3 +- src/mo_basis/mos_in_r.irp.f | 10 ++-- src/mo_guess/h_core_guess_routine.irp.f | 4 +- src/mu_of_r/mu_of_r_conditions.irp.f | 67 ++++++++++++++++++++++++ src/utils/util.irp.f | 12 +++++ 7 files changed, 91 insertions(+), 9 deletions(-) diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg index 342acae9..926255e0 100644 --- a/config/gfortran_debug.cfg +++ b/config/gfortran_debug.cfg @@ -51,7 +51,8 @@ FCFLAGS : -Ofast # -g : Extra debugging information # [DEBUG] -FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan +#FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan +FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow -finit-real=nan # OpenMP flags ################# diff --git a/src/basis_correction/print_routine.irp.f b/src/basis_correction/print_routine.irp.f index 05fbbf60..3d497552 100644 --- a/src/basis_correction/print_routine.irp.f +++ b/src/basis_correction/print_routine.irp.f @@ -75,6 +75,7 @@ subroutine print_basis_correction 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)') 'mu_average_trans_corr, state ',istate,' = ',mu_average_trans_corr(istate) enddo end diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index 81863c3a..fa3ec73c 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -79,7 +79,7 @@ END_DOC integer :: m integer :: i,j - mos_grad_in_r_array = 0.d0 + mos_grad_in_r_array_tranp = 0.d0 do i = 1, n_points_final_grid do j = 1, mo_num do m = 1, 3 @@ -126,4 +126,3 @@ enddo END_PROVIDER - diff --git a/src/mo_basis/mos_in_r.irp.f b/src/mo_basis/mos_in_r.irp.f index 049db8aa..ffcd5a9a 100644 --- a/src/mo_basis/mos_in_r.irp.f +++ b/src/mo_basis/mos_in_r.irp.f @@ -12,18 +12,18 @@ subroutine give_all_mos_and_grad_at_r(r,mos_array,mos_grad_array) implicit none double precision, intent(in) :: r(3) double precision, intent(out) :: mos_array(mo_num) - double precision, intent(out) :: mos_grad_array(mo_num,3) + double precision, intent(out) :: mos_grad_array(3,mo_num) integer :: i,j,k - double precision :: aos_array(ao_num),aos_grad_array(ao_num,3) + double precision :: aos_array(ao_num),aos_grad_array(3,ao_num) call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) mos_array=0d0 mos_grad_array=0d0 do j = 1, mo_num do k=1, ao_num mos_array(j) += mo_coef(k,j)*aos_array(k) - mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1) - mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2) - mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3) + mos_grad_array(1,j) += mo_coef(k,j)*aos_grad_array(1,k) + mos_grad_array(2,j) += mo_coef(k,j)*aos_grad_array(2,k) + mos_grad_array(3,j) += mo_coef(k,j)*aos_grad_array(3,k) enddo enddo end diff --git a/src/mo_guess/h_core_guess_routine.irp.f b/src/mo_guess/h_core_guess_routine.irp.f index 8fc3f6f2..9c5e76ff 100644 --- a/src/mo_guess/h_core_guess_routine.irp.f +++ b/src/mo_guess/h_core_guess_routine.irp.f @@ -9,5 +9,7 @@ subroutine hcore_guess size(mo_one_e_integrals,1), & size(mo_one_e_integrals,2),label,1,.false.) call save_mos - SOFT_TOUCH mo_coef mo_label +! SOFT_TOUCH mo_coef mo_label + TOUCH mo_coef mo_label + print*,'mo_one_e_integrals(1,1) = ',mo_one_e_integrals(1,1) 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 05a2a4d7..b35a2a8a 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -28,6 +28,8 @@ mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) + else if(mu_of_r_potential.EQ."transcorr")then + mu_of_r_prov(ipoint,istate) = mu_of_r_transcorr(ipoint,istate) else print*,'you requested the following mu_of_r_potential' print*,mu_of_r_potential @@ -124,6 +126,47 @@ print*,'Time to provide mu_of_r_psi_cas = ',wall1-wall0 END_PROVIDER + BEGIN_PROVIDER [double precision, mu_of_r_transcorr, (n_points_final_grid,N_states) ] + implicit none + BEGIN_DOC + ! mu(r) computed with a wave function developped in an active space + ! + ! corresponds to \sqrt(pi) * (W(0) + 1/4)/3 + ! + ! !!!!!! 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 :: wall0,wall1,f_psi,on_top,w_psi,sqpi + print*,'providing mu_of_r_transcorr ...' + call wall_time(wall0) + sqpi = dsqrt(dacos(-1.d0)) + + provide f_psi_cas_ab + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) & + !$OMP SHARED (n_points_final_grid,mu_of_r_transcorr,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states) + 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) + 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 + w_psi = f_psi / on_top + endif + mu_of_r_transcorr(ipoint,istate) = (0.25d0 + w_psi) * sqpi / 3.d0 + enddo + enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide mu_of_r_transcorr = ',wall1-wall0 + END_PROVIDER + BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)] implicit none @@ -148,3 +191,27 @@ mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate) enddo END_PROVIDER + + BEGIN_PROVIDER [double precision, mu_average_trans_corr, (N_states)] + implicit none + BEGIN_DOC + ! average value of mu(r) weighted with the total one-e density and divised by the number of electrons + ! + ! !!!!!! 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 + mu_average_trans_corr = 0.d0 + do istate = 1, N_states + 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_transcorr(ipoint,istate).gt.1.d+09)cycle + mu_average_trans_corr(istate) += mu_of_r_transcorr(ipoint,istate) * weight * density + enddo + mu_average_trans_corr(istate) = mu_average_trans_corr(istate) / elec_num_grid_becke(istate) + enddo + END_PROVIDER diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f index 1b01a1ec..cfb42fd1 100644 --- a/src/utils/util.irp.f +++ b/src/utils/util.irp.f @@ -1,3 +1,15 @@ +double precision function derf_mu_x(mu,x) + implicit none + include 'utils/constants.include.F' + double precision, intent(in) :: mu,x + if(dabs(x).gt.1.d-6)then + derf_mu_x = derf(mu * x)/x + else + derf_mu_x = inv_sq_pi * 2.d0 * mu * (1.d0 - mu*mu*x*x/3.d0) + endif +end + + double precision function binom_func(i,j) implicit none BEGIN_DOC