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_real
8 : c ********************************************************
9 : c calculates the value of the periodic part of the
10 : c wavefunction at the given real-grid point p(:)
11 : c Y.Mokrousov 16.8.6
12 : c ********************************************************
13 : CONTAINS
14 0 : SUBROUTINE wann_real(
15 : > p,n,na,iv,iflag,bkpt,iband,
16 : > n3d,nmzxyd,n2d,ntypsd,lmaxd,jmtd,
17 0 : > natd,ntypd,nmzd,nop,nop2,mrot,tau,invtab,
18 : > nq3,nvac,invs,z1,delz,nmz,nmzxy,nq2,
19 0 : > lmax,rmsh,jri,pos,ngopr,ntypsy,nvd,
20 0 : > omtil,amat,bmat,nlod,llod,nlo,llo,
21 0 : > ff,gg,flo,acof,bcof,ccof,zMat,
22 0 : > nv,k1,k2,k3,lmd,nbasfcn,l_ss,qss,jspin,addnoco,
23 : < xdnout)
24 :
25 : USE m_types
26 : USE m_ylm
27 : USE m_constants
28 : use m_juDFT
29 :
30 : IMPLICIT NONE
31 :
32 : TYPE(t_mat),INTENT(IN) :: zMat
33 :
34 : C .. Scalar Arguments ..
35 : INTEGER, INTENT (IN) :: n3d,nmzxyd,n2d,ntypsd,llod,nlod,iband
36 : INTEGER, INTENT (IN) :: lmaxd,jmtd,ntypd,natd,nmzd
37 : INTEGER, INTENT (IN) :: iflag,n,na,iv,lmd,nv,nvd,nbasfcn
38 : INTEGER, INTENT (IN) :: nq3,nvac,nmz,nmzxy,nq2,nop,nop2
39 : LOGICAL, INTENT (IN) :: invs,l_ss
40 : REAL, INTENT (IN) :: z1,delz,omtil,bkpt(3),qss(3)
41 : INTEGER, INTENT (IN) :: jspin,addnoco
42 : COMPLEX, INTENT (OUT):: xdnout
43 :
44 : C ..
45 : C .. Array Arguments ..
46 : INTEGER, INTENT (IN) :: ngopr(natd),ntypsy(natd),jri(ntypd)
47 : INTEGER, INTENT (IN) :: lmax(ntypd),mrot(3,3,nop),invtab(nop)
48 : INTEGER, INTENT (IN) :: nlo(ntypd),llo(nlod,ntypd)
49 : REAL, INTENT (IN) :: amat(3,3),bmat(3,3),pos(3,natd)
50 : REAL, INTENT (IN) :: rmsh(jmtd,ntypd),tau(3,nop)
51 : INTEGER, INTENT (IN) :: k1(nvd),k2(nvd),k3(nvd)
52 : COMPLEX, INTENT (IN) :: ccof(-llod:llod,nlod,natd)
53 : COMPLEX, INTENT (IN) :: acof(0:lmd,natd)
54 : COMPLEX, INTENT (IN) :: bcof(0:lmd,natd)
55 : REAL, INTENT (IN) :: ff(ntypd,jmtd,2,0:lmaxd)
56 : REAL, INTENT (IN) :: gg(ntypd,jmtd,2,0:lmaxd)
57 : REAL, INTENT (IN) :: flo(ntypd,jmtd,2,nlod)
58 : REAL, INTENT (INOUT) :: p(3)
59 :
60 : C ..
61 : C .. Local Scalars ..
62 : REAL delta,sx,xx1,xx2,rrr,phi,const,arg,tpi,arg1
63 : INTEGER i,j,jp3,jr,k,l,nd,nopa,ivac,lm,m,gzi,kk
64 : INTEGER kk1,kk2,kk3
65 : COMPLEX const2,s,xd1,xd2,const3
66 : C ..
67 : C .. Local Arrays ..
68 0 : COMPLEX sf2(n2d),sf3(n3d),ylm((lmaxd+1)**2)
69 : REAL rcc(3),x(3),rcc2(3)
70 : REAL bqpt(3)
71 :
72 0 : call timestart("wann_real")
73 0 : tpi = 2 * pimach()
74 0 : const = 1./(sqrt(omtil))
75 :
76 : c..define the factor e^{-ikr}
77 0 : rcc2=matmul(bmat,p)/tpi_const
78 :
79 0 : bqpt = 0.0
80 : ! if(l_ss.and.jspin.eq.1) then
81 : ! bqpt = -qss/2.0
82 : ! elseif(l_ss.and.jspin.eq.2) then
83 : ! bqpt = +qss/2.0
84 : ! endif
85 :
86 : arg = -tpi*( (bkpt(1)+bqpt(1))*rcc2(1)
87 : > + (bkpt(2)+bqpt(2))*rcc2(2)
88 0 : > + (bkpt(3)+bqpt(3))*rcc2(3) )
89 :
90 0 : arg1 = tpi*( bkpt(1)*rcc2(1) + bkpt(2)*rcc2(2) + bkpt(3)*rcc2(3) )
91 0 : const2 = cmplx(cos(arg),sin(arg))
92 : const3 = cmplx(cos(arg1),sin(arg1))
93 : c write (oUnit,*) 'bkpt,const2,const3=',bkpt(:),const2,const3
94 :
95 0 : ivac=iv
96 :
97 0 : IF (iflag.EQ.0) GO TO 20
98 0 : IF (iflag.EQ.1) GO TO 40
99 : c ---> interstitial part
100 0 : rcc=matmul(bmat,p)/tpi_const
101 0 : xdnout = cmplx(0.,0.)
102 : c write (oUnit,*) 'nv,nvd=',nv,nvd
103 0 : IF (zMat%l_real) THEN
104 0 : DO k = 1,nv
105 : c write (oUnit,*) 'k1,k2,k3=',k1(k),k2(k),k3(k)
106 : c write (oUnit,*) 'z(k,iband)=', z(k,iband)
107 0 : arg = tpi * ((k1(k))*rcc(1)+(k2(k))*rcc(2)+(k3(k))*rcc(3))
108 : xdnout = xdnout + zMat%data_r(k+addnoco,iband)*
109 0 : + cmplx(cos(arg),sin(arg))*const
110 : IF (((abs(p(1)-2.2).le.0.0001).and.(abs(p(2)).le.0.0001))
111 0 : & .or.((abs(p(2)-2.2).le.0.0001).and.(abs(p(1)).le.0.0001)))then
112 : c write (oUnit,*) 'p(i)=',p(1:2)
113 : c write (oUnit,*) 'G=',k1(k),k2(k),k3(k)
114 : c write (oUnit,*) 'z(k,iband)=',z(k,iband)
115 : c write (oUnit,*) 'val=',z(k,iband)*cmplx(cos(arg),sin(arg))
116 : ENDIF
117 : END DO
118 : ELSE
119 0 : DO k = 1,nv
120 : c write (oUnit,*) 'k1,k2,k3=',k1(k),k2(k),k3(k)
121 : c write (oUnit,*) 'z(k,iband)=', z(k,iband)
122 0 : arg = tpi * ((k1(k))*rcc(1)+(k2(k))*rcc(2)+(k3(k))*rcc(3))
123 : xdnout = xdnout + zMat%data_c(k+addnoco,iband)*
124 0 : + cmplx(cos(arg),sin(arg))*const
125 : IF (((abs(p(1)-2.2).le.0.0001).and.(abs(p(2)).le.0.0001))
126 0 : & .or.((abs(p(2)-2.2).le.0.0001).and.(abs(p(1)).le.0.0001)))then
127 : c write (oUnit,*) 'p(i)=',p(1:2)
128 : c write (oUnit,*) 'G=',k1(k),k2(k),k3(k)
129 : c write (oUnit,*) 'z(k,iband)=',z(k,iband)
130 : c write (oUnit,*) 'val=',z(k,iband)*cmplx(cos(arg),sin(arg))
131 : ENDIF
132 : END DO
133 : END IF
134 : c write (oUnit,*) 'ir:p(i)',p(:)
135 0 : call timestop("wann_real")
136 0 : RETURN
137 : c ---> vacuum part
138 : 20 CONTINUE
139 0 : xdnout = cmplx(0.,0.)
140 0 : call timestop("wann_real")
141 0 : return
142 : c ----> m.t. part
143 : 40 CONTINUE
144 :
145 :
146 0 : nd = ntypsy(na)
147 0 : nopa = ngopr(na)
148 0 : nopa=1
149 :
150 :
151 0 : sx = 0.0
152 0 : DO 50 i = 1,3
153 0 : x(i) = p(i) - pos(i,na)
154 0 : sx = sx + x(i)*x(i)
155 0 : 50 CONTINUE
156 0 : sx = sqrt(sx)
157 : IF (nopa.NE.1) THEN
158 : c... switch to internal units
159 : rcc=matmul(bmat,p)/tpi_const
160 : c... rotate into representative
161 : DO 70 i = 1,3
162 : p(i) = 0.
163 : DO 60 j = 1,3
164 : p(i) = p(i) + mrot(i,j,nopa)*rcc(j)
165 :
166 : 60 CONTINUE
167 : 70 CONTINUE
168 : c... switch back to cartesian units
169 : x=matmul(amat,p)/tpi_const
170 : END IF
171 0 : DO 80 j = jri(n),2,-1
172 0 : IF (sx.GE.rmsh(j,n)) GO TO 90
173 0 : 80 CONTINUE
174 0 : 90 jr = j
175 : CALL ylm4(
176 : > lmax(n),x,
177 0 : < ylm)
178 0 : xd1 = cmplx(0.,0.)
179 0 : xd2 = cmplx(0.,0.)
180 0 : DO l = 0,lmax(n)
181 : c if (p(1).eq.0. .and. p(2).eq.0. .and. p(3).eq.0)then
182 : c write (oUnit,*) 'ff(l,300)=',ff(1,300,1,l)
183 : c write (oUnit,*) 'ff(l,300)=',ff(1,300,2,l)
184 : c write (oUnit,*) 'gg(l,300)=',gg(1,300,1,l)
185 : c write (oUnit,*) 'gg(l,300)=',gg(1,300,2,l)
186 : c endif
187 0 : DO 110 m = -l,l
188 0 : lm = l*(l+1)+m
189 0 : s = ylm(lm+1)*(ImagUnit)**l
190 : c if (p(1).eq.0. .and. p(2).eq.0. .and. p(3).eq.0)then
191 : c write (oUnit,*) 'acof=',acof(lm,1)
192 : c write (oUnit,*) 'bcof=',bcof(lm,1)
193 : c endif
194 : xd1 = xd1 + (acof(lm,na)*cmplx(ff(n,jr,1,l),0.)+
195 : + bcof(lm,na)*cmplx(gg(n,jr,1,l),0.))*s/
196 0 : / (rmsh(jr,n))
197 : c / (rmsh(jr,n)*rmsh(jr,n))
198 0 : IF (jr.EQ.1) GO TO 110
199 : xd2 = xd2 + (acof(lm,na)*cmplx(ff(n,jr+1,1,l),0.)+
200 : + bcof(lm,na)*cmplx(gg(n,jr+1,1,l),0.))*s/
201 0 : / (rmsh(jr+1,n))
202 : c / (rmsh(jr+1,n)*rmsh(jr+1,n))
203 0 : 110 CONTINUE
204 : ENDDO
205 : c..contributions from the local orbitals
206 0 : IF (nlo(n).GE.1) THEN
207 0 : DO l = 1,nlo(n)
208 0 : DO 111 m = -llo(l,n),llo(l,n)
209 0 : lm = llo(l,n)*(llo(l,n)+1)+m
210 0 : s = ylm(lm+1)*(ImagUnit)**llo(l,n)
211 : xd1 = xd1 + ccof(m,l,na)*flo(n,jr,1,l)*s/
212 0 : / (rmsh(jr,n))
213 0 : IF (jr.EQ.1) GO TO 111
214 : xd2 = xd2 + ccof(m,l,na)*flo(n,jr+1,1,l)*s/
215 0 : / (rmsh(jr+1,n))
216 0 : 111 CONTINUE
217 : ENDDO
218 : ENDIF
219 0 : IF (jr.EQ.1) THEN
220 0 : xdnout = xd1
221 : ELSE
222 : xdnout = xd1 + (xd2-xd1) *
223 0 : + (sx-rmsh(jr,n)) / (rmsh(jr+1,n)-rmsh(jr,n))
224 :
225 : END IF
226 0 : xdnout = xdnout*const2
227 : c write (oUnit,*) 'mt:p(i)',p(:)
228 : 8000 FORMAT (2f10.6)
229 : c
230 0 : call timestop("wann_real")
231 0 : RETURN
232 : END SUBROUTINE wann_real
233 : END MODULE m_wann_real
|
