BEGIN_PROVIDER [ real, ao_axis_power_p, (-2:ao_power_max,3,nucl_num) ]
 implicit none

 BEGIN_DOC
! Evaluation of power of x, y, z at the current point for each
! nucleus. Negative power -> 0.
 END_DOC

 integer :: i,k,l
 do i=1,nucl_num
   do l=1,3
    ao_axis_power_p(-2,l,i) = 0.
    ao_axis_power_p(-1,l,i) = 0.
    ao_axis_power_p(0,l,i) = 0.
    ao_axis_power_p(0,l,i) = 1.
    do k=1,ao_power_max_nucl(i,l)
     ao_axis_power_p(k,l,i) = point_nucl_dist_vec(i,l)*ao_axis_power_p(k-1,l,i)
    enddo
   enddo
 enddo
END_PROVIDER

BEGIN_PROVIDER [ real, ao_axis_p, (ao_num) ]
 implicit none

 BEGIN_DOC
! Cartesian polynomial part of the atomic orbitals.
 END_DOC
 integer :: i

  do i=1,ao_num
    ao_axis_p(i)     &
        = ao_axis_power_p( ao_power(i,1) , 1 , ao_nucl(i) ) &
        * ao_axis_power_p( ao_power(i,2) , 2 , ao_nucl(i) ) &
        * ao_axis_power_p( ao_power(i,3) , 3 , ao_nucl(i) )
  enddo

END_PROVIDER

BEGIN_PROVIDER [ real, ao_axis_grad_p, (ao_num,3) ]
 implicit none

 BEGIN_DOC
! Gradients of the cartesian polynomial part of the atomic orbitals.
 END_DOC
 integer :: i, l
 real:: real_of_int(-1:10)
 data real_of_int /0.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10./

  do i=1,ao_num
     ao_axis_grad_p(i,1) = real_of_int(ao_power(i,1)) &
                              * ao_axis_power_p( ao_power(i,1)-1, 1 , ao_nucl(i) ) &
                              * ao_axis_power_p( ao_power(i,2)  , 2 , ao_nucl(i) ) &
                              * ao_axis_power_p( ao_power(i,3)  , 3 , ao_nucl(i) )
  enddo

  do i=1,ao_num
     ao_axis_grad_p(i,2) = real_of_int(ao_power(i,2)) &
                              * ao_axis_power_p( ao_power(i,1)  , 1 , ao_nucl(i) ) &
                              * ao_axis_power_p( ao_power(i,2)-1, 2 , ao_nucl(i) ) &
                              * ao_axis_power_p( ao_power(i,3)  , 3 , ao_nucl(i) )
  enddo

  do i=1,ao_num
     ao_axis_grad_p(i,3) = real_of_int(ao_power(i,3)) &
                              * ao_axis_power_p( ao_power(i,1)  , 1 , ao_nucl(i) ) &
                              * ao_axis_power_p( ao_power(i,2)  , 2 , ao_nucl(i) ) &
                              * ao_axis_power_p( ao_power(i,3)-1, 3 , ao_nucl(i) )
  enddo

END_PROVIDER

BEGIN_PROVIDER [ real, ao_axis_lapl_p, (ao_num) ]
 implicit none
 BEGIN_DOC
! Laplacian of the cartesian atomic orbitals 
 END_DOC
 integer :: i, j, l

 do i=1,ao_num
     real:: real_of_int(-2:10)
     data real_of_int /0.,0.,0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10./

      ao_axis_lapl_p(i)                                           &
             = real_of_int(ao_power(i,1))                       &
             * real_of_int(ao_power(i,1)-1)                     &
             * ao_axis_power_p( ao_power(i,1)-2, 1 , ao_nucl(i) ) &
             * ao_axis_power_p( ao_power(i,2)  , 2 , ao_nucl(i) ) &
             * ao_axis_power_p( ao_power(i,3)  , 3 , ao_nucl(i) ) &
             + real_of_int(ao_power(i,2))                       &
             * real_of_int(ao_power(i,2)-1)                     &
             * ao_axis_power_p( ao_power(i,1)  , 1 , ao_nucl(i) ) &
             * ao_axis_power_p( ao_power(i,2)-2, 2 , ao_nucl(i) ) &
             * ao_axis_power_p( ao_power(i,3)  , 3 , ao_nucl(i) ) &
             + real_of_int(ao_power(i,3))                       &
             * real_of_int(ao_power(i,3)-1)                     &
             * ao_axis_power_p( ao_power(i,1)  , 1 , ao_nucl(i) ) &
             * ao_axis_power_p( ao_power(i,2)  , 2 , ao_nucl(i) ) &
             * ao_axis_power_p( ao_power(i,3)-2, 3 , ao_nucl(i) )
   enddo

END_PROVIDER