Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor into solve_dirac_eigenproblem() #14

Merged
merged 13 commits into from
Oct 12, 2023
123 changes: 73 additions & 50 deletions src/dirac.f90
Original file line number Diff line number Diff line change
Expand Up @@ -18,12 +18,78 @@ module dirac
use iso_c_binding, only: c_double, c_int
implicit none
private
public solve_dirac, csolve_dirac
public solve_dirac, csolve_dirac, solve_dirac_eigenproblem

contains

subroutine solve_dirac(Z, p, xiq, wtq, xe, eps, energies, Etot, V, DOFs)
subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_gj, &
xq, xq1, wtq_gj, V, Z, Vin, D, S, H, lam, rho0, rho1, accurate_eigensolver, fullc, &
ib, in, idx, lam_tmp, uq, wtq, xin, xiq, focc, focc_idx, eng, xq2, E_dirac_shift)
logical, intent(in) :: accurate_eigensolver
integer, intent(in) :: Lmin, Lmax, Z, Nb, Nq
real(dp), intent(in) :: alpha_j(Lmin:), alpha(Lmin:), Vin(:,:), xq(:,:)
real(dp), intent(in) :: xe(:), xiq_gj(:,Lmin:), wtq_gj(:,Lmin:)
real(dp), intent(in) :: wtq(:), xin(:)
real(dp), intent(in) :: xiq(:), focc(:,Lmin:), xq2(:,:)
real(dp), intent(inout) :: xq1(:,:), fullc(:), uq(:,:), rho0(:,:), rho1(:,:)
integer, intent(in) :: ib(:,:), in(:,:), focc_idx(:,Lmin:)
real(dp), intent(inout) :: D(:,:), S(:,:), H(:,:), lam(:), lam_tmp(:), eng(:)
real(dp), intent(out) :: V(:,:)
integer, intent(out) :: idx
real(dp), intent(in) :: E_dirac_shift
integer :: kappa, i
idx = 0
do kappa = Lmin, Lmax
if (kappa == 0) cycle
!print *, "Calculating kappa =", kappa
if (alpha_j(kappa) > -1) then
call get_quad_pts(xe(:2), xiq_gj(:, kappa), xq1)
call proj_fn(Nq-1, xe(:2), xiq_gj(:,-1), wtq_gj(:,-1), xiq_gj(:, kappa), Vin, V(:,:1))
V(:,1) = V(:,1) - Z/xq1(:,1) - E_dirac_shift
V(:,2:) = Vin(:,2:) - Z/xq(:,2:) - E_dirac_shift
else
V = Vin - Z/xq - E_dirac_shift
endif

call assemble_radial_dirac_SH(V, kappa, xin, xe, ib, xiq, wtq, &
xiq_gj(:, kappa), wtq_gj(:, kappa), alpha(kappa), alpha_j(kappa), c, S, H)

! One can enforce symmetry using the following lines, it helps but
! we still need two seperate eigensolves for 1e-8 accuracy:
!H = (H + transpose(H))/2
!S = (S + transpose(S))/2

if (accurate_eigensolver) then
call eigh(H, S, lam)
call eigh(H, S, lam_tmp, D)
else
call eigh(H, S, lam, D)
end if

do i = 1, size(focc,1)
if (focc(i,kappa) < tiny(1._dp)) cycle

call c2fullc2(in, ib, D(:Nb,i), fullc)
call fe2quad(xe, xin, xiq, in, fullc, uq)
rho0 = rho0 - focc(i,kappa)*uq**2 * xq**alpha_j(kappa)
call fe2quad(xe, xin, xiq_gj(:,-1), in, fullc, uq)
rho1(:,1) = rho1(:,1) - focc(i,kappa)*uq(:,1)**2 * xq2(:,1)**alpha_j(kappa)

call c2fullc2(in, ib, D(Nb+1:,i), fullc)
call fe2quad(xe, xin, xiq, in, fullc, uq)
rho0 = rho0 - focc(i,kappa)*uq**2 * xq**alpha_j(kappa)
call fe2quad(xe, xin, xiq_gj(:,-1), in, fullc, uq)
rho1(:,1) = rho1(:,1) - focc(i,kappa)*uq(:,1)**2 * xq2(:,1)**alpha_j(kappa)

idx = idx + 1
eng(focc_idx(i,kappa)) = sqrt(lam(i)) - c**2 + E_dirac_shift
end do

end do
end subroutine


subroutine solve_dirac(Z, p, xiq, wtq, xe, eps, energies, Etot, V, DOFs)
integer, intent(in) :: Z, p
real(dp), intent(in) :: xe(:) ! element coordinates
real(dp), intent(in) :: xiq(:) ! quadrature points
Expand Down Expand Up @@ -165,6 +231,7 @@ subroutine Ffunc(x, y, eng)
! Converge Vee+Vxc only (the other components are constant)
real(dp), intent(in) :: x(:)
real(dp), intent(out) :: y(:), eng(:)
real(dp) :: E_dirac_shift
integer :: idx
logical :: accurate_eigensolver
accurate_eigensolver = .true.
Expand All @@ -176,54 +243,10 @@ subroutine Ffunc(x, y, eng)
Vee = 0
rho0 = 0
rho1 = 0
do kappa = Lmin, Lmax
if (kappa == 0) cycle
!print *, "Calculating kappa =", kappa
if (alpha_j(kappa) > -1) then
call get_quad_pts(xe(:2), xiq_gj(:, kappa), xq1)
call proj_fn(Nq-1, xe(:2), xiq_gj(:,-1), wtq_gj(:,-1), xiq_gj(:, kappa), Vin, V(:,:1))
V(:,1) = V(:,1) - Z/xq1(:,1)
V(:,2:) = Vin(:,2:) - Z/xq(:,2:)
else
V = Vin - Z/xq
endif

call assemble_radial_dirac_SH(V, kappa, xin, xe, ib, xiq, wtq, &
xiq_gj(:, kappa), wtq_gj(:, kappa), alpha(kappa), alpha_j(kappa), c, S, H)

! One can enforce symmetry using the following lines, it helps but
! we still need two seperate eigensolves for 1e-8 accuracy:
!H = (H + transpose(H))/2
!S = (S + transpose(S))/2

if (accurate_eigensolver) then
call eigh(H, S, lam)
call eigh(H, S, lam_tmp, D)
else
call eigh(H, S, lam, D)
end if

do i = 1, size(focc,1)
if (focc(i,kappa) < tiny(1._dp)) cycle

call c2fullc2(in, ib, D(:Nb,i), fullc)
call fe2quad(xe, xin, xiq, in, fullc, uq)
rho0 = rho0 - focc(i,kappa)*uq**2 * xq**alpha_j(kappa)
call fe2quad(xe, xin, xiq_gj(:,-1), in, fullc, uq)
rho1(:,1) = rho1(:,1) - focc(i,kappa)*uq(:,1)**2 * xq2(:,1)**alpha_j(kappa)

call c2fullc2(in, ib, D(Nb+1:,i), fullc)
call fe2quad(xe, xin, xiq, in, fullc, uq)
rho0 = rho0 - focc(i,kappa)*uq**2 * xq**alpha_j(kappa)
call fe2quad(xe, xin, xiq_gj(:,-1), in, fullc, uq)
rho1(:,1) = rho1(:,1) - focc(i,kappa)*uq(:,1)**2 * xq2(:,1)**alpha_j(kappa)

idx = idx + 1
eng(focc_idx(i,kappa)) = sqrt(lam(i)) - c**2
end do

end do

E_dirac_shift = 0
call solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_gj, &
xq, xq1, wtq_gj, V, Z, Vin, D, S, H, lam, rho0, rho1, accurate_eigensolver, fullc, &
ib, in, idx, lam_tmp, uq, wtq, xin, xiq, focc, focc_idx, eng, xq2, E_dirac_shift)
if ( .not. (size(eng) == idx) ) then
error stop 'Size mismatch in energy array'
end if
Expand Down
55 changes: 32 additions & 23 deletions src/schroed_dirac_solver.f90
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ module schroed_dirac_solver
use constants, only: pi
use mesh, only: meshexp
use schroed_glob, only: solve_schroed
use dirac, only: solve_dirac
use dirac, only: solve_dirac, solve_dirac_eigenproblem
use feutils, only: define_connect, get_quad_pts, get_parent_quad_pts_wts, &
get_parent_nodes, phih, dphih, phih_array, dphih_array, c2fullc2, &
fe2quad, fe2quad_core
Expand All @@ -23,16 +23,17 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
integer, intent(out) :: DOFs
real(dp), allocatable, intent(out) :: lam(:), eigfn(:,:,:)
real(dp), intent(out) :: xq(Nq, Ne)
real(dp), allocatable :: xe(:), xiq(:), wtq(:), V(:,:), xin(:), &
real(dp), allocatable :: xe(:), xiq(:), wtq(:), V(:,:), Vin(:,:), xin(:), &
H(:,:), S(:), S2(:,:), DSQ(:), phihq(:,:), dphihq(:,:), xiq_lob(:), wtq_lob(:), &
D(:,:), lam2(:), xiq1(:), wtq1(:), xq1(:,:), fullc(:), uq(:,:), rho(:,:)
integer, allocatable :: in(:,:), ib(:,:)
D(:,:), lam2(:), xiq1(:), wtq1(:), xq1(:,:), fullc(:), uq(:,:), rho(:,:), &
focc(:,:), lam_tmp(:), rho1(:,:), xiq_gj(:,:), wtq_gj(:,:)
integer, allocatable :: in(:,:), ib(:,:), focc_idx(:,:)
real(dp) :: rmin
integer :: al, i, j, l, k, ind, kappa, Nb, Nn
integer :: al, i, j, l, k, ind, kappa, Nb, Nn, idx, Lmin2, Lmax
real(dp) :: E_dirac_shift
rmin = 0
allocate(xe(Ne+1), xiq(Nq), xiq1(Nq), wtq(Nq), wtq1(Nq), V(Nq, Ne), xiq_lob(p+1), wtq_lob(p+1))
allocate(phihq(Nq, p+1), dphihq(Nq, p+1), in(p+1, Ne), ib(p+1, Ne), xin(p+1), xq1(Nq, Ne))
allocate(phihq(Nq, p+1), dphihq(Nq, p+1), in(p+1, Ne), ib(p+1, Ne), xin(p+1), xq1(Nq,1))
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

1 is correct here, consistent with DFT. Inside solve_dirac_eigenproblem only the first element is used. This should eventually be converted to just a 1D array.


HaoZeke marked this conversation as resolved.
Show resolved Hide resolved
xe = meshexp(rmin, rmax, a, Ne)
call get_parent_quad_pts_wts(1, Nq, xiq, wtq)
Expand Down Expand Up @@ -71,11 +72,15 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
allocate(eigfn(Nq, Ne, 28))
end if

allocate(Vin(size(xq,1), size(xq,2)))

if (potential_type == 0) then
V = -Z/xq
Vin = 0
E_dirac_shift = 500
else
V = xq**2/2
Vin = xq**2/2
V = Vin
E_dirac_shift = 1000
end if
if (dirac_int == 1) then
Expand Down Expand Up @@ -110,16 +115,24 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
end do
end do
else
do kappa = -6, 5
Lmin2 = -6
Lmax = 5
allocate(xiq_gj(size(xiq1),Lmin:Lmax))
allocate(wtq_gj(size(wtq1),Lmin:Lmax))
allocate(rho1(Nq,Ne))

! Initialize focc and focc_idx
allocate(focc(max(Lmax,abs(Lmin2))+1,Lmin2:Lmax))
allocate(focc_idx(max(Lmax,abs(Lmin2))+1,Lmin2:Lmax))
focc = 0
focc_idx = 0
do kappa = Lmin2, Lmax
if (kappa == 0) cycle
if (alpha_j(kappa) > -1) then
call gauss_jacobi_gw(Nq, 0.0_dp, alpha_j(kappa), xiq1, wtq1)
call get_quad_pts(xe(:2), xiq1, xq1)
V(:,:1) = -Z/xq1(:,:1)-E_dirac_shift
xiq_gj(:,kappa) = xiq1
wtq_gj(:,kappa) = wtq1
end if
call assemble_radial_dirac_SH(V, kappa, xin, xe, ib, xiq, wtq, xiq1, wtq1, &
alpha(kappa), alpha_j(kappa), c, S2, H)
call eigh(H, S2, lam2, D)
if (kappa < 0) then
l = -kappa-1
else
Expand All @@ -133,18 +146,14 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
if (kappa == -1) then
ind = ind + 1
end if
lam(ind) = sqrt(lam2(k)) - c**2
lam(ind) = lam(ind) + E_dirac_shift

call c2fullc2(in, ib, D(:Nb,k), fullc)
call fe2quad(xe, xin, xiq, in, fullc, uq)
eigfn(:,:,ind) = uq * xq**(alpha_j(kappa)/2)
if (alpha_j(kappa) > -1) then
call fe2quad(xe, xin, xiq1, in, fullc, uq)
eigfn(:,1,ind) = uq(:,1) * xq1(:,1)**(alpha_j(kappa)/2)
end if
focc_idx(k,kappa) = ind
focc(k,kappa) = 1
end do
end do
call solve_dirac_eigenproblem(Nb, Nq, Lmin2, Lmax, alpha, alpha_j, xe, xiq_gj, &
xq, xq1, wtq_gj, V, Z, Vin, D, S2, H, lam2, rho, rho1, .false., fullc, &
ib, in, idx, lam_tmp, uq, wtq, xin, xiq, focc, focc_idx, lam, xq, &
E_dirac_shift)
end if
end subroutine total_energy

Expand Down
6 changes: 3 additions & 3 deletions test/test_harmonic_dirac.f90
Original file line number Diff line number Diff line change
Expand Up @@ -145,8 +145,8 @@ program test_harmonic_dirac
Ne = p_or_Ne
end if

Lmax=6
Lmin=-7
Lmax=5
Lmin=-6

allocate(alpha(Lmin:Lmax), alpha_j(Lmin:Lmax))
do kappa = Lmin, Lmax
Expand Down Expand Up @@ -184,7 +184,7 @@ program test_harmonic_dirac
!end if
p = i
allocate(xq(Nq, Ne))
call total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
call total_energy(0, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
c, potential_type, Lmin, alpha_j, alpha, energies, eigfn, xq)
Etot = sum(energies)
allocate(energies_ref(size(energies)))
Expand Down
Loading