Skip to content
Snippets Groups Projects
Commit 9ab80078 authored by Manuel Dias's avatar Manuel Dias
Browse files

KKRimp: added option to switch back to old rllsll

parent e22e4c3d
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
......@@ -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
......@@ -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
......
......@@ -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
......
......@@ -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, &
......
......@@ -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
......
......@@ -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))
......
......@@ -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')"
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment