Merge pull request #41 from QuantumPackage/dev-stable

Dev stable

docs/source/users_guide/qpsh.rst

@@ -8,7 +8,7 @@ qpsh
 
 
 :command:`qpsh` is the |qp| shell. It is a Bash shell with all the
-required evironment variables loaded, a modified prompt, and the
+required environment variables loaded, a modified prompt, and the
 :ref:`qp` command.
 
 

ocaml/Input_ao_basis.ml

@@ -58,26 +58,36 @@ end = struct
   ;;
 
   let read_ao_prim_num () =
+    if Ezfio.has_ao_basis_ao_prim_num () then
       Ezfio.get_ao_basis_ao_prim_num ()
       |> Ezfio.flattened_ezfio
       |> Array.map AO_prim_number.of_int
+    else
+      [||]
   ;;
 
   let read_ao_prim_num_max () =
+    if Ezfio.has_ao_basis_ao_prim_num () then
       Ezfio.get_ao_basis_ao_prim_num ()
       |> Ezfio.flattened_ezfio
       |> Array.fold_left (fun x y -> if x>y then x else y) 0
       |> AO_prim_number.of_int
+    else
+      AO_prim_number.of_int 0
   ;;
 
   let read_ao_nucl () =
+    if Ezfio.has_ao_basis_ao_nucl () then
       let nmax = Nucl_number.get_max () in
       Ezfio.get_ao_basis_ao_nucl ()
       |> Ezfio.flattened_ezfio
       |> Array.map (fun x-> Nucl_number.of_int ~max:nmax x)
+    else
+      [||]
   ;;
 
   let read_ao_power () =
+    if Ezfio.has_ao_basis_ao_power () then
       let x = Ezfio.get_ao_basis_ao_power () in
       let dim = x.Ezfio.dim.(0) in
       let data = Ezfio.flattened_ezfio x in
@@ -92,18 +102,26 @@ end = struct
           result.(i-1) <- result.(i-1)^"z"^(string_of_int data.(2*dim+i-1));
       done;
       Array.map Angmom.Xyz.of_string result
+    else
+      [||]
   ;;
 
   let read_ao_coef () =
+    if Ezfio.has_ao_basis_ao_coef () then
       Ezfio.get_ao_basis_ao_coef ()
       |> Ezfio.flattened_ezfio
       |> Array.map AO_coef.of_float
+    else
+      [||]
   ;;
 
   let read_ao_expo () =
+    if Ezfio.has_ao_basis_ao_expo () then
       Ezfio.get_ao_basis_ao_expo ()
       |> Ezfio.flattened_ezfio
       |> Array.map AO_expo.of_float
+    else
+      [||]
   ;;
 
   let read_ao_cartesian () =

ocaml/qp_run.ml

@@ -132,6 +132,7 @@ let run slave ?prefix exe ezfio_file =
   (** Run executable *)
   let prefix =
     match prefix with
+    | Some "gdb" -> "gdb --args "
     | Some x -> x^" "
     | None -> ""
   and exe =

plugins/local/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f

@@ -155,6 +155,7 @@ subroutine run_stochastic_cipsi
   call pt2_alloc(pt2_data_err, N_states)
   call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection
   call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm)
+  call print_summary_tc(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2)
   call pt2_dealloc(pt2_data)
   call pt2_dealloc(pt2_data_err)
 

plugins/local/jastrow/EZFIO.cfg

@@ -1,7 +1,25 @@
+[log_jpsi]
+type: logical
+doc: If |true|, the Jpsi is taken as log(1+psi_cor)
+interface: ezfio,provider,ocaml
+default: False
+
+
+[mu_of_r_tc]
+type: character*(32)
+doc: type of the mu(r): [ Standard | Erfmu | Erfmugauss ]
+interface: ezfio,provider,ocaml
+default: Standard
+
+[mu_of_r_av]
+type: logical
+doc: If |true|, take the second formula for mu(r)
+interface: ezfio,provider,ocaml
+default: False
 
 [j2e_type]
 type: character*(32)
-doc: type of the 2e-Jastrow: [ None | Mu | Mu_Nu | Mur | Boys | Boys_Handy | Qmckl ]
+doc: type of the 2e-Jastrow: [ None | Mu | Mugauss | Mu_Nu | Mur | Murgauss | Bump | Boys | Boys_Handy | Qmckl ]
 interface: ezfio,provider,ocaml
 default: Mu
 

plugins/local/non_h_ints_mu/NEED

@@ -4,3 +4,4 @@ jastrow
 ao_tc_eff_map
 bi_ortho_mos
 trexio
+mu_of_r

plugins/local/non_h_ints_mu/deb_deriv_mu.irp.f

@@ -0,0 +1,28 @@
+program test_j_mu_of_r
+ implicit none
+ double precision :: x,mu_min,dmu,mu_max, mu, mu_p, mu_m
+ double precision :: j_simple,j_p, j_m,numeric_d_mu,d_dx_mu
+ double precision :: accu
+ integer :: npt,i
+ npt = 1000
+ mu_min = 0.3d0
+ mu_max = 10.d0
+ dmu = (mu_max - mu_min)/dble(npt)
+ x = 0.7d0
+ mu = mu_min
+ do i = 1, npt
+  call get_deriv_mu_j12(x,mu,d_dx_mu)
+  mu_p = mu + dmu
+  mu_m = mu - dmu
+  j_p = j_simple(x,mu_p)
+  j_m = j_simple(x,mu_m)
+  numeric_d_mu = 0.5d0 * (j_p - j_m)/dmu
+  print*,mu
+  print*,numeric_d_mu,d_dx_mu,dabs(d_dx_mu-numeric_d_mu)
+  accu += dabs(d_dx_mu-numeric_d_mu)
+  mu += dmu
+ enddo
+ accu *= dmu
+ print*,'accu = ',accu
+end
+

plugins/local/non_h_ints_mu/deb_j_bump.irp.f

@@ -0,0 +1,98 @@
+program test_j_mu_of_r
+ implicit none
+! call routine_test_mu_of_r
+ call routine_test_mu_of_r_tot
+end
+
+subroutine routine_test_mu_of_r_tot
+ implicit none
+ integer :: ipoint,k 
+ double precision :: r2(3), weight, dr, r1(3), r1bis(3)
+ double precision :: accu_grad(3)
+ double precision :: jast,grad_jast_mu_r1(3),j_bump
+ double precision :: jast_p,jast_m,num_grad_jast_mu_r1(3)
+
+ dr = 0.00001d0
+ r2 = 0.d0
+ r2(1) =  0.5d0
+ r2(2) = -0.1d0
+ r2(3) =  1.0d0
+ accu_grad = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r1(1:3) = final_grid_points(1:3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+!   call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast_mu_r1)
+   call get_grad_j_bump_mu_of_r(r1,r2,grad_jast_mu_r1)
+   double precision :: norm,error
+   norm = 0.D0
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    jast_p = j_bump(r1bis,r2,a_boys)
+
+    r1bis= r1
+    r1bis(k) -= dr
+    jast_m = j_bump(r1bis,r2,a_boys)
+
+    num_grad_jast_mu_r1(k) = (jast_p - jast_m)/(2.d0* dr)
+    norm += num_grad_jast_mu_r1(k)*num_grad_jast_mu_r1(k)
+   enddo
+   error = 0.d0
+   do k = 1, 3
+    error += dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
+   enddo
+   error *= 0.33333333d0
+   norm = dsqrt(norm)
+   if(norm.gt.1.d-05)then
+    if(dabs(error/norm).gt.dr)then
+     print*,'/////'
+     print*,error,norm
+     print*,grad_jast_mu_r1
+     print*,num_grad_jast_mu_r1
+    endif
+   endif
+   do k = 1,3
+    accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end
+
+subroutine routine_test_mu_of_r
+ implicit none
+ integer :: ipoint,k 
+ double precision :: weight, dr, r1(3), r1bis(3),accu_grad(3),num_grad_mu_r1(3)
+ double precision :: mu_r1,dm_r1, mu_der_r1(3), grad_dm_r1(3)
+ double precision :: mu_der_rp(3), grad_dm_rp(3),mu_rp
+ double precision :: mu_der_rm(3), grad_dm_rm(3),mu_rm
+
+ dr = 0.0001d0
+ accu_grad = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r1(1:3) = final_grid_points(1:3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+   call grad_mu_of_r_mean_field(r1,mu_r1,dm_r1, mu_der_r1, grad_dm_r1) 
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    call grad_mu_of_r_mean_field(r1bis,mu_rp, dm_r1, mu_der_rp, grad_dm_r1) 
+
+    r1bis= r1
+    r1bis(k) -= dr
+    call grad_mu_of_r_mean_field(r1bis,mu_rm, dm_r1, mu_der_rm, grad_dm_r1) 
+    num_grad_mu_r1(k) = (mu_rp - mu_rm)/(2.d0* dr)
+!    print*,jast_mu_r1_p,jast_mu_r1_m
+   enddo
+   print*,'/////'
+   print*,mu_der_r1
+   print*,num_grad_mu_r1
+   do k = 1,3
+    accu_grad(k) += weight * dabs(mu_der_r1(k) - num_grad_mu_r1(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end

plugins/local/non_h_ints_mu/deb_j_gauss.irp.f

@@ -0,0 +1,62 @@
+program test_j_mu_of_r
+ implicit none
+! call routine_test_mu_of_r
+ call routine_test_mu_of_r_tot
+end
+
+subroutine routine_test_mu_of_r_tot
+ implicit none
+ integer :: ipoint,k 
+ double precision :: r2(3), weight, dr, r1(3), r1bis(3)
+ double precision :: accu_grad(3)
+ double precision :: jast,grad_jast(3),j_bump,j12_mu
+ double precision :: jast_p,jast_m,num_grad_jast(3)
+
+ dr = 0.00001d0
+ r2 = 0.d0
+ r2(1) =  0.5d0
+ r2(2) = -0.1d0
+ r2(3) =  1.0d0
+ accu_grad = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r1(1:3) = final_grid_points(1:3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+   call grad1_j12_mu(r1, r2, grad_jast)
+   grad_jast = - grad_jast
+   double precision :: norm,error
+   norm = 0.D0
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    jast_p = j12_mu(r1bis, r2)
+
+    r1bis= r1
+    r1bis(k) -= dr
+    jast_m = j12_mu(r1bis, r2)
+
+    num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr)
+    norm += num_grad_jast(k)*num_grad_jast(k)
+   enddo
+   error = 0.d0
+   do k = 1, 3
+    error += dabs(grad_jast(k) - num_grad_jast(k))
+   enddo
+   error *= 0.33333333d0
+   norm = dsqrt(norm)
+   if(norm.gt.1.d-05)then
+    if(dabs(error/norm).gt.dr)then
+     print*,'/////'
+     print*,error,norm
+     print*,grad_jast
+     print*,num_grad_jast
+    endif
+   endif
+   do k = 1,3
+    accu_grad(k) += weight * dabs(grad_jast(k) - num_grad_jast(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end
+

plugins/local/non_h_ints_mu/deb_j_mu_of_r.irp.f

@@ -0,0 +1,97 @@
+program test_j_mu_of_r
+ implicit none
+! call routine_test_mu_of_r
+ call routine_test_mu_of_r_tot
+end
+
+subroutine routine_test_mu_of_r_tot
+ implicit none
+ integer :: ipoint,k 
+ double precision :: r2(3), weight, dr, r1(3), r1bis(3)
+ double precision :: accu_grad(3)
+ double precision :: jast,grad_jast_mu_r1(3)
+ double precision :: jast_p,jast_m,num_grad_jast_mu_r1(3)
+
+ dr = 0.000001d0
+ r2 = 0.d0
+ r2(1) =  0.5d0
+ r2(2) = -0.1d0
+ r2(3) =  1.0d0
+ accu_grad = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r1(1:3) = final_grid_points(1:3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+   call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast_mu_r1)
+   double precision :: norm,error
+   norm = 0.D0
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    call get_j_sum_mu_of_r(r1bis,r2,jast_p)
+
+    r1bis= r1
+    r1bis(k) -= dr
+    call get_j_sum_mu_of_r(r1bis,r2,jast_m)
+
+    num_grad_jast_mu_r1(k) = (jast_p - jast_m)/(2.d0* dr)
+    norm += num_grad_jast_mu_r1(k)*num_grad_jast_mu_r1(k)
+   enddo
+   error = 0.d0
+   do k = 1, 3
+    error += dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
+   enddo
+   error *= 0.33333333d0
+   norm = dsqrt(norm)
+   if(norm.gt.1.d-05)then
+    if(dabs(error/norm).gt.10.d0*dr)then
+     print*,'/////'
+     print*,error,norm,dabs(error/norm)
+     print*,grad_jast_mu_r1
+     print*,num_grad_jast_mu_r1
+    endif
+   endif
+   do k = 1,3
+    accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end
+
+subroutine routine_test_mu_of_r
+ implicit none
+ integer :: ipoint,k 
+ double precision :: weight, dr, r1(3), r1bis(3),accu_grad(3),num_grad_mu_r1(3)
+ double precision :: mu_r1,dm_r1, mu_der_r1(3), grad_dm_r1(3)
+ double precision :: mu_der_rp(3), grad_dm_rp(3),mu_rp
+ double precision :: mu_der_rm(3), grad_dm_rm(3),mu_rm
+
+ dr = 0.0001d0
+ accu_grad = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r1(1:3) = final_grid_points(1:3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+   call grad_mu_of_r_mean_field(r1,mu_r1,dm_r1, mu_der_r1, grad_dm_r1) 
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    call grad_mu_of_r_mean_field(r1bis,mu_rp, dm_r1, mu_der_rp, grad_dm_r1) 
+
+    r1bis= r1
+    r1bis(k) -= dr
+    call grad_mu_of_r_mean_field(r1bis,mu_rm, dm_r1, mu_der_rm, grad_dm_r1) 
+    num_grad_mu_r1(k) = (mu_rp - mu_rm)/(2.d0* dr)
+!    print*,jast_mu_r1_p,jast_mu_r1_m
+   enddo
+   print*,'/////'
+   print*,mu_der_r1
+   print*,num_grad_mu_r1
+   do k = 1,3
+    accu_grad(k) += weight * dabs(mu_der_r1(k) - num_grad_mu_r1(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end

plugins/local/non_h_ints_mu/deb_j_psi.irp.f

@@ -0,0 +1,131 @@
+program test_j_mu_of_r
+ implicit none
+ call routine_deb_j_psi
+! call routine_deb_denom
+end
+
+subroutine routine_deb_j_psi
+ implicit none
+ integer :: ipoint,k 
+ double precision :: r2(3), weight, dr, r1(3), r1bis(3)
+ double precision :: accu_grad(3)
+ double precision :: jast,grad_jast(3),j_bump,jastrow_psi,grad_jast_bis(3)
+ double precision :: jast_p,jast_m,num_grad_jast(3)
+
+ dr = 0.00001d0
+ r2 = 0.d0
+ r2(1) =  0.5d0
+ r2(2) = -0.1d0
+ r2(3) =  1.0d0
+ accu_grad = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r1(1:3) = final_grid_points(1:3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+   call get_grad_r1_jastrow_psi(r1,r2,grad_jast,jast)
+!   grad_jast = - grad_jast
+   double precision :: norm,error
+   norm = 0.D0
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    call get_grad_r1_jastrow_psi(r1bis,r2,grad_jast_bis,jast_p)
+
+    r1bis= r1
+    r1bis(k) -= dr
+    call get_grad_r1_jastrow_psi(r1bis,r2,grad_jast_bis,jast_m)
+
+    num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr)
+    norm += num_grad_jast(k)*num_grad_jast(k)
+   enddo
+   error = 0.d0
+   do k = 1, 3
+    error += dabs(grad_jast(k) - num_grad_jast(k))
+   enddo
+   error *= 0.33333333d0
+   norm = dsqrt(norm)
+   if(norm.gt.1.d-05)then
+    if(dabs(error/norm).gt.dr)then
+     print*,'/////'
+     print*,error,norm
+     print*,grad_jast
+     print*,num_grad_jast
+    endif
+   endif
+   do k = 1,3
+    accu_grad(k) += weight * dabs(grad_jast(k) - num_grad_jast(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end
+
+
+subroutine routine_deb_denom
+ implicit none
+ integer :: ipoint,k,i,j
+ double precision :: r2(3), weight, dr, r1(3), r1bis(3)
+ double precision :: accu_grad(3)
+ double precision :: jast,grad_jast(3),j_bump,jastrow_psi,grad_jast_bis(3)
+ double precision :: jast_p,jast_m,num_grad_jast(3)
+
+ dr = 0.00001d0
+ r2 = 0.d0
+ r2(1) =  0.5d0
+ r2(2) = -0.1d0
+ r2(3) =  1.0d0
+ double precision, allocatable :: mos_array_r1(:), mos_array_r2(:)
+ double precision, allocatable :: mos_grad_array_r1(:,:),mos_grad_array_r2(:,:)
+ allocate(mos_array_r1(mo_num), mos_array_r2(mo_num))
+ allocate(mos_grad_array_r1(3,mo_num), mos_grad_array_r2(3,mo_num))
+ do i = 1, 1
+  do j = 1, 1
+   accu_grad = 0.d0
+   call give_all_mos_and_grad_at_r(r2,mos_array_r2,mos_grad_array_r2)
+    do ipoint = 1, n_points_final_grid
+     r1(1:3) = final_grid_points(1:3,ipoint)
+     weight = final_weight_at_r_vector(ipoint)
+     call give_all_mos_and_grad_at_r(r1,mos_array_r1,mos_grad_array_r1)
+     call denom_jpsi(i,j,a_boys, mos_array_r1,mos_grad_array_r1,mos_array_r2,jast, grad_jast)
+     double precision :: norm,error
+     norm = 0.D0
+     do k = 1, 3
+      r1bis= r1
+      r1bis(k) += dr
+      call give_all_mos_and_grad_at_r(r1bis,mos_array_r1,mos_grad_array_r1)
+      call denom_jpsi(i,j,a_boys, mos_array_r1,mos_grad_array_r1,mos_array_r2,jast_p, grad_jast_bis)
+  
+      r1bis= r1
+      r1bis(k) -= dr
+      call give_all_mos_and_grad_at_r(r1bis,mos_array_r1,mos_grad_array_r1)
+      call denom_jpsi(i,j,a_boys, mos_array_r1,mos_grad_array_r1,mos_array_r2,jast_m, grad_jast_bis)
+  
+      num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr)
+      norm += num_grad_jast(k)*num_grad_jast(k)
+     enddo
+     error = 0.d0
+     do k = 1, 3
+      error += dabs(grad_jast(k) - num_grad_jast(k))
+     enddo
+     error *= 0.33333333d0
+     norm = dsqrt(norm)
+     if(norm.gt.1.d-05)then
+      if(dabs(error/norm).gt.dr)then
+       print*,'/////'
+       print*,error,norm
+       print*,grad_jast
+       print*,num_grad_jast
+      endif
+     endif
+     do k = 1,3
+      accu_grad(k) += weight * dabs(grad_jast(k) - num_grad_jast(k))
+     enddo
+    enddo
+   print*,'i,j = ',i,j
+   print*,'accu_grad = '
+   print*, accu_grad
+  enddo
+ enddo
+
+end
+

plugins/local/non_h_ints_mu/j_bump.irp.f

@@ -0,0 +1,90 @@
+double precision function wigner_radius(rho)
+ implicit none
+ include 'constants.include.F'
+ double precision, intent(in) :: rho
+ wigner_radius = 4.d0 * pi * rho * 0.333333333333d0
+ wigner_radius = wigner_radius**(-0.3333333d0)
+end
+
+double precision function j_bump(r1,r2,a)
+ implicit none
+ include 'constants.include.F'
+ double precision, intent(in) :: r1(3),r2(3),a
+ double precision :: inv_a,factor,x_scaled,scalar
+ double precision :: r12
+ r12   = (r1(1) - r2(1))*(r1(1) - r2(1))
+ r12  += (r1(2) - r2(2))*(r1(2) - r2(2))
+ r12  += (r1(3) - r2(3))*(r1(3) - r2(3))
+ r12 = dsqrt(r12)
+ inv_a = 1.d0/a
+ x_scaled = r12*inv_a*inv_sq_pi
+ x_scaled*= x_scaled
+ j_bump = 0.5d0 * (r12-a) * dexp(-x_scaled)
+end
+
+subroutine get_grad_j_bump(x,a,grad)
+ implicit none
+ BEGIN_DOC
+ ! gradient of the Jastrow with a bump
+ !
+ ! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2]
+ !
+ ! d/dx j(x,a) = 1/(2 pi a^2) * exp[-(x/(a*sqrt(pi)))^2] * (pi a^2 + 2 a x - 2x^2)
+ END_DOC
+  include 'constants.include.F'
+ double precision, intent(in)  :: x,a
+ double precision, intent(out) :: grad
+ double precision :: inv_a,factor,x_scaled,scalar
+ inv_a = 1.d0/a
+ factor = 0.5d0*inv_pi*inv_a*inv_a 
+ x_scaled = x*inv_a*inv_sq_pi
+ x_scaled*= x_scaled
+ grad = factor * dexp(-x_scaled) * (pi*a*a + 2.d0 * a*x - 2.d0*x*x)
+end
+
+subroutine get_d_da_j_bump(x,a,d_da)
+ implicit none
+ BEGIN_DOC
+ ! Derivative with respect by to the parameter "a" of the Jastrow with a bump
+ !
+ ! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2]
+ !
+ ! d/da j(x,a) = - 1/(pi*a^3) * exp[-(x/(a*sqrt(pi)))^2] * (-2 x^3 + 2 a x^2 + pi a^x3)
+ END_DOC
+ include 'constants.include.F'
+ double precision, intent(in)  :: x,a
+ double precision, intent(out) :: d_da
+ double precision :: factor, inv_a,x_scaled,scalar
+ inv_a = 1.d0/a
+ factor = inv_a*inv_a*inv_a*inv_pi
+ x_scaled = x*inv_a*inv_sq_pi
+ x_scaled*= x_scaled
+ d_da = factor * dexp(-x_scaled) * (-2.d0 * x*x*x + 2.d0*x*x*a+pi*a*a*a)
+end
+
+subroutine get_grad_j_bump_mu_of_r(r1,r2,grad_j_bump)
+ implicit none
+ BEGIN_DOC
+ ! d/dx1 j(x,a(r1,r2)) where j(x,a) is the Jastrow with a bump
+ !
+ ! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2]
+ !
+ ! a(r1,r2) = [rho(r1) a(r1) + rho(r2) a(r2)]/[rho(r1) + rho(r2)]
+ !
+ ! d/dx1 j(x,a) = d/dx1 j(x,a(r1,r2))
+ END_DOC
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: grad_j_bump(3)
+ double precision :: r12,r12_vec(3),grad_scal,inv_r12
+ r12_vec = r1 - r2
+ r12   = (r1(1) - r2(1))*(r1(1) - r2(1))
+ r12  += (r1(2) - r2(2))*(r1(2) - r2(2))
+ r12  += (r1(3) - r2(3))*(r1(3) - r2(3))
+ r12 = dsqrt(r12)
+ call get_grad_j_bump(r12,a_boys,grad_scal)
+ if(r12.lt.1.d-10)then
+  grad_j_bump = 0.d0
+ else 
+  grad_j_bump = grad_scal * r12_vec/r12
+ endif
+end

plugins/local/non_h_ints_mu/jast_deriv.irp.f

@@ -31,7 +31,7 @@
   grad1_u12_squared_num = 0.d0
 
   if( ((j2e_type .eq. "Mu") .and. (env_type .eq. "None")) .or. &
-       (j2e_type .eq. "Mur") ) then
+       (j2e_type .eq. "Mur").or.(j2e_type .eq. "Mugauss") .or. (j2e_type .eq. "Murgauss")) then
 
     !$OMP PARALLEL                                                                                    &
     !$OMP DEFAULT (NONE)                                                                              &

plugins/local/non_h_ints_mu/jast_deriv_mu_of_r.irp.f

@@ -0,0 +1,306 @@
+subroutine get_j_sum_mu_of_r(r1,r2,jast)
+ implicit none
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: jast
+ double precision :: mu_r1, dm_r1, grad_mu_r1(3), grad_dm_r1(3), j_mu_r1
+ double precision :: mu_r2, dm_r2, grad_mu_r2(3), grad_dm_r2(3), j_mu_r2
+ double precision :: j12_mu_input,mu_tot,r12,j_simple
+ jast = 0.d0
+ if(murho_type==0)then
+! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)]
+  call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) 
+  call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) 
+  j_mu_r1 = j12_mu_input(r1, r2, mu_r1)
+  j_mu_r2 = j12_mu_input(r1, r2, mu_r2)
+  if(dm_r1 + dm_r2.lt.1.d-7)return
+  jast = (dm_r1 * j_mu_r1 + dm_r2 * j_mu_r2) / (dm_r1 + dm_r2)
+ else if(murho_type==1)then
+! J(r1,r2) = j(0.5 * (mu(r1)+mu(r2)),r12), MU(r1,r2) = 0.5 *(mu(r1)+mu(r2))
+  call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) 
+  call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) 
+  mu_tot = 0.5d0 * (mu_r1 + mu_r2)
+  jast = j12_mu_input(r1, r2, mu_tot)
+ else if(murho_type==2)then
+! MU(r1,r2) = (rho(1) * mu(r1)+ rho(2) * mu(r2))/(rho(1)+rho(2))
+! J(r1,r2) = j(MU(r1,r2),r12)
+  call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) 
+  call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) 
+  double precision :: mu_tmp, dm_tot, dm_tot_inv
+  dm_tot = dm_r1**a_boys + dm_r2**a_boys  ! rho(1)**alpha+rho(2)**alpha
+  if(dm_tot.lt.1.d-12)then
+   dm_tot_inv = 1.d+12
+  else
+   dm_tot_inv = 1.d0/dm_tot
+  endif
+  mu_tmp = dm_r1**a_boys * mu_r1 + dm_r2**a_boys * mu_r2 !rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2)
+  mu_tot = nu_erf * mu_tmp*dm_tot_inv ! 
+  r12  = (r1(1) - r2(1)) * (r1(1) - r2(1))
+  r12 += (r1(2) - r2(2)) * (r1(2) - r2(2))
+  r12 += (r1(3) - r2(3)) * (r1(3) - r2(3))
+  r12 = dsqrt(r12)
+  jast = j_simple(r12,mu_tot)
+ endif
+
+end
+
+subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
+ implicit none
+  include 'constants.include.F'
+ BEGIN_DOC
+ END_DOC
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: jast, grad_jast(3)
+ jast = 0.d0
+ grad_jast = 0.d0
+ double precision :: num, denom, grad_num(3), grad_denom(3)
+ double precision :: j_r1, grad_j_r1(3),j_r2, grad_j_r2(3)
+ double precision :: dm_r1, grad_dm_r1(3), grad_jmu_r2(3)
+ double precision :: dm_r2, grad_dm_r2(3),mu_r2, grad_mu_r2(3),mu_r1
+ double precision :: j12_mu_input,r12,grad_mu_r1(3),grad_jmu_r1(3)
+ double precision :: mu_tot,dx,dy,dz,r12_vec(3),d_dmu_j,d_dr12_j
+
+ dx = r1(1) - r2(1)
+ dy = r1(2) - r2(2)
+ dz = r1(3) - r2(3)
+
+ r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+ if(r12.gt.1.d-10)then
+  r12_vec(1) = dx
+  r12_vec(2) = dy
+  r12_vec(3) = dz
+  r12_vec *= 1.d0/r12
+  ! r12_vec = grad_r1 (r12)
+ else
+  r12 = 1.d-10
+  r12_vec = 0.d0
+ endif
+
+ if(murho_type==0)then
+! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)]
+! 
+!          = num(r1,r2) / denom(r1,r2)
+!
+! d/dx1 J(r1,r2) = [denom(r1,r2) X d/dx1 num(r1,r2) - num(r1,r2) X d/dx1 denom(r1,r2) ] / denom(r1,r2)^2
+!
+! d/dx1 num(r1,r2) =  j(mu(r1),r12)*d/dx1 rho(r1) + rho(r1) * d/dx1 j(mu(r1),r12) 
+!                   + rho(r2) d/dx1 j(mu(r2),r12)
+! d/dx1 denom(r1,r2) = d/dx1 rho(r1)
+  call grad_j_mu_of_r_1(r1,r2,j_r1, grad_j_r1,dm_r1, grad_dm_r1)
+  call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) 
+  j_r2 = j12_mu_input(r1, r2, mu_r2) ! j(mu(r2),r1,r2)
+  num = dm_r1 * j_r1 + dm_r2 * j_r2
+  denom = dm_r1 + dm_r2
+  if(denom.lt.1.d-7)return
+  jast = num / denom
+ 
+  grad_denom = grad_dm_r1
+  call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2,r12)
+  grad_num =  j_r1 * grad_dm_r1 + dm_r1 * grad_j_r1 + dm_r2 * grad_jmu_r2
+  grad_jast = (grad_num * denom - num * grad_denom)/(denom*denom)
+ else if(murho_type==1)then
+! J(r1,r2) = j(0.5 * (mu(r1)+mu(r2)),r12), MU(r1,r2) = 0.5 *(mu(r1)+mu(r2))
+!
+! d/dx1 J(r1,r2) = d/dx1 j(MU(r1,r2),r12)|MU=cst 
+!                + d/dMU [j(MU,r12)]
+!                x d/d(mu(r1)) MU(r1,r2)
+!                x d/dx1 mu(r1)
+!                = 0.5 * (1 - erf(MU(r1,r2) *r12))/r12 * (x1 - x2) == grad_jmu_r1
+!                + e^{-(r12*MU(r1,r2))^2}/(2 sqrt(pi) * MU(r1,r2)^2) 
+!                x 0.5 
+!                x d/dx1 mu(r1)
+ call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1)
+ call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2)
+ mu_tot = 0.5d0 * (mu_r1 + mu_r2)
+ call grad_j12_mu_input(r1, r2, mu_tot, grad_jmu_r1,r12)
+ grad_jast = grad_jmu_r1
+ grad_jast+= dexp(-r12*mu_tot*r12*mu_tot) * inv_sq_pi_2 /(mu_tot* mu_tot) * 0.5d0 * grad_mu_r1
+ else if(murho_type==2)then
+! MU(r1,r2) = beta * (rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2))/(rho(1)**alpha+rho(2)**alpha)
+! J(r1,r2) = j(MU(r1,r2),r12)
+!
+! d/dx1 J(r1,r2) = d/dx1 j(MU(r1,r2),r12)|MU=cst 
+!                + d/dMU [j(MU,r12)] 
+!                x d/d(mu(r1)) MU(r1,r2)
+!                x d/dx1 mu(r1)
+!                = 0.5 * (1 - erf(MU(r1,r2) *r12))/r12 * (x1 - x2) == grad_jmu_r1
+!                + 0.5 e^{-(r12*MU(r1,r2))^2}/(2 sqrt(pi) * MU(r1,r2)^2) 
+!                x d/dx1 MU(r1,r2)
+! with d/dx1 MU(r1,r2) = beta * {[mu(1) d/dx1 [rho(1)**alpha] + rho(1)**alpha * d/dx1 mu(1)](rho(1)**alpha+rho(2)**alpha) 
+!                       - MU(1,2) d/dx1 [rho(1)]**alpha}/(rho(1)**alpha+rho(2)**alpha)^2
+! d/dx1 [rho(1)]**alpha = alpha [rho(1)]**(alpha-1) d/dx1 rho(1)
+!                        
+ call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1)
+ call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2)
+ double precision :: dm_tot,dm_tot_inv,grad_mu_tot(3),mu_tmp,grad_dm_r1_alpha(3),d_dx_j
+ dm_tot = dm_r1**a_boys + dm_r2**a_boys  ! rho(1)**alpha+rho(2)**alpha
+ grad_dm_r1_alpha = a_boys * dm_r1**(a_boys-1) * grad_dm_r1
+ if(dm_tot.lt.1.d-12)then
+  dm_tot_inv = 1.d+12
+ else
+  dm_tot_inv = 1.d0/dm_tot
+ endif
+ mu_tmp = dm_r1**a_boys * mu_r1 + dm_r2**a_boys * mu_r2 !rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2)
+ mu_tot = nu_erf * mu_tmp*dm_tot_inv ! 
+ grad_mu_tot = ( mu_r1 * grad_dm_r1_alpha + dm_r1**a_boys * grad_mu_r1 ) * dm_tot & 
+              -  mu_tmp * grad_dm_r1_alpha
+ grad_mu_tot *= dm_tot_inv * dm_tot_inv * nu_erf
+ call get_deriv_r12_j12(r12,mu_tot,d_dr12_j) ! d/dr12 j(MU(r1,r2,r12)
+ ! d/dx1 j(MU(r1,r2),r12) | MU(r1,r2) = cst
+ ! d/dr12 j(MU(r1,r2,r12) x d/dx1 r12
+ grad_jmu_r1 = d_dr12_j * r12_vec 
+! call grad_j12_mu_input(r1, r2, mu_tot, grad_jmu_r1,r12)
+ grad_jast = grad_jmu_r1
+ ! d/dMU j(MU(r1,r2),r12)
+ call get_deriv_mu_j12(r12,mu_tot,d_dmu_j)
+ grad_jast+= d_dmu_j * grad_mu_tot
+ else if(murho_type==-1)then
+! J(r1,r2) = 0.5 * [j(mu(r1),r12) + j(mu(r2),r12)] 
+!
+! d/dx1 J(r1,r2)  = 0.5 * (d/dx1 j(mu(r1),r12) + d/dx1 j(mu(r2),r12))
+  call grad_j_mu_of_r_1(r1,r2,j_r1, grad_j_r1,dm_r1, grad_dm_r1)
+  call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) 
+  j_r2 = j12_mu_input(r1, r2, mu_r2) ! j(mu(r2),r1,r2)
+  call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2,r12)
+  jast = 0.5d0 * (j_r1 + j_r2)
+  grad_jast = 0.5d0 * (grad_j_r1 + grad_jmu_r2)
+  
+ endif
+
+end
+
+subroutine grad_j_mu_of_r_1(r1,r2,jast, grad_jast, dm_r1, grad_dm_r1)
+ implicit none
+  include 'constants.include.F'
+ BEGIN_DOC
+! grad_r1 of j(mu(r1),r12)
+  !
+  !
+  ! d/dx1 j(mu(r1),r12) = exp(-(mu(r1)*r12)**2) /(2 *sqrt(pi) * mu(r1)**2 ) d/dx1 mu(r1) 
+  !                     + d/dx1 j(mu(r1),r12) 
+  !
+  ! with 
+  !
+  !           j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2}
+  !
+  ! and d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2)
+  !
+  !     d/d mu j(mu,r12) = e^{-(r12*mu)^2}/(2 sqrt(pi) * mu^2)
+  !
+  ! here mu(r1) is obtained by MU MEAN FIELD 
+ END_DOC
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: jast, grad_jast(3),dm_r1, grad_dm_r1(3)
+ double precision :: dx, dy, dz,  r12, mu_der(3)
+ double precision :: mu_tmp, tmp, grad(3), mu_val
+ jast = 0.d0
+ grad = 0.d0
+
+ dx  = r1(1) - r2(1)
+ dy  = r1(2) - r2(2)
+ dz  = r1(3) - r2(3)
+ r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+ ! get mu(r1) == mu_val and its gradient d/dx1 mu(r1) == mu_der 
+ call grad_mu_of_r_mean_field(r1,mu_val, dm_r1, mu_der, grad_dm_r1) 
+ mu_tmp  = mu_val * r12
+ ! evalulation of the jastrow j(mu(r1),r12)
+ jast = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_val
+
+ ! tmp = exp(-(mu(r1)*r12)**2) /(2 *sqrt(pi) * mu(r1)**2 )
+ tmp     = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
+ ! grad = 
+ grad(1) = tmp * mu_der(1)
+ grad(2) = tmp * mu_der(2)
+ grad(3) = tmp * mu_der(3)
+
+ if(r12 .lt. 1d-10) return
+ tmp     = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 ! d/dx1 j(mu(r1),r12) 
+ grad(1) = grad(1) + tmp * dx
+ grad(2) = grad(2) + tmp * dy
+ grad(3) = grad(3) + tmp * dz
+
+ grad_jast = grad 
+end
+
+! ---
+
+double precision function j12_mu_input(r1, r2, mu)
+
+  BEGIN_DOC
+  ! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2}
+  END_DOC
+  include 'constants.include.F'
+
+  implicit none
+  double precision, intent(in) :: r1(3), r2(3), mu
+  double precision             :: mu_tmp, r12
+
+    r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+               + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+               + (r1(3) - r2(3)) * (r1(3) - r2(3)) )
+    mu_tmp = mu * r12
+
+    j12_mu_input = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu
+
+end
+
+subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu,r12)
+ implicit none
+ BEGIN_DOC
+ ! grad_jmu = d/dx1 j(mu,r12) assuming mu=cst(r1)
+ !
+ !          = 0.5/r_12 * (x_1 - x_2) * [1 - erf(mu*r12)]
+ END_DOC
+ double precision, intent(in) :: r1(3), r2(3), mu
+ double precision, intent(out):: grad_jmu(3),r12
+ double precision             :: mu_tmp, dx, dy, dz, grad(3), tmp
+ grad_jmu = 0.d0
+ dx  = r1(1) - r2(1)
+ dy  = r1(2) - r2(2)
+ dz  = r1(3) - r2(3)
+ r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+ if(r12 .lt. 1d-10) return
+ mu_tmp  = mu * r12
+ tmp     = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 ! d/dx1 j(mu(r1),r12) 
+ grad(1) = tmp * dx
+ grad(2) = tmp * dy
+ grad(3) = tmp * dz
+
+ grad_jmu = grad 
+end
+
+subroutine j12_and_grad_j12_mu_input(r1, r2, mu, jmu, grad_jmu)
+ implicit none
+ include 'constants.include.F'
+ BEGIN_DOC
+ ! jmu = j(mu,r12) 
+ ! grad_jmu = d/dx1 j(mu,r12) assuming mu=cst(r1)
+ !
+ !          = 0.5/r_12 * (x_1 - x_2) * [1 - erf(mu*r12)]
+ END_DOC
+ double precision, intent(in) :: r1(3), r2(3), mu
+ double precision, intent(out):: grad_jmu(3), jmu
+ double precision             :: mu_tmp, r12, dx, dy, dz, grad(3), tmp
+ double precision :: erfc_mur12,inv_mu
+ inv_mu = 1.d0/mu
+
+ grad_jmu = 0.d0 ! initialization when r12 --> 0
+ jmu = - inv_sq_pi_2 * inv_mu ! initialization when r12 --> 0
+
+ dx  = r1(1) - r2(1)
+ dy  = r1(2) - r2(2)
+ dz  = r1(3) - r2(3)
+ r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+ if(r12 .lt. 1d-10) return
+ erfc_mur12 = (1.d0 - derf(mu_tmp))
+ mu_tmp  = mu * r12
+ tmp     = 0.5d0 * erfc_mur12  / r12 ! d/dx1 j(mu(r1),r12) 
+ grad(1) = tmp * dx
+ grad(2) = tmp * dy
+ grad(3) = tmp * dz
+
+ grad_jmu = grad 
+
+ jmu= 0.5d0 * r12 * erfc_mur12 - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) * inv_mu
+
+
+end

plugins/local/non_h_ints_mu/jast_deriv_utils.irp.f

@@ -1,8 +1,73 @@
+subroutine get_deriv_r12_j12(x,mu,d_dx_j)
+ implicit none
+  include 'constants.include.F'
+ BEGIN_DOC
+ ! d/dr12 j(mu,r12)
+ END_DOC
+ double precision, intent(in)  :: x,mu
+ double precision, intent(out) :: d_dx_j
+ 
+ d_dx_j = 0.d0
+ if(x .lt. 1d-10) return
+ if(j2e_type .eq. "Mu" .or. j2e_type .eq. "Mur") then
+  d_dx_j = 0.5d0 * (1.d0 - derf(mu * x))
+ else if(j2e_type .eq. "Mugauss" .or. j2e_type .eq. "Murgauss" ) then
+  double precision :: x_tmp
+  x_tmp = mu * x
+  ! gradient of j(mu,x)
+  d_dx_j = 0.5d0 * (1.d0 - derf(x_tmp))
+
+  ! gradient of gaussian additional term
+  x_tmp *= alpha_mu_gauss
+  x_tmp *= x_tmp 
+  d_dx_j += -0.5d0 * mu * c_mu_gauss * x * dexp(-x_tmp)
+ else
+   print *, ' Error in get_deriv_r12_j12: Unknown j2e_type = ', j2e_type
+   stop
+ endif
+end
+
+
+subroutine get_deriv_mu_j12(x,mu,d_d_mu)
+ implicit none
+ BEGIN_DOC
+ ! d/dmu j(mu,r12)
+ END_DOC
+  include 'constants.include.F'
+ double precision, intent(in)  :: x,mu
+ double precision, intent(out) :: d_d_mu
+ double precision :: x_tmp,inv_mu_2,inv_alpha_2
+ 
+ d_d_mu = 0.d0
+ if(x .lt. 1d-10) return
+ x_tmp = x*mu
+ if(mu.lt.1.d-10) return
+ inv_mu_2 = mu*mu
+ inv_mu_2 = 1.d0/inv_mu_2
+ if(j2e_type .eq. "Mu" .or. j2e_type .eq. "Mur") then
+  ! e^{-(r12*mu)^2}/(2 sqrt(pi) * mu^2)
+  d_d_mu = dexp(-x_tmp*x_tmp) * inv_sq_pi_2 * inv_mu_2
+ else if(j2e_type .eq. "Mugauss"  .or. j2e_type .eq. "Murgauss" ) then
+  d_d_mu  = dexp(-x_tmp*x_tmp) * inv_sq_pi_2 * inv_mu_2
+  inv_alpha_2 = 1.d0/alpha_mu_gauss
+  inv_alpha_2 *= inv_alpha_2
+  x_tmp *= alpha_mu_gauss
+  x_tmp *= x_tmp
+  d_d_mu += -0.25d0 * c_mu_gauss*inv_alpha_2*dexp(-x_tmp) * (1.d0 + 2.d0 * x_tmp) * inv_mu_2
+ else
+   print *, ' Error in get_deriv_r12_j12: Unknown j2e_type = ', j2e_type
+   stop
+ endif
+end
+
 
 ! ---
 
 double precision function j12_mu(r1, r2)
 
+  BEGIN_DOC
+  ! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2}
+  END_DOC
   include 'constants.include.F'
 
   implicit none
@@ -18,6 +83,18 @@ double precision function j12_mu(r1, r2)
 
     j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf
 
+  else if(j2e_type .eq. "Mugauss") then
+
+    r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+               + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+               + (r1(3) - r2(3)) * (r1(3) - r2(3)) )
+    double precision :: r12_tmp
+    r12_tmp = mu_erf * r12
+
+    j12_mu = 0.5d0 * r12 * (1.d0 - derf(r12_tmp)) - inv_sq_pi_2 * dexp(-r12_tmp*r12_tmp) / mu_erf
+    r12_tmp *= alpha_mu_gauss
+    j12_mu += 0.25d0 * c_mu_gauss / (alpha_mu_gauss*alpha_mu_gauss*mu_erf) * dexp(-r12_tmp*r12_tmp)
+
   else
 
     print *, ' Error in j12_mu: Unknown j2e_type = ', j2e_type
@@ -57,7 +134,7 @@ subroutine grad1_j12_mu(r1, r2, grad)
 
   grad = 0.d0
 
-  if(j2e_type .eq. "Mu") then
+  if(j2e_type .eq. "Mu".or.j2e_type .eq. "Mugauss") then
 
     dx = r1(1) - r2(1)
     dy = r1(2) - r2(2)
@@ -66,31 +143,42 @@ subroutine grad1_j12_mu(r1, r2, grad)
     r12 = dsqrt(dx * dx + dy * dy + dz * dz)
     if(r12 .lt. 1d-10) return
 
-    tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
+    call get_deriv_r12_j12(r12,mu_erf,tmp)
+!    tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
 
     grad(1) = tmp * dx
     grad(2) = tmp * dy
     grad(3) = tmp * dz
+    grad *= 1.d0/r12
 
-  elseif(j2e_type .eq. "Mur") then
+  elseif(j2e_type .eq. "Mur" .or. j2e_type .eq. "Murgauss") then
+   double precision :: jast
+   call grad_j_sum_mu_of_r(r1,r2,jast,grad)
  
+  elseif(j2e_type .eq. "Bump") then
+   double precision ::grad_jast(3)
+   call get_grad_j_bump_mu_of_r(r1,r2,grad_jast)
      dx = r1(1) - r2(1)
      dy = r1(2) - r2(2)
      dz = r1(3) - r2(3)
+
      r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+     if(r12 .lt. 1d-10) then
+       grad(1) = 0.d0
+       grad(2) = 0.d0
+       grad(3) = 0.d0
+       return
+     endif                                                                                                                                                                            
 
-    call mu_r_val_and_grad(r1, r2, mu_val, mu_der)
+     tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
-    mu_tmp  = mu_val * r12
+
-    tmp     = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
+     grad(1) = 0.5d0 * tmp * dx
-    grad(1) = tmp * mu_der(1)
+     grad(2) = 0.5d0 * tmp * dy
-    grad(2) = tmp * mu_der(2)
+     grad(3) = 0.5d0 * tmp * dz
-    grad(3) = tmp * mu_der(3)
+     grad(1) += 0.5d0 * grad_jast(1)
+     grad(2) += 0.5d0 * grad_jast(2)
+     grad(3) += 0.5d0 * grad_jast(3)
 
-    if(r12 .lt. 1d-10) return
-    tmp     = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
-    grad(1) = grad(1) + tmp * dx
-    grad(2) = grad(2) + tmp * dy
-    grad(3) = grad(3) + tmp * dz
 
   else
 
@@ -369,7 +457,18 @@ end
 ! ---
 
 subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
-
+ BEGIN_DOC
+! various flavours of mu(r1,r2) 
+! depends on essentially the density and other related quantities 
+!
+! change the variable "murho_type" to change type
+!
+!  murho_type == -1 :: mu(r1,r2) = (rho(r1) mu_mf(r1) + rho(r2) mu_mf(r2))/[rho(r1)+rho(r2)]
+!
+!             ==  0 :: mu(r1,r2) = (sqrt(rho(r1)) mu_mf(r1) + sqrt(rho(r2)) mu_mf(r2))/[sqrt(rho(r1))+sqrt(rho(r2))]
+!             
+!             == -2 :: mu(r1,r2) = 0.5(mu_mf(r1) + mu_mf(r2))
+ END_DOC
   implicit none
   double precision, intent(in)  :: r1(3), r2(3)
   double precision, intent(out) :: mu_val, mu_der(3)
@@ -379,11 +478,50 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
   double precision :: rho1, grad_rho1(3),rho2,rho_tot,inv_rho_tot
   double precision :: f_rho1, f_rho2, d_drho_f_rho1
   double precision :: d_dx1_f_rho1(3),d_dx_rho_f_rho(3),nume
+  double precision :: mu_mf_r1, dm_r1, grad_mu_mf_r1(3), grad_dm_r1(3)
+  double precision :: mu_mf_r2, dm_r2, grad_mu_mf_r2(3), grad_dm_r2(3)
+
+  double precision :: num, denom, grad_denom(3), grad_num(3)
+  double precision :: dsqrt_dm_r1
 
   PROVIDE murho_type
   PROVIDE mu_r_ct mu_erf
 
-  if(murho_type .eq. 1) then
+  if(murho_type .eq. 0) then 
+   call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1)
+   call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2)
+   dsqrt_dm_r1 = dsqrt(dm_r1)
+   denom = (dsqrt_dm_r1 + dsqrt(dm_r2) )
+   if(denom.lt.1.d-7)then
+    mu_val = 1.d+10
+    mu_der = 0.d0
+    return
+   endif
+   num = (dsqrt(dm_r1) * mu_mf_r1 + dsqrt(dm_r2) * mu_mf_r2) 
+   mu_val = num / denom
+   grad_denom = grad_dm_r1/dsqrt_dm_r1
+   grad_num = dsqrt(dm_r1) * grad_mu_mf_r1 + mu_mf_r1 * grad_dm_r1
+   mu_der = (grad_num * denom - num * grad_denom)/(denom*denom)
+  else if(murho_type .eq. -1) then
+   call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1)
+   call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2)
+   denom = (dm_r1 + dm_r2 )
+   if(denom.lt.1.d-7)then
+    mu_val = 1.d+10
+    mu_der = 0.d0
+    return
+   endif
+   num = (dm_r1 * mu_mf_r1 + dm_r2 * mu_mf_r2) 
+   mu_val = num / denom
+   grad_denom = grad_dm_r1
+   grad_num = dm_r1 * grad_mu_mf_r1 + mu_mf_r1 * grad_dm_r1
+   mu_der = (grad_num * denom - num * grad_denom)/(denom*denom)
+  else if(murho_type .eq. -2) then
+   call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1)
+   call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2)
+   mu_val = 0.5d0 * (mu_mf_r1 + mu_mf_r2) 
+   mu_der = 0.5d0 * grad_mu_mf_r1
+  else if(murho_type .eq. 1) then
 
     !
     ! r = 0.5 (r1 + r2)

plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f

@@ -35,6 +35,10 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res)
 
   if( (j2e_type .eq. "Mu")        .or. &
       (j2e_type .eq. "Mur")       .or. &
+      (j2e_type .eq. "Jpsi")   .or. &
+      (j2e_type .eq. "Mugauss")   .or. &
+      (j2e_type .eq. "Murgauss")  .or. &
+      (j2e_type .eq. "Bump")      .or. &
       (j2e_type .eq. "Boys") ) then
 
     if(env_type .eq. "None") then
@@ -206,10 +210,9 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
       gradz(jpoint) = tmp * dz
     enddo
 
-  elseif(j2e_type .eq. "Mur") then
+  else if(j2e_type .eq. "Mugauss") then
 
-    ! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2) 
+    !  d/dx1 j(mu,r12) = 0.5 * [(1 - erf(mu * r12)) / r12 - mu*c*r12*exp(-(mu*alpha*r12)^2] * (x1 - x2)
-    !                        + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2)
 
     do jpoint = 1, n_points_extra_final_grid ! r2 
 
@@ -220,15 +223,8 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
       dx = r1(1) - r2(1)
       dy = r1(2) - r2(2)
       dz = r1(3) - r2(3)
+
       r12 = dsqrt(dx * dx + dy * dy + dz * dz)
-
-      call mu_r_val_and_grad(r1, r2, mu_val, mu_der)
-      mu_tmp  = mu_val * r12
-      tmp     = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
-      gradx(jpoint) = tmp * mu_der(1)
-      grady(jpoint) = tmp * mu_der(2)
-      gradz(jpoint) = tmp * mu_der(3)
-
       if(r12 .lt. 1d-10) then
         gradx(jpoint) = 0.d0 
         grady(jpoint) = 0.d0 
@@ -236,11 +232,70 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
         cycle
       endif
 
-      tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
+      double precision :: r12_tmp
+      r12_tmp = mu_erf * r12
+      ! gradient of j(mu,r12)
+      tmp = 0.5d0 * (1.d0 - derf(r12_tmp)) / r12
+      ! gradient of gaussian additional term
+      r12_tmp *= alpha_mu_gauss
+      r12_tmp *= r12_tmp 
+      tmp += -0.5d0 * mu_erf * c_mu_gauss * r12 * dexp(-r12_tmp)/r12
 
-      gradx(jpoint) = gradx(jpoint) + tmp * dx
+      gradx(jpoint) = tmp * dx
-      grady(jpoint) = grady(jpoint) + tmp * dy
+      grady(jpoint) = tmp * dy
-      gradz(jpoint) = gradz(jpoint) + tmp * dz
+      gradz(jpoint) = tmp * dz
+    enddo
+
+  elseif(j2e_type .eq. "Mur".or.j2e_type .eq. "Murgauss") then
+
+    ! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2) 
+    !                        + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2)
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+      double precision :: jast, grad_jast(3)
+      call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
+      gradx(jpoint) = grad_jast(1)
+      grady(jpoint) = grad_jast(2)
+      gradz(jpoint) = grad_jast(3)
+    enddo
+  elseif(j2e_type .eq. "Bump") then
+
+    ! d/dx1 jbump(r1,r2)
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+      call get_grad_j_bump_mu_of_r(r1,r2,grad_jast)
+
+      dx = r1(1) - r2(1)
+      dy = r1(2) - r2(2)
+      dz = r1(3) - r2(3)
+
+      r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+      if(r12 .lt. 1d-10) then
+        gradx(jpoint) = 0.d0 
+        grady(jpoint) = 0.d0 
+        gradz(jpoint) = 0.d0 
+        cycle
+      endif
+
+      tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
+
+      gradx(jpoint) = 0.5d0 * tmp * dx
+      grady(jpoint) = 0.5d0 * tmp * dy
+      gradz(jpoint) = 0.5d0 * tmp * dz
+      gradx(jpoint) += 0.5d0 * grad_jast(1)
+      grady(jpoint) += 0.5d0 * grad_jast(2)
+      gradz(jpoint) += 0.5d0 * grad_jast(3)
+!      gradx(jpoint) =  grad_jast(1)
+!      grady(jpoint) =  grad_jast(2)
+!      gradz(jpoint) =  grad_jast(3)
     enddo
 
   elseif(j2e_type .eq. "Boys") then
@@ -363,6 +418,17 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
       enddo ! i_nucl
     enddo ! jpoint
 
+  elseif(j2e_type .eq. "Jpsi") then
+   double precision :: grad_j_psi_r1(3),jast_psi
+   do jpoint = 1, n_points_extra_final_grid ! r2 
+     r2(1) = final_grid_points_extra(1,jpoint)
+     r2(2) = final_grid_points_extra(2,jpoint)
+     r2(3) = final_grid_points_extra(3,jpoint)
+     call get_grad_r1_jastrow_psi(r1,r2,grad_j_psi_r1,jast_psi)
+     gradx(jpoint) = grad_j_psi_r1(1)
+     grady(jpoint) = grad_j_psi_r1(2)
+     gradz(jpoint) = grad_j_psi_r1(3)
+   enddo
   else
 
     print *, ' Error in grad1_j12_r1_seq: Unknown j2e_type = ', j2e_type
@@ -667,7 +733,11 @@ subroutine get_grad1_u12_2e_r1_seq(ipoint, n_grid2, resx, resy, resz)
   r1(3) = final_grid_points(3,ipoint)
 
   if( (j2e_type .eq. "Mu")        .or. &
+      (j2e_type .eq. "Mugauss")   .or. &
       (j2e_type .eq. "Mur")       .or. &
+      (j2e_type .eq. "Jpsi")       .or. &
+      (j2e_type .eq. "Murgauss")  .or. &
+      (j2e_type .eq. "Bump")      .or. &
       (j2e_type .eq. "Boys") ) then
 
     if(env_type .eq. "None") then
@@ -786,6 +856,9 @@ subroutine get_u12_2e_r1_seq(ipoint, n_grid2, res)
 
   if( (j2e_type .eq. "Mu")        .or. &
       (j2e_type .eq. "Mur")       .or. &
+      (j2e_type .eq. "Mugauss")  .or. &
+      (j2e_type .eq. "Murgauss")  .or. &
+      (j2e_type .eq. "Mugauss")   .or. &
       (j2e_type .eq. "Boys") ) then
 
     if(env_type .eq. "None") then

plugins/local/non_h_ints_mu/jastrow_psi.irp.f

@@ -0,0 +1,124 @@
+BEGIN_PROVIDER [ double precision, c_ij_ab_jastrow, (mo_num, mo_num, elec_alpha_num, elec_beta_num)]
+ implicit none
+ integer :: iunit, getUnitAndOpen
+ c_ij_ab_jastrow = 0.d0
+ iunit = getUnitAndOpen(trim(ezfio_work_dir)//'c_ij_ab', 'R')                                                     
+ read(iunit) c_ij_ab_jastrow
+ close(iunit)
+ print*,'c_ij_ab_jastrow = '
+ integer :: i,j,a,b
+ do i = 1, elec_beta_num ! r2
+  do j = 1, elec_alpha_num ! r1
+   do a = elec_beta_num+1, mo_num ! r2
+    do b = elec_alpha_num+1, mo_num ! r1 
+!     print*,b,a,j,i
+     print*,c_ij_ab_jastrow(b,a,j,i),b,a,j,i
+     if(dabs(c_ij_ab_jastrow(b,a,j,i)).lt.1.d-12)then
+      c_ij_ab_jastrow(b,a,j,i) = 0.d0
+     endif
+    enddo
+   enddo
+  enddo
+ enddo
+END_PROVIDER 
+
+double precision function jastrow_psi(r1,r2)
+ implicit none
+ double precision, intent(in) :: r1(3), r2(3)
+ integer :: i,j,a,b
+ double precision, allocatable :: mos_array_r1(:), mos_array_r2(:)
+ allocate(mos_array_r1(mo_num), mos_array_r2(mo_num))
+ call give_all_mos_at_r(r1,mos_array_r1)
+ call give_all_mos_at_r(r2,mos_array_r2)
+ double precision :: eps,coef, numerator,denominator
+ double precision :: phi_i_phi_j
+ eps = a_boys
+ jastrow_psi= 0.d0
+ do i = 1, elec_beta_num ! r1
+  do j = 1, elec_alpha_num ! r2
+   phi_i_phi_j = mos_array_r1(i) * mos_array_r2(j) + eps
+   denominator = 1.d0/phi_i_phi_j
+   do a = elec_beta_num+1, mo_num ! r1
+    do b = elec_alpha_num+1, mo_num ! r2 
+     coef = c_ij_ab_jastrow(b,a,j,i)
+     numerator = mos_array_r2(b) * mos_array_r1(a)
+     jastrow_psi += coef * numerator*denominator
+    enddo
+   enddo
+  enddo
+ enddo
+end
+
+subroutine get_grad_r1_jastrow_psi(r1,r2,grad_j_psi_r1,jast)
+ implicit none
+ double precision, intent(in) :: r1(3), r2(3)
+ double precision, intent(out):: grad_j_psi_r1(3),jast
+ integer :: i,j,a,b
+ double precision, allocatable :: mos_array_r1(:), mos_array_r2(:)
+ double precision, allocatable :: mos_grad_array_r1(:,:),mos_grad_array_r2(:,:)
+ double precision :: num_j, denom_j, num_j_grad(3), denom_j_grad(3),delta,coef
+ double precision :: inv_denom_j
+ allocate(mos_array_r1(mo_num), mos_array_r2(mo_num))
+ allocate(mos_grad_array_r1(3,mo_num), mos_grad_array_r2(3,mo_num))
+ delta = a_boys
+ call give_all_mos_and_grad_at_r(r1,mos_array_r1,mos_grad_array_r1)
+ call give_all_mos_and_grad_at_r(r2,mos_array_r2,mos_grad_array_r2)
+ grad_j_psi_r1 = 0.d0
+ jast = 0.d0
+ do i = 1, elec_beta_num ! r1
+  do j = 1, elec_alpha_num ! r2
+   call denom_jpsi(i,j,delta,mos_array_r1,mos_grad_array_r1,mos_array_r2,denom_j, denom_j_grad)
+   inv_denom_j = 1.d0/denom_j
+   do a = elec_beta_num+1, mo_num ! r1
+    do b = elec_alpha_num+1, mo_num ! r2 
+     call numerator_psi(a,b,mos_array_r1,mos_grad_array_r1,mos_array_r2,num_j, num_j_grad)
+     coef = c_ij_ab_jastrow(b,a,j,i)
+     jast += coef * num_j * inv_denom_j
+     grad_j_psi_r1 += coef * (num_j_grad * denom_j - num_j * denom_j_grad) * inv_denom_j * inv_denom_j
+    enddo
+   enddo
+  enddo
+ enddo
+ if(jast.lt.-1.d0.or.dabs(jast).gt.1.d0)then
+  print*,'pb ! '
+  print*,jast
+  print*,dsqrt(r1(1)**2+r1(2)**2+r1(3)**2),dsqrt(r2(1)**2+r2(2)**2+r2(3)**2)
+  print*,r1
+!  print*,mos_array_r1(1:2)
+  print*,r2
+!  print*,mos_array_r2(1:2)
+  stop
+ endif
+ if(log_jpsi)then
+  grad_j_psi_r1 = grad_j_psi_r1/(1.d0 + jast)
+ endif
+
+end
+
+
+subroutine denom_jpsi(i,j,delta,mos_array_r1,mos_grad_array_r1,mos_array_r2,denom, grad_denom)
+ implicit none
+ integer, intent(in)           :: i,j
+ double precision, intent(in)  :: mos_array_r1(mo_num),mos_grad_array_r1(3,mo_num),mos_array_r2(mo_num),delta
+ double precision, intent(out) :: denom, grad_denom(3)
+ double precision :: coef,phi_i_phi_j,inv_phi_i_phi_j,inv_phi_i_phi_j_2
+ phi_i_phi_j = mos_array_r1(i) * mos_array_r2(j)
+ if(phi_i_phi_j /= 0.d0)then
+  inv_phi_i_phi_j   = 1.d0/phi_i_phi_j
+  inv_phi_i_phi_j_2 = 1.d0/(phi_i_phi_j * phi_i_phi_j)
+ else
+  inv_phi_i_phi_j   = huge(1.0)
+  inv_phi_i_phi_j_2 = huge(1.d0)
+ endif
+ denom         = phi_i_phi_j + delta * inv_phi_i_phi_j
+ grad_denom(:) = (1.d0 - delta*inv_phi_i_phi_j_2) * mos_array_r2(j) * mos_grad_array_r1(:,i)
+end
+
+subroutine numerator_psi(a,b,mos_array_r1,mos_grad_array_r1,mos_array_r2,num, grad_num)
+ implicit none
+ integer, intent(in)           :: a,b
+ double precision, intent(in)  :: mos_array_r1(mo_num),mos_grad_array_r1(3,mo_num),mos_array_r2(mo_num)
+ double precision, intent(out) :: num, grad_num(3)
+ num = mos_array_r1(a) * mos_array_r2(b)
+ grad_num(:) = mos_array_r2(b) * mos_grad_array_r1(:,a)
+end

plugins/local/non_h_ints_mu/mu_of_r.irp.f

@@ -0,0 +1,43 @@
+
+subroutine grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm)
+ implicit none
+ BEGIN_DOC
+ ! returns the value and gradients of the mu(r) mean field, together with the HF density and its gradients. 
+ END_DOC
+ include 'constants.include.F'
+ double precision, intent(in) :: r(3)
+ double precision, intent(out):: grad_mu_mf(3), grad_dm(3)
+ double precision, intent(out):: mu_mf, dm
+ double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3),grad_dm_a(3), grad_dm_b(3)
+ double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b
+ 
+ double precision :: dist
+ call get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b)
+ 
+ dm = dm_a + dm_b
+ grad_dm(1:3) = grad_dm_a(1:3) + grad_dm_b(1:3)
+
+ if(dabs(two_bod_dens).lt.1.d-10)then
+  mu_mf = 1.d+10
+  grad_mu_mf = 0.d0
+ else
+  if(mu_of_r_tc=="Erfmu")then
+   mu_mf           = 0.3333333333d0 * sqpi * (f_mf_ab/two_bod_dens + 0.25d0)
+   grad_mu_mf(1:3) = 0.3333333333d0 * sqpi * (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& 
+                                    /(two_bod_dens*two_bod_dens)
+  else if(mu_of_r_tc=="Standard")then
+   mu_mf = 0.5d0 * sqpi * f_mf_ab/two_bod_dens
+   grad_mu_mf(1:3) = 0.5d0 * sqpi * (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& 
+                                    /(two_bod_dens*two_bod_dens)
+  else if(mu_of_r_tc=="Erfmugauss")then
+   mu_mf = (f_mf_ab/two_bod_dens + 0.25d0)/c_mu_gauss_tot 
+   grad_mu_mf(1:3) = 1.d0/c_mu_gauss_tot* (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& 
+                                    /(two_bod_dens*two_bod_dens)
+  else 
+   print*,'Wrong value for mu_of_r_tc !'
+   stop
+  endif
+ endif 
+
+end
+

plugins/local/non_h_ints_mu/mu_of_r_mean_field.irp.f

@@ -57,6 +57,9 @@ end
 
 subroutine get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b)
  implicit none
+ BEGIN_DOC
+ ! gradient of mu(r) mean field, together with the gradient of the one- and two-body HF density. 
+ END_DOC
  double precision, intent(in)  :: r(3)
  double precision, intent(out) :: f_mf_ab, two_bod_dens
  double precision, intent(out) :: grad_two_bod_dens(3), grad_f_mf_ab(3)
@@ -146,26 +149,18 @@ subroutine mu_of_r_mean_field(r,mu_mf, dm)
  endif
 end
 
-subroutine grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm)
+subroutine mu_of_r_mean_field_tc(r,mu_mf, dm)
  implicit none
   include 'constants.include.F'
  double precision, intent(in) :: r(3)
- double precision, intent(out):: grad_mu_mf(3), grad_dm(3)
  double precision, intent(out):: mu_mf, dm
- double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3),grad_dm_a(3), grad_dm_b(3)
  double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b
- call get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b)
+ call get_f_mf_ab(r,f_mf_ab,two_bod_dens, dm_a, dm_b)
-
  dm = dm_a + dm_b
- grad_dm(1:3) = grad_dm_a(1:3) + grad_dm_b(1:3)
-
  if(dabs(two_bod_dens).lt.1.d-10)then
   mu_mf = 1.d+10
-  grad_mu_mf = 0.d0
  else
-  mu_mf = 0.5d0 * sqpi * f_mf_ab/two_bod_dens
+  mu_mf = 0.3333333333d0 * sqpi * (f_mf_ab/two_bod_dens + 0.25d0)
-  grad_mu_mf(1:3) = 0.5d0 * sqpi * (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& 
-                                   /(two_bod_dens*two_bod_dens)
  endif
-
 end
+

plugins/local/non_h_ints_mu/plot_j_gauss.irp.f

@@ -0,0 +1,59 @@
+program plot_j_gauss
+ implicit none
+ double precision :: xmin, xmax, x, dx
+ double precision :: mu_min, mu_max, mu, d_mu
+ double precision :: pot_j_gauss,j_mu_simple,j_gauss_simple,pot_j_mu
+ double precision, allocatable :: mu_tab(:),j_mu(:),j_mu_gauss(:)
+ double precision, allocatable :: w_mu(:), w_mu_gauss(:)
+
+ character*(128) :: output
+ integer :: getUnitAndOpen
+ integer :: i_unit_output_wee_gauss,i_unit_output_wee_mu
+ integer :: i_unit_output_j_gauss,i_unit_output_j_mu
+ output=trim(ezfio_filename)//'.w_ee_mu_gauss'
+ i_unit_output_wee_gauss = getUnitAndOpen(output,'w')
+ output=trim(ezfio_filename)//'.w_ee_mu'
+ i_unit_output_wee_mu = getUnitAndOpen(output,'w')
+ output=trim(ezfio_filename)//'.j_mu_gauss'
+ i_unit_output_j_gauss = getUnitAndOpen(output,'w')
+ output=trim(ezfio_filename)//'.j_mu'
+ i_unit_output_j_mu = getUnitAndOpen(output,'w')
+
+ integer :: npt, i, j, n_mu
+ n_mu = 3
+ allocate(mu_tab(n_mu),j_mu(n_mu),j_mu_gauss(n_mu),w_mu(n_mu), w_mu_gauss(n_mu))
+ mu_min = 0.5d0
+ mu_max = 2.d0
+ d_mu = (mu_max - mu_min)/dble(n_mu)
+ mu = mu_min
+ do i = 1, n_mu
+  mu_tab(i) = mu
+  print*,'mu = ',mu
+  mu += d_mu
+ enddo
+ mu_tab(1) = 0.9d0
+ mu_tab(2) = 0.95d0
+ mu_tab(3) = 1.d0
+
+ xmin = 0.01d0
+ xmax = 10.d0
+ npt = 1000
+ dx = (xmax - xmin)/dble(npt)
+ x = xmin
+ do i = 1, npt
+  do j = 1, n_mu
+   mu = mu_tab(j)
+   w_mu_gauss(j) = pot_j_gauss(x,mu)
+   w_mu(j) = pot_j_mu(x,mu)
+   j_mu(j) = j_mu_simple(x,mu)
+   j_mu_gauss(j) = j_gauss_simple(x,mu) + j_mu(j)
+  enddo
+  write(i_unit_output_wee_gauss,'(100(F16.10,X))')x,w_mu_gauss(:)
+  write(i_unit_output_wee_mu,'(100(F16.10,X))')x,w_mu(:)
+  write(i_unit_output_j_gauss,'(100(F16.10,X))')x,j_mu_gauss(:)
+  write(i_unit_output_j_mu,'(100(F16.10,X))')x,j_mu(:)
+  x += dx
+ enddo
+
+
+end

plugins/local/non_h_ints_mu/plot_mo.irp.f

@@ -0,0 +1,19 @@
+program plot_mo
+ implicit none
+ integer :: i,npt
+ double precision :: xmin,xmax,dx,r(3)
+ double precision,allocatable :: mos_array(:)
+ allocate(mos_array(mo_num))
+ npt = 10000
+ xmin =0.d0
+ xmax =10.d0
+ dx=(xmax-xmin)/dble(npt)
+ r=0.d0
+ r(1) = xmin
+ do i = 1, npt
+  call give_all_mos_at_r(r,mos_array)
+  write(33,'(100(F16.10,X))')r(1),mos_array(1),mos_array(2),mos_array(3)
+  r(1) += dx
+ enddo
+
+end

plugins/local/non_h_ints_mu/plot_mu_of_r.irp.f

@@ -16,15 +16,16 @@ subroutine routine_print
  integer :: ipoint,nx,i
  double precision :: xmax,xmin,r(3),dx,sigma
  double precision :: mu_val, mu_der(3),dm_a,dm_b,grad,grad_dm_a(3), grad_dm_b(3)
- xmax =  5.D0
+ xmax =  3.9D0
- xmin = -5.D0
+ xmin = -3.9D0
  nx = 10000
  dx = (xmax - xmin)/dble(nx)
  r = 0.d0
  r(1) = xmin
  do ipoint = 1, nx
-  call mu_r_val_and_grad(r, r, mu_val, mu_der)
+!  call mu_r_val_and_grad(r, r, mu_val, mu_der)
-  call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)
+  call grad_mu_of_r_mean_field(r,mu_val, dm_a, mu_der, grad_dm_a)
+!  call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)
   sigma = 0.d0
   do i = 1,3
    sigma += grad_dm_a(i)**2
@@ -32,7 +33,8 @@ subroutine routine_print
   sigma=dsqrt(sigma)
   grad = mu_der(1)**2 + mu_der(2)**2 + mu_der(3)**2 
   grad = dsqrt(grad)
-  write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a+dm_b,grad,sigma/dm_a
+  print*,r(1),mu_val
+  write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a,grad,sigma/dm_a
   r(1) += dx
  enddo
 end

plugins/local/non_h_ints_mu/pot_j_gauss.irp.f

@@ -0,0 +1,146 @@
+double precision function j_simple(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ double precision :: j_mu_simple,j_gauss_simple
+ if(j2e_type .eq. "Mu".or.j2e_type .eq. "Mur") then
+  j_simple = j_mu_simple(x,mu)
+ else if(j2e_type .eq. "Mugauss".or.j2e_type .eq. "Murgauss") then
+  j_simple = j_gauss_simple(x,mu) + j_mu_simple(x,mu)
+ endif
+end
+
+
+double precision function j_mu_simple(x,mu)
+ implicit none
+ double precision, intent(in):: x,mu
+ include 'constants.include.F'
+ BEGIN_DOC
+! j_mu(mu,x) = 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2)
+ END_DOC
+ j_mu_simple = 0.5d0 * x * (1.D0 - derf(mu*x)) - 0.5d0 * inv_sq_pi/mu *  dexp(-x*mu*x*mu)
+
+end
+
+double precision function j_gauss_simple(x,mu)
+ implicit none
+ double precision, intent(in):: x,mu
+ include 'constants.include.F'
+ BEGIN_DOC
+! j_mu(mu,x) = c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2)
+!      with c = 27/(8 sqrt(pi)), alpha=3/2
+ END_DOC
+ double precision :: x_tmp
+ x_tmp = alpha_mu_gauss * mu * x
+ j_gauss_simple = 0.25d0 * c_mu_gauss / (alpha_mu_gauss*alpha_mu_gauss*mu) * dexp(-x_tmp*x_tmp)
+
+end
+
+double precision function j_mu_deriv(x,mu)
+ implicit none
+ BEGIN_DOC
+! d/dx j_mu(mu,x) = d/dx 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2)
+!              = 0.5*(1 - erf(mu x))
+ END_DOC
+ include 'constants.include.F'
+ double precision, intent(in) :: x,mu
+ j_mu_deriv = 0.5d0 * (1.d0 - derf(mu*x))
+end
+
+double precision function j_mu_deriv_2(x,mu)
+ implicit none
+ BEGIN_DOC
+! d^2/dx^2 j_mu(mu,x) = d^2/dx^2 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2)
+!                  = -mu/sqrt(pi) * exp(-(mu x)^2)
+ END_DOC
+ include 'constants.include.F'
+ double precision, intent(in) :: x,mu
+ j_mu_deriv_2 = - mu * inv_sq_pi * dexp(-x*mu*x*mu)
+end
+
+double precision function j_gauss_deriv(x,mu)
+ implicit none
+ include 'constants.include.F'
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+! d/dx j_gauss(mu,x) = d/dx c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2)
+!      with c = 27/(8 sqrt(pi)), alpha=3/2
+!                    = -0.5 * mu *  c * x * exp(-(alpha * mu * x)^2)
+ END_DOC
+ double precision :: x_tmp
+ x_tmp = alpha_mu_gauss * mu * x
+ j_gauss_deriv = -0.5d0 * mu * c_mu_gauss * x * exp(-x_tmp*x_tmp)
+end
+
+double precision function j_gauss_deriv_2(x,mu)
+ implicit none
+ include 'constants.include.F'
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+! d/dx j_gauss(mu,x) = d/dx c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2)
+!      with c = 27/(8 sqrt(pi)), alpha=3/2
+!                    = 0.5 * mu *  c * exp(-(alpha * mu * x)^2) * (2 (alpha*mu*x)^2 - 1)
+ END_DOC
+ double precision :: x_tmp
+ x_tmp = alpha_mu_gauss * mu * x
+ x_tmp = x_tmp * x_tmp
+ j_gauss_deriv_2 = 0.5d0 * mu * c_mu_gauss * exp(-x_tmp) * (2.d0*x_tmp - 1.d0)
+end
+
+double precision function j_erf_gauss_deriv(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+! d/dx (j_gauss(mu,x)+j_mu(mu,x)) 
+ END_DOC
+ double precision :: j_gauss_deriv,j_mu_deriv
+ j_erf_gauss_deriv = j_gauss_deriv(x,mu)+j_mu_deriv(x,mu)
+end
+
+double precision function j_erf_gauss_deriv_2(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+! d^2/dx^2 (j_gauss(mu,x)+j_mu(mu,x)) 
+ END_DOC
+ double precision :: j_gauss_deriv_2,j_mu_deriv_2
+ j_erf_gauss_deriv_2 = j_gauss_deriv_2(x,mu)+j_mu_deriv_2(x,mu)
+end
+
+
+double precision function pot_j_gauss(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+ ! effective scalar potential associated with the erf_gauss correlation factor
+ ! 
+ ! 1/x( 1 - 2 * d/dx j_erf_gauss(x,mu)) - d^2/dx^2 j_erf_gauss(x,mu)) - d/dx d/dx (j_erf_gauss(x,mu))^2
+ END_DOC
+ double precision :: j_erf_gauss_deriv_2,j_erf_gauss_deriv
+ double precision :: deriv_1, deriv_2
+ pot_j_gauss = 0.d0
+ if(x.ne.0.d0)then
+  deriv_1 = j_erf_gauss_deriv(x,mu)
+  deriv_2 = j_erf_gauss_deriv_2(x,mu)
+  pot_j_gauss = 1.d0/x * (1.d0 - 2.d0 * deriv_1) - deriv_1 * deriv_1 - deriv_2
+ endif
+
+end
+
+double precision function pot_j_mu(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+ ! effective scalar potential associated with the correlation factor
+ ! 
+ ! 1/x( 1 - 2 * d/dx j_erf(x,mu)) - d^2/dx^2 j_erf(x,mu)) - d/dx d/dx (j_erf(x,mu))^2
+ END_DOC
+ double precision :: j_mu_deriv_2,j_mu_deriv
+ double precision :: deriv_1, deriv_2
+ pot_j_mu = 0.d0
+ if(x.ne.0.d0)then
+  deriv_1 = j_mu_deriv(x,mu)
+  deriv_2 = j_mu_deriv_2(x,mu)
+  pot_j_mu= 1.d0/x * (1.d0 - 2.d0 * deriv_1) - deriv_1 * deriv_1 - deriv_2
+ endif
+
+end

plugins/local/non_h_ints_mu/print_jastrow_psi.irp.f

@@ -0,0 +1,15 @@
+program print_j_psi
+ implicit none
+ integer :: i,j,a,b
+ do i = 1, elec_beta_num ! r2
+  do j = 1, elec_alpha_num ! r1
+   do a = elec_beta_num+1, mo_num ! r2
+    do b = elec_alpha_num+1, mo_num ! r1 
+     print*,b,a,j,i
+     print*,c_ij_ab_jastrow(b,a,j,i)
+    enddo
+   enddo
+  enddo
+ enddo
+
+end

plugins/local/non_h_ints_mu/test_mu_of_r_tc.irp.f

@@ -0,0 +1,157 @@
+program test_mu_of_r_tc
+  implicit none
+  BEGIN_DOC
+! TODO
+  END_DOC
+  ! You specify that you want to avoid any contribution from 
+  ! orbitals coming from core 
+ call test_grad_f_mean_field
+ call test_grad_mu_mf
+ call plot_mu_of_r_mf
+end
+
+
+subroutine test_grad_f_mean_field
+ implicit none
+ integer :: i_point,k
+ double precision :: weight,r(3)
+ double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3)
+ double precision :: grad_dm_a(3), grad_dm_b(3)
+ double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b
+
+ double precision :: num_grad_f_mf_ab(3), num_grad_two_bod_dens(3)
+ double precision :: num_grad_dm_a(3), num_grad_dm_b(3)
+ double precision :: f_mf_ab_p,f_mf_ab_m
+ double precision :: two_bod_dens_p, two_bod_dens_m
+ double precision :: dm_a_p, dm_a_m
+ double precision :: dm_b_p, dm_b_m
+ double precision :: rbis(3), dr
+ double precision :: accu_grad_f_mf_ab(3),accu_grad_two_bod_dens(3)
+ double precision :: accu_grad_dm_a(3),accu_grad_dm_b(3)
+ double precision :: accu_f_mf_ab, accu_two_bod_dens, accu_dm_a, accu_dm_b
+ dr = 0.00001d0
+ accu_f_mf_ab = 0.d0 
+ accu_two_bod_dens = 0.d0 
+ accu_dm_a = 0.d0 
+ accu_dm_b = 0.d0
+
+ accu_grad_f_mf_ab = 0.d0
+ accu_grad_two_bod_dens = 0.d0
+ accu_grad_dm_a = 0.d0
+ accu_grad_dm_b = 0.d0
+ do i_point = 1, n_points_final_grid
+  r(1:3)   = final_grid_points(1:3,i_point)
+  weight = final_weight_at_r_vector(i_point)
+  call get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b)
+  call get_f_mf_ab(r,f_mf_ab_p,two_bod_dens_p, dm_a_p, dm_b_p)
+  accu_f_mf_ab += weight * dabs(f_mf_ab - f_mf_ab_p)
+  accu_two_bod_dens += weight * dabs(two_bod_dens - two_bod_dens_p)
+  accu_dm_a += weight*dabs(dm_a - dm_a_p)
+  accu_dm_b += weight*dabs(dm_b - dm_b_p)
+  do k = 1, 3
+   rbis = r
+   rbis(k) += dr
+   call get_f_mf_ab(rbis,f_mf_ab_p,two_bod_dens_p, dm_a_p, dm_b_p)
+   rbis = r
+   rbis(k) -= dr
+   call get_f_mf_ab(rbis,f_mf_ab_m,two_bod_dens_m, dm_a_m, dm_b_m)
+   num_grad_f_mf_ab(k) = (f_mf_ab_p - f_mf_ab_m)/(2.d0*dr)
+   num_grad_two_bod_dens(k) = (two_bod_dens_p - two_bod_dens_m)/(2.d0*dr)
+   num_grad_dm_a(k) = (dm_a_p - dm_a_m)/(2.d0*dr)
+   num_grad_dm_b(k) = (dm_b_p - dm_b_m)/(2.d0*dr)
+  enddo
+  do k = 1, 3
+   accu_grad_f_mf_ab(k) += weight * dabs(grad_f_mf_ab(k) - num_grad_f_mf_ab(k))
+   accu_grad_two_bod_dens(k) += weight * dabs(grad_two_bod_dens(k) - num_grad_two_bod_dens(k))
+   accu_grad_dm_a(k) += weight * dabs(grad_dm_a(k) - num_grad_dm_a(k))
+   accu_grad_dm_b(k) += weight * dabs(grad_dm_b(k) - num_grad_dm_b(k))
+  enddo
+ enddo
+ print*,'accu_f_mf_ab = ',accu_f_mf_ab
+ print*,'accu_two_bod_dens = ',accu_two_bod_dens
+ print*,'accu_dm_a = ',accu_dm_a
+ print*,'accu_dm_b = ',accu_dm_b
+ print*,'accu_grad_f_mf_ab = '
+ print*,accu_grad_f_mf_ab
+ print*,'accu_grad_two_bod_dens = '
+ print*,accu_grad_two_bod_dens
+ print*,'accu_dm_a = '
+ print*,accu_grad_dm_a
+ print*,'accu_dm_b = '
+ print*,accu_grad_dm_b
+
+end
+
+subroutine test_grad_mu_mf
+ implicit none
+ integer :: i_point,k
+ double precision :: weight,r(3),rbis(3)
+ double precision :: mu_mf, dm,grad_mu_mf(3), grad_dm(3)
+ double precision :: mu_mf_p, mu_mf_m, dm_m, dm_p, num_grad_mu_mf(3),dr, num_grad_dm(3)
+ double precision :: accu_mu, accu_dm, accu_grad_dm(3), accu_grad_mu_mf(3)
+ dr = 0.00001d0
+ accu_grad_mu_mf = 0.d0
+ accu_mu = 0.d0
+ accu_grad_dm = 0.d0
+ accu_dm = 0.d0
+ do i_point = 1, n_points_final_grid
+  r(1:3)   = final_grid_points(1:3,i_point)
+  weight = final_weight_at_r_vector(i_point)
+  call grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm)
+  call mu_of_r_mean_field(r,mu_mf_p, dm_p)
+  accu_mu += weight*dabs(mu_mf_p - mu_mf)
+  accu_dm += weight*dabs(dm_p - dm)
+  do k = 1, 3
+   rbis = r
+   rbis(k) += dr
+   call mu_of_r_mean_field(rbis,mu_mf_p, dm_p)
+   rbis = r
+   rbis(k) -= dr
+   call mu_of_r_mean_field(rbis,mu_mf_m, dm_m)
+
+   num_grad_mu_mf(k) = (mu_mf_p - mu_mf_m)/(2.d0*dr)
+   num_grad_dm(k) = (dm_p - dm_m)/(2.d0*dr)
+  enddo
+  do k = 1, 3
+   accu_grad_dm(k)+= weight *dabs(num_grad_dm(k) - grad_dm(k))
+   accu_grad_mu_mf(k)+= weight *dabs(num_grad_mu_mf(k) - grad_mu_mf(k))
+  enddo
+ enddo
+ print*,'accu_mu = ',accu_mu
+ print*,'accu_dm = ',accu_dm
+ print*,'accu_grad_dm = '
+ print*, accu_grad_dm
+ print*,'accu_grad_mu_mf = '
+ print*, accu_grad_mu_mf
+
+end
+
+subroutine plot_mu_of_r_mf
+ implicit none
+  include 'constants.include.F'
+ integer :: ipoint,npoint
+ double precision :: dx,r(3),xmax,xmin
+ double precision :: accu_mu,accu_nelec,mu_mf, dm,mu_mf_tc
+ character*(128) :: output
+ integer :: i_unit_output,getUnitAndOpen
+ output=trim(ezfio_filename)//'.mu_mf'
+ i_unit_output = getUnitAndOpen(output,'w')
+ xmax = 5.D0
+ xmin = 0.d0
+ npoint = 10000
+ dx = (xmax - xmin)/dble(npoint)
+ r = 0.d0
+ r(1) = xmin
+ accu_mu = 0.d0
+ accu_nelec = 0.d0
+ do ipoint = 1, npoint
+  call mu_of_r_mean_field(r,mu_mf, dm)
+  call mu_of_r_mean_field_tc(r,mu_mf_tc, dm)
+  write(i_unit_output,'(100(F16.10,X))')r(1),mu_mf,mu_mf_tc,dm
+  accu_mu    += mu_mf * dm * r(1)**2*dx*4.D0*pi
+  accu_nelec +=         dm * r(1)**2*dx*4.D0*pi
+  r(1) += dx
+ enddo
+ print*,'nelec      = ',accu_nelec
+ print*,'mu average = ',accu_mu/accu_nelec
+end

plugins/local/tc_int/uninstall

@@ -9,5 +9,5 @@ then
 fi
 
 rm -rf ${PWD}/CuTC
-rm ${QP_ROOT}/lib/libcutcint.so
+rm -f ${QP_ROOT}/lib/libcutcint.so
 

plugins/local/tc_keywords/EZFIO.cfg

@@ -230,7 +230,7 @@ default: 70
 type: character*(32)
 doc: approach used to evaluate TC integrals [ analytic | numeric | semi-analytic ]
 interface: ezfio,ocaml,provider
-default: semi-analytic
+default: numeric
 
 [minimize_lr_angles]
 type: logical

scripts/import_champ_jastrow.py

@@ -45,7 +45,7 @@ if __name__ == '__main__':
   jastrow_file = sys.argv[2]
   jastrow = import_jastrow(jastrow_file)
   print (jastrow)
-  ezfio.set_jastrow_jast_type("Qmckl")
+  ezfio.set_jastrow_j2e_type("Qmckl")
   ezfio.set_jastrow_jast_qmckl_type_nucl_num(jastrow['type_num'])
   charges = ezfio.get_nuclei_nucl_charge()
   types = {}

scripts/qp_import_trexio.py

@@ -261,13 +261,10 @@ def write_ezfio(trexio_filename, filename):
     except:
         cartesian = True
 
-    if not cartesian:
-        raise TypeError('Only cartesian TREXIO files can be converted')
-
     ao_num = trexio.read_ao_num(trexio_file)
     ezfio.set_ao_basis_ao_num(ao_num)
 
-    if shell_num > 0:
+    if cartesian and shell_num > 0:
         ao_shell    = trexio.read_ao_shell(trexio_file)
         at = [ nucl_index[i]+1 for i in ao_shell ]
         ezfio.set_ao_basis_ao_nucl(at)
@@ -330,7 +327,7 @@ def write_ezfio(trexio_filename, filename):
         print("OK")
 
     else:
-        print("None")
+        print("None: integrals should be also imported using qp run import_trexio_integrals")
 
 
     #                _

src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f

@@ -111,8 +111,9 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a
   complex*16                   :: accu, P_center(3)
   complex*16                   :: d(0:n_pt_in)
 
-  complex*16                   :: V_n_e_cosgtos
+  complex*16, external         :: V_n_e_cosgtos
-  complex*16                   :: crint
+  complex*16, external         :: crint_2
+  complex*16, external         :: crint_sum_2
 
   if ( (A_center(1)/=B_center(1)) .or. (A_center(2)/=B_center(2)) .or. (A_center(3)/=B_center(3)) .or. &
        (A_center(1)/=C_center(1)) .or. (A_center(2)/=C_center(2)) .or. (A_center(3)/=C_center(3)) ) then
@@ -158,27 +159,27 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a
 
   n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) )
   if(n_pt == 0) then
-    NAI_pol_mult_cosgtos = coeff * crint(0, const)
+    NAI_pol_mult_cosgtos = coeff * crint_2(0, const)
     return
   endif
 
-  call give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
+  call give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, &
-                                         , power_A, power_B, C_center, n_pt_in, d, n_pt_out)
+                                          power_A, power_B, C_center, n_pt_in, d, n_pt_out)
 
   if(n_pt_out < 0) then
     NAI_pol_mult_cosgtos = (0.d0, 0.d0)
     return
   endif
 
-  accu = (0.d0, 0.d0)
+  !accu = (0.d0, 0.d0)
-  do i = 0, n_pt_out, 2
+  !do i = 0, n_pt_out, 2
-    accu += crint(shiftr(i, 1), const) * d(i)
+  !  accu += crint_2(shiftr(i, 1), const) * d(i)
-
+  !enddo
-!    print *, shiftr(i, 1), real(const), real(d(i)), real(crint(shiftr(i, 1), const))
+  accu = crint_sum_2(n_pt_out, const, d)
-  enddo
   NAI_pol_mult_cosgtos = accu * coeff
 
-end function NAI_pol_mult_cosgtos
+  return
+end
 
 ! ---
 
@@ -312,7 +313,7 @@ subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
     d(i) = d1(i)
   enddo
 
-end subroutine give_cpolynomial_mult_center_one_e
+end
 
 ! ---
 
@@ -405,7 +406,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_
 
   endif
 
-end subroutine I_x1_pol_mult_one_e_cosgtos
+end
 
 ! ---
 
@@ -467,7 +468,7 @@ recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim)
 
   endif
 
-end subroutine I_x2_pol_mult_one_e_cosgtos
+end
 
 ! ---
 
@@ -502,7 +503,7 @@ complex*16 function V_n_e_cosgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta)
                   * V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1)
   endif
 
-end function V_n_e_cosgtos
+end
 
 ! ---
 
@@ -529,7 +530,7 @@ complex*16 function V_r_cosgtos(n, alpha)
     V_r_cosgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1)
   endif
 
-end function V_r_cosgtos
+end
 
 ! ---
 

src/ao_two_e_ints/cholesky.irp.f

@@ -73,14 +73,13 @@ END_PROVIDER
    integer, external              :: getUnitAndOpen
    integer                        :: iunit, ierr
 
-   ndim8 = ao_num*ao_num*1_8
+   ndim8 = ao_num*ao_num*1_8+1
    double precision :: wall0,wall1
 
-   type(c_ptr)                    :: c_pointer(2)
+   type(mmap_type) :: map
-   integer                        :: fd(2)
 
    PROVIDE nproc ao_cholesky_threshold do_direct_integrals qp_max_mem
-   PROVIDE nucl_coord ao_two_e_integral_schwartz
+   PROVIDE nucl_coord
    call set_multiple_levels_omp(.False.)
 
    call wall_time(wall0)
@@ -143,19 +142,21 @@ END_PROVIDER
 
      if (do_direct_integrals) then
        !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21)
-       do i8=ndim8,1,-1
+       do i8=ndim8-1,1,-1
          D(i8) = ao_two_e_integral(addr1(i8), addr2(i8),              &
              addr1(i8), addr2(i8))
        enddo
        !$OMP END PARALLEL DO
      else
        !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i8) SCHEDULE(dynamic,21)
-       do i8=ndim8,1,-1
+       do i8=ndim8-1,1,-1
          D(i8) = get_ao_two_e_integral(addr1(i8), addr1(i8),          &
              addr2(i8), addr2(i8), ao_integrals_map)
        enddo
        !$OMP END PARALLEL DO
      endif
+     ! Just to guarentee termination
+     D(ndim8) = 0.d0
 
      D_sorted(:) = -D(:)
      call dsort_noidx_big(D_sorted,ndim8)
@@ -179,14 +180,9 @@ END_PROVIDER
      if (elec_num > 10) then
        rank_max = min(np,20*elec_num*elec_num)
      endif
-     call mmap(trim(ezfio_work_dir)//'cholesky_ao_tmp', (/ ndim8, rank_max /), 8, fd(1), .False., .True., c_pointer(1))
-     call c_f_pointer(c_pointer(1), L, (/ ndim8, rank_max /))
-
-     ! Deleting the file while it is open makes the file invisible on the filesystem,
-     ! and automatically deleted, even if the program crashes
-     iunit = getUnitAndOpen(trim(ezfio_work_dir)//'cholesky_ao_tmp', 'R')
-     close(iunit,status='delete')
 
+     call mmap_create_d('', (/ ndim8, rank_max /), .False., .True., map)
+     L => map%d2
 
      ! 3.
      N = 0
@@ -205,6 +201,7 @@ END_PROVIDER
        i = i+1
 
 
+
        block_size = max(N,24)
 
        ! Determine nq so that Delta fits in memory
@@ -314,9 +311,10 @@ END_PROVIDER
        ! g.
 
        iblock = 0
+
        do j=1,nq
 
-         if ( (Qmax <= Dmin).or.(N+j*1_8 > ndim8) ) exit
+         if ( (Qmax < Dmin).or.(N+j*1_8 > ndim8) ) exit
 
          ! i.
          rank = N+j
@@ -476,7 +474,7 @@ END_PROVIDER
      enddo
      !$OMP END PARALLEL DO
 
-     call munmap( (/ ndim8, rank_max /), 8, fd(1), c_pointer(1) )
+     call mmap_destroy(map)
 
      cholesky_ao_num = rank
 

src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f

@@ -35,11 +35,9 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
 
   if(ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024) then
 
-    !print *, ' with shwartz acc '
     ao_two_e_integral_cosgtos = ao_2e_cosgtos_schwartz_accel(i, j, k, l)
 
   else
-    !print *, ' without shwartz acc '
 
     dim1 = n_pt_max_integrals
 
@@ -51,7 +49,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
     ao_two_e_integral_cosgtos = 0.d0
 
     if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then
-      !print *, ' not the same center'
 
       do p = 1, 3
         I_power(p)  = ao_power(i,p)
@@ -72,72 +69,22 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
           coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j)
           expo2 = ao_expo_ord_transp_cosgtos(q,j) 
 
-          call give_explicit_cpoly_and_cgaussian( P1_new, P1_center, pp1, fact_p1, iorder_p1               &
+          call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, &
-                                                , expo1, expo2, I_power, J_power, I_center, J_center, dim1 )
+                                                 expo1, expo2, I_power, J_power, I_center, J_center, dim1)
           p1_inv = (1.d0,0.d0) / pp1
 
-          call give_explicit_cpoly_and_cgaussian( P2_new, P2_center, pp2, fact_p2, iorder_p2                      &
+          call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, &
-                                                , conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1 )
+                                                 conjg(expo1), expo2, I_power, J_power, I_center, J_center, dim1)
           p2_inv = (1.d0,0.d0) / pp2
 
-          call give_explicit_cpoly_and_cgaussian( P3_new, P3_center, pp3, fact_p3, iorder_p3                      &
+          call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, &
-                                                , expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1 )
+                                                 expo1, conjg(expo2), I_power, J_power, I_center, J_center, dim1)
           p3_inv = (1.d0,0.d0) / pp3
 
-          call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4                             &
+          call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, &
-                                                , conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1 )
+                                                 conjg(expo1), conjg(expo2), I_power, J_power, I_center, J_center, dim1)
           p4_inv = (1.d0,0.d0) / pp4
 
-          !integer :: ii
-          !do ii = 1, 3
-          !  print *, 'fact_p1', fact_p1
-          !  print *, 'fact_p2', fact_p2
-          !  print *, 'fact_p3', fact_p3
-          !  print *, 'fact_p4', fact_p4
-          !  !print *, pp1, p1_inv
-          !  !print *, pp2, p2_inv
-          !  !print *, pp3, p3_inv
-          !  !print *, pp4, p4_inv
-          !enddo 
-          !  if( abs(aimag(P1_center(ii))) .gt. 0.d0 ) then
-          !    print *, ' P_1 is complex !!'
-          !    print *, P1_center
-          !    print *, expo1, expo2
-          !    print *, conjg(expo1), conjg(expo2)
-          !    stop
-          !  endif
-          !  if( abs(aimag(P2_center(ii))) .gt. 0.d0 ) then
-          !    print *, ' P_2 is complex !!'
-          !    print *, P2_center
-          !    print *, ' old expos:'
-          !    print *, expo1, expo2
-          !    print *, conjg(expo1), conjg(expo2)
-          !    print *, ' new expo:'
-          !    print *, pp2, p2_inv
-          !    print *, ' factor:'
-          !    print *, fact_p2
-          !    print *, ' old centers:'
-          !    print *, I_center, J_center
-          !    print *, ' powers:'
-          !    print *, I_power, J_power
-          !    stop
-          !  endif
-          !  if( abs(aimag(P3_center(ii))) .gt. 0.d0 ) then
-          !    print *, ' P_3 is complex !!'
-          !    print *, P3_center
-          !    print *, expo1, expo2
-          !    print *, conjg(expo1), conjg(expo2)
-          !    stop
-          !  endif
-          !  if( abs(aimag(P4_center(ii))) .gt. 0.d0 ) then
-          !    print *, ' P_4 is complex !!'
-          !    print *, P4_center
-          !    print *, expo1, expo2
-          !    print *, conjg(expo1), conjg(expo2)
-          !    stop
-          !  endif
-          !enddo
-
           do r = 1, ao_prim_num(k)
             coef3 = coef2 * ao_coef_norm_ord_transp_cosgtos(r,k)
             expo3 = ao_expo_ord_transp_cosgtos(r,k) 
@@ -146,66 +93,40 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
               coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l)
               expo4 = ao_expo_ord_transp_cosgtos(s,l) 
 
-              call give_explicit_cpoly_and_cgaussian( Q1_new, Q1_center, qq1, fact_q1, iorder_q1               &
+              call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, &
-                                                    , expo3, expo4, K_power, L_power, K_center, L_center, dim1 )
+                                                     expo3, expo4, K_power, L_power, K_center, L_center, dim1)
               q1_inv = (1.d0,0.d0) / qq1
 
-              call give_explicit_cpoly_and_cgaussian( Q2_new, Q2_center, qq2, fact_q2, iorder_q2               &
+              call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, &
-                                                    , conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1 )
+                                                     conjg(expo3), expo4, K_power, L_power, K_center, L_center, dim1)
               q2_inv = (1.d0,0.d0) / qq2
 
-              !do ii = 1, 3
+              integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
-              !  !print *, qq1, q1_inv
+                                                                   Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
-              !  !print *, qq2, q2_inv
-              !  print *, 'fact_q1', fact_q1
-              !  print *, 'fact_q2', fact_q2
-              !enddo 
-              !  if( abs(aimag(Q1_center(ii))) .gt. 0.d0 ) then
-              !    print *, ' Q_1 is complex !!'
-              !    print *, Q1_center
-              !    print *, expo3, expo4
-              !    print *, conjg(expo3), conjg(expo4)
-              !    stop
-              !  endif
-              !  if( abs(aimag(Q2_center(ii))) .gt. 0.d0 ) then
-              !    print *, ' Q_2 is complex !!'
-              !    print *, Q2_center
-              !    print *, expo3, expo4
-              !    print *, conjg(expo3), conjg(expo4)
-              !    stop
-              !  endif
-              !enddo
 
+              integral2 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
+                                                                   Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
 
-              integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 &
+              integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
-                                                                  , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 )
+                                                                   Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
 
-              integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 &
+              integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
-                                                                  , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 )
+                                                                   Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
 
-              integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 &
+              integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
-                                                                  , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 )
+                                                                   Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
 
-              integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 &
+              integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
-                                                                  , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 )
+                                                                   Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
 
-              integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 &
+              integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
-                                                                  , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 )
+                                                                   Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1)
 
-              integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 &
+              integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
-                                                                  , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 )
+                                                                   Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2)
-
-              integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 &
-                                                                  , Q1_new, Q1_center, fact_q1, qq1, q1_inv, iorder_q1 )
-
-              integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 &
-                                                                  , Q2_new, Q2_center, fact_q2, qq2, q2_inv, iorder_q2 )
 
               integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
 
-              !integral_tot = integral1
-              !print*, integral_tot
-
               ao_two_e_integral_cosgtos = ao_two_e_integral_cosgtos + coef4 * 2.d0 * real(integral_tot)
             enddo ! s
           enddo ! r
@@ -213,7 +134,6 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
       enddo ! p
 
     else
-      !print *, ' the same center'
 
       do p = 1, 3
         I_power(p) = ao_power(i,p)
@@ -290,7 +210,7 @@ double precision function ao_two_e_integral_cosgtos(i, j, k, l)
     endif
   endif
 
-end function ao_two_e_integral_cosgtos
+end
 
 ! ---
 
@@ -326,8 +246,8 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
   double precision              :: thr
   double precision              :: schwartz_ij
 
-  complex*16                    :: ERI_cosgtos
+  complex*16, external          :: ERI_cosgtos
-  complex*16                    :: general_primitive_integral_cosgtos
+  complex*16, external          :: general_primitive_integral_cosgtos
 
   ao_2e_cosgtos_schwartz_accel = 0.d0
 
@@ -366,45 +286,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
         coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l)
         expo2 = ao_expo_ord_transp_cosgtos(s,l) 
 
-        call give_explicit_cpoly_and_cgaussian( P1_new, P1_center, pp1, fact_p1, iorder_p1               &
+        call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, &
-                                              , expo1, expo2, K_power, L_power, K_center, L_center, dim1 )
+                                               expo1, expo2, K_power, L_power, K_center, L_center, dim1)
         p1_inv = (1.d0,0.d0) / pp1
 
-        call give_explicit_cpoly_and_cgaussian( P2_new, P2_center, pp2, fact_p2, iorder_p2                      &
+        call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, &
-                                              , conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1 )
+                                               conjg(expo1), expo2, K_power, L_power, K_center, L_center, dim1)
         p2_inv = (1.d0,0.d0) / pp2
 
-        call give_explicit_cpoly_and_cgaussian( P3_new, P3_center, pp3, fact_p3, iorder_p3                      &
+        call give_explicit_cpoly_and_cgaussian(P3_new, P3_center, pp3, fact_p3, iorder_p3, &
-                                              , expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1 )
+                                               expo1, conjg(expo2), K_power, L_power, K_center, L_center, dim1)
         p3_inv = (1.d0,0.d0) / pp3
 
-        call give_explicit_cpoly_and_cgaussian( P4_new, P4_center, pp4, fact_p4, iorder_p4                             &
+        call give_explicit_cpoly_and_cgaussian(P4_new, P4_center, pp4, fact_p4, iorder_p4, &
-                                              , conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1 )
+                                               conjg(expo1), conjg(expo2), K_power, L_power, K_center, L_center, dim1)
         p4_inv = (1.d0,0.d0) / pp4
 
-        integral1 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 &
+        integral1 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
-                                                            , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 )
+                                                             P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
 
-        integral2 = general_primitive_integral_cosgtos( dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 &
+        integral2 = general_primitive_integral_cosgtos(dim1, P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1, &
-                                                            , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 )
+                                                             P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
 
-        integral3 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 &
+        integral3 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
-                                                            , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 )
+                                                             P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
 
-        integral4 = general_primitive_integral_cosgtos( dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 &
+        integral4 = general_primitive_integral_cosgtos(dim1, P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2, &
-                                                            , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 )
+                                                             P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
 
-        integral5 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 &
+        integral5 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
-                                                            , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 )
+                                                             P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
 
-        integral6 = general_primitive_integral_cosgtos( dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3 &
+        integral6 = general_primitive_integral_cosgtos(dim1, P3_new, P3_center, fact_p3, pp3, p3_inv, iorder_p3, &
-                                                            , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 )
+                                                             P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
 
-        integral7 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 &
+        integral7 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
-                                                            , P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1 )
+                                                             P1_new, P1_center, fact_p1, pp1, p1_inv, iorder_p1)
 
-        integral8 = general_primitive_integral_cosgtos( dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4 &
+        integral8 = general_primitive_integral_cosgtos(dim1, P4_new, P4_center, fact_p4, pp4, p4_inv, iorder_p4, &
-                                                            , P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2 )
+                                                             P2_new, P2_center, fact_p2, pp2, p2_inv, iorder_p2)
 
         integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
 
@@ -544,41 +464,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
         coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(s,l) * ao_coef_norm_ord_transp_cosgtos(s,l)
         expo2 = ao_expo_ord_transp_cosgtos(s,l)
 
-        integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2                     &
+        integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2,                     &
-                               , K_power(1), L_power(1), K_power(1), L_power(1) &
+                                K_power(1), L_power(1), K_power(1), L_power(1), &
-                               , K_power(2), L_power(2), K_power(2), L_power(2) &
+                                K_power(2), L_power(2), K_power(2), L_power(2), &
-                               , K_power(3), L_power(3), K_power(3), L_power(3) )
+                                K_power(3), L_power(3), K_power(3), L_power(3))
-        integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2              &
-                               , K_power(1), L_power(1), K_power(1), L_power(1) &
-                               , K_power(2), L_power(2), K_power(2), L_power(2) &
-                               , K_power(3), L_power(3), K_power(3), L_power(3) )
 
-        integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2              &
+        integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2,              &
-                               , K_power(1), L_power(1), K_power(1), L_power(1) &
+                                K_power(1), L_power(1), K_power(1), L_power(1), &
-                               , K_power(2), L_power(2), K_power(2), L_power(2) &
+                                K_power(2), L_power(2), K_power(2), L_power(2), &
-                               , K_power(3), L_power(3), K_power(3), L_power(3) )
+                                K_power(3), L_power(3), K_power(3), L_power(3))
-        integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2       &
-                               , K_power(1), L_power(1), K_power(1), L_power(1) &
-                               , K_power(2), L_power(2), K_power(2), L_power(2) &
-                               , K_power(3), L_power(3), K_power(3), L_power(3) )
 
-        integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2              &
+        integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2,              &
-                               , K_power(1), L_power(1), K_power(1), L_power(1) &
+                                K_power(1), L_power(1), K_power(1), L_power(1), &
-                               , K_power(2), L_power(2), K_power(2), L_power(2) &
+                                K_power(2), L_power(2), K_power(2), L_power(2), &
-                               , K_power(3), L_power(3), K_power(3), L_power(3) )
+                                K_power(3), L_power(3), K_power(3), L_power(3))
-        integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2       &
-                               , K_power(1), L_power(1), K_power(1), L_power(1) &
-                               , K_power(2), L_power(2), K_power(2), L_power(2) &
-                               , K_power(3), L_power(3), K_power(3), L_power(3) )
 
-        integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2       &
+        integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2,       &
-                               , K_power(1), L_power(1), K_power(1), L_power(1) &
+                                K_power(1), L_power(1), K_power(1), L_power(1), &
-                               , K_power(2), L_power(2), K_power(2), L_power(2) &
+                                K_power(2), L_power(2), K_power(2), L_power(2), &
-                               , K_power(3), L_power(3), K_power(3), L_power(3) )
+                                K_power(3), L_power(3), K_power(3), L_power(3))
-        integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 &
+
-                               , K_power(1), L_power(1), K_power(1), L_power(1)  &
+        integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2,              &
-                               , K_power(2), L_power(2), K_power(2), L_power(2)  &
+                                K_power(1), L_power(1), K_power(1), L_power(1), &
-                               , K_power(3), L_power(3), K_power(3), L_power(3)  )
+                                K_power(2), L_power(2), K_power(2), L_power(2), &
+                                K_power(3), L_power(3), K_power(3), L_power(3))
+
+        integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2,       &
+                                K_power(1), L_power(1), K_power(1), L_power(1), &
+                                K_power(2), L_power(2), K_power(2), L_power(2), &
+                                K_power(3), L_power(3), K_power(3), L_power(3))
+
+        integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2,       &
+                                K_power(1), L_power(1), K_power(1), L_power(1), &
+                                K_power(2), L_power(2), K_power(2), L_power(2), &
+                                K_power(3), L_power(3), K_power(3), L_power(3))
+
+        integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, &
+                                K_power(1), L_power(1), K_power(1), L_power(1),  &
+                                K_power(2), L_power(2), K_power(2), L_power(2),  &
+                                K_power(3), L_power(3), K_power(3), L_power(3))
 
         integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
 
@@ -598,45 +522,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
         coef2 = coef1 * ao_coef_norm_ord_transp_cosgtos(q,j)
         expo2 = ao_expo_ord_transp_cosgtos(q,j)
 
-        integral1 = ERI_cosgtos( expo1, expo2, expo1, expo2                     &
+        integral1 = ERI_cosgtos(expo1, expo2, expo1, expo2,                     &
-                               , I_power(1), J_power(1), I_power(1), J_power(1) &
+                                I_power(1), J_power(1), I_power(1), J_power(1), &
-                               , I_power(2), J_power(2), I_power(2), J_power(2) &
+                                I_power(2), J_power(2), I_power(2), J_power(2), &
-                               , I_power(3), J_power(3), I_power(3), J_power(3) )
+                                I_power(3), J_power(3), I_power(3), J_power(3))
 
-        integral2 = ERI_cosgtos( expo1, expo2, conjg(expo1), expo2              &
+        integral2 = ERI_cosgtos(expo1, expo2, conjg(expo1), expo2,              &
-                               , I_power(1), J_power(1), I_power(1), J_power(1) &
+                                I_power(1), J_power(1), I_power(1), J_power(1), &
-                               , I_power(2), J_power(2), I_power(2), J_power(2) &
+                                I_power(2), J_power(2), I_power(2), J_power(2), &
-                               , I_power(3), J_power(3), I_power(3), J_power(3) )
+                                I_power(3), J_power(3), I_power(3), J_power(3))
 
-        integral3 = ERI_cosgtos( conjg(expo1), expo2, expo1, expo2              &
+        integral3 = ERI_cosgtos(conjg(expo1), expo2, expo1, expo2,              &
-                               , I_power(1), J_power(1), I_power(1), J_power(1) &
+                                I_power(1), J_power(1), I_power(1), J_power(1), &
-                               , I_power(2), J_power(2), I_power(2), J_power(2) &
+                                I_power(2), J_power(2), I_power(2), J_power(2), &
-                               , I_power(3), J_power(3), I_power(3), J_power(3) )
+                                I_power(3), J_power(3), I_power(3), J_power(3))
 
-        integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo1), expo2       &
+        integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo1), expo2,       &
-                               , I_power(1), J_power(1), I_power(1), J_power(1) &
+                                I_power(1), J_power(1), I_power(1), J_power(1), &
-                               , I_power(2), J_power(2), I_power(2), J_power(2) &
+                                I_power(2), J_power(2), I_power(2), J_power(2), &
-                               , I_power(3), J_power(3), I_power(3), J_power(3) )
+                                I_power(3), J_power(3), I_power(3), J_power(3))
 
-        integral5 = ERI_cosgtos( expo1, conjg(expo2), expo1, expo2              &
+        integral5 = ERI_cosgtos(expo1, conjg(expo2), expo1, expo2,              &
-                               , I_power(1), J_power(1), I_power(1), J_power(1) &
+                                I_power(1), J_power(1), I_power(1), J_power(1), &
-                               , I_power(2), J_power(2), I_power(2), J_power(2) &
+                                I_power(2), J_power(2), I_power(2), J_power(2), &
-                               , I_power(3), J_power(3), I_power(3), J_power(3) )
+                                I_power(3), J_power(3), I_power(3), J_power(3))
 
-        integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo1), expo2       &
+        integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo1), expo2,       &
-                               , I_power(1), J_power(1), I_power(1), J_power(1) &
+                                I_power(1), J_power(1), I_power(1), J_power(1), &
-                               , I_power(2), J_power(2), I_power(2), J_power(2) &
+                                I_power(2), J_power(2), I_power(2), J_power(2), &
-                               , I_power(3), J_power(3), I_power(3), J_power(3) )
+                                I_power(3), J_power(3), I_power(3), J_power(3))
 
-        integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo1, expo2       &
+        integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo1, expo2,       &
-                               , I_power(1), J_power(1), I_power(1), J_power(1) &
+                                I_power(1), J_power(1), I_power(1), J_power(1), &
-                               , I_power(2), J_power(2), I_power(2), J_power(2) &
+                                I_power(2), J_power(2), I_power(2), J_power(2), &
-                               , I_power(3), J_power(3), I_power(3), J_power(3) )
+                                I_power(3), J_power(3), I_power(3), J_power(3))
 
-        integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo1), expo2 &
+        integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo1), expo2, &
-                               , I_power(1), J_power(1), I_power(1), J_power(1)  &
+                                I_power(1), J_power(1), I_power(1), J_power(1),  &
-                               , I_power(2), J_power(2), I_power(2), J_power(2)  &
+                                I_power(2), J_power(2), I_power(2), J_power(2),  &
-                               , I_power(3), J_power(3), I_power(3), J_power(3)  )
+                                I_power(3), J_power(3), I_power(3), J_power(3))
 
         integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
 
@@ -655,45 +579,45 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
             coef4 = coef3 * ao_coef_norm_ord_transp_cosgtos(s,l)
             expo4 = ao_expo_ord_transp_cosgtos(s,l)
 
-            integral1 = ERI_cosgtos( expo1, expo2, expo3, expo4                     &
+            integral1 = ERI_cosgtos(expo1, expo2, expo3, expo4,                     &
-                                   , I_power(1), J_power(1), K_power(1), L_power(1) &
+                                    I_power(1), J_power(1), K_power(1), L_power(1), &
-                                   , I_power(2), J_power(2), K_power(2), L_power(2) &
+                                    I_power(2), J_power(2), K_power(2), L_power(2), &
-                                   , I_power(3), J_power(3), K_power(3), L_power(3) )
+                                    I_power(3), J_power(3), K_power(3), L_power(3))
 
-            integral2 = ERI_cosgtos( expo1, expo2, conjg(expo3), expo4              &
+            integral2 = ERI_cosgtos(expo1, expo2, conjg(expo3), expo4,              &
-                                   , I_power(1), J_power(1), K_power(1), L_power(1) &
+                                    I_power(1), J_power(1), K_power(1), L_power(1), &
-                                   , I_power(2), J_power(2), K_power(2), L_power(2) &
+                                    I_power(2), J_power(2), K_power(2), L_power(2), &
-                                   , I_power(3), J_power(3), K_power(3), L_power(3) )
+                                    I_power(3), J_power(3), K_power(3), L_power(3))
 
-            integral3 = ERI_cosgtos( conjg(expo1), expo2, expo3, expo4              &
+            integral3 = ERI_cosgtos(conjg(expo1), expo2, expo3, expo4,              &
-                                   , I_power(1), J_power(1), K_power(1), L_power(1) &
+                                    I_power(1), J_power(1), K_power(1), L_power(1), &
-                                   , I_power(2), J_power(2), K_power(2), L_power(2) &
+                                    I_power(2), J_power(2), K_power(2), L_power(2), &
-                                   , I_power(3), J_power(3), K_power(3), L_power(3) )
+                                    I_power(3), J_power(3), K_power(3), L_power(3))
 
-            integral4 = ERI_cosgtos( conjg(expo1), expo2, conjg(expo3), expo4       &
+            integral4 = ERI_cosgtos(conjg(expo1), expo2, conjg(expo3), expo4,       &
-                                   , I_power(1), J_power(1), K_power(1), L_power(1) &
+                                    I_power(1), J_power(1), K_power(1), L_power(1), &
-                                   , I_power(2), J_power(2), K_power(2), L_power(2) &
+                                    I_power(2), J_power(2), K_power(2), L_power(2), &
-                                   , I_power(3), J_power(3), K_power(3), L_power(3) )
+                                    I_power(3), J_power(3), K_power(3), L_power(3))
 
-            integral5 = ERI_cosgtos( expo1, conjg(expo2), expo3, expo4              &
+            integral5 = ERI_cosgtos(expo1, conjg(expo2), expo3, expo4,              &
-                                   , I_power(1), J_power(1), K_power(1), L_power(1) &
+                                    I_power(1), J_power(1), K_power(1), L_power(1), &
-                                   , I_power(2), J_power(2), K_power(2), L_power(2) &
+                                    I_power(2), J_power(2), K_power(2), L_power(2), &
-                                   , I_power(3), J_power(3), K_power(3), L_power(3) )
+                                    I_power(3), J_power(3), K_power(3), L_power(3))
 
-            integral6 = ERI_cosgtos( expo1, conjg(expo2), conjg(expo3), expo4       &
+            integral6 = ERI_cosgtos(expo1, conjg(expo2), conjg(expo3), expo4,       &
-                                   , I_power(1), J_power(1), K_power(1), L_power(1) &
+                                    I_power(1), J_power(1), K_power(1), L_power(1), &
-                                   , I_power(2), J_power(2), K_power(2), L_power(2) &
+                                    I_power(2), J_power(2), K_power(2), L_power(2), &
-                                   , I_power(3), J_power(3), K_power(3), L_power(3) )
+                                    I_power(3), J_power(3), K_power(3), L_power(3))
 
-            integral7 = ERI_cosgtos( conjg(expo1), conjg(expo2), expo3, expo4       &
+            integral7 = ERI_cosgtos(conjg(expo1), conjg(expo2), expo3, expo4,       &
-                                   , I_power(1), J_power(1), K_power(1), L_power(1) &
+                                    I_power(1), J_power(1), K_power(1), L_power(1), &
-                                   , I_power(2), J_power(2), K_power(2), L_power(2) &
+                                    I_power(2), J_power(2), K_power(2), L_power(2), &
-                                   , I_power(3), J_power(3), K_power(3), L_power(3) )
+                                    I_power(3), J_power(3), K_power(3), L_power(3))
 
-            integral8 = ERI_cosgtos( conjg(expo1), conjg(expo2), conjg(expo3), expo4 &
+            integral8 = ERI_cosgtos(conjg(expo1), conjg(expo2), conjg(expo3), expo4, &
-                                   , I_power(1), J_power(1), K_power(1), L_power(1)  &
+                                    I_power(1), J_power(1), K_power(1), L_power(1),  &
-                                   , I_power(2), J_power(2), K_power(2), L_power(2)  &
+                                    I_power(2), J_power(2), K_power(2), L_power(2),  &
-                                   , I_power(3), J_power(3), K_power(3), L_power(3)  )
+                                    I_power(3), J_power(3), K_power(3), L_power(3))
 
             integral_tot = integral1 + integral2 + integral3 + integral4 + integral5 + integral6 + integral7 + integral8
 
@@ -707,7 +631,7 @@ double precision function ao_2e_cosgtos_schwartz_accel(i, j, k, l)
 
   deallocate(schwartz_kl)
 
-end function ao_2e_cosgtos_schwartz_accel
+end
 
 ! ---
 
@@ -739,8 +663,8 @@ END_PROVIDER
 
 ! ---
 
-complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fact_p, p, p_inv, iorder_p &
+complex*16 function general_primitive_integral_cosgtos(dim, P_new, P_center, fact_p, p, p_inv, iorder_p, &
-                                                           , Q_new, Q_center, fact_q, q, q_inv, iorder_q )
+                                                            Q_new, Q_center, fact_q, q, q_inv, iorder_q)
 
   BEGIN_DOC
   !
@@ -765,7 +689,7 @@ complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fa
   complex*16             :: dx(0:max_dim), Ix_pol(0:max_dim), dy(0:max_dim), Iy_pol(0:max_dim), dz(0:max_dim), Iz_pol(0:max_dim)
   complex*16             :: d1(0:max_dim), d_poly(0:max_dim)
 
-  complex*16             :: crint_sum
+  complex*16             :: crint_sum_2
 
 
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol
@@ -912,13 +836,11 @@ complex*16 function general_primitive_integral_cosgtos( dim, P_new, P_center, fa
   !DIR$ FORCEINLINE
   call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out)
 
-  accu = crint_sum(n_pt_out, const, d1)
+  accu = crint_sum_2(n_pt_out, const, d1)
-!  print *, n_pt_out, real(d1(0:n_pt_out))
-!  print *, real(accu)
 
   general_primitive_integral_cosgtos = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / sq_ppq
 
-end function general_primitive_integral_cosgtos
+end
 
 ! ---
 
@@ -994,7 +916,7 @@ complex*16 function ERI_cosgtos(alpha, beta, delta, gama, a_x, b_x, c_x, d_x, a_
 
   ERI_cosgtos = I_f * coeff
 
-end function ERI_cosgtos
+end
 
 ! ---
 
@@ -1076,7 +998,7 @@ subroutine integrale_new_cosgtos(I_f, a_x, b_x, c_x, d_x, a_y, b_y, c_y, d_y, a_
     I_f += gauleg_w(i, j) * t1(i)
   enddo
 
-end subroutine integrale_new_cosgtos
+end
 
 ! ---
 
@@ -1123,7 +1045,7 @@ recursive subroutine I_x1_new_cosgtos(a, c, B_10, B_01, B_00, res, n_pt)
 
   endif
 
-end subroutine I_x1_new_cosgtos
+end
 
 ! ---
 
@@ -1163,7 +1085,7 @@ recursive subroutine I_x2_new_cosgtos(c, B_10, B_01, B_00, res, n_pt)
 
   endif
 
-end subroutine I_x2_new_cosgtos
+end
 
 ! ---
 
@@ -1234,7 +1156,7 @@ subroutine give_cpolynom_mult_center_x( P_center, Q_center, a_x, d_x, p, q, n_pt
     return
   endif
 
-end subroutine give_cpolynom_mult_center_x
+end
 
 ! ---
 
@@ -1277,7 +1199,7 @@ subroutine I_x1_pol_mult_cosgtos(a, c, B_10, B_01, B_00, C_00, D_00, d, nd, n_pt
 
   endif
 
-end subroutine I_x1_pol_mult_cosgtos
+end
 
 ! ---
 
@@ -1370,7 +1292,7 @@ recursive subroutine I_x1_pol_mult_recurs_cosgtos(a, c, B_10, B_01, B_00, C_00,
   !DIR$ FORCEINLINE
   call multiply_cpoly(Y, ny, C_00, 2, d, nd)
 
-end subroutine I_x1_pol_mult_recurs_cosgtos
+end
 
 ! ---
 
@@ -1426,7 +1348,7 @@ recursive subroutine I_x1_pol_mult_a1_cosgtos(c,B_10,B_01,B_00,C_00,D_00,d,nd,n_
   !DIR$ FORCEINLINE
   call multiply_cpoly(Y, ny, C_00, 2, d, nd)
 
-end subroutine I_x1_pol_mult_a1_cosgtos
+end
 
 ! ---
 
@@ -1490,7 +1412,7 @@ recursive subroutine I_x1_pol_mult_a2_cosgtos(c, B_10, B_01, B_00, C_00, D_00, d
   !DIR$ FORCEINLINE
   call multiply_cpoly(Y, ny, C_00, 2, d, nd)
 
-end subroutine I_x1_pol_mult_a2_cosgtos
+end
 
 ! ---
 
@@ -1575,7 +1497,7 @@ recursive subroutine I_x2_pol_mult_cosgtos(c, B_10, B_01, B_00, C_00, D_00, d, n
 
   end select
 
-end subroutine I_x2_pol_mult_cosgtos
+end
 
 ! ---
 

src/ao_two_e_ints/two_e_integrals.irp.f

@@ -1,6 +1,22 @@
 
 ! ---
 
+logical function do_schwartz_accel(i,j,k,l)
+ implicit none
+ BEGIN_DOC
+ ! If true, use Schwatrz to accelerate direct integral calculation
+ END_DOC
+  integer, intent(in)            :: i, j, k, l
+  if (do_ao_cholesky) then
+    do_schwartz_accel = .False.
+  else
+    do_schwartz_accel = (ao_prim_num(i) * ao_prim_num(j) * &
+                         ao_prim_num(k) * ao_prim_num(l) > 1024 )
+  endif
+
+end function
+
+
 double precision function ao_two_e_integral(i, j, k, l)
 
   BEGIN_DOC
@@ -25,6 +41,7 @@ double precision function ao_two_e_integral(i, j, k, l)
   double precision, external     :: ao_two_e_integral_cosgtos
   double precision, external     :: ao_two_e_integral_schwartz_accel
 
+  logical, external :: do_schwartz_accel
 
   if(use_cosgtos) then
     !print *, ' use_cosgtos for ao_two_e_integral ?', use_cosgtos
@@ -35,7 +52,7 @@ double precision function ao_two_e_integral(i, j, k, l)
 
     ao_two_e_integral = ao_two_e_integral_erf(i, j, k, l)
 
-  else if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
+  else if (do_schwartz_accel(i,j,k,l)) then
 
        ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
 

src/ccsd/ccsd_space_orb_sub.irp.f

@@ -154,14 +154,14 @@ subroutine run_ccsd_space_orb
 
     allocate(all_err(nO*nV+nO*nO*nV*(nV*1_8),cc_diis_depth), all_t(nO*nV+nO*nO*nV*(nV*1_8),cc_diis_depth))
     !$OMP PARALLEL PRIVATE(i,j) DEFAULT(SHARED)
+    !$OMP DO COLLAPSE(2)
     do j=1,cc_diis_depth
-      !$OMP DO
       do i=1, size(all_err,1)
         all_err(i,j) = 0d0
         all_t(i,j)   = 0d0
       enddo
-      !$OMP END DO NOWAIT
     enddo
+    !$OMP END DO NOWAIT
     !$OMP END PARALLEL
   endif
 
@@ -237,6 +237,7 @@ subroutine run_ccsd_space_orb
       call update_t2(nO,nV,cc_space_f_o,cc_space_f_v,r2%f,t2%f)
     else
       print*,'Unkown cc_method_method: '//cc_update_method
+      call abort
     endif
 
     call update_tau_space(nO,nV,t1%f,t1,t2,tau)

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -223,12 +223,11 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       exit
     endif
 
-    if (itermax > 4) then
+    if (disk_based_davidson) then
-      itermax = itermax - 1
-    else if (m==1.and.disk_based_davidson) then
       m=0
       disk_based = .True.
-      itermax = 6
+    else if (itermax > 4) then
+      itermax = itermax - 1
     else
       nproc_target = nproc_target - 1
     endif
@@ -267,14 +266,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
 
   if (disk_based) then
     ! Create memory-mapped files for W and S
-    type(c_ptr) :: ptr_w, ptr_s
+    type(mmap_type) :: map_s, map_w
-    integer :: fd_s, fd_w
+
-    call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),&
+    call mmap_create_d('', (/ 1_8*sze, 1_8*N_st_diag*itermax /), .False., .True., map_w)
-        8, fd_w, .False., .True., ptr_w)
+    call mmap_create_s('', (/ 1_8*sze, 1_8*N_st_diag*itermax /), .False., .True., map_s)
-    call mmap(trim(ezfio_work_dir)//'davidson_s', (/int(sze,8),int(N_st_diag*itermax,8)/),&
+    w => map_w%d2
-        4, fd_s, .False., .True., ptr_s)
+    s => map_s%s2
-    call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/))
-    call c_f_pointer(ptr_s, s, (/sze,N_st_diag*itermax/))
   else
     allocate(W(sze,N_st_diag*itermax), S(sze,N_st_diag*itermax))
   endif
@@ -755,13 +752,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
 
   if (disk_based)then
     ! Remove temp files
-    integer, external :: getUnitAndOpen
+    call mmap_destroy(map_w)
-    call munmap( (/int(sze,8),int(N_st_diag*itermax,8)/), 8, fd_w, ptr_w )
+    call mmap_destroy(map_s)
-    fd_w = getUnitAndOpen(trim(ezfio_work_dir)//'davidson_w','r')
-    close(fd_w,status='delete')
-    call munmap( (/int(sze,8),int(N_st_diag*itermax,8)/), 4, fd_s, ptr_s )
-    fd_s = getUnitAndOpen(trim(ezfio_work_dir)//'davidson_s','r')
-    close(fd_s,status='delete')
   else
     deallocate(W,S)
   endif
@@ -774,6 +766,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       lambda                                                         &
       )
   FREE nthreads_davidson
+
 end
 
 

src/davidson/diagonalize_ci.irp.f

@@ -330,6 +330,10 @@ END_PROVIDER
      deallocate(eigenvectors,eigenvalues)
    endif
 
+!    ! Dominant determinants for each states
+!    call print_dominant_det(psi_det,CI_eigenvectors,N_det,N_states,N_int)
+!    call wf_overlap(psi_det,psi_coef,N_states,N_det,psi_det,CI_eigenvectors,N_states,N_det)
+
 END_PROVIDER
 
 subroutine diagonalize_CI

src/davidson/u0_h_u0.irp.f

@@ -179,11 +179,13 @@ subroutine H_u_0_nstates_openmp_work_$N_int(v_t,u_t,N_st,sze,istart,iend,ishift,
 !
 !  compute_singles = (mem+rss > qp_max_mem)
 !
-!  if (.not.compute_singles) then
-!    provide singles_beta_csc
-!  endif
    compute_singles=.True.
 
+   if (.not.compute_singles) then
+     provide singles_alpha_csc singles_beta_csc
+   endif
+
+
 
   maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
   allocate(idx0(maxab))
@@ -287,8 +289,7 @@ compute_singles=.True.
 
       tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
 
-!---
+      if (compute_singles) then
-!      if (compute_singles) then
 
         l_a = psi_bilinear_matrix_columns_loc(lcol)
         ASSERT (l_a <= N_det)
@@ -311,69 +312,67 @@ compute_singles=.True.
             buffer, idx, tmp_det(1,1), j,                              &
             singles_a, n_singles_a )
 
-!-----
+      else
-!      else
+
-!
+ ! Search for singles
-! ! Search for singles
+
-!
+ ! Right boundary
-!call cpu_time(time0)
+        l_a = psi_bilinear_matrix_columns_loc(lcol+1)-1
-! ! Right boundary
+        ASSERT (l_a <= N_det)
-!        l_a = psi_bilinear_matrix_columns_loc(lcol+1)-1
+        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
-!        ASSERT (l_a <= N_det)
+          lrow = psi_bilinear_matrix_rows(l_a)
-!        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
+          ASSERT (lrow <= N_det_alpha_unique)
-!          lrow = psi_bilinear_matrix_rows(l_a)
+
-!          ASSERT (lrow <= N_det_alpha_unique)
+          left = singles_alpha_csc_idx(krow)
-!
+          right_max = -1_8
-!          left = singles_alpha_csc_idx(krow)
+          right = singles_alpha_csc_idx(krow+1)
-!          right_max = -1_8
+          do while (right-left>0_8)
-!          right = singles_alpha_csc_idx(krow+1)
+            k8 = shiftr(right+left,1)
-!          do while (right-left>0_8)
+            if (singles_alpha_csc(k8) > lrow) then
-!            k8 = shiftr(right+left,1)
+              right = k8
-!            if (singles_alpha_csc(k8) > lrow) then
+            else if (singles_alpha_csc(k8) < lrow) then
-!              right = k8
+              left = k8 + 1_8
-!            else if (singles_alpha_csc(k8) < lrow) then
+            else
-!              left = k8 + 1_8
+              right_max = k8+1_8
-!            else
+              exit
-!              right_max = k8+1_8
+            endif
-!              exit
+          enddo
-!            endif
+          if (right_max > 0_8) exit
-!          enddo
+          l_a = l_a-1
-!          if (right_max > 0_8) exit
+        enddo
-!          l_a = l_a-1
+        if (right_max < 0_8) right_max = singles_alpha_csc_idx(krow)
-!        enddo
+
-!        if (right_max < 0_8) right_max = singles_alpha_csc_idx(krow)
+ ! Search
-!
+        n_singles_a = 0
-! ! Search
+        l_a = psi_bilinear_matrix_columns_loc(lcol)
-!        n_singles_a = 0
+        ASSERT (l_a <= N_det)
-!        l_a = psi_bilinear_matrix_columns_loc(lcol)
+
-!        ASSERT (l_a <= N_det)
+        last_found = singles_alpha_csc_idx(krow)
-!
+        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
-!        last_found = singles_alpha_csc_idx(krow)
+          lrow = psi_bilinear_matrix_rows(l_a)
-!        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
+          ASSERT (lrow <= N_det_alpha_unique)
-!          lrow = psi_bilinear_matrix_rows(l_a)
+
-!          ASSERT (lrow <= N_det_alpha_unique)
+          left = last_found
-!
+          right = right_max
-!          left = last_found
+          do while (right-left>0_8)
-!          right = right_max
+            k8 = shiftr(right+left,1)
-!          do while (right-left>0_8)
+            if (singles_alpha_csc(k8) > lrow) then
-!            k8 = shiftr(right+left,1)
+              right = k8
-!            if (singles_alpha_csc(k8) > lrow) then
+            else if (singles_alpha_csc(k8) < lrow) then
-!              right = k8
+              left = k8 + 1_8
-!            else if (singles_alpha_csc(k8) < lrow) then
+            else
-!              left = k8 + 1_8
+              n_singles_a += 1
-!            else
+              singles_a(n_singles_a) = l_a
-!              n_singles_a += 1
+              last_found = k8+1_8
-!              singles_a(n_singles_a) = l_a
+              exit
-!              last_found = k8+1_8
+            endif
-!              exit
+          enddo
-!            endif
+          l_a = l_a+1
-!          enddo
+        enddo
-!          l_a = l_a+1
+        j = j-1
-!        enddo
+
-!        j = j-1
+      endif
-!
+
-!      endif
-!-----
 
       ! Loop over alpha singles
       ! -----------------------

src/davidson/u0_hs2_u0.irp.f

@@ -218,11 +218,14 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
 !
 !  compute_singles = (mem+rss > qp_max_mem)
 !
-!  if (.not.compute_singles) then
-!    provide singles_beta_csc
-!  endif
   compute_singles=.True.
 
+  if (.not.compute_singles) then
+    provide singles_alpha_csc
+    provide singles_beta_csc
+  endif
+
+
 
   maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
   allocate(idx0(maxab))
@@ -314,6 +317,7 @@ compute_singles=.True.
           singles_b(n_singles_b) = singles_beta_csc(k8)
         enddo
       endif
+
     endif
     kcol_prev = kcol
 
@@ -326,8 +330,7 @@ compute_singles=.True.
 
       tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
 
-!---
+      if (compute_singles) then
-!      if (compute_singles) then
 
         l_a = psi_bilinear_matrix_columns_loc(lcol)
         ASSERT (l_a <= N_det)
@@ -352,69 +355,66 @@ compute_singles=.True.
             buffer, idx, tmp_det(1,1), j,                              &
             singles_a, n_singles_a )
 
-!-----
+      else
-!      else
+
-!
+ ! Search for singles
-! ! Search for singles
+
-!
+ ! Right boundary
-!call cpu_time(time0)
+        l_a = psi_bilinear_matrix_columns_loc(lcol+1)-1
-! ! Right boundary
+        ASSERT (l_a <= N_det)
-!        l_a = psi_bilinear_matrix_columns_loc(lcol+1)-1
+        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
-!        ASSERT (l_a <= N_det)
+          lrow = psi_bilinear_matrix_rows(l_a)
-!        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
+          ASSERT (lrow <= N_det_alpha_unique)
-!          lrow = psi_bilinear_matrix_rows(l_a)
+
-!          ASSERT (lrow <= N_det_alpha_unique)
+          left = singles_alpha_csc_idx(krow)
-!
+          right_max = -1_8
-!          left = singles_alpha_csc_idx(krow)
+          right = singles_alpha_csc_idx(krow+1)
-!          right_max = -1_8
+          do while (right-left>0_8)
-!          right = singles_alpha_csc_idx(krow+1)
+            k8 = shiftr(right+left,1)
-!          do while (right-left>0_8)
+            if (singles_alpha_csc(k8) > lrow) then
-!            k8 = shiftr(right+left,1)
+              right = k8
-!            if (singles_alpha_csc(k8) > lrow) then
+            else if (singles_alpha_csc(k8) < lrow) then
-!              right = k8
+              left = k8 + 1_8
-!            else if (singles_alpha_csc(k8) < lrow) then
+            else
-!              left = k8 + 1_8
+              right_max = k8+1_8
-!            else
+              exit
-!              right_max = k8+1_8
+            endif
-!              exit
+          enddo
-!            endif
+          if (right_max > 0_8) exit
-!          enddo
+          l_a = l_a-1
-!          if (right_max > 0_8) exit
+        enddo
-!          l_a = l_a-1
+        if (right_max < 0_8) right_max = singles_alpha_csc_idx(krow)
-!        enddo
+
-!        if (right_max < 0_8) right_max = singles_alpha_csc_idx(krow)
+ ! Search
-!
+        n_singles_a = 0
-! ! Search
+        l_a = psi_bilinear_matrix_columns_loc(lcol)
-!        n_singles_a = 0
+        ASSERT (l_a <= N_det)
-!        l_a = psi_bilinear_matrix_columns_loc(lcol)
+
-!        ASSERT (l_a <= N_det)
+        last_found = singles_alpha_csc_idx(krow)
-!
+        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
-!        last_found = singles_alpha_csc_idx(krow)
+          lrow = psi_bilinear_matrix_rows(l_a)
-!        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
+          ASSERT (lrow <= N_det_alpha_unique)
-!          lrow = psi_bilinear_matrix_rows(l_a)
+
-!          ASSERT (lrow <= N_det_alpha_unique)
+          left = last_found
-!
+          right = right_max
-!          left = last_found
+          do while (right-left>0_8)
-!          right = right_max
+            k8 = shiftr(right+left,1)
-!          do while (right-left>0_8)
+            if (singles_alpha_csc(k8) > lrow) then
-!            k8 = shiftr(right+left,1)
+              right = k8
-!            if (singles_alpha_csc(k8) > lrow) then
+            else if (singles_alpha_csc(k8) < lrow) then
-!              right = k8
+              left = k8 + 1_8
-!            else if (singles_alpha_csc(k8) < lrow) then
+            else
-!              left = k8 + 1_8
+              n_singles_a += 1
-!            else
+              singles_a(n_singles_a) = l_a
-!              n_singles_a += 1
+              last_found = k8+1_8
-!              singles_a(n_singles_a) = l_a
+              exit
-!              last_found = k8+1_8
+            endif
-!              exit
+          enddo
-!            endif
+          l_a = l_a+1
-!          enddo
+        enddo
-!          l_a = l_a+1
+        j = j-1
-!        enddo
+
-!        j = j-1
+      endif
-!
-!      endif
-!-----
 
       ! Loop over alpha singles
       ! -----------------------

src/davidson_keywords/EZFIO.cfg

@@ -48,7 +48,7 @@ default: false
 
 [distributed_davidson]
 type: logical
-doc: If |true|, use the distributed algorithm
+doc: If |true|, use the distributed algorithm. If you plan to run multi-node calculations, set this to true before running.
-default: True
+default: False
 interface: ezfio,provider,ocaml
 

src/determinants/generate_cas_space.irp.f

@@ -33,7 +33,7 @@ subroutine generate_cas_space
   print *, 'CAS(', n_alpha_act+n_beta_act, ', ', n_act_orb, ')'
   print *, ''
 
-  n_det_alpha_unique = binom_int(n_act_orb, n_alpha_act)
+  n_det_alpha_unique = int(binom_int(n_act_orb, n_alpha_act),4)
   TOUCH n_det_alpha_unique
 
   n = n_alpha_act
@@ -56,7 +56,7 @@ subroutine generate_cas_space
     u = ior(t1,t2)
   enddo
 
-  n_det_beta_unique = binom_int(n_act_orb, n_beta_act)
+  n_det_beta_unique = int(binom_int(n_act_orb, n_beta_act),4)
   TOUCH n_det_beta_unique
 
   n = n_beta_act

src/determinants/ref_bitmask.irp.f

@@ -30,31 +30,30 @@
     ref_bitmask_energy += mo_one_e_integrals(occ(i,1),occ(i,1)) + mo_one_e_integrals(occ(i,2),occ(i,2))
     ref_bitmask_kinetic_energy += mo_kinetic_integrals(occ(i,1),occ(i,1)) + mo_kinetic_integrals(occ(i,2),occ(i,2))
     ref_bitmask_n_e_energy += mo_integrals_n_e(occ(i,1),occ(i,1)) + mo_integrals_n_e(occ(i,2),occ(i,2))
+    do j = i+1, elec_alpha_num
+      ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(j,1),occ(i,1))
+      ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(j,1),occ(i,1))
+    enddo
+    do j= 1, elec_alpha_num
+      ref_bitmask_two_e_energy += mo_two_e_integrals_jj(occ(j,1),occ(i,2))
+      ref_bitmask_energy += mo_two_e_integrals_jj(occ(j,1),occ(i,2))
+    enddo
+    do j = i+1, elec_beta_num
+      ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(j,2),occ(i,2))
+      ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(j,2),occ(i,2))
+    enddo
   enddo
 
   do i = elec_beta_num+1,elec_alpha_num
     ref_bitmask_energy += mo_one_e_integrals(occ(i,1),occ(i,1))
     ref_bitmask_kinetic_energy += mo_kinetic_integrals(occ(i,1),occ(i,1))
     ref_bitmask_n_e_energy += mo_integrals_n_e(occ(i,1),occ(i,1))
-  enddo
+    do j = i+1, elec_alpha_num
-
+      ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(j,1),occ(i,1))
-  do j= 1, elec_alpha_num
+      ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(j,1),occ(i,1))
-    do i = j+1, elec_alpha_num
-      ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(i,1),occ(j,1))
-      ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(i,1),occ(j,1))
     enddo
   enddo
 
-  do j= 1, elec_beta_num
-    do i = j+1, elec_beta_num
-      ref_bitmask_two_e_energy += mo_two_e_integrals_jj_anti(occ(i,2),occ(j,2))
-      ref_bitmask_energy += mo_two_e_integrals_jj_anti(occ(i,2),occ(j,2))
-    enddo
-    do i= 1, elec_alpha_num
-      ref_bitmask_two_e_energy += mo_two_e_integrals_jj(occ(i,1),occ(j,2))
-      ref_bitmask_energy += mo_two_e_integrals_jj(occ(i,1),occ(j,2))
-    enddo
-  enddo
   ref_bitmask_one_e_energy = ref_bitmask_kinetic_energy +   ref_bitmask_n_e_energy
 
  ref_bitmask_energy_ab = 0.d0

src/determinants/spindeterminants.irp.f

@@ -910,6 +910,8 @@ subroutine copy_psi_bilinear_to_psi(psi, isize)
 end
 
 
+use mmap_module
+
  BEGIN_PROVIDER [ integer*8, singles_alpha_csc_idx, (N_det_alpha_unique+1) ]
 &BEGIN_PROVIDER [ integer*8, singles_alpha_csc_size ]
   implicit none
@@ -927,10 +929,9 @@ end
 
   !$OMP PARALLEL DEFAULT(NONE)                                 &
       !$OMP   SHARED(N_det_alpha_unique, psi_det_alpha_unique, &
-      !$OMP          idx0, N_int, singles_alpha_csc,                 &
+      !$OMP          idx0, N_int, singles_alpha_csc_idx)       &
-      !$OMP          elec_alpha_num, mo_num, singles_alpha_csc_idx)      &
       !$OMP   PRIVATE(i,s,j)
-  allocate (s(elec_alpha_num * (mo_num-elec_alpha_num) ))
+  allocate (s(N_det_alpha_unique))
   !$OMP DO SCHEDULE(static,64)
   do i=1, N_det_alpha_unique
     call get_all_spin_singles(                                       &
@@ -966,7 +967,7 @@ BEGIN_PROVIDER [ integer, singles_alpha_csc, (singles_alpha_csc_size) ]
   !$OMP PARALLEL DO DEFAULT(NONE)                                    &
       !$OMP   SHARED(N_det_alpha_unique, psi_det_alpha_unique,       &
       !$OMP          idx0, N_int, singles_alpha_csc, singles_alpha_csc_idx)&
-      !$OMP   PRIVATE(i,k) SCHEDULE(static,1)
+      !$OMP   PRIVATE(i,k) SCHEDULE(static)
   do i=1, N_det_alpha_unique
     call get_all_spin_singles(                                       &
         psi_det_alpha_unique, idx0, psi_det_alpha_unique(1,i), N_int,&
@@ -978,7 +979,36 @@ BEGIN_PROVIDER [ integer, singles_alpha_csc, (singles_alpha_csc_size) ]
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ type(mmap_type), singles_alpha_csc_map ]
+  implicit none
+  BEGIN_DOC
+  ! Indices of all single excitations
+  END_DOC
+  integer                        :: i, k
+  integer, allocatable           :: idx0(:)
 
+  call mmap_create_i('', (/ 1_8*singles_alpha_csc_size /), &
+                           .False., .False., singles_alpha_csc_map)
+
+  allocate (idx0(N_det_alpha_unique))
+  do i=1, N_det_alpha_unique
+    idx0(i) = i
+  enddo
+
+  !$OMP PARALLEL DO DEFAULT(NONE)                                    &
+      !$OMP   SHARED(N_det_alpha_unique, psi_det_alpha_unique,       &
+      !$OMP          idx0, N_int, singles_alpha_csc_map, singles_alpha_csc_idx)&
+      !$OMP   PRIVATE(i,k) SCHEDULE(static)
+  do i=1, N_det_alpha_unique
+    call get_all_spin_singles(                                       &
+        psi_det_alpha_unique, idx0, psi_det_alpha_unique(1,i), N_int, N_det_alpha_unique, &
+        singles_alpha_csc_map%i1(singles_alpha_csc_idx(i):singles_alpha_csc_idx(i)+N_det_alpha_unique-1),&
+        k)
+  enddo
+  !$OMP END PARALLEL DO
+  deallocate(idx0)
+
+END_PROVIDER
 
 
  BEGIN_PROVIDER [ integer*8, singles_beta_csc_idx, (N_det_beta_unique+1) ]
@@ -998,11 +1028,10 @@ END_PROVIDER
 
   !$OMP PARALLEL DEFAULT(NONE)                                     &
       !$OMP   SHARED(N_det_beta_unique, psi_det_beta_unique,       &
-      !$OMP          idx0, N_int, singles_beta_csc,                 &
+      !$OMP          idx0, N_int, singles_beta_csc_idx)            &
-      !$OMP          elec_beta_num, mo_num, singles_beta_csc_idx)      &
       !$OMP   PRIVATE(i,s,j)
-  allocate (s(elec_beta_num*(mo_num-elec_beta_num)))
+  allocate (s(N_det_beta_unique))
-  !$OMP DO SCHEDULE(static,1)
+  !$OMP DO SCHEDULE(static)
   do i=1, N_det_beta_unique
     call get_all_spin_singles(                                       &
         psi_det_beta_unique, idx0, psi_det_beta_unique(1,i), N_int,&
@@ -1037,7 +1066,7 @@ BEGIN_PROVIDER [ integer, singles_beta_csc, (singles_beta_csc_size) ]
   !$OMP PARALLEL DO DEFAULT(NONE)                                    &
       !$OMP   SHARED(N_det_beta_unique, psi_det_beta_unique,       &
       !$OMP          idx0, N_int, singles_beta_csc, singles_beta_csc_idx)&
-      !$OMP   PRIVATE(i,k) SCHEDULE(static,64)
+      !$OMP   PRIVATE(i,k) SCHEDULE(static)
   do i=1, N_det_beta_unique
     call get_all_spin_singles(                                       &
         psi_det_beta_unique, idx0, psi_det_beta_unique(1,i), N_int,&
@@ -1049,6 +1078,37 @@ BEGIN_PROVIDER [ integer, singles_beta_csc, (singles_beta_csc_size) ]
 
 END_PROVIDER
 
+BEGIN_PROVIDER [ type(mmap_type), singles_beta_csc_map ]
+  implicit none
+  BEGIN_DOC
+  ! Indices of all single excitations
+  END_DOC
+  integer                        :: i, k
+  integer, allocatable           :: idx0(:)
+
+  call mmap_create_i('', (/ 1_8*singles_beta_csc_size /), &
+                           .False., .False., singles_beta_csc_map)
+
+  allocate (idx0(N_det_beta_unique))
+  do i=1, N_det_beta_unique
+    idx0(i) = i
+  enddo
+
+  !$OMP PARALLEL DO DEFAULT(NONE)                                    &
+      !$OMP   SHARED(N_det_beta_unique, psi_det_beta_unique,       &
+      !$OMP          idx0, N_int, singles_beta_csc_map, singles_beta_csc_idx)&
+      !$OMP   PRIVATE(i,k) SCHEDULE(static)
+  do i=1, N_det_beta_unique
+    call get_all_spin_singles(                                       &
+        psi_det_beta_unique, idx0, psi_det_beta_unique(1,i), N_int, N_det_beta_unique, &
+        singles_beta_csc_map%i1(singles_beta_csc_idx(i):singles_beta_csc_idx(i)+N_det_beta_unique-1),&
+        k)
+  enddo
+  !$OMP END PARALLEL DO
+  deallocate(idx0)
+
+END_PROVIDER
+
 
 
 
@@ -1111,16 +1171,16 @@ subroutine get_all_spin_singles_1(buffer, idx, spindet, size_buffer, singles, n_
   integer                        :: i
   integer(bit_kind)              :: v
   integer                        :: degree
-  integer                        :: add_single(0:64) = (/ 0, 0, 1, 0, 0, (0, i=1,60) /)
   include 'utils/constants.include.F'
 
-  n_singles = 1
+  n_singles = 0
   do i=1,size_buffer
     degree = popcnt(xor( spindet, buffer(i) ))
+    if (degree == 2) then
+      n_singles = n_singles+1
       singles(n_singles) = idx(i)
-    n_singles = n_singles+add_single(degree)
+    endif
   enddo
-  n_singles = n_singles-1
 
 end
 
@@ -1142,15 +1202,15 @@ subroutine get_all_spin_doubles_1(buffer, idx, spindet, size_buffer, doubles, n_
   integer                        :: i
   include 'utils/constants.include.F'
   integer                        :: degree
-  integer                        :: add_double(0:64) = (/ 0, 0, 0, 0, 1, (0, i=1,60) /)
 
-  n_doubles = 1
+  n_doubles = 0
   do i=1,size_buffer
     degree = popcnt(xor( spindet, buffer(i) ))
+    if (degree == 4) then
+      n_doubles = n_doubles+1
       doubles(n_doubles) = idx(i)
-    n_doubles = n_doubles+add_double(degree)
+    endif
   enddo
-  n_doubles = n_doubles-1
 
 end
 
@@ -1181,8 +1241,8 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu
   integer(bit_kind)              :: xorvec($N_int)
   integer                        :: degree
 
-  n_singles = 1
+  n_singles = 0
-  n_doubles = 1
+  n_doubles = 0
   do i=1,size_buffer
 
     do k=1,$N_int
@@ -1196,16 +1256,14 @@ subroutine get_all_spin_singles_and_doubles_$N_int(buffer, idx, spindet, size_bu
     enddo
 
     if ( degree == 4 ) then
-      doubles(n_doubles) = idx(i)
       n_doubles = n_doubles+1
+      doubles(n_doubles) = idx(i)
     else if ( degree == 2 ) then
-      singles(n_singles) = idx(i)
       n_singles = n_singles+1
+      singles(n_singles) = idx(i)
     endif
 
   enddo
-  n_singles = n_singles-1
-  n_doubles = n_doubles-1
 
 end
 
@@ -1230,7 +1288,7 @@ subroutine get_all_spin_singles_$N_int(buffer, idx, spindet, size_buffer, single
   integer(bit_kind)              :: xorvec($N_int)
   integer                        :: degree
 
-  n_singles = 1
+  n_singles = 0
   do i=1,size_buffer
 
     do k=1,$N_int
@@ -1247,11 +1305,10 @@ subroutine get_all_spin_singles_$N_int(buffer, idx, spindet, size_buffer, single
       cycle
     endif
 
-    singles(n_singles) = idx(i)
     n_singles = n_singles+1
+    singles(n_singles) = idx(i)
 
   enddo
-  n_singles = n_singles-1
 
 end
 
@@ -1275,7 +1332,7 @@ subroutine get_all_spin_doubles_$N_int(buffer, idx, spindet, size_buffer, double
   include 'utils/constants.include.F'
   integer(bit_kind)              :: xorvec($N_int)
 
-  n_doubles = 1
+  n_doubles = 0
   do i=1,size_buffer
 
     do k=1,$N_int
@@ -1292,13 +1349,11 @@ subroutine get_all_spin_doubles_$N_int(buffer, idx, spindet, size_buffer, double
       cycle
     endif
 
-    doubles(n_doubles) = idx(i)
     n_doubles = n_doubles+1
+    doubles(n_doubles) = idx(i)
 
   enddo
 
-  n_doubles = n_doubles-1
-
 end
 
 SUBST [ N_int ]

src/ezfio_files/ezfio.irp.f

@@ -60,3 +60,16 @@ BEGIN_PROVIDER [ character*(1024), ezfio_work_dir ]
  ezfio_work_dir = trim(ezfio_filename)//'/work/'
 END_PROVIDER
 
+BEGIN_PROVIDER [ character*(1024), ezfio_work_dir_pid ]
+ use c_functions
+ implicit none
+ BEGIN_DOC
+ ! EZFIO/work/pid_
+ END_DOC
+ character*(32)                 :: pid_str
+ integer :: getpid
+
+ write(pid_str,*) getpid()
+ ezfio_work_dir_pid = trim(ezfio_work_dir)//'/'//trim(pid_str)//'_'
+END_PROVIDER
+

src/gpu/gpu_module.F90

@@ -143,7 +143,7 @@ module gpu
       b, ldb, c, ldc) bind(C, name='gpu_dgeam')
       import
       type(c_ptr), value, intent(in)        :: handle
-      character(c_char), intent(in), value  :: transa, transb
+      character(c_char), intent(in)         :: transa, transb
       integer(c_int64_t), intent(in), value :: m, n, lda, ldb, ldc
       real(c_double), intent(in)            :: alpha, beta
       type(c_ptr), value :: a, b, c
@@ -153,7 +153,7 @@ module gpu
       b, ldb, c, ldc) bind(C, name='gpu_sgeam')
       import
       type(c_ptr), value, intent(in)        :: handle
-      character(c_char), intent(in), value  :: transa, transb
+      character(c_char), intent(in)         :: transa, transb
       integer(c_int64_t), intent(in), value :: m, n, lda, ldb, ldc
       real(c_float), intent(in)             :: alpha, beta
       real(c_float) :: a, b, c
@@ -194,7 +194,7 @@ module gpu
       b, ldb, beta, c, ldc) bind(C, name='gpu_sgemm')
       import
       type(c_ptr), value, intent(in)        :: handle
-      character(c_char), intent(in), value  :: transa, transb
+      character(c_char), intent(in)         :: transa, transb
       integer(c_int64_t), intent(in), value :: m, n, k, lda, ldb, ldc
       real(c_float), intent(in)             :: alpha, beta
       real(c_float) :: a, b, c
@@ -268,8 +268,12 @@ module gpu
       implicit none
       type(gpu_double1), intent(inout) :: ptr
       integer, intent(in) :: s
+      integer*8 :: s_8, n
 
-      call gpu_allocate_c(ptr%c, s*8_8)
+      s_8 = s
+      n = s_8 * 8_8
+
+      call gpu_allocate_c(ptr%c, n)
       call c_f_pointer(ptr%c, ptr%f, (/ s /))
     end subroutine
 
@@ -277,8 +281,13 @@ module gpu
       implicit none
       type(gpu_double2), intent(inout) :: ptr
       integer, intent(in) :: s1, s2
+      integer*8 :: s1_8, s2_8, n
 
-      call gpu_allocate_c(ptr%c, s1*s2*8_8)
+      s1_8 = s1
+      s2_8 = s2
+      n = s1_8 * s2_8 * 8_8
+
+      call gpu_allocate_c(ptr%c, n)
       call c_f_pointer(ptr%c, ptr%f, (/ s1, s2 /))
     end subroutine
 
@@ -286,8 +295,14 @@ module gpu
       implicit none
       type(gpu_double3), intent(inout) :: ptr
       integer, intent(in) :: s1, s2, s3
+      integer*8 :: s1_8, s2_8, s3_8, n
 
-      call gpu_allocate_c(ptr%c, s1*s2*s3*8_8)
+      s1_8 = s1
+      s2_8 = s2
+      s3_8 = s3
+      n = s1_8 * s2_8 * s3_8 * 8_8
+
+      call gpu_allocate_c(ptr%c, n)
       call c_f_pointer(ptr%c, ptr%f, (/ s1, s2, s3 /))
     end subroutine
 
@@ -295,8 +310,15 @@ module gpu
       implicit none
       type(gpu_double4), intent(inout) :: ptr
       integer, intent(in) :: s1, s2, s3, s4
+      integer*8 :: s1_8, s2_8, s3_8, s4_8, n
 
-      call gpu_allocate_c(ptr%c, s1*s2*s3*s4*8_8)
+      s1_8 = s1
+      s2_8 = s2
+      s3_8 = s3
+      s4_8 = s4
+      n = s1_8 * s2_8 * s3_8 * s4_8 * 8_8
+
+      call gpu_allocate_c(ptr%c, n)
       call c_f_pointer(ptr%c, ptr%f, (/ s1, s2, s3, s4 /))
     end subroutine
 
@@ -304,8 +326,16 @@ module gpu
       implicit none
       type(gpu_double5), intent(inout) :: ptr
       integer, intent(in) :: s1, s2, s3, s4, s5
+      integer*8 :: s1_8, s2_8, s3_8, s4_8, s5_8, n
 
-      call gpu_allocate_c(ptr%c, s1*s2*s3*s4*s5*8_8)
+      s1_8 = s1
+      s2_8 = s2
+      s3_8 = s3
+      s4_8 = s4
+      s5_8 = s5
+      n = s1_8 * s2_8 * s3_8 * s4_8 * s5_8 * 8_8
+
+      call gpu_allocate_c(ptr%c, n)
       call c_f_pointer(ptr%c, ptr%f, (/ s1, s2, s3, s4, s5 /))
     end subroutine
 
@@ -313,8 +343,17 @@ module gpu
       implicit none
       type(gpu_double6), intent(inout) :: ptr
       integer, intent(in) :: s1, s2, s3, s4, s5, s6
+      integer*8 :: s1_8, s2_8, s3_8, s4_8, s5_8, s6_8, n
 
-      call gpu_allocate_c(ptr%c, s1*s2*s3*s4*s5*s6*8_8)
+      s1_8 = s1
+      s2_8 = s2
+      s3_8 = s3
+      s4_8 = s4
+      s5_8 = s5
+      s6_8 = s6
+      n = s1_8 * s2_8 * s3_8 * s4_8 * s5_8 * s6_8 * 8_8
+
+      call gpu_allocate_c(ptr%c, n)
       call c_f_pointer(ptr%c, ptr%f, (/ s1, s2, s3, s4, s5, s6 /))
     end subroutine
 

src/hartree_fock/deb_ao_2e_int.irp.f

@@ -0,0 +1,188 @@
+
+program deb_ao_2e_int
+
+  !call check_ao_two_e_integral_cosgtos()
+  call check_crint1()
+  !call check_crint2()
+
+end
+
+! ---
+
+subroutine check_ao_two_e_integral_cosgtos()
+
+  implicit none
+
+  integer                    :: i, j, k, l
+  double precision           :: tmp1, tmp2
+  double precision           :: acc, nrm, dif
+
+  double precision, external :: ao_two_e_integral
+  double precision, external :: ao_two_e_integral_cosgtos
+
+  acc = 0.d0
+  nrm = 0.d0
+
+  i = 1
+  j = 6
+  k = 1
+  l = 16
+!  do i = 1, ao_num
+!    do k = 1, ao_num
+!      do j = 1, ao_num
+!        do l = 1, ao_num
+
+          tmp1 = ao_two_e_integral        (i, j, k, l)
+          tmp2 = ao_two_e_integral_cosgtos(i, j, k, l)
+
+          dif  = dabs(tmp1 - tmp2)
+          if(dif .gt. 1d-12) then
+            print*, ' error on:', i, j, k, l
+            print*, tmp1, tmp2, dif
+            stop
+          endif
+
+!          acc += dif
+!          nrm += dabs(tmp1)
+!        enddo
+!      enddo
+!    enddo
+!  enddo
+
+  print *, ' acc (%) = ', dif * 100.d0 / nrm
+
+end
+
+! ---
+
+subroutine check_crint1()
+
+  implicit none
+  integer                    :: i, n, i_rho
+  double precision           :: dif_thr
+  double precision           :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im
+  complex*16                 :: rho_test(1:10) = (/ (1d-12, 0.d0),  &
+                                                    (+1d-9, +1d-6), &
+                                                    (-1d-6, -1d-5), &
+                                                    (+1d-3, -1d-2), &
+                                                    (-1d-1, +1d-1), &
+                                                    (+1d-0, +1d-1), &
+                                                    (-1d+1, +1d+1), &
+                                                    (+1d+2, +1d+1), &
+                                                    (-1d+3, +1d+2), &
+                                                    (+1d+4, +1d+4) /)
+  complex*16                 :: rho
+  complex*16                 :: int_an, int_nm
+  double precision, external :: rint
+  complex*16, external       :: crint_1, crint_2, crint_quad
+
+  n = 10
+  dif_thr = 1d-7
+
+  do i_rho = 8, 10
+  !do i_rho = 7, 7
+  
+    !rho = (-10.d0, 0.1d0)
+    !rho = (+10.d0, 0.1d0)
+    rho = rho_test(i_rho)
+    print*, "rho = ", real(rho), aimag(rho)
+
+    acc_re = 0.d0
+    nrm_re = 0.d0
+    acc_im = 0.d0
+    nrm_im = 0.d0
+    do i = 0, n
+      !int_an = crint_1   (i, rho)
+      int_an = crint_2   (i, rho)
+      int_nm = crint_quad(i, rho)
+      
+      dif_re = dabs(real(int_an) - real(int_nm))
+      dif_im = dabs(aimag(int_an) - aimag(int_nm))
+
+      if((dif_re .gt. dif_thr) .or. (dif_im .gt. dif_thr)) then
+        print*, ' error on i =', i
+        print*, real(int_an), real(int_nm), dif_re
+        print*, aimag(int_an), aimag(int_nm), dif_im
+        !print*, rint(i, real(rho))
+        print*, crint_1(i, rho)
+        !print*, crint_2(i, rho)
+        stop
+      endif
+      acc_re += dif_re
+      nrm_re += dabs(real(int_nm))
+      acc_im += dif_im
+      nrm_im += dabs(aimag(int_nm))
+    enddo
+  
+    print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re+1d-15)
+    print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im+1d-15)
+  enddo
+
+end
+
+! ---
+
+subroutine check_crint2()
+
+  implicit none
+
+  integer              :: i, n, i_rho
+  double precision     :: dif_thr
+  double precision     :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im
+  complex*16           :: rho_test(1:10) = (/ (1d-12, 0.d0),  &
+                                              (+1d-9, +1d-6), &
+                                              (-1d-6, -1d-5), &
+                                              (+1d-3, -1d-2), &
+                                              (-1d-1, +1d-1), &
+                                              (+1d-0, +1d-1), &
+                                              (-1d+1, +1d+1), &
+                                              (+1d+2, +1d+1), &
+                                              (-1d+3, +1d+2), &
+                                              (+1d+4, +1d+4) /)
+  complex*16           :: rho
+  complex*16           :: int_an, int_nm
+  complex*16, external :: crint_1, crint_2
+
+  n = 30
+  dif_thr = 1d-12
+
+  do i_rho = 1, 10
+    rho = rho_test(i_rho)
+    print*, "rho = ", real(rho), aimag(rho)
+
+    acc_re = 0.d0
+    nrm_re = 0.d0
+    acc_im = 0.d0
+    nrm_im = 0.d0
+    do i = 0, n
+      int_an = crint_1(i, rho)
+      int_nm = crint_2(i, rho)
+      
+      dif_re = dabs(real(int_an) - real(int_nm))
+      !if(dif_re .gt. dif_thr) then
+      !  print*, ' error in real part:', i
+      !  print*, real(int_an), real(int_nm), dif_re
+      !  stop
+      !endif
+      acc_re += dif_re
+      nrm_re += dabs(real(int_nm))
+  
+      dif_im = dabs(aimag(int_an) - aimag(int_nm))
+      !if(dif_im .gt. dif_thr) then
+      !  print*, ' error in imag part:', i
+      !  print*, aimag(int_an), aimag(int_nm), dif_im
+      !  stop
+      !endif
+      acc_im += dif_im
+      nrm_im += dabs(aimag(int_nm))
+    enddo
+  
+    print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re+1d-15)
+    print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im+1d-15)
+  enddo
+
+end
+
+! ---
+
+

src/hartree_fock/hf_energy.irp.f

@@ -19,16 +19,41 @@ END_PROVIDER
  ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components.
  END_DOC
  integer :: i,j
- HF_energy = nuclear_repulsion
+ double precision :: tmp1, tmp2
+ HF_energy = 0.d0
  HF_two_electron_energy = 0.d0
  HF_one_electron_energy = 0.d0
  do j=1,ao_num
    do i=1,ao_num
-    HF_two_electron_energy += 0.5d0 * ( ao_two_e_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) &
+    tmp1 = 0.5d0 * ( ao_two_e_integral_alpha(i,j) * SCF_density_matrix_ao_alpha(i,j) &
                     +ao_two_e_integral_beta (i,j) * SCF_density_matrix_ao_beta (i,j) )
-    HF_one_electron_energy += ao_one_e_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) )
+    tmp2 = ao_one_e_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) )
+    HF_two_electron_energy += tmp1
+    HF_one_electron_energy += tmp2
+    HF_energy += tmp1 + tmp2
+   enddo
+ enddo
+ HF_energy += nuclear_repulsion
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, HF_kinetic_energy]
+&BEGIN_PROVIDER [ double precision, HF_n_e_energy]
+ implicit none
+ BEGIN_DOC
+ ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components.
+ END_DOC
+ integer :: i,j
+ double precision :: tmp1, tmp2
+ HF_n_e_energy = 0.d0
+ HF_kinetic_energy = 0.d0
+ do j=1,ao_num
+   do i=1,ao_num
+    tmp1 = ao_integrals_n_e(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) )
+    tmp2 = ao_kinetic_integrals(i,j) * (SCF_density_matrix_ao_alpha(i,j) + SCF_density_matrix_ao_beta (i,j) )
+    HF_n_e_energy += tmp1
+    HF_kinetic_energy += tmp2
    enddo
  enddo
- HF_energy += HF_two_electron_energy + HF_one_electron_energy
 END_PROVIDER
 

src/hartree_fock/print_scf_int.irp.f

@@ -1,114 +0,0 @@
-
-program print_scf_int
-
-  call main()
-
-end
-
-subroutine main()
-
-  implicit none
-  integer ::  i, j, k, l
-
-  print *, " Hcore:"
-  do j = 1, ao_num
-    do i = 1, ao_num
-      print *, i, j, ao_one_e_integrals(i,j)
-    enddo
-  enddo
-
-  print *, " P:"
-  do j = 1, ao_num
-    do i = 1, ao_num
-      print *, i, j, SCF_density_matrix_ao_alpha(i,j)
-    enddo
-  enddo
-
-
-  double precision :: integ, density_a, density_b, density
-  double precision :: J_scf(ao_num, ao_num)
-  double precision :: K_scf(ao_num, ao_num)
-
-
-  double precision, external :: get_ao_two_e_integral
-  PROVIDE ao_integrals_map
-
-  print *, " J:"
-  !do j = 1, ao_num
-  !  do l = 1, ao_num
-  !    do i = 1, ao_num
-  !      do k = 1, ao_num
-  !        !  < 1:k, 2:l | 1:i, 2:j > 
-  !        print *, '< k l | i j >', k, l, i, j
-  !        print *, get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
-  !      enddo
-  !    enddo
-  !  enddo
-  !enddo
-
-  !do k = 1, ao_num
-  !  do i = 1, ao_num
-  !    do j = 1, ao_num
-  !      do l = 1, ao_num
-  !        !  ( 1:k, 1:i | 2:l, 2:j ) 
-  !        print *, '(k i | l j)', k, i, l, j
-  !        print *, get_ao_two_e_integral(l, j, k, i, ao_integrals_map)
-  !      enddo
-  !    enddo
-  !    print *, ''
-  !  enddo
-  !enddo
-
-  J_scf = 0.d0
-  K_scf = 0.d0
-  do i = 1, ao_num
-    do k = 1, ao_num
-      do j = 1, ao_num
-        do l = 1, ao_num
-
-          density_a = SCF_density_matrix_ao_alpha(l,j)
-          density_b = SCF_density_matrix_ao_beta (l,j)
-          density   = density_a + density_b
-
-          integ = get_ao_two_e_integral(l, j, k, i, ao_integrals_map)
-          J_scf(k,i) += density   * integ
-          integ = get_ao_two_e_integral(l, i, k, j, ao_integrals_map)
-          K_scf(k,i) -= density_a * integ
-        enddo
-      enddo
-    enddo
-  enddo
-
-  print *, 'J x P'
-  do i = 1, ao_num
-    do k = 1, ao_num
-      print *, k, i, J_scf(k,i)
-    enddo
-  enddo
-
-  print *, ''
-  print *, 'K x P'
-  do i = 1, ao_num
-    do k = 1, ao_num
-      print *, k, i, K_scf(k,i)
-    enddo
-  enddo
-
-  print *, ''
-  print *, 'F in AO'
-  do i = 1, ao_num
-    do k = 1, ao_num
-      print *, k, i, Fock_matrix_ao(k,i)
-    enddo
-  enddo
-
-  print *, ''
-  print *, 'F in MO'
-  do i = 1, ao_num
-    do k = 1, ao_num
-      print *, k, i, 2.d0 * Fock_matrix_mo_alpha(k,i)
-    enddo
-  enddo
-
-end
-

src/mo_basis/mos.irp.f

@@ -277,7 +277,7 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
       T, ao_num,                                                     &
       0.d0, A_mo, size(A_mo,1))
 
-  call restore_symmetry(mo_num,mo_num,A_mo,size(A_mo,1),1.d-12)
+  call restore_symmetry(mo_num,mo_num,A_mo,size(A_mo,1),1.d-15)
   deallocate(T)
 end
 

src/mo_one_e_ints/mo_one_e_ints.irp.f

@@ -18,6 +18,6 @@ BEGIN_PROVIDER [ double precision, mo_one_e_integrals,(mo_num,mo_num)]
         call ezfio_set_mo_one_e_ints_mo_one_e_integrals(mo_one_e_integrals)
        print *,  'MO one-e integrals written to disk'
   ENDIF
-  call nullify_small_elements(mo_num,mo_num,mo_one_e_integrals,size(mo_one_e_integrals,1),1.d-10)
+  call nullify_small_elements(mo_num,mo_num,mo_one_e_integrals,size(mo_one_e_integrals,1),1.d-15)
 
 END_PROVIDER

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -70,6 +70,10 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
     else
       call add_integrals_to_map(full_ijkl_bitmask_4)
     endif
+    double precision, external     :: map_mb
+    print*,'Molecular integrals provided:'
+    print*,' Size of MO map           ', map_mb(mo_integrals_map) ,'MB'
+    print*,' Number of MO integrals: ',  mo_map_size
   endif
 
   call wall_time(wall_2)
@@ -78,10 +82,6 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
   integer*8                      :: get_mo_map_size, mo_map_size
   mo_map_size = get_mo_map_size()
 
-  double precision, external     :: map_mb
-  print*,'Molecular integrals provided:'
-  print*,' Size of MO map           ', map_mb(mo_integrals_map) ,'MB'
-  print*,' Number of MO integrals: ',  mo_map_size
   print*,' cpu  time :',cpu_2 - cpu_1, 's'
   print*,' wall time :',wall_2 - wall_1, 's  ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
 

src/mol_properties/print_properties.irp.f

@@ -100,7 +100,7 @@ subroutine print_transition_dipole_moment
       dip_str = d_x**2 + d_y**2 + d_z**2
       d = multi_s_dipole_moment(istate,jstate)
       f = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
-      write(*,'(I4,I4,A4,I3,6(F12.6))') (istate-1), (jstate-1), '  ->', (istate-1), d_x, d_y, d_z, d, dip_str, f
+      write(*,'(I4,I4,A4,I3,6(F12.6))') (jstate -1) * (2*N_states-jstate)/2 + istate - jstate, (jstate-1), '  ->', (istate-1), d_x, d_y, d_z, d, dip_str, f
     enddo
   enddo
 
@@ -117,7 +117,7 @@ subroutine print_transition_dipole_moment
       dip_str = d_x**2 + d_y**2 + d_z**2
       f = 2d0/3d0 * d * d * dabs(ci_energy_no_diag(istate) - ci_energy_no_diag(jstate))
       d = multi_s_dipole_moment(istate,jstate) * au_to_D
-      write(*,'(I4,I4,A4,I3,6(F12.6))') (istate-1), (jstate-1), '  ->', (istate-1), d_x, d_y, d_z, d, dip_str, f
+      write(*,'(I4,I4,A4,I3,6(F12.6))') (jstate -1) * (2*N_states-jstate)/2 + istate - jstate, (jstate-1), '  ->', (istate-1), d_x, d_y, d_z, d, dip_str, f
     enddo
   enddo
   print*,'=============================================='
@@ -181,10 +181,9 @@ subroutine print_oscillator_strength
       ! Mixed gauge
       f_m = 2d0/3d0 * d * v
 
-      write(*,'(A19,I3,A9,F10.6,A5,F7.1,A10,F9.6,A6,F9.6,A6,F9.6,A8,F7.3)') '   #  Transition n.', (istate-1), ': Excit.=', dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*ha_to_ev, &
+      write(*,'(A19,I3,A9,F10.6,A5,F7.1,A10,F9.6,A6,F9.6,A6,F9.6,A8,F7.3)') '   #  Transition n.', (jstate -1) * (2*N_states-jstate)/2 + istate - jstate, ': Excit.=', dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*ha_to_ev, &
       ' eV ( ',dabs((ci_energy_no_diag(istate) - ci_energy_no_diag(jstate)))*Ha_to_nm,' nm), f_l=',f_l, ', f_v=', f_v, ', f_m=', f_m, ', <S^2>=', s2_values(istate)
       !write(*,'(I4,I4,A4,I3,A6,F6.1,A6,F6.1)') (istate-1), (jstate-1), '  ->', (istate-1), ', %T1=', percent_exc(2,istate), ', %T2=',percent_exc(3,istate)
-  
     enddo
   enddo
 

src/mu_of_r/test_proj_op.irp.f

@@ -12,15 +12,12 @@ program projected_operators
   mu_of_r_potential = "cas_full"
   touch mu_of_r_potential 
   print*,'Using Valence Only functions'
-!  call test_f_HF_valence_ab
+  call test_f_HF_valence_ab
-!  call routine_full_mos
+  call routine_full_mos
-!   call test_f_ii_valence_ab
+  call test_f_ii_valence_ab
-!   call test_f_ia_valence_ab
+  call test_f_ia_valence_ab
-!  call test_f_ii_ia_aa_valence_ab
+  call test_f_ii_ia_aa_valence_ab
-! call test
+  call test
-!  call test_f_mean_field
-! call test_grad_f_mean_field
- call test_grad_mu_mf
 end
 
 
@@ -39,138 +36,3 @@ subroutine test
 
 end
 
-subroutine test_f_mean_field
- implicit none
- integer :: i_point
- double precision :: weight,r(3)
- double precision :: ref_f, new_f, accu_f
- double precision :: ref_two_dens, new_two_dens, accu_two_dens, dm_a, dm_b
- accu_f = 0.d0
- accu_two_dens = 0.d0
- do i_point = 1, n_points_final_grid
-  r(1:3)   = final_grid_points(1:3,i_point)
-  weight = final_weight_at_r_vector(i_point)
-  call get_f_mf_ab(r,new_f,new_two_dens, dm_a, dm_b)
-  call f_HF_valence_ab(r,r,ref_f,ref_two_dens)
-  accu_f += weight * dabs(new_f- ref_f)
-  accu_two_dens += weight * dabs(new_two_dens - ref_two_dens)
- enddo
- print*,'accu_f        = ',accu_f
- print*,'accu_two_dens = ',accu_two_dens
-
-end
-
-subroutine test_grad_f_mean_field
- implicit none
- integer :: i_point,k
- double precision :: weight,r(3)
- double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3)
- double precision :: grad_dm_a(3), grad_dm_b(3)
- double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b
-
- double precision :: num_grad_f_mf_ab(3), num_grad_two_bod_dens(3)
- double precision :: num_grad_dm_a(3), num_grad_dm_b(3)
- double precision :: f_mf_ab_p,f_mf_ab_m
- double precision :: two_bod_dens_p, two_bod_dens_m
- double precision :: dm_a_p, dm_a_m
- double precision :: dm_b_p, dm_b_m
- double precision :: rbis(3), dr
- double precision :: accu_grad_f_mf_ab(3),accu_grad_two_bod_dens(3)
- double precision :: accu_grad_dm_a(3),accu_grad_dm_b(3)
- double precision :: accu_f_mf_ab, accu_two_bod_dens, accu_dm_a, accu_dm_b
- dr = 0.00001d0
- accu_f_mf_ab = 0.d0 
- accu_two_bod_dens = 0.d0 
- accu_dm_a = 0.d0 
- accu_dm_b = 0.d0
-
- accu_grad_f_mf_ab = 0.d0
- accu_grad_two_bod_dens = 0.d0
- accu_grad_dm_a = 0.d0
- accu_grad_dm_b = 0.d0
- do i_point = 1, n_points_final_grid
-  r(1:3)   = final_grid_points(1:3,i_point)
-  weight = final_weight_at_r_vector(i_point)
-  call get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b)
-  call get_f_mf_ab(r,f_mf_ab_p,two_bod_dens_p, dm_a_p, dm_b_p)
-  accu_f_mf_ab += weight * dabs(f_mf_ab - f_mf_ab_p)
-  accu_two_bod_dens += weight * dabs(two_bod_dens - two_bod_dens_p)
-  accu_dm_a += weight*dabs(dm_a - dm_a_p)
-  accu_dm_b += weight*dabs(dm_b - dm_b_p)
-  do k = 1, 3
-   rbis = r
-   rbis(k) += dr
-   call get_f_mf_ab(rbis,f_mf_ab_p,two_bod_dens_p, dm_a_p, dm_b_p)
-   rbis = r
-   rbis(k) -= dr
-   call get_f_mf_ab(rbis,f_mf_ab_m,two_bod_dens_m, dm_a_m, dm_b_m)
-   num_grad_f_mf_ab(k) = (f_mf_ab_p - f_mf_ab_m)/(2.d0*dr)
-   num_grad_two_bod_dens(k) = (two_bod_dens_p - two_bod_dens_m)/(2.d0*dr)
-   num_grad_dm_a(k) = (dm_a_p - dm_a_m)/(2.d0*dr)
-   num_grad_dm_b(k) = (dm_b_p - dm_b_m)/(2.d0*dr)
-  enddo
-  do k = 1, 3
-   accu_grad_f_mf_ab(k) += weight * dabs(grad_f_mf_ab(k) - num_grad_f_mf_ab(k))
-   accu_grad_two_bod_dens(k) += weight * dabs(grad_two_bod_dens(k) - num_grad_two_bod_dens(k))
-   accu_grad_dm_a(k) += weight * dabs(grad_dm_a(k) - num_grad_dm_a(k))
-   accu_grad_dm_b(k) += weight * dabs(grad_dm_b(k) - num_grad_dm_b(k))
-  enddo
- enddo
- print*,'accu_f_mf_ab = ',accu_f_mf_ab
- print*,'accu_two_bod_dens = ',accu_two_bod_dens
- print*,'accu_dm_a = ',accu_dm_a
- print*,'accu_dm_b = ',accu_dm_b
- print*,'accu_grad_f_mf_ab = '
- print*,accu_grad_f_mf_ab
- print*,'accu_grad_two_bod_dens = '
- print*,accu_grad_two_bod_dens
- print*,'accu_dm_a = '
- print*,accu_grad_dm_a
- print*,'accu_dm_b = '
- print*,accu_grad_dm_b
-
-end
-
-subroutine test_grad_mu_mf
- implicit none
- integer :: i_point,k
- double precision :: weight,r(3),rbis(3)
- double precision :: mu_mf, dm,grad_mu_mf(3), grad_dm(3)
- double precision :: mu_mf_p, mu_mf_m, dm_m, dm_p, num_grad_mu_mf(3),dr, num_grad_dm(3)
- double precision :: accu_mu, accu_dm, accu_grad_dm(3), accu_grad_mu_mf(3)
- dr = 0.00001d0
- accu_grad_mu_mf = 0.d0
- accu_mu = 0.d0
- accu_grad_dm = 0.d0
- accu_dm = 0.d0
- do i_point = 1, n_points_final_grid
-  r(1:3)   = final_grid_points(1:3,i_point)
-  weight = final_weight_at_r_vector(i_point)
-  call grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm)
-  call mu_of_r_mean_field(r,mu_mf_p, dm_p)
-  accu_mu += weight*dabs(mu_mf_p - mu_mf)
-  accu_dm += weight*dabs(dm_p - dm)
-  do k = 1, 3
-   rbis = r
-   rbis(k) += dr
-   call mu_of_r_mean_field(rbis,mu_mf_p, dm_p)
-   rbis = r
-   rbis(k) -= dr
-   call mu_of_r_mean_field(rbis,mu_mf_m, dm_m)
-
-   num_grad_mu_mf(k) = (mu_mf_p - mu_mf_m)/(2.d0*dr)
-   num_grad_dm(k) = (dm_p - dm_m)/(2.d0*dr)
-  enddo
-  do k = 1, 3
-   accu_grad_dm(k)+= weight *dabs(num_grad_dm(k) - grad_dm(k))
-   accu_grad_mu_mf(k)+= weight *dabs(num_grad_mu_mf(k) - grad_mu_mf(k))
-  enddo
- enddo
- print*,'accu_mu = ',accu_mu
- print*,'accu_dm = ',accu_dm
- print*,'accu_grad_dm = '
- print*, accu_grad_dm
- print*,'accu_grad_mu_mf = '
- print*, accu_grad_mu_mf
-
-end

src/scf_utils/EZFIO.cfg

@@ -45,6 +45,12 @@ type: double precision
 doc: Calculated HF energy
 interface: ezfio
 
+[do_mom]
+type: logical
+doc: If true, this will run a MOM calculation. The overlap will be computed at each step with respect to the initial MOs. After an initial Hartree-Fock calculation, the guess can be created by swapping molecular orbitals through the qp run swap_mos command. 
+interface: ezfio,provider,ocaml
+default: False
+
 [frozen_orb_scf]
 type: logical
 doc: If true, leave untouched all the orbitals defined as core and optimize all the orbitals defined as active with qp_set_mo_class

src/scf_utils/fock_matrix.irp.f

@@ -253,17 +253,18 @@ BEGIN_PROVIDER [ double precision, SCF_energy ]
  BEGIN_DOC
  ! Hartree-Fock energy
  END_DOC
- SCF_energy = nuclear_repulsion
-
  integer                        :: i,j
+
+ SCF_energy = 0.d0
+
  do j=1,ao_num
    do i=1,ao_num
-     SCF_energy += 0.5d0 * (                                          &
+     SCF_energy +=  &
          (ao_one_e_integrals(i,j) + Fock_matrix_ao_alpha(i,j) ) *  SCF_density_matrix_ao_alpha(i,j) +&
-         (ao_one_e_integrals(i,j) + Fock_matrix_ao_beta (i,j) ) *  SCF_density_matrix_ao_beta (i,j) )
+         (ao_one_e_integrals(i,j) + Fock_matrix_ao_beta (i,j) ) *  SCF_density_matrix_ao_beta (i,j) 
    enddo
  enddo
- SCF_energy += extra_e_contrib_density
+ SCF_energy = 0.5d0 * SCF_energy + extra_e_contrib_density + nuclear_repulsion 
 
 END_PROVIDER
 

src/scf_utils/reorder_mo_max_overlap.irp.f

@@ -0,0 +1,96 @@
+ subroutine reorder_mo_max_overlap
+  implicit none
+   BEGIN_DOC
+   ! routines that compute the projection of each MO of the current `mo_coef` on the space spanned by the occupied orbitals of `mo_coef_begin_iteration`
+   END_DOC
+   integer :: i,j,k,l
+   double precision, allocatable :: overlap(:,:)
+   double precision, allocatable :: proj(:)
+   integer, allocatable :: iorder(:)
+   double precision, allocatable :: mo_coef_tmp(:,:)
+   double precision, allocatable :: tmp(:,:)
+   allocate(overlap(mo_num,mo_num),proj(mo_num),iorder(mo_num),mo_coef_tmp(ao_num,mo_num),tmp(mo_num,ao_num))
+   
+   overlap(:,:) = 0d0
+   mo_coef_tmp(:,:) = 0d0
+   proj(:) = 0d0
+   iorder(:) = 0d0
+   tmp(:,:) = 0d0
+   
+   ! These matrix products compute the overlap bewteen the initial and the current MOs
+   call dgemm('T','N', mo_num, ao_num, ao_num, 1.d0, &
+        mo_coef_begin_iteration, size(mo_coef_begin_iteration,1), &
+        ao_overlap, size(ao_overlap,1), 0.d0, &
+        tmp, size(tmp,1))
+
+   call dgemm('N','N', mo_num, mo_num, ao_num, 1.d0, &
+        tmp, size(tmp,1), &
+        mo_coef, size(mo_coef, 1), 0.d0, &
+        overlap, size(overlap,1) )
+   
+   
+   ! for each orbital compute the best overlap
+   do i = 1, mo_num
+      iorder(i) = i ! initialize the iorder list as we need it to sort later
+      do j = 1, elec_alpha_num
+         proj(i) += overlap(j,i)*overlap(j,i) ! compute the projection of current orbital i on the occupied space of the initial orbitals
+      enddo
+      proj(i) = dsqrt(proj(i))
+   enddo
+   ! sort the list of projection to find the mos with the largest overlap
+   call dsort(proj(:),iorder(:),mo_num)
+   ! reorder orbitals according to projection
+   do i=1,mo_num
+      mo_coef_tmp(:,i) = mo_coef(:,iorder(mo_num+1-i)) 
+   enddo
+
+   ! update the orbitals
+   mo_coef(:,:) =  mo_coef_tmp(:,:)
+   
+   ! if the determinant is open-shell we need to make sure that the singly occupied orbital correspond to the initial ones
+   if (elec_alpha_num > elec_beta_num) then
+      double precision, allocatable :: overlap_alpha(:,:)
+      double precision, allocatable :: proj_alpha(:)
+      integer, allocatable :: iorder_alpha(:)
+      allocate(overlap_alpha(mo_num,elec_alpha_num),proj_alpha(elec_alpha_num),iorder_alpha(elec_alpha_num))
+      overlap_alpha(:,:) = 0d0
+      mo_coef_tmp(:,:) = 0d0
+      proj_alpha(:) = 0d0
+      iorder_alpha(:) = 0d0
+      tmp(:,:) = 0d0 
+      ! These matrix products compute the overlap bewteen the initial and the current MOs
+      call dgemm('T','N', mo_num, ao_num, ao_num, 1.d0, &
+           mo_coef_begin_iteration, size(mo_coef_begin_iteration,1), &
+           ao_overlap, size(ao_overlap,1), 0.d0, &
+           tmp, size(tmp,1))
+
+      call dgemm('N','N', mo_num, elec_alpha_num, ao_num, 1.d0, &
+           tmp, size(tmp,1), &
+           mo_coef, size(mo_coef, 1), 0.d0, &
+           overlap_alpha, size(overlap_alpha,1) )
+
+      do i = 1, elec_alpha_num
+         iorder_alpha(i) = i ! initialize the iorder list as we need it to sort later
+         do j = 1, elec_beta_num
+            proj_alpha(i) += overlap_alpha(j,i)*overlap_alpha(j,i) ! compute the projection of current orbital i on the beta occupied space of the initial orbitals
+         enddo
+         proj_alpha(i) = dsqrt(proj_alpha(i))
+      enddo
+      
+      ! sort the list of projection to find the mos with the largest overlap
+      call dsort(proj_alpha(:),iorder_alpha(:),elec_alpha_num)
+      ! reorder orbitals according to projection
+      do i=1,elec_alpha_num
+         mo_coef_tmp(:,i) = mo_coef(:,iorder_alpha(elec_alpha_num+1-i)) 
+      enddo
+      do i=1,elec_alpha_num
+         mo_coef(:,i) = mo_coef_tmp(:,i) 
+      enddo
+
+      deallocate(overlap_alpha, proj_alpha, iorder_alpha)
+   endif
+   
+   deallocate(overlap, proj, iorder, mo_coef_tmp, tmp)
+
+ end
+

src/scf_utils/roothaan_hall_scf.irp.f

@@ -51,6 +51,11 @@ END_DOC
 !
   PROVIDE FPS_SPF_matrix_AO Fock_matrix_AO 
 
+  ! Initialize MO to run IMOM
+  if(do_mom)then
+     call initialize_mo_coef_begin_iteration
+  endif
+  
   converged = .False.
   do while ( .not.converged .and. (iteration_SCF < n_it_SCF_max) )
 
@@ -89,15 +94,16 @@ END_DOC
       TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta
 
    endif
-
     MO_coef = eigenvectors_Fock_matrix_MO
+    if(do_mom)then
+       call reorder_mo_max_overlap
+    endif
     if(frozen_orb_scf)then
       call reorder_core_orb
       call initialize_mo_coef_begin_iteration
     endif
 
     TOUCH MO_coef
-
 !   Calculate error vectors
 
     max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_MO))
@@ -106,13 +112,14 @@ END_DOC
 
     energy_SCF = SCF_energy
     Delta_Energy_SCF = energy_SCF - energy_SCF_previous
-    if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) then
+    if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0).and.(.not.do_mom) ) then
       Fock_matrix_AO(1:ao_num,1:ao_num) = Fock_matrix_DIIS (1:ao_num,1:ao_num,index_dim_DIIS)
       Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0
       Fock_matrix_AO_beta  = Fock_matrix_AO*0.5d0
       TOUCH Fock_matrix_AO_alpha Fock_matrix_AO_beta
     endif
 
+    if (.not.do_mom) then
        double precision :: level_shift_save
        level_shift_save = level_shift
        mo_coef_save(1:ao_num,1:mo_num) = mo_coef(1:ao_num,1:mo_num)
@@ -125,6 +132,9 @@ END_DOC
           endif
           TOUCH mo_coef level_shift
           mo_coef(1:ao_num,1:mo_num) = eigenvectors_Fock_matrix_MO(1:ao_num,1:mo_num)
+          if(do_mom)then
+             call reorder_mo_max_overlap
+          endif
           if(frozen_orb_scf)then
              call reorder_core_orb
              call initialize_mo_coef_begin_iteration
@@ -141,6 +151,7 @@ END_DOC
        enddo
        level_shift = level_shift * 0.5d0
        SOFT_TOUCH level_shift
+    endif
     energy_SCF_previous = energy_SCF
 
     converged = ( (max_error_DIIS <= threshold_DIIS_nonzero) .and. &
@@ -228,6 +239,9 @@ END_DOC
     i = j+1
   enddo
 
+  if(do_mom)then
+     call reorder_mo_max_overlap
+  endif
 
   call save_mos
 

src/tools/print_energy.irp.f

@@ -16,4 +16,20 @@ subroutine run
  implicit none
  call print_mol_properties
  print *, psi_energy + nuclear_repulsion
+ call print_energy_components
+ print *, 'E(HF) = ', HF_energy
+ print *, 'E(CI) = ', psi_energy + nuclear_repulsion
+ print *, ''
+ print *, 'E_kin(CI) = ', ref_bitmask_kinetic_energy
+ print *, 'E_kin(HF) = ', HF_kinetic_energy
+ print *, ''
+ print *, 'E_ne (CI) = ', ref_bitmask_n_e_energy
+ print *, 'E_ne (HF) = ', HF_n_e_energy
+ print *, ''
+ print *, 'E_1e (CI) = ', ref_bitmask_one_e_energy
+ print *, 'E_1e (HF) = ', HF_one_electron_energy
+ print *, ''
+ print *, 'E_2e (CI) = ', ref_bitmask_two_e_energy
+ print *, 'E_2e (HF) = ', HF_two_electron_energy
+
 end

src/tools/truncate_wf.irp.f

@@ -56,9 +56,14 @@ subroutine routine_s2
  double precision :: accu(N_states)
 
  print *,  'Weights of the CFG'
- do i=1,N_det
+ integer :: step
+
+ step = max(1,N_det/100)
+ do i=1,N_det-1,step
    print *,  i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
  enddo
+ i=N_det
+ print *,  i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
  print*, 'Min weight of the configuration?'
  read(5,*) wmin
 

src/trexio/import_trexio_integrals.irp.f

@@ -46,6 +46,8 @@ subroutine run(f)
   double precision, allocatable :: tmp(:,:,:)
   integer*8 :: offset, icount
 
+  integer :: k_num
+
   integer, external :: getUnitAndOpen
 
   if (trexio_has_nucleus_repulsion(f) == TREXIO_SUCCESS) then
@@ -163,7 +165,8 @@ subroutine run(f)
 
           deallocate(Vi, V, tmp)
           print *, 'Cholesky AO integrals read from TREXIO file'
-      endif
+
+      else
 
           rc = trexio_has_ao_2e_int_eri(f)
           if (rc /= TREXIO_HAS_NOT) then
@@ -204,6 +207,7 @@ subroutine run(f)
               deallocate(buffer_i, buffer_values, Vi, V)
               print *, 'AO integrals read from TREXIO file'
           endif
+      endif
   else
       print *, 'AO integrals not found in TREXIO file'
   endif
@@ -270,7 +274,8 @@ subroutine run(f)
 
           deallocate(Vi, V, tmp)
           print *, 'Cholesky MO integrals read from TREXIO file'
-      endif
+
+      else
 
         rc = trexio_has_mo_2e_int_eri(f)
         if (rc /= TREXIO_HAS_NOT) then
@@ -310,6 +315,8 @@ subroutine run(f)
             print *, 'MO integrals read from TREXIO file'
         endif
 
+      endif
+
   else
       print *, 'MO integrals not found in TREXIO file'
   endif

src/utils/cgtos_utils.irp.f

@@ -56,7 +56,7 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde
   call multiply_cpoly(P_a(0), a, P_b(0), b, P_new(0), n_new)
   iorder = a + b
 
-end subroutine give_explicit_cpoly_and_cgaussian_x
+end
 
 ! ---
 
@@ -141,7 +141,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder,
   !DIR$ FORCEINLINE
   call multiply_cpoly(P_a(0,3), a(3), P_b(0,3), b(3), P_new(0,3), n_new)
 
-end subroutine give_explicit_cpoly_and_cgaussian
+end
 
 ! ---
 
@@ -249,7 +249,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
   xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv
   xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv
 
-end subroutine cgaussian_product
+end
 
 ! ---
 
@@ -290,7 +290,7 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp)
   k  = zexp(-k)
   xp = (a*xa + b*xb) * p_inv
 
-end subroutine cgaussian_product_x
+end
 
 ! ---
 
@@ -338,7 +338,7 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd)
     exit
   enddo
 
-end subroutine multiply_cpoly
+end
 
 ! ---
 
@@ -373,7 +373,7 @@ subroutine add_cpoly(b, nb, c, nc, d, nd)
     if(nd < 0) exit
   enddo
 
-end subroutine add_cpoly
+end
 
 ! ---
 
@@ -413,7 +413,7 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd)
 
   endif
 
-end subroutine add_cpoly_multiply
+end
 
 ! ---
 
@@ -475,7 +475,7 @@ subroutine recentered_cpoly2(P_A, x_A, x_P, a, P_B, x_B, x_Q, b)
     P_B(i) = binom_func(b,b-i) * pows_b(b-i)
   enddo
 
-end subroutine recentered_cpoly2
+end
 
 ! ---
 
@@ -514,267 +514,7 @@ complex*16 function Fc_integral(n, inv_sq_p)
 
   Fc_integral = sqpi * 0.5d0**n * inv_sq_p**dble(n+1) * fact(n) / fact(shiftr(n, 1))
 
-end function Fc_integral
+end
-
-! ---
-
-complex*16 function crint(n, rho)
-
-  implicit none
-  include 'constants.include.F'
-
-  integer,    intent(in)  :: n
-  complex*16, intent(in)  :: rho
-
-  integer                 :: i, mmax
-  double precision        :: rho_mod, rho_re, rho_im 
-  double precision        :: sq_rho_re, sq_rho_im
-  double precision        :: n_tmp
-  complex*16              :: sq_rho, rho_inv, rho_exp
-
-  complex*16              :: crint_smallz, cpx_erf
-
-  rho_re  = REAL (rho)
-  rho_im  = AIMAG(rho)
-  rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
-
-  if(rho_mod < 10.d0) then
-    ! small z
-
-    if(rho_mod .lt. 1.d-10) then
-      crint = 1.d0 / dble(n + n + 1)
-    else
-      crint = crint_smallz(n, rho)
-    endif
-
-  else
-    ! large z
-
-    if(rho_mod .gt. 40.d0) then
-
-      n_tmp = dble(n) + 0.5d0
-      crint = 0.5d0 * gamma(n_tmp) / (rho**n_tmp)
-
-    else
-
-      ! get \sqrt(rho)
-      sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
-      sq_rho_im = 0.5d0 * rho_im / sq_rho_re 
-      sq_rho    = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
-
-      rho_exp = 0.5d0 * zexp(-rho)
-      rho_inv = (1.d0, 0.d0) / rho
-
-      crint = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho
-      mmax = n
-      if(mmax .gt. 0) then
-        do i = 0, mmax-1
-          crint = ((dble(i) + 0.5d0) * crint - rho_exp) * rho_inv
-        enddo
-      endif
-
-      ! ***
-      
-    endif
-
-  endif
-
-!  print *, n, real(rho), real(crint)
-
-end function crint
-
-! ---
-
-complex*16 function crint_sum(n_pt_out, rho, d1)
-
-  implicit none
-  include 'constants.include.F'
-
-  integer,    intent(in)  :: n_pt_out
-  complex*16, intent(in)  :: rho, d1(0:n_pt_out)
-
-  integer                 :: n, i, mmax
-  double precision        :: rho_mod, rho_re, rho_im 
-  double precision        :: sq_rho_re, sq_rho_im
-  complex*16              :: sq_rho, F0
-  complex*16              :: rho_tmp, rho_inv, rho_exp
-  complex*16, allocatable :: Fm(:)
-
-  complex*16              :: crint_smallz, cpx_erf
-
-  rho_re  = REAL (rho)
-  rho_im  = AIMAG(rho)
-  rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
-
-  if(rho_mod < 10.d0) then
-    ! small z
-
-    if(rho_mod .lt. 1.d-10) then
-
-!      print *, ' 111'
-!      print *, ' rho = ', rho
-
-      crint_sum = d1(0)
-!      print *, 0, 1
-
-      do i = 2, n_pt_out, 2
-
-        n = shiftr(i, 1)
-        crint_sum = crint_sum + d1(i) / dble(n+n+1)
-
-!        print *, n, 1.d0 / dble(n+n+1)
-      enddo
-
-      ! ***
-
-    else
-
-!      print *, ' 222'
-!      print *, ' rho = ', real(rho)
-!      if(abs(aimag(rho)) .gt. 1d-15) then
-!        print *, ' complex rho', rho
-!        stop
-!      endif
-
-      crint_sum = d1(0) * crint_smallz(0, rho)
-
-!      print *, 0, real(d1(0)), real(crint_smallz(0, rho))
-!      if(abs(aimag(d1(0))) .gt. 1d-15) then
-!        print *, ' complex d1(0)', d1(0)
-!        stop
-!      endif
-
-      do i = 2, n_pt_out, 2
-        n = shiftr(i, 1)
-        crint_sum = crint_sum + d1(i) * crint_smallz(n, rho)
-
-!        print *, n, real(d1(i)), real(crint_smallz(n, rho))
-!        if(abs(aimag(d1(i))) .gt. 1d-15) then
-!          print *, ' complex d1(i)', i, d1(i)
-!          stop
-!        endif
-
-      enddo
-
-!      print *, 'sum = ', real(crint_sum)
-!      if(abs(aimag(crint_sum)) .gt. 1d-15) then
-!        print *, ' complex crint_sum', crint_sum
-!        stop
-!      endif
-
-      ! ***
-
-    endif
-
-  else
-    ! large z
-
-    if(rho_mod .gt. 40.d0) then
-
-!      print *, ' 333'
-!      print *, ' rho = ', rho
-
-      rho_inv   = (1.d0, 0.d0) / rho
-      rho_tmp   = 0.5d0 * sqpi * zsqrt(rho_inv)
-      crint_sum = rho_tmp * d1(0)
-!      print *, 0, rho_tmp
-
-      do i = 2, n_pt_out, 2
-        n = shiftr(i, 1)
-        rho_tmp   = rho_tmp * (dble(n) + 0.5d0) * rho_inv
-        crint_sum = crint_sum + rho_tmp * d1(i)
-!        print *, n, rho_tmp
-      enddo
-
-      ! ***
-
-    else
-
-!      print *, ' 444'
-!      print *, ' rho = ', rho
-
-      ! get \sqrt(rho)
-      sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
-      sq_rho_im = 0.5d0 * rho_im / sq_rho_re 
-      sq_rho    = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
-      !sq_rho = zsqrt(rho)
-
-
-      F0        = 0.5d0 * sqpi * cpx_erf(sq_rho_re, sq_rho_im) / sq_rho
-      crint_sum = F0 * d1(0)
-!      print *, 0, F0
-
-      rho_exp = 0.5d0 * zexp(-rho)
-      rho_inv = (1.d0, 0.d0) / rho
-
-      mmax = shiftr(n_pt_out, 1)
-      if(mmax .gt. 0) then
-
-        allocate( Fm(mmax) )
-        Fm(1:mmax) = (0.d0, 0.d0)
-
-        do n = 0, mmax-1
-          F0      = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv
-          Fm(n+1) = F0
-!          print *, n, F0
-        enddo
-    
-        do i = 2, n_pt_out, 2
-          n = shiftr(i, 1)
-          crint_sum = crint_sum + Fm(n) * d1(i)
-        enddo
-        deallocate(Fm)
-      endif
-
-      ! ***
-      
-    endif
-
-  endif
-
-end function crint_sum
-
-! ---
-
-complex*16 function crint_smallz(n, rho)
-
-  BEGIN_DOC
-  ! Standard version of rint
-  END_DOC
-
-  implicit none
-  integer,    intent(in)      :: n
-  complex*16, intent(in)      :: rho
-
-  integer,          parameter :: kmax = 40
-  double precision, parameter :: eps = 1.d-13
-
-  integer                     :: k
-  double precision            :: delta_mod
-  complex*16                  :: rho_k, ct, delta_k
-
-  ct           = 0.5d0 * zexp(-rho) * gamma(dble(n) + 0.5d0)
-  rho_k        = (1.d0, 0.d0)
-  crint_smallz = ct * rho_k / gamma(dble(n) + 1.5d0)
-
-  do k = 1, kmax
-
-    rho_k        = rho_k * rho
-    delta_k      = ct * rho_k / gamma(dble(n+k) + 1.5d0)
-    crint_smallz = crint_smallz + delta_k
-
-    delta_mod = dsqrt(REAL(delta_k)*REAL(delta_k) + AIMAG(delta_k)*AIMAG(delta_k))
-    if(delta_mod .lt. eps) return
-  enddo
-
-  if(delta_mod > eps) then
-    write(*,*) ' pb in crint_smallz !'
-    write(*,*) ' n, rho    = ', n, rho
-    write(*,*) ' delta_mod = ', delta_mod
-    stop 1
-  endif
-
-end function crint_smallz
 
 ! ---
 

src/utils/constants.include.F

@@ -9,6 +9,9 @@ double precision, parameter :: pi_5_2 =  34.9868366552d0
 double precision, parameter :: dfour_pi =  4.d0*dacos(-1.d0)
 double precision, parameter :: dtwo_pi =  2.d0*dacos(-1.d0)
 double precision, parameter :: inv_sq_pi =  1.d0/dsqrt(dacos(-1.d0))
+double precision, parameter :: c_mu_gauss =  27.d0/(8.d0*dsqrt(dacos(-1.d0)))
+double precision, parameter :: c_mu_gauss_tot =  1.5d0*27.d0/(8.d0*dsqrt(dacos(-1.d0)))+3.d0/dsqrt(dacos(-1.d0))
+double precision, parameter :: alpha_mu_gauss =  1.5d0
 double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
 double precision, parameter :: thresh = 1.d-15
 double precision, parameter :: cx_lda = -0.73855876638202234d0

src/utils/cpx_boys.irp.f

@@ -0,0 +1,543 @@
+
+! ---
+
+complex*16 function crint_1(n, rho)
+
+  implicit none
+  include 'constants.include.F'
+
+  integer,    intent(in) :: n
+  complex*16, intent(in) :: rho
+
+  integer                :: i, mmax
+  double precision       :: rho_mod, rho_re, rho_im 
+  double precision       :: sq_rho_re, sq_rho_im
+  double precision       :: n_tmp
+  complex*16             :: sq_rho, rho_inv, rho_exp
+
+  complex*16             :: crint_smallz, cpx_erf_1
+  complex*16             :: cpx_erf_2
+
+  rho_re  = real (rho)
+  rho_im  = aimag(rho)
+  rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
+
+  if(rho_mod < 10.d0) then
+    ! small z
+    if(rho_mod .lt. 1.d-15) then
+      crint_1 = 1.d0 / dble(n + n + 1)
+    else
+      crint_1 = crint_smallz(n, rho)
+    endif
+
+  else
+    ! large z
+
+    if(rho_mod .gt. 40.d0) then
+
+      n_tmp = dble(n) + 0.5d0
+      crint_1 = 0.5d0 * gamma(n_tmp) / (rho**n_tmp)
+
+    else
+
+      ! get \sqrt(rho)
+      sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
+      sq_rho_im = 0.5d0 * rho_im / sq_rho_re 
+      sq_rho    = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
+
+      rho_exp = 0.5d0 * zexp(-rho)
+      rho_inv = (1.d0, 0.d0) / rho
+
+      !print*, sq_rho_re, sq_rho_im
+      !print*, cpx_erf_1(sq_rho_re, sq_rho_im)
+      !print*, cpx_erf_2(sq_rho_re, sq_rho_im)
+
+      crint_1 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho
+      mmax = n
+      if(mmax .gt. 0) then
+        do i = 0, mmax-1
+          crint_1 = ((dble(i) + 0.5d0) * crint_1 - rho_exp) * rho_inv
+        enddo
+      endif
+
+    endif
+
+  endif
+
+end
+
+! ---
+
+complex*16 function crint_quad(n, rho)
+
+  implicit none
+
+  integer,    intent(in)  :: n
+  complex*16, intent(in)  :: rho
+
+  integer                 :: i_quad, n_quad
+  double precision        :: tmp_inv, tmp
+
+  n_quad = 1000000000
+  tmp_inv = 1.d0 / dble(n_quad)
+
+  !crint_quad = 0.5d0 * zexp(-rho)
+  !do i_quad = 1, n_quad - 1
+  !  tmp = tmp_inv * dble(i_quad)
+  !  tmp = tmp * tmp
+  !  crint_quad += zexp(-rho*tmp) * tmp**n
+  !enddo
+  !crint_quad = crint_quad * tmp_inv
+
+  !crint_quad = 0.5d0 * zexp(-rho)
+  !do i_quad = 1, n_quad - 1
+  !  tmp = tmp_inv * dble(i_quad)
+  !  crint_quad += zexp(-rho*tmp) * tmp**n / dsqrt(tmp)
+  !enddo
+  !crint_quad = crint_quad * 0.5d0 * tmp_inv
+
+  ! Composite Boole's Rule
+  crint_quad = 7.d0 * zexp(-rho)
+  do i_quad = 1, n_quad - 1
+    tmp = tmp_inv * dble(i_quad)
+    tmp = tmp * tmp
+    if(modulo(i_quad, 4) .eq. 0) then
+      crint_quad += 14.d0 * zexp(-rho*tmp) * tmp**n
+    else if(modulo(i_quad, 2) .eq. 0) then
+      crint_quad += 12.d0 * zexp(-rho*tmp) * tmp**n
+    else
+      crint_quad += 32.d0 * zexp(-rho*tmp) * tmp**n
+    endif
+  enddo
+  crint_quad = crint_quad * 2.d0 * tmp_inv / 45.d0
+
+  ! Composite Simpson's 3/8 rule
+  !crint_quad = zexp(-rho)
+  !do i_quad = 1, n_quad - 1
+  !  tmp = tmp_inv * dble(i_quad)
+  !  tmp = tmp * tmp
+  !  if(modulo(i_quad, 3) .eq. 0) then
+  !    crint_quad += 2.d0 * zexp(-rho*tmp) * tmp**n
+  !  else
+  !    crint_quad += 3.d0 * zexp(-rho*tmp) * tmp**n
+  !  endif
+  !enddo
+  !crint_quad = crint_quad * 3.d0 * tmp_inv / 8.d0
+
+end
+
+! ---
+
+complex*16 function crint_sum_1(n_pt_out, rho, d1)
+
+  implicit none
+  include 'constants.include.F'
+
+  integer,    intent(in)  :: n_pt_out
+  complex*16, intent(in)  :: rho, d1(0:n_pt_out)
+
+  integer                 :: n, i, mmax
+  double precision        :: rho_mod, rho_re, rho_im 
+  double precision        :: sq_rho_re, sq_rho_im
+  complex*16              :: sq_rho, F0
+  complex*16              :: rho_tmp, rho_inv, rho_exp
+  complex*16, allocatable :: Fm(:)
+
+  complex*16              :: crint_smallz, cpx_erf_1
+
+
+  rho_re  = real (rho)
+  rho_im  = aimag(rho)
+  rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
+
+!  ! debug
+!  double precision        :: d1_real(0:n_pt_out)
+!  double precision        :: rint_sum
+!  do i = 0, n_pt_out
+!    d1_real(i) = real(d1(i))
+!  enddo
+!  crint_sum_1 = rint_sum(n_pt_out, rho_re, d1_real)
+!  return
+
+  if(rho_mod < 10.d0) then
+    ! small z
+
+    if(rho_mod .lt. 1.d-15) then
+
+      crint_sum_1 = d1(0)
+      do i = 2, n_pt_out, 2
+        n = shiftr(i, 1)
+        crint_sum_1 = crint_sum_1 + d1(i) / dble(n+n+1)
+      enddo
+
+    else
+
+      crint_sum_1 = d1(0) * crint_smallz(0, rho)
+      do i = 2, n_pt_out, 2
+        n = shiftr(i, 1)
+        crint_sum_1 = crint_sum_1 + d1(i) * crint_smallz(n, rho)
+      enddo
+
+    endif
+
+  else
+    ! large z
+
+    if(rho_mod .gt. 40.d0) then
+
+      rho_inv = (1.d0, 0.d0) / rho
+      rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv)
+
+      crint_sum_1 = rho_tmp * d1(0)
+      do i = 2, n_pt_out, 2
+        n = shiftr(i, 1)
+        rho_tmp   = rho_tmp * (dble(n) + 0.5d0) * rho_inv
+        crint_sum_1 = crint_sum_1 + rho_tmp * d1(i)
+      enddo
+
+    else
+
+      ! get \sqrt(rho)
+      sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
+      sq_rho_im = 0.5d0 * rho_im / sq_rho_re 
+      sq_rho    = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
+
+      F0        = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho
+      crint_sum_1 = F0 * d1(0)
+
+      rho_exp = 0.5d0 * zexp(-rho)
+      rho_inv = (1.d0, 0.d0) / rho
+
+      mmax = shiftr(n_pt_out, 1)
+      if(mmax .gt. 0) then
+
+        allocate(Fm(mmax))
+        Fm(1:mmax) = (0.d0, 0.d0)
+
+        do n = 0, mmax-1
+          F0      = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv
+          Fm(n+1) = F0
+        enddo
+    
+        do i = 2, n_pt_out, 2
+          n = shiftr(i, 1)
+          crint_sum_1 = crint_sum_1 + Fm(n) * d1(i)
+        enddo
+
+        deallocate(Fm)
+      endif
+
+    endif ! rho_mod
+  endif ! rho_mod
+
+end
+
+! ---
+
+complex*16 function crint_smallz(n, rho)
+
+  BEGIN_DOC
+  ! Standard version of rint
+  END_DOC
+
+  implicit none
+  integer,    intent(in)      :: n
+  complex*16, intent(in)      :: rho
+
+  integer,          parameter :: kmax = 40
+  double precision, parameter :: eps = 1.d-13
+
+  integer                     :: k
+  double precision            :: delta_mod
+  complex*16                  :: rho_k, ct, delta_k
+
+  ct           = 0.5d0 * zexp(-rho) * gamma(dble(n) + 0.5d0)
+  crint_smallz = ct / gamma(dble(n) + 1.5d0)
+
+  rho_k = (1.d0, 0.d0)
+  do k = 1, kmax
+
+    rho_k        = rho_k * rho
+    delta_k      = ct * rho_k / gamma(dble(n+k) + 1.5d0)
+    crint_smallz = crint_smallz + delta_k
+
+    delta_mod = dsqrt(real(delta_k)*real(delta_k) + aimag(delta_k)*aimag(delta_k))
+    if(delta_mod .lt. eps) return
+  enddo
+
+  if(delta_mod > eps) then
+    write(*,*) ' pb in crint_smallz !'
+    write(*,*) ' n, rho    = ', n, rho
+    write(*,*) ' delta_mod = ', delta_mod
+    !stop 1
+  endif
+
+end
+
+! ---
+
+complex*16 function crint_2(n, rho)
+
+  implicit none
+
+  integer,    intent(in) :: n
+  complex*16, intent(in) :: rho
+
+  double precision :: tmp
+  complex*16 :: rho2
+  complex*16 :: vals(0:n)
+  complex*16, external :: crint_smallz
+
+  if(abs(rho) < 10.d0) then
+
+    if(abs(rho) .lt. 1d-6) then
+      tmp = 2.d0 * dble(n)
+      rho2 = rho * rho
+      crint_2 = 1.d0 / (tmp + 1.d0)         &
+              - rho / (tmp + 3.d0)          &
+              + 0.5d0 * rho2 / (tmp + 5.d0) &
+              - 0.16666666666666666d0 * rho * rho2 / (tmp + 7.d0)
+    else
+      crint_2 = crint_smallz(n, rho)
+    endif
+
+  else
+
+    if(real(rho) .ge. 0.d0) then
+      call zboysfun(n, rho, vals)
+      crint_2 = vals(n)
+    else
+      call zboysfunnrp(n, rho, vals)
+      crint_2 = vals(n) * zexp(-rho)
+    endif
+
+  endif
+
+  return
+end
+
+! ---
+
+subroutine zboysfun(n_max, x, vals)
+
+  BEGIN_DOC
+  !
+  ! Computes values of the Boys function for n = 0, 1, ..., n_max
+  ! for a complex valued argument
+  !
+  ! Input: x --- argument, complex*16, Re(x) >= 0
+  ! Output: vals  --- values of the Boys function, n = 0, 1, ..., n_max
+  !
+  END_DOC
+
+  implicit none
+
+  integer,    intent(in)  :: n_max
+  complex*16, intent(in)  :: x
+  complex*16, intent(out) :: vals(0:n_max)
+
+  integer                 :: n
+  complex*16              :: yy, x_inv
+
+  call zboysfun00(x, vals(0))
+
+  yy = 0.5d0 * zexp(-x)
+  x_inv = (1.d0, 0.d0) / x
+  do n = 1, n_max
+    vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - yy) * x_inv
+  enddo
+
+  return
+end
+
+! ---
+
+subroutine zboysfunnrp(n_max, x, vals)
+
+  BEGIN_DOC
+  !
+  ! Computes values of e^z F(n,z) for n = 0, 1, ..., n_max
+  ! (where F(n,z) are the Boys functions)
+  ! for a complex valued argument WITH NEGATIVE REAL PART
+  !
+  ! Input: x  --- argument, complex *16 Re(x)<=0
+  ! Output: vals  --- values of e^z F(n,z), n = 0, 1, ..., n_max
+  !
+  END_DOC
+
+  implicit none
+
+  integer,    intent(in)  :: n_max
+  complex*16, intent(in)  :: x
+  complex*16, intent(out) :: vals(0:n_max)
+
+  integer                 :: n
+  complex*16              :: x_inv
+
+  call zboysfun00nrp(x, vals(0))
+
+  x_inv = (1.d0, 0.d0) / x
+  do n = 1, n_max
+    vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - 0.5d0) * x_inv
+  enddo
+
+  return
+end
+
+! ---
+
+complex*16 function crint_sum_2(n_pt_out, rho, d1)
+
+  implicit none
+
+  integer,    intent(in)  :: n_pt_out
+  complex*16, intent(in)  :: rho, d1(0:n_pt_out)
+                          
+  integer                 :: n, i
+  integer                 :: n_max
+
+  complex*16, allocatable :: vals(:)
+
+  !complex*16, external    :: crint_2
+  !crint_sum_2 = (0.d0, 0.d0)
+  !do i = 0, n_pt_out, 2
+  !  n = shiftr(i, 1)
+  !  crint_sum_2 = crint_sum_2 + d1(i) * crint_2(n, rho)
+  !enddo
+
+  n_max = shiftr(n_pt_out, 1)
+
+  allocate(vals(0:n_max))
+  call crint_2_vec(n_max, rho, vals)
+
+  crint_sum_2 = d1(0) * vals(0)
+  do i = 2, n_pt_out, 2
+    n = shiftr(i, 1)
+    crint_sum_2 += d1(i) * vals(n)
+  enddo
+
+  deallocate(vals)
+
+  return
+end
+
+! ---
+
+subroutine crint_2_vec(n_max, rho, vals)
+
+  implicit none
+
+  integer,    intent(in)  :: n_max
+  complex*16, intent(in)  :: rho
+  complex*16, intent(out) :: vals(0:n_max)
+
+  integer                 :: n
+  double precision        :: tmp, abs_rho
+  complex*16              :: rho2, rho3, erho
+
+
+  abs_rho = abs(rho)
+
+  if(abs_rho < 10.d0) then
+
+    if(abs_rho .lt. 1d-6) then
+
+      ! use finite expansion for very small rho
+
+      ! rho^2 / 2
+      rho2 = 0.5d0 * rho * rho
+      ! rho^3 / 6
+      rho3 = 0.3333333333333333d0 * rho * rho2
+
+      vals(0) = 1.d0 - 0.3333333333333333d0 * rho + 0.2d0 * rho2 - 0.14285714285714285d0 * rho3
+      do n = 1, n_max
+        tmp = 2.d0 * dble(n)
+        vals(n) = 1.d0 / (tmp + 1.d0) - rho  / (tmp + 3.d0) &
+                + rho2 / (tmp + 5.d0) - rho3 / (tmp + 7.d0)
+      enddo
+
+    else
+
+      call crint_smallz_vec(n_max, rho, vals)
+
+    endif
+
+  else
+
+    if(real(rho) .ge. 0.d0) then
+
+      call zboysfun(n_max, rho, vals)
+
+    else
+
+      call zboysfunnrp(n_max, rho, vals)
+      erho = zexp(-rho)
+      do n = 0, n_max
+        vals(n) = vals(n) * erho
+      enddo
+
+    endif
+
+  endif
+
+  return
+end
+
+! ---
+
+subroutine crint_smallz_vec(n_max, rho, vals)
+
+  BEGIN_DOC
+  ! Standard version of rint
+  END_DOC
+
+  implicit none
+  integer,    intent(in)      :: n_max
+  complex*16, intent(in)      :: rho
+  complex*16, intent(out)     :: vals(0:n_max)
+
+  integer,          parameter :: kmax = 40
+  double precision, parameter :: eps = 1.d-13
+
+  integer                     :: k, n
+  complex*16                  :: ct, delta_k
+  complex*16                  :: rhoe
+  complex*16, allocatable     :: rho_k(:)
+
+
+  allocate(rho_k(0:kmax))
+
+  rho_k(0) = (1.d0, 0.d0)
+  do k = 1, kmax
+    rho_k(k) = rho_k(k-1) * rho
+  enddo
+
+  rhoe = 0.5d0 * zexp(-rho)
+
+  do n = 0, n_max
+
+    ct = rhoe * gamma(dble(n) + 0.5d0)
+    vals(n) = ct / gamma(dble(n) + 1.5d0)
+  
+    do k = 1, kmax
+      delta_k = ct * rho_k(k) / gamma(dble(n+k) + 1.5d0)
+      vals(n) += delta_k
+      if(abs(delta_k) .lt. eps) then
+        exit
+      endif
+    enddo
+  
+    !if(abs(delta_k) > eps) then
+    !  write(*,*) ' pb in crint_smallz_vec !'
+    !  write(*,*) ' n, rho    = ', n, rho
+    !  write(*,*) ' |delta_k| = ', abs(delta_k)
+    !endif
+  enddo
+
+  deallocate(rho_k)
+
+  return
+end
+
+! ---
+

src/utils/cpx_erf.irp.f

@@ -1,7 +1,7 @@
 
 ! ---
 
-complex*16 function cpx_erf(x, y)
+complex*16 function cpx_erf_1(x, y)
 
   BEGIN_DOC
   !
@@ -25,7 +25,7 @@ complex*16 function cpx_erf(x, y)
 
   if(yabs .lt. 1.d-15) then
 
-    cpx_erf = (1.d0, 0.d0) * derf(x)
+    cpx_erf_1 = (1.d0, 0.d0) * derf(x)
     return
 
   else
@@ -38,12 +38,12 @@ complex*16 function cpx_erf(x, y)
   endif
 
   if(y .gt. 0.d0) then
-    cpx_erf = erf_tot
+    cpx_erf_1 = erf_tot
   else
-    cpx_erf = CONJG(erf_tot)
+    cpx_erf_1 = conjg(erf_tot)
   endif
 
-end function cpx_erf
+end
 
 ! ---
 
@@ -70,7 +70,7 @@ complex*16 function erf_E(x, yabs)
 
   endif
 
-end function erf_E
+end
 
 ! ---
 
@@ -109,7 +109,7 @@ double precision function erf_F(x, yabs)
 
   endif 
 
-end function erf_F
+end
 
 ! ---
 
@@ -149,7 +149,7 @@ complex*16 function erf_G(x, yabs)
 
   enddo
 
-end function erf_G
+end
 
 ! ---
 
@@ -186,7 +186,7 @@ complex*16 function erf_H(x, yabs)
       tmp2  = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1
       erf_H = erf_H + tmp2
 
-      tmp_mod = dsqrt(REAL(tmp2)*REAL(tmp2) + AIMAG(tmp2)*AIMAG(tmp2))
+      tmp_mod = dsqrt(real(tmp2)*real(tmp2) + aimag(tmp2)*aimag(tmp2))
       if(tmp_mod .lt. 1d-15) exit
     enddo
     erf_H = ct * erf_H
@@ -197,8 +197,394 @@ complex*16 function erf_H(x, yabs)
 
   endif 
 
-end function erf_H
+end
 
 ! ---
 
+complex*16 function cpx_erf_2(x, y)
+
+  BEGIN_DOC
+  !
+  ! compute erf(z) for z = x + i y
+  !
+  ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021)
+  ! https://doi.org/10.1063/5.0062444
+  !
+  END_DOC
+
+  implicit none
+  
+  double precision, intent(in) :: x, y
+
+  double precision             :: yabs
+  complex*16                   :: z
+
+  yabs = dabs(y)
+
+  if(yabs .lt. 1.d-15) then
+
+    cpx_erf_2 = (1.d0, 0.d0) * derf(x)
+    return
+
+  else
+
+    z = x + (0.d0, 1.d0) * y
+
+    if(x .ge. 0.d0) then
+      call zboysfun00(z, cpx_erf_2)
+    else
+      call zboysfun00nrp(z, cpx_erf_2)
+      cpx_erf_2 = cpx_erf_2 * zexp(-z)
+    endif
+
+  endif
+
+  return
+end
+
+! ---
+
+subroutine zboysfun00(z, val)
+
+  BEGIN_DOC
+  !
+  ! Computes values of the Boys function for n=0 
+  ! for a complex valued argument
+  !
+  ! Input: z  --- argument, complex*16, Real(z) >= 0
+  ! Output: val  --- value of the Boys function n=0
+  !
+  ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021)
+  ! https://doi.org/10.1063/5.0062444
+  !
+  END_DOC
+
+  implicit none
+
+  double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, &
+                                                     0.249999999999993161d0, &
+                                                    -0.374999999999766599d0, &
+                                                     0.937499999992027020d0, &
+                                                    -3.28124999972738868d0,  &
+                                                     14.7656249906697030d0,  &
+                                                    -81.2109371803307752d0 /)
+
+  double precision, parameter :: taylcoef(0:10) = (/ 1.0d0,                    &
+                                                    -0.333333333333333333d0,   &
+                                                     0.1d0,                    &
+                                                    -0.238095238095238095d-01, &
+                                                     0.462962962962962963d-02, &
+                                                    -0.757575757575757576d-03, &
+                                                     0.106837606837606838d-03, & 
+                                                    -0.132275132275132275d-04, &
+                                                     1.458916900093370682d-06, &
+                                                    -1.450385222315046877d-07, &
+                                                     1.3122532963802805073d-08 /)
+
+  double precision, parameter :: sqpio2 = 0.886226925452758014d0
+
+  double precision, parameter :: pp(1:22) = (/ 0.001477878263796956477d0, &
+                                               0.013317276413725817441d0, &
+                                               0.037063591452052541530d0, &
+                                               0.072752512422882761543d0, &
+                                               0.120236941228785688896d0, &
+                                               0.179574293958937717967d0, &
+                                               0.253534046984087292596d0, &
+                                               0.350388652780721927513d0, &
+                                               0.482109575931276669313d0, &
+                                               0.663028993158374107103d0, &
+                                               0.911814736856590885929d0, &
+                                               1.2539502287919293d0,      &
+                                               1.7244634233573395d0,      &
+                                               2.3715248262781863d0,      &
+                                               3.2613796996078355d0,      &
+                                               4.485130169059591d0,       &
+                                               6.168062135122484d0,       &
+                                               8.48247187231787d0,        &
+                                               11.665305486296793d0,      &
+                                               16.042417132288328d0,      &
+                                               22.06192951814709d0,       &
+                                               30.340112094708307d0 /)
+
+  double precision, parameter :: ff(1:22) = (/ 0.0866431027201416556d0, &
+                                               0.0857720608434394764d0, &
+                                               0.0839350436829178814d0, &
+                                               0.0809661970413229146d0, &
+                                               0.0769089548492978618d0, &
+                                               0.0731552078711821626d0, &
+                                               0.0726950035163157228d0, &
+                                               0.0752842556089304050d0, &
+                                               0.0770943953645196145d0, &
+                                               0.0754250625677530441d0, &
+                                               0.0689686192650315305d0, &
+                                               0.05744480422143023d0,   &
+                                               0.04208199434694545d0,   &
+                                               0.025838539448223282d0,  &
+                                               0.012445024157255563d0,  &
+                                               0.004292541592599837d0,  &
+                                               0.0009354342987735969d0, &
+                                               0.10840885466502504d-03, &
+                                               5.271867966761674d-06,   &
+                                               7.765974039750418d-08,   &
+                                               2.2138172422680093d-10,  &
+                                               6.594161760037707d-14 /)
+
+  complex*16, intent(in)  :: z
+  complex*16, intent(out) :: val
+
+  integer    :: k
+  complex*16 :: z1, zz, y
+
+  zz = zexp(-z)
+
+  if(abs(z) .ge. 100.0d0) then
+
+    ! large |z|
+    z1 = 1.0d0 / zsqrt(z)
+    y  = 1.0d0 / z
+    val = asymcoef(7)
+    do k = 6, 1, -1
+      val = val * y + asymcoef(k)
+    enddo
+    val = zz * val * y + z1 * sqpio2
+
+  else if(abs(z) .le. 0.35d0) then
+
+    ! small |z|
+    val = taylcoef(10) * (1.d0, 0.d0)
+    do k = 9, 0, -1
+      val = val * z + taylcoef(k)
+    enddo
+
+  else
+
+    ! intermediate |z|
+    val = sqpio2 / zsqrt(z) - 0.5d0 * zz * sum(ff(1:22)/(z+pp(1:22)))
+    !val = (0.d0, 0.d0)
+    !do k = 1, 22
+    !  val += ff(k) / (z + pp(k))
+    !enddo
+    !val = sqpio2 / zsqrt(z) - 0.5d0 * zz * val
+
+  endif
+
+  return
+end
+
+! ---
+
+subroutine zboysfun00nrp(z, val)
+
+  BEGIN_DOC
+  !
+  ! Computes values of the exp(z) F(0,z)
+  ! (where F(0,z) is the Boys function)
+  ! for a complex valued argument with Real(z)<=0
+  !
+  ! Input: z  --- argument, complex*16, !!! Real(z)<=0 !!!
+  ! Output: val  --- value of the function !!! exp(z) F(0,z) !!!, where F(0,z) is the Boys function
+  !
+  ! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021)
+  ! https://doi.org/10.1063/5.0062444
+  !
+  END_DOC
+
+  implicit none
+
+  double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, &
+                                                     0.249999999999993161d0, &
+                                                    -0.374999999999766599d0, &
+                                                     0.937499999992027020d0, &
+                                                    -3.28124999972738868d0,  &
+                                                     14.7656249906697030d0,  &
+                                                    -81.2109371803307752d0 /)
+
+  double precision, parameter :: taylcoef(0:10) = (/ 1.0d0,                    &
+                                                    -0.333333333333333333d0,   &
+                                                     0.1d0,                    &
+                                                    -0.238095238095238095d-01, &
+                                                     0.462962962962962963d-02, &
+                                                    -0.757575757575757576d-03, &
+                                                     0.106837606837606838d-03, & 
+                                                    -0.132275132275132275d-04, &
+                                                     1.458916900093370682d-06, &
+                                                    -1.450385222315046877d-07, &
+                                                     1.3122532963802805073d-08 /)
+
+  double precision, parameter :: tol     = 1.0d-03
+  double precision, parameter :: sqpio2  = 0.886226925452758014d0 ! sqrt(pi)/2
+  double precision, parameter :: pi      = 3.14159265358979324d0
+  double precision, parameter :: etmax   = 25.7903399171930621d0
+  double precision, parameter :: etmax1  = 26.7903399171930621d0
+  complex*16, parameter :: ima = (0.d0, 1.d0)
+
+  double precision, parameter :: pp(1:16) = (/ 0.005299532504175031d0, &
+                                               0.0277124884633837d0,   &
+                                               0.06718439880608407d0,  &
+                                               0.12229779582249845d0,  &
+                                               0.19106187779867811d0,  &
+                                               0.27099161117138637d0,  &
+                                               0.35919822461037054d0,  &
+                                               0.45249374508118123d0,  &
+                                               0.5475062549188188d0,   &
+                                               0.6408017753896295d0,   &
+                                               0.7290083888286136d0,   &
+                                               0.8089381222013219d0,   &
+                                               0.8777022041775016d0,   &
+                                               0.9328156011939159d0,   &
+                                               0.9722875115366163d0,   &
+                                               0.994700467495825d0 /)
+
+  double precision, parameter :: ww(1:16) = (/ 0.013576229705876844d0, &
+                                               0.03112676196932382d0,  &
+                                               0.04757925584124612d0,  &
+                                               0.062314485627766904d0, &
+                                               0.07479799440828848d0,  &
+                                               0.08457825969750153d0,  &
+                                               0.09130170752246194d0,  &
+                                               0.0947253052275344d0,   &
+                                               0.0947253052275344d0,   &
+                                               0.09130170752246194d0,  &
+                                               0.08457825969750153d0,  &
+                                               0.07479799440828848d0,  &
+                                               0.062314485627766904d0, &
+                                               0.04757925584124612d0,  &
+                                               0.03112676196932382d0,  &
+                                               0.013576229705876844d0 /)
+
+  double precision, parameter :: qq (1:16) = (/ 0.0007243228510223928d0, &
+                                                0.01980651726441906d0,   &
+                                                0.11641097769229371d0,   &
+                                                0.38573968881461146d0,   &
+                                                0.9414671037609641d0,    &
+                                                1.8939510935716377d0,    &
+                                                3.3275564293459383d0,    &
+                                                5.280587297262129d0,     &
+                                                7.730992222360452d0,     &
+                                                10.590207725831563d0,    &
+                                                13.706359477128965d0,    &
+                                                16.876705473663804d0,    &
+                                                19.867876155236257d0,    &
+                                                22.441333930203022d0,    &
+                                                24.380717439613566d0,    &
+                                                25.51771075067431d0 /)
+
+
+  double precision, parameter :: qq1 (1:16) = (/ 0.0007524078957852004d0,&
+                                                 0.020574499281252233d0, &
+                                                 0.12092472113522865d0,  &
+                                                 0.40069643967765295d0,  &
+                                                 0.9779717449089211d0,   &
+                                                 1.9673875468969015d0,   &
+                                                 3.4565797939091802d0,   &
+                                                 5.485337886599723d0,    &
+                                                 8.030755321535683d0,    &
+                                                 11.000834641174064d0,   &
+                                                 14.237812708111456d0,   &
+                                                 17.531086359214406d0,   &
+                                                 20.6382373144543d0,     &
+                                                 23.31147887603379d0,    &
+                                                 25.326060444703632d0,   &
+                                                 26.507139770710722d0 /)
+
+  double precision, parameter :: uu(1:16) = (/ 0.9992759394074501d0,      &
+                                               0.9803883431758104d0,      &
+                                               0.8901093330366746d0,      &
+                                               0.6799475005849274d0,      &
+                                               0.3900551639790145d0,      &
+                                               0.15047608763371934d0,     &
+                                               0.0358806749968974d0,      &
+                                               0.005089440900100864d0,    &
+                                               0.00043900830706867264d0,  &
+                                               0.000025161192619824898d0, &
+                                               1.1153308427285078d-6,     &
+                                               4.68317018372038d-8,       &
+                                               2.3522908467181876d-9,     &
+                                               1.7941242138648815d-10,    &
+                                               2.5798173021885247d-11,    &
+                                               8.27559122014575d-12 /)
+
+
+  double precision, parameter :: uu1(1:16) = (/ 0.999247875092057d0,       &
+                                                0.979635711599488d0,       &
+                                                0.8861006617341018d0,      &
+                                                0.6698533710831932d0,      &
+                                                0.3760730980014839d0,      &
+                                                0.13982165701683388d0,     &
+                                                0.031537442321301304d0,    &
+                                                0.004147133581658446d0,    &
+                                                0.0003253024081883165d0,   &
+                                                0.000016687766678889653d0, &
+                                                6.555359391864376d-7,      &
+                                                2.4341421258295026d-8,     &
+                                                1.0887481200652014d-9,     &
+                                                7.51542178140961d-11,      &
+                                                1.002378402152542d-11,     &
+                                                3.0767730761654096d-12 /)
+
+  complex*16, intent(in)  :: z
+  complex*16, intent(out) :: val
+
+  integer    :: k
+  complex*16 :: z1, zz, y, zsum, tmp, zt, q, p
+
+  zz = zexp(z)
+
+  if(abs(z) .ge. 100.0d0) then
+    ! large |z|
+    z1 = 1.0d0 / zsqrt(z)
+    y  = 1.0d0 / z
+    val = asymcoef(7)
+    do k = 6, 1, -1
+      val = val * y + asymcoef(k)
+    enddo
+    val = val * y + z1 * sqpio2 * zz
+    return
+  endif
+
+  if(abs(z) .le. 0.35d0) then
+    ! small |z|
+    val = taylcoef(10) * (1.d0, 0.d0)
+    do k = 9, 0, -1
+      val = val * z + taylcoef(k)
+    enddo
+    val = val * zz
+    return
+  endif
+
+  if(abs(etmax+z) .ge. 0.5d0) then
+    ! intermediate |z|
+    zsum = (0.d0, 0.d0)
+    do k = 1, 16
+      if(abs(z + qq(k)) .ge. tol) then 
+        zsum = zsum + ww(k) * (zz - uu(k)) / (qq(k) + z)
+      else
+        q = z + qq(k)
+        p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0
+        zsum = zsum + ww(k) * p *zz
+      endif
+    enddo
+    zt = ima * sqrt(z / etmax)
+    tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt))
+    val = sqrt(etmax) * zsum / sqrt(pi) + zz * tmp / sqrt(pi*z)
+  else
+    zsum = (0.d0, 0.d0)
+    do k = 1, 16
+      if(abs(z + qq1(k)) .ge. tol) then 
+        zsum = zsum + ww(k) * (zz - uu1(k)) / (qq1(k) + z)
+      else
+        q = z + qq1(k)
+        p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0
+        zsum = zsum + ww(k) * p * zz
+      endif
+    enddo
+    zt = ima * zsqrt(z / etmax1)
+    tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt))
+    val = dsqrt(etmax1) * zsum / dsqrt(pi) + zz * tmp / zsqrt(pi*z)
+  endif
+
+  return
+end
+
+! ---
 

src/utils/mmap.f90

@@ -2,6 +2,34 @@ module mmap_module
 
   use iso_c_binding
 
+  type mmap_type
+    type(c_ptr)     :: ptr          ! Pointer to the data
+    character*(128) :: filename     ! Name of the file
+    integer*8       :: length       ! Size of the array in bytes
+    integer         :: fd           ! File descriptor
+
+    ! Pointers to data
+    integer, pointer  :: i1(:)
+    integer, pointer  :: i2(:,:)
+    integer, pointer  :: i3(:,:,:)
+    integer, pointer  :: i4(:,:,:,:)
+
+    integer*8, pointer :: i81(:)
+    integer*8, pointer :: i82(:,:)
+    integer*8, pointer :: i83(:,:,:)
+    integer*8, pointer :: i84(:,:,:,:)
+
+    double precision, pointer :: d1(:)
+    double precision, pointer :: d2(:,:)
+    double precision, pointer :: d3(:,:,:)
+    double precision, pointer :: d4(:,:,:,:)
+
+    real, pointer :: s1(:)
+    real, pointer :: s2(:,:)
+    real, pointer :: s3(:,:,:)
+    real, pointer :: s4(:,:,:,:)
+  end type mmap_type
+
   interface
 
     ! File descriptors
@@ -106,6 +134,200 @@ module mmap_module
       call c_msync_fortran( length, fd_, map)
   end subroutine
 
+
+  ! Functions for the mmap_type
+
+  subroutine mmap_create(filename, shape, bytes, read_only, single_node, map)
+      implicit none
+      character*(*), intent(in)      :: filename   ! Name of the mapped file
+      integer*8, intent(in)          :: shape(:)   ! Shape of the array to map
+      integer, intent(in)            :: bytes      ! Number of bytes per element
+      logical, intent(in)            :: read_only  ! If true, mmap is read-only
+      logical, intent(in)            :: single_node! If true, mmap is on a single node
+      type(mmap_type), intent(out)   :: map        ! mmap
+
+      integer :: i
+      logical :: temporary
+
+      temporary = ( trim(filename) == '' )
+
+      if (.not.temporary) then
+        map%filename = filename
+      else
+        call getenv('EZFIO_FILE', map%filename)
+        map%filename = trim(map%filename) // '/work/tmpfile'
+      endif
+
+      map%length = int(bytes,8)
+      do i=1,size(shape)
+        map%length = map%length * shape(i)
+      enddo
+      call mmap(map%filename, &
+                shape,        &
+                bytes,        &
+                map%fd,       &
+                read_only,    &
+                single_node,  &
+                map%ptr)
+
+      if (temporary) then
+        ! Deleting the file while it is open makes the file invisible on the filesystem,
+        ! and automatically deleted, even if the program crashes
+        open(UNIT=47, FILE=trim(map%filename), STATUS='OLD')
+        close(47,STATUS='DELETE')
+      endif
+
+      map%d1 => NULL()
+      map%d2 => NULL()
+      map%d3 => NULL()
+      map%d4 => NULL()
+      map%s1 => NULL()
+      map%s2 => NULL()
+      map%s3 => NULL()
+      map%s4 => NULL()
+      map%i1 => NULL()
+      map%i2 => NULL()
+      map%i3 => NULL()
+      map%i4 => NULL()
+      map%i81 => NULL()
+      map%i82 => NULL()
+      map%i83 => NULL()
+      map%i84 => NULL()
+
+  end
+
+  subroutine mmap_create_d(filename, shape, read_only, single_node, map)
+      implicit none
+      character*(*), intent(in)      :: filename   ! Name of the mapped file
+      integer*8, intent(in)          :: shape(:)   ! Shape of the array to map
+      logical, intent(in)            :: read_only  ! If true, mmap is read-only
+      logical, intent(in)            :: single_node! If true, mmap is on a single node
+      type(mmap_type), intent(out)   :: map        ! mmap
+
+      call mmap_create(filename, shape, 8, read_only, single_node, map)
+
+      select case (size(shape))
+        case (1)
+          call c_f_pointer(map%ptr, map%d1, shape)
+        case (2)
+          call c_f_pointer(map%ptr, map%d2, shape)
+        case (3)
+          call c_f_pointer(map%ptr, map%d3, shape)
+        case (4)
+          call c_f_pointer(map%ptr, map%d4, shape)
+        case default
+          stop 'mmap: dimension not implemented'
+      end select
+  end subroutine
+
+  subroutine mmap_create_s(filename, shape, read_only, single_node, map)
+      implicit none
+      character*(*), intent(in)      :: filename   ! Name of the mapped file
+      integer*8, intent(in)          :: shape(:)   ! Shape of the array to map
+      logical, intent(in)            :: read_only  ! If true, mmap is read-only
+      logical, intent(in)            :: single_node! If true, mmap is on a single node
+      type(mmap_type), intent(out)   :: map        ! mmap
+
+      call mmap_create(filename, shape, 4, read_only, single_node, map)
+
+      select case (size(shape))
+        case (1)
+          call c_f_pointer(map%ptr, map%s1, shape)
+        case (2)
+          call c_f_pointer(map%ptr, map%s2, shape)
+        case (3)
+          call c_f_pointer(map%ptr, map%s3, shape)
+        case (4)
+          call c_f_pointer(map%ptr, map%s4, shape)
+        case default
+          stop 'mmap: dimension not implemented'
+      end select
+  end subroutine
+
+  subroutine mmap_create_i(filename, shape, read_only, single_node, map)
+      implicit none
+      character*(*), intent(in)      :: filename   ! Name of the mapped file
+      integer*8, intent(in)          :: shape(:)   ! Shape of the array to map
+      logical, intent(in)            :: read_only  ! If true, mmap is read-only
+      logical, intent(in)            :: single_node! If true, mmap is on a single node
+      type(mmap_type), intent(out)   :: map        ! mmap
+
+      call mmap_create(filename, shape, 4, read_only, single_node, map)
+
+      select case (size(shape))
+        case (1)
+          call c_f_pointer(map%ptr, map%i1, shape)
+        case (2)
+          call c_f_pointer(map%ptr, map%i2, shape)
+        case (3)
+          call c_f_pointer(map%ptr, map%i3, shape)
+        case (4)
+          call c_f_pointer(map%ptr, map%i4, shape)
+        case default
+          stop 'mmap: dimension not implemented'
+      end select
+  end subroutine
+
+  subroutine mmap_create_i8(filename, shape, read_only, single_node, map)
+      implicit none
+      character*(*), intent(in)      :: filename   ! Name of the mapped file
+      integer*8, intent(in)          :: shape(:)   ! Shape of the array to map
+      logical, intent(in)            :: read_only  ! If true, mmap is read-only
+      logical, intent(in)            :: single_node! If true, mmap is on a single node
+      type(mmap_type), intent(out)   :: map        ! mmap
+
+      call mmap_create(filename, shape, 8, read_only, single_node, map)
+
+      select case (size(shape))
+        case (1)
+          call c_f_pointer(map%ptr, map%i81, shape)
+        case (2)
+          call c_f_pointer(map%ptr, map%i82, shape)
+        case (3)
+          call c_f_pointer(map%ptr, map%i83, shape)
+        case (4)
+          call c_f_pointer(map%ptr, map%i84, shape)
+        case default
+          stop 'mmap: dimension not implemented'
+      end select
+  end subroutine
+
+  subroutine mmap_destroy(map)
+      implicit none
+      type(mmap_type), intent(inout)    :: map
+
+      call c_munmap_fortran(map%length, map%fd, map%ptr)
+
+      map%ptr      = C_NULL_PTR
+      map%filename = ''
+      map%length   = 0
+      map%fd       = 0
+      map%s1 => NULL()
+      map%s2 => NULL()
+      map%s3 => NULL()
+      map%s4 => NULL()
+      map%d1 => NULL()
+      map%d2 => NULL()
+      map%d3 => NULL()
+      map%d4 => NULL()
+      map%i1 => NULL()
+      map%i2 => NULL()
+      map%i3 => NULL()
+      map%i4 => NULL()
+      map%i81 => NULL()
+      map%i82 => NULL()
+      map%i83 => NULL()
+      map%i84 => NULL()
+  end subroutine
+
+
+  subroutine mmap_sync(map)
+      implicit none
+      type(mmap_type), intent(inout)    :: map
+
+      call c_msync_fortran(map%length, map%fd, map%ptr)
+  end subroutine
+
 end module mmap_module
 
 

src/utils_cc/diis.irp.f

@@ -53,10 +53,10 @@ subroutine diis_cc(all_err,all_t,sze,m,iter,t)
   !$OMP END PARALLEL
   
   do i = 1, m_iter
-    B(i,m_iter+1) = -1
+    B(i,m_iter+1) = -1.d0
   enddo
   do j = 1, m_iter
-    B(m_iter+1,j) = -1
+    B(m_iter+1,j) = -1.d0
   enddo
   ! Debug
   !print*,'B'
@@ -493,7 +493,7 @@ subroutine update_t_ccsd_diis_v3(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err,all_t
   do i = 1, nO*nV
     tmp(i) = t1(i)
   enddo
-  !$OMP END DO NOWAIT
+  !$OMP END DO
   !$OMP DO
   do i = 1, nO*nO*nV*nV
     tmp(i+nO*nV) = t2(i)
@@ -515,7 +515,7 @@ subroutine update_t_ccsd_diis_v3(nO,nV,nb_iter,f_o,f_v,r1,r2,t1,t2,all_err,all_t
   do i = 1, nO*nV
     t1(i) = tmp(i)
   enddo
-  !$OMP END DO NOWAIT
+  !$OMP END DO
   !$OMP DO
   do i = 1, nO*nO*nV*nV
     t2(i) = tmp(i+nO*nV)