diff --git a/source/KKRhost/main1b.F90 b/source/KKRhost/main1b.F90
index f7dda9d8074ad2279d2ea3b31b679368387dc543..96d1c33f3a4addfa16e1430b69e4d59644d23f69 100644
--- a/source/KKRhost/main1b.F90
+++ b/source/KKRhost/main1b.F90
@@ -37,7 +37,7 @@ contains
use :: mod_datatypes, only: dp, sp
use :: mod_runoptions, only: calc_exchange_couplings, formatted_file, set_gmat_to_zero, use_Chebychev_solver, &
use_qdos, use_readcpa, write_deci_tmat, write_gmat_plain, write_green_host, write_green_imp, write_kkrimp_input, &
- write_pkkr_input, write_pkkr_operators, write_rhoq_input, write_gmat_ascii, decouple_spin_cheby, write_double_precision
+ write_pkkr_input, write_pkkr_operators, write_rhoq_input, write_gmat_ascii, decouple_spin_cheby, write_double_precision, use_gmat_unity
use mod_constants, only: czero, cone, pi
use mod_operators_for_fscode, only: operators_for_fscode
use mod_getscratch, only: opendafile
@@ -662,11 +662,19 @@ contains
#ifdef CPP_TIMING
call timing_start('main1b - kloopz')
#endif
- call kloopz1_qdos(eryd, gmatll, ins, alat, ie, igf, nshell, naez, nofks(nmesh), volbz(nmesh), bzkp(1,1,nmesh), volcub(1,nmesh), cls, nacls, naclsmax, ncls, rr, rbasis, &
- ezoa, atom, rcls, icc, ginp, ideci, lefttinvll(1,1,1,1,ie), righttinvll(1,1,1,1,ie), vacflag, nlbasis, nrbasis, factl, natomimp, nsymat, dsymll, ratom, rrot, nsh1, &
- nsh2, ijtabsym, ijtabsh, icheck, invmod, refpot, trefll, tsst, msst, cfctor, cfctorinv, crel, rc, rrel, srrel, irrel, nrrel, drotq, symunitary, kmrot, natyp, ncpa, &
- icpa, itcpamax, cpatol, noq, iqat, itoq, conc, iprint, icpaflag, ispin, nspindd, tqdos, iqdosrun, & ! qdos
- dtrefll, dtmatll, dginp, lly_grtr(ie,ispin), tracet(ie,ispin), lly) ! LLY Lloyd
+ if (.not. use_gmat_unity) then
+ call kloopz1_qdos(eryd, gmatll, ins, alat, ie, igf, nshell, naez, nofks(nmesh), volbz(nmesh), bzkp(1,1,nmesh), volcub(1,nmesh), cls, nacls, naclsmax, ncls, rr, rbasis, &
+ ezoa, atom, rcls, icc, ginp, ideci, lefttinvll(1,1,1,1,ie), righttinvll(1,1,1,1,ie), vacflag, nlbasis, nrbasis, factl, natomimp, nsymat, dsymll, ratom, rrot, nsh1, &
+ nsh2, ijtabsym, ijtabsh, icheck, invmod, refpot, trefll, tsst, msst, cfctor, cfctorinv, crel, rc, rrel, srrel, irrel, nrrel, drotq, symunitary, kmrot, natyp, ncpa, &
+ icpa, itcpamax, cpatol, noq, iqat, itoq, conc, iprint, icpaflag, ispin, nspindd, tqdos, iqdosrun, & ! qdos
+ dtrefll, dtmatll, dginp, lly_grtr(ie,ispin), tracet(ie,ispin), lly) ! LLY Lloyd
+ else
+ ! set gmatll artificially to the unity matrix
+ if (myrank==master .and. ie_num==1 .and. ispin==1) write(*,'(/A/)') ' !!! setting GMATLL to unity matrix !!!'
+ do lm1=1, lmmaxd
+ gmatll(lm1, lm1, :) = cone
+ end do
+ end if !.not. use_gmat_unity
#ifdef CPP_TIMING
call timing_pause('main1b - kloopz')
diff --git a/source/KKRhost/rhooutnew.f90 b/source/KKRhost/rhooutnew.f90
index 7bd9eee288c31de05c7808154a01391ff3e1dbbf..e5fbe1153a4001ff3a1250782cef7cae721a52b5 100644
--- a/source/KKRhost/rhooutnew.f90
+++ b/source/KKRhost/rhooutnew.f90
@@ -22,7 +22,7 @@ contains
rll, sll, ull, rllleft, sllleft, ullleft, 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, decouple_spin_cheby
+ use :: mod_runoptions, only: calc_gmat_lm_full, use_ldau, decouple_spin_cheby, calc_onsite_only
use :: mod_profiling, only: memocc
use :: global_variables, only: lmmaxd, ncleb, ntotd, nfund, korbit
use :: mod_datatypes, only: dp
@@ -73,7 +73,7 @@ contains
integer :: ir, jspin, lm1, lm2, lm3, m1, l1, j, ifun
integer :: i_stat, i_all
real (kind=dp) :: c0ll
- complex (kind=dp) :: cltdf
+ complex (kind=dp) :: cltdf, alpha
integer, dimension (4) :: lmshift1
integer, dimension (4) :: lmshift2
complex (kind=dp), dimension (lmmaxd, lmmaxd, 3) :: loperator
@@ -130,6 +130,12 @@ contains
! big component of Dirac spinor
do ir = 1, irmdnew
+
+ ! 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 (calc_onsite_only) alpha = czero
+
do lm1 = 1, lmmaxd
do lm2 = 1, lmmaxd
qnsi(lm1, lm2) = sllleft(lm1, lm2, ir)
@@ -144,7 +150,7 @@ contains
pnsi(lm1, lm2) = rllleft(lm1, lm2, ir)
end do
end do
- call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, cone, pnsi, lmmaxd, gmatll, lmmaxd, czero, qnsi, lmmaxd)
+ call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, alpha, pnsi, lmmaxd, gmatll, lmmaxd, czero, qnsi, lmmaxd)
do lm1 = 1, lmmaxd
do lm2 = 1, lmmaxd
pnsi(lm1, lm2) = rll(lm1, lm2, ir)
@@ -172,7 +178,7 @@ contains
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)
+ call zgemm('N', 'T', lmmaxd, lmmaxd, lmmaxd, alpha, 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/KKRimp/rhooutnew.f90 b/source/KKRimp/rhooutnew.f90
index c56048b1407171ca74800ad1841d184249b4c545..7c87351439fbd3dbe2b4b3b68f61a93d17289116 100644
--- a/source/KKRimp/rhooutnew.f90
+++ b/source/KKRimp/rhooutnew.f90
@@ -54,7 +54,7 @@ subroutine rhooutnew(gauntcoeff, df, gmatin, ek, cellnew, wavefunction, rho2nsc,
complex (kind=dp) :: gflle_part(lmsize,lmsize) ! lda+u
! ..
! .. Local Scalars ..
- complex (kind=dp) :: cltdf
+ complex (kind=dp) :: cltdf, alpha
real (kind=dp) :: c0ll
integer :: ifun, ir, j, l1, lm1, lm2, lm3, m1
integer :: ispin, jspin, ierr
@@ -134,8 +134,13 @@ subroutine rhooutnew(gauntcoeff, df, gmatin, ek, cellnew, wavefunction, rho2nsc,
rlltemp = 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,cone,rlltemp, lmsize,gmat,lmsize,ek,qnsi,lmsize)
+ call zgemm('n','t',lmsize,lmsize,lmsize, alpha,rlltemp, lmsize,gmat,lmsize,ek,qnsi,lmsize)
rlltemp(:,:) = wavefunction%rll(1:lmsize,1:lmsize,ir,1)
@@ -172,7 +177,7 @@ subroutine rhooutnew(gauntcoeff, df, gmatin, ek, cellnew, wavefunction, rho2nsc,
end if
! changed the second mode to transpose - bauer
- call zgemm('n','t',lmsize,lmsize,lmsize,cone,rlltemp, &
+ 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
diff --git a/source/common/runoptions.F90 b/source/common/runoptions.F90
index e4e73da644e21930d0ac32cc2defc1c5d9a4640a..23cddb18c1b9ada2d3d439c1b63c798f63ecfd41 100644
--- a/source/common/runoptions.F90
+++ b/source/common/runoptions.F90
@@ -247,6 +247,8 @@ module mod_runoptions
call set_runoption(write_angles_alliter , '<write_angles_alliter>')
call set_runoption(write_tmat_all , '<write_tmat_all>')
call set_runoption(write_double_precision , '<write_double_precision>')
+ call set_runoption(calc_onsite_only , '<calc_onsite_only>')
+ call set_runoption(use_gmat_unity , '<use_gmat_unity>')
end subroutine read_runoptions
@@ -792,6 +794,8 @@ module mod_runoptions
call mpi_bcast(write_angles_alliter , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(write_tmat_all , 1, mpi_logical, master, mpi_comm_world, ierr)
call mpi_bcast(write_double_precision , 1, mpi_logical, master, mpi_comm_world, ierr)
+ call mpi_bcast(calc_onsite_only , 1, mpi_logical, master, mpi_comm_world, ierr)
+ call mpi_bcast(use_gmat_unity , 1, mpi_logical, master, mpi_comm_world, ierr)
end subroutine bcast_runoptions
#endif
@@ -909,6 +913,8 @@ module mod_runoptions
write(iounit, '(A35,1x,1L,3x,A)') '<write_angles_alliter>=', write_angles_alliter, "write out noncollinear angles for all iterations"
write(iounit, '(A35,1x,1L,3x,A)') '<write_tmat_all>=', write_tmat_all, "write out the tmat for all energies and all atoms"
write(iounit, '(A35,1x,1L,3x,A)') '<write_double_precision>=', write_double_precision, "write out kkrflex files in double precision"
+ write(iounit, '(A35,1x,1L,3x,A)') '<calc_onsite_only>=', calc_onsite_only, "calculate not the full Green function for the density but take the onsite part alone"
+ write(iounit, '(A35,1x,1L,3x,A)') '<use_gmat_unity>=', use_gmat_unity, "set the structural GF to the unity matrix for test purposes"
end subroutine print_runoptions