fleur.F90 21.6 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 8 9 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 39 40 41 42
  IMPLICIT NONE
CONTAINS
  SUBROUTINE fleur_execute(mpi_comm)

    !     ***************************************************************
    !
    !     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

    USE m_types
    USE m_constants
    USE m_fleur_init
    USE m_optional
43
    USE m_cdn_io
44
    USE m_mixing_history
45
    USE m_qfix
46 47 48 49 50 51 52 53 54 55 56 57 58
    USE m_vgen
    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
59
    USE m_rdmft
Gregor Michalicek's avatar
Gregor Michalicek committed
60
    USE m_io_hybrid
61 62 63 64 65
    USE m_wann_optional
    USE m_wannier
    USE m_bs_comfort
    USE m_dwigner
    USE m_ylm
Matthias Redies's avatar
Matthias Redies committed
66
    USE m_metagga
67
#ifdef CPP_MPI
68
    USE m_mpi_bc_potden
69
#endif
70
    USE m_eig66_io
71
    USE m_chase_diag
72
    USE m_writeBasis
Matthias Redies's avatar
Matthias Redies committed
73
    !$ USE omp_lib
74 75
    IMPLICIT NONE

76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
    INTEGER, INTENT(IN)             :: mpi_comm

    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
    TYPE(t_oneD)                    :: oneD
    TYPE(t_mpi)                     :: mpi
    TYPE(t_coreSpecInput)           :: coreSpecInput
    TYPE(t_wann)                    :: wann
    TYPE(t_potden)                  :: vTot, vx, vCoul, vTemp
Matthias Redies's avatar
Matthias Redies committed
100
    TYPE(t_potden)                  :: inDen, outDen, EnergyDen
101 102
    CLASS(t_xcpot),     ALLOCATABLE :: xcpot
    CLASS(t_forcetheo), ALLOCATABLE :: forcetheo
103
    TYPE(t_greensf)                 :: gOnsite
104
    TYPE(t_greensf)                 :: gIntersite
105
    TYPE(t_hub1ham)                 :: hub1
106 107

    ! local scalars
Matthias Redies's avatar
Matthias Redies committed
108
    INTEGER :: eig_id,archiveType, num_threads
109
    INTEGER :: iter,iterHF
Henning Janssen's avatar
Henning Janssen committed
110
    LOGICAL :: l_opti,l_cont,l_qfix,l_real
111
    REAL    :: fix
112
#ifdef CPP_MPI
113
    INCLUDE 'mpif.h'
114
    INTEGER :: ierr(2),n
115
#endif
Gregor Michalicek's avatar
Gregor Michalicek committed
116

117
    mpi%mpi_comm = mpi_comm
118

119
    CALL timestart("Initialization")
120
    CALL fleur_init(mpi,input,field,DIMENSION,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,sliceplot,&
121
                    banddos,obsolete,enpara,xcpot,results,kpts,hybrid,oneD,coreSpecInput,wann,hub1,l_opti)
122
    CALL timestop("Initialization")
Daniel Wortmann's avatar
Daniel Wortmann committed
123

124 125
    IF ( ( input%preconditioning_param /= 0 ) .AND. oneD%odi%d1 ) THEN
      CALL juDFT_error('Currently no preconditioner for 1D calculations', calledby = 'fleur')
126
    END IF
127

128 129
    IF (l_opti) CALL optional(mpi,atoms,sphhar,vacuum,dimension,&
                              stars,input,sym,cell,sliceplot,obsolete,xcpot,noco,oneD)
130 131 132

    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')
133
       CALL wann_optional(input,kpts,atoms,sym,cell,oneD,noco,wann)
134
    END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
135
  
136

137 138
    iter     = 0
    iterHF   = 0
139
    hub1%iter  = 0
Henning Janssen's avatar
Henning Janssen committed
140
    hub1%l_runthisiter = .FALSE.
141
    l_cont = (iter < input%itmax)
142
    
143 144
    IF (mpi%irank.EQ.0) CALL openXMLElementNoAttributes('scfLoop')

145
    ! Initialize and load inDen density (start)
146
    CALL inDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
147 148
    archiveType = CDN_ARCHIVE_TYPE_CDN1_const
    IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
149 150
    IF(mpi%irank.EQ.0) THEN
       CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
151
                        0,results%ef,l_qfix,inDen)
152
       CALL timestart("Qfix")
153
       CALL qfix(mpi,stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false.,fix)
154 155
       CALL timestop("Qfix")
       CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
156
                         0,-1.0,results%ef,.FALSE.,inDen)
157
    END IF
158
    ! Initialize and load inDen density (end)
159

160
    ! Initialize potentials (start)
161 162 163 164
    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)
165
    ! Initialize potentials (end)
166

167
    ! Initialize Green's function (start)
168
    CALL gOnsite%init(input,lmaxU_const,atoms,noco)
169
    ! Initialize Green's function (end)
170

171
    ! Open/allocate eigenvector storage (start)
172
    l_real=sym%invs.AND..NOT.noco%l_noco.AND..NOT.(noco%l_soc.AND.atoms%n_u+atoms%n_hia>0)
173
    eig_id=open_eig(mpi%mpi_comm,DIMENSION%nbasfcn,DIMENSION%neigd,kpts%nkpt,input%jspins,&
Frank Freimuth's avatar
Frank Freimuth committed
174
                    noco%l_noco,.NOT.INPUT%eig66(1),l_real,noco%l_soc,INPUT%eig66(1),mpi%n_size)
175 176

#ifdef CPP_CHASE
177
    CALL init_chase(mpi,dimension,input,atoms,kpts,noco,.AND..NOT.(noco%l_soc.AND.atoms%n_u+atoms%n_hia>0))
178
#endif
179
    ! Open/allocate eigenvector storage (end)
180

181
    scfloop:DO WHILE (l_cont)
182

183
       iter = iter + 1
184
       IF(hub1%l_runthisiter) hub1%iter = hub1%iter + 1
185
       IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber      '/),&
186
                                                       (/iter,inden%iter/), RESHAPE((/19,13,5,5/),(/2,2/)))
Gregor Michalicek's avatar
Gregor Michalicek committed
187

188 189 190
!!$       !+t3e
!!$       IF (input%alpha.LT.10.0) THEN
!!$
191
!!$          IF (iter.GT.1) THEN
192 193 194
!!$             input%alpha = input%alpha - NINT(input%alpha)
!!$          END IF

195
       !CALL resetIterationDependentTimers()
196 197
       CALL timestart("Iteration")
       IF (mpi%irank.EQ.0) THEN
198 199
          WRITE (6,FMT=8100) iter
8100      FORMAT (/,10x,'   iter=  ',i5)
200
       ENDIF !mpi%irank.eq.0
201
       input%total = .TRUE.
202

Matthias Redies's avatar
Matthias Redies committed
203
#ifdef CPP_CHASE
204
       CALL chase_distance(results%last_distance)
Matthias Redies's avatar
Matthias Redies committed
205
#endif
206

207
#ifdef CPP_MPI
208
       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,inDen)
209 210
#endif

211 212
       dimension%neigd2 = dimension%neigd
       IF (noco%l_soc) dimension%neigd2 = dimension%neigd*2
213 214

       !HF
Matthias Redies's avatar
Matthias Redies committed
215
       !$ num_threads = omp_get_max_threads()
Matthias Redies's avatar
Matthias Redies committed
216
       !$ call omp_set_num_threads(1)
217 218 219
       IF (hybrid%l_hybrid) THEN
          SELECT TYPE(xcpot)
          TYPE IS(t_xcpot_inbuild)
220 221
             CALL calc_hybrid(eig_id,hybrid,kpts,atoms,input,DIMENSION,mpi,noco,&
                              cell,oneD,enpara,results,sym,xcpot,vTot,iter,iterHF)
222
          END SELECT
223
          IF(hybrid%l_calhf) THEN
224
             call mixing_history_reset(mpi)
225 226
             iter = 0
          END IF
227
       ENDIF
Gregor Michalicek's avatar
Gregor Michalicek committed
228 229 230 231
       !RDMFT
       IF(input%l_rdmft) THEN
          CALL open_hybrid_io1(DIMENSION,sym%invs)
       END IF
Frank Freimuth's avatar
Frank Freimuth committed
232 233 234
       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
Matthias Redies's avatar
Matthias Redies committed
235
       !$ call omp_set_num_threads(num_threads)
236

237 238 239 240 241 242 243 244 245 246 247 248 249 250
       !#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

       CALL timestart("generation of potential")
251
       CALL vgen(hybrid,field,input,xcpot,DIMENSION,atoms,sphhar,stars,vacuum,sym,&
Matthias Redies's avatar
Matthias Redies committed
252
                 obsolete,cell,oneD,sliceplot,mpi,results,noco,EnergyDen,inDen,vTot,vx,vCoul)
253 254
       CALL timestop("generation of potential")

255
#ifdef CPP_MPI
256
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
257
#endif
258
       CALL forcetheo%start(vtot,mpi%irank==0)
Daniel Wortmann's avatar
Daniel Wortmann committed
259
       forcetheoloop:DO WHILE(forcetheo%next_job(iter==input%itmax,atoms,noco))
260

Henning Janssen's avatar
Henning Janssen committed
261
          CALL timestart("gen. of hamil. and diag. (total)")
262 263
          CALL timestart("eigen")
          vTemp = vTot
Henning Janssen's avatar
Henning Janssen committed
264 265
          vTemp%mmpMat = 0.0 !To avoid errors later on (When ldaUAdjEnpara is T the density
                             !is carried over after vgen)
266
          CALL timestart("Updating energy parameters")
267
          CALL enpara%update(mpi,atoms,vacuum,input,vToT)
268
          CALL timestop("Updating energy parameters")
Frank Freimuth's avatar
Frank Freimuth committed
269 270
          IF(.not.input%eig66(1))THEN
            CALL eigen(mpi,stars,sphhar,atoms,xcpot,sym,kpts,DIMENSION,vacuum,input,&
271
                     cell,enpara,banddos,noco,oneD,hybrid,iter,eig_id,results,inDen,vTemp,vx,gOnsite,hub1)
Frank Freimuth's avatar
Frank Freimuth committed
272
          ENDIF             
273 274 275
          vTot%mmpMat = vTemp%mmpMat
!!$          eig_idList(pc) = eig_id
          CALL timestop("eigen")
276 277

          ! add all contributions to total energy
278
#ifdef CPP_MPI
279
          ! send all result of local total energies to the r
280 281 282 283 284 285 286 287 288 289 290
          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
291
          END IF
292 293
#endif

294
          ! WRITE(6,fmt='(A)') 'Starting 2nd variation ...'
Frank Freimuth's avatar
Frank Freimuth committed
295
          IF (noco%l_soc.AND..NOT.noco%l_noco.AND..NOT.INPUT%eig66(1)) &
296
             CALL eigenso(eig_id,mpi,DIMENSION,stars,vacuum,atoms,sphhar,&
297
                          obsolete,sym,cell,noco,input,kpts, oneD,vTot,enpara,results,hub1)
Henning Janssen's avatar
Henning Janssen committed
298
          CALL timestop("gen. of hamil. and diag. (total)")
299 300

#ifdef CPP_MPI
301
          CALL MPI_BARRIER(mpi%mpi_comm,ierr)
302 303
#endif

304 305
          ! fermi level and occupancies
          IF (noco%l_soc.AND.(.NOT.noco%l_noco)) DIMENSION%neigd = 2*DIMENSION%neigd
306 307 308 309 310 311 312 313 314 315 316

	  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

317
          !IF ((mpi%irank.EQ.0)) THEN
318 319
             CALL timestart("determination of fermi energy")

320
             IF (noco%l_soc.AND.(.NOT.noco%l_noco)) THEN
321
                input%zelec = input%zelec*2
Henning Janssen's avatar
Henning Janssen committed
322
                CALL fermie(eig_id,mpi,kpts,input,noco,enpara%epara_min,cell,results)
323 324 325 326
                results%seigscv = results%seigscv/2
                results%ts = results%ts/2
                input%zelec = input%zelec/2
             ELSE
Henning Janssen's avatar
Henning Janssen committed
327
                CALL fermie(eig_id,mpi,kpts,input,noco,enpara%epara_min,cell,results)
328 329
             ENDIF
             CALL timestop("determination of fermi energy")
330

331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349
!!$          !+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
350

351
          !ENDIF
352
#ifdef CPP_MPI
353 354
          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)
355 356
#endif

357
          IF (forcetheo%eval(eig_id,DIMENSION,atoms,kpts,sym,cell,noco,input,mpi,oneD,enpara,vToT,results)) THEN
358
             IF (noco%l_soc.AND.(.NOT.noco%l_noco)) DIMENSION%neigd=DIMENSION%neigd/2
359 360
             CYCLE forcetheoloop
          ENDIF
361 362

          
Daniel Wortmann's avatar
Daniel Wortmann committed
363 364 365 366 367 368 369
          !+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
370

371
          ! charge density generation
372
          CALL timestart("generation of new charge density (total)")
373
          CALL outDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
374
          outDen%iter = inDen%iter
375 376 377
          CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum, &
                      dimension,kpts,atoms,sphhar,stars,sym,&
                      enpara,cell,noco,vTot,results,oneD,coreSpecInput,&
Henning Janssen's avatar
Henning Janssen committed
378
                      archiveType,xcpot,outDen,EnergyDen,gOnsite,hub1)
Henning Janssen's avatar
Henning Janssen committed
379
          outDen%mmpMat(:,:,atoms%n_u+1:atoms%n_u+atoms%n_hia,:) = inDen%mmpMat(:,:,atoms%n_u+1:atoms%n_u+atoms%n_hia,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
380 381 382 383
          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,&
Gregor Michalicek's avatar
Gregor Michalicek committed
384
                              sphhar,sym,field,vTot,vCoul,oneD,noco,xcpot,hybrid,results,coreSpecInput,archiveType,outDen)
Gregor Michalicek's avatar
Gregor Michalicek committed
385 386
             END SELECT
          END IF
387

388 389
          IF (noco%l_soc.AND.(.NOT.noco%l_noco)) DIMENSION%neigd=DIMENSION%neigd/2

390
#ifdef CPP_MPI
391
          CALL MPI_BCAST(enpara%evac,SIZE(enpara%evac),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
392 393 394
          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)
395

396 397 398 399 400
          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)
401
                ENDIF
402 403 404
             ENDDO
             IF (noco%l_constr) THEN
                CALL MPI_BCAST(noco%b_con,SIZE(noco%b_con),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
405
             ENDIF
406
          ENDIF
407
#endif
408
          CALL timestop("generation of new charge density (total)")
409

410 411 412 413 414 415 416 417 418 419 420 421 422
             
!!$             !----> 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,&
!!$                     obsolete,cell,oneD,sliceplot,mpi, results,noco,outDen,inDenRot,vTot,vx,vCoul)
!!$                CALL timestop("generation of potential (total)")
!!$
!!$                CALL potdis(stars,vacuum,atoms,sphhar, input,cell,sym)
!!$             END IF
             
423
             ! total energy
424
             CALL timestart('determination of total energy')
Daniel Wortmann's avatar
Daniel Wortmann committed
425
             CALL totale(mpi,atoms,sphhar,stars,vacuum,DIMENSION,sym,input,noco,cell,oneD,&
426
                         xcpot,hybrid,vTot,vCoul,iter,inDen,results)
427
             CALL timestop('determination of total energy')
428
          IF (hybrid%l_hybrid) CALL close_eig(eig_id)
429 430 431 432

       END DO forcetheoloop

       CALL forcetheo%postprocess()
433

434
       CALL enpara%mix(mpi,atoms,vacuum,input,vTot%mt(:,0,:,:),vtot%vacz)
435
       field2 = field
436
       ! mix input and output densities
437 438
       CALL mix_charge(field2,DIMENSION,mpi,(iter==input%itmax.OR.judft_was_argument("-mix_io")),&
            stars,atoms,sphhar,vacuum,input,&
Henning Janssen's avatar
Henning Janssen committed
439
            sym,cell,noco,oneD,archiveType,xcpot,iter,inDen,outDen,results,hub1%l_runthisiter)
440
       
441
       IF(mpi%irank == 0) THEN
442
         WRITE (6,FMT=8130) iter
443
8130     FORMAT (/,5x,'******* it=',i3,'  is completed********',/,/)
444
         WRITE(*,*) "Iteration:",iter," Distance:",results%last_distance
445 446 447 448 449 450 451 452
         !Write out information if a hubbard 1 Iteration was performed
         IF(hub1%l_runthisiter)  THEN
            WRITE(*,*) "Hubbard 1 Iteration: ", hub1%iter," Distance: ", results%last_mmpMatdistance
            WRITE(6,*) "nmmp occupation distance: ", results%last_occdistance
            WRITE(6,*) "nmmp element distance: ", results%last_mmpMatdistance
            WRITE(6,FMT=8140) hub1%iter
8140        FORMAT (/,5x,'******* Hubbard 1 it=',i3,'  is completed********',/,/)
         ENDIF
453
         CALL timestop("Iteration")
454
       END IF ! mpi%irank.EQ.0
455
          
456
#ifdef CPP_MPI
457
       CALL MPI_BCAST(results%last_distance,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
458 459
       CALL MPI_BCAST(results%last_occdistance,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
       CALL MPI_BCAST(results%last_mmpMatdistance,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
460
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
461
#endif
462 463
       CALL priv_geo_end(mpi)

464
       l_cont = .TRUE.
Gregor Michalicek's avatar
Gregor Michalicek committed
465 466 467 468 469 470 471 472
       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
473 474 475
          IF (hybrid%l_subvxc) THEN
             results%te_hfex%valence = 0
          END IF
476 477
       ELSE
          l_cont = l_cont.AND.(iter < input%itmax)
Henning Janssen's avatar
Henning Janssen committed
478
         ! MetaGGAs need a at least 2 iterations
479 480
          l_cont = l_cont.AND.((input%mindistance<=results%last_distance).OR.input%l_f & 
                               .OR. (xcpot%exc_is_MetaGGA() .and. iter == 1))
Henning Janssen's avatar
Henning Janssen committed
481
          !If we have converged run hia if the density matrix has not converged 
Henning Janssen's avatar
Henning Janssen committed
482
          IF(atoms%n_hia>0) THEN
Henning Janssen's avatar
Henning Janssen committed
483
             hub1%l_runthisiter = .NOT.l_cont.AND.(iter < input%itmax)&
484 485
                                  .AND.(input%minoccDistance<=results%last_occdistance&
                                  .OR.input%minmatDistance<=results%last_mmpMatdistance)
Henning Janssen's avatar
Henning Janssen committed
486
             !Run after first overall iteration to generate a starting density matrix
Henning Janssen's avatar
Fixes  
Henning Janssen committed
487
             hub1%l_runthisiter = hub1%l_runthisiter.OR.(iter==5.AND.(hub1%iter == 0&
Henning Janssen's avatar
Henning Janssen committed
488
                                  .AND.ALL(vTot%mmpMat(:,:,atoms%n_u+1:atoms%n_u+atoms%n_hia,:).EQ.0.0)))
Henning Janssen's avatar
Henning Janssen committed
489 490 491
             !Prevent that the scf loop terminates
             l_cont = l_cont.OR.hub1%l_runthisiter
          ENDIF
492
          CALL check_time_for_next_iteration(iter,l_cont)
493
       END IF
494

495
       !CALL writeTimesXML()
496

Henning Janssen's avatar
Henning Janssen committed
497 498
       IF (mpi%irank.EQ.0) THEN
          IF (isCurrentXMLElement("iteration")) CALL closeXMLElement('iteration')
499 500
       END IF

501
    END DO scfloop ! DO WHILE (l_cont)
502

503
    CALL add_usage_data("Iterations",iter)
504 505

    IF (mpi%irank.EQ.0) CALL closeXMLElement('scfLoop')
506 507

    CALL close_eig(eig_id)
508
    CALL juDFT_end("all done",mpi%irank)
509
    
510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531
  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
            OPEN(2,file='qfix')
            WRITE(2,*)"F"
            CLOSE(2)
            PRINT *,"qfix set to F"
         ENDIF
532
         call mixing_history_reset(mpi)
533 534 535
      ENDIF
      CALL juDFT_end(" GEO new inp.xml created ! ",mpi%irank)
    END SUBROUTINE priv_geo_end
536
    
537 538
  END SUBROUTINE fleur_execute
END MODULE m_fleur