Commit 687daedf authored by Gregor Michalicek's avatar Gregor Michalicek

Merge branch 'develop' of iffgit.fz-juelich.de:fleur/fleur into develop

parents e498c60c 03299042
......@@ -47,8 +47,6 @@ pages:
paths:
- public
only:
- schedules
- triggers
- web
doxygen:
......@@ -65,8 +63,6 @@ doxygen:
- build/fleur_MPI
- build/inpgen
only:
- schedules
- triggers
- web
......@@ -179,4 +175,4 @@ gfortran-coverage:
url: https://fleur.iffgit.fz-juelich.de/fleur/coverage_html
only:
- web
- schedules
\ No newline at end of file
......@@ -3,16 +3,41 @@ Welcome to the source code of FLEUR
Please note that the documentation of the
code can be found at the [FLEUR Homepage]
(http://www.flapw.de/pm/index.php?n=User-Documentation.Documentation).
(http://www.flapw.de/).
For further instructions on Installation/Usage,
please check the [FLEUR Homepage]
(http://www.flapw.de/pm/index.php?n=User-Documentation.Documentation).
(http://www.flapw.de/).
The rest of this document summarizes only the
structure of the FLEUR source code. Did
we mention to check the Homepage?
## Bugs in FLEUR
You might experience bugs in FLEUR :-).
If you find a bug you should:
A) Report this bug by generating an Issue. Please describe in
detail the relevant input and what happens. You should consider using
the bug-template for your issue as this will help you providing us with
the relevant information.
or/and
B) Provide a bugfix. If the bug is only present in the development branch/ is due
to a new feature under development simply commit your fix to the development branch.
If you are fixing a bug in a release-version, please:
* check out the git release branch: ```git checkout --track origin/release```
* create a bugfix branch: ```git checkout -b bugfix_SOME_NAME_HERE```
* apply your changes, test them and commit them
* push your bugfix branch to the server: ``` git push -u origin bugfix_SOME_NAME_HERE```
* create a merge request on the gitlab to have you bugfix merged with the release branch
* check out the develop branch: ```git checkout develop```
* merge your fix into the develop branch: ```git merge bugfix_SOME_NAME_HERE```
## Structure of the source code
The source of FLEUR is organized in several
subdirectories. Some of them collect code
specific for particular features, others code
......
......@@ -101,7 +101,7 @@ CONTAINS
qtot, qistot
call tmp_potden%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
call init_mt_grid(input%jspins, atoms, sphhar, xcpot, sym)
call init_mt_grid(input%jspins, atoms, sphhar, xcpot%needs_grad(), sym)
do n_atm =1,atoms%ntype
call mt_from_grid(atoms, sphhar, n_atm, input%jspins, mt(:,:,n_atm), &
tmp_potden%mt(:,0:,n_atm,:))
......@@ -112,8 +112,8 @@ CONTAINS
enddo
call finish_mt_grid()
call init_pw_grid(xcpot, stars, sym, cell)
call pw_from_grid(xcpot, stars, .False., is, tmp_potden%pw)
call init_pw_grid(xcpot%needs_grad(), stars, sym, cell)
call pw_from_grid(xcpot%needs_grad(), stars, .False., is, tmp_potden%pw)
call finish_pw_grid()
call integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, tmp_potden, &
......
......@@ -205,13 +205,13 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
IF (.NOT.((jspin.EQ.2).AND.noco%l_noco)) THEN
! valence density in the interstitial region
CALL pwden(stars,kpts,banddos,oneD,input,mpi,noco,cell,atoms,sym,ikpt,&
jspin,lapw,noccbd,we,eig,den,results,force%f_b8,zMat,dos)
jspin,lapw,noccbd,ev_list,we,eig,den,results,force%f_b8,zMat,dos)
! charge of each valence state in this k-point of the SBZ in the layer interstitial region of the film
IF (l_dosNdir.AND.PRESENT(slab)) CALL q_int_sl(jspin,ikpt,stars,atoms,sym,cell,noccbd,lapw,slab,oneD,zMat)
IF (l_dosNdir.AND.PRESENT(slab)) CALL q_int_sl(jspin,ikpt,stars,atoms,sym,cell,noccbd,ev_list,lapw,slab,oneD,zMat)
! valence density in the vacuum region
IF (input%film) THEN
CALL vacden(vacuum,dimension,stars,oneD, kpts,input,sym,cell,atoms,noco,banddos,&
gVacMap,we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,lapw,enpara%evac,eig,den,zMat,dos)
gVacMap,we,ikpt,jspin,vTot%vacz(:,:,jspin),noccbd,ev_list,lapw,enpara%evac,eig,den,zMat,dos)
END IF
END IF
IF (input%film) CALL regCharges%sumBandsVac(vacuum,dos,noccbd,ikpt,jsp_start,jsp_end,eig,we)
......@@ -227,18 +227,18 @@ SUBROUTINE cdnval(eig_id, mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,st
! perform Brillouin zone integration and summation over the
! bands in order to determine the energy parameters for each atom and angular momentum
CALL eparas(ispin,atoms,noccbd,mpi,ikpt,noccbd,we,eig,&
CALL eparas(ispin,atoms,noccbd,ev_list,mpi,ikpt,noccbd,we,eig,&
skip_t,cdnvalJob%l_evp,eigVecCoeffs,usdus,regCharges,dos,banddos%l_mcd,mcd)
IF (noco%l_mperp.AND.(ispin==jsp_end)) CALL qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos)
IF (noco%l_mperp.AND.(ispin==jsp_end)) CALL qal_21(dimension,atoms,input,noccbd,ev_list,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos)
! layer charge of each valence state in this k-point of the SBZ from the mt-sphere region of the film
IF (l_dosNdir) THEN
IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
IF (PRESENT(slab)) CALL q_mt_sl(ispin,atoms,noccbd,ev_list,ikpt,noccbd,skip_t,noccbd,eigVecCoeffs,usdus,slab)
IF (banddos%l_orb.AND.ANY((/banddos%alpha,banddos%beta,banddos%gamma/).NE.0.0)) THEN
CALL abcrot2(atoms,banddos,noccbd,eigVecCoeffs,ispin) ! rotate ab-coeffs
IF (PRESENT(orbcomp)) CALL orb_comp(ispin,ikpt,noccbd,atoms,noccbd,usdus,eigVecCoeffs,orbcomp)
IF (PRESENT(orbcomp)) CALL orb_comp(ispin,ikpt,noccbd,ev_list,atoms,noccbd,usdus,eigVecCoeffs,orbcomp)
END IF
ENDIF
CALL calcDenCoeffs(atoms,sphhar,sym,we,noccbd,eigVecCoeffs,ispin,denCoeffs)
......
......@@ -23,7 +23,7 @@ MODULE m_eparas
!***********************************************************************
!
CONTAINS
SUBROUTINE eparas(jsp,atoms,noccbd, mpi,ikpt,ne,we,eig,skip_t,l_evp,eigVecCoeffs,&
SUBROUTINE eparas(jsp,atoms,noccbd,ev_list,mpi,ikpt,ne,we,eig,skip_t,l_evp,eigVecCoeffs,&
usdus,regCharges,dos,l_mcd,mcd)
USE m_types
IMPLICIT NONE
......@@ -39,6 +39,7 @@ CONTAINS
INTEGER, INTENT (IN) :: noccbd,jsp
INTEGER, INTENT (IN) :: ne,ikpt ,skip_t
LOGICAL, INTENT (IN) :: l_mcd,l_evp
INTEGER, INTENT (IN) :: ev_list(noccbd)
! ..
! .. Array Arguments ..
REAL, INTENT (IN) :: eig(:)!(dimension%neigd),
......@@ -98,7 +99,7 @@ CONTAINS
DO icore = 1, mcd%ncore(n)
DO ipol = 1, 3
index = 3*(n-1) + ipol
mcd%mcd(index,icore,i,ikpt,jsp)=mcd%mcd(index,icore,i,ikpt,jsp) + fac*(&
mcd%mcd(index,icore,ev_list(i),ikpt,jsp)=mcd%mcd(index,icore,ev_list(i),ikpt,jsp) + fac*(&
suma * CONJG(mcd%m_mcd(icore,lm+1,index,1))*mcd%m_mcd(icore,lm+1,index,1) +&
sumb * CONJG(mcd%m_mcd(icore,lm+1,index,2))*mcd%m_mcd(icore,lm+1,index,2) +&
sumab* CONJG(mcd%m_mcd(icore,lm+1,index,2))*mcd%m_mcd(icore,lm+1,index,1) +&
......@@ -107,7 +108,7 @@ CONTAINS
ENDDO
ENDIF ! end MCD
ENDDO
dos%qal(l,n,i,ikpt,jsp) = (suma+sumb*usdus%ddn(l,n,jsp))/atoms%neq(n)
dos%qal(l,n,ev_list(i),ikpt,jsp) = (suma+sumb*usdus%ddn(l,n,jsp))/atoms%neq(n)
ENDDO
nt1 = nt1 + atoms%neq(n)
ENDDO
......@@ -120,8 +121,8 @@ CONTAINS
DO l = 0,3
DO n = 1,atoms%ntype
DO i = (skip_t+1),noccbd
regCharges%ener(l,n,jsp) = regCharges%ener(l,n,jsp) + dos%qal(l,n,i,ikpt,jsp)*we(i)*eig(i)
regCharges%sqal(l,n,jsp) = regCharges%sqal(l,n,jsp) + dos%qal(l,n,i,ikpt,jsp)*we(i)
regCharges%ener(l,n,jsp) = regCharges%ener(l,n,jsp) + dos%qal(l,n,ev_list(i),ikpt,jsp)*we(i)*eig(i)
regCharges%sqal(l,n,jsp) = regCharges%sqal(l,n,jsp) + dos%qal(l,n,ev_list(i),ikpt,jsp)*we(i)
ENDDO
ENDDO
ENDDO
......@@ -174,7 +175,7 @@ CONTAINS
! llo > 3 used for unoccupied states only
IF( l .GT. 3 ) CYCLE
DO i = 1,ne
dos%qal(l,ntyp,i,ikpt,jsp)= dos%qal(l,ntyp,i,ikpt,jsp) + ( 1.0/atoms%neq(ntyp) )* (&
dos%qal(l,ntyp,ev_list(i),ikpt,jsp)= dos%qal(l,ntyp,ev_list(i),ikpt,jsp) + ( 1.0/atoms%neq(ntyp) )* (&
qaclo(i,lo,ntyp)*usdus%uulon(lo,ntyp,jsp)+qbclo(i,lo,ntyp)*usdus%dulon(lo,ntyp,jsp) )
END DO
DO lop = 1,atoms%nlo(ntyp)
......@@ -182,7 +183,7 @@ CONTAINS
DO i = 1,ne
regCharges%enerlo(lo,ntyp,jsp) = regCharges%enerlo(lo,ntyp,jsp) +qlo(i,lop,lo,ntyp)*we(i)*eig(i)
regCharges%sqlo(lo,ntyp,jsp) = regCharges%sqlo(lo,ntyp,jsp) + qlo(i,lop,lo,ntyp)*we(i)
dos%qal(l,ntyp,i,ikpt,jsp)= dos%qal(l,ntyp,i,ikpt,jsp) + ( 1.0/atoms%neq(ntyp) ) *&
dos%qal(l,ntyp,ev_list(i),ikpt,jsp)= dos%qal(l,ntyp,ev_list(i),ikpt,jsp) + ( 1.0/atoms%neq(ntyp) ) *&
qlo(i,lop,lo,ntyp)*usdus%uloulopn(lop,lo,ntyp,jsp)
ENDDO
ENDIF
......
......@@ -36,7 +36,7 @@ CONTAINS
! .. Local Scalars ..
INTEGER iri
REAL b_xavh,scale,b_con_outx,b_con_outy,mx,my,mz,&
& alphh,betah,mz_tmp,mx_mix,my_mix,mz_mix
& alphh,betah,mz_tmp,mx_mix,my_mix,mz_mix,absmag
REAL rho11,rho22, alphdiff
COMPLEX rho21
! ..
......@@ -56,11 +56,12 @@ CONTAINS
mx = 2*REAL(qa21(itype))
my = 2*AIMAG(qa21(itype))
mz = chmom(itype,1) - chmom(itype,2)
WRITE (6,8025) mx,my
absmag=SQRT(mx*mx+my*my+mz*mz)
WRITE (6,8025) mx,my,mz,absmag
!---> determine the polar angles of the moment vector in the local frame
CALL pol_angle(mx,my,mz,betah,alphh)
WRITE (6,8026) betah,alphh
8025 FORMAT(2x,'--> local frame: ','mx=',f9.5,' my=',f9.5)
8025 FORMAT(2x,'--> local frame: ','mx=',f9.5,' my=',f9.5,' mz=',f9.5,' |m|=',f9.5)
8026 FORMAT(2x,'-->',10x,' delta beta=',f9.5,&
& ' delta alpha=',f9.5)
......
......@@ -7,7 +7,7 @@
MODULE m_pwden
CONTAINS
SUBROUTINE pwden(stars,kpts,banddos,oneD, input,mpi,noco,cell,atoms,sym, &
ikpt,jspin,lapw,ne,we,eig,den,results,f_b8,zMat,dos)
ikpt,jspin,lapw,ne,ev_list,we,eig,den,results,f_b8,zMat,dos)
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! In this subroutine the star function expansion coefficients of
! the plane wave charge density is determined.
......@@ -98,6 +98,7 @@ CONTAINS
REAL,INTENT(IN) :: we(:) !(nobd)
REAL,INTENT(IN) :: eig(:)!(dimension%neigd)
INTEGER, INTENT(IN) :: ev_list(ne)
!-----> BASIS FUNCTION INFORMATION
INTEGER,INTENT(IN):: ne
!-----> CHARGE DENSITY INFORMATION
......@@ -408,7 +409,7 @@ CONTAINS
length_zfft(1) = stars%kq1_fft
length_zfft(2) = stars%kq2_fft
length_zfft(3) = stars%kq3_fft
call fft_interface(3,length_zfft,zfft,forw)
call fft_interface(3,length_zfft,zfft,forw,iv1d(1:lapw%nv(jspin),jspin))
psir = real(zfft)
psii = aimag(zfft)
!--------------------------------
......@@ -442,7 +443,7 @@ CONTAINS
length_zfft(1) = stars%kq1_fft
length_zfft(2) = stars%kq2_fft
length_zfft(3) = stars%kq3_fft
call fft_interface(3,length_zfft,zfft,forw)
call fft_interface(3,length_zfft,zfft,forw,iv1d(1:lapw%nv(jspin),jspin))
kpsir = real(zfft)
kpsii = aimag(zfft)
!--------------------------------
......@@ -494,7 +495,7 @@ CONTAINS
ENDDO
DO istr = 1,stars%ng3_fft
CALL pwint(stars,atoms,sym, oneD,cell,istr,x)
dos%qis(nu,ikpt,1) = dos%qis(nu,ikpt,1) + REAL(cwk(istr)*x)/cell%omtil/REAL(ifftq3)
dos%qis(ev_list(nu),ikpt,1) = dos%qis(ev_list(nu),ikpt,1) + REAL(cwk(istr)*x)/cell%omtil/REAL(ifftq3)
ENDDO
cwk=0.0
......@@ -504,7 +505,7 @@ CONTAINS
ENDDO
DO istr = 1,stars%ng3_fft
CALL pwint(stars,atoms,sym, oneD,cell, istr, x)
dos%qis(nu,ikpt,input%jspins) = dos%qis(nu,ikpt,input%jspins) + REAL(cwk(istr)*x)/cell%omtil/REAL(ifftq3)
dos%qis(ev_list(nu),ikpt,input%jspins) = dos%qis(ev_list(nu),ikpt,input%jspins) + REAL(cwk(istr)*x)/cell%omtil/REAL(ifftq3)
ENDDO
ENDIF
ELSE
......
MODULE m_qintsl
USE m_juDFT
CONTAINS
SUBROUTINE q_int_sl(isp,ikpt,stars,atoms,sym,cell,ne,lapw,slab,oneD,zMat)
SUBROUTINE q_int_sl(isp,ikpt,stars,atoms,sym,cell,ne,ev_list,lapw,slab,oneD,zMat)
! *******************************************************
! calculate the charge of the En(k) state
! in the interstitial region of each leyer
......@@ -24,6 +24,9 @@ CONTAINS
!
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: ne,isp,ikpt
INTEGER, INTENT (IN) :: ev_list(ne)
! ..
! .. Local Scalars ..
REAL q1,zsl1,zsl2,qi,volsli,volintsli
......@@ -105,7 +108,7 @@ CONTAINS
DO j = 1,stars%ng3
qi = qi + z_z(j)*stfunint(j,i)
ENDDO
slab%qintsl(i,n,ikpt,isp) = qi
slab%qintsl(i,ev_list(n),ikpt,isp) = qi
ENDDO ! over vacuum%layers
ENDDO ! over states
......
......@@ -8,7 +8,7 @@ CONTAINS
!
!***********************************************************************
!
SUBROUTINE q_mt_sl(jsp,atoms,nobd,ikpt,ne,skip_t,noccbd,eigVecCoeffs,usdus,slab)
SUBROUTINE q_mt_sl(jsp,atoms,nobd,ev_list,ikpt,ne,skip_t,noccbd,eigVecCoeffs,usdus,slab)
USE m_types_setup
USE m_types_usdus
USE m_types_cdnval, ONLY: t_eigVecCoeffs, t_slab
......@@ -21,6 +21,9 @@ CONTAINS
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: nobd,jsp
INTEGER, INTENT (IN) :: ne,ikpt ,skip_t,noccbd
INTEGER, INTENT (IN) :: ev_list(nobd)
! ..
! .. Local Scalars ..
INTEGER i,l,lo ,natom,nn,ntyp,nt1,nt2,m
......@@ -136,7 +139,7 @@ CONTAINS
DO ntyp = 1,atoms%ntype
qq = qq + qmttot(ntyp,i)*slab%nmtsl(ntyp,nl)
ENDDO
slab%qmtsl(nl,i,ikpt,jsp) = qq
slab%qmtsl(nl,ev_list(i),ikpt,jsp) = qq
ENDDO
ENDDO
! DO ntyp = 1,ntype
......
......@@ -5,7 +5,7 @@ MODULE m_qal21
!***********************************************************************
!
CONTAINS
SUBROUTINE qal_21(dimension,atoms,input,noccbd,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos)
SUBROUTINE qal_21(dimension,atoms,input,noccbd,ev_list,noco,eigVecCoeffs,denCoeffsOffdiag,ikpt,dos)
USE m_types_setup
USE m_types_dos
......@@ -25,6 +25,8 @@ CONTAINS
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: noccbd,ikpt
INTEGER, INTENT (IN) :: ev_list(noccbd)
! .. Local Scalars ..
INTEGER i,l,lo,lop ,natom,nn,ntyp
INTEGER nt1,nt2,lm,n,ll1,ipol,icore,index,m
......@@ -152,10 +154,10 @@ CONTAINS
state : DO i = 1, noccbd
lls : DO l = 0,3
CALL rot_den_mat(noco%alph(n),noco%beta(n),&
dos%qal(l,n,i,ikpt,1),dos%qal(l,n,i,ikpt,2),qal21(l,n,i))
dos%qal(l,n,ev_list(i),ikpt,1),dos%qal(l,n,ev_list(i),ikpt,2),qal21(l,n,i))
IF (.FALSE.) THEN
IF (n==1) WRITE(*,'(3i3,4f10.5)') l,n,i,qal21(l,n,i),dos%qal(l,n,i,ikpt,:)
q_loc(1,1) = dos%qal(l,n,i,ikpt,1); q_loc(2,2) = dos%qal(l,n,i,ikpt,2)
IF (n==1) WRITE(*,'(3i3,4f10.5)') l,n,i,qal21(l,n,i),dos%qal(l,n,ev_list(i),ikpt,:)
q_loc(1,1) = dos%qal(l,n,ev_list(i),ikpt,1); q_loc(2,2) = dos%qal(l,n,ev_list(i),ikpt,2)
q_loc(1,2) = qal21(l,n,i); q_loc(2,1) = CONJG(q_loc(1,2))
q_hlp = MATMUL( TRANSPOSE( CONJG(chi) ) ,q_loc)
q_loc = MATMUL(q_hlp,chi)
......
......@@ -6,7 +6,7 @@ MODULE m_vacden
! *************************************************************
CONTAINS
SUBROUTINE vacden(vacuum,DIMENSION,stars,oneD,kpts,input,sym,cell,atoms,noco,banddos,&
gVacMap,we,ikpt,jspin,vz,ne,lapw,evac,eig,den,zMat,dos)
gVacMap,we,ikpt,jspin,vz,ne,ev_list,lapw,evac,eig,den,zMat,dos)
!***********************************************************************
! ****** change vacden(....,q) for vacuum density of states shz Jan.96
......@@ -70,6 +70,7 @@ CONTAINS
INTEGER,PARAMETER :: n2max=13
REAL,PARAMETER :: emax=2.0/hartree_to_ev_const
! .. Array Arguments ..
INTEGER, INTENT(IN) :: ev_list(ne)
REAL, INTENT(IN) :: evac(2,input%jspins)
REAL, INTENT(IN) :: we(DIMENSION%neigd)
REAL :: vz(vacuum%nmzd,2) ! Note this breaks the INTENT(IN) from cdnval. It may be read from a file in this subroutine.
......@@ -530,7 +531,7 @@ CONTAINS
ab = ab + we(n)*CONJG(ac(l,n,jspin))*bc(l,n,jspin)
ba = ba + we(n)*CONJG(bc(l,n,jspin))*ac(l,n,jspin)
qout = REAL(CONJG(ac(l,n,jspin))*ac(l,n,jspin)+tei(l,jspin)*CONJG(bc(l,n,jspin))*bc(l,n,jspin))
dos%qvac(n,ivac,ikpt,jspin) = dos%qvac(n,ivac,ikpt,jspin) + qout*cell%area
dos%qvac(ev_list(n),ivac,ikpt,jspin) = dos%qvac(ev_list(n),ivac,ikpt,jspin) + qout*cell%area
END DO
aae=-aa*vacuum%tworkf*2/3
bbe=-bb*vacuum%tworkf*2/3
......@@ -578,7 +579,7 @@ CONTAINS
qout = REAL(CONJG(ac_1(l,m,n,ispin))*ac_1(l,m,n,ispin) +&
tei_1(l,m,ispin)*CONJG(bc_1(l,m,n,ispin))*&
bc_1(l,m,n,ispin))
dos%qvac(n,ivac,ikpt,ispin) = dos%qvac(n,ivac,ikpt,ispin)+qout*cell%area
dos%qvac(ev_list(n),ivac,ikpt,ispin) = dos%qvac(ev_list(n),ivac,ikpt,ispin)+qout*cell%area
END DO
DO jz = 1,vacuum%nmz
ui = u_1(jz,l,m,ispin)
......@@ -599,7 +600,7 @@ CONTAINS
ab=ab + we(n)*CONJG(ac(l,n,ispin))*bc(l,n,ispin)
ba=ba + we(n)*CONJG(bc(l,n,ispin))*ac(l,n,ispin)
qout = REAL(CONJG(ac(l,n,ispin))*ac(l,n,ispin)+tei(l,ispin)*CONJG(bc(l,n,ispin))*bc(l,n,ispin))
dos%qvac(n,ivac,ikpt,ispin) = dos%qvac(n,ivac,ikpt,ispin) + qout*cell%area
dos%qvac(ev_list(n),ivac,ikpt,ispin) = dos%qvac(ev_list(n),ivac,ikpt,ispin) + qout*cell%area
END DO
DO jz = 1,vacuum%nmz
ui = u(jz,l,ispin)
......@@ -624,7 +625,7 @@ CONTAINS
ba = ba + we(n)*CONJG(bc_1(l,m,n,jspin))*ac_1(l,m,n,jspin)
qout = REAL(CONJG(ac_1(l,m,n,jspin))*ac_1(l,m,n,jspin)+&
tei_1(l,m,jspin)*CONJG(bc_1(l,m,n,jspin))*bc_1(l,m,n,jspin))
dos%qvac(n,ivac,ikpt,jspin) = dos%qvac(n,ivac,ikpt,jspin)+qout*cell%area
dos%qvac(ev_list(n),ivac,ikpt,jspin) = dos%qvac(ev_list(n),ivac,ikpt,jspin)+qout*cell%area
END DO
DO jz = 1,vacuum%nmz
ui = u_1(jz,l,m,jspin)
......@@ -645,7 +646,7 @@ CONTAINS
ab = ab + we(n)*CONJG(ac(l,n,jspin))*bc(l,n,jspin)
ba = ba + we(n)*CONJG(bc(l,n,jspin))*ac(l,n,jspin)
qout = REAL(CONJG(ac(l,n,jspin))*ac(l,n,jspin)+tei(l,jspin)*CONJG(bc(l,n,jspin))*bc(l,n,jspin))
dos%qvac(n,ivac,ikpt,jspin) = dos%qvac(n,ivac,ikpt,jspin) + qout*cell%area
dos%qvac(ev_list(n),ivac,ikpt,jspin) = dos%qvac(ev_list(n),ivac,ikpt,jspin) + qout*cell%area
END DO
DO jz = 1,vacuum%nmz
ui = u(jz,l,jspin)
......@@ -695,7 +696,7 @@ CONTAINS
ll = ll+1
END DO
CALL qsf(vacuum%delz,yy,RESULT,ll-1,0)
dos%qvlay(n,jj,ivac,ikpt,jspin) = dos%qvlay(n,jj,ivac,ikpt,jspin) + RESULT(1)
dos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) = dos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) + RESULT(1)
ELSE
ui = u(vacuum%izlay(jj,1),l,jspin)
uei = ue(vacuum%izlay(jj,1),l,jspin)
......@@ -706,7 +707,7 @@ CONTAINS
bb*dduei*dduei+(ab+ba)*ddui*dduei+&
2*aae*ui*ddui+2*bbe*uei*dduei+&
(abe+bae)*(ui*dduei+uei*ddui))
dos%qvlay(n,jj,ivac,ikpt,jspin) = dos%qvlay(n,jj,ivac,ikpt,jspin) +yy (1)
dos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) = dos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) +yy (1)
END IF
END DO
END DO
......@@ -746,11 +747,11 @@ CONTAINS
ll = ll+1
END DO
CALL qsf(vacuum%delz,yy,RESULT,ll-1,0)
dos%qvlay(n,jj,ivac,ikpt,ispin) = dos%qvlay(n,jj,ivac,ikpt,ispin) + RESULT(1)
dos%qvlay(ev_list(n),jj,ivac,ikpt,ispin) = dos%qvlay(ev_list(n),jj,ivac,ikpt,ispin) + RESULT(1)
ELSE
ui = u(vacuum%izlay(jj,1),l,ispin)
uei = ue(vacuum%izlay(jj,1),l,ispin)
dos%qvlay(n,jj,ivac,ikpt,ispin) = dos%qvlay(n,jj,ivac,ikpt,ispin) + REAL(&
dos%qvlay(ev_list(n),jj,ivac,ikpt,ispin) = dos%qvlay(ev_list(n),jj,ivac,ikpt,ispin) + REAL(&
aa*ui*ui+bb*uei*uei+(ab+ba)*ui*uei)
END IF
......@@ -806,13 +807,13 @@ CONTAINS
ll = ll+1
END DO
CALL qsf(vacuum%delz,yy,RESULT,ll-1,0)
dos%qvlay(n,jj,ivac,ikpt,jspin) = dos%qvlay(n,jj,ivac,ikpt,jspin) + RESULT(1)
dos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) = dos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) + RESULT(1)
ELSE
ui = u(vacuum%izlay(jj,1),l,jspin)
uei = ue(vacuum%izlay(jj,1),l,jspin)
uj = u(vacuum%izlay(jj,1),l1,jspin)
uej = ue(vacuum%izlay(jj,1),l1,jspin)
dos%qvlay(n,jj,ivac,ikpt,jspin) = REAL((aa*ui*uj + bb*uei*uej+ab*uei*uj+ba*ui**uej)*factorx*factory)
dos%qvlay(ev_list(n),jj,ivac,ikpt,jspin) = REAL((aa*ui*uj + bb*uei*uej+ab*uei*uj+ba*ui**uej)*factorx*factory)
END IF
END DO
END DO
......@@ -1205,9 +1206,9 @@ CONTAINS
uej = ue(vacuum%izlay(jj,1),l1,jspin)
t1 = aa*ui*uj + bb*uei*uej +ba*ui*uej + ab*uei*uj
IF (ind2.GE.2.AND.ind2.LE.vacuum%nstars) &