fleur.F90 40.2 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 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 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
MODULE m_fleur

   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
      !
      ! flapw --+-- dimens      ... dimensioning routine(s)
      !         +-- inped - rw_inp
      !         +-- setup       ... readin & setup
      !         +-- optional    ... plots, startdensity etc.
      !        +-+
      !        | +-- vgen       ... construct potential
      !        | +-- eigen      ... calculate eigenvalues & vectors
      !        |(+-- eigenso)       ... spin orbit coupling
      !        | +-- fermie     ... calculate fermi energy
      !        |(+-- force0)        ... initialize force
      !        | +-- cdngen     ... construct charge density
      !        |(+-- vgen)          ... output potential
      !        |(+-- potdis)        ... potential distance
      !        | +-- totale     ... calculate total energy
      !        | +-- mix        ... mix densities
      !        +-+
      !----------------------------------------
      !
      ! I/O file list
      !        2... time.info ...
      !        5... inp       ...
      !        6... out       ...
      !        8... pottot    ...
      !        9... nrp       ...
      !       11... potcoul   ...
      !       14... wkf2      ...
      !       15... QGpsi,CLASS,NLAindx,SYMOPS,LATTC,gwa,abcoeff,ecore
      !       16... inf       ...
      !       17... cdnc      ...
      !       18... plotin    ...
      !       19... cdnout    ...
      !       20... cdn_slice ...
      !       22... dirofmag  ...
      !       24... nocoinf   ...
      !       25... potmat    ...
      !       26... rhomat_inp & rhomat_out
      !       28... tmat      ...
      !       29... tlst      ...
      !       38... tmas      ...
      !       40... enpara    ...
      !       41... kpts      ...
      !       43... force.dat ...
      !       45... inp_new   ...
      !       51... stars     ...
      !       57... broyd     ...
      !       59... broyd.x   ...
      !       66... eig       ...
      !       67... eig.soc   ...
      !       71... cdn1      ...
      !       79... cdn9      ...
      !       85... dosinp    ...
      !       86... vacdos    ...
      !       87... tmp_vacdos ...
      !       88... vacwave   ...
      !       92... CBS_eigenvalues.1/2
      !       93... CBS_states
      !       99... qfix      ...
      !      101... Smat...
      !      102... Tmat..
      !      110... scratch unit for GW
      !      113... qpts      ...
      !      114... jenerg    ...
      !      115... jconst    ...
      !      116... qptsinfo  ...
      !      117... shells    ...
      !      118... MCinp     ...
      !      300... matelsso  (for ss+soc)
      !      301+XXX matelXXX (for ss+soc)
      !
      !----------------------------------------

      USE m_types
      USE m_constants
      USE m_fleur_init
      USE m_pldngen
      USE m_optional
      USE m_vgen
      USE m_rhodirgen
      USE m_writexcstuff
      USE m_vmatgen
      USE m_icorrkeys
      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_jcoff
!      USE m_jcoff2
!      USE m_ssomat
136

137 138
      USE m_wann_optional
      USE m_wannier
139
      USE m_bs_comfort
140

141 142 143 144 145 146 147 148
!      USE m_mixedbasis
!      USE m_coulomb
      USE m_gen_map
      USE m_dwigner
      USE m_gen_bz
!      USE m_generate_pntgpt
!      USE m_rotate_eig
      USE m_ylm
149
#ifdef CPP_MPI
150
      USE m_mpi_bc_all,  ONLY : mpi_bc_all
151
#endif
152 153 154 155 156 157 158 159 160 161 162 163 164 165
      USE m_eig66_io,   ONLY : open_eig, close_eig

      IMPLICIT NONE

      INTEGER,INTENT(IN) :: mpi_comm

      !     scalar bmu changed to array bmu(ntypd) to enable calculation of
      !     both ferromagnetic and antiferromagnetic CASE - r.pentcheva Jan'96
      !     another possibility for the antiferromagnetic case ist to start
      !     from the ferromagnetic charge density and flip the magnetic moment
      !     of certain atoms: logical variable lflip and integer nflip(ntype)
      !                                 r.pentcheva,kfa,Feb'96

      !     Types, these variables contain a lot of data!
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186
      TYPE(t_input)        :: input
      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_xcpot)        :: xcpot
      TYPE(t_results)      :: results
      TYPE(t_jij)          :: jij
      TYPE(t_kpts)         :: kpts
      TYPE(t_hybrid)       :: hybrid
      TYPE(t_oneD)         :: oneD
      TYPE(t_mpi)          :: mpi
      TYPE(t_wann)         :: wann
187 188

      !     .. Local Scalars ..
189 190
      INTEGER              :: eig_id
      INTEGER              :: i,it,ithf,jspin,n,pc
191
      LOGICAL              :: stop80,reap,l_endit,l_opti,l_cont,l_real
192
      !--- J<
193 194 195
      INTEGER              :: phn
      REAL, PARAMETER      :: tol = 1.e-8
      INTEGER              :: qcount ,imt,i_J,j_J
196 197 198
      !--- J>

      !     .. Local Arrays ..
199
      CHARACTER(8)         :: name(10)
200 201

      !     HF/hybrid-functionals/EXX
202 203 204 205
      INTEGER              ::  l,j,m1,m2,isym,iisym
      INTEGER              ::  ok
      LOGICAL              ::  l_restart
      COMPLEX              ::  cdum
206
#ifdef CPP_MPI
207
      include 'mpif.h'
208
      integer              :: ierr(2)
209
#endif
210
      INTEGER, ALLOCATABLE :: eig_idList(:)
211
      mpi%mpi_comm = mpi_comm
212
         
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237
      CALL timestart("Initialization")
      CALL fleur_init(mpi,input,dimension,atoms,sphhar,cell,stars,sym,noco,vacuum,&
                      sliceplot,banddos,obsolete,enpara,xcpot,results,jij,kpts,hybrid,&
                      oneD,wann,l_opti)
      CALL timestop("Initialization")

      IF (l_opti) THEN
         IF (sliceplot%iplot .AND. (mpi%irank==0) ) THEN
            IF (noco%l_noco) THEN
               CALL pldngen(sym,stars,atoms,sphhar,vacuum,&
                            cell,input,noco,oneD,sliceplot)
            END IF
         END IF
         CALL OPTIONAL(mpi,atoms,sphhar,vacuum,dimension,&
                       stars,input,sym,cell,sliceplot,obsolete,&
                       xcpot,noco,oneD)
      END IF

      IF (sliceplot%iplot)   CALL juDFT_end("density plot o.k.",mpi%irank)
      IF (input%strho)       CALL juDFT_end("starting density generated",mpi%irank)
      IF (input%swsp)        CALL juDFT_end("spin polarised density generated",mpi%irank)
      IF (input%lflip)       CALL juDFT_end("magnetic moments flipped",mpi%irank)
      IF (obsolete%l_f2u)    CALL juDFT_end("conversion to unformatted",mpi%irank)
      IF (obsolete%l_u2f)    CALL juDFT_end("conversion to formatted",mpi%irank)
      IF (input%l_bmt)       CALL juDFT_end('"cdnbmt" written',mpi%irank)
238

239
!+Wannier
240 241 242
      input%l_wann = .FALSE.
      INQUIRE (file='wann_inp',exist=input%l_wann)
      IF (input%l_wann.AND.(mpi%irank==0).AND.(.NOT.wann%l_bs_comf)) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
243
         IF(mpi%isize.NE.1) CALL juDFT_error('No Wannier+MPI at the moment',calledby = 'fleur')
244 245 246
         CALL wann_optional(input,atoms,sym,cell,oneD,noco,wann)
      END IF
      IF (wann%l_gwf) input%itmax = 1
247 248 249
      l_real = (sym%invs).AND.(.NOT.noco%l_noco) ! Is this right? I took it from eigen.
!-Wannier

250
      ALLOCATE (eig_idList(wann%nparampts))
251

252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274
      l_restart = .TRUE.

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

      DO WHILE (l_cont)

         it = it + 1
         !+t3e

         IF (input%alpha.LT.10.0) THEN

            IF (it.GT.1) THEN
               obsolete%pot8 = .FALSE.
               input%alpha = input%alpha - NINT(input%alpha)
            END IF

            CALL resetIterationDependentTimers()
            CALL timestart("Iteration")
            IF (mpi%irank.EQ.0) THEN
275
               WRITE(*,"(a,i3)",advance="no") "Iteration:",it
276 277 278
               !-t3e
               WRITE (6,FMT=8100) it
               WRITE (16,FMT=8100) it
279
8100               FORMAT (/,10x,'   it=    ',i5)
280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301

               IF (.NOT.obsolete%pot8) THEN

                  !      ----> potential generator
                  !
                  !---> pk non-collinear
                  !--->        reload the density matrix from file rhomat_in
                  !--->        calculate spin-up and -down density for USE in the
                  !--->        potential generator and store the direction of
                  !--->        magnetization on file dirofmag
                  IF (noco%l_noco) THEN
                     CALL timestart("gen. spin-up and -down density")
                     CALL rhodirgen(dimension,sym,stars,atoms,sphhar,&
                                    vacuum,22,cell,input,oneD)
                     CALL timestop("gen. spin-up and -down density")
                  END IF
                  !---> pk non-collinear

                  reap = .NOT.obsolete%disp
                  input%total = .TRUE.
               END IF!(obsolete%pot8)
            END IF !mpi%irank.eq.0
302
#ifdef CPP_MPI
303
            CALL MPI_BCAST(input%total,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
304
#endif
305

306 307 308 309 310 311
            !--- J<
            IF(jij%l_jenerg) GOTO 234

            jij%alph1(:)=noco%alph(:)
            stop80= .FALSE.
            IF ((noco%l_soc.AND.noco%l_ss)) THEN
312 313
               IF ((.NOT.wann%l_gwf).AND.&
                   ((jij%l_J).OR.(jij%nqpt/=1).OR.(jij%nmagn/=1).OR.(jij%phnd/=1))) THEN
314 315 316 317 318 319 320 321 322 323
                  CALL juDFT_error("fleur: J-loop with ss+soc", calledby ="fleur")
               END IF
            END IF

            DO qcount=1,jij%nqpt

               IF (jij%l_J) THEN
                  noco%qss(:)=jij%qj(:,qcount)
                  jij%qn = ( noco%qss(1)**2 + noco%qss(2)**2 + noco%qss(3)**2 )
               END IF
324
               IF ((input%l_wann.OR.jij%l_J).AND.(mpi%irank.EQ.0)) THEN
325 326 327 328 329
                  WRITE(6,*) 'qss=(',noco%qss(1),',',noco%qss(2),',',noco%qss(3),')'
                  CALL timestart("Q-point for J_ij(total)")
               END IF

               !HF
330
#if defined(CPP_HF) && (2==1)
331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
               hybrid%l_subvxc = (xcpot%icorr == icorr_hf .OR.xcpot%icorr == icorr_pbe0.OR. &
                                  xcpot%icorr == icorr_hse.OR.xcpot%icorr == icorr_vhse)

               IF (noco%l_soc) THEN
                  dimension%neigd2 = dimension%neigd*2
               ELSE
                  dimension%neigd2 = dimension%neigd
               END IF

               IF (.NOT.ALLOCATED(results%w_iks)) THEN
                  ALLOCATE (results%w_iks(dimension%neigd2,kpts%nkpt,dimension%jspd))
               END IF

               IF((xcpot%icorr == icorr_hf .OR.xcpot%icorr == icorr_pbe0.OR.&
                   xcpot%icorr == icorr_hse.OR.xcpot%icorr == icorr_exx .OR.&
                   xcpot%icorr == icorr_vhse).AND.it == 1) THEN

                  CALL juDFT_WARN("Hybrid functionals not working in this version")
                  CALL timestart("generation of mixedbasis and coulombmatrix")

                  IF (mpi%irank == 0) THEN
                     WRITE(*,'(/A)',advance='no') ' calculation of mixedbasis...'
                  END IF

                  IF (it==1) THEN
                     eig_id=open_eig(mpi%mpi_comm,dimension%nbasfcn,dimension%neigd,kpts%nkpt(1),&
357 358
                                             dimension%jspd,atoms%lmaxd,atoms%nlod,atoms%ntype,atoms%nlotot,&
                                             noco%l_noco,.FALSE.,.FALSE.)
359 360 361 362
                      !               CALL open_eig(mpi_comm,
                      !     >              nbasfcn,neigd,nkpt(1),jspd,lmaxd,nlod,ntypd,nlotot,
                      !     >                 l_noco,.false.,.true.)
                      !DW TODO! eig_id has to be adjusted here
363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381
                  END IF

                  CALL mixedbasis(atoms,kpts,obsolete,&
                                  sphhar,dimension,input,&
                                  enpara,cell,vacuum,sym,&
                                  oneD,stars,xcpot,hybrid,&
                                  eig_id,mpi,l_restart)

                  IF (irank == 0) WRITE(*,'(A)')'...done'

                  IF (mpi%irank == 0) THEN
                     WRITE(*,'(A)',advance='no') ' calculation of coulomb matrix ...'
                  END IF

                  CALL coulombmatrix(mpi,obsolete,atoms,kpts,&
                                     cell,sym,hybrid,xcpot,l_restart,&
                                     oneD)

                  IF (irank == 0) WRITE(*,'(A)')'...done'
382
#ifdef CPP_MPI
383 384
                  CALL MPI_BCAST(hybrid%maxbasm1,1,MPI_INTEGER4,0,mpi%mpi_comm,ierr(1))
                  CALL MPI_BCAST(hybrid%radshmin,1,MPI_REAL8,0,mpi%mpi_comm,ierr(1))
385
#endif
386
                  CALL timestop("generation of mixedbasis and coulombmatrix")
387

388 389 390
                  IF (noco%l_soc) THEN
                     input%zelec = input%zelec * 2
                  END IF
391

392 393 394
                  IF (mpi%irank == 0) THEN
                     WRITE(*,'(A)',advance='no') ' start fermie....'
                  END IF
395

396 397
                  CALL fermie(eig_id,mpi,dimension,kpts,obsolete,atoms,&
                              input,noco,results,jij,cell)
398

399 400 401
                  IF (noco%l_soc) THEN
                     input%zelec = input%zelec / 2
                  END IF
402

403
                  IF (mpi%irank == 0) WRITE(*,'(A)') '...done'
404

405
               ELSE IF (it == 1) THEN ! allocate some dummy arrays
406
#endif
407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428
               IF (it==1) THEN ! temporary until HF is not excluded by #if any more
                  IF ( noco%l_soc ) THEN
                     dimension%neigd2 = dimension%neigd*2
                  ELSE
                     dimension%neigd2 = dimension%neigd
                  END IF
                  kpts%nkptf = 0; hybrid%maxindx = 0; hybrid%gptmd = 0; hybrid%maxgptm = 0
                  hybrid%maxgptm1 = 0; hybrid%maxgptm2 = 0; hybrid%maxindxm1 = 0; hybrid%maxindxm2 = 0
                  hybrid%maxlcutm1 = 0; hybrid%maxlcutm2 = 0; hybrid%maxindxp1 = 0; hybrid%maxindxp2 = 0
                  ALLOCATE(hybrid%gptm(0,0),hybrid%ngptm(0),hybrid%pgptm(0,0),hybrid%ngptm1(0))
                  ALLOCATE(hybrid%pgptm1(0,0),hybrid%ngptm2(0),hybrid%pgptm2(0,0),hybrid%basm1(0,0,0,0))
                  ALLOCATE(hybrid%basm2(0,0,0,0),hybrid%nindxm1(0,0),hybrid%nindxm2(0,0))
               END IF ! first iteration hybrids
               !HF

               DO pc = 1, wann%nparampts
!---> gwf
                  IF (wann%l_sgwf.or.wann%l_ms) THEN
                     noco%qss(:) = wann%param_vec(:,pc)
                     jij%qn = (noco%qss(1)**2 + noco%qss(2)**2 + noco%qss(3)**2)
                     noco%alph(:) = wann%param_alpha(:,pc)
                  ELSE IF (wann%l_socgwf) THEN
429 430
                     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)
431 432 433 434 435
                  END IF
!---< gwf

                  IF (.NOT.obsolete%pot8) THEN
                     CALL timestart("generation of potential")
436
                     IF (mpi%irank==0) WRITE(*,"(a)",advance="no") " * Potential generation "
437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465
                     CALL vgen(reap,input,xcpot,dimension,&
                               atoms,sphhar,stars,vacuum,&
                               sym,obsolete,cell,&
                               oneD,sliceplot,mpi,&
                               results,noco)

                     IF (mpi%irank.EQ.0) THEN
                        CALL timestop("generation of potential")

                        !---> pk non-collinear
                        !--->          generate the four component matrix potential from spin up
                        !--->          and down potentials and direction of the magnetic field
                        IF (noco%l_noco) THEN
                           CALL timestart("generation of potential-matrix")
                           CALL vmatgen(stars,&
                                        atoms,sphhar,vacuum,&
                                        sym,input,oneD,8,22,26)
                           CALL timestop("generation of potential-matrix")
                        END IF
                        !---> end pk non-collinear
                        !---> do some output for the tddft calculations:
                        IF (input%gw /= 0) THEN
                           CALL write_xcstuff(sphhar,atoms,dimension,sym,&
                                              stars,vacuum,input)
                        END IF
                     END IF ! mpi%irank.eq.0

                  !+t3e
                  END IF ! .not.obsolete%pot8
466
#ifdef CPP_MPI
467
                  CALL MPI_BARRIER(mpi%mpi_comm,ierr)
468 469
#endif

470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485
                  !          ----> eigenvalues and eigenfunctions
                  !
                  !--- J<
                  IF (jij%l_disp) THEN
                     jij%mtypes = 1
                     jij%nmagn = 1
                     jij%phnd = 1
                  END IF
                  i_J = 1
                  DO imt = 1, jij%mtypes
                     DO j_J = i_J, jij%nmagn
                        DO phn = 1, jij%phnd

                           input%eigvar(1)= .TRUE.
                           input%eigvar(2)= .TRUE.
                           input%eigvar(3)= .TRUE.
486
                            
487 488 489 490 491 492 493 494 495 496 497 498 499 500 501
                           input%eigvar(2)= input%eigvar(2).AND.( noco%l_soc.AND.(.NOT.noco%l_noco))
                           ! eigvar(1/2)= 1st/2nd var. ; eigvar(3)= calc density,etc

                           IF (noco%l_soc.AND.(.NOT.noco%l_noco)) THEN
                              input%evonly(1)= .FALSE.
                              input%evonly(2)= input%eonly
                           ELSE
                              input%evonly(1)= input%eonly
                              input%evonly(2)= .FALSE.
                           ENDIF

                           IF (input%eigvar(1).OR.input%eigvar(2)) THEN

                              IF (jij%l_J) THEN
                                 input%tkb=0.
502
#ifdef CPP_NEVER
503 504
                                 CALL jcoff(i_J,j_J,phn,mpi,atoms,atoms,&
                                            noco,jij)
505
#endif
506 507 508 509 510 511
                              END IF
                              IF (input%eigvar(1)) THEN

                                 CALL timestart("generation of hamiltonian and diagonalization (total)")
                                 ! WRITE(6,fmt='(A)') 'Starting 1st variation ...'
                                 CALL timestart("eigen")
512
                                 IF (mpi%irank==0) WRITE(*,"(a)",advance="no") "* Eigenvalue problem "
513 514 515
                                 CALL eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,&
                                            sym,kpts,dimension,vacuum,input,cell,enpara,banddos,noco,jij,oneD,hybrid,&
                                            it,eig_id, results)
516
                                 eig_idList(pc) = eig_id
517 518 519 520 521 522 523
                                 CALL timestop("eigen")
                                 !
                                 !                   add all contributions to total energy
                                 !
                                 IF(xcpot%icorr == icorr_pbe0.OR.xcpot%icorr == icorr_hse .OR. &
                                    xcpot%icorr == icorr_hf  .OR.xcpot%icorr == icorr_vhse) THEN
                                    DEALLOCATE( results%w_iks )
524
#ifdef CPP_MPI
525 526 527 528 529 530 531 532 533 534 535 536
                                    ! 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
537
#endif
538
                                 END IF ! xcpot%icorr == any hybrid
539

540 541 542 543
                                 ! for orbital decomposed calculations all kpts are n
                                 ! so transform the eig file such that each kpt is ro
                                 ! to its symmetry equivalent ones
                                 IF (banddos%dos.AND.banddos%ndir == -3) THEN
544
#ifdef CPP_NEVER
545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563
                                    ! This code section has to be reimplemented to enable
                                    ! calculations of the orbital decomposed DOS with
                                    ! eigenvectors up to this point only calculated for the
                                    ! IBZ. Here we reconstruct the eigenvectors for the full
                                    ! BZ from this subset.

                                    ! generate pointer from gpts at one kpt to a diffe
                                    ALLOCATE(kpts%pntgptd(3))
                                    CALL generate_pntgpt(dimension,obsolete,input,&
                                                         noco,kpts,cell,sym)

                                    ! rotate the eigenvectors to the symmetryequivalen
                                    CALL rotate_eig(kpts,dimension,atoms,&
                                                    sym,mpi)

                                    DEALLOCATE(kpts%pntgptd, kpts%pntgpt)
                                    ! this change is sufficient to modify fermie and c
                                    ! the enlarged kpt mesh
                                    kpts%nkpt = kpts%nkptf
564
#endif
565 566 567 568 569 570 571 572 573 574
                                 END IF
                              END IF ! (input%eigvar(1))
                              IF (input%eigvar(2)) THEN
                                 ! RS: open unit for SOC vectors for GW
                                 IF (noco%l_soc.AND.input%gw.EQ.2) THEN
                                    WRITE(6,'(A)') 'RS: open SOCVEC unit 4649'
                                    OPEN(4649,file='SOCVEC',form='unformatted')
                                 END IF

                                 ! WRITE(6,fmt='(A)') 'Starting 2nd variation ...'
575
                                 IF (mpi%irank==0) WRITE(*,"(a)",advance="no") "* SOC "
576 577 578 579 580 581 582 583 584 585 586 587 588 589
                                 CALL eigenso(eig_id,mpi,dimension,stars,vacuum,atoms,sphhar,&
                                              obsolete,sym,cell,noco,input,kpts, oneD)

                                 IF (noco%l_soc.AND.input%gw.EQ.2) THEN
                                    CLOSE(4649)
                                    INQUIRE(1014,opened=l_endit)
                                    IF (l_endit) CLOSE(1014)
                                    INQUIRE(667,opened=l_endit)
                                    IF (l_endit) CLOSE(667)
                                    CALL juDFT_end("GW+SOC finished",mpi%irank)
                                 END IF
                              END IF ! (input%eigvar(2))

                              CALL timestop("generation of hamiltonian and diagonalization (total)")
590 591

#ifdef CPP_MPI
592
                              CALL MPI_BARRIER(mpi%mpi_comm,ierr)
593 594
#endif

595
                           END IF ! (input%eigvar(1).OR.input%eigvar(2))
596

597
                           IF (input%eigvar(3).AND.noco%l_soc.AND.noco%l_ss) THEN
598
#ifdef CPP_NEVER
599 600
                              CALL ssomat(eig_id, mpi,dimension,stars,vacuum,atoms,sphhar,&
                                          sym,cell,noco,input,obsolete,kpts,oneD,MPI_DOUBLE_PRECISION)
601
#endif
602 603 604 605 606 607 608
                              stop80= .TRUE.
                           END IF

                           !-t3e
                           !
                           !              ----> fermi level and occupancies
                           !
609
                           IF (mpi%irank==0) WRITE(*,"(a)",advance="no") "* Fermi-energy "
610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638
                           IF (input%eigvar(3).AND.(.NOT.(noco%l_soc.AND.noco%l_ss))) THEN

                              IF (jij%l_J) THEN
                                 CALL timestart("determination of fermi energy")
                                 ALLOCATE (results%w_iks(dimension%neigd,kpts%nkpt,dimension%jspd))
                                 CALL fermie(eig_id, mpi,kpts,obsolete,input,&
                                             noco,enpara%epara_min,jij,cell,results)
                                 DEALLOCATE (results%w_iks)
                                 CALL timestop("determination of fermi energy")
                              END IF

                              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).AND.(.NOT.jij%l_J)) 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,&
                                                input,noco,enpara%epara_min,jij,cell,results)
                                    results%seigscv = results%seigscv/2
                                    results%ts = results%ts/2
                                    input%zelec = input%zelec/2
                                 ELSE
                                    CALL fermie(eig_id,mpi,kpts,obsolete,&
                                                input,noco,enpara%epara_min,jij,cell,results)
                                 END IF
                                 CALL timestop("determination of fermi energy")
639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657
!+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
658 659 660 661 662 663 664
                              END IF

                              IF (input%eonly) THEN
                                 CALL close_eig(eig_id)

                                 IF (.NOT.jij%l_J) THEN
                                    DEALLOCATE( results%w_iks )
665
#ifdef CPP_MPI
666
                                    CALL MPI_BARRIER(mpi%mpi_comm,ierr)
667
#endif
668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691
                                    IF (mpi%irank==0) THEN
                                       WRITE (*,fmt='(A)') 'eigenvalues written, program stops'
                                    END IF
                                    stop80= .TRUE.
                                 END IF
                              END IF ! input%eonly

                           END IF ! (input%eigvar(3).and..not.(noco%l_soc.and.noco%l_ss))

                           IF(jij%l_J) THEN
                              IF (.NOT. input%eonly) THEN
                                 DEALLOCATE ( results%w_iks )
                              END IF
                              IF (((i_J.EQ.j_J)).OR.(sym%invs.AND.(jij%qn.GT.tol))) GOTO 33
                           END IF

                        END DO !phn

33                      CONTINUE
                     END DO !j_J
                     i_J = i_J + jij%nmagtype(imt)

                  END DO !imt

692
                  IF ((input%l_wann.OR.jij%l_J).AND.(mpi%irank.EQ.0)) THEN
693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713
                     CALL timestop("Q-point for J_ij(total)")
                  END IF

               END DO !pc

            END DO !qcount

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

234         CONTINUE

            IF (mpi%irank.EQ.0) THEN
               IF (jij%l_J) THEN
                  IF (.NOT.jij%l_disp) THEN
                     REWIND(113)
                     REWIND(114)
714
#ifdef CPP_NEVER
715
                     CALL jcoff2(atoms,sym,cell,jij,input)
716
#endif
717 718 719 720 721 722 723 724
                  END IF
                  CLOSE(113)
                  CLOSE(114)
               END IF
            END IF

            IF (.NOT.jij%l_J) THEN
               !--- J>
725
#ifdef CPP_MPI
726 727
               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)
728
#endif
729 730 731 732 733 734 735 736
               !-t3e
               !
               !        ----> initialise force_old
               !
               CALL force_0(results)
               !
               !        ----> charge density

737
!+Wannier functions
738 739
               input%l_wann = .FALSE.
               INQUIRE (file='wann_inp',exist=input%l_wann)
740
               IF ((input%l_wann).AND.(.NOT.wann%l_bs_comf)) THEN
741 742 743
                  CALL wannier(DIMENSION,mpi,input,sym,atoms,stars,vacuum,sphhar,oneD,&
                               wann,noco,cell,enpara,banddos,sliceplot,results,&
                               eig_idList,l_real,kpts%nkpt)
744
               END IF
745
               IF (wann%l_gwf) CALL juDFT_error("provide wann_inp if l_gwf=T", calledby = "fleur")
746
!-Wannier
747 748

               CALL timestart("generation of new charge density (total)")
749
               IF (mpi%irank==0) WRITE(*,"(a)",advance="no") "* New Charge "
750 751 752 753 754 755 756 757 758 759 760 761 762 763 764
               CALL cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
                           dimension,kpts,atoms,sphhar,stars,sym,obsolete,&
                           enpara,cell,noco,jij,results,oneD)
               !
               ! the w_iks are needed for hybrid functionals so do not
               ! deallocate them in that case

               IF (xcpot%icorr /= icorr_pbe0.AND.xcpot%icorr /= icorr_hse.AND. &
                   xcpot%icorr /= icorr_hf  .AND.xcpot%icorr /= icorr_exx.AND. &
                   xcpot%icorr /= icorr_vhse ) THEN
                  DEALLOCATE ( results%w_iks )
               END IF

               IF (noco%l_soc.AND.(.NOT.noco%l_noco)) dimension%neigd=dimension%neigd/2
               !+t3e
765
#ifdef CPP_MPI
766 767 768
               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)
Daniel Wortmann's avatar
Daniel Wortmann committed
769
                
770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826
               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)
                     END IF
                  END DO
                  IF (noco%l_constr) THEN
                     CALL MPI_BCAST(noco%b_con,SIZE(noco%b_con),MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
                  END IF
               END IF
#endif
               CALL timestop("generation of new charge density (total)")
               IF (mpi%irank.EQ.0) THEN
                  !-t3e

                  IF (banddos%ndir.GT.0) THEN
                     CALL juDFT_end("NDIR",mpi%irank)
                  END IF

                  !          ----> output potential and potential difference
                  IF (obsolete%disp) THEN
                     reap = .FALSE.
                     input%total = .FALSE.
                     CALL timestart("generation of potential (total)")
                     CALL vgen(reap,input,xcpot,dimension, atoms,sphhar,stars,vacuum,sym,&
                               obsolete,cell,oneD,sliceplot,mpi, results,noco)
                     CALL timestop("generation of potential (total)")

                     CALL potdis(stars,vacuum,atoms,sphhar, input,cell,sym)
                  END IF
                  !
                  !i         ----> total energy
                  !


                  CALL timestart('determination of total energy')
                  CALL totale(atoms,sphhar,stars,vacuum,dimension,&
                              sym,input,noco,cell,oneD,xcpot,hybrid,it,results)
                  CALL timestop('determination of total energy')

                  ! in case of parallel processing, the total energy calculation is done
                  ! only for irank.eq.0, since no parallelization is required here. once
                  ! a force calculation is applied, however, the irank.eq.0 process is
                  ! led into a MPI_FINALIZE after convergence, while the other processes
                  ! are not, resulting in fleur not terminating despite having finished
                  ! the calculation. the next 7 lines correct that issue.
                  ! (other files subject to this correction: geo.F, force_w.F)
                  ! Schlipf/Klueppelberg Jun 2012
#ifdef CPP_MPI
               ELSE IF (input%l_f) THEN ! forces, but mpi%irank.ne.0
                  !This does not work, you can not call MPI_BCAST within a
                  !else part of irank==0 as PE=0 will not call this!
                  !CALL MPI_BCAST(lconv,1,MPI_LOGICAL,0,mpi_comm,ierr)
                  !IF (lconv) THEN
                  !  CALL MPI_FINALIZE(ierr)
                  !END IF
827
#endif
828 829 830 831 832 833
               END IF ! mpi%irank.EQ.0

               !Close file if not a hybrid calculation
               IF (xcpot%icorr /= icorr_pbe0.AND.xcpot%icorr /= icorr_hse .AND. &
                   xcpot%icorr /= icorr_hf  .AND.xcpot%icorr /= icorr_exx .AND. &
                   xcpot%icorr /= icorr_vhse) THEN
834 835 836
                  DO pc = 1, wann%nparampts
                     CALL close_eig(eig_idList(pc))
                  END DO
837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853
               END IF

            END IF !(if not jij%l_J)

         ELSE

            input%alpha = input%alpha - 10.0

         END IF !(if input%alpha <10.)

         IF (.NOT.jij%l_J) THEN
            IF (mpi%irank.EQ.0) THEN
               !-t3e
               !
               !          ----> mix input and output densities
               !
               CALL timestart("mixing")
854
               IF (mpi%irank==0) WRITE(*,"(a)",advance="no") "* Mixing distance: "
855 856 857 858 859
               CALL mix(stars,atoms,sphhar,vacuum,input,sym,cell,it,noco,oneD,hybrid,results)

               CALL timestop("mixing")
               WRITE (6,FMT=8130) it
               WRITE (16,FMT=8130) it
860 861 862 863
               IF (mpi%irank==0) THEN
                  WRITE(*,"(f11.7)",advance='no') results%last_distance
                  WRITE(*,*)
               ENDIF
864 865 866 867 868 869
8130           FORMAT (/,5x,'******* it=',i3,'  is completed********',/,/)
               CALL timestop("Iteration")
               !+t3e
            END IF ! mpi%irank.EQ.0

            ! hybrid functionals - ddist is needed on all processes
870
#ifdef CPP_MPI
871
            CALL MPI_BCAST(hybrid%ddist,dimension%jspd,MPI_REAL8,0,mpi%mpi_comm,ierr)
872 873
#endif

874 875 876 877
                !--- J<
         ELSE
         END IF !(if not jij%l_J)
         !--- J>
878 879

#ifdef CPP_MPI
880 881
         CALL MPI_BCAST(results%last_distance,1,MPI_DOUBLE_PRECISION,0,mpi%mpi_comm,ierr)
         CALL MPI_BARRIER(mpi%mpi_comm,ierr)
882
#endif
883 884 885 886 887 888 889 890 891 892
         !-t3e
         ! Delete the broyden files after the fifth iteration
         ! in the case of a HF or hybrid functional calculation
         IF ((it.EQ.5).AND.(xcpot%icorr.EQ.icorr_hf .OR.xcpot%icorr.EQ.icorr_pbe0.OR. &
                            xcpot%icorr.EQ.icorr_hse.OR.xcpot%icorr.EQ.icorr_vhse).AND. &
                           input%imix.LE.10) THEN
            CALL system('rm -f broyd*')
         END IF
         !+fo
         CALL priv_geo_end(mpi)
893
       
894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923
         !-fo
         IF (hybrid%l_calhf) ithf = ithf + 1
         IF (xcpot%icorr == icorr_hf .OR.xcpot%icorr == icorr_pbe0.OR.&
             xcpot%icorr == icorr_hse.OR.xcpot%icorr == icorr_vhse) THEN
            l_cont = ( ithf < input%itmax )
            results%te_hfex%core    = 0
            results%te_hfex%valence = 0
         ELSE
            l_cont = ( it < input%itmax )
         END IF

         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

      END DO ! DO WHILE (l_cont)

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

      CALL juDFT_end("all done",mpi%irank)

   CONTAINS

      SUBROUTINE priv_geo_end(mpi)

Daniel Wortmann's avatar
Daniel Wortmann committed
924
         TYPE(t_mpi),INTENT(IN)::mpi
925
         LOGICAL :: l_exist
926

927
         !Check if a new input was generated
928
         INQUIRE (file='inp_new', exist=l_exist)
929
         IF (l_exist) THEN
930
            CALL juDFT_end(" GEO new inp created ! ", mpi%irank)
931
         END IF
932

933
         !check for inp.xml
934
         INQUIRE (file='inp_new.xml', exist=l_exist)
935 936 937 938 939 940 941 942 943
         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)
944
               PRINT *,"qfix set to F"
945 946 947 948
            ENDIF
            INQUIRE(file='broyd',exist=l_exist)
            IF (l_exist) THEN
               CALL system('rm broyd')
949
               PRINT *,"broyd file deleted"
950 951 952 953
            ENDIF
            INQUIRE(file='broyd.7',exist=l_exist)
            IF (l_exist) THEN
               CALL system('rm broyd.7')
954 955 956 957 958
               PRINT *,"broyd.7 file deleted"
            END IF
         END IF
         CALL juDFT_end(" GEO new inp.xml created ! ", mpi%irank)
      END SUBROUTINE priv_geo_end
959
       
960 961
   END SUBROUTINE
END MODULE