LCOV - code coverage report
Current view: top level - dos - evaldos.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 194 253 76.7 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             :       MODULE m_evaldos
       2             :       CONTAINS
       3           7 :       SUBROUTINE evaldos(eig_id,input,banddos,vacuum,kpts,atoms,sym,noco,oneD,cell,results,dos,&
       4             :                          dimension,efermiarg,bandgap,l_mcd,mcd,slab,orbcomp)
       5             : !----------------------------------------------------------------------
       6             : !
       7             : !     vk: k-vectors
       8             : !     wk weight of k-point (not used)
       9             : !     nevk no of eigenvalues
      10             : !     ev eigenvalue
      11             : !     qal partial charges
      12             : !             partial charge interstitial: qal(lmax*ntype+1...
      13             : !             partial charge vacuum      : qal(lmax*ntype+2...
      14             : !     qlay,qstar read in vacuum charges
      15             : !     qval partial charges in the vacuum
      16             : !             qval(m*nstars,neigd,nkptd):charge in m'th layer
      17             : !             qval(m*nstars+nstars,... ):star-resolved charge of that layer
      18             : !             qval(layers*nstars+....  ):same for second vacuum
      19             : !     ntb=max(nevk)
      20             : !
      21             : !----------------------------------------------------------------------
      22             :       USE m_triang
      23             :       USE m_maketetra
      24             :       USE m_tetrados
      25             :       USE m_dosbin
      26             :       USE m_ptdos
      27             :       USE m_smooth
      28             :       USE m_types
      29             :       USE m_constants
      30             :       USE m_cdn_io
      31             :       IMPLICIT NONE
      32             :       INTEGER,INTENT(IN)             :: eig_id
      33             :       TYPE(t_dimension),INTENT(IN)   :: dimension
      34             :       TYPE(t_oneD),INTENT(IN)        :: oneD
      35             :       TYPE(t_banddos),INTENT(IN)     :: banddos
      36             :       TYPE(t_input),INTENT(IN)       :: input
      37             :       TYPE(t_vacuum),INTENT(IN)      :: vacuum
      38             :       TYPE(t_noco),INTENT(IN)        :: noco
      39             :       TYPE(t_sym),INTENT(IN)         :: sym
      40             :       TYPE(t_cell),INTENT(IN)        :: cell
      41             :       TYPE(t_results),INTENT(IN)     :: results
      42             :       TYPE(t_dos),INTENT(IN)         :: dos
      43             :       TYPE(t_mcd),INTENT(IN)         :: mcd
      44             :       TYPE(t_slab),INTENT(IN)        :: slab
      45             :       TYPE(t_orbcomp),INTENT(IN)     :: orbcomp
      46             :       TYPE(t_kpts),INTENT(IN)        :: kpts
      47             :       TYPE(t_atoms),INTENT(IN)       :: atoms
      48             : 
      49             :       REAL,    INTENT(IN) :: efermiarg, bandgap
      50             :       LOGICAL, INTENT(IN) :: l_mcd 
      51             :  
      52             : !    locals
      53             :       INTEGER, PARAMETER ::  lmax= 4, ned = 1301
      54             :       INTEGER  i,s,v,index,jspin,k,l,l1,l2,ln,n,nl,ntb,ntria,ntetra
      55             :       INTEGER  icore,qdim,n_orb,ncored
      56             :       REAL     as,de,efermi,emax,emin,qmt,sigma,totdos,efermiPrev
      57             :       REAL     e_up,e_lo,e_test1,e_test2,fac,sumwei,dk,eFermiCorrection
      58             :       LOGICAL  l_tria,l_orbcomp,l_error
      59             : 
      60          14 :       INTEGER  itria(3,2*kpts%nkpt),itetra(4,6*kpts%nkpt)
      61          28 :       REAL     voltet(6*kpts%nkpt),kx(kpts%nkpt),vkr(3,kpts%nkpt)
      62          21 :       REAL     ev(dimension%neigd,kpts%nkpt),e(ned),gpart(ned,atoms%ntype),atr(2*kpts%nkpt)
      63          14 :       REAL     e_grid(ned+1),spect(ned,3*atoms%ntype),ferwe(dimension%neigd,kpts%nkpt)
      64          14 :       REAL,    ALLOCATABLE :: qal(:,:,:),qval(:,:,:),qlay(:,:,:),g(:,:)
      65           7 :       REAL,    ALLOCATABLE :: mcd_local(:,:,:)
      66             :       REAL,    ALLOCATABLE :: qvac(:,:)
      67             :       CHARACTER(len=2) :: spin12(2),ch_mcd(3)
      68             :       CHARACTER(len=8) :: chntype*2,chform*19
      69             :       DATA spin12/'.1' , '.2'/
      70             :       DATA ch_mcd/'.+' , '.-' , '.0'/
      71             : 
      72           7 :       ncored =  MAX(0,MAXVAL(mcd%ncore))
      73           7 :       qdim = lmax*atoms%ntype+3
      74           7 :       l_orbcomp = banddos%l_orb
      75           7 :       IF (banddos%ndir.EQ.-3) THEN
      76           0 :         qdim = 2*slab%nsld 
      77           0 :         n_orb = 0
      78           0 :         IF (banddos%l_orb) THEN
      79           0 :            n_orb = banddos%orbCompAtom
      80           0 :            WRITE (*,*) 'DOS: orbcomp',n_orb
      81           0 :            qdim = 23
      82             :         END IF
      83             :       ENDIF
      84             :       ALLOCATE( qal(qdim,dimension%neigd,kpts%nkpt),&
      85             :      &          qval(vacuum%nstars*vacuum%layers*vacuum%nvac,dimension%neigd,kpts%nkpt),&
      86           7 :      &          qlay(dimension%neigd,vacuum%layerd,2))
      87           7 :       IF (l_mcd) THEN 
      88           1 :          ALLOCATE(mcd_local(3*atoms%ntype*ncored,dimension%neigd,kpts%nkpt) )
      89             :       ELSE
      90           6 :          ALLOCATE(mcd_local(0,0,0))
      91             :       ENDIF
      92             : !
      93             : ! scale energies
      94           7 :       sigma = banddos%sig_dos*hartree_to_ev_const
      95           7 :       emin =min(banddos%e1_dos*hartree_to_ev_const,banddos%e2_dos*hartree_to_ev_const)
      96           7 :       emax =max(banddos%e1_dos*hartree_to_ev_const,banddos%e2_dos*hartree_to_ev_const)
      97           7 :       efermi = efermiarg*hartree_to_ev_const
      98             :  
      99           7 :       WRITE (6,'(a)') 'DOS-Output is generated!'
     100             : 
     101           7 :       IF ( NINT((emax - emin)/sigma) > ned ) THEN
     102           0 :         WRITE(6,*) 'sig_dos too small for DOS smoothing:'   
     103           0 :         WRITE(6,*) 'Reduce energy window or enlarge banddos%sig_dos!'
     104           0 :         WRITE(6,*) 'For now: setting sigma to zero !'
     105           0 :         sigma = 0.0
     106             :       ENDIF
     107             : 
     108           7 :       WRITE (6,*) 'sigma=   ' , sigma
     109           7 :       WRITE (6,*) 'emax=   ' , emax
     110           7 :       WRITE (6,*) 'emin=   ' , emin
     111           7 :       WRITE (6,*) 'ef_inp=   ' , efermi
     112             : !
     113             : !     create energy grid
     114           7 :       emax = emax - efermi
     115           7 :       emin = emin - efermi
     116           7 :       de = (emax-emin)/(ned-1)
     117        9114 :       DO i=1,ned
     118        9114 :          e(i) = emin + (i-1)*de
     119             :       ENDDO
     120             :  
     121           7 :       IF ( l_mcd ) THEN ! create an energy grid for mcd-spectra
     122           1 :         e_lo =  9.9*10.0**9 
     123           1 :         e_up = -9.9*10.0**9     
     124           3 :         DO jspin = 1,input%jspins
     125           7 :           DO n = 1,atoms%ntype
     126          46 :             DO icore = 1 , mcd%ncore(n)
     127          40 :               e_lo = min(mcd%e_mcd(n,jspin,icore),e_lo)
     128          44 :               e_up = max(mcd%e_mcd(n,jspin,icore),e_up)
     129             :             ENDDO
     130             :           ENDDO
     131             :         ENDDO
     132           1 :         e_lo = e_lo*hartree_to_ev_const - efermi - emax 
     133           1 :         e_up = e_up*hartree_to_ev_const - efermi
     134           1 :         de = (e_up-e_lo)/(ned-1)
     135        1302 :         DO i=1,ned
     136        1301 :           e_grid(i) = e_lo + (i-1)*de
     137        1302 :           spect(i,:) = 0.0
     138             :         ENDDO
     139           1 :         e_grid(ned+1) = e_lo + ned*de
     140             :       ENDIF
     141             : 
     142          15 :       DO jspin = 1,input%jspins
     143           8 :          ntb = 0
     144         148 :          DO k = 1,kpts%nkpt
     145             : 
     146         140 :             qal(:,:,k) = 0.0
     147         140 :             qval(:,:,k) = 0.0
     148             : 
     149         140 :             ntb = max(ntb,results%neig(k,jspin))
     150         140 :             IF (l_mcd) mcd_local(:,:,k) = RESHAPE(mcd%mcd(:,1:ncored,:,k,jspin),(/3*atoms%ntype*ncored,dimension%neigd/))
     151         140 :             IF (.NOT.l_orbcomp) THEN
     152         140 :                qal(1:lmax*atoms%ntype,:,k)=reshape(dos%qal(0:,:,:,k,jspin),(/lmax*atoms%ntype,size(dos%qal,3)/))
     153         140 :                qal(lmax*atoms%ntype+2,:,k)=dos%qvac(:,1,k,jspin) ! vacuum 1
     154         140 :                qal(lmax*atoms%ntype+3,:,k)=dos%qvac(:,2,k,jspin) ! vacuum 2
     155         140 :                qal(lmax*atoms%ntype+1,:,k)=dos%qis(:,k,jspin)    ! interstitial
     156             :             ELSE
     157           0 :                IF (n_orb == 0) THEN
     158           0 :                   qal(1:slab%nsld,:,k)             = slab%qintsl(:,:,k,jspin)
     159           0 :                   qal(slab%nsld+1:2*slab%nsld,:,k) = slab%qmtsl(:,:,k,jspin)
     160             :                ELSE
     161           0 :                   DO i = 1, 23
     162           0 :                      DO l = 1, results%neig(k,jspin)
     163           0 :                         qal(i,l,k) = orbcomp%comp(l,i,n_orb,k,jspin)*orbcomp%qmtp(l,n_orb,k,jspin)/10000.
     164             :                      END DO
     165           0 :                      DO l = results%neig(k,jspin)+1, dimension%neigd
     166           0 :                         qal(i,l,k) = 0.0
     167             :                      END DO
     168             :                   END DO
     169             :                END IF
     170             :             END IF
     171             : !
     172             : !     set vacuum partial charge zero, if bulk calculation
     173             : !     otherwise, write vacuum charge in correct arrays
     174             : !
     175         140 :             IF ((.NOT.input%film).AND.(banddos%ndir.NE.-3)) THEN
     176        8600 :                DO n = 1,dimension%neigd
     177        4160 :                   qal(lmax*atoms%ntype+2,n,k) = 0.0
     178        4300 :                   qal(lmax*atoms%ntype+3,n,k) = 0.0
     179             :                ENDDO
     180           0 :             ELSEIF ( banddos%vacdos .and. input%film ) THEN
     181           0 :                DO i = 1,results%neig(k,jspin)
     182           0 :                   DO v = 1,vacuum%nvac
     183           0 :                      DO l = 1,vacuum%layers
     184           0 :                         index = (l-1)*vacuum%nstars + (v-1)*(vacuum%nstars*vacuum%layers) + 1
     185           0 :                         qval(index,i,k) = qlay(i,l,v)
     186           0 :                         DO s = 1,vacuum%nstars - 1
     187           0 :                            qval(index+s,i,k) = real(dos%qstars(s,i,l,v,k,jspin))
     188             :                         ENDDO
     189             :                      ENDDO
     190             :                   ENDDO
     191             :                ENDDO
     192             :             ENDIF
     193             : !
     194             : !     calculate interstitial dos if not noco
     195             : !     in the noco case, qis has been calculated in pwden and is read in from tmp_dos
     196             : !
     197         140 :             IF ((.NOT.noco%l_noco).AND.(banddos%ndir.NE.-3)) THEN
     198        4300 :                DO i = 1 , dimension%neigd
     199        4160 :                   qal(lmax*atoms%ntype+1,i,k) = 1.
     200       12480 :                   DO nl = 1 , atoms%ntype
     201        8320 :                      l1 = lmax*(nl-1) + 1
     202        8320 :                      l2 = lmax*nl
     203        8320 :                      qmt=0.0
     204       41600 :                      DO l = l1 , l2
     205       41600 :                         qmt = qmt + qal(l,i,k)*atoms%neq(nl)
     206             :                      ENDDO
     207       12480 :                      qal(lmax*atoms%ntype+1,i,k) = qal(lmax*atoms%ntype+1,i,k) - qmt
     208             :                   ENDDO
     209             :                   qal(lmax*atoms%ntype+1,i,k) = qal(lmax*atoms%ntype+1,i,k)&
     210        4300 :                        -qal(lmax*atoms%ntype+2,i,k)*(3-vacuum%nvac) -qal(lmax*atoms%ntype+3,i,k)*(vacuum%nvac-1) 
     211             :                ENDDO
     212             :             ENDIF
     213             : !
     214             : !---- >     convert eigenvalues to ev and shift them by efermi
     215             : !
     216        4300 :             DO i = 1 , results%neig(k,jspin)
     217        4300 :                ev(i,k) = results%eig(i,k,jspin)*hartree_to_ev_const - efermi
     218             :             ENDDO
     219         148 :             DO i = results%neig(k,jspin) + 1, dimension%neigd
     220         140 :                ev(i,k) = 9.9e+99
     221             :             ENDDO
     222             : !
     223             : !
     224             :          ENDDO                                                 ! end of k-point loop
     225             : !
     226             : !     calculate the triangles!
     227             : !
     228           8 :          IF ( jspin.EQ.1 ) THEN
     229           7 :            l_tria=.true.
     230           7 :            IF (input%film .AND. .NOT.oneD%odi%d1) THEN
     231           0 :              CALL triang(kpts%bk,kpts%nkpt,itria,ntria,atr,as,l_tria)
     232           0 :              IF (sym%invs) THEN
     233           0 :                IF (abs(sym%nop2*as-0.5).GT.0.000001) l_tria=.false.
     234             :              ELSE
     235           0 :                IF (abs(sym%nop2*as-1.0).GT.0.000001) l_tria=.false.
     236             :              ENDIF
     237           0 :              write(*,*) as,sym%nop2,l_tria
     238             : !             l_tria=.true.
     239             :            ELSE
     240           7 :              IF (input%l_inpXML) THEN
     241           3 :                 IF (input%tria) THEN
     242           1 :                    ntetra = kpts%ntet
     243          45 :                    DO i = 1, ntetra
     244          44 :                       itetra(1:4,i) = kpts%ntetra(1:4,i)
     245          45 :                       voltet(i) = kpts%voltet(i) / ntetra
     246             :                    END DO
     247           1 :                    l_tria = input%tria
     248           1 :                    GOTO 67
     249             :                 ELSE
     250             :                    GOTO 66
     251             :                 END IF
     252             :              END IF
     253           4 :              OPEN (41,file='kpts',FORM='formatted',STATUS='old')
     254          88 :              DO i = 1, kpts%nkpt+1
     255          84 :                 READ (41,*,END=66,ERR=66)
     256             :              ENDDO
     257           6 :              READ (41,'(i5)',END=66,ERR=66) ntetra
     258          92 :              READ (41,'(4(4i6,4x))') ((itetra(i,k),i=1,4),k=1,ntetra)
     259           2 :              READ (41,'(4f20.13)') (voltet(k),k=1,ntetra)
     260           2 :              CLOSE(41)
     261           2 :              voltet(1:ntetra) = voltet(1:ntetra) / ntetra
     262             :              l_tria=.true.
     263           2 :              GOTO 67
     264             :  66          CONTINUE                       ! no tetrahedron-information of file
     265           4 :              CALL triang(kpts%bk,kpts%nkpt,itria,ntria,atr,as,l_tria)
     266           4 :              l_tria=.true.
     267             : ! YM: tetrahedrons is not the way in 1D
     268           4 :              IF (oneD%odi%d1) as = 0.0         
     269           4 :              IF (sym%invs) THEN
     270           3 :                IF (abs(sym%nop2*as-1.0).GT.0.000001) l_tria=.false.
     271             :              ELSE
     272           1 :                IF (abs(sym%nop2*as-0.5).GT.0.000001) l_tria=.false.
     273             :              ENDIF
     274             : 
     275           4 :              IF (l_tria) THEN
     276             :                CALL make_tetra(kpts%nkpt,kpts%bk,ntria,itria,atr,&
     277           0 :                     ntetra,itetra,voltet)
     278             :              ELSE
     279           4 :                WRITE (6,*) 'no tetrahedron method with these k-points!'
     280           4 :                WRITE (6,*) sym%nop2,as
     281             :              ENDIF
     282             :  67          CONTINUE                       ! tetrahedron-information read or created
     283             :            ENDIF
     284             :          ENDIF
     285             : !
     286           8 :         IF ( .not.l_mcd ) THEN
     287           6 :          ALLOCATE (g(ned,qdim))
     288             :         ELSE
     289           2 :          ALLOCATE (g(ned,3*atoms%ntype*ncored))
     290             :         ENDIF
     291             : !
     292           8 :          IF ( l_tria.and.(.not.l_mcd).and.(banddos%ndir.NE.-3) ) THEN
     293             : !
     294             : !     DOS calculation: use triangular method!!
     295             : !
     296           3 :             IF ( input%film ) THEN
     297             : !             CALL ptdos(
     298             : !    >                  emin,emax,jspins,ned,qdim,neigd,
     299             : !    >                  ntria,as,atr,2*nkpt,itria,nkpt,ev,qal,e,
     300             : !    <                  g)
     301             :               CALL ptdos(emin,emax,input%jspins,ned,qdim,ntb,ntria,as,&
     302             :                         atr,2*kpts%nkpt,itria,kpts%nkpt,ev(1:ntb,1:kpts%nkpt),&
     303           0 :                         qal(:,1:ntb,1:kpts%nkpt),e, g)
     304             :             ELSE
     305           3 :               write(*,*) efermi
     306             :               CALL tetra_dos(lmax,atoms%ntype,dimension%neigd,ned,ntetra,kpts%nkpt,&
     307           3 :                             itetra,efermi,voltet,e,results%neig(:,jspin), ev,qal, g)
     308           3 :               IF (input%jspins.EQ.1) g(:,:) = 2 * g(:,:)
     309             :             ENDIF
     310             :          ELSE
     311             : !
     312             : !     DOS calculation: use histogram method
     313             : !
     314           5 :             IF ( .not.l_mcd ) THEN
     315             :             CALL dos_bin(input%jspins,qdim,ned,emin,emax,dimension%neigd,kpts%nkpt,&
     316           3 :                  results%neig(:,jspin),kpts%wtkpt(1:kpts%nkpt),ev,qal, g)
     317             :             ELSE
     318             :             CALL dos_bin(input%jspins,3*atoms%ntype*ncored,ned,emin,emax,ntb,kpts%nkpt,&
     319           2 :                  results%neig(:,jspin),kpts%wtkpt(1:kpts%nkpt),ev(1:ntb,1:kpts%nkpt), mcd_local(1:3*atoms%ntype*ncored,1:ntb,1:kpts%nkpt), g)
     320             :             ENDIF
     321             :          ENDIF
     322             : !
     323             : !---- >     smoothening
     324             : !
     325           8 :          IF ( .not.l_mcd ) THEN
     326           6 :             IF ( sigma.GT.0.0 ) THEN
     327          72 :               DO ln = 1 , qdim
     328          72 :                 CALL smooth(e,g(1,ln),sigma,ned)
     329             :               ENDDO
     330             :             ENDIF
     331             :  
     332             : !*** sum up for all atoms
     333             :  
     334           6 :          IF (banddos%ndir.NE.-3) THEN
     335          18 :             DO l = 1 , atoms%ntype
     336          12 :                l1 = lmax*(l-1) + 1
     337          12 :                l2 = lmax*l
     338       15630 :                DO i = 1 , ned
     339       15612 :                   gpart(i,l) = 0.0
     340       78072 :                   DO nl = l1 , l2
     341       78060 :                      gpart(i,l) = gpart(i,l) + g(i,nl)
     342             :                   ENDDO
     343             :                ENDDO
     344             :             ENDDO
     345           0 :          ELSEIF (n_orb == 0) THEN
     346           0 :             DO l = 1, slab%nsld
     347           0 :                nl = slab%nsld+l
     348           0 :                DO i = 1 , ned
     349           0 :                   gpart(i,l) = g(i,l) + g(i,nl)
     350             :                ENDDO
     351             :             ENDDO
     352             :          ENDIF
     353             :     
     354             : !**** write out DOS
     355           6 :          OPEN (18,FILE='DOS'//spin12(jspin))
     356             : 
     357        7812 :          DO i = 1 , ned
     358        7806 :            totdos = 0.0
     359        7812 :            IF (banddos%ndir.NE.-3) THEN
     360       23418 :              DO nl = 1 , atoms%ntype
     361       23418 :                 totdos = totdos + gpart(i,nl)*atoms%neq(nl)
     362             :              ENDDO
     363             :              totdos = totdos + g(i,lmax*atoms%ntype+1) + g(i,lmax*atoms%ntype+2) *&
     364        7806 :                   (3 - vacuum%nvac) + g(i,lmax*atoms%ntype+3)*(vacuum%nvac - 1)
     365        7806 :              IF (atoms%ntype < 20) THEN
     366        7806 :                 WRITE (18,99001)  e(i),totdos,g(i,lmax*atoms%ntype+1), &
     367        7806 :                      g(i,lmax*atoms%ntype+2),g(i,lmax*atoms%ntype+3),&
     368       15612 :                      (gpart(i,l),l=1,atoms%ntype), (g(i,l),l=1,atoms%ntype*lmax)
     369             :              ELSE
     370           0 :              WRITE (18,99001)  e(i),totdos,g(i,lmax*atoms%ntype+1), &
     371           0 :                   g(i,lmax*atoms%ntype+2),g(i,lmax*atoms%ntype+3), (gpart(i,l),l=1,atoms%ntype)
     372             :           ENDIF
     373           0 :        ELSEIF (n_orb == 0) THEN
     374           0 :           DO nl = 1, slab%nsld
     375           0 :              totdos = totdos + gpart(i,nl)
     376             :           ENDDO
     377           0 :           WRITE (18,99001)  e(i),totdos,(gpart(i,nl),nl=1,slab%nsld), (g(i,l),l=1,2*slab%nsld)
     378             :            ELSE
     379           0 :              DO nl = 1 , 23
     380           0 :                 totdos = totdos + g(i,nl)
     381             :              ENDDO
     382           0 :              WRITE (18,99001)  e(i),totdos,(g(i,l),l=1,23)
     383             :            ENDIF
     384             :          ENDDO
     385           6 :          CLOSE (18)
     386             : 
     387             :          ELSE
     388           2 :            write(*,'(4f15.8)') ((mcd%e_mcd(n,jspin,i),n=1,atoms%ntype),i=1,ncored)
     389           2 :            write(*,*)
     390           2 :            write(*,'(4f15.8)') (g(800,n),n=1,3*atoms%ntype*ncored)
     391           2 :            write(*,*)
     392           2 :            write(*,'(4f15.8)') (mcd_local(n,10,8),n=1,3*atoms%ntype*ncored)
     393           6 :            DO n = 1,atoms%ntype
     394       10414 :              DO l = 1 , ned
     395       57248 :                DO icore = 1 , mcd%ncore(n)
     396   135361244 :                  DO i = 1 , ned-1
     397    67704040 :                    IF (e(i).GT.0) THEN     ! take unoccupied part only
     398     9939640 :                    e_test1 = -e(i) - efermi +mcd%e_mcd(n,jspin,icore)*hartree_to_ev_const
     399     9939640 :                    e_test2 = -e(i+1)-efermi +mcd%e_mcd(n,jspin,icore)*hartree_to_ev_const
     400     9939640 :                    IF ((e_test2.LE.e_grid(l)).AND. (e_test1.GT.e_grid(l))) THEN
     401       26836 :                      fac = (e_grid(l)-e_test1)/(e_test2-e_test1)
     402      107344 :                      DO k = 3*(n-1)+1,3*(n-1)+3
     403             :                        spect(l,k) = spect(l,k)+ g(i,3*atoms%ntype*(icore-1)+k)&
     404      107344 :                            *(1.-fac) + fac * g(i+1,3*atoms%ntype*(icore-1)+k)
     405             :                      ENDDO
     406             :                    ENDIF
     407             :                    ENDIF
     408             :                  ENDDO
     409             :                ENDDO
     410             :              ENDDO
     411             :            ENDDO
     412           2 :            CLOSE (18) 
     413             :          ENDIF
     414           8 :          DEALLOCATE (g)
     415             : !         
     416             : !------------------------------------------------------------------------------
     417             : !     now calculate the VACOS
     418             : !------------------------------------------------------------------------------
     419             :             
     420           8 :          IF ( banddos%vacdos .and. input%film ) THEN
     421           0 :             ALLOCATE(g(ned,vacuum%nstars*vacuum%layers*vacuum%nvac))
     422             : !            CALL ptdos(
     423             : !     >                 emin,emax,jspins,ned,nstars*nvac*layers,neigd,
     424             : !     >                 ntria,as,atr,2*nkpt,itria,nkptd,ev,qval,e,
     425             : !     <                 g)
     426             :             CALL ptdos(emin,emax,input%jspins,ned,vacuum%nstars*vacuum%nvac*vacuum%layers,ntb,ntria&
     427           0 :                 ,as,atr,2*kpts%nkpt,itria,kpts%nkpt,ev(1:ntb,1:kpts%nkpt), qval(:,1:ntb,1:kpts%nkpt),e,g)
     428             :             
     429             : !---- >     smoothening
     430           0 :             IF ( sigma.GT.0.0 ) THEN
     431           0 :                DO ln = 1 , vacuum%nstars*vacuum%nvac*vacuum%layers
     432           0 :                   CALL smooth(e,g(1,ln),sigma,ned)
     433             :                ENDDO
     434             :             ENDIF
     435             :             
     436             : !     write VACDOS
     437             :             
     438           0 :             OPEN (18,FILE='VACDOS'//spin12(jspin))
     439             : !            WRITE (18,'(i2,25(2x,i3))') Layers , (Zlay(l),l=1,Layers)
     440           0 :             DO i = 1 , ned
     441           0 :              WRITE (18,99001) e(i) , (g(i,l),l=1,vacuum%Layers*vacuum%Nstars*vacuum%Nvac)
     442             :             ENDDO
     443           0 :             CLOSE (18)
     444           0 :             DEALLOCATE(g)
     445             :          ENDIF
     446             : !
     447             : !------------------------------------------------------------------------------
     448             : !     for bandstructures
     449             : !------------------------------------------------------------------------------
     450             : 
     451          15 :          IF (banddos%ndir == -4) THEN
     452           3 :             eFermiCorrection = 0.0
     453           3 :             IF(bandgap.LT.(8.0*input%tkb*hartree_to_ev_const)) THEN
     454           3 :                CALL readPrevEFermi(eFermiPrev,l_error)
     455           3 :                IF(.NOT.l_error) THEN
     456           3 :                   WRITE(*,*) 'Fermi energy is automatically corrected in bands.* files.'
     457           3 :                   WRITE(*,*) 'It is consistent with last calculated density!'
     458           3 :                   WRITE(*,*) 'No manual correction (e.g. in band.gnu file) required.'
     459           3 :                   eFermiCorrection = (eFermiPrev-efermiarg)*hartree_to_ev_const
     460             :                ELSE
     461           0 :                   WRITE(*,*) 'Fermi energy in bands.* files may not be consistent with last density.'
     462           0 :                   WRITE(*,*) 'Please correct it manually (e.g. in band.gnu file).'
     463             :                END IF
     464             :             END IF
     465             : 
     466           3 :             OPEN (18,FILE='bands'//spin12(jspin))
     467           3 :             ntb = minval(results%neig(:,jspin))    
     468           3 :             kx(1) = 0.0
     469          12 :             vkr(:,1)=matmul(kpts%bk(:,1),cell%bmat)
     470          60 :             DO k = 2, kpts%nkpt
     471             :               
     472          57 :                vkr(:,k)=matmul(kpts%bk(:,k),cell%bmat)
     473             :                dk = (vkr(1,k)-vkr(1,k-1))**2 + (vkr(2,k)-vkr(2,k-1) )**2 + &
     474          57 :                     (vkr(3,k)-vkr(3,k-1))**2
     475          60 :                kx(k) = kx(k-1) + sqrt(dk)
     476             :             ENDDO
     477         111 :             DO i = 1, ntb
     478        2217 :                DO k = 1, kpts%nkpt
     479        1134 :                   write(18,'(2f15.9)') kx(k),ev(i,k)-eFermiCorrection
     480             :                ENDDO
     481             :             ENDDO
     482           3 :             CLOSE (18)
     483             :          ENDIF
     484             : 
     485             :       ENDDO
     486             : !         
     487             : !------------------------------------------------------------------------------
     488             : !     for MCD calculations ...
     489             : !------------------------------------------------------------------------------
     490             : 
     491           7 :       IF (l_mcd) THEN
     492           1 :         WRITE (chntype,'(i2)') atoms%ntype+1
     493           1 :         chform = '('//chntype//'f15.8)'
     494           1 :         IF ( sigma.GT.0.0 ) THEN
     495             :            IF ( l_mcd ) THEN
     496           7 :              DO ln = 1 , 3*atoms%ntype
     497           1 :                CALL smooth(e_grid,spect(1,ln),sigma,ned)
     498             :              ENDDO
     499             :            ENDIF
     500             :         ENDIF
     501           7 :         DO l = 1,3
     502           3 :           OPEN (18,FILE='MCD_SPEC'//ch_mcd(l))
     503        3906 :           DO i = 1 , ned
     504        3906 :           WRITE (18,FMT=chform) e_grid(i),(spect(i,3*(n-1)+l),n=1,atoms%ntype)
     505             :           ENDDO
     506           4 :           CLOSE (18)
     507             :         ENDDO
     508             :       ENDIF
     509             : 
     510           7 :       DEALLOCATE(qal,qval,qlay)
     511           7 :       IF (l_mcd) DEALLOCATE( mcd_local )
     512             : 99001 FORMAT (f10.5,110(1x,e10.3))
     513             : 
     514          14 :       END SUBROUTINE evaldos
     515             :       END MODULE m_evaldos 

Generated by: LCOV version 1.13