eigen.F90 28.1 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_eigen
8
  USE m_juDFT
9 10
CONTAINS
  SUBROUTINE eigen(mpi,stars,sphhar,atoms,obsolete,xcpot,&
11
       sym,kpts,DIMENSION, vacuum, input, cell, enpara_in,banddos, noco,jij, oneD,hybrid,&
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
       it,eig_id,results)
    !*********************************************************************
    !     sets up and solves the eigenvalue problem for a basis of lapws.
    !
    ! nv,   nvd     ... actual length & dimension of EV without LO's
    ! nmat, nbasfcn                                   including LO's
    !        g. bihlmayer '96
    !**********************************************************************
    USE m_constants, ONLY : pi_const,sfp_const
    USE m_types
    USE m_lodpot
    USE m_tlmplm
    USE m_tlmplm_store
    USE m_apws
    USE m_hsmt
    USE m_hsint
    USE m_hsvac
    USE m_od_hsvac
    USE m_usetup
31
    USE m_pot_io
32
    USE m_eigen_diag
33 34
    USE m_eigen_hf_init
    USE m_eigen_hf_setup
35 36 37 38 39 40 41 42 43 44 45 46 47 48
#ifdef CPP_NOTIMPLEMENTED
    USE m_symm_hf,  ONLY : symm_hf_nkpt_EIBZ
    USE m_gen_bz
    USE m_gen_wavf
    USE m_hsfock
    USE m_read_core
    USE m_subvxc
    USE m_gweig
    USE m_gw_qsgw
    USE m_checkolap
#endif
    USE m_hsefunctional
    USE m_hybridmix    , ONLY: amix_pbe0,amix_hf
    USE m_util
49
    USE m_icorrkeys
50
    USE m_eig66_io, ONLY : open_eig, write_eig, close_eig,read_eig
51
    USE m_xmlOutput
52 53 54
#ifdef CPP_MPI
    USE m_mpi_bc_pot
#endif
55 56 57 58 59

    IMPLICIT NONE
    TYPE(t_results),INTENT(INOUT):: results
    TYPE(t_xcpot),INTENT(IN)     :: xcpot
    TYPE(t_mpi),INTENT(IN)       :: mpi
60
    TYPE(t_dimension),INTENT(IN) :: DIMENSION
61
    TYPE(t_oneD),INTENT(IN)      :: oneD
62
    TYPE(t_hybrid),INTENT(INOUT) :: hybrid
63 64 65 66 67
    TYPE(t_enpara),INTENT(INOUT) :: enpara_in
    TYPE(t_obsolete),INTENT(IN)  :: obsolete
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_vacuum),INTENT(IN)    :: vacuum
    TYPE(t_noco),INTENT(IN)      :: noco
68
    TYPE(t_banddos),INTENT(IN)   :: banddos
69
    TYPE(t_jij),INTENT(IN)       :: jij
Daniel Wortmann's avatar
Daniel Wortmann committed
70
    TYPE(t_sym),INTENT(IN)       :: sym  
71 72 73 74 75 76 77 78 79 80 81
    TYPE(t_stars),INTENT(IN)     :: stars
    TYPE(t_cell),INTENT(IN)      :: cell
    TYPE(t_kpts),INTENT(IN)      :: kpts
    TYPE(t_sphhar),INTENT(IN)    :: sphhar
    TYPE(t_atoms),INTENT(INOUT)  :: atoms !in u_setup n_u might be modified

#ifdef CPP_MPI
    INCLUDE 'mpif.h'
#endif
    !     ..
    !     .. Scalar Arguments ..
82
    INTEGER,INTENT(IN) :: it
83 84 85 86 87 88
    INTEGER,INTENT(INOUT):: eig_id
    !     ..
    !-odim
    !+odim
    !     ..
    !     .. Local Scalars ..
89
    INTEGER jsp,nk,nred,ne_all,n_u_in,ne_found
90 91 92
    INTEGER iter,ne,matsize  ,nrec,lh0
    INTEGER nspins,isp,l,i,j,err,gwc
    INTEGER mlotot,mlolotot,mlot_d,mlolot_d,nlot_d
Daniel Wortmann's avatar
Daniel Wortmann committed
93
    LOGICAL l_wu,lcal_qsgw,l_file,l_real,l_zref
94
    REAL evac_sv(DIMENSION%jspd)
95
    INTEGER ::eig_id_hf=-1
96
    
97 98 99 100
    !     ..
    !     .. Local Arrays ..
    INTEGER, PARAMETER :: lmaxb=3
    INTEGER, ALLOCATABLE :: matind(:,:),kveclo(:)
101
    INTEGER, ALLOCATABLE :: nv2(:)
102
    REAL,    ALLOCATABLE :: bkpt(:)
103
    REAL,    ALLOCATABLE :: eig(:)
104

105 106
    COMPLEX, ALLOCATABLE :: vpw(:,:),vzxy(:,:,:,:)
    COMPLEX, ALLOCATABLE :: vpwtot(:,:)
107
    REAL,    ALLOCATABLE :: vz(:,:,:),vr(:,:,:,:)
108 109 110
    REAL,    ALLOCATABLE :: vrtot(:,:,:,:)

    COMPLEX, ALLOCATABLE :: vs_mmp(:,:,:,:)
111 112 113
    TYPE(t_tlmplm)  :: td
    TYPE(t_usdus)   :: ud
    TYPE(t_lapw)    :: lapw
114
    TYPE(t_enpara)  :: enpara
115 116
    TYPE(t_zMat)    :: zMat
    TYPE(t_hamOvlp) :: hamOvlp
117
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
118
    INTEGER fh,nn,n
119 120 121 122 123 124 125 126
    INTEGER ierr(3)

    !
    !     .. variables for HF or hybrid functional calculation ..
    !
    !      - scalar -
#ifdef CPP_NEVER
    INTEGER, INTENT(IN)     ::  maxlcutm,maxindxm,maxbasm
127 128
    INTEGER, INTENT(IN)     ::  maxindxp
    INTEGER, INTENT(IN)     ::  bands
129
    !     - arrays -
Daniel Wortmann's avatar
Daniel Wortmann committed
130 131 132
    INTEGER, INTENT(IN)     ::  nindxm(0:maxlcutm,atoms%ntype)
    INTEGER, INTENT(IN)     ::  lcutm(atoms%ntype)
    REAL   , INTENT(IN)     ::  basm(atoms%jmtd,maxindxm,0:maxlcutm,atoms%ntype)
133 134 135 136 137 138 139 140 141 142
#endif
    !     - local scalar -
    INTEGER                 ::  itype,ispin,isym,iisym
    INTEGER                 ::  indx,ic
    INTEGER                 ::  ll,lm,l1,l2
    INTEGER                 ::  lmaxcd
    INTEGER                 ::  maxindxc,mnobd
    INTEGER                 ::  maxfac
    INTEGER                 ::  maxbands
    !     - local arrays -
143 144 145
    TYPE(t_hybdat)   :: hybdat
    INTEGER                 ::  comm(kpts%nkpt),irank2(kpts%nkpt),isize2(kpts%nkpt)
  
146 147
#ifdef CPP_NEVER
    INTEGER                 ::  nobd(kpts%nkptf)
Daniel Wortmann's avatar
Daniel Wortmann committed
148
    INTEGER                 ::  lmaxc(atoms%ntype)
149
    INTEGER                 ::  g(3)
Daniel Wortmann's avatar
Daniel Wortmann committed
150
    INTEGER                 ::  nindxp(0:maxlcutm,atoms%ntype)
151
    INTEGER , ALLOCATABLE   ::  nkpt_EIBZ(:)
152 153 154
    INTEGER , ALLOCATABLE   ::  nindxc(:,:)
    INTEGER , ALLOCATABLE   ::  kveclo_eig(:,:)
    INTEGER , ALLOCATABLE   ::  nbasm(:)
Daniel Wortmann's avatar
Daniel Wortmann committed
155 156 157
    REAL                    ::  el_eig(0:atoms%lmaxd,atoms%ntype), ello_eig(atoms%nlod,atoms%ntype),rarr(3)
    REAL                    ::  bas1_MT(hybrid%maxindx,0:atoms%lmaxd,atoms%ntype)
    REAL                    ::  drbas1_MT(hybrid%maxindx,0:atoms%lmaxd,atoms%ntype)
158 159 160 161 162 163 164
    REAL,    ALLOCATABLE    ::  eig_c(:,:,:)
    REAL,    ALLOCATABLE    ::  core1(:,:,:,:),core2(:,:,:,:)
    REAL,    ALLOCATABLE    ::  gauntarr(:,:,:,:,:,:)
    REAL,    ALLOCATABLE    ::  sfac(:),fac(:)
    REAL,    ALLOCATABLE    ::  prodm(:,:,:,:)
    TYPE(PRODTYPE),ALLOCATABLE :: prod(:,:,:)
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
165
    INTEGER                 ::  ne_eig(kpts%nkpt),nbands(kpts%nkpt)
166
    REAL,    ALLOCATABLE    ::  eig_irr(:,:),vr0(:,:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
167 168
    REAL                    ::  bas1(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype)
    REAL                    ::  bas2(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype)
169 170

#ifdef CPP_MPI
Daniel Wortmann's avatar
Daniel Wortmann committed
171
    INTEGER   :: sndreqd,sndreq(mpi%isize*kpts%nkpt)
172 173 174
#endif
    !
    !
175
    ! --> Allocate
176
    !
177 178 179 180 181 182
    ALLOCATE ( ud%uloulopn(atoms%nlod,atoms%nlod,atoms%ntype,DIMENSION%jspd),nv2(DIMENSION%jspd) )
    ALLOCATE ( ud%ddn(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd),eig(DIMENSION%neigd),bkpt(3) )
    ALLOCATE ( ud%us(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd),ud%uds(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd) )
    ALLOCATE ( ud%dus(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd),ud%duds(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd))
    ALLOCATE ( ud%ulos(atoms%nlod,atoms%ntype,DIMENSION%jspd),ud%dulos(atoms%nlod,atoms%ntype,DIMENSION%jspd) )
    ALLOCATE ( ud%uulon(atoms%nlod,atoms%ntype,DIMENSION%jspd),ud%dulon(atoms%nlod,atoms%ntype,DIMENSION%jspd) )
Daniel Wortmann's avatar
Daniel Wortmann committed
183 184
   ! ALLOCATE ( enpara%ello(atoms%nlod,atoms%ntype,dimension%jspd) )
   ! ALLOCATE ( enpara%el(0:atoms%lmaxd,atoms%ntype,dimension%jspd),enpara%evac(2,dimension%jspd) )
185
    ALLOCATE ( lapw%k1(DIMENSION%nvd,DIMENSION%jspd),lapw%k2(DIMENSION%nvd,DIMENSION%jspd),lapw%k3(DIMENSION%nvd,DIMENSION%jspd),lapw%rk(DIMENSION%nvd,DIMENSION%jspd) )
186 187 188 189 190 191
    !
    ! --> some parameters first
    !
    !     determine the total number of lo's : nlotot
    !
    mlotot = 0 ; mlolotot = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
192 193 194
    DO nn = 1, atoms%ntype
       mlotot = mlotot + atoms%nlo(nn)
       mlolotot = mlolotot + atoms%nlo(nn)*(atoms%nlo(nn)+1)/2
195 196 197 198 199 200
    ENDDO
    nlot_d = atoms%nlotot !max(atoms%nlotot,1)
    ALLOCATE ( kveclo(nlot_d) )
    !     ..
    nbands     = 0
    bas1 = 0 ; bas2 = 0
201
    hybrid%l_hybrid   = (&
202 203 204 205 206
         xcpot%icorr == icorr_pbe0 .OR.&
         xcpot%icorr == icorr_hse  .OR.&
         xcpot%icorr == icorr_vhse .OR.&
         xcpot%icorr == icorr_hf   .OR.&
         xcpot%icorr == icorr_exx)
207
    l_real=sym%invs.AND..NOT.noco%l_noco
208
    IF (noco%l_soc.AND.l_real.AND.hybrid%l_hybrid ) THEN
209
       CALL juDFT_error('hybrid functional + SOC + inv.symmetry is not tested', calledby='eigen')
210
    END IF
211

212 213 214 215 216 217 218 219 220 221 222 223 224 225
    !
    !  if gw = 1 or 2, we are in the first or second run of a GW  calculation
    !  if gw = 1 we just proceed as normal (one round),
    !  if gw = 2 it's the second run: write out the eigenfunctions and
    !  the matrix elements with the xc-potential (needs two rounds)
    !  if gw = 3 energy-independet hermitian self-energy is read in from file
    !            spex.qsgw, transformed to the APW basis, and SCF is performed
    !
    gwc = 1
    fh = 0
    !
    ! look, if WU diagonalisation
    !
    IF (it.LT.input%isec1) THEN
226
       IF (mpi%irank.EQ.0) WRITE (6,FMT=8110) it,input%isec1
227
8110   FORMAT (' IT=',i4,'  ISEC1=',i4,' standard diagonalization')
228
       l_wu = .FALSE.
229
    ELSE
230
       IF (mpi%irank.EQ.0) WRITE (6,FMT=8120) it,input%isec1
231
8120   FORMAT (' IT=',i4,'  ISEC1=',i4,' reduced diagonalization')
232
       l_wu = .TRUE.
233 234 235 236
    END IF
    !
    ! load potential from file pottot (=unit 8)
    !
237 238 239 240 241
    ALLOCATE ( vpw(stars%ng3,DIMENSION%jspd),vzxy(vacuum%nmzxyd,oneD%odi%n2d-1,2,DIMENSION%jspd) )
    ALLOCATE ( vz(vacuum%nmzd,2,4), vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd) )
    ALLOCATE ( vr0(atoms%jmtd,atoms%ntype,DIMENSION%jspd) ) ; vr0 = 0
    IF (input%gw.EQ.2) THEN
       ALLOCATE ( vpwtot(stars%ng3,DIMENSION%jspd), vrtot(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,DIMENSION%jspd) )
242 243 244 245 246
       IF ( mpi%irank == 0 ) WRITE(6,'(A/A/A/A)')&
            &  'Info: vxc matrix elements for GW will be calculated in gw_vxc',&
            &  'Info: and stored in "vxc", the values obtained from the',&
            &  'Info: original implementation are saved to "vxc.old".'
    ENDIF
247 248 249 250 251 252 253 254 255 256

    IF (mpi%irank.EQ.0) THEN
       CALL readPotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_TOT_const,&
                          iter,vr,vpw,vz,vzxy)
    END IF
#ifdef CPP_MPI
    CALL mpi_bc_pot(mpi,stars,sphhar,atoms,input,vacuum,&
                    iter,vr,vpw,vz,vzxy)
#endif

257
999 CONTINUE
258
    IF (mpi%irank.EQ.0) CALL openXMLElementFormPoly('iteration',(/'numberForCurrentRun','overallNumber      '/),(/it,iter/),&
259
                                                    RESHAPE((/19,13,5,5/),(/2,2/)))
260 261 262 263

    !
    ! some modifications for gw-calculations
    !
264
    IF (input%gw.EQ.2.AND.gwc.EQ.1) THEN
265 266 267 268 269 270 271 272 273 274 275
       vrtot(:,:,:,:)  = vr  ! store potential for subroutine gw_vxc
       vpwtot(:,:) = vpw !
    ENDIF

    IF (gwc==1) THEN
       vr0(:,:,:) = vr(:,0,:,:)
       lh0 = 1
    ELSE IF (gwc==2) THEN
       lh0 = 0                         ! for a input%gw-calculation, we
                                       ! now evaluate matrix elements
       DO jsp = 1,input%jspins               ! with the coulomb potential
Daniel Wortmann's avatar
Daniel Wortmann committed
276 277 278
          DO nn = 1,atoms%ntype                ! but with explicit kinetic energy
             DO j = 1,atoms%jri(nn)
                vr(j,0,nn,jsp) = vr(j,0,nn,jsp)-vr0(j,nn,jsp)*sfp_const/atoms%rmsh(j,nn)
279 280 281 282 283 284
             ENDDO
          ENDDO
       ENDDO
    ENDIF

    INQUIRE(file='fleur.qsgw',EXIST=lcal_qsgw)
285
    lcal_qsgw = .NOT. lcal_qsgw
286 287 288 289 290 291 292 293 294

    !
    ! set energy parameters (normally to that, what we read in)
    !
    IF (gwc /= 2) THEN
       CALL lodpot(mpi,atoms,sphhar,obsolete,vacuum,&
            input, vr,vz, enpara_in, enpara)
    ENDIF
    !
295
   
296
    CALL eigen_hf_init(hybrid,kpts,sym,atoms,input,dimension,hybdat,irank2,isize2)
297 298 299 300

    !---> set up and solve the eigenvalue problem
    !---> loop over energy windows

Daniel Wortmann's avatar
Daniel Wortmann committed
301
!check if z-reflection trick can be used
302

303
    l_zref=(sym%zrfs.AND.(SUM(ABS(kpts%bk(3,:kpts%nkpt))).LT.1e-9).AND..NOT.noco%l_noco) 
304

Daniel Wortmann's avatar
Daniel Wortmann committed
305 306

#if ( defined(CPP_MPI))
307
    IF (mpi%n_size > 1) l_zref = .FALSE.
308
    IF ( hybrid%l_calhf ) THEN
309
       CALL judft_error("BUG parallelization in HF case must be fixed")
Daniel Wortmann's avatar
Daniel Wortmann committed
310 311
       !n_start  = 1
       !n_stride = 1
312 313
    END IF
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
314 315
    !Count number of matrix columns on this PE
    n=0
316
    DO i=1+mpi%n_rank,DIMENSION%nbasfcn,mpi%n_size
Daniel Wortmann's avatar
Daniel Wortmann committed
317
       n=n+1
318
    ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
319
    IF (mpi%n_size>1) THEN
320
       matsize = DIMENSION%nbasfcn * n
Daniel Wortmann's avatar
Daniel Wortmann committed
321
    ELSE
322 323 324 325 326
       IF (MOD(DIMENSION%nbasfcn,2)==0 ) THEN !same formula, division by 2 different to avail int overflow
          matsize = (DIMENSION%nbasfcn+1)*(DIMENSION%nbasfcn/2)
       ELSE
          matsize = ((DIMENSION%nbasfcn+1)/2)*DIMENSION%nbasfcn
       ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
327
    ENDIF
328 329
    IF (matsize<2) CALL judft_error("Wrong size of matrix",calledby="eigen",hint="Your basis might be too large or the parallelization fail or ??")
    ne = MAX(5,DIMENSION%neigd)
Daniel Wortmann's avatar
Daniel Wortmann committed
330

331
    IF (hybrid%l_hybrid.OR.hybrid%l_calhf) THEN
332
       eig_id_hf=eig_id
333
    ENDIF
334
    eig_id=open_eig(&
335 336 337 338
         mpi%mpi_comm,DIMENSION%nbasfcn,DIMENSION%neigd,kpts%nkpt,DIMENSION%jspd,atoms%lmaxd,&
         atoms%nlod,atoms%ntype,atoms%nlotot,noco%l_noco,.TRUE.,l_real,noco%l_soc,.FALSE.,&
         mpi%n_size,layers=vacuum%layers,nstars=vacuum%nstars,ncored=DIMENSION%nstd,&
         nsld=atoms%nat,nat=atoms%nat,l_dos=banddos%dos.OR.input%cdinf,l_mcd=banddos%l_mcd,&
339
         l_orb=banddos%l_orb)
340 341

    IF (l_real) THEN
342
       ALLOCATE ( hamOvlp%a_r(matsize), stat = err )
343
    ELSE
344
       ALLOCATE ( hamOvlp%a_c(matsize), stat = err )
345
    ENDIF
346 347 348 349 350
    IF (err.NE.0) THEN
       WRITE (*,*) 'eigen: an error occured during allocation of'
       WRITE (*,*) 'the Hamilton Matrix: ',err,'  size: ',matsize
       CALL juDFT_error("eigen: Error during allocation of Hamilton" //"matrix",calledby ="eigen")
    ENDIF
351
    IF (l_real) THEN
352
       ALLOCATE ( hamOvlp%b_r(matsize), stat = err )
353
    ELSE
354
       ALLOCATE ( hamOvlp%b_c(matsize), stat = err )
355
    ENDIF
356

357 358 359 360 361
    IF (err.NE.0) THEN
       WRITE (*,*) 'eigen: an error occured during allocation of'
       WRITE (*,*) 'the overlap Matrix: ',err,'  size: ',matsize
       CALL juDFT_error("eigen: Error during allocation of overlap " //"matrix",calledby ="eigen")
    ENDIF
362 363 364 365

    hamOvlp%l_real = l_real
    hamOvlp%matsize = matsize

366
    ALLOCATE (  matind(DIMENSION%nbasfcn,2) )
367 368 369 370 371 372
    !
    !--->    loop over spins
    nspins = input%jspins
    IF (noco%l_noco) nspins = 1
    !
    !        Append information about file eig to gwa
373
    IF(input%gw.EQ.2.AND.gwc.EQ.1) THEN
374 375 376 377 378
       IF ( mpi%irank == 0 ) THEN
          OPEN(15,file='gwa',status='old',form='unformatted')
          READ(15)
          READ(15)
          READ(15)
379 380
          WRITE(15) mpi%n_start,mpi%n_stride,mpi%n_rank,mpi%n_size,DIMENSION%nvd,&
               &                 DIMENSION%nbasfcn,atoms%nlotot
381 382 383 384 385
          CLOSE(15)
       END IF
    ENDIF
    !  ..
    !  LDA+U
386
    n_u_in=atoms%n_u
387 388
    IF ((atoms%n_u.GT.0)) THEN
       ALLOCATE( vs_mmp(-lmaxb:lmaxb,-lmaxb:lmaxb,atoms%n_u,input%jspins) )
389
       CALL u_setup(sym,atoms,lmaxb,sphhar,input, enpara%el0(0:,:,:),vr,mpi, vs_mmp,results)
390 391 392 393 394 395 396
    ELSE
       ALLOCATE( vs_mmp(-lmaxb:-lmaxb,-lmaxb:-lmaxb,1,2) )
    ENDIF
    !
    !--->    loop over k-points: each can be a separate task

    DO jsp = 1,nspins
397 398
       CALL eigen_HF_setup(hybrid,input,sym,kpts,dimension,atoms,mpi,noco,cell,oneD,results,jsp,eig_id_hf,&
         hybdat,irank2,it,vr0)  
399 400 401 402 403 404 405

       !
       !--->       set up k-point independent t(l'm',lm) matrices
       !
       CALL timestart("tlmplm")
       err=0
       j = 1 ; IF (noco%l_noco) j = 2
406 407 408 409 410 411 412
       ALLOCATE(td%tuu(0:DIMENSION%lmplmd,atoms%ntype,j),stat=err)
       ALLOCATE(td%tud(0:DIMENSION%lmplmd,atoms%ntype,j),stat=err)
       ALLOCATE(td%tdd(0:DIMENSION%lmplmd,atoms%ntype,j),stat=err)
       ALLOCATE(td%tdu(0:DIMENSION%lmplmd,atoms%ntype,j),stat=err)
       mlot_d = MAX(mlotot,1) ; mlolot_d = MAX(mlolotot,1)
       ALLOCATE(td%tdulo(0:DIMENSION%lmd,-atoms%llod:atoms%llod,mlot_d,j),stat=err)
       ALLOCATE(td%tuulo(0:DIMENSION%lmd,-atoms%llod:atoms%llod,mlot_d,j),stat=err)
413
       ALLOCATE(td%tuloulo(-atoms%llod:atoms%llod,-atoms%llod:atoms%llod,mlolot_d,j), stat=err)
414
       ALLOCATE(td%ind(0:DIMENSION%lmd,0:DIMENSION%lmd,atoms%ntype,j),stat=err )
415 416 417 418 419
       IF (err.NE.0) THEN
          WRITE (*,*) 'eigen: an error occured during allocation of'
          WRITE (*,*) 'the tlmplm%tuu, tlmplm%tdd etc.: ',err,'  size: ',mlotot
          CALL juDFT_error("eigen: Error during allocation of tlmplm, tdd  etc.",calledby ="eigen")
       ENDIF
420 421
       CALL tlmplm(sphhar,atoms,DIMENSION,enpara, jsp,1,mpi, vr(1,0,1,jsp),gwc,lh0,input, td,ud)
       IF (input%l_f) CALL write_tlmplm(td,vs_mmp,atoms%n_u>0,1,jsp,input%jspins)
422 423 424 425 426 427 428 429 430
       CALL timestop("tlmplm")

       !---> pk non-collinear
       !--->       call tlmplm again for the second spin direction in
       !--->       each MT, because the t-matrices are needed for both
       !--->       spins at once in hsmt
       IF (noco%l_noco) THEN
          isp = 2
          CALL timestart("tlmplm")
431 432
          CALL tlmplm(sphhar,atoms,DIMENSION,enpara,isp,isp,mpi, vr(1,0,1,isp),gwc,lh0,input, td,ud)
          IF (input%l_f) CALL write_tlmplm(td,vs_mmp,atoms%n_u>0,2,2,input%jspins)
433 434 435 436
          CALL timestop("tlmplm")
       ENDIF
       !

437
#ifdef CPP_MPI
438 439
       ! check that all sending operations are completed
       IF ( hybrid%l_calhf ) CALL MPI_WAITALL(sndreqd,sndreq,MPI_STATUSES_IGNORE,ierr)
440
#endif
441

Daniel Wortmann's avatar
Daniel Wortmann committed
442
       k_loop:DO nk = mpi%n_start,kpts%nkpt,mpi%n_stride
443
#if defined(CPP_MPI)&&defined(CPP_NEVER)
444 445 446 447 448 449
          IF ( hybrid%l_calhf ) THEN
             ! jump to next k-point if this process is not present in communicator
             IF ( comm(nk) == MPI_COMM_NULL ) CYCLE
             ! allocate buffer for communication of the results
             IF ( irank2(nk) /= 0 ) CALL work_dist_reserve_buffer( nbands(nk) )
          END IF
450
#endif
451 452

          nrec =  kpts%nkpt*(jsp-1) + nk
Daniel Wortmann's avatar
Daniel Wortmann committed
453
          nrec = mpi%n_size*(nrec-1) + mpi%n_rank + 1
454 455 456
          !
          !--->         set up lapw list
          !
457
          CALL timestart("Setup of LAPW")
458
          lapw%rk = 0 ; lapw%k1 = 0 ; lapw%k2 = 0 ; lapw%k3 = 0
459
          CALL apws(DIMENSION,input,noco, kpts,nk,cell,l_zref, mpi%n_size,jsp, bkpt,lapw,matind,nred)
460

461
          CALL timestop("Setup of LAPW")
462 463 464 465 466 467 468 469
          IF (noco%l_noco) THEN
             !--->         the file potmat contains the 2x2 matrix-potential in
             !--->         the interstitial region and the vacuum
             OPEN (25,FILE='potmat',FORM='unformatted', STATUS='old')
          ENDIF
          !
          !--->         set up interstitial hamiltonian and overlap matrices
          !
470
          CALL timestart("Interstitial Hamiltonian&Overlap")
Daniel Wortmann's avatar
Daniel Wortmann committed
471
          CALL hsint(input,noco,jij,stars, vpw(:,jsp),lapw,jsp, mpi%n_size,mpi%n_rank,kpts%bk(:,nk),cell,atoms,l_real,hamOvlp)
472

473
          CALL timestop("Interstitial Hamiltonian&Overlap")
474 475 476
          !
          !--->         update with sphere terms
          !
477 478 479
          IF (.NOT.l_wu) THEN
             CALL timestart("MT Hamiltonian&Overlap")
             CALL hsmt(DIMENSION,atoms,sphhar,sym,enpara, mpi%SUB_COMM,mpi%n_size,mpi%n_rank,jsp,input,mpi,&
480
                  lmaxb,gwc, noco,cell, lapw, bkpt,vr, vs_mmp, oneD,ud, kveclo,td,l_real,hamOvlp)
481
             CALL timestop("MT Hamiltonian&Overlap")
482 483 484
          ENDIF
          !
#ifdef CPP_NOTIMPLEMENTED
485
          IF( hybrid%l_hybrid ) THEN
486

487
             CALL hsfock(nk,atoms,lcutm,obsolete,lapw, DIMENSION,kpts,jsp,input,hybrid,maxbasm,&
488
                  maxindxp,maxlcutm,maxindxm,nindxm, basm,bas1,bas2,bas1_MT,drbas1_MT,ne_eig,eig_irr,&
Daniel Wortmann's avatar
Daniel Wortmann committed
489
                  mpi%n_size,sym,cell, noco,noco,oneD, nbasp,nbasm, results,results,it,nbands(nk),maxbands,nobd,&
490 491 492 493 494 495
                  mnobd,xcpot, core1,core2,nindxc,maxindxc,lmaxc, lmaxcd, kveclo_eig,maxfac,fac,sfac,gauntarr,&
                  nindxp,prod,prodm,gwc, mpi,irank2(nk),isize2(nk),comm(nk), a)

             IF ( irank2(nk) /= 0 ) CYCLE

             IF( hybrid%l_subvxc ) THEN
496
                CALL subvxc(lapw,kpts(:,nk),obsolete,DIMENSION, input,jsp,atoms, hybrid,matsize,enpara%el0,enpara%ello0,&
497 498 499 500
                     sym, nlot_d,kveclo, cell,sphhar, stars,stars, xcpot,mpi, irank2(nk),vacuum,&
                     oneD, vr(:,:,:,jsp),vpw(:,jsp), a)
             END IF

501
          END IF ! hybrid%l_hybrid
502 503 504 505
#endif
          !
          !--->         update with vacuum terms
          !
506
          CALL timestart("Vacuum Hamiltonian&Overlap")
507
          IF (input%film .AND. .NOT.oneD%odi%d1) THEN
508
             CALL hsvac(vacuum,stars,DIMENSION, atoms, jsp,input,vzxy(1,1,1,jsp),vz,enpara%evac0,cell, &
Daniel Wortmann's avatar
Daniel Wortmann committed
509
                  bkpt,lapw,sym, noco,jij, mpi%n_size,mpi%n_rank,nv2,l_real,hamOvlp)
510
          ELSEIF (oneD%odi%d1) THEN
511
             CALL od_hsvac(vacuum,stars,DIMENSION, oneD,atoms, jsp,input,vzxy(1,1,1,jsp),vz, &
512
                  enpara%evac0,cell, bkpt,lapw, oneD%odi%M,oneD%odi%mb,oneD%odi%m_cyl,oneD%odi%n2d, &
Daniel Wortmann's avatar
Daniel Wortmann committed
513
                  mpi%n_size,mpi%n_rank,sym,noco,jij,nv2,l_real,hamOvlp)
514
          END IF
515
          CALL timestop("Vacuum Hamiltonian&Overlap")
516 517

#ifdef CPP_NOTIMPLEMENTED
518
          IF ( input%gw.EQ.3.OR.(input%gw.EQ.2.AND.gwc.EQ.1.AND..NOT.lcal_qsgw)) THEN
519 520

             CALL gw_qsgw ( lcal_qsgw, b,cell,sym,atoms,&
521
                  jsp,DIMENSION,lapw, nk,kpts, matsize,oneD%tau,noco, a )
522 523 524 525 526


          END IF

          IF (gwc==2) THEN
527
             CALL gw_eig(eig_id,nk,kpts,atoms,DIMENSION,neigd,sym,&
528 529 530 531 532 533 534 535 536 537
                  kveclo,cell, ud%us(0,1,jsp),ud%dus(0,1,jsp),ud%uds(0,1,jsp),&
                  ud%duds(0,1,jsp),ud%ddn(0,1,jsp),ud%ulos(1,1,jsp),ud%uulon(1,1,jsp),ud%dulon(1,1,jsp),&
                  ud%dulos(1,1,jsp),nrec,noco,jsp,matsize,a,sphhar,stars,stars,&
                  vrtot(1,0,1,jsp),vpwtot,vr,vpw,vs_mmp(-lmaxb,-lmaxb,1,jsp),lmaxb,oneD)
             CYCLE k_loop
          ENDIF
#endif
          IF (noco%l_noco) CLOSE (25)

          !write overlap matrix b to direct access file olap
538 539 540
          INQUIRE(file='olap',exist=l_file)
          IF (l_file) THEN
             IF (l_real) THEN
541
                OPEN(88,file='olap',form='unformatted',access='direct', recl=matsize*8)
542
                WRITE(88,rec=nrec) hamOvlp%b_r
543
                CLOSE(88)
544
             ELSE
545
                OPEN(88,file='olap',form='unformatted',access='direct', recl=matsize*16)
546
                WRITE(88,rec=nrec) hamOvlp%b_c
547
                CLOSE(88)
548 549
             ENDIF
          ENDIF
550

551
       
552
          CALL eigen_diag(jsp,eig_id,it,atoms,DIMENSION,matsize,mpi, mpi%n_rank,mpi%n_size,ne,nk,lapw,input,&
Daniel Wortmann's avatar
Daniel Wortmann committed
553
               nred,mpi%sub_comm, sym,l_zref,matind,kveclo, noco,cell,bkpt,enpara%el0,jij,l_wu,&
554
               oneD,td,ud, eig,ne_found,hamOvlp,zMat)
555
          
556 557 558 559
          !
          !--->         output results
          !
          CALL timestart("EV output")
560
          ne_all=ne_found
561 562
#if defined(CPP_MPI)
          !Collect number of all eigenvalues
Daniel Wortmann's avatar
Daniel Wortmann committed
563
          CALL MPI_ALLREDUCE(ne_found,ne_all,1,MPI_INTEGER,MPI_SUM, mpi%sub_comm,ierr)
564
          ne_all=MIN(DIMENSION%neigd,ne_all)
565 566
#endif
          !jij%eig_l = 0.0 ! need not be used, if hdf-file is present
567 568 569
          IF (.NOT.l_real) THEN
             IF (.NOT.jij%l_J) THEN
                zMat%z_c(:lapw%nmat,:ne_found) = CONJG(zMat%z_c(:lapw%nmat,:ne_found))
570
             ELSE
571
                zMat%z_c(:lapw%nmat,:ne_found) = CMPLX(0.0,0.0)
572
             ENDIF
573
          ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
574
	  zmat%nbands=ne_found
575
          CALL write_eig(eig_id, nk,jsp,ne_found,ne_all,lapw%nv(jsp),lapw%nmat,&
576 577
                  lapw%k1(:lapw%nv(jsp),jsp),lapw%k2 (:lapw%nv(jsp),jsp),lapw%k3(:lapw%nv(jsp),jsp),&
                  bkpt, kpts%wtkpt(nk),eig(:ne_found),enpara%el0(0:,:,jsp), enpara%ello0(:,:,jsp),enpara%evac0(:,jsp),&
578
                  atoms%nlotot,kveclo,mpi%n_size,mpi%n_rank,zMat)
579
          IF (noco%l_noco) THEN
580
             CALL write_eig(eig_id, nk,2,ne_found,ne_all,lapw%nv(2),lapw%nmat,&
581
                  lapw%k1(:lapw%nv(2),2),lapw%k2 (:lapw%nv(2),2),lapw%k3(:lapw%nv(2),2),&
582
                  bkpt, kpts%wtkpt(nk),eig(:ne_found),enpara%el0(0:,:,2), enpara%ello0(:,:,2),enpara%evac0(:,2),&
583 584
                  atoms%nlotot,kveclo)
          ENDIF
585 586 587 588
#if defined(CPP_MPI)
          !RMA synchronization
          CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
#endif
589

590
#if defined(CPP_MPI)&&defined(CPP_NEVER)
591 592 593 594 595 596 597 598 599 600 601 602 603 604 605
          IF ( hybrid%l_calhf ) THEN
             IF ( isize2(nk) == 1 ) THEN
                WRITE(*,'(a,i6,a,i6,a)') 'HF: kpt ', nk, ' was done by rank ', mpi%irank, '.'
             ELSE
                WRITE(*,'(a,i6,a,i6,a,i6,a)')&
                     'HF: kpt ', nk, ' was done by rank ', mpi%irank, ' and ', isize2(nk)-1, ' more.'
             END IF
             !                ELSE
             !                  WRITE(*,'(a,i6,a,i6,a)')        '    kpt ', nk, ' was done by rank ', irank, '.'
          END IF
          !#             else
          !                WRITE (*,*) 'pe: ',irank,' wrote ',nrec
#             endif
          CALL timestop("EV output")
          !#ifdef CPP_MPI
606
          IF (l_real) THEN
607
             DEALLOCATE ( zMat%z_r )
608
          ELSE
609
             DEALLOCATE ( zMat%z_c )
610
ENDIF
611 612 613 614 615 616 617 618 619 620 621 622 623
          !
       END DO  k_loop

       DEALLOCATE (td%tuu,td%tud,td%tdu,td%tdd)
       DEALLOCATE (td%ind,td%tuulo,td%tdulo)
       DEALLOCATE (td%tuloulo)
#ifdef CPP_NEVER
       IF ( hybrid%l_calhf ) THEN
          DEALLOCATE ( eig_irr,kveclo_eig )
       END IF
#endif
    END DO ! spin loop ends
    DEALLOCATE( vs_mmp )
624
    DEALLOCATE (matind)
625 626 627 628 629
    IF (l_real) THEN
       DEALLOCATE(hamOvlp%a_r,hamOvlp%b_r)
    ELSE
       DEALLOCATE(hamOvlp%a_c,hamOvlp%b_c)
    ENDIF
630 631 632 633 634 635
#ifdef CPP_NEVER
    IF( hybrid%l_calhf ) THEN
       DEALLOCATE( fac,sfac,gauntarr )
       DEALLOCATE( nindxc,core1,core2,nbasm,eig_c )
    END IF
#endif
636
#if defined(CPP_MPI)&&defined(CPP_NEVER)
637
    IF ( hybrid%l_calhf ) DEALLOCATE (nkpt_EIBZ)
638
#endif
639

640
    IF ( input%gw.EQ.2.AND.(gwc==1) )  THEN        ! go for another round
641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664
       !
       !       Generate input file abcoeff for subsequent GW calculation
       !       28.10.2003 Arno Schindlmayr
       !
       IF ( mpi%irank == 0 ) THEN
          WRITE(6,'(A)') 'Info: Write out vxc for GW and vxc.old.'
          WRITE(6,'(A)') 'Info: Write out abcoeff for GW.'
          WRITE(6,'(A)') 'Info: Write out radfun for gw_vxc and GW.'
       END IF
       OPEN (12,file='vxc.old',form='formatted',status='unknown') ! contains vxc from gw_eig
       OPEN (13,file='vxc',form='formatted',status='unknown')     ! contains vxc from gw_vxc
       OPEN (1013,file='vxcfull',form='unformatted',status='unknown')
       INQUIRE(file='fleur.qsgw',exist=l_file)
       IF(l_file) THEN
          WRITE(6,'(A)') 'Info: Write file qsgw for GW.'
          OPEN(1014,file='qsgw',form='unformatted')
       ENDIF
       OPEN (15,file='abcoeff',form='unformatted',status='unknown', action='write')
       OPEN (14,file='radfun',form='unformatted',status='unknown')
       WRITE(14) atoms%jri(1:atoms%ntype)
       OPEN (16,file='latharm',form='unformatted',status='unknown')
       WRITE(16) sphhar%nlhd,sphhar%memd
       l = 0
       DO i = 1,atoms%ntype
665
          j = atoms%ntypsy(SUM(atoms%neq(:i-1))+1)
666 667 668 669 670 671 672 673 674 675 676 677 678
          WRITE(16) sphhar%nlh(j),sphhar%llh(:sphhar%nlh(j),j),sphhar%nmem(:sphhar%nlh(j),j),&
               sphhar%mlh(:sphhar%memd,:sphhar%nlh(j),j),sphhar%clnu(:sphhar%memd,:sphhar%nlh(j),j)
          DO j = 1,atoms%neq(i)
             l = l + 1
             IF(atoms%invsat(l).EQ.2) THEN
                WRITE(16) -atoms%ngopr(sym%invsatnr(l))
             ELSE
                WRITE(16)  atoms%ngopr(l)
             ENDIF
          ENDDO
       ENDDO
       CLOSE (16)
       gwc=2
679 680 681 682 683 684 685 686 687

       IF (mpi%irank.EQ.0) THEN
          CALL readPotential(stars,vacuum,atoms,sphhar,input,sym,POT_ARCHIVE_TYPE_COUL_const,&
                             iter,vr,vpw,vz,vzxy)
       END IF
#ifdef CPP_MPI
       CALL mpi_bc_pot(mpi,stars,sphhar,atoms,input,vacuum,&
                       iter,vr,vpw,vz,vzxy)
#endif
688
       GOTO 999
689
    ELSE IF ( input%gw.EQ.2.AND.(gwc==2) )  THEN
690 691 692 693 694 695 696 697 698 699
       CLOSE (12)
       CLOSE (13)
       CLOSE (1013)
       CLOSE (14)
       CLOSE (15)
       IF(.NOT.noco%l_soc)  THEN
          INQUIRE(1014,opened=l_file)
          IF(l_file) CLOSE(1014)
          INQUIRE(667,opened=l_file)
          IF(l_file) CLOSE(667)
700
          CALL juDFT_end("GW finished",mpi%irank)
701 702 703 704
       ENDIF
    ENDIF

    !     hf: write out radial potential vr0
705
    IF (hybrid%l_hybrid.OR.hybrid%l_calhf) THEN
706 707
       OPEN(unit=120,file='vr0',form='unformatted')
       DO isp=1,DIMENSION%jspd
Daniel Wortmann's avatar
Daniel Wortmann committed
708
          DO nn=1,atoms%ntype
709
             DO i=1,atoms%jmtd
Daniel Wortmann's avatar
Daniel Wortmann committed
710
                WRITE(120) vr0(i,nn,isp)
711 712 713 714 715 716 717 718 719
             END DO
          END DO
       END DO
       CLOSE(120)
    ENDIF

    DEALLOCATE ( vpw,vzxy,vz,vr,vr0 )

#ifdef CPP_MPI
720
    CALL MPI_BARRIER(mpi%MPI_COMM,ierr)
721
#endif
722
    IF (hybrid%l_hybrid.OR.hybrid%l_calhf) CALL close_eig(eig_id_hf)
723
    atoms%n_u=n_u_in
724 725


726
    IF( input%jspins .EQ. 1 .AND. hybrid%l_hybrid ) THEN
727 728 729
       results%te_hfex%valence = 2*results%te_hfex%valence
       results%te_hfex%core    = 2*results%te_hfex%core
    END IF
730 731
    enpara_in%epara_min = MINVAL(enpara%el0)
    enpara_in%epara_min = MIN(MINVAL(enpara%ello0),enpara_in%epara_min)
732
!    enpara_in=enpara
733 734
  END SUBROUTINE eigen
END MODULE m_eigen