Merge branch 'dev-stable' of github.com:QuantumPackage/qp2 into dev-stable

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
+ if(mu_of_r_potential.EQ."hf".or. &
-   print*, ''
+    mu_of_r_potential.EQ."hf_old".or.&
-   print*,'Using a HF-like two-body density to define mu(r)'
+    mu_of_r_potential.EQ."hf_sparse".or.&
-   print*,'This assumes that HF is a qualitative representation of the wave function '
+    mu_of_r_potential.EQ."proj")then
-   print*,'********************************************'
+     print*, ''
-   print*,'Functionals more suited for weak correlation'
+     print*,'Using a HF-like two-body density to define mu(r)'
-   print*,'********************************************'
+     print*,'This assumes that HF is a qualitative representation of the wave function '
-   print*,'+) LDA Ecmd functional     : purely based on the UEG (JCP,149,194301,1-15 (2018)) '
+     print*,'********************************************'
-   do istate = 1, N_states
+     print*,'Functionals more suited for weak correlation'
-    write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD LDA           , state ',istate,' = ',ecmd_lda_mu_of_r(istate)
+     print*,'********************************************'
-   enddo
+     print*,'+) LDA Ecmd functional     : purely based on the UEG (JCP,149,194301,1-15 (2018)) '
-   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
-   do istate = 1, N_states
+      write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD LDA           , state ',istate,' = ',ecmd_lda_mu_of_r(istate)
-    write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-UEG       , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
+     enddo
-   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
+  else if(mu_of_r_potential.EQ."cas_full".or. &
-   print*, ''
+          mu_of_r_potential.EQ."cas_truncated".or. &
-   print*,'Using a CAS-like two-body density to define mu(r)'
+          mu_of_r_potential.EQ."pure_act") then
-   print*,'This assumes that the CAS is a qualitative representation of the wave function '
+     print*, ''
-   print*,'********************************************'
+     print*,'Using a CAS-like two-body density to define mu(r)'
-   print*,'Functionals more suited for weak correlation'
+     print*,'This assumes that the CAS is a qualitative representation of the wave function '
-   print*,'********************************************'
+     print*,'********************************************'
-   print*,'+) LDA Ecmd functional     : purely based on the UEG (JCP,149,194301,1-15 (2018)) '
+     print*,'Functionals more suited for weak correlation'
-   do istate = 1, N_states
+     print*,'********************************************'
-    write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD LDA           , state ',istate,' = ',ecmd_lda_mu_of_r(istate)
+     print*,'+) LDA Ecmd functional     : purely based on the UEG (JCP,149,194301,1-15 (2018)) '
-   enddo
+     do istate = 1, N_states
-   print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))'
+      write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD LDA           , state ',istate,' = ',ecmd_lda_mu_of_r(istate)
-   do istate = 1, N_states
+     enddo
-    write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-UEG       , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
+     print*,'+) PBE-UEG Ecmd functional : PBE at mu=0, UEG ontop pair density at large mu (JPCL, 10, 2931-2937 (2019))'
-   enddo
+     do istate = 1, N_states
-   print*,''
+      write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-UEG       , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
-   print*,'********************************************'
+     enddo
-   print*,'********************************************'
+     print*,''
-   print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) '
+     print*,'********************************************'
-   print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization'
+     print*,'********************************************'
-   do istate = 1, N_states
+     print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) '
-    write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-OT        , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate)
+     print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization'
-   enddo
+     do istate = 1, N_states
-   print*,''
+      write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-OT        , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate)
-   print*,'********************************************'
+     enddo
-   print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)'
+     print*,''
-   print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION'
+     print*,'********************************************'
-   do istate = 1, N_states
+     print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)'
-    write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD SU-PBE-OT     , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)
+     print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION'
-   enddo
+     do istate = 1, N_states
-   print*,''
+      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*,''

src/ao_two_e_ints/cholesky.irp.f

@@ -178,7 +178,7 @@ END_PROVIDER
      rank_max = np
      ! Avoid too large arrays when there are many electrons
      if (elec_num > 10) then
-       rank_max = min(np,20*elec_num*elec_num)
+       rank_max = min(np,25*elec_num*elec_num)
      endif
 
      call mmap_create_d('', (/ ndim8, rank_max /), .False., .True., map)

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
 

src/bitmask/bitmasks_routines.irp.f

@@ -283,33 +283,16 @@ subroutine print_det_one_dimension(string,Nint)
 
 end
 
-logical function is_integer_in_string(bite,string,Nint)
+logical function is_integer_in_string(orb,bitmask,Nint)
-  use bitmasks
+ use bitmasks
  implicit none
- integer, intent(in) :: bite,Nint
+ BEGIN_DOC
- integer(bit_kind), intent(in) :: string(Nint)
+! Checks is the orbital orb is set to 1 in the bit string
- integer(bit_kind) :: string_bite(Nint)
+ END_DOC
- integer   :: i,itot,itot_and
+ integer, intent(in) :: orb, Nint
- character*(2048)                :: output(1)
+ integer(bit_kind), intent(in) :: bitmask(Nint)
- string_bite = 0_bit_kind
+ integer :: j, k
- call set_bit_to_integer(bite,string_bite,Nint)
+ k = ishft(orb-1,-bit_kind_shift)+1
- itot = 0
+ j = orb-ishft(k-1,bit_kind_shift)-1
- itot_and = 0
+ is_integer_in_string = iand(bitmask(k), ibset(0_bit_kind, j)) /= 0_bit_kind
- is_integer_in_string = .False.
-!print*,''
-!print*,''
-!print*,'bite = ',bite
-!call bitstring_to_str( output(1), string_bite, Nint )
-! print *,  trim(output(1))
-!call bitstring_to_str( output(1), string, Nint )
-! print *,  trim(output(1))
- do i = 1, Nint
-  itot += popcnt(string(i))
-  itot_and += popcnt(ior(string(i),string_bite(i)))
- enddo
-!print*,'itot,itot_and',itot,itot_and
- if(itot == itot_and)then
-  is_integer_in_string = .True. 
- endif
-!pause 
 end

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

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
+  integer          :: i
-  double precision :: tmp_array(ao_num), r(3)
 
   !$OMP PARALLEL DO               &
   !$OMP DEFAULT (NONE)            &
-  !$OMP PRIVATE (i,r,tmp_array,j) & 
+  !$OMP PRIVATE (i) &
-  !$OMP SHARED(aos_in_r_array,n_points_final_grid,ao_num,final_grid_points)
+  !$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)
+    call give_all_aos_at_r(final_grid_points(1,i), aos_in_r_array(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
   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 :: r(3)
-  double precision :: aos_grad_array(3,ao_num)
+  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)
+    call give_all_aos_and_grad_at_r(final_grid_points(1,i),aos_array,aos_grad_array)
-    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)
     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, allocatable :: aos_lapl_array(:,:), aos_grad_array(:,:), aos_array(:)
- double precision :: aos_grad_array(3,ao_num)
+
- double precision :: aos_lapl_array(3,ao_num)
+ !$OMP PARALLEL        &
- !$OMP PARALLEL DO &
  !$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)
+  call give_all_aos_and_grad_and_lapl_at_r(final_grid_points(1,i),aos_array,aos_grad_array,aos_lapl_array)
-  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)
   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
+ integer :: i
- double precision :: aos_array(ao_num), r(3)
  !$OMP PARALLEL DO &
- !$OMP DEFAULT (NONE)  &
+ !$OMP DEFAULT (NONE)  PRIVATE (i) &
- !$OMP PRIVATE (i,r,aos_array,j) & 
+ !$OMP SHARED(aos_in_r_array_extra,n_points_extra_final_grid,final_grid_points_extra)
- !$OMP SHARED(aos_in_r_array_extra,n_points_extra_final_grid,ao_num,final_grid_points_extra)
  do i = 1, n_points_extra_final_grid
-  r(1) = final_grid_points_extra(1,i)
+  call give_all_aos_at_r(final_grid_points_extra(1,i),aos_in_r_array_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
  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, allocatable :: aos_array(:), aos_grad_array(:,:)
-  double precision :: aos_grad_array(3,ao_num)
 
-  !$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)
+    call give_all_aos_and_grad_at_r(final_grid_points_extra(1,i), aos_array, aos_grad_array)
-    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)
     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
 
 ! ---
 

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
+ integer :: i
- double precision :: mos_array(mo_num), r(3)
+ !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE (i) & 
- !$OMP PARALLEL DO &
- !$OMP DEFAULT (NONE)  &
- !$OMP PRIVATE (i,r,mos_array,j) & 
  !$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)
+  call give_all_mos_at_r(final_grid_points(1,i),mos_in_r_array_omp(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
  enddo
  !$OMP END PARALLEL DO
  END_PROVIDER

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( 2.d0*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
+

src/two_body_rdm/act_2_rdm.irp.f

@@ -145,6 +145,7 @@
  print*,''
  print*,'Providing act_2_rdm_spin_trace_mo '
  character*(128) :: name_file
+ PROVIDE all_mo_integrals
  name_file = 'act_2_rdm_spin_trace_mo'
  ispin = 4
  act_2_rdm_spin_trace_mo = 0.d0

src/two_rdm_routines/davidson_like_2rdm.irp.f

@@ -13,7 +13,7 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
    END_DOC
    integer, intent(in)             :: N_st,sze
    integer, intent(in)             :: dim1,norb,list_orb(norb),ispin
-   double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
+   double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st)
    double precision, intent(in)    :: u_0(sze,N_st)
 
    integer                        :: k
@@ -50,7 +50,7 @@ subroutine orb_range_2_rdm_openmp_work(big_array,dim1,norb,list_orb,ispin,u_t,N_
    END_DOC
    integer, intent(in)            :: N_st,sze,istart,iend,ishift,istep
    integer, intent(in)             :: dim1,norb,list_orb(norb),ispin
-   double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
+   double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st)
    double precision, intent(in)   :: u_t(N_st,N_det)
 
    integer :: k
@@ -91,7 +91,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
    integer, intent(in)            :: N_st,sze,istart,iend,ishift,istep
    double precision, intent(in)   :: u_t(N_st,N_det)
    integer, intent(in)            :: dim1,norb,list_orb(norb),ispin
-   double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
+   double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1,N_st)
 
    integer(omp_lock_kind)          :: lock_2rdm
    integer                        :: i,j,k,l
@@ -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
@@ -154,6 +155,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
    ! Prepare the array of all alpha single excitations
    ! -------------------------------------------------
 
+   double precision, allocatable :: big_array_local(:,:,:,:,:)
+
    PROVIDE N_int nthreads_davidson elec_alpha_num
     !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson)      &
             !$OMP   SHARED(psi_bilinear_matrix_rows, N_det,lock_2rdm,&
@@ -173,7 +176,7 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
             !$OMP          buffer, doubles, n_doubles,               &
             !$OMP          tmp_det2, idx, l, kcol_prev,              &
             !$OMP          singles_a, n_singles_a, singles_b,        &
-            !$OMP          n_singles_b, nkeys, keys, values)
+            !$OMP          n_singles_b, nkeys, keys, values, big_array_local)
 
    ! Alpha/Beta double excitations
    ! =============================
@@ -184,6 +187,8 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
        singles_b(maxab),                                             &
        doubles(maxab),                                               &
        idx(maxab))
+   allocate( big_array_local(N_states,dim1, dim1, dim1, dim1) )
+   big_array_local(:,:,:,:,:) = 0.d0
 
    kcol_prev=-1
 
@@ -191,8 +196,9 @@ 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
+!print *, 'aa', k_a, '/', iend
 
      krow = psi_bilinear_matrix_rows(k_a)
      ASSERT (krow <= N_det_alpha_unique)
@@ -254,33 +260,36 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
          do l= 1, N_states
            c_1(l) = u_t(l,l_a) * u_t(l,k_a)
          enddo
-         if(alpha_beta)then
+!         if(alpha_beta)then
-          ! only ONE contribution
+!          ! only ONE contribution
-          if (nkeys+1 .ge. sze_buff) then
+!          if (nkeys+1 .ge. sze_buff) then
-            call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+!            call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
-            nkeys = 0
+!            nkeys = 0
-          endif
+!          endif
-         else if (spin_trace)then
+!         else if (spin_trace)then
-          ! TWO contributions
+!          ! TWO contributions
           if (nkeys+2 .ge. sze_buff) then
-            call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+!            call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+            call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
             nkeys = 0
           endif
-         endif
+!         endif
          call orb_range_off_diag_double_to_all_states_ab_dm_buffer(tmp_det,tmp_det2,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
      enddo
 
    enddo
-     !$OMP END DO
+     !$OMP END DO NOWAIT
+!     call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+     call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
+     nkeys = 0
 
-     !$OMP DO SCHEDULE(dynamic,64)
+     !$OMP DO SCHEDULE(dynamic)
    do k_a=istart+ishift,iend,istep
+!print *, 'ab', k_a, '/', iend
 
 
      ! Single and double alpha exitations
@@ -331,36 +340,39 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
      ! ----------------------------------
 
      tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
-     do i=1,n_singles_a
+     if(alpha_beta.or.spin_trace.or.alpha_alpha)then
-       l_a = singles_a(i)
+       do i=1,n_singles_a
-       ASSERT (l_a <= N_det)
+         l_a = singles_a(i)
+         ASSERT (l_a <= N_det)
 
-       lrow = psi_bilinear_matrix_rows(l_a)
+         lrow = psi_bilinear_matrix_rows(l_a)
-       ASSERT (lrow <= N_det_alpha_unique)
+         ASSERT (lrow <= N_det_alpha_unique)
+
+         tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
+         do l= 1, N_states
+           c_1(l) = u_t(l,l_a) * u_t(l,k_a)
+         enddo
+
+         ! 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)
+           call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
+           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)
+           call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
+           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)
 
-       tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
-       do l= 1, N_states
-         c_1(l) = u_t(l,l_a) * u_t(l,k_a)
        enddo
-       if(alpha_beta.or.spin_trace.or.alpha_alpha)then
+     endif
-       ! 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
      ! ----------------------------------
@@ -377,14 +389,15 @@ 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 (nkeys+4 .ge. sze_buff) then
-          call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+!          call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+          call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
           nkeys = 0
         endif
         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)
+!     call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
-     nkeys = 0
+!     nkeys = 0
 
 
      ! Single and double beta excitations
@@ -432,35 +445,39 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
      ! ----------------------------------
 
      tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
-     do i=1,n_singles_b
+     if(alpha_beta.or.spin_trace.or.beta_beta)then
-       l_b = singles_b(i)
+      do i=1,n_singles_b
-       ASSERT (l_b <= N_det)
+        l_b = singles_b(i)
+        ASSERT (l_b <= N_det)
 
-       lcol = psi_bilinear_matrix_transp_columns(l_b)
+        lcol = psi_bilinear_matrix_transp_columns(l_b)
-       ASSERT (lcol <= N_det_beta_unique)
+        ASSERT (lcol <= N_det_beta_unique)
+
+        tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
+        l_a = psi_bilinear_matrix_transp_order(l_b)
+        do l= 1, N_states
+          c_1(l) = u_t(l,l_a) * u_t(l,k_a)
+        enddo
 
-       tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
-       l_a = psi_bilinear_matrix_transp_order(l_b)
-       do l= 1, N_states
-         c_1(l) = u_t(l,l_a) * u_t(l,k_a)
-       enddo
-       if(alpha_beta.or.spin_trace.or.beta_beta)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)
+!          call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+          call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
           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 beta /beta  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)
+!          call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+          call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
           nkeys = 0
         endif
         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
-     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
 
      ! Compute Hij for all beta doubles
      ! ----------------------------------
@@ -478,7 +495,8 @@ 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 (nkeys+4 .ge. sze_buff) then
-          call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+!          call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+          call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
           nkeys = 0
         endif
         call orb_range_off_diag_double_to_all_states_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
@@ -487,8 +505,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)
+!     call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
-     nkeys = 0
+!     nkeys = 0
 
 
      ! Diagonal contribution
@@ -514,16 +532,28 @@ subroutine orb_range_2_rdm_openmp_work_$N_int(big_array,dim1,norb,list_orb,ispin
        c_1(l) = u_t(l,k_a) * u_t(l,k_a)
      enddo
 
-     call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+     if (nkeys+elec_alpha_num*elec_alpha_num .ge. sze_buff) then
-     nkeys = 0
+!       call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+       call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
+       nkeys = 0
+     endif
+
      call orb_range_diag_to_all_states_2_rdm_dm_buffer(tmp_det,c_1,N_states,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-     call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+
+!     call update_keys_values_n_states(keys,values,nkeys,dim1,n_st,big_array,lock_2rdm)
+     call update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
      nkeys = 0
 
    end do
-    !$OMP END DO
+    !$OMP END DO NOWAIT
    deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values)
-    !$OMP END PARALLEL
+   !$OMP CRITICAL
+   do i=1,N_states
+     big_array(:,:,:,:,i) = big_array(:,:,:,:,i) + big_array_local(i,:,:,:,:)
+   enddo
+   !$OMP END CRITICAL
+   deallocate(big_array_local)
+   !$OMP END PARALLEL
 
 end
 
@@ -550,22 +580,66 @@ 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 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
+
+subroutine update_keys_values_n_states_local(keys,values,nkeys,dim1,n_st,big_array_local)
+ use omp_lib
+ implicit none
+ integer, intent(in) :: n_st,nkeys,dim1
+ integer, intent(in) :: keys(4,nkeys)
+ double precision, intent(in) :: values(n_st,nkeys)
+ double precision, intent(inout) :: big_array_local(n_st,dim1,dim1,dim1,dim1)
+
+ integer :: istate
+ integer :: i,h1,h2,p1,p2
+
  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_local(istate,h1,h2,p1,p2) = big_array_local(istate,h1,h2,p1,p2) + values(istate,i)
-    big_array(h1,h2,p1,p2,istate) += values(istate,i)
    enddo
  enddo
- call omp_unset_lock(lock_2rdm)
 
 end