Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_fleur
7 : !! Legacy comment from before reformatting:
8 : !!
9 : !! Based on flapw7 (c.l.fu, m.weinert, e.wimmer):
10 : !! full potential linearized augmented plane wave method for thin
11 : !! films and superlattices (version 7 ---- general symmetry)
12 : !! symmetry part --- e.wimmer
13 : !! potential generator --- c.l.fu,r.podloucky
14 : !! matrix elements --- m.weinert
15 : !! charge density --- c.l.fu
16 : !! c.l.fu 1987
17 : !!
18 : !! 2nd variation diagon. --- r.-q. wu 1992
19 : !! forces a la Yu et al --- r.podloucky 1995
20 : !! full relativistic core --- a.shick 1996
21 : !! broyden mixing --- r.pentcheva 1996
22 : !! gga (pw91, pbe) --- t.asada 1997
23 : !! local orbitals --- p.kurz 1997
24 : !! automatic symmetry --- w.hofer 1997
25 : !! core tails & start --- r.abt 1998
26 : !! spin orbit coupling --- a.shick,x.nie 1998
27 : !! non-colinear magnet. --- p.kurz 1999
28 : !! one-dimensional --- y.mokrousov 2002
29 : !! exchange parameters --- m.lezaic 2004
30 : !! g.bihlmayer, s.bluegel 1999
31 :
32 : IMPLICIT NONE
33 : CONTAINS
34 160 : SUBROUTINE fleur_execute(fmpi, fi, sphhar, stars, nococonv, forcetheo, enpara, results, &
35 : xcpot, wann, hybdat, mpdata)
36 : !! This routine is the main program of the FLEUR code.
37 :
38 : USE m_types
39 : USE m_types_forcetheo_extended
40 : USE m_constants
41 : USE m_optional
42 : USE m_cdn_io
43 : USE m_mixing_history
44 : USE m_qfix
45 : USE m_vgen
46 : USE m_vgen_coulomb
47 : USE m_writexcstuff
48 : USE m_eigen
49 : USE m_eigenso
50 : USE m_fermie
51 : USE m_cdngen
52 : USE m_totale
53 : USE m_potdis
54 : USE m_mix
55 : USE m_xmlOutput
56 : USE m_juDFT_time
57 : USE m_calc_hybrid
58 : USE m_rdmft
59 : USE m_io_hybrid
60 : USE m_wann_optional
61 : USE m_wannier
62 : USE m_bs_comfort
63 : USE m_dwigner
64 : USE m_ylm
65 : USE m_metagga
66 : USE m_plot
67 : USE m_usetup
68 : USE m_hubbard1_setup
69 : USE m_writeCFOutput
70 : USE m_mpi_bc_tool
71 : USE m_eig66_io
72 : USE m_chase_diag
73 : USE m_writeBasis
74 : USE m_RelaxSpinAxisMagn
75 : USE m_dfpt
76 :
77 :
78 : !$ USE omp_lib
79 :
80 : TYPE(t_mpi), INTENT(INOUT) :: fmpi
81 : TYPE(t_fleurinput), INTENT(IN) :: fi
82 : CLASS(t_xcpot), INTENT(IN) :: xcpot
83 : TYPE(t_sphhar), INTENT(IN) :: sphhar
84 : TYPE(t_stars), INTENT(IN) :: stars
85 : TYPE(t_nococonv), INTENT(INOUT) :: nococonv
86 : TYPE(t_results), INTENT(INOUT) :: results
87 : TYPE(t_wann), INTENT(INOUT) :: wann
88 :
89 : CLASS(t_forcetheo), INTENT(INOUT) :: forcetheo
90 : TYPE(t_enpara), INTENT(INOUT) :: enpara
91 : TYPE(t_hybdat), INTENT(INOUT) :: hybdat
92 : TYPE(t_mpdata), INTENT(INOUT) :: mpdata
93 :
94 : ! TODO: This is just fi%input with neig=2*neig and should be refactored out.
95 : TYPE(t_input) :: input_soc
96 :
97 480 : TYPE(t_field) :: field2
98 320 : TYPE(t_potden) :: vTot, vx, vCoul, vxc, exc
99 160 : TYPE(t_potden) :: inDen, outDen, EnergyDen, sliceDen
100 160 : TYPE(t_hub1data) :: hub1data
101 :
102 160 : TYPE(t_greensf), ALLOCATABLE :: greensFunction(:)
103 160 : TYPE(t_log_message) :: log
104 :
105 : INTEGER :: eig_id, archiveType, num_threads
106 : INTEGER :: iter, iterHF, i, n, i_gf
107 : INTEGER :: wannierspin
108 : LOGICAL :: l_opti, l_cont, l_qfix, l_real, l_olap, l_error, l_dummy
109 : LOGICAL :: l_forceTheorem, l_lastIter, l_exist
110 : REAL :: fix, sfscale, rdummy, tempDistance
111 : REAL :: mmpmatDistancePrev, occDistancePrev
112 : INTEGER :: tempI, tempK, tempJSP
113 :
114 : #ifdef CPP_MPI
115 : INTEGER :: ierr
116 : #endif
117 :
118 : ! Check, whether we already have a suitable density file and if not,
119 : ! generate a starting density.
120 : CALL optional(fmpi, fi%atoms, sphhar, fi%vacuum, stars, fi%input, &
121 160 : fi%sym, fi%cell, fi%sliceplot, xcpot, fi%noco)
122 :
123 154 : IF (fi%input%l_wann .AND. (fmpi%irank == 0) .AND. (.NOT. wann%l_bs_comf)) THEN
124 : ! TODO: If this warning is commented out, can it be erased?
125 : !IF(fmpi%isize.NE.1) CALL juDFT_error('No Wannier+MPI at the moment',calledby = 'fleur')
126 2 : CALL wann_optional(fmpi, fi%input, fi%kpts, fi%atoms, fi%sym, fi%cell, fi%noco, wann)
127 : END IF
128 :
129 154 : iter = 0
130 154 : iterHF = 0
131 154 : l_cont = (iter < fi%input%itmax)
132 :
133 : ! Read in last Hubbard 1 distances
134 154 : l_error = .TRUE.
135 154 : IF (fi%atoms%n_hia>0 .AND. fmpi%irank==0) CALL readPrevmmpDistances(mmpmatDistancePrev,occDistancePrev,l_error)
136 :
137 154 : CALL hub1data%init(fi%atoms, fi%input, fi%hub1inp, fmpi, mmpmatDistancePrev, occDistancePrev, l_error)
138 154 : CALL hub1data%mpi_bc(fmpi%mpi_comm)
139 :
140 154 : IF(fi%atoms%n_hia>0 .AND. .NOT.l_error .AND. .NOT.fi%hub1inp%l_forceHIAiteration) THEN
141 : !Set the current HIA distance to the read in value
142 : !Prevents too many HIA iterations after restart
143 0 : results%last_mmpmatDistance = mmpmatDistancePrev
144 0 : results%last_occDistance = occDistancePrev
145 : END IF
146 154 : CALL mpi_bc(results%last_mmpmatDistance,0,fmpi%mpi_comm)
147 154 : CALL mpi_bc(results%last_occDistance,0,fmpi%mpi_comm)
148 :
149 154 : IF (fmpi%irank==0) CALL openXMLElementNoAttributes('scfLoop')
150 :
151 : ! Initialize (and load) the input density inDen
152 154 : CALL timestart("Load inDen")
153 :
154 : ! Warning on strange choice of switches before starting density is generated.
155 400 : IF (fi%input%l_onlyMtStDen .AND. .NOT. ANY(fi%noco%l_unrestrictMT)) THEN
156 0 : CALL juDFT_warn("l_onlyMtStDen='T' and l_mtNocoPot='F' makes no sense.", calledby='types_input')
157 : END IF
158 :
159 154 : CALL inDen%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN)
160 :
161 154 : archiveType = CDN_ARCHIVE_TYPE_CDN1_const
162 154 : IF (fi%noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
163 400 : IF (ANY(fi%noco%l_unrestrictMT)) archiveType = CDN_ARCHIVE_TYPE_FFN_const
164 :
165 154 : IF (fmpi%irank==0) CALL readDensity(stars, fi%noco, fi%vacuum, fi%atoms, fi%cell, sphhar, &
166 : fi%input, fi%sym, archiveType, CDN_INPUT_DEN_const, 0, &
167 77 : results%ef, results%last_distance, l_qfix, inDen)
168 154 : call mpi_bc(results%last_distance, 0, fmpi%mpi_comm)
169 :
170 : !IF (fi%noco%l_alignMT .AND. fmpi%irank .EQ. 0) THEN
171 : !CALL initRelax(fi%noco, nococonv, fi%atoms, fi%input, fi%vacuum, sphhar, stars, fi%sym, fi%cell, inDen)
172 : !CALL doRelax(fi%vacuum, sphhar, stars, fi%sym, fi%cell, fi%noco, nococonv, fi%input, fi%atoms, inDen)
173 : !END IF
174 :
175 154 : CALL timestart("Qfix main")
176 154 : CALL qfix(fmpi, stars,nococonv, fi%atoms, fi%sym, fi%vacuum, sphhar, fi%input, fi%cell, inDen, fi%noco%l_noco, .FALSE., .FALSE., .FALSE., fix)
177 : !CALL magMoms(fi%input,fi%atoms,fi%noco,nococonv,den=inDen)
178 154 : CALL timestop("Qfix main")
179 :
180 154 : IF (fmpi%irank==0) THEN
181 : CALL writeDensity(stars, fi%noco, fi%vacuum, fi%atoms, fi%cell, sphhar, fi%input, fi%sym, archiveType, CDN_INPUT_DEN_const, &
182 77 : 0, -1.0, results%ef, results%last_mmpmatDistance, results%last_occDistance, .FALSE., inDen)
183 : END IF
184 :
185 418 : IF (ANY(fi%noco%l_alignMT)) CALL toLocalSpinFrame(fmpi, fi%vacuum, sphhar, stars, fi%sym, fi%cell, &
186 2 : fi%noco, nococonv, fi%input, fi%atoms, .TRUE., inDen, .TRUE.)
187 154 : CALL timestop("Load inDen")
188 :
189 : ! Initialize potentials
190 154 : CALL timestart("Initialize potentials")
191 154 : CALL vTot%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT)
192 154 : CALL vCoul%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTCOUL)
193 154 : CALL vx%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT)
194 154 : CALL vxc%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT)
195 154 : CALL exc%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_POTTOT)
196 154 : CALL timestop("Initialize potentials")
197 :
198 : ! Initialize Green's function
199 154 : CALL timestart("Initialize GF")
200 1592 : ALLOCATE (greensFunction(MAX(1, fi%gfinp%n)))
201 154 : IF (fi%gfinp%n > 0) THEN
202 1048 : DO i_gf = 1, fi%gfinp%n
203 1048 : CALL greensFunction(i_gf)%init(fi%gfinp%elem(i_gf), fi%gfinp, fi%atoms, fi%input, nkpt=fi%kpts%nkpt)
204 : END DO
205 : END IF
206 154 : CALL timestop("Initialize GF")
207 :
208 : ! Open/allocate eigenvector storage
209 154 : CALL timestart("Open/allocate eigenvector storage")
210 154 : IF (fi%noco%l_soc .AND. fi%input%l_wann) THEN
211 : ! Weed up and down spinor components for SOC MLWFs.
212 : ! When jspins=1 Fleur usually writes only the up-spinor into the eig-file.
213 : ! Make sure we always get up and down spinors when SOC=true.
214 0 : wannierspin = 2
215 : ELSE
216 154 : wannierspin = fi%input%jspins
217 : END IF
218 :
219 154 : l_olap = fi%hybinp%l_hybrid .OR. fi%input%l_rdmft
220 : eig_id = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, wannierspin, fi%noco%l_noco, &
221 154 : .NOT.fi%INPUT%eig66(1), fi%input%l_real, fi%noco%l_soc, fi%INPUT%eig66(1), l_olap, fmpi%n_size)
222 :
223 : ! The hybrid implementation doesn't use spinors. In the case of SOC calculations, a separate eig file has to be created
224 : ! It contains those eigenvalues/vectors found in the 'eigen' subroutine
225 154 : IF (fi%noco%l_soc .AND. fi%hybinp%l_hybrid) THEN
226 : hybdat%eig_id = open_eig(fmpi%mpi_comm, lapw_dim_nbasfcn, fi%input%neig, fi%kpts%nkpt, wannierspin, fi%noco%l_noco, &
227 0 : .NOT.fi%INPUT%eig66(1), fi%input%l_real, .FALSE., fi%INPUT%eig66(1), l_olap, fmpi%n_size)
228 : ELSE
229 154 : hybdat%eig_id = eig_id
230 : ENDIF
231 :
232 : ! TODO: Isn't this comment kind of lost here?
233 : ! Rotate cdn to local frame if specified.
234 :
235 : #ifdef CPP_CHASE
236 : CALL init_chase(fmpi, fi%input, fi%atoms, fi%kpts, fi%noco, l_real)
237 : #endif
238 :
239 154 : CALL timestop("Open/allocate eigenvector storage")
240 :
241 : ! Perform some pre-scf-loop checks (start)
242 154 : CALL timestart("Pre-scf sanity check")
243 154 : IF (fi%input%gw .EQ. 2) THEN
244 0 : IF (ABS(results%last_distance).GT.fi%input%mindistance) THEN
245 0 : CALL juDFT_warn('Performing spex="2" step without a selfconsistent density!')
246 : END IF
247 : END IF
248 154 : CALL timestop("Pre-scf sanity check")
249 :
250 : ! Start the scf loop.
251 154 : l_lastIter = .FALSE.
252 818 : scfloop: DO WHILE (l_cont)
253 686 : iter = iter + 1
254 686 : l_lastIter = l_lastIter.OR.(iter.EQ.fi%input%itmax)
255 686 : hub1data%overallIteration = hub1data%overallIteration + 1
256 :
257 686 : IF (fmpi%irank==0) CALL openXMLElementFormPoly('iteration', (/'numberForCurrentRun', 'overallNumber '/), &
258 1715 : (/iter, inden%iter/), RESHAPE((/19, 13, 5, 5/), (/2, 2/)))
259 :
260 : ! TODO: What is commented out here and should it perhaps be removed?
261 :
262 : ! !$ !+t3e
263 : ! !$ IF (fi%input%alpha.LT.10.0) THEN
264 : ! !$
265 : ! !$ IF (iter.GT.1) THEN
266 : ! !$ fi%input%alpha = fi%input%alpha - NINT(fi%input%alpha)
267 : ! !$ END IF
268 :
269 : !CALL resetIterationDependentTimers()
270 :
271 686 : CALL timestart("Iteration")
272 686 : IF (fmpi%irank==0) THEN
273 343 : WRITE (oUnit, FMT=8100) iter
274 : 8100 FORMAT(/, 10x, ' iter= ', i5)
275 : END IF !fmpi%irank==0
276 :
277 : #ifdef CPP_CHASE
278 : CALL chase_distance(results%last_distance)
279 : #endif
280 :
281 686 : CALL inDen%distribute(fmpi%mpi_comm)
282 686 : CALL nococonv%mpi_bc(fmpi%mpi_comm)
283 :
284 : ! Plot the input density if specified
285 686 : IF (fi%sliceplot%iplot .NE. 0) THEN
286 0 : IF (.NOT.fi%sliceplot%slice) THEN
287 : CALL makeplots(stars, fi%atoms, sphhar, fi%vacuum, fi%input, fmpi, fi%sym, fi%cell, &
288 0 : fi%noco, nococonv, inDen, PLOT_INPDEN, fi%sliceplot)
289 : ELSE
290 0 : CALL sliceDen%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN)
291 0 : IF (fmpi%irank .EQ. 0) CALL readDensity(stars, fi%noco, fi%vacuum, fi%atoms, fi%cell, sphhar, &
292 : fi%input, fi%sym, CDN_ARCHIVE_TYPE_CDN_const, &
293 0 : CDN_INPUT_DEN_const, 0, rdummy, tempDistance, l_dummy, sliceDen, inFilename='cdn_slice')
294 0 : CALL sliceden%distribute(fmpi%mpi_comm)
295 0 : call nococonv%mpi_bc(fmpi%mpi_comm)
296 : CALL makeplots(stars, fi%atoms, sphhar, fi%vacuum, fi%input, fmpi, fi%sym, fi%cell, &
297 0 : fi%noco, nococonv, sliceDen, PLOT_INPDEN, fi%sliceplot)
298 : END IF
299 0 : IF ((fmpi%irank .EQ. 0) .AND. (fi%sliceplot%iplot .EQ. 2)) THEN
300 0 : CALL juDFT_end("Stopped self consistency loop after plots have been generated.")
301 : END IF
302 : END IF
303 :
304 : !HF
305 686 : IF (fi%hybinp%l_hybrid) THEN
306 186 : hybdat%l_calhf = (results%last_distance >= 0.0) .AND. (results%last_distance < fi%input%minDistance)
307 : SELECT TYPE (xcpot)
308 : TYPE IS (t_xcpot_inbuild)
309 : CALL calc_hybrid(fi, mpdata, hybdat, fmpi, nococonv, stars, enpara, &
310 186 : hybdat%results, xcpot, vTot, iter, iterHF)
311 : END SELECT
312 :
313 : #ifdef CPP_MPI
314 186 : CALL MPI_Barrier(fmpi%mpi_comm, ierr)
315 : #endif
316 186 : IF (hybdat%l_calhf) THEN
317 12 : CALL mixing_history_reset(fmpi)
318 12 : iter = 0
319 : END IF
320 : END IF
321 :
322 : ! TODO: What is commented out here and should it perhaps be removed?
323 :
324 : ! !$ DO pc = 1, wann%nparampts
325 : ! !$ !---> gwf
326 : ! !$ IF (wann%l_sgwf.OR.wann%l_ms) THEN
327 : ! !$ fi%noco%qss(:) = wann%param_vec(:,pc)
328 : ! !$ fi%noco%alph(:) = wann%param_alpha(:,pc)
329 : ! !$ ELSE IF (wann%l_socgwf) THEN
330 : ! !$ IF(wann%l_dim(2)) fi%noco%phi = tpi_const * wann%param_vec(2,pc)
331 : ! !$ IF(wann%l_dim(3)) fi%noco%theta = tpi_const * wann%param_vec(3,pc)
332 : ! !$ END IF
333 : !---< gwf
334 :
335 : ! Optionally scale up the magnetization density before the potential calculation.
336 1862 : IF (ANY(fi%noco%l_unrestrictMT).AND.fi%noco%l_scaleMag) THEN
337 2 : sfscale = fi%noco%mag_scale
338 2 : CALL inDen%SpinsToChargeAndMagnetisation()
339 736796 : inDen%mt(:, 0:, :, 2:4) = sfscale*inDen%mt(:, 0:, :, 2:4)
340 12290 : inDen%pw(:, 2:3) = sfscale*inDen%pw(:, 2:3)
341 14 : inDen%vac(:, :, :, 2:3) = sfscale*inDen%vac(:, :, :, 2:3)
342 2 : CALL inDen%ChargeAndMagnetisationToSpins()
343 : END IF
344 :
345 686 : CALL timestart("generation of potential")
346 : CALL vgen(hybdat, fi%field, fi%input, xcpot, fi%atoms, sphhar, stars, fi%vacuum, fi%sym, fi%juphon, &
347 686 : fi%cell, fi%sliceplot, fmpi, results, fi%noco, nococonv, EnergyDen, inDen, vTot, vx, vCoul, vxc, exc)
348 686 : CALL timestop("generation of potential")
349 :
350 : ! Scale the magnetization back.
351 1862 : IF (ANY(fi%noco%l_unrestrictMT).AND.fi%noco%l_scaleMag) THEN
352 2 : CALL inDen%SpinsToChargeAndMagnetisation()
353 736796 : inDen%mt(:, 0:, :, 2:4) = inDen%mt(:, 0:, :, 2:4)/sfscale
354 12290 : inDen%pw(:, 2:3) = inDen%pw(:, 2:3)/sfscale
355 14 : inDen%vac(:, :, :, 2:3) = inDen%vac(:, :, :, 2:3)/sfscale
356 2 : CALL inDen%ChargeAndMagnetisationToSpins()
357 : END IF
358 :
359 : ! Set up HIA stuff.
360 686 : IF (hub1data%l_runthisiter .AND. fi%atoms%n_hia > 0) THEN
361 0 : DO i_gf = 1, fi%gfinp%n
362 0 : CALL greensFunction(i_gf)%mpi_bc(fmpi%mpi_comm)
363 : END DO
364 0 : IF (ALL(greensFunction(fi%gfinp%hiaElem)%l_calc)) THEN
365 0 : hub1data%iter = hub1data%iter + 1
366 : CALL hubbard1_setup(fi%atoms, fi%cell, fi%gfinp, fi%hub1inp, fi%input, fmpi, fi%noco, fi%kpts, sphhar, fi%sym, nococonv, vTot, &
367 0 : greensFunction(fi%gfinp%hiaElem), hub1data, results, inDen)
368 : ELSE
369 0 : IF (fmpi%irank .EQ. 0) WRITE (*, *) 'Not all Greens Functions available: Running additional iteration'
370 0 : hub1data%l_runthisiter = .FALSE. !To prevent problems in mixing later on
371 : END IF
372 : END IF
373 :
374 : ! Set up DFT+U stuff.
375 686 : IF (fi%atoms%n_u + fi%atoms%n_hia > 0) THEN
376 64 : CALL u_setup(fi%atoms, fi%input, fi%noco, fmpi, hub1data, inDen, vTot, results)
377 : END IF
378 :
379 : #ifdef CPP_MPI
380 686 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
381 : #endif
382 :
383 686 : CALL forcetheo%start(vtot, fmpi%irank==0)
384 1350 : forcetheoloop: DO WHILE (forcetheo%next_job(fmpi,l_lastIter, fi%atoms, fi%noco, nococonv))
385 :
386 686 : CALL timestart("H generation and diagonalization (total)")
387 :
388 686 : CALL timestart("eigen")
389 :
390 686 : CALL timestart("Updating energy parameters")
391 686 : CALL enpara%update(fmpi, fi%atoms, fi%vacuum, fi%input, vToT, hub1data)
392 686 : CALL timestop("Updating energy parameters")
393 :
394 686 : IF (.NOT. fi%input%eig66(1)) THEN
395 : CALL eigen(fi, fmpi, stars, sphhar, xcpot, forcetheo, enpara, nococonv, mpdata, &
396 686 : hybdat, iter, eig_id, results, inDen, vToT, vx, hub1data)
397 : END IF
398 : ! TODO: What is commented out here and should it perhaps be removed?
399 : ! !$ eig_idList(pc) = eig_id
400 686 : CALL timestop("eigen")
401 :
402 : ! Add all HF contributions to the total energy
403 686 : IF( fi%input%jspins .EQ. 1 .AND. fi%hybinp%l_hybrid ) THEN
404 132 : hybdat%results%te_hfex%valence = 2*hybdat%results%te_hfex%valence
405 132 : IF(hybdat%l_calhf) hybdat%results%te_hfex%core = 2*hybdat%results%te_hfex%core
406 : END IF
407 : ! Send all result of local total energies to the r ! TODO: Is half the comment missing?
408 686 : IF (fi%hybinp%l_hybrid .AND. hybdat%l_calhf) THEN
409 12 : results%te_hfex=hybdat%results%te_hfex
410 : #ifdef CPP_MPI
411 12 : CALL fmpi%set_root_comm()
412 12 : IF (fmpi%n_rank==0) THEN
413 6 : IF (fmpi%irank==0) THEN
414 6 : CALL MPI_Reduce(MPI_IN_PLACE, results%te_hfex%valence, 1, MPI_REAL8, MPI_SUM, 0, fmpi%root_comm, ierr)
415 : ELSE
416 0 : CALL MPI_Reduce(results%te_hfex%valence, MPI_IN_PLACE, 1, MPI_REAL8, MPI_SUM, 0, fmpi%root_comm, ierr)
417 : END IF
418 : END IF
419 : #endif
420 : END IF
421 :
422 686 : CALL timestart("2nd variation SOC")
423 686 : IF (fi%noco%l_soc .AND. .NOT. fi%noco%l_noco .AND. .NOT. fi%INPUT%eig66(1)) THEN
424 68 : IF (fi%hybinp%l_hybrid) hybdat%results = results !Store old eigenvalues for later call to fermie
425 68 : CALL eigenso(eig_id, fmpi, stars, sphhar, nococonv, vTot, enpara, results, fi%hub1inp, hub1data,fi)
426 : ENDIF
427 686 : CALL timestop("2nd variation SOC")
428 686 : CALL timestop("H generation and diagonalization (total)")
429 :
430 : #ifdef CPP_MPI
431 686 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
432 : #endif
433 :
434 : ! Fermi level and occupancies
435 686 : input_soc = fi%input
436 686 : IF (fi%noco%l_soc .AND. (.NOT. fi%noco%l_noco)) THEN
437 68 : input_soc = fi%input
438 68 : input_soc%neig = 2*fi%input%neig
439 : END IF
440 :
441 686 : IF (fi%input%gw>0) THEN
442 0 : IF (fmpi%irank==0) THEN
443 : CALL writeBasis(input_soc, fi%noco, nococonv, fi%kpts, fi%atoms, fi%sym, fi%cell, enpara, &
444 0 : hub1data, vTot, vCoul, vx, fmpi, results, eig_id, sphhar, stars, fi%vacuum)
445 : END IF
446 0 : IF (fi%input%gw==2) THEN
447 0 : CALL juDFT_end("SPEX data written. Fleur ends.", fmpi%irank)
448 : END IF
449 : END IF
450 :
451 686 : CALL timestart("determination of fermi energy")
452 :
453 686 : IF (fi%noco%l_soc .AND. (.NOT. fi%noco%l_noco)) THEN
454 68 : input_soc%zelec = fi%input%zelec*2
455 68 : IF (fi%hybinp%l_hybrid) &
456 0 : CALL fermie(hybdat%eig_id, fmpi, fi%kpts, fi%input, fi%noco, enpara%epara_min, fi%cell, hybdat%results)
457 68 : CALL fermie(eig_id, fmpi, fi%kpts, input_soc, fi%noco, enpara%epara_min, fi%cell, results)
458 :
459 68 : results%seigv = results%seigv / 2.0
460 68 : results%ts = results%ts / 2.0
461 : ELSE
462 618 : CALL fermie(eig_id, fmpi, fi%kpts, fi%input, fi%noco, enpara%epara_min, fi%cell, results)
463 618 : IF (fi%hybinp%l_hybrid) hybdat%results = results
464 : ENDIF
465 : #ifdef CPP_MPI
466 686 : CALL MPI_BCAST(results%ef, 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
467 2744 : CALL MPI_BCAST(results%w_iks, SIZE(results%w_iks), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
468 686 : if (fi%hybinp%l_hybrid) then
469 186 : CALL MPI_BCAST(hybdat%results%ef, 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
470 744 : CALL MPI_BCAST(hybdat%results%w_iks, SIZE(hybdat%results%w_iks), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
471 : endif
472 : #endif
473 686 : CALL timestop("determination of fermi energy")
474 :
475 : ! TODO: What is commented out here and should it perhaps be removed?
476 : ! !$ !+Wannier
477 : ! !$ IF(wann%l_bs_comf)THEN
478 : ! !$ IF(pc.EQ.1) THEN
479 : ! !$ OPEN(777,file='out_eig.1')
480 : ! !$ OPEN(778,file='out_eig.2')
481 : ! !$ OPEN(779,file='out_eig.1_diag')
482 : ! !$ OPEN(780,file='out_eig.2_diag')
483 : ! !$ END IF
484 : ! !$
485 : ! !$ CALL bs_comfort(eig_id,fi%input,fi%noco,fi%kpts%nkpt,pc)
486 : ! !$
487 : ! !$ IF(pc.EQ.wann%nparampts)THEN
488 : ! !$ CLOSE(777)
489 : ! !$ CLOSE(778)
490 : ! !$ CLOSE(779)
491 : ! !$ CLOSE(780)
492 : ! !$ END IF
493 : ! !$ END IF
494 : ! !$ !-Wannier
495 :
496 : !ENDIF
497 :
498 :
499 686 : IF (forcetheo%eval(eig_id, fi%atoms, fi%kpts, fi%sym, fi%cell, fi%noco, nococonv, input_soc, fmpi, enpara, vToT, results)) THEN
500 : CYCLE forcetheoloop
501 : END IF
502 :
503 686 : CALL timestart("Wannier")
504 686 : IF ((fi%input%l_wann) .AND. (.NOT. wann%l_bs_comf)) THEN
505 : CALL wannier(fmpi, input_soc, fi%kpts, fi%sym, fi%atoms, stars, fi%vacuum, sphhar, &
506 : wann, fi%noco, nococonv, fi%cell, enpara, fi%banddos, fi%sliceplot, vTot, results, &
507 8 : (/eig_id/), (fi%sym%invs) .AND. (.NOT. fi%noco%l_soc) .AND. (.NOT. fi%noco%l_noco), fi%kpts%nkpt)
508 : END IF
509 682 : CALL timestop("Wannier")
510 :
511 : ! Check if the greensFunction have to be calculated
512 682 : IF (fi%gfinp%n > 0) THEN
513 1060 : DO i_gf = 1, fi%gfinp%n
514 : ! Either the set distance has been reached (or is negative)
515 : greensFunction(i_gf)%l_calc = (results%last_distance >= 0.0 .AND. &
516 : results%last_distance < fi%gfinp%minCalcDistance) &
517 : .OR. fi%gfinp%minCalcDistance < 0.0 & !No minCalcDistance distance set
518 1018 : .OR. iter == fi%input%itmax !Maximum iteration reached
519 : ! or we are in the first iteration for Hubbard 1
520 1060 : IF (fi%atoms%n_hia > 0) THEN
521 : greensFunction(i_gf)%l_calc = greensFunction(i_gf)%l_calc .OR. (iter == 1 .AND. (hub1data%iter == 0 &
522 0 : .AND. ALL(ABS(vTot%mmpMat(:, :, fi%atoms%n_u + 1:fi%atoms%n_u + fi%atoms%n_hia, :)) .LT. 1e-12)))
523 : END IF
524 : END DO
525 : END IF
526 :
527 : ! charge density generation
528 682 : CALL timestart("generation of new charge density (total)")
529 682 : CALL outDen%init(stars, fi%atoms, sphhar, fi%vacuum, fi%noco, fi%input%jspins, POTDEN_TYPE_DEN)
530 682 : outDen%iter = inDen%iter
531 : CALL cdngen(eig_id, fmpi, input_soc, fi%banddos, fi%sliceplot, fi%vacuum, &
532 : fi%kpts, fi%atoms, sphhar, stars, fi%sym, fi%juphon, fi%gfinp, fi%hub1inp, &
533 : enpara, fi%cell, fi%noco, nococonv, vTot, results, fi%corespecinput, &
534 682 : archiveType, xcpot, outDen, EnergyDen, greensFunction, hub1data,vxc,exc)
535 : ! The density matrix for DFT+Hubbard1 only changes in hubbard1_setup and is kept constant otherwise
536 1790 : outDen%mmpMat(:, :, fi%atoms%n_u + 1:fi%atoms%n_u + fi%atoms%n_hia, :) = inDen%mmpMat(:, :, fi%atoms%n_u + 1:fi%atoms%n_u + fi%atoms%n_hia, :)
537 :
538 670 : IF (fi%sliceplot%iplot/=0) THEN
539 : ! Plot the full output density.
540 : CALL makeplots(stars, fi%atoms, sphhar, fi%vacuum, fi%input, fmpi, fi%sym, fi%cell, &
541 0 : fi%noco, nococonv, outDen, PLOT_OUTDEN_Y_CORE, fi%sliceplot)
542 :
543 0 : IF ((fi%sliceplot%iplot .NE. 0) .AND. (fmpi%irank .EQ. 0) .AND. (fi%sliceplot%iplot .LT. 64) .AND. (MODULO(fi%sliceplot%iplot, 2) .NE. 1)) THEN
544 0 : CALL juDFT_end("Stopped self consistency loop after plots have been generated.")
545 : END IF
546 : END IF
547 :
548 670 : IF (fi%input%l_rdmft) THEN
549 : SELECT TYPE (xcpot)
550 : TYPE IS (t_xcpot_inbuild)
551 : CALL rdmft(eig_id, fmpi, fi, enpara, stars, &
552 : sphhar, vTot, vCoul, nococonv, xcpot, mpdata, hybdat, &
553 0 : results, archiveType, outDen)
554 : END SELECT
555 : END IF
556 :
557 : #ifdef CPP_MPI
558 670 : CALL MPI_BCAST(enpara%evac, SIZE(enpara%evac), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
559 670 : CALL MPI_BCAST(enpara%evac0, SIZE(enpara%evac0), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
560 2680 : CALL MPI_BCAST(enpara%el0, SIZE(enpara%el0), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
561 2680 : CALL MPI_BCAST(enpara%ello0, SIZE(enpara%ello0), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
562 :
563 : ! TODO: This might be broken at the moment. Fix!
564 : !IF (fi%noco%l_noco) THEN
565 : ! DO n = 1, fi%atoms%ntype
566 : ! IF (fi%noco%l_relax(n)) THEN
567 : ! CALL MPI_BCAST(nococonv%alph(n), 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
568 : ! CALL MPI_BCAST(nococonv%beta(n), 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
569 : ! END IF
570 : ! END DO
571 : ! IF (any(fi%noco%l_constrained)) THEN
572 : ! CALL MPI_BCAST(nococonv%b_con, SIZE(nococonv%b_con), MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
573 : ! END IF
574 : !END IF
575 : #endif
576 670 : CALL timestop("generation of new charge density (total)")
577 :
578 670 : IF (fi%juPhon%l_dfpt) THEN
579 : ! Sideline the actual scf loop for a phonon calculation.
580 : ! It is assumed that the density was converged beforehand.
581 0 : CALL timestop("Iteration")
582 0 : CALL timestart("juPhon DFPT")
583 0 : CALL dfpt(fi, sphhar, stars, nococonv, fi%kpts, fmpi, results, enpara, outDen, vTot, vxc, eig_id, xcpot, hybdat, mpdata, forcetheo)
584 0 : CALL timestop("juPhon DFPT")
585 : END IF
586 :
587 : !CRYSTAL FIELD OUTPUT
588 3010 : IF(ANY(fi%atoms%l_outputCFpot(:)).OR.ANY(fi%atoms%l_outputCFcdn(:))) THEN
589 2 : CALL hub1data%mpi_bc(fmpi%mpi_comm)
590 2 : CALL writeCFOutput(fi,stars,hybdat,sphhar,xcpot,EnergyDen,outDen,hub1data,nococonv,enpara,fmpi)
591 2 : CALL juDFT_end("Crystal Field Output written",fmpi%irank)
592 : END IF
593 :
594 : ! TODO: What is commented out here and should it perhaps be removed?
595 : ! !$ !----> output potential and potential difference
596 : ! !$ IF (disp) THEN
597 : ! !$ reap = .FALSE.
598 : ! !$ CALL timestart("generation of potential (total)")
599 : ! !$ CALL vgen(fi%hybinp,reap,fi%input,xcpot, fi%atoms,sphhar,stars,fi%vacuum,fi%sym,&
600 : ! !$ fi%cell, fi%sliceplot,fmpi, results,fi%noco,outDen,inDenRot,vTot,vx,vCoul)
601 : ! !$ CALL timestop("generation of potential (total)")
602 : ! !$
603 : ! !$ CALL potdis(stars,fi%vacuum,fi%atoms,sphhar, fi%input,fi%cell,fi%sym)
604 : ! !$ END IF
605 :
606 : ! total energy
607 :
608 : ! Rotating from local MT frame in global frame for mixing
609 : ! TODO: Should this be done before the total energy calculation already?
610 668 : CALL toGlobalSpinFrame(fi%noco, nococonv, fi%vacuum, sphhar, stars, fi%sym, fi%cell, fi%input, fi%atoms, inDen, fmpi)
611 668 : CALL toGlobalSpinFrame(fi%noco, nococonv, fi%vacuum, sphhar, stars, fi%sym, fi%cell, fi%input, fi%atoms, outDen, fmpi, .TRUE.)
612 668 : CALL timestart('determination of total energy')
613 : CALL totale(fmpi, fi%atoms, sphhar, stars, fi%vacuum, fi%sym, fi%input, fi%noco, fi%cell, &
614 668 : xcpot, hybdat, vTot, vCoul, iter, inDen, results)
615 664 : CALL timestop('determination of total energy')
616 : END DO forcetheoloop
617 :
618 664 : CALL forcetheo%postprocess()
619 :
620 664 : CALL enpara%mix(fmpi%mpi_comm, fi%atoms, fi%vacuum, fi%input, vTot)
621 664 : field2 = fi%field
622 :
623 : ! mix input and output densities
624 : CALL mix_charge(field2, fmpi, (iter == fi%input%itmax .OR. judft_was_argument("-mix_io")), stars, &
625 : fi%atoms, sphhar, fi%vacuum, fi%input, fi%sym, fi%juphon, fi%cell, fi%noco, nococonv, &
626 800 : archiveType, xcpot, iter, inDen, outDen, results, hub1data%l_runthisiter, fi%sliceplot)
627 :
628 : ! Rotating to the local MT frame
629 : CALL toLocalSpinFrame(fmpi, fi%vacuum, sphhar, stars, fi%sym, fi%cell, fi%noco, &
630 664 : nococonv, fi%input, fi%atoms, .TRUE., inDen, .TRUE.)
631 :
632 664 : IF (fmpi%irank==0) THEN
633 332 : WRITE (oUnit, FMT=8130) iter
634 : 8130 FORMAT(/, 5x, '******* it=', i3, ' is completed********', /,/)
635 332 : call log%add("Interation",int2str(iter))
636 332 : IF (fi%hybinp%l_hybrid) THEN
637 93 : WRITE (*, *) "Iteration:", iter, " Distance:", results%last_distance, " hyb distance:", hybdat%results%last_distance
638 93 : call log%add("Hybrid-distance",float2str(hybdat%results%last_distance))
639 : ELSE
640 239 : WRITE (*, *) "Iteration:", iter, " Distance:", results%last_distance
641 : ENDIF
642 332 : call log%add("Distance",float2str(results%last_distance))
643 332 : call log%report(logmode_info)
644 : END IF ! fmpi%irank.EQ.0
645 664 : CALL timestop("Iteration")
646 :
647 : #ifdef CPP_MPI
648 664 : CALL MPI_BCAST(results%last_distance, 1, MPI_DOUBLE_PRECISION, 0, fmpi%mpi_comm, ierr)
649 664 : CALL MPI_BARRIER(fmpi%mpi_comm, ierr)
650 : #endif
651 :
652 :
653 : l_cont = .TRUE.
654 664 : IF (fi%hybinp%l_hybrid) THEN
655 186 : IF (hybdat%l_calhf) THEN
656 12 : l_cont = l_cont .AND. (iterHF < fi%input%itmax)
657 12 : l_cont = l_cont .AND. (fi%input%mindistance <= results%last_distance)
658 12 : CALL check_time_for_next_iteration(iterHF, l_cont)
659 : ELSE
660 174 : l_cont = l_cont .AND. (iter < 50) ! Security stop for non-converging nested PBE calculations
661 : END IF
662 :
663 186 : IF (hybdat%l_subvxc) THEN
664 98 : results%te_hfex%valence = 0
665 : END IF
666 478 : ELSE IF (fi%atoms%n_hia > 0) THEN
667 0 : l_cont = l_cont .AND. (iter < fi%input%itmax) !The SCF cycle reached the maximum iteration
668 0 : l_cont = l_cont .AND. ((fi%input%mindistance <= results%last_distance) .OR. fi%input%l_f)
669 : !If we have converged run hia if the density matrix has not converged
670 : hub1data%l_runthisiter = .NOT. l_cont .AND. (fi%hub1inp%minoccDistance <= results%last_occdistance &
671 : .OR. results%last_occdistance <= 0.0 .OR. results%last_mmpMatdistance <= 0.0 &
672 0 : .OR. fi%hub1inp%minmatDistance <= results%last_mmpMatdistance)
673 : !Run after first overall iteration to generate a starting density matrix
674 : hub1data%l_runthisiter = hub1data%l_runthisiter .OR. (iter == 1 .AND. (hub1data%iter == 0 &
675 0 : .AND. ALL(ABS(vTot%mmpMat(:, :, fi%atoms%n_u + 1:fi%atoms%n_u + fi%atoms%n_hia, :)) .LT. 1e-12)))
676 0 : hub1data%l_runthisiter = hub1data%l_runthisiter .AND. (iter < fi%input%itmax)
677 0 : hub1data%l_runthisiter = hub1data%l_runthisiter .AND. (hub1data%iter < fi%hub1inp%itmax)
678 : !Prevent that the scf loop terminates
679 0 : l_cont = l_cont .OR. hub1data%l_runthisiter
680 0 : CALL check_time_for_next_iteration(hub1data%overallIteration, l_cont)
681 : ELSE
682 478 : l_cont = l_cont .AND. (iter < fi%input%itmax)
683 : ! MetaGGAs need a at least 2 iterations
684 : l_cont = l_cont .AND. ((fi%input%mindistance <= results%last_distance) .OR. fi%input%l_f &
685 478 : .OR. (xcpot%exc_is_MetaGGA() .and. iter == 1))
686 478 : CALL check_time_for_next_iteration(iter, l_cont)
687 : END IF
688 :
689 : ! Add extra iteration for force theorem if necessary
690 664 : l_forceTheorem = .FALSE.
691 : SELECT TYPE(forcetheo)
692 : TYPE IS(t_forcetheo_mae)
693 : l_forceTheorem = .TRUE.
694 : TYPE IS(t_forcetheo_dmi)
695 : l_forceTheorem = .TRUE.
696 : TYPE IS(t_forcetheo_jij)
697 : l_forceTheorem = .TRUE.
698 : TYPE IS(t_forcetheo_ssdisp)
699 : l_forceTheorem = .TRUE.
700 : END SELECT
701 :
702 0 : IF(l_forceTheorem.AND..NOT.l_cont) THEN
703 0 : IF(.NOT.l_lastIter) THEN
704 0 : l_lastIter = .TRUE.
705 0 : l_cont = .TRUE.
706 : END IF
707 0 : ELSE IF(l_forceTheorem.AND.l_lastIter) THEN
708 0 : l_cont = .FALSE.
709 : END IF
710 :
711 : ! IF file JUDFT_NO_MORE_ITERATIONS is present in the directory, don't do any more iterations
712 664 : IF(fmpi%irank.EQ.0) THEN
713 332 : l_exist = .FALSE.
714 332 : INQUIRE (file='JUDFT_NO_MORE_ITERATIONS', exist=l_exist)
715 332 : IF (l_exist) l_cont = .FALSE.
716 : END IF
717 : #ifdef CPP_MPI
718 664 : CALL MPI_BCAST(l_cont,1,MPI_LOGICAL,0,fmpi%mpi_comm,ierr)
719 : #endif
720 :
721 : ! TODO: What is commented out here and should it perhaps be removed?
722 : !CALL writeTimesXML()
723 :
724 664 : IF (fmpi%irank .EQ. 0) THEN
725 332 : IF (isCurrentXMLElement("iteration")) CALL closeXMLElement('iteration')
726 : END IF
727 :
728 664 : IF ((fi%sliceplot%iplot .NE. 0)) THEN
729 : ! Plot the mixed density
730 : CALL makeplots(stars, fi%atoms, sphhar, fi%vacuum, fi%input, fmpi, fi%sym, &
731 0 : fi%cell, fi%noco, nococonv, inDen, PLOT_MIXDEN_Y_CORE, fi%sliceplot)
732 : ! Plot the mixed valence density
733 : !CALL makeplots(fi%sym,stars,fi%vacuum,fi%atoms,sphhar,fi%input,fi%cell, fi%noco,fi%sliceplot,inDen,PLOT_MIXDEN_N_CORE)
734 : !CALL makeplots(stars, fi%atoms, sphhar, fi%vacuum, fi%input, fi%sym, fi%cell, fi%noco, inDen, PLOT_OUTDEN_N_CORE, fi%sliceplot)
735 : END IF
736 :
737 : ! Break SCF loop if Plots were generated in ongoing run (iplot/=0). This needs to happen here, as the mixed density
738 : ! is the last plottable t_potden to appear in the scf loop and with no mixed density written out (so it is quasi
739 : ! post-process).
740 :
741 2014 : IF ((fi%sliceplot%iplot .NE. 0) .AND. (fmpi%irank .EQ. 0)) THEN
742 0 : CALL juDFT_end("Stopped self consistency loop after plots have been generated.")
743 : END IF
744 :
745 : END DO scfloop ! DO WHILE (l_cont)
746 :
747 132 : CALL add_usage_data("Iterations", iter)
748 :
749 132 : IF (fmpi%irank .EQ. 0) CALL closeXMLElement('scfLoop')
750 :
751 132 : CALL close_eig(eig_id)
752 132 : IF (fi%noco%l_soc .AND. fi%hybinp%l_hybrid) CALL close_eig(hybdat%eig_id)
753 132 : CALL juDFT_end("all done", fmpi%irank)
754 480 : END SUBROUTINE fleur_execute
755 850 : END MODULE m_fleur
|
