Fix LDA+U

- fix initialization of gfint
- implement QBOUND_LDAU to control how well LDAU potential should be converged

source/KKRimp/calcwldau.f

@@ -15,7 +15,7 @@
 !>           
 !-------------------------------------------------------------------------------
       SUBROUTINE CALCWLDAU(
-     >     NSPIN,NATOM,LMAXD,IRMD,LMAXATOM,DENSITY,STRMIX,
+     >     NSPIN,NATOM,LMAXD,IRMD,LMAXATOM,DENSITY,STRMIX, QBOUND_LDAU, 
      X     LDAU)
 
 C **********************************************************************
@@ -72,6 +72,7 @@ C Dummy arguments
 C
       INTEGER NATOM,NSPIN
       INTEGER LMAXATOM(:)
+      REAL*8 QBOUND_LDAU
       TYPE(LDAU_TYPE),ALLOCATABLE :: LDAU(:) ! lda+u variables dimension: (NATOM)
       TYPE(DENSITY_TYPE) :: DENSITY(:)
 
@@ -289,6 +290,10 @@ C Calculate rms error for wldau
 
 C Apply damping to the interaction matrix WLDAU (mixing):
       XMIX =STRMIX
+      if (QBOUND_LDAU > DSQRT(RMSTOT)) then
+        write(*,*) 'LDAU mixing: reached QBOUND_LDAU',QBOUND_LDAU
+        xmix = 0.0d0
+      end if
       write(*,*) 'Mixing old-new wldau with xmix=',xmix
       CALL WMIX(XMIX,LDAU(IAT)%WLDAU(:,:,:),WLDAU_OLD(:,:,:),MMAX,NSPIN)
       

source/KKRimp/config.f90

@@ -162,6 +162,12 @@ subroutine config_read(config)
                write(*,*) 'error in config ', keyword1
                stop
             end if
+         case ('QBOUND_LDAU')
+            read(string1,*,iostat=ios2) keyword1, config%qbound_ldau
+            if (ios2 /= 0) then
+               write(*,*) 'error in config ', keyword1
+               stop
+            end if
          case ('IMIX')
             read(string1,*,iostat=ios2) keyword1, config%IMIX
             if (ios2 /= 0) then

source/KKRimp/energyloop.F90

@@ -351,6 +351,8 @@ do iatom=1,natom
   density(iatom)%ncharge=0.0D0
   density(iatom)%rho2ns_integrated=(0.0D0,0.0D0)
   density(iatom)%rho2ns_integrated_scattering=(0.0D0,0.0D0)
+  density(iatom)%gfint=(0.0D0,0.0D0)
+  density(iatom)%gflle=(0.0D0,0.0D0)
 
   density(iatom)%orbitalmom=(0.0D0,0.0D0)
   density(iatom)%orbitalmom_lm=(0.0D0,0.0D0)

source/KKRimp/kkrflex.F90

@@ -355,7 +355,6 @@ if ( config_runflag('LLYsimple') ) then
     open(192837, file='kkrflex_llyfac', form='formatted', iostat=ierror)
     if(ierror/=0) stop 'Error: File kkrflex_llyfac not found, needed for LLYsimple option'
     read(192837, *) llyfac
-    !llyfac = 0.993769945E+00
     close(192837)
     write(*,*) 'Renormalize weights with factor:',llyfac
   end if
@@ -858,7 +857,7 @@ end if
       if (mod(itscf,config%ITDBRY).eq.0) mixldau = config%mixfac                           ! lda+u
     endif                                                                                  ! lda+u
     if ( config_testflag('freezeldau').or.config_runflag('freezeldau') ) mixldau = 0.d0    ! lda+u
-    call calcwldau(nspin,natom,lmaxd,cell(1)%nrmaxd,lmaxatom,density,mixldau,ldau)         ! lda+u
+    call calcwldau(nspin,natom,lmaxd,cell(1)%nrmaxd,lmaxatom,density,mixldau,config%qbound_ldau,ldau)         ! lda+u
 !-------------------------------------------------------------------                       ! lda+u
 
 

source/KKRimp/type_config.f90

@@ -76,6 +76,7 @@ integer                        ::  spinmixbound=99999
 real(kind=8)                   ::  mixfac   = 0.1      		!! mixing factor 
 real(kind=8)                   ::  fcm      = 2.0         	!! 
 real(kind=8)                   ::  qbound   = 1d-8     		!! 
+real(kind=8)                   ::  qbound_ldau   = 1d-4    	!! qbound for mixing of ldau potential (see calcwldau)
 integer                        ::  itdbry   = 40          	!! number of iterations which are used for the broyden mixing
 ! -------------------------
 ! lattice relaxation