IffGit has a new shared runner for building Docker images in GitLab CI. Visit https://iffgit.fz-juelich.de/examples/ci-docker-in-docker for more details.

fleur.F90 17.3 KB
Newer Older
1

2
3
4
5
6
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------
7
MODULE m_fleur
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
43
  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
44
    USE m_cdn_io
45
    USE m_broyd_io
46
    USE m_qfix
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
    USE m_vgen
    USE m_rhodirgen
    USE m_writexcstuff
    USE m_vmatgen
    USE m_eigen
    USE m_eigenso
    USE m_fermie
    USE m_force0
    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_wann_optional
    USE m_wannier
    USE m_bs_comfort
    USE m_gen_map
    USE m_dwigner
    USE m_ylm
68
#ifdef CPP_MPI
69
    USE m_mpi_bc_all,  ONLY : mpi_bc_all
70
    USE m_mpi_bc_potden
71
#endif
72
73
74
75
76
77
78
    USE m_eig66_io,   ONLY : open_eig, close_eig
    IMPLICIT NONE

    INTEGER,INTENT(IN) :: mpi_comm

    !     Types, these variables contain a lot of data!
    TYPE(t_input)    :: input
79
    TYPE(t_field)    :: field
80
81
82
83
84
85
86
87
88
89
90
    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
91
    TYPE(t_enpara)   :: enpara
92
93
94
95
96
97
    TYPE(t_xcpot)    :: xcpot
    TYPE(t_results)  :: results
    TYPE(t_kpts)     :: kpts
    TYPE(t_hybrid)   :: hybrid
    TYPE(t_oneD)     :: oneD
    TYPE(t_mpi)      :: mpi
98
    TYPE(t_coreSpecInput) :: coreSpecInput
99
    TYPE(t_wann)     :: wann
100
    TYPE(t_potden)   :: vTot,vx,vCoul,vTemp
101
    TYPE(t_potden)   :: inDen, outDen
102
    CLASS(t_forcetheo),ALLOCATABLE:: forcetheo
103
104

    !     .. Local Scalars ..
105
    INTEGER:: eig_id, archiveType
106
107
    INTEGER:: n,it,ithf
    LOGICAL:: reap,l_opti,l_cont,l_qfix, l_wann_inp
108
    REAL   :: fermiEnergyTemp, fix
109
#ifdef CPP_MPI
110
111
    INCLUDE 'mpif.h'
    INTEGER:: ierr(2)
112
#endif
Gregor Michalicek's avatar
Gregor Michalicek committed
113

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


123

124
125
126
127
    IF (l_opti) &
         CALL OPTIONAL(mpi,atoms,sphhar,vacuum,DIMENSION,&
         stars,input,sym,cell,sliceplot,obsolete,xcpot,noco,oneD)
 
128
129

    !+Wannier
130
131
    INQUIRE (file='wann_inp',exist=l_wann_inp)
    input%l_wann = input%l_wann.OR.l_wann_inp
132
133
    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')
134
       CALL wann_optional(input,kpts,atoms,sym,cell,oneD,noco,wann)
135
136
137
    END IF
    IF (wann%l_gwf) input%itmax = 1

138
    !-Wannier
139
140
141
142
143


    it     = 0
    ithf   = 0
    l_cont = ( it < input%itmax )
144
    
145
146
147
    results%last_distance = -1.0
    IF (mpi%irank.EQ.0) CALL openXMLElementNoAttributes('scfLoop')

148
    ! Initialize and load inDen density (start)
149
    CALL inDen%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
150
 
151
152
153
    archiveType = CDN_ARCHIVE_TYPE_CDN1_const
    IF (noco%l_noco) archiveType = CDN_ARCHIVE_TYPE_NOCO_const

154
155
    IF(mpi%irank.EQ.0) THEN
       CALL readDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
156
            0,fermiEnergyTemp,l_qfix,inDen)
157
       CALL timestart("Qfix")
158
       CALL qfix(stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.false.,fix)
159
160
       CALL timestop("Qfix")
       CALL writeDensity(stars,vacuum,atoms,cell,sphhar,input,sym,oneD,archiveType,CDN_INPUT_DEN_const,&
161
            0,-1.0,0.0,.FALSE.,inDen)
162
    END IF
163
    ! Initialize and load inDen density (end)
164

165
    ! Initialize potentials (start)
166
167
168
169
    CALL vTot%init(stars,atoms,sphhar,vacuum,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTTOT)
    CALL vCoul%init(stars,atoms,sphhar,vacuum,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTCOUL)
    CALL vx%init(stars,atoms,sphhar,vacuum,DIMENSION%jspd,.FALSE.,POTDEN_TYPE_POTCOUL)
    CALL vTemp%init(stars,atoms,sphhar,vacuum,DIMENSION%jspd,noco%l_noco,POTDEN_TYPE_POTTOT)
170
    ! Initialize potentials (end)
171

172
    scfloop:DO WHILE (l_cont)
173
174

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

178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
!!$       !+t3e
!!$       IF (input%alpha.LT.10.0) THEN
!!$
!!$          IF (it.GT.1) THEN
!!$             input%alpha = input%alpha - NINT(input%alpha)
!!$          END IF

       CALL resetIterationDependentTimers()
       CALL timestart("Iteration")
       IF (mpi%irank.EQ.0) THEN
          !-t3e
          WRITE (6,FMT=8100) it
          WRITE (16,FMT=8100) it
8100      FORMAT (/,10x,'   it=    ',i5)


          !      ----> potential generator
          !
          reap=.NOT.obsolete%disp
          input%total = .TRUE.
       ENDIF !mpi%irank.eq.0
199
#ifdef CPP_MPI
200
       CALL MPI_BCAST(input%total,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
201
#endif
202

203
#ifdef CPP_MPI
204
       CALL mpi_bc_potden(mpi,stars,sphhar,atoms,input,vacuum,oneD,noco,inDen)
205
206
#endif

Daniel Wortmann's avatar
Daniel Wortmann committed
207

208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
       IF ( noco%l_soc ) THEN
          dimension%neigd2 = dimension%neigd*2
       ELSE
          dimension%neigd2 = dimension%neigd
       END IF


       !HF
       IF (hybrid%l_hybrid) CALL  calc_hybrid(hybrid,kpts,atoms,input,DIMENSION,mpi,noco,&
            cell,vacuum,oneD,banddos,results,sym,xcpot,vTot,it)
       !#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")
232
       CALL vgen(hybrid,field,input,xcpot,DIMENSION, atoms,sphhar,stars,vacuum,&
233
            sym,obsolete,cell, oneD,sliceplot,mpi ,results,noco,inDen,vTot,vx,vCoul)
234
235
       CALL timestop("generation of potential")

236

Daniel Wortmann's avatar
Daniel Wortmann committed
237
238


239
#ifdef CPP_MPI
240
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
241
242
#endif

243
244
       CALL forcetheo%start()

245
       forcetheoloop:DO WHILE(forcetheo%next_job(it==input%itmax,noco))
246

247

248
249
250
          CALL timestart("generation of hamiltonian and diagonalization (total)")
          CALL timestart("eigen")
          vTemp = vTot
251
          CALL enpara%update(mpi,atoms,vacuum,input,vToT)
252
          CALL eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,&
253
               sym,kpts,DIMENSION,vacuum,input,cell,enpara,banddos,noco,oneD,hybrid,&
254
               it,eig_id,results,inDen,vTemp,vx)
255
256
257
258
259
260
          vTot%mmpMat = vTemp%mmpMat
!!$          eig_idList(pc) = eig_id
          CALL timestop("eigen")
          !
          !                   add all contributions to total energy
          !
261
#ifdef CPP_MPI
262
263
264
265
266
267
268
269
270
271
272
273
274
          ! send all result of local total energies to the r
          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))
             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%valence,MPI_IN_PLACE,&
                  1,MPI_REAL8,MPI_SUM,0, mpi%mpi_comm,ierr(1))
             CALL MPI_Reduce(results%te_hfex%core,MPI_IN_PLACE,&
                  1,MPI_REAL8,MPI_SUM,0, mpi%mpi_comm,ierr(1))
          ENDIF
          !                                  END IF
275
276
#endif

Daniel Wortmann's avatar
Daniel Wortmann committed
277

278

279
280
281
          ! WRITE(6,fmt='(A)') 'Starting 2nd variation ...'
          IF (noco%l_soc.AND..NOT.noco%l_noco) &
               CALL eigenso(eig_id,mpi,DIMENSION,stars,vacuum,atoms,sphhar,&
282
               obsolete,sym,cell,noco,input,kpts, oneD,vTot,enpara)
283
          CALL timestop("generation of hamiltonian and diagonalization (total)")
284
285

#ifdef CPP_MPI
286
          CALL MPI_BARRIER(mpi%mpi_comm,ierr)
287
288
#endif

289
290
291
292
293
294
295
296
297
298
299
300
          !
          !              ----> fermi level and occupancies

          IF ( noco%l_soc .AND. (.NOT. noco%l_noco) ) DIMENSION%neigd = 2*DIMENSION%neigd
          IF( .NOT. ALLOCATED(results%w_iks) )&
               ALLOCATE ( results%w_iks(DIMENSION%neigd,kpts%nkpt,DIMENSION%jspd) )
          IF ( (mpi%irank.EQ.0)) THEN
             CALL timestart("determination of fermi energy")

             IF ( noco%l_soc .AND. (.NOT. noco%l_noco) ) THEN
                input%zelec = input%zelec*2
                CALL fermie(eig_id,mpi,kpts,obsolete,&
301
                     input,noco,enpara%epara_min,cell,results)
302
303
304
305
306
                results%seigscv = results%seigscv/2
                results%ts = results%ts/2
                input%zelec = input%zelec/2
             ELSE
                CALL fermie(eig_id,mpi,kpts,obsolete,&
307
                     input,noco,enpara%epara_min,cell,results)
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
             ENDIF
             CALL timestop("determination of fermi energy")
!!$             
!!$          !+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
          ENDIF
331
#ifdef CPP_MPI
332
333
          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)
334
335
336
#endif


337
          IF (forcetheo%eval(eig_id,DIMENSION,atoms,kpts,sym,&
338
339
340
341
               cell,noco, input,mpi, oneD,enpara,vToT,results)) THEN
             IF ( noco%l_soc .AND. (.NOT. noco%l_noco) ) DIMENSION%neigd=DIMENSION%neigd/2
             CYCLE forcetheoloop
          ENDIF
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356


          CALL force_0(results)!        ----> initialise force_old
          !        ----> charge density
          
!!$          !+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_idList,(sym%invs).AND.(.NOT.noco%l_soc).AND.(.NOT.noco%l_noco),kpts%nkpt)
!!$          END IF
!!$          IF (wann%l_gwf) CALL juDFT_error("provide wann_inp if l_gwf=T", calledby = "fleur")
!!$          !-Wannier

          CALL timestart("generation of new charge density (total)")
357
          CALL outDen%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
358
359
360
          outDen%iter = inDen%iter
          CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
               DIMENSION,kpts,atoms,sphhar,stars,sym,obsolete,&
361
               enpara,cell,noco,vTot,results,oneD,coreSpecInput,&
362
363
364
365
               archiveType,outDen)

          IF ( noco%l_soc .AND. (.NOT. noco%l_noco) ) DIMENSION%neigd=DIMENSION%neigd/2
          !+t3e
366
#ifdef CPP_MPI
367
368
369
370
371
372
373
374
375
          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)
          
          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)
376
                ENDIF
377
378
379
             ENDDO
             IF (noco%l_constr) THEN
                CALL MPI_BCAST(noco%b_con,SIZE(noco%b_con),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
380
             ENDIF
381
          ENDIF
382
#endif
383
          CALL timestop("generation of new charge density (total)")
384
385
          IF (mpi%irank.EQ.0) THEN
             !-t3e
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
             
!!$             !----> 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
             
             !----> total energy
             CALL timestart('determination of total energy')
             CALL totale(atoms,sphhar,stars,vacuum,DIMENSION,&
                  sym,input,noco,cell,oneD,xcpot,hybrid,vTot,vCoul,it,inDen,results)
403

404
405
406
407
408
409
410
411
             CALL timestop('determination of total energy')
          ENDIF ! mpi%irank.EQ.0
          IF ( hybrid%l_hybrid ) CALL close_eig(eig_id)

       END DO forcetheoloop

       CALL forcetheo%postprocess()
       
412
       CALL enpara%mix(mpi,atoms,vacuum,input,vTot%mt(:,0,:,:),vtot%vacz)
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
       IF (mpi%irank.EQ.0) THEN
          !          ----> mix input and output densities
          CALL timestart("mixing")
          CALL mix(stars,atoms,sphhar,vacuum,input,sym,cell,noco,oneD,hybrid,archiveType,inDen,outDen,results)
          CALL timestop("mixing")
          
          WRITE (6,FMT=8130) it
          WRITE (16,FMT=8130) it
8130      FORMAT (/,5x,'******* it=',i3,'  is completed********',/,/)
          WRITE(*,*) "Iteration:",it," Distance:",results%last_distance
          CALL timestop("Iteration")
          !+t3e
       ENDIF ! mpi%irank.EQ.0
       
          
428
#ifdef CPP_MPI
429
430
       CALL MPI_BCAST(results%last_distance,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
431
#endif
432
433
       CALL priv_geo_end(mpi)

434
435
436
437
438
439
440
441
!!$       IF ( hybrid%l_calhf ) ithf = ithf + 1
!!$    IF ( hybrid%l_subvxc ) THEN
!!$       l_cont = ( ithf < input%itmax )
!!$       results%te_hfex%core    = 0
!!$       results%te_hfex%valence = 0
!!$    ELSE
       l_cont = ( it < input%itmax )
!!$    END IF
442
443
444
445
446
447
448
449
450
       CALL writeTimesXML()
       CALL check_time_for_next_iteration(it,l_cont)

       l_cont = l_cont.AND.((input%mindistance<=results%last_distance).OR.input%l_f)

       IF ((mpi%irank.EQ.0).AND.(isCurrentXMLElement("iteration"))) THEN
          CALL closeXMLElement('iteration')
       END IF

451
    END DO scfloop ! DO WHILE (l_cont)
452
453
454

    IF (mpi%irank.EQ.0) CALL closeXMLElement('scfLoop')
    CALL juDFT_end("all done",mpi%irank)
455
    
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
  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
478
         CALL resetBroydenHistory()
479
480
481
      ENDIF
      CALL juDFT_end(" GEO new inp.xml created ! ",mpi%irank)
    END SUBROUTINE priv_geo_end
482
    
483
484
  END SUBROUTINE fleur_execute
END MODULE m_fleur