dftatom
dftatom copied to clipboard
Expose P, Q in orbitals
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. Even better, simply create:
d%orbitalsPQ(:, 1, i) = P
d%orbitalsPQ(:, 2, i) = Q
This can be done in the last iteration.