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_radflo
8 : USE m_juDFT
9 : CONTAINS
10 4928 : SUBROUTINE radflo(atoms, ntyp,jsp,ello,vr, f,g,fmpi, usdus,&
11 4928 : uuilon,duilon,ulouilopn,flo,lout_all)
12 : !
13 : !***********************************************************************
14 : ! generates the scalar relativistic wavefunctions (flo) needed for the
15 : ! local orbitals at atom type n for angular momentum l.
16 : ! the values of the function and the radial derivative on the sphere
17 : ! (ulos,dulos) boundaries, the overlap with the other radial functions
18 : ! (uulon,dulon) and between the different local orbitals (uloulopn) are
19 : ! also calculated.
20 : !
21 : ! ulos : the value of the radial function of a local orbital
22 : ! at the muffin tin radius
23 : ! dulos : the value of the radial derivative of the radial function
24 : ! function of a local orbital at the muffin tin radius
25 : ! uulon : overlap integral between the radial functions of a local
26 : ! obital and the flapw radial function with the same l
27 : ! dulon : overlap integral between the radial functions of a local
28 : ! obital and the energy derivative of the flapw radial
29 : ! function with the same l
30 : ! uloulopn: overlap integral between the radial functions of two
31 : ! different local orbitals
32 : ! (only needed if they have the same l)
33 : ! l_dulo : whether we use a dot u as local orbital (see setlomap.F)
34 : !
35 : ! p.kurz jul. 1996 gb 2001
36 : !
37 : ! ulo_der : specifies the order of the energy derivative to be used
38 : ! (0, 1, 2, ... for primitive, first, second derivatives etc.)
39 : ! uuilon : overlap integral between the radial functions of the
40 : ! integral (multiplied by ulo_der) of a local orbital and the
41 : ! flapw radial function with the same l
42 : ! duilon : overlap integral between the radial functions of the
43 : ! integral of a local orbital and the energy derivative of the
44 : ! flapw radial function with the same l
45 : ! ulouilopn: overlap integral between the radial functions of the
46 : ! integral of a local orbital and another local orbital with
47 : ! the same l.
48 : !
49 : ! C. Friedrich Feb. 2005
50 : !***********************************************************************
51 : !
52 : USE m_intgr, ONLY : intgr0
53 : USE m_constants
54 : USE m_radsra
55 : USE m_radsrdn
56 : USE m_differ
57 : USE m_types
58 : IMPLICIT NONE
59 : TYPE(t_usdus),INTENT(INOUT):: usdus !lo part is calculated here
60 : TYPE(t_mpi),INTENT(IN) :: fmpi
61 : TYPE(t_atoms),INTENT(IN) :: atoms
62 : ! ..
63 : ! .. Scalar Arguments ..
64 : INTEGER, INTENT (IN) :: jsp,ntyp
65 : LOGICAL, INTENT (IN), OPTIONAL :: lout_all
66 : ! ..
67 : ! .. Array Arguments ..
68 : REAL, INTENT (IN) :: ello(atoms%nlod,atoms%ntype),vr(atoms%jmtd)
69 : REAL, INTENT (IN) :: f(atoms%jmtd,2,0:atoms%lmaxd),g(atoms%jmtd,2,0:atoms%lmaxd)
70 : REAL, INTENT (OUT):: uuilon(atoms%nlod,atoms%ntype),duilon(atoms%nlod,atoms%ntype)
71 : REAL, INTENT (OUT):: ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)
72 : REAL, INTENT (OUT):: flo(atoms%jmtd,2,atoms%nlod)
73 : ! ..
74 : ! .. Local Scalars ..
75 : INTEGER i,j,k,l,ilo,jlo,nodelo,noded ,ierr,msh
76 : REAL rn,t1,t2,d ,rr,fn,fl,fj,e,uds,duds,ddn,c,rmt
77 : LOGICAL ofdiag, loutput
78 : ! ..
79 : ! .. Local Arrays ..
80 4928 : REAL dulo(atoms%jmtd),ulo(atoms%jmtd),glo(atoms%jmtd,2),filo(atoms%jmtd,2,atoms%nlod)
81 4928 : REAL, ALLOCATABLE :: f_rel(:,:),vrd(:)
82 4928 : REAL help(atoms%nlod+2,atoms%nlod+2)
83 : ! ..
84 4928 : c = c_light(1.0)
85 : !
86 : IF ( PRESENT(lout_all) ) THEN
87 : loutput = ( fmpi%irank == 0 ) .OR. lout_all
88 : ELSE
89 : loutput = ( fmpi%irank == 0 )
90 : END IF
91 4928 : !$ loutput=.false.
92 : IF (loutput) THEN
93 : WRITE (oUnit,FMT=8000)
94 : END IF
95 4928 : ofdiag = .FALSE.
96 : !---> calculate the radial wavefunction with the appropriate
97 : !---> energy parameter (ello)
98 13880 : DO ilo = 1,atoms%nlo(ntyp)
99 : CALL radsra(ello(ilo,ntyp),atoms%llo(ilo,ntyp),vr(:),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c,&
100 8952 : usdus%ulos(ilo,ntyp,jsp),usdus%dulos(ilo,ntyp,jsp),nodelo,flo(:,1,ilo), flo(:,2,ilo))
101 : !
102 : !+apw+lo
103 8952 : IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN
104 36 : IF (atoms%ulo_der(ilo,ntyp).LE.8) THEN
105 : !---> calculate orthogonal energy derivative at e
106 36 : i = atoms%ulo_der(ilo,ntyp)
107 36 : IF(atoms%l_dulo(ilo,ntyp)) i=1
108 : CALL radsrdn(ello(ilo,ntyp),atoms%llo(ilo,ntyp),vr(:),atoms%rmsh(1,ntyp), atoms%dx(ntyp),&
109 36 : atoms%jri(ntyp),c, uds,duds,ddn,noded,glo,filo(:,:,ilo), flo(:,:,ilo),usdus%dulos(ilo,ntyp,jsp),i)
110 27288 : DO i=1,atoms%jri(ntyp)
111 27252 : flo(i,1,ilo) = glo(i,1)
112 27288 : flo(i,2,ilo) = glo(i,2)
113 : ENDDO
114 36 : nodelo = noded
115 36 : ddn = SQRT(ddn)
116 36 : IF(atoms%l_dulo(ilo,ntyp)) ddn=1.0
117 54612 : flo (:,:,ilo) = flo (:,:,ilo)/ddn ! Normalize ulo (flo) if APW+lo is not used
118 54612 : filo(:,:,ilo) = filo(:,:,ilo)/ddn ! and scale its integral (filo) accordingly
119 36 : usdus%dulos(ilo,ntyp,jsp) = duds/ddn ! (setabc1lo and slomat assume <flo|flo>=1)
120 36 : usdus%ulos (ilo,ntyp,jsp) = uds/ddn !
121 : ELSE
122 : !
123 : ! test:
124 : !
125 : ! set up core-mesh
126 0 : d = EXP(atoms%dx(ntyp))
127 0 : rmt = atoms%rmsh(1,ntyp)
128 0 : DO i = 1, atoms%jri(ntyp) - 1
129 0 : rmt = rmt * d
130 : ENDDO
131 0 : rn = rmt
132 0 : msh = atoms%jri(ntyp)
133 0 : DO WHILE (rn < rmt + 20.0)
134 0 : msh = msh + 1
135 0 : rn = rn*d
136 : ENDDO
137 0 : rn = atoms%rmsh(1,ntyp)*( d**(msh-1) )
138 0 : ALLOCATE ( f_rel(msh,2),vrd(msh) )
139 :
140 : ! extend core potential (linear with slope t1 / a.u.)
141 :
142 0 : DO j = 1, atoms%jri(ntyp)
143 0 : vrd(j) = vr(j)
144 : ENDDO
145 0 : t1=0.125
146 0 : t2 = vrd(atoms%jri(ntyp))/rmt - rmt*t1
147 0 : rr = rmt
148 0 : DO j = atoms%jri(ntyp) + 1, msh
149 0 : rr = d*rr
150 0 : vrd(j) = rr*( t2 + rr*t1 )
151 : ENDDO
152 0 : e = ello(ilo,ntyp)
153 0 : fn = 6.0 ; fl = 1.0 ; fj = 0.5
154 : CALL differ(fn,fl,fj,c,82.0,atoms%dx(ntyp),atoms%rmsh(1,ntyp),&
155 0 : rn,d,msh,vrd, e, f_rel(1,1),f_rel(1,2),ierr)
156 :
157 0 : f_rel(:,1) = 2 * f_rel(:,1)
158 0 : f_rel(:,2) = 2 * f_rel(:,2)
159 0 : rn = atoms%rmsh(1,ntyp)
160 0 : DO i = 1, atoms%jri(ntyp)
161 0 : rn = rn * d
162 0 : WRITE(123,'(5f20.15)') rn,f_rel(i,1),f_rel(i,2), f(i,1,1),f(i,2,1)
163 : ENDDO
164 :
165 0 : STOP
166 : ENDIF
167 :
168 : ENDIF
169 : !-apw+lo
170 : !
171 : !---> calculate the overlap between these fcn. and the radial functions
172 : !---> of the flapw basis with the same l
173 7130792 : DO i = 1,atoms%jri(ntyp)
174 7121840 : ulo(i) = f(i,1,atoms%llo(ilo,ntyp))*flo(i,1,ilo) + f(i,2,atoms%llo(ilo,ntyp))*flo(i,2,ilo)
175 7130792 : dulo(i) = g(i,1,atoms%llo(ilo,ntyp))*flo(i,1,ilo) + g(i,2,atoms%llo(ilo,ntyp))*flo(i,2,ilo)
176 : END DO
177 8952 : CALL intgr0(ulo, atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%uulon(ilo,ntyp,jsp))
178 8952 : CALL intgr0(dulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%dulon(ilo,ntyp,jsp))
179 8952 : IF (atoms%l_dulo(ilo,ntyp)) usdus%dulon(ilo,ntyp,jsp) = 0.0
180 : IF (loutput) THEN
181 : WRITE (oUnit,FMT=8010) ilo,atoms%llo(ilo,ntyp),ello(ilo,ntyp),&
182 : usdus%ulos(ilo,ntyp,jsp),usdus%dulos(ilo,ntyp,jsp),nodelo,usdus%uulon(ilo,ntyp,jsp),&
183 : usdus%dulon(ilo,ntyp,jsp)
184 : END IF
185 : !
186 : !---> case LO = energy derivative (ulo_der>=1):
187 : !---> calculate the overlap between the LO-integral (filo) and the radial functions
188 8952 : IF(atoms%ulo_der(ilo,ntyp).GE.1) THEN
189 27288 : DO i=1,atoms%jri(ntyp)
190 27252 : ulo(i) = f(i,1,atoms%llo(ilo,ntyp))*filo(i,1,ilo) + f(i,2,atoms%llo(ilo,ntyp))*filo(i,2,ilo)
191 27288 : dulo(i) = g(i,1,atoms%llo(ilo,ntyp))*filo(i,1,ilo) + g(i,2,atoms%llo(ilo,ntyp))*filo(i,2,ilo)
192 : ENDDO
193 36 : CALL intgr0(ulo, atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),uuilon(ilo,ntyp))
194 36 : CALL intgr0(dulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),duilon(ilo,ntyp))
195 : ELSE
196 8916 : uuilon(ilo,ntyp)=0
197 8916 : duilon(ilo,ntyp)=0
198 : ENDIF
199 : !---> calculate overlap between radial fcn. of different local
200 : !---> orbitals (only if both have the same l)
201 : ! uloulopn(ilo,ilo,ntyp) = 1.0
202 : ! DO jlo = 1, (ilo-1)
203 27000 : DO jlo = 1, ilo
204 22072 : IF (atoms%llo(ilo,ntyp).EQ.atoms%llo(jlo,ntyp)) THEN
205 7158080 : DO i = 1,atoms%jri(ntyp)
206 7158080 : ulo(i) = flo(i,1,ilo)*flo(i,1,jlo) + flo(i,2,ilo)*flo(i,2,jlo)
207 : END DO
208 8988 : CALL intgr0(ulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%uloulopn(ilo,jlo,ntyp,jsp))
209 8988 : usdus%uloulopn(jlo,ilo,ntyp,jsp) = usdus%uloulopn(ilo,jlo,ntyp,jsp)
210 8988 : ofdiag = .TRUE.
211 : ELSE
212 4132 : usdus%uloulopn(ilo,jlo,ntyp,jsp) = 0.0
213 4132 : usdus%uloulopn(jlo,ilo,ntyp,jsp) = 0.0
214 : END IF
215 : END DO
216 : END DO
217 : !
218 : !---> case: one of LOs = energy derivative (ulo_der>=1):
219 : !---> calculate overlap between LOs and integrals of LOs
220 13880 : DO ilo = 1,atoms%nlo(ntyp)
221 31168 : DO jlo = 1,atoms%nlo(ntyp)
222 26240 : IF(atoms%ulo_der(jlo,ntyp).GE.1.AND.atoms%llo(ilo,ntyp).EQ.atoms%llo(jlo,ntyp)) THEN
223 27288 : DO i = 1,atoms%jri(ntyp)
224 27288 : ulo(i) = flo(i,1,ilo)*filo(i,1,jlo) + flo(i,2,ilo)*filo(i,2,jlo)
225 : ENDDO
226 36 : CALL intgr0(ulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),ulouilopn(ilo,jlo,ntyp))
227 : ELSE
228 17252 : ulouilopn(ilo,jlo,ntyp)=0.0
229 : ENDIF
230 : ENDDO
231 : ENDDO
232 : !
233 : IF ( (ofdiag).AND.(loutput) ) THEN
234 : WRITE (oUnit,FMT=*) 'overlap matrix between different local orbitals'
235 : WRITE (oUnit,FMT=8020) (i,i=1,atoms%nlo(ntyp))
236 : DO ilo = 1,atoms%nlo(ntyp)
237 : WRITE (oUnit,FMT='(i3,40e13.6)') ilo, (usdus%uloulopn(ilo,jlo,ntyp,jsp),jlo=1,atoms%nlo(ntyp))
238 :
239 : END DO
240 : END IF
241 : !
242 : !
243 : ! Diagonalize overlap matrix of normalized MT functions
244 : ! to check linear dependencies
245 : ! help is overlap matrix of all MT functions (flapw & LO)
246 : IF (.FALSE.) THEN
247 : DO l=0,MAXVAL(atoms%llo(1:atoms%nlo(ntyp),ntyp))
248 : IF(ALL(atoms%llo(1:atoms%nlo(ntyp),ntyp).NE.l)) CYCLE
249 : help(1,1)=1.0
250 : help(1,2)=0.0
251 : help(2,1)=0.0
252 : help(2,2)=1.0
253 : DO i=1,atoms%jri(ntyp)
254 : ulo(i)=g(i,1,l)**2+g(i,2,l)**2
255 : ENDDO
256 : CALL intgr0(ulo,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),ddn)
257 : ddn=SQRT(ddn)
258 : j=2
259 : k=2
260 : DO ilo=1,atoms%nlo(ntyp)
261 : IF(atoms%llo(ilo,ntyp).EQ.l) THEN
262 : j=j+1
263 : help(1,j)=usdus%uulon(ilo,ntyp,jsp)
264 : help(2,j)=usdus%dulon(ilo,ntyp,jsp)/ddn
265 : k=2
266 : DO jlo=1,ilo
267 : IF(atoms%llo(jlo,ntyp).EQ.l) THEN
268 : k=k+1
269 : help(k,j)=usdus%uloulopn(ilo,jlo,ntyp,jsp)
270 : ENDIF
271 : ENDDO
272 : ENDIF
273 : ENDDO
274 : IF ( loutput ) THEN
275 : WRITE(oUnit,'(A,I2)')&
276 : & 'Overlap matrix of normalized MT functions '//&
277 : & 'and its eigenvalues for l =',l
278 : DO i=1,j
279 : WRITE(oUnit,'(30F11.7)') (help(k,i),k=1,i)
280 : ENDDO
281 : END IF
282 : CALL dsyev('N','U',j,help,atoms%nlod+2,ulo,dulo, i)
283 : IF(i/=0) CALL juDFT_error("ssyev failed.",calledby ="radflo")
284 : IF ( loutput ) THEN
285 : WRITE(oUnit,'(30F11.7)') (ulo(i),i=1,j)
286 : END IF
287 : ENDDO
288 : ENDIF ! .false.
289 : !
290 :
291 : 8000 FORMAT (/,t20,'radial function for local orbitals',/,t2,'lo',t6,&
292 : & 'l',t11,'energy',t29,'value',t42,'derivative',t56,'nodes',&
293 : & t63,'ovlp with u',t78,'ovlp with udot')
294 : 8010 FORMAT (i3,i3,f10.5,5x,1p,2e16.7,i5,2e16.7)
295 : 8020 FORMAT (' lo',i9,19i15)
296 4928 : RETURN
297 : END SUBROUTINE radflo
298 : END MODULE m_radflo
|
