BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [ integer, ao_ortho_canonical_num ]
  implicit none
  BEGIN_DOC
! matrix of the coefficients of the mos generated by the 
! orthonormalization by the S^{-1/2} canonical transformation of the aos
! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
  END_DOC
  integer :: i
  ao_ortho_canonical_coef(:,:) = 0.d0
  if (ao_cartesian) then
    do i=1,ao_num
      ao_ortho_canonical_coef(i,i) = 1.d0
    enddo
  else
    integer, external              :: ao_power_index
    integer                        :: ibegin,j,k
    ! Assume order provided by ao_power_index
    i = 1
    do while (i <= ao_num)
      select case ( ao_l(i) )
        case (0)
          ao_ortho_canonical_coef(i,i) = 1.d0
          i += 1
        BEGIN_TEMPLATE
        case ($SHELL)
          if (ao_power(i,1) == $SHELL) then
            do k=1,size(cart_to_sphe_$SHELL,2)
              do j=1,size(cart_to_sphe_$SHELL,1)
                ao_ortho_canonical_coef(i+j-1,i+k-1) = cart_to_sphe_$SHELL(j,k)
              enddo
            enddo
            i += size(cart_to_sphe_$SHELL,1)
          endif
        SUBST [ SHELL ]
          1;;
          2;;
          3;;
          4;;
          5;;
          6;;
          7;;
          8;;
          9;;
        END_TEMPLATE
        case default
          stop 'Error in sphe_to_cart'
      end select
    enddo
  endif
  call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,ao_ortho_canonical_num)
END_PROVIDER

BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)]
  implicit none
  BEGIN_DOC
! overlap matrix of the ao_ortho_canonical.
! Expected to be the Identity
  END_DOC
  integer                        :: i,j,k,l
  double precision               :: c
  do j=1, ao_ortho_canonical_num
    do i=1, ao_ortho_canonical_num
      ao_ortho_canonical_overlap(i,j) = 0.d0
    enddo
  enddo
  do j=1, ao_ortho_canonical_num
    do k=1, ao_num
      c = 0.d0
      do l=1, ao_num
        c +=  ao_ortho_canonical_coef(l,j) * ao_overlap(l,k)
      enddo
      do i=1, ao_ortho_canonical_num
        ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(k,i) * c
      enddo
    enddo
  enddo
END_PROVIDER