BEGIN_PROVIDER [ integer, nucl_num ]
 implicit none

 BEGIN_DOC
! Number of nuclei
 END_DOC

 nucl_num = -1
 call get_nuclei_nucl_num(nucl_num)
 if (nucl_num <= 0) then
  call abrt(irp_here,'Number of nuclei should be > 0')
 endif

END_PROVIDER

BEGIN_PROVIDER  [ real, nucl_charge, (nucl_num) ]
 implicit none

 BEGIN_DOC
! Nuclear charge
 END_DOC

 nucl_charge = -1.d0
 call get_nuclei_nucl_charge(nucl_charge)

 integer :: i
 do i=1,nucl_num
   if (nucl_charge(i) <= 0.) then
    call abrt(irp_here,'Nuclear charges should be > 0')
   endif
 enddo
END_PROVIDER

BEGIN_PROVIDER [ real, nucl_coord,  (nucl_num,3) ]
 implicit none

 BEGIN_DOC
! Nuclear coordinates
 END_DOC

 nucl_coord = 0.
 call get_nuclei_nucl_coord(nucl_coord)

END_PROVIDER


 BEGIN_PROVIDER [ real, nucl_dist_2, (nucl_num,nucl_num) ]
&BEGIN_PROVIDER [ real, nucl_dist_vec, (nucl_num,nucl_num,3) ]
 implicit none
 BEGIN_DOC
! nucl_dist_2   : Nucleus-nucleus distances squared
! nucl_dist_vec : Nucleus-nucleus distances vectors
 END_DOC

 integer :: ie1, ie2, l

 do ie2 = 1,nucl_num
  do ie1 = 1,nucl_num
   nucl_dist_2(ie1,ie2) = 0.d0
  enddo
 enddo

 do l=1,3
  do ie2 = 1,nucl_num
   do ie1 = 1,nucl_num
    nucl_dist_vec(ie1,ie2,l) = nucl_coord(ie1,l) - nucl_coord(ie2,l)
   enddo
   do ie1 = 1,nucl_num
    nucl_dist_2(ie1,ie2) = nucl_dist_2(ie1,ie2) &
      + nucl_dist_vec(ie1,ie2,l)*nucl_dist_vec(ie1,ie2,l) 
   enddo
  enddo
 enddo
END_PROVIDER

BEGIN_PROVIDER [ real, nucl_dist, (nucl_num,nucl_num) ]
 implicit none
 BEGIN_DOC
! Nucleus-nucleus distances
 END_DOC
 integer :: ie1, ie2
 do ie2 = 1,nucl_num
  do ie1 = 1,nucl_num
   nucl_dist(ie1,ie2) = sqrt(nucl_dist_2(ie1,ie2))
  enddo
 enddo
END_PROVIDER