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_wann_amn
8 : use m_juDFT
9 : contains
10 0 : subroutine wann_amn(
11 0 : > chi,nslibd,nwfsd,ntypd,nlod,llod,llo,nlo,
12 0 : > lmaxd,jmtd,lmd,neq,natd,ikpt,nbnd,
13 0 : > rmsh,rmt,jri,dx,lmax,
14 0 : > us,dus,uds,duds,flo,
15 0 : > ff,gg,acof,bcof,ccof,l_nocosoc,jspin,
16 : & l_amn2_in,
17 0 : < amn,
18 : > bkpt)
19 : c*********************************************************************
20 : c...here: twf - trial wf Y. Mokrousov June'06
21 : c.. and August'06
22 :
23 : c...For the moment calculates only the overlaps inside the spheres.
24 :
25 : c...As an input file the routine reads the file 'proj', where
26 : c...it reads the following parameters:
27 : c...nwfs: total number of wfs
28 : c...For each atom (not atom type) the following values
29 : c...are specified:
30 : c.. lwf, mrwf - specify the angular part of the twf (see tutorial)
31 : c.. rwf - specifies the radial index of the twf
32 : c.. alpha,beta,gamma - the euler angles of the twf
33 : c.. zona - the diffusivity of the twf
34 : c.. regio - 2*Rmt for instance, will calculate the radial
35 : c.. integrals of the twf and the wavefunction in the
36 : c.. sphere around the atom center with the radius of
37 : c.. 2*Rmt. Outside the MT the integrals with the
38 : c.. radial part of the planewave (in terms of bessels)
39 : c.. will be calculated without seeing other MTs, vacuum
40 : c.. boundaries whatsoever. But for the moment works
41 : c.. only for 1*Rmt (so regio=1.0)
42 : c*********************************************************************
43 : c l_amn2_in=.true. => Variant of standard wann_amn:
44 : c Projections are defined in file pro2.1/pro2.2
45 : c Includes the possibility to define displacements
46 : c of the MT-localized
47 : c trial orbitals onto atoms in neighboring unit cells. Otherwise
48 : c equal to standard case. May
49 : c be useful for covalent systems: The proj
50 : c file specifies the localized
51 : c trial orbital for one partner of the covalent bond
52 : c while proj2 does
53 : c the same for the second. In many cases the second partner will be
54 : c located in a neighboring unit cell, however.
55 : c Frank Freimuth, March 2007
56 : c*************************************************************************
57 : c Inclusion of Noco&&soc; tidied up version
58 : c Frank Freimuth
59 : c*********************************************************************
60 : use m_intgr, only : intgr3
61 : use m_constants, only: pimach, ImagUnit
62 : use m_dwigner
63 : use m_wann_tlmw
64 : use m_wann_rad_twf
65 : use m_eulerrot
66 : implicit none
67 :
68 : integer, intent (in) :: nwfsd,nslibd,ntypd,nlod,llod
69 : integer, intent (in) :: jmtd,lmaxd,lmd,natd,ikpt,nbnd
70 : integer, intent (in) :: jri(ntypd),nlo(ntypd),llo(nlod,ntypd)
71 : integer, intent (in) :: neq(ntypd),lmax(ntypd),jspin
72 : real, intent (in) :: rmsh(jmtd,ntypd),rmt(ntypd),dx(ntypd)
73 : real, intent (in) :: us(0:lmaxd,ntypd),dus(0:lmaxd,ntypd)
74 : real, intent (in) :: uds(0:lmaxd,ntypd),duds(0:lmaxd,ntypd)
75 : real, intent (in) :: ff(ntypd,jmtd,2,0:lmaxd)
76 : real, intent (in) :: gg(ntypd,jmtd,2,0:lmaxd)
77 : real, intent (in) :: flo(ntypd,jmtd,2,nlod)
78 : complex, intent (in) :: acof(nslibd,0:lmd,natd)
79 : complex, intent (in) :: bcof(nslibd,0:lmd,natd)
80 : complex, intent (in) :: ccof(-llod:llod,nslibd,nlod,natd)
81 : logical, intent (in) :: l_nocosoc
82 : logical, intent (inout) :: l_amn2_in
83 : complex, intent (inout) :: amn(nbnd,nwfsd)
84 : real,intent(in),optional:: bkpt(3)
85 : complex, intent (in) :: chi
86 : c...local
87 : integer :: nwf,nwfs,nat,j,ntyp,ne,l,m,lm,iatom,i,mp,lo
88 0 : integer :: ind(nwfsd),ntp(natd),banddummy
89 0 : integer :: lwf(nwfsd),mrwf(nwfsd),rwf(nwfsd),spi(nwfsd)
90 0 : real :: alpha(nwfsd),beta(nwfsd),gamma(nwfsd)
91 0 : real :: jwf(nwfsd),jmwf(nwfsd),pi
92 0 : real,allocatable :: posshifts(:,:)
93 0 : real,allocatable :: weights(:)
94 0 : real :: zona(nwfsd),regio(nwfsd),amx(3,3,nwfsd)
95 : real :: rr,vl,vld,r0,vlo
96 0 : real :: rads(nwfsd,0:3,jmtd,2),vlpr(jmtd),vlprd(jmtd)
97 0 : complex :: tlmwf(0:3,-3:3,nwfsd),tlmwft(0:3,-3:3,nwfsd)
98 : logical :: l_oldproj,l_amn2,l_file
99 0 : complex :: d_wgn(-3:3,-3:3,1:3),wign(-3:3,-3:3,3,nwfsd)
100 : complex :: factor
101 : real :: arg
102 : real :: mrott(3,3),bmatt(3,3),imx(3,3)
103 : character(len=2) :: spin12(0:2)
104 : data spin12/' ', '.1', '.2'/
105 : character(len=6) :: filename
106 :
107 0 : l_amn2=l_amn2_in
108 0 : if(.not.l_amn2_in)then
109 0 : inquire(file='pro2.1',exist=l_amn2_in)
110 : endif
111 0 : pi=pimach()
112 : c..generates an array giving the atom type for each atom
113 0 : ntp(:) = 0
114 : iatom = 0
115 0 : do ntyp = 1,ntypd
116 0 : do nat = 1,neq(ntyp)
117 0 : iatom = iatom + 1
118 0 : ntp(iatom) = ntyp
119 : enddo
120 : enddo
121 :
122 0 : if(l_amn2)then
123 0 : do j=jspin,1,-1
124 0 : inquire(file=trim('pro2'//spin12(j)),exist=l_file)
125 0 : if(l_file)then
126 0 : filename='pro2'//spin12(j)
127 0 : exit
128 : endif
129 : enddo
130 : else
131 : c..reading the proj.1 / proj.2 / proj file
132 0 : do j=jspin,0,-1
133 0 : inquire(file=trim('proj'//spin12(j)),exist=l_file)
134 0 : if(l_file)then
135 0 : filename='proj'//spin12(j)
136 0 : exit
137 : endif
138 : enddo
139 : endif
140 :
141 0 : if(l_file)then
142 0 : open (203,file=trim(filename),status='old')
143 0 : rewind (203)
144 : else
145 0 : CALL juDFT_error("no proj/proj.1/proj.2",calledby="wann_amn")
146 : endif
147 :
148 0 : if(l_amn2)then
149 0 : allocate(posshifts(3,nwfsd))
150 0 : allocate(weights(nwfsd))
151 : endif
152 :
153 0 : if(l_nocosoc)then
154 0 : read (203,*)nwfs,banddummy,l_oldproj
155 0 : if(.not.l_oldproj)then
156 0 : do nwf=1,nwfs
157 0 : read(203,*)ind(nwf),lwf(nwf),jwf(nwf),jmwf(nwf),rwf(nwf)
158 0 : read(203,*)alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
159 0 : if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
160 : enddo !nwf
161 : else
162 0 : do nwf = 1,nwfs
163 : read (203,*)
164 0 : & ind(nwf),lwf(nwf),mrwf(nwf),rwf(nwf),spi(nwf)
165 : read (203,*)
166 0 : & alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
167 0 : if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
168 : enddo !nwf
169 : endif !oldproj
170 : else
171 0 : read (203,*) nwfs
172 0 : do nwf = 1,nwfs
173 : read (203,*)
174 0 : & ind(nwf),lwf(nwf),mrwf(nwf),rwf(nwf)
175 : read (203,*)
176 0 : & alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
177 0 : if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
178 : enddo !nwf
179 : endif !l_nocosoc
180 0 : rewind (203)
181 0 : close (203)
182 :
183 0 : if (ikpt.eq.1) then
184 0 : write (6,*) 'Number of trial WFs:',nwfs
185 0 : write (6,*)
186 0 : do nwf = 1,nwfs
187 0 : write (6,*) 'Twfs N:',nwf,' Atom N:',ind(nwf)
188 0 : write (6,*) 'l=',lwf(nwf),' mr=',mrwf(nwf),' r=',rwf(nwf)
189 0 : write (6,*) 'zona=',zona(nwf),' region=',regio(nwf),'*Rmt'
190 0 : write (6,*) 'alpha=',alpha(nwf),
191 0 : & ' beta=',beta(nwf),' gamma=',gamma(nwf)
192 0 : write (6,*)
193 : enddo
194 : endif
195 :
196 : c..generating the radial twf function
197 :
198 0 : rads(:,:,:,:) = 0.
199 :
200 : call wann_rad_twf(
201 : > nwfs,jmtd,natd,ind,rwf,zona,regio,
202 : > us,dus,uds,duds,ff,gg,lmaxd,ikpt,
203 : > ntypd,ntp,jri,rmsh,dx,
204 : > nlod,flo,llo,nlo,
205 0 : < rads)
206 :
207 0 : open (100,file='rads')
208 0 : do i = 1,jmtd
209 : c write (100,'(i3,2x,4f10.6)') i,rads(1,0:3,i,1)
210 0 : write (100,'(f10.6,2x,4f10.6)') rmsh(i,1),rads(1,0:3,i,1)
211 : enddo
212 0 : close(100)
213 :
214 : c..now generate the coefficients in the expansion in lattice
215 : c..harmonics of the angular part of the twf
216 :
217 0 : tlmwft(:,:,:) = cmplx(0.,0.)
218 0 : tlmwf(:,:,:) = cmplx(0.,0.)
219 :
220 0 : if(l_nocosoc.and..not.l_oldproj)then
221 0 : call soc_tlmw(nwfs,lwf,jwf,jmwf,jspin,tlmwf)
222 : else
223 : call wann_tlmw(
224 : > nwfs,lwf,mrwf,
225 : > l_nocosoc,spi,jspin,
226 0 : < tlmwf)
227 : endif
228 :
229 0 : call eulerrot(nwfs,alpha,beta,gamma,amx)
230 :
231 0 : imx(:,:) = 0. ; imx(1,1) = 1. ; imx(2,2) = 1. ; imx(3,3) = 1.
232 :
233 : c..performing the wigner rotation
234 : c..These rotations are specified in the proj file in terms of the
235 : c..euler angles. The wave functions inside the muffin-tins are represen-
236 : c..-ted in terms of the global-frame-lattice-harmonics, therefore,
237 : c..the wigner transformation has to be applied to tlmwf, resulting
238 : c..in the tlmwft matrix, which is suitable for further use
239 :
240 : c..given the euler angles, the following procedure has to be
241 : c..performed in order to match the two local frames:
242 : c..1. Rotation around the z-axis by alpha
243 : c..2. Rotation around the x-axis by beta
244 : c..3. Rotation around the z-axis by gamma again.
245 0 : call d_wigner(nwfs,amx,imx,3,wign)
246 :
247 : c..now we transform the tlmwf coefficients
248 :
249 0 : do nwf = 1,nwfs
250 0 : tlmwft(0,:,nwf) = tlmwf(0,:,nwf)
251 0 : do l = 1,3
252 0 : do m = -l,l
253 0 : do mp = -l,l
254 : tlmwft(l,m,nwf) = tlmwft(l,m,nwf) +
255 0 : + wign(mp,m,l,nwf)*tlmwf(l,mp,nwf)
256 : enddo
257 : enddo
258 : enddo
259 : enddo
260 :
261 : c..calculating the amn matrix
262 :
263 0 : vlpr(:) = 0. ; vlprd(:) = 0.
264 :
265 : c...sum by wfs, each of them is localized at a certain mt
266 0 : do nwf = 1,nwfs
267 0 : if(l_amn2)then
268 0 : arg=-bkpt(1)*posshifts(1,nwf)
269 0 : arg=arg-bkpt(2)*posshifts(2,nwf)
270 0 : arg=arg-bkpt(3)*posshifts(3,nwf)
271 0 : arg=2*pi*arg
272 0 : factor=cmplx(cos(arg),sin(arg))*weights(nwf)
273 : else
274 : factor=1.0
275 : endif
276 0 : factor = factor*chi
277 :
278 0 : nat = ind(nwf)
279 0 : ntyp = ntp(nat)
280 : c...sum by bands
281 0 : do ne = 1,nslibd
282 : c...sum by l,m
283 0 : do l = 0,min(lmax(ntyp),3)
284 0 : do j = 1,jri(ntyp)
285 : vlpr(j) = ff(ntyp,j,1,l)*rads(nwf,l,j,1)+
286 0 : + ff(ntyp,j,2,l)*rads(nwf,l,j,2)
287 : vlprd(j) = gg(ntyp,j,1,l)*rads(nwf,l,j,1) +
288 0 : + gg(ntyp,j,2,l)*rads(nwf,l,j,2)
289 0 : if (rwf(nwf).gt.0)then
290 0 : vlpr(j) = vlpr(j)*rmsh(j,ntyp)
291 0 : vlprd(j) = vlprd(j)*rmsh(j,ntyp)
292 : endif
293 :
294 : enddo
295 :
296 : c..these integrations are not necessary if rads is the lin.comb.
297 : c..of the u_l and \dot{u}_l, but are necessary for other ways of
298 : c..constructing the radial part, therefore, we do it anyway
299 :
300 0 : r0 = rmsh(1,ntyp)
301 0 : call intgr3(vlpr,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vl)
302 0 : call intgr3(vlprd,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vld)
303 0 : do m = -l,l
304 0 : lm = l*(l+1) + m
305 : amn(ne,nwf) = amn(ne,nwf) +
306 : + tlmwft(l,m,nwf)*conjg(( acof(ne,lm,nat)*vl +
307 0 : + bcof(ne,lm,nat)*vld )*(ImagUnit)**l)*factor
308 :
309 : enddo
310 : enddo
311 :
312 : c..local orbitals
313 0 : if (nlo(ntyp).ge.1) then
314 0 : do lo = 1,nlo(ntyp)
315 0 : l = llo(lo,ntyp)
316 0 : do j = 1,jri(ntyp)
317 : vlpr(j) = flo(ntyp,j,1,lo)*rads(nwf,l,j,1) +
318 0 : + flo(ntyp,j,2,lo)*rads(nwf,l,j,2)
319 0 : if (rwf(nwf).gt.0)then
320 0 : vlpr(j) = vlpr(j)*rmsh(j,ntyp)
321 : endif
322 : enddo
323 :
324 :
325 0 : r0 = rmsh(1,ntyp)
326 0 : call intgr3(vlpr,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vlo)
327 0 : do m = -l,l
328 : amn(ne,nwf) = amn(ne,nwf) +
329 : + tlmwft(l,m,nwf)*conjg((ccof(m,ne,lo,nat)*vlo)*
330 0 : * (ImagUnit)**l)*factor
331 : enddo
332 : enddo
333 : endif
334 :
335 : enddo
336 : enddo
337 :
338 0 : return
339 :
340 0 : end subroutine wann_amn
341 :
342 : c*********************************************************************
343 : c Calculate the expansion of the angular part of the trial orbital
344 : c in terms of spherical harmonics.
345 : c Version for Soc: The trial orbital is a spinor.
346 : c Frank Freimuth, June 2007
347 : c*********************************************************************
348 0 : subroutine soc_tlmw(nwfs,lwf,jwf,jmwf,jspin,tlmwf)
349 : use m_clebsch
350 : implicit none
351 :
352 : integer, intent (in) :: nwfs
353 : integer, intent (in) :: lwf(nwfs)
354 : real, intent (in) :: jwf(nwfs),jmwf(nwfs)
355 : integer, intent (in) :: jspin
356 : complex, intent (out) :: tlmwf(0:3,-3:3,nwfs)
357 :
358 : integer nwf,l,m
359 : complex tlm(0:3,-3:3,1:7)
360 : real j,jm
361 :
362 :
363 0 : do nwf=1,nwfs
364 0 : l=lwf(nwf)
365 0 : IF(l<0) CALL juDFT_error("not yet implemented",calledby
366 0 : + ="wann_amn")
367 0 : j=jwf(nwf)
368 0 : jm=jmwf(nwf)
369 0 : IF(j<0) CALL juDFT_error("jwf",calledby ="wann_amn")
370 0 : IF( ABS(jm) - j >1e-10) CALL juDFT_error("jmwf",calledby
371 0 : + ="wann_amn")
372 0 : if( abs( l + 0.5 -j ).gt. 1e-10 .and.
373 : & abs( l - 0.5 -j ).gt. 1e-10 )
374 : & CALL juDFT_error ('regula trianguli violata',
375 0 : & calledby="wann_amn")
376 0 : tlmwf(0:3,-3:3,nwf)=cmplx(0.0,0.0)
377 0 : do m=-l,l
378 0 : tlmwf(l,m,nwf)=clebsch(real(l),0.5,real(m),1.5-jspin,j,jm)
379 : enddo
380 : enddo
381 :
382 0 : end subroutine soc_tlmw
383 :
384 : end module m_wann_amn
|
