From 70745cbeaaf59900c3ce4e1df042ef88ff1ecb11 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 23 May 2024 14:45:33 +0200 Subject: [PATCH] added sparse cholesky mu_of_r --- external/irpf90 | 2 +- .../basis_correction/basis_correction.irp.f | 4 -- .../basis_correction/print_routine.irp.f | 2 +- .../basis_correction/test_chol_bas.irp.f | 2 +- src/mu_of_r/f_hf_cholesky.irp.f | 52 ++++++++++++------- src/mu_of_r/mu_of_r_conditions.irp.f | 44 ++++++++++++++-- 6 files changed, 76 insertions(+), 30 deletions(-) diff --git a/external/irpf90 b/external/irpf90 index 0007f72f..4ab1b175 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102 +Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6 diff --git a/plugins/local/basis_correction/basis_correction.irp.f b/plugins/local/basis_correction/basis_correction.irp.f index a7ea7244..f17b5d5b 100644 --- a/plugins/local/basis_correction/basis_correction.irp.f +++ b/plugins/local/basis_correction/basis_correction.irp.f @@ -7,10 +7,6 @@ program basis_correction touch read_wf no_core_density = .True. touch no_core_density - if(io_mo_two_e_integrals .ne. "Read")then - provide ao_two_e_integrals_in_map - endif - provide mo_two_e_integrals_in_map call print_basis_correction end diff --git a/plugins/local/basis_correction/print_routine.irp.f b/plugins/local/basis_correction/print_routine.irp.f index 96faba30..b3b38673 100644 --- a/plugins/local/basis_correction/print_routine.irp.f +++ b/plugins/local/basis_correction/print_routine.irp.f @@ -22,7 +22,7 @@ subroutine print_basis_correction print*, '****************************************' print*, '****************************************' print*, 'mu_of_r_potential = ',mu_of_r_potential - if(mu_of_r_potential.EQ."hf")then + 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 ' diff --git a/plugins/local/basis_correction/test_chol_bas.irp.f b/plugins/local/basis_correction/test_chol_bas.irp.f index ae47ec09..076d888c 100644 --- a/plugins/local/basis_correction/test_chol_bas.irp.f +++ b/plugins/local/basis_correction/test_chol_bas.irp.f @@ -12,7 +12,7 @@ subroutine test do ipoint = 1, n_points_final_grid weight = final_weight_at_r_vector(ipoint) ! accu += dabs(mu_of_r_hf(ipoint) - mu_of_r_hf_old(ipoint)) * weight - accu += dabs(f_hf_sparse_cholesky(ipoint) - f_hf_cholesky(ipoint)) * weight + accu += dabs(f_hf_cholesky_sparse(ipoint) - f_hf_cholesky(ipoint)) * weight enddo print*,'accu = ',accu end diff --git a/src/mu_of_r/f_hf_cholesky.irp.f b/src/mu_of_r/f_hf_cholesky.irp.f index b937addf..17f0229a 100644 --- a/src/mu_of_r/f_hf_cholesky.irp.f +++ b/src/mu_of_r/f_hf_cholesky.irp.f @@ -189,7 +189,7 @@ BEGIN_PROVIDER [ double precision, f_hf_cholesky, (n_points_final_grid)] endif END_PROVIDER -BEGIN_PROVIDER [ double precision, f_hf_sparse_cholesky, (n_points_final_grid)] +BEGIN_PROVIDER [ double precision, f_hf_cholesky_sparse, (n_points_final_grid)] implicit none integer :: ipoint,m,mm,i,ii,p !!f(R) = \sum_{I} \sum_{J} Phi_I(R) Phi_J(R) V_IJ @@ -198,43 +198,55 @@ BEGIN_PROVIDER [ double precision, f_hf_sparse_cholesky, (n_points_final_grid)] !! = \sum_A V_AR G_AR !! V_AR = \sum_{I}Phi_IR V_AI = \sum_{I}Phi^t_RI V_AI double precision :: u_dot_v,wall0,wall1,accu_1, accu_2,mo_i_r1,mo_b_r1 + double precision :: thresh_1,thresh_2 + double precision, allocatable :: accu_vec(:) + thresh_2 = ao_cholesky_threshold * 100.d0 + thresh_1 = dsqrt(thresh_2) + provide cholesky_mo_transp if(elec_alpha_num == elec_beta_num)then call wall_time(wall0) - !$OMP PARALLEL DO & - !$OMP DEFAULT (NONE) & - !$OMP PRIVATE (accu_1,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) & - !$OMP ShARED (n_occ_val_orb_for_hf,list_valence_orb_for_hf,list_basis,mos_in_r_array_omp) & - !$OMP ShARED (cholesky_mo_num,f_hf_sparse_cholesky,n_points_final_grid,cholesky_mo,n_basis_orb) + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (accu_vec,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) & + !$OMP ShARED (n_occ_val_orb_for_hf,list_valence_orb_for_hf,list_basis,mos_in_r_array_omp,thresh_1,thresh_2) & + !$OMP ShARED (cholesky_mo_num,f_hf_cholesky_sparse,n_points_final_grid,cholesky_mo_transp,n_basis_orb) + allocate(accu_vec(cholesky_mo_num)) + !$OMP DO do ipoint = 1, n_points_final_grid - f_hf_sparse_cholesky(ipoint) = 0.d0 - do p = 1, cholesky_mo_num - accu_1 = 0.d0 + f_hf_cholesky_sparse(ipoint) = 0.d0 + accu_vec = 0.d0 do ii = 1, n_occ_val_orb_for_hf(1) i = list_valence_orb_for_hf(ii,1) mo_i_r1 = mos_in_r_array_omp(i,ipoint) + if(dabs(mo_i_r1).lt.thresh_1)cycle do mm = 1, n_basis_orb ! electron 1 m = list_basis(mm) mo_b_r1 = mos_in_r_array_omp(m,ipoint) - accu_1 += mo_i_r1 * mo_b_r1 * cholesky_mo(m,i,p) + if(dabs(mo_i_r1*mo_b_r1).lt.thresh_2)cycle + do p = 1, cholesky_mo_num + accu_vec(p) += mo_i_r1 * mo_b_r1 * cholesky_mo_transp(p,m,i) + enddo enddo enddo - f_hf_sparse_cholesky(ipoint) += accu_1 * accu_1 - enddo - f_hf_sparse_cholesky(ipoint) *= 2.D0 + do p = 1, cholesky_mo_num + f_hf_cholesky_sparse(ipoint) += accu_vec(p) * accu_vec(p) + enddo + f_hf_cholesky_sparse(ipoint) *= 2.D0 enddo - !$OMP END PARALLEL DO + !$OMP END DO + deallocate(accu_vec) + !$OMP END PARALLEL call wall_time(wall1) - print*,'Time to provide f_hf_sparse_cholesky = ',wall1-wall0 + print*,'Time to provide f_hf_cholesky_sparse = ',wall1-wall0 else call wall_time(wall0) !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (accu_2,accu_1,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) & !$OMP ShARED (n_occ_val_orb_for_hf,list_valence_orb_for_hf,list_basis,mos_in_r_array_omp) & - !$OMP ShARED (cholesky_mo_num,f_hf_sparse_cholesky,n_points_final_grid,cholesky_mo,n_basis_orb) + !$OMP ShARED (cholesky_mo_num,f_hf_cholesky_sparse,n_points_final_grid,cholesky_mo,n_basis_orb) do ipoint = 1, n_points_final_grid - f_hf_sparse_cholesky(ipoint) = 0.d0 + f_hf_cholesky_sparse(ipoint) = 0.d0 do p = 1, cholesky_mo_num accu_2 = 0.d0 do ii = 1, n_occ_val_orb_for_hf(2) @@ -256,13 +268,13 @@ BEGIN_PROVIDER [ double precision, f_hf_sparse_cholesky, (n_points_final_grid)] accu_1 += mo_i_r1 * mo_b_r1 * cholesky_mo(m,i,p) enddo enddo - f_hf_sparse_cholesky(ipoint) += accu_1 * accu_2 + f_hf_cholesky_sparse(ipoint) += accu_1 * accu_2 enddo - f_hf_sparse_cholesky(ipoint) *= 2.D0 + f_hf_cholesky_sparse(ipoint) *= 2.D0 enddo !$OMP END PARALLEL DO call wall_time(wall1) - print*,'Time to provide f_hf_sparse_cholesky = ',wall1-wall0 + print*,'Time to provide f_hf_cholesky_sparse = ',wall1-wall0 endif 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 5b4d4b83..f2bb7145 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -13,7 +13,6 @@ integer :: ipoint,istate double precision :: wall0,wall1 print*,'providing mu_of_r ...' -! PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals call wall_time(wall0) if (read_mu_of_r) then @@ -26,6 +25,10 @@ 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) + else if(mu_of_r_potential.EQ."hf_old")then + mu_of_r_prov(ipoint,istate) = mu_of_r_hf_old(ipoint) + else if(mu_of_r_potential.EQ."hf_sparse")then + mu_of_r_prov(ipoint,istate) = mu_of_r_hf_sparse(ipoint) 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) else @@ -61,11 +64,10 @@ END_DOC integer :: ipoint double precision :: wall0,wall1,f_hf,on_top,w_hf,sqpi - PROVIDE f_hf_cholesky on_top_hf_grid print*,'providing mu_of_r_hf ...' call wall_time(wall0) + PROVIDE f_hf_cholesky on_top_hf_grid sqpi = dsqrt(dacos(-1.d0)) - provide f_psi_hf_ab !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & @@ -85,6 +87,42 @@ print*,'Time to provide mu_of_r_hf = ',wall1-wall0 END_PROVIDER + BEGIN_PROVIDER [double precision, mu_of_r_hf_sparse, (n_points_final_grid) ] + implicit none + BEGIN_DOC + ! mu(r) computed with a HF wave function (assumes that HF MOs are stored in the EZFIO) + ! + ! corresponds to Eq. (37) of J. Chem. Phys. 149, 194301 (2018) but for \Psi^B = HF^B + ! + ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals + ! + ! in the two-body density matrix are excluded + END_DOC + integer :: ipoint + double precision :: wall0,wall1,f_hf,on_top,w_hf,sqpi + print*,'providing mu_of_r_hf_sparse ...' + call wall_time(wall0) + sqpi = dsqrt(dacos(-1.d0)) + PROVIDE f_hf_cholesky_sparse on_top_hf_grid + !$OMP PARALLEL DO & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (ipoint,f_hf,on_top,w_hf) & + !$OMP ShARED (n_points_final_grid,mu_of_r_hf_sparse,f_hf_cholesky_sparse,on_top_hf_grid,sqpi) + do ipoint = 1, n_points_final_grid + f_hf = f_hf_cholesky_sparse(ipoint) + on_top = on_top_hf_grid(ipoint) + if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then + w_hf = 1.d+10 + else + w_hf = f_hf / on_top + endif + mu_of_r_hf_sparse(ipoint) = w_hf * sqpi * 0.5d0 + enddo + !$OMP END PARALLEL DO + call wall_time(wall1) + print*,'Time to provide mu_of_r_hf_sparse = ',wall1-wall0 + END_PROVIDER + BEGIN_PROVIDER [double precision, mu_of_r_hf_old, (n_points_final_grid) ] implicit none BEGIN_DOC