Line data Source code
1 : MODULE m_fermie
2 : USE m_juDFT
3 : #ifdef CPP_MPI
4 : use mpi
5 : #endif
6 : !-----------------------------------------------------------------------
7 : ! determines the fermi energy by
8 : ! gaussian-integration method c.l.fu
9 : ! triangular method (or tetrahedrons)
10 : ! or fermi-function p.kurz
11 : !----------------------------------------------------------------------
12 : CONTAINS
13 686 : SUBROUTINE fermie(eig_id,fmpi,kpts,input,noco,e_min,cell,results,l_output_param)
14 :
15 : !---------------------------------------------------f--------------------
16 : !
17 : ! a fist (T=0) approximation to the fermi-energy is determined
18 : ! by:
19 : ! zelec = sum { spindg * we }
20 : ! e=<ef
21 : !
22 : ! TREE STRUCTURE: fermie----sort
23 : ! ----fergwt
24 : ! ----fertri---triang
25 : ! |--dosint
26 : ! |--dosef -- trisrt
27 : ! +--doswt
28 : ! ----ferhis---ef_newton
29 : !
30 : !-----------------------------------------------------------------------
31 :
32 : USE m_types
33 : USE m_constants
34 : USE m_eig66_io, ONLY : read_eig,write_eig
35 : USE m_sort
36 : USE m_fertri
37 : USE m_ferhis
38 : USE m_fergwt
39 : USE m_fertetra
40 : USE m_xmlOutput
41 :
42 : IMPLICIT NONE
43 :
44 : TYPE(t_results), INTENT(INOUT) :: results
45 : TYPE(t_mpi), INTENT(IN) :: fmpi
46 : TYPE(t_input), INTENT(IN) :: input
47 : TYPE(t_noco), INTENT(IN) :: noco
48 : TYPE(t_cell), INTENT(IN) :: cell
49 : TYPE(t_kpts), INTENT(IN) :: kpts
50 : ! ..
51 : ! .. Scalar Arguments ..
52 : INTEGER, INTENT (IN) :: eig_id
53 : REAL,INTENT(IN) :: e_min
54 : ! .. Logical Arguments ..
55 : LOGICAL,INTENT(IN),OPTIONAL :: l_output_param
56 : ! ..
57 : ! .. Array Arguments ..
58 : !REAL, INTENT (OUT):: w(:,:,:) !(input%neig,kpts%nkpt,dimension%jspd)
59 : ! ..
60 : ! .. Local Scalars ..
61 : REAL del ,spindg,ssc ,ws,zc,weight,efermi,seigv
62 : INTEGER i,idummy,j,jsp,k,l,n,nbands,nstef,nv,nmat,nspins
63 : INTEGER n_help,m_spins,mspin,sslice(2)
64 : LOGICAL :: l_output
65 : ! ..
66 : ! .. Local Arrays ..
67 : !
68 8918 : INTEGER :: idxeig(SIZE(results%w_iks)),idxjsp(SIZE(results%w_iks)),idxkpt(SIZE(results%w_iks)),INDEX(SIZE(results%w_iks))
69 5488 : REAL :: e(SIZE(results%w_iks)),we(SIZE(results%w_iks))
70 : CHARACTER(LEN=20) :: attributes(5)
71 :
72 : !--- J constants
73 : !--- J constants
74 :
75 : #ifdef CPP_MPI
76 : INTEGER, PARAMETER :: comm = MPI_COMM_SELF
77 : INTEGER*4 :: nv_mpi(2)
78 : INTEGER ierr
79 : #endif
80 :
81 : ! ..
82 : ! ..
83 : !***********************************************************************
84 : ! ABBREVIATIONS
85 : !
86 : ! eig : array of eigenvalues
87 : ! wtkpt : list of the weights of each k-point (from inp-file)
88 : ! e : linear list of the eigenvalues
89 : ! we : list of weights of the eigenvalues in e
90 : ! zelec : number of electrons
91 : ! spindg : spindegeneracy (2 in nonmagnetic calculations)
92 : ! seigv : weighted sum of the occupied valence eigenvalues
93 : ! ts : entropy contribution to the free energy
94 : !
95 : !***********************************************************************
96 : ! .. Data statements ..
97 : DATA del/1.0e-6/
98 :
99 : ! initiliaze e
100 394570 : e = 0
101 :
102 : ! Logical that controls the output
103 686 : IF (.NOT.present(l_output_param)) THEN
104 686 : l_output = .true.
105 : ELSE
106 0 : l_output = l_output_param
107 : ENDIF
108 :
109 686 : IF ( fmpi%irank == 0 .and. l_output) WRITE (oUnit,FMT=8000)
110 : 8000 FORMAT (/,/,1x,'fermi energy and band-weighting factors:')
111 : !
112 : !---> READ IN EIGENVALUES
113 : !
114 686 : spindg = 2.0/REAL(input%jspins)
115 686 : n = 0
116 686 : ssc = 0.0
117 686 : n_help = 0
118 : !
119 : !---> pk non-collinear
120 686 : IF (noco%l_noco) THEN
121 192 : nspins = 1
122 : ELSE
123 494 : nspins = input%jspins
124 : ENDIF
125 : !---> pk non-collinear
126 : !
127 686 : IF (fmpi%irank == 0. .and. .NOT.judft_was_argument("-minimalOutput") .and. l_output) CALL openXMLElementNoAttributes('eigenvalues')
128 :
129 1576 : DO jsp = 1,nspins
130 11388 : DO k = 1,kpts%nkpt
131 9812 : IF (fmpi%irank == 0) THEN
132 4906 : if(input%eig66(1))CALL read_eig(eig_id,k,jsp,neig=results%neig(k,jsp),eig=results%eig(:,k,jsp))
133 4906 : if (l_output) then
134 4906 : WRITE (oUnit,'(a2,3f10.5,a12,f12.6)') 'at',kpts%bk(:,k), " k-weight:", kpts%wtkpt(k)
135 4906 : WRITE (oUnit,'(i5,a14)') results%neig(k,jsp),' eigenvalues :'
136 185836 : WRITE (oUnit,'(8f12.6)') (results%eig(i,k,jsp),i=1,results%neig(k,jsp))
137 : end if
138 4906 : IF(.NOT.judft_was_argument("-minimalOutput") .and. l_output) THEN
139 29436 : attributes = ''
140 4906 : WRITE(attributes(1),'(i0)') jsp
141 4906 : WRITE(attributes(2),'(i0)') k
142 4906 : WRITE(attributes(3),'(f15.8)') kpts%bk(1,k)
143 4906 : WRITE(attributes(4),'(f15.8)') kpts%bk(2,k)
144 4906 : WRITE(attributes(5),'(f15.8)') kpts%bk(3,k)
145 29436 : CALL writeXMLElementPoly('eigenvaluesAt',(/'spin','ikpt','k_x ','k_y ','k_z '/),attributes,results%eig(1:results%neig(k,jsp),k,jsp))
146 : END IF
147 : END IF
148 : #ifdef CPP_MPI
149 10702 : CALL MPI_BARRIER(fmpi%mpi_comm,ierr)
150 : #endif
151 : END DO
152 : ENDDO
153 : !finished reading of eigenvalues
154 686 : IF (fmpi%irank == 0 .and. .NOT.judft_was_argument("-minimalOutput") .and. l_output) CALL closeXMLElement('eigenvalues')
155 686 : IF (fmpi%irank == 0) THEN
156 :
157 343 : IF (ABS(input%fixed_moment)<1E-6) THEN
158 : !this is a standard calculation
159 : m_spins=1
160 : else
161 : !total moment is fixed
162 0 : m_spins=2
163 : END IF
164 :
165 343 : results%seigv = 0.0e0
166 686 : do mspin=1,m_spins
167 343 : IF (m_spins == 1) THEN
168 1029 : sslice = (/1,nspins/)
169 : ELSE
170 0 : sslice = (/mspin,mspin/)
171 0 : nspins = 1
172 : ENDIF
173 343 : n = 0
174 788 : DO jsp = sslice(1),sslice(2)
175 : !Generate a list of energies
176 5694 : DO k = 1,kpts%nkpt
177 : !
178 : !---> STORE EIGENVALUES AND WEIGHTS IN A LINEAR LIST. AND MEMORIZE
179 : !---> CONECTION TO THE ORIGINAL ARRAYS
180 : !
181 185836 : DO j = 1,results%neig(k,jsp)
182 180930 : e(n+j) = results%eig(j,k,jsp)
183 180930 : we(n+j) = kpts%wtkpt(k)
184 180930 : idxeig(n+j) = j+n_help
185 180930 : idxkpt(n+j) = k
186 185836 : idxjsp(n+j) = jsp
187 : END DO
188 : !---> COUNT THE NUMBER OF EIGENVALUES
189 5351 : n = n + results%neig(k,jsp)
190 : END DO
191 : END DO
192 :
193 343 : CALL sort(index(:n),e)
194 :
195 : ! Check if no deep eigenvalue is found
196 181273 : IF (e_min-MINVAL(e(1:n))>1.0) THEN
197 0 : if (l_output) then
198 0 : WRITE(oUnit,*) 'WARNING: Too low eigenvalue detected:'
199 0 : WRITE(oUnit,*) 'min E=', MINVAL(e(1:n)),' min(enpara)=',e_min
200 : endif
201 : CALL juDFT_warn("Too low eigenvalue detected",calledby="fermi", &
202 : hint ="If the lowest eigenvalue is more than 1Htr below "//&
203 : "the lowest energy parameter, you probably have picked up"//&
204 0 : " a ghoststate")
205 : END IF
206 : !
207 : !---> DETERMINE EF BY SUMMING WEIGHTS
208 : !
209 343 : weight = input%zelec/spindg
210 343 : seigv=0.0
211 343 : IF(m_spins /= 1) weight = weight/2.0 -(mspin-1.5)*input%fixed_moment
212 343 : ws = 0.0e0
213 343 : l = 0
214 75807 : DO WHILE ((ws+del).LT.weight)
215 75464 : l = l + 1
216 75464 : IF (l.GT.n) THEN
217 0 : IF ( fmpi%irank == 0 .and. l_output) THEN
218 0 : WRITE (oUnit,FMT=8010) n,ws,weight
219 : END IF
220 0 : CALL juDFT_error("Not enough wavefunctions",calledby="fermie")
221 : 8010 FORMAT (/,10x,'error: not enough wavefunctions.',i10,2d20.10)
222 : END IF
223 75464 : ws = ws + we(INDEX(l))
224 75464 : seigv =seigv + e(INDEX(l))*we(INDEX(l))*spindg
225 : !IF (l_output) WRITE (oUnit,FMT='(2f10.7)') e(index(l)),we(index(l))
226 : END DO
227 343 : results%ef = -100000.0
228 343 : IF(l.GT.0) THEN
229 343 : results%ef = e(INDEX(l))
230 : END IF
231 343 : nstef = l
232 343 : zc = input%zelec
233 343 : IF(m_spins /= 1) THEN
234 0 : zc = zc/2.0-(mspin-1.5)*input%fixed_moment
235 0 : idxjsp = 1 !assume single spin in following calculations
236 0 : if (l_output) then
237 0 : IF (mspin == 1) THEN
238 0 : WRITE(oUnit,*) "Fixed total moment calculation"
239 0 : WRITE(oUnit,*) "Moment:",input%fixed_moment
240 0 : write(oUnit,*) "First Spin:"
241 : ELSE
242 0 : WRITE(oUnit,*) "Second Spin:"
243 : ENDIF
244 : end if
245 : ENDIF
246 :
247 343 : IF ( fmpi%irank == 0 .and. l_output) WRITE (oUnit,FMT=8020) results%ef,nstef,seigv,ws,ssc
248 :
249 : !+po
250 343 : results%ts = 0.0
251 : !-po
252 189764 : results%w_iks(:,:,sslice(1):sslice(2)) = 0.0
253 343 : results%bandgap = 0.0
254 343 : IF(input%bz_integration==BZINT_METHOD_HIST) THEN
255 : CALL ferhis(input,kpts,fmpi,index,idxeig,idxkpt,idxjsp,nspins, n,&
256 : nstef,ws,spindg,weight,e,results%neig(:,sslice(1):sslice(2)),&
257 318 : we, noco,cell,results%ef,results%seigv,results%w_iks(:,:,sslice(1):sslice(2)),results,l_output)
258 25 : ELSE IF (input%bz_integration==BZINT_METHOD_GAUSS) THEN
259 : CALL fergwt(kpts,input,fmpi,results%neig(:,sslice(1):sslice(2)), results%eig(:,:,sslice(1):sslice(2)),&
260 0 : results%ef,results%w_iks(:,:,sslice(1):sslice(2)),results%seigv,l_output)
261 25 : ELSE IF (input%bz_integration==BZINT_METHOD_TRIA) THEN
262 : CALL fertri(input,noco,kpts,fmpi%irank, results%neig(:,sslice(1):sslice(2)),nspins,zc,results%eig(:,:,sslice(1):sslice(2)),spindg,&
263 2 : results%ef,results%seigv,results%w_iks(:,:,sslice(1):sslice(2)),l_output)
264 23 : ELSE IF (input%bz_integration==BZINT_METHOD_TETRA) THEN
265 : CALL fertetra(input,noco,kpts,fmpi,results%neig(:,sslice(1):sslice(2)), results%eig(:,:,sslice(1):sslice(2)),&
266 23 : results%ef,results%w_iks(:,:,sslice(1):sslice(2)),results%seigv,l_output)
267 : ENDIF
268 :
269 343 : IF (mspin == 2 .AND. l_output) THEN
270 0 : WRITE(oUnit,*) "Different Fermi-energies for both spins:"
271 0 : WRITE(oUnit,"(a,f0.3,a,f0.4,a,f0.4,a,f0.4)") "Fixed Moment:" &
272 0 : ,input%fixed_moment," Difference(EF):",efermi," - ",results%ef,"="&
273 0 : ,efermi-results%ef
274 : ENDIF
275 686 : efermi = results%ef
276 : enddo
277 :
278 343 : IF (m_spins == 2) nspins = 2
279 :
280 343 : IF (l_output) THEN
281 2058 : attributes = ''
282 343 : WRITE(attributes(1),'(f20.10)') results%ef
283 343 : WRITE(attributes(2),'(a)') 'Htr'
284 1029 : IF (fmpi%irank.EQ.0) CALL writeXMLElement('FermiEnergy',(/'value','units'/),attributes(1:2))
285 : END IF
286 : ENDIF
287 :
288 686 : RETURN
289 : 8020 FORMAT (/,'FERMIE:',/,&
290 : & 10x,'first approx. to ef (T=0) :',f11.6,' htr',&
291 : & ' (energy of the highest occ. eigenvalue)',/,&
292 : & 10x,'number of occ. states (T=0) :',i11,/,&
293 : & 10x,'first approx. to seigv (T=0) :',f11.6,' htr',/,&
294 : & 10x,'sum of weights of occ. states :',f11.6,/,&
295 : & 10x,'sum of semicore charge :',f11.6,' e',/)
296 : END SUBROUTINE fermie
297 : END MODULE m_fermie
|
