anchor.diff:         if (rho_mrcc(i,s) == 0.d0) then
            rho_mrcc(i,s) = 1.d-32
          endif
  +
  +       ! f is such that f.\tilde{c_i} = c_i
          f = psi_non_ref_coef(i,s) / rho_mrcc(i,s)
  +
          ! Avoid numerical instabilities
          f = min(f,2.d0)
          f = max(f,-2.d0)
  +
          norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s)
          rho_mrcc(i,s) = f
        enddo
  +     ! norm now contains the norm of |T.Psi_0>
  +     ! rho_mrcc now contains the f factors
   
        f = 1.d0/norm
  -     norm = 0.d0
  -     do i=1,N_det_non_ref
  -       norm = norm + psi_non_ref_coef(i,s)*psi_non_ref_coef(i,s)
  +     ! f now contains 1/ <T.Psi_0|T.Psi_0>
  +
  +     norm = 1.d0
  +     do i=1,N_det_ref
  +       norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s)
        enddo
  +     ! norm now contains <Psi_SD|Psi_SD>
        f = dsqrt(f*norm)
  -
  -     print *,  'norm of |T Psi_0> = ', norm*f*f
  +     ! f normalises T.Psi_0 such that (1+T)|Psi> is normalized
  +
  +     norm = norm*f
  +     print *,  'norm of |T Psi_0> = ', dsqrt(norm)
  +
  +     do i=1,N_det_ref
  +       norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
  +     enddo
   
        do i=1,N_det_non_ref
          rho_mrcc(i,s) = rho_mrcc(i,s) * f
        enddo
  +     ! rho_mrcc now contains the product of the scaling factors and the
  +     ! normalization constant
   
      end do
      
  -   print *, "done"
   END_PROVIDER
   
   
  
positive.code:        
         if(res < 1d-8) exit
       end do
       
       norm = 0.d0
       do i=1,N_det_non_ref
         norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)
       enddo
       do i=1,N_det_ref
         norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
       enddo
       
       print *, k, "res : ", res, "norm : ", sqrt(norm)
       
  
negative.code:   double precision :: HIi, phase
  
    if(lambda_type == 0) then
      call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int)
      get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase
      get_dij_index = get_dij_index * rho_mrcc(i,s) 
    else
      call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi)
      get_dij_index = HIi * lambda_mrcc(s, i)
  