fleur.F90 24.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
   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
59
   USE m_calc_hybrid
60
   USE m_rdmft
Matthias Redies's avatar
Matthias Redies committed
61
   USE m_io_hybinp
62 63 64 65 66 67 68
   USE m_wann_optional
   USE m_wannier
   USE m_bs_comfort
   USE m_dwigner
   USE m_ylm
   USE m_metagga
   USE m_plot
69
   USE m_hubbard1_setup
Robin Hilgers's avatar
Robin Hilgers committed
70

Robin Hilgers's avatar
Robin Hilgers committed
71 72


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

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

84
   INTEGER, INTENT(IN)             :: mpi_comm
85 86 87

    TYPE(t_input)                   :: input
    TYPE(t_field)                   :: field, field2
88

89 90 91 92 93
    TYPE(t_atoms) ,TARGET           :: atoms
    TYPE(t_sphhar),TARGET           :: sphhar
    TYPE(t_cell) ,TARGET            :: cell
    TYPE(t_stars),TARGET            :: stars
    TYPE(t_sym),TARGET              :: sym
94 95 96 97 98 99 100
    TYPE(t_noco)                    :: noco
    TYPE(t_vacuum)                  :: vacuum
    TYPE(t_sliceplot)               :: sliceplot
    TYPE(t_banddos)                 :: banddos
    TYPE(t_enpara)                  :: enpara
    TYPE(t_results)                 :: results
    TYPE(t_kpts)                    :: kpts
Matthias Redies's avatar
Matthias Redies committed
101
    TYPE(t_mpinp)                   :: mpinp
Matthias Redies's avatar
Matthias Redies committed
102
    TYPE(t_hybinp)                  :: hybinp
103
    TYPE(t_hybdat)                  :: hybdat
104
    TYPE(t_mpdata)                 :: mpdata
105 106 107 108
    TYPE(t_oneD)                    :: oneD
    TYPE(t_mpi)                     :: mpi
    TYPE(t_coreSpecInput)           :: coreSpecInput
    TYPE(t_wann)                    :: wann
109
    TYPE(t_potden)                  :: vTot, vx, vCoul, vTemp, vxcForPlotting
Matthias Redies's avatar
Matthias Redies committed
110
    TYPE(t_potden)                  :: inDen, outDen, EnergyDen
111

112 113
    CLASS(t_xcpot),     ALLOCATABLE :: xcpot
    CLASS(t_forcetheo), ALLOCATABLE :: forcetheo
114
    TYPE(t_greensf)                 :: gOnsite
115
    TYPE(t_greensf)                 :: gIntersite
116
    TYPE(t_hub1ham)                 :: hub1
117 118

    ! local scalars
Matthias Redies's avatar
Matthias Redies committed
119
    INTEGER :: eig_id,archiveType, num_threads
120
    INTEGER :: iter,iterHF,i,n
Frank Freimuth's avatar
Frank Freimuth committed
121
    INTEGER :: wannierspin
Daniel Wortmann's avatar
Daniel Wortmann committed
122
    LOGICAL :: l_opti,l_cont,l_qfix,l_real
123
    REAL    :: fix
124

125
#ifdef CPP_MPI
126
    INCLUDE 'mpif.h'
127
    INTEGER :: ierr(2)
128
#endif
129
    REAL, ALLOCATABLE :: flh(:,:),flh2(:,:)
130
    COMPLEX, ALLOCATABLE :: flm(:,:)
Gregor Michalicek's avatar
Gregor Michalicek committed
131

132
    mpi%mpi_comm = mpi_comm
133

134
    CALL timestart("Initialization")
135
    CALL fleur_init(mpi,input,field,atoms,sphhar,cell,stars,sym,noco,vacuum,forcetheo,sliceplot,&
Matthias Redies's avatar
Matthias Redies committed
136
                    banddos,enpara,xcpot,results,kpts,mpinp,hybinp,oneD,coreSpecInput,hub1,wann)
137
    CALL timestop("Initialization")
Daniel Wortmann's avatar
Daniel Wortmann committed
138

139 140
    IF ( ( input%preconditioning_param /= 0 ) .AND. oneD%odi%d1 ) THEN
      CALL juDFT_error('Currently no preconditioner for 1D calculations', calledby = 'fleur')
141
    END IF
142

143
    CALL optional(mpi,atoms,sphhar,vacuum,&
Daniel Wortmann's avatar
Daniel Wortmann committed
144
                              stars,input,sym,cell,sliceplot,xcpot,noco,oneD)
145 146 147

    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')
148
       CALL wann_optional(input,kpts,atoms,sym,cell,oneD,noco,wann)
149
    END IF
150

151 152
    iter     = 0
    iterHF   = 0
153
    hub1%iter  = 0
Henning Janssen's avatar
Henning Janssen committed
154
    hub1%l_runthisiter = .FALSE.
155
    l_cont = (iter < input%itmax)
156

157 158
    IF (mpi%irank.EQ.0) CALL openXMLElementNoAttributes('scfLoop')

159
    ! Initialize and load inDen density (start)
160
    CALL inDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
161

162 163
    archiveType = CDN_ARCHIVE_TYPE_CDN1_const
    IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const
164
    IF(mpi%irank.EQ.0) THEN
165
       CALL readDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
166
                        0,results%ef,l_qfix,inDen)
167
       CALL timestart("Qfix")
168
       CALL qfix(mpi,stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false.,fix)
169
       CALL timestop("Qfix")
170
       CALL writeDensity(stars,noco,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
171
                         0,-1.0,results%ef,.FALSE.,inDen)
172
    END IF
173

174
    IF ((sliceplot%iplot.NE.0 ).AND.(mpi%irank==0) ) THEN
175
       CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
176 177
                      noco, inDen, PLOT_INPDEN, sliceplot)
    END IF
178

179
    ! Initialize and load inDen density (end)
180

181
    ! Initialize potentials (start)
182 183 184 185
    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)
186
    ! Initialize potentials (end)
187

188
    ! Initialize Green's function (start)
189
    CALL gOnsite%init(input,lmaxU_const,atoms,noco)
190
    ! Initialize Green's function (end)
191

192
    ! Open/allocate eigenvector storage (start)
193
    l_real=sym%invs.AND..NOT.noco%l_noco.AND..NOT.(noco%l_soc.AND.atoms%n_u+atoms%n_hia>0)
Frank Freimuth's avatar
Frank Freimuth committed
194
    if(noco%l_soc.and.input%l_wann)then
195 196 197
    	 !! 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.
Frank Freimuth's avatar
Frank Freimuth committed
198 199
       wannierspin=2
    else
200
       wannierspin = input%jspins
Frank Freimuth's avatar
Frank Freimuth committed
201
    endif
202

203 204
    eig_id=open_eig(mpi%mpi_comm,lapw_dim_nbasfcn,input%neig,kpts%nkpt,wannierspin,&
                    noco%l_noco,.true.,l_real,noco%l_soc,.false.,mpi%n_size)
205 206

#ifdef CPP_CHASE
207
    CALL init_chase(mpi,input,atoms,kpts,noco,l_real)
208
#endif
209
    ! Open/allocate eigenvector storage (end)
210

211
    scfloop:DO WHILE (l_cont)
212

213
       iter = iter + 1
214 215
       IF(hub1%l_runthisiter.AND.atoms%n_hia>0) THEN
          hub1%iter = hub1%iter + 1
216
          CALL hubbard1_setup(atoms,input,sym,mpi,noco,vTot,gOnsite,hub1,results,inDen)
217
       ENDIF
218
       IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber      '/),&
219
                                                       (/iter,inden%iter/), RESHAPE((/19,13,5,5/),(/2,2/)))
Gregor Michalicek's avatar
Gregor Michalicek committed
220

221 222 223
!!$       !+t3e
!!$       IF (input%alpha.LT.10.0) THEN
!!$
224
!!$          IF (iter.GT.1) THEN
225 226 227
!!$             input%alpha = input%alpha - NINT(input%alpha)
!!$          END IF

228
       !CALL resetIterationDependentTimers()
229 230
       CALL timestart("Iteration")
       IF (mpi%irank.EQ.0) THEN
231 232
          WRITE (6,FMT=8100) iter
8100      FORMAT (/,10x,'   iter=  ',i5)
233
       ENDIF !mpi%irank.eq.0
234
       input%total = .TRUE.
235

Matthias Redies's avatar
Matthias Redies committed
236
#ifdef CPP_CHASE
237
       CALL chase_distance(results%last_distance)
Matthias Redies's avatar
Matthias Redies committed
238
#endif
239

240
#ifdef CPP_MPI
241
       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,inDen)
242 243
#endif

244

245
       !HF
Matthias Redies's avatar
Matthias Redies committed
246
       IF (hybinp%l_hybrid) THEN
247 248
          SELECT TYPE(xcpot)
          TYPE IS(t_xcpot_inbuild)
249
             CALL calc_hybrid(eig_id,mpdata,hybinp,hybdat,kpts,atoms,input,mpi,noco,&
Matthias Redies's avatar
Matthias Redies committed
250
                              cell,oneD,enpara,results,sym,xcpot,vTot,iterHF)
251
          END SELECT
252
          IF(hybdat%l_calhf) THEN
253
             call mixing_history_reset(mpi)
254 255
             iter = 0
          END IF
256
       ENDIF
Gregor Michalicek's avatar
Gregor Michalicek committed
257 258
       !RDMFT
       IF(input%l_rdmft) THEN
Matthias Redies's avatar
Matthias Redies committed
259
          CALL open_hybinp_io1(sym%invs)
Gregor Michalicek's avatar
Gregor Michalicek committed
260
       END IF
261

262
       !IF(.not.input%eig66(1))THEN
263
          CALL reset_eig(eig_id,noco%l_soc) ! This has to be placed after the calc_hybrid call but before eigen
264
       !END IF
265

266 267 268 269 270 271 272 273 274 275 276 277 278
       !#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

279
!START Rot For Testing (HIGHLY EXPERIMENTAL ROUTINE)
280 281
    !   IF(.FALSE.)CALL rotateMagnetToSpinAxis(vacuum,sphhar,stars&
    !          ,sym,oneD,cell,noco,input,atoms,inDen)
282
!END Rot For Testing (HIGHLY EXPERIMENTAL ROUTINE)
283
       CALL timestart("generation of potential")
284
       CALL vgen(hybdat,field,input,xcpot,atoms,sphhar,stars,vacuum,sym,&
Daniel Wortmann's avatar
Daniel Wortmann committed
285
                 cell,oneD,sliceplot,mpi,results,noco,EnergyDen,inDen,vTot,vx,vCoul)
286 287
       CALL timestop("generation of potential")

288
       IF ((sliceplot%iplot.NE.0 ).AND.(mpi%irank==0) ) THEN
289
          CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
290
                         noco, vTot, PLOT_POT_TOT, sliceplot)
291 292 293
          !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)
294
       END IF
295

296
#ifdef CPP_MPI
297
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
298
#endif
299
       CALL forcetheo%start(vtot,mpi%irank==0)
Daniel Wortmann's avatar
Daniel Wortmann committed
300
       forcetheoloop:DO WHILE(forcetheo%next_job(iter==input%itmax,atoms,noco))
301

Matthias Redies's avatar
Matthias Redies committed
302
          CALL timestart("gen. of hamil. and diag. (total)")
303 304
          CALL timestart("eigen")
          vTemp = vTot
Henning Janssen's avatar
Henning Janssen committed
305 306
          vTemp%mmpMat = 0.0 !To avoid errors later on (When ldaUAdjEnpara is T the density
                             !is carried over after vgen)
307
          CALL timestart("Updating energy parameters")
Daniel Wortmann's avatar
Daniel Wortmann committed
308
          CALL enpara%update(mpi%mpi_comm,atoms,vacuum,input,vToT)
309
          CALL timestop("Updating energy parameters")
310 311
          !IF(.not.input%eig66(1))THEN
            CALL eigen(mpi,stars,sphhar,atoms,xcpot,sym,kpts,vacuum,input,&
312 313
                       cell,enpara,banddos,noco,oneD,mpdata,hybinp,hybdat,&
                       iter,eig_id,results,inDen,vTemp,vx,hub1)
314
          !ENDIF
315 316 317
          vTot%mmpMat = vTemp%mmpMat
!!$          eig_idList(pc) = eig_id
          CALL timestop("eigen")
318 319

          ! add all contributions to total energy
320
#ifdef CPP_MPI
321
          ! send all result of local total energies to the r
322
          IF (hybinp%l_hybrid.AND.hybdat%l_calhf) THEN
323 324 325 326 327 328 329 330 331 332
             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
333
          END IF
334 335
#endif

336
          ! WRITE(6,fmt='(A)') 'Starting 2nd variation ...'
337
          IF (noco%l_soc.AND..NOT.noco%l_noco) &
338
             CALL eigenso(eig_id,mpi,stars,vacuum,atoms,sphhar,&
339
                          sym,cell,noco,input,kpts, oneD,vTot,enpara,results,hub1)
Matthias Redies's avatar
Matthias Redies committed
340
          CALL timestop("gen. of hamil. and diag. (total)")
341 342

#ifdef CPP_MPI
343
          CALL MPI_BARRIER(mpi%mpi_comm,ierr)
344 345
#endif

346
          ! fermi level and occupancies
347
          IF (noco%l_soc.AND.(.NOT.noco%l_noco)) input%neig = 2*input%neig
348 349 350

	  IF (input%gw.GT.0) THEN
	    IF (mpi%irank.EQ.0) THEN
351
	       CALL writeBasis(input,noco,kpts,atoms,sym,cell,enpara,vTot,vCoul,vx,mpi,&
352 353 354 355 356 357 358
		  	     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

359
          !IF ((mpi%irank.EQ.0)) THEN
360 361
             CALL timestart("determination of fermi energy")

362
             IF (noco%l_soc.AND.(.NOT.noco%l_noco)) THEN
363
                input%zelec = input%zelec*2
Daniel Wortmann's avatar
Daniel Wortmann committed
364
                CALL fermie(eig_id,mpi,kpts,input,noco,enpara%epara_min,cell,results)
365 366 367 368
                results%seigscv = results%seigscv/2
                results%ts = results%ts/2
                input%zelec = input%zelec/2
             ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
369
                CALL fermie(eig_id,mpi,kpts,input,noco,enpara%epara_min,cell,results)
370 371
             ENDIF
             CALL timestop("determination of fermi energy")
372

373 374 375 376 377 378 379 380 381
!!$          !+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
!!$
382
!!$             CALL bs_comfort(eig_id,input,noco,kpts%nkpt,pc)
383 384 385 386 387 388 389 390 391
!!$
!!$             IF(pc.EQ.wann%nparampts)THEN
!!$                CLOSE(777)
!!$                CLOSE(778)
!!$                CLOSE(779)
!!$                CLOSE(780)
!!$             END IF
!!$          END IF
!!$          !-Wannier
392

393
          !ENDIF
394
#ifdef CPP_MPI
395 396
          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)
397 398
#endif

399 400
          IF (forcetheo%eval(eig_id,atoms,kpts,sym,cell,noco,input,mpi,oneD,enpara,vToT,results)) THEN
             IF (noco%l_soc.AND.(.NOT.noco%l_noco)) input%neig=input%neig/2
401 402
             CYCLE forcetheoloop
          ENDIF
403

404

Daniel Wortmann's avatar
Daniel Wortmann committed
405 406
          !+Wannier functions
          IF ((input%l_wann).AND.(.NOT.wann%l_bs_comf)) THEN
407
             CALL wannier(mpi,input,kpts,sym,atoms,stars,vacuum,sphhar,oneD,&
Daniel Wortmann's avatar
Daniel Wortmann committed
408 409 410 411
                  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
412

413
          ! charge density generation
414
          CALL timestart("generation of new charge density (total)")
415
          CALL outDen%init(stars,atoms,sphhar,vacuum,noco,input%jspins,POTDEN_TYPE_DEN)
416
          outDen%iter = inDen%iter
417
          CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum, &
418
                      kpts,atoms,sphhar,stars,sym,&
419
                      enpara,cell,noco,vTot,results,oneD,coreSpecInput,&
Henning Janssen's avatar
Henning Janssen committed
420
                      archiveType,xcpot,outDen,EnergyDen,gOnsite,hub1)
421
          !The density matrix for DFT+Hubbard1 only changes in hubbard1_setup and is kept constant otherwise
Henning Janssen's avatar
Henning Janssen committed
422
          outDen%mmpMat(:,:,atoms%n_u+1:atoms%n_u+atoms%n_hia,:) = inDen%mmpMat(:,:,atoms%n_u+1:atoms%n_u+atoms%n_hia,:)
423
          IF ((sliceplot%iplot.NE.0 ).AND.(mpi%irank==0) ) THEN
Robin Hilgers's avatar
Robin Hilgers committed
424
!               CDN including core charge
425
               ! CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, &
Robin Hilgers's avatar
Robin Hilgers committed
426
!                               cell, noco, outDen, PLOT_OUTDEN_Y_CORE, sliceplot)
427
!!               CDN subtracted by core charge
428
               ! CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, &
Robin Hilgers's avatar
Robin Hilgers committed
429
!                               cell, noco, outDen, PLOT_OUTDEN_N_CORE, sliceplot)
430
          END IF
431

Gregor Michalicek's avatar
Gregor Michalicek committed
432 433 434
          IF (input%l_rdmft) THEN
             SELECT TYPE(xcpot)
                TYPE IS(t_xcpot_inbuild)
435
                   CALL rdmft(eig_id,mpi,input,kpts,banddos,sliceplot,cell,atoms,enpara,stars,vacuum,&
436 437
                              sphhar,sym,field,vTot,vCoul,oneD,noco,xcpot,mpdata,hybinp,hybdat,&
                              results,coreSpecInput,archiveType,outDen)
Gregor Michalicek's avatar
Gregor Michalicek committed
438 439
             END SELECT
          END IF
440

441
          IF (noco%l_soc.AND.(.NOT.noco%l_noco)) input%neig=input%neig/2
442

443
#ifdef CPP_MPI
444
          CALL MPI_BCAST(enpara%evac,SIZE(enpara%evac),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
445 446 447
          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)
448

449 450 451 452 453
          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)
454
                ENDIF
455 456 457
             ENDDO
             IF (noco%l_constr) THEN
                CALL MPI_BCAST(noco%b_con,SIZE(noco%b_con),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
458
             ENDIF
459
          ENDIF
460
#endif
461
          CALL timestop("generation of new charge density (total)")
462
!START Rot For Testing (HIGHLY EXPERIMENTAL ROUTINE)
463
!IF (.FALSE.) CALL rotateMagnetFromSpinAxis(noco,vacuum,sphhar,stars,sym,oneD,cell,input,atoms,outDen,inDen)
464
!END Rot For Testing (HIGHLY EXPERIMENTAL ROUTINE)
465
!!$             !----> output potential and potential difference
Daniel Wortmann's avatar
Daniel Wortmann committed
466
!!$             IF (disp) THEN
467 468 469
!!$                reap = .FALSE.
!!$                input%total = .FALSE.
!!$                CALL timestart("generation of potential (total)")
Matthias Redies's avatar
Matthias Redies committed
470
!!$                CALL vgen(hybinp,reap,input,xcpot, atoms,sphhar,stars,vacuum,sym,&
Daniel Wortmann's avatar
Daniel Wortmann committed
471
!!$                     cell,oneD,sliceplot,mpi, results,noco,outDen,inDenRot,vTot,vx,vCoul)
472 473 474 475
!!$                CALL timestop("generation of potential (total)")
!!$
!!$                CALL potdis(stars,vacuum,atoms,sphhar, input,cell,sym)
!!$             END IF
476

477
             ! total energy
478
             CALL timestart('determination of total energy')
479
             CALL totale(mpi,atoms,sphhar,stars,vacuum,sym,input,noco,cell,oneD,&
480
                         xcpot,hybdat,vTot,vCoul,iter,inDen,results)
481
             CALL timestop('determination of total energy')
Matthias Redies's avatar
Matthias Redies committed
482
          IF (hybinp%l_hybrid) CALL close_eig(eig_id)
483 484 485 486

       END DO forcetheoloop

       CALL forcetheo%postprocess()
487

488
       CALL enpara%mix(mpi%mpi_comm,atoms,vacuum,input,vTot%mt(:,0,:,:),vtot%vacz)
489
       field2 = field
490
       ! mix input and output densities
491
       CALL mix_charge(field2,mpi,(iter==input%itmax.OR.judft_was_argument("-mix_io")),&
492
            stars,atoms,sphhar,vacuum,input,&
Henning Janssen's avatar
Henning Janssen committed
493
            sym,cell,noco,oneD,archiveType,xcpot,iter,inDen,outDen,results,hub1%l_runthisiter)
494 495
!Plots of mixed density
       IF ((sliceplot%iplot.NE.0 ).AND.(mpi%irank==0) ) THEN
496
!               CDN including core charge
497 498
!                CALL makeplots(stars, atoms, sphhar, vacuum, input, oneD, sym, &
!                               cell, noco, outDen, PLOT_MIXDEN_Y_CORE, sliceplot)
499
!!               CDN subtracted by core charge
500 501 502
!                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)
503
       END IF
504

505
       IF(mpi%irank == 0) THEN
506
         WRITE (6,FMT=8130) iter
507
8130     FORMAT (/,5x,'******* it=',i3,'  is completed********',/,/)
508
         WRITE(*,*) "Iteration:",iter," Distance:",results%last_distance
509 510 511 512 513 514 515 516
         !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
517
         CALL timestop("Iteration")
518
       END IF ! mpi%irank.EQ.0
519

520
#ifdef CPP_MPI
521
       CALL MPI_BCAST(results%last_distance,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
522 523
       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)
524
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
525
#endif
526 527
       CALL priv_geo_end(mpi)

528
       l_cont = .TRUE.
Matthias Redies's avatar
Matthias Redies committed
529
       IF (hybinp%l_hybrid) THEN
530
          IF(hybdat%l_calhf) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
531 532 533 534 535 536
             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
537
          IF (hybdat%l_subvxc) THEN
538 539
             results%te_hfex%valence = 0
          END IF
540 541
       ELSE
          l_cont = l_cont.AND.(iter < input%itmax)
542
          ! MetaGGAs need a at least 2 iterations
543
          l_cont = l_cont.AND.((input%mindistance<=results%last_distance).OR.input%l_f &
544
                               .OR. (xcpot%exc_is_MetaGGA() .and. iter == 1))
Henning Janssen's avatar
Henning Janssen committed
545
          !If we have converged run hia if the density matrix has not converged
Henning Janssen's avatar
Henning Janssen committed
546
          IF(atoms%n_hia>0) THEN
547
             hub1%l_runthisiter = .NOT.l_cont.AND.(input%minoccDistance<=results%last_occdistance&
548
                                  .OR.input%minmatDistance<=results%last_mmpMatdistance)
Henning Janssen's avatar
Henning Janssen committed
549
             !Run after first overall iteration to generate a starting density matrix
550
             hub1%l_runthisiter = hub1%l_runthisiter.OR.(iter==1.AND.(hub1%iter == 0&
Henning Janssen's avatar
Henning Janssen committed
551
                                  .AND.ALL(vTot%mmpMat(:,:,atoms%n_u+1:atoms%n_u+atoms%n_hia,:).EQ.0.0)))
552
             hub1%l_runthisiter = hub1%l_runthisiter.AND.(iter < input%itmax)
Henning Janssen's avatar
Henning Janssen committed
553 554 555
             !Prevent that the scf loop terminates
             l_cont = l_cont.OR.hub1%l_runthisiter
          ENDIF
556
          CALL check_time_for_next_iteration(iter,l_cont)
557
       END IF
558

559
       !CALL writeTimesXML()
560

561 562
       IF (mpi%irank.EQ.0) THEN
          IF (isCurrentXMLElement("iteration")) CALL closeXMLElement('iteration')
563 564
       END IF

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

570
    END DO scfloop ! DO WHILE (l_cont)
571

572
    CALL add_usage_data("Iterations",iter)
573 574

    IF (mpi%irank.EQ.0) CALL closeXMLElement('scfLoop')
575 576

    CALL close_eig(eig_id)
577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594
    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
595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610
            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)
611
         ENDIF
612
         call mixing_history_reset(mpi)
613 614 615
      ENDIF
      CALL juDFT_end(" GEO new inp.xml created ! ",mpi%irank)
    END SUBROUTINE priv_geo_end
616

617 618
  END SUBROUTINE fleur_execute
END MODULE m_fleur