From c97279c1456e96736d7689b25c0bab8db41f4efb Mon Sep 17 00:00:00 2001
From: Philipp Ruessmann <p.ruessmann@fz-juelich.de>
Date: Thu, 2 Dec 2021 09:30:43 +0100
Subject: [PATCH] Activate use_gmat_unity and calc_onsite_only options in
 KKRhost and calc_onsite_only testflag for KKRimp

---
 source/KKRhost/main1b.F90    | 20 ++++++++++++++------
 source/KKRhost/rhooutnew.f90 | 14 ++++++++++----
 source/KKRimp/rhooutnew.f90  | 11 ++++++++---
 source/common/runoptions.F90 |  6 ++++++
 4 files changed, 38 insertions(+), 13 deletions(-)

diff --git a/source/KKRhost/main1b.F90 b/source/KKRhost/main1b.F90
index f7dda9d80..96d1c33f3 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 7bd9eee28..e5fbe1153 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 c56048b14..7c8735143 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 e4e73da64..23cddb18c 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
 
-- 
GitLab