cdnval.F90 26.8 KB
Newer Older
1 2 3 4
MODULE m_cdnval
  use m_juDFT
CONTAINS
  SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,sliceplot,noco, input,banddos,cell,atoms,enpara,stars,&
5
                    vacuum,dimension,sphhar,sym,obsolete,vTot,oneD,coreSpecInput,den,regCharges,results,&
6
                    moments)
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
    !
    !     ***********************************************************
    !         this subroutin is a modified version of cdnval.F.
    !         it calculates a layer charge distribution and an orbital
    !         composition for each state in a film.
    !         this information is written in the  'ek_orbcomp'  file
    !                                    Yu.Koroteev  01.12.2003
    !     ***********************************************************
    !     flapw7 valence density generator
    !                                         c.l.fu
    !     zelec used to calculate ne - 6.12.95 r.pentcheva
    !
    !     changed subroutine to allow parallel writing of vacdos&dosinp
    !     used temporary direct access file 84,tmp_dos to store data used
    !     in cdninf
    !     call of cdninf changed, sympsi is called from cdnval now
    !     look for 'ifdef CPP_MPI' -blocks!               d.wortmann 6.5.99
    !
    !******** ABBREVIATIONS ************************************************
    !     nbands   : number of bands in the energy window
    !     noccbd   : number of occupied bands
    !     slice    : when set to .true. the charge density of a enery range
    !                (slice) or a single state is calculated
    !     e1s,e2s  : (used if slice) energy range for the slice. if both
    !                are set to 0.0 the charge density of the band nr. nnne
    !                is calculated
    !     pallst   : (used if slice) if set to .true. bands above the
    !                Fermi-Energy are taken into account
    !     nnne     : (used if slice) number of the band of which the charge
    !                density is calculated if e1s and e2s are 0.0
    !     kk       : (used if slice) if set to 0 all k-points are used to
    !                calculate the charge density of the slice, otherwise
    !                only k-points kk is taken into account
    !     nslibd   : number of bands in slice
    !     ener     : band energy averaged over all bands and k-points,
    !                wheighted with the l-like charge of each atom type
    !     sqal     : l-like charge of each atom type. sum over all k-points
    !                and bands
    !***********************************************************************
    !
Gregor Michalicek's avatar
Gregor Michalicek committed
47
    USE m_constants
48
    USE m_eig66_io,ONLY: write_dos
49
    USE m_genMTBasis
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
    USE m_rhomt
    USE m_rhonmt
    USE m_rhomtlo
    USE m_rhonmtlo
    USE m_mcdinit
    USE m_sympsi
    USE m_eparas      ! energy parameters and partial charges
    USE m_qal21       ! off-diagonal part of partial charges
    USE m_abcof
    USE m_nmat        ! calculate density matrix for LDA + U
    USE m_vacden
    USE m_nstm3
    USE m_pwden
    USE m_forcea8
    USE m_forcea12
    USE m_forcea21
66
    USE m_checkdopall
67 68 69 70
    USE m_int21       ! integrate (spin) off-diagonal radial functions
    USE m_int21lo     ! -"- for u_lo
    USE m_rhomt21     ! calculate (spin) off-diagonal MT-density coeff's
    USE m_rhonmt21    ! -"-                       non-MT-density coeff's
71
    USE m_cdnmt       ! calculate the density and orbital moments etc.
72 73 74 75 76 77 78 79
    USE m_orbmom      ! coeffd for orbital moments
    USE m_qmtsl       ! These subroutines divide the input%film into vacuum%layers
    USE m_qintsl      ! (slabs) and intergate the DOS in these vacuum%layers
    USE m_orbcomp     ! calculate corbital composition (like p_x,p_y,p_z)
    USE m_Ekwritesl   ! and write to file.
    USE m_abcrot2
    USE m_doswrite
    USE m_cdnread, ONLY : cdn_read0, cdn_read
80 81 82
    USE m_corespec, only : l_cs    ! calculation of core spectra (EELS)
    USE m_corespec_io, only : corespec_init
    USE m_corespec_eval, only : corespec_gaunt,corespec_rme,corespec_dos,corespec_ddscs
83
#ifdef CPP_MPI
84
    USE m_mpi_col_den ! collect density data from parallel nodes
85 86
#endif
    USE m_types
87
    USE m_xmlOutput
88
    IMPLICIT NONE
89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110
    TYPE(t_results),       INTENT(INOUT) :: results
    TYPE(t_mpi),           INTENT(IN)    :: mpi
    TYPE(t_dimension),     INTENT(IN)    :: dimension
    TYPE(t_oneD),          INTENT(IN)    :: oneD
    TYPE(t_enpara),        INTENT(INOUT) :: enpara
    TYPE(t_obsolete),      INTENT(IN)    :: obsolete
    TYPE(t_banddos),       INTENT(IN)    :: banddos
    TYPE(t_sliceplot),     INTENT(IN)    :: sliceplot
    TYPE(t_input),         INTENT(IN)    :: input
    TYPE(t_vacuum),        INTENT(IN)    :: vacuum
    TYPE(t_noco),          INTENT(IN)    :: noco
    TYPE(t_sym),           INTENT(IN)    :: sym
    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(IN)    :: atoms
    TYPE(t_coreSpecInput), INTENT(IN)    :: coreSpecInput
    TYPE(t_potden),        INTENT(IN)    :: vTot
    TYPE(t_potden),        INTENT(INOUT) :: den
    TYPE(t_regionCharges), INTENT(INOUT) :: regCharges
    TYPE(t_moments),       INTENT(INOUT) :: moments
111

112
    !     .. Scalar Arguments ..
113
    INTEGER, INTENT(IN)              :: eig_id,jspin
114 115 116

#ifdef CPP_MPI
    INCLUDE 'mpif.h'
117
    LOGICAL :: mpi_flag, mpi_status
118 119
#endif
    !     .. Local Scalars ..
120
    INTEGER :: llpd,ikpt,jsp_start,jsp_end,ispin
121
    INTEGER :: i,ie,iv,ivac,j,k,l,n,ilo,isp,nbands,noccbd,nslibd
122
    INTEGER :: skip_t,skip_tt, nkpt_extended, ikptStart, ikptIncrement
123
    INTEGER :: n_start,n_end,noccbd_l,nbasfcn
124
    LOGICAL :: l_fmpl,l_evp,l_orbcomprot,l_real, l_write
125
    !     ...Local Arrays ..
126 127 128
    INTEGER :: noccbd_in(kpts%nkpt)
    INTEGER :: nStart_in(kpts%nkpt)
    INTEGER :: nEnd_in(kpts%nkpt)
129
    REAL    :: eig(dimension%neigd)
130

Daniel Wortmann's avatar
Daniel Wortmann committed
131
    INTEGER, ALLOCATABLE :: gvac1d(:),gvac2d(:)
132
    INTEGER, ALLOCATABLE :: jsym(:),ksym(:)
133
    REAL,    ALLOCATABLE :: we(:)
134
    REAL,    ALLOCATABLE :: f(:,:,:,:),g(:,:,:,:),flo(:,:,:,:) ! radial functions
135

136
    TYPE (t_lapw)             :: lapw
137 138 139
    TYPE (t_orb)              :: orb
    TYPE (t_denCoeffs)        :: denCoeffs
    TYPE (t_denCoeffsOffdiag) :: denCoeffsOffdiag
140
    TYPE (t_force)            :: force
141
    TYPE (t_slab)             :: slab
142
    TYPE (t_eigVecCoeffs)     :: eigVecCoeffs
143
    TYPE (t_mcd)              :: mcd
144 145
    TYPE (t_usdus)            :: usdus
    TYPE (t_zMat)             :: zMat
146
    TYPE (t_orbcomp)          :: orbcomp
147

148
    l_real = sym%invs.AND.(.NOT.noco%l_soc).AND.(.NOT.noco%l_noco)
Gregor Michalicek's avatar
Gregor Michalicek committed
149

150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
    llpd=(atoms%lmaxd*(atoms%lmaxd+3))/2
    !---> l_fmpl is meant as a switch to to a plot of the full magnet.
    !---> density without the atomic sphere approximation for the magnet.
    !---> density. It is not completely implemented (lo's missing).
    l_fmpl = .false.
    IF (noco%l_mperp) THEN
       !--->    when the off-diag. part of the desinsity matrix, i.e. m_x and
       !--->    m_y, is calculated inside the muffin-tins (l_mperp = T), cdnval
       !--->    is called only once. therefore, several spin loops have been
       !--->    added. if l_mperp = F, these loops run only from jspin - jspin.
       jsp_start = 1
       jsp_end   = 2
    ELSE
       jsp_start = jspin
       jsp_end   = jspin
    ENDIF
166
    !---> if l_mperp = F, these variables are only needed for one spin
167 168 169
    !---> at a time, otherwise for both spins:
    ALLOCATE ( f(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )      ! Deallocation before mpi_col_den
    ALLOCATE ( g(atoms%jmtd,2,0:atoms%lmaxd,jsp_start:jsp_end) )
170
    ALLOCATE ( flo(atoms%jmtd,2,atoms%nlod,dimension%jspd) )
171 172
    ALLOCATE ( jsym(dimension%neigd),ksym(dimension%neigd) )
    ALLOCATE ( gvac1d(dimension%nv2d),gvac2d(dimension%nv2d) )
173

174
    ! --> Initializations
Gregor Michalicek's avatar
Gregor Michalicek committed
175

176
    CALL usdus%init(atoms,input%jspins)
177 178
    CALL denCoeffs%init(atoms,sphhar,jsp_start,jsp_end)
    CALL denCoeffsOffdiag%init(atoms,noco,sphhar,l_fmpl)
179
    CALL force%init1(input,atoms)
180
    CALL orb%init(atoms,noco,jsp_start,jsp_end)
181
    CALL mcd%init1(banddos,dimension,input,atoms)
182 183
    CALL slab%init(banddos,dimension,atoms,cell)
    CALL orbcomp%init(banddos,dimension,atoms)
184

185
    IF ((l_fmpl).AND.(.not.noco%l_mperp)) CALL juDFT_error("for fmpl set noco%l_mperp = T!" ,calledby ="cdnval")
186
    IF ((banddos%ndir.EQ.-3).AND.banddos%dos.AND.oneD%odi%d1) CALL juDFT_error("layer-resolved feature does not work with 1D",calledby ="cdnval")
187

188
! calculation of core spectra (EELS) initializations -start-
189
    CALL corespec_init(input,atoms,coreSpecInput)
Gregor Michalicek's avatar
Gregor Michalicek committed
190
    IF(l_cs.AND.(mpi%isize.NE.1)) CALL juDFT_error('EELS + MPI not implemented', calledby = 'cdnval')
191 192
    IF(l_cs.AND.jspin.EQ.1) CALL corespec_gaunt()
! calculation of core spectra (EELS) initializations -end-
193

194 195 196
    IF (mpi%irank==0) THEN
       WRITE (6,FMT=8000) jspin
       WRITE (16,FMT=8000) jspin
197
       CALL openXMLElementPoly('mtCharges',(/'spin'/),(/jspin/))
198
    END IF
199 200 201 202
8000 FORMAT (/,/,10x,'valence density: spin=',i2)

    skip_tt = dot_product(enpara%skiplo(:atoms%ntype,jspin),atoms%neq(:atoms%ntype))
    IF (noco%l_soc.OR.noco%l_noco)  skip_tt = 2 * skip_tt
203

204 205
    l_write = input%cdinf.AND.mpi%irank==0

206
    DO n = 1,atoms%ntype
207 208 209

       DO ispin = jsp_start, jsp_end
          CALL genMTBasis(atoms,enpara,vTot,mpi,n,ispin,l_write,usdus,f(:,:,0:,ispin),g(:,:,0:,ispin),flo(:,:,:,ispin))
210
       END DO
211 212 213 214 215 216 217 218 219 220

       IF (noco%l_mperp) THEN
          DO l = 0,atoms%lmax(n)
             CALL int_21(f,g,atoms,n,l,denCoeffsOffdiag)
          END DO
          DO ilo = 1, atoms%nlo(n)
             CALL int_21lo(f,g,atoms,n,flo,ilo,denCoeffsOffdiag)
          END DO
       END IF

221
       IF (banddos%l_mcd) CALL mcd_init(atoms,input,dimension,vTot%mt(:,0,:,:),g,f,mcd,n,jspin)
222 223 224

       IF(l_cs) CALL corespec_rme(atoms,input,n,dimension%nstd,&
                                  input%jspins,jspin,results%ef,&
225
                                  dimension%msh,vTot%mt(:,0,:,:),f,g)
226
    END DO
227
    DEALLOCATE (f,g,flo)
228

229 230 231 232
    ALLOCATE (we(dimension%neigd))

    ! For k-point paralelization:
    l_evp = .FALSE.
233
    IF (kpts%nkpt < mpi%isize) THEN
234 235 236 237
       l_evp = .TRUE.
       nkpt_extended = kpts%nkpt
       ikptStart = 1
       ikptIncrement = 1
238
    ELSE
239 240 241 242
       ! the number of iterations is adjusted to the number of MPI processes to synchronize RMA operations
       nkpt_extended = (kpts%nkpt / mpi%isize + 1) * mpi%isize
       ikptStart = mpi%irank + 1
       ikptIncrement = mpi%isize
243
    END IF
244

245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
    ! determine bands to be used for each k point, MPI process
    DO ikpt = ikptStart, kpts%nkpt, ikptIncrement
       noccbd_in(ikpt) = 0
       DO i = 1,MERGE(results%neig(ikpt,1),results%neig(ikpt,jspin),noco%l_noco)
          we(i) = MERGE(results%w_iks(i,ikpt,1),results%w_iks(i,ikpt,jspin),noco%l_noco)
          IF ((we(i).GE.1.e-8).OR.input%pallst) THEN
             noccbd_in(ikpt) = noccbd_in(ikpt) + 1
          END IF
       END DO
       IF (banddos%dos) noccbd_in(ikpt) = MERGE(results%neig(ikpt,1),results%neig(ikpt,jspin),noco%l_noco)
       IF (l_evp) THEN
          noccbd_l = CEILING(real(noccbd_in(ikpt)) / mpi%isize)
          nStart_in(ikpt) = mpi%irank*noccbd_l + 1
          nEnd_in(ikpt)   = min((mpi%irank+1)*noccbd_l, noccbd_in(ikpt))

          noccbd_in(ikpt) = nEnd_in(ikpt) - nStart_in(ikpt) + 1
          IF (noccbd_in(ikpt).LT.1) noccbd_in(ikpt) = 0
       ELSE
          nStart_in(ikpt) = 1
          nEnd_in(ikpt)   = noccbd_in(ikpt)
       END IF
    END DO

268
    DO ikpt = ikptStart, nkpt_extended, ikptIncrement
269 270 271 272 273 274 275 276 277 278
       IF (ikpt.GT.kpts%nkpt) THEN
#ifdef CPP_MPI
          ! Synchronizes the RMA operations
          CALL MPI_BARRIER(mpi%mpi_comm,ie)
#endif
          EXIT
       END IF

       we=0.0
       !--->    determine number of occupied bands and set weights (we)
279
       DO i = nStart_in(ikpt),nEnd_in(ikpt)
280
          we(i) = MERGE(results%w_iks(i,ikpt,1),results%w_iks(i,ikpt,jspin),noco%l_noco)
281
          IF (.NOT.((we(i).GE.1.e-8).OR.input%pallst)) we(i)=0.0
282 283 284 285 286 287 288
       END DO
       ! uncomment this so that cdinf plots works for all states
       ! noccbd = neigd

       ! -> Gu test: distribute ev's among the processors...
       CALL lapw%init(input,noco, kpts,atoms,sym,ikpt,cell,.false., mpi)
       skip_t = skip_tt
289 290 291 292 293 294
       noccbd = noccbd_in(ikpt)
       n_start = nStart_in(ikpt)
       n_end = nEnd_in(ikpt)
       IF (l_evp) THEN
          IF(noccbd.GT.0) THEN
             we(1:noccbd) = we(nStart_in(ikpt):nEnd_in(ikpt))
295
          END IF
296

297
          IF (n_start > skip_tt) skip_t = 0
298 299 300
          IF (n_end <= skip_tt) skip_t = noccbd
          IF ((n_start <= skip_tt).AND.(n_end > skip_tt)) skip_t = mod(skip_tt,noccbd)
       END IF
301

302
       nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
303
       CALL zMat%init(l_real,nbasfcn,noccbd)
304 305 306
       CALL cdn_read(eig_id,dimension%nvd,dimension%jspd,mpi%irank,mpi%isize,&
                     ikpt,jspin,zmat%nbasfcn,noco%l_ss,noco%l_noco,&
                     noccbd,n_start,n_end,nbands,eig,zMat)
307
#ifdef CPP_MPI
308 309
       ! Synchronizes the RMA operations
       CALL MPI_BARRIER(mpi%mpi_comm,ie)
310
#endif
311 312 313 314 315 316

       IF (noco%l_noco) THEN
          eig(1:noccbd) = results%eig(n_start:n_end,ikpt,1)
       ELSE
          eig(1:noccbd) = results%eig(n_start:n_end,ikpt,jspin)
       END IF
317 318 319 320 321 322

       IF (vacuum%nstm.EQ.3.AND.input%film) THEN
          CALL nstm3(sym,atoms,vacuum,stars,ikpt,lapw%nv(jspin),input,jspin,kpts,&
                     cell,kpts%wtkpt(ikpt),lapw%k1(:,jspin),lapw%k2(:,jspin),&
                     enpara%evac0(1,jspin),vTot%vacz(:,:,jspin),gvac1d,gvac2d)
       END IF
323

324 325 326 327 328
       !--->    if slice, only a certain bands are taken into account
       !--->    in order to do this the coresponding eigenvalues, eigenvectors
       !--->    and weights have to be copied to the beginning of the arrays
       !--->    eig, z and we and the number of occupied bands (noccbd) has to
       !--->    changed
329
       IF (sliceplot%slice.AND.noccbd.GT.0) THEN
330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348
          IF (mpi%irank==0) WRITE (16,FMT=*) 'NNNE',sliceplot%nnne
          IF (mpi%irank==0) WRITE (16,FMT=*) 'sliceplot%kk',sliceplot%kk
          nslibd = 0
          IF (input%pallst) we(:nbands) = kpts%wtkpt(ikpt)
          IF (sliceplot%kk.EQ.0) THEN
             IF (mpi%irank==0) THEN
                WRITE (16,FMT='(a)') 'ALL K-POINTS ARE TAKEN IN SLICE'
                WRITE (16,FMT='(a,i2)') ' sliceplot%slice: k-point nr.',ikpt
             END IF
             DO i = 1,nbands
                IF (eig(i).GE.sliceplot%e1s .AND. eig(i).LE.sliceplot%e2s) THEN
                   nslibd = nslibd + 1
                   eig(nslibd) = eig(i)
                   we(nslibd) = we(i)
                   if (zmat%l_real) THEN
                      zMat%z_r(:,nslibd) = zMat%z_r(:,i)
                   else
                      zMat%z_c(:,nslibd) = zMat%z_c(:,i)
                   endif
349
                END IF
350 351 352 353 354 355 356
             END DO
             IF (mpi%irank==0) WRITE (16,'(a,i3)') ' eigenvalues in sliceplot%slice:',nslibd
          ELSE IF (sliceplot%kk.EQ.ikpt) THEN
             IF (mpi%irank==0) WRITE (16,FMT='(a,i2)') ' sliceplot%slice: k-point nr.',ikpt
             IF ((sliceplot%e1s.EQ.0.0) .AND. (sliceplot%e2s.EQ.0.0)) THEN
                IF (mpi%irank==0) WRITE (16,FMT='(a,i5,f10.5)') 'slice: eigenvalue nr.',&
                     sliceplot%nnne,eig(sliceplot%nnne)
357 358 359
                nslibd = 1
                eig(1) = eig(sliceplot%nnne)
                we(1) = we(sliceplot%nnne)
360
                if (zmat%l_real) Then
361
                   zMat%z_r(:,1) = zMat%z_r(:,sliceplot%nnne)
362
                else
363
                   zMat%z_c(:,1) = zMat%z_c(:,sliceplot%nnne)
364 365
                endif
             ELSE
366 367 368 369 370
                DO i = 1,nbands
                   IF (eig(i).GE.sliceplot%e1s .AND. eig(i).LE.sliceplot%e2s) THEN
                      nslibd = nslibd + 1
                      eig(nslibd) = eig(i)
                      we(nslibd) = we(i)
371
                      if (zmat%l_real) THEN
372
                         zMat%z_r(:,nslibd) = zMat%z_r(:,i)
373
                      else
374
                         zMat%z_c(:,nslibd) = zMat%z_c(:,i)
375
                      endif
376 377
                   END IF
                END DO
378
                IF (mpi%irank==0) WRITE (16,FMT='(a,i3)')' eigenvalues in sliceplot%slice:',nslibd
379 380
             END IF
          END IF
381 382 383
          noccbd = nslibd
       END IF ! sliceplot%slice

384 385
       IF (noccbd.EQ.0) GO TO 199

386 387 388 389 390 391 392 393 394
       !--->    in normal iterations the charge density of the unoccupied
       !--->    does not need to be calculated (in pwden, vacden and abcof)
       IF (banddos%dos.AND. .NOT.(l_evp.AND.(mpi%isize.GT.1)) ) THEN
          noccbd=nbands
       END IF
       !     ----> add in spin-doubling factor
       we(:noccbd) = 2.0 * we(:noccbd) / input%jspins

       !--->    valence density in the interstitial and vacuum region
395 396
       !--->    has to be called only once (if jspin=1) in the non-collinear case

397 398 399 400 401 402 403
       !     ----> valence density in the interstitial region
       IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
          CALL timestart("cdnval: pwden")
          CALL pwden(stars,kpts,banddos,oneD,input,mpi,noco,cell,atoms,sym,ikpt,&
                     jspin,lapw,noccbd,we,eig,den,regCharges%qis,results,force%f_b8,zMat)
          CALL timestop("cdnval: pwden")
       END IF
404

405 406 407
       !--->    charge of each valence state in this k-point of the SBZ
       !--->    in the layer interstitial region of the film
       IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
408
          IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
409
             CALL q_int_sl(jspin,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat)
410
          END IF
411
       END IF
412

413 414 415 416 417 418 419 420
       !--->    valence density in the vacuum region
       IF (input%film) THEN
          IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
             CALL timestart("cdnval: vacden")
             CALL vacden(vacuum,dimension,stars,oneD, kpts,input, cell,atoms,noco,banddos,&
                         gvac1d,gvac2d, we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,lapw,enpara%evac0,eig,&
                         den,regCharges%qvac,regCharges%qvlay,regCharges%qstars,zMat)
             CALL timestop("cdnval: vacden")
421
          END IF
422 423 424 425 426 427
          !--->       perform Brillouin zone integration and summation over the
          !--->       bands in order to determine the vacuum energy parameters.
          DO ispin = jsp_start,jsp_end
             DO ivac = 1,vacuum%nvac
                regCharges%pvac(ivac,ispin)=regCharges%pvac(ivac,ispin)+dot_product(eig(:noccbd)*regCharges%qvac(:noccbd,ivac,ikpt,ispin),we(:noccbd))
                regCharges%svac(ivac,ispin)=regCharges%svac(ivac,ispin)+dot_product(regCharges%qvac(:noccbd,ivac,ikpt,ispin),we(:noccbd))
428
             END DO
429 430
          END DO
       END IF
431

432 433
       !--->    valence density in the atomic spheres
       CALL eigVecCoeffs%init(dimension,atoms,noco,jspin,noccbd)
434

435 436 437 438 439 440 441 442 443 444 445 446 447
       DO ispin = jsp_start,jsp_end
          IF (input%l_f) THEN
             CALL force%init2(noccbd,input,atoms)
          END IF
          CALL timestart("cdnval: abcof")
          CALL abcof(input,atoms,sym,cell,lapw,noccbd,usdus,noco,ispin,oneD,&
                     eigVecCoeffs%acof(:,0:,:,ispin),eigVecCoeffs%bcof(:,0:,:,ispin),&
                     eigVecCoeffs%ccof(-atoms%llod:,:,:,:,ispin),zMat,eig,force)
          CALL timestop("cdnval: abcof")

          IF (atoms%n_u.GT.0) THEN
             CALL n_mat(atoms,sym,noccbd,usdus,ispin,we,eigVecCoeffs,den%mmpMat(:,:,:,jspin))
          END IF
448

449 450 451 452 453 454
          !--->       perform Brillouin zone integration and summation over the
          !--->       bands in order to determine the energy parameters for each
          !--->       atom and angular momentum
          IF (.not.sliceplot%slice) THEN
             CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,&
                         skip_t,l_evp,eigVecCoeffs,usdus,regCharges,mcd,banddos%l_mcd)
455

456 457
             IF (noco%l_mperp.AND.(ispin == jsp_end)) THEN
                CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,regCharges)
458
             END IF
459
          END IF
460

461 462 463 464
          !--->    layer charge of each valence state in this k-point of the SBZ
          !--->    from the mt-sphere region of the film
          IF (banddos%dos.AND.(banddos%ndir.EQ.-3))  THEN
             CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
465

466 467 468
             INQUIRE (file='orbcomprot',exist=l_orbcomprot)
             IF (l_orbcomprot) THEN                           ! rotate ab-coeffs
                CALL abcrot2(atoms,noccbd,eigVecCoeffs,ispin)
469
             END IF
470

471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502
             CALL orb_comp(ispin,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp)
          END IF
          !-new
          !--->          set up coefficients for the spherical and
          CALL timestart("cdnval: rhomt")
          CALL rhomt(atoms,we,noccbd,eigVecCoeffs,denCoeffs,ispin)
          CALL timestop("cdnval: rhomt")

          IF (noco%l_soc) CALL orbmom(atoms,noccbd,we,ispin,eigVecCoeffs,orb)
          !--->          non-spherical m.t. density
          CALL timestart("cdnval: rhonmt")
          CALL rhonmt(atoms,sphhar,we,noccbd,sym,eigVecCoeffs,denCoeffs,ispin)
          CALL timestop("cdnval: rhonmt")

          !--->          set up coefficients of the local orbitals and the
          !--->          flapw - lo cross terms for the spherical and
          !--->          non-spherical mt density
          CALL timestart("cdnval: rho(n)mtlo")
          CALL rhomtlo(atoms,noccbd,we,eigVecCoeffs,denCoeffs,ispin)

          CALL rhonmtlo(atoms,sphhar,noccbd,we,eigVecCoeffs,denCoeffs,ispin)
          CALL timestop("cdnval: rho(n)mtlo")

          IF (input%l_f) THEN
             CALL timestart("cdnval: force_a12/21")
             IF (.not.input%l_useapw) THEN
                CALL force_a12(atoms,noccbd,sym,dimension,cell,oneD,&
                               we,ispin,noccbd,usdus,eigVecCoeffs,force,results)
             ENDIF
             CALL force_a21(input,atoms,dimension,noccbd,sym,oneD,cell,we,ispin,&
                            enpara%el0(0:,:,ispin),noccbd,eig,usdus,eigVecCoeffs,force,results)
             CALL timestop("cdnval: force_a12/21")
503 504
          END IF

505 506 507 508 509
          IF(l_cs) THEN
             CALL corespec_dos(atoms,usdus,ispin,dimension%lmd,kpts%nkpt,ikpt,dimension%neigd,&
                               noccbd,results%ef,banddos%sig_dos,eig,we,eigVecCoeffs)
          END IF
       END DO !--->    end loop over ispin
510

511 512 513 514
       IF (noco%l_mperp) THEN
          CALL rhomt21(atoms,we,noccbd,eigVecCoeffs,denCoeffsOffdiag)
          IF (l_fmpl) CALL rhonmt21(atoms,sphhar,we,noccbd,sym,eigVecCoeffs,denCoeffsOffdiag)
       END IF
515

516 517 518 519 520 521
199    CONTINUE
       IF ((banddos%dos .OR. banddos%vacdos .OR. input%cdinf)  ) THEN
          CALL timestart("cdnval: write_info")
          !--->    calculate charge distribution of each state (l-character ...)
          !--->    and write the information to the files dosinp and vacdos
          !--->    for dos and bandstructure plots
522

523 524 525
          !--dw    parallel writing of vacdos,dosinp....
          !        write data to direct access file first, write to formated file later by PE 0 only!
          !--dw    since z is no longer an argument of cdninf sympsi has to be called here!
526

527 528 529
          IF (banddos%ndir.GT.0) THEN
             CALL sympsi(lapw%bkpt,lapw%nv(jspin),lapw%k1(:,jspin),lapw%k2(:,jspin),&
                         lapw%k3(:,jspin),sym,dimension,nbands,cell,eig,noco, ksym,jsym,zMat)
530 531
          END IF

532 533 534 535 536 537 538
          !--dw   now write k-point data to tmp_dos
          CALL write_dos(eig_id,ikpt,jspin,regCharges%qal(:,:,:,jspin),regCharges%qvac(:,:,ikpt,jspin),regCharges%qis(:,ikpt,jspin),&
                         regCharges%qvlay(:,:,:,ikpt,jspin),regCharges%qstars,ksym,jsym,mcd%mcd,slab%qintsl,&
                         slab%qmtsl(:,:),orbcomp%qmtp(:,:),orbcomp%comp)

          CALL timestop("cdnval: write_info")
          !-new_sl
539
       END IF
540
    END DO !---> end of k-point loop
541
    DEALLOCATE (we)
542
    !+t3e
543
#ifdef CPP_MPI
544 545 546
    CALL timestart("cdnval: mpi_col_den")
    DO ispin = jsp_start,jsp_end
       CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,&
547
                        input,noco,l_fmpl,ispin,llpd, den%vacxy(1,1,1,ispin),&
548 549
                        den%vacz(1,1,ispin),den%pw(1,ispin),regCharges,&
                        results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin))
550 551
    END DO
    CALL timestop("cdnval: mpi_col_den")
552
#endif
553

554 555 556 557 558 559
    IF (mpi%irank==0) THEN
       CALL cdnmt(dimension%jspd,atoms,sphhar,llpd,noco,l_fmpl,jsp_start,jsp_end,&
                  enpara%el0,enpara%ello0,vTot%mt(:,0,:,:),denCoeffs,&
                  usdus,orb,denCoeffsOffdiag,moments,den%mt)

       IF(l_cs) CALL corespec_ddscs(jspin,input%jspins)
560

561 562
       IF (((jspin.eq.input%jspins).OR.noco%l_mperp) .AND. (banddos%dos.or.banddos%vacdos.or.input%cdinf) ) THEN
          CALL timestart("cdnval: dos")
563
          CALL doswrite(eig_id,dimension,kpts,atoms,vacuum,input,banddos,&
564
                        sliceplot,noco,sym,cell,mcd,results,slab%nsld,oneD)
565
          IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
566
             CALL Ek_write_sl(eig_id,dimension,kpts,atoms,vacuum,input,jspin,sym,cell,slab)
567
          END IF
568
          CALL timestop("cdnval: dos")
569 570 571
       END IF

       DO ispin = jsp_start,jsp_end
572 573
          IF (.NOT.sliceplot%slice) THEN
             DO n=1,atoms%ntype
574 575
                enpara%el1(0:3,n,ispin)=regCharges%ener(0:3,n,ispin)/regCharges%sqal(0:3,n,ispin)
                IF (atoms%nlo(n)>0) enpara%ello1(:atoms%nlo(n),n,ispin)=regCharges%enerlo(:atoms%nlo(n),n,ispin)/regCharges%sqlo(:atoms%nlo(n),n,ispin)
576
             END DO
577
             IF (input%film) enpara%evac1(:vacuum%nvac,ispin)=regCharges%pvac(:vacuum%nvac,ispin)/regCharges%svac(:vacuum%nvac,ispin)
578 579 580 581 582 583 584
          END IF

          !--->      check continuity of charge density
          IF (input%cdinf) THEN
             CALL timestart("cdnval: cdninf-stuff")
             WRITE (6,FMT=8210) ispin
8210         FORMAT (/,5x,'check continuity of cdn for spin=',i2)
585
             CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,den,ispin)
586 587 588 589 590 591
             CALL timestop("cdnval: cdninf-stuff")
          END IF
          !+for
          !--->      forces of equ. A8 of Yu et al.
          IF ((input%l_f)) THEN
             CALL timestart("cdnval: force_a8")
592
             CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),den%mt,force,results)
593 594
             CALL timestop("cdnval: force_a8")
          END IF
595
          !-for
596
       END DO ! end of loop ispin = jsp_start,jsp_end
597
       CALL closeXMLElement('mtCharges')
598 599 600 601 602 603 604 605 606 607

       IF(vacuum%nvac.EQ.1) THEN
          den%vacz(:,2,:) = den%vacz(:,1,:)
          IF (sym%invs) THEN
             den%vacxy(:,:,2,:) = CONJG(den%vacxy(:,:,1,:))
          ELSE
             den%vacxy(:,:,2,:) = den%vacxy(:,:,1,:)
          END IF
       END IF

608
    END IF ! end of (mpi%irank==0)
609 610 611
#ifdef CPP_MPI
    CALL MPI_BARRIER(mpi%mpi_comm,ie)
#endif
612

613 614 615 616
    IF ((jsp_end.EQ.input%jspins)) THEN
       IF ((banddos%dos.OR.banddos%vacdos).AND.(banddos%ndir/=-2))  CALL juDFT_end("DOS OK",mpi%irank)
       IF (vacuum%nstm.EQ.3)  CALL juDFT_end("VACWAVE OK",mpi%irank)
    END IF
617

618 619
  END SUBROUTINE cdnval
END MODULE m_cdnval