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 8 : subroutine wann_amn(
11 8 : > chi,nslibd,nwfsd,ntypd,nlod,llod,llo,nlo,
12 8 : > lmaxd,jmtd,lmd,neq,natd,ikpt,nbnd,
13 8 : > rmsh,rmt,jri,dx,lmax,
14 8 : > us,dus,uds,duds,flo,
15 8 : > ff,gg,acof,bcof,ccof,l_nocosoc,jspin,
16 : & l_amn2_in,
17 8 : < 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
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 8 : integer :: ind(nwfsd),ntp(natd),banddummy
89 8 : integer :: lwf(nwfsd),mrwf(nwfsd),rwf(nwfsd),spi(nwfsd)
90 8 : real :: alpha(nwfsd),beta(nwfsd),gamma(nwfsd)
91 8 : real :: jwf(nwfsd),jmwf(nwfsd),pi
92 8 : real,allocatable :: posshifts(:,:)
93 8 : real,allocatable :: weights(:)
94 8 : real :: zona(nwfsd),regio(nwfsd),amx(3,3,nwfsd)
95 : real :: rr,vl,vld,r0,vlo
96 8 : real :: rads(nwfsd,0:3,jmtd,2),vlpr(jmtd),vlprd(jmtd)
97 8 : complex :: tlmwf(0:3,-3:3,nwfsd),tlmwft(0:3,-3:3,nwfsd)
98 : logical :: l_oldproj,l_amn2,l_file
99 8 : 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 8 : call timestart("wann_amn")
108 72 : spi(:) = 0
109 :
110 8 : l_amn2=l_amn2_in
111 8 : if(.not.l_amn2_in)then
112 8 : inquire(file='pro2.1',exist=l_amn2_in)
113 : endif
114 8 : pi=pimach()
115 : c..generates an array giving the atom type for each atom
116 8 : call timestart("gen ntp")
117 24 : ntp(:) = 0
118 : iatom = 0
119 16 : do ntyp = 1,ntypd
120 32 : do nat = 1,neq(ntyp)
121 16 : iatom = iatom + 1
122 24 : ntp(iatom) = ntyp
123 : enddo
124 : enddo
125 8 : call timestop("gen ntp")
126 :
127 8 : call timestart("inquire pro2")
128 8 : if(l_amn2)then
129 0 : do j=jspin,1,-1
130 0 : inquire(file=trim('pro2'//spin12(j)),exist=l_file)
131 0 : if(l_file)then
132 0 : filename='pro2'//spin12(j)
133 0 : exit
134 : endif
135 : enddo
136 : else
137 : c..reading the proj.1 / proj.2 / proj file
138 16 : do j=jspin,0,-1
139 16 : inquire(file=trim('proj'//spin12(j)),exist=l_file)
140 16 : if(l_file)then
141 8 : filename='proj'//spin12(j)
142 8 : exit
143 : endif
144 : enddo
145 : endif
146 8 : call timestop("inquire pro2")
147 :
148 8 : if(l_file)then
149 8 : open (203,file=trim(filename),status='old')
150 8 : rewind (203)
151 : else
152 0 : CALL juDFT_error("no proj/proj.1/proj.2",calledby="wann_amn")
153 : endif
154 :
155 8 : if(l_amn2)then
156 0 : allocate(posshifts(3,nwfsd))
157 0 : allocate(weights(nwfsd))
158 : endif
159 :
160 8 : call timestart("read proj 203")
161 8 : if(l_nocosoc)then
162 0 : read (203,*)nwfs,banddummy,l_oldproj
163 0 : if(.not.l_oldproj)then
164 0 : do nwf=1,nwfs
165 0 : read(203,*)ind(nwf),lwf(nwf),jwf(nwf),jmwf(nwf),rwf(nwf)
166 0 : read(203,*)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 : else
170 0 : do nwf = 1,nwfs
171 : read (203,*)
172 0 : & ind(nwf),lwf(nwf),mrwf(nwf),rwf(nwf),spi(nwf)
173 : read (203,*)
174 0 : & alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
175 0 : if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
176 : enddo !nwf
177 : endif !oldproj
178 : else
179 8 : read (203,*) nwfs
180 72 : do nwf = 1,nwfs
181 : read (203,*)
182 64 : & ind(nwf),lwf(nwf),mrwf(nwf),rwf(nwf)
183 : read (203,*)
184 64 : & alpha(nwf),beta(nwf),gamma(nwf),zona(nwf),regio(nwf)
185 72 : if(l_amn2) read(203,*) posshifts(1:3,nwf),weights(nwf)
186 : enddo !nwf
187 : endif !l_nocosoc
188 8 : rewind (203)
189 8 : close (203)
190 8 : call timestop("read proj 203")
191 :
192 8 : call timestart("output trail WFs")
193 8 : if (ikpt.eq.1) then
194 1 : write (oUnit,*) 'Number of trial WFs:',nwfs
195 1 : write (oUnit,*)
196 9 : do nwf = 1,nwfs
197 8 : write (oUnit,*) 'Twfs N:',nwf,' Atom N:',ind(nwf)
198 8 : write (oUnit,*) 'l=',lwf(nwf),' mr=',mrwf(nwf),' r=',rwf(nwf)
199 8 : write (oUnit,*) 'zona=',zona(nwf),' region=',regio(nwf),'*Rmt'
200 8 : write (oUnit,*) 'alpha=',alpha(nwf),
201 16 : & ' beta=',beta(nwf),' gamma=',gamma(nwf)
202 9 : write (oUnit,*)
203 : enddo
204 : endif
205 8 : call timestop("output trail WFs")
206 :
207 : c..generating the radial twf function
208 :
209 278856 : rads(:,:,:,:) = 0.
210 :
211 : call wann_rad_twf(
212 : > nwfs,jmtd,natd,ind,rwf,zona,regio,
213 : > us,dus,uds,duds,ff,gg,lmaxd,ikpt,
214 : > ntypd,ntp,jri,rmsh,dx,
215 : > nlod,flo,llo,nlo,
216 8 : < rads)
217 :
218 8 : call timestart("write rads")
219 8 : open (100,file='rads')
220 3776 : do i = 1,jmtd
221 : c write (100,'(i3,2x,4f10.6)') i,rads(1,0:3,i,1)
222 3776 : write (100,'(f10.6,2x,4f10.6)') rmsh(i,1),rads(1,0:3,i,1)
223 : enddo
224 8 : close(100)
225 8 : call timestop("write rads")
226 :
227 : c..now generate the coefficients in the expansion in lattice
228 : c..harmonics of the angular part of the twf
229 8 : call timestart("generate coeffs in latham")
230 2312 : tlmwft(:,:,:) = cmplx(0.,0.)
231 2312 : tlmwf(:,:,:) = cmplx(0.,0.)
232 :
233 8 : if(l_nocosoc.and..not.l_oldproj)then
234 0 : call soc_tlmw(nwfs,lwf,jwf,jmwf,jspin,tlmwf)
235 : else
236 : call wann_tlmw(
237 : > nwfs,lwf,mrwf,
238 : > l_nocosoc,spi,jspin,
239 8 : < tlmwf)
240 : endif
241 :
242 8 : call eulerrot(nwfs,alpha,beta,gamma,amx)
243 :
244 8 : imx(:,:) = 0. ; imx(1,1) = 1. ; imx(2,2) = 1. ; imx(3,3) = 1.
245 :
246 : c..performing the wigner rotation
247 : c..These rotations are specified in the proj file in terms of the
248 : c..euler angles. The wave functions inside the muffin-tins are represen-
249 : c..-ted in terms of the global-frame-lattice-harmonics, therefore,
250 : c..the wigner transformation has to be applied to tlmwf, resulting
251 : c..in the tlmwft matrix, which is suitable for further use
252 :
253 : c..given the euler angles, the following procedure has to be
254 : c..performed in order to match the two local frames:
255 : c..1. Rotation around the z-axis by alpha
256 : c..2. Rotation around the x-axis by beta
257 : c..3. Rotation around the z-axis by gamma again.
258 8 : call d_wigner(nwfs,amx,imx,3,wign)
259 :
260 : c..now we transform the tlmwf coefficients
261 :
262 72 : do nwf = 1,nwfs
263 512 : tlmwft(0,:,nwf) = tlmwf(0,:,nwf)
264 264 : do l = 1,3
265 1216 : do m = -l,l
266 6464 : do mp = -l,l
267 : tlmwft(l,m,nwf) = tlmwft(l,m,nwf) +
268 6272 : + wign(mp,m,l,nwf)*tlmwf(l,mp,nwf)
269 : enddo
270 : enddo
271 : enddo
272 : enddo
273 8 : call timestop("generate coeffs in latham")
274 :
275 : c..calculating the amn matrix
276 :
277 7544 : vlpr(:) = 0. ; vlprd(:) = 0.
278 :
279 : c...sum by wfs, each of them is localized at a certain mt
280 8 : call timestart("sum by wfs")
281 72 : do nwf = 1,nwfs
282 64 : if(l_amn2.and.present(bkpt))then
283 0 : arg=-bkpt(1)*posshifts(1,nwf)
284 0 : arg=arg-bkpt(2)*posshifts(2,nwf)
285 0 : arg=arg-bkpt(3)*posshifts(3,nwf)
286 0 : arg=2*pi*arg
287 0 : factor=cmplx(cos(arg),sin(arg))*weights(nwf)
288 : else
289 : factor=1.0
290 : endif
291 64 : factor = factor*chi
292 :
293 64 : nat = ind(nwf)
294 64 : ntyp = ntp(nat)
295 : c...sum by bands
296 584 : do ne = 1,nslibd
297 : c...sum by l,m
298 2560 : do l = 0,min(lmax(ntyp),3)
299 966656 : do j = 1,jri(ntyp)
300 : vlpr(j) = ff(ntyp,j,1,l)*rads(nwf,l,j,1)+
301 964608 : + ff(ntyp,j,2,l)*rads(nwf,l,j,2)
302 : vlprd(j) = gg(ntyp,j,1,l)*rads(nwf,l,j,1) +
303 964608 : + gg(ntyp,j,2,l)*rads(nwf,l,j,2)
304 966656 : if (rwf(nwf).gt.0)then
305 0 : vlpr(j) = vlpr(j)*rmsh(j,ntyp)
306 0 : vlprd(j) = vlprd(j)*rmsh(j,ntyp)
307 : endif
308 :
309 : enddo
310 :
311 : c..these integrations are not necessary if rads is the lin.comb.
312 : c..of the u_l and \dot{u}_l, but are necessary for other ways of
313 : c..constructing the radial part, therefore, we do it anyway
314 :
315 2048 : r0 = rmsh(1,ntyp)
316 2048 : call intgr3(vlpr,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vl)
317 2048 : call intgr3(vlprd,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vld)
318 10752 : do m = -l,l
319 8192 : lm = l*(l+1) + m
320 : amn(ne,nwf) = amn(ne,nwf) +
321 : + tlmwft(l,m,nwf)*conjg(( acof(ne,lm,nat)*vl +
322 10240 : + bcof(ne,lm,nat)*vld )*(ImagUnit)**l)*factor
323 :
324 : enddo
325 : enddo
326 :
327 : c..local orbitals
328 576 : if (nlo(ntyp).ge.1) then
329 0 : do lo = 1,nlo(ntyp)
330 0 : l = llo(lo,ntyp)
331 0 : do j = 1,jri(ntyp)
332 : vlpr(j) = flo(ntyp,j,1,lo)*rads(nwf,l,j,1) +
333 0 : + flo(ntyp,j,2,lo)*rads(nwf,l,j,2)
334 0 : if (rwf(nwf).gt.0)then
335 0 : vlpr(j) = vlpr(j)*rmsh(j,ntyp)
336 : endif
337 : enddo
338 :
339 :
340 0 : r0 = rmsh(1,ntyp)
341 0 : call intgr3(vlpr,rmsh(1,ntyp),dx(ntyp),jri(ntyp),vlo)
342 0 : do m = -l,l
343 : amn(ne,nwf) = amn(ne,nwf) +
344 : + tlmwft(l,m,nwf)*conjg((ccof(m,ne,lo,nat)*vlo)*
345 0 : * (ImagUnit)**l)*factor
346 : enddo
347 : enddo
348 : endif
349 :
350 : enddo
351 : enddo
352 8 : call timestop("sum by wfs")
353 8 : call timestop("wann_amn")
354 8 : return
355 :
356 8 : end subroutine wann_amn
357 :
358 : c*********************************************************************
359 : c Calculate the expansion of the angular part of the trial orbital
360 : c in terms of spherical harmonics.
361 : c Version for Soc: The trial orbital is a spinor.
362 : c Frank Freimuth, June 2007
363 : c*********************************************************************
364 0 : subroutine soc_tlmw(nwfs,lwf,jwf,jmwf,jspin,tlmwf)
365 : use m_clebsch
366 : implicit none
367 :
368 : integer, intent (in) :: nwfs
369 : integer, intent (in) :: lwf(nwfs)
370 : real, intent (in) :: jwf(nwfs),jmwf(nwfs)
371 : integer, intent (in) :: jspin
372 : complex, intent (out) :: tlmwf(0:3,-3:3,nwfs)
373 :
374 : integer nwf,l,m
375 : complex tlm(0:3,-3:3,1:7)
376 : real j,jm
377 :
378 0 : call timestart("soc_tlmw")
379 0 : do nwf=1,nwfs
380 0 : l=lwf(nwf)
381 0 : IF(l<0) CALL juDFT_error("not yet implemented",calledby
382 0 : + ="wann_amn")
383 0 : j=jwf(nwf)
384 0 : jm=jmwf(nwf)
385 0 : IF(j<0) CALL juDFT_error("jwf",calledby ="wann_amn")
386 0 : IF( ABS(jm) - j >1e-10) CALL juDFT_error("jmwf",calledby
387 0 : + ="wann_amn")
388 0 : if( abs( l + 0.5 -j ).gt. 1e-10 .and.
389 : & abs( l - 0.5 -j ).gt. 1e-10 )
390 : & CALL juDFT_error ('regula trianguli violata',
391 0 : & calledby="wann_amn")
392 0 : tlmwf(0:3,-3:3,nwf)=cmplx(0.0,0.0)
393 0 : do m=-l,l
394 0 : tlmwf(l,m,nwf)=clebsch(real(l),0.5,real(m),1.5-jspin,j,jm)
395 : enddo
396 : enddo
397 :
398 0 : call timestop("soc_tlmw")
399 0 : end subroutine soc_tlmw
400 :
401 : end module m_wann_amn
|
