Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 : MODULE m_fertri
7 : !
8 : ! calculates fermi energy and weights using triangular (tetrahedron) method
9 : !
10 : USE m_juDFT
11 : USE m_types
12 : USE m_constants
13 : USE m_tetraef
14 : USE m_dosef
15 : USE m_dosint
16 : USE m_doswt
17 : USE m_xmlOutput
18 :
19 : IMPLICIT NONE
20 :
21 : CONTAINS
22 :
23 2 : SUBROUTINE fertri(input,noco,kpts,irank,ne,jspins,zc,eig,sfac,&
24 2 : ef,seigv,w,l_output)
25 :
26 : TYPE(t_input), INTENT(IN) :: input
27 : TYPE(t_noco), INTENT(IN) :: noco
28 : TYPE(t_kpts), INTENT(IN) :: kpts
29 : INTEGER, INTENT(IN) :: jspins,irank
30 : REAL, INTENT(IN) :: zc,sfac
31 : INTEGER, INTENT(IN) :: ne(:,:)!(nkpt,jspins)
32 : REAL, INTENT(OUT) :: seigv
33 : REAL, INTENT(INOUT) :: ef
34 : REAL, INTENT(OUT) :: w(:,:,:) !(neig,nkpt,jspins)
35 : REAL, INTENT(INOUT) :: eig(:,:,:)!(neig,nkpt,jspins)
36 : LOGICAL, INTENT(IN) :: l_output
37 :
38 : REAL :: chmom,ct,del,dez,ei,emax,emin,s,s1,workf
39 : REAL :: lb,ub,e_set,seigvTemp
40 : INTEGER :: i,ic,j,jsp,k,neig,jj
41 : INTEGER :: nemax(2)
42 : CHARACTER(LEN=20) :: attributes(2)
43 : REAL, PARAMETER :: de = 5.0e-3 !Step for initial search
44 :
45 2 : IF (irank == 0 .and. l_output) THEN
46 2 : WRITE (oUnit,FMT=8000)
47 : 8000 FORMAT (/,/,10x,'linear triangular method')
48 : END IF
49 764 : w=0.0
50 : !
51 : !---> clear w and set eig=-9999.9
52 2 : e_set = -9999.9
53 2 : IF (.NOT.input%film) e_set = 1.0e10
54 4 : DO jsp = 1,jspins
55 2 : nemax(jsp) = 0.0
56 44 : DO k = 1,kpts%nkpt
57 40 : nemax(jsp) = max0(nemax(jsp),ne(k,jsp))
58 760 : DO i = 1,ne(k,jsp)
59 760 : w(i,k,jsp) = 0.
60 : ENDDO
61 42 : DO i = ne(k,jsp)+1,size(w,1)
62 0 : w(i,k,jsp) = 0.
63 40 : eig(i,k,jsp) = e_set
64 : ENDDO
65 : ENDDO
66 : ENDDO
67 : !
68 : ! sfac = 2.0/real(jspins)
69 :
70 2 : IF(.not.input%film) THEN
71 764 : lb = MINVAL(eig) - 0.01
72 2 : ub = ef + 0.2
73 2 : CALL tetra_ef(kpts,jspins,lb,ub,eig,zc,sfac,ef,w,l_output)
74 : ELSE
75 0 : IF (irank == 0) THEN
76 0 : WRITE (oUnit,FMT=*) 'ef_hist=',ef
77 : END IF
78 :
79 0 : ei = ef
80 : !jr emin = -9999.9
81 0 : emin = +9999.9
82 0 : emax = -emin
83 : !Find initial boundaries for the bisection method
84 0 : ic = 1
85 0 : DO WHILE (emin.GT.emax)
86 0 : ic = ic + 1
87 :
88 0 : CALL dosint(ei,nemax,jspins,kpts,sfac,eig,ct)
89 :
90 0 : IF (irank == 0 .and. l_output) WRITE (oUnit,FMT=*) 'ct=',ct
91 :
92 0 : IF (ct.LT.zc) THEN ! ei < ef
93 0 : emin = ei
94 0 : ei = ei + de
95 0 : ELSEIF (ct.GT.zc) THEN ! ei > ef
96 0 : emax = ei
97 0 : ei = ei - de
98 : ENDIF
99 :
100 0 : IF (irank==0 .AND. ic.GT.100) THEN
101 0 : WRITE (oUnit,FMT=8050) ei,ef,emin,emax,ct,zc
102 : 8050 FORMAT (/,/,10x,'error fertri: initial guess of ef off by 25 mry',&
103 : ' ei,ef,emin,emax,ct,zc',/,10x,6e16.7,/,10x,&
104 : 'check number of bands')
105 0 : CALL juDFT_error("initial guess of ef off by 25 mry",calledby="fertri")
106 : END IF
107 :
108 : ENDDO
109 :
110 0 : IF (ct.NE.zc) THEN
111 0 : IF (irank == 0 .and. l_output) WRITE (oUnit,FMT=*) '2nd dosint'
112 : !---> refine ef to a value of 5 mry * (2**-20)
113 0 : iterate : DO i = 1, 40
114 0 : ei = 0.5* (emin+emax)
115 :
116 0 : CALL dosint(ei,nemax,jspins,kpts,sfac,eig,ct)
117 :
118 0 : IF (irank == 0) WRITE (oUnit,FMT=*) 'i=',i,', ct=',ct
119 0 : IF ( ct == zc ) THEN
120 : EXIT iterate
121 0 : ELSEIF ( ct > zc ) THEN
122 0 : emax = ei
123 : ELSE
124 0 : emin = ei
125 : ENDIF
126 : ENDDO iterate
127 : ENDIF
128 0 : ef = ei
129 0 : del = emax - emin
130 0 : dez = zc - ct
131 0 : workf = -hartree_to_ev_const*ef
132 0 : IF (irank == 0 .and. l_output) THEN
133 0 : WRITE (oUnit,FMT=8030) ef,workf,del,dez
134 : 8030 FORMAT(/,10x,'fermi energy=',f10.5,' har',/,10x,'work function='&
135 : ,f10.5,' ev',/,10x,'uncertainity in energy and weights=',&
136 : 2e16.6)
137 : END IF
138 : !
139 : !---> obtain dos at ef
140 : !
141 0 : CALL dosef(ei,nemax,jspins,kpts,sfac,eig,l_output)
142 : !
143 : !---> obtain weights needed for integration
144 : !
145 0 : CALL doswt(ei,nemax,jspins,kpts,eig,w)
146 :
147 : ENDIF ! .NOT.input%film
148 :
149 : !
150 : !---> write weights
151 : !
152 : ! DO jsp = 1,jspins
153 : ! neig = nemax(jsp)
154 : ! DO i = 1,neig
155 : ! DO k = 1,kpts%nkpt
156 : ! WRITE (oUnit,FMT=*) 'w(',i,',',k,',',jsp,')=',w(i,k,jsp)
157 : ! ENDDO
158 : ! ENDDO
159 : ! ENDDO
160 :
161 : !find degenerate states
162 4 : DO jsp=1,jspins
163 44 : DO k=1,kpts%nkpt
164 40 : i=1
165 546 : DO while(i<nemax(jsp))
166 : j=1
167 680 : do while (abs(eig(i,k,jsp)-eig(i+j,k,jsp))<1E-9)
168 196 : j=j+1
169 680 : if (i+j>nemax(jsp)) exit
170 : ENDDO
171 504 : if (j>1) THEN
172 176 : j=j-1
173 : !Make sure all weights are equal
174 920 : w(i:i+j,k,jsp)=sum(w(i:i+j,k,jsp))/(j+1)
175 176 : i=i+j
176 : endif
177 504 : i=i+1
178 : enddo
179 : enddo
180 : enddo
181 : !
182 : !---> obtain sum of weights and valence eigenvalues
183 : !
184 :
185 2 : s1 = 0.
186 2 : seigv = 0.
187 4 : DO jsp = 1,jspins
188 2 : s = 0.
189 2 : neig = nemax(jsp)
190 38 : DO i = 1,neig
191 758 : DO k = 1,kpts%nkpt
192 720 : s = s + w(i,k,jsp)
193 756 : seigv = seigv + w(i,k,jsp)*eig(i,k,jsp)
194 : ENDDO
195 : ENDDO
196 4 : s1 = s1 + s
197 : ENDDO
198 2 : seigv = sfac*seigv
199 2 : chmom = s1 - jspins*s
200 :
201 2 : seigvTemp = seigv
202 2 : IF (noco%l_soc .AND. (.NOT. noco%l_noco)) THEN
203 0 : seigvTemp = seigvTemp / 2.0
204 : END IF
205 :
206 2 : IF (irank == 0 .and. l_output) THEN
207 6 : attributes = ''
208 2 : WRITE(attributes(1),'(f20.10)') seigvTemp
209 2 : WRITE(attributes(2),'(a)') 'Htr'
210 6 : CALL writeXMLElement('sumValenceSingleParticleEnergies',(/'value','units'/),attributes)
211 2 : WRITE (oUnit,FMT=8040) seigvTemp,s1,chmom
212 : 8040 FORMAT (/,10x,'sum of valence eigenvalues=',f20.10,5x,&
213 : 'sum of weights=',f10.6,/,10x,'moment=',f12.6)
214 : END IF
215 :
216 2 : END SUBROUTINE fertri
217 : END MODULE m_fertri
|
