From dbe78402c45932d5e1f03a8afa59914c78c7553a Mon Sep 17 00:00:00 2001
From: Rudolf Zeller <ru.zeller@fz-juelich.de>
Date: Sat, 9 May 2020 12:07:52 +0200
Subject: [PATCH] Major change of Chebyshev solver to improve stability and
accuracy
---
source/KKRhost/rhooutnew.f90 | 28 +-
source/KKRhost/rhovalnew.F90 | 61 +++-
source/KKRhost/tmat_newsolver.F90 | 20 +-
.../create_newmesh.f90 | 6 +-
.../rll_global_solutions.F90 | 9 +-
.../sll_global_solutions.F90 | 268 ++++++++++++++++--
source/common/runoptions.F90 | 6 +
7 files changed, 347 insertions(+), 51 deletions(-)
diff --git a/source/KKRhost/rhooutnew.f90 b/source/KKRhost/rhooutnew.f90
index bdc7bf6eb..293f00d39 100644
--- a/source/KKRhost/rhooutnew.f90
+++ b/source/KKRhost/rhooutnew.f90
@@ -19,7 +19,7 @@ contains
!> Deprecated: False
!> Calculation of the density for the new solver
subroutine rhooutnew(nsra, lmax, gmatll, ek, lmpot, df, npan_tot, ncheb, cleb, icleb, iend, irmdnew, thetasnew, ifunm, imt1, lmsp, &
- rll, rllleft, sllleft, cden, cdenlm, cdenns, rho2nsc, corbital, gflle_part, rpan_intervall, ipan_intervall, nspin)
+ rll, ull, rllleft, sllleft, cden, cdenlm, cdenns, rho2nsc, corbital, gflle_part, rpan_intervall, ipan_intervall, nspin)
use :: mod_constants, only: cone,czero,pi
use :: mod_runoptions, only: calc_gmat_lm_full, use_ldau, set_cheby_nosoc
@@ -51,6 +51,7 @@ contains
real (kind=dp), dimension (0:ntotd), intent (in) :: rpan_intervall
real (kind=dp), dimension (ntotd*(ncheb+1), nfund), intent (in) :: thetasnew
complex (kind=dp), dimension (lmmaxd, lmmaxd), intent (in) :: gmatll !! GMATLL=diagonal elements of the G matrix (system) Note that SLL is not needed for calculation of density, only needed for calculation of Green function
+ complex (kind=dp), dimension (nsra*lmmaxd, lmmaxd, irmdnew), intent (in) :: ull
complex (kind=dp), dimension (nsra*lmmaxd, lmmaxd, irmdnew), intent (in) :: rll
complex (kind=dp), dimension (nsra*lmmaxd, lmmaxd, irmdnew), intent (in) :: rllleft
complex (kind=dp), dimension (nsra*lmmaxd, lmmaxd, irmdnew), intent (in) :: sllleft
@@ -131,18 +132,24 @@ contains
do lm2 = 1, lmmaxd
qnsi(lm1, lm2) = sllleft(lm1, lm2, ir)
! PNSI(LM1,LM2)=RLL(LM1,LM2,IR)
+ pnsi(lm1, lm2) = ull(lm1, lm2, ir)
+ end do
+ end do
+! call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, gmatll, lmmaxd, ek, qnsi, lmmaxd)
+ call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, ek, pnsi, lmmaxd, qnsi, lmmaxd, czero, wr(1,1,ir), lmmaxd)
+ do lm1 = 1, lmmaxd
+ do lm2 = 1, lmmaxd
pnsi(lm1, lm2) = rllleft(lm1, lm2, ir)
end do
end do
- ! CALL ZGEMM('N','N',lmmaxd,lmmaxd,lmmaxd,CONE,PNSI,
- ! + lmmaxd,GMATLL,lmmaxd,EK,QNSI,lmmaxd)
- call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, gmatll, lmmaxd, ek, qnsi, lmmaxd)
+ call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, gmatll, lmmaxd, czero, qnsi, lmmaxd)
do lm1 = 1, lmmaxd
do lm2 = 1, lmmaxd
pnsi(lm1, lm2) = rll(lm1, lm2, ir)
end do
end do
- call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, qnsi, lmmaxd, czero, wr(1,1,ir), lmmaxd)
+! call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, qnsi, lmmaxd, czero, wr(1,1,ir), lmmaxd)
+ call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, qnsi, lmmaxd, cone, wr(1,1,ir), lmmaxd)
! small component of Dirac spinor
if (nsra==2) then
@@ -151,12 +158,19 @@ contains
! QNSI(LM1,LM2)=SLLLEFT(LM1+lmmaxd,LM2,IR)
qnsi(lm1, lm2) = -sllleft(lm1+lmmaxd, lm2, ir)
! PNSI(LM1,LM2)=RLLLEFT(LM1+lmmaxd,LM2,IR)
- pnsi(lm1, lm2) = -rllleft(lm1+lmmaxd, lm2, ir)
+ pnsi(lm1, lm2) = ull(lm1+lmmaxd, lm2, ir)
end do
end do
! CALL ZGEMM('N','N',lmmaxd,lmmaxd,lmmaxd,CONE,PNSI,
! + lmmaxd,GMATLL,lmmaxd,EK,QNSI,lmmaxd)
- call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, gmatll, lmmaxd, ek, qnsi, lmmaxd)
+! call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, gmatll, lmmaxd, ek, qnsi, lmmaxd)
+ call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, ek, pnsi, lmmaxd, qnsi, lmmaxd, cone, wr(1,1,ir), lmmaxd)
+ do lm1 = 1, lmmaxd
+ do lm2 = 1, lmmaxd
+ pnsi(lm1, lm2) = -rllleft(lm1+lmmaxd, lm2, ir)
+ end do
+ end do
+ call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, gmatll, lmmaxd, czero, qnsi, lmmaxd)
do lm1 = 1, lmmaxd
do lm2 = 1, lmmaxd
pnsi(lm1, lm2) = rll(lm1+lmmaxd, lm2, ir)
diff --git a/source/KKRhost/rhovalnew.F90 b/source/KKRhost/rhovalnew.F90
index 24dbae836..eae7166e2 100644
--- a/source/KKRhost/rhovalnew.F90
+++ b/source/KKRhost/rhovalnew.F90
@@ -119,7 +119,7 @@ contains
real (kind=dp) :: thetanew, phinew
real (kind=dp) :: totmoment
- real (kind=dp) :: totxymoment
+ real (kind=dp) :: totxymoment, rsum
complex (kind=dp) :: ek
complex (kind=dp) :: df
complex (kind=dp) :: eryd
@@ -172,7 +172,10 @@ contains
complex (kind=dp), dimension (:, :, :, :), allocatable :: vnspll
complex (kind=dp), dimension (:, :, :, :), allocatable :: r2orbc
complex (kind=dp), dimension (:, :, :, :), allocatable :: vnspll1
+ complex (kind=dp), dimension (:, :, :, :), allocatable :: ull
+ complex (kind=dp), dimension (:, :, :, :), allocatable :: ullleft
complex (kind=dp), dimension (:, :, :, :), allocatable :: rllleft
+
complex (kind=dp), dimension (:, :, :, :), allocatable :: sllleft
complex (kind=dp), dimension (:, :, :, :), allocatable :: r2nefc_loop
complex (kind=dp), dimension (:, :, :, :), allocatable :: rho2nsc_loop
@@ -281,6 +284,12 @@ contains
allocate (sll(nsra*lmmaxd,lmmaxd,irmdnew,0:nth-1), stat=i_stat)
call memocc(i_stat, product(shape(sll))*kind(sll), 'SLL', 'RHOVALNEW')
sll(:,:,:,:) = czero
+ allocate (ull(nsra*lmmaxd,lmmaxd,irmdnew,0:nth-1), stat=i_stat)
+ call memocc(i_stat, product(shape(ull))*kind(ull), 'ULL', 'RHOVALNEW')
+ ull(:,:,:,:) = czero
+ allocate (ullleft(nsra*lmmaxd,lmmaxd,irmdnew,0:nth-1), stat=i_stat)
+ call memocc(i_stat, product(shape(ullleft))*kind(ullleft), 'ULLLEFT', 'RHOVALNEW')
+ ullleft(:,:,:,:) = czero
allocate (rllleft(nsra*lmmaxd,lmmaxd,irmdnew,0:nth-1), stat=i_stat)
call memocc(i_stat, product(shape(rllleft))*kind(rllleft), 'RLLLEFT', 'RHOVALNEW')
rllleft(:,:,:,:) = czero
@@ -523,8 +532,10 @@ contains
! no irregular solutions SLL are needed in self-consistent iterations
! because the density depends only on RLL, RLLLEFT and SLLLEFT
if (.not.set_cheby_nospeedup .and. .not. (calc_exchange_couplings .or. write_pkkr_operators)) then
- call rll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith), rll(:,:,:,ith), tmatll, ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, jlk_index, &
+ call rll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith), ull(:,:,:,ith), rll(:,:,:,ith), tmatll, ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, jlk_index, &
hlk(:,:,ith), jlk(:,:,ith), hlk2(:,:,ith), jlk2(:,:,ith), gmatprefactor, '1', use_sratrick, alphall)
+! call rll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith), rll(:,:,:,ith), tmatll, ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, jlk_index, &
+! hlk(:,:,ith), jlk(:,:,ith), hlk2(:,:,ith), jlk2(:,:,ith), gmatprefactor, '1', use_sratrick, alphall)
else
call rllsll(rpan_intervall, rnew, vnspll(:,:,:,ith), rll(:,:,:,ith), sll(:,:,:,ith), tmatll, ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, jlk_index, &
hlk(:,:,ith), jlk(:,:,ith), hlk2(:,:,ith), jlk2(:,:,ith), gmatprefactor, '1', '1', '0', use_sratrick, alphall)
@@ -543,7 +554,7 @@ contains
! read/recalc wavefunctions left contruct the TRANSPOSE spin-orbit coupling hamiltonian and add to potential
if ( .not. set_cheby_nosoc) then
call spinorbit_ham(lmax, lmmax0d, vins, rnew, eryd, zat, cvlight, socscale, nspin, lmpotd, theta, phi, ipan_intervall, rpan_intervall, npan_tot, ncheb, irmdnew, nrmaxd, &
- vnspll0, vnspll1(:,:,:,ith), 'transpose')
+ vnspll0, vnspll1(:,:,:,ith), 'rll')
else
vnspll1(:,:,:,ith) = vnspll0(:,:,:)
end if
@@ -579,6 +590,7 @@ contains
end if
! calculate the tmat and wavefunctions
+ ullleft(:, :, :, ith) = czero
rllleft(:, :, :, ith) = czero
sllleft(:, :, :, ith) = czero
@@ -588,8 +600,8 @@ contains
alphall = czero
! faster calculation of RLLLEFT and SLLLEFT.
if (.not.set_cheby_nospeedup .and. .not. (calc_exchange_couplings .or. write_pkkr_operators)) then
- call rll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith), rllleft(:,:,:,ith), tmattemp, ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, &
- jlk_index, hlk2(:,:,ith), jlk2(:,:,ith), hlk(:,:,ith), jlk(:,:,ith), gmatprefactor, '1', use_sratrick, alphall)
+ call rll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith),ullleft(:,:,:,ith), rllleft(:,:,:,ith), tmattemp, ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, &
+ jlk_index, hlk2(:,:,ith), jlk2(:,:,ith), hlk(:,:,ith), jlk(:,:,ith), gmatprefactor, '1', use_sratrick, alphall)
call sll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith), sllleft(:,:,:,ith), ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, jlk_index, &
hlk2(:,:,ith), jlk2(:,:,ith), hlk(:,:,ith), jlk(:,:,ith), gmatprefactor, '1', use_sratrick)
else
@@ -628,9 +640,12 @@ contains
gmatll(lm1, lm2, ie) = gmat0(lm1, lm2)
end do
end do
+
! calculate density
+
call rhooutnew(nsra, lmax, gmatll(1,1,ie), ek, lmpotd, df, npan_tot, ncheb, cleb, icleb, iend, irmdnew, thetasnew, ifunm, imt1, lmsp, rll(:,:,:,ith), & ! SLL(:,:,:,ith), commented out since sll is not used in rhooutnew
- rllleft(:,:,:,ith), sllleft(:,:,:,ith), cden(:,:,:,ith), cdenlm(:,:,:,ith), cdenns(:,:,ith), rho2nsc_loop(:,:,:,ie), 0, gflle(:,:,ie,iq), rpan_intervall, ipan_intervall, nspin/(nspin-korbit))
+ ull(:,:,:,ith), rllleft(:,:,:,ith), sllleft(:,:,:,ith), cden(:,:,:,ith), cdenlm(:,:,:,ith), cdenns(:,:,ith), rho2nsc_loop(:,:,:,ie), 0, gflle(:,:,ie,iq), rpan_intervall, ipan_intervall, nspin/(nspin-korbit))
+
do jspin = 1, nspin/(nspin-korbit)*(1+korbit)
do lm1 = 0, lmax
@@ -677,7 +692,7 @@ contains
!------------------------------------------------------------------------------
if (ie==ielast .and. ldorhoef) then
call rhooutnew(nsra, lmax, gmatll(1,1,ie), ek, lmpotd, cone, npan_tot, ncheb, cleb, icleb, iend, irmdnew, thetasnew, ifunm, imt1, lmsp, rll(:,:,:,ith), & ! SLL(:,:,:,ith), ! commented out since sll is not used in rhooutnew
- rllleft(:,:,:,ith), sllleft(:,:,:,ith), cden(:,:,:,ith), cdenlm(:,:,:,ith), cdenns(:,:,ith), r2nefc_loop(:,:,:,ith), 0, gflle_part(:,:,ith), rpan_intervall, &
+ ull(:,:,:,ith), rllleft(:,:,:,ith), sllleft(:,:,:,ith), cden(:,:,:,ith), cdenlm(:,:,:,ith), cdenns(:,:,ith), r2nefc_loop(:,:,:,ith), 0, gflle_part(:,:,ith), rpan_intervall, &
ipan_intervall, nspin/(nspin-korbit))
end if
@@ -687,7 +702,7 @@ contains
if (.not. set_cheby_nosoc) then
do iorb = 1, 3
call rhooutnew(nsra, lmax, gmatll(1,1,ie), ek, lmpotd, cone, npan_tot, ncheb, cleb, icleb, iend, irmdnew, thetasnew, ifunm, imt1, lmsp, rll(:,:,:,ith), & ! SLL(:,:,:,ith), ! commented out since sll is not used in rhooutnew
- rllleft(:,:,:,ith), sllleft(:,:,:,ith), cden(:,:,:,ith), cdenlm(:,:,:,ith), cdenns(:,:,ith), r2orbc(:,:,:,ith), iorb, gflle_part(:,:,ith), rpan_intervall, ipan_intervall, nspin)
+ ull(:,:,:,ith), rllleft(:,:,:,ith), sllleft(:,:,:,ith), cden(:,:,:,ith), cdenlm(:,:,:,ith), cdenns(:,:,ith), r2orbc(:,:,:,ith), iorb, gflle_part(:,:,ith), rpan_intervall, ipan_intervall, nspin)
do jspin = 1, nspin*(1+korbit)
if (jspin<=2) then
do lm1 = 0, lmax
@@ -949,6 +964,30 @@ contains
allocate (rhonewtemp(irws,lmpotd), stat=i_stat)
call memocc(i_stat, product(shape(rhonewtemp))*kind(rhonewtemp), 'RHONEWTEMP', 'RHOVALNEW')
+ open (131,file='rho2',form='formatted')
+ write(131,fmt='(A,I8)') '# ',irmdnew
+ do jspin = 1, nspin/(nspin-korbit)*(1+korbit)
+ do lm1 = 1, lmpotd
+ rsum = 0.0d0
+ do ir = 1, irmdnew
+ rsum = rsum + abs(aimag(rho2nsc(ir, lm1, jspin)))
+ end do
+ if(rsum.gt.1.d-10) then
+ write(131,fmt='(A,I8)') '# ',lm1
+ do ir = 1, irmdnew
+ if(abs(thetasnew(ir,1)).gt.1.d-6) then
+ write(131,*) rnew(ir),aimag(rho2nsc(ir, lm1, jspin))/rnew(ir)**2,aimag(rho2nsc(ir, lm1, jspin))/rnew(ir)**2*thetasnew(ir,1)
+ else
+ write(131,*) rnew(ir),aimag(rho2nsc(ir, lm1, jspin))/rnew(ir)**2,aimag(rho2nsc(ir, lm1, jspin))/rnew(ir)**2*3.544907701811032d0
+ end if
+ end do
+ end if
+ end do
+ end do
+ lm1 =999999
+ write(131,fmt='(A,I8)') '# ',lm1
+ close (131)
+
do jspin = 1, nspin/(nspin-korbit)*(1+korbit)
rhotemp = czero
rhonewtemp = czero
@@ -1074,12 +1113,18 @@ contains
i_all = -product(shape(rll))*kind(rll)
deallocate (rll, stat=i_stat)
call memocc(i_stat, i_all, 'RLL', 'RHOVALNEW')
+ i_all = -product(shape(ull))*kind(ull)
+ deallocate (ull, stat=i_stat)
+ call memocc(i_stat, i_all, 'ULL', 'RHOVALNEW')
i_all = -product(shape(sll))*kind(sll)
deallocate (sll, stat=i_stat)
call memocc(i_stat, i_all, 'SLL', 'RHOVALNEW')
i_all = -product(shape(rllleft))*kind(rllleft)
deallocate (rllleft, stat=i_stat)
call memocc(i_stat, i_all, 'RLLLEFT', 'RHOVALNEW')
+ i_all = -product(shape(ullleft))*kind(ullleft)
+ deallocate (ullleft, stat=i_stat)
+ call memocc(i_stat, i_all, 'ULLLEFT', 'RHOVALNEW')
i_all = -product(shape(sllleft))*kind(sllleft)
deallocate (sllleft, stat=i_stat)
call memocc(i_stat, i_all, 'SLLLEFT', 'RHOVALNEW')
diff --git a/source/KKRhost/tmat_newsolver.F90 b/source/KKRhost/tmat_newsolver.F90
index 1fe1d7e1e..1e293bd45 100644
--- a/source/KKRhost/tmat_newsolver.F90
+++ b/source/KKRhost/tmat_newsolver.F90
@@ -121,6 +121,7 @@ contains
complex (kind=dp), dimension (:, :, :), allocatable :: hlk2
complex (kind=dp), dimension (:, :, :), allocatable :: jlk2
complex (kind=dp), dimension (:, :, :), allocatable :: vnspll0
+ complex (kind=dp), dimension (:, :, :, :), allocatable :: ull !! regular solution of radial equation
complex (kind=dp), dimension (:, :, :, :), allocatable :: rll !! regular solution of radial equation
complex (kind=dp), dimension (:, :, :, :), allocatable :: sll !! irregular solution of radial equation
complex (kind=dp), dimension (:, :, :, :), allocatable :: vnspll
@@ -186,7 +187,7 @@ contains
! .. allocate and initialize arrays
call allocate_locals_tmat_newsolver(1, irmdnew, lmpot, nspin/(nspin-korbit), vins, aux, ipiv, tmat0, tmatll, alpha0, dtmatll, alphall, dalphall, jlk_index, nsra, lmmaxd, nth, lmax, vnspll, &
- vnspll0, vnspll1, hlk, jlk, hlk2, jlk2, tmatsph, rll, sll, rllleft, sllleft)
+ vnspll0, vnspll1, hlk, jlk, hlk2, jlk2, tmatsph, ull, rll, sll, rllleft, sllleft)
vins(1:irmdnew, 1:lmpot, 1) = vinsnew(1:irmdnew, 1:lmpot, ipot)
if (.not.set_cheby_nosoc) vins(1:irmdnew, 1:lmpot, nspin) = vinsnew(1:irmdnew, 1:lmpot, ipot+nspin-1)
@@ -415,6 +416,7 @@ contains
else
use_fullgmat = 1
end if
+!RZ
call rllsllsourceterms(nsra, nvec, eryd, rnew, irmdnew, nrmaxd, lmax, lmmaxd, use_fullgmat, jlk_index, hlk(:,:,ith), jlk(:,:,ith), hlk2(:,:,ith), jlk2(:,:,ith), gmatprefactor)
#ifdef CPP_OMP
@@ -494,8 +496,10 @@ contains
! no irregular solutions are needed in self-consistent iterations
! because the t-matrix depends only on RLL
if (.not. set_cheby_nospeedup .and. .not. (calc_exchange_couplings .or. write_pkkr_operators) .and. .not.calc_wronskian) then
- call rll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith), rll(:,:,:,ith), tmat0(:,:), ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, &
- jlk_index, hlk(:,:,ith), jlk(:,:,ith), hlk2(:,:,ith), jlk2(:,:,ith), gmatprefactor, '1', use_sratrick, alpha0(:,:))
+! call rll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith), ull(:,:,:,ith), rll(:,:,:,ith), tmat0(:,:), ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, &
+! jlk_index, hlk(:,:,ith), jlk(:,:,ith), hlk2(:,:,ith), jlk2(:,:,ith), gmatprefactor, '1', use_sratrick, alpha0(:,:))
+ call rll_global_solutions(rpan_intervall, rnew, vnspll(:,:,:,ith), ull(:,:,:,ith), rll(:,:,:,ith), tmatll, ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, jlk_index, &
+ hlk(:,:,ith), jlk(:,:,ith), hlk2(:,:,ith), jlk2(:,:,ith), gmatprefactor, '1', use_sratrick, alphall)
else
call rllsll(rpan_intervall, rnew, vnspll(:,:,:,ith), rll(:,:,:,ith), sll(:,:,:,ith), tmat0(:,:), ncheb, npan_tot, lmmaxd, nvec*lmmaxd, nsra*(1+korbit)*(lmax+1), irmdnew, nsra, &
jlk_index, hlk(:,:,ith), jlk(:,:,ith), hlk2(:,:,ith), jlk2(:,:,ith), gmatprefactor, '1', '1', '0', use_sratrick, alpha0(:,:))
@@ -826,7 +830,7 @@ contains
! deallocate arrays
call allocate_locals_tmat_newsolver(-1,irmdnew,lmpot,nspin,vins,aux,ipiv,tmat0, &
tmatll,alpha0,dtmatll,alphall,dalphall,jlk_index,nsra,lmmaxd,nth,lmax,vnspll,&
- vnspll0,vnspll1,hlk,jlk,hlk2,jlk2,tmatsph,rll,sll,rllleft,sllleft)
+ vnspll0,vnspll1,hlk,jlk,hlk2,jlk2,tmatsph,ull,rll,sll,rllleft,sllleft)
end subroutine tmat_newsolver
@@ -842,7 +846,7 @@ contains
!-------------------------------------------------------------------------------
subroutine allocate_locals_tmat_newsolver(allocmode,irmdnew,lmpot,nspin,vins,aux, &
ipiv,tmat0,tmatll,alpha0,dtmatll,alphall,dalphall,jlk_index,nsra,lmmaxd,nth, &
- lmax,vnspll,vnspll0,vnspll1,hlk,jlk,hlk2,jlk2,tmatsph,rll,sll,rllleft,sllleft)
+ lmax,vnspll,vnspll0,vnspll1,hlk,jlk,hlk2,jlk2,tmatsph,ull,rll,sll,rllleft,sllleft)
use :: mod_datatypes, only: dp
use :: mod_runoptions, only: calc_exchange_couplings, write_rhoq_input, calc_wronskian
use :: mod_constants, only: czero
@@ -873,6 +877,7 @@ contains
complex (kind=dp), allocatable, dimension (:, :, :), intent (inout) :: jlk2
complex (kind=dp), allocatable, dimension (:, :, :), intent (inout) :: hlk2
complex (kind=dp), allocatable, dimension (:, :), intent (inout) :: tmatsph
+ complex (kind=dp), allocatable, dimension (:, :, :, :), intent (inout) :: ull
complex (kind=dp), allocatable, dimension (:, :, :, :), intent (inout) :: rll
complex (kind=dp), allocatable, dimension (:, :, :, :), intent (inout) :: sll
complex (kind=dp), allocatable, dimension (:, :, :, :), intent (inout) :: rllleft
@@ -916,6 +921,9 @@ contains
tmatsph = czero
! Regular and irregular wavefunctions
+ allocate (ull(nsra*lmmaxd,lmmaxd,irmdnew,0:nth-1), stat=i_stat)
+ call memocc(i_stat, product(shape(ull))*kind(ull), 'ULL', 'allocate_locals_tmat_newsolver')
+ ull = czero
allocate (rll(nsra*lmmaxd,lmmaxd,irmdnew,0:nth-1), stat=i_stat)
call memocc(i_stat, product(shape(rll))*kind(rll), 'RLL', 'allocate_locals_tmat_newsolver')
rll = czero
@@ -997,6 +1005,8 @@ contains
deallocate (tmatsph, stat=i_stat)
call memocc(i_stat, -product(shape(tmatsph))*kind(tmatsph), 'TMATSPH', 'allocate_locals_tmat_newsolver')
+ deallocate (ull, stat=i_stat)
+ call memocc(i_stat, -product(shape(ull))*kind(ull), 'ULL', 'allocate_locals_tmat_newsolver')
deallocate (rll, stat=i_stat)
call memocc(i_stat, -product(shape(rll))*kind(rll), 'RLL', 'allocate_locals_tmat_newsolver')
deallocate (sll, stat=i_stat)
diff --git a/source/common/radial_solver_Chebychev/create_newmesh.f90 b/source/common/radial_solver_Chebychev/create_newmesh.f90
index 7c5f961dd..8ce134ce8 100644
--- a/source/common/radial_solver_Chebychev/create_newmesh.f90
+++ b/source/common/radial_solver_Chebychev/create_newmesh.f90
@@ -71,7 +71,8 @@ contains
real (kind=dp), dimension (ntotd*(ncheb+1), nfund, ncelld), intent (inout), optional :: thetasnew !! interpolated shape function in Chebychev radial mesh
! .. Local variables
- real (kind=dp), parameter :: fac=2e0_dp
+! real (kind=dp), parameter :: fac=2e0_dp
+ real (kind=dp), parameter :: fac=1.2e0_dp
integer :: npan_inst, i_stat, i_all
integer :: i1, ir2, ip, icell
integer :: imin, imax, iminnew, imaxnew, lm1
@@ -178,7 +179,10 @@ contains
call interpolspline(rmesh(imin:imax,i1), rnew(iminnew:imaxnew,i1), thetasin(imin-ircut(1,i1):imax-ircut(1,i1),lm1,icell), thetasnew(iminnew:imaxnew,lm1,icell), imax-imin+1, &
imaxnew-iminnew+1)
end do
+
end do
+
+
end if ! present(ntcell)
end do ! i1
diff --git a/source/common/radial_solver_Chebychev/rll_global_solutions.F90 b/source/common/radial_solver_Chebychev/rll_global_solutions.F90
index 2e114b1d2..e64d55514 100644
--- a/source/common/radial_solver_Chebychev/rll_global_solutions.F90
+++ b/source/common/radial_solver_Chebychev/rll_global_solutions.F90
@@ -20,7 +20,7 @@ contains
!> Deprecated: False
!> Wrapper for the calculation of the regular solutions for the impurity code `KKRimp`
!-------------------------------------------------------------------------------
- subroutine rll_global_solutions(rpanbound,rmesh,vll,rll,tllp,ncheb,npan,lmsize, &
+ subroutine rll_global_solutions(rpanbound,rmesh,vll,ull,rll,tllp,ncheb,npan,lmsize, &
lmsize2,lbessel,nrmax,nvec,jlk_index,hlk,jlk,hlk2,jlk2,gmatprefactor,cmoderll, &
use_sratrick1,alphaget) ! LLY
! ************************************************************************
@@ -61,12 +61,12 @@ contains
! cmoderll ="1" : op( )=identity
! cmoderll ="T" : op( )=transpose in L
-
+
+ complex (kind=dp) :: ull(lmsize2, lmsize, nrmax)
complex (kind=dp) :: rll(lmsize2, lmsize, nrmax), & ! reg. fredholm sol.
tllp(lmsize, lmsize), & ! t-matrix
vll(lmsize*nvec, lmsize*nvec, nrmax) ! potential term in 5.7
! bauer, phd
- complex (kind=dp), allocatable :: ull(:, :, :) ! reg. volterra sol.
complex (kind=dp), allocatable :: work(:, :), allp(:, :, :), bllp(:, :, :), & ! eq. 5.9, 5.10 for reg. sol
mrnvy(:, :, :), mrnvz(:, :, :), & ! ***************
@@ -140,8 +140,6 @@ contains
allocate (work(lmsize,lmsize))
allocate (allp(lmsize,lmsize,0:npan), bllp(lmsize,lmsize,0:npan))
- allocate (ull(lmsize2,lmsize,nrmax))
-
if (idotime==1) call timing_start('local')
#ifdef CPP_HYBRID
@@ -240,6 +238,7 @@ contains
alphaget(lm1, lm2) = allp(lm1, lm2, npan) ! LLY
end do ! LLY
end do ! LLY
+
! calculation of the t-matrix ! calc t-matrix tll = bll*alpha^-1
call zgemm('n','n',lmsize,lmsize,lmsize,cone/gmatprefactor,bllp(1,1,npan), &
lmsize,allp(1,1,npan),lmsize,czero,tllp,lmsize)
diff --git a/source/common/radial_solver_Chebychev/sll_global_solutions.F90 b/source/common/radial_solver_Chebychev/sll_global_solutions.F90
index cca49096c..27dfa572e 100644
--- a/source/common/radial_solver_Chebychev/sll_global_solutions.F90
+++ b/source/common/radial_solver_Chebychev/sll_global_solutions.F90
@@ -33,6 +33,7 @@ contains
use :: mod_constants, only: czero, cone
use :: mod_datatypes, only: dp
+ use :: mod_runoptions, only: use_cheby_quadprec
use :: mod_chebint, only: chebint
use :: mod_sll_local_solutions, only: sll_local_solutions
@@ -48,7 +49,9 @@ contains
! running indices
integer :: lm1, lm2
- integer :: icheb, ipan, mn
+ integer :: icheb, ipan, mn, nm
+ integer :: info, ipiv(lmsize)
+ integer :: iter_beta, niter_beta
! source terms
complex (kind=dp) :: gmatprefactor ! prefactor of green function
@@ -61,15 +64,26 @@ contains
! cmodesll ="1" : op( )=identity
! cmodesll ="T" : op( )=transpose in L
+
+ complex*32, allocatable :: qcllp(:, :, :), qdllp(:, :, :)
+ complex*32, allocatable :: qmihvy(:, :), qmihvz(:, :), qmijvy(:, :), qmijvz(:, :)
+ complex*32, allocatable :: qyif(:, :, :)
+ complex*32, allocatable :: qcllptemp(:, :), qdllptemp(:, :)
+ complex*32, allocatable :: qsll(:, :)
+ complex*32, allocatable :: qcone, qczero
+ complex*32, allocatable :: qbetainv(:, :), qbetainv_save(:, :)
complex (kind=dp) :: sll(lmsize2, lmsize, nrmax), & ! irr. volterra sol.
vll(lmsize*nvec, lmsize*nvec, nrmax) ! potential term in 5.7
! bauer, phd
+ real (kind=dp) :: dllpmax,dllpval
+ complex (kind=dp), allocatable :: cllptemp(:, :), dllptemp(:, :)
complex (kind=dp), allocatable :: work(:, :), cllp(:, :, :), dllp(:, :, :), mihvy(:, :, :), mihvz(:, :, :), &
mijvy(:, :, :), mijvz(:, :, :) ! ***************
- complex (kind=dp), allocatable :: yif(:, :, :, :), & ! source terms (different array
- zif(:, :, :, :)
+ complex (kind=dp), allocatable :: yif(:, :, :, :), zif(:, :, :, :)
+ complex (kind=dp), allocatable :: betainv(:, :), betainv_save(:, :)
+
! chebyshev arrays
real (kind=dp) :: c1(0:ncheb, 0:ncheb), rpanbound(0:npan), drpan2
real (kind=dp) :: cslc1(0:ncheb, 0:ncheb), & ! Integration matrix from left ( C*S_L*C^-1 in eq. 5.53)
@@ -117,6 +131,14 @@ contains
#endif
if (idotime==1) call timing_start('sll-glob')
+ allocate (mihvy(lmsize,lmsize,npan), mihvz(lmsize,lmsize,npan))
+ allocate (mijvy(lmsize,lmsize,npan), mijvz(lmsize,lmsize,npan))
+ allocate (yif(lmsize2,lmsize,0:ncheb,npan))
+ allocate (zif(lmsize2,lmsize,0:ncheb,npan))
+ allocate (betainv(lmsize,lmsize), betainv_save(lmsize,lmsize))
+ allocate (cllp(lmsize,lmsize,0:npan), dllp(lmsize,lmsize,0:npan))
+ allocate (cllptemp(lmsize,lmsize), dllptemp(lmsize,lmsize))
+ allocate (work(lmsize,lmsize))
do ipan = 1, npan
do icheb = 0, ncheb
@@ -127,14 +149,6 @@ contains
call chebint(cslc1, csrc1, slc1sum, c1, ncheb)
- allocate (mihvy(lmsize,lmsize,npan), mihvz(lmsize,lmsize,npan))
- allocate (mijvy(lmsize,lmsize,npan), mijvz(lmsize,lmsize,npan))
- allocate (yif(lmsize2,lmsize,0:ncheb,npan))
- allocate (zif(lmsize2,lmsize,0:ncheb,npan))
-
- allocate (work(lmsize,lmsize))
- allocate (cllp(lmsize,lmsize,0:npan), dllp(lmsize,lmsize,0:npan))
-
if (idotime==1) call timing_start('local')
#ifdef CPP_HYBRID
@@ -172,36 +186,151 @@ contains
if (idotime==1) call timing_start('afterlocal')
! ***********************************************************************
- ! calculate A(M), B(M), C(M), D(M)
- ! according to 5.17-5.18 (regular solution) of Bauer PhD
- ! C,D are calculated accordingly for the irregular solution
- ! (starting condition: A(0) = 1, B(0) = 0, C(MMAX) = 0 and D(MMAX) = 1)
+ ! calculate C(M), D(M)
+ ! (starting condition: C(MMAX) = 0 and D(MMAX) = 1)
! ***********************************************************************
- ! irregular
- do lm2 = 1, lmsize
- do lm1 = 1, lmsize
- dllp(lm1, lm2, npan) = czero
- cllp(lm1, lm2, npan) = czero
- end do
- end do
+ niter_beta = 3
+ cllp(:, :, npan) = czero
+ dllp(:, :, npan) = czero
do lm1 = 1, lmsize
dllp(lm1, lm1, npan) = cone
end do
do ipan = npan, 1, -1
- call zcopy(lmsize*lmsize, cllp(1,1,ipan), 1, cllp(1,1,ipan-1), 1)
- call zcopy(lmsize*lmsize, dllp(1,1,ipan), 1, dllp(1,1,ipan-1), 1)
+
+ cllp(:, :, ipan-1) = cllp(:, :, ipan)
+ dllp(:, :, ipan-1) = dllp(:, :, ipan)
+
call zgemm('n', 'n', lmsize, lmsize, lmsize, cone, mihvz(1,1,ipan), lmsize, cllp(1,1,ipan), lmsize, cone, cllp(1,1,ipan-1), lmsize)
call zgemm('n', 'n', lmsize, lmsize, lmsize, cone, mihvy(1,1,ipan), lmsize, dllp(1,1,ipan), lmsize, cone, cllp(1,1,ipan-1), lmsize)
call zgemm('n', 'n', lmsize, lmsize, lmsize, -cone, mijvz(1,1,ipan), lmsize, cllp(1,1,ipan), lmsize, cone, dllp(1,1,ipan-1), lmsize)
call zgemm('n', 'n', lmsize, lmsize, lmsize, -cone, mijvy(1,1,ipan), lmsize, dllp(1,1,ipan), lmsize, cone, dllp(1,1,ipan-1), lmsize)
+
end do
+ ! ***********************************************************************
+ ! invert beta = dllp(:, :, 0)
+ ! ***********************************************************************
+
+ betainv = dllp(:, :, 0)
+ call zgetrf(lmsize, lmsize, betainv, lmsize, ipiv, info)
+ call zgetri(lmsize, betainv, lmsize, ipiv, work, lmsize*lmsize, info)
+
+ if(use_cheby_quadprec) then
+ allocate (qcone, qczero)
+ qcone = (1.q0,0.0q0)
+ qczero = (0.q0,0.0q0)
+ allocate (qmihvy(lmsize,lmsize), qmihvz(lmsize,lmsize))
+ allocate (qmijvy(lmsize,lmsize), qmijvz(lmsize,lmsize))
+ allocate (qyif(lmsize2,lmsize,0:ncheb))
+ allocate (qbetainv(lmsize,lmsize), qbetainv_save(lmsize,lmsize))
+ allocate (qsll(lmsize2,lmsize))
+ allocate (qcllp(lmsize,lmsize,0:npan), qdllp(lmsize,lmsize,0:npan))
+ allocate (qcllptemp(lmsize,lmsize), qdllptemp(lmsize,lmsize))
+ qbetainv = betainv
+ end if
+
+ do iter_beta = 1, niter_beta
+
+ if(.not.use_cheby_quadprec) then
+
+ dllp(:, :, npan) = betainv
+ cllp(:, :, npan) = czero
+
+ do lm2 = 1, lmsize
+ dllp(lm2, lm2, npan) = betainv(lm2,lm2) - cone
+ end do
+
+ do ipan = npan, 1, -1
+
+ cllp(:, :, ipan-1) = cllp(:, :, ipan) + mihvy(:, :, ipan)
+ dllp(:, :, ipan-1) = dllp(:, :, ipan) - mijvy(:, :, ipan)
+
+ call zgemm('n', 'n', lmsize, lmsize, lmsize, cone, mihvz(1,1,ipan), lmsize, cllp(1,1,ipan), lmsize, cone, cllp(1,1,ipan-1), lmsize)
+ call zgemm('n', 'n', lmsize, lmsize, lmsize, cone, mihvy(1,1,ipan), lmsize, dllp(1,1,ipan), lmsize, cone, cllp(1,1,ipan-1), lmsize)
+ call zgemm('n', 'n', lmsize, lmsize, lmsize, -cone, mijvz(1,1,ipan), lmsize, cllp(1,1,ipan), lmsize, cone, dllp(1,1,ipan-1), lmsize)
+ call zgemm('n', 'n', lmsize, lmsize, lmsize, -cone, mijvy(1,1,ipan), lmsize, dllp(1,1,ipan), lmsize, cone, dllp(1,1,ipan-1), lmsize)
+
+ end do
+
+ betainv_save = betainv
+
+ call zgemm('n', 'n', lmsize, lmsize, lmsize, -cone, betainv_save, lmsize, dllp(1,1,0), lmsize, cone, betainv, lmsize)
+
+ dllpmax = 0.d0
+ do lm1 = 1,lmsize
+ do lm2 = 1,lmsize
+ dllpval = dllp(lm1,lm2,0)
+ if(lm1.ne.lm2.and.abs(dllpval).gt.dllpmax) dllpmax = abs(dllpval)
+ if(lm1.eq.lm2.and.abs(dllpval-cone).gt.dllpmax) dllpmax = abs(dllpval-cone)
+ end do
+ end do
+
+ else
+
+ qdllp(:, :, npan) = qbetainv
+ qcllp(:, :, npan) = qczero
+
+ do lm2 = 1, lmsize
+ qdllp(lm2, lm2, npan) = qbetainv(lm2,lm2) - qcone
+ end do
+
+ do ipan = npan, 1, -1
+
+ qmihvz(:, :) = mihvz(:, :, ipan)
+ qmihvy(:, :) = mihvy(:, :, ipan)
+ qmijvz(:, :) = mijvz(:, :, ipan)
+ qmijvy(:, :) = mijvy(:, :, ipan)
+
+ qcllp(:, :, ipan-1) = qcllp(:, :, ipan) + qmihvy(:, :)
+ qdllp(:, :, ipan-1) = qdllp(:, :, ipan) - qmijvy(:, :)
+
+ call cqgemm(lmsize,lmsize,lmsize,qcone,qmihvz,lmsize,qcllp(1,1,ipan),lmsize,qcone,qcllp(1,1,ipan-1),lmsize)
+ call cqgemm(lmsize,lmsize,lmsize,qcone,qmihvy,lmsize,qdllp(1,1,ipan),lmsize,qcone,qcllp(1,1,ipan-1),lmsize)
+ call cqgemm(lmsize,lmsize,lmsize,-qcone,qmijvz,lmsize,qcllp(1,1,ipan),lmsize,qcone,qdllp(1,1,ipan-1),lmsize)
+ call cqgemm(lmsize,lmsize,lmsize,-qcone,qmijvy,lmsize,qdllp(1,1,ipan),lmsize,qcone,qdllp(1,1,ipan-1),lmsize)
+
+ end do
+
+ qbetainv_save = qbetainv
+
+ call cqgemm(lmsize, lmsize, lmsize, -qcone, qbetainv_save, lmsize, qdllp(1,1,0), lmsize, qcone, qbetainv, lmsize)
+ dllpmax = 0.0d0
+ do lm1 = 1,lmsize
+ do lm2 = 1,lmsize
+ dllpval = qdllp(lm1,lm2,0)
+ if(abs(dllpval).gt.dllpmax) dllpmax = abs(dllpval)
+ end do
+ end do
+ end if
+
+ write(6,*) 'dllpmax',dllpmax,'iter_beta',iter_beta
+
+ end do ! niter_beta
+
! ***********************************************************************
! determine the irregular solution sll by using 5.14
! ***********************************************************************
+
+ if(.not.use_cheby_quadprec) then
+
+ do ipan = 0, npan
+ do lm1 = 1, lmsize
+ dllp(lm1,lm1,ipan) = dllp(lm1,lm1,ipan) + cone
+ end do
+ end do
+
+ do ipan = 1, npan
+ cllptemp(:, :) = cllp(:, :, ipan)
+ dllptemp(:, :) = dllp(:, :, ipan)
+ cllp(:, :,ipan) = cllptemp(:, :)*(cone+cone)
+ dllp(:, :,ipan) = dllptemp(:, :)*(cone+cone)
+ call zgemm('n', 'n', lmsize, lmsize, lmsize, -cone, cllptemp, lmsize, dllp(1,1,0), lmsize, cone, cllp(1,1,ipan), lmsize)
+ call zgemm('n', 'n', lmsize, lmsize, lmsize, -cone, dllptemp, lmsize, dllp(1,1,0), lmsize, cone, dllp(1,1,ipan), lmsize)
+ end do
+
do ipan = 1, npan
do icheb = 0, ncheb
mn = ipan*ncheb + ipan - icheb
@@ -209,7 +338,38 @@ contains
call zgemm('n', 'n', lmsize2, lmsize, lmsize, cone, yif(1,1,icheb,ipan), lmsize2, dllp(1,1,ipan), lmsize, cone, sll(1,1,mn), lmsize2)
end do
end do
+
+ else
+
+ do ipan = 0, npan
+ do lm1 = 1, lmsize
+ qdllp(lm1,lm1,ipan) = qdllp(lm1,lm1,ipan) + qcone
+ end do
+ end do
+ do ipan = 1, npan
+ qcllptemp(:, :) = qcllp(:, :, ipan)
+ qdllptemp(:, :) = qdllp(:, :, ipan)
+ qcllp(:, :,ipan) = qcllptemp(:, :)*(qcone + qcone)
+ qdllp(:, :,ipan) = qdllptemp(:, :)*(qcone + qcone)
+ call cqgemm(lmsize, lmsize, lmsize, -qcone, qcllptemp, lmsize, qdllp(1,1,0), lmsize, qcone, qcllp(1,1,ipan), lmsize)
+ call cqgemm(lmsize, lmsize, lmsize, -qcone, qdllptemp, lmsize, qdllp(1,1,0), lmsize, qcone, qdllp(1,1,ipan), lmsize)
+ end do
+
+ cllp = qcllp
+ do ipan = 1, npan
+ qyif(:,:,:) = yif(:,:,:,ipan)
+ do icheb = 0, ncheb
+ mn = ipan*ncheb + ipan - icheb
+ call zgemm('n', 'n', lmsize2, lmsize, lmsize, cone, zif(1,1,icheb,ipan), lmsize2, cllp(1,1,ipan), lmsize, czero, sll(1,1,mn), lmsize2)
+ call cqgemm(lmsize2, lmsize, lmsize, qcone, qyif(1,1,icheb), lmsize2, qdllp(1,1,ipan), lmsize, qczero, qsll, lmsize2)
+
+ sll(:,:,mn) = sll(:,:,mn) + qsll(:,:)
+
+ end do
+ end do
+ end if
+
if (idotime==1) call timing_stop('afterlocal')
if (idotime==1) call timing_start('endstuff')
@@ -218,8 +378,66 @@ contains
if (idotime==1) call timing_stop('checknan')
if (idotime==1) call timing_stop('sll-glob')
- deallocate (work, cllp, dllp, mihvy, mihvz, mijvy, mijvz, yif, zif, stat=ierror)
+ deallocate (work, betainv, betainv_save, cllp, dllp, cllptemp, dllptemp, mihvy, mihvz, mijvy, mijvz, yif, zif, stat=ierror)
+ if (ierror/=0) stop '[sll-glob] ERROR in deallocating arrays'
+ if(use_cheby_quadprec) deallocate (qbetainv, qbetainv_save, qcllp, qdllp, qcllptemp, qdllptemp, qmihvy, qmihvz, qmijvy, qmijvz, qyif, qsll, stat=ierror)
if (ierror/=0) stop '[sll-glob] ERROR in deallocating arrays'
end subroutine sll_global_solutions
end module mod_sll_global_solutions
+
+ SUBROUTINE CQGEMM (M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC )
+ IMPLICIT NONE
+! .. Scalar Arguments ..
+ INTEGER M, N, K, LDA, LDB, LDC
+ COMPLEX*32 ALPHA, BETA
+! .. Array Arguments ..
+ COMPLEX*32 A( LDA, * ), B( LDB, * ), C( LDC, * )
+! ..
+! .. Local Scalars ..
+ COMPLEX*32 TEMP
+ INTEGER I, J, L
+! .. Parameters ..
+ COMPLEX*32 ONE
+ PARAMETER ( ONE = ( 1.0Q+0, 0.0Q+0 ) )
+ COMPLEX*32 ZERO
+ PARAMETER ( ZERO = ( 0.0Q+0, 0.0Q+0 ) )
+! ..
+!
+ IF( ALPHA.EQ.ZERO )THEN
+ IF( BETA.EQ.ZERO )THEN
+ DO 20, J = 1, N
+ DO 10, I = 1, M
+ C( I, J ) = ZERO
+ 10 CONTINUE
+ 20 CONTINUE
+ ELSE
+ DO 40, J = 1, N
+ DO 30, I = 1, M
+ C( I, J ) = BETA*C( I, J )
+ 30 CONTINUE
+ 40 CONTINUE
+ END IF
+ RETURN
+ END IF
+ DO 90, J = 1, N
+ IF( BETA.EQ.ZERO )THEN
+ DO 50, I = 1, M
+ C( I, J ) = ZERO
+ 50 CONTINUE
+ ELSE IF( BETA.NE.ONE )THEN
+ DO 60, I = 1, M
+ C( I, J ) = BETA*C( I, J )
+ 60 CONTINUE
+ END IF
+ DO 80, L = 1, K
+ IF( B( L, J ).NE.ZERO )THEN
+ TEMP = ALPHA*B( L, J )
+ DO 70, I = 1, M
+ C( I, J ) = C( I, J ) + TEMP*A( I, L )
+ 70 CONTINUE
+ END IF
+ 80 CONTINUE
+ 90 CONTINUE
+ RETURN
+ END
diff --git a/source/common/runoptions.F90 b/source/common/runoptions.F90
index 064a1fb7a..e5a9f44d2 100644
--- a/source/common/runoptions.F90
+++ b/source/common/runoptions.F90
@@ -20,6 +20,7 @@ module mod_runoptions
logical :: calc_GF_Efermi = .false. !!calculation of cluster Green function at E Fermi (former: 'GF-EF')
logical :: set_cheby_nospeedup = .false. !!always calculate irregular solution in Chebychev solver (even if not needed) (former: 'norllsll')
logical :: set_cheby_nosoc = .false. !!decouple matrices in Chebychev solver neglecting SOC (former: 'NOSOC')
+ logical :: use_cheby_quadprec = .false. !!use quadruple precision in Chebychev solver for irregular solution (former: 'quadprec')
logical :: calc_complex_bandstructure = .false. !!complex band structure (former: 'COMPLEX')
logical :: calc_exchange_couplings = .false. !!calculate magnetic exchange coupling parameters (former: 'XCPL')
logical :: calc_exchange_couplings_energy = .false. !!write energy-resolved Jij-files also if npol/=0 (former: 'Jijenerg')
@@ -199,6 +200,7 @@ module mod_runoptions
call set_runoption(print_refpot , '<print_refpot>' , '<REFPOT>' )
call set_runoption(write_lloyd_file , '<write_lloyd_file>' , '<llyfile>')
call set_runoption(symmetrize_potential_madelung , '<symmetrize_potential_madelung>' , '<potsymm>' )
+ call set_runoption(use_cheby_quadprec , '<use_cheby_quadprec>' , '<quadprec>')
call set_runoption(set_cheby_nosoc , '<set_cheby_nosoc>' , '<NOSOC>' )
call set_runoption(set_empty_system , '<set_empty_system>' , '<zeropot>' )
call set_runoption(write_tmat_file , '<write_tmat_file>' , '<tmatfile>')
@@ -513,6 +515,9 @@ module mod_runoptions
else if (keyword == 'POTSYMM ') then
symmetrize_potential_madelung = .true.
write (1337, *) " Enable runoption 'symmetrize_potential_madelung'"
+ else if (keyword == 'QUADPREC') then
+ use_cheby_quadprec = .true.
+ write (1337, *) " Enable runoption 'use_cheby_quadprec'"
else if (keyword == 'NOSOC ') then
set_cheby_nosoc = .true.
write (1337, *) " Enable runoption 'set_cheby_nosoc'"
@@ -736,6 +741,7 @@ module mod_runoptions
call mpi_bcast(print_refpot , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(write_lloyd_file , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(symmetrize_potential_madelung , 1, mpi_logical, master, mpi_comm_world, ierr)
+ call mpi_bcast(use_cheby_quadprec , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(set_cheby_nosoc , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(set_empty_system , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(write_tmat_file , 1, mpi_logical, master, mpi_comm_world, ierr)
--
GitLab