anchor.diff:    integer :: s,i,j
     double precision, external :: get_dij_index
     print *, "computing amplitudes..."
  +  !$OMP PARALLEL DEFAULT(shared) PRIVATE(s,i,j)
     do s=1, N_states
  +    !$OMP DO
       do i=1, N_det_non_ref
         do j=1, N_det_ref
           dij(j, i, s) = get_dij_index(j, i, s, N_int)
         end do
       end do
  -  end do
  +    !$OMP END DO 
  +  end do
  +  !$OMP END PARALLEL
     print *, "done computing amplitudes"
   END_PROVIDER
   
  
positive.code:        if (rho_mrcc(i,s) == 0.d0) then
           rho_mrcc(i,s) = 1.d-32
         endif
         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
  
       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)
       enddo
       f = dsqrt(f*norm)
  
       print *,  'norm of |T Psi_0> = ', norm*f*f
  
       do i=1,N_det_non_ref
         rho_mrcc(i,s) = rho_mrcc(i,s) * f
       enddo
  
     end do
     
     print *, "done"
  END_PROVIDER
  
  
  
negative.code:        
         if(res < 1d-8) exit
       end do
       ! rho_mrcc now contains A.X
  
       norm = 0.d0
       do i=1,N_det_non_ref
         norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)
       enddo
       ! Norm now contains the norm of A.X
  
       do i=1,N_det_ref
         norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
       enddo
       ! Norm now contains the norm of Psi + A.X
       
       print *, k, "res : ", res, "norm : ", sqrt(norm)
       
  