Merge remote-tracking branch 'origin/cleaning_dft' into dev

src/ao_one_e_ints/ao_one_e_ints.irp.f

@@ -11,9 +11,6 @@
   ELSE
         ao_one_e_integrals = ao_integrals_n_e + ao_kinetic_integrals
 
-        IF (DO_PSEUDO) THEN
-              ao_one_e_integrals += ao_pseudo_integrals
-        ENDIF
   ENDIF
 
   DO j = 1, ao_num

src/ao_one_e_ints/pot_ao_ints.irp.f

@@ -76,6 +76,12 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
     !$OMP END DO
     !$OMP END PARALLEL
   endif
+
+  IF (DO_PSEUDO) THEN                                                                                               
+        ao_integrals_n_e += ao_pseudo_integrals
+  ENDIF
+
+
   if (write_ao_integrals_n_e) then
     call ezfio_set_ao_one_e_ints_ao_integrals_n_e(ao_integrals_n_e)
     print *,  'AO N-e integrals written to disk'

src/dft_utils_func/ecmd_pbe_general.irp.f

@@ -26,6 +26,14 @@ subroutine ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,
   pi = 4.d0 * datan(1.d0)
 
   eps_c_md_on_top_PBE = 0.d0
+  ! convertion from (alpha,beta) formalism to (closed, open) formalism for the density 
+  call rho_ab_to_rho_oc(rho_a,rho_b,rhoo,rhoc)
+  if(rhoc.lt.1.d-10)then
+   return
+  else if(on_top/(rhoc**2) .lt. 1.d-6)then
+   return
+  endif
+
   grad_rho_a_2 = 0.d0
   grad_rho_b_2 = 0.d0
   grad_rho_a_b = 0.d0
@@ -34,8 +42,7 @@ subroutine ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,
    grad_rho_b_2 += grad_rho_b(m)*grad_rho_b(m)
    grad_rho_a_b += grad_rho_a(m)*grad_rho_b(m)
   enddo
-  ! convertion from (alpha,beta) formalism to (closed, open) formalism
-  call rho_ab_to_rho_oc(rho_a,rho_b,rhoo,rhoc)
+  ! same same for gradients : convertion from (alpha,beta) formalism to (closed, open) formalism
   call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,sigmaoo,sigmacc,sigmaco)
 
   ! usual PBE correlation energy using the density, spin polarization and density gradients for alpha/beta electrons

src/mo_one_e_ints/mo_one_e_ints.irp.f

@@ -12,10 +12,6 @@ BEGIN_PROVIDER [ double precision, mo_one_e_integrals,(mo_num,mo_num)]
   ELSE
       mo_one_e_integrals  = mo_integrals_n_e + mo_kinetic_integrals
 
-      IF (DO_PSEUDO) THEN
-            mo_one_e_integrals  += mo_pseudo_integrals
-      ENDIF
-
   ENDIF
 
   IF (write_mo_one_e_integrals) THEN

src/mu_of_r/f_psi_i_a_v_utils.irp.f

@@ -311,3 +311,34 @@ BEGIN_PROVIDER [double precision, two_e_int_ii_f, (n_basis_orb,n_basis_orb,n_ina
  enddo
 END_PROVIDER 
 
+
+subroutine give_mu_of_r_cas(r,istate,mu_of_r,f_psi,n2_psi)
+ implicit none
+ BEGIN_DOC
+! returns mu(r), f_psi, n2_psi for a general cas wave function
+ END_DOC
+ integer, intent(in) :: istate
+ double precision, intent(in)  :: r(3)
+ double precision, intent(out) :: mu_of_r,f_psi,n2_psi
+ double precision :: f_ii_val_ab,two_bod_dens_ii
+ double precision :: f_ia_val_ab,two_bod_dens_ia
+ double precision :: f_aa_val_ab,two_bod_dens_aa
+ double precision :: sqpi,w_psi
+ sqpi = dsqrt(dacos(-1.d0))
+ ! inactive-inactive part of f_psi(r1,r2)
+ call give_f_ii_val_ab(r,r,f_ii_val_ab,two_bod_dens_ii)
+ ! inactive-active part of f_psi(r1,r2)
+ call give_f_ia_val_ab(r,r,f_ia_val_ab,two_bod_dens_ia,istate)
+ ! active-active part of f_psi(r1,r2)
+ call give_f_aa_val_ab(r,r,f_aa_val_ab,two_bod_dens_aa,istate)
+
+ f_psi  = f_ii_val_ab     + f_ia_val_ab     + f_aa_val_ab 
+ n2_psi = two_bod_dens_ii + two_bod_dens_ia + two_bod_dens_aa
+ if(n2_psi.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * n2_psi.lt.0.d0)then
+   w_psi   = 1.d+10
+ else 
+   w_psi  = f_psi /  n2_psi
+ endif
+ mu_of_r  = w_psi * sqpi * 0.5d0
+ 
+end

src/two_body_rdm/act_2_rdm.irp.f

@@ -25,6 +25,7 @@
  print*,'Providing act_2_rdm_ab_mo '
  ispin = 3 
  act_2_rdm_ab_mo = 0.d0
+ provide mo_two_e_integrals_in_map
  call wall_time(wall_1)
  if(read_two_body_rdm_ab)then
   print*,'Reading act_2_rdm_ab_mo from disk ...'

src/two_rdm_routines/davidson_like_2rdm.irp.f

@@ -19,6 +19,7 @@ subroutine orb_range_2_rdm_openmp(big_array,dim1,norb,list_orb,ispin,u_0,N_st,sz
    integer                        :: k
    double precision, allocatable  :: u_t(:,:)
    !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
+   PROVIDE mo_two_e_integrals_in_map
    allocate(u_t(N_st,N_det))
    do k=1,N_st
      call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)

src/utils/prim_in_r.irp.f

@@ -0,0 +1,41 @@
+double precision function primitive_value_explicit(power_prim,center_prim,alpha,r)
+ implicit none
+ BEGIN_DOC
+! Evaluates at "r" a primitive of type : 
+! (x - center_prim(1))**power_prim(1)  (y - center_prim(2))**power_prim(2) * (z - center_prim(3))**power_prim(3)
+!
+! exp(-alpha * [(x - center_prim(1))**2 + (y - center_prim(2))**2 + (z - center_prim(3))**2] )
+ END_DOC
+ integer, intent(in)          :: power_prim(3)
+ double precision, intent(in) :: center_prim(3),alpha
+ double precision, intent(in) :: r(3)
+
+ double precision :: dx,dy,dz,r2
+ dx = (r(1) - center_prim(1))
+ dy = (r(2) - center_prim(2))
+ dz = (r(3) - center_prim(3))
+ r2 = dx*dx + dy*dy + dz*dz
+ dx = dx**power_prim(1)
+ dy = dy**power_prim(2)
+ dz = dz**power_prim(3)
+
+ primitive_value_explicit = dexp(-alpha*r2) * dx * dy * dz
+
+end
+
+double precision function give_pol_in_r(r,pol,center, alpha,iorder, max_dim)
+ double precision    :: r(3), center(3), alpha,pol(0:max_dim,3)
+ integer, intent(in) :: iorder(3), max_dim
+ integer :: i,m
+ double precision :: gauss(3), x
+ gauss  = 0.d0
+
+ do m = 1, 3
+  x = r(m) - center(m)
+  do i = 0, iorder(m)
+   gauss(m) += pol(i,m) * dexp(-alpha *x**2 ) * x**i 
+  enddo
+ enddo
+ give_pol_in_r = gauss(1) * gauss(2) * gauss(3)
+
+end