Line data Source code
1 : MODULE m_wann_nocoplot
2 : USE m_juDFT
3 : c ++++++++++++++++++++++++++++++++++++++++++++++++
4 : c + If noco=t, transform wave functions within +
5 : c + mu-th MT back to the global frame to plot +
6 : c + to plot WFs in a meaningful way. +
7 : c + +
8 : c ++++++++++++++++++++++++++++++++++++++++++++++++
9 : CONTAINS
10 0 : SUBROUTINE wann_nocoplot(atoms,slice,nnne,amat,bmat
11 : > ,nkpts,film,natd,ntypd,jmtd
12 0 : > ,ntype,neq,pos,jri,rmsh,alph,beta
13 0 : > ,nqpts,qss,z1,zatom)
14 : ! *****************************************************
15 : USE m_types
16 : USE m_xsf_io
17 : USE m_wann_gwf_tools, ONLY : get_index_q
18 : USE m_constants
19 :
20 : IMPLICIT NONE
21 :
22 : TYPE(t_atoms),INTENT(IN):: atoms
23 :
24 : LOGICAL, INTENT(IN) :: slice,film
25 : INTEGER, INTENT(IN) :: nnne,nkpts,ntypd,jmtd,natd
26 : INTEGER, INTENT(IN) :: ntype,neq(ntypd)
27 : INTEGER, INTENT(IN) :: jri(ntypd),nqpts
28 : REAL, INTENT(IN) :: amat(3,3),bmat(3,3),qss(3,nqpts)
29 : REAL, INTENT(IN) :: rmsh(jmtd,ntypd)
30 : REAL, INTENT(IN) :: pos(3,natd),z1,zatom(:)
31 : REAL, INTENT(IN) :: alph(ntypd),beta(ntypd)
32 :
33 :
34 : LOGICAL :: twodim,xsf,cartesian
35 : INTEGER :: nbmin,nbmax,nplot,nplo,nbn
36 : INTEGER :: nslibd,nkqpts
37 : INTEGER :: ix,iy,iz,ikpt,jspin
38 : INTEGER :: ii1,ii2,ii3
39 : INTEGER :: i1,i2,i3,iintsp
40 : INTEGER :: gx,gy,gz,help_kpt
41 : INTEGER :: na,nt,nq,iqpt
42 : INTEGER :: grid(3),count_mt,count_int,count_vac
43 : REAL :: vec1(3),vec2(3),vec3(3),zero(3)
44 : REAL :: pt(3),point(3),rcc2(3)
45 : REAL :: s,arg,pi,u_r,u_i
46 : COMPLEX :: phasfac
47 : COMPLEX :: U(2,2)
48 : COMPLEX :: wf_local(2),wf_global,xdnout
49 : CHARACTER(len=30) :: filename
50 : CHARACTER(len=20) :: vandername,name1,name2,name3
51 :
52 : NAMELIST /plot/twodim,cartesian,vec1,vec2,vec3,grid,zero,filename
53 :
54 : intrinsic real,aimag,conjg,exp,cmplx
55 :
56 0 : call timestart("wann_nocoplot")
57 0 : pi = pimach()
58 0 : nkqpts = nkpts*nqpts
59 :
60 0 : INQUIRE(file ="plot_inp",exist=twodim)
61 0 : IF(.NOT.twodim) THEN
62 : CALL juDFT_error("Need the plot_inp from RNK generation!"
63 0 : > ,calledby="wann_nocoplot")
64 : ENDIF
65 :
66 : !<-- Open the plot_inp file for input
67 0 : OPEN (18,file='plot_inp')
68 0 : READ(18,'(i2,5x,l1)') nplot,xsf
69 :
70 0 : IF (nplot.ge.2)
71 : & CALL juDFT_error
72 : + ("plots one by one, please, this is not charge density"
73 0 : + ,calledby="wann_nocoplot")
74 :
75 : ! the defaults
76 0 : twodim = .TRUE.;cartesian=.TRUE.;grid=(/100,100,100/)
77 0 : vec1 = (/0.,0.,0./);vec2=(/0.,0.,0./);vec3=(/0.,0.,0./)
78 0 : zero = (/0.,0.,0./);filename="default"
79 0 : READ(18,plot)
80 0 : IF (twodim.AND.ANY(grid(1:2)<1))
81 : + CALL juDFT_error("Illegal grid size in plot",calledby
82 0 : + ="wann_nocoplot")
83 0 : IF (.NOT.twodim.AND.ANY(grid<1))
84 : + CALL juDFT_error("Illegal grid size in plot",calledby
85 0 : + ="wann_nocoplot")
86 0 : IF (twodim) grid(3) = 1
87 : !calculate cartesian coordinates if needed
88 0 : IF (.NOT.cartesian) THEN
89 0 : vec1=matmul(amat,vec1)
90 0 : vec2=matmul(amat,vec2)
91 0 : vec3=matmul(amat,vec3)
92 0 : zero=matmul(amat,zero)
93 : ENDIF
94 0 : IF (filename =="default") WRITE(filename,'(a,i2)') "plot",1
95 0 : CLOSE(18)
96 :
97 : ! loop over k-points
98 0 : DO ikpt=1,nkqpts
99 0 : count_mt=0; count_int=0; count_vac=0
100 0 : iqpt = get_index_q(ikpt,nkpts)
101 0 : write(*,*)'kq',ikpt,' q',iqpt,' qz',qss(3,iqpt)
102 :
103 : ! open the old RNK.ikpt.jspin files
104 0 : DO jspin=1,2 ! local frame axis (inside MT)
105 0 : WRITE(vandername,202) ikpt,jspin
106 0 : OPEN(400+jspin,file=vandername)
107 0 : READ(400+jspin,7)gx,gy,gz,help_kpt,nslibd
108 : ENDDO
109 :
110 0 : IF(.not.xsf) THEN
111 : ! open the new UNK.ikpt.iintsp files
112 0 : DO iintsp=1,2 ! global frame axis (INT and vacuum)
113 0 : WRITE(vandername,201) ikpt,iintsp
114 0 : OPEN(500+iintsp,file=vandername,status='unknown')
115 0 : WRITE(500+iintsp,7)grid(1),grid(2),grid(3),ikpt,nslibd
116 : ENDDO
117 : ENDIF
118 :
119 : ! loop over all bands
120 0 : nbmin=1
121 0 : nbmax=nslibd
122 0 : bands:DO nbn = nbmin,nbmax
123 :
124 0 : IF (xsf) THEN
125 0 : do jspin=1,2
126 0 : write (name1,22) ikpt,nbn,jspin
127 : 22 format (i5.5,'.',i3.3,'.real.',i1,'.xsf')
128 0 : write (name2,23) ikpt,nbn,jspin
129 : 23 format (i5.5,'.',i3.3,'.imag.',i1,'.xsf')
130 0 : write (name3,24) ikpt,nbn,jspin
131 : 24 format (i5.5,'.',i3.3,'.absv.',i1,'.xsf')
132 0 : OPEN(600+jspin,file=name1)
133 0 : CALL xsf_WRITE_atoms(600+jspin,atoms,film,amat)
134 0 : OPEN(602+jspin,file=name2)
135 0 : CALL xsf_WRITE_atoms(602+jspin,atoms,film,amat)
136 0 : OPEN(604+jspin,file=name3)
137 0 : CALL xsf_WRITE_atoms(604+jspin,atoms,film,amat)
138 : CALL xsf_WRITE_header(600+jspin,twodim,filename,(vec1),
139 : & (vec2),(vec3),zero
140 0 : $ ,grid)
141 : CALL xsf_WRITE_header(602+jspin,twodim,filename,(vec1),
142 : & (vec2),(vec3),zero
143 0 : $ ,grid)
144 : CALL xsf_WRITE_header(604+jspin,twodim,filename,(vec1),
145 : & (vec2),(vec3),zero
146 0 : $ ,grid)
147 : enddo
148 : ENDIF
149 :
150 : ! loop over real space grid points
151 0 : DO iz = 0,grid(3)-1
152 0 : DO iy = 0,grid(2)-1
153 0 : xloop:DO ix = 0,grid(1)-1
154 : point = zero+vec1*REAL(ix)/grid(1)+vec2*REAL(iy)
155 0 : $ /grid(2)
156 0 : IF (.NOT.twodim) point = point+vec3*REAL(iz)/grid(3)
157 :
158 : ! read old RNK information at grid point
159 0 : DO jspin=1,2
160 0 : READ(400+jspin,8)u_r,u_i
161 0 : wf_local(jspin) = cmplx(u_r,u_i)
162 : ENDDO
163 :
164 : ! is point in MT?
165 0 : ii1 = 3
166 0 : ii2 = 3
167 0 : ii3 = 3
168 0 : IF (film ) ii3 = 0
169 :
170 0 : DO i1 = -ii1,ii1
171 0 : DO i2 = -ii2,ii2
172 0 : DO i3 = -ii3,ii3
173 0 : pt = point+MATMUL(amat,(/i1,i2,i3/))
174 0 : na = 0
175 0 : DO nt = 1,ntype
176 0 : DO nq = 1,neq(nt)
177 0 : na = na + 1
178 0 : s = SQRT(dot_PRODUCT(pos(:,na)-pt,pos(:,na)-pt))
179 0 : IF (s<rmsh(jri(nt),nt)) THEN
180 0 : count_mt=count_mt+1
181 : ! we are inside the MT with alph(nt),beta(nt)
182 : ! set up transformation local -> global
183 0 : U(1,1) = exp(-ImagUnit*alph(nt)/2.)*cos(beta(nt)/2.)
184 0 : U(1,2) = -exp(-ImagUnit*alph(nt)/2.)*sin(beta(nt)/2.)
185 0 : U(2,1) = exp( ImagUnit*alph(nt)/2.)*sin(beta(nt)/2.)
186 0 : U(2,2) = exp( ImagUnit*alph(nt)/2.)*cos(beta(nt)/2.)
187 :
188 : ! transform wfs to global frame
189 0 : DO iintsp=1,2
190 0 : pt = matmul(bmat,rcc2)/tpi_const
191 : arg = -pi*real(2*iintsp-3)*( qss(1,iqpt)*rcc2(1)
192 : > +qss(2,iqpt)*rcc2(2)
193 0 : > +qss(3,iqpt)*rcc2(3))
194 0 : phasfac = cmplx(cos(arg),sin(arg))
195 :
196 : wf_global= phasfac*( U(iintsp,1)*wf_local(1)
197 0 : > + U(iintsp,2)*wf_local(2) )
198 :
199 0 : if(xsf) THEN
200 0 : xdnout=wf_global
201 0 : WRITE(600+iintsp,*) real(xdnout)
202 0 : WRITE(602+iintsp,*) aimag(xdnout)
203 0 : WRITE(604+iintsp,*) real(xdnout*conjg(xdnout))
204 : ELSE
205 0 : WRITE(500+iintsp,8)real(wf_global),
206 0 : > aimag(wf_global)
207 : ENDIF
208 : ENDDO
209 :
210 0 : CYCLE xloop
211 : ENDIF
212 : ENDDO
213 : ENDDO !nt
214 : ENDDO
215 : ENDDO
216 : ENDDO !i1
217 :
218 : ! VACUUM region
219 0 : IF (SQRT((pt(1))**2+(pt(2))**2)>=z1)THEN
220 0 : count_vac=count_vac+1
221 0 : DO iintsp=1,2
222 0 : phasfac=cmplx(1.,0.)
223 0 : wf_global = phasfac*wf_local(iintsp)
224 :
225 0 : if(xsf) THEN
226 0 : xdnout=wf_global
227 0 : WRITE(600+iintsp,*) real(xdnout)
228 0 : WRITE(602+iintsp,*) aimag(xdnout)
229 0 : WRITE(604+iintsp,*) real(xdnout*conjg(xdnout))
230 : ELSE
231 0 : WRITE(500+iintsp,8)real(wf_global),
232 0 : > aimag(wf_global)
233 : ENDIF
234 : ENDDO
235 :
236 : CYCLE xloop
237 : ENDIF
238 :
239 : ! if we are here, point is in INTERSTITIAL
240 : ! therefore just copy wfs
241 0 : count_int = count_int+1
242 0 : DO iintsp=1,2
243 0 : phasfac=cmplx(1.,0.)
244 0 : wf_global = phasfac*wf_local(iintsp)
245 :
246 0 : if(xsf) THEN
247 0 : xdnout=wf_global
248 0 : WRITE(600+iintsp,*) real(xdnout)
249 0 : WRITE(602+iintsp,*) aimag(xdnout)
250 0 : WRITE(604+iintsp,*) real(xdnout*conjg(xdnout))
251 : ELSE
252 0 : WRITE(500+iintsp,8)real(wf_global),
253 0 : > aimag(wf_global)
254 : ENDIF
255 : ENDDO
256 :
257 : ENDDO xloop
258 : ENDDO
259 : ENDDO !z-loop
260 :
261 0 : IF (xsf) THEN
262 0 : DO iintsp=1,2
263 0 : CALL xsf_WRITE_endblock(600+iintsp,twodim)
264 0 : CALL xsf_WRITE_endblock(602+iintsp,twodim)
265 0 : CALL xsf_WRITE_endblock(604+iintsp,twodim)
266 0 : CLOSE(600+iintsp); CLOSE(602+iintsp); CLOSE(604+iintsp)
267 : ENDDO
268 : ENDIF
269 :
270 : ENDDO bands
271 :
272 0 : IF(.not.xsf) THEN
273 : ! close new UNK.ikpt.iintsp files
274 0 : DO iintsp=1,2
275 0 : CLOSE(500+iintsp)
276 : ENDDO
277 : ENDIF
278 :
279 : ! delete RNK.ikpt.jspin files
280 0 : DO jspin=1,2
281 0 : IF(xsf) THEN
282 0 : CLOSE(400+jspin)
283 : ELSE
284 0 : CLOSE(400+jspin,status='delete')
285 : ENDIF
286 : ENDDO
287 :
288 0 : write(*,*)count_mt,count_int,count_vac
289 :
290 : ENDDO ! end ikpt loop
291 :
292 : 201 FORMAT ('UNK',i5.5,'.',i1)
293 : 202 FORMAT ('RNK',i5.5,'.',i1)
294 : 7 FORMAT (5i4)
295 : 8 FORMAT (f20.12,1x,f20.12)
296 :
297 0 : call timestop("wann_nocoplot")
298 0 : END SUBROUTINE wann_nocoplot
299 : !------------------------------------------
300 :
301 : END MODULE m_wann_nocoplot
|