LCOV - code coverage report
Current view: top level - cdn - vacden.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 143 695 20.6 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_vacden
       2             :   USE m_juDFT
       3             :   !     *************************************************************
       4             :   !     determines the 2-d star function expansion coefficients of
       5             :   !     vacuum charge density. speed up by r. wu 1992
       6             :   !     *************************************************************
       7             : CONTAINS
       8          28 :   SUBROUTINE vacden(vacuum,DIMENSION,stars,oneD,kpts,input,sym,cell,atoms,noco,banddos,&
       9          28 :                     gVacMap,we,ikpt,jspin,vz,ne,lapw,evac,eig,den,zMat,dos)
      10             : 
      11             :     !***********************************************************************
      12             :     !     ****** change vacden(....,q) for vacuum density of states shz Jan.96
      13             :     !     ****** change vacden(......,dos%qstars) for starcoefficients, shz. Jan.99
      14             :     !     ****** changed for fleur dw
      15             :     !     In non-collinear calculations the density becomes a hermitian 2x2
      16             :     !     matrix. This subroutine generates this density matrix in the 
      17             :     !     vacuum region. The diagonal elements of this matrix (n_11 & n_22)
      18             :     !     are store in den%vacz and den%vacxy, while the real and imaginary part
      19             :     !     of the off-diagonal element are stored in den%vacz(:,:,3:4) and den%vacxy(:,:,:,3). 
      20             :     !
      21             :     !     Philipp Kurz 99/07
      22             :     !***********************************************************************
      23             : 
      24             :     !******** ABBREVIATIONS ************************************************
      25             :     !     qvac     : vacuum charge of each eigenstate, needed in in cdnval
      26             :     !                to determine the vacuum energy parameters
      27             :     !     vz       : non-warping part of the vacuum potential (matrix)
      28             :     !                collinear    : 2. index = ivac (# of vaccum)
      29             :     !                non-collinear: 2. index = ipot (comp. of pot. matr.)
      30             :     !     den%vacz : non-warping part of the vacuum density matrix,
      31             :     !                diagonal elements n_11 and n_22
      32             :     !     den%vacxy: warping part of the vacuum density matrix,
      33             :     !                diagonal elements n_11 and n_22
      34             :     !     den%vacz(:,:,3:4): non-warping part of the vacuum density matrix,
      35             :     !                off-diagonal elements n_21 (real part in (:,:,3), imaginary part in (:,:,4))
      36             :     !     den%vacxy(:,:,:,3): warping part of the vacuum density matrix,
      37             :     !                off-diagonal elements n_21
      38             :     !***********************************************************************
      39             :     !
      40             :     USE m_constants
      41             :     USE m_grdchlh
      42             :     USE m_qsf
      43             :     USE m_cylbes
      44             :     USE m_dcylbs
      45             :     USE m_od_abvac
      46             :     USE m_vacuz
      47             :     USE m_vacudz
      48             :     USE m_types
      49             :     IMPLICIT NONE
      50             :     TYPE(t_lapw),INTENT(INOUT)    :: lapw !for some reason the second spin data is reset in noco case
      51             :     TYPE(t_dimension),INTENT(IN)  :: DIMENSION
      52             :     TYPE(t_oneD),INTENT(IN)       :: oneD
      53             :     TYPE(t_banddos),INTENT(IN)    :: banddos
      54             :     TYPE(t_input),INTENT(IN)      :: input
      55             :     TYPE(t_vacuum),INTENT(IN)     :: vacuum
      56             :     TYPE(t_noco),INTENT(IN)       :: noco
      57             :     TYPE(t_stars),INTENT(IN)      :: stars
      58             :     TYPE(t_sym),INTENT(IN)        :: sym
      59             :     TYPE(t_cell),INTENT(IN)       :: cell
      60             :     TYPE(t_kpts),INTENT(IN)       :: kpts
      61             :     TYPE(t_atoms),INTENT(IN)      :: atoms
      62             :     TYPE(t_mat),INTENT(IN)        :: zMat
      63             :     TYPE(t_gVacMap),INTENT(IN)    :: gVacMap
      64             :     TYPE(t_potden),INTENT(INOUT)  :: den
      65             :     TYPE(t_dos),   INTENT(INOUT)  :: dos
      66             :     !     .. Scalar Arguments ..
      67             :     INTEGER, INTENT (IN) :: jspin      
      68             :     INTEGER, INTENT (IN) :: ne    
      69             :     INTEGER, INTENT (IN) :: ikpt
      70             :     INTEGER,PARAMETER    :: n2max=13
      71             :     REAL,PARAMETER        :: emax=2.0/hartree_to_ev_const
      72             :     !     .. Array Arguments ..
      73             :     REAL,    INTENT(IN)     :: evac(2,input%jspins)
      74             :     REAL,    INTENT(IN)     :: we(DIMENSION%neigd)
      75             :     REAL                    :: vz(vacuum%nmzd,2) ! Note this breaks the INTENT(IN) from cdnval. It may be read from a file in this subroutine.
      76             :     !     STM-Arguments
      77             :     REAL,    INTENT (IN)    :: eig(DIMENSION%neigd)
      78             :     !     local STM variables
      79          56 :     INTEGER nv2(input%jspins)
      80         112 :     INTEGER kvac1(DIMENSION%nv2d,input%jspins),kvac2(DIMENSION%nv2d,input%jspins),map2(DIMENSION%nvd,input%jspins)
      81          56 :     INTEGER kvac3(DIMENSION%nv2d,input%jspins),map1(DIMENSION%nvd,input%jspins)
      82          56 :     INTEGER mapg2k(DIMENSION%nv2d)
      83             :     !     .. Local Scalars ..
      84             :     COMPLEX aa,ab,av,ba,bb,bv,t1,aae,bbe,abe,bae,aaee,bbee,abee,baee,&
      85             :          &     factorx,factory,c_1,aa_1,ab_1,ba_1,bb_1,ic,av_1,bv_1,d,tempCmplx
      86             : 
      87             :     REAL arg,const,ddui,dduj,dduei,dduej,eps,ev,evacp,phs,phsp,qout,&
      88             :          &     scale,sign,uei,uej,ui,uj,wronk,zks,RESULT(1),ui2,uei2,&
      89             :          &     k_diff,k_d1,k_d2,ui_1,uei_1,uj_1,uej_1,wronk_1
      90             : 
      91             :     INTEGER i,ii,i1,i2,i3,ig3,ig3p,ik,ind2,ind2p,istar ,m1,m3,&
      92             :          &        ivac,j,jj,jz,k,l,ll,l1,n,n2,ispin,kspin,jsp_start,jsp_end,&
      93             :          &        ipot,ie,imz,isp_start,isp_end,&
      94             :          &        ind1,ind1p,irec2,irec3,m
      95             :     !
      96             :     !     .. Local Arrays ..
      97             :     REAL qssbti(3,2),qssbtii,vz0(2)
      98          56 :     REAL bess(-oneD%odi%mb:oneD%odi%mb),dbss(-oneD%odi%mb:oneD%odi%mb)
      99          56 :     COMPLEX, ALLOCATABLE :: ac(:,:,:),bc(:,:,:)
     100         140 :     REAL,    ALLOCATABLE :: dt(:),dte(:),du(:),ddu(:,:),due(:)
     101         140 :     REAL,    ALLOCATABLE :: ddue(:,:),t(:),te(:),tei(:,:),dummy(:)
     102          84 :     REAL,    ALLOCATABLE :: u(:,:,:),ue(:,:,:),v(:),yy(:)
     103             :     !-odim
     104          56 :     COMPLEX, ALLOCATABLE :: ac_1(:,:,:,:),bc_1(:,:,:,:)
     105          56 :     REAL,    ALLOCATABLE :: dt_1(:,:),dte_1(:,:)
     106          84 :     REAL,    ALLOCATABLE :: du_1(:,:),ddu_1(:,:,:),due_1(:,:)
     107          56 :     REAL,    ALLOCATABLE :: ddue_1(:,:,:),t_1(:,:)
     108          56 :     REAL,    ALLOCATABLE :: tei_1(:,:,:),te_1(:,:)
     109          56 :     REAL,    ALLOCATABLE :: u_1(:,:,:,:),ue_1(:,:,:,:)
     110             :     !+odim
     111             :     !     ..
     112             :     !     ..
     113             : 
     114             :     !     *******************************************************************************
     115             :     !
     116             : 
     117             :     !
     118             :     !    layers: no. of layers to be calculated (in vertical direction with z-values as given by izlay)
     119             :     !    izlay : defines vertical position of layers in delz (=0.1 a.u.) units from begining of vacuum region
     120             :     !    vacdos: =T: calculate vacuum dos in layers as given by the above
     121             :     !    integ : =T: integrate in vertical position between izlay(layer,1)..izlay(layer,2)
     122             :     !    nstm  : 0: s-Tip, 1: p_z-Tip, 2: d_z^2-Tip (following Chen's derivative rule) ->rhzgrd.f is used 
     123             :     !                 to calculate derivatives 
     124             :     !    tworkf: Workfunction of Tip (in hartree units) is needed for d_z^2-Orbital)
     125             :     !    starcoeff: =T: star coefficients are calculated at values of izlay for 0th (=q) to nstars-1. star 
     126             :     !                (dos%qstars(1..nstars-1))
     127             :     !    nstars: number of star functions to be used (0th star is given by value of q=charge integrated in 2D) 
     128             :     !
     129             :     !    further possibility: (readin of locx, locy has to be implemented in flapw7.f or they have to be set explicitly)
     130             :     !
     131             :     !     locx and locy can be used to calculate local DOS at a certain vertical position z (or integrated in z)
     132             :     !     within a restricted area of the 2D unit cell, the corners of this area is given by locx and locy
     133             :     !     they are defined in internal coordinates, i.e. \vec{r}_1=locx(1)*\vec{a}_1+locy(1)*\vec{a}_2
     134             :     !                                                    \vec{r}_2=locx(2)*\vec{a}_1+locy(2)*\vec{a}_2
     135             :     !                 \vec{a}_1,2 are the 2D lattice vectors           
     136             :     !  
     137             :     !     **************************************************************************************************
     138             : 
     139          28 :     CALL timestart("vacden")
     140             : 
     141             :     ALLOCATE ( ac(DIMENSION%nv2d,DIMENSION%neigd,input%jspins),bc(DIMENSION%nv2d,DIMENSION%neigd,input%jspins),dt(DIMENSION%nv2d),&
     142             :          &           dte(DIMENSION%nv2d),du(vacuum%nmzd),ddu(vacuum%nmzd,DIMENSION%nv2d),due(vacuum%nmzd),&
     143             :          &           ddue(vacuum%nmzd,DIMENSION%nv2d),t(DIMENSION%nv2d),te(DIMENSION%nv2d),&
     144             :          &           tei(DIMENSION%nv2d,input%jspins),u(vacuum%nmzd,DIMENSION%nv2d,input%jspins),ue(vacuum%nmzd,DIMENSION%nv2d,input%jspins),&
     145          28 :          &           v(3),yy(vacuum%nmzd))
     146          28 :     IF (oneD%odi%d1) THEN
     147             :        ALLOCATE (      ac_1(DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb,DIMENSION%neigd,input%jspins),&
     148             :             &                  bc_1(DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb,DIMENSION%neigd,input%jspins),&
     149             :             &                  dt_1(DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb),&
     150             :             &                 dte_1(DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb),&
     151             :             &                  du_1(vacuum%nmzd,-oneD%odi%mb:oneD%odi%mb),&
     152             :             &            ddu_1(vacuum%nmzd,DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb),&
     153             :             &                 due_1(vacuum%nmzd,-oneD%odi%mb:oneD%odi%mb),&
     154             :             &           ddue_1(vacuum%nmzd,DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb),&
     155             :             &                   t_1(DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb),&
     156             :             &                  te_1(DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb),&
     157             :             &                 tei_1(DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb,input%jspins),&
     158             :             &              u_1(vacuum%nmzd,DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb,input%jspins),&
     159           0 :             &             ue_1(vacuum%nmzd,DIMENSION%nv2d,-oneD%odi%mb:oneD%odi%mb,input%jspins) )
     160             :     END IF ! oneD%odi%d1
     161             :     !
     162             : 
     163          28 :     vz0(:) = vz(vacuum%nmz,:)
     164          28 :     eps=0.01
     165          28 :     ic = CMPLX(0.,1.)
     166             :     !    ------------------
     167             :    
     168             :     !     -----> set up mapping arrays
     169          28 :     IF (noco%l_ss) THEN
     170             :        jsp_start = 1
     171             :        jsp_end   = 2
     172             :     ELSE
     173          28 :        jsp_start = jspin
     174          28 :        jsp_end   = jspin
     175             :     ENDIF
     176          56 :     DO ispin = jsp_start,jsp_end
     177          56 :        IF (oneD%odi%d1) THEN
     178           0 :           n2 = 0
     179           0 :           k_loop:DO  k = 1,lapw%nv(ispin)
     180           0 :              DO  j = 1,n2
     181           0 :                 IF (lapw%k3(k,ispin).EQ.kvac3(j,ispin)) THEN
     182           0 :                    map1(k,ispin) = j
     183           0 :                    CYCLE k_loop
     184             :                 END IF
     185             :              ENDDO
     186           0 :              n2 = n2 + 1
     187           0 :              IF (n2>DIMENSION%nv2d)  CALL juDFT_error("vacden0",calledby ="vacden")
     188           0 :              kvac3(n2,ispin) =  lapw%k3(k,ispin)
     189           0 :              map1(k,ispin) = n2
     190             :           ENDDO k_loop
     191           0 :           nv2(ispin) = n2
     192             :        ELSE
     193          28 :           n2 = 0
     194       11486 :           k_loop2:DO  k = 1,lapw%nv(ispin)
     195      380566 :              DO  j = 1,n2
     196      194918 :                 IF ( lapw%gvec(1,k,ispin).EQ.kvac1(j,ispin) .AND.&
     197        1094 :                      lapw%gvec(2,k,ispin).EQ.kvac2(j,ispin) ) THEN
     198       10364 :                    map2(k,ispin) = j
     199       10364 :                    CYCLE k_loop2
     200             :                 END IF
     201             :              ENDDO
     202        1094 :              n2 = n2 + 1
     203        1094 :              IF (n2>DIMENSION%nv2d)  CALL juDFT_error("vacden0","vacden")
     204        1094 :              kvac1(n2,ispin) = lapw%gvec(1,k,ispin)
     205        1094 :              kvac2(n2,ispin) = lapw%gvec(2,k,ispin)
     206        1122 :              map2(k,ispin) = n2
     207             :           ENDDO k_loop2
     208          28 :           nv2(ispin) = n2
     209             :        END IF
     210             :     ENDDO
     211          28 :     IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) THEN
     212           0 :        lapw%nv(2)  = lapw%nv(1)
     213           0 :        nv2(2) = nv2(1)
     214           0 :        DO k = 1,nv2(1)
     215           0 :           kvac1(k,2) = kvac1(k,1)
     216           0 :           kvac2(k,2) = kvac2(k,1)
     217           0 :           IF(oneD%odi%d1) kvac3(k,2) =  kvac3(k,1)
     218             :        ENDDO
     219           0 :        DO k = 1,lapw%nv(1)
     220           0 :           lapw%k3(k,2) = lapw%k3(k,1)
     221           0 :           map2(k,2) = map2(k,1)
     222           0 :           IF(oneD%odi%d1)THEN
     223           0 :              lapw%k1(k,2) = lapw%k1(k,1)
     224           0 :              lapw%k2(k,2) = lapw%k2(k,1)
     225           0 :              map1(k,2) = map1(k,1)
     226             :           ENDIF
     227             :        ENDDO
     228             :     ENDIF
     229             : 
     230             :     !+dw
     231             :     !    if tunneling current should be calculated we need to write out 
     232             :     !     info on electronic structure: --> mapping from kvac to gvac by mapg2k
     233             :     !                                             shz, Jan.99
     234          28 :     IF (vacuum%nstm.EQ.3) THEN
     235           0 :        DO j=1, n2max
     236           0 :           mapg2k(j)=j
     237           0 :           DO i=1, nv2(jspin)
     238           0 :              IF (kvac1(i,jspin).EQ.gVacMap%gvac1d(j).AND.kvac2(i,jspin).EQ.gVacMap%gvac2d(j)) mapg2k(j)=i
     239             :           END DO
     240             :        END DO
     241             :     END IF
     242             :     !
     243             :     !-dw
     244             :    
     245             : 
     246          28 :     wronk = 2.0
     247          28 :     const = 1.0 / ( SQRT(cell%omtil)*wronk )
     248             :     !-odim
     249          28 :     IF (oneD%odi%d1) THEN
     250           0 :        ac_1(:,:,:,:) = CMPLX(0.0,0.0)
     251           0 :        bc_1(:,:,:,:) = CMPLX(0.0,0.0)
     252             :     END IF
     253             :     !+odim
     254          56 :     DO  ivac = 1,vacuum%nvac
     255          76 :        ac(:,:,:) = CMPLX(0.0,0.0)
     256          76 :        bc(:,:,:) = CMPLX(0.0,0.0)
     257          28 :        sign = 3. - 2.*ivac
     258             :       
     259          28 :        IF (noco%l_noco) THEN
     260             :           !--->    In a non-collinear calculation vacden is only called once.
     261             :           !--->    Thus, the vaccum wavefunctions and the A- and B-coeff. (ac bc)
     262             :           !--->    have to be calculated for both spins on that call.
     263             :           !--->       setup the spin-spiral q-vector
     264           0 :           qssbti(1,1) = - noco%qss(1)/2
     265           0 :           qssbti(2,1) = - noco%qss(2)/2
     266           0 :           qssbti(1,2) = + noco%qss(1)/2
     267           0 :           qssbti(2,2) = + noco%qss(2)/2
     268           0 :           qssbti(3,1) = - noco%qss(3)/2
     269           0 :           qssbti(3,2) = + noco%qss(3)/2
     270           0 :           DO ispin = 1,input%jspins
     271             :              !     -----> set up vacuum wave functions
     272           0 :              IF (oneD%odi%d1) THEN
     273             :                 CALL od_abvac(&
     274             :                      cell,vacuum,DIMENSION,stars,&
     275             :                      oneD,qssbti(3,ispin),&
     276             :                      oneD%odi%n2d,&
     277             :                      wronk,evacp,lapw%bkpt,oneD%odi%M,oneD%odi%mb,&
     278             :                      vz(1,ispin),kvac3(1,ispin),nv2(ispin),&
     279             :                      t_1(1,-oneD%odi%mb),dt_1(1,-oneD%odi%mb),u_1(1,1,-oneD%odi%mb,ispin),&
     280             :                      te_1(1,-oneD%odi%mb),dte_1(1,-oneD%odi%mb),&
     281             :                      tei_1(1,-oneD%odi%mb,ispin),&
     282           0 :                      ue_1(1,1,-oneD%odi%mb,ispin))
     283           0 :                 DO k = 1,lapw%nv(ispin)
     284           0 :                    kspin = (lapw%nv(1)+atoms%nlotot)*(ispin-1) + k
     285           0 :                    l = map1(k,ispin) 
     286           0 :                    irec3 = stars%ig(lapw%gvec(1,k,ispin),lapw%gvec(2,k,ispin),lapw%gvec(3,k,ispin))
     287           0 :                    IF (irec3.NE.0) THEN
     288           0 :                       irec2 = stars%ig2(irec3)
     289           0 :                       zks = stars%sk2(irec2)*cell%z1
     290           0 :                       arg = stars%phi2(irec2)
     291           0 :                       CALL cylbes(oneD%odi%mb,zks,bess)
     292           0 :                       CALL dcylbs(oneD%odi%mb,zks,bess,dbss)
     293           0 :                       DO m = -oneD%odi%mb,oneD%odi%mb
     294             :                          wronk_1 = t_1(l,m)*dte_1(l,m) -&
     295           0 :                               te_1(l,m)*dt_1(l,m)
     296             :                          av_1 = EXP(-CMPLX(0.0,m*arg))*(ic**m)*&
     297             :                               CMPLX(dte_1(l,m)*bess(m) -&
     298             :                               te_1(l,m)*stars%sk2(irec2)*dbss(m),0.0)/&
     299           0 :                               ((wronk_1)*SQRT(cell%omtil))
     300             :                          bv_1 = EXP(-CMPLX(0.0,m*arg))*(ic**m)*&
     301             :                               CMPLX(-dt_1(l,m)*bess(m) +&
     302             :                               t_1(l,m)*stars%sk2(irec2)*dbss(m),0.0)/&
     303           0 :                               ((wronk_1)*SQRT(cell%omtil))
     304           0 :                          IF (zmat%l_real) THEN
     305           0 :                             ac_1(l,m,:ne,ispin) = ac_1(l,m,:ne,ispin) + zMat%data_r(kspin,:ne)*av_1
     306           0 :                             bc_1(l,m,:ne,ispin) = bc_1(l,m,:ne,ispin) + zMat%data_r(kspin,:ne)*bv_1
     307             :                          ELSE
     308           0 :                             ac_1(l,m,:ne,ispin) = ac_1(l,m,:ne,ispin) + zMat%data_c(kspin,:ne)*av_1
     309           0 :                             bc_1(l,m,:ne,ispin) = bc_1(l,m,:ne,ispin) + zMat%data_c(kspin,:ne)*bv_1
     310             :                          END IF
     311             :                       END DO      ! -mb:mb
     312             :                    END IF
     313             :                 END DO
     314             :              ELSE ! 1-dimensional
     315           0 :                 vz0(ispin) = vz(vacuum%nmz,ispin)
     316           0 :                 evacp = evac(ivac,ispin)
     317           0 :                 DO ik = 1,nv2(ispin)
     318           0 :                    v(1) = lapw%bkpt(1) + kvac1(ik,ispin) + qssbti(1,ispin)
     319           0 :                    v(2) = lapw%bkpt(2) + kvac2(ik,ispin) + qssbti(2,ispin)
     320           0 :                    v(3) = 0.
     321           0 :                    ev = evacp - 0.5*DOT_PRODUCT(v,MATMUL(v,cell%bbmat))
     322             :                    CALL vacuz(ev,vz(1,ispin),vz0(ispin),vacuum%nmz,vacuum%delz,t(ik),&
     323           0 :                         dt(ik),u(1,ik,ispin))
     324             :                    CALL vacudz(ev,vz(1,ispin),vz0(ispin),vacuum%nmz,vacuum%delz,te(ik),&
     325             :                         dte(ik),tei(ik,ispin),ue(1,ik,ispin),dt(ik),&
     326           0 :                         u(1,ik,ispin))
     327           0 :                    scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
     328           0 :                    te(ik) = scale*te(ik)
     329           0 :                    dte(ik) = scale*dte(ik)
     330           0 :                    tei(ik,ispin) = scale*tei(ik,ispin)
     331           0 :                    DO j = 1,vacuum%nmz
     332           0 :                       ue(j,ik,ispin) = scale*ue(j,ik,ispin)
     333             :                    ENDDO
     334             :                 ENDDO
     335             :                 !     -----> construct a and b coefficients
     336           0 :                 DO k = 1,lapw%nv(ispin)
     337             :                    !--->          the coefficients of the spin-down basis functions are
     338             :                    !--->          stored in the second half of the eigenvector
     339           0 :                    kspin = (lapw%nv(1)+atoms%nlotot)*(ispin-1) + k
     340           0 :                    l = map2(k,ispin)
     341           0 :                    zks = lapw%k3(k,ispin)*cell%bmat(3,3)*sign
     342           0 :                    arg = zks*cell%z1
     343           0 :                    c_1 = CMPLX(COS(arg),SIN(arg)) * const
     344           0 :                    av = -c_1 * CMPLX( dte(l),zks*te(l) ) 
     345           0 :                    bv =  c_1 * CMPLX(  dt(l),zks* t(l) ) 
     346             :                    !     -----> loop over basis functions
     347           0 :                    IF (zmat%l_real) THEN
     348           0 :                       ac(l,:ne,ispin) = ac(l,:ne,ispin) + zMat%data_r(kspin,:ne)*av
     349           0 :                       bc(l,:ne,ispin) = bc(l,:ne,ispin) + zMat%data_r(kspin,:ne)*bv
     350             :                    ELSE
     351           0 :                       ac(l,:ne,ispin) = ac(l,:ne,ispin) + zMat%data_c(kspin,:ne)*av
     352           0 :                       bc(l,:ne,ispin) = bc(l,:ne,ispin) + zMat%data_c(kspin,:ne)*bv
     353             :                    ENDIF
     354             :                 ENDDO
     355             :                 !--->       end of spin loop
     356             :              ENDIF
     357             :              !---Y       end of geometry 1d/film
     358             :           ENDDO
     359             :           !--->       output for testing
     360             :           !            DO k = 1,10
     361             :           !               DO n = 1,5
     362             :           !                  DO ispin = 1,jspins
     363             :           !                     write(*,9000)k,n,ispin,ac(k,n,ispin),bc(k,n,ispin)
     364             :           !                  ENDDO
     365             :           !               ENDDO
     366             :           !            ENDDO
     367             :           ! 9000       FORMAT('k=',i3,' ie=',i3,' isp=',i3,
     368             :           !     +             ' ac= (',e12.6,',',e12.6,')',
     369             :           !     +             ' bc= (',e12.6,',',e12.6,')')
     370             :        ELSE
     371             :           !     -----> set up vacuum wave functions
     372          28 :           IF (oneD%odi%d1) THEN
     373           0 :              qssbtii = 0.
     374           0 :              evacp = evac(ivac,jspin)
     375             :              CALL od_abvac(&
     376             :                   &           cell,vacuum,DIMENSION,stars,&
     377             :                   &           oneD,qssbtii,&
     378             :                   &           oneD%odi%n2d,&
     379             :                   &           wronk,evacp,lapw%bkpt,oneD%odi%M,oneD%odi%mb,&
     380             :                   &           vz(1,ivac),kvac3(1,jspin),nv2(jspin),&
     381             :                   &           t_1(1,-oneD%odi%mb),dt_1(1,-oneD%odi%mb),u_1(1,1,-oneD%odi%mb,jspin),&
     382             :                   &           te_1(1,-oneD%odi%mb),dte_1(1,-oneD%odi%mb),&
     383             :                   &           tei_1(1,-oneD%odi%mb,jspin),&
     384           0 :                   &           ue_1(1,1,-oneD%odi%mb,jspin))
     385           0 :              DO k = 1,lapw%nv(jspin)
     386           0 :                 l = map1(k,jspin)
     387           0 :                 irec3 = stars%ig(lapw%gvec(1,k,jspin),lapw%gvec(2,k,jspin),lapw%gvec(3,k,jspin))
     388           0 :                 IF (irec3.NE.0) THEN
     389           0 :                    irec2 = stars%ig2(irec3)
     390           0 :                    zks = stars%sk2(irec2)*cell%z1
     391           0 :                    arg = stars%phi2(irec2)
     392           0 :                    CALL cylbes(oneD%odi%mb,zks,bess)
     393           0 :                    CALL dcylbs(oneD%odi%mb,zks,bess,dbss)
     394           0 :                    DO m = -oneD%odi%mb,oneD%odi%mb
     395           0 :                       wronk_1 = t_1(l,m)*dte_1(l,m) -te_1(l,m)*dt_1(l,m)
     396             :                       av_1 = EXP(-CMPLX(0.0,m*arg))*(ic**m)*&
     397             :                            CMPLX(dte_1(l,m)*bess(m) -&
     398             :                            te_1(l,m)*stars%sk2(irec2)*dbss(m),0.0)/&
     399           0 :                            ((wronk_1)*SQRT(cell%omtil))
     400             :                       bv_1 = EXP(-CMPLX(0.0,m*arg))*(ic**m)*&
     401             :                            CMPLX(-dt_1(l,m)*bess(m) +&
     402             :                            t_1(l,m)*stars%sk2(irec2)*dbss(m),0.0)/&
     403           0 :                            ((wronk_1)*SQRT(cell%omtil))
     404           0 :                       IF (zmat%l_real) THEN
     405           0 :                          ac_1(l,m,:ne,jspin) = ac_1(l,m,:ne,jspin) + zMat%data_r(k,:ne)*av_1
     406           0 :                          bc_1(l,m,:ne,jspin) = bc_1(l,m,:ne,jspin) + zMat%data_r(k,:ne)*bv_1
     407             :                       ELSE
     408           0 :                          ac_1(l,m,:ne,jspin) = ac_1(l,m,:ne,jspin) + zMat%data_c(k,:ne)*av_1
     409           0 :                          bc_1(l,m,:ne,jspin) = bc_1(l,m,:ne,jspin) + zMat%data_c(k,:ne)*bv_1
     410             :                       ENDIF
     411             :                    END DO      ! -mb:mb
     412             :                 END IF
     413             :              END DO         ! k = 1,lapw%nv
     414             :           ELSE     !oneD%odi%d1
     415          28 :              evacp = evac(ivac,jspin)
     416        1122 :              DO ik = 1,nv2(jspin)
     417        1094 :                 v(1) = lapw%bkpt(1) + kvac1(ik,jspin)
     418        1094 :                 v(2) = lapw%bkpt(2) + kvac2(ik,jspin)
     419        1094 :                 v(3) = 0.
     420        1094 :                 ev = evacp - 0.5*DOT_PRODUCT(v,MATMUL(v,cell%bbmat))
     421        1094 :                 CALL vacuz(ev,vz(1,ivac),vz0(ivac),vacuum%nmz,vacuum%delz,t(ik),dt(ik),u(1,ik,jspin))
     422             :                 CALL vacudz(ev,vz(1,ivac),vz0(ivac),vacuum%nmz,vacuum%delz,te(ik),&
     423             :                      &              dte(ik),tei(ik,jspin),ue(1,ik,jspin),dt(ik),&
     424        1094 :                      &              u(1,ik,jspin))
     425        1094 :                 scale = wronk/ (te(ik)*dt(ik)-dte(ik)*t(ik))
     426        1094 :                 te(ik) = scale*te(ik)
     427        1094 :                 dte(ik) = scale*dte(ik)
     428        1094 :                 tei(ik,jspin) = scale*tei(ik,jspin)
     429      274622 :                 DO j = 1,vacuum%nmz
     430      274594 :                    ue(j,ik,jspin) = scale*ue(j,ik,jspin)
     431             :                 ENDDO
     432             :              ENDDO
     433             :              !     -----> construct a and b coefficients
     434       11486 :              DO k = 1,lapw%nv(jspin)
     435       11458 :                 l = map2(k,jspin)
     436       11458 :                 zks = lapw%k3(k,jspin)*cell%bmat(3,3)*sign
     437       11458 :                 arg = zks*cell%z1
     438       11458 :                 c_1 = CMPLX(COS(arg),SIN(arg)) * const
     439       11458 :                 av = -c_1 * CMPLX( dte(l),zks*te(l) ) 
     440       11458 :                 bv =  c_1 * CMPLX(  dt(l),zks* t(l) ) 
     441             :                 !     -----> loop over basis functions
     442       11486 :                 IF (zmat%l_real) THEN
     443        9386 :                    ac(l,:ne,jspin) = ac(l,:ne,jspin) + zMat%data_r(k,:ne)*av
     444       55559 :                    bc(l,:ne,jspin) = bc(l,:ne,jspin) + zMat%data_r(k,:ne)*bv
     445             :                 ELSE
     446        2072 :                    ac(l,:ne,jspin) = ac(l,:ne,jspin) + zMat%data_c(k,:ne)*av
     447       20720 :                    bc(l,:ne,jspin) = bc(l,:ne,jspin) + zMat%data_c(k,:ne)*bv
     448             :                 ENDIF
     449             :              ENDDO
     450             :           END IF ! D1
     451             :        ENDIF
     452             :        !
     453             :        !   ----> calculate first and second derivative of u,ue 
     454             :        !        in order to simulate p_z or d_z^2 Tip in Chen's model , shz. 97
     455             :        ! 
     456          28 :        IF (vacuum%nstm.GT.0) THEN
     457           0 :           DO  ik = 1,nv2(jspin)
     458             :              !               CALL rhzgrd(nmz,delz,u(1,ik,jspin),4,du,ddu(1,ik))
     459             :              !               CALL rhzgrd(nmz,delz,ue(1,ik,jspin),4,due,ddue(1,ik))   
     460             : 
     461           0 :              ALLOCATE ( dummy(vacuum%nmz) )
     462             :              CALL grdchlh(&
     463             :                   0,1,vacuum%nmz,vacuum%delz,dummy,u(1,ik,jspin),4,&
     464           0 :                   du,ddu(1,ik))
     465             :              CALL grdchlh(&
     466             :                   0,1,vacuum%nmz,vacuum%delz,dummy,ue(1,ik,jspin),4,&
     467           0 :                   due,ddue(1,ik))
     468           0 :              DEALLOCATE ( dummy )
     469             : 
     470           0 :              IF (vacuum%nstm.EQ.1) THEN
     471           0 :                 u(:vacuum%nmz,ik,jspin)=du(:vacuum%nmz)
     472           0 :                 ue(:vacuum%nmz,ik,jspin)=due(:vacuum%nmz)
     473             :              END IF
     474             :           ENDDO
     475             :        END IF
     476             : 
     477             :        !+dw
     478             : 
     479             :        !
     480             :        !       --> to calculate Tunneling Current between two systems
     481             :        !           within Bardeens Approach one needs ac(l,n), bc(l,n);
     482             :        !           they are written to the file vacwave 
     483             :        !                           IF nstm=3
     484             :        !                              tworkf is then the fermi energy (in hartree)
     485             :        !
     486          28 :        IF (vacuum%nstm.EQ.3) THEN
     487             : #ifdef CPP_MPI
     488           0 :           CALL judft_error("nstm==3 does not work in parallel",calledby="vacden")
     489             : #else
     490             :           i=0
     491             :           DO n = 1, ne
     492             :              IF (ABS(eig(n)-vacuum%tworkf).LE.emax) i=i+1
     493             :           END DO
     494             :           WRITE (87,FMT=990) lapw%bkpt(1),lapw%bkpt(2), i, n2max
     495             :           DO n = 1, ne
     496             :              IF (ABS(eig(n)-vacuum%tworkf).LE.emax) THEN
     497             :                 WRITE (87,FMT=1000) eig(n)
     498             :                 DO j=1,n2max
     499             :                    WRITE (87,FMT=1010) ac(mapg2k(j),n,jspin),&
     500             :                         bc(mapg2k(j),n,jspin)
     501             :                 END DO
     502             :              END IF
     503             :           END DO
     504             : #endif
     505             : 
     506             : 
     507             : 
     508             :        END IF
     509             : 990    FORMAT(2(f8.4,1x),i3,1x,i3)
     510             : 1000   FORMAT(e12.4)
     511             : 1010   FORMAT(2(2e20.8,1x))
     512             :        !
     513             :        !        ------------------------------------------------------------
     514             :        !-dw
     515             :        !
     516             :        !---->   non-warping part of the density (g=0 star)
     517             :        !
     518          28 :        IF (vacuum%nstm.EQ.2) THEN
     519             :           !
     520             :           !  ----> d_z^2-Tip needs: |d^2(psi)/dz^2 - kappa^2/3 psi|^2
     521             :           !
     522           0 :           DO l = 1,nv2(jspin)
     523           0 :              aa = 0.0
     524           0 :              bb = 0.0
     525           0 :              ba = 0.0
     526           0 :              ab = 0.0
     527           0 :              DO n = 1,ne
     528           0 :                 aa = aa + we(n)*CONJG(ac(l,n,jspin))*ac(l,n,jspin)
     529           0 :                 bb = bb + we(n)*CONJG(bc(l,n,jspin))*bc(l,n,jspin)
     530           0 :                 ab = ab + we(n)*CONJG(ac(l,n,jspin))*bc(l,n,jspin)
     531           0 :                 ba = ba + we(n)*CONJG(bc(l,n,jspin))*ac(l,n,jspin)
     532           0 :                 qout = REAL(CONJG(ac(l,n,jspin))*ac(l,n,jspin)+tei(l,jspin)*CONJG(bc(l,n,jspin))*bc(l,n,jspin))
     533           0 :                 dos%qvac(n,ivac,ikpt,jspin) = dos%qvac(n,ivac,ikpt,jspin) + qout*cell%area
     534             :              END DO
     535           0 :              aae=-aa*vacuum%tworkf*2/3
     536           0 :              bbe=-bb*vacuum%tworkf*2/3
     537           0 :              abe=-ab*vacuum%tworkf*2/3
     538           0 :              bae=-ba*vacuum%tworkf*2/3
     539           0 :              aaee=aa*vacuum%tworkf*vacuum%tworkf*4/9
     540           0 :              bbee=bb*vacuum%tworkf*vacuum%tworkf*4/9
     541           0 :              abee=ab*vacuum%tworkf*vacuum%tworkf*4/9
     542           0 :              baee=ba*vacuum%tworkf*vacuum%tworkf*4/9  
     543           0 :              DO  jz = 1,vacuum%nmz
     544           0 :                 ui = u(jz,l,jspin)
     545           0 :                 uei = ue(jz,l,jspin)
     546           0 :                 ddui = ddu(jz,l) 
     547           0 :                 dduei = ddue(jz,l)
     548             :                 den%vacz(jz,ivac,jspin) = den%vacz(jz,ivac,jspin) +&
     549             :                      REAL(aaee*ui*ui+bbee*uei*uei+&
     550             :                      (abee+baee)*ui*uei+aa*ddui*ddui+&
     551             :                      bb*dduei*dduei+(ab+ba)*ddui*dduei+&
     552             :                      2*aae*ui*ddui+2*bbe*uei*dduei+&
     553           0 :                      (abe+bae)*(ui*dduei+uei*ddui))
     554             : 
     555             :              ENDDO
     556             :           END DO
     557             :           !
     558             :           !    -----> s-Tip: |psi|^2 and p-Tip: |d(psi)/dz|^2 
     559             :           !
     560             :        ELSE
     561          28 :           IF (noco%l_noco) THEN
     562             :              !--->          diagonal elements of the density matrix, n_11 and n_22
     563             :              !--->          the non-warping part of n_21 is calculated together with
     564             :              !--->          the warping part of n_21
     565           0 :              DO ispin = 1,input%jspins
     566           0 :                 IF (oneD%odi%d1) THEN
     567           0 :                    DO  l = 1,nv2(ispin)
     568           0 :                       DO  m = -oneD%odi%mb,oneD%odi%mb
     569           0 :                          aa = CMPLX(0.0,0.0)
     570           0 :                          bb = CMPLX(0.0,0.0)
     571           0 :                          ba = CMPLX(0.0,0.0)
     572           0 :                          ab = CMPLX(0.0,0.0)
     573           0 :                          DO  n = 1,ne
     574           0 :                             aa = aa + we(n)*CONJG(ac_1(l,m,n,ispin))*ac_1(l,m,n,ispin)
     575           0 :                             bb = bb + we(n)*CONJG(bc_1(l,m,n,ispin))*bc_1(l,m,n,ispin)
     576           0 :                             ab = ab + we(n)*CONJG(ac_1(l,m,n,ispin))*bc_1(l,m,n,ispin)
     577           0 :                             ba = ba + we(n)*CONJG(bc_1(l,m,n,ispin))*ac_1(l,m,n,ispin)
     578             :                             qout = REAL(CONJG(ac_1(l,m,n,ispin))*ac_1(l,m,n,ispin) +&
     579             :                                  tei_1(l,m,ispin)*CONJG(bc_1(l,m,n,ispin))*&
     580           0 :                                  bc_1(l,m,n,ispin))
     581           0 :                             dos%qvac(n,ivac,ikpt,ispin) = dos%qvac(n,ivac,ikpt,ispin)+qout*cell%area
     582             :                          END DO
     583           0 :                          DO  jz = 1,vacuum%nmz
     584           0 :                             ui = u_1(jz,l,m,ispin)
     585           0 :                             uei = ue_1(jz,l,m,ispin)
     586           0 :                             den%vacz(jz,ivac,ispin) = den%vacz(jz,ivac,ispin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
     587             :                          END DO
     588             :                       END DO
     589             :                    END DO
     590             :                 ELSE 
     591           0 :                    DO l = 1,nv2(ispin)
     592           0 :                       aa = 0.0
     593           0 :                       bb = 0.0
     594           0 :                       ba = 0.0
     595           0 :                       ab = 0.0
     596           0 :                       DO n = 1,ne
     597           0 :                          aa=aa + we(n)*CONJG(ac(l,n,ispin))*ac(l,n,ispin)
     598           0 :                          bb=bb + we(n)*CONJG(bc(l,n,ispin))*bc(l,n,ispin)
     599           0 :                          ab=ab + we(n)*CONJG(ac(l,n,ispin))*bc(l,n,ispin)
     600           0 :                          ba=ba + we(n)*CONJG(bc(l,n,ispin))*ac(l,n,ispin)
     601           0 :                          qout = REAL(CONJG(ac(l,n,ispin))*ac(l,n,ispin)+tei(l,ispin)*CONJG(bc(l,n,ispin))*bc(l,n,ispin))
     602           0 :                          dos%qvac(n,ivac,ikpt,ispin) = dos%qvac(n,ivac,ikpt,ispin) + qout*cell%area
     603             :                       END DO
     604           0 :                       DO jz = 1,vacuum%nmz
     605           0 :                          ui = u(jz,l,ispin)
     606           0 :                          uei = ue(jz,l,ispin)
     607           0 :                          den%vacz(jz,ivac,ispin) = den%vacz(jz,ivac,ispin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
     608             :                       ENDDO
     609             :                    ENDDO
     610             :                 END IF ! one-dimensional
     611             :              ENDDO
     612             :           ELSE
     613          28 :              IF (oneD%odi%d1) THEN
     614           0 :                 DO  l = 1,nv2(jspin)
     615           0 :                    DO  m = -oneD%odi%mb,oneD%odi%mb
     616           0 :                       aa = CMPLX(0.0,0.0)
     617           0 :                       bb = CMPLX(0.0,0.0)
     618           0 :                       ba = CMPLX(0.0,0.0)
     619           0 :                       ab = CMPLX(0.0,0.0)
     620           0 :                       DO  n = 1,ne
     621           0 :                          aa = aa + we(n)*CONJG(ac_1(l,m,n,jspin))*ac_1(l,m,n,jspin)
     622           0 :                          bb = bb + we(n)*CONJG(bc_1(l,m,n,jspin))*bc_1(l,m,n,jspin)
     623           0 :                          ab = ab + we(n)*CONJG(ac_1(l,m,n,jspin))*bc_1(l,m,n,jspin)
     624           0 :                          ba = ba + we(n)*CONJG(bc_1(l,m,n,jspin))*ac_1(l,m,n,jspin)
     625             :                          qout = REAL(CONJG(ac_1(l,m,n,jspin))*ac_1(l,m,n,jspin)+&
     626           0 :                               tei_1(l,m,jspin)*CONJG(bc_1(l,m,n,jspin))*bc_1(l,m,n,jspin))
     627           0 :                          dos%qvac(n,ivac,ikpt,jspin) = dos%qvac(n,ivac,ikpt,jspin)+qout*cell%area
     628             :                       END DO
     629           0 :                       DO  jz = 1,vacuum%nmz
     630           0 :                          ui = u_1(jz,l,m,jspin)
     631           0 :                          uei = ue_1(jz,l,m,jspin)
     632           0 :                          den%vacz(jz,ivac,jspin) = den%vacz(jz,ivac,jspin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
     633             :                       END DO
     634             :                    END DO
     635             :                 END DO
     636             :              ELSE          ! D1
     637        1122 :                 DO l = 1,nv2(jspin)
     638        1094 :                    aa = CMPLX(0.0,0.0)
     639        1094 :                    bb = CMPLX(0.0,0.0)
     640        1094 :                    ba = CMPLX(0.0,0.0)
     641        1094 :                    ab = CMPLX(0.0,0.0)
     642        7539 :                    DO n = 1,ne
     643        6445 :                       aa = aa + we(n)*CONJG(ac(l,n,jspin))*ac(l,n,jspin)
     644        6445 :                       bb = bb + we(n)*CONJG(bc(l,n,jspin))*bc(l,n,jspin)
     645        6445 :                       ab = ab + we(n)*CONJG(ac(l,n,jspin))*bc(l,n,jspin)
     646        6445 :                       ba = ba + we(n)*CONJG(bc(l,n,jspin))*ac(l,n,jspin)
     647        6445 :                       qout = REAL(CONJG(ac(l,n,jspin))*ac(l,n,jspin)+tei(l,jspin)*CONJG(bc(l,n,jspin))*bc(l,n,jspin))
     648        7539 :                       dos%qvac(n,ivac,ikpt,jspin) = dos%qvac(n,ivac,ikpt,jspin) + qout*cell%area
     649             :                    END DO
     650      274622 :                    DO  jz = 1,vacuum%nmz
     651      273500 :                       ui = u(jz,l,jspin)
     652      273500 :                       uei = ue(jz,l,jspin)
     653      274594 :                       den%vacz(jz,ivac,jspin) = den%vacz(jz,ivac,jspin) +REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
     654             :                    ENDDO
     655             :                 END DO
     656             :              END IF ! D1
     657             :           ENDIF
     658             :        END IF
     659             :        !
     660             :        !     ****************** change for vacuum density of states shz Jan.96 ***
     661             :        !
     662          28 :        IF (banddos%vacdos) THEN                       
     663             :           !
     664             :           !  ----> d_z^2-Tip needs: |d^2(psi)/dz^2 - kappa^2/3 psi|^2
     665             :           !
     666           0 :           IF (vacuum%nstm.EQ.2) THEN
     667           0 :              DO l=1,nv2(jspin)
     668           0 :                 DO n = 1,ne
     669           0 :                    aa = CONJG(ac(l,n,jspin))*ac(l,n,jspin)
     670           0 :                    bb = CONJG(bc(l,n,jspin))*bc(l,n,jspin)
     671           0 :                    ab = CONJG(ac(l,n,jspin))*bc(l,n,jspin)
     672           0 :                    ba = CONJG(bc(l,n,jspin))*ac(l,n,jspin)
     673           0 :                    aae = -vacuum%tworkf*aa*2/3
     674           0 :                    bbe = -vacuum%tworkf*bb*2/3
     675           0 :                    abe = -vacuum%tworkf*ab*2/3
     676           0 :                    bae = -vacuum%tworkf*ba*2/3
     677           0 :                    aaee = aa*vacuum%tworkf*vacuum%tworkf*4/9
     678           0 :                    bbee = bb*vacuum%tworkf*vacuum%tworkf*4/9
     679           0 :                    abee = ab*vacuum%tworkf*vacuum%tworkf*4/9
     680           0 :                    baee = ba*vacuum%tworkf*vacuum%tworkf*4/9
     681           0 :                    DO jj = 1,vacuum%layers
     682             :                       !
     683             :                       !     ----> either integrated LDOS(z1,z2) or LDOS(z1)
     684             :                       !     
     685           0 :                       IF (input%integ) THEN
     686           0 :                          ll = 1
     687           0 :                          DO ii = vacuum%izlay(jj,1),vacuum%izlay(jj,2)
     688           0 :                             ui = u(ii,l,jspin)
     689           0 :                             uei = ue(ii,l,jspin)
     690           0 :                             ddui = ddu(ii,l)
     691           0 :                             dduei = ddue(ii,l)
     692             :                             yy(ll) = REAL(aaee*ui*ui+bbee*uei*uei+(abee+baee)*ui*uei+aa*ddui*ddui+bb*&
     693             :                                  dduei*dduei+(ab+ba)*ddui*dduei+2*aae*ui*ddui+2*bbe*uei*dduei+&
     694           0 :                                  (abe+bae)*(ui*dduei+uei*ddui))*cell%area
     695           0 :                             ll = ll+1
     696             :                          END DO
     697           0 :                          CALL qsf(vacuum%delz,yy,RESULT,ll-1,0)
     698           0 :                          dos%qvlay(n,jj,ivac,ikpt,jspin) = dos%qvlay(n,jj,ivac,ikpt,jspin) + RESULT(1)
     699             :                       ELSE
     700           0 :                          ui = u(vacuum%izlay(jj,1),l,jspin)       
     701           0 :                          uei = ue(vacuum%izlay(jj,1),l,jspin)
     702           0 :                          ddui = ddu(vacuum%izlay(jj,1),l)
     703           0 :                          dduei = ddue(vacuum%izlay(jj,1),l)
     704             :                          yy(1) = REAL(aaee*ui*ui+bbee*uei*uei+&
     705             :                               (abee+baee)*ui*uei+aa*ddui*ddui+&
     706             :                               bb*dduei*dduei+(ab+ba)*ddui*dduei+&
     707             :                               2*aae*ui*ddui+2*bbe*uei*dduei+&
     708           0 :                               (abe+bae)*(ui*dduei+uei*ddui))
     709           0 :                          dos%qvlay(n,jj,ivac,ikpt,jspin) = dos%qvlay(n,jj,ivac,ikpt,jspin) +yy (1)
     710             :                       END IF
     711             :                    END DO
     712             :                 END DO
     713             :              END DO
     714             :              !     
     715             :              !     ----> s-Tip = calculate LDOS and(!) p_z-Tip (since u->du/dz, ue->due/dz)
     716             :              !     
     717             :           ELSE 
     718           0 :              IF (ABS(vacuum%locx(1)-vacuum%locx(2)).LE.eps) THEN
     719             :                 !     
     720             :                 !     ----> integrated over 2D-unit cell
     721             :                 !
     722           0 :                 IF (noco%l_noco) THEN
     723           0 :                    isp_start = 1
     724           0 :                    isp_end   = input%jspins
     725             :                 ELSE
     726           0 :                    isp_start = jspin
     727           0 :                    isp_end   = jspin
     728             :                 ENDIF
     729           0 :                 DO ispin = isp_start, isp_end
     730           0 :                    DO l=1,nv2(ispin)
     731           0 :                       DO n = 1,ne
     732           0 :                          aa = CONJG(ac(l,n,ispin))*ac(l,n,ispin)
     733           0 :                          bb = CONJG(bc(l,n,ispin))*bc(l,n,ispin)
     734           0 :                          ab = CONJG(ac(l,n,ispin))*bc(l,n,ispin)
     735           0 :                          ba = CONJG(bc(l,n,ispin))*ac(l,n,ispin) 
     736           0 :                          DO jj = 1,vacuum%layers
     737             :                             !     
     738             :                             !     ---> either integrated (z1,z2) or slice (z1)
     739             :                             !     
     740           0 :                             IF (input%integ) THEN
     741           0 :                                ll = 1
     742           0 :                                DO ii = vacuum%izlay(jj,1),vacuum%izlay(jj,2)
     743           0 :                                   ui = u(ii,l,ispin)
     744           0 :                                   uei = ue(ii,l,ispin)
     745           0 :                                   yy(ll) = REAL(aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
     746           0 :                                   ll = ll+1
     747             :                                END DO
     748           0 :                                CALL qsf(vacuum%delz,yy,RESULT,ll-1,0)
     749           0 :                                dos%qvlay(n,jj,ivac,ikpt,ispin) = dos%qvlay(n,jj,ivac,ikpt,ispin) + RESULT(1)
     750             :                             ELSE
     751           0 :                                ui = u(vacuum%izlay(jj,1),l,ispin)       
     752           0 :                                uei = ue(vacuum%izlay(jj,1),l,ispin)
     753             :                                dos%qvlay(n,jj,ivac,ikpt,ispin) = dos%qvlay(n,jj,ivac,ikpt,ispin) + REAL(&
     754           0 :                                     aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
     755             : 
     756             :                             END IF
     757             :                          END DO
     758             :                       END DO
     759             :                    END DO
     760             :                 ENDDO
     761             :              ELSE
     762             :                 !     
     763             :                 !     ----> if LDOS should be calculated over restricted area of the 2D-unit cell
     764             :                 !     lower left corner: (locx(1), locy(1))   }  in internal
     765             :                 !     upper right corner: (locx(2), locy(2))  }  coordinates
     766             :                 !     
     767           0 :                 DO l=1, nv2(jspin)
     768           0 :                    DO l1=1, nv2(jspin)
     769           0 :                       IF (kvac1(l,jspin).EQ.kvac1(l1,jspin)) THEN
     770           0 :                          factorx = CMPLX((vacuum%locx(2)-vacuum%locx(1)), 0.)
     771             :                       ELSE
     772           0 :                          k_diff=tpi_const*(kvac1(l,jspin)-kvac1(l1,jspin))
     773           0 :                          k_d1 = k_diff*vacuum%locx(1)
     774           0 :                          k_d2 = k_diff*vacuum%locx(2)
     775             :                          factorx=( CMPLX( COS(k_d2), SIN(k_d2)) -&
     776             :                               CMPLX( COS(k_d1), SIN(k_d1)) ) /&
     777           0 :                               CMPLX( 0.,k_diff )
     778             :                       END IF
     779           0 :                       IF (kvac2(l,jspin).EQ.kvac2(l1,jspin)) THEN
     780           0 :                          factory = CMPLX((vacuum%locy(2)-vacuum%locy(1)), 0.)
     781             :                       ELSE
     782           0 :                          k_diff=tpi_const*(kvac2(l,jspin)-kvac2(l1,jspin))
     783           0 :                          k_d1 = k_diff*vacuum%locy(1)
     784           0 :                          k_d2 = k_diff*vacuum%locy(2)
     785             :                          factory=( CMPLX( COS(k_d2), SIN(k_d2)) -&
     786             :                               CMPLX( COS(k_d1), SIN(k_d1)) ) /&
     787           0 :                               CMPLX( 0.,k_diff )
     788             :                       END IF
     789           0 :                       DO n=1, ne
     790           0 :                          aa = CONJG(ac(l1,n,jspin))*ac(l,n,jspin)
     791           0 :                          bb = CONJG(bc(l1,n,jspin))*bc(l,n,jspin)
     792           0 :                          ab = CONJG(ac(l1,n,jspin))*bc(l,n,jspin)
     793           0 :                          ba = CONJG(bc(l1,n,jspin))*ac(l,n,jspin)   
     794           0 :                          DO jj = 1,vacuum%layers
     795             :                             !     
     796             :                             !     ---> either integrated (z1,z2) or slice (z1)
     797             :                             !     
     798           0 :                             IF (input%integ) THEN
     799           0 :                                ll = 1
     800           0 :                                DO ii = vacuum%izlay(jj,1), vacuum%izlay(jj,2)
     801           0 :                                   ui = u(ii,l,jspin)
     802           0 :                                   uei = ue(ii,l,jspin)
     803           0 :                                   uj = u(ii,l1,jspin)
     804           0 :                                   uej = ue(ii,l1,jspin)
     805           0 :                                   yy(ll) = REAL((aa*ui*uj+bb*uei*uej+ab*uei*uj+ba*ui*uej)*factorx*factory)
     806           0 :                                   ll = ll+1
     807             :                                END DO
     808           0 :                                CALL qsf(vacuum%delz,yy,RESULT,ll-1,0)
     809           0 :                                dos%qvlay(n,jj,ivac,ikpt,jspin) = dos%qvlay(n,jj,ivac,ikpt,jspin) + RESULT(1)
     810             :                             ELSE
     811           0 :                                ui = u(vacuum%izlay(jj,1),l,jspin)       
     812           0 :                                uei = ue(vacuum%izlay(jj,1),l,jspin)
     813           0 :                                uj = u(vacuum%izlay(jj,1),l1,jspin)
     814           0 :                                uej = ue(vacuum%izlay(jj,1),l1,jspin)
     815           0 :                                dos%qvlay(n,jj,ivac,ikpt,jspin) = REAL((aa*ui*uj + bb*uei*uej+ab*uei*uj+ba*ui**uej)*factorx*factory)
     816             :                             END IF
     817             :                          END DO
     818             :                       END DO
     819             :                    END DO
     820             :                 END DO
     821             :              END IF
     822             :           END IF
     823             :        END IF
     824             : 
     825             :        !
     826             :        !     **********************************************************************
     827             :        !
     828             :        !--->    warping part of the density (g.ne.0 stars)
     829             :        !
     830             :        !   ---> d_z^2-Tip
     831             :        !
     832          28 :        IF (vacuum%nstm.EQ.2) THEN
     833           0 :           DO l = 1,nv2(jspin)
     834           0 :              DO  l1 = 1,l - 1
     835           0 :                 i1 = kvac1(l,jspin) - kvac1(l1,jspin)
     836           0 :                 i2 = kvac2(l,jspin) - kvac2(l1,jspin)
     837           0 :                 i3 = 0
     838           0 :                 IF (iabs(i1).GT.stars%mx1) CYCLE
     839           0 :                 IF (iabs(i2).GT.stars%mx2) CYCLE
     840           0 :                 ig3 = stars%ig(i1,i2,i3)
     841           0 :                 IF (ig3.EQ.0)  CYCLE
     842           0 :                 phs = stars%rgphs(i1,i2,i3)
     843           0 :                 ig3p = stars%ig(-i1,-i2,i3)
     844           0 :                 phsp = stars%rgphs(-i1,-i2,i3)
     845           0 :                 ind2 = stars%ig2(ig3)
     846           0 :                 ind2p = stars%ig2(ig3p)
     847           0 :                 aa = 0.0
     848           0 :                 bb = 0.0
     849           0 :                 ba = 0.0
     850           0 :                 ab = 0.0
     851           0 :                 DO n = 1,ne
     852           0 :                    aa = aa + we(n)*CONJG(ac(l1,n,jspin))*ac(l,n,jspin)
     853           0 :                    bb = bb + we(n)*CONJG(bc(l1,n,jspin))*bc(l,n,jspin)
     854           0 :                    ab = ab + we(n)*CONJG(ac(l1,n,jspin))*bc(l,n,jspin)
     855           0 :                    ba = ba + we(n)*CONJG(bc(l1,n,jspin))*ac(l,n,jspin)
     856             :                 END DO
     857           0 :                 aae=-aa*2/3*vacuum%tworkf
     858           0 :                 bbe=-bb*2/3*vacuum%tworkf
     859           0 :                 abe=-ab*2/3*vacuum%tworkf
     860           0 :                 bae=-ba*2/3*vacuum%tworkf
     861           0 :                 aaee=aa*4/9*vacuum%tworkf*vacuum%tworkf
     862           0 :                 bbee=bb*4/9*vacuum%tworkf*vacuum%tworkf
     863           0 :                 abee=ab*4/9*vacuum%tworkf*vacuum%tworkf
     864           0 :                 baee=ba*4/9*vacuum%tworkf*vacuum%tworkf
     865           0 :                 DO  jz = 1,vacuum%nmzxy
     866           0 :                    ui = u(jz,l,jspin)
     867           0 :                    uj = u(jz,l1,jspin)
     868           0 :                    ddui = ddu(jz,l)
     869           0 :                    dduj = ddu(jz,l1)
     870           0 :                    uei = ue(jz,l,jspin)
     871           0 :                    uej = ue(jz,l1,jspin)
     872           0 :                    dduei = ddue(jz,l)
     873           0 :                    dduej = ddue(jz,l1)
     874             :                    t1=aaee*ui*uj+bbee*uei*uej+baee*ui*uej+abee*uei*uj&
     875             :                         + aae*(ui*dduj+uj*ddui)+bbe*(uei*dduej+uej*dduei)&
     876             :                         + abe*(ui*dduej+uj*dduei)+bae*(ddui*uej+dduj*uei)&
     877             :                         + aa*ddui*dduj+bb*dduei*dduej+ba*ddui*dduej&
     878           0 :                         + ab*dduei*dduj  
     879           0 :                    den%vacxy(jz,ind2-1,ivac,jspin) = den%vacxy(jz,ind2-1,ivac,jspin) + t1*phs/stars%nstr2(ind2)
     880           0 :                    den%vacxy(jz,ind2p-1,ivac,jspin)= den%vacxy(jz,ind2p-1,ivac,jspin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
     881             :                 ENDDO
     882             :              ENDDO
     883             :           END DO
     884             :           !
     885             :           ! ---> s-Tip and p_z-Tip
     886             :           !
     887             :        ELSE
     888             :           !=============================================================
     889             :           !           continuation of vacden....
     890             :           !=============================================================
     891          28 :           IF (noco%l_noco) THEN
     892             :              !--->       diagonal elements of the density matrix, n_11 and n_22
     893           0 :              DO ispin = 1,input%jspins
     894           0 :                 IF (oneD%odi%d1) THEN
     895           0 :                    DO l = 1,nv2(ispin)
     896           0 :                       DO m = -oneD%odi%mb,oneD%odi%mb
     897           0 :                          lprimee: DO l1 = 1,l-1
     898           0 :                             mprimee: DO m1 = -oneD%odi%mb,m-1
     899           0 :                                i3 = kvac3(l,ispin) - kvac3(l1,ispin)
     900           0 :                                m3 = m-m1
     901           0 :                                IF (m3.EQ.0 .AND. i3.EQ.0) CYCLE mprimee
     902           0 :                                IF (iabs(m3).GT.oneD%odi%M) CYCLE mprimee
     903           0 :                                IF (iabs(i3).GT.stars%mx3) CYCLE lprimee
     904           0 :                                ind1 = oneD%odi%ig(i3,m3)
     905           0 :                                ind1p = oneD%odi%ig(-i3,-m3)
     906           0 :                                IF (ind1.NE.0 .OR. ind1p.NE.0) THEN
     907           0 :                                   aa = CMPLX(0.,0.)
     908           0 :                                   bb = CMPLX(0.,0.)
     909           0 :                                   ba = CMPLX(0.,0.)
     910           0 :                                   ab = CMPLX(0.,0.)
     911           0 :                                   DO n = 1,ne 
     912           0 :                                      aa=aa+we(n)*CONJG(ac_1(l1,m1,n,ispin))*ac_1(l,m,n,ispin)
     913           0 :                                      bb=bb+we(n)*CONJG(bc_1(l1,m1,n,ispin))*bc_1(l,m,n,ispin)
     914           0 :                                      ab=ab+we(n)*CONJG(ac_1(l1,m1,n,ispin))*bc_1(l,m,n,ispin)
     915           0 :                                      ba=ba+we(n)*CONJG(bc_1(l1,m1,n,ispin))*ac_1(l,m,n,ispin)
     916             :                                   END DO
     917           0 :                                   xys1: DO jz = 1,vacuum%nmzxy
     918           0 :                                      ui = u_1(jz,l,m,ispin)
     919           0 :                                      uj = u_1(jz,l1,m1,ispin)
     920           0 :                                      uei = ue_1(jz,l,m,ispin)
     921           0 :                                      uej = ue_1(jz,l1,m1,ispin)
     922           0 :                                      t1 = aa*ui*uj + bb*uei*uej + ba*ui*uej + ab*uei*uj
     923           0 :                                      IF (ind1.NE.0) THEN
     924           0 :                                         den%vacxy(jz,ind1-1,ivac,ispin) = den%vacxy(jz,ind1-1,ivac,ispin) + t1/ oneD%odi%nst2(ind1)
     925             :                                      END IF
     926           0 :                                      IF (ind1p.NE.0) THEN
     927           0 :                                         den%vacxy(jz,ind1p-1,ivac,ispin) = den%vacxy(jz,ind1p-1,ivac,ispin) + CONJG(t1)/ oneD%odi%nst2(ind1p)
     928             :                                      END IF
     929             : 
     930             :                                   END DO xys1
     931             :                                END IF   ! ind1 and ind1p =0
     932             :                             END DO mprimee
     933             :                          END DO lprimee
     934             :                       END DO  ! m
     935             :                    END DO   ! l
     936             :                 ELSE
     937           0 :                    DO l = 1,nv2(ispin)
     938           0 :                       DO  l1 = 1,l - 1
     939           0 :                          i1 = kvac1(l,ispin) - kvac1(l1,ispin)
     940           0 :                          i2 = kvac2(l,ispin) - kvac2(l1,ispin)
     941           0 :                          i3 = 0
     942           0 :                          IF (iabs(i1).GT.stars%mx1) CYCLE
     943           0 :                          IF (iabs(i2).GT.stars%mx2) CYCLE
     944           0 :                          ig3 = stars%ig(i1,i2,i3)
     945           0 :                          IF (ig3.EQ.0)  CYCLE
     946           0 :                          phs = stars%rgphs(i1,i2,i3)
     947           0 :                          ig3p = stars%ig(-i1,-i2,i3)
     948           0 :                          phsp = stars%rgphs(-i1,-i2,i3)
     949           0 :                          ind2 = stars%ig2(ig3)
     950           0 :                          ind2p = stars%ig2(ig3p)
     951           0 :                          aa = 0.0
     952           0 :                          bb = 0.0
     953           0 :                          ba = 0.0
     954           0 :                          ab = 0.0
     955           0 :                          DO n = 1,ne
     956           0 :                             aa=aa+we(n)*CONJG(ac(l1,n,ispin))*ac(l,n,ispin)
     957           0 :                             bb=bb+we(n)*CONJG(bc(l1,n,ispin))*bc(l,n,ispin)
     958           0 :                             ab=ab+we(n)*CONJG(ac(l1,n,ispin))*bc(l,n,ispin)
     959           0 :                             ba=ba+we(n)*CONJG(bc(l1,n,ispin))*ac(l,n,ispin)
     960             :                          END DO
     961           0 :                          DO jz = 1,vacuum%nmzxy
     962           0 :                             ui = u(jz,l,ispin)
     963           0 :                             uj = u(jz,l1,ispin)
     964           0 :                             uei = ue(jz,l,ispin)
     965           0 :                             uej = ue(jz,l1,ispin)
     966           0 :                             t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
     967           0 :                             den%vacxy(jz,ind2-1,ivac,ispin) = den%vacxy(jz,ind2-1,ivac,ispin) + t1*phs/stars%nstr2(ind2)
     968           0 :                             den%vacxy(jz,ind2p-1,ivac,ispin) = den%vacxy(jz,ind2p-1,ivac,ispin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
     969             :                          ENDDO
     970             :                       ENDDO
     971             :                    ENDDO
     972             :                 END IF ! oneD%odi%d1
     973             :              END DO
     974             :              !--->          off-diagonal element of the density matrix, n_21
     975           0 :              IF (oneD%odi%d1) THEN
     976           0 :                 DO l = 1,nv2(1)
     977           0 :                    DO m = -oneD%odi%mb,oneD%odi%mb
     978           0 :                       lprimea: DO l1 = 1,nv2(2)
     979           0 :                          mprimea: DO m1 = -oneD%odi%mb,oneD%odi%mb
     980           0 :                             i3 = kvac3(l,1) - kvac3(l1,2)
     981           0 :                             m3 = m-m1
     982           0 :                             IF (iabs(m3).GT.oneD%odi%M) CYCLE mprimea
     983           0 :                             IF (iabs(i3).GT.stars%mx3) CYCLE lprimea
     984           0 :                             ind1 = oneD%odi%ig(i3,m3)
     985           0 :                             IF (ind1.NE.0) THEN
     986           0 :                                IF (m3.EQ.0 .AND. i3.EQ.0) THEN
     987           0 :                                   aa = CMPLX(0.,0.)
     988           0 :                                   bb = CMPLX(0.,0.)
     989           0 :                                   ba = CMPLX(0.,0.)
     990           0 :                                   ab = CMPLX(0.,0.)
     991           0 :                                   DO n = 1,ne
     992           0 :                                      aa=aa+we(n)*CONJG(ac_1(l1,m1,n,2))* ac_1(l,m,n,1)
     993           0 :                                      bb=bb+we(n)*CONJG(bc_1(l1,m1,n,2))* bc_1(l,m,n,1)
     994           0 :                                      ab=ab+we(n)*CONJG(ac_1(l1,m1,n,2))* bc_1(l,m,n,1)
     995           0 :                                      ba=ba+we(n)*CONJG(bc_1(l1,m1,n,2))* ac_1(l,m,n,1)
     996             :                                   END DO
     997           0 :                                   xys3: DO jz = 1,vacuum%nmzxy
     998           0 :                                      ui = u_1(jz,l,m,1)
     999           0 :                                      uj = u_1(jz,l1,m1,2)
    1000           0 :                                      uei = ue_1(jz,l,m,1)
    1001           0 :                                      uej = ue_1(jz,l1,m1,2)
    1002           0 :                                      tempCmplx = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
    1003           0 :                                      den%vacz(jz,ivac,3) = den%vacz(jz,ivac,3) + REAL(tempCmplx)
    1004           0 :                                      den%vacz(jz,ivac,4) = den%vacz(jz,ivac,4) + AIMAG(tempCmplx)
    1005             :                                   END DO xys3
    1006             :                                ELSE ! the warped part (ind1 > 1)
    1007           0 :                                   aa = CMPLX(0.,0.)
    1008           0 :                                   bb = CMPLX(0.,0.)
    1009           0 :                                   ba = CMPLX(0.,0.)
    1010           0 :                                   ab = CMPLX(0.,0.)
    1011           0 :                                   DO n = 1,ne
    1012           0 :                                      aa=aa+we(n)*CONJG(ac_1(l1,m1,n,2))* ac_1(l,m,n,1)
    1013           0 :                                      bb=bb+we(n)*CONJG(bc_1(l1,m1,n,2))* bc_1(l,m,n,1)
    1014           0 :                                      ab=ab+we(n)*CONJG(ac_1(l1,m1,n,2))* bc_1(l,m,n,1)
    1015           0 :                                      ba=ba+we(n)*CONJG(bc_1(l1,m1,n,2))* ac_1(l,m,n,1)
    1016             :                                   END DO
    1017           0 :                                   xys2: DO jz = 1,vacuum%nmzxy
    1018           0 :                                      ui = u_1(jz,l,m,1)
    1019           0 :                                      uj = u_1(jz,l1,m1,2)
    1020           0 :                                      uei = ue_1(jz,l,m,1)
    1021           0 :                                      uej = ue_1(jz,l1,m1,2)
    1022           0 :                                      t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
    1023           0 :                                      den%vacxy(jz,ind1-1,ivac,3) = den%vacxy(jz,ind1-1,ivac,3) + t1/oneD%odi%nst2(ind1)
    1024             :                                   END DO xys2
    1025             :                                END IF  ! the non-warped (ind1 = 1)
    1026             :                             END IF   ! ind1.ne.0
    1027             :                          END DO mprimea
    1028             :                       END DO lprimea
    1029             :                    END DO  ! m
    1030             :                 END DO   ! l
    1031             :              ELSE ! oneD%odi%d1
    1032           0 :                 DO l = 1,nv2(1)
    1033           0 :                    DO  l1 = 1,nv2(2)
    1034           0 :                       i1 = kvac1(l,1) - kvac1(l1,2)
    1035           0 :                       i2 = kvac2(l,1) - kvac2(l1,2)
    1036           0 :                       i3 = 0
    1037             :                       !--->                treat only the warping part
    1038           0 :                       IF (iabs(i1).GT.stars%mx1) CYCLE
    1039           0 :                       IF (iabs(i2).GT.stars%mx2) CYCLE
    1040           0 :                       ig3 = stars%ig(i1,i2,i3)
    1041           0 :                       IF (ig3.EQ.0)  CYCLE
    1042           0 :                       phs = stars%rgphs(i1,i2,i3)
    1043           0 :                       ind2 = stars%ig2(ig3)
    1044           0 :                       IF ( ind2.EQ.1) THEN
    1045             :                          !--->                non-warping part (1st star G=0)
    1046           0 :                          aa = 0.0
    1047           0 :                          bb = 0.0
    1048           0 :                          ba = 0.0
    1049           0 :                          ab = 0.0
    1050           0 :                          DO ie = 1,ne
    1051           0 :                             aa=aa+we(ie)*CONJG(ac(l1,ie,2))*ac(l,ie,1)
    1052           0 :                             bb=bb+we(ie)*CONJG(bc(l1,ie,2))*bc(l,ie,1)
    1053           0 :                             ab=ab+we(ie)*CONJG(ac(l1,ie,2))*bc(l,ie,1)
    1054           0 :                             ba=ba+we(ie)*CONJG(bc(l1,ie,2))*ac(l,ie,1)
    1055             :                          END DO
    1056           0 :                          DO jz = 1,vacuum%nmz
    1057           0 :                             ui = u(jz,l,1)
    1058           0 :                             ui2 = u(jz,l1,2)
    1059           0 :                             uei = ue(jz,l,1)
    1060           0 :                             uei2 = ue(jz,l1,2)
    1061           0 :                             tempCmplx = aa*ui2*ui + bb*uei2*uei + ab*ui2*uei + ba*uei2*ui
    1062           0 :                             den%vacz(jz,ivac,3) = den%vacz(jz,ivac,3) + REAL(tempCmplx)
    1063           0 :                             den%vacz(jz,ivac,4) = den%vacz(jz,ivac,4) + AIMAG(tempCmplx)
    1064             :                          ENDDO
    1065             :                       ELSE
    1066             :                          !--->                warping part
    1067           0 :                          aa = 0.0
    1068           0 :                          bb = 0.0
    1069           0 :                          ba = 0.0
    1070           0 :                          ab = 0.0
    1071           0 :                          DO ie = 1,ne
    1072           0 :                             aa=aa + we(ie)*CONJG(ac(l1,ie,2))*ac(l,ie,1)
    1073           0 :                             bb=bb + we(ie)*CONJG(bc(l1,ie,2))*bc(l,ie,1)
    1074           0 :                             ab=ab + we(ie)*CONJG(ac(l1,ie,2))*bc(l,ie,1)
    1075           0 :                             ba=ba + we(ie)*CONJG(bc(l1,ie,2))*ac(l,ie,1)
    1076             :                          END DO
    1077           0 :                          DO  jz = 1,vacuum%nmzxy
    1078           0 :                             ui = u(jz,l,1)
    1079           0 :                             uj = u(jz,l1,2)
    1080           0 :                             uei = ue(jz,l,1)
    1081           0 :                             uej = ue(jz,l1,2)
    1082           0 :                             t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
    1083           0 :                             den%vacxy(jz,ind2-1,ivac,3) = den%vacxy(jz, ind2-1,ivac,3) + t1*phs/stars%nstr2(ind2)
    1084             :                          ENDDO
    1085             :                       ENDIF
    1086             :                    ENDDO
    1087             :                 END DO
    1088             :              END IF ! oneD%odi%d1
    1089             :           ELSE                                ! collinear part
    1090             : 
    1091          28 :              IF (oneD%odi%d1) THEN
    1092           0 :                 DO l = 1,nv2(jspin)
    1093           0 :                    DO m = -oneD%odi%mb,oneD%odi%mb
    1094           0 :                       lprime: DO l1 = 1,l-1
    1095           0 :                          mprime: DO m1 = -oneD%odi%mb,m-1
    1096           0 :                             i3 = kvac3(l,jspin) - kvac3(l1,jspin)
    1097           0 :                             m3 = m-m1
    1098           0 :                             IF (m3.EQ.0 .AND. i3.EQ.0) CYCLE mprime
    1099           0 :                             IF (iabs(m3).GT.oneD%odi%M) CYCLE mprime
    1100           0 :                             IF (iabs(i3).GT.stars%mx3) CYCLE lprime
    1101           0 :                             ind1 = oneD%odi%ig(i3,m3)
    1102           0 :                             ind1p = oneD%odi%ig(-i3,-m3)
    1103           0 :                             IF (ind1.NE.0 .OR. ind1p.NE.0) THEN
    1104           0 :                                aa = CMPLX(0.,0.)
    1105           0 :                                bb = CMPLX(0.,0.)
    1106           0 :                                ba = CMPLX(0.,0.)
    1107           0 :                                ab = CMPLX(0.,0.)
    1108           0 :                                DO n = 1,ne
    1109           0 :                                   aa=aa+we(n)*CONJG(ac_1(l1,m1,n,jspin))* ac_1(l,m,n,jspin)
    1110           0 :                                   bb=bb+we(n)*CONJG(bc_1(l1,m1,n,jspin))* bc_1(l,m,n,jspin)
    1111           0 :                                   ab=ab+we(n)*CONJG(ac_1(l1,m1,n,jspin))* bc_1(l,m,n,jspin)
    1112           0 :                                   ba=ba+we(n)*CONJG(bc_1(l1,m1,n,jspin))* ac_1(l,m,n,jspin)
    1113             :                                END DO
    1114           0 :                                xys: DO jz = 1,vacuum%nmzxy
    1115           0 :                                   ui = u_1(jz,l,m,jspin)
    1116           0 :                                   uj = u_1(jz,l1,m1,jspin)
    1117           0 :                                   uei = ue_1(jz,l,m,jspin)
    1118           0 :                                   uej = ue_1(jz,l1,m1,jspin)
    1119           0 :                                   t1 = aa*ui*uj + bb*uei*uej + ba*ui*uej + ab*uei*uj
    1120           0 :                                   IF (ind1.NE.0) THEN
    1121           0 :                                      den%vacxy(jz,ind1-1,ivac,jspin) = den%vacxy(jz,ind1-1,ivac,jspin) + t1/ oneD%odi%nst2(ind1)
    1122             :                                   END IF
    1123           0 :                                   IF (ind1p.NE.0) THEN
    1124           0 :                                      den%vacxy(jz,ind1p-1,ivac,jspin) = den%vacxy(jz,ind1p-1,ivac,jspin) + CONJG(t1)/ oneD%odi%nst2(ind1p)
    1125             :                                   END IF
    1126             : 
    1127             :                                END DO xys
    1128             :                             END IF   ! ind1 and ind1p =0
    1129             :                          END DO mprime
    1130             :                       END DO lprime
    1131             :                    END DO  ! m
    1132             :                 END DO   ! l
    1133             : 
    1134             :              ELSE         !D1
    1135             : 
    1136        1122 :                 DO l = 1,nv2(jspin)
    1137       21478 :                    DO  l1 = 1,l - 1
    1138       21450 :                       i1 = kvac1(l,jspin) - kvac1(l1,jspin)
    1139       21450 :                       i2 = kvac2(l,jspin) - kvac2(l1,jspin)
    1140       21450 :                       i3 = 0
    1141       21450 :                       IF (iabs(i1).GT.stars%mx1) CYCLE
    1142       21450 :                       IF (iabs(i2).GT.stars%mx2) CYCLE
    1143       21450 :                       ig3 = stars%ig(i1,i2,i3)
    1144       21450 :                       IF (ig3.EQ.0)  CYCLE
    1145       21450 :                       phs = stars%rgphs(i1,i2,i3)
    1146       21450 :                       ig3p = stars%ig(-i1,-i2,i3)
    1147       21450 :                       phsp = stars%rgphs(-i1,-i2,i3)
    1148       21450 :                       ind2 = stars%ig2(ig3)
    1149       21450 :                       ind2p = stars%ig2(ig3p)
    1150       21450 :                       aa = 0.0
    1151       21450 :                       bb = 0.0
    1152       21450 :                       ba = 0.0
    1153       21450 :                       ab = 0.0
    1154      146153 :                       DO n = 1,ne
    1155      124703 :                          aa=aa+we(n)*CONJG(ac(l1,n,jspin))*ac(l,n,jspin)
    1156      124703 :                          bb=bb+we(n)*CONJG(bc(l1,n,jspin))*bc(l,n,jspin)
    1157      124703 :                          ab=ab+we(n)*CONJG(ac(l1,n,jspin))*bc(l,n,jspin)
    1158      146153 :                          ba=ba+we(n)*CONJG(bc(l1,n,jspin))*ac(l,n,jspin)
    1159             :                       END DO
    1160     2167544 :                       DO  jz = 1,vacuum%nmzxy
    1161     2145000 :                          ui = u(jz,l,jspin)
    1162     2145000 :                          uj = u(jz,l1,jspin)
    1163     2145000 :                          uei = ue(jz,l,jspin)
    1164     2145000 :                          uej = ue(jz,l1,jspin)
    1165     2145000 :                          t1 = aa*ui*uj+bb*uei*uej+ba*ui*uej+ab*uei*uj
    1166     2145000 :                          den%vacxy(jz,ind2-1,ivac,jspin) = den%vacxy(jz,ind2-1, ivac,jspin) + t1*phs/stars%nstr2(ind2)
    1167     2166450 :                          den%vacxy(jz,ind2p-1,ivac,jspin) = den%vacxy(jz,ind2p-1, ivac,jspin) + CONJG(t1)*phsp/stars%nstr2(ind2p)
    1168             :                       ENDDO
    1169             :                    ENDDO
    1170             :                 END DO
    1171             :              END IF ! D1
    1172             :           ENDIF
    1173             :        END IF
    1174             :        !=============================================================
    1175             :        !
    1176             :        !       calculate 1. to nstars. starcoefficient for each k and energy eigenvalue 
    1177             :        !           to dos%qstars(ne,layer,ivac,ikpt) if starcoeff=T (the star coefficient values are written to vacdos)
    1178             :        !
    1179          56 :        IF (vacuum%starcoeff .AND. banddos%vacdos) THEN
    1180           0 :           DO  n=1,ne
    1181           0 :              DO l = 1,nv2(jspin)
    1182           0 :                 DO  l1 = 1,l - 1
    1183           0 :                    i1 = kvac1(l,jspin) - kvac1(l1,jspin)
    1184           0 :                    i2 = kvac2(l,jspin) - kvac2(l1,jspin)
    1185           0 :                    i3 = 0
    1186           0 :                    IF (iabs(i1).GT.stars%mx1) CYCLE
    1187           0 :                    IF (iabs(i2).GT.stars%mx2) CYCLE
    1188           0 :                    ig3 = stars%ig(i1,i2,i3)
    1189           0 :                    IF (ig3.EQ.0)  CYCLE
    1190           0 :                    ind2 = stars%ig2(ig3)
    1191           0 :                    ig3p = stars%ig(-i1,-i2,i3)
    1192           0 :                    ind2p = stars%ig2(ig3p)
    1193           0 :                    IF ((ind2.GE.2.AND.ind2.LE.vacuum%nstars).OR.&
    1194           0 :                         (ind2p.GE.2.AND.ind2p.LE.vacuum%nstars)) THEN
    1195           0 :                       phs = stars%rgphs(i1,i2,i3)
    1196           0 :                       phsp = stars%rgphs(-i1,-i2,i3)
    1197           0 :                       aa = CONJG(ac(l1,n,jspin))*ac(l,n,jspin)
    1198           0 :                       bb = CONJG(bc(l1,n,jspin))*bc(l,n,jspin)
    1199           0 :                       ab = CONJG(ac(l1,n,jspin))*bc(l,n,jspin)
    1200           0 :                       ba = CONJG(bc(l1,n,jspin))*ac(l,n,jspin)
    1201           0 :                       DO jj = 1,vacuum%layers
    1202           0 :                          ui = u(vacuum%izlay(jj,1),l,jspin)
    1203           0 :                          uj = u(vacuum%izlay(jj,1),l1,jspin)
    1204           0 :                          uei = ue(vacuum%izlay(jj,1),l,jspin)
    1205           0 :                          uej = ue(vacuum%izlay(jj,1),l1,jspin)
    1206           0 :                          t1 = aa*ui*uj + bb*uei*uej +ba*ui*uej + ab*uei*uj
    1207           0 :                          IF (ind2.GE.2.AND.ind2.LE.vacuum%nstars) &
    1208           0 :                               dos%qstars(ind2-1,n,jj,ivac,ikpt,jspin) = dos%qstars(ind2-1,n,jj,ivac,ikpt,jspin)+ t1*phs/stars%nstr2(ind2)
    1209           0 :                          IF (ind2p.GE.2.AND.ind2p.LE.vacuum%nstars) &
    1210           0 :                               dos%qstars(ind2p-1,n,jj,ivac,ikpt,jspin) = dos%qstars(ind2p-1,n,jj,ivac,ikpt,jspin) +CONJG(t1)*phs/stars%nstr2(ind2p)
    1211             :                       END DO
    1212             :                    END IF
    1213             :                 ENDDO
    1214             :              END DO
    1215             :           ENDDO
    1216             :        END IF
    1217             :     ENDDO
    1218          28 :     DEALLOCATE (ac,bc,dt,dte,du,ddu,due,ddue,t,te,tei,u,ue,v,yy )
    1219             : 
    1220          28 :     IF (oneD%odi%d1) THEN
    1221           0 :        DEALLOCATE (ac_1,bc_1,dt_1,dte_1,du_1,ddu_1,due_1,ddue_1)
    1222           0 :        DEALLOCATE (t_1,te_1,tei_1,u_1,ue_1)
    1223             :     END IF ! oneD%odi%d1
    1224             : 
    1225          28 :     IF(vacuum%nvac.EQ.1) THEN
    1226          28 :        den%vacz(:,2,:) = den%vacz(:,1,:)
    1227          28 :        IF (sym%invs) THEN
    1228          28 :           den%vacxy(:,:,2,:) = CONJG(den%vacxy(:,:,1,:))
    1229             :        ELSE
    1230           0 :           den%vacxy(:,:,2,:) = den%vacxy(:,:,1,:)
    1231             :        END IF
    1232             :     END IF
    1233             : 
    1234          28 :     CALL timestop("vacden")
    1235             : 
    1236          84 :   END SUBROUTINE vacden
    1237             : END MODULE m_vacden

Generated by: LCOV version 1.13