hubbard1_setup.F90 22.8 KB
Newer Older
1 2 3 4 5 6
MODULE m_hubbard1_setup

   USE m_juDFT

   IMPLICIT NONE

7
   LOGICAL, PARAMETER :: l_setupdebug = .TRUE.  !Enable/Disable Debug outputs like dependency of occupation on chemical potential shift 
8 9 10
   CHARACTER(len=30), PARAMETER :: main_folder = "Hubbard1"


11

12 13
   CONTAINS

14
   SUBROUTINE hubbard1_setup(atoms,hub1,sym,mpi,noco,input,usdus,den,pot,gdft,results)
15 16 17

      USE m_types
      USE m_constants
18
      USE m_uj2f
19 20
      USE m_umtx
      USE m_vmmp
21
      USE m_vmmp21
22
      USE m_nmat_rot
23 24
      USE m_mudc
      USE m_denmat_dist
25 26 27
      USE m_gfcalc
      USE m_hubbard1_io
      USE m_add_selfen
28 29 30
#ifdef CPP_EDSOLVER
      USE EDsolver, only: EDsolver_from_cfg
#endif
Henning Janssen's avatar
Henning Janssen committed
31 32 33
#ifdef CPP_MPI
      INCLUDE "mpif.h"
#endif
34

35

36
      TYPE(t_atoms),    INTENT(IN)     :: atoms
37
      TYPE(t_hub1ham),  INTENT(INOUT)  :: hub1
38 39
      TYPE(t_sym),      INTENT(IN)     :: sym
      TYPE(t_mpi),      INTENT(IN)     :: mpi
40
      TYPE(t_noco),     INTENT(IN)     :: noco
41
      TYPE(t_input),    INTENT(IN)     :: input
Henning Janssen's avatar
Henning Janssen committed
42
      TYPE(t_usdus),    INTENT(IN)     :: usdus
43
      TYPE(t_potden),   INTENT(INOUT)  :: den
44
      TYPE(t_potden),   INTENT(INOUT)  :: pot
45
      TYPE(t_results),  INTENT(INOUT)  :: results
46
      TYPE(t_greensf),  INTENT(IN)     :: gdft !green's function calculated from the Kohn-Sham system
47
      
Henning Janssen's avatar
Henning Janssen committed
48 49 50
#ifdef CPP_MPI
      EXTERNAL MPI_BCAST
#endif
51

52
      INTEGER i_hia,nType,l,n_occ,ispin,m,iz,k,j,i_exc,i
53
      INTEGER io_error,ierr
54
      INTEGER indStart,indEnd
55
      REAL    mu_dc,e_lda_hia,exc,e_off
56

Henning Janssen's avatar
Henning Janssen committed
57
      CHARACTER(len=300) :: cwd,path,folder,xPath
58
      CHARACTER(len=8)   :: l_type*2,l_form*9
59

60
      TYPE(t_greensf)    :: gu 
61

62 63
      REAL     f0(atoms%n_hia,input%jspins),f2(atoms%n_hia,input%jspins)
      REAL     f4(atoms%n_hia,input%jspins),f6(atoms%n_hia,input%jspins)
64
      REAL     u(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,&
65
               -lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_hia)
66
      REAL     zero(atoms%n_hia)
67

68
      COMPLEX  mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_hia,3)
69
      COMPLEX  n_mmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_hia,input%jspins)
70 71
      COMPLEX  selfen(atoms%n_hia,2*(2*lmaxU_const+1),2*(2*lmaxU_const+1),2*gdft%nz)
      COMPLEX  e(2*gdft%nz)
72
      REAL     n_l(atoms%n_hia,input%jspins)
73
      LOGICAL  l_selfenexist,l_exist,l_linkedsolver,l_ccfexist,l_bathexist
74

75 76 77 78
      !To avoid confusing structure later on
#ifdef CPP_EDSOLVER
      l_linkedsolver = .TRUE.
#else
Henning Janssen's avatar
Henning Janssen committed
79
      l_linkedsolver = .FALSE.
80 81
#endif

82
      e = 0.0
Henning Janssen's avatar
Henning Janssen committed
83
      e_lda_hia = 0.0
84

85 86 87
      !Positions of the DFT+HIA elements in all DFT+U related arrays
      indStart = atoms%n_u+1
      indEnd   = atoms%n_u+atoms%n_hia 
88

89

90 91
      !TODO: We don't need to calculate the green's function in every iteration 
      !(what is an appropriate cutoff for the distance under which we calculate the greens function)
92

93
      IF(ANY(den%mmpMat(:,:,indStart:indEnd,:).NE.0.0).OR.hub1%l_runthisiter) THEN
94
         !Get the slater integrals from the U and J parameters
95
         CALL uj2f(input%jspins,atoms%lda_u(indStart:indEnd),atoms%n_hia,f0,f2,f4,f6)
96 97 98 99
         f0(:,1) = (f0(:,1) + f0(:,input%jspins) ) / 2
         f2(:,1) = (f2(:,1) + f2(:,input%jspins) ) / 2
         f4(:,1) = (f4(:,1) + f4(:,input%jspins) ) / 2
         f6(:,1) = (f6(:,1) + f6(:,input%jspins) ) / 2
Henning Janssen's avatar
Henning Janssen committed
100
         CALL umtx(atoms%lda_u(indStart:indEnd),atoms%n_hia,f0(:,1),f2(:,1),f4(:,1),f6(:,1),u)
101

102
         ! check for possible rotation of n_mmp
Henning Janssen's avatar
Fixes  
Henning Janssen committed
103
         n_mmp = den%mmpMat(:,:,indStart:indEnd,1:input%jspins)
104 105 106 107
         zero = 0.0
         CALL nmat_rot(atoms%lda_u(indStart:indEnd)%phi,atoms%lda_u(indStart:indEnd)%theta,zero,3,atoms%n_hia,input%jspins,atoms%lda_u(indStart:indEnd)%l,n_mmp)

         CALL v_mmp(sym,atoms,atoms%lda_u(indStart:indEnd),atoms%n_hia,input%jspins,input%l_dftspinpol,n_mmp,&
108
         u,f0,f2,pot%mmpMat(:,:,indStart:indEnd,:),e_lda_hia)
109

110 111 112 113 114 115
         IF(noco%l_mperp) THEN
             IF(ANY(atoms%lda_u(1:atoms%n_u)%phi.NE.0.0).OR.ANY(atoms%lda_u(1:atoms%n_u)%theta.NE.0.0)) CALL juDFT_error("vmmp21+Rot not implemented", calledby="u_setup")
             CALL v_mmp_21(atoms%lda_u(indStart:indEnd),atoms%n_u,den%mmpMat(:,:,indStart:indEnd,3),u,f0,f2,pot%mmpMat(:,:,indStart:indEnd,3),e_off)
             e_lda_hia = e_lda_hia + e_off
         ENDIF

116
         IF(hub1%l_runthisiter.AND.(ANY(gdft%gmmpMat(:,:,:,:,:,:).NE.0.0)).AND.mpi%irank.EQ.0) THEN 
117 118 119
            !The onsite green's function was calculated but the solver 
            !was not yet run
            !--> write out the configuration for the hubbard 1 solver 
120
            CALL gu%init(input,lmaxU_const,atoms,noco,nz_in=gdft%nz, e_in=gdft%e,de_in=gdft%de,matsub_in=gdft%nmatsub)
121 122 123

            !Get the working directory
            CALL get_environment_variable('PWD',cwd)
124
            path = TRIM(ADJUSTL(cwd)) // "/" // TRIM(ADJUSTL(main_folder))
125
            CALL SYSTEM('mkdir -p ' // TRIM(ADJUSTL(path)))
126 127
            !Remove everything from the last iteration (Good Idea??)
            CALL SYSTEM('rm -rf ' // TRIM(ADJUSTL(path)) // "/*")
128

129
            DO i_hia = 1, atoms%n_hia
130 131
               nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
               l = atoms%lda_u(atoms%n_u+i_hia)%l
132

Henning Janssen's avatar
Henning Janssen committed
133 134 135
               CALL hubbard1_path(atoms,i_hia,folder)
               WRITE(xPath,*) TRIM(ADJUSTL(cwd)) // "/" // TRIM(ADJUSTL(folder)) 
               CALL SYSTEM('mkdir -p ' // TRIM(ADJUSTL(xPath)))
136 137

               !For the first iteration we can fix the occupation and magnetic moments in the inp.xml file
138
               IF(hub1%iter.EQ.1.AND.ALL(den%mmpMat(:,:,indStart:indEnd,:).EQ.0.0)) THEN
139 140 141
                  n_l(i_hia,:) = hub1%init_occ(i_hia)/input%jspins
                  DO i_exc = 1, hub1%n_exc_given(i_hia)
                     hub1%mag_mom(i_hia,i_exc) = hub1%init_mom(i_hia,i_exc)
Henning Janssen's avatar
Henning Janssen committed
142
                  ENDDO
143 144 145 146 147 148
               ELSE
                  !calculate the occupation of the correlated shell
                  CALL occmtx(gdft,l,nType,atoms,sym,input,mmpMat(:,:,i_hia,:))
                  n_l(i_hia,:) = 0.0
                  DO ispin = 1, input%jspins
                     DO m = -l, l
149
                        n_l(i_hia,ispin) = n_l(i_hia,ispin) + REAL(mmpMat(m,m,i_hia,ispin))
150
                     ENDDO
151
                  ENDDO
152
               ENDIF
153

Henning Janssen's avatar
Henning Janssen committed
154
               !Initial Information (We are already on irank 0)
155
               WRITE(6,*)
156
               WRITE(6,9010) nType
157
               WRITE(6,"(A)") "Everything related to the solver (e.g. mu_dc) is given in eV"
158
               WRITE(6,*)
159 160 161
               WRITE(6,"(A)") "Occupation from DFT-Green's function:"
               WRITE(6,9020) 'spin-up','spin-dn'
               WRITE(6,9030) n_l(i_hia,:)
162

163
               !--------------------------------------------------------------------------
164
               ! Calculate the chemical potential for the solver 
165 166 167
               ! This is equal to the double-counting correction used in DFT+U
               !--------------------------------------------------------------------------
               ! V_FLL = U (n - 1/2) - J (n - 1) / 2
168
               ! V_AMF = U n/2 + 2l/[2(2l+1)] (U-J) n
169 170
               !--------------------------------------------------------------------------
               CALL mudc(atoms%lda_u(atoms%n_u+i_hia),n_l(i_hia,:),mu_dc,input%jspins)
171

172
               !Check wether the hubbard 1 solver was run:(old version)
173 174
               INQUIRE(file=TRIM(ADJUSTL(xPath)) // "se.atom",exist=l_selfenexist)

175
               IF(mpi%irank.EQ.0.AND..NOT.l_selfenexist) THEN
176 177
                  !Nearest Integer occupation
                  n_occ = ANINT(SUM(n_l(i_hia,:)))
178 179 180 181 182 183 184 185 186 187 188 189 190 191
                  IF(l_linkedsolver) THEN
                     !-------------------------------------------------------
                     ! Check for additional input files
                     !-------------------------------------------------------
                     !Is a crystal field matrix present in the work directory
                     INQUIRE(file=TRIM(ADJUSTL(cwd)) // "/" // TRIM(ADJUSTL(cfg_file_ccf)),exist=l_ccfexist)
                     IF(l_ccfexist) CALL read_ccfmat(TRIM(ADJUSTL(cwd)),hub1%ccfmat(i_hia,-l:l,-l:l),l)
                     !Is a bath parameter file present 
                     INQUIRE(file=TRIM(ADJUSTL(cwd)) // "/" // TRIM(ADJUSTL(cfg_file_bath)),exist=l_bathexist)
                     !Copy the bath file to the Hubbard 1 solver if its present
                     IF(l_bathexist) CALL SYSTEM('cp ' // TRIM(ADJUSTL(cfg_file_bath)) // ' ' // TRIM(ADJUSTL(xPath)))
                     !-------------------------------------------------------
                     ! Write the main config files
                     !-------------------------------------------------------
Henning Janssen's avatar
Fixes  
Henning Janssen committed
192 193
                     CALL hubbard1_input(xPath,i_hia,l,f0(i_hia,1),f2(i_hia,1),f4(i_hia,1),f6(i_hia,1),hub1,mu_dc,n_occ,&
                                        l_bathexist,hub1%iter==1.AND.ALL(pot%mmpmat(:,:,indStart:indEnd,:).EQ.0.0),.true.)
194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
                  ELSE
                     !If no Solver is linked we assume that the old solver is used and we write out some additional files
                     !Crystal field matrix (old version)
                     IF(hub1%ccf(i_hia).NE.0.0) THEN
                        CALL write_ccfmat(xPath,hub1%ccfmat(i_hia,-l:l,-l:l),l)
                     ENDIF
                     !Energy contour (old version)
                     IF(gdft%mode.NE.3) THEN
                        OPEN(unit=1337,file=TRIM(ADJUSTL(xPath)) // "contour.dat",status="replace",action="write")
                        DO iz = 1, gdft%nz
                           WRITE(1337,"(2f14.8)") REAL(gdft%e(iz))*hartree_to_ev_const, AIMAG(gdft%e(iz))*hartree_to_ev_const
                        ENDDO
                        CLOSE(unit= 1337)
                     ENDIF 
                     !-------------------------------------------------------
                     ! Write the main config files (old version)
                     !-------------------------------------------------------  
Henning Janssen's avatar
Fixes  
Henning Janssen committed
211
                     CALL hubbard1_input(xPath,i_hia,l,f0(i_hia,1),f2(i_hia,1),f4(i_hia,1),f6(i_hia,1),hub1,mu_dc,n_occ,.false.,.false.,.false.)
212
                  ENDIF
213 214
               ENDIF
            ENDDO
215

216 217 218 219
            !------------------------------------------------------------
            ! This loop runs the solver if it is available
            ! If not the program terminates here
            !------------------------------------------------------------
220
            DO i_hia = 1, atoms%n_hia 
221 222
               nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
               l = atoms%lda_u(atoms%n_u+i_hia)%l
Henning Janssen's avatar
Henning Janssen committed
223 224 225

               CALL hubbard1_path(atoms,i_hia,folder)
               WRITE(xPath,*) TRIM(ADJUSTL(cwd)) // "/" // TRIM(ADJUSTL(folder)) 
226
               INQUIRE(file=TRIM(ADJUSTL(xPath)) // "se.atom",exist=l_selfenexist)
Henning Janssen's avatar
Henning Janssen committed
227 228
               IF(mpi%irank.EQ.0) THEN
                  IF(.NOT.l_selfenexist) THEN
229
                     IF(.NOT.l_linkedsolver) CALL juDFT_END("Hubbard1 input has been written into Hubbard1/ (No Solver linked)",mpi%irank)
Henning Janssen's avatar
Henning Janssen committed
230
                     CALL timestart("Hubbard 1: EDsolver")
231
                     !We have to change into the Hubbard1 directory so that the solver routines can read the config
232 233 234
                     CALL CHDIR(TRIM(ADJUSTL(xPath)))
                     !Set up the energy points (z and z*)
                     e(1:gdft%nz) = gdft%e(1:gdft%nz)*hartree_to_ev_const  
Henning Janssen's avatar
Henning Janssen committed
235 236
                     e(gdft%nz+1:2*gdft%nz) = conjg(gdft%e(1:gdft%nz))*hartree_to_ev_const 
#ifdef CPP_EDSOLVER                 
237
                     CALL EDsolver_from_cfg(2*(2*l+1),2*gdft%nz,e,selfen(i_hia,:,:,1:2*gdft%nz),1)
Henning Janssen's avatar
Henning Janssen committed
238
#endif
Henning Janssen's avatar
Henning Janssen committed
239 240 241 242 243 244 245 246
                     !The solver is given everything in eV by default, so we need to convert back to htr
                     selfen(i_hia,:,:,:) = selfen(i_hia,:,:,:)/hartree_to_ev_const
                     CALL CHDIR(TRIM(ADJUSTL(cwd)))
                     CALL timestop("Hubbard 1: EDsolver")
                  ELSE               
                     !If there is no linked solver library we read in the selfenergy here
                     CALL read_selfen(xPath,selfen(i_hia,1:2*(2*l+1),1:2*(2*l+1),1:gdft%nz),gdft%nz,2*(2*l+1),.false.)
                  ENDIF
247 248 249 250
               ENDIF
            ENDDO

            IF(l_selfenexist.OR.l_linkedsolver) THEN
Henning Janssen's avatar
Henning Janssen committed
251 252 253
               !-------------------------------------------
               ! Postprocess selfenergy
               !-------------------------------------------
254 255 256
               DO i_hia = 1, atoms%n_hia
                  nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
                  l = atoms%lda_u(atoms%n_u+i_hia)%l
257
                  DO iz = 1, 2*gdft%nz
Henning Janssen's avatar
Henning Janssen committed
258 259 260
                     !---------------------------------------------
                     ! The order of spins is reversed in the Solver
                     !---------------------------------------------
261
                     CALL swapSpin(selfen(i_hia,:,:,iz),2*l+1)
Henning Janssen's avatar
Henning Janssen committed
262 263 264 265
                     ! The DFT green's function also includes the previous DFT+U correction
                     ! This is removed by substracting it from the selfenergy
                     !-------------------------------------------------------------
                     CALL removeU(selfen(i_hia,:,:,iz),l,input%jspins,pot%mmpMat(:,:,atoms%n_u+i_hia,:))
266 267
                  ENDDO
               ENDDO
268 269 270 271 272 273 274 275 276
               !----------------------------------------------------------------------
               ! Solution of the Dyson Equation
               !----------------------------------------------------------------------
               ! G(z)^(-1) = G_0(z)^(-1) - mu - Sigma(z)
               !----------------------------------------------------------------------
               ! Sigma(z) is the self-energy from the impurity solver 
               ! We introduce an additional chemical potential mu, which is determined
               ! so that the occupation of the correlated orbital does not change
               !----------------------------------------------------------------------
277
               CALL timestart("Hubbard 1: Add Selfenenergy")
Henning Janssen's avatar
Henning Janssen committed
278
               CALL add_selfen(gdft,gu,selfen,atoms,noco,hub1,sym,input,results%ef,n_l,mu_dc/hartree_to_ev_const,mmpMat)
279
               CALL timestop("Hubbard 1: Add Selfenenergy")
Henning Janssen's avatar
Henning Janssen committed
280 281 282 283
               IF(l_setupdebug) THEN
                  DO i_hia = 1, atoms%n_hia
                     nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
                     l = atoms%lda_u(atoms%n_u+i_hia)%l
284 285
                     CALL gfDOS(gdft,l,nType,800+i_hia+hub1%iter,atoms,input,results%ef)
                     CALL gfDOS(gu,l,nType,900+i_hia+hub1%iter,atoms,input,results%ef)
Henning Janssen's avatar
Henning Janssen committed
286
                     CALL writeSelfenElement(selfen(i_hia,:,:,1:gdft%nz),gdft%e,results%ef,gdft%nz,2*l+1)
Henning Janssen's avatar
Henning Janssen committed
287 288
                  ENDDO
               ENDIF
289 290 291
               !----------------------------------------------------------------------
               ! Calculate the distance and update the density matrix 
               !----------------------------------------------------------------------
Henning Janssen's avatar
Fixes  
Henning Janssen committed
292
               CALL n_mmp_dist(den%mmpMat(:,:,indStart:indEnd,:),mmpMat,atoms%n_hia,results,input)
293
               den%mmpMat(:,:,indStart:indEnd,:) = mmpMat(:,:,:,1:MERGE(3,input%jspins,noco%l_mperp)) !For now LDA+U in FLEUR ignores spin offdiagonal elements
294 295 296 297 298 299 300 301
               !----------------------------------------------------------------------
               ! Calculate DFT+U potential correction
               !----------------------------------------------------------------------
               ! check for possible rotation of n_mmp
               n_mmp = den%mmpMat(:,:,indStart:indEnd,:)
               zero = 0.0
               CALL nmat_rot(atoms%lda_u(indStart:indEnd)%phi,atoms%lda_u(indStart:indEnd)%theta,zero,3,atoms%n_hia,input%jspins,atoms%lda_u(indStart:indEnd)%l,n_mmp)
               CALL v_mmp(sym,atoms,atoms%lda_u(indStart:indEnd),atoms%n_hia,input%jspins,input%l_dftspinpol,n_mmp,&
Henning Janssen's avatar
Henning Janssen committed
302
                     u,f0,f2,pot%mmpMat(:,:,indStart:indEnd,:), e_lda_hia)
303 304 305 306 307 308

               IF(noco%l_mperp) THEN
                  IF(ANY(atoms%lda_u(1:atoms%n_u)%phi.NE.0.0).OR.ANY(atoms%lda_u(1:atoms%n_u)%theta.NE.0.0)) CALL juDFT_error("vmmp21+Rot not implemented", calledby="u_setup")
                  CALL v_mmp_21(atoms%lda_u(indStart:indEnd),atoms%n_u,den%mmpMat(:,:,indStart:indEnd,3),u,f0,f2,pot%mmpMat(:,:,indStart:indEnd,3),e_off)
                  e_lda_hia = e_lda_hia + e_off
               ENDIF
Henning Janssen's avatar
Henning Janssen committed
309
               results%e_ldau = MERGE(results%e_ldau,0.0,atoms%n_u>0) + e_lda_hia 
310
            ENDIF
311 312
         ELSE 
            !The solver does not need to be run so we just add the current energy correction from LDA+HIA 
313
            results%e_ldau = MERGE(results%e_ldau,0.0,atoms%n_u>0) + e_lda_hia 
314 315 316 317 318
         ENDIF
         !Write out the density matrix and potential matrix (compare u_setup.f90)
         IF (mpi%irank.EQ.0) THEN
            DO ispin = 1,input%jspins
               WRITE (6,'(a7,i3)') 'spin #',ispin
319
               DO i_hia = indStart, indEnd
320 321 322 323 324
                  nType = atoms%lda_u(i_hia)%atomType
                  l = atoms%lda_u(i_hia)%l
                  WRITE (l_type,'(i2)') 2*(2*l+1)
                  l_form = '('//l_type//'f12.7)'
                  WRITE (6,'(a20,i3)') 'n-matrix for atom # ',nType
325
                  WRITE (6,l_form) ((n_mmp(k,j,i_hia,ispin),k=-l,l),j=-l,l)
326 327 328 329 330 331 332
                  WRITE (6,'(a20,i3)') 'V-matrix for atom # ',nType
                  IF (atoms%lda_u(i_hia)%l_amf) THEN
                     WRITE (6,*) 'using the around-mean-field limit '
                  ELSE
                     WRITE (6,*) 'using the atomic limit of LDA+U '
                  ENDIF
                  WRITE (6,l_form) ((pot%mmpMat(k,j,i_hia,ispin),k=-l,l),j=-l,l)
333
               END DO
334 335
            END DO
            WRITE (6,*) results%e_ldau
336
         ENDIF
337
      ELSE IF(mpi%irank.NE.0) THEN
338
         results%e_ldau = MERGE(results%e_ldau,0.0,atoms%n_u>0)
339
         pot%mmpMat(:,:,indStart:indEnd,:) = CMPLX(0.0,0.0)
340 341
         !If we are on a different mpi%irank and no solver is linked we need to call juDFT_end here if the solver was not run
         !kind of a weird workaround (replace with something better)
342 343 344 345 346 347 348 349 350 351
         IF(.NOT.l_linkedsolver) THEN
            CALL get_environment_variable('PWD',cwd)
            DO i_hia = atoms%n_u+1, atoms%n_u+atoms%n_hia
               CALL hubbard1_path(atoms,i_hia,folder)
               WRITE(xPath,*) TRIM(ADJUSTL(cwd)) // "/" // TRIM(ADJUSTL(folder)) 
               INQUIRE(file=TRIM(ADJUSTL(xPath)) // "se.atom",exist=l_selfenexist)
               IF(.NOT.l_selfenexist) CALL juDFT_END("Hubbard1 input has been written into Hubbard1/ (No Solver linked)",mpi%irank)
            ENDDO
            !If we are here the solver was run and we go to MPI_BCAST
         ENDIF
352 353 354
      ELSE
         !occupation matrix is zero and LDA+Hubbard 1 shouldn't be run yet
         !There is nothing to be done yet just set the potential correction to 0
355
         WRITE(*,*) "No density matrix and GF found -> skipping LDA+HIA"
356
         pot%mmpMat(:,:,indStart:indEnd,:) = CMPLX(0.0,0.0)
357
         results%e_ldau = MERGE(results%e_ldau,0.0,atoms%n_u>0)
358
      ENDIF
359
#ifdef CPP_MPI
360 361
      !Broadcast both the potential and the density matrix here
      CALL MPI_BCAST(pot%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,indStart:indEnd,:),&
Henning Janssen's avatar
Henning Janssen committed
362
                     49*atoms%n_hia*MERGE(3,input%jspins,input%l_gfmperp),MPI_DOUBLE_COMPLEX,0,mpi%mpi_comm,ierr)
363
      CALL MPI_BCAST(den%mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,indStart:indEnd,:),&
Henning Janssen's avatar
Henning Janssen committed
364
                     49*atoms%n_hia*MERGE(3,input%jspins,input%l_gfmperp),MPI_DOUBLE_COMPLEX,0,mpi%mpi_comm,ierr)
Henning Janssen's avatar
Henning Janssen committed
365
      CALL MPI_BCAST(results%e_ldau,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
366
#endif 
367
      
368 369 370 371 372
      !FORMAT Statements:
9010  FORMAT("Setup for Hubbard 1 solver for atom ", I3, ": ")
9020  FORMAT(TR8,A7,TR3,A7)
9030  FORMAT(TR7,f8.4,TR2,f8.4)

373 374
   END SUBROUTINE hubbard1_setup

375 376 377 378 379 380 381 382
   SUBROUTINE swapSpin(mat,ns)

      IMPLICIT NONE

      COMPLEX,       INTENT(INOUT) :: mat(2*ns,2*ns)
      INTEGER,       INTENT(IN)    :: ns

      COMPLEX tmp(ns,ns)
Henning Janssen's avatar
Henning Janssen committed
383
      INTEGER i,j
384 385 386 387 388 389 390

      !Spin-diagonal
      tmp = mat(1:ns,1:ns)
      mat(1:ns,1:ns) = mat(ns+1:2*ns,ns+1:2*ns)
      mat(ns+1:2*ns,ns+1:2*ns) = tmp
      !Spin-offdiagonal
      tmp = mat(ns+1:2*ns,1:ns)
Henning Janssen's avatar
Henning Janssen committed
391 392 393 394 395 396 397 398 399 400 401
      DO i = 1, ns
         DO j = 1, ns
            mat(ns+j,i) = tmp(ns-i+1,ns-j+1)
         ENDDO
      ENDDO
      tmp = mat(1:ns,ns+1:2*ns)
      DO i = 1, ns
         DO j = 1, ns
            mat(j,ns+i) = tmp(ns-i+1,ns-j+1)
         ENDDO
      ENDDO
Henning Janssen's avatar
Henning Janssen committed
402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430
   END SUBROUTINE swapSpin


   SUBROUTINE removeU(mat,l,jspins,vmmp)

      USE m_constants

      IMPLICIT NONE
      COMPLEX,       INTENT(INOUT) :: mat(2*(2*l+1),2*(2*l+1))
      INTEGER,       INTENT(IN)    :: l
      INTEGER,       INTENT(IN)    :: jspins
      COMPLEX,       INTENT(IN)     :: vmmp(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,jspins)

      INTEGER ns,i,j,m,mp,ispin

      ns = 2*l+1

      DO i = 1, ns 
         DO j = 1, ns 
            m = i-1-l
            mp = j-1-l
            DO ispin = 1, jspins
               mat(i+(ispin-1)*ns,j+(ispin-1)*ns) = mat(i+(ispin-1)*ns,j+(ispin-1)*ns) - vmmp(m,mp,ispin)
               IF(jspins.EQ.1) mat(i+ns,j+ns) = mat(i+ns,j+ns) - vmmp(-m,-mp,ispin)
            ENDDO
         ENDDO
      ENDDO

   END SUBROUTINE removeU
431 432


Henning Janssen's avatar
Henning Janssen committed
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457
   SUBROUTINE hubbard1_path(atoms,i_hia,xPath)

      !Defines the folder structure
      ! The Solver is run in the subdirectories
      ! Hubbard1/ if only one Hubbard1 prodcedure is run 
      ! Hubbard1/atom_label_l if there are more

      USE m_types

      IMPLICIT NONE 

      TYPE(t_atoms),       INTENT(IN)  :: atoms 
      INTEGER,             INTENT(IN)  :: i_hia
      CHARACTER(len=300),  INTENT(OUT) :: xPath

      CHARACTER(len=300) :: folder,fmt
      CHARACTER(len=1)   :: l_name(0:3)
      INTEGER nType,l

      l_name(0:3) = (/"s","p","d","f"/)

      nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
      l = atoms%lda_u(atoms%n_u+i_hia)%l
      xPath = TRIM(ADJUSTL(main_folder))
      IF(atoms%n_hia>1) THEN
458 459
         WRITE(fmt,'("(A",I2.2,",A1,A1,A1)")') LEN(TRIM(ADJUSTL(atoms%label(nType))))
         WRITE(folder,fmt) TRIM(ADJUSTL(atoms%label(nType))),"_",l_name(l),"/"
Henning Janssen's avatar
Henning Janssen committed
460 461 462 463 464 465 466
      ELSE
         folder=""
      ENDIF
      xPath = TRIM(ADJUSTL(xPath)) // "/" // folder

   END SUBROUTINE hubbard1_path

Henning Janssen's avatar
Henning Janssen committed
467
   SUBROUTINE writeSelfenElement(selfen,e,ef,nz,ns)
468

Henning Janssen's avatar
Henning Janssen committed
469 470
      USE m_constants
      
471 472
      IMPLICIT NONE

Henning Janssen's avatar
Henning Janssen committed
473
      INTEGER,       INTENT(IN)  :: nz,ns
Henning Janssen's avatar
Henning Janssen committed
474
      REAL,          INTENT(IN)  :: ef
Henning Janssen's avatar
Henning Janssen committed
475
      COMPLEX,       INTENT(IN)  :: selfen(2*ns,2*ns,nz)
476 477
      COMPLEX,       INTENT(IN)  :: e(nz)

Henning Janssen's avatar
Henning Janssen committed
478
      INTEGER iz,i,io_error
479
      CHARACTER(len=300) file
Henning Janssen's avatar
Henning Janssen committed
480
      COMPLEX up, down
481

Henning Janssen's avatar
Henning Janssen committed
482
      OPEN(unit=3456,file="selfen",status="replace",iostat = io_error)
483 484

      DO iz = 1, nz
Henning Janssen's avatar
Henning Janssen committed
485 486 487 488 489 490 491
         up = 0.0 
         down = 0.0
         DO i = 1, ns 
            up = up + selfen(i,i,iz)
            down = down + selfen(i+ns,i+ns,iz)
         ENDDO
         WRITE(3456,"(5f14.8)") REAL(e(iz))*hartree_to_ev_const, -AIMAG(up), -AIMAG(down), REAL(up), REAL(down)
492 493 494 495 496 497 498
      ENDDO

      CLOSE(unit=3456)


   END SUBROUTINE writeSelfenElement

Henning Janssen's avatar
Henning Janssen committed
499
END MODULE m_hubbard1_setup