cdnval.F90 26 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
    REAL    :: eig(dimension%neigd)
127

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

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

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

147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
    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
163
    !---> if l_mperp = F, these variables are only needed for one spin
164 165 166
    !---> 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) )
167
    ALLOCATE ( flo(atoms%jmtd,2,atoms%nlod,dimension%jspd) )
168 169
    ALLOCATE ( jsym(dimension%neigd),ksym(dimension%neigd) )
    ALLOCATE ( gvac1d(dimension%nv2d),gvac2d(dimension%nv2d) )
170

171
    ! --> Initializations
Gregor Michalicek's avatar
Gregor Michalicek committed
172

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

182
    IF ((l_fmpl).AND.(.not.noco%l_mperp)) CALL juDFT_error("for fmpl set noco%l_mperp = T!" ,calledby ="cdnval")
183
    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")
184

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

191 192 193
    IF (mpi%irank==0) THEN
       WRITE (6,FMT=8000) jspin
       WRITE (16,FMT=8000) jspin
194
       CALL openXMLElementPoly('mtCharges',(/'spin'/),(/jspin/))
195
    END IF
196 197 198 199
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
200

201 202
    l_write = input%cdinf.AND.mpi%irank==0

203
    DO n = 1,atoms%ntype
204 205 206

       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))
207
       END DO
208 209 210 211 212 213 214 215 216 217

       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

218
       IF (banddos%l_mcd) CALL mcd_init(atoms,input,dimension,vTot%mt(:,0,:,:),g,f,mcd,n,jspin)
219 220 221

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

226 227 228 229
    ALLOCATE (we(dimension%neigd))

    ! For k-point paralelization:
    l_evp = .FALSE.
230
    IF (kpts%nkpt < mpi%isize) THEN
231 232 233 234
       l_evp = .TRUE.
       nkpt_extended = kpts%nkpt
       ikptStart = 1
       ikptIncrement = 1
235
    ELSE
236 237 238 239
       ! 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
240
    END IF
241

242
    DO ikpt = ikptStart, nkpt_extended, ikptIncrement
243 244 245 246 247 248 249 250 251 252 253 254
       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)
       noccbd = 0
       DO i = 1,dimension%neigd ! nbands
255
          we(i) = MERGE(results%w_iks(i,ikpt,1),results%w_iks(i,ikpt,jspin),noco%l_noco)
256 257 258 259 260 261 262 263 264 265 266 267
          IF ((we(i).GE.1.e-8).OR.input%pallst) THEN
             noccbd = noccbd + 1
          ELSE
             we(i)=0.0
          END IF
       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
268
       IF (banddos%dos) noccbd = dimension%neigd
269 270 271 272 273 274 275 276 277 278 279
       IF (l_evp.AND.(mpi%isize.GT.1)) THEN
          noccbd_l = CEILING(real(noccbd) / mpi%isize)
          n_start = mpi%irank*noccbd_l + 1
          n_end   = min( (mpi%irank+1)*noccbd_l , noccbd )

          noccbd = n_end - n_start + 1
          IF (noccbd<1) THEN
             noccbd=0
          ELSE
             we(1:noccbd) = we(n_start:n_end)
          END IF
280

281
          IF (n_start > skip_tt) skip_t = 0
282 283 284 285
          IF (n_end <= skip_tt) skip_t = noccbd
          IF ((n_start <= skip_tt).AND.(n_end > skip_tt)) skip_t = mod(skip_tt,noccbd)
       ELSE
          n_start = 1
286 287
          noccbd_l = noccbd
          n_end    = noccbd
288
       END IF
289

290 291 292 293 294
       nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
       CALL zMat%init(l_real,nbasfcn,dimension%neigd)
       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)
295
#ifdef CPP_MPI
296 297
       ! Synchronizes the RMA operations
       CALL MPI_BARRIER(mpi%mpi_comm,ie)
298
#endif
299 300 301 302 303 304 305 306 307
       !IF (l_evp.AND.(isize.GT.1)) THEN
       !  eig(1:noccbd) = eig(n_start:n_end)
       !ENDIF

       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
308

309 310 311 312 313
       !--->    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
314
       IF (sliceplot%slice.AND.noccbd.GT.0) THEN
315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333
          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
334
                END IF
335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350
             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)
                nslibd = nslibd + 1
                eig(nslibd) = eig(sliceplot%nnne)
                we(nslibd) = we(sliceplot%nnne)
                if (zmat%l_real) Then
                   zMat%z_r(:,nslibd) = zMat%z_r(:,sliceplot%nnne)
                else
                   zMat%z_c(:,nslibd) = zMat%z_c(:,sliceplot%nnne)
                endif
             ELSE
351 352 353 354 355
                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)
356
                      if (zmat%l_real) THEN
357
                         zMat%z_r(:,nslibd) = zMat%z_r(:,i)
358
                      else
359
                         zMat%z_c(:,nslibd) = zMat%z_c(:,i)
360
                      endif
361 362
                   END IF
                END DO
363
                IF (mpi%irank==0) WRITE (16,FMT='(a,i3)')' eigenvalues in sliceplot%slice:',nslibd
364 365
             END IF
          END IF
366 367 368
          noccbd = nslibd
       END IF ! sliceplot%slice

369 370
       IF (noccbd.EQ.0) GO TO 199

371 372 373 374 375 376 377 378 379
       !--->    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
380 381
       !--->    has to be called only once (if jspin=1) in the non-collinear case

382 383 384 385 386 387 388
       !     ----> 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
389

390 391 392
       !--->    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
393
          IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
394
             CALL q_int_sl(jspin,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat)
395
          END IF
396
       END IF
397

398 399 400 401 402 403 404 405
       !--->    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")
406
          END IF
407 408 409 410 411 412
          !--->       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))
413
             END DO
414 415
          END DO
       END IF
416

417 418
       !--->    valence density in the atomic spheres
       CALL eigVecCoeffs%init(dimension,atoms,noco,jspin,noccbd)
419

420 421 422 423 424 425 426 427 428 429 430 431 432
       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
433

434 435 436 437 438 439
          !--->       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)
440

441 442
             IF (noco%l_mperp.AND.(ispin == jsp_end)) THEN
                CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,regCharges)
443
             END IF
444
          END IF
445

446 447 448 449
          !--->    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)
450

451 452 453
             INQUIRE (file='orbcomprot',exist=l_orbcomprot)
             IF (l_orbcomprot) THEN                           ! rotate ab-coeffs
                CALL abcrot2(atoms,noccbd,eigVecCoeffs,ispin)
454
             END IF
455

456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487
             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")
488 489
          END IF

490 491 492 493 494
          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
495

496 497 498 499
       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
500

501 502 503 504 505 506
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
507

508 509 510
          !--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!
511

512 513 514
          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)
515 516
          END IF

517 518 519 520 521 522 523
          !--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
524
       END IF
525
    END DO !---> end of k-point loop
526
    DEALLOCATE (we)
527
    !+t3e
528
#ifdef CPP_MPI
529 530 531
    CALL timestart("cdnval: mpi_col_den")
    DO ispin = jsp_start,jsp_end
       CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,&
532
                        input,noco,l_fmpl,ispin,llpd, den%vacxy(1,1,1,ispin),&
533 534
                        den%vacz(1,1,ispin),den%pw(1,ispin),regCharges,&
                        results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin))
535 536
    END DO
    CALL timestop("cdnval: mpi_col_den")
537
#endif
538

539 540 541 542 543 544
    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)
545

546 547
       IF (((jspin.eq.input%jspins).OR.noco%l_mperp) .AND. (banddos%dos.or.banddos%vacdos.or.input%cdinf) ) THEN
          CALL timestart("cdnval: dos")
548
          CALL doswrite(eig_id,dimension,kpts,atoms,vacuum,input,banddos,&
549
                        sliceplot,noco,sym,cell,mcd,results,slab%nsld,oneD)
550
          IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
551
             CALL Ek_write_sl(eig_id,dimension,kpts,atoms,vacuum,input,jspin,sym,cell,slab)
552
          END IF
553
          CALL timestop("cdnval: dos")
554 555 556
       END IF

       DO ispin = jsp_start,jsp_end
557 558
          IF (.NOT.sliceplot%slice) THEN
             DO n=1,atoms%ntype
559 560
                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)
561
             END DO
562
             IF (input%film) enpara%evac1(:vacuum%nvac,ispin)=regCharges%pvac(:vacuum%nvac,ispin)/regCharges%svac(:vacuum%nvac,ispin)
563 564 565 566 567 568 569
          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)
570
             CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,den,ispin)
571 572 573 574 575 576
             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")
577
             CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),den%mt,force,results)
578 579
             CALL timestop("cdnval: force_a8")
          END IF
580
          !-for
581
       END DO ! end of loop ispin = jsp_start,jsp_end
582
       CALL closeXMLElement('mtCharges')
583 584 585 586 587 588 589 590 591 592

       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

593
    END IF ! end of (mpi%irank==0)
594 595 596
#ifdef CPP_MPI
    CALL MPI_BARRIER(mpi%mpi_comm,ie)
#endif
597

598 599 600 601
    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
602

603 604
  END SUBROUTINE cdnval
END MODULE m_cdnval