Expose P, Q in self-consistent DFT orbitals

[certik]

We should also make the following change:

diff --git a/src/dft.f90 b/src/dft.f90
index 94db0ee..58a51c6 100644
--- a/src/dft.f90
+++ b/src/dft.f90
@@ -36,7 +36,7 @@ subroutine V2rho(d)
 ! Calculates rho from V by solving Kohn-Sham equations
 type(dft_data_t), intent(inout) :: d

-real(dp), dimension(size(d%R)) :: P, Q, Y
+real(dp), dimension(size(d%R)) :: P, Q, Y2
 integer :: converged, i, n, l, relat
 real(dp) :: Ein, Emin_init, Emax_init

@@ -71,12 +71,12 @@ do i = 1, size(d%no)
         call stop_error("V2rho: Radial eigen problem didn't converge")
     end if
     if (relat == 0) then
-        Y = P / d%R
+        Y2 = (P / d%R)**2
     else
-        Y = sqrt(P**2 + Q**2) / d%R
+        Y2 = (P**2 + Q**2) / d%R**2
     end if
-    d%rho = d%rho + d%fo(i) * Y**2
-    d%orbitals(:, i) = Y
+    d%rho = d%rho + d%fo(i) * Y2
+    d%orbitals(:, i) = Y2
 end do
 d%rho = d%rho / (4*pi)
 end subroutine

and create orbitalsP and orbitalsQ instead of orbitals that would return P and Q directly, so that the user has access to the original wavefunctions.