find_enpara.f90 10.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_find_enpara
  USE m_judft
  IMPLICIT NONE
  PRIVATE
  CHARACTER(len=1),PARAMETER,DIMENSION(0:9):: ch=(/'s','p','d','f','g','h','i','j','k','l'/)
  PUBLIC:: find_enpara

CONTAINS
 
  !> Function to determine the energy parameter given the quantum number and the potential
  !! Different schemes are implemented. Nqn (main quantum number) is used as a switch.
  !! This code was previously in lodpot.f
  REAL FUNCTION find_enpara(lo,l,n,jsp,nqn,atoms,mpi,vr)RESULT(e)
    USE m_types_setup
    USE m_types_mpi
    USE m_radsra
    USE m_differ
    USE m_xmlOutput
    USE m_constants
    IMPLICIT NONE
    LOGICAL,INTENT(IN):: lo
    INTEGER,INTENT(IN):: l,n,nqn,jsp
    TYPE(t_atoms),INTENT(IN)::atoms
    TYPE(t_mpi),INTENT(IN)  ::mpi
    REAL,INTENT(IN):: vr(:)

    IF (nqn>0) e=priv_method1(lo,l,n,jsp,nqn,atoms,mpi,vr)
    IF (nqn<0) e=priv_method2(lo,l,n,jsp,nqn,atoms,mpi,vr)
  END FUNCTION find_enpara


  REAL FUNCTION priv_method1(lo,l,n,jsp,nqn,atoms,mpi,vr)RESULT(e)
    USE m_types_setup
    USE m_types_mpi
    USE m_radsra
    USE m_differ
    USE m_xmlOutput
    USE m_constants
    IMPLICIT NONE
    LOGICAL,INTENT(IN):: lo
    INTEGER,INTENT(IN):: l,n,nqn,jsp
    TYPE(t_atoms),INTENT(IN)::atoms
    TYPE(t_mpi),INTENT(IN)  ::mpi
    REAL,INTENT(IN):: vr(:)


    INTEGER j,ilo,i
    INTEGER nodeu,node,ierr,msh
    REAL   e_up,e_lo,lnd,e1
    REAL   d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c
    LOGICAL start
    !     ..
    !     .. Local Arrays .. 
    REAL, ALLOCATABLE :: f(:,:),vrd(:)
    CHARACTER(LEN=20)    :: attributes(6)
    c=c_light(1.0)

    !Core potential setup done for each n,l now 
    d = EXP(atoms%dx(n))
    ! set up core-mesh
    rn = atoms%rmt(n)
    msh = atoms%jri(n)
    DO WHILE (rn < atoms%rmt(n) + 20.0)
       msh = msh + 1
       rn = rn*d
    ENDDO
    rn = atoms%rmsh(1,n)*( d**(msh-1) )
    ALLOCATE ( f(msh,2),vrd(msh) )
75
76
    f = 0.0
    vrd = 0.0
77
78
79
80
81
82
83
84
85
86
87
    ! extend core potential (linear with slope t1 / a.u.)
    vrd(:atoms%jri(n))=vr(:atoms%jri(n))
    t1=0.125
    t2 = vrd(atoms%jri(n))/atoms%rmt(n) - atoms%rmt(n)*t1
    rr = atoms%rmt(n)
    DO j = atoms%jri(n) + 1, msh
       rr = d*rr
       vrd(j) = rr*( t2 + rr*t1 )
    ENDDO

    node = nqn - (l+1)
88
    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")
89
90
91
92
93
94
95
96
    e = 0.0 
    ! determine upper edge
    nodeu = -1 ; start = .TRUE.
    DO WHILE ( nodeu <= node ) 
       CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
            atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
       IF  ( ( nodeu > node ) .AND. start ) THEN
          e = e - 1.0
97
          IF (e<-1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
98
99
100
          nodeu = -1
       ELSE
          e = e + 0.01
101
          IF (e>1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
102
103
104
105
106
107
108
109
110
111
112
113
          start = .FALSE.
       ENDIF
    ENDDO

    e_up = e
    IF (node /= 0) THEN
       ! determine lower edge
       nodeu = node + 1
       DO WHILE ( nodeu >= node ) 
          CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
               atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
          e = e - 0.01
114
          IF (e<-1E10) CALL judft_error("Determination of energy parameters did not converge",hint="Perhaps your potential is broken?")
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
       ENDDO
       e_lo = e
    ELSE
       e_lo = -9.99 
    ENDIF
    ! calculate core
    e  = (e_up+e_lo)/2
    fn = REAL(nqn) ; fl = REAL(l) ; fj = fl + 0.5
    CALL differ(fn,fl,fj,c,atoms%zatom(n),atoms%dx(n),atoms%rmsh(1,n),&
         rn,d,msh,vrd, e, f(:,1),f(:,2),ierr)
    IF (lo.AND.l>0) THEN
       e1  = (e_up+e_lo)/2
       fn = REAL(nqn) ; fl = REAL(l) ; fj = fl-0.5
       CALL differ(fn,fl,fj,c,atoms%zatom(n),atoms%dx(n),atoms%rmsh(1,n),&
            rn,d,msh,vrd, e1, f(:,1),f(:,2),ierr)
       e = (2.0*e + e1 ) / 3.0      
    ENDIF
 

    IF (mpi%irank  == 0) THEN
       attributes = ''
       WRITE(attributes(1),'(i0)') n
       WRITE(attributes(2),'(i0)') jsp
       WRITE(attributes(3),'(i0,a1)') nqn, ch(l)
       WRITE(attributes(4),'(f8.2)') e_lo
       WRITE(attributes(5),'(f8.2)') e_up
       WRITE(attributes(6),'(f16.10)') e
       IF (lo) THEN
          CALL writeXMLElementForm('loAtomicEP',(/'atomType     ','spin         ','branch       ',&
               'branchLowest ','branchHighest','value        '/),&
               attributes,RESHAPE((/10,4,6,12,13,5,6,1,3,8,8,16/),(/6,2/)))
       ELSE
          CALL writeXMLElementForm('atomicEP',(/'atomType     ','spin         ','branch       ',&
               'branchLowest ','branchHighest','value        '/),&
               attributes,RESHAPE((/12,4,6,12,13,5,6,1,3,8,8,16/),(/6,2/)))
       ENDIF
Uliana Alekseeva's avatar
Uliana Alekseeva committed
151
       WRITE(6,'(a6,i5,i2,a1,a12,f6.2,a3,f6.2,a13,f8.4)') '  Atom',n,nqn,ch(l),' branch from',&
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
            e_lo, ' to',e_up,' htr. ; e_l =',e
    ENDIF
  END FUNCTION priv_method1

  REAL FUNCTION priv_method2(lo,l,n,jsp,nqn,atoms,mpi,vr)RESULT(e)
    USE m_types_setup
    USE m_types_mpi
    USE m_radsra
    USE m_differ
    USE m_xmlOutput
    USE m_constants
    IMPLICIT NONE
    LOGICAL,INTENT(IN):: lo
    INTEGER,INTENT(IN):: l,n,nqn,jsp
    TYPE(t_atoms),INTENT(IN)::atoms
    TYPE(t_mpi),INTENT(IN)  ::mpi
    REAL,INTENT(IN):: vr(:)

    INTEGER j,ilo,i
    INTEGER nodeu,node,ierr,msh
    REAL   e_up,e_lo,lnd,e_up_temp,e_lo_temp,large_e_step
    REAL   d,rn,fl,fn,fj,t2,rr,t1,ldmt,us,dus,c
    LOGICAL start
    !     ..
    !     .. Local Arrays .. 

    REAL, ALLOCATABLE :: f(:,:),vrd(:)
    CHARACTER(LEN=20)    :: attributes(6)
    
    c=c_light(1.0)

    !Core potential setup done for each n,l now 
    d = EXP(atoms%dx(n))
    ! set up core-mesh
    rn = atoms%rmt(n)
    msh = atoms%jri(n)
    DO WHILE (rn < atoms%rmt(n) + 20.0)
       msh = msh + 1
       rn = rn*d
    ENDDO
    rn = atoms%rmsh(1,n)*( d**(msh-1) )
    ALLOCATE ( f(msh,2),vrd(msh) )
194
195
    f = 0.0
    vrd = 0.0
196
197
198
199
200
201
202
203
204
205
206
    ! extend core potential (linear with slope t1 / a.u.)
    vrd(:atoms%jri(n))=vr(:atoms%jri(n))
    t1=0.125
    t2 = vrd(atoms%jri(n))/atoms%rmt(n) - atoms%rmt(n)*t1
    rr = atoms%rmt(n)
    DO j = atoms%jri(n) + 1, msh
       rr = d*rr
       vrd(j) = rr*( t2 + rr*t1 )
    ENDDO
    ! search for branches
    node = ABS(nqn) - (l+1)
207
    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")
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
    e = 0.0 ! The initial value of e is arbitrary.
    large_e_step = 5.0 ! 5.0 Htr steps for coarse energy searches

    ! determine upper band edge
    ! Step 1: Coarse search for the band edge
    CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
         atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
    DO WHILE ( nodeu > node )
       e = e - large_e_step
       CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
            atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
    END DO
    DO WHILE ( nodeu <= node )
       e = e + large_e_step
       CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
            atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
    END DO
    e_up_temp = e
    e_lo_temp = e - large_e_step
    ! Step 2: Fine band edge determination by bisection search
    DO WHILE ((e_up_temp - e_lo_temp) > 1e-2)
       e = (e_up_temp + e_lo_temp) / 2.0
       CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
            atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
       IF (nodeu > node) THEN
          e_up_temp = e
       ELSE
          e_lo_temp = e
       END IF
    END DO
    e_up = (e_up_temp + e_lo_temp) / 2.0
    e    = e_up

    ! determine lower band edge
    IF (node == 0) THEN
       e_lo = -49.99
    ELSE
       ! Step 1: Coarse search for the band edge
       nodeu = node
       DO WHILE ( nodeu >= node )
          e = e - large_e_step
          CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
               atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
       ENDDO
       e_up_temp = e + large_e_step
       e_lo_temp = e
       ! Step 2: Fine band edge determination by bisection search
       DO WHILE ((e_up_temp - e_lo_temp) > 1e-2)
          e = (e_up_temp + e_lo_temp) / 2.0
          CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
               atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
          IF (nodeu < node) THEN
             e_lo_temp = e
          ELSE
             e_up_temp = e
          END IF
       END DO
       e_lo = (e_up_temp + e_lo_temp) / 2.0
    END IF


    ! determince notches by intersection
    ldmt= -99.0 !ldmt = logarithmic derivative @ MT boundary
    lnd = -l-1
    DO WHILE ( ABS(ldmt-lnd) .GE. 1E-07) 
       e = (e_up+e_lo)/2
       CALL radsra(e,l,vr(:),atoms%rmsh(1,n),&
            atoms%dx(n),atoms%jri(n),c, us,dus,nodeu,f(:,1),f(:,2))
       ldmt = dus/us
       IF( ldmt .GT. lnd) THEN
          e_lo = e
       ELSE IF( ldmt .LT. lnd ) THEN
          e_up = e
          e_lo = e_lo
       END IF
    END DO

    IF (mpi%irank == 0) THEN
       attributes = ''
       WRITE(attributes(1),'(i0)') n
       WRITE(attributes(2),'(i0)') jsp
       WRITE(attributes(3),'(i0,a1)') ABS(nqn), ch(l)
       WRITE(attributes(4),'(f16.10)') ldmt
       WRITE(attributes(5),'(f16.10)') e
       IF (lo) THEN
          CALL writeXMLElementForm('heloAtomicEP',(/'atomType      ','spin          ','branch        ',&
               'logDerivMT    ','value         '/),&
               attributes(1:5),reshape((/8,4,6,12,5+17,6,1,3,16,16/),(/5,2/)))
       ELSE
          CALL writeXMLElementForm('heAtomicEP',(/'atomType      ','spin          ','branch        ',&
               'logDerivMT    ','value         '/),&
               attributes(1:5),reshape((/10,4,6,12,5+17,6,1,3,16,16/),(/5,2/)))
       ENDIF
       WRITE (6,'(a7,i3,i2,a1,a12,f7.2,a4,f7.2,a5)') "  Atom ",n,nqn,ch(l)," branch, D = ",&
            ldmt, " at ",e," htr."
    ENDIF
  END FUNCTION priv_method2
  
END MODULE m_find_enpara