program print_density_on_grid
  implicit none
  BEGIN_DOC
! Print sum of alpha and beta electron density on the x axis
  END_DOC

  double precision :: r(3) ! Array for spatial coordinates 'r(x, y, z)' (a.u.)
  double precision :: dm_a(N_states), dm_b(N_states) ! densities at r
  double precision :: n_at_r
  double precision :: xlength, dx
  integer :: nx, ipoint

  r(1) = nucl_coord(1,1)
  r(2) = nucl_coord(1,2)
  r(3) = nucl_coord(1,3)

  nx = 500
  xlength = 7.d0
  dx =  xlength/dble(nx)
  r(1)+=-xlength*0.5d0

  do ipoint=0, nx
   call dm_dft_alpha_beta_at_r(r, dm_a, dm_b)
   n_at_r= dm_a(1) + dm_b(1)
   print *, r(1), n_at_r
   r(1)+=dx
  enddo

end