BEGIN_PROVIDER [ double precision, alpha_coef, (0:1) ]
 implicit none
 BEGIN_DOC
 ! SavCar-JCP-23
 END_DOC

 double precision :: num, den

 num = 0.319820d0 + mu_erf * (1.063846d0 + mu_erf)
 den = 0.487806d0 + mu_erf * (1.375439d0 + mu_erf)
 alpha_coef(0) = num/den

 num = 0.113074d0 + mu_erf * (0.638308d0 + mu_erf)
 den = 0.122652d0 + mu_erf * (0.674813d0 + mu_erf)
 alpha_coef(1) = num/den
END_PROVIDER


BEGIN_PROVIDER [ double precision, alpha_coef_r, (0:1,N_states) ]
 implicit none
 BEGIN_DOC
 ! SavCar-JCP-23, corrected for small mu
 END_DOC

 integer :: istate
 double precision :: num, den
 double precision :: a0, b0, delta_E, E

 do istate = 1,N_states
   delta_E = (energy_mu(istate) + correction_alpha_1(istate)) - energy_mu0(istate)
   a0 = delta_E * correction_mu0(istate) + 2.d0/dsqrt(dacos(-1.d0))*(delta_E - correction_mu0(istate)) * mu_erf
   b0 = correction_mu0(istate)*correction_mu0(istate)

   num = 0.319820d0 + mu_erf * (1.063846d0 + mu_erf)
   den = 0.487806d0 + mu_erf * (1.375439d0 + mu_erf)
   alpha_coef_r(0,istate) = (a0 + mu_erf*num)/(b0 + mu_erf*den)


  ! delta_E = (energy_mu + correction_alpha_1) - energy_mu0
  ! a0 = delta_E * correction_mu0 + 2.d0/dsqrt(dacos(-1.d0))*(delta_E - correction_mu0) * mu_erf
  ! b0 = correction_mu0*correction_mu0

   num = 0.113074d0 + mu_erf * (0.638308d0 + mu_erf)
   den = 0.122652d0 + mu_erf * (0.674813d0 + mu_erf)
   alpha_coef_r(1,istate) = (a0 + mu_erf*num)/(b0 + mu_erf*den)
 enddo

END_PROVIDER

BEGIN_PROVIDER [ double precision, energy_mu0, (N_states) ]
 implicit none
 BEGIN_DOC
 ! E(mu=0)
 END_DOC

 integer :: l, istate
 double precision, external :: ddot

 do istate = 1,N_states
   energy_mu0(istate) = 0.d0
   do l=1,mo_num
     energy_mu0(istate) += ddot(mo_num, one_e_dm_mo_alpha(1,l,istate), 1, mo_one_e_integrals(1,l), 1)
     energy_mu0(istate) += ddot(mo_num, one_e_dm_mo_beta(1,l,istate), 1, mo_one_e_integrals(1,l), 1)
   enddo
   energy_mu0(istate) = energy_mu0(istate) + nuclear_repulsion
 enddo

END_PROVIDER

BEGIN_PROVIDER [ double precision, correction_alpha_0_r, (N_states) ]
 implicit none
 BEGIN_DOC
 ! alpha_0r(mu) * <Psi(mu)|W_bar_mu|Psi(mu)>
 END_DOC

 integer :: k,l, istate
 double precision :: c0
 double precision, external :: ddot

 do istate = 1,N_states
   c0 = 0.5d0 * alpha_coef_r(0,istate)
   correction_alpha_0_r(istate) = 0.d0
   do l=1,mo_num
     do k=1,mo_num
       correction_alpha_0_r += c0 * ddot (mo_num*mo_num, full_occ_2_rdm_spin_trace_mo(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1)
     enddo
   enddo
 enddo
END_PROVIDER

BEGIN_PROVIDER [ double precision, correction_alpha_1_r, (N_states) ]
 implicit none
 BEGIN_DOC
 ! alpha_0(mu) * <Psi(mu)|W_bar_mu|Psi(mu)>_s + alpha_1(mu) * <Psi(mu)|W_bar_mu|Psi(mu)>_t
 END_DOC

 double precision :: c1, c0

 integer :: k,l, istate
 double precision, external :: ddot

 do istate = 1,N_states
   c0 = 0.5d0 * alpha_coef_r(0,istate)
   c1 = 0.5d0 * alpha_coef_r(1,istate)

   correction_alpha_1_r(istate) = 0.d0
   do l=1,mo_num
     do k=1,mo_num
       correction_alpha_1_r(istate) += c0 * ddot (mo_num*mo_num, two_e_dm_mo_singlet(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1)
       correction_alpha_1_r(istate) += c1 * ddot (mo_num*mo_num, two_e_dm_mo_triplet(1,1,k,l,istate), 1, W_bar_mu(1,1,k,l), 1)
     enddo
   enddo
 enddo
END_PROVIDER