diff --git a/src/dirac.f90 b/src/dirac.f90 index 0fc99c0..f9aa6a1 100644 --- a/src/dirac.f90 +++ b/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 diff --git a/src/fe.f90 b/src/fe.f90 index e08afd3..ad8bc86 100644 --- a/src/fe.f90 +++ b/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 diff --git a/src/lapack.f90 b/src/lapack.f90 index 75fee7d..e2f64cc 100644 --- a/src/lapack.f90 +++ b/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 diff --git a/src/linalg.f90 b/src/linalg.f90 index a12ee3b..e5175ad 100644 --- a/src/linalg.f90 +++ b/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,7 +13,7 @@ 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. @@ -21,11 +21,14 @@ subroutine deigh_generalized(Am, Bm, lam, 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 diff --git a/src/schroed_dirac_solver.f90 b/src/schroed_dirac_solver.f90 index 376a4ef..f51f87c 100644 --- a/src/schroed_dirac_solver.f90 +++ b/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,7 +138,7 @@ 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) @@ -146,10 +146,12 @@ subroutine total_energy(Z, rmax, Ne, a, p, Nq, DOFs, alpha_int, dirac_int, & 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, & diff --git a/test/test_coulomb_dirac.f90 b/test/test_coulomb_dirac.f90 index 5d614a5..cd04ed4 100644 --- a/test/test_coulomb_dirac.f90 +++ b/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)