added sparse cholesky mu_of_r

external/irpf90

@@ -1 +1 @@
-Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102
+Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6

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
 

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 '

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

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 PARALLEL DEFAULT(NONE)                                      &
-  !$OMP DEFAULT (NONE)  &
+  !$OMP PRIVATE (accu_vec,ipoint,p,ii,i,mm,m,mo_i_r1,mo_b_r1) & 
-  !$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,thresh_1,thresh_2) & 
-  !$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_cholesky_sparse,n_points_final_grid,cholesky_mo_transp,n_basis_orb) 
-  !$OMP ShARED (cholesky_mo_num,f_hf_sparse_cholesky,n_points_final_grid,cholesky_mo,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
+    f_hf_cholesky_sparse(ipoint) = 0.d0
-    do p = 1, cholesky_mo_num
+     accu_vec = 0.d0
-     accu_1 = 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
-     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
-  !$OMP END PARALLEL DO
+    f_hf_cholesky_sparse(ipoint) *= 2.D0
+   enddo
+  !$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 
 

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