Line data Source code
1 : MODULE m_hubbard1_setup
2 :
3 : USE m_juDFT
4 : USE m_types
5 : USE m_constants
6 : USE m_uj2f
7 : USE m_doubleCounting
8 : USE m_hubbard1Distance
9 : USE m_polangle
10 : USE m_hubbard1_io
11 : USE m_types_selfen
12 : USE m_add_selfen
13 : USE m_mpi_bc_tool
14 : USE m_greensf_io
15 : USE m_rotMMPmat
16 : use m_xmlOutput
17 : #ifdef CPP_EDSOLVER
18 : USE EDsolver, only: EDsolver_from_cfg
19 : #endif
20 : #ifdef CPP_MPI
21 : use mpi
22 : #endif
23 : IMPLICIT NONE
24 :
25 : CHARACTER(len=30), PARAMETER :: hubbard1CalcFolder = "Hubbard1"
26 : CHARACTER(len=30), PARAMETER :: hubbard1Outfile = "out"
27 :
28 : CHARACTER(len=30), PARAMETER :: cfg_file_ccf = "somefile"
29 : CHARACTER(len=30), PARAMETER :: cfg_file_bath = "somefile"
30 :
31 : CONTAINS
32 :
33 0 : SUBROUTINE hubbard1_setup(atoms,cell,gfinp,hub1inp,input,fmpi,noco,kpts,sphhar,sym,nococonv,pot,gdft,hub1data,results,den)
34 :
35 : TYPE(t_atoms), INTENT(IN) :: atoms
36 : TYPE(t_cell), INTENT(IN) :: cell
37 : TYPE(t_gfinp), INTENT(IN) :: gfinp
38 : TYPE(t_hub1inp), INTENT(IN) :: hub1inp
39 : TYPE(t_input), INTENT(IN) :: input
40 : TYPE(t_mpi), INTENT(IN) :: fmpi
41 : TYPE(t_noco), INTENT(IN) :: noco
42 : TYPE(t_kpts), INTENT(IN) :: kpts
43 : type(t_sphhar), intent(in) :: sphhar
44 : type(t_sym), intent(in) :: sym
45 : TYPE(t_nococonv), INTENT(IN) :: nococonv
46 : TYPE(t_potden), INTENT(IN) :: pot
47 : TYPE(t_greensf), INTENT(IN) :: gdft(:) !green's function calculated from the Kohn-Sham system
48 : TYPE(t_hub1data), INTENT(INOUT) :: hub1data
49 : TYPE(t_results), INTENT(INOUT) :: results
50 : TYPE(t_potden), INTENT(INOUT) :: den
51 :
52 : LOGICAL, PARAMETER :: l_mix = .FALSE.
53 :
54 : INTEGER :: i_hia,nType,l,occDFT_INT,ispin,m,i_exc,n,i_u
55 : INTEGER :: io_error,ierr
56 : INTEGER :: indStart,indEnd
57 : INTEGER :: hubbardioUnit
58 : INTEGER :: n_hia_task,extra,i_hia_start,i_hia_end
59 : REAL :: U,J,mx,my,mz,alpha_mix,beta,alpha
60 : COMPLEX :: offdtrace
61 : LOGICAL :: l_firstIT_HIA,l_ccfexist,l_bathexist,l_amf
62 :
63 : CHARACTER(len=300) :: folder
64 0 : TYPE(t_greensf),ALLOCATABLE :: gdft_rot(:)
65 0 : TYPE(t_greensf),ALLOCATABLE :: gu(:)
66 0 : TYPE(t_selfen), ALLOCATABLE :: selfen(:)
67 :
68 : #ifdef CPP_HDF
69 : INTEGER(HID_T) :: greensf_fileID
70 : #endif
71 :
72 0 : REAL :: mu_dc(input%jspins)
73 0 : REAL :: f0(atoms%n_hia),f2(atoms%n_hia)
74 0 : REAL :: f4(atoms%n_hia),f6(atoms%n_hia)
75 0 : REAL :: diffalpha(atoms%n_hia), diffbeta(atoms%n_hia)
76 0 : REAL :: occDFT(atoms%n_hia,input%jspins)
77 0 : COMPLEX :: mmpMat(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_hia,3)
78 : COMPLEX :: dcpot(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,3)
79 0 : COMPLEX, ALLOCATABLE :: e(:)
80 0 : COMPLEX, ALLOCATABLE :: ctmp(:)
81 : character(len=30) :: attributes(7)
82 :
83 : !Check if the EDsolver library is linked
84 : #ifndef CPP_EDSOLVER
85 0 : CALL juDFT_error("No solver linked for Hubbard 1", hint="Link the edsolver library",calledby="hubbard1_setup")
86 : #endif
87 :
88 0 : ALLOCATE(gdft_rot(atoms%n_hia))
89 0 : DO i_hia = 1, atoms%n_hia
90 0 : CALL gdft_rot(i_hia)%init(gdft(i_hia)%elem,gfinp,atoms,input,contour_in=gdft(i_hia)%contour)
91 0 : gdft_rot(i_hia) = gdft(i_hia)
92 : ENDDO
93 :
94 :
95 0 : IF(fmpi%irank.EQ.0) THEN
96 0 : write(attributes(1),'(i0)') hub1data%iter
97 0 : call openXMLElement('hubbard1Iteration',['number'], attributes(:1))
98 : !-------------------------------------------
99 : ! Create the Input for the Hubbard 1 Solver
100 : !-------------------------------------------
101 0 : CALL SYSTEM('mkdir -p ' // TRIM(ADJUSTL(hubbard1CalcFolder)))
102 :
103 : !Positions of the DFT+HIA elements in all DFT+U related arrays
104 0 : indStart = atoms%n_u+1
105 0 : indEnd = atoms%n_u+atoms%n_hia
106 :
107 : ! calculate slater integrals from u and j
108 0 : CALL uj2f(input%jspins,atoms%lda_u(indStart:indEnd),atoms%n_hia,f0,f2,f4,f6)
109 :
110 0 : DO i_hia = 1, atoms%n_hia
111 :
112 0 : l = atoms%lda_u(atoms%n_u+i_hia)%l
113 0 : U = atoms%lda_u(atoms%n_u+i_hia)%U
114 0 : J = atoms%lda_u(atoms%n_u+i_hia)%J
115 0 : l_amf = atoms%lda_u(atoms%n_u+i_hia)%l_amf
116 0 : nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
117 :
118 0 : IF(ALL(ABS(gdft(i_hia)%gmmpMat).LT.1e-12)) THEN
119 0 : CALL juDFT_error("Hubbard-1 has no DFT greensf available",calledby="hubbard1_setup")
120 : ENDIF
121 :
122 : !Create Subfolder (if there are multiple Hubbard 1 procedures)
123 0 : CALL hubbard1_path(atoms,i_hia,folder)
124 0 : CALL SYSTEM('mkdir -p ' // TRIM(ADJUSTL(folder)))
125 :
126 : !-------------------------------------------------------
127 : ! Calculate the DFT occupation of the correlated shell
128 : !-------------------------------------------------------
129 0 : mmpMat(:,:,i_hia,:) = gdft(i_hia)%occupationMatrix(gfinp,input,atoms,noco,nococonv)
130 :
131 : !For the first iteration we can fix the occupation and magnetic moments in the inp.xml file
132 0 : l_firstIT_HIA = hub1data%iter.EQ.1 .AND.ALL(ABS(den%mmpMat(:,:,indStart:indEnd,:)).LT.1e-12)
133 : IF(l_firstIT_HIA) THEN
134 0 : IF(hub1inp%init_occ(i_hia) > -9e98) THEN
135 0 : occDFT(i_hia,1) = MIN(2.0*l+1.0,hub1inp%init_occ(i_hia))
136 0 : IF(input%jspins.EQ.2) occDFT(i_hia,2) = MAX(0.0,hub1inp%init_occ(i_hia)-2.0*l-1.0)
137 : ELSE
138 0 : occDFT(i_hia,:) = 0.0
139 0 : DO ispin = 1, input%jspins
140 0 : DO m = -l, l
141 0 : occDFT(i_hia,ispin) = occDFT(i_hia,ispin) + REAL(mmpMat(m,m,i_hia,ispin))
142 : ENDDO
143 : ENDDO
144 : ENDIF
145 :
146 : !For the first iteration we just consider the diagonal elements of the density matrix
147 0 : mmpMat(:,:,i_hia,:) = cmplx_0
148 0 : do ispin = 1, input%jspins
149 0 : do m = -l, l
150 0 : mmpMat(m,m,i_hia,ispin) = occDFT(i_hia,ispin)/real(2*l+1)
151 : enddo
152 : enddo
153 :
154 0 : DO i_exc = 1, hub1inp%n_exc(i_hia)
155 0 : IF(hub1inp%init_mom(i_hia,i_exc) > -9e98) THEN
156 0 : hub1data%mag_mom(i_hia,i_exc) = hub1inp%init_mom(i_hia,i_exc)
157 : ENDIF
158 : ENDDO
159 : ELSE
160 0 : occDFT(i_hia,:) = 0.0
161 0 : DO ispin = 1, input%jspins
162 0 : DO m = -l, l
163 0 : occDFT(i_hia,ispin) = occDFT(i_hia,ispin) + REAL(mmpMat(m,m,i_hia,ispin))
164 : ENDDO
165 : ENDDO
166 : ENDIF
167 :
168 0 : IF (noco%l_unrestrictMT(nType).OR.noco%l_spinoffd_ldau(ntype) .and. gfinp%l_mperp) then
169 : !Calculate local magnetization vector from greens function
170 0 : mz = occDFT(i_hia,1) - occDFT(i_hia,2)
171 0 : offdtrace = cmplx_0
172 0 : DO m = -l, l
173 0 : offdtrace = offdtrace + mmpMat(m,m,i_hia,3)
174 : ENDDO
175 0 : mx = 2.0 * REAL(offdtrace)
176 0 : my = 2.0 * AIMAG(offdtrace)
177 0 : CALL pol_angle(mx,my,mz,beta,alpha)
178 :
179 0 : diffbeta(i_hia) = nococonv%beta(nType) - beta
180 0 : diffalpha(i_hia) = nococonv%alph(nType) - alpha
181 :
182 : !TODO: ROTATE GREENS FUNCTION SPINFRAME TO ALIGN WITH CURRENT MAGNETIZATION
183 0 : CALL gdft_rot(i_hia)%rotate_euler_angles(atoms,diffalpha(i_hia),diffbeta(i_hia),0.0,spin_rotation=.TRUE.)
184 : ENDIF
185 :
186 : !Nearest Integer occupation
187 0 : occDFT_INT = ANINT(SUM(occDFT(i_hia,:)))
188 :
189 : !Initial Information (We are already on irank 0)
190 0 : WRITE(oUnit,9010) nType
191 : 9010 FORMAT(/,"Setup for Hubbard 1 solver for atom ", I3, ": ")
192 0 : WRITE(oUnit,"(A)") "Everything related to the solver (e.g. mu_dc) is given in eV"
193 0 : WRITE(oUnit,"(/,A)") "Occupation from DFT-Green's function:"
194 0 : WRITE(oUnit,9020) 'spin-up','spin-dn'
195 : 9020 FORMAT(TR8,A7,TR3,A7)
196 0 : WRITE(oUnit,9030) occDFT(i_hia,:)
197 : 9030 FORMAT(TR7,f8.4,TR2,f8.4)
198 :
199 : !--------------------------------------------------------------------------
200 : ! Calculate the chemical potential for the solver
201 : ! This is equal to the double-counting correction used in DFT+U
202 : !--------------------------------------------------------------------------
203 : ! V_FLL = U (n - 1/2) - J (n - 1) / 2
204 : ! V_AMF = U n/2 + 2l/[2(2l+1)] (U-J) n
205 : !--------------------------------------------------------------------------
206 : IF(l_mix) alpha_mix = doubleCountingMixFactor(mmpMat(:,:,i_hia,:), l, SUM(occDFT(i_hia,:)))
207 : dcpot = doubleCountingPot(mmpMat(:,:,i_hia,:),atoms%lda_u(atoms%n_u+i_hia), U,J,any(noco%l_unrestrictMT).OR.input%ldauSpinoffd,l_mix,hub1data%l_performSpinavg,&
208 0 : alpha_mix, l_write=fmpi%irank==0)
209 :
210 0 : mu_dc = 0.0
211 0 : DO ispin = 1, input%jspins
212 0 : DO m=-l,l
213 0 : mu_dc(ispin) = mu_dc(ispin) + REAL(dcpot(m,m,ispin))/REAL(2*l+1)
214 : ENDDO
215 : ENDDO
216 :
217 : !-------------------------------------------------------
218 : ! Check for additional input files
219 : !-------------------------------------------------------
220 : !Is a crystal field matrix present in the work directory (overwrites the calculated matrix)
221 0 : INQUIRE(file=TRIM(ADJUSTL(cfg_file_ccf)),exist=l_ccfexist)
222 0 : IF(l_ccfexist) CALL read_ccfmat(hub1data%ccfmat(i_hia,-l:l,-l:l),l)
223 : !Is a bath parameter file present
224 0 : INQUIRE(file=TRIM(ADJUSTL(cfg_file_bath)),exist=l_bathexist)
225 : !Copy the bath file to the Hubbard 1 solver if its present
226 0 : IF(l_bathexist) CALL SYSTEM('cp ' // TRIM(ADJUSTL(cfg_file_bath)) // ' ' // TRIM(ADJUSTL(folder)))
227 :
228 : !Create XML output for the solver parameters
229 0 : write(attributes(1), '(i0)') nType
230 0 : write(attributes(2), '(i0)') l
231 0 : write(attributes(3), '(f14.8)') mu_dc(1)
232 0 : write(attributes(4), '(i0)') occDFT_INT
233 0 : write(attributes(5), '(a)') 'eV'
234 0 : call openXMLElement('solverParameters', ['atomType ', 'l ', 'chemPot ', 'occupation ', 'energy_unit'], attributes(:5))
235 :
236 0 : write(attributes(1), '(f14.8)') f0(i_hia)
237 0 : write(attributes(2), '(f14.8)') f2(i_hia)
238 0 : write(attributes(3), '(f14.8)') f4(i_hia)
239 0 : write(attributes(4), '(f14.8)') f6(i_hia)
240 0 : call writeXMLElementNoAttributes('slaterParameters',attributes(:4))
241 :
242 0 : do i_exc = 1, hub1inp%n_exc(i_hia)
243 0 : write(attributes(1), '(i0)') hub1inp%exc_l(i_hia,i_exc)
244 0 : write(attributes(2), '(f14.8)') hub1inp%exc(i_hia,i_exc)
245 0 : write(attributes(3), '(f14.8)') hub1data%mag_mom(i_hia,i_exc)
246 0 : call writeXMLElement('exchange',['l ', 'J ', 'moment'],attributes(:3))
247 : enddo
248 0 : write(attributes(1), '(f14.8)') hub1data%xi(i_hia)
249 0 : call writeXMLElement('socParameter',['value'],attributes(:1))
250 :
251 0 : IF(ABS(hub1inp%ccf(i_hia)).GT.1e-12) THEN
252 0 : write(attributes(1), '(f14.8)') hub1inp%ccf(i_hia)
253 : CALL writeXMLElementMatrixFormPoly('crystalField',['factor'],attributes(:1),&
254 : reshape([6,14],[1,2]),&
255 0 : hub1data%ccfmat(i_hia,-l:l,-l:l)*hartree_to_ev_const*hub1inp%ccf(i_hia))
256 : ENDIF
257 :
258 0 : call closeXMLElement('solverParameters')
259 : !-------------------------------------------------------
260 : ! Write the main config files
261 : !-------------------------------------------------------
262 : CALL write_hubbard1_input(folder,i_hia,l,f0(i_hia),f2(i_hia),f4(i_hia),f6(i_hia),&
263 0 : hub1inp,hub1data,mu_dc(1),occDFT_INT,l_bathexist)
264 : ENDDO
265 : ENDIF !fmpi%irank == 0
266 :
267 0 : IF(fmpi%irank.EQ.0) THEN
268 0 : WRITE(*,*) "Calculating new density matrix ..."
269 : ENDIF
270 :
271 : !Argument order different because occDFT is not allocatable
272 0 : CALL mpi_bc(0,fmpi%mpi_comm,occDFT)
273 :
274 : !Initializations
275 0 : ALLOCATE(gu(atoms%n_hia))
276 0 : ALLOCATE(selfen(atoms%n_hia))
277 0 : DO i_hia = 1, atoms%n_hia
278 0 : CALL gu(i_hia)%init(gdft(i_hia)%elem,gfinp,atoms,input,contour_in=gdft(i_hia)%contour)
279 0 : CALL gdft_rot(i_hia)%mpi_bc(fmpi%mpi_comm)
280 : CALL selfen(i_hia)%init(lmaxU_const,gdft(i_hia)%contour%nz,input%jspins,&
281 0 : noco%l_noco.AND.(noco%l_soc.OR.gfinp%l_mperp).OR.hub1inp%l_fullmatch)
282 : ENDDO
283 :
284 : #ifdef CPP_MPI
285 : !distribute the individual hubbard1 elements over the ranks
286 0 : n_hia_task = FLOOR(REAL(atoms%n_hia)/(fmpi%isize))
287 0 : extra = atoms%n_hia - n_hia_task*fmpi%isize
288 0 : i_hia_start = fmpi%irank*n_hia_task + 1 + extra
289 0 : i_hia_end =(fmpi%irank+1)*n_hia_task + extra
290 0 : IF(fmpi%irank < extra) THEN
291 0 : i_hia_start = i_hia_start - (extra - fmpi%irank)
292 0 : i_hia_end = i_hia_end - (extra - fmpi%irank - 1)
293 : ENDIF
294 : #else
295 : i_hia_start = 1
296 : i_hia_end = atoms%n_hia
297 : #endif
298 :
299 : #ifdef CPP_MPI
300 : !Make sure that the ranks are synchronized
301 0 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr)
302 : #endif
303 :
304 0 : mmpMat = cmplx_0
305 : !------------------------------------------------------------
306 : ! This loop runs the solver
307 : !------------------------------------------------------------
308 0 : DO i_hia = i_hia_start, i_hia_end
309 :
310 0 : IF(i_hia > atoms%n_hia .OR. i_hia < 1) CYCLE
311 :
312 0 : nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
313 0 : l = atoms%lda_u(atoms%n_u+i_hia)%l
314 :
315 0 : CALL hubbard1_path(atoms,i_hia,folder)
316 :
317 0 : ALLOCATE(e(gdft(i_hia)%contour%nz),source=cmplx_0)
318 :
319 0 : CALL timestart("Hubbard 1: EDsolver")
320 : !We have to change into the Hubbard1 directory so that the solver routines can read the config
321 0 : CALL CHDIR(TRIM(ADJUSTL(folder)))
322 : #ifdef CPP_EDSOLVER
323 : !Open the output file for the solver
324 : hubbardioUnit = 4000+i_hia
325 : OPEN(unit=hubbardioUnit, file=TRIM(ADJUSTL(hubbard1Outfile)),&
326 : status="replace", action="write", iostat=io_error)
327 : IF(io_error/=0) CALL juDFT_error("Error in opening EDsolver out file",calledby="hubbard1_setup")
328 : e = gdft(i_hia)%contour%e*hartree_to_ev_const
329 : CALL EDsolver_from_cfg(2*(2*l+1),gdft(i_hia)%contour%nz,e,selfen(i_hia)%data(:,:,:,1),1,hubbardioUnit)
330 : !---------------------------------------------------
331 : ! Calculate selfenergy on lower contour explicitly
332 : ! Mainly out of paranoia :D
333 : ! No rediagonalization (last argument switches this)
334 : !---------------------------------------------------
335 : e = conjg(gdft(i_hia)%contour%e)*hartree_to_ev_const
336 : CALL EDsolver_from_cfg(2*(2*l+1),gdft(i_hia)%contour%nz,e,selfen(i_hia)%data(:,:,:,2),0,hubbardioUnit)
337 : CLOSE(hubbardioUnit, iostat=io_error)
338 : IF(io_error/=0) CALL juDFT_error("Error in closing EDsolver out file",calledby="hubbard1_setup")
339 : #endif
340 0 : IF(atoms%n_hia>1) THEN
341 0 : CALL CHDIR("../../")
342 : ELSE
343 0 : CALL CHDIR("../")
344 : ENDIF
345 0 : CALL timestop("Hubbard 1: EDsolver")
346 :
347 0 : DEALLOCATE(e)
348 :
349 : !-------------------------------------------
350 : ! Postprocess selfenergy
351 : !-------------------------------------------
352 0 : CALL selfen(i_hia)%postProcess(noco,nococonv,nType,l,input%jspins,pot%mmpMat(:,:,atoms%n_u+i_hia,:))
353 :
354 : !----------------------------------------------------------------------
355 : ! Solution of the Dyson Equation
356 : !----------------------------------------------------------------------
357 : ! G(z)^(-1) = G_0(z)^(-1) - mu - Sigma(z)
358 : !----------------------------------------------------------------------
359 : ! Sigma(z) is the self-energy from the impurity solver
360 : ! We introduce an additional chemical potential mu, which is determined
361 : ! so that the occupation of the correlated orbital does not change
362 : !----------------------------------------------------------------------
363 : #ifdef CPP_DEBUG
364 : OPEN(unit=1337, file=TRIM(ADJUSTL(folder)) // 'mu',&
365 : status="replace", action="write", iostat=io_error)
366 : #endif
367 :
368 0 : CALL timestart("Hubbard 1: Add Selfenergy")
369 : CALL add_selfen(gdft_rot(i_hia),selfen(i_hia),gfinp,input,atoms,noco,nococonv,&
370 0 : occDFT(i_hia,:),gu(i_hia),mmpMat(:,:,i_hia,:))
371 0 : CALL timestop("Hubbard 1: Add Selfenergy")
372 :
373 : #ifdef CPP_DEBUG
374 : CLOSE(unit=1337)
375 : #endif
376 :
377 : ENDDO
378 :
379 0 : IF(fmpi%irank.EQ.0) THEN
380 0 : WRITE(oUnit,*)
381 0 : WRITE(oUnit,'(A)') "Calculated mu to match Self-energy to DFT-GF"
382 : ENDIF
383 : !Collect the impurity Green's Function
384 0 : DO i_hia = 1, atoms%n_hia
385 0 : CALL gu(i_hia)%collect(fmpi%mpi_comm)
386 0 : CALL selfen(i_hia)%collect(fmpi%mpi_comm)
387 0 : IF(fmpi%irank.EQ.0) THEN
388 : !We found the chemical potential to within the desired accuracy
389 0 : WRITE(oUnit,*) 'i_hia: ',i_hia, " muMatch = ", selfen(i_hia)%muMatch(:)
390 : ENDIF
391 : ENDDO
392 :
393 : #ifdef CPP_MPI
394 : !Collect the density matrix to rank 0
395 0 : n = SIZE(mmpMat)
396 0 : ALLOCATE(ctmp(n))
397 0 : CALL MPI_REDUCE(mmpMat,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,0,MPI_COMM_WORLD,ierr)
398 0 : IF(fmpi%irank.EQ.0) CALL zcopy(n,ctmp,1,mmpMat,1)
399 0 : DEALLOCATE(ctmp)
400 : #endif
401 :
402 : !--------------------------------------------------------------------
403 : ! Calculate Distances from last density matrix and update den%mmpmat
404 : !--------------------------------------------------------------------
405 0 : results%last_mmpMatdistance = 0.0
406 0 : results%last_occdistance = 0.0
407 :
408 0 : IF(fmpi%irank.EQ.0) THEN
409 0 : DO i_hia = 1, atoms%n_hia
410 0 : nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
411 0 : l = atoms%lda_u(atoms%n_u+i_hia)%l
412 0 : IF (noco%l_unrestrictMT(nType) .OR. noco%l_spinoffd_ldau(ntype) .and. gfinp%l_mperp) then
413 : !TODO: ROTATE GREENS FUNCTION AND mmpmat SPINFRAME TO LOCAL FRAME
414 0 : CALL gu(i_hia)%rotate_euler_angles(atoms,-diffalpha(i_hia),-diffbeta(i_hia),0.0,spin_rotation=.TRUE.)
415 0 : mmpMat(:,:,i_hia,:) = rotMMPmat(mmpMat(:,:,i_hia,:),-diffalpha(i_hia),-diffbeta(i_hia),0.0,l,spin_rotation=.TRUE.)
416 : ENDIF
417 :
418 0 : CALL hubbard1Distance(den%mmpMat(:,:,atoms%n_u+i_hia,:),mmpMat(:,:,i_hia,:),results)
419 0 : DO ispin = 1, MERGE(3,input%jspins,gfinp%l_mperp)
420 0 : den%mmpMat(-lmaxU_const:,-lmaxU_const:,atoms%n_u+i_hia,ispin) = mmpMat(-lmaxU_const:,-lmaxU_const:,i_hia,ispin)
421 : ENDDO
422 : ENDDO
423 : ENDIF
424 :
425 : #ifdef CPP_HDF
426 0 : IF(fmpi%irank.EQ.0) THEN
427 : !------------------------------
428 : !Write out DFT Green's Function
429 : !------------------------------
430 0 : CALL timestart("Hubbard 1: IO/Write")
431 0 : CALL openGreensFFile(greensf_fileID, input, gfinp, atoms, sym, cell, kpts,sphhar, inFilename="greensf_DFT.hdf")
432 : CALL writeGreensFData(greensf_fileID, input, gfinp, atoms, nococonv, noco, cell,&
433 0 : GREENSF_HUBBARD_CONST, gdft, mmpmat)
434 0 : CALL closeGreensFFile(greensf_fileID)
435 :
436 : !-------------------------------------
437 : !Write out correlated Green's Function
438 : !-------------------------------------
439 0 : CALL openGreensFFile(greensf_fileID, input, gfinp, atoms, sym, cell, kpts,sphhar, inFilename="greensf_IMP.hdf")
440 : CALL writeGreensFData(greensf_fileID, input, gfinp, atoms, nococonv, noco, cell,&
441 0 : GREENSF_HUBBARD_CONST, gu, mmpmat,selfen=selfen)
442 0 : CALL closeGreensFFile(greensf_fileID)
443 0 : CALL timestop("Hubbard 1: IO/Write")
444 : ENDIF
445 : #endif
446 :
447 : !Broadcast the density matrix
448 0 : CALL mpi_bc(den%mmpMat,0,fmpi%mpi_comm)
449 0 : CALL mpi_bc(results%last_occdistance,0,fmpi%mpi_comm)
450 0 : CALL mpi_bc(results%last_mmpMatdistance,0,fmpi%mpi_comm)
451 :
452 0 : IF(fmpi%irank.EQ.0) THEN
453 0 : WRITE(*,*) "Hubbard 1 Iteration: ", hub1data%iter
454 0 : WRITE(*,*) "Distances: "
455 0 : WRITE(*,*) "-----------------------------------------------------"
456 0 : WRITE(*,*) "Occupation Distance: " , results%last_occdistance
457 0 : WRITE(*,*) "Element Distance: " , results%last_mmpMatdistance
458 0 : WRITE(*,*) "-----------------------------------------------------"
459 0 : WRITE(oUnit,*) "nmmp occupation distance: ", results%last_occdistance
460 0 : WRITE(oUnit,*) "nmmp element distance: ", results%last_mmpMatdistance
461 0 : WRITE(oUnit,FMT=8140) hub1data%iter
462 : 8140 FORMAT (/,5x,'******* Hubbard 1 it=',i3,' is completed********',/,/)
463 0 : IF(hub1inp%l_correctEtot.AND..NOT.hub1inp%l_dftspinpol) THEN
464 0 : IF(results%last_mmpMatdistance <= hub1inp%minMatDistance .AND. &
465 : results%last_occdistance <= hub1inp%minOccDistance) THEN
466 :
467 0 : WRITE(*,*) 'Density matrix converged: switching off Spin averaging in DFT'
468 0 : hub1data%l_performSpinavg = .FALSE.
469 : ENDIF
470 : ENDIF
471 0 : CALL openXMLElementNoAttributes('hubbard1DensityMatrix')
472 0 : DO ispin = 1, SIZE(den%mmpMat,4)
473 0 : DO i_hia = 1, atoms%n_hia
474 0 : i_u = atoms%n_u + i_hia
475 0 : attributes = ''
476 0 : WRITE(attributes(1),'(i0)') ispin
477 0 : WRITE(attributes(2),'(i0)') atoms%lda_u(i_u)%atomType
478 0 : WRITE(attributes(3),'(i0)') i_hia
479 0 : WRITE(attributes(4),'(i0)') atoms%lda_u(i_u)%l
480 0 : WRITE(attributes(5),'(f15.8)') atoms%lda_u(i_u)%u
481 0 : WRITE(attributes(6),'(f15.8)') atoms%lda_u(i_u)%j
482 0 : WRITE(attributes(7),'(2f15.8)') selfen(i_hia)%muMatch(:)
483 : CALL writeXMLElementMatrixPoly('densityMatrixFor',&
484 : ['spin ','atomType','hiaIndex','l ','U ','J ','muMatch '],&
485 0 : attributes,den%mmpMat(-l:l,-l:l,i_u,ispin))
486 : END DO
487 : END DO
488 0 : CALL closeXMLElement('hubbard1DensityMatrix')
489 0 : write(attributes(1),'(f14.8)') results%last_occdistance
490 0 : write(attributes(2),'(f14.8)') results%last_mmpMatdistance
491 0 : call writeXMLElement('hubbard1Distance',['occupationDistance','elementDistance '], attributes(:2))
492 0 : call closeXMLElement('hubbard1Iteration')
493 : ENDIF
494 :
495 0 : CALL mpi_bc(hub1data%l_performSpinavg,0,fmpi%mpi_comm)
496 :
497 0 : END SUBROUTINE hubbard1_setup
498 :
499 0 : SUBROUTINE hubbard1_path(atoms,i_hia,xPath)
500 :
501 : !Defines the folder structure
502 : ! The Solver is run in the subdirectories
503 : ! Hubbard1/ if only one Hubbard1 procedure is run
504 : ! Hubbard1/atom_label_l if there are more
505 :
506 : TYPE(t_atoms), INTENT(IN) :: atoms
507 : INTEGER, INTENT(IN) :: i_hia
508 : CHARACTER(len=300), INTENT(OUT) :: xPath
509 :
510 : CHARACTER(len=300) :: folder,fmt
511 : CHARACTER(len=1),PARAMETER :: spdfg(0:4) = ['s','p','d','f','g']
512 : INTEGER nType,l
513 :
514 0 : nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
515 0 : l = atoms%lda_u(atoms%n_u+i_hia)%l
516 0 : xPath = TRIM(ADJUSTL(hubbard1CalcFolder))
517 0 : IF(atoms%n_hia>1) THEN
518 0 : WRITE(fmt,'("(A",I2.2,",A1,A1,A1)")') LEN(TRIM(ADJUSTL(atoms%label(nType))))
519 0 : WRITE(folder,fmt) TRIM(ADJUSTL(atoms%label(nType))),"_",spdfg(l),"/"
520 : ELSE
521 0 : folder=""
522 : ENDIF
523 0 : xPath = TRIM(ADJUSTL(xPath)) // "/" // folder
524 :
525 0 : END SUBROUTINE hubbard1_path
526 :
527 : END MODULE m_hubbard1_setup
|
