Line data Source code
1 : MODULE m_efnewton
2 : use m_juDFT
3 : CONTAINS
4 134 : RECURSIVE SUBROUTINE ef_newton(
5 : > n,irank,
6 134 : > inkem,nocst,index,tkb,e,
7 134 : X w_near_ef,ef,we,recursion_level)
8 :
9 : c***********************************************************************
10 : c
11 : c This subroutine adjusts the Fermi-Energy to obtain charge
12 : c neutrality. This is done by solving
13 : c
14 : c ( sum ff(e) * we ) - w_near_ef = 0
15 : c e
16 : c
17 : c using the Newton-Method.
18 : c Here the sum is over the eigenvalues between ef-8kt and ef+8kt,
19 : c ff is the Fermi-Function, we is the weight of the according
20 : c k-point and w_near_ef is the weight required between ef-8kt
21 : c and ef+8kt to obtain neutrality.
22 : c
23 : c***********************************************************************
24 :
25 : USE m_constants
26 :
27 : IMPLICIT NONE
28 :
29 : C .. Scalar Arguments ..
30 : INTEGER, INTENT (IN) :: n
31 : INTEGER, INTENT (IN) :: inkem,nocst,irank
32 : REAL, INTENT (IN) :: tkb
33 : REAL, INTENT (INOUT) :: ef,w_near_ef
34 : INTEGER,INTENT(IN),OPTIONAL:: recursion_level
35 : C ..
36 : C .. Array Arguments ..
37 : INTEGER, INTENT (IN) :: index(n)
38 : REAL, INTENT (IN) :: e(:) !(nkptd*neigd*jspd)
39 : REAL, INTENT (INOUT) :: we(:) !(nkptd*neigd*jspd)
40 : C ..
41 : C .. Local Scalars ..
42 : REAL dff,expo,sdff,sff,efmin,efmax
43 : INTEGER icnt,idim,rec_level
44 : LOGICAL efminok,efmaxok
45 :
46 : C ..
47 : C .. Local Arrays ..
48 134 : REAL ff(size(e))
49 : C ..
50 : C .. Parameters ..
51 : REAL,PARAMETER:: eps=1.0e-10
52 : C ..
53 : c***********************************************************************
54 : c ABBREVIATIONS
55 : c
56 : c e : linear list of the eigenvalues within the highest
57 : c energy-window
58 : c we : list of weights of the eigenvalues in e
59 : c ef : fermi energy determined to obtain charge neutrality
60 : c tkb : value of temperature (kt) broadening around fermi
61 : c energy in htr units
62 : c index : index list, e(index(n)) is the list of the
63 : c eigenvalues in ascending order
64 : c ikem : number of eigenvalues below ef-8kt
65 : c nocst : number of eigenvalues below ef+8kt
66 : c w_near_ef : weight (charge/spindg) between ef-8kt and ef+8kt
67 : c needed to obtain charge neutrality
68 : c
69 : c***********************************************************************
70 :
71 134 : IF (PRESENT(recursion_level)) THEN
72 0 : rec_level=recursion_level+1
73 0 : IF (rec_level>20) CALL juDFT_error
74 : + ("Determination of fermi-level did not converge",hint=
75 : + 'change temp. broad., set gauss = T, or use finer k mesh',
76 0 : + calledby="ef_newton")
77 : ELSE
78 134 : rec_level=0
79 134 : IF ( irank == 0 ) WRITE (oUnit,FMT='(/,5x,''EF_NEWTON: '',
80 134 : +''Adjust Fermi-Energy by Newton-Method.'',/)')
81 : ENDIF
82 : c
83 : c ---> NEWTON CYCLE
84 : c
85 134 : efminok= .false.
86 134 : efmaxok= .false.
87 422 : DO icnt = 1,50
88 422 : sff = 0.0
89 422 : sdff = 0.0
90 6960 : DO idim = inkem + 1,nocst
91 : c
92 : c ---> COMPUTE FERMI-FUNCTION ff AND THE DERIVATIVE dff
93 : c
94 6538 : expo= -ABS(e(index(idim))-ef)/tkb
95 6538 : expo= EXP(expo)
96 6538 : IF (e(index(idim))<ef) THEN
97 2556 : ff(idim) = 1./ (expo+1.)
98 : ELSE
99 3982 : ff(idim)= expo/ (expo+1.)
100 : ENDIF
101 6538 : dff= (1.+expo)**2
102 6538 : dff= expo/(dff*tkb)
103 : c
104 : c ---> MULTIPLY WITH THE WEIGHT
105 : c
106 6538 : ff(idim) = we(index(idim))*ff(idim)
107 6538 : dff = we(index(idim))*dff
108 : c
109 : c ---> AND SUM IT UP
110 : c
111 6538 : sff = sff + ff(idim)
112 6960 : sdff = sdff + dff
113 : END DO
114 422 : sff = sff - w_near_ef
115 422 : IF (abs(sff).LT.eps) THEN
116 : !Converged, so do some output and return
117 134 : w_near_ef = sff + w_near_ef
118 134 : IF (irank == 0) WRITE (oUnit,FMT=8010) icnt,sff,-sff/sdff
119 :
120 1870 : DO idim = inkem + 1,nocst
121 1870 : we(index(idim)) = ff(idim)
122 : END DO
123 :
124 : 8000 FORMAT (15x,'ef_newton failed after :',i3,'iterations.',/,
125 : + 15x,'The error in the weight is : ',e12.5,/,
126 : + 15x,'The error in ef is : ',e12.5,' htr',/)
127 : 8010 FORMAT (15x,'Number of iterations needed : ',i3,/,
128 : + 15x,'The error in the weight is : ',e12.5,/,
129 : + 15x,'The error in ef is : ',e12.5,' htr',/)
130 134 : RETURN
131 : ENDIF
132 :
133 288 : IF (abs(sdff).LT.1e-29) THEN
134 0 : if (irank==0) THEN
135 0 : write(oUnit,*) "Instability in determination of fermi-level,"
136 0 : write(oUnit,*) "doubled temperature broading to continue"
137 0 : write(*,*) "Instability in determination of fermi-level,"
138 0 : write(*,*) "doubled temperature broading to continue"
139 : ENDIF
140 : CALL ef_newton(
141 : > n,irank,
142 : > inkem,nocst,index,tkb*2.0,e,
143 0 : X w_near_ef,ef,we,rec_level)
144 0 : RETURN
145 : ENDIF
146 288 : IF (sff > 0.) THEN
147 158 : efmax= ef
148 158 : efmaxok= .true.
149 : ELSE
150 130 : efmin= ef
151 130 : efminok= .true.
152 : ENDIF
153 288 : ef = ef - sff/sdff
154 : ! if the Newton method is not stable, nested intervals are used
155 288 : IF ( efminok .and. efmaxok ) THEN
156 53 : IF ( (ef<efmin) .or. (ef>efmax) ) THEN
157 2 : ef= (efmin+efmax)/2.
158 : ENDIF
159 : ENDIF
160 : END DO
161 : c
162 : c--- > NOT CONVERGED AFTER 50 ITERATIONS
163 : c
164 0 : IF ( irank == 0 ) WRITE (oUnit,FMT=8000) icnt,sff,-sff/sdff
165 0 : ef=ef+0.001
166 : CALL ef_newton(
167 : > n,irank,
168 : > inkem,nocst,index,tkb*2.0,e,
169 0 : X w_near_ef,ef,we,rec_level)
170 :
171 :
172 : END SUBROUTINE ef_newton
173 : END MODULE m_efnewton
|
