From 7f4e9e499151d27245c9d01d8d1814e0f9f20830 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Thu, 21 Sep 2023 14:55:16 -0600 Subject: [PATCH 01/13] Create solve_dirac_eigenproblem() function --- src/dirac.f90 | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 65 insertions(+), 1 deletion(-) diff --git a/src/dirac.f90 b/src/dirac.f90 index 4c22964..f39ff85 100644 --- a/src/dirac.f90 +++ b/src/dirac.f90 @@ -22,8 +22,72 @@ module dirac 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) +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(:,:), V(:,:), fullc(:), uq(:,:), rho0(:,:), rho1(:,:) +integer, intent(in) :: ib(:,:), in(:,:), focc_idx(:,Lmin:) +real(dp), intent(inout) :: D(:,:), S(:,:), H(:,:), lam(:), lam_tmp(:), eng(:) +integer, intent(out) :: idx +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) + 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 +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 From 0aeb587aad170ac17ddc94a16ddb2e72d306d496 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Thu, 21 Sep 2023 15:08:18 -0600 Subject: [PATCH 02/13] Call it --- src/dirac.f90 | 51 +++------------------------------------------------ 1 file changed, 3 insertions(+), 48 deletions(-) diff --git a/src/dirac.f90 b/src/dirac.f90 index f39ff85..21c9845 100644 --- a/src/dirac.f90 +++ b/src/dirac.f90 @@ -240,54 +240,9 @@ 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 - + 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) if ( .not. (size(eng) == idx) ) then error stop 'Size mismatch in energy array' end if From 58a793dafbd0c3f9492debc2d62b4d685d169129 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Tue, 26 Sep 2023 15:02:19 -0600 Subject: [PATCH 03/13] Use it in schroed_dirac_solver.f90 --- src/schroed_dirac_solver.f90 | 39 ++++-------------------------------- 1 file changed, 4 insertions(+), 35 deletions(-) diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index d84c6a6..2fe689b 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -110,41 +110,10 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & end do end do else - do kappa = -6, 5 - 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 - 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 - l = kappa - end if - do k = 1, 7-l - ind = (kappa+1)*(kappa+1)+(k-1)*(2*kappa+1)+(k-1)*(k-2) - if (kappa < 0) then - ind = ind + 2*k-2+(4*k-2)*(-1-kappa) - end if - 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 - end do - end do + ! TODO: generate focc + call solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_gj, & + xq, xq1, wtq_gj, V, Z, Vin, D, S, H, lam2, rho0, rho1, .false., fullc, & + ib, in, idx, lam_tmp, uq, wtq, xin, xiq, focc, focc_idx, lam, xq2) end if end subroutine total_energy From a5997afa2555ca4d80005af821c965c518b46196 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Tue, 26 Sep 2023 16:13:14 -0600 Subject: [PATCH 04/13] Get it to compile further --- src/dirac.f90 | 2 +- src/schroed_dirac_solver.f90 | 21 ++++++++++++--------- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/dirac.f90 b/src/dirac.f90 index 21c9845..00c1714 100644 --- a/src/dirac.f90 +++ b/src/dirac.f90 @@ -18,7 +18,7 @@ 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 diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index 2fe689b..7bb84ed 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -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 @@ -23,12 +23,13 @@ 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(:,:) + 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)) @@ -110,10 +111,12 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & end do end do else - ! TODO: generate focc - call solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_gj, & - xq, xq1, wtq_gj, V, Z, Vin, D, S, H, lam2, rho0, rho1, .false., fullc, & - ib, in, idx, lam_tmp, uq, wtq, xin, xiq, focc, focc_idx, lam, xq2) + ! TODO: allocate focc, focc_idx + Lmin2 = -6 + Lmax = 5 + call solve_dirac_eigenproblem(Nb, Nq, Lmin2, Lmax, alpha, alpha_j, xe, xiq1, & + xq, xq1, wtq1, V, Z, Vin, D, S, H, lam2, rho, rho1, .false., fullc, & + ib, in, idx, lam_tmp, uq, wtq, xin, xiq, focc, focc_idx, lam, xq) end if end subroutine total_energy From d3d822dabf9e9144979338923fb8486e156bc28e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Thu, 28 Sep 2023 14:43:12 -0600 Subject: [PATCH 05/13] Make it compile --- src/schroed_dirac_solver.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index 7bb84ed..ef26b94 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -26,8 +26,8 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & 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(:,:), & - focc(:,:), lam_tmp(:), rho1(:,:) - integer, allocatable :: in(:,:), ib(:,:), focc_idx(:) + 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, idx, Lmin2, Lmax real(dp) :: E_dirac_shift @@ -114,8 +114,10 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & ! TODO: allocate focc, focc_idx Lmin2 = -6 Lmax = 5 - call solve_dirac_eigenproblem(Nb, Nq, Lmin2, Lmax, alpha, alpha_j, xe, xiq1, & - xq, xq1, wtq1, V, Z, Vin, D, S, H, lam2, rho, rho1, .false., fullc, & + allocate(xiq_gj(size(xe),Lmin:Lmax)) + allocate(wtq_gj(size(xe),Lmin:Lmax)) + 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) end if end subroutine total_energy From c635884e7b72f6ab611ffb13fdaa3c4679936f07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Thu, 28 Sep 2023 14:49:49 -0600 Subject: [PATCH 06/13] Work on xq --- src/schroed_dirac_solver.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index ef26b94..f3cbf00 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -33,7 +33,7 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & 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)) xe = meshexp(rmin, rmax, a, Ne) call get_parent_quad_pts_wts(1, Nq, xiq, wtq) From 514faed43add9b918c18b600378e1566d9d19119 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Wed, 11 Oct 2023 11:52:03 +0200 Subject: [PATCH 07/13] Initialize xiq_gj and w corectly --- src/schroed_dirac_solver.f90 | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index f3cbf00..803ee98 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -114,8 +114,18 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & ! TODO: allocate focc, focc_idx Lmin2 = -6 Lmax = 5 - allocate(xiq_gj(size(xe),Lmin:Lmax)) - allocate(wtq_gj(size(xe),Lmin:Lmax)) + allocate(xiq_gj(size(xiq1),Lmin:Lmax)) + allocate(wtq_gj(size(wtq1),Lmin:Lmax)) + do kappa = -6, 5 + if (kappa == 0) cycle + if (alpha_j(kappa) > -1) then + call gauss_jacobi_gw(Nq, 0.0_dp, alpha_j(kappa), xiq1, wtq1) + xiq_gj(:,kappa) = xiq1 + wtq_gj(:,kappa) = wtq1 + end if + end do + allocate(Vin(size(xq,1), size(xq,2))) + Vin = 0 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) From 193fc6f23d2c687dffeb9b7b95fac2f3834f2ab6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Wed, 11 Oct 2023 12:03:57 +0200 Subject: [PATCH 08/13] Initialize focc --- src/schroed_dirac_solver.f90 | 29 +++++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index 803ee98..fef8533 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -111,21 +111,42 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & end do end do else - ! TODO: allocate focc, focc_idx Lmin2 = -6 Lmax = 5 + allocate(Vin(size(xq,1), size(xq,2))) + Vin = 0 allocate(xiq_gj(size(xiq1),Lmin:Lmax)) allocate(wtq_gj(size(wtq1),Lmin:Lmax)) - do kappa = -6, 5 + + ! 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) xiq_gj(:,kappa) = xiq1 wtq_gj(:,kappa) = wtq1 end if + if (kappa < 0) then + l = -kappa-1 + else + l = kappa + end if + do k = 1, 7-l + ind = (kappa+1)*(kappa+1)+(k-1)*(2*kappa+1)+(k-1)*(k-2) + if (kappa < 0) then + ind = ind + 2*k-2+(4*k-2)*(-1-kappa) + end if + if (kappa == -1) then + ind = ind + 1 + end if + focc_idx(k,kappa) = ind + focc(k,kappa) = 1 + end do end do - allocate(Vin(size(xq,1), size(xq,2))) - Vin = 0 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) From 871abe81393b945741ada51a5f9f02ef23cf67b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Wed, 11 Oct 2023 12:12:53 +0200 Subject: [PATCH 09/13] Allocate rho1 --- src/schroed_dirac_solver.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index fef8533..b95c3a0 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -117,6 +117,7 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & Vin = 0 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)) From c26c0f2466c26b3ac667578efb3e6552de3566f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Wed, 11 Oct 2023 13:37:54 +0200 Subject: [PATCH 10/13] Implement E_dirac_shift --- src/dirac.f90 | 15 +++++++++------ src/schroed_dirac_solver.f90 | 3 ++- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/src/dirac.f90 b/src/dirac.f90 index 00c1714..4084a0c 100644 --- a/src/dirac.f90 +++ b/src/dirac.f90 @@ -24,7 +24,7 @@ module dirac 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) + 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(:,:) @@ -35,6 +35,7 @@ subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_ integer, intent(in) :: ib(:,:), in(:,:), focc_idx(:,Lmin:) real(dp), intent(inout) :: D(:,:), S(:,:), H(:,:), lam(:), lam_tmp(:), eng(:) integer, intent(out) :: idx +real(dp), intent(in) :: E_dirac_shift integer :: kappa, i idx = 0 do kappa = Lmin, Lmax @@ -43,10 +44,10 @@ subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_ 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:) + 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 + V = Vin - Z/xq - E_dirac_shift endif call assemble_radial_dirac_SH(V, kappa, xin, xe, ib, xiq, wtq, & @@ -80,7 +81,7 @@ subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_ 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 + eng(focc_idx(i,kappa)) = sqrt(lam(i)) - c**2 + E_dirac_shift end do end do @@ -229,6 +230,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. @@ -240,9 +242,10 @@ subroutine Ffunc(x, y, eng) Vee = 0 rho0 = 0 rho1 = 0 + 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) + 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 diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index b95c3a0..6e4c906 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -150,7 +150,8 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & 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) + ib, in, idx, lam_tmp, uq, wtq, xin, xiq, focc, focc_idx, lam, xq, & + E_dirac_shift) end if end subroutine total_energy From 954761caffacd255184e41bbbcf1c0f008154cbc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Wed, 11 Oct 2023 14:11:52 +0200 Subject: [PATCH 11/13] Use Vin --- src/dirac.f90 | 5 +++-- src/schroed_dirac_solver.f90 | 7 ++++--- test/test_harmonic_dirac.f90 | 2 +- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/dirac.f90 b/src/dirac.f90 index 4084a0c..0fc99c0 100644 --- a/src/dirac.f90 +++ b/src/dirac.f90 @@ -31,9 +31,10 @@ subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_ 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(:,:), V(:,:), fullc(:), uq(:,:), rho0(:,:), rho1(:,:) +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 @@ -45,7 +46,7 @@ subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_ 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 + V(:,2:) = Vin(:,2:) - Z/xq(:,2:) - E_dirac_shift else V = Vin - Z/xq - E_dirac_shift endif diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index 6e4c906..624d935 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -72,11 +72,14 @@ 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 E_dirac_shift = 1000 end if if (dirac_int == 1) then @@ -113,8 +116,6 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & else Lmin2 = -6 Lmax = 5 - allocate(Vin(size(xq,1), size(xq,2))) - Vin = 0 allocate(xiq_gj(size(xiq1),Lmin:Lmax)) allocate(wtq_gj(size(wtq1),Lmin:Lmax)) allocate(rho1(Nq,Ne)) diff --git a/test/test_harmonic_dirac.f90 b/test/test_harmonic_dirac.f90 index cdde4c7..5b9870a 100644 --- a/test/test_harmonic_dirac.f90 +++ b/test/test_harmonic_dirac.f90 @@ -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))) From bab8d122604c31e0d4b9366023a142d546e4a0a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Wed, 11 Oct 2023 15:00:22 +0200 Subject: [PATCH 12/13] Get harmonic_dirac working --- test/test_harmonic_dirac.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_harmonic_dirac.f90 b/test/test_harmonic_dirac.f90 index 5b9870a..50dfc6d 100644 --- a/test/test_harmonic_dirac.f90 +++ b/test/test_harmonic_dirac.f90 @@ -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 From 00a29141f83dfb0b940cfc951507df854be231f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ond=C5=99ej=20=C4=8Cert=C3=ADk?= Date: Wed, 11 Oct 2023 15:03:19 +0200 Subject: [PATCH 13/13] Workaround for Schroedinger --- src/schroed_dirac_solver.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index 624d935..376a4ef 100644 --- a/src/schroed_dirac_solver.f90 +++ b/src/schroed_dirac_solver.f90 @@ -80,6 +80,7 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & E_dirac_shift = 500 else Vin = xq**2/2 + V = Vin E_dirac_shift = 1000 end if if (dirac_int == 1) then