cleaning in dft

src/README.rst

@@ -1,3 +1,51 @@
 ==========================
 The core modules of the QP
 ==========================
+
+*** How are handled the DFT functionals in QP2 ?
+================================================
+    The Exchange and Correlation energies/potentials can be accessed by the following providers 
+    energy_x
+    energy_c
+    potential_x_alpha_ao
+    potential_c_alpha_ao
+    potential_x_beta_ao
+    potential_c_beta_ao
+
+    These providers are automatically linked to the providers of the actual exchange/correlation energies of a given functional 
+    through the character keywords 
+    "exchange_functional"
+    "correlation_functional"
+
+    All the providers for the available functionals are in the folder "functionals", with one file "my_functional.irp.f" per functional.  
+
+    Ex : if "exchange_functional" == "sr_pbe", then energy_x will contain the exchange correlation functional defined in "functiona/sr_pbe.irp.f", which corresponds to the short-range PBE functional (at the value mu_erf for the range separation parameter) 
+
+
+*** How are handled the DFT functionals in QP2 ?
+================================================
+
+    Creating a new functional and propagating it through the whole QP2 programs is easy as all dependencies are handled by a script. 
+
+    To do so, let us assume that the name of your functional is "my_func". 
+    Then you just have to create the file "my_func.irp.f" in the folder "functional" which shoud contain 
+
+    +) if you're adding an exchange functional, then create the provider "energy_x_my_func" 
+
+    +) if you're adding a correlation functional, create the provider "energy_c_my_func" 
+    
+    +) if you want to add the echange potentials, create the providers "potential_x_alpha_ao_my_func", "potential_x_beta_ao_my_func"  which are the exchange potentials on the AO basis for the alpha/beta electrons 
+
+    +) if you want to add the correlation potentials, create the providers "potential_c_alpha_ao_my_func", "potential_c_beta_ao_my_func"  which are the correlation potentials on the AO basis for the alpha/beta electrons 
+    
+    That's all :) 
+   
+    Then, when running whatever DFT calculation or accessing/using the providers: 
+    energy_x
+    energy_c
+    potential_x_alpha_ao
+    potential_c_alpha_ao
+    potential_x_beta_ao
+    potential_c_beta_ao
+ 
+    if exchange_functional = mu_func, then you will automatically have access to what you need, such as kohn sham orbital optimization and so on ... 

src/dft_utils_one_e/utils.irp.f

@@ -1,58 +1,32 @@
 
-subroutine GGA_sr_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, &
+subroutine GGA_sr_type_functionals(mu,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, &
                                 ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &
                                 ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
  implicit none
  BEGIN_DOC
  ! routine that helps in building the x/c potentials on the AO basis for a GGA functional with a short-range interaction
  END_DOC
- double precision, intent(in)  :: r(3),rho_a(N_states),rho_b(N_states),grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
+ double precision, intent(in)  :: mu,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
- double precision, intent(out) :: ex(N_states),vx_rho_a(N_states),vx_rho_b(N_states),vx_grad_rho_a_2(N_states),vx_grad_rho_b_2(N_states),vx_grad_rho_a_b(N_states)
+ double precision, intent(out) :: ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b
- double precision, intent(out) :: ec(N_states),vc_rho_a(N_states),vc_rho_b(N_states),vc_grad_rho_a_2(N_states),vc_grad_rho_b_2(N_states),vc_grad_rho_a_b(N_states)
+ double precision, intent(out) :: ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b
- integer          :: istate
- double precision :: r2(3),dr2(3), local_potential,r12,dx2,mu
- do istate = 1, N_states
-  call ex_pbe_sr(mu_erf_dft,rho_a(istate),rho_b(istate),grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),ex(istate),vx_rho_a(istate),vx_rho_b(istate),vx_grad_rho_a_2(istate),vx_grad_rho_b_2(istate),vx_grad_rho_a_b(istate))
-
-  double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
-  ! convertion from (alpha,beta) formalism to (closed, open) formalism
-   call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
-   call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
-
-   call ec_pbe_sr(mu_erf_dft,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec(istate),vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo)
-
-   call v_rho_oc_to_v_rho_ab(vrhoo,vrhoc,vc_rho_a(istate),vc_rho_b(istate))
-   call v_grad_rho_oc_to_v_grad_rho_ab(vsigmaoo,vsigmacc,vsigmaco,vc_grad_rho_a_2(istate),vc_grad_rho_b_2(istate),vc_grad_rho_a_b(istate))
- enddo
-end
-
-
-subroutine GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, &
-                                ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &
-                                ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b )
- implicit none
- BEGIN_DOC
- ! routine that helps in building the x/c potentials on the AO basis for a GGA functional
- END_DOC
- double precision, intent(in)  :: r(3),rho_a(N_states),rho_b(N_states),grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
- double precision, intent(out) :: ex(N_states),vx_rho_a(N_states),vx_rho_b(N_states),vx_grad_rho_a_2(N_states),vx_grad_rho_b_2(N_states),vx_grad_rho_a_b(N_states)
- double precision, intent(out) :: ec(N_states),vc_rho_a(N_states),vc_rho_b(N_states),vc_grad_rho_a_2(N_states),vc_grad_rho_b_2(N_states),vc_grad_rho_a_b(N_states)
  integer          :: istate
  double precision :: r2(3),dr2(3), local_potential,r12,dx2
- double precision :: mu_local
- mu_local = 1.d-9
- do istate = 1, N_states
-   call ex_pbe_sr(mu_local,rho_a(istate),rho_b(istate),grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),ex(istate),vx_rho_a(istate),vx_rho_b(istate),vx_grad_rho_a_2(istate),vx_grad_rho_b_2(istate),vx_grad_rho_a_b(istate))
 
   double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo
+
+   ! exhange energy and potentials 
+   call ex_pbe_sr(mu,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b)
+
    ! convertion from (alpha,beta) formalism to (closed, open) formalism
-   call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
+   call rho_ab_to_rho_oc(rho_a,rho_b,rhoo,rhoc)
-   call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
+   call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,sigmaoo,sigmacc,sigmaco)
 
-   call ec_pbe_sr(mu_local,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec(istate),vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo)
+   ! correlation energy and potentials 
+   call ec_pbe_sr(mu,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,ec,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo)
 
-   call v_rho_oc_to_v_rho_ab(vrhoo,vrhoc,vc_rho_a(istate),vc_rho_b(istate))
+   ! convertion from (closed, open) formalism to (alpha,beta) formalism 
-   call v_grad_rho_oc_to_v_grad_rho_ab(vsigmaoo,vsigmacc,vsigmaco,vc_grad_rho_a_2(istate),vc_grad_rho_b_2(istate),vc_grad_rho_a_b(istate))
+   call v_rho_oc_to_v_rho_ab(vrhoo,vrhoc,vc_rho_a,vc_rho_b)
- enddo
+   call v_grad_rho_oc_to_v_grad_rho_ab(vsigmaoo,vsigmacc,vsigmaco,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b)
 end
 
+

src/functionals/pbe.irp.f

@@ -6,7 +6,6 @@
 ! exchange/correlation energy with the short range pbe functional
  END_DOC
  integer :: istate,i,j,m
- double precision :: r(3)
  double precision :: mu,weight
  double precision, allocatable :: ex(:), ec(:)
  double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
@@ -22,9 +21,6 @@
  energy_x_pbe = 0.d0
  do istate = 1, N_states
   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)
    weight = final_weight_at_r_vector(i)
    rho_a(istate) =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
    rho_b(istate) =  one_e_dm_and_grad_beta_in_r(4,i,istate)
@@ -40,7 +36,7 @@
    enddo
 
                              ! inputs
-   call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
+   call GGA_sr_type_functionals(0.d0,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
                              ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &  ! outputs correlation
                              ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b  )
    energy_x_pbe += ex * weight
@@ -56,7 +52,6 @@ BEGIN_PROVIDER[double precision, energy_c_pbe, (N_states) ]
 ! exchange/correlation energy with the short range pbe functional
  END_DOC
  integer :: istate,i,j,m
- double precision :: r(3)
  double precision :: mu,weight
  double precision, allocatable :: ex(:), ec(:)
  double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
@@ -72,9 +67,6 @@ BEGIN_PROVIDER[double precision, energy_c_pbe, (N_states) ]
  energy_c_pbe = 0.d0
  do istate = 1, N_states
   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)
    weight = final_weight_at_r_vector(i)
    rho_a(istate) =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
    rho_b(istate) =  one_e_dm_and_grad_beta_in_r(4,i,istate)
@@ -90,7 +82,7 @@ BEGIN_PROVIDER[double precision, energy_c_pbe, (N_states) ]
    enddo
 
                              ! inputs
-   call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
+   call GGA_sr_type_functionals(0.d0,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
                              ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &  ! outputs correlation
                              ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b  )
    energy_c_pbe += ec * weight
@@ -159,7 +151,6 @@ END_PROVIDER
 ! aos_vxc_alpha_pbe_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
  END_DOC
  integer :: istate,i,j,m
- double precision :: r(3)
  double precision :: mu,weight
  double precision, allocatable :: ex(:), ec(:)
  double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
@@ -179,9 +170,6 @@ END_PROVIDER
 
  do istate = 1, N_states
   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)
    weight = final_weight_at_r_vector(i)
    rho_a(istate) =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
    rho_b(istate) =  one_e_dm_and_grad_beta_in_r(4,i,istate)
@@ -197,7 +185,7 @@ END_PROVIDER
    enddo
 
                              ! inputs
-   call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
+   call GGA_sr_type_functionals(0.d0,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
                              ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &  ! outputs correlation
                              ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b  )
    vx_rho_a(istate) *= weight
@@ -325,63 +313,54 @@ END_PROVIDER
 ! aos_vxc_alpha_pbe_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
  END_DOC
  integer :: istate,i,j,m
- double precision :: r(3)
  double precision :: mu,weight
- double precision, allocatable :: ex(:), ec(:)
+ double precision :: ex, ec
- double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
+ double precision :: rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
- double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
+ double precision :: contrib_grad_xa(3),contrib_grad_xb(3),contrib_grad_ca(3),contrib_grad_cb(3)
- double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
+ double precision :: vc_rho_a, vc_rho_b, vx_rho_a, vx_rho_b
- double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
+ double precision :: vx_grad_rho_a_2, vx_grad_rho_b_2, vx_grad_rho_a_b, vc_grad_rho_a_2, vc_grad_rho_b_2, vc_grad_rho_a_b
- allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
- allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
- allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
- allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
- allocate(contrib_grad_xa(3,N_states),contrib_grad_xb(3,N_states),contrib_grad_ca(3,N_states),contrib_grad_cb(3,N_states))
 
   aos_dvxc_alpha_pbe_w = 0.d0
   aos_dvxc_beta_pbe_w  = 0.d0
 
  do istate = 1, N_states
   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)
    weight = final_weight_at_r_vector(i)
-   rho_a(istate) =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
+   rho_a =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
-   rho_b(istate) =  one_e_dm_and_grad_beta_in_r(4,i,istate)
+   rho_b =  one_e_dm_and_grad_beta_in_r(4,i,istate)
-   grad_rho_a(1:3,istate) =  one_e_dm_and_grad_alpha_in_r(1:3,i,istate)
+   grad_rho_a(1:3) =  one_e_dm_and_grad_alpha_in_r(1:3,i,istate)
-   grad_rho_b(1:3,istate) =  one_e_dm_and_grad_beta_in_r(1:3,i,istate)
+   grad_rho_b(1:3) =  one_e_dm_and_grad_beta_in_r(1:3,i,istate)
    grad_rho_a_2 = 0.d0
    grad_rho_b_2 = 0.d0
    grad_rho_a_b = 0.d0
    do m = 1, 3
-    grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
+    grad_rho_a_2 += grad_rho_a(m) * grad_rho_a(m)
-    grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
+    grad_rho_b_2 += grad_rho_b(m) * grad_rho_b(m)
-    grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
+    grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m)
    enddo
 
-                             ! inputs
+   ! call exc_sr_pbe 
-   call GGA_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
+   call GGA_sr_type_functionals(0.d0,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,         &  ! inputs
-                             ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &  ! outputs correlation
+                             ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &  ! outputs exchange
-                             ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b  )
+                             ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b  )  ! outputs correlation
-   vx_rho_a(istate) *= weight
+   vx_rho_a *= weight
-   vc_rho_a(istate) *= weight
+   vc_rho_a *= weight
-   vx_rho_b(istate) *= weight
+   vx_rho_b *= weight
-   vc_rho_b(istate) *= weight
+   vc_rho_b *= weight
    do m= 1,3
-    contrib_grad_ca(m,istate) = weight * (2.d0 * vc_grad_rho_a_2(istate) *  grad_rho_a(m,istate) + vc_grad_rho_a_b(istate)  * grad_rho_b(m,istate))
+    contrib_grad_ca(m) = weight * (2.d0 * vc_grad_rho_a_2 *  grad_rho_a(m) + vc_grad_rho_a_b  * grad_rho_b(m))
-    contrib_grad_xa(m,istate) = weight * (2.d0 * vx_grad_rho_a_2(istate) *  grad_rho_a(m,istate) + vx_grad_rho_a_b(istate)  * grad_rho_b(m,istate))
+    contrib_grad_xa(m) = weight * (2.d0 * vx_grad_rho_a_2 *  grad_rho_a(m) + vx_grad_rho_a_b  * grad_rho_b(m))
-    contrib_grad_cb(m,istate) = weight * (2.d0 * vc_grad_rho_b_2(istate) *  grad_rho_b(m,istate) + vc_grad_rho_a_b(istate)  * grad_rho_a(m,istate))
+    contrib_grad_cb(m) = weight * (2.d0 * vc_grad_rho_b_2 *  grad_rho_b(m) + vc_grad_rho_a_b  * grad_rho_a(m))
-    contrib_grad_xb(m,istate) = weight * (2.d0 * vx_grad_rho_b_2(istate) *  grad_rho_b(m,istate) + vx_grad_rho_a_b(istate)  * grad_rho_a(m,istate))
+    contrib_grad_xb(m) = weight * (2.d0 * vx_grad_rho_b_2 *  grad_rho_b(m) + vx_grad_rho_a_b  * grad_rho_a(m))
    enddo
    do j = 1, ao_num
-    aos_vxc_alpha_pbe_w(j,i,istate) = ( vc_rho_a(istate) + vx_rho_a(istate) ) * aos_in_r_array(j,i)
+    aos_vxc_alpha_pbe_w(j,i,istate) = ( vc_rho_a + vx_rho_a ) * aos_in_r_array(j,i)
-    aos_vxc_beta_pbe_w (j,i,istate) = ( vc_rho_b(istate) + vx_rho_b(istate) ) * aos_in_r_array(j,i)
+    aos_vxc_beta_pbe_w (j,i,istate) = ( vc_rho_b + vx_rho_b ) * aos_in_r_array(j,i)
    enddo
    do j = 1, ao_num
     do m = 1,3
-     aos_dvxc_alpha_pbe_w(j,i,istate) += ( contrib_grad_ca(m,istate) + contrib_grad_xa(m,istate) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
+     aos_dvxc_alpha_pbe_w(j,i,istate) += ( contrib_grad_ca(m) + contrib_grad_xa(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
-     aos_dvxc_beta_pbe_w (j,i,istate) += ( contrib_grad_cb(m,istate) + contrib_grad_xb(m,istate) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
+     aos_dvxc_beta_pbe_w (j,i,istate) += ( contrib_grad_cb(m) + contrib_grad_xb(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
     enddo
    enddo
   enddo

src/functionals/sr_pbe.irp.f

@@ -3,55 +3,95 @@
 &BEGIN_PROVIDER[double precision, energy_c_sr_pbe, (N_states) ]
  implicit none
  BEGIN_DOC
+ ! exchange / correlation energies  with the short-range version Perdew-Burke-Ernzerhof GGA functional 
+ !
+ ! defined in Chem. Phys.329, 276 (2006)
+ END_DOC 
+ BEGIN_DOC
 ! exchange/correlation energy with the short range pbe functional
  END_DOC
  integer :: istate,i,j,m
- double precision :: r(3)
  double precision :: mu,weight
- double precision, allocatable :: ex(:), ec(:)
+ double precision, allocatable :: ex, ec
- double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
+ double precision :: rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
- double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
+ double precision :: vc_rho_a, vc_rho_b, vx_rho_a, vx_rho_b
- double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
+ double precision :: vx_grad_rho_a_2, vx_grad_rho_b_2, vx_grad_rho_a_b, vc_grad_rho_a_2, vc_grad_rho_b_2, vc_grad_rho_a_b
- double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
- allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
- allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
 
 
- allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
- allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
  energy_x_sr_pbe = 0.d0
  energy_c_sr_pbe = 0.d0
  do istate = 1, N_states
   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)
    weight = final_weight_at_r_vector(i)
-   rho_a(istate) =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
+   rho_a =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
-   rho_b(istate) =  one_e_dm_and_grad_beta_in_r(4,i,istate)
+   rho_b =  one_e_dm_and_grad_beta_in_r(4,i,istate)
-   grad_rho_a(1:3,istate) =  one_e_dm_and_grad_alpha_in_r(1:3,i,istate)
+   grad_rho_a(1:3) =  one_e_dm_and_grad_alpha_in_r(1:3,i,istate)
-   grad_rho_b(1:3,istate) =  one_e_dm_and_grad_beta_in_r(1:3,i,istate)
+   grad_rho_b(1:3) =  one_e_dm_and_grad_beta_in_r(1:3,i,istate)
    grad_rho_a_2 = 0.d0
    grad_rho_b_2 = 0.d0
    grad_rho_a_b = 0.d0
    do m = 1, 3
-    grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
+    grad_rho_a_2 += grad_rho_a(m) * grad_rho_a(m)
-    grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
+    grad_rho_b_2 += grad_rho_b(m) * grad_rho_b(m)
-    grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
+    grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m)
    enddo
 
                              ! inputs
-   call GGA_sr_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
+   call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
                              ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &  ! outputs correlation
                              ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b  )
-   energy_x_sr_pbe += ex * weight
+   energy_x_sr_pbe(istate) += ex * weight
-   energy_c_sr_pbe += ec * weight
+   energy_c_sr_pbe(istate) += ec * weight
   enddo
  enddo
 
 
 END_PROVIDER
 
+ BEGIN_PROVIDER [double precision, potential_x_alpha_ao_sr_pbe,(ao_num,ao_num,N_states)]
+&BEGIN_PROVIDER [double precision, potential_x_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
+&BEGIN_PROVIDER [double precision, potential_c_alpha_ao_sr_pbe,(ao_num,ao_num,N_states)]
+&BEGIN_PROVIDER [double precision, potential_c_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
+   implicit none
+ BEGIN_DOC
+ ! exchange / correlation potential for alpha / beta electrons  with the short-range version Perdew-Burke-Ernzerhof GGA functional 
+ !
+ ! defined in Chem. Phys.329, 276 (2006)
+ END_DOC 
+   integer :: i,j,istate
+   do istate = 1, n_states 
+    do i = 1, ao_num
+     do j = 1, ao_num
+      potential_x_alpha_ao_sr_pbe(j,i,istate) = pot_sr_scal_x_alpha_ao_pbe(j,i,istate) + pot_sr_grad_x_alpha_ao_pbe(j,i,istate) + pot_sr_grad_x_alpha_ao_pbe(i,j,istate)
+      potential_x_beta_ao_sr_pbe(j,i,istate) = pot_sr_scal_x_beta_ao_pbe(j,i,istate) + pot_sr_grad_x_beta_ao_pbe(j,i,istate) + pot_sr_grad_x_beta_ao_pbe(i,j,istate)
+
+      potential_c_alpha_ao_sr_pbe(j,i,istate) = pot_sr_scal_c_alpha_ao_pbe(j,i,istate) + pot_sr_grad_c_alpha_ao_pbe(j,i,istate) + pot_sr_grad_c_alpha_ao_pbe(i,j,istate)
+      potential_c_beta_ao_sr_pbe(j,i,istate) = pot_sr_scal_c_beta_ao_pbe(j,i,istate) + pot_sr_grad_c_beta_ao_pbe(j,i,istate) + pot_sr_grad_c_beta_ao_pbe(i,j,istate)
+     enddo
+    enddo
+   enddo
+
+END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, potential_xc_alpha_ao_sr_pbe,(ao_num,ao_num,N_states)]
+&BEGIN_PROVIDER [double precision, potential_xc_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
+   implicit none
+ BEGIN_DOC
+ ! exchange / correlation potential for alpha / beta electrons  with the Perdew-Burke-Ernzerhof GGA functional 
+ END_DOC 
+   integer :: i,j,istate
+   do istate = 1, n_states 
+    do i = 1, ao_num
+     do j = 1, ao_num
+      potential_xc_alpha_ao_sr_pbe(j,i,istate) = pot_sr_scal_xc_alpha_ao_pbe(j,i,istate) + pot_sr_grad_xc_alpha_ao_pbe(j,i,istate) + pot_sr_grad_xc_alpha_ao_pbe(i,j,istate)
+      potential_xc_beta_ao_sr_pbe(j,i,istate)  = pot_sr_scal_xc_beta_ao_pbe(j,i,istate)  + pot_sr_grad_xc_beta_ao_pbe(j,i,istate)  + pot_sr_grad_xc_beta_ao_pbe(i,j,istate)
+     enddo
+    enddo
+   enddo
+
+END_PROVIDER 
+
+
 
  BEGIN_PROVIDER[double precision, aos_sr_vc_alpha_pbe_w  , (ao_num,n_points_final_grid,N_states)]
 &BEGIN_PROVIDER[double precision, aos_sr_vc_beta_pbe_w   , (ao_num,n_points_final_grid,N_states)]
@@ -63,72 +103,64 @@ END_PROVIDER
 &BEGIN_PROVIDER[double precision, aos_dsr_vx_beta_pbe_w   ,  (ao_num,n_points_final_grid,N_states)]
  implicit none
  BEGIN_DOC
+! intermediates to compute the sr_pbe potentials 
+! 
 ! aos_sr_vxc_alpha_pbe_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
  END_DOC
  integer :: istate,i,j,m
- double precision :: r(3)
  double precision :: mu,weight
- double precision, allocatable :: ex(:), ec(:)
+ double precision :: ex, ec
- double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
+ double precision :: rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
- double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
+ double precision :: contrib_grad_xa(3),contrib_grad_xb(3),contrib_grad_ca(3),contrib_grad_cb(3)
- double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
+ double precision :: vc_rho_a, vc_rho_b, vx_rho_a, vx_rho_b
- double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
+ double precision :: vx_grad_rho_a_2, vx_grad_rho_b_2, vx_grad_rho_a_b, vc_grad_rho_a_2, vc_grad_rho_b_2, vc_grad_rho_a_b
- allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
- allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
-
-
- allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
- allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
- allocate(contrib_grad_xa(3,N_states),contrib_grad_xb(3,N_states),contrib_grad_ca(3,N_states),contrib_grad_cb(3,N_states))
  aos_dsr_vc_alpha_pbe_w= 0.d0
  aos_dsr_vc_beta_pbe_w = 0.d0
  aos_dsr_vx_alpha_pbe_w= 0.d0
  aos_dsr_vx_beta_pbe_w = 0.d0
  do istate = 1, N_states
   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)
    weight = final_weight_at_r_vector(i)
-   rho_a(istate) =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
+
-   rho_b(istate) =  one_e_dm_and_grad_beta_in_r(4,i,istate)
+   rho_a =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
-   grad_rho_a(1:3,istate) =  one_e_dm_and_grad_alpha_in_r(1:3,i,istate)
+   rho_b =  one_e_dm_and_grad_beta_in_r(4,i,istate)
-   grad_rho_b(1:3,istate) =  one_e_dm_and_grad_beta_in_r(1:3,i,istate)
+   grad_rho_a(1:3) =  one_e_dm_and_grad_alpha_in_r(1:3,i,istate)
+   grad_rho_b(1:3) =  one_e_dm_and_grad_beta_in_r(1:3,i,istate)
    grad_rho_a_2 = 0.d0
    grad_rho_b_2 = 0.d0
    grad_rho_a_b = 0.d0
    do m = 1, 3
-    grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
+    grad_rho_a_2 += grad_rho_a(m) * grad_rho_a(m)
-    grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
+    grad_rho_b_2 += grad_rho_b(m) * grad_rho_b(m)
-    grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
+    grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m)
    enddo
 
                              ! inputs
-   call GGA_sr_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
+   call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
                              ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &  ! outputs correlation
                              ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b  )
-   vx_rho_a(istate) *= weight
+   vx_rho_a *= weight
-   vc_rho_a(istate) *= weight
+   vc_rho_a *= weight
-   vx_rho_b(istate) *= weight
+   vx_rho_b *= weight
-   vc_rho_b(istate) *= weight
+   vc_rho_b *= weight
    do m= 1,3
-    contrib_grad_ca(m,istate) = weight * (2.d0 * vc_grad_rho_a_2(istate) *  grad_rho_a(m,istate) + vc_grad_rho_a_b(istate)  * grad_rho_b(m,istate))
+    contrib_grad_ca(m) = weight * (2.d0 * vc_grad_rho_a_2 *  grad_rho_a(m) + vc_grad_rho_a_b  * grad_rho_b(m) )
-    contrib_grad_xa(m,istate) = weight * (2.d0 * vx_grad_rho_a_2(istate) *  grad_rho_a(m,istate) + vx_grad_rho_a_b(istate)  * grad_rho_b(m,istate))
+    contrib_grad_xa(m) = weight * (2.d0 * vx_grad_rho_a_2 *  grad_rho_a(m) + vx_grad_rho_a_b  * grad_rho_b(m) )
-    contrib_grad_cb(m,istate) = weight * (2.d0 * vc_grad_rho_b_2(istate) *  grad_rho_b(m,istate) + vc_grad_rho_a_b(istate)  * grad_rho_a(m,istate))
+    contrib_grad_cb(m) = weight * (2.d0 * vc_grad_rho_b_2 *  grad_rho_b(m) + vc_grad_rho_a_b  * grad_rho_a(m) )
-    contrib_grad_xb(m,istate) = weight * (2.d0 * vx_grad_rho_b_2(istate) *  grad_rho_b(m,istate) + vx_grad_rho_a_b(istate)  * grad_rho_a(m,istate))
+    contrib_grad_xb(m) = weight * (2.d0 * vx_grad_rho_b_2 *  grad_rho_b(m) + vx_grad_rho_a_b  * grad_rho_a(m) )
    enddo
    do j = 1, ao_num
-    aos_sr_vc_alpha_pbe_w(j,i,istate) = vc_rho_a(istate) * aos_in_r_array(j,i)
+    aos_sr_vc_alpha_pbe_w(j,i,istate) = vc_rho_a * aos_in_r_array(j,i)
-    aos_sr_vc_beta_pbe_w (j,i,istate) = vc_rho_b(istate) * aos_in_r_array(j,i)
+    aos_sr_vc_beta_pbe_w (j,i,istate) = vc_rho_b * aos_in_r_array(j,i)
-    aos_sr_vx_alpha_pbe_w(j,i,istate) = vx_rho_a(istate) * aos_in_r_array(j,i)
+    aos_sr_vx_alpha_pbe_w(j,i,istate) = vx_rho_a * aos_in_r_array(j,i)
-    aos_sr_vx_beta_pbe_w (j,i,istate) = vx_rho_b(istate) * aos_in_r_array(j,i)
+    aos_sr_vx_beta_pbe_w (j,i,istate) = vx_rho_b * aos_in_r_array(j,i)
    enddo
    do j = 1, ao_num
     do m = 1,3
-     aos_dsr_vc_alpha_pbe_w(j,i,istate) += contrib_grad_ca(m,istate) * aos_grad_in_r_array_transp_xyz(m,j,i)
+     aos_dsr_vc_alpha_pbe_w(j,i,istate) += contrib_grad_ca(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
-     aos_dsr_vc_beta_pbe_w (j,i,istate) += contrib_grad_cb(m,istate) * aos_grad_in_r_array_transp_xyz(m,j,i)
+     aos_dsr_vc_beta_pbe_w (j,i,istate) += contrib_grad_cb(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
-     aos_dsr_vx_alpha_pbe_w(j,i,istate) += contrib_grad_xa(m,istate) * aos_grad_in_r_array_transp_xyz(m,j,i)
+     aos_dsr_vx_alpha_pbe_w(j,i,istate) += contrib_grad_xa(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
-     aos_dsr_vx_beta_pbe_w (j,i,istate) += contrib_grad_xb(m,istate) * aos_grad_in_r_array_transp_xyz(m,j,i)
+     aos_dsr_vx_beta_pbe_w (j,i,istate) += contrib_grad_xb(m) * aos_grad_in_r_array_transp_xyz(m,j,i)
     enddo
    enddo
   enddo
@@ -142,6 +174,8 @@ END_PROVIDER
 &BEGIN_PROVIDER [double precision, pot_sr_scal_x_beta_ao_pbe, (ao_num,ao_num,N_states)]
 &BEGIN_PROVIDER [double precision, pot_sr_scal_c_beta_ao_pbe, (ao_num,ao_num,N_states)]
  implicit none
+! intermediates to compute the sr_pbe potentials 
+! 
  integer                        :: istate
    BEGIN_DOC
    ! intermediate quantity for the calculation of the vxc potentials for the GGA functionals  related to the scalar part of the potential 
@@ -220,29 +254,6 @@ END_PROVIDER
    
  call wall_time(wall_2)
 
-END_PROVIDER
-
- BEGIN_PROVIDER [double precision, potential_x_alpha_ao_sr_pbe,(ao_num,ao_num,N_states)]
-&BEGIN_PROVIDER [double precision, potential_x_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
-&BEGIN_PROVIDER [double precision, potential_c_alpha_ao_sr_pbe,(ao_num,ao_num,N_states)]
-&BEGIN_PROVIDER [double precision, potential_c_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
-   implicit none
- BEGIN_DOC
- ! exchange / correlation potential for alpha / beta electrons  with the Perdew-Burke-Ernzerhof GGA functional 
- END_DOC 
-   integer :: i,j,istate
-   do istate = 1, n_states 
-    do i = 1, ao_num
-     do j = 1, ao_num
-      potential_x_alpha_ao_sr_pbe(j,i,istate) = pot_sr_scal_x_alpha_ao_pbe(j,i,istate) + pot_sr_grad_x_alpha_ao_pbe(j,i,istate) + pot_sr_grad_x_alpha_ao_pbe(i,j,istate)
-      potential_x_beta_ao_sr_pbe(j,i,istate) = pot_sr_scal_x_beta_ao_pbe(j,i,istate) + pot_sr_grad_x_beta_ao_pbe(j,i,istate) + pot_sr_grad_x_beta_ao_pbe(i,j,istate)
-
-      potential_c_alpha_ao_sr_pbe(j,i,istate) = pot_sr_scal_c_alpha_ao_pbe(j,i,istate) + pot_sr_grad_c_alpha_ao_pbe(j,i,istate) + pot_sr_grad_c_alpha_ao_pbe(i,j,istate)
-      potential_c_beta_ao_sr_pbe(j,i,istate) = pot_sr_scal_c_beta_ao_pbe(j,i,istate) + pot_sr_grad_c_beta_ao_pbe(j,i,istate) + pot_sr_grad_c_beta_ao_pbe(i,j,istate)
-     enddo
-    enddo
-   enddo
-
 END_PROVIDER
 
 
@@ -255,65 +266,54 @@ END_PROVIDER
 ! aos_sr_vxc_alpha_pbe_w(j,i) = ao_i(r_j) * (v^x_alpha(r_j) + v^c_alpha(r_j)) * W(r_j)
  END_DOC
  integer :: istate,i,j,m
- double precision :: r(3)
  double precision :: mu,weight
- double precision, allocatable :: ex(:), ec(:)
+ double precision :: ex, ec
- double precision, allocatable :: rho_a(:),rho_b(:),grad_rho_a(:,:),grad_rho_b(:,:),grad_rho_a_2(:),grad_rho_b_2(:),grad_rho_a_b(:)
+ double precision :: rho_a,rho_b,grad_rho_a(3),grad_rho_b(3),grad_rho_a_2,grad_rho_b_2,grad_rho_a_b
- double precision, allocatable :: contrib_grad_xa(:,:),contrib_grad_xb(:,:),contrib_grad_ca(:,:),contrib_grad_cb(:,:)
+ double precision :: contrib_grad_xa(3),contrib_grad_xb(3),contrib_grad_ca(3),contrib_grad_cb(3)
- double precision, allocatable :: vc_rho_a(:), vc_rho_b(:), vx_rho_a(:), vx_rho_b(:)
+ double precision :: vc_rho_a, vc_rho_b, vx_rho_a, vx_rho_b
- double precision, allocatable :: vx_grad_rho_a_2(:), vx_grad_rho_b_2(:), vx_grad_rho_a_b(:), vc_grad_rho_a_2(:), vc_grad_rho_b_2(:), vc_grad_rho_a_b(:)
+ double precision :: vx_grad_rho_a_2, vx_grad_rho_b_2, vx_grad_rho_a_b, vc_grad_rho_a_2, vc_grad_rho_b_2, vc_grad_rho_a_b
- allocate(vc_rho_a(N_states), vc_rho_b(N_states), vx_rho_a(N_states), vx_rho_b(N_states))
- allocate(vx_grad_rho_a_2(N_states), vx_grad_rho_b_2(N_states), vx_grad_rho_a_b(N_states), vc_grad_rho_a_2(N_states), vc_grad_rho_b_2(N_states), vc_grad_rho_a_b(N_states))
-
-
- allocate(rho_a(N_states), rho_b(N_states),grad_rho_a(3,N_states),grad_rho_b(3,N_states))
- allocate(grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states), ex(N_states), ec(N_states))
- allocate(contrib_grad_xa(3,N_states),contrib_grad_xb(3,N_states),contrib_grad_ca(3,N_states),contrib_grad_cb(3,N_states))
 
  aos_dsr_vxc_alpha_pbe_w = 0.d0
  aos_dsr_vxc_beta_pbe_w = 0.d0
 
  do istate = 1, N_states
   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)
    weight = final_weight_at_r_vector(i)
-   rho_a(istate) =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
+   rho_a =  one_e_dm_and_grad_alpha_in_r(4,i,istate)
-   rho_b(istate) =  one_e_dm_and_grad_beta_in_r(4,i,istate)
+   rho_b =  one_e_dm_and_grad_beta_in_r(4,i,istate)
-   grad_rho_a(1:3,istate) =  one_e_dm_and_grad_alpha_in_r(1:3,i,istate)
+   grad_rho_a(1:3) =  one_e_dm_and_grad_alpha_in_r(1:3,i,istate)
-   grad_rho_b(1:3,istate) =  one_e_dm_and_grad_beta_in_r(1:3,i,istate)
+   grad_rho_b(1:3) =  one_e_dm_and_grad_beta_in_r(1:3,i,istate)
    grad_rho_a_2 = 0.d0
    grad_rho_b_2 = 0.d0
    grad_rho_a_b = 0.d0
    do m = 1, 3
-    grad_rho_a_2(istate) += grad_rho_a(m,istate) * grad_rho_a(m,istate)
+    grad_rho_a_2 += grad_rho_a(m) * grad_rho_a(m)
-    grad_rho_b_2(istate) += grad_rho_b(m,istate) * grad_rho_b(m,istate)
+    grad_rho_b_2 += grad_rho_b(m) * grad_rho_b(m)
-    grad_rho_a_b(istate) += grad_rho_a(m,istate) * grad_rho_b(m,istate)
+    grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m)
    enddo
 
                              ! inputs
-   call GGA_sr_type_functionals(r,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
+   call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,                 &  ! outputs exchange
                              ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, &  ! outputs correlation
                              ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b  )
-   vx_rho_a(istate) *= weight
+   vx_rho_a *= weight
-   vc_rho_a(istate) *= weight
+   vc_rho_a *= weight
-   vx_rho_b(istate) *= weight
+   vx_rho_b *= weight
-   vc_rho_b(istate) *= weight
+   vc_rho_b *= weight
    do m= 1,3
-    contrib_grad_ca(m,istate) = weight * (2.d0 * vc_grad_rho_a_2(istate) *  grad_rho_a(m,istate) + vc_grad_rho_a_b(istate)  * grad_rho_b(m,istate))
+    contrib_grad_ca(m) = weight * (2.d0 * vc_grad_rho_a_2 *  grad_rho_a(m) + vc_grad_rho_a_b  * grad_rho_b(m) )
-    contrib_grad_xa(m,istate) = weight * (2.d0 * vx_grad_rho_a_2(istate) *  grad_rho_a(m,istate) + vx_grad_rho_a_b(istate)  * grad_rho_b(m,istate))
+    contrib_grad_xa(m) = weight * (2.d0 * vx_grad_rho_a_2 *  grad_rho_a(m) + vx_grad_rho_a_b  * grad_rho_b(m) )
-    contrib_grad_cb(m,istate) = weight * (2.d0 * vc_grad_rho_b_2(istate) *  grad_rho_b(m,istate) + vc_grad_rho_a_b(istate)  * grad_rho_a(m,istate))
+    contrib_grad_cb(m) = weight * (2.d0 * vc_grad_rho_b_2 *  grad_rho_b(m) + vc_grad_rho_a_b  * grad_rho_a(m) )
-    contrib_grad_xb(m,istate) = weight * (2.d0 * vx_grad_rho_b_2(istate) *  grad_rho_b(m,istate) + vx_grad_rho_a_b(istate)  * grad_rho_a(m,istate))
+    contrib_grad_xb(m) = weight * (2.d0 * vx_grad_rho_b_2 *  grad_rho_b(m) + vx_grad_rho_a_b  * grad_rho_a(m) )
    enddo
    do j = 1, ao_num
-    aos_sr_vxc_alpha_pbe_w(j,i,istate) = ( vc_rho_a(istate) + vx_rho_a(istate) ) * aos_in_r_array(j,i)
+    aos_sr_vxc_alpha_pbe_w(j,i,istate) = ( vc_rho_a + vx_rho_a ) * aos_in_r_array(j,i)
-    aos_sr_vxc_beta_pbe_w (j,i,istate) = ( vc_rho_b(istate) + vx_rho_b(istate) ) * aos_in_r_array(j,i)
+    aos_sr_vxc_beta_pbe_w (j,i,istate) = ( vc_rho_b + vx_rho_b ) * aos_in_r_array(j,i)
    enddo
    do j = 1, ao_num
     do m = 1,3
-     aos_dsr_vxc_alpha_pbe_w(j,i,istate) += ( contrib_grad_ca(m,istate) + contrib_grad_xa(m,istate) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
+     aos_dsr_vxc_alpha_pbe_w(j,i,istate) += ( contrib_grad_ca(m) + contrib_grad_xa(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
-     aos_dsr_vxc_beta_pbe_w (j,i,istate) += ( contrib_grad_cb(m,istate) + contrib_grad_xb(m,istate) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
+     aos_dsr_vxc_beta_pbe_w (j,i,istate) += ( contrib_grad_cb(m) + contrib_grad_xb(m) ) * aos_grad_in_r_array_transp_xyz(m,j,i)
     enddo
    enddo
   enddo
@@ -378,20 +378,3 @@ END_PROVIDER
 
 END_PROVIDER
 
- BEGIN_PROVIDER [double precision, potential_xc_alpha_ao_sr_pbe,(ao_num,ao_num,N_states)]
-&BEGIN_PROVIDER [double precision, potential_xc_beta_ao_sr_pbe,(ao_num,ao_num,N_states)]
-   implicit none
- BEGIN_DOC
- ! exchange / correlation potential for alpha / beta electrons  with the Perdew-Burke-Ernzerhof GGA functional 
- END_DOC 
-   integer :: i,j,istate
-   do istate = 1, n_states 
-    do i = 1, ao_num
-     do j = 1, ao_num
-      potential_xc_alpha_ao_sr_pbe(j,i,istate) = pot_sr_scal_xc_alpha_ao_pbe(j,i,istate) + pot_sr_grad_xc_alpha_ao_pbe(j,i,istate) + pot_sr_grad_xc_alpha_ao_pbe(i,j,istate)
-      potential_xc_beta_ao_sr_pbe(j,i,istate)  = pot_sr_scal_xc_beta_ao_pbe(j,i,istate)  + pot_sr_grad_xc_beta_ao_pbe(j,i,istate)  + pot_sr_grad_xc_beta_ao_pbe(i,j,istate)
-     enddo
-    enddo
-   enddo
-
-END_PROVIDER 

src/hartree_fock/10.hf.bats

@@ -11,7 +11,7 @@ function run() {
   qp edit --check
   qp reset --mos
   qp run scf 
-  qp set_frozen_core 
+#  qp set_frozen_core 
   energy="$(ezfio get hartree_fock energy)"
   eq $energy $2 $thresh
 }