Coulomb benchmark

src/dirac.f90

@@ -37,7 +37,7 @@ subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_
 real(dp), intent(out) :: V(:,:)
 integer, intent(out) :: idx
 real(dp), intent(in) :: E_dirac_shift
-integer :: kappa, i
+integer :: kappa, i, nlam
 idx = 0
 do kappa = Lmin, Lmax
     if (kappa == 0) cycle
@@ -61,25 +61,27 @@ subroutine solve_dirac_eigenproblem(Nb, Nq, Lmin, Lmax, alpha, alpha_j, xe, xiq_
 
     if (accurate_eigensolver) then
         call eigh(H, S, lam)
-        call eigh(H, S, lam_tmp, D)
+        call eigh(H, S, lam_tmp, D, 7)
     else
-        call eigh(H, S, lam, D)
+        nlam = count(focc(:,kappa) > tiny(1._dp))
+        call eigh(H, S, lam, D, nlam)
     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)
+        !print *, "i=", i
+
+        !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
@@ -199,6 +201,7 @@ subroutine solve_dirac(Z, p, xiq, wtq, xe, eps, energies, Etot, V, DOFs)
     rho0(Nq,Ne), rho1(Nq,Ne))
 
 nband = count(focc > 0)
+!print *, "nband=", nband
 scf_max_iter = 100
 scf_alpha = 0.4_dp
 scf_L2_eps = 1e-4_dp

src/fe.f90

@@ -367,20 +367,6 @@ subroutine assemble_radial_dirac_SH(V, kappa, xin, xe, ib, xiq, wtq, xiq1, wtq1,
         end do
     end do
 end do
-
-!print *, "Checking symmetry"
-do j = 1, 2*Nb
-    do i = 1, j-1
-        if (abs(H(i,j) - H(j,i)) > 1e-8_dp) then
-            if (abs(H(i,j)-H(j,i)) / max(abs(H(i,j)), abs(H(j,i))) &
-                    > 1e-8_dp) then
-                print *, i, j, H(i,j)-H(j,i), H(i,j), H(j,i)
-                error stop "H not symmetric"
-            end if
-        end if
-        if (abs(S(i,j)-S(j,i)) > 1e-12_dp) error stop "S not symmetric"
-   end do
-end do
 end subroutine
 
 

src/lapack.f90

@@ -77,6 +77,18 @@ SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, &
     REAL(dp)           A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
     END SUBROUTINE
 
+    SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, &
+                       VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, &
+                       LWORK, IWORK, IFAIL, INFO )
+    import :: dp
+    CHARACTER          JOBZ, RANGE, UPLO
+    INTEGER            IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
+    REAL(dp)           ABSTOL, VL, VU
+    INTEGER            IFAIL( * ), IWORK( * )
+    REAL(dp)           A( LDA, * ), B( LDB, * ), W( * ), WORK( * ), &
+                       Z( LDZ, * )
+    END SUBROUTINE
+
     REAL(dp) FUNCTION DLAMCH( CMACH )
     import :: dp
     CHARACTER          CMACH

src/linalg.f90

@@ -1,6 +1,6 @@
 module linalg
   use types, only: dp
-  use lapack, only: dsyevd, dsygvd, dgesv
+  use lapack, only: dsyevd, dsygvd, dgesv, dsygvx
   implicit none
   private
   public eigh, solve
@@ -13,19 +13,22 @@ module linalg
 
 contains
 
-  subroutine deigh_generalized(Am, Bm, lam, c)
+  subroutine deigh_generalized(Am, Bm, lam, c, nlam)
     ! solves generalized eigen value problem for all eigenvalues and eigenvectors
     ! Am must by symmetric, Bm symmetric positive definite.
     ! Only the lower triangular part of Am and Bm is used.
     real(dp), intent(in) :: Am(:,:)   ! LHS matrix: Am c = lam Bm c
     real(dp), intent(in) :: Bm(:,:)   ! RHS matrix: Am c = lam Bm c
     real(dp), intent(out) :: lam(:)   ! eigenvalues: Am c = lam Bm c
     real(dp), intent(out) :: c(:,:)   ! eigenvectors: Am c = lam Bm c; c(i,j) = ith component of jth vec.
+    integer, intent(in) :: nlam
     integer :: n
     ! lapack variables
     integer :: lwork, liwork, info
-    integer, allocatable :: iwork(:)
-    real(dp), allocatable :: Bmt(:,:), work(:)
+    integer, allocatable :: iwork(:), ifail(:)
+    real(dp), allocatable :: Bmt(:,:), work(:), Z(:,:), Amt(:,:)
+    integer :: il, iu, M
+    real(dp) :: abstol
 
     ! solve
     n = size(Am,1)
@@ -34,9 +37,22 @@ subroutine deigh_generalized(Am, Bm, lam, c)
     call assert_shape(c, [n, n], "eigh", "c")
     lwork = 1 + 6*n + 2*n**2
     liwork = 3 + 5*n
-    allocate(Bmt(n,n), work(lwork), iwork(liwork))
-    c = Am; Bmt = Bm  ! Bmt temporaries overwritten by dsygvd
-    call dsygvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,iwork,liwork,info)
+    allocate(work(lwork), iwork(liwork))
+    allocate(ifail(n))
+    !call dsygvd(1,'V','L',n,c,n,Bmt,n,lam,work,lwork,iwork,liwork,info)
+    il = 1
+    iu = nlam
+    M = iu-il+1
+    allocate(z(n,M))
+    abstol = 1e-3_dp
+    call dsygvx(1,'V','I','L',n,Am,n,Bm,n, &
+        0._dp, 0._dp, il, iu, abstol, M, lam, c, n, work, &
+        lwork, iwork, ifail, info)
+    !SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, &
+    !                   LWORK, IWORK, LIWORK, INFO )
+    !SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, &
+    !                   VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, &
+    !                   LWORK, IWORK, IFAIL, INFO )
     if (info /= 0) then
        print *, "dsygvd returned info =", info
        if (info < 0) then

src/schroed_dirac_solver.f90

@@ -65,7 +65,7 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
     allocate(H(DOFS, DOFS), S(DOFS), DSQ(DOFS), uq(Nq,Ne))
     allocate(D(DOFS, DOFS), lam2(DOFS), rho(Nq,Ne), fullc(Nn))
     if (dirac_int == 1) then
-        allocate(lam(47))
+        allocate(lam(49))
         allocate(eigfn(Nq, Ne, 49))
     else
         allocate(lam(28))
@@ -115,18 +115,18 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
             end do
         end do
     else
-        Lmin2 = -6
-        Lmax = 5
+        Lmin2 = -7
+        Lmax = 6
         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))
+        allocate(focc(max(Lmax,abs(Lmin2))+2,Lmin2:Lmax))
+        allocate(focc_idx(max(Lmax,abs(Lmin2))+2,Lmin2:Lmax))
         focc = 0
         focc_idx = 0
-        do kappa = Lmin2, Lmax
+        end: 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)
@@ -138,18 +138,20 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
             else
                 l = kappa
             end if
-            do k = 1, 7-l
+            do k = 1, 8-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
+                if (ind > 49) cycle
+                !print *, kappa, l, k, ind
                 focc_idx(k,kappa) = ind
                 focc(k,kappa) = 1
             end do
-        end do
+        end do end
         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, &

test/test_coulomb_dirac.f90

@@ -49,8 +49,8 @@ program test_coulomb_dirac
 rmax = 50
 a = 200
 Ne = 4
-Nq = 53
-p = 26
+Nq = 20
+p = 1
 
 i = 7
 !call run_convergence_potential(0, 1, i, &
@@ -72,8 +72,8 @@ program test_coulomb_dirac
     Ne = p_or_Ne
 end if
 
-Lmax=5
-Lmin=-6
+Lmax=6
+Lmin=-7
 
 allocate(alpha(Lmin:Lmax), alpha_j(Lmin:Lmax))
 do kappa = Lmin, Lmax
@@ -109,7 +109,8 @@ program test_coulomb_dirac
 !if (dirac_int == 1 .and. i > 23) then
 !    exit
 !end if
-p = i
+p = 18
+Nq = p + 1
 allocate(xq(Nq, Ne))
 call total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, &
     c, potential_type, Lmin, alpha_j, alpha, energies, eigfn, xq)
@@ -143,7 +144,7 @@ program test_coulomb_dirac
 err = abs(Etot - Etot_ref)
 print "(f20.12, f20.12, es10.2)", Etot, Etot_ref, err
 if ( .not. (err < 5e-9_dp)) then
-   error stop 'assert failed'
+!   error stop 'assert failed'
 end if
 
 print *, "Eigenvalues:"
@@ -152,18 +153,10 @@ program test_coulomb_dirac
     err = abs(energies(i) - energies_ref(i))
     print "(i4, f20.12, f20.12, es10.2)", i, energies(i), energies_ref(i), err
     if ( .not. (err < 5e-9_dp)) then
-       error stop 'assert failed'
+!       error stop 'assert failed'
     end if
 end do
 
-print *, "Eigenfunctions saved in data_coulomb_dirac.txt"
-open(newunit=u, file="data_coulomb_dirac.txt", status="replace")
-write(u, *) xq
-do i = 1, size(energies)
-    write(u, *) eigfn(:,:,i)
-end do
-close(u)
-
 contains
 
     real(dp) function E_nl(c, n, l, Z, relat)