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
  double precision, intent(in) :: r1(3), r2(3)
  double precision             :: mu_tmp, r12

  if(j2e_type .eq. "Mu") 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)) )
    mu_tmp = mu_erf * r12

    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
    stop

  endif ! j2e_type

  return
end

! ---

subroutine grad1_j12_mu(r1, r2, grad)

  BEGIN_DOC
  !
  !  gradient of j(mu(r1,r2),r12) form of jastrow. 
  !
  ! if mu(r1,r2) = cst ---> 
  !
  !  d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2)
  !
  ! if mu(r1,r2) /= cst ---> 
  !
  ! 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)
  !
  END_DOC

  include 'constants.include.F'

  implicit none
  double precision, intent(in)  :: r1(3), r2(3)
  double precision, intent(out) :: grad(3)
  double precision              :: dx, dy, dz, r12, tmp
  double precision              :: mu_val, mu_tmp, mu_der(3)

  grad = 0.d0

  if(j2e_type .eq. "Mu".or.j2e_type .eq. "Mugauss") then

    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

    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" .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                                                                                                                                                                            

     tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12

     grad(1) = 0.5d0 * tmp * dx
     grad(2) = 0.5d0 * tmp * dy
     grad(3) = 0.5d0 * tmp * dz
     grad(1) += 0.5d0 * grad_jast(1)
     grad(2) += 0.5d0 * grad_jast(2)
     grad(3) += 0.5d0 * grad_jast(3)


  else

    print *, ' Error in grad1_j12_mu: Unknown j2e_type = ', j2e_type
    stop

  endif ! j2e_type

  grad = -grad

  return
end

! ---

double precision function env_nucl(r)

  implicit none
  double precision, intent(in) :: r(3)
  integer                      :: i
  double precision             :: a, d, e, x, y, z

  if(env_type .eq. "Sum_Slat") then

    env_nucl = 1.d0
    do i = 1, nucl_num
      a = env_expo(i)
      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
      env_nucl = env_nucl - env_coef(i) * dexp(-a*dsqrt(d))
    enddo

  elseif(env_type .eq. "Prod_Gauss") then

    env_nucl = 1.d0
    do i = 1, nucl_num
      a = env_expo(i)
      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
      e = 1.d0 - dexp(-a*d)
      env_nucl = env_nucl * e
    enddo

  elseif(env_type .eq. "Sum_Gauss") then

    env_nucl = 1.d0
    do i = 1, nucl_num
      a = env_expo(i)
      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
      env_nucl = env_nucl - env_coef(i) * dexp(-a*d)
    enddo

  elseif(env_type .eq. "Sum_Quartic") then

    env_nucl = 1.d0
    do i = 1, nucl_num
      a = env_expo(i)
      x = r(1) - nucl_coord(i,1)
      y = r(2) - nucl_coord(i,2)
      z = r(3) - nucl_coord(i,3)
      d = x*x + y*y + z*z
      env_nucl = env_nucl - env_coef(i) * dexp(-a*d*d)
    enddo

  else

    print *, ' Error in env_nucl: Unknown env_type = ', env_type
    stop

  endif

  return
end

! ---

double precision function env_nucl_square(r)

  implicit none
  double precision, intent(in) :: r(3)
  integer                      :: i
  double precision             :: a, d, e, x, y, z

  if(env_type .eq. "Sum_Slat") then

    env_nucl_square = 1.d0
    do i = 1, nucl_num
      a = env_expo(i)
      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
      env_nucl_square = env_nucl_square - env_coef(i) * dexp(-a*dsqrt(d))
    enddo
    env_nucl_square = env_nucl_square * env_nucl_square

  elseif(env_type .eq. "Prod_Gauss") then

    env_nucl_square = 1.d0
    do i = 1, nucl_num
      a = env_expo(i)
      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
      e = 1.d0 - dexp(-a*d)
      env_nucl_square = env_nucl_square * e
    enddo
    env_nucl_square = env_nucl_square * env_nucl_square

  elseif(env_type .eq. "Sum_Gauss") then

    env_nucl_square = 1.d0
    do i = 1, nucl_num
      a = env_expo(i)
      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
      env_nucl_square = env_nucl_square - env_coef(i) * dexp(-a*d)
    enddo
    env_nucl_square = env_nucl_square * env_nucl_square

  elseif(env_type .eq. "Sum_Quartic") then

    env_nucl_square = 1.d0
    do i = 1, nucl_num
      a = env_expo(i)
      x = r(1) - nucl_coord(i,1)
      y = r(2) - nucl_coord(i,2)
      z = r(3) - nucl_coord(i,3)
      d = x*x + y*y + z*z
      env_nucl_square = env_nucl_square - env_coef(i) * dexp(-a*d*d)
    enddo
    env_nucl_square = env_nucl_square * env_nucl_square

  else

    print *, ' Error in env_nucl_square: Unknown env_type = ', env_type
    stop

  endif

  return
end

! ---

subroutine grad1_env_nucl(r, grad)

  implicit none
  double precision, intent(in)  :: r(3)
  double precision, intent(out) :: grad(3)
  integer                       :: ipoint, i, j, phase
  double precision              :: x, y, z, dx, dy, dz
  double precision              :: a, d, e
  double precision              :: fact_x, fact_y, fact_z
  double precision              :: ax_der, ay_der, az_der, a_expo

  if(env_type .eq. "Sum_Slat") then

    fact_x = 0.d0
    fact_y = 0.d0
    fact_z = 0.d0
    do i = 1, nucl_num
      a = env_expo(i)
      x = r(1) - nucl_coord(i,1)
      y = r(2) - nucl_coord(i,2)
      z = r(3) - nucl_coord(i,3)
      d = dsqrt(x*x + y*y + z*z)
      e = a * env_coef(i) * dexp(-a*d) / d

      fact_x += e * x
      fact_y += e * y
      fact_z += e * z
    enddo

    grad(1) = fact_x
    grad(2) = fact_y
    grad(3) = fact_z

  elseif(env_type .eq. "Prod_Gauss") then

    x = r(1)
    y = r(2)
    z = r(3)

    fact_x = 0.d0
    fact_y = 0.d0
    fact_z = 0.d0
    do i = 1, List_env1s_size

      phase  = 0
      a_expo = 0.d0
      ax_der = 0.d0
      ay_der = 0.d0
      az_der = 0.d0
      do j = 1, nucl_num
        a  = dble(List_env1s(j,i)) * env_expo(j)
        dx = x - nucl_coord(j,1)
        dy = y - nucl_coord(j,2)
        dz = z - nucl_coord(j,3)

        phase  += List_env1s(j,i)
        a_expo += a * (dx*dx + dy*dy + dz*dz)
        ax_der += a * dx
        ay_der += a * dy
        az_der += a * dz
      enddo
      e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo)

      fact_x += e * ax_der
      fact_y += e * ay_der
      fact_z += e * az_der
    enddo

    grad(1) = fact_x
    grad(2) = fact_y
    grad(3) = fact_z

  elseif(env_type .eq. "Sum_Gauss") then

    fact_x = 0.d0
    fact_y = 0.d0
    fact_z = 0.d0
    do i = 1, nucl_num
      a = env_expo(i)
      x = r(1) - nucl_coord(i,1)
      y = r(2) - nucl_coord(i,2)
      z = r(3) - nucl_coord(i,3)
      d = x*x + y*y + z*z
      e = a * env_coef(i) * dexp(-a*d)

      fact_x += e * x
      fact_y += e * y
      fact_z += e * z
    enddo

    grad(1) = 2.d0 * fact_x
    grad(2) = 2.d0 * fact_y
    grad(3) = 2.d0 * fact_z

  elseif(env_type .eq. "Sum_Quartic") then

    fact_x = 0.d0
    fact_y = 0.d0
    fact_z = 0.d0
    do i = 1, nucl_num
      a = env_expo(i)
      x = r(1) - nucl_coord(i,1)
      y = r(2) - nucl_coord(i,2)
      z = r(3) - nucl_coord(i,3)
      d = x*x + y*y + z*z
      e = a * env_coef(i) * d * dexp(-a*d*d)

      fact_x += e * x
      fact_y += e * y
      fact_z += e * z
    enddo

    grad(1) = 4.d0 * fact_x
    grad(2) = 4.d0 * fact_y
    grad(3) = 4.d0 * fact_z

  else

    print *, ' Error in grad1_env_nucl: Unknown env_type = ', env_type
    stop

  endif

  return
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)
  double precision              :: r(3)
  double precision              :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1)
  double precision              :: dm_tot, tmp1, tmp2, tmp3
  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. 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)
    !
    ! mu[rho(r)] = alpha sqrt(rho) + mu0 exp(-rho)
    !
    ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx
    ! d mu[rho(r)] / dx  = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx)
    !

    r(1) = 0.5d0 * (r1(1) + r2(1))
    r(2) = 0.5d0 * (r1(2) + r2(2))
    r(3) = 0.5d0 * (r1(3) + r2(3))

    call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b)

    dm_tot = dm_a(1) + dm_b(1)
    tmp1   = dsqrt(dm_tot)
    tmp2   = mu_erf * dexp(-dm_tot)

    mu_val = mu_r_ct * tmp1 + tmp2

    mu_der = 0.d0
    if(dm_tot .lt. 1d-7) return

    tmp3      = 0.25d0 * mu_r_ct / tmp1 - 0.5d0 * tmp2
    mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1))
    mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1))
    mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1))

  elseif(murho_type .eq. 2) then

    !
    ! r = 0.5 (r1 + r2)
    !
    ! mu[rho(r)] = alpha rho + mu0 exp(-rho)
    !
    ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx
    ! d mu[rho(r)] / dx  = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx)
    !

    r(1) = 0.5d0 * (r1(1) + r2(1))
    r(2) = 0.5d0 * (r1(2) + r2(2))
    r(3) = 0.5d0 * (r1(3) + r2(3))

    call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b)

    dm_tot = dm_a(1) + dm_b(1)
    tmp2   = mu_erf * dexp(-dm_tot)

    mu_val = mu_r_ct * dm_tot + tmp2

    tmp3      = 0.5d0 * (mu_r_ct - tmp2)
    mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1))
    mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1))
    mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1))

  elseif(murho_type .eq. 3) then

    ! mu(r1,r2) = {rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]} / RHO
    !
    ! RHO = rho(r1) + rho(r2)
    !
    ! f[rho] = alpha rho^beta + mu0 exp(-rho)
    !
    ! d/dx1 mu(r1,r2) = 1/RHO^2 * {RHO * d/dx1 (rho(r1) f[rho(r1)]) 
    !                              - d/dx1 rho(r1) * [rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]] }
    !
    ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) - mu0 exp(-rho(r1))] (d rho(r1) / dx1)
    !
    ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1)
     
    !!!!!!!!! rho1,rho2,rho1+rho2
    call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1)
    rho_tot = rho1 + rho2
    if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
    inv_rho_tot = 1.d0/rho_tot
    ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf * exp(-rho)
    call get_all_f_rho(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2)
    d_dx1_f_rho1(1:3)   = d_drho_f_rho1 * grad_rho1(1:3)
    d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3)
    nume   = rho1 * f_rho1 + rho2 * f_rho2
    mu_val = nume * inv_rho_tot
    mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume)

  elseif(murho_type .eq. 4) then

    ! mu(r1,r2) = {rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]} / RHO
    !
    ! RHO = rho(r1) + rho(r2)
    !
    ! f[rho] = alpha rho^beta + mu0 
    !
    ! d/dx1 mu(r1,r2) = 1/RHO^2 * {RHO * d/dx1 (rho(r1) f[rho(r1)]) 
    !                              - d/dx1 rho(r1) * [rho(r1) f[rho(r1)] + rho(r2) f[rho(r2)]] }
    !
    ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) ] (d rho(r1) / dx1)
    !
    ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1)
     
    !!!!!!!!! rho1,rho2,rho1+rho2
    call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1)
    rho_tot = rho1 + rho2
!    if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
    if(rho_tot.lt.1.d-10)then
     mu_val = mu_erf 
     mu_der = 0.d0
     return
    endif
    
    if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
    inv_rho_tot = 1.d0/rho_tot
    ! f(rho) = mu_r_ct * rho**beta_rho_power + mu_erf 
    call get_all_f_rho_simple(rho1,rho2,mu_r_ct,mu_erf,beta_rho_power,f_rho1,d_drho_f_rho1,f_rho2)
    d_dx1_f_rho1(1:3)   = d_drho_f_rho1 * grad_rho1(1:3)
    d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3)
    nume   = rho1 * f_rho1 + rho2 * f_rho2
    mu_val = nume * inv_rho_tot
    mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume)

  elseif(murho_type .eq. 5) then

    ! mu(r1,r2) = 1/2 * (f[rho(r1)] + f[rho(r2)]} 
    !
    ! f[rho] = alpha rho^beta + mu0 
    !
    ! d/dx1 mu(r1,r2) = 1/2 * d/dx1 (rho(r1) f[rho(r1)])
    !                   
    ! d/dx1 f[rho(r1)] = [0.5 alpha / sqrt(rho(r1)) ] (d rho(r1) / dx1)
    !
    ! d/dx1 (rho(r1) f[rho(r1)] = rho(r1) * d/dx1 f[rho(r1)] + f[rho(r1)] * d/dx1 rho(r1)
    !!!!!!!!! rho1,rho2,rho1+rho2
    call get_all_rho_grad_rho(r1,r2,rho1,rho2,grad_rho1)
    rho_tot = rho1 + rho2
!    if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
    if(rho_tot.lt.1.d-10)then
     mu_val = mu_erf 
     mu_der = 0.d0
     return
    endif
    
    if(rho_tot.lt.1.d-10)rho_tot = 1.d-10
    inv_rho_tot = 1.d0/rho_tot
    ! f(rho) = (mu_r_ct* rho)**beta_rho_power * erf(zeta_erf_mu_of_r * rho) + mu_eff * (1 - erf(zeta_erf_mu_of_r*rho))
    call get_all_f_rho_erf(rho1,rho2,mu_r_ct,beta_rho_power,mu_erf,zeta_erf_mu_of_r,f_rho1,d_drho_f_rho1,f_rho2)
    d_dx1_f_rho1(1:3)   = d_drho_f_rho1 * grad_rho1(1:3)
    d_dx_rho_f_rho(1:3) = rho1 * d_dx1_f_rho1(1:3) + f_rho1 * grad_rho1(1:3)
    nume   = rho1 * f_rho1 + rho2 * f_rho2
    mu_val = nume * inv_rho_tot
    mu_der(1:3) = inv_rho_tot*inv_rho_tot * (rho_tot * d_dx_rho_f_rho(1:3) - grad_rho1(1:3) * nume)
     
  else

    print *, ' Error in mu_r_val_and_grad: Unknown env_type = ', env_type
    stop

  endif

  return
end

! ---

subroutine grad1_env_nucl_square_num(r1, grad)

  implicit none
  double precision, intent(in)  :: r1(3)
  double precision, intent(out) :: grad(3)
  double precision              :: r(3), eps, tmp_eps, vp, vm
  double precision, external    :: env_nucl_square

  eps     = 1d-5
  tmp_eps = 0.5d0 / eps

  r(1:3) = r1(1:3)

  r(1) = r(1) + eps
  vp   = env_nucl_square(r)
  r(1) = r(1) - 2.d0 * eps
  vm   = env_nucl_square(r)
  r(1) = r(1) + eps
  grad(1) = tmp_eps * (vp - vm)

  r(2) = r(2) + eps
  vp   = env_nucl_square(r)
  r(2) = r(2) - 2.d0 * eps
  vm   = env_nucl_square(r)
  r(2) = r(2) + eps
  grad(2) = tmp_eps * (vp - vm)

  r(3) = r(3) + eps
  vp   = env_nucl_square(r)
  r(3) = r(3) - 2.d0 * eps
  vm   = env_nucl_square(r)
  r(3) = r(3) + eps
  grad(3) = tmp_eps * (vp - vm)
  
  return
end

! ---

subroutine grad1_j12_mu_square_num(r1, r2, grad)

  include 'constants.include.F'

  implicit none
  double precision, intent(in)  :: r1(3), r2(3)
  double precision, intent(out) :: grad(3)
  double precision              :: r(3)
  double precision              :: eps, tmp_eps, vp, vm
  double precision, external    :: j12_mu_square

  eps     = 1d-5
  tmp_eps = 0.5d0 / eps

  r(1:3) = r1(1:3)

  r(1)    = r(1) + eps
  vp      = j12_mu_square(r, r2)
  r(1)    = r(1) - 2.d0 * eps
  vm      = j12_mu_square(r, r2)
  r(1)    = r(1) + eps
  grad(1) = tmp_eps * (vp - vm)

  r(2)    = r(2) + eps
  vp      = j12_mu_square(r, r2)
  r(2)    = r(2) - 2.d0 * eps
  vm      = j12_mu_square(r, r2)
  r(2)    = r(2) + eps
  grad(2) = tmp_eps * (vp - vm)

  r(3)    = r(3) + eps
  vp      = j12_mu_square(r, r2)
  r(3)    = r(3) - 2.d0 * eps
  vm      = j12_mu_square(r, r2)
  r(3)    = r(3) + eps
  grad(3) = tmp_eps * (vp - vm)

  return
end

! ---

double precision function j12_mu_square(r1, r2)

  implicit none
  double precision, intent(in) :: r1(3), r2(3)
  double precision, external   :: j12_mu

  j12_mu_square = j12_mu(r1, r2) * j12_mu(r1, r2)

  return
end

! ---

subroutine f_mu_and_deriv_mu(rho, alpha, mu0, beta, f_mu, d_drho_f_mu)

  BEGIN_DOC
  ! function giving mu as a function of rho
  !
  ! f_mu = alpha * rho**beta + mu0 * exp(-rho)
  !
  ! and its derivative with respect to rho d_drho_f_mu
  END_DOC

  implicit none
  double precision, intent(in)  :: rho, alpha, mu0, beta
  double precision, intent(out) :: f_mu, d_drho_f_mu

  f_mu = alpha * (rho)**beta + mu0 * dexp(-rho)
  d_drho_f_mu = alpha * beta * rho**(beta-1.d0) - mu0 * dexp(-rho)

end

! ---

subroutine get_all_rho_grad_rho(r1, r2, rho1, rho2, grad_rho1)

  BEGIN_DOC
  ! returns the density in r1,r2 and grad_rho at r1
  END_DOC

  implicit none
  double precision, intent(in)  :: r1(3), r2(3)
  double precision, intent(out) :: grad_rho1(3), rho1, rho2
  double precision              :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1)

  call density_and_grad_alpha_beta(r1, dm_a, dm_b, grad_dm_a, grad_dm_b)
  rho1 = dm_a(1) + dm_b(1)
  grad_rho1(1:3) = grad_dm_a(1:3,1) + grad_dm_b(1:3,1)
  call density_and_grad_alpha_beta(r2, dm_a, dm_b, grad_dm_a, grad_dm_b)
  rho2 = dm_a(1) + dm_b(1)

end

! ---

subroutine get_all_f_rho(rho1, rho2, alpha, mu0, beta, f_rho1, d_drho_f_rho1, f_rho2)

  BEGIN_DOC
  ! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1))
  END_DOC

  implicit none
  double precision, intent(in)  :: rho1, rho2, alpha, mu0, beta
  double precision, intent(out) :: f_rho1, d_drho_f_rho1, f_rho2
  double precision              :: tmp

  call f_mu_and_deriv_mu(rho1, alpha, mu0, beta, f_rho1, d_drho_f_rho1)
  call f_mu_and_deriv_mu(rho2, alpha, mu0, beta, f_rho2, tmp)

end

! ---

subroutine get_all_f_rho_simple(rho1,rho2,alpha,mu0,beta,f_rho1,d_drho_f_rho1,f_rho2)

  BEGIN_DOC
  ! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1))
  END_DOC

  implicit none
  double precision, intent(in)  :: rho1, rho2, alpha, mu0, beta
  double precision, intent(out) :: f_rho1, d_drho_f_rho1, f_rho2
  double precision              :: tmp

  if(rho1.lt.1.d-10) then
    f_rho1 = 0.d0
    d_drho_f_rho1 = 0.d0
  else
    call f_mu_and_deriv_mu_simple(rho1, alpha, mu0, beta, f_rho1, d_drho_f_rho1)
  endif

  if(rho2.lt.1.d-10)then
    f_rho2 = 0.d0
  else
    call f_mu_and_deriv_mu_simple(rho2, alpha, mu0, beta, f_rho2, tmp)
  endif

end

! ---

subroutine f_mu_and_deriv_mu_simple(rho, alpha, mu0, beta, f_mu, d_drho_f_mu)

  BEGIN_DOC
  ! function giving mu as a function of rho
  !
  ! f_mu = alpha * rho**beta + mu0 
  !
  ! and its derivative with respect to rho d_drho_f_mu
  END_DOC

  implicit none
  double precision, intent(in)  :: rho, alpha, mu0, beta
  double precision, intent(out) :: f_mu, d_drho_f_mu

  f_mu = alpha**beta * (rho)**beta + mu0 
  d_drho_f_mu = alpha**beta * beta * rho**(beta-1.d0) 

end

! ---

subroutine f_mu_and_deriv_mu_erf(rho,alpha,zeta,mu0,beta,f_mu,d_drho_f_mu)

  include 'constants.include.F'

  BEGIN_DOC
  ! function giving mu as a function of rho
  !
  ! f_mu = (alpha * rho)**zeta * erf(beta * rho) + mu0 * (1 - erf(beta*rho))
  !
  ! and its derivative with respect to rho d_drho_f_mu
  !
  ! d_drho_f_mu = 2 beta/sqrt(pi) * exp(-(beta*rho)**2) * ( (alpha*rho)**zeta - mu0) 
  !               + alpha * zeta * (alpha *rho)**(zeta-1)  * erf(beta*rho)
  END_DOC

  implicit none
  double precision, intent(in)  :: rho, alpha, mu0, beta, zeta
  double precision, intent(out) :: f_mu, d_drho_f_mu

  f_mu = (alpha * rho)**zeta * derf(beta * rho) + mu0 * (1.d0 - derf(beta*rho))
  d_drho_f_mu = 2.d0  * beta * inv_sq_pi * dexp(-(beta*rho)**2) * ( (alpha*rho)**zeta - mu0) & 
              + alpha * zeta * (alpha *rho)**(zeta-1)  * derf(beta*rho)

end

! ---

subroutine get_all_f_rho_erf(rho1, rho2, alpha, zeta, mu0, beta, f_rho1, d_drho_f_rho1, f_rho2)

  BEGIN_DOC
  ! returns the values f(mu(r1)), f(mu(r2)) and d/drho(1) f(mu(r1))
  ! with f_mu = (alpha * rho)**zeta * erf(beta * rho) + mu0 * (1 - erf(beta*rho))
  END_DOC

  implicit none
  double precision, intent(in)  :: rho1, rho2, alpha, mu0, beta, zeta
  double precision, intent(out) :: f_rho1, d_drho_f_rho1, f_rho2
  double precision              :: tmp

  if(rho1 .lt. 1.d-10) then
    f_rho1 = mu_erf
    d_drho_f_rho1 = 0.d0
  else
    call f_mu_and_deriv_mu_erf(rho1, alpha, zeta, mu0, beta, f_rho1, d_drho_f_rho1)
  endif

  if(rho2 .lt. 1.d-10)then
    f_rho2 = mu_erf
  else
    call f_mu_and_deriv_mu_erf(rho2, alpha, zeta, mu0, beta, f_rho2, tmp)
  endif

end

! ---