fleur.F90 17.5 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
    USE m_vgen
    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
61
    USE m_rdmft
62 63 64 65 66 67
    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_potden
70
#endif
71
    USE m_eig66_io,   ONLY : open_eig, close_eig
72
    USE m_chase_diag
73 74 75 76 77 78
    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
    TYPE(t_results)  :: results
    TYPE(t_kpts)     :: kpts
    TYPE(t_hybrid)   :: hybrid
    TYPE(t_oneD)     :: oneD
    TYPE(t_mpi)      :: mpi
97
    TYPE(t_coreSpecInput) :: coreSpecInput
98
    TYPE(t_wann)     :: wann
99
    TYPE(t_potden)   :: vTot,vx,vCoul,vTemp
100
    TYPE(t_potden)   :: inDen, outDen
101
    CLASS(t_xcpot),ALLOCATABLE    :: xcpot
102
    CLASS(t_forcetheo),ALLOCATABLE:: forcetheo
103 104

    !     .. Local Scalars ..
105 106 107
    INTEGER:: eig_id, archiveType
    INTEGER:: n,it,ithf
    LOGICAL:: 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
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
    IF (l_opti) CALL optional(mpi,atoms,sphhar,vacuum,dimension,&
                              stars,input,sym,cell,sliceplot,obsolete,xcpot,noco,oneD)
124
 
125 126

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

135
    !-Wannier
136

137
#ifdef CPP_CHASE
138
    CALL init_chase(mpi,dimension,atoms,kpts,noco,sym%invs.AND..NOT.noco%l_noco)
139
#endif
140 141 142 143

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

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

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

164
    ! Initialize potentials (start)
165 166 167 168
    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)
169
    ! Initialize potentials (end)
170

171
    scfloop:DO WHILE (l_cont)
172 173

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

177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196
!!$       !+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
          !
          input%total = .TRUE.
       ENDIF !mpi%irank.eq.0
197
#ifdef CPP_MPI
198
       CALL MPI_BCAST(input%total,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
199
#endif
200

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

Daniel Wortmann's avatar
Daniel Wortmann committed
205

206 207 208 209 210 211 212 213
       IF ( noco%l_soc ) THEN
          dimension%neigd2 = dimension%neigd*2
       ELSE
          dimension%neigd2 = dimension%neigd
       END IF


       !HF
214 215 216
       IF (hybrid%l_hybrid) THEN
          SELECT TYPE(xcpot)
          TYPE IS(t_xcpot_inbuild)
217
             CALL calc_hybrid(hybrid,kpts,atoms,input,DIMENSION,mpi,noco,&
218
                              cell,oneD,results,sym,xcpot,vTot,it)
219 220
          END SELECT
       ENDIF
221 222 223 224 225 226 227 228 229 230 231 232 233 234
       !#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")
235
       CALL vgen(hybrid,field,input,xcpot,DIMENSION, atoms,sphhar,stars,vacuum,&
236
            sym,obsolete,cell, oneD,sliceplot,mpi ,results,noco,inDen,vTot,vx,vCoul)
237 238
       CALL timestop("generation of potential")

239

Daniel Wortmann's avatar
Daniel Wortmann committed
240 241


242
#ifdef CPP_MPI
243
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
244 245
#endif

246 247
       CALL forcetheo%start()

248
       forcetheoloop:DO WHILE(forcetheo%next_job(it==input%itmax,noco))
249

250

251 252 253
          CALL timestart("generation of hamiltonian and diagonalization (total)")
          CALL timestart("eigen")
          vTemp = vTot
254
          CALL enpara%update(mpi,atoms,vacuum,input,vToT)
255
          CALL eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,&
256
                     sym,kpts,DIMENSION,vacuum,input,cell,enpara,banddos,noco,oneD,hybrid,&
257
                     it,eig_id,results,inDen,vTemp,vx)
258 259 260 261 262 263
          vTot%mmpMat = vTemp%mmpMat
!!$          eig_idList(pc) = eig_id
          CALL timestop("eigen")
          !
          !                   add all contributions to total energy
          !
264
#ifdef CPP_MPI
265 266 267 268 269 270 271 272 273 274 275 276 277
          ! 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
278 279
#endif

Daniel Wortmann's avatar
Daniel Wortmann committed
280

281

282 283 284
          ! 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,&
285
               obsolete,sym,cell,noco,input,kpts, oneD,vTot,enpara,results)
286
          CALL timestop("generation of hamiltonian and diagonalization (total)")
287 288

#ifdef CPP_MPI
289
          CALL MPI_BARRIER(mpi%mpi_comm,ierr)
290 291
#endif

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

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


338
          IF (forcetheo%eval(eig_id,DIMENSION,atoms,kpts,sym,&
339 340 341 342
               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
343 344 345 346 347 348 349 350 351 352 353 354 355 356 357


          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)")
358
          CALL outDen%init(stars,atoms,sphhar,vacuum,input%jspins,noco%l_noco,POTDEN_TYPE_DEN)
359 360
          outDen%iter = inDen%iter
          CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
361
               DIMENSION,kpts,atoms,sphhar,stars,sym,&
362
               enpara,cell,noco,vTot,results,oneD,coreSpecInput,&
363 364
               archiveType,outDen)

365
          IF (.FALSE.) CALL rdmft(eig_id,mpi,input,kpts,banddos,cell,atoms,enpara,stars,vacuum,dimension,&
366
                                  sphhar,sym,field,vTot,oneD,noco,results)
367

368 369
          IF ( noco%l_soc .AND. (.NOT. noco%l_noco) ) DIMENSION%neigd=DIMENSION%neigd/2
          !+t3e
370
#ifdef CPP_MPI
371 372 373 374 375 376 377 378 379
          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)
380
                ENDIF
381 382 383
             ENDDO
             IF (noco%l_constr) THEN
                CALL MPI_BCAST(noco%b_con,SIZE(noco%b_con),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
384
             ENDIF
385
          ENDIF
386
#endif
387
          CALL timestop("generation of new charge density (total)")
388 389
          IF (mpi%irank.EQ.0) THEN
             !-t3e
390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406
             
!!$             !----> 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)
407

408 409 410 411 412 413 414 415
             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()
       
416
       CALL enpara%mix(mpi,atoms,vacuum,input,vTot%mt(:,0,:,:),vtot%vacz)
417 418 419 420 421 422 423 424 425 426 427 428 429 430 431
       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
       
          
432
#ifdef CPP_MPI
433 434
       CALL MPI_BCAST(results%last_distance,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
       CALL MPI_BARRIER(mpi%mpi_comm,ierr)
435
#endif
436 437
       CALL priv_geo_end(mpi)

438 439 440 441 442 443 444 445
!!$       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
446 447 448 449 450 451 452 453 454
       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

455
    END DO scfloop ! DO WHILE (l_cont)
456 457 458

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