From ed1253f62956c10866ea83c707cf532b7f3d4d67 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 25 Feb 2025 18:09:57 +0100 Subject: [PATCH] Sort indices for faster access in RDM --- external/ezfio | 2 +- .../basis_correction/print_routine.irp.f | 97 ++++++++------- src/ao_two_e_ints/two_e_integrals.irp.f | 1 + src/cas_based_on_top/on_top_cas_prov.irp.f | 13 +- src/dft_utils_in_r/ao_in_r.irp.f | 112 ++++++++---------- src/dft_utils_in_r/mo_in_r.irp.f | 15 +-- src/mu_of_r/mu_of_r_conditions.irp.f | 85 ++++++++++++- src/two_rdm_routines/davidson_like_2rdm.irp.f | 88 ++++++++------ 8 files changed, 248 insertions(+), 165 deletions(-) diff --git a/external/ezfio b/external/ezfio index dba01c4f..d02132ea 160000 --- a/external/ezfio +++ b/external/ezfio @@ -1 +1 @@ -Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3 +Subproject commit d02132ea79217c16fd24242e8f8b8a6c3ff68091 diff --git a/plugins/local/basis_correction/print_routine.irp.f b/plugins/local/basis_correction/print_routine.irp.f index b3b38673..8879fd5d 100644 --- a/plugins/local/basis_correction/print_routine.irp.f +++ b/plugins/local/basis_correction/print_routine.irp.f @@ -22,53 +22,58 @@ subroutine print_basis_correction print*, '****************************************' print*, '****************************************' print*, 'mu_of_r_potential = ',mu_of_r_potential - if(mu_of_r_potential.EQ."hf".or.mu_of_r_potential.EQ."hf_old".or.mu_of_r_potential.EQ."hf_sparse")then - print*, '' - print*,'Using a HF-like two-body density to define mu(r)' - print*,'This assumes that HF is a qualitative representation of the wave function ' - print*,'********************************************' - print*,'Functionals more suited for weak correlation' - print*,'********************************************' - print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) ' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate) - enddo - print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) - enddo + if(mu_of_r_potential.EQ."hf".or. & + mu_of_r_potential.EQ."hf_old".or.& + mu_of_r_potential.EQ."hf_sparse".or.& + mu_of_r_potential.EQ."proj")then + print*, '' + print*,'Using a HF-like two-body density to define mu(r)' + print*,'This assumes that HF is a qualitative representation of the wave function ' + print*,'********************************************' + print*,'Functionals more suited for weak correlation' + print*,'********************************************' + print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) ' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate) + enddo + print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) + enddo - else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then - print*, '' - print*,'Using a CAS-like two-body density to define mu(r)' - print*,'This assumes that the CAS is a qualitative representation of the wave function ' - print*,'********************************************' - print*,'Functionals more suited for weak correlation' - print*,'********************************************' - print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) ' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate) - enddo - print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))' - do istate = 1, N_states - write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) - enddo - print*,'' - print*,'********************************************' - print*,'********************************************' - 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 : 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) - enddo - print*,'' + else if(mu_of_r_potential.EQ."cas_full".or. & + mu_of_r_potential.EQ."cas_truncated".or. & + mu_of_r_potential.EQ."pure_act") then + print*, '' + print*,'Using a CAS-like two-body density to define mu(r)' + print*,'This assumes that the CAS is a qualitative representation of the wave function ' + print*,'********************************************' + print*,'Functionals more suited for weak correlation' + print*,'********************************************' + print*,'+) LDA Ecmd functional : purely based on the UEG (JCP,149,194301,1-15 (2018)) ' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD LDA , state ',istate,' = ',ecmd_lda_mu_of_r(istate) + enddo + print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))' + do istate = 1, N_states + write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) + enddo + print*,'' + print*,'********************************************' + print*,'********************************************' + 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 : 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) + enddo + print*,'' endif print*,'' diff --git a/src/ao_two_e_ints/two_e_integrals.irp.f b/src/ao_two_e_ints/two_e_integrals.irp.f index fb376ce1..1cb7617e 100644 --- a/src/ao_two_e_ints/two_e_integrals.irp.f +++ b/src/ao_two_e_ints/two_e_integrals.irp.f @@ -54,6 +54,7 @@ double precision function ao_two_e_integral(i, j, k, l) else if (use_only_lr) then ao_two_e_integral = ao_two_e_integral_erf(i, j, k, l) + return else if (do_schwartz_accel(i,j,k,l)) then diff --git a/src/cas_based_on_top/on_top_cas_prov.irp.f b/src/cas_based_on_top/on_top_cas_prov.irp.f index dd93ed40..9f9a1f06 100644 --- a/src/cas_based_on_top/on_top_cas_prov.irp.f +++ b/src/cas_based_on_top/on_top_cas_prov.irp.f @@ -15,14 +15,17 @@ pure_act_on_top_of_r = 0.d0 do l = 1, n_act_orb phi_l = act_mos_in_r_array(l,ipoint) + if (dabs(phi_l) < 1.d-12) cycle do k = 1, n_act_orb - phi_k = act_mos_in_r_array(k,ipoint) + phi_k = act_mos_in_r_array(k,ipoint) * phi_l + if (dabs(phi_k) < 1.d-12) cycle do j = 1, n_act_orb - phi_j = act_mos_in_r_array(j,ipoint) + phi_j = act_mos_in_r_array(j,ipoint) * phi_k + if (dabs(phi_j) < 1.d-12) cycle do i = 1, n_act_orb - phi_i = act_mos_in_r_array(i,ipoint) - ! 1 2 1 2 - pure_act_on_top_of_r += act_2_rdm_ab_mo(i,j,k,l,istate) * phi_i * phi_j * phi_k * phi_l + phi_i = act_mos_in_r_array(i,ipoint) * phi_j + ! 1 2 1 2 + pure_act_on_top_of_r = pure_act_on_top_of_r + act_2_rdm_ab_mo(i,j,k,l,istate) * phi_i !* phi_j * phi_k * phi_l enddo enddo enddo diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index c8822776..ffd8c9d5 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -8,21 +8,14 @@ BEGIN_PROVIDER[double precision, aos_in_r_array, (ao_num,n_points_final_grid)] END_DOC implicit none - integer :: i, j - double precision :: tmp_array(ao_num), r(3) + integer :: i !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,tmp_array,j) & - !$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points) + !$OMP PRIVATE (i) & + !$OMP SHARED(aos_in_r_array,n_points_final_grid,final_grid_points) do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_aos_at_r(r, tmp_array) - do j = 1, ao_num - aos_in_r_array(j,i) = tmp_array(j) - enddo + call give_all_aos_at_r(final_grid_points(1,i), aos_in_r_array(1,i)) enddo !$OMP END PARALLEL DO @@ -42,7 +35,7 @@ BEGIN_PROVIDER[double precision, aos_in_r_array_transp, (n_points_final_grid,ao_ do i = 1, n_points_final_grid do j = 1, ao_num - aos_in_r_array_transp(i,j) = aos_in_r_array(j,i) + aos_in_r_array_transp(i,j) = aos_in_r_array(j,i) enddo enddo @@ -62,27 +55,29 @@ BEGIN_PROVIDER[double precision, aos_grad_in_r_array, (ao_num,n_points_final_gri implicit none integer :: i, j, m - double precision :: aos_array(ao_num), r(3) - double precision :: aos_grad_array(3,ao_num) + double precision :: r(3) + double precision, allocatable :: aos_grad_array(:,:), aos_array(:) - !$OMP PARALLEL DO & + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) & + !$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) & !$OMP SHARED(aos_grad_in_r_array,n_points_final_grid,ao_num,final_grid_points) + allocate(aos_grad_array(3,ao_num), aos_array(ao_num)) + + !$OMP DO do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) + call give_all_aos_and_grad_at_r(final_grid_points(1,i),aos_array,aos_grad_array) do m = 1, 3 do j = 1, ao_num aos_grad_in_r_array(j,i,m) = aos_grad_array(m,j) enddo enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO + deallocate(aos_grad_array,aos_array) + !$OMP END PARALLEL -END_PROVIDER +END_PROVIDER ! --- @@ -116,7 +111,7 @@ END_PROVIDER enddo enddo enddo - END_PROVIDER + END_PROVIDER BEGIN_PROVIDER [double precision, aos_lapl_in_r_array, (3,ao_num,n_points_final_grid)] implicit none @@ -126,32 +121,32 @@ END_PROVIDER ! k = 1 : x, k= 2, y, k 3, z END_DOC integer :: i,j,m - double precision :: aos_array(ao_num), r(3) - double precision :: aos_grad_array(3,ao_num) - double precision :: aos_lapl_array(3,ao_num) - !$OMP PARALLEL DO & + double precision, allocatable :: aos_lapl_array(:,:), aos_grad_array(:,:), aos_array(:) + + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,aos_array,aos_grad_array,aos_lapl_array,j,m) & + !$OMP PRIVATE (i,aos_array,aos_grad_array,aos_lapl_array,j,m) & !$OMP SHARED(aos_lapl_in_r_array,n_points_final_grid,ao_num,final_grid_points) + allocate( aos_array(ao_num), aos_grad_array(3,ao_num), aos_lapl_array(3,ao_num)) + !$OMP DO do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array) + call give_all_aos_and_grad_and_lapl_at_r(final_grid_points(1,i),aos_array,aos_grad_array,aos_lapl_array) do j = 1, ao_num do m = 1, 3 aos_lapl_in_r_array(m,j,i) = aos_lapl_array(m,j) enddo enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO + deallocate( aos_array, aos_grad_array, aos_lapl_array) + !$OMP END PARALLEL END_PROVIDER BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_bis, (n_points_final_grid,ao_num,3)] implicit none BEGIN_DOC -! Transposed gradients -! +! Transposed gradients +! END_DOC integer :: i,j,m double precision :: aos_array(ao_num), r(3) @@ -169,8 +164,8 @@ END_PROVIDER BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_3, (3,n_points_final_grid,ao_num)] implicit none BEGIN_DOC -! Transposed gradients -! +! Transposed gradients +! END_DOC integer :: i,j,m double precision :: aos_array(ao_num), r(3) @@ -187,22 +182,14 @@ END_PROVIDER BEGIN_PROVIDER[double precision, aos_in_r_array_extra, (ao_num,n_points_extra_final_grid)] implicit none BEGIN_DOC - ! aos_in_r_array_extra(i,j) = value of the ith ao on the jth grid point of the EXTRA grid + ! aos_in_r_array_extra(i,j) = value of the ith ao on the jth grid point of the EXTRA grid END_DOC - integer :: i,j - double precision :: aos_array(ao_num), r(3) + integer :: i !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,aos_array,j) & - !$OMP SHARED(aos_in_r_array_extra,n_points_extra_final_grid,ao_num,final_grid_points_extra) + !$OMP DEFAULT (NONE) PRIVATE (i) & + !$OMP SHARED(aos_in_r_array_extra,n_points_extra_final_grid,final_grid_points_extra) do i = 1, n_points_extra_final_grid - r(1) = final_grid_points_extra(1,i) - r(2) = final_grid_points_extra(2,i) - r(3) = final_grid_points_extra(3,i) - call give_all_aos_at_r(r,aos_array) - do j = 1, ao_num - aos_in_r_array_extra(j,i) = aos_array(j) - enddo + call give_all_aos_at_r(final_grid_points_extra(1,i),aos_in_r_array_extra(1,i)) enddo !$OMP END PARALLEL DO @@ -214,9 +201,9 @@ END_PROVIDER BEGIN_PROVIDER[double precision, aos_in_r_array_extra_transp, (n_points_extra_final_grid,ao_num)] BEGIN_DOC - ! aos_in_r_array_extra_transp(i,j) = value of the jth ao on the ith grid point of the EXTRA grid + ! aos_in_r_array_extra_transp(i,j) = value of the jth ao on the ith grid point of the EXTRA grid END_DOC - + implicit none integer :: i, j double precision :: aos_array(ao_num), r(3) @@ -235,27 +222,28 @@ BEGIN_PROVIDER[double precision, aos_grad_in_r_array_extra, (ao_num,n_points_ext implicit none integer :: i, j, m - double precision :: aos_array(ao_num), r(3) - double precision :: aos_grad_array(3,ao_num) + double precision, allocatable :: aos_array(:), aos_grad_array(:,:) - !$OMP PARALLEL DO & + + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,j,m,r,aos_array,aos_grad_array) & + !$OMP PRIVATE (i,j,m,aos_array,aos_grad_array) & !$OMP SHARED(aos_grad_in_r_array_extra,n_points_extra_final_grid,ao_num,final_grid_points_extra) + allocate(aos_array(ao_num), aos_grad_array(3,ao_num)) + !$OMP DO do i = 1, n_points_extra_final_grid - r(1) = final_grid_points_extra(1,i) - r(2) = final_grid_points_extra(2,i) - r(3) = final_grid_points_extra(3,i) - call give_all_aos_and_grad_at_r(r, aos_array, aos_grad_array) + call give_all_aos_and_grad_at_r(final_grid_points_extra(1,i), aos_array, aos_grad_array) do m = 1, 3 do j = 1, ao_num aos_grad_in_r_array_extra(j,i,m) = aos_grad_array(m,j) enddo enddo enddo - !$OMP END PARALLEL DO + !$OMP END DO + deallocate(aos_array,aos_grad_array) + !$OMP END PARALLEL -END_PROVIDER +END_PROVIDER ! --- 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 623de4f8..f3cc5559 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -21,20 +21,11 @@ BEGIN_DOC ! mos_in_r_array(i,j) = value of the ith mo on the jth grid point END_DOC - integer :: i,j - double precision :: mos_array(mo_num), r(3) - !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (i,r,mos_array,j) & + integer :: i + !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE (i) & !$OMP SHARED(mos_in_r_array_omp,n_points_final_grid,mo_num,final_grid_points) do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_mos_at_r(r,mos_array) - do j = 1, mo_num - mos_in_r_array_omp(j,i) = mos_array(j) - enddo + call give_all_mos_at_r(final_grid_points(1,i),mos_in_r_array_omp(1,i)) enddo !$OMP END PARALLEL DO END_PROVIDER 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 88dad8c3..3ad52a05 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -22,22 +22,32 @@ endif do istate = 1, N_states - do ipoint = 1, n_points_final_grid if(mu_of_r_potential.EQ."hf")then - mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) + enddo else if(mu_of_r_potential.EQ."hf_old")then - mu_of_r_prov(ipoint,istate) = mu_of_r_hf_old(ipoint) + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_hf_old(ipoint) + enddo else if(mu_of_r_potential.EQ."hf_sparse")then - mu_of_r_prov(ipoint,istate) = mu_of_r_hf_sparse(ipoint) + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_hf_sparse(ipoint) + enddo else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then - mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) + enddo + else if(mu_of_r_potential.EQ."proj")then + do ipoint = 1, n_points_final_grid + mu_of_r_prov(ipoint,istate) = mu_of_r_projector_mo(ipoint) + enddo else print*,'you requested the following mu_of_r_potential' print*,mu_of_r_potential print*,'which does not correspond to any of the options for such keyword' stop endif - enddo enddo if (write_mu_of_r) then @@ -225,3 +235,66 @@ enddo END_PROVIDER + +BEGIN_PROVIDER [double precision, mu_of_r_projector_mo, (n_points_final_grid) ] + implicit none + BEGIN_DOC + ! mu(r) computed with the projector onto the atomic basis + ! P_B(\mathbf{r},\mathbf{r}') = \sum_{ij} | + ! \chi_{i} \rangle \left[S^{-1}\right]_{ij} \langle \chi_{j} | + ! \] where $i$ and $j$ denote all atomic orbitals. + END_DOC + + double precision, parameter :: factor = dsqrt(2.d0*dacos(-1.d0)) + double precision, allocatable :: tmp(:,:) + integer :: ipoint + + + do ipoint=1,n_points_final_grid + mu_of_r_projector_mo(ipoint) = 0.d0 + integer :: i,j + do j=1,n_inact_act_orb + i = list_inact_act(j) + mu_of_r_projector_mo(ipoint) = mu_of_r_projector_mo(ipoint) + & + mos_in_r_array_omp(i,ipoint) * mos_in_r_array_omp(i,ipoint) + enddo + do j=1,n_virt_orb + i = list_virt(j) + mu_of_r_projector_mo(ipoint) = mu_of_r_projector_mo(ipoint) + & + mos_in_r_array_omp(i,ipoint) * mos_in_r_array_omp(i,ipoint) + enddo + enddo + + do ipoint=1,n_points_final_grid + ! epsilon + mu_of_r_projector_mo(ipoint) = 1.d0/(2.d0*dacos(-1.d0) * mu_of_r_projector_mo(ipoint)**(2.d0/3.d0)) + ! mu + mu_of_r_projector_mo(ipoint) = 1.d0/dsqrt( mu_of_r_projector_mo(ipoint) ) + enddo +END_PROVIDER + + + +BEGIN_PROVIDER [double precision, mu_average_proj, (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 + ! + ! !!!!!! 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 + do istate = 1, N_states + mu_average_proj(istate) = 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) + mu_average_proj(istate) += mu_of_r_projector_mo(ipoint) * weight * density + enddo + mu_average_proj(istate) = mu_average_proj(istate) / elec_num_grid_becke(istate) + enddo +END_PROVIDER + diff --git a/src/two_rdm_routines/davidson_like_2rdm.irp.f b/src/two_rdm_routines/davidson_like_2rdm.irp.f index 09436663..f0b40459 100644 --- a/src/two_rdm_routines/davidson_like_2rdm.irp.f +++ b/src/two_rdm_routines/davidson_like_2rdm.irp.f @@ -139,6 +139,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) sze_buff = 6 * norb + elec_alpha_num * elec_alpha_num * 60 + sze_buff = sze_buff*100 list_orb_reverse = -1000 do i = 1, norb list_orb_reverse(list_orb(i)) = i @@ -191,7 +192,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin ASSERT (istart > 0) ASSERT (istep > 0) - !$OMP DO SCHEDULE(dynamic,64) + !$OMP DO SCHEDULE(dynamic) do k_a=istart+ishift,iend,istep krow = psi_bilinear_matrix_rows(k_a) @@ -272,14 +273,14 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin enddo endif - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 enddo +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 enddo !$OMP END DO - !$OMP DO SCHEDULE(dynamic,64) + !$OMP DO SCHEDULE(dynamic) do k_a=istart+ishift,iend,istep @@ -343,24 +344,24 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin c_1(l) = u_t(l,l_a) * u_t(l,k_a) enddo if(alpha_beta.or.spin_trace.or.alpha_alpha)then - ! increment the alpha/beta part for single excitations - if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - ! increment the alpha/alpha part for single excitations - if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 - endif - call orb_range_off_diag_single_to_all_states_aa_dm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + ! increment the alpha/beta part for single excitations + if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then + call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_ab_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + ! increment the alpha/alpha part for single excitations + if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then + call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) + nkeys = 0 + endif + call orb_range_off_diag_single_to_all_states_aa_dm_buffer(tmp_det,tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) endif enddo - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 ! Compute Hij for all alpha doubles ! ---------------------------------- @@ -383,8 +384,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin call orb_range_off_diag_double_to_all_states_aa_dm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) enddo endif - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 ! Single and double beta excitations @@ -459,8 +460,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin call orb_range_off_diag_single_to_all_states_bb_dm_buffer(tmp_det, tmp_det2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) endif enddo - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 ! Compute Hij for all beta doubles ! ---------------------------------- @@ -487,8 +488,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin enddo endif - call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) - nkeys = 0 +! call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm) +! nkeys = 0 ! Diagonal contribution @@ -550,22 +551,43 @@ subroutine update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,loc integer :: istate integer :: i,h1,h2,p1,p2 - call omp_set_lock(lock_2rdm) + integer, allocatable :: iorder(:) + integer*8, allocatable :: to_sort(:) + + allocate(iorder(nkeys)) + do i=1,nkeys + iorder(i) = i + enddo + + ! If the lock is already taken, sort the keys while waiting for a faster access + if (.not.omp_test_lock(lock_2rdm)) then + allocate(to_sort(nkeys)) + do i=1,nkeys + h1 = keys(1,iorder(i)) + h2 = keys(2,iorder(i))-1 + p1 = keys(3,iorder(i))-1 + p2 = keys(4,iorder(i))-1 + to_sort(i) = int(h1,8) + int(dim1,8)*(int(h2,8) + int(dim1,8)*(int(p1,8) + int(dim1,8)*int(p2,8))) + enddo + call i8sort(to_sort, iorder, nkeys) + deallocate(to_sort) + call omp_set_lock(lock_2rdm) + endif ! print*,'*************' ! print*,'updating' ! print*,'nkeys',nkeys - do i = 1, nkeys - h1 = keys(1,i) - h2 = keys(2,i) - p1 = keys(3,i) - p2 = keys(4,i) - do istate = 1, N_st -! print*,h1,h2,p1,p2,values(istate,i) - big_array(h1,h2,p1,p2,istate) += values(istate,i) + do istate = 1, N_st + do i = 1, nkeys + h1 = keys(1,iorder(i)) + h2 = keys(2,iorder(i)) + p1 = keys(3,iorder(i)) + p2 = keys(4,iorder(i)) + big_array(h1,h2,p1,p2,istate) = big_array(h1,h2,p1,p2,istate) + values(istate,iorder(i)) enddo enddo call omp_unset_lock(lock_2rdm) + deallocate(iorder) end