From 3b75503bc189527aaaef0f0d03934ec1ef544978 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Fri, 7 Aug 2020 11:47:51 +0200
Subject: [PATCH 01/21] minor changes

---
 src/determinants/psi_energy_mono_elec.irp.f | 30 +++++++++++++++++++++
 src/utils/integration.irp.f                 | 10 ++++++-
 2 files changed, 39 insertions(+), 1 deletion(-)

diff --git a/src/determinants/psi_energy_mono_elec.irp.f b/src/determinants/psi_energy_mono_elec.irp.f
index 74e69160..fc0ffba9 100644
--- a/src/determinants/psi_energy_mono_elec.irp.f
+++ b/src/determinants/psi_energy_mono_elec.irp.f
@@ -27,3 +27,33 @@
   psi_energy_h_core(i) = psi_energy_h_core(i) * accu
  enddo
 END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, v_ne_psi_energy, (N_states) ]
+  implicit none
+  integer :: i
+  integer :: j,k
+  double precision :: tmp(mo_num,mo_num),mono_ints(mo_num,mo_num)
+  BEGIN_DOC
+! v_ne_psi_energy = $\langle \Psi | v_ne |\Psi \rangle$
+!
+! computed using the :c:data:`one_e_dm_mo_alpha` +
+! :c:data:`one_e_dm_mo_beta` and :c:data:`mo_one_e_integrals`
+  END_DOC
+  v_ne_psi_energy = 0.d0
+  do i = 1, N_states
+   do j = 1, mo_num
+    do k = 1, mo_num
+     v_ne_psi_energy(i) += mo_integrals_n_e(k,j) * (one_e_dm_mo_alpha(k,j,i) + one_e_dm_mo_beta(k,j,i))
+    enddo
+   enddo
+  enddo
+ double precision :: accu
+ do i = 1, N_states
+  accu = 0.d0
+  do j = 1, mo_num
+   accu += one_e_dm_mo_alpha(j,j,i) + one_e_dm_mo_beta(j,j,i)
+  enddo
+  accu = (elec_alpha_num + elec_beta_num ) / accu
+  v_ne_psi_energy(i) = v_ne_psi_energy(i) * accu
+ enddo
+END_PROVIDER
diff --git a/src/utils/integration.irp.f b/src/utils/integration.irp.f
index c907e425..da8120a1 100644
--- a/src/utils/integration.irp.f
+++ b/src/utils/integration.irp.f
@@ -30,7 +30,11 @@ subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alph
   ab = alpha * beta
   d_AB = (A_center - B_center) * (A_center - B_center)
   P_center = (alpha * A_center + beta * B_center) * p_inv
-  fact_k = exp(-ab*p_inv * d_AB)
+  if(dabs(ab*p_inv * d_AB).lt.50.d0)then
+   fact_k = exp(-ab*p_inv * d_AB)
+  else
+   fact_k = 0.d0
+  endif
 
   ! Recenter the polynomials P_a and P_b on x
   !DIR$ FORCEINLINE
@@ -78,6 +82,10 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
   !DIR$ FORCEINLINE
   call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center)
   if (fact_k < thresh) then
+    P_center = 0.d0
+    p = 1.d-10
+    P_new = 0.d0
+    iorder = -1
     fact_k = 0.d0
     return
   endif

From 03445e1a6ec05e02b8c5268b2db67e3e3c6a447a Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Fri, 11 Sep 2020 12:53:06 +0200
Subject: [PATCH 02/21] added hcore_guess in tools

---
 src/tools/hcore_guess.irp.f | 3 +++
 1 file changed, 3 insertions(+)
 create mode 100644 src/tools/hcore_guess.irp.f

diff --git a/src/tools/hcore_guess.irp.f b/src/tools/hcore_guess.irp.f
new file mode 100644
index 00000000..87d0cb7d
--- /dev/null
+++ b/src/tools/hcore_guess.irp.f
@@ -0,0 +1,3 @@
+program hcore_guess_prog
+ call hcore_guess
+end

From ee267e27e9b6027ea76a05e2b43271cc8906fd55 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Mon, 21 Sep 2020 15:38:26 +0200
Subject: [PATCH 03/21] minor modifs

---
 config/gfortran_debug.cfg                |  3 +-
 src/basis_correction/print_routine.irp.f |  1 +
 src/dft_utils_in_r/mo_in_r.irp.f         |  3 +-
 src/mo_basis/mos_in_r.irp.f              | 10 ++--
 src/mo_guess/h_core_guess_routine.irp.f  |  4 +-
 src/mu_of_r/mu_of_r_conditions.irp.f     | 67 ++++++++++++++++++++++++
 src/utils/util.irp.f                     | 12 +++++
 7 files changed, 91 insertions(+), 9 deletions(-)

diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg
index 342acae9..926255e0 100644
--- a/config/gfortran_debug.cfg
+++ b/config/gfortran_debug.cfg
@@ -51,7 +51,8 @@ FCFLAGS : -Ofast
 # -g           : Extra debugging information
 #
 [DEBUG]
-FCFLAGS : -g -msse4.2  -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized  -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
+#FCFLAGS : -g -msse4.2  -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized  -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
+FCFLAGS : -g -msse4.2  -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized  -fbacktrace -ffpe-trap=zero,overflow -finit-real=nan
 
 # OpenMP flags
 #################
diff --git a/src/basis_correction/print_routine.irp.f b/src/basis_correction/print_routine.irp.f
index 05fbbf60..3d497552 100644
--- a/src/basis_correction/print_routine.irp.f
+++ b/src/basis_correction/print_routine.irp.f
@@ -75,6 +75,7 @@ subroutine print_basis_correction
   print*,'**************'
   do istate = 1, N_states
     write(*, '(A29,X,I3,X,A3,X,F16.10)') '  Average mu(r)      , state ',istate,' = ',mu_average_prov(istate)
+    write(*, '(A29,X,I3,X,A3,X,F16.10)') 'mu_average_trans_corr, state ',istate,' = ',mu_average_trans_corr(istate)
   enddo
 
 end
diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f
index 81863c3a..fa3ec73c 100644
--- a/src/dft_utils_in_r/mo_in_r.irp.f
+++ b/src/dft_utils_in_r/mo_in_r.irp.f
@@ -79,7 +79,7 @@
  END_DOC
  integer :: m
  integer  :: i,j
- mos_grad_in_r_array = 0.d0
+ mos_grad_in_r_array_tranp = 0.d0
  do i = 1, n_points_final_grid
   do j = 1, mo_num
    do m = 1, 3
@@ -126,4 +126,3 @@
  enddo
  END_PROVIDER
 
-
diff --git a/src/mo_basis/mos_in_r.irp.f b/src/mo_basis/mos_in_r.irp.f
index 049db8aa..ffcd5a9a 100644
--- a/src/mo_basis/mos_in_r.irp.f
+++ b/src/mo_basis/mos_in_r.irp.f
@@ -12,18 +12,18 @@ subroutine give_all_mos_and_grad_at_r(r,mos_array,mos_grad_array)
  implicit none
  double precision, intent(in) :: r(3)
  double precision, intent(out) :: mos_array(mo_num)
- double precision, intent(out) :: mos_grad_array(mo_num,3)
+ double precision, intent(out) :: mos_grad_array(3,mo_num)
  integer :: i,j,k
- double precision :: aos_array(ao_num),aos_grad_array(ao_num,3)
+ double precision :: aos_array(ao_num),aos_grad_array(3,ao_num)
  call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
  mos_array=0d0
  mos_grad_array=0d0
  do j = 1, mo_num
   do k=1, ao_num
    mos_array(j) += mo_coef(k,j)*aos_array(k)
-   mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1)
-   mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2)
-   mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3)
+   mos_grad_array(1,j) += mo_coef(k,j)*aos_grad_array(1,k)
+   mos_grad_array(2,j) += mo_coef(k,j)*aos_grad_array(2,k)
+   mos_grad_array(3,j) += mo_coef(k,j)*aos_grad_array(3,k)
   enddo
  enddo
 end
diff --git a/src/mo_guess/h_core_guess_routine.irp.f b/src/mo_guess/h_core_guess_routine.irp.f
index 8fc3f6f2..9c5e76ff 100644
--- a/src/mo_guess/h_core_guess_routine.irp.f
+++ b/src/mo_guess/h_core_guess_routine.irp.f
@@ -9,5 +9,7 @@ subroutine hcore_guess
                                      size(mo_one_e_integrals,1),  &
                                      size(mo_one_e_integrals,2),label,1,.false.)
   call save_mos
-  SOFT_TOUCH mo_coef mo_label
+!  SOFT_TOUCH mo_coef mo_label
+  TOUCH mo_coef mo_label
+  print*,'mo_one_e_integrals(1,1) = ',mo_one_e_integrals(1,1)
 end
diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f
index 05a2a4d7..b35a2a8a 100644
--- a/src/mu_of_r/mu_of_r_conditions.irp.f
+++ b/src/mu_of_r/mu_of_r_conditions.irp.f
@@ -28,6 +28,8 @@
     mu_of_r_prov(ipoint,istate) =  mu_of_r_hf(ipoint)
    else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then
     mu_of_r_prov(ipoint,istate) =  mu_of_r_psi_cas(ipoint,istate)
+   else if(mu_of_r_potential.EQ."transcorr")then
+    mu_of_r_prov(ipoint,istate) = mu_of_r_transcorr(ipoint,istate)
    else 
     print*,'you requested the following mu_of_r_potential'
     print*,mu_of_r_potential
@@ -124,6 +126,47 @@
  print*,'Time to provide mu_of_r_psi_cas = ',wall1-wall0
  END_PROVIDER 
 
+ BEGIN_PROVIDER [double precision, mu_of_r_transcorr, (n_points_final_grid,N_states) ]
+ implicit none 
+ BEGIN_DOC
+ ! mu(r) computed with a wave function developped in an active space
+ !
+ ! corresponds to \sqrt(pi) * (W(0) + 1/4)/3
+ !
+ ! !!!!!! 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 :: wall0,wall1,f_psi,on_top,w_psi,sqpi
+ print*,'providing mu_of_r_transcorr ...'
+ call wall_time(wall0)
+ sqpi = dsqrt(dacos(-1.d0))
+
+ provide f_psi_cas_ab
+ !$OMP PARALLEL DO &
+ !$OMP DEFAULT (NONE)  &
+ !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) & 
+ !$OMP SHARED (n_points_final_grid,mu_of_r_transcorr,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states) 
+ do istate = 1, N_states
+  do ipoint = 1, n_points_final_grid
+   f_psi  = f_psi_cas_ab(ipoint,istate) 
+   on_top = on_top_cas_mu_r(ipoint,istate)
+   if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then
+     w_psi   = 1.d+10
+   else 
+     w_psi  = f_psi /  on_top
+   endif
+   mu_of_r_transcorr(ipoint,istate) = (0.25d0 + w_psi) * sqpi / 3.d0
+  enddo
+ enddo
+ !$OMP END PARALLEL DO
+ call wall_time(wall1)
+ print*,'Time to provide mu_of_r_transcorr = ',wall1-wall0
+ END_PROVIDER 
+
 
  BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)]
  implicit none
@@ -148,3 +191,27 @@
   mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate)
  enddo
  END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, mu_average_trans_corr, (N_states)]
+ implicit none
+ BEGIN_DOC
+ ! average value of mu(r) weighted with the total one-e density and divised 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
+ mu_average_trans_corr = 0.d0
+ do istate = 1, N_states
+  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)
+   if(mu_of_r_transcorr(ipoint,istate).gt.1.d+09)cycle
+   mu_average_trans_corr(istate) += mu_of_r_transcorr(ipoint,istate) * weight * density
+  enddo
+  mu_average_trans_corr(istate) = mu_average_trans_corr(istate) / elec_num_grid_becke(istate)
+ enddo
+ END_PROVIDER 
diff --git a/src/utils/util.irp.f b/src/utils/util.irp.f
index 1b01a1ec..cfb42fd1 100644
--- a/src/utils/util.irp.f
+++ b/src/utils/util.irp.f
@@ -1,3 +1,15 @@
+double precision function derf_mu_x(mu,x)
+ implicit none
+ include 'utils/constants.include.F'
+ double precision, intent(in) :: mu,x
+  if(dabs(x).gt.1.d-6)then
+   derf_mu_x = derf(mu * x)/x
+  else
+   derf_mu_x =  inv_sq_pi * 2.d0 * mu * (1.d0 - mu*mu*x*x/3.d0)                                                           
+  endif                      
+end    
+
+
 double precision function binom_func(i,j)
   implicit none
   BEGIN_DOC

From 01b1ee327392efbab00af9ff0badef73f9f91ca9 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Wed, 7 Oct 2020 11:03:23 +0200
Subject: [PATCH 04/21] fixed bug in laplacians

---
 src/dft_utils_in_r/ao_in_r.irp.f | 39 --------------------------------
 src/dft_utils_in_r/mo_in_r.irp.f | 21 +++++++++++++++--
 src/mo_basis/mos_in_r.irp.f      | 16 ++++++-------
 3 files changed, 27 insertions(+), 49 deletions(-)

diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f
index 4b1526dd..2892d2e2 100644
--- a/src/dft_utils_in_r/ao_in_r.irp.f
+++ b/src/dft_utils_in_r/ao_in_r.irp.f
@@ -121,42 +121,3 @@
  END_PROVIDER
 
 
- BEGIN_PROVIDER[double precision, aos_lapl_in_r_array_transp, (n_points_final_grid,ao_num,3)]
- implicit none
- !
- ! aos_lapl_in_r_array_transp(i,j,k)   = value of the kth component of the laplacian of jth ao on the ith grid point
- !
- ! k = 1 : x, k= 2, y, k  3, z
- integer :: i,j,m
- do m = 1, 3
-  do i = 1, n_points_final_grid
-   do j = 1, ao_num
-    aos_lapl_in_r_array_transp(i,j,m) = aos_lapl_in_r_array(j,i,m) 
-   enddo
-  enddo
- enddo
- END_PROVIDER
-
-
- BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)]
-&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)]
- implicit none
- BEGIN_DOC
- ! aos_in_r_array_per_atom(i,j,k)        = value of the ith ao on the jth grid point attached on the kth atom 
- END_DOC
- integer :: i,j,k
- double precision :: aos_array(ao_num), r(3)
- do k = 1, nucl_num
-  do i = 1, n_pts_per_atom(k)
-   r(1) = final_grid_points_per_atom(1,i,k)
-   r(2) = final_grid_points_per_atom(2,i,k)
-   r(3) = final_grid_points_per_atom(3,i,k)
-   call give_all_aos_at_r(r,aos_array)
-   do j = 1, ao_num
-    aos_in_r_array_per_atom(j,i,k) = aos_array(j)
-    aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j)
-   enddo
-  enddo
- enddo
- END_PROVIDER
-
diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f
index fa3ec73c..d4afcf76 100644
--- a/src/dft_utils_in_r/mo_in_r.irp.f
+++ b/src/dft_utils_in_r/mo_in_r.irp.f
@@ -115,8 +115,6 @@
  BEGIN_DOC
  ! mos_lapl_in_r_array(i,j,k)          = value of the kth component of the laplacian of ith mo on the jth grid point
  !
- ! mos_lapl_in_r_array_transp(i,j,k)   = value of the kth component of the laplacian of jth mo on the ith grid point
- !
  ! k = 1 : x, k= 2, y, k  3, z
  END_DOC
  integer :: m
@@ -126,3 +124,22 @@
  enddo
  END_PROVIDER
 
+
+ BEGIN_PROVIDER[double precision, mos_lapl_in_r_array_tranp,(3,mo_num,n_points_final_grid)]
+ implicit none
+ BEGIN_DOC
+ ! mos_lapl_in_r_array_transp(i,j,k)   = value of the kth component of the laplient of jth mo on the ith grid point
+ !
+ ! k = 1 : x, k= 2, y, k  3, z
+ END_DOC
+ integer :: m
+ integer  :: i,j
+ mos_lapl_in_r_array_tranp = 0.d0
+ do i = 1, n_points_final_grid
+  do j = 1, mo_num
+   do m = 1, 3
+     mos_lapl_in_r_array_tranp(m,j,i) = mos_lapl_in_r_array(j,i,m)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER
diff --git a/src/mo_basis/mos_in_r.irp.f b/src/mo_basis/mos_in_r.irp.f
index ffcd5a9a..ee2795d0 100644
--- a/src/mo_basis/mos_in_r.irp.f
+++ b/src/mo_basis/mos_in_r.irp.f
@@ -33,9 +33,9 @@ subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_la
  implicit none
  double precision, intent(in) :: r(3)
  double precision, intent(out) :: mos_array(mo_num)
- double precision, intent(out) :: mos_grad_array(mo_num,3),mos_lapl_array(mo_num,3)
+ double precision, intent(out) :: mos_grad_array(3,mo_num),mos_lapl_array(3,mo_num)
  integer :: i,j,k
- double precision :: aos_array(ao_num),aos_grad_array(ao_num,3),aos_lapl_array(ao_num,3)
+ double precision :: aos_array(ao_num),aos_grad_array(3,ao_num),aos_lapl_array(3,ao_num)
  call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
  mos_array      = 0.d0
  mos_grad_array = 0.d0
@@ -43,12 +43,12 @@ subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_la
  do j = 1, mo_num
   do k=1, ao_num
    mos_array(j)        += mo_coef(k,j) * aos_array(k)
-   mos_grad_array(j,1) += mo_coef(k,j) * aos_grad_array(k,1)
-   mos_grad_array(j,2) += mo_coef(k,j) * aos_grad_array(k,2)
-   mos_grad_array(j,3) += mo_coef(k,j) * aos_grad_array(k,3)
-   mos_lapl_array(j,1) += mo_coef(k,j) * aos_lapl_array(k,1)
-   mos_lapl_array(j,2) += mo_coef(k,j) * aos_lapl_array(k,2)
-   mos_lapl_array(j,3) += mo_coef(k,j) * aos_lapl_array(k,3)
+   mos_grad_array(1,j) += mo_coef(k,j) * aos_grad_array(1,k)
+   mos_grad_array(2,j) += mo_coef(k,j) * aos_grad_array(2,k)
+   mos_grad_array(3,j) += mo_coef(k,j) * aos_grad_array(3,k)
+   mos_lapl_array(1,j) += mo_coef(k,j) * aos_lapl_array(1,k)
+   mos_lapl_array(2,j) += mo_coef(k,j) * aos_lapl_array(2,k)
+   mos_lapl_array(3,j) += mo_coef(k,j) * aos_lapl_array(3,k)
   enddo
  enddo
 end

From c545098af083a423e2b4f6ce7597c1d9ce566c38 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Mon, 19 Oct 2020 12:20:14 +0200
Subject: [PATCH 05/21] added tools/huckel_guess.irp.f

---
 src/tools/huckel_guess.irp.f | 5 +++++
 1 file changed, 5 insertions(+)
 create mode 100644 src/tools/huckel_guess.irp.f

diff --git a/src/tools/huckel_guess.irp.f b/src/tools/huckel_guess.irp.f
new file mode 100644
index 00000000..5ec37df4
--- /dev/null
+++ b/src/tools/huckel_guess.irp.f
@@ -0,0 +1,5 @@
+program pouet
+ implicit none 
+ call huckel_guess
+
+end

From 2c33bca408aaca9003d2af93c8ba63721cefcb88 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Wed, 28 Oct 2020 11:24:32 +0100
Subject: [PATCH 06/21] minor changes

---
 src/dft_utils_in_r/ao_in_r.irp.f | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f
index 2892d2e2..6abda550 100644
--- a/src/dft_utils_in_r/ao_in_r.irp.f
+++ b/src/dft_utils_in_r/ao_in_r.irp.f
@@ -120,4 +120,18 @@
  !$OMP END PARALLEL DO
  END_PROVIDER
 
+ BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_bis, (n_points_final_grid,ao_num,3)]
+ implicit none
+ integer :: i,j,m
+ double precision :: aos_array(ao_num), r(3)
+ double precision :: aos_grad_array(3,ao_num)
+ do m = 1, 3
+  do j = 1, ao_num
+   do i = 1, n_points_final_grid
+    aos_grad_in_r_array_transp_bis(i,j,m) = aos_grad_in_r_array(j,i,m)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER
+
 

From bd8fcdb03030af73d5f4bef74a6ccd0a8899702d Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Mon, 2 Nov 2020 17:24:35 +0100
Subject: [PATCH 07/21] removed stupid bug in utils/one_e_integration.irp.f

---
 src/ao_one_e_ints/ao_overlap.irp.f |  7 +++++++
 src/utils/one_e_integration.irp.f  | 22 +++++++++++-----------
 2 files changed, 18 insertions(+), 11 deletions(-)

diff --git a/src/ao_one_e_ints/ao_overlap.irp.f b/src/ao_one_e_ints/ao_overlap.irp.f
index 11c95e42..d9061d67 100644
--- a/src/ao_one_e_ints/ao_overlap.irp.f
+++ b/src/ao_one_e_ints/ao_overlap.irp.f
@@ -54,6 +54,13 @@
         call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
         c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)
         ao_overlap(i,j) += c * overlap
+        if(isnan(ao_overlap(i,j)))then
+         print*,'i,j',i,j
+         print*,'l,n',l,n
+         print*,'c,overlap',c,overlap
+         print*,overlap_x,overlap_y,overlap_z
+         stop
+        endif
         ao_overlap_x(i,j) += c * overlap_x
         ao_overlap_y(i,j) += c * overlap_y
         ao_overlap_z(i,j) += c * overlap_z
diff --git a/src/utils/one_e_integration.irp.f b/src/utils/one_e_integration.irp.f
index 97eef89d..e7b3e051 100644
--- a/src/utils/one_e_integration.irp.f
+++ b/src/utils/one_e_integration.irp.f
@@ -15,10 +15,10 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
   call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,&
       beta,power_A,power_B,A_center,B_center,dim)
 
-! if(fact_p.lt.0.000001d0)then
-!   overlap_gaussian_x = 0.d0
-!   return
-! endif
+  if(fact_p.lt.1.d-20)then
+    overlap_gaussian_x = 0.d0
+    return
+  endif
 
   overlap_gaussian_x = 0.d0
   integer                        :: i
@@ -53,13 +53,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
   integer                        :: iorder_p(3)
 
   call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
-!  if(fact_p.lt.1d-20)then
-!    overlap_x = 0.d0
-!    overlap_y = 0.d0
-!    overlap_z = 0.d0
-!    overlap = 0.d0
-!    return
-!  endif
+   if(fact_p.lt.1d-20)then
+     overlap_x = 0.d0
+     overlap_y = 0.d0
+     overlap_z = 0.d0
+     overlap = 0.d0
+     return
+   endif
   integer                        :: nmax
   double precision               :: F_integral
   nmax = maxval(iorder_p)

From bdbb08207956051bb427051f184449985628fd23 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Sun, 8 Nov 2020 17:53:37 +0100
Subject: [PATCH 08/21] fixed a little bun in src/utils/one_e_integration.irp.f

---
 src/utils/one_e_integration.irp.f | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/utils/one_e_integration.irp.f b/src/utils/one_e_integration.irp.f
index 97eef89d..2c597104 100644
--- a/src/utils/one_e_integration.irp.f
+++ b/src/utils/one_e_integration.irp.f
@@ -53,13 +53,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
   integer                        :: iorder_p(3)
 
   call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
-!  if(fact_p.lt.1d-20)then
-!    overlap_x = 0.d0
-!    overlap_y = 0.d0
-!    overlap_z = 0.d0
-!    overlap = 0.d0
-!    return
-!  endif
+   if(fact_p.lt.1d-20)then
+     overlap_x = 1.d-10
+     overlap_y = 1.d-10
+     overlap_z = 1.d-10
+     overlap = 1.d-10
+     return
+   endif
   integer                        :: nmax
   double precision               :: F_integral
   nmax = maxval(iorder_p)

From 6652da23f57dbb84f490913346ad31c13e628f68 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Sun, 8 Nov 2020 18:26:54 +0100
Subject: [PATCH 09/21] removed symmetry bug in a-b-2-RDM : warning it has no
 impacts on any quantity as long as all the elements of the two-rdm were used

---
 src/two_rdm_routines/update_rdm.irp.f         | 20 +++++++++++++++++--
 .../update_state_av_rdm.irp.f                 | 16 +++++++++++++--
 2 files changed, 32 insertions(+), 4 deletions(-)

diff --git a/src/two_rdm_routines/update_rdm.irp.f b/src/two_rdm_routines/update_rdm.irp.f
index 4d74280e..662b3ee9 100644
--- a/src/two_rdm_routines/update_rdm.irp.f
+++ b/src/two_rdm_routines/update_rdm.irp.f
@@ -356,12 +356,20 @@
       h2 = list_orb_reverse(h2)
       nkeys += 1
       do istate = 1, N_st
-       values(istate,nkeys) = c_1(istate) * phase
+       values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
       enddo
       keys(1,nkeys) = h1
       keys(2,nkeys) = h2
       keys(3,nkeys) = p1
       keys(4,nkeys) = h2
+      nkeys += 1
+      do istate = 1, N_st
+       values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
+      enddo
+      keys(1,nkeys) = h2
+      keys(2,nkeys) = h1
+      keys(3,nkeys) = h2
+      keys(4,nkeys) = p1
      enddo 
    else 
     ! Mono beta
@@ -377,12 +385,20 @@
       h2 = list_orb_reverse(h2)
       nkeys += 1
       do istate = 1, N_st
-       values(istate,nkeys) = c_1(istate) * phase
+       values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
       enddo
       keys(1,nkeys) = h1
       keys(2,nkeys) = h2
       keys(3,nkeys) = p1
       keys(4,nkeys) = h2
+      nkeys += 1
+      do istate = 1, N_st
+       values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
+      enddo
+      keys(1,nkeys) = h2
+      keys(2,nkeys) = h1
+      keys(3,nkeys) = h2
+      keys(4,nkeys) = p1
      enddo 
    endif
   else if(spin_trace)then
diff --git a/src/two_rdm_routines/update_state_av_rdm.irp.f b/src/two_rdm_routines/update_state_av_rdm.irp.f
index 35024331..bf7b6148 100644
--- a/src/two_rdm_routines/update_state_av_rdm.irp.f
+++ b/src/two_rdm_routines/update_state_av_rdm.irp.f
@@ -327,11 +327,17 @@
       if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
       h2 = list_orb_reverse(h2)
       nkeys += 1
-      values(nkeys) = c_1 * phase
+      values(nkeys) = 0.5d0 * c_1 * phase
       keys(1,nkeys) = h1
       keys(2,nkeys) = h2
       keys(3,nkeys) = p1
       keys(4,nkeys) = h2
+      nkeys += 1
+      values(nkeys) = 0.5d0 * c_1 * phase
+      keys(1,nkeys) = h2
+      keys(2,nkeys) = h1
+      keys(3,nkeys) = h2
+      keys(4,nkeys) = p1
      enddo 
    else 
     ! Mono beta
@@ -346,11 +352,17 @@
       if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
       h2 = list_orb_reverse(h2)
       nkeys += 1
-      values(nkeys) = c_1 * phase
+      values(nkeys) = 0.5d0 * c_1 * phase
       keys(1,nkeys) = h1
       keys(2,nkeys) = h2
       keys(3,nkeys) = p1
       keys(4,nkeys) = h2
+      nkeys += 1
+      values(nkeys) = 0.5d0 * c_1 * phase
+      keys(1,nkeys) = h2
+      keys(2,nkeys) = h1
+      keys(3,nkeys) = h2
+      keys(4,nkeys) = p1
      enddo 
    endif
   else if(spin_trace)then

From aa0c44959c69d0355e35c09d2798b3c053122305 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Mon, 9 Nov 2020 18:30:10 +0100
Subject: [PATCH 10/21] removed stupid dependency in basis_set_correction

---
 src/basis_correction/print_routine.irp.f | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/basis_correction/print_routine.irp.f b/src/basis_correction/print_routine.irp.f
index 3d497552..05fbbf60 100644
--- a/src/basis_correction/print_routine.irp.f
+++ b/src/basis_correction/print_routine.irp.f
@@ -75,7 +75,6 @@ subroutine print_basis_correction
   print*,'**************'
   do istate = 1, N_states
     write(*, '(A29,X,I3,X,A3,X,F16.10)') '  Average mu(r)      , state ',istate,' = ',mu_average_prov(istate)
-    write(*, '(A29,X,I3,X,A3,X,F16.10)') 'mu_average_trans_corr, state ',istate,' = ',mu_average_trans_corr(istate)
   enddo
 
 end

From 5b8580fe2dc993a469f3376d4e01ee7b31b7a552 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Sat, 2 Jan 2021 15:40:03 +0100
Subject: [PATCH 11/21] some cleaning

---
 src/dft_utils_in_r/ao_in_r.irp.f        |  4 ++
 src/mo_guess/h_core_guess_routine.irp.f |  2 -
 src/mu_of_r/mu_of_r_conditions.irp.f    | 66 -------------------------
 3 files changed, 4 insertions(+), 68 deletions(-)

diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f
index 6abda550..a98afd5d 100644
--- a/src/dft_utils_in_r/ao_in_r.irp.f
+++ b/src/dft_utils_in_r/ao_in_r.irp.f
@@ -122,6 +122,10 @@
 
  BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_bis, (n_points_final_grid,ao_num,3)]
  implicit none
+ BEGIN_DOC
+! Transposed gradients 
+! 
+ END_DOC
  integer :: i,j,m
  double precision :: aos_array(ao_num), r(3)
  double precision :: aos_grad_array(3,ao_num)
diff --git a/src/mo_guess/h_core_guess_routine.irp.f b/src/mo_guess/h_core_guess_routine.irp.f
index 9c5e76ff..246dfef2 100644
--- a/src/mo_guess/h_core_guess_routine.irp.f
+++ b/src/mo_guess/h_core_guess_routine.irp.f
@@ -9,7 +9,5 @@ subroutine hcore_guess
                                      size(mo_one_e_integrals,1),  &
                                      size(mo_one_e_integrals,2),label,1,.false.)
   call save_mos
-!  SOFT_TOUCH mo_coef mo_label
   TOUCH mo_coef mo_label
-  print*,'mo_one_e_integrals(1,1) = ',mo_one_e_integrals(1,1)
 end
diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f
index b35a2a8a..148c65b3 100644
--- a/src/mu_of_r/mu_of_r_conditions.irp.f
+++ b/src/mu_of_r/mu_of_r_conditions.irp.f
@@ -28,8 +28,6 @@
     mu_of_r_prov(ipoint,istate) =  mu_of_r_hf(ipoint)
    else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then
     mu_of_r_prov(ipoint,istate) =  mu_of_r_psi_cas(ipoint,istate)
-   else if(mu_of_r_potential.EQ."transcorr")then
-    mu_of_r_prov(ipoint,istate) = mu_of_r_transcorr(ipoint,istate)
    else 
     print*,'you requested the following mu_of_r_potential'
     print*,mu_of_r_potential
@@ -126,47 +124,6 @@
  print*,'Time to provide mu_of_r_psi_cas = ',wall1-wall0
  END_PROVIDER 
 
- BEGIN_PROVIDER [double precision, mu_of_r_transcorr, (n_points_final_grid,N_states) ]
- implicit none 
- BEGIN_DOC
- ! mu(r) computed with a wave function developped in an active space
- !
- ! corresponds to \sqrt(pi) * (W(0) + 1/4)/3
- !
- ! !!!!!! 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 :: wall0,wall1,f_psi,on_top,w_psi,sqpi
- print*,'providing mu_of_r_transcorr ...'
- call wall_time(wall0)
- sqpi = dsqrt(dacos(-1.d0))
-
- provide f_psi_cas_ab
- !$OMP PARALLEL DO &
- !$OMP DEFAULT (NONE)  &
- !$OMP PRIVATE (ipoint,f_psi,on_top,w_psi,istate) & 
- !$OMP SHARED (n_points_final_grid,mu_of_r_transcorr,f_psi_cas_ab,on_top_cas_mu_r,sqpi,N_states) 
- do istate = 1, N_states
-  do ipoint = 1, n_points_final_grid
-   f_psi  = f_psi_cas_ab(ipoint,istate) 
-   on_top = on_top_cas_mu_r(ipoint,istate)
-   if(on_top.le.1.d-12.or.f_psi.le.0.d0.or.f_psi * on_top.lt.0.d0)then
-     w_psi   = 1.d+10
-   else 
-     w_psi  = f_psi /  on_top
-   endif
-   mu_of_r_transcorr(ipoint,istate) = (0.25d0 + w_psi) * sqpi / 3.d0
-  enddo
- enddo
- !$OMP END PARALLEL DO
- call wall_time(wall1)
- print*,'Time to provide mu_of_r_transcorr = ',wall1-wall0
- END_PROVIDER 
-
 
  BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)]
  implicit none
@@ -192,26 +149,3 @@
  enddo
  END_PROVIDER 
 
- BEGIN_PROVIDER [double precision, mu_average_trans_corr, (N_states)]
- implicit none
- BEGIN_DOC
- ! average value of mu(r) weighted with the total one-e density and divised 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
- mu_average_trans_corr = 0.d0
- do istate = 1, N_states
-  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)
-   if(mu_of_r_transcorr(ipoint,istate).gt.1.d+09)cycle
-   mu_average_trans_corr(istate) += mu_of_r_transcorr(ipoint,istate) * weight * density
-  enddo
-  mu_average_trans_corr(istate) = mu_average_trans_corr(istate) / elec_num_grid_becke(istate)
- enddo
- END_PROVIDER 

From 970f846a4d4615d7cb73853f9c8b7d2b1e2802c9 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Thu, 7 Jan 2021 14:20:11 +0100
Subject: [PATCH 12/21] changed the read/write two_rdm

---
 src/two_body_rdm/act_2_rdm.irp.f  | 24 ++++++++++++++++--------
 src/two_body_rdm/io_two_rdm.irp.f | 29 +++++++++++++++++++++++++++++
 2 files changed, 45 insertions(+), 8 deletions(-)
 create mode 100644 src/two_body_rdm/io_two_rdm.irp.f

diff --git a/src/two_body_rdm/act_2_rdm.irp.f b/src/two_body_rdm/act_2_rdm.irp.f
index e3265572..41b28aea 100644
--- a/src/two_body_rdm/act_2_rdm.irp.f
+++ b/src/two_body_rdm/act_2_rdm.irp.f
@@ -22,6 +22,8 @@
  END_DOC 
  integer :: ispin
  double precision :: wall_1, wall_2
+ character*(128) :: name_file 
+ name_file = 'act_2_rdm_ab_mo'
  ! condition for alpha/beta spin
  print*,''
  print*,'Providing act_2_rdm_ab_mo '
@@ -31,13 +33,13 @@
  call wall_time(wall_1)
  if(read_two_body_rdm_ab)then
   print*,'Reading act_2_rdm_ab_mo from disk ...'
-  call ezfio_get_two_body_rdm_two_rdm_ab_disk(act_2_rdm_ab_mo)
+  call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_ab_mo,name_file)
  else 
   call orb_range_2_rdm_openmp(act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  endif
  if(write_two_body_rdm_ab)then
   print*,'Writing act_2_rdm_ab_mo on disk ...'
-  call ezfio_set_two_body_rdm_two_rdm_ab_disk(act_2_rdm_ab_mo)
+  call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_ab_mo,name_file)
   call ezfio_set_two_body_rdm_io_two_body_rdm_ab("Read")
  endif
  call wall_time(wall_2)
@@ -63,18 +65,20 @@
  ! condition for alpha/beta spin
  print*,''
  print*,'Providing act_2_rdm_aa_mo '
+ character*(128) :: name_file 
+ name_file = 'act_2_rdm_aa_mo'
  ispin = 1 
  act_2_rdm_aa_mo = 0.d0
  call wall_time(wall_1)
  if(read_two_body_rdm_aa)then
   print*,'Reading act_2_rdm_aa_mo from disk ...'
-  call ezfio_get_two_body_rdm_two_rdm_aa_disk(act_2_rdm_aa_mo)
+  call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_aa_mo,name_file)
  else 
   call orb_range_2_rdm_openmp(act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  endif
  if(write_two_body_rdm_aa)then
   print*,'Writing act_2_rdm_aa_mo on disk ...'
-  call ezfio_set_two_body_rdm_two_rdm_aa_disk(act_2_rdm_aa_mo)
+  call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_aa_mo,name_file)
   call ezfio_set_two_body_rdm_io_two_body_rdm_aa("Read")
  endif
 
@@ -101,18 +105,20 @@
  ! condition for beta/beta spin
  print*,''
  print*,'Providing act_2_rdm_bb_mo '
+ character*(128) :: name_file 
+ name_file = 'act_2_rdm_bb_mo'
  ispin = 2 
  act_2_rdm_bb_mo = 0.d0
  call wall_time(wall_1)
  if(read_two_body_rdm_bb)then
   print*,'Reading act_2_rdm_bb_mo from disk ...'
-  call ezfio_get_two_body_rdm_two_rdm_bb_disk(act_2_rdm_bb_mo)
+  call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_bb_mo,name_file)
  else 
   call orb_range_2_rdm_openmp(act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  endif
  if(write_two_body_rdm_bb)then
   print*,'Writing act_2_rdm_bb_mo on disk ...'
-  call ezfio_set_two_body_rdm_two_rdm_bb_disk(act_2_rdm_bb_mo)
+  call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_bb_mo,name_file)
   call ezfio_set_two_body_rdm_io_two_body_rdm_bb("Read")
  endif
 
@@ -138,18 +144,20 @@
  ! condition for beta/beta spin
  print*,''
  print*,'Providing act_2_rdm_spin_trace_mo '
+ character*(128) :: name_file 
+ name_file = 'act_2_rdm_spin_trace_mo'
  ispin = 4 
  act_2_rdm_spin_trace_mo = 0.d0
  call wall_time(wall_1)
  if(read_two_body_rdm_spin_trace)then
   print*,'Reading act_2_rdm_spin_trace_mo from disk ...'
-  call ezfio_get_two_body_rdm_two_rdm_spin_trace_disk(act_2_rdm_spin_trace_mo)
+  call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_spin_trace_mo,name_file)
  else 
   call orb_range_2_rdm_openmp(act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
  endif
  if(write_two_body_rdm_spin_trace)then
   print*,'Writing act_2_rdm_spin_trace_mo on disk ...'
-  call ezfio_set_two_body_rdm_two_rdm_spin_trace_disk(act_2_rdm_spin_trace_mo)
+  call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_spin_trace_mo,name_file)
   call ezfio_set_two_body_rdm_io_two_body_rdm_spin_trace("Read")
  endif
 
diff --git a/src/two_body_rdm/io_two_rdm.irp.f b/src/two_body_rdm/io_two_rdm.irp.f
new file mode 100644
index 00000000..f7008ca9
--- /dev/null
+++ b/src/two_body_rdm/io_two_rdm.irp.f
@@ -0,0 +1,29 @@
+subroutine write_array_two_rdm(n_orb,nstates,array_tmp,name_file)
+ implicit none
+ integer, intent(in) :: n_orb,nstates
+ character*(128),  intent(in) :: name_file 
+ double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,nstates)
+
+ character*(128)                :: output
+ integer                        :: i_unit_output,getUnitAndOpen
+ PROVIDE ezfio_filename                                                                                                  
+ output=trim(ezfio_filename)//'/work/'//trim(name_file)
+ i_unit_output = getUnitAndOpen(output,'W')
+ write(i_unit_output)array_tmp
+ close(unit=i_unit_output)
+end
+
+subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file)
+ implicit none
+ character*(128)                :: output
+ integer                        :: i_unit_output,getUnitAndOpen
+ integer, intent(in) :: n_orb,nstates
+ character*(128),  intent(in)  :: name_file 
+ double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,N_states)
+ PROVIDE ezfio_filename                                                                                                  
+ output=trim(ezfio_filename)//'/work/'//trim(name_file)
+ i_unit_output = getUnitAndOpen(output,'R')
+ read(i_unit_output)array_tmp
+ close(unit=i_unit_output)
+end
+

From 7874e4f0373addc220e153d4fee362469629f850 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Thu, 7 Jan 2021 14:24:44 +0100
Subject: [PATCH 13/21] removed useless arrays for two_rdm

---
 src/two_body_rdm/EZFIO.cfg | 24 ------------------------
 1 file changed, 24 deletions(-)

diff --git a/src/two_body_rdm/EZFIO.cfg b/src/two_body_rdm/EZFIO.cfg
index 4ca39d73..d73f704d 100644
--- a/src/two_body_rdm/EZFIO.cfg
+++ b/src/two_body_rdm/EZFIO.cfg
@@ -1,45 +1,21 @@
-[two_rdm_ab_disk]  
-type: double precision
-doc: active part of the two body rdm alpha/beta stored on disk
-interface: ezfio
-size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
-
 [io_two_body_rdm_ab]
 type: Disk_access
 doc: Read/Write the active part of the two-body rdm for alpha/beta electrons from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
 default: None
 
-[two_rdm_aa_disk]  
-type: double precision
-doc: active part of the two body rdm alpha/alpha stored on disk
-interface: ezfio
-size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
-
 [io_two_body_rdm_aa]
 type: Disk_access
 doc: Read/Write the active part of the two-body rdm for alpha/alpha electrons from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
 default: None
 
-[two_rdm_bb_disk]  
-type: double precision
-doc: active part of the two body rdm beta/beta stored on disk
-interface: ezfio
-size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
-
 [io_two_body_rdm_bb]
 type: Disk_access
 doc: Read/Write the active part of the two-body rdm for beta/beta electrons from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
 default: None
 
-[two_rdm_spin_trace_disk]  
-type: double precision
-doc: active part of the two body rdm spin trace stored on disk
-interface: ezfio
-size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
-
 [io_two_body_rdm_spin_trace]
 type: Disk_access
 doc: Read/Write the active part of the two-body rdm for spin trace electrons from/to disk [ Write | Read | None ]

From e5da77577fdaa29f2990256d3a245e035f2feb60 Mon Sep 17 00:00:00 2001
From: bpradines <barth371@hotmail.com>
Date: Thu, 14 Jan 2021 11:05:03 +0100
Subject: [PATCH 14/21] removed stupid dependency on the AO integrals in
 basis_correction

---
 src/basis_correction/basis_correction.irp.f | 18 +++---------------
 1 file changed, 3 insertions(+), 15 deletions(-)

diff --git a/src/basis_correction/basis_correction.irp.f b/src/basis_correction/basis_correction.irp.f
index 8d3264bc..a7ea7244 100644
--- a/src/basis_correction/basis_correction.irp.f
+++ b/src/basis_correction/basis_correction.irp.f
@@ -7,22 +7,10 @@ program basis_correction
   touch read_wf
   no_core_density = .True.
   touch no_core_density
-  provide ao_two_e_integrals_in_map
+  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
-!  call print_e_b
 end
 
-subroutine print_e_b
- implicit none
-  print *, 'Hello world'
-  print*,'ecmd_lda_mu_of_r             = ',ecmd_lda_mu_of_r
-  print*,'ecmd_pbe_ueg_mu_of_r         = ',ecmd_pbe_ueg_mu_of_r
-  print*,'ecmd_pbe_ueg_eff_xi_mu_of_r  = ',ecmd_pbe_ueg_eff_xi_mu_of_r
-  print*,''
-  print*,'psi_energy + E^B_LDA         = ',psi_energy + ecmd_lda_mu_of_r
-  print*,'psi_energy + E^B_PBE_UEG     = ',psi_energy + ecmd_pbe_ueg_mu_of_r
-  print*,'psi_energy + E^B_PBE_UEG_Xi  = ',psi_energy + ecmd_pbe_ueg_eff_xi_mu_of_r
-  print*,''
-  print*,'mu_average_prov              = ',mu_average_prov
-end

From acd9192aa4393f41d1c81cea20248d8b6d22083b Mon Sep 17 00:00:00 2001
From: bpradines <barth371@hotmail.com>
Date: Tue, 19 Jan 2021 19:01:52 +0100
Subject: [PATCH 15/21] fixed bug in weird dependence of mu(r)

---
 src/mu_of_r/f_hf_utils.irp.f        | 1 +
 src/mu_of_r/f_psi_i_a_v_utils.irp.f | 3 +++
 2 files changed, 4 insertions(+)

diff --git a/src/mu_of_r/f_hf_utils.irp.f b/src/mu_of_r/f_hf_utils.irp.f
index b89dda18..8480a288 100644
--- a/src/mu_of_r/f_hf_utils.irp.f
+++ b/src/mu_of_r/f_hf_utils.irp.f
@@ -9,6 +9,7 @@ BEGIN_PROVIDER [double precision, two_e_int_hf_f, (n_basis_orb,n_basis_orb,n_max
  END_DOC
  integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
  double precision :: get_two_e_integral
+ PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
  do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1 
   m = list_valence_orb_for_hf(orb_m,1)
   do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2 
diff --git a/src/mu_of_r/f_psi_i_a_v_utils.irp.f b/src/mu_of_r/f_psi_i_a_v_utils.irp.f
index aed054ae..427da199 100644
--- a/src/mu_of_r/f_psi_i_a_v_utils.irp.f
+++ b/src/mu_of_r/f_psi_i_a_v_utils.irp.f
@@ -235,6 +235,7 @@ BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act
  END_DOC
  integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
  double precision :: integrals_array(mo_num,mo_num),get_two_e_integral
+ PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
  do orb_m = 1, n_act_orb ! electron 1 
   m = list_act(orb_m)
   do orb_n = 1, n_act_orb ! electron 2 
@@ -264,6 +265,7 @@ BEGIN_PROVIDER [double precision, two_e_int_ia_f, (n_basis_orb,n_basis_orb,n_ina
  END_DOC
  integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
  double precision :: integrals_array(mo_num,mo_num),get_two_e_integral
+ PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
  do orb_m = 1, n_act_orb ! electron 1 
   m = list_act(orb_m)
   do orb_n = 1, n_inact_orb ! electron 2 
@@ -293,6 +295,7 @@ BEGIN_PROVIDER [double precision, two_e_int_ii_f, (n_basis_orb,n_basis_orb,n_ina
  END_DOC
  integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
  double precision :: get_two_e_integral,integrals_array(mo_num,mo_num)
+ PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals 
  do orb_m = 1, n_inact_orb ! electron 1 
   m = list_inact(orb_m)
   do orb_n = 1, n_inact_orb ! electron 2 

From d6d8b0586c60cd3184e9e419b476ddaf5f7219d1 Mon Sep 17 00:00:00 2001
From: bpradines <barth371@hotmail.com>
Date: Fri, 26 Feb 2021 10:18:42 +0100
Subject: [PATCH 16/21] added basis_correction/print_su_pbe_ot.irp.f

---
 src/basis_correction/print_su_pbe_ot.irp.f | 28 ++++++++++++++++++++++
 1 file changed, 28 insertions(+)
 create mode 100644 src/basis_correction/print_su_pbe_ot.irp.f

diff --git a/src/basis_correction/print_su_pbe_ot.irp.f b/src/basis_correction/print_su_pbe_ot.irp.f
new file mode 100644
index 00000000..49f90ade
--- /dev/null
+++ b/src/basis_correction/print_su_pbe_ot.irp.f
@@ -0,0 +1,28 @@
+program basis_corr_su_pbe_ot
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  read_wf = .True.
+  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_su_pbe_ot
+ 
+end
+
+subroutine print_su_pbe_ot
+ implicit none
+ integer :: istate 
+ 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
+ do istate = 1, N_states
+  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
+
+end

From 3e551e153742c6effd0680ef97620816a382b5f6 Mon Sep 17 00:00:00 2001
From: bpradines <barth371@hotmail.com>
Date: Fri, 26 Feb 2021 10:48:26 +0100
Subject: [PATCH 17/21] added dipole moments in determinants

---
 src/determinants/dipole_moments.irp.f | 65 +++++++++++++++++++++++++++
 1 file changed, 65 insertions(+)
 create mode 100644 src/determinants/dipole_moments.irp.f

diff --git a/src/determinants/dipole_moments.irp.f b/src/determinants/dipole_moments.irp.f
new file mode 100644
index 00000000..8a5f1a2d
--- /dev/null
+++ b/src/determinants/dipole_moments.irp.f
@@ -0,0 +1,65 @@
+ BEGIN_PROVIDER [double precision, z_dipole_moment, (N_states)]
+&BEGIN_PROVIDER [double precision, y_dipole_moment, (N_states)]
+&BEGIN_PROVIDER [double precision, x_dipole_moment, (N_states)]
+ implicit none
+ BEGIN_DOC
+ ! blablabla 
+ END_DOC
+ integer :: ipoint,istate,i,j 
+ double precision :: weight, r(3)  
+ double precision :: cpu0,cpu1,nuclei_part_z,nuclei_part_y,nuclei_part_x
+
+ call cpu_time(cpu0)
+ z_dipole_moment = 0.d0
+ y_dipole_moment = 0.d0
+ x_dipole_moment = 0.d0
+ do istate = 1, N_states
+  do i = 1, mo_num  
+   z_dipole_moment(istate) += -(one_e_dm_mo_alpha(i,i,istate)+one_e_dm_mo_beta(i,i,istate)) *  mo_dipole_z(i,i)
+   y_dipole_moment(istate) += -(one_e_dm_mo_alpha(i,i,istate)+one_e_dm_mo_beta(i,i,istate)) *  mo_dipole_y(i,i)
+   x_dipole_moment(istate) += -(one_e_dm_mo_alpha(i,i,istate)+one_e_dm_mo_beta(i,i,istate)) *  mo_dipole_x(i,i)
+   do j = i+1, mo_num  
+    z_dipole_moment(istate) += - 2.d0 * (one_e_dm_mo_alpha(j,i,istate)+one_e_dm_mo_beta(j,i,istate)) *  mo_dipole_z(j,i) 
+    y_dipole_moment(istate) += - 2.d0 * (one_e_dm_mo_alpha(j,i,istate)+one_e_dm_mo_beta(j,i,istate)) *  mo_dipole_y(j,i) 
+    x_dipole_moment(istate) += - 2.d0 * (one_e_dm_mo_alpha(j,i,istate)+one_e_dm_mo_beta(j,i,istate)) *  mo_dipole_x(j,i) 
+   enddo
+  enddo
+ enddo
+ 
+ print*,'electron part for z_dipole = ',z_dipole_moment
+ print*,'electron part for y_dipole = ',y_dipole_moment
+ print*,'electron part for x_dipole = ',x_dipole_moment
+ 
+ nuclei_part_z = 0.d0
+ nuclei_part_y = 0.d0
+ nuclei_part_x = 0.d0
+ do i = 1,nucl_num 
+  nuclei_part_z += nucl_charge(i) * nucl_coord(i,3) 
+  nuclei_part_y += nucl_charge(i) * nucl_coord(i,2) 
+  nuclei_part_x += nucl_charge(i) * nucl_coord(i,1) 
+ enddo
+ print*,'nuclei   part for z_dipole = ',nuclei_part_z
+ print*,'nuclei   part for y_dipole = ',nuclei_part_y
+ print*,'nuclei   part for x_dipole = ',nuclei_part_x
+
+ do istate = 1, N_states
+  z_dipole_moment(istate) += nuclei_part_z
+  y_dipole_moment(istate) += nuclei_part_y
+  x_dipole_moment(istate) += nuclei_part_x
+ enddo
+
+ call cpu_time(cpu1)
+ print*,'Time to provide the dipole moment :',cpu1-cpu0
+END_PROVIDER
+
+
+
+
+ subroutine print_z_dipole_moment_only
+ implicit none
+  print*, ''
+  print*, ''
+  print*,  '****************************************'
+  print*, 'z_dipole_moment                  = ',z_dipole_moment
+  print*,  '****************************************'
+ end

From c1669181b2e66b2a4f2c24090001ad06d53a2152 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Mon, 15 Mar 2021 17:40:37 +0100
Subject: [PATCH 18/21] minor changes

---
 src/mo_two_e_erf_ints/core_quantities_erf.irp.f |  2 +-
 src/tools/print_dipole.irp.f                    |  5 +++++
 src/tools/print_var_energy.irp.f                | 11 +++++++++++
 3 files changed, 17 insertions(+), 1 deletion(-)
 create mode 100644 src/tools/print_dipole.irp.f
 create mode 100644 src/tools/print_var_energy.irp.f

diff --git a/src/mo_two_e_erf_ints/core_quantities_erf.irp.f b/src/mo_two_e_erf_ints/core_quantities_erf.irp.f
index 3cd68205..43d77904 100644
--- a/src/mo_two_e_erf_ints/core_quantities_erf.irp.f
+++ b/src/mo_two_e_erf_ints/core_quantities_erf.irp.f
@@ -7,7 +7,7 @@ BEGIN_PROVIDER [double precision, core_energy_erf]
  core_energy_erf = 0.d0
  do i = 1, n_core_orb
   j = list_core(i)
-  core_energy_erf += 2.d0 * mo_one_e_integrals(j,j) + mo_two_e_int_erf_jj(j,j)
+  core_energy_erf +=  mo_two_e_int_erf_jj(j,j)
   do k = i+1, n_core_orb
    l = list_core(k)
    core_energy_erf += 2.d0 * (2.d0 * mo_two_e_int_erf_jj(j,l) - mo_two_e_int_erf_jj_exchange(j,l))
diff --git a/src/tools/print_dipole.irp.f b/src/tools/print_dipole.irp.f
new file mode 100644
index 00000000..8351308e
--- /dev/null
+++ b/src/tools/print_dipole.irp.f
@@ -0,0 +1,5 @@
+program print_dipole
+ implicit none
+ call print_z_dipole_moment_only
+
+end
diff --git a/src/tools/print_var_energy.irp.f b/src/tools/print_var_energy.irp.f
new file mode 100644
index 00000000..bab57d8c
--- /dev/null
+++ b/src/tools/print_var_energy.irp.f
@@ -0,0 +1,11 @@
+program print_var_energy
+ implicit none
+ read_wf = .True.
+ touch read_wf
+ call routine
+end
+
+subroutine routine
+ implicit none
+ print*,'psi_energy = ',psi_energy
+end

From 1935cc845e9a6ad11c8181e5afb45a821a2a76f4 Mon Sep 17 00:00:00 2001
From: bpradines <barth371@hotmail.com>
Date: Tue, 30 Mar 2021 14:58:21 +0200
Subject: [PATCH 19/21] removed nan in density_for_dft when no beta electrons

---
 src/density_for_dft/density_for_dft.irp.f | 20 +++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/src/density_for_dft/density_for_dft.irp.f b/src/density_for_dft/density_for_dft.irp.f
index ee70cd38..364cedd8 100644
--- a/src/density_for_dft/density_for_dft.irp.f
+++ b/src/density_for_dft/density_for_dft.irp.f
@@ -41,7 +41,11 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_alpha_for_dft, (mo_num,mo_num, N_s
      elec_alpha_valence(istate) += one_e_dm_mo_alpha_for_dft(i,i,istate)
     enddo
     elec_alpha_valence(istate) = elec_alpha_frozen_num/elec_alpha_valence(istate)
-    one_e_dm_mo_alpha_for_dft(:,:,istate) = one_e_dm_mo_alpha_for_dft(:,:,istate) * elec_alpha_valence(istate)
+    if( dabs(elec_alpha_valence(istate)) .lt.1.d-12)then
+     one_e_dm_mo_alpha_for_dft = 0.d0
+    else
+     one_e_dm_mo_alpha_for_dft(:,:,istate) = one_e_dm_mo_alpha_for_dft(:,:,istate) * elec_alpha_valence(istate)
+    endif
    enddo 
    
   endif
@@ -55,6 +59,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st
 ! density matrix for beta  electrons in the MO basis used for all DFT calculations based on the density
  END_DOC
  double precision :: delta_beta(mo_num,mo_num,N_states)
+ one_e_dm_mo_beta_for_dft = 0.d0
  if(density_for_dft .EQ. "damping_rs_dft")then
   delta_beta = one_e_dm_mo_beta - data_one_e_dm_beta_mo
   one_e_dm_mo_beta_for_dft = data_one_e_dm_beta_mo + damping_for_rs_dft * delta_beta
@@ -73,6 +78,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st
   one_e_dm_mo_beta_for_dft(:,:,1) = one_e_dm_mo_beta_average(:,:)
  endif
 
+
  if(no_core_density)then
   integer :: ii,i,j
   do ii = 1, n_core_orb
@@ -82,17 +88,21 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st
     one_e_dm_mo_beta_for_dft(i,j,:) = 0.d0
    enddo
   enddo
+  double precision :: elec_beta_valence(N_states),elec_beta_frozen_num
+  integer :: istate 
   if(normalize_dm)then
-   double precision :: elec_beta_valence(N_states),elec_beta_frozen_num
    elec_beta_frozen_num = elec_beta_num - n_core_orb 
    elec_beta_valence = 0.d0
-   integer :: istate 
    do istate = 1, N_states
     do i = 1, mo_num
      elec_beta_valence(istate) += one_e_dm_mo_beta_for_dft(i,i,istate)
     enddo
-    elec_beta_valence(istate) = elec_beta_frozen_num/elec_beta_valence(istate)
-    one_e_dm_mo_beta_for_dft(:,:,istate) = one_e_dm_mo_beta_for_dft(:,:,istate) * elec_beta_valence(istate)
+    if(dabs(elec_beta_valence(istate)).lt.1.d-12)then
+     one_e_dm_mo_beta_for_dft = 0.d0
+    else
+     elec_beta_valence(istate) = elec_beta_frozen_num/elec_beta_valence(istate)
+     one_e_dm_mo_beta_for_dft(:,:,istate) = one_e_dm_mo_beta_for_dft(:,:,istate) * elec_beta_valence(istate)
+    endif
    enddo 
   endif
  endif

From 4a6f7a3a928fa634fc1a5384eb43494c59f3d2bf Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Thu, 8 Apr 2021 20:37:17 +0200
Subject: [PATCH 20/21] added some stuffs for foboscf

---
 src/bitmask/modify_bitmasks.irp.f         |  8 ++--
 src/cas_based_on_top/example.irp.f        | 16 ++++---
 src/dft_utils_in_r/dm_in_r_routines.irp.f | 56 +++++++++++++++++++++++
 src/scf_utils/EZFIO.cfg                   |  7 +++
 src/scf_utils/diagonalize_fock.irp.f      | 21 +++++++++
 src/scf_utils/fock_matrix.irp.f           | 21 +++++++++
 6 files changed, 118 insertions(+), 11 deletions(-)

diff --git a/src/bitmask/modify_bitmasks.irp.f b/src/bitmask/modify_bitmasks.irp.f
index 834be6c8..ed291adf 100644
--- a/src/bitmask/modify_bitmasks.irp.f
+++ b/src/bitmask/modify_bitmasks.irp.f
@@ -143,10 +143,10 @@ subroutine print_generators_bitmasks_holes
      key_tmp(j,1) = generators_bitmask(j,1,i)
      key_tmp(j,2) = generators_bitmask(j,2,i)
    enddo
-   print*,''
-   print*,'index hole  = ',i
-   call print_det(key_tmp,N_int)
-   print*,''
+!   print*,''
+!   print*,'index hole  = ',i
+!   call print_det(key_tmp,N_int)
+!   print*,''
  enddo
  deallocate(key_tmp)
 
diff --git a/src/cas_based_on_top/example.irp.f b/src/cas_based_on_top/example.irp.f
index f00956e3..2f709495 100644
--- a/src/cas_based_on_top/example.irp.f
+++ b/src/cas_based_on_top/example.irp.f
@@ -5,37 +5,39 @@ subroutine write_on_top_in_real_space
 ! This routines is a simple example of how to plot the on-top pair density on a simple 1D grid 
  END_DOC
  double precision :: zmax,dz,r(3),on_top_in_r,total_density,zcenter,dist
+ double precision :: core_dens, inact_dens,act_dens(2,1)
 
  integer :: nz,i,istate
  character*(128)                :: output
  integer                        :: i_unit_output,getUnitAndOpen
  PROVIDE ezfio_filename
  output=trim(ezfio_filename)//'.on_top'
- print*,'output = ',trim(output)                                                                                         
+ print*,'output = ',trim(output)
  i_unit_output = getUnitAndOpen(output,'w')
 
 
- zmax = 2.0d0
+ zmax = 5.0d0
  print*,'nucl_coord(1,3) = ',nucl_coord(1,3)
  print*,'nucl_coord(2,3) = ',nucl_coord(2,3)
  dist = dabs(nucl_coord(1,3) - nucl_coord(2,3))
- zmax += dist 
+ zmax += dist
  zcenter = (nucl_coord(1,3) + nucl_coord(2,3))*0.5d0
  print*,'zcenter = ',zcenter
  print*,'zmax    = ',zmax
  nz = 1000
  dz = zmax / dble(nz)
- r(:) = 0.d0 
- r(3) = zcenter -zmax * 0.5d0 
+ r(:) = 0.d0
+ r(3) = zcenter -zmax * 0.5d0                                                                                                                                          
  print*,'r(3)    = ',r(3)
  istate = 1
+
+  write(i_unit_output,*)" z,    on-top(z),    n(z) "
  do i = 1, nz
   call give_on_top_in_r_one_state(r,istate,on_top_in_r)
-  call give_cas_density_in_r(r,total_density)
+  call give_cas_density_in_r(core_dens,inact_dens,act_dens,total_density,r)
   write(i_unit_output,*)r(3),on_top_in_r,total_density
   r(3) += dz
  enddo
 
-
 end
 
diff --git a/src/dft_utils_in_r/dm_in_r_routines.irp.f b/src/dft_utils_in_r/dm_in_r_routines.irp.f
index 6fa99e22..9991289c 100644
--- a/src/dft_utils_in_r/dm_in_r_routines.irp.f
+++ b/src/dft_utils_in_r/dm_in_r_routines.irp.f
@@ -110,6 +110,62 @@ end
    grad_dm_b *= 2.d0
  end
 
+ subroutine density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)
+ implicit none
+ BEGIN_DOC
+! input:
+!
+! * r(1) ==> r(1) = x, r(2) = y, r(3) = z
+!
+! output:
+!
+! * dm_a = alpha density evaluated at r
+! * dm_b = beta  density evaluated at r
+! * grad_dm_a(1) = X gradient of the alpha density evaluated in r
+! * grad_dm_a(1) = X gradient of the beta  density evaluated in r
+!
+ END_DOC
+ double precision, intent(in)  :: r(3)
+ double precision, intent(out) :: dm_a(N_states),dm_b(N_states)
+ double precision, intent(out) :: grad_dm_a(3,N_states),grad_dm_b(3,N_states)
+ double precision              :: grad_aos_array(3,ao_num)
+ integer :: i,j,istate
+ double precision  :: aos_array(ao_num),aos_array_bis(ao_num),u_dot_v
+ double precision  :: aos_grad_array(ao_num,3), aos_grad_array_bis(ao_num,3)
+
+ call give_all_aos_and_grad_at_r(r,aos_array,grad_aos_array)
+ do i = 1, ao_num
+  do j = 1, 3
+   aos_grad_array(i,j) =  grad_aos_array(j,i)
+  enddo
+ enddo
+
+ do istate = 1, N_states
+   ! alpha density
+   ! aos_array_bis = \rho_ao * aos_array
+   call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
+   dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
+
+   ! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
+   grad_dm_a(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
+   grad_dm_a(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
+   grad_dm_a(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
+   ! aos_grad_array_bis = \rho_ao * aos_grad_array
+
+   ! beta density
+   call dsymv('U',ao_num,1.d0,one_e_dm_beta_ao_for_dft(1,1,istate),size(one_e_dm_beta_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
+   dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
+
+   ! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
+   grad_dm_b(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
+   grad_dm_b(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
+   grad_dm_b(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
+   ! aos_grad_array_bis = \rho_ao * aos_grad_array
+ enddo
+   grad_dm_a *= 2.d0
+   grad_dm_b *= 2.d0
+ end
+
 
 
  subroutine density_and_grad_lapl_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, grad_dm_a, grad_dm_b, lapl_dm_a, lapl_dm_b, aos_array, grad_aos_array, lapl_aos_array)
diff --git a/src/scf_utils/EZFIO.cfg b/src/scf_utils/EZFIO.cfg
index 4a56a35b..694590ec 100644
--- a/src/scf_utils/EZFIO.cfg
+++ b/src/scf_utils/EZFIO.cfg
@@ -51,3 +51,10 @@ doc: If true, leave untouched all the orbitals defined as core and optimize all
 interface: ezfio,provider,ocaml
 default: False
 
+
+[no_oa_or_av_opt]
+type: logical
+doc: If true, you set to zero all Fock elements between the orbital set to active and all the other orbitals 
+interface: ezfio,provider,ocaml
+default: False
+
diff --git a/src/scf_utils/diagonalize_fock.irp.f b/src/scf_utils/diagonalize_fock.irp.f
index d501278f..5188581a 100644
--- a/src/scf_utils/diagonalize_fock.irp.f
+++ b/src/scf_utils/diagonalize_fock.irp.f
@@ -31,6 +31,27 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
       enddo
      enddo
    endif
+   if(no_oa_or_av_opt)then
+     do i = 1, n_act_orb
+       iorb = list_act(i)
+       do j = 1, n_inact_orb
+         jorb = list_inact(j)
+         F(iorb,jorb) = 0.d0
+         F(jorb,iorb) = 0.d0
+       enddo
+       do j = 1, n_virt_orb
+         jorb = list_virt(j)
+         F(iorb,jorb) = 0.d0
+         F(jorb,iorb) = 0.d0
+       enddo
+       do j = 1, n_core_orb
+         jorb = list_core(j)
+         F(iorb,jorb) = 0.d0
+         F(jorb,iorb) = 0.d0                                                                                                                 
+       enddo
+     enddo
+   endif
+
 
    ! Insert level shift here
    do i = elec_beta_num+1, elec_alpha_num
diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f
index fc9eaadd..61633d3b 100644
--- a/src/scf_utils/fock_matrix.irp.f
+++ b/src/scf_utils/fock_matrix.irp.f
@@ -92,6 +92,27 @@
      enddo
    endif
 
+   if(no_oa_or_av_opt)then
+     do i = 1, n_act_orb
+       iorb = list_act(i)
+       do j = 1, n_inact_orb
+         jorb = list_inact(j)
+         Fock_matrix_mo(iorb,jorb) = 0.d0
+         Fock_matrix_mo(jorb,iorb) = 0.d0
+       enddo
+       do j = 1, n_virt_orb
+         jorb = list_virt(j)
+         Fock_matrix_mo(iorb,jorb) = 0.d0
+         Fock_matrix_mo(jorb,iorb) = 0.d0
+       enddo
+       do j = 1, n_core_orb
+         jorb = list_core(j)
+         Fock_matrix_mo(iorb,jorb) = 0.d0
+         Fock_matrix_mo(jorb,iorb) = 0.d0                                                                                                                 
+       enddo
+     enddo
+   endif
+
 END_PROVIDER
 
 

From 6ecb8d1be6c5a7a7183e5306d7363e79466c5d8c Mon Sep 17 00:00:00 2001
From: Gilles Frison <gfrison@master.lct.jussieu.fr>
Date: Thu, 15 Apr 2021 15:46:41 +0200
Subject: [PATCH 21/21] added opam_installer_no_usr_bin.sh

---
 scripts/opam_installer_no_usr_bin.sh | 365 +++++++++++++++++++++++++++
 1 file changed, 365 insertions(+)
 create mode 100755 scripts/opam_installer_no_usr_bin.sh

diff --git a/scripts/opam_installer_no_usr_bin.sh b/scripts/opam_installer_no_usr_bin.sh
new file mode 100755
index 00000000..0f7219d1
--- /dev/null
+++ b/scripts/opam_installer_no_usr_bin.sh
@@ -0,0 +1,365 @@
+#!/bin/sh
+
+set -ue
+
+# (c) Copyright Fabrice Le Fessant INRIA/OCamlPro 2013
+# (c) Copyright Louis Gesbert OCamlPro 2014-2017
+
+VERSION='2.0.7'
+DEV_VERSION='2.1.0~beta2'
+DEFAULT_BINDIR=/usr/local/bin
+
+bin_sha512() {
+  case "$OPAM_BIN" in
+    opam-2.0.6-arm64-linux)     echo "d2b3d92fd5fae7f053702b53ddbc7c224fcfbfc9b232247ba4e40cbf1cda28f160d8c14fde87aebeebfd2545e13265c0ee4a47e292f035767fb944b1b8ff5c90";;
+    opam-2.0.6-armhf-linux)     echo "a42a7ad8c1afdb20ac5746934306576e6364f5453b176ccd42a3e5a116a5db05c2758cec31800ffab11411290cf671f9eee3f299df48c7ceca8e4d7e33dfedc8";;
+    opam-2.0.6-i686-linux)      echo "6c0d965f89a2026ead3120e217d12b2df7426740d54bc94e2c46faaeff5893081e68aac162621bfa694ab597a18be28165f10cdda1217a4d73653789a9928b64";;
+    opam-2.0.6-x86_64-linux)    echo "2b9d4a99aa28a193c88c7c6f6265203bd3cfeef98929d6f5cfce4b52cd9ddbd7be7eddc1d3d9c440f81d65074dd7851b8d29cd397fb06d2cfccffb54d3cdcc6a";;
+    opam-2.0.6-x86_64-macos)    echo "cf02546b22ca91b1d97a3657b970b34d4acf4dc745696b7200ff185d25ebb5914ea8b6a94b503eb8c999634de6fdb944998a970105cd6a4c6df538c262b48b7f";;
+    opam-2.0.6-x86_64-openbsd)  echo "2f58b3d4902d4c3fb823d251a50e034f9101b0c5a3827725876bb3bcb6c013c4f54138054d82abba0a9e917675275e26f05b98630cf7116c465d2110756f1309";;
+
+    opam-2.0.7-arm64-linux)     echo "0dd4d80496545f684af39dc5b4b28867bc19a74186577c38bd2a8934d871c2cbcdb9891bfd41c080b5f12d5a3c8801e203df8a76d55e1e22fe80d31447402e46";;
+    opam-2.0.7-armhf-linux)     echo "ea691bc9565acc1207dea3dfb89192b1865b5b5809efe804a329f39878640fb19771edcb05c5699f8e914e88e3155f31132b845c54b0095bedd3952d336bae0b";;
+    opam-2.0.7-i686-linux)      echo "5fa8fb9664d36ead5760e7e1c337f6ae7b0fd4be5089ddfb50ae74028deec30893b1f4dee040402bc3f15da197ba89a45c7d626ecf6e5be80d176f43526c4bad";;
+    opam-2.0.7-x86_64-linux)    echo "da75b0cb5ad50f95d31857a7d72f1836132a1fa1cdbfdedf684342b798e7107b4add4c74c05d5ce44881309fa1e57707538dbcda874e7f74b269b1bb204f3ae3";;
+    opam-2.0.7-x86_64-macos)    echo "de1194c8e97e53956e5e47502c28881bbf26d1beaac4f33a43a922b8ca7ce97725533cfaf65a33fc0e183eab5a95e9ecd2e20f72faeaec333dc3850b79b5fe8a";;
+    opam-2.0.7-x86_64-openbsd)  echo "b253809c4388847e1a33b5c4f1f5d72bef79a2f0c43b19ef65b40d0c10341aa0bee4a4b1f3a9ab70eb026e4cc220a63cfc56a18c035b6b0297c92f2bdb7f9a78";;
+
+    opam-2.1.0-alpha-arm64-linux)     echo "1bf0acfa64aa01c3244e65eed60eef1caaa6de53aa8b32dd0d2446f91905a1e41591f53cd350e85b2b9f5edba9b137d723c32949115623e9753e77b707bb25b0";;
+    opam-2.1.0-alpha-armhf-linux)     echo "87c12a422bd14a0d10a94ddaaa46de23700e3b89810a0c06232eff8d96b37c2fd43dcb5a8da5a2004aa8040d1b93293209f1ff1aab865ffd150364e24c87c716";;
+    opam-2.1.0-alpha-i686-linux)      echo "b8369da6d4795a461ff1b49e687b027325d4e90bc8f19517e52a94ee3be167c4faaaf33bd0b3536be552d2add54865d0e33933acaa674f2e1a17249b022738af";;
+    opam-2.1.0-alpha-x86_64-linux)    echo "2e22747829fb0bada3a74a23f5e0ff2228520d647fc4fe08a1ce76f3cb357cc7240f7b45e422c5f4b8eafe832ae3a8973ecbd4814ae0e8ce1096bcff39482020";;
+    opam-2.1.0-alpha-x86_64-macos)    echo "c440e8ae1970fa7533e6e1b96ba3e3dd65b04432d41bc57ce4c768ed9b4229954546d59ec06f3d4ee49cbe00bb4bfd0b3f509d6d9a27de2db17725e097a61c86";;
+    opam-2.1.0-alpha-x86_64-openbsd)  echo "d87afe99fee541a1c6fae30b72653db7a5ea2abdec3fa3b2b480daddf3fcd8d4096e2a40458310755faec3722119f29ed981ffbfa65142e618f99b70572f892f";;
+
+    opam-2.1.0-alpha2-arm64-linux)    echo "b67520bb2a6c59f800da100278d74e58f2bbf66924f94643023dc46b97b16f17a30de95e439c6f9b032bd555c062ddba325f3e5169cac186615b959a8c434788";;
+    opam-2.1.0-alpha2-armhf-linux)    echo "9a6312eb54d6c9c2036ca90f7816789c27c23f1b1d325cd69d27a910cdd8760b82f19c9e9b61b5b6214818f1f40f8b4d2ef081acb43f0dad68c976986a7c6a45";;
+    opam-2.1.0-alpha2-i686-linux)     echo "0dc07f236405777ad74d58fcc6cb6c3247e7dfc31408df4a199599077d5cb41ec86895f1d0c5eaa2a9c70842a2a998226674f986ba0044c82896c073ac90b209";;
+    opam-2.1.0-alpha2-x86_64-linux)   echo "21509e8abd8463f4e18a55398f690700772e25f0ddb9f3fd7644e2f9a9a89ebbf5c09efbeceafe4a0ab5015d0d03b2f29506be514aae813a2f3dac7dd01261f3";;
+    opam-2.1.0-alpha2-x86_64-macos)   echo "1c1bd26621eebb5bf3783dec80d5555aa5ff02dcbf43eb44398798e6162c1964bc1964e3980391ea115e5c068c1bb66960f8ebdd91bc4f0bac844f3a61433f1e";;
+    opam-2.1.0-alpha2-x86_64-openbsd) echo "941f3e306bc36e8e44e4245ca5e635b04e0a54f33439d55d41875ced47384cad8c222b649027d3c4eacc3c2c569cf5006c872763b19c490d9b289c9cfe4f491a";;
+
+    opam-2.1.0-alpha3-arm64-linux)     echo "ad906bb2ab764a92fabdf0b906310c5034bf5daf0ebfb2529e9b87661ddbf8fd14f51dee5ce75b4fd4bb5789e29c7be71063f1ebcc92e92333be12aa62efdff9";;
+    opam-2.1.0-alpha3-armhf-linux)     echo "2a7022c1f5dbc855a0d067f29677b13253dccbc9792b8170fa72a743802bbcd6e41ce7512c4845091af0f73b8ba7573038ec53ea9aaf74be04367ac1767e7220";;
+    opam-2.1.0-alpha3-i686-linux)      echo "6f2fce0c45ae700e7a1b32d0a24988645c9aed3afc45998c8fbe70e97a65e3ba5d824069914a892bb3f9b1336383cfd492c28678ff16db5cada863da924b07d8";;
+    opam-2.1.0-alpha3-x86_64-linux)    echo "1d219dbf670e1550bf71c28e586d14f1d8af2605f0e13bea2f11ad52a7f176bd9a89637e44a91a024f0088db1b2aba8dc3207bc81fa930580e54f4031255c178";;
+    opam-2.1.0-alpha3-x86_64-macos)    echo "93edb6c1151f8f5bd017f230ffd9277f6ad943e3f5032ea000c37f012738fb3ab4b4add172e1f624c37e6564963fef0716b876b0113c8e43f5943d77bbbc173c";;
+    opam-2.1.0-alpha3-x86_64-openbsd)  echo "0e3b3761e877c57f5b333aacb70c86bf60f50eecdca6e9e1a552e3d666cea034d8873f3a87e585a5970b1aef7e540adb18c71e0e8fd8794843dd5d1d421a87ec";;
+
+    opam-2.1.0-beta-arm64-linux)     echo "954670c74ea8244b440756e4f7755bd2b5548ab67428ce577c4c507fc33c8d00eb73c4d7b59ccb0ef800f4465b5c704573c63486b78a23e9568f3751bf9aef78";;
+    opam-2.1.0-beta-armhf-linux)     echo "cc666f2c6b1ac07d1bc8a035c6b3a9455794b51a827c54bb92786ae1a75c6c55839d3f48b378508f42a66ac887fdc68f7628a67e2826813cb6df048c906755ca";;
+    opam-2.1.0-beta-i686-linux)      echo "66ac48b298741f753ca868be362851ccd9bf84fd8772d18f3307e99cf72c8c68ac9fa17bf2d610d7f3b5dc6209eb8371bf0e10b363e963fc6c31d70e5938017f";;
+    opam-2.1.0-beta-x86_64-linux)    echo "e316f1b5f1c668affba6c2819f692c28776e131a17fb64b2c0e23f8a3b7d456575a8109fcdcb9babfad13bc33c17fa619cbb4a48ca6198765f86296b7e611f24";;
+    opam-2.1.0-beta-x86_64-macos)    echo "acb29b7c64df314c6629e14f6d8f079504d39b7fd3104867fd22df3395ccfea9f1014a3a87dff9c12bf03ca451e9ee2918b9d9d8f17ce1a6d7de0c0649452fa9";;
+    opam-2.1.0-beta-x86_64-openbsd)  echo "ff9fa1ee0ae7e54b4e18999cf5ea9b899c0b4039b744a950e96221e3e86c21eaa50904bdbc836ff8103f7713506d0de3d32ec77b169561e0cd694bfeea812cae";;
+
+    opam-2.1.0-beta2-arm64-linux)     echo "a58ba3ebb4431d3cabfe96b806c9897205153e8a546ebe74f0229982758d140b4fcbcea421db70589b1eb3080dc86534522a3cba0330ce82e0898a60048d51ba";;
+    opam-2.1.0-beta2-armhf-linux)     echo "fc4e6b753ce6368f75a0d3005f4b21ce9606599d21607a67015db55a38b6ef473b4205f5b128c5808189feed8ae58f93bd79348988be7c5007ae1b39307a5cd0";;
+    opam-2.1.0-beta2-i686-linux)      echo "a376a6e0e1e2b08ea4d0a5c1c38487e67984bef2e89f978536dd08283f945f74dd31ee287bc68d91690603ba0fa657e91ff0d30bea217743f79ed99d2390eba5";;
+    opam-2.1.0-beta2-x86_64-linux)    echo "12c5e2b0087ed389fa12fdb0e1f6f7dc0b3df3f95c59e8bc576279b7780921d47bbc4ebcba6caddde30f4fb1cc9e4a873cc8a6aef80fcc48a878aba69be7af44";;
+    opam-2.1.0-beta2-x86_64-macos)    echo "4acc12672a2e3ad7e78540634edcae2e7e84860057b86a56b1cdf7eacf8d97957aaa864f571d6fb8f61ee8280f8a4ed73df7881d91a22c9d8c2d73e8a558f61d";;
+    opam-2.1.0-beta2-x86_64-openbsd)  echo "84d7d409220c72e3ed7e6acdd7cce3b5a208f2966d232648a57a48641ab8ce4fa58e94e40b7176201455d82260e6c501a6ba4a30b1426a552f8d09cfd027ddde";;
+
+    *) echo "no sha";;
+  esac
+}
+
+usage() {
+    echo "opam binary installer v.$VERSION"
+    echo "Downloads and installs a pre-compiled binary of opam $VERSION to the system."
+    echo "This can also be used to switch between opam versions"
+    echo
+    echo "Options:"
+    echo "    --dev                  Install the latest alpha or beta instead: $DEV_VERSION"
+    echo "    --no-backup            Don't attempt to backup the current opam root"
+    echo "    --backup               Force the backup the current opam root (even if it"
+    echo "                           is from the 2.0 branch already)"
+    echo "    --fresh                Create the opam $VERSION root from scratch"
+    echo "    --restore   VERSION    Restore a backed up opam binary and root"
+    echo
+    echo "The default is to backup if the current version of opam is 1.*, or when"
+    echo "using '--fresh' or '--dev'"
+}
+
+RESTORE=
+NOBACKUP=
+FRESH=
+DOWNLOAD_ONLY=
+
+while [ $# -gt 0 ]; do
+    case "$1" in
+        --dev)
+            VERSION=$DEV_VERSION
+            if [ -z "$NOBACKUP" ]; then NOBACKUP=0; fi;;
+        --restore)
+            if [ $# -lt 2 ]; then echo "Option $1 requires an argument"; exit 2; fi
+            shift;
+            RESTORE=$1;;
+        --no-backup)
+            NOBACKUP=1;;
+        --backup)
+            NOBACKUP=0;;
+        --fresh)
+            FRESH=1;;
+        --download-only)
+            DOWNLOAD_ONLY=1;;
+        --help|-h)
+            usage; exit 0;;
+        *)
+            usage; exit 2;;
+    esac
+    shift
+done
+
+
+TMP=${TMPDIR:-/tmp}
+
+ARCH=$(uname -m || echo unknown)
+case "$ARCH" in
+    x86|i?86) ARCH="i686";;
+    x86_64|amd64) ARCH="x86_64";;
+    ppc|powerpc|ppcle) ARCH="ppc";;
+    aarch64_be|aarch64|armv8b|armv8l) ARCH="arm64";;
+    armv5*|armv6*|earmv6*|armv7*|earmv7*) ARCH="armhf";;
+    *) ARCH=$(echo "$ARCH" | awk '{print tolower($0)}')
+esac
+
+OS=$( (uname -s || echo unknown) | awk '{print tolower($0)}')
+
+if [ "$OS" = "darwin" ] ; then
+  OS=macos
+fi
+
+TAG=$(echo "$VERSION" | tr '~' '-')
+
+OPAM_BIN_URL_BASE='https://github.com/ocaml/opam/releases/download/'
+OPAM_BIN="opam-${TAG}-${ARCH}-${OS}"
+OPAM_BIN_URL="${OPAM_BIN_URL_BASE}${TAG}/${OPAM_BIN}"
+
+download() {
+    if command -v wget >/dev/null; then wget -q -O "$@"
+    else curl -s -L -o "$@"
+    fi
+}
+
+check_sha512() {
+    OPAM_BIN_LOC="$1"
+    if command -v openssl > /dev/null; then
+        sha512_devnull="cf83e1357eefb8bdf1542850d66d8007d620e4050b5715dc83f4a921d36ce9ce47d0d13c5d85f2b0ff8318d2877eec2f63b931bd47417a81a538327af927da3e"
+        sha512_check=`openssl sha512 2>&1 < /dev/null | cut -f 2 -d ' '`
+        if [ "x$sha512_devnull" = "x$sha512_check" ]; then
+            sha512=`openssl sha512 "$OPAM_BIN_LOC" 2> /dev/null | cut -f 2 -d ' '`
+            check=`bin_sha512`
+            test "x$sha512" = "x$check"
+        else
+            echo "openssl 512 option not handled, binary integrity check can't be performed."
+            return 0
+        fi
+    else
+        echo "openssl not found, binary integrity check can't be performed."
+        return 0
+    fi
+}
+
+download_and_check() {
+    OPAM_BIN_LOC="$1"
+    echo "## Downloading opam $VERSION for $OS on $ARCH..."
+
+    if ! download "$OPAM_BIN_LOC" "$OPAM_BIN_URL"; then
+        echo "There may not yet be a binary release for your architecture or OS, sorry."
+        echo "See https://github.com/ocaml/opam/releases/tag/$TAG for pre-compiled binaries,"
+        echo "or run 'make cold' from https://github.com/ocaml/opam/archive/$TAG.tar.gz"
+        echo "to build from scratch"
+        exit 10
+    else
+        if check_sha512 "$OPAM_BIN_LOC"; then
+            echo "## Downloaded."
+        else
+            echo "Checksum mismatch, a problem occurred during download."
+            exit 10
+        fi
+    fi
+}
+
+DOWNLOAD_ONLY=${DOWNLOAD_ONLY:-0}
+
+if [ $DOWNLOAD_ONLY -eq 1 ]; then
+    OPAM_BIN_LOC="$PWD/$OPAM_BIN"
+    if [ -e "$OPAM_BIN_LOC" ]; then
+        echo "Found opam binary in $OPAM_BIN_LOC ..."
+        if check_sha512 "$OPAM_BIN_LOC" ; then
+            echo "... with matching sha512"
+            exit 0;
+        else
+            echo "... with mismatching sha512, download the good one."
+        fi
+    fi
+    download_and_check "$OPAM_BIN_LOC"
+    exit 0;
+fi
+
+EXISTING_OPAM=$(command -v opam || echo)
+EXISTING_OPAMV=
+if [ -n "$EXISTING_OPAM" ]; then
+   EXISTING_OPAMV=$("$EXISTING_OPAM" --version || echo "unknown")
+fi
+
+FRESH=${FRESH:-0}
+
+OPAMROOT=${OPAMROOT:-$HOME/.opam}
+
+if [ ! -d "$OPAMROOT" ]; then FRESH=1; fi
+
+if [ -z "$NOBACKUP" ] && [ ! "$FRESH" = 1 ] && [ -z "$RESTORE" ]; then
+    case "$EXISTING_OPAMV" in
+        2.*) NOBACKUP=1;;
+        *) NOBACKUP=0;;
+    esac
+fi
+
+xsudo() {
+    local CMD=$1; shift
+    local DST
+    for DST in "$@"; do : ; done
+
+    local DSTDIR=$(dirname "$DST")
+    if [ ! -w "$DSTDIR" ]; then
+        echo "Write access to $DSTDIR required, using 'sudo'."
+        echo "Command: $CMD $@"
+        if [ "$CMD" = "install" ]; then
+            sudo "$CMD" -g 0 -o root "$@"
+        else
+            sudo "$CMD" "$@"
+        fi
+    else
+        "$CMD" "$@"
+    fi
+}
+
+if [ -n "$RESTORE" ]; then
+    OPAM=$(command -v opam)
+    OPAMV=$("$OPAM" --version)
+    OPAM_BAK="$OPAM.$RESTORE"
+    OPAMROOT_BAK="$OPAMROOT.$RESTORE"
+    if [ ! -e "$OPAM_BAK" ] || [ ! -d "$OPAMROOT_BAK" ]; then
+        echo "No backup of opam $RESTORE was found"
+        exit 1
+    fi
+    if [ "$NOBACKUP" = 1 ]; then
+        printf "## This will clear $OPAM and $OPAMROOT. Continue ? [Y/n] "
+        read R
+        case "$R" in
+            ""|"y"|"Y"|"yes")
+                xsudo rm -f "$OPAM"
+                rm -rf "$OPAMROOT";;
+            *) exit 1
+        esac
+    else
+        xsudo mv "$OPAM" "$OPAM.$OPAMV"
+        mv "$OPAMROOT" "$OPAMROOT.$OPAMV"
+    fi
+    xsudo mv "$OPAM_BAK" "$OPAM"
+    mv "$OPAMROOT_BAK" "$OPAMROOT"
+    printf "## Opam $RESTORE and its root were restored."
+    if [ "$NOBACKUP" = 1 ]; then echo
+    else echo " Opam $OPAMV was backed up."
+    fi
+    exit 0
+fi
+
+#if [ -e "$TMP/$OPAM_BIN" ] && ! check_sha512 "$TMP/$OPAM_BIN" || [ ! -e "$TMP/$OPAM_BIN" ]; then
+    download_and_check "$TMP/$OPAM_BIN"
+#else
+#    echo "## Using already downloaded \"$TMP/$OPAM_BIN\""
+#fi
+
+if [ -n "$EXISTING_OPAM" ]; then
+    DEFAULT_BINDIR=$(dirname "$EXISTING_OPAM")
+fi
+
+while true; do
+    printf "## Where should it be installed ? [$DEFAULT_BINDIR] "
+    read BINDIR
+    if [ -z "$BINDIR" ]; then BINDIR="$DEFAULT_BINDIR"; fi
+
+    if [ -d "$BINDIR" ]; then break
+    else
+        printf "## $BINDIR does not exist. Create ? [Y/n] "
+        read R
+        case "$R" in
+            ""|"y"|"Y"|"yes")
+            xsudo mkdir -p $BINDIR
+            break;;
+        esac
+    fi
+done
+
+#if [ -e "$EXISTING_OPAM" ]; then
+#    if [ "$NOBACKUP" = 1 ]; then
+#        xsudo rm -f "$EXISTING_OPAM"
+#    else
+#        xsudo mv "$EXISTING_OPAM" "$EXISTING_OPAM.$EXISTING_OPAMV"
+#        echo "## $EXISTING_OPAM backed up as $(basename $EXISTING_OPAM).$EXISTING_OPAMV"
+#    fi
+#fi
+
+if [ -d "$OPAMROOT" ]; then
+    if [ "$FRESH" = 1 ]; then
+        if [ "$NOBACKUP" = 1 ]; then
+            printf "## This will clear $OPAMROOT. Continue ? [Y/n] "
+            read R
+            case "$R" in
+                ""|"y"|"Y"|"yes")
+                    rm -rf "$OPAMROOT";;
+                *) exit 1
+            esac
+        else
+            mv "$OPAMROOT" "$OPAMROOT.$EXISTING_OPAMV"
+            echo "## $OPAMROOT backed up as $(basename $OPAMROOT).$EXISTING_OPAMV"
+        fi
+        echo "## opam $VERSION installed. Please run 'opam init' to get started"
+    elif [ ! "$NOBACKUP" = 1 ]; then
+        echo "## Backing up $OPAMROOT to $(basename $OPAMROOT).$EXISTING_OPAMV (this may take a while)"
+        if [ -e "$OPAMROOT.$EXISTING_OPAMV" ]; then
+            echo "ERROR: there is already a backup at $OPAMROOT.$EXISTING_OPAMV"
+            echo "Please move it away or run with --no-backup"
+        fi
+        FREE=$(df -k "$OPAMROOT" | awk 'NR>1 {print $4}')
+        NEEDED=$(du -sk "$OPAMROOT" | awk '{print $1}')
+        if ! [ $NEEDED -lt $FREE ]; then
+            echo "Error: not enough free space to backup. You can retry with --no-backup,"
+            echo "--fresh, or remove '$OPAMROOT'"
+            exit 1
+        fi
+        cp -a "$OPAMROOT" "$OPAMROOT.$EXISTING_OPAMV"
+        echo "## $OPAMROOT backed up as $(basename $OPAMROOT).$EXISTING_OPAMV"
+    fi
+    rm -f "$OPAMROOT"/repo/*/*.tar.gz*
+fi
+
+xsudo install -m 755 "$TMP/$OPAM_BIN" "$BINDIR/opam"
+echo "## opam $VERSION installed to $BINDIR"
+
+if [ ! "$FRESH" = 1 ]; then
+    echo "## Converting the opam root format & updating"
+    "$BINDIR/opam" init --reinit -ni
+fi
+
+WHICH=$(command -v opam || echo notfound)
+
+case "$WHICH" in
+    "$BINDIR/opam") ;;
+    notfound) echo "## Remember to add $BINDIR to your PATH";;
+    *)
+        echo "## WARNING: 'opam' command found in PATH does not match the installed one:"
+        echo "   - Installed: '$BINDIR/opam'"
+        echo "   - Found:     '$WHICH'"
+        echo "Make sure to remove the second or fix your PATH to use the new opam"
+        echo
+esac
+
+if [ ! "$NOBACKUP" = 1 ]; then
+    echo "## Run this script again with '--restore $EXISTING_OPAMV' to revert."
+fi
+
+rm -f $TMP/$OPAM_BIN