fleur.F90 22.8 KB
Newer Older
1
2
3
4
5
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
6
MODULE m_fleur
7
   IMPLICIT NONE
8
CONTAINS
9
   SUBROUTINE fleur_execute(mpi_comm)
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

    !     ***************************************************************
    !
    !     based on flapw7 (c.l.fu, m.weinert, e.wimmer):
    !     full potential linearized augmented plane wave method for thin
    !     films and superlattices (version 7 ---- general symmetry)
    !     symmetry part       ---  e.wimmer
    !     potential generator ---  c.l.fu,r.podloucky
    !     matrix elements     ---  m.weinert
    !     charge density      ---  c.l.fu
    !                                c.l.fu        1987
    !     2nd variation diagon.  --- r.-q. wu      1992
    !     forces a la Yu et al   --- r.podloucky   1995
    !     full relativistic core --- a.shick       1996
    !     broyden mixing         --- r.pentcheva   1996
    !     gga (pw91, pbe)        --- t.asada       1997
    !     local orbitals         --- p.kurz        1997
    !     automatic symmetry     --- w.hofer       1997
    !     core tails & start     --- r.abt         1998
    !     spin orbit coupling    --- a.shick,x.nie 1998
    !     non-colinear magnet.   --- p.kurz        1999
    !     one-dimensional        --- y.mokrousov   2002
    !     exchange parameters    --- m.lezaic      2004
    !
    !                       g.bihlmayer, s.bluegel 1999
    !     ***************************************************************
    !----------------------------------------
    ! this routine is the main PROGRAM

39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
   USE m_types
   USE m_constants
   USE m_fleur_init
   USE m_optional
   USE m_cdn_io
   USE m_mixing_history
   USE m_qfix
   USE m_vgen
   USE m_vgen_coulomb
   USE m_writexcstuff
   USE m_vmatgen
   USE m_eigen
   USE m_eigenso
   USE m_fermie
   USE m_cdngen
   USE m_totale
   USE m_potdis
   USE m_mix
   USE m_xmlOutput
   USE m_juDFT_time
   USE m_calc_hybrid
   USE m_rdmft
   USE m_io_hybrid
   USE m_wann_optional
   USE m_wannier
   USE m_bs_comfort
   USE m_dwigner
   USE m_ylm
   USE m_metagga
   USE m_plot
Robin Hilgers's avatar
Robin Hilgers committed
69

Robin Hilgers's avatar
Robin Hilgers committed
70
71


72
#ifdef CPP_MPI
73
   USE m_mpi_bc_potden
74
#endif
75
76
77
   USE m_eig66_io
   USE m_chase_diag
   USE m_writeBasis
78

Robin Hilgers's avatar
Robin Hilgers committed
79
   USE m_alignSpinAxisMagn
80
81
   !$ USE omp_lib
   IMPLICIT NONE
82

83
   INTEGER, INTENT(IN)             :: mpi_comm
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101

    TYPE(t_input)                   :: input
    TYPE(t_field)                   :: field, field2
    TYPE(t_dimension)               :: DIMENSION
    TYPE(t_atoms)                   :: atoms
    TYPE(t_sphhar)                  :: sphhar
    TYPE(t_cell)                    :: cell
    TYPE(t_stars)                   :: stars
    TYPE(t_sym)                     :: sym
    TYPE(t_noco)                    :: noco
    TYPE(t_vacuum)                  :: vacuum
    TYPE(t_sliceplot)               :: sliceplot
    TYPE(t_banddos)                 :: banddos
    TYPE(t_obsolete)                :: obsolete
    TYPE(t_enpara)                  :: enpara
    TYPE(t_results)                 :: results
    TYPE(t_kpts)                    :: kpts
    TYPE(t_hybrid)                  :: hybrid
Matthias Redies's avatar
Matthias Redies committed
102
    TYPE(t_mpbasis)                 :: mpbasis
103
104
105
106
    TYPE(t_oneD)                    :: oneD
    TYPE(t_mpi)                     :: mpi
    TYPE(t_coreSpecInput)           :: coreSpecInput
    TYPE(t_wann)                    :: wann
107
    TYPE(t_potden)                  :: vTot, vx, vCoul, vTemp, vxcForPlotting
Matthias Redies's avatar
Matthias Redies committed
108
    TYPE(t_potden)                  :: inDen, outDen, EnergyDen
109

110
111
112
113
    CLASS(t_xcpot),     ALLOCATABLE :: xcpot
    CLASS(t_forcetheo), ALLOCATABLE :: forcetheo

    ! local scalars
Matthias Redies's avatar
Matthias Redies committed
114
    INTEGER :: eig_id,archiveType, num_threads
115
    INTEGER :: iter,iterHF,i
Frank Freimuth's avatar
Frank Freimuth committed
116
    INTEGER :: wannierspin
Daniel Wortmann's avatar
Daniel Wortmann committed
117
    LOGICAL :: l_opti,l_cont,l_qfix,l_real
118
    REAL    :: fix
119

120
#ifdef CPP_MPI
121
    INCLUDE 'mpif.h'
122
    INTEGER :: ierr(2),n
123
#endif
124
125
    REAL, ALLOCATABLE :: flh(:,:),flh2(:,:)
    COMPLEX, ALLOCATABLE :: flm(:,:) 
Gregor Michalicek's avatar
Gregor Michalicek committed
126

127
    mpi%mpi_comm = mpi_comm
128

129
    CALL timestart("Initialization")
130
    CALL fleur_init(mpi,input,field,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,sliceplot,&
131
                    banddos,obsolete,enpara,xcpot,results,kpts,mpbasis,hybrid,oneD,coreSpecInput,wann,l_opti)
132
    CALL timestop("Initialization")
133
    
134
135
    IF ( ( input%preconditioning_param /= 0 ) .AND. oneD%odi%d1 ) THEN
      CALL juDFT_error('Currently no preconditioner for 1D calculations', calledby = 'fleur')
136
    END IF
137

138
139
    IF (l_opti) CALL optional(mpi,atoms,sphhar,vacuum,dimension,&
                              stars,input,sym,cell,sliceplot,obsolete,xcpot,noco,oneD)
140
141
142

    IF (input%l_wann.AND.(mpi%irank==0).AND.(.NOT.wann%l_bs_comf)) THEN
       IF(mpi%isize.NE.1) CALL juDFT_error('No Wannier+MPI at the moment',calledby = 'fleur')
143
       CALL wann_optional(input,kpts,atoms,sym,cell,oneD,noco,wann)
144
    END IF
145

146
147
148
    iter     = 0
    iterHF   = 0
    l_cont = (iter < input%itmax)
149

150
151
    IF (mpi%irank.EQ.0) CALL openXMLElementNoAttributes('scfLoop')

152
    ! Initialize and load inDen density (start)
153
    CALL inDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
154
    
155
156
    archiveType = CDN_ARCHIVE_TYPE_CDN1_const
    IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
157
    IF(mpi%irank.EQ.0) THEN
158
       CALL readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
159
                        0,results%ef,l_qfix,inDen)
160
       CALL timestart("Qfix")
161
       CALL qfix(mpi,stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false.,fix)
162
       CALL timestop("Qfix")
163
       CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
164
                         0,-1.0,results%ef,.FALSE.,inDen)
165
    END IF
166
    
167
    IF ((sliceplot%iplot.NE.0 ).AND.(mpi%irank==0) ) THEN
168
169
       CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
                      noco, inDen, PLOT_INPDEN, sliceplot) 
170
    END IF 
171

172
    ! Initialize and load inDen density (end)
173

174
    ! Initialize potentials (start)
175
176
177
178
    CALL vTot%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTTOT)
    CALL vCoul%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTCOUL)
    CALL vx%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTCOUL)
    CALL vTemp%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_POTTOT)
179
    ! Initialize potentials (end)
180

181
    ! Open/allocate eigenvector storage (start)
182
    l_real=sym%invs.AND..NOT.noco%l_noco
Frank Freimuth's avatar
Frank Freimuth committed
183
184
185
186
187
188
189
190
191
192
    if(noco%l_soc.and.input%l_wann)then
    !! Weed up and down spinor components for SOC MLWFs.
    !! When jspins=1 Fleur usually writes only the up-spinor into the eig-file.
    !! Make sure we always get up and down spinors when SOC=true.
       wannierspin=2
    else
       wannierspin = input%jspins       
    endif
    
    eig_id=open_eig(mpi%mpi_comm,DIMENSION%nbasfcn,DIMENSION%neigd,kpts%nkpt,wannierspin,&
Frank Freimuth's avatar
Frank Freimuth committed
193
                    noco%l_noco,.NOT.INPUT%eig66(1),l_real,noco%l_soc,INPUT%eig66(1),mpi%n_size)
194
195

#ifdef CPP_CHASE
Daniel Wortmann's avatar
Daniel Wortmann committed
196
    CALL init_chase(mpi,dimension,input,atoms,kpts,noco,sym%invs.AND..NOT.noco%l_noco)
197
#endif
198
    ! Open/allocate eigenvector storage (end)
199

200
    scfloop:DO WHILE (l_cont)
201

202
       iter = iter + 1
203
       IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber      '/),&
204
                                                       (/iter,inden%iter/), RESHAPE((/19,13,5,5/),(/2,2/)))
Gregor Michalicek's avatar
Gregor Michalicek committed
205

206
207
208
!!$       !+t3e
!!$       IF (input%alpha.LT.10.0) THEN
!!$
209
!!$          IF (iter.GT.1) THEN
210
211
212
!!$             input%alpha = input%alpha - NINT(input%alpha)
!!$          END IF

213
       !CALL resetIterationDependentTimers()
214
215
       CALL timestart("Iteration")
       IF (mpi%irank.EQ.0) THEN
216
217
          WRITE (6,FMT=8100) iter
8100      FORMAT (/,10x,'   iter=  ',i5)
218
       ENDIF !mpi%irank.eq.0
219
       input%total = .TRUE.
220

Matthias Redies's avatar
Matthias Redies committed
221
#ifdef CPP_CHASE
222
       CALL chase_distance(results%last_distance)
Matthias Redies's avatar
Matthias Redies committed
223
#endif
224

225
#ifdef CPP_MPI
226
       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,inDen)
227
228
#endif

229
230
       dimension%neigd2 = dimension%neigd
       IF (noco%l_soc) dimension%neigd2 = dimension%neigd*2
231
232

       !HF
233
234
235
       IF (hybrid%l_hybrid) THEN
          SELECT TYPE(xcpot)
          TYPE IS(t_xcpot_inbuild)
Matthias Redies's avatar
Matthias Redies committed
236
             CALL calc_hybrid(eig_id,mpbasis,hybrid,kpts,atoms,input,DIMENSION,mpi,noco,&
Matthias Redies's avatar
Matthias Redies committed
237
                              cell,oneD,enpara,results,sym,xcpot,vTot,iterHF)
238
          END SELECT
239
          IF(hybrid%l_calhf) THEN
240
             call mixing_history_reset(mpi)
241
242
             iter = 0
          END IF
243
       ENDIF
Gregor Michalicek's avatar
Gregor Michalicek committed
244
245
246
247
       !RDMFT
       IF(input%l_rdmft) THEN
          CALL open_hybrid_io1(DIMENSION,sym%invs)
       END IF
248

Frank Freimuth's avatar
Frank Freimuth committed
249
250
251
       IF(.not.input%eig66(1))THEN
          CALL reset_eig(eig_id,noco%l_soc) ! This has to be placed after the calc_hybrid call but before eigen
       END IF
252

253
254
255
256
257
258
259
260
261
262
263
264
265
       !#endif

!!$             DO pc = 1, wann%nparampts
!!$                !---> gwf
!!$                IF (wann%l_sgwf.OR.wann%l_ms) THEN
!!$                   noco%qss(:) = wann%param_vec(:,pc)
!!$                   noco%alph(:) = wann%param_alpha(:,pc)
!!$                ELSE IF (wann%l_socgwf) THEN
!!$                   IF(wann%l_dim(2)) noco%phi   = tpi_const * wann%param_vec(2,pc)
!!$                   IF(wann%l_dim(3)) noco%theta = tpi_const * wann%param_vec(3,pc)
!!$                END IF
       !---< gwf

266
!START Rot For Testing (HIGHLY EXPERIMENTAL ROUTINE)
267
268
    !   IF(.FALSE.)CALL rotateMagnetToSpinAxis(vacuum,sphhar,stars&
    !          ,sym,oneD,cell,noco,input,atoms,inDen)
269
!END Rot For Testing (HIGHLY EXPERIMENTAL ROUTINE)
270
       CALL timestart("generation of potential")
271
       CALL vgen(hybrid,field,input,xcpot,DIMENSION,atoms,sphhar,stars,vacuum,sym,&
272
                 obsolete,cell,oneD,sliceplot,mpi,results,noco,EnergyDen,inDen,vTot,vx,vCoul)
273
274
       CALL timestop("generation of potential")

275
276
277
       IF ((sliceplot%iplot.NE.0 ).AND.(mpi%irank==0) ) THEN            
          CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
                         noco, vTot, PLOT_POT_TOT, sliceplot)      
278
279
280
          !CALL makeplots(sym,stars,vacuum,atoms,sphhar,input,cell,oneD,noco,sliceplot,vCoul,PLOT_POT_COU)
          !CALL subPotDen(vxcForPlotting,vTot,vCoul)
          !CALL makeplots(sym,stars,vacuum,atoms,sphhar,input,cell,oneD,noco,sliceplot,vxcForPlotting,PLOT_POT_VXC)
281
       END IF 
282

283
#ifdef CPP_MPI
284
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
285
286
#endif

287
       CALL forcetheo%start(vtot,mpi%irank==0)
Daniel Wortmann's avatar
Daniel Wortmann committed
288
       forcetheoloop:DO WHILE(forcetheo%next_job(iter==input%itmax,atoms,noco))
289

Matthias Redies's avatar
Matthias Redies committed
290
          CALL timestart("gen. of hamil. and diag. (total)")
291
292
          CALL timestart("eigen")
          vTemp = vTot
293
          CALL timestart("Updating energy parameters")
294
          CALL enpara%update(mpi,atoms,vacuum,input,vToT)
295
          CALL timestop("Updating energy parameters")
Frank Freimuth's avatar
Frank Freimuth committed
296
297
          IF(.not.input%eig66(1))THEN
            CALL eigen(mpi,stars,sphhar,atoms,xcpot,sym,kpts,DIMENSION,vacuum,input,&
298
                     cell,enpara,banddos,noco,oneD,mpbasis,hybrid,iter,eig_id,results,inDen,vTemp,vx)
299
          ENDIF
300
301
302
          vTot%mmpMat = vTemp%mmpMat
!!$          eig_idList(pc) = eig_id
          CALL timestop("eigen")
303
304

          ! add all contributions to total energy
305
#ifdef CPP_MPI
306
          ! send all result of local total energies to the r
307
308
309
310
311
312
313
314
315
316
317
          IF (hybrid%l_hybrid.AND.hybrid%l_calhf) THEN
             IF (mpi%irank==0) THEN
                CALL MPI_Reduce(MPI_IN_PLACE,results%te_hfex%core,1,MPI_REAL8,MPI_SUM,0,mpi%mpi_comm,ierr(1))
             ELSE
                CALL MPI_Reduce(results%te_hfex%core,MPI_IN_PLACE,1,MPI_REAL8,MPI_SUM,0, mpi%mpi_comm,ierr(1))
             END IF
             IF (mpi%irank==0) THEN
                CALL MPI_Reduce(MPI_IN_PLACE,results%te_hfex%valence,1,MPI_REAL8,MPI_SUM,0,mpi%mpi_comm,ierr(1))
             ELSE
                CALL MPI_Reduce(results%te_hfex%valence,MPI_IN_PLACE,1,MPI_REAL8,MPI_SUM,0, mpi%mpi_comm,ierr(1))
             END IF
318
          END IF
319
320
#endif

321
          ! WRITE(6,fmt='(A)') 'Starting 2nd variation ...'
Frank Freimuth's avatar
Frank Freimuth committed
322
          IF (noco%l_soc.AND..NOT.noco%l_noco.AND..NOT.INPUT%eig66(1)) &
323
324
             CALL eigenso(eig_id,mpi,DIMENSION,stars,vacuum,atoms,sphhar,&
                          obsolete,sym,cell,noco,input,kpts, oneD,vTot,enpara,results)
Matthias Redies's avatar
Matthias Redies committed
325
          CALL timestop("gen. of hamil. and diag. (total)")
326
327

#ifdef CPP_MPI
328
          CALL MPI_BARRIER(mpi%mpi_comm,ierr)
329
330
#endif

331
332
          ! fermi level and occupancies
          IF (noco%l_soc.AND.(.NOT.noco%l_noco)) DIMENSION%neigd = 2*DIMENSION%neigd
333
334
335
336
337
338
339
340
341
342
343

	  IF (input%gw.GT.0) THEN
	    IF (mpi%irank.EQ.0) THEN
	       CALL writeBasis(input,noco,kpts,atoms,sym,cell,enpara,vTot,vCoul,vx,mpi,DIMENSION,&
		  	     results,eig_id,oneD,sphhar,stars,vacuum)
	    END IF
	    IF (input%gw.EQ.2) THEN
	       CALL juDFT_end("GW data written. Fleur ends.",mpi%irank)
	    END IF
	  END IF

344
          !IF ((mpi%irank.EQ.0)) THEN
345
346
             CALL timestart("determination of fermi energy")

347
             IF (noco%l_soc.AND.(.NOT.noco%l_noco)) THEN
348
                input%zelec = input%zelec*2
349
                CALL fermie(eig_id,mpi,kpts,input,noco,enpara%epara_min,cell,results)
350
351
352
353
                results%seigscv = results%seigscv/2
                results%ts = results%ts/2
                input%zelec = input%zelec/2
             ELSE
354
                CALL fermie(eig_id,mpi,kpts,input,noco,enpara%epara_min,cell,results)
355
356
             ENDIF
             CALL timestop("determination of fermi energy")
357

358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
!!$          !+Wannier
!!$          IF(wann%l_bs_comf)THEN
!!$             IF(pc.EQ.1) THEN
!!$                OPEN(777,file='out_eig.1')
!!$                OPEN(778,file='out_eig.2')
!!$                OPEN(779,file='out_eig.1_diag')
!!$                OPEN(780,file='out_eig.2_diag')
!!$             END IF
!!$
!!$             CALL bs_comfort(eig_id,DIMENSION,input,noco,kpts%nkpt,pc)
!!$
!!$             IF(pc.EQ.wann%nparampts)THEN
!!$                CLOSE(777)
!!$                CLOSE(778)
!!$                CLOSE(779)
!!$                CLOSE(780)
!!$             END IF
!!$          END IF
!!$          !-Wannier
377

378
          !ENDIF
379
#ifdef CPP_MPI
380
381
          CALL MPI_BCAST(results%ef,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
          CALL MPI_BCAST(results%w_iks,SIZE(results%w_iks),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
382
383
#endif

384
          IF (forcetheo%eval(eig_id,DIMENSION,atoms,kpts,sym,cell,noco,input,mpi,oneD,enpara,vToT,results)) THEN
385
             IF (noco%l_soc.AND.(.NOT.noco%l_noco)) DIMENSION%neigd=DIMENSION%neigd/2
386
387
             CYCLE forcetheoloop
          ENDIF
388

389

Daniel Wortmann's avatar
Daniel Wortmann committed
390
391
392
393
394
395
396
          !+Wannier functions
          IF ((input%l_wann).AND.(.NOT.wann%l_bs_comf)) THEN
             CALL wannier(DIMENSION,mpi,input,kpts,sym,atoms,stars,vacuum,sphhar,oneD,&
                  wann,noco,cell,enpara,banddos,sliceplot,vTot,results,&
                  (/eig_id/),(sym%invs).AND.(.NOT.noco%l_soc).AND.(.NOT.noco%l_noco),kpts%nkpt)
          END IF
          !-Wannier
397

398
          ! charge density generation
399
          CALL timestart("generation of new charge density (total)")
400
          CALL outDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
401
          outDen%iter = inDen%iter
402
403
404
405
          CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum, &
                      dimension,kpts,atoms,sphhar,stars,sym,&
                      enpara,cell,noco,vTot,results,oneD,coreSpecInput,&
                      archiveType,xcpot,outDen,EnergyDen)
406
          IF ((sliceplot%iplot.NE.0 ).AND.(mpi%irank==0) ) THEN        
Robin Hilgers's avatar
Robin Hilgers committed
407
!               CDN including core charge
408
               ! CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, &
Robin Hilgers's avatar
Robin Hilgers committed
409
!                               cell, noco, outDen, PLOT_OUTDEN_Y_CORE, sliceplot)
410
!!               CDN subtracted by core charge
411
               ! CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, &
Robin Hilgers's avatar
Robin Hilgers committed
412
!                               cell, noco, outDen, PLOT_OUTDEN_N_CORE, sliceplot)
413
          END IF 
414

Gregor Michalicek's avatar
Gregor Michalicek committed
415
416
417
418
          IF (input%l_rdmft) THEN
             SELECT TYPE(xcpot)
                TYPE IS(t_xcpot_inbuild)
                   CALL rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars,vacuum,dimension,&
Matthias Redies's avatar
Matthias Redies committed
419
                              sphhar,sym,field,vTot,vCoul,oneD,noco,xcpot,mpbasis,hybrid,results,coreSpecInput,archiveType,outDen)
Gregor Michalicek's avatar
Gregor Michalicek committed
420
421
             END SELECT
          END IF
422

423
424
          IF (noco%l_soc.AND.(.NOT.noco%l_noco)) DIMENSION%neigd=DIMENSION%neigd/2

425
#ifdef CPP_MPI
426
          CALL MPI_BCAST(enpara%evac,SIZE(enpara%evac),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
427
428
429
          CALL MPI_BCAST(enpara%evac0,SIZE(enpara%evac0),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
          CALL MPI_BCAST(enpara%el0,SIZE(enpara%el0),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
          CALL MPI_BCAST(enpara%ello0,SIZE(enpara%ello0),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
430

431
432
433
434
435
          IF (noco%l_noco) THEN
             DO n= 1,atoms%ntype
                IF (noco%l_relax(n)) THEN
                   CALL MPI_BCAST(noco%alph(n),1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
                   CALL MPI_BCAST(noco%beta(n),1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
436
                ENDIF
437
438
439
             ENDDO
             IF (noco%l_constr) THEN
                CALL MPI_BCAST(noco%b_con,SIZE(noco%b_con),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
440
             ENDIF
441
          ENDIF
442
#endif
443
          CALL timestop("generation of new charge density (total)")
444
!START Rot For Testing (HIGHLY EXPERIMENTAL ROUTINE)
445
!IF (.FALSE.) CALL rotateMagnetFromSpinAxis(noco,vacuum,sphhar,stars,sym,oneD,cell,input,atoms,outDen,inDen)
446
!END Rot For Testing (HIGHLY EXPERIMENTAL ROUTINE)
447
448
449
450
451
452
!!$             !----> output potential and potential difference
!!$             IF (obsolete%disp) THEN
!!$                reap = .FALSE.
!!$                input%total = .FALSE.
!!$                CALL timestart("generation of potential (total)")
!!$                CALL vgen(hybrid,reap,input,xcpot,DIMENSION, atoms,sphhar,stars,vacuum,sym,&
453
!!$                     obsolete,cell,oneD,sliceplot,mpi, results,noco,outDen,inDenRot,vTot,vx,vCoul)
454
455
456
457
!!$                CALL timestop("generation of potential (total)")
!!$
!!$                CALL potdis(stars,vacuum,atoms,sphhar, input,cell,sym)
!!$             END IF
458

459
             ! total energy
460
             CALL timestart('determination of total energy')
Daniel Wortmann's avatar
Daniel Wortmann committed
461
             CALL totale(mpi,atoms,sphhar,stars,vacuum,DIMENSION,sym,input,noco,cell,oneD,&
462
                         xcpot,hybrid,vTot,vCoul,iter,inDen,results)
463
             CALL timestop('determination of total energy')
464
          IF (hybrid%l_hybrid) CALL close_eig(eig_id)
465
466
467
468

       END DO forcetheoloop

       CALL forcetheo%postprocess()
469

470
       CALL enpara%mix(mpi,atoms,vacuum,input,vTot%mt(:,0,:,:),vtot%vacz)
471
       field2 = field
472
473

       ! mix input and output densities
474
475
       CALL mix_charge(field2,DIMENSION,mpi,(iter==input%itmax.OR.judft_was_argument("-mix_io")),&
            stars,atoms,sphhar,vacuum,input,&
Matthias Redies's avatar
Matthias Redies committed
476
            sym,cell,noco,oneD,archiveType,xcpot,iter,inDen,outDen,results)
477
!Plots of mixed density       
478
479
       IF ((sliceplot%iplot.NE.0 ).AND.(mpi%irank==0) ) THEN        
!               CDN including core charge
480
481
!                CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, &
!                               cell, noco, outDen, PLOT_MIXDEN_Y_CORE, sliceplot)
482
!!               CDN subtracted by core charge
483
484
485
!                CALL makeplots(sym,stars,vacuum,atoms,sphhar,input,cell,oneD,noco,sliceplot,inDen,PLOT_MIXDEN_N_CORE)
!                CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, &
!                               cell, noco, outDen, PLOT_OUTDEN_N_CORE, sliceplot)
Matthias Redies's avatar
Matthias Redies committed
486
       END IF 
487

488
       IF(mpi%irank == 0) THEN
489
         WRITE (6,FMT=8130) iter
490
8130     FORMAT (/,5x,'******* it=',i3,'  is completed********',/,/)
491
         WRITE(*,"(A,I3,A,ES15.8)") "Iteration: ",iter," Distance: ",results%last_distance
492
         CALL timestop("Iteration")
493
       END IF ! mpi%irank.EQ.0
494

495
#ifdef CPP_MPI
496
497
       CALL MPI_BCAST(results%last_distance,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
498
#endif
499
500
       CALL priv_geo_end(mpi)

501
       l_cont = .TRUE.
Gregor Michalicek's avatar
Gregor Michalicek committed
502
503
504
505
506
507
508
509
       IF (hybrid%l_hybrid) THEN
          IF(hybrid%l_calhf) THEN
             l_cont = l_cont.AND.(iterHF < input%itmax)
             l_cont = l_cont.AND.(input%mindistance<=results%last_distance)
             CALL check_time_for_next_iteration(iterHF,l_cont)
          ELSE
             l_cont = l_cont.AND.(iter < 50) ! Security stop for non-converging nested PBE calculations
          END IF
510
511
512
          IF (hybrid%l_subvxc) THEN
             results%te_hfex%valence = 0
          END IF
513
514
       ELSE
          l_cont = l_cont.AND.(iter < input%itmax)
515
          ! MetaGGAs need a at least 2 iterations
516
          l_cont = l_cont.AND.((input%mindistance<=results%last_distance).OR.input%l_f &
517
                               .OR. (xcpot%exc_is_MetaGGA() .and. iter == 1))
518
          CALL check_time_for_next_iteration(iter,l_cont)
519
       END IF
520

521
       !CALL writeTimesXML()
522

523
524
       IF (mpi%irank.EQ.0) THEN
          IF (isCurrentXMLElement("iteration")) CALL closeXMLElement('iteration')
525
526
       END IF

Robin Hilgers's avatar
Robin Hilgers committed
527
  !Break SCF loop if Plots were generated in ongoing run (iplot=/=0).
528
529
       IF(sliceplot%iplot.NE.0) THEN
          CALL juDFT_end("Stopped self consistency loop after plots have been generated.")
530
531
       END IF

532
    END DO scfloop ! DO WHILE (l_cont)
533

534
    CALL add_usage_data("Iterations",iter)
535
536

    IF (mpi%irank.EQ.0) CALL closeXMLElement('scfLoop')
537
538
539

    CALL close_eig(eig_id)

540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
    CALL juDFT_end("all done",mpi%irank)
  CONTAINS
    SUBROUTINE priv_geo_end(mpi)
      TYPE(t_mpi),INTENT(IN)::mpi
      LOGICAL :: l_exist
      !Check if a new input was generated
      INQUIRE (file='inp_new',exist=l_exist)
      IF (l_exist) THEN
         CALL juDFT_end(" GEO new inp created ! ",mpi%irank)
      END IF
      !check for inp.xml
      INQUIRE (file='inp_new.xml',exist=l_exist)
      IF (.NOT.l_exist) RETURN
      IF (mpi%irank==0) THEN
         CALL system('mv inp.xml inp_old.xml')
         CALL system('mv inp_new.xml inp.xml')
         INQUIRE (file='qfix',exist=l_exist)
         IF (l_exist) THEN
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
            CALL juDFT_end(" GEO new inp created ! ",mpi%irank)
         END IF
         !check for inp.xml
         INQUIRE (file='inp_new.xml',exist=l_exist)
         IF (.NOT.l_exist) RETURN
         IF (mpi%irank==0) THEN
            CALL system('mv inp.xml inp_old.xml')
            CALL system('mv inp_new.xml inp.xml')
            INQUIRE (file='qfix',exist=l_exist)
            IF (l_exist) THEN
               OPEN(2,file='qfix')
               WRITE(2,*)"F"
               CLOSE(2)
               PRINT *,"qfix set to F"
            ENDIF
            CALL mixing_history_reset(mpi)
574
         ENDIF
575
         call mixing_history_reset(mpi)
576
577
578
      ENDIF
      CALL juDFT_end(" GEO new inp.xml created ! ",mpi%irank)
    END SUBROUTINE priv_geo_end
579

580
581
  END SUBROUTINE fleur_execute
END MODULE m_fleur