LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2025-04-25 17:54:44 +02:00
Code Releases Activity

Sort indices for faster access in RDM

Browse Source
This commit is contained in:
Anthony Scemama 2025-02-25 18:09:57 +01:00
parent d20ac13c4f
commit ed1253f629
8 changed files with 248 additions and 165 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
external/ezfio vendored

@ -1 +1 @@
Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3
Subproject commit d02132ea79217c16fd24242e8f8b8a6c3ff68091

9
plugins/local/basis_correction/print_routine.irp.f
View File

@ -22,7 +22,10 @@ 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
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 '
@ -38,7 +41,9 @@ subroutine print_basis_correction
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
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 '

1
src/ao_two_e_ints/two_e_integrals.irp.f
View File

@ -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

11
src/cas_based_on_top/on_top_cas_prov.irp.f
View File

@ -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)
phi_i = act_mos_in_r_array(i,ipoint) * phi_j
! 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
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

88
src/dft_utils_in_r/ao_in_r.irp.f
View File

@ -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
@ -62,25 +55,27 @@ 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 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
@ -126,25 +121,25 @@ 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)]
@ -189,20 +184,12 @@ END_PROVIDER
BEGIN_DOC
! 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
@ -235,25 +222,26 @@ 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

15
src/dft_utils_in_r/mo_in_r.irp.f
View File

@ -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

77
src/mu_of_r/mu_of_r_conditions.irp.f
View File

@ -22,15 +22,26 @@
endif
do istate = 1, N_states
do ipoint = 1, n_points_final_grid
if(mu_of_r_potential.EQ."hf")then
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
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
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
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
@ -38,7 +49,6 @@
stop
endif
enddo
enddo
if (write_mu_of_r) then
print*,'Writing mu(r) on disk ...'
@ -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

60
src/two_rdm_routines/davidson_like_2rdm.irp.f
View File

@ -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
@ -359,8 +360,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
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
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 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