From 9ab80078b061f748972999b1e7d436f66d2e73c4 Mon Sep 17 00:00:00 2001
From: Manuel dos Santos Dias <m.dos.santos.dias@fz-juelich.de>
Date: Mon, 6 Dec 2021 12:44:05 +0100
Subject: [PATCH] KKRimp: added option to switch back to old rllsll
---
cmake/source_list_KKRimp.txt | 4 +
source/KKRimp/calctmat_bauernew.f90 | 141 +++++++++++++++----------
source/KKRimp/energyloop.F90 | 4 +-
source/KKRimp/rhooutnew.f90 | 157 +++++++++++++++++++---------
source/KKRimp/rhoval_new.f90 | 9 +-
source/KKRimp/type_wavefunction.f90 | 3 +-
source/KKRimp/wavefunctodisc.f90 | 7 +-
source/common/runoptions.F90 | 4 +
8 files changed, 214 insertions(+), 115 deletions(-)
diff --git a/cmake/source_list_KKRimp.txt b/cmake/source_list_KKRimp.txt
index 0b77fbb70..d0f10c2a1 100644
--- a/cmake/source_list_KKRimp.txt
+++ b/cmake/source_list_KKRimp.txt
@@ -47,8 +47,12 @@ add_executable(
source/common/radial_solver_Chebychev/cheb.f90
source/common/radial_solver_Chebychev/cheb2oldgrid.f90
source/common/radial_solver_Chebychev/chebint.f90
+ source/common/radial_solver_Chebychev/rll_global_solutions.F90
+ source/common/radial_solver_Chebychev/rll_local_solutions.F90
source/common/radial_solver_Chebychev/rllsll.F90
source/common/radial_solver_Chebychev/rllsllsourceterms.f90
+ source/common/radial_solver_Chebychev/sll_global_solutions.F90
+ source/common/radial_solver_Chebychev/sll_local_solutions.F90
source/common/radial_solver_Chebychev/create_newmesh.f90
source/common/radial_solver_Chebychev/vllmat.f90
source/common/radial_solver_Chebychev/vllmatsra.f90
diff --git a/source/KKRimp/calctmat_bauernew.f90 b/source/KKRimp/calctmat_bauernew.f90
index b0c3d136d..c90985134 100644
--- a/source/KKRimp/calctmat_bauernew.f90
+++ b/source/KKRimp/calctmat_bauernew.f90
@@ -35,6 +35,8 @@ use mod_interpolpot, only: interpolpot
use mod_vllmat, only: vllmat
use mod_vllmatsra, only: vllmatsra
use mod_rllsll, only: rllsll
+use mod_rll_global_solutions, only: rll_global_solutions
+use mod_sll_global_solutions, only: sll_global_solutions ! MdSD: TEST
use mod_calccouplingconstants, only: calccouplingdeltat
use mod_config, only: config_testflag
use mod_spinorbit_ham, only: spinorbit_ham
@@ -399,12 +401,12 @@ end if
!#######################################################
if ( .not. config_testflag('nosph') .and. nsra/=5 ) then
- allocate(tmatsph(nspin*(lmaxatom+1)), stat=istat)
- tmatsph = czero
- if(istat/=0) stop 'Error allocating tmatsph in calctmat_bauernew'
- call calcsph(nsra, cellnew%nrmaxnew, cellnew%nrmaxnew, lmaxatom, nspin/(2-kspinorbit), zatom, eryd, &
- lmpot, lmsize, cellnew%rmeshnew, vins, cellnew%ncheb, cellnew%npan_tot, cellnew%rpan_intervall, jlk_index, &
- hlk, jlk, hlk2, jlk2, gmatprefactor, tmatsph, tmattemp, use_sratrick)
+ allocate(tmatsph(nspin*(lmaxatom+1)), stat=istat)
+ if(istat/=0) stop 'Error allocating tmatsph in calctmat_bauernew'
+ tmatsph = czero
+ call calcsph(nsra, cellnew%nrmaxnew, cellnew%nrmaxnew, lmaxatom, nspin/(2-kspinorbit), zatom, eryd, &
+ lmpot, lmsize, cellnew%rmeshnew, vins, cellnew%ncheb, cellnew%npan_tot, cellnew%rpan_intervall, jlk_index, &
+ hlk, jlk, hlk2, jlk2, gmatprefactor, tmatsph, tmattemp, use_sratrick, .true.)
end if
if ( config_testflag('writesourceterms')) then
@@ -445,12 +447,15 @@ end if
if (.not. allocated(wavefunction%SLL)) then
allocate (wavefunction%SLL(lmsize2,lmsize,cellnew%nrmaxnew,1),&
+ wavefunction%ULL(lmsize2,lmsize,cellnew%nrmaxnew,1),&
wavefunction%RLL(lmsize2,lmsize,cellnew%nrmaxnew,1))
- wavefunction%SLL = czero
- wavefunction%RLL = czero
+! MdSD: redundant
+! wavefunction%SLL = czero
+! wavefunction%RLL = czero
end if
wavefunction%rll=(0.0D0,0.0D0)
+wavefunction%ull=(0.0D0,0.0D0)
wavefunction%sll=(0.0D0,0.0D0)
! might be deleted in the future
@@ -472,18 +477,28 @@ if (nsra==4) then
jlk2= jlk2 / sqrt( (1.0D0,0.0D0)+(eryd)/cvlight**2)
end if
-
!#######################################################!
! calculate the right-hand side solution of the single-site wave functions
!#######################################################!
tmat%tmat = czero
-call rllsll(cellnew%rpan_intervall, cellnew%rmeshnew, Vpotll, wavefunction%RLL(:,:,:,1), wavefunction%SLL(:,:,:,1), &
- tmat%tmat, cellnew%ncheb, cellnew%npan_tot, lmsize, lmsize2, nsra*(1+kspinorbit)*(lmaxatom+1), cellnew%nrmaxnew, nsra, &
- jlk_index, hlk, jlk, hlk2, jlk2, GMATPREFACTOR, '1', '1', '0', use_sratrick, tmattemp)
-
+if ( config_testflag('use_rllsll') ) then
+ call rllsll(cellnew%rpan_intervall, cellnew%rmeshnew, Vpotll, wavefunction%RLL(:,:,:,1), wavefunction%SLL(:,:,:,1), &
+ tmat%tmat, cellnew%ncheb, cellnew%npan_tot, lmsize, lmsize2, nsra*(1+kspinorbit)*(lmaxatom+1), cellnew%nrmaxnew, nsra, &
+ jlk_index, hlk, jlk, hlk2, jlk2, GMATPREFACTOR, '1', '1', '0', use_sratrick, tmattemp)
+ ! MdSD: if using the old rllsll this is needed for rhooutnew
+ wavefunction%ULL(:,:,:,1)=wavefunction%RLL(:,:,:,1)
+else
+ call rll_global_solutions(cellnew%rpan_intervall, cellnew%rmeshnew, Vpotll, wavefunction%ULL(:,:,:,1), wavefunction%RLL(:,:,:,1), &
+ tmat%tmat, cellnew%ncheb, cellnew%npan_tot, lmsize, lmsize2, nsra*(1+kspinorbit)*(lmaxatom+1), cellnew%nrmaxnew, nsra, &
+ jlk_index, hlk, jlk, hlk2, jlk2, gmatprefactor, '1', use_sratrick, tmattemp)
+ call sll_global_solutions(cellnew%rpan_intervall, cellnew%rmeshnew, Vpotll, wavefunction%SLL(:,:,:,1), &
+ cellnew%ncheb, cellnew%npan_tot, lmsize, lmsize2, nsra*(1+kspinorbit)*(lmaxatom+1), cellnew%nrmaxnew, nsra, &
+ jlk_index, hlk, jlk, hlk2, jlk2, gmatprefactor, '1', use_sratrick)
+end if
if (nsra==2) then ! for nummerical reasons a factor of cvlight has been added to the equations
! which needs to be removed now
wavefunction%RLL(lmsize+1:,:,:,1)=wavefunction%RLL(lmsize+1:,:,:,1)/cvlight
+ wavefunction%ULL(lmsize+1:,:,:,1)=wavefunction%ULL(lmsize+1:,:,:,1)/cvlight
wavefunction%SLL(lmsize+1:,:,:,1)=wavefunction%SLL(lmsize+1:,:,:,1)/cvlight
end if
@@ -496,6 +511,8 @@ if (nsra==3) then
do ir=1,cellnew%nrmaxnew
call SINGLE_TRANSFORM(lmaxatom,(lmaxatom+1)**2,wavefunction%RLL(1:lmsize,1:lmsize,ir,1),'REL>RLM')
call SINGLE_TRANSFORM(lmaxatom,(lmaxatom+1)**2,wavefunction%RLL(lmsize+1:2*lmsize,1:lmsize,ir,1),'REL>RLM')
+ call SINGLE_TRANSFORM(lmaxatom,(lmaxatom+1)**2,wavefunction%ULL(1:lmsize,1:lmsize,ir,1),'REL>RLM')
+ call SINGLE_TRANSFORM(lmaxatom,(lmaxatom+1)**2,wavefunction%ULL(lmsize+1:2*lmsize,1:lmsize,ir,1),'REL>RLM')
call SINGLE_TRANSFORM(lmaxatom,(lmaxatom+1)**2,wavefunction%SLL(1:lmsize,1:lmsize,ir,1),'REL>RLM')
call SINGLE_TRANSFORM(lmaxatom,(lmaxatom+1)**2,wavefunction%SLL(lmsize+1:2*lmsize,1:lmsize,ir,1),'REL>RLM')
end do
@@ -505,6 +522,7 @@ end if
if ( config_testflag('conjgtest')) then
call conjugate4(wavefunction%RLL)
+ call conjugate4(wavefunction%ULL)
call conjugate4(wavefunction%SLL)
end if
@@ -512,10 +530,10 @@ call timing_stop('---rll call---')
if ( config_testflag('kappamutest')) then
call RLL_TRANSFORM(wavefunction%RLL(:,:,:,1),lmaxatom,'REL>RLM')
+ call RLL_TRANSFORM(wavefunction%ULL(:,:,:,1),lmaxatom,'REL>RLM')
call RLL_TRANSFORM(wavefunction%SLL(:,:,:,1),lmaxatom,'REL>RLM')
end if
-
!#######################################################
! In case the option 'nosph' is not set. The output t-matrix
! just contains the non-sph part of the t-matrix. Thus,
@@ -527,8 +545,6 @@ if ( .not. config_testflag('nosph') .or. nsra==5 ) then
end do
end if
-
-
if (config_testflag('write_tmat_all')) then
write (filename, '(A,I0.3,A,I0.3,A)') 'tmat_atom_', iatom, '_energ_', ie, '.dat'
open (888888, file=trim(filename), form='formatted')
@@ -536,6 +552,25 @@ if (config_testflag('write_tmat_all')) then
write (888888, '(2ES25.16)') tmat%tmat(:, :)
close (888888)
end if
+if ( config_testflag('tmatdebug') ) then
+ do lm1=1,lmsize
+ do lm2=1,lmsize
+ write(4000,'(50000E25.14)') wavefunction%RLL(lm2,lm1,:,1)
+ write(4001,'(50000E25.14)') wavefunction%SLL(lm2,lm1,:,1)
+ end do
+ end do
+end if
+
+if (wavefunction%nvec==2) then
+ if ( config_testflag('tmatdebug') ) then
+ do lm1=1,lmsize
+ do lm2=lmsize+1,2*lmsize
+ write(4010,'(50000E25.14)') wavefunction%SLL(lm2,lm1,:,1)
+ write(4011,'(50000E25.14)') wavefunction%RLL(lm2,lm1,:,1)
+ end do
+ end do
+ end if
+end if
!#######################################################
! If spin-orbit coupling is used the left solution of the
@@ -546,8 +581,9 @@ if ((kspinorbit==1).and.calcleft) then
if (.not. allocated(wavefunction%SLLleft)) then
allocate (wavefunction%SLLleft(lmsize2,lmsize,cellnew%nrmaxnew,1),&
wavefunction%RLLleft(lmsize2,lmsize,cellnew%nrmaxnew,1))
- wavefunction%SLLleft = czero
- wavefunction%RLLleft = czero
+ ! MdSD: redundant
+ ! wavefunction%SLLleft = czero
+ ! wavefunction%RLLleft = czero
end if
wavefunction%SLLleft=(0.0D0,0.0D0)
@@ -560,28 +596,46 @@ if ((kspinorbit==1).and.calcleft) then
hlk2 = czero
call rllsllsourceterms(nsra, wavefunction%nvec, eryd, cellnew%rmeshnew, cellnew%nrmaxnew, cellnew%nrmaxnew, &
lmaxatom, lmsize, use_fullgmat, jlk_index, hlk, jlk, hlk2, jlk2, GMATPREFACTOR)
-
+
+ ! MdSD: I put this here just in case it's needed
+ ! if (nsra==4) then
+ ! hlk = hlk / sqrt( (1.0D0,0.0D0)+(eryd)/cvlight**2)
+ ! jlk = jlk / sqrt( (1.0D0,0.0D0)+(eryd)/cvlight**2)
+ ! hlk2= hlk2 / sqrt( (1.0D0,0.0D0)+(eryd)/cvlight**2)
+ ! jlk2= jlk2 / sqrt( (1.0D0,0.0D0)+(eryd)/cvlight**2)
+ ! end if
if ( .not. config_testflag('nosph') .and. nsra/=5 ) then
- call calcsph(nsra, cellnew%nrmaxnew, cellnew%nrmaxnew, lmaxatom, nspin/(2-kspinorbit), zatom, eryd, &
- lmpot, lmsize, cellnew%rmeshnew, vins, cellnew%ncheb, cellnew%npan_tot, cellnew%rpan_intervall, jlk_index, &
- hlk2, jlk2, hlk, jlk, gmatprefactor, tmatsph, tmattemp, use_sratrick)
+ call calcsph(nsra, cellnew%nrmaxnew, cellnew%nrmaxnew, lmaxatom, nspin/(2-kspinorbit), zatom, eryd, &
+ lmpot, lmsize, cellnew%rmeshnew, vins, cellnew%ncheb, cellnew%npan_tot, cellnew%rpan_intervall, jlk_index, &
+ hlk2, jlk2, hlk, jlk, gmatprefactor, tmatsph, tmattemp, use_sratrick, .true.)
end if
- call rllsll(cellnew%rpan_intervall, cellnew%rmeshnew, Vpotll2, wavefunction%RLLleft(:,:,:,1), wavefunction%SLLleft(:,:,:,1), &
- tmattemp, cellnew%ncheb, cellnew%npan_tot, lmsize, lmsize2, nsra*(1+kspinorbit)*(lmaxatom+1), cellnew%nrmaxnew, nsra, &
- ! ------------> watch out here changed the order for left and right solution <-----------
- jlk_index, hlk2, jlk2, hlk, jlk, GMATPREFACTOR, '1', '1', '0', use_sratrick, tmattemp)
+ if ( config_testflag('use_rllsll') ) then
+ call rllsll(cellnew%rpan_intervall, cellnew%rmeshnew, Vpotll2, wavefunction%RLLleft(:,:,:,1), wavefunction%SLLleft(:,:,:,1), &
+ tmattemp, cellnew%ncheb, cellnew%npan_tot, lmsize, lmsize2, nsra*(1+kspinorbit)*(lmaxatom+1), cellnew%nrmaxnew, nsra, &
+ ! ------------> watch out here changed the order for left and right solution <-----------
+ jlk_index, hlk2, jlk2, hlk, jlk, GMATPREFACTOR, '1', '1', '0', use_sratrick, tmattemp)
+ else
+ ! MdSD: here SLLleft is being used to save memory, as ULLleft is never needed
+ call rll_global_solutions(cellnew%rpan_intervall, cellnew%rmeshnew, Vpotll2, wavefunction%SLLleft(:,:,:,1), wavefunction%RLLleft(:,:,:,1), &
+ tmattemp, cellnew%ncheb, cellnew%npan_tot, lmsize, lmsize2, nsra*(1+kspinorbit)*(lmaxatom+1), cellnew%nrmaxnew, nsra, &
+ jlk_index, hlk2, jlk2, hlk, jlk, gmatprefactor, '1', use_sratrick, tmattemp)
+ wavefunction%SLLleft=(0.0D0,0.0D0)
+ call sll_global_solutions(cellnew%rpan_intervall, cellnew%rmeshnew, Vpotll2, wavefunction%SLLleft(:,:,:,1), &
+ cellnew%ncheb, cellnew%npan_tot, lmsize, lmsize2, nsra*(1+kspinorbit)*(lmaxatom+1), cellnew%nrmaxnew, nsra, &
+ jlk_index, hlk2, jlk2, hlk, jlk, gmatprefactor, '1', use_sratrick)
+ end if
if (nsra==2) then
- wavefunction%RLLleft(lmsize+1:,:,:,1) = wavefunction%RLLleft(lmsize+1:,:,:,1)/(cvlight)
- wavefunction%SLLleft(lmsize+1:,:,:,1) = wavefunction%SLLleft(lmsize+1:,:,:,1)/(cvlight)
+ wavefunction%RLLleft(lmsize+1:,:,:,1) = wavefunction%RLLleft(lmsize+1:,:,:,1)/cvlight
+ wavefunction%SLLleft(lmsize+1:,:,:,1) = wavefunction%SLLleft(lmsize+1:,:,:,1)/cvlight
end if
if ( config_testflag('tmatdebug') ) then
do lm1=1,lmsize
do lm2=1,lmsize
- write(4100,'(50000E25.14)') wavefunction%rllleft(lm2,lm1,:,1)
- write(4101,'(50000E25.14)') wavefunction%sllleft(lm2,lm1,:,1)
+ write(4100,'(50000E25.14)') wavefunction%RLLleft(lm2,lm1,:,1)
+ write(4101,'(50000E25.14)') wavefunction%SLLleft(lm2,lm1,:,1)
end do
end do
end if
@@ -590,8 +644,8 @@ if ((kspinorbit==1).and.calcleft) then
if ( config_testflag('tmatdebug') ) then
do lm1=1,lmsize
do lm2=lmsize+1,2*lmsize
- write(4110,'(50000E25.14)') wavefunction%rllleft(lm2,lm1,:,1)
- write(4111,'(50000E25.14)') wavefunction%sllleft(lm2,lm1,:,1)
+ write(4110,'(50000E25.14)') wavefunction%RLLleft(lm2,lm1,:,1)
+ write(4111,'(50000E25.14)') wavefunction%SLLleft(lm2,lm1,:,1)
end do
end do
end if
@@ -599,25 +653,6 @@ if ((kspinorbit==1).and.calcleft) then
end if !(kspinorbit==1)
-if ( config_testflag('tmatdebug') ) then
- do lm1=1,lmsize
- do lm2=1,lmsize
- write(4000,'(50000E25.14)') wavefunction%rll(lm2,lm1,:,1)
- write(4001,'(50000E25.14)') wavefunction%sll(lm2,lm1,:,1)
- end do
- end do
-end if
-
-if (wavefunction%nvec==2) then
- if ( config_testflag('tmatdebug') ) then
- do lm1=1,lmsize
- do lm2=lmsize+1,2*lmsize
- write(4010,'(50000E25.14)') wavefunction%rll(lm2,lm1,:,1)
- write(4011,'(50000E25.14)') wavefunction%sll(lm2,lm1,:,1)
- end do
- end do
- end if
-end if
!#######################################################
! calculation of Jij's by a Lichtenstein-like approach
@@ -629,7 +664,7 @@ if (.not. allocated (tmat%deltaT_Jij)) then
end if
if (config%calcJijmat==1) then
-call calccouplingdeltat(wavefunction,tmat%deltaT_Jij,cellnew,gauntcoeff(lmaxatom),theta,phi,lmmax,lmsize,lmaxatom,lmpot,cellnew%nrmaxnew)
+ call calccouplingdeltat(wavefunction,tmat%deltaT_Jij,cellnew,gauntcoeff(lmaxatom),theta,phi,lmmax,lmsize,lmaxatom,lmpot,cellnew%nrmaxnew)
end if
!#######################################################
@@ -668,7 +703,7 @@ if ( .not. config_testflag('nosph') .and. nsra/=5 ) then
end if
-
end subroutine calctmat_bauernew
+
end module mod_calctmat_bauernew
diff --git a/source/KKRimp/energyloop.F90 b/source/KKRimp/energyloop.F90
index 554ed0cfd..d72b18a00 100644
--- a/source/KKRimp/energyloop.F90
+++ b/source/KKRimp/energyloop.F90
@@ -554,7 +554,7 @@ do ie=mpi_iebounds(1,my_rank), mpi_iebounds(2,my_rank)
if ( config_testflag('rlltodisc') ) then
call wavefunctodisc_write(wavefunction(iatom,ispin),cellnew(iatom),iatom,ispin,my_rank)
else
- deallocate(wavefunction(iatom,ispin)%rll, wavefunction(iatom,ispin)%sll)
+ deallocate(wavefunction(iatom,ispin)%rll, wavefunction(iatom,ispin)%ull, wavefunction(iatom,ispin)%sll)
if ((config%kspinorbit==1).and.calcleft) then
deallocate(wavefunction(iatom,ispin)%rllleft, wavefunction(iatom,ispin)%sllleft)
end if
@@ -695,7 +695,7 @@ do ie=mpi_iebounds(1,my_rank), mpi_iebounds(2,my_rank)
if ( wavefunction(iatom,ispin)%deallocate==1 ) then
- deallocate(wavefunction(iatom,ispin)%rll, wavefunction(iatom,ispin)%sll)
+ deallocate(wavefunction(iatom,ispin)%rll, wavefunction(iatom,ispin)%ull, wavefunction(iatom,ispin)%sll)
if (config%kspinorbit==1) then
deallocate(wavefunction(iatom,ispin)%rllleft, wavefunction(iatom,ispin)%sllleft)
end if
diff --git a/source/KKRimp/rhooutnew.f90 b/source/KKRimp/rhooutnew.f90
index 7c8735143..2ada7ef10 100644
--- a/source/KKRimp/rhooutnew.f90
+++ b/source/KKRimp/rhooutnew.f90
@@ -48,7 +48,7 @@ subroutine rhooutnew(gauntcoeff, df, gmatin, ek, cellnew, wavefunction, rho2nsc,
complex (kind=dp) :: cdenns(irmd,nspinden)
complex (kind=dp) :: gmat(lmsize,lmsize)
complex (kind=dp) :: gmatin(lmsize,lmsize)
- complex (kind=dp) :: qnsi(lmsize,lmsize), rlltemp(lmsize,lmsize)
+ complex (kind=dp) :: pnsi(lmsize,lmsize), qnsi(lmsize,lmsize)
complex (kind=dp) :: cdenlm(irmd,lmmaxatom,nspinden) ! lm-dos
complex (kind=dp) :: rho2nsc(irmd,lmpotd,nspinden)
complex (kind=dp) :: gflle_part(lmsize,lmsize) ! lda+u
@@ -122,75 +122,132 @@ subroutine rhooutnew(gauntcoeff, df, gmatin, ek, cellnew, wavefunction, rho2nsc,
! WATCH OUT CHECK IF A FACTOR OF M_0 needs to be added into the Greensfunction
!
! #########################################################3
- if (allocated(wavefunction%sllleft)) then
- qnsi(:,:) = wavefunction%sllleft(1:lmsize,1:lmsize,ir,1)
- else
- qnsi(:,:) = wavefunction%sll(1:lmsize,1:lmsize,ir,1)
- end if
-
- if (allocated(wavefunction%sllleft)) then
- rlltemp = wavefunction%rllleft(1:lmsize,1:lmsize,ir,1)
- else
- rlltemp = wavefunction%rll(1:lmsize,1:lmsize,ir,1)
- end if
+ ! if (allocated(wavefunction%sllleft)) then
+ ! qnsi(:,:) = wavefunction%sllleft(1:lmsize,1:lmsize,ir,1)
+ ! else
+ ! qnsi(:,:) = wavefunction%sll(1:lmsize,1:lmsize,ir,1)
+ ! end if
+
+ ! if (allocated(wavefunction%sllleft)) then
+ ! pnsi = wavefunction%rllleft(1:lmsize,1:lmsize,ir,1)
+ ! else
+ ! pnsi = wavefunction%rll(1:lmsize,1:lmsize,ir,1)
+ ! end if
! this is the prefactor for the gmatll*rllleft term in the first zgemm
! if the onsite densit is calculated alone we set this to zero
alpha = cone
if (config_testflag('calc_onsite_only')) alpha = czero
- ! changed the second mode to transpose - bauer
- call zgemm('n','t',lmsize,lmsize,lmsize, alpha,rlltemp, lmsize,gmat,lmsize,ek,qnsi,lmsize)
+ ! ! changed the second mode to transpose - bauer
+ ! call zgemm('n','t',lmsize,lmsize,lmsize, alpha,pnsi, lmsize,gmat,lmsize,ek,qnsi,lmsize)
- rlltemp(:,:) = wavefunction%rll(1:lmsize,1:lmsize,ir,1)
+ ! pnsi(:,:) = wavefunction%rll(1:lmsize,1:lmsize,ir,1)
- call zgemm('n','t',lmsize,lmsize,lmsize,cone,rlltemp, lmsize,qnsi,lmsize,czero,wr(1,1,ir),lmsize)
+ ! call zgemm('n','t',lmsize,lmsize,lmsize,cone,pnsi, lmsize,qnsi,lmsize,czero,wr(1,1,ir),lmsize)
+
+
+ ! ------------------------------------------------------------------
+ ! MdSD: now coded like in KKRhost
+ ! ------------------------------------------------------------------
+ ! S' irregular left wfn
+ if (allocated(wavefunction%sllleft)) then
+ qnsi(1:lmsize,1:lmsize) = wavefunction%sllleft(1:lmsize,1:lmsize,ir,1)
+ else
+ qnsi(:,:) = wavefunction%sll(1:lmsize,1:lmsize,ir,1)
+ end if
+ ! R regular right wfn (redefined to U)
+ pnsi(1:lmsize,1:lmsize) = wavefunction%ull(1:lmsize,1:lmsize,ir,1)
+ ! onsite contribution sqrt(E)*R*S'
+ call zgemm('N','T',lmsize,lmsize,lmsize,ek,pnsi,lmsize,qnsi,lmsize,czero,wr(:,:,ir),lmsize)
+ ! ------------------------------------------------------------------
+ ! R' regular left wfn
+ if (allocated(wavefunction%sllleft)) then
+ pnsi(1:lmsize,1:lmsize) = wavefunction%rllleft(1:lmsize,1:lmsize,ir,1)
+ else
+ pnsi(1:lmsize,1:lmsize) = wavefunction%rll(1:lmsize,1:lmsize,ir,1)
+ end if
+ ! MdSD: note that this transpose is followed by another transpose in the next zgemm
+ call zgemm('N','T',lmsize,lmsize,lmsize,alpha,pnsi,lmsize,gmat,lmsize,czero,qnsi,lmsize)
+ ! R regular right wfn (not redefined)
+ pnsi(1:lmsize,1:lmsize) = wavefunction%rll(1:lmsize,1:lmsize,ir,1)
+ ! backscattering contribution R*G*R' is added to onsite contribution
+ call zgemm('N','T',lmsize,lmsize,lmsize,cone,pnsi,lmsize,qnsi,lmsize,cone,wr(:,:,ir),lmsize)
+ ! ------------------------------------------------------------------
if (nsra.eq.2 .and. (.not. config_testflag('nosmallcomp')) ) then
- if (allocated(wavefunction%sllleft)) then
- qnsi(:,:) = -wavefunction%sllleft(lmsize+1:2*lmsize,1:lmsize,ir,1) ! attention to the
- ! additional minus sign
- ! ##########################################################################################
- ! Drittler assumes that for the left solution, is given by the right solution with an
- ! additional minus sign. This minus sign is contained inside the equations to calculate
- ! the electronic density. While calculating the left solution, the minus sign is already
- ! included in the left solution. To make calculations consistant a factor of -1 is included
- ! which cancels out by the routines of Drittler
- ! ##########################################################################################
+ ! if (allocated(wavefunction%sllleft)) then
+ ! qnsi(:,:) = -wavefunction%sllleft(lmsize+1:2*lmsize,1:lmsize,ir,1) ! attention to the
+ ! ! additional minus sign
+ ! ! ##########################################################################################
+ ! ! Drittler assumes that for the left solution, is given by the right solution with an
+ ! ! additional minus sign. This minus sign is contained inside the equations to calculate
+ ! ! the electronic density. While calculating the left solution, the minus sign is already
+ ! ! included in the left solution. To make calculations consistant a factor of -1 is included
+ ! ! which cancels out by the routines of Drittler
+ ! ! ##########################################################################################
+ ! else
+ ! qnsi(:,:) = wavefunction%sll(lmsize+1:2*lmsize,1:lmsize,ir,1)
+ ! end if
+ ! if (allocated(wavefunction%rllleft)) then
+ ! pnsi = -wavefunction%rllleft(lmsize+1:2*lmsize,1:lmsize,ir,1) ! attention to the
+ ! ! additional minus sign
+ ! ! ##########################################################################################
+ ! ! Drittler assumes that for the left solution, is given by the right solution with an
+ ! ! additional minus sign. This minus sign is contained inside the equations to calculate
+ ! ! the electronic density. While calculating the left solution, the minus sign is already
+ ! ! included in the left solution. To make calculations consistant a factor of -1 is included
+ ! ! which cancels out by the routines of Drittler
+ ! ! ##########################################################################################
+ ! else
+ ! pnsi = wavefunction%rll(lmsize+1:2*lmsize,1:lmsize,ir,1)
+ ! end if
+
+ ! ! changed the second mode to transpose - bauer
+ ! call zgemm('n','t',lmsize,lmsize,lmsize,gmat,pnsi, &
+ ! lmsize,gmat,lmsize,ek,qnsi,lmsize)
+
+ ! pnsi = wavefunction%rll(lmsize+1:2*lmsize,1:lmsize,ir,1)!/cvlight
+
+ ! call zgemm('n','t',lmsize,lmsize,lmsize,cone,pnsi, &
+ ! lmsize,qnsi,lmsize,cone,wr(1,1,ir),lmsize)
+
+
+ ! ------------------------------------------------------------------
+ ! MdSD: now coded like in KKRhost
+ ! ------------------------------------------------------------------
+ ! S' irregular left wfn
+ if (allocated(wavefunction%sllleft)) then ! this minus sign is explained above
+ qnsi(1:lmsize,1:lmsize) = -wavefunction%sllleft(lmsize+1:2*lmsize,1:lmsize,ir,1)
else
- qnsi(:,:) = wavefunction%sll(lmsize+1:2*lmsize,1:lmsize,ir,1)
+ qnsi(1:lmsize,1:lmsize) = wavefunction%sll(lmsize+1:2*lmsize,1:lmsize,ir,1)
end if
- if (allocated(wavefunction%rllleft)) then
- rlltemp = -wavefunction%rllleft(lmsize+1:2*lmsize,1:lmsize,ir,1) ! attention to the
- ! additional minus sign
- ! ##########################################################################################
- ! Drittler assumes that for the left solution, is given by the right solution with an
- ! additional minus sign. This minus sign is contained inside the equations to calculate
- ! the electronic density. While calculating the left solution, the minus sign is already
- ! included in the left solution. To make calculations consistant a factor of -1 is included
- ! which cancels out by the routines of Drittler
- ! ##########################################################################################
+ ! R regular right wfn (redefined to U)
+ pnsi(1:lmsize,1:lmsize) = wavefunction%ull(lmsize+1:2*lmsize,1:lmsize,ir,1)
+ ! onsite contribution sqrt(E)*R*S'
+ call zgemm('N','T',lmsize,lmsize,lmsize,ek,pnsi,lmsize,qnsi,lmsize,cone,wr(:,:,ir),lmsize)
+ ! ------------------------------------------------------------------
+ ! R' regular left wfn
+ if (allocated(wavefunction%sllleft)) then ! this minus sign is explained above
+ pnsi(1:lmsize,1:lmsize) = -wavefunction%rllleft(lmsize+1:2*lmsize,1:lmsize,ir,1)
else
- rlltemp = wavefunction%rll(lmsize+1:2*lmsize,1:lmsize,ir,1)
+ pnsi(1:lmsize,1:lmsize) = wavefunction%rll(lmsize+1:2*lmsize,1:lmsize,ir,1)
end if
-
- ! changed the second mode to transpose - bauer
- call zgemm('n','t',lmsize,lmsize,lmsize,gmat,rlltemp, &
- lmsize,gmat,lmsize,ek,qnsi,lmsize)
-
- rlltemp = wavefunction%rll(lmsize+1:2*lmsize,1:lmsize,ir,1)!/cvlight
-
- call zgemm('n','t',lmsize,lmsize,lmsize,cone,rlltemp, &
- lmsize,qnsi,lmsize,cone,wr(1,1,ir),lmsize)
-
+ ! MdSD: note that this transpose is followed by another transpose in the next zgemm
+ call zgemm('N','T',lmsize,lmsize,lmsize,alpha,pnsi,lmsize,gmat,lmsize,czero,qnsi,lmsize)
+ ! R regular right wfn (not redefined)
+ pnsi(1:lmsize,1:lmsize) = wavefunction%rll(lmsize+1:2*lmsize,1:lmsize, ir,1)
+ ! backscattering contribution R*G*R' is added to onsite contribution
+ call zgemm('N','T',lmsize,lmsize,lmsize,cone,pnsi,lmsize,qnsi,lmsize,cone,wr(:,:,ir),lmsize)
+ ! ------------------------------------------------------------------
end if
if (corbital/=0) then
call zgemm('n','n',lmsize,lmsize,lmsize,cone,loperator(:,:,corbital), &
- lmsize,wr(:,:,ir),lmsize,czero,rlltemp,lmsize)
- wr(:,:,ir) = rlltemp
+ lmsize,wr(:,:,ir),lmsize,czero,pnsi,lmsize)
+ wr(1:lmsize,1:lmsize,ir) = pnsi(1:lmsize,1:lmsize)
end if
end do !ir
diff --git a/source/KKRimp/rhoval_new.f90 b/source/KKRimp/rhoval_new.f90
index 5669dd313..d389ec484 100644
--- a/source/KKRimp/rhoval_new.f90
+++ b/source/KKRimp/rhoval_new.f90
@@ -120,14 +120,7 @@ if ( .not. allocated(wavefunction%rll)) then
stop '[rhoval_new] rll not allocated'
end if
-
-
-
-
- imt1=cellnew%ipan_intervall(cellnew%npan_log+cellnew%npan_eq)+1
-
-
-
+imt1=cellnew%ipan_intervall(cellnew%npan_log+cellnew%npan_eq)+1
call rhooutnew(gauntcoeff,df,gmatll,ek,cellnew,wavefunction,rho2ns_complex(:,:,:), &
config%nsra, &
diff --git a/source/KKRimp/type_wavefunction.f90 b/source/KKRimp/type_wavefunction.f90
index 9fa4614b4..29b73e0df 100644
--- a/source/KKRimp/type_wavefunction.f90
+++ b/source/KKRimp/type_wavefunction.f90
@@ -11,7 +11,8 @@ module type_wavefunction
integer :: lmsize,lmsize2,nrmaxnew
integer :: nvec
- double complex,allocatable :: sll(:,:,:,:), rll(:,:,:,:)
+ ! MdSD: sll never seems to be used, better skip their allocation
+ double complex,allocatable :: sll(:,:,:,:), rll(:,:,:,:), ull(:,:,:,:)
double complex,allocatable :: sllleft(:,:,:,:), rllleft(:,:,:,:)
integer :: deallocate = 0
diff --git a/source/KKRimp/wavefunctodisc.f90 b/source/KKRimp/wavefunctodisc.f90
index 08d7a0b51..8f25475df 100644
--- a/source/KKRimp/wavefunctodisc.f90
+++ b/source/KKRimp/wavefunctodisc.f90
@@ -45,14 +45,17 @@ contains
recl=wlength*4*wavefunction%lmsize2*wavefunction%lmsize*cellnew%nrmaxnew,form='unformatted')
open(unit=2343366,file='temp_sllleft_'//trim(ctemp)//'.txt',access='direct',&
recl=wlength*4*wavefunction%lmsize2*wavefunction%lmsize*cellnew%nrmaxnew,form='unformatted')
+ open(unit=2343367,file='temp_ull_' //trim(ctemp)//'.txt',access='direct',&
+ recl=wlength*4*wavefunction%lmsize2*wavefunction%lmsize*cellnew%nrmaxnew,form='unformatted')
end if
call setlabel(iatom,ispin,ilabel)
if ( allocated(wavefunction%rll) ) then
write(2343363,rec=ilabel) wavefunction%rll
+ write(2343367,rec=ilabel) wavefunction%ull
wavefunction%rll_saved=1
- deallocate(wavefunction%rll)
+ deallocate(wavefunction%rll,wavefunction%ull)
end if
if ( allocated(wavefunction%sll) ) then
@@ -100,7 +103,9 @@ contains
if ( .not. allocated(wavefunction%rll) .and. wavefunction%rll_saved==1) then
allocate(wavefunction%rll(wavefunction%lmsize2,wavefunction%lmsize,cellnew%nrmaxnew,1))
+ allocate(wavefunction%ull(wavefunction%lmsize2,wavefunction%lmsize,cellnew%nrmaxnew,1))
read(2343363,rec=ilabel) wavefunction%rll
+ read(2343367,rec=ilabel) wavefunction%ull
end if
if ( .not. allocated(wavefunction%sll) .and. wavefunction%sll_saved==1) then
allocate(wavefunction%sll(wavefunction%lmsize2,wavefunction%lmsize,cellnew%nrmaxnew,1))
diff --git a/source/common/runoptions.F90 b/source/common/runoptions.F90
index 23cddb18c..bb874a85c 100644
--- a/source/common/runoptions.F90
+++ b/source/common/runoptions.F90
@@ -66,6 +66,7 @@ module mod_runoptions
logical :: torque_operator_onlySph = .false. !! for torque operator: include only the spherically symmetric part (former: 'ONLYSPH')
logical :: use_BdG = .false. !! use Bogoliubov-de-Gennes Formalism (former: 'useBdG')
logical :: use_Chebychev_solver = .false. !! use the Chebychev solver (former: 'NEWSOSOL')
+ logical :: use_rllsll = .false. !! switch to previous approach to compute wavefunctions in Chebyshev solver
logical :: use_cond_LB = .false. !! perform calculation of conductance in Landauer-Büttiker formalism (former: 'CONDUCT')
logical :: use_cont = .false. !! no usage of embedding points. NEMB is set to 0. (former: 'CONT')
logical :: use_deci_onebulk= .false. !! in case of decimation: use same bulk on right and left. Speeds up calculations. (former: 'ONEBULK')
@@ -155,6 +156,7 @@ module mod_runoptions
call set_runoption(write_gmat_ascii , '<write_gmat_ascii>' , '<gmatasci>')
call set_runoption(use_rigid_Efermi , '<use_rigid_Efermi>' , '<rigid-ef>')
call set_runoption(use_Chebychev_solver , '<use_Chebychev_solver>' , '<NEWSOSOL>')
+ call set_runoption(use_rllsll , '<use_rllsll>')
call set_runoption(write_pkkr_input , '<write_pkkr_input>' , '<FERMIOUT>')
call set_runoption(calc_complex_bandstructure , '<calc_complex_bandstructure>' , '<COMPLEX>' )
call set_runoption(write_pkkr_operators , '<write_pkkr_operators>' , '<OPERATOR>')
@@ -702,6 +704,7 @@ module mod_runoptions
call mpi_bcast(write_gmat_ascii , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(use_rigid_Efermi , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(use_Chebychev_solver , 1, mpi_logical, master, mpi_comm_world, ierr)
+ call mpi_bcast(use_rllsll , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(write_pkkr_input , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(calc_complex_bandstructure , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(write_pkkr_operators , 1, mpi_logical, master, mpi_comm_world, ierr)
@@ -861,6 +864,7 @@ module mod_runoptions
write(iounit, '(A35,1x,1L,3x,A)') '<torque_operator_onlySph>=', torque_operator_onlySph, "for torque operator: include only the spherically symmetric part (former: 'ONLYSPH')"
write(iounit, '(A35,1x,1L,3x,A)') '<use_BdG>=', use_BdG, "use Bogoliubov-de-Gennes Formalism (former: 'useBdG')"
write(iounit, '(A35,1x,1L,3x,A)') '<use_Chebychev_solver>=', use_Chebychev_solver, "use the Chebychev solver (former: 'NEWSOSOL')"
+ write(iounit, '(A35,1x,1L,3x,A)') '<use_rllsll>=', use_rllsll, "switch to previous approach to compute wavefunctions in Chebyshev solver"
write(iounit, '(A35,1x,1L,3x,A)') '<use_cond_LB>=', use_cond_LB, "perform calculation of conductance in Landauer-Büttiker formalism (former: 'CONDUCT')"
write(iounit, '(A35,1x,1L,3x,A)') '<use_cont>=', use_cont, "no usage of embedding points. NEMB is set to 0. (former: 'CONT')"
write(iounit, '(A35,1x,1L,3x,A)') '<use_deci_onebulk>=', use_deci_onebulk, "in case of decimation: use same bulk on right and left. Speeds up calculations. (former: 'ONEBULK')"
--
GitLab