cdnval.F90 27.5 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_size,i_rec,n_rank,n_start,n_end,noccbd_l,nbasfcn
124
    LOGICAL :: l_fmpl,l_evp,l_orbcomprot,l_real, l_write
125
    !     ...Local Arrays ..
126 127
    INTEGER :: n_bands(0:dimension%neigd)
    REAL    :: eig(dimension%neigd)
128

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

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

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

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

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

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

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

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

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

199
    CALL cdn_read0(eig_id,mpi%irank,mpi%isize,jspin,dimension%jspd,noco%l_noco,n_bands,n_size)
200
#ifdef CPP_MPI
201
    ! Synchronizes the RMA operations
202 203
    CALL MPI_BARRIER(mpi%mpi_comm,ie) 
#endif
204 205 206

    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
207

208 209
    l_write = input%cdinf.AND.mpi%irank==0

210
    DO n = 1,atoms%ntype
211 212 213

       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))
214
       END DO
215 216 217 218 219 220 221 222 223 224

       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

225
       IF (banddos%l_mcd) CALL mcd_init(atoms,input,dimension,vTot%mt(:,0,:,:),g,f,mcd,n,jspin)
226 227 228

       IF(l_cs) CALL corespec_rme(atoms,input,n,dimension%nstd,&
                                  input%jspins,jspin,results%ef,&
229
                                  dimension%msh,vTot%mt(:,0,:,:),f,g)
230
    END DO
231
    DEALLOCATE (f,g,flo)
232

233 234 235 236 237
    ALLOCATE (we(dimension%neigd))
    n_rank = 0

    ! For k-point paralelization:
    l_evp = .FALSE.
238
    IF (kpts%nkpt < mpi%isize) THEN
239 240 241 242
       l_evp = .TRUE.
       nkpt_extended = kpts%nkpt
       ikptStart = 1
       ikptIncrement = 1
243
    ELSE
244 245 246 247
       ! 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
248
    END IF
249

250 251 252
    DO ikpt = ikptStart, nkpt_extended, ikptIncrement
       i_rec = ikpt
       IF (ikpt.LE.kpts%nkpt) THEN
253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
          !-t3e
          we=0.0
          !--->    determine number of occupied bands and set weights (we)
          noccbd = 0
          DO i = 1,dimension%neigd ! nbands
             we(i) = results%w_iks(n_bands(n_rank)+i,ikpt,jspin)
             IF (noco%l_noco) we(i) = results%w_iks(i,ikpt,1)
             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
268 269

          !
270
          ! -> Gu test: distribute ev's among the processors...
271
          !
Daniel Wortmann's avatar
Daniel Wortmann committed
272
          CALL lapw%init(input,noco, kpts,atoms,sym,ikpt,cell,.false., mpi)
273 274
          skip_t = skip_tt
          IF (l_evp.AND.(mpi%isize.GT.1)) THEN
Gregor Michalicek's avatar
Gregor Michalicek committed
275 276 277 278 279
             noccbd = MERGE(n_bands(1),noccbd,banddos%dos)
             noccbd_l = CEILING(real(noccbd) / mpi%isize)
             n_start = mpi%irank*noccbd_l + 1
             n_end   = min( (mpi%irank+1)*noccbd_l , noccbd )

280 281 282 283 284 285
             noccbd = n_end - n_start + 1
             IF (noccbd<1) THEN
                noccbd=0
             ELSE
                we(1:noccbd) = we(n_start:n_end)
             END IF
286 287 288 289

             IF (n_start > skip_tt) skip_t  = 0
             IF (n_end <= skip_tt) skip_t = noccbd
             IF ((n_start <= skip_tt).AND.(n_end > skip_tt)) skip_t = mod(skip_tt,noccbd)
290
          ELSE
291 292 293 294 295 296 297 298 299 300
             n_start = 1
             IF (banddos%dos) THEN
                noccbd_l = n_bands(1)
                n_end    = n_bands(1)
                noccbd   = n_bands(1)
             ELSE
                noccbd_l = noccbd
                n_end    = noccbd
             END IF
          END IF
301

Gregor Michalicek's avatar
Gregor Michalicek committed
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,dimension%neigd)
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
          ! Synchronizes the RMA operations
309
          CALL MPI_BARRIER(mpi%mpi_comm,ie)
310
#endif
311 312 313 314 315
          !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
316
             CALL nstm3(sym,atoms,vacuum,stars,ikpt,lapw%nv(jspin),input,jspin,kpts,&
317
                        cell,kpts%wtkpt(ikpt),lapw%k1(:,jspin),lapw%k2(:,jspin),&
318
                        enpara%evac0(1,jspin),vTot%vacz(:,:,jspin),gvac1d,gvac2d)
319
          END IF
320

321
          IF (noccbd.EQ.0) GO TO 199
322

323 324 325 326 327 328 329 330 331
          !--->    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
          IF (sliceplot%slice) THEN
             IF (mpi%irank==0) WRITE (16,FMT=*) 'NNNE',sliceplot%nnne
             IF (mpi%irank==0) WRITE (16,FMT=*) 'sliceplot%kk',sliceplot%kk
             nslibd = 0
332
             IF (input%pallst) we(:nbands) = kpts%wtkpt(ikpt)
333 334 335 336 337 338 339 340 341 342
             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)
343
                      if (zmat%l_real) THEN
344
                         zMat%z_r(:,nslibd) = zMat%z_r(:,i)
345
                      else
346
                         zMat%z_c(:,nslibd) = zMat%z_c(:,i)
347
                      endif
348 349 350 351 352 353 354
                   END IF
                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.',&
355
                        sliceplot%nnne,eig(sliceplot%nnne)
356 357 358
                   nslibd = nslibd + 1
                   eig(nslibd) = eig(sliceplot%nnne)
                   we(nslibd) = we(sliceplot%nnne)
359
                   if (zmat%l_real) Then
360
                      zMat%z_r(:,nslibd) = zMat%z_r(:,sliceplot%nnne)
361
                   else
362
                      zMat%z_c(:,nslibd) = zMat%z_c(:,sliceplot%nnne)
363
                   endif
364 365 366
                ELSE
                   DO i = 1,nbands
                      IF (eig(i).GE.sliceplot%e1s .AND. eig(i).LE.sliceplot%e2s) THEN
367
                         nslibd = nslibd + 1
368 369
                         eig(nslibd) = eig(i)
                         we(nslibd) = we(i)
370
                         if (zmat%l_real) THEN
371
                            zMat%z_r(:,nslibd) = zMat%z_r(:,i)
372
                         else
373
                            zMat%z_c(:,nslibd) = zMat%z_c(:,i)
374
                         endif
375
                      END IF
376 377 378 379 380 381 382 383 384 385 386 387 388 389
                   END DO
                   IF (mpi%irank==0) WRITE (16,FMT='(a,i3)')' eigenvalues in sliceplot%slice:',nslibd
                END IF
             END IF
             noccbd = nslibd
             IF (nslibd.EQ.0) GO TO 199 !200
          END IF ! sliceplot%slice

          !--->    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
390
          we(:noccbd) = 2.0 * we(:noccbd) / input%jspins
391 392 393 394 395 396 397 398

          !---> pk non-collinear
          !--->    valence density in the interstitial and vacuum region
          !--->    has to be called only once (if jspin=1) in the non-collinear
          !--->    case
          !     ----> valence density in the interstitial region
          IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
             CALL timestart("cdnval: pwden")
399
             CALL pwden(stars,kpts,banddos,oneD,input,mpi,noco,cell,atoms,sym,ikpt,&
400
                        jspin,lapw,noccbd,we,eig,den,regCharges%qis,results,force%f_b8,zMat)
401 402 403 404 405 406
             CALL timestop("cdnval: pwden")
          END IF
          !+new
          !--->    charge of each valence state in this k-point of the SBZ
          !--->    in the layer interstitial region of the film
          !
407
          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 411 412 413 414 415 416
             END IF
          END IF
          !-new c
          !--->    valence density in the vacuum region
          IF (input%film) THEN
             IF (.NOT.((jspin.EQ.2) .AND. noco%l_noco)) THEN
                CALL timestart("cdnval: vacden")
417
                CALL vacden(vacuum,dimension,stars,oneD, kpts,input, cell,atoms,noco,banddos,&
418
                            gvac1d,gvac2d, we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,lapw,enpara%evac0,eig,&
419
                            den,regCharges%qvac,regCharges%qvlay,regCharges%qstars,zMat)
420 421 422 423 424 425
                CALL timestop("cdnval: vacden")
             END IF
             !--->       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
426 427
                   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 429 430
                END DO
             END DO
          END IF
431

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

          DO ispin = jsp_start,jsp_end
             IF (input%l_f) THEN
437
                CALL force%init2(noccbd,input,atoms)
438
             END IF
439 440 441 442 443
             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")
444

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

449 450 451 452
             !--->       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
453
                CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,&
454
                            skip_t,l_evp,eigVecCoeffs,usdus,regCharges,mcd,banddos%l_mcd)
455 456

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

461 462 463
             !--->    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
464
                CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
465 466 467

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

471
                CALL orb_comp(ispin,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp)
472 473 474 475
             END IF
             !-new
             !--->          set up coefficients for the spherical and
             CALL timestart("cdnval: rhomt")
476
             CALL rhomt(atoms,we,noccbd,eigVecCoeffs,denCoeffs,ispin)
477
             CALL timestop("cdnval: rhomt")
478 479

             IF (noco%l_soc) CALL orbmom(atoms,noccbd,we,ispin,eigVecCoeffs,orb)
480 481
             !--->          non-spherical m.t. density
             CALL timestart("cdnval: rhonmt")
482
             CALL rhonmt(atoms,sphhar,we,noccbd,sym,eigVecCoeffs,denCoeffs,ispin)
483 484 485 486 487 488
             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")
489
             CALL rhomtlo(atoms,noccbd,we,eigVecCoeffs,denCoeffs,ispin)
490

491
             CALL rhonmtlo(atoms,sphhar,noccbd,we,eigVecCoeffs,denCoeffs,ispin)
492 493 494 495
             CALL timestop("cdnval: rho(n)mtlo")

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

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

          IF (noco%l_mperp) THEN
512
             CALL rhomt21(atoms,we,noccbd,eigVecCoeffs,denCoeffsOffdiag)
513
             IF (l_fmpl) CALL rhonmt21(atoms,sphhar,we,noccbd,sym,eigVecCoeffs,denCoeffsOffdiag)
514 515 516 517 518 519 520 521 522 523 524 525
          END IF

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

             !--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
             IF (banddos%ndir.GT.0) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
528
                CALL sympsi(lapw%bkpt,lapw%nv(jspin),lapw%k1(:,jspin),lapw%k2(:,jspin),&
529
                            lapw%k3(:,jspin),sym,dimension,nbands,cell,eig,noco, ksym,jsym,zMat)
530
             END IF
531

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

537 538 539 540 541
             CALL timestop("cdnval: write_info")
             !-new_sl
          END IF

          !--->  end of loop over PE's
542
        ELSE !(ikpt.LE.nkpt)
543 544 545 546
#ifdef CPP_MPI
          ! Synchronizes the RMA operations
          CALL MPI_BARRIER(mpi%mpi_comm,ie)
#endif
547
       END IF
548
    END DO !---> end of k-point loop
549
    DEALLOCATE (we)
550
    !+t3e
551
#ifdef CPP_MPI
552 553 554
    CALL timestart("cdnval: mpi_col_den")
    DO ispin = jsp_start,jsp_end
       CALL mpi_col_den(mpi,sphhar,atoms,oneD,stars,vacuum,&
555
                        input,noco,l_fmpl,ispin,llpd, den%vacxy(1,1,1,ispin),&
556 557
                        den%vacz(1,1,ispin),den%pw(1,ispin),regCharges,&
                        results,denCoeffs,orb,denCoeffsOffdiag,den,den%mmpMat(:,:,:,jspin))
558 559
    END DO
    CALL timestop("cdnval: mpi_col_den")
560
#endif
561 562 563

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

564 565 566
    IF (((jspin.eq.input%jspins).OR.noco%l_mperp) .AND. (banddos%dos.or.banddos%vacdos.or.input%cdinf) ) THEN
       CALL timestart("cdnval: dos")
       IF (mpi%irank==0) THEN
567
          CALL doswrite(eig_id,dimension,kpts,atoms,vacuum,input,banddos,&
568
                        sliceplot,noco,sym,cell,mcd,results,slab%nsld,oneD)
569
          IF (banddos%dos.AND.(banddos%ndir.EQ.-3)) THEN
570
             CALL Ek_write_sl(eig_id,dimension,kpts,atoms,vacuum,input,jspin,sym,cell,slab)
571 572
          END IF
       END IF
573
#ifdef CPP_MPI                
574
       CALL MPI_BARRIER(mpi%mpi_comm,ie)
575
#endif
576 577 578 579
       CALL timestop("cdnval: dos")
    END IF

    IF (mpi%irank==0) THEN
580
       CALL cdnmt(dimension%jspd,atoms,sphhar,llpd,noco,l_fmpl,jsp_start,jsp_end,&
581
                  enpara%el0,enpara%ello0,vTot%mt(:,0,:,:),denCoeffs,&
582
                  usdus,orb,denCoeffsOffdiag,moments,den%mt)
583 584

       DO ispin = jsp_start,jsp_end
585 586
          IF (.NOT.sliceplot%slice) THEN
             DO n=1,atoms%ntype
587 588
                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)
589
             END DO
590
             IF (input%film) enpara%evac1(:vacuum%nvac,ispin)=regCharges%pvac(:vacuum%nvac,ispin)/regCharges%svac(:vacuum%nvac,ispin)
591 592 593 594 595 596 597
          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)
598
             CALL checkDOPAll(input,dimension,sphhar,stars,atoms,sym,vacuum,oneD,cell,den,ispin)
599 600 601 602 603 604
             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")
605
             CALL force_a8(input,atoms,sphhar,ispin,vTot%mt(:,:,:,ispin),den%mt,force,results)
606 607
             CALL timestop("cdnval: force_a8")
          END IF
608
          !-for
609
       END DO ! end of loop ispin = jsp_start,jsp_end
610
       CALL closeXMLElement('mtCharges')
611 612 613 614 615 616 617 618 619 620

       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

621
    END IF ! end of (mpi%irank==0)
622

623 624 625 626
    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
627

628 629
  END SUBROUTINE cdnval
END MODULE m_cdnval