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 :
7 : MODULE m_find_enpara
8 : USE m_judft
9 : IMPLICIT NONE
10 : PRIVATE
11 : PUBLIC:: find_enpara
12 :
13 : CONTAINS
14 :
15 : !> Function to determine the energy parameter given the quantum number and the potential
16 : !! Different schemes are implemented. Nqn (main quantum number) is used as a switch.
17 : !! This code was previously in lodpot.f
18 7205 : REAL FUNCTION find_enpara(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up,l_scalar_relativi)RESULT(e)
19 : USE m_types_setup
20 : USE m_radsra
21 : USE m_differ
22 : USE m_constants
23 : IMPLICIT NONE
24 : LOGICAL,INTENT(IN):: lo
25 : INTEGER,INTENT(IN):: l,n,nqn,jsp
26 : REAL,INTENT(OUT) :: e_lo,e_up
27 : TYPE(t_atoms),INTENT(IN)::atoms
28 : REAL,INTENT(IN):: vr(:)
29 : LOGICAL, OPTIONAL, INTENT(IN) :: l_scalar_relativi
30 :
31 7205 : IF(PRESENT(l_scalar_relativi)) THEN
32 2994 : e = priv_scalar_relativi(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up)
33 2994 : RETURN
34 : END IF
35 :
36 4211 : IF (nqn>0) e=priv_method1(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up)
37 4211 : IF (nqn<0) e=priv_method2(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up)
38 : END FUNCTION find_enpara
39 :
40 :
41 4207 : REAL FUNCTION priv_method1(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up)RESULT(e)
42 : USE m_types_setup
43 : USE m_radsra
44 : USE m_differ
45 : USE m_constants
46 : IMPLICIT NONE
47 : LOGICAL,INTENT(IN):: lo
48 : INTEGER,INTENT(IN):: l,n,nqn,jsp
49 : REAL,INTENT(OUT) :: e_lo,e_up
50 : TYPE(t_atoms),INTENT(IN)::atoms
51 : REAL,INTENT(IN):: vr(:)
52 :
53 :
54 : INTEGER j,ilo,i
55 : INTEGER nodeu,node,ierr,msh
56 : REAL lnd,e1, e_up_local, e_lo_local
57 : REAL d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c
58 : ! ..
59 : ! .. Local Arrays ..
60 4207 : REAL, ALLOCATABLE :: f(:,:),vrd(:)
61 : CHARACTER(LEN=20) :: attributes(6)
62 4207 : c=c_light(1.0)
63 : !Core potential setup done for each n,l now
64 4207 : d = EXP(atoms%dx(n))
65 : ! set up core-mesh
66 4207 : rn = atoms%rmt(n)
67 4207 : msh = atoms%jri(n)
68 578967 : DO WHILE (rn < atoms%rmt(n) + 20.0)
69 574760 : msh = msh + 1
70 574760 : rn = rn*d
71 : ENDDO
72 4207 : rn = atoms%rmsh(1,n)*( d**(msh-1) )
73 21035 : ALLOCATE ( f(msh,2),vrd(msh) )
74 7134107 : f = 0.0
75 3564950 : vrd = 0.0
76 : ! extend core potential (linear with slope t1 / a.u.)
77 2990190 : vrd(:atoms%jri(n))=vr(:atoms%jri(n))
78 4207 : t1=0.125
79 4207 : t2 = vrd(atoms%jri(n))/atoms%rmt(n) - atoms%rmt(n)*t1
80 4207 : rr = atoms%rmt(n)
81 578967 : DO j = atoms%jri(n) + 1, msh
82 574760 : rr = d*rr
83 578967 : vrd(j) = rr*( t2 + rr*t1 )
84 : ENDDO
85 :
86 4207 : node = nqn - (l+1)
87 4207 : IF (node<0) CALL judft_error("Error in setup of energy-parameters",hint="This could e.g. happen if you try to use 1p-states")
88 4207 : e = 0.0
89 : ! determine upper edge
90 4207 : nodeu = -1
91 :
92 11955 : DO WHILE ( nodeu <= node )
93 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
94 7748 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
95 11955 : IF ( nodeu.LE.node ) THEN
96 3541 : e = e + 10.0
97 3541 : IF (e>1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
98 : ENDIF
99 : ENDDO
100 :
101 4207 : e_up_local = e
102 :
103 12621 : DO WHILE (nodeu.GT.node)
104 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
105 8414 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
106 12621 : IF ( nodeu.GT.node ) THEN
107 4207 : e = e - 10.0
108 4207 : IF (e<-1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
109 : ENDIF
110 : ENDDO
111 :
112 4207 : e_lo_local = e
113 :
114 88347 : DO WHILE ( (e_up_local-e_lo_local).GT.1.0e-5)
115 84140 : e = (e_up_local+e_lo_local) / 2.0
116 84140 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),c,us,dus,nodeu,f(:,1),f(:,2))
117 88347 : IF (nodeu.GT.node) THEN
118 : e_up_local = e
119 : ELSE
120 40808 : e_lo_local = e
121 : END IF
122 : END DO
123 :
124 4207 : e_up = (e_up_local+e_lo_local) / 2.0
125 :
126 4207 : IF (node /= 0) THEN
127 : ! determine lower edge
128 :
129 2461 : e_up_local = e
130 :
131 8886 : DO WHILE (nodeu.GE.node)
132 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
133 6425 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
134 6425 : e = e - 10.0
135 : END DO
136 :
137 2461 : e_lo_local = e
138 :
139 54961 : DO WHILE ( (e_up_local-e_lo_local).GT.1.0e-5)
140 52500 : e = (e_up_local+e_lo_local) / 2.0
141 52500 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),c,us,dus,nodeu,f(:,1),f(:,2))
142 54961 : IF (nodeu.GE.node) THEN
143 : e_up_local = e
144 : ELSE
145 26751 : e_lo_local = e
146 : END IF
147 : END DO
148 :
149 2461 : e_lo = (e_up_local+e_lo_local) / 2.0
150 : ELSE
151 1746 : e_lo = -9.99
152 : ENDIF
153 : ! calculate core
154 4207 : e = (e_up+e_lo)/2
155 :
156 4207 : fn = REAL(nqn) ; fl = REAL(l) ; fj = fl + 0.5
157 : CALL differ(fn,fl,fj,c,atoms%zatom(n),atoms%dx(n),atoms%rmsh(1,n),&
158 4207 : rn,d,msh,vrd, e, f(:,1),f(:,2),ierr)
159 4207 : IF (lo.AND.l>0) THEN
160 521 : e1 = (e_up+e_lo)/2
161 521 : fn = REAL(nqn) ; fl = REAL(l) ; fj = fl-0.5
162 : CALL differ(fn,fl,fj,c,atoms%zatom(n),atoms%dx(n),atoms%rmsh(1,n),&
163 521 : rn,d,msh,vrd, e1, f(:,1),f(:,2),ierr)
164 521 : e = (2.0*e + e1 ) / 3.0
165 : ENDIF
166 4207 : END FUNCTION priv_method1
167 :
168 4 : REAL FUNCTION priv_method2(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up)RESULT(e)
169 : USE m_types_setup
170 : USE m_radsra
171 : USE m_differ
172 : USE m_constants
173 : IMPLICIT NONE
174 : LOGICAL,INTENT(IN):: lo
175 : INTEGER,INTENT(IN):: l,n,nqn,jsp
176 : REAL,INTENT(OUT) :: e_lo,e_up
177 : TYPE(t_atoms),INTENT(IN)::atoms
178 : REAL,INTENT(IN):: vr(:)
179 :
180 : INTEGER j,ilo,i
181 : INTEGER nodeu,node,ierr,msh
182 : REAL lnd,e_up_temp,e_lo_temp,large_e_step
183 : REAL d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c
184 : ! ..
185 : ! .. Local Arrays ..
186 :
187 4 : REAL, ALLOCATABLE :: f(:,:),vrd(:)
188 : CHARACTER(LEN=20) :: attributes(6)
189 :
190 4 : c=c_light(1.0)
191 :
192 : !Core potential setup done for each n,l now
193 4 : d = EXP(atoms%dx(n))
194 : ! set up core-mesh
195 4 : rn = atoms%rmt(n)
196 4 : msh = atoms%jri(n)
197 576 : DO WHILE (rn < atoms%rmt(n) + 20.0)
198 572 : msh = msh + 1
199 572 : rn = rn*d
200 : ENDDO
201 4 : rn = atoms%rmsh(1,n)*( d**(msh-1) )
202 20 : ALLOCATE ( f(msh,2),vrd(msh) )
203 7212 : f = 0.0
204 3604 : vrd = 0.0
205 : ! extend core potential (linear with slope t1 / a.u.)
206 3032 : vrd(:atoms%jri(n))=vr(:atoms%jri(n))
207 4 : t1=0.125
208 4 : t2 = vrd(atoms%jri(n))/atoms%rmt(n) - atoms%rmt(n)*t1
209 4 : rr = atoms%rmt(n)
210 576 : DO j = atoms%jri(n) + 1, msh
211 572 : rr = d*rr
212 576 : vrd(j) = rr*( t2 + rr*t1 )
213 : ENDDO
214 : ! search for branches
215 4 : node = ABS(nqn) - (l+1)
216 4 : IF (node<0) CALL judft_error("Error in setup of energy-parameters",hint="This could e.g. happen if you try to use 1p-states")
217 4 : e = 0.0 ! The initial value of e is arbitrary.
218 4 : large_e_step = 5.0 ! 5.0 Htr steps for coarse energy searches
219 :
220 : ! determine upper band edge
221 : ! Step 1: Coarse search for the band edge
222 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
223 4 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
224 4 : DO WHILE ( nodeu > node )
225 0 : e = e - large_e_step
226 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
227 4 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
228 : END DO
229 8 : DO WHILE ( nodeu <= node )
230 4 : e = e + large_e_step
231 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
232 8 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
233 : END DO
234 4 : e_up_temp = e
235 4 : e_lo_temp = e - large_e_step
236 : ! Step 2: Fine band edge determination by bisection search
237 40 : DO WHILE ((e_up_temp - e_lo_temp) > 1e-2)
238 36 : e = (e_up_temp + e_lo_temp) / 2.0
239 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
240 36 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
241 40 : IF (nodeu > node) THEN
242 : e_up_temp = e
243 : ELSE
244 20 : e_lo_temp = e
245 : END IF
246 : END DO
247 4 : e_up = (e_up_temp + e_lo_temp) / 2.0
248 4 : e = e_up
249 :
250 : ! determine lower band edge
251 4 : IF (node == 0) THEN
252 0 : e_lo = -49.99
253 : ELSE
254 : ! Step 1: Coarse search for the band edge
255 4 : nodeu = node
256 8 : DO WHILE ( nodeu >= node )
257 4 : e = e - large_e_step
258 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
259 8 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
260 : ENDDO
261 4 : e_up_temp = e + large_e_step
262 4 : e_lo_temp = e
263 : ! Step 2: Fine band edge determination by bisection search
264 40 : DO WHILE ((e_up_temp - e_lo_temp) > 1e-2)
265 36 : e = (e_up_temp + e_lo_temp) / 2.0
266 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
267 36 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
268 40 : IF (nodeu < node) THEN
269 : e_lo_temp = e
270 : ELSE
271 17 : e_up_temp = e
272 : END IF
273 : END DO
274 4 : e_lo = (e_up_temp + e_lo_temp) / 2.0
275 : END IF
276 :
277 :
278 : ! determince notches by intersection
279 4 : ldmt= -99.0 !ldmt = logarithmic derivative @ MT boundary
280 4 : lnd = -l-1
281 101 : DO WHILE ( ABS(ldmt-lnd) .GE. 1E-07)
282 97 : e = (e_up+e_lo)/2
283 : CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
284 97 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
285 97 : ldmt = dus/us
286 101 : IF( ldmt .GT. lnd) THEN
287 51 : e_lo = e
288 46 : ELSE IF( ldmt .LT. lnd ) THEN
289 46 : e_up = e
290 : e_lo = e_lo
291 : END IF
292 : END DO
293 :
294 4 : END FUNCTION priv_method2
295 :
296 2994 : REAL FUNCTION priv_scalar_relativi(lo,l,n,jsp,nqn,atoms,vr,e_lo,e_up)RESULT(e)
297 : USE m_types_setup
298 : USE m_radsra
299 : USE m_differ
300 : USE m_constants
301 : IMPLICIT NONE
302 : LOGICAL,INTENT(IN):: lo
303 : INTEGER,INTENT(IN):: l,n,nqn,jsp
304 : REAL,INTENT(OUT) :: e_lo,e_up
305 : TYPE(t_atoms),INTENT(IN)::atoms
306 : REAL,INTENT(IN):: vr(:)
307 :
308 :
309 : INTEGER j,ilo,i
310 : INTEGER nodeu,node,ierr,msh
311 : REAL lnd,e1
312 : REAL d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c
313 : ! ..
314 : ! .. Local Arrays ..
315 2994 : REAL, ALLOCATABLE :: f(:,:),vrd(:)
316 : CHARACTER(LEN=20) :: attributes(6)
317 2994 : c=c_light(1.0)
318 :
319 : !Core potential setup done for each n,l now
320 2994 : d = EXP(atoms%dx(n))
321 : ! set up core-mesh
322 2994 : rn = atoms%rmt(n)
323 2994 : msh = atoms%jri(n)
324 263752 : DO WHILE (rn < atoms%rmt(n) + 8.0)
325 260758 : msh = msh + 1
326 260758 : rn = rn*d
327 : ENDDO
328 2994 : rn = atoms%rmsh(1,n)*( d**(msh-1) )
329 14970 : ALLOCATE ( f(msh,2),vrd(msh) )
330 4734042 : f = 0.0
331 2365524 : vrd = 0.0
332 : ! extend core potential (linear with slope t1 / a.u.)
333 2104766 : vrd(:atoms%jri(n))=vr(:atoms%jri(n))
334 2994 : t1=0.125
335 2994 : t2 = vrd(atoms%jri(n))/atoms%rmt(n) - atoms%rmt(n)*t1
336 2994 : rr = atoms%rmt(n)
337 263752 : DO j = atoms%jri(n) + 1, msh
338 260758 : rr = d*rr
339 263752 : vrd(j) = rr*( t2 + rr*t1 )
340 : ENDDO
341 :
342 2994 : node = nqn - (l+1)
343 2994 : IF (node<0) CALL judft_error("Error in setup of energy-parameters",hint="This could e.g. happen if you try to use 1p-states")
344 2994 : e_up = 0.0
345 : ! determine upper edge
346 2994 : nodeu = -1
347 5988 : DO WHILE ( nodeu.LE.node )
348 : CALL radsra(e_up,l,vr(:),atoms%rmsh(1,n),&
349 2994 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
350 5988 : IF (nodeu.LE.node) THEN
351 0 : e_up = e_up + 5.0
352 : END IF
353 : ENDDO
354 :
355 2994 : e_lo = e_up - 50.0
356 :
357 5994 : DO WHILE ( nodeu.GT.node )
358 : CALL radsra(e_lo,l,vr(:),atoms%rmsh(1,n),&
359 3000 : atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
360 5994 : IF (nodeu.GT.node) THEN
361 6 : e_lo = e_lo - 50.0
362 : END IF
363 : ENDDO
364 :
365 2994 : e_lo = e_lo - 0.1
366 119766 : DO WHILE ( (e_up-e_lo).GT.1.0e-10)
367 116772 : e = (e_up+e_lo) / 2.0
368 116772 : CALL radsra(e,l,vrd,atoms%rmsh(1,n),atoms%dx(n),msh,c,us,dus,nodeu,f(:,1),f(:,2))
369 119766 : IF (nodeu.GT.node) THEN
370 56650 : e_up = e
371 : ELSE
372 60122 : e_lo = e
373 : END IF
374 : END DO
375 2994 : e = (e_up+e_lo) / 2.0
376 :
377 2994 : END FUNCTION priv_scalar_relativi
378 :
379 :
380 : END MODULE m_find_enpara
|
