BEGIN_PROVIDER [ double precision, Energy_mu ]
  implicit none
  BEGIN_DOC
  ! E mu
  ! Eq.(26-30)
  END_DOC

  integer                        :: i
  energy_mu = E_nucl
  !DIR$ VECTOR ALIGNED
  !DIR$ LOOP COUNT(200)
  do i=1,elec_num
   energy_mu += E_kin_elec(i) 
  enddo
  energy_mu += Eff_pot_mu + eff_pot_deriv_mu + E_nucl_elec - three_body_mu 

  energy_mu_min = min(energy_mu_min,energy_mu)
  energy_mu_max = max(energy_mu_max,energy_mu)
  SOFT_TOUCH energy_mu_min energy_mu_max
END_PROVIDER

BEGIN_PROVIDER [double precision, E_nucl_elec]
 implicit none
 !TODO

 integer :: i,j
 E_nucl_elec = 0.d0
 do i = 1, elec_num
!  E_nucl_elec += E_pot_elec_one(i) + E_pot_elec_two(i)
  E_nucl_elec += E_pot_elec_one(i) 
 enddo
  E_nucl_elec_min = min(E_nucl_elec_min,E_nucl_elec)
  E_nucl_elec_max = max(E_nucl_elec_max,E_nucl_elec)
END_PROVIDER 


 BEGIN_PROVIDER [double precision, Eff_pot_mu_elec, (elec_num)]
&BEGIN_PROVIDER [double precision, Eff_pot_mu_elec_simple, (elec_num)]
 implicit none
 include '../constants.F'
 BEGIN_DOC
 ! Eq.(32)
 END_DOC
 integer :: i,j
 double precision :: rij, mu
 mu = mu_erf
 Eff_pot_mu_elec = 0.d0
 do i=1,elec_num
   !DIR$ VECTOR ALIGNED
   !DIR$ LOOP COUNT(50)
   do j=1,elec_num
     rij = elec_dist(j,i)
     if(i==j)cycle
     Eff_pot_mu_elec(i) = Eff_pot_mu_elec(i) + 0.5d0 * derf(mu * rij) * elec_dist_inv(j,i)
     Eff_pot_mu_elec(i) = Eff_pot_mu_elec(i) + 0.5d0 * mu/dsqpi * dexp(-mu*mu*rij*rij)
     Eff_pot_mu_elec_simple(i) = Eff_pot_mu_elec(i)
     Eff_pot_mu_elec(i) = Eff_pot_mu_elec(i) + 0.5d0 * (-  0.25d0 * (1.d0 - derf(mu*rij))**2.d0 )
                                             
   enddo
 enddo

END_PROVIDER 

 BEGIN_PROVIDER [double precision, Eff_pot_mu ]
 implicit none
 include '../constants.F'
 BEGIN_DOC
 ! Eq.(32) 
 END_DOC
 integer :: i
 Eff_pot_mu = 0.d0
 do i=1,elec_num
  Eff_pot_mu += eff_pot_mu_elec(i) 
 enddo
  Eff_pot_mu_min = min(Eff_pot_mu_min,Eff_pot_mu)
  Eff_pot_mu_max = max(Eff_pot_mu_max,Eff_pot_mu)
  SOFT_TOUCH Eff_pot_mu_min Eff_pot_mu_max

END_PROVIDER 

BEGIN_PROVIDER [double precision, Eff_pot_mu_simple ]
 implicit none
 BEGIN_DOC
 ! Eq.(32)
 END_DOC
 include '../constants.F'
 integer :: i
 Eff_pot_mu_simple = 0.d0
 do i=1,elec_num
  Eff_pot_mu_simple += Eff_pot_mu_elec_simple(i) 
 enddo
  Eff_pot_mu_simple_min = min(Eff_pot_mu_simple_min,Eff_pot_mu_simple)
  Eff_pot_mu_simple_max = max(Eff_pot_mu_simple_max,Eff_pot_mu_simple)
  SOFT_TOUCH Eff_pot_mu_simple_min Eff_pot_mu_simple_max

END_PROVIDER 

BEGIN_PROVIDER [double precision, eff_pot_deriv_mu_elec, (elec_num) ]
 implicit none
 BEGIN_DOC
 ! Eq.(33)
 END_DOC
 integer :: i,j
 double precision :: rij, mu
 mu = mu_erf
 eff_pot_deriv_mu_elec = 0.d0
 do i = 1, elec_num
  do j = 1, elec_num
   if(i==j)cycle
   rij = elec_dist(i,j)
   eff_pot_deriv_mu_elec(i) += 0.5d0 * ( derf(mu * rij) - 1.d0 ) * elec_dist_inv(j,i)  &
                                     * (  - elec_dist_vec_x(j,i) * psidet_grad_lapl(1,i)        &
                                          - elec_dist_vec_y(j,i) * psidet_grad_lapl(2,i)        &
                                          - elec_dist_vec_z(j,i) * psidet_grad_lapl(3,i) ) * psidet_inv
  enddo
 enddo

END_PROVIDER 


BEGIN_PROVIDER [double precision, three_body_mu ]
 implicit none
 BEGIN_DOC
 ! Eq.(30)
 END_DOC
 integer :: i,j,k
 three_body_mu = 0.d0
 do i = 1, elec_num
  do j = i+1, elec_num
   do k = j+1, elec_num
    three_body_mu += grad_j_mu_x(i,j) * grad_j_mu_x(i,k) 
    three_body_mu += grad_j_mu_y(i,j) * grad_j_mu_y(i,k) 
    three_body_mu += grad_j_mu_z(i,j) * grad_j_mu_z(i,k) 

    three_body_mu += grad_j_mu_x(j,i) * grad_j_mu_x(j,k) 
    three_body_mu += grad_j_mu_y(j,i) * grad_j_mu_y(j,k) 
    three_body_mu += grad_j_mu_z(j,i) * grad_j_mu_z(j,k) 

    three_body_mu += grad_j_mu_x(k,i) * grad_j_mu_x(k,j) 
    three_body_mu += grad_j_mu_y(k,i) * grad_j_mu_y(k,j) 
    three_body_mu += grad_j_mu_z(k,i) * grad_j_mu_z(k,j) 
   enddo
  enddo
 enddo
 three_body_mu_min = min(three_body_mu_min,three_body_mu)
 three_body_mu_max = max(three_body_mu_max,three_body_mu)
 SOFT_TOUCH  three_body_mu_min three_body_mu_max
END_PROVIDER 

BEGIN_PROVIDER [double precision, eff_pot_deriv_mu]
 implicit none
 BEGIN_DOC
 ! Eq.(33)
 END_DOC
 integer :: i
 eff_pot_deriv_mu = 0.d0
 do i = 1, elec_num
  eff_pot_deriv_mu += eff_pot_deriv_mu_elec(i)
 enddo
  eff_pot_deriv_mu_min = min(eff_pot_deriv_mu_min,eff_pot_deriv_mu)
  eff_pot_deriv_mu_max = max(eff_pot_deriv_mu_max,eff_pot_deriv_mu)
  SOFT_TOUCH eff_pot_deriv_mu_min eff_pot_deriv_mu_max

END_PROVIDER