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_plot
8 : use m_juDFT
9 : ! +++++++++++++++++++++++++++++++++++++++++++++++++
10 : ! plots the periodic part of the wavefunction at the given k-point
11 : ! and specified bands. Needed for real-space
12 : ! visualization of the Wannier Functions
13 : !
14 : ! if twodim = .false. a 3-D plot with nplot layers in z-direction
15 : ! is constructed; the 3x3 matrix gives the 3 vectors of the cell ..
16 : ! .gustav
17 : !
18 : ! Changed the input/output for easier use.
19 : ! This subroutine uses the file plot_inp for input.
20 : ! based on the routine by: Juelich, 21.1.06 D.Wortmann
21 : ! Y.Mokrousov 16.8.6 Hamburg
22 : !
23 : ! First of all, if xsf then the files are written in the following
24 : ! format: k.bnd.real.jsp.xsf and k.bnd.imag.jsp.xsf
25 : ! for reading with XCrysDen,where
26 : ! k is the kpoint, bnd is the band index
27 : ! imag/real stand for the real and imaginary parts of wavefn
28 : ! jsp stands for spin, for all the k-points and bands
29 : ! If the xsf is false, than the files are written in the
30 : ! format suitable for the Vanderbilt's code:
31 : ! UNK0000k.jsp, having all the bands in this file, real and imag
32 : !
33 : ! Then, if xsf=true, then:
34 : !
35 : ! If slice=T, the kk specifies the k-point and the nnne
36 : ! specifies the number of the band, the wavefunction is
37 : ! written to the file kk.nnne.real.jsp.xsf and kk.nnne.imag.jsp.xsf
38 : ! +++++++++++++++++++++++++++++++++++++++++++++++++
39 : CONTAINS
40 0 : SUBROUTINE wann_plot(
41 : > vacuum,stars,cell,atoms,
42 : > nv2d,jspin,n3d,nmzxyd,n2d,ntypsd,
43 0 : > ntype,lmaxd,jmtd,ntypd,natd,nmzd,neq,nq3,nvac,
44 : > nmz,nmzxy,nq2,nop,nop2,volint,film,slice,symor,
45 0 : > invs,invs2,z1,delz,ngopr,ntypsy,jri,pos,zatom,
46 0 : > lmax,mrot,tau,rmsh,invtab,amat,bmat,bbmat,ikpt,nnne,kk,
47 0 : > nvd,nlod,llod,nv,lmd,bkpt,omtil,nlo,llo,k1,k2,k3,evac,vz,
48 0 : > nslibd,nbasfcn,neigd,ff,gg,flo,acof,bcof,ccof,zMat,
49 : > k1d,k2d,k3d,ig,ig2,sk2,phi2,l_noco,l_ss,qss,addnoco,
50 : > index_kq,l_sgwf)
51 : ! *****************************************************
52 : USE m_constants
53 : USE m_types
54 : USE m_wann_real
55 : USE m_xsf_io
56 : USE m_wann_plot_vac
57 : USE m_wann_2dvacabcof
58 :
59 : IMPLICIT NONE
60 :
61 :
62 :
63 : TYPE(t_vacuum),INTENT(IN) :: vacuum
64 : TYPE(t_stars),INTENT(IN) :: stars
65 : TYPE(t_cell),INTENT(IN) :: cell
66 : TYPE(t_atoms),INTENT(IN) :: atoms
67 : TYPE(t_mat),INTENT(IN) :: zMat
68 :
69 : ! .. Scalar Arguments ..
70 : INTEGER, INTENT (IN) :: n3d,nmzxyd,n2d,ntypsd,ikpt,jspin,nv2d
71 : INTEGER, INTENT (IN) :: lmaxd,jmtd,ntypd,natd,nmzd
72 : INTEGER, INTENT (IN) :: nq3,nvac,nmz,nmzxy,nq2,nop,nop2,ntype
73 : INTEGER, INTENT (IN) :: nvd,nv,lmd,llod,nlod
74 : INTEGER, INTENT (IN) :: nslibd,nbasfcn,neigd
75 : INTEGER, INTENT (IN) :: nnne,kk,addnoco,index_kq
76 : LOGICAL, INTENT (IN) :: symor,invs,slice,invs2,film
77 : LOGICAL, INTENT (IN) :: l_noco,l_ss,l_sgwf
78 : REAL, INTENT (IN) :: z1,delz,volint,omtil
79 : REAL, INTENT (IN) :: qss(3)
80 : real, intent (in) :: vz(nmzd,2)
81 : ! ..
82 : ! .. Array Arguments ..
83 : INTEGER, INTENT (IN) :: ngopr(natd),ntypsy(natd),lmax(ntypd)
84 : INTEGER, INTENT (IN) :: jri(ntypd),neq(ntypd),mrot(3,3,nop)
85 : INTEGER, INTENT (IN) :: invtab(nop),nlo(ntypd),llo(nlod,ntypd)
86 : INTEGER, INTENT (IN) :: k1(nvd),k2(nvd),k3(nvd)
87 : REAL, INTENT (IN) :: zatom(:),amat(3,3),bmat(3,3),pos(3,natd)
88 : REAL, INTENT (IN) :: rmsh(jmtd,ntypd),tau(3,nop),bkpt(3)
89 : REAL, INTENT (IN) :: ff(ntypd,jmtd,2,0:lmaxd),bbmat(3,3)
90 : REAL, INTENT (IN) :: gg(ntypd,jmtd,2,0:lmaxd),evac(2)
91 : REAL, INTENT (IN) :: flo(ntypd,jmtd,2,nlod)
92 : COMPLEX, INTENT (IN) :: ccof(-llod:llod,nslibd,nlod,natd)
93 : COMPLEX, INTENT (IN) :: acof(nslibd,0:lmd,natd)
94 : COMPLEX, INTENT (IN) :: bcof(nslibd,0:lmd,natd)
95 : integer, intent (in) :: k1d,k2d,k3d,ig2(n3d)
96 : integer, intent (in) :: ig(-k1d:k1d,-k2d:k2d,-k3d:k3d)
97 : real, intent (in) :: sk2(n2d),phi2(n2d)
98 :
99 : ! ..
100 : ! .. Local Scalars ..
101 : integer n2,k,j,l
102 : REAL :: s
103 : COMPLEX :: xdnout
104 : INTEGER :: i,i1,i2,i3,ii3,ix,iy,iz,na,nplo,nbn
105 : INTEGER :: nbmin,nbmax,n
106 : INTEGER :: nplot,nq,nt,jm,iter,ii1,ii2
107 : CHARACTER*8 :: dop,iop
108 : LOGICAL :: twodim
109 : ! ..
110 : ! .. Local Arrays ..
111 : REAL :: pt(3),vec1(3),vec2(3),vec3(3),zero(3)
112 : INTEGER :: grid(3)
113 : LOGICAL :: cartesian,xsf
114 : REAL :: rhocc(jmtd)
115 : REAL :: point(3),post(3,natd)
116 : CHARACTER(len=30):: filename
117 : CHARACTER(len=20):: name1,name2,name3
118 : CHARACTER(len=10):: vandername
119 : CHARACTER*8 :: name(10)
120 : c ..basis wavefunctions in the vacuum
121 0 : complex, allocatable :: ac(:,:,:),bc(:,:,:)
122 : complex, allocatable :: ac_1(:,:,:),bc_1(:,:,:)
123 0 : real, allocatable :: u(:,:,:),ue(:,:,:)
124 : real, allocatable :: u_1(:,:,:),ue_1(:,:,:)
125 :
126 : integer ik,nv2,ivac,jvac
127 : complex factor
128 : real fas
129 :
130 :
131 : NAMELIST /plot/twodim,cartesian,vec1,vec2,vec3,grid,zero,filename
132 :
133 : intrinsic real,aimag,conjg
134 : c external dotirp
135 : c real dotirp !module now
136 :
137 0 : call timestart("wann_plot")
138 0 : IF (slice) THEN
139 0 : nbmin = nnne
140 0 : nbmax = nnne
141 : ELSE
142 : nbmin = 1
143 : nbmax = nslibd
144 : ENDIF
145 :
146 : !write(*,*)'in wann_plot'
147 : !write(*,*)'jspin:',jspin
148 : !write(*,*)'addnoco:',addnoco
149 : !write(*,*)'l_noco,l_ss: ',l_noco,l_ss
150 : !write(*,*)'qss:',qss
151 :
152 0 : INQUIRE(file ="plot_inp",exist= twodim)
153 0 : IF (.NOT.twodim) THEN !no input file exists, create a template and
154 : !exit
155 0 : OPEN(20,file ="plot_inp")
156 0 : WRITE(20,'(i2,a5,l1)') 1,",xsf=",.false.
157 : c WRITE(20,*) "&PLOT twodim=t,cartesian=t"
158 : c WRITE(20,*) " vec1(1)=10.0 vec2(2)=10.0"
159 : c WRITE(20,*) " filename='plot1' /"
160 0 : WRITE(20,*) "&PLOT twodim=f,cartesian=f"
161 0 : WRITE(20,*) " vec1(1)=1.0 vec1(2)=0.0 vec1(3)=0.0 "
162 0 : WRITE(20,*) " vec2(1)=0.0 vec2(2)=1.0 vec2(3)=0.0 "
163 0 : WRITE(20,*) " vec3(1)=0.0 vec3(2)=0.0 vec3(3)=1.0 "
164 0 : WRITE(20,*) " grid(1)=30 grid(2)=30 grid(3)=30 "
165 0 : WRITE(20,*) " zero(1)=0.0 zero(2)=0.0 zero(3)=0.0 "
166 0 : WRITE(20,*) " filename ='plot2' /"
167 0 : CLOSE(20)
168 0 : WRITE(*,*) "No plot_inp file found. Created a template"
169 : CALL juDFT_error("Missing input for plot; modify plot_inp"
170 0 : + ,calledby ="wann_plot")
171 : ENDIF
172 :
173 : c make preparations for plotting in vacuum
174 :
175 0 : if(film)then
176 : allocate ( ac(nv2d,nslibd,2),bc(nv2d,nslibd,2),
177 0 : + u(nmzd,nv2d,nvac),ue(nmzd,nv2d,nvac))
178 : call wann_2dvacabcof(
179 : > nv2d,nslibd,nvac,nmzd,nmz,omtil,vz,nv,bkpt,z1,
180 : > nvd,k1,k2,k3,evac,bbmat,delz,bmat,nbasfcn,neigd,zMat,
181 0 : < ac,bc,u,ue,addnoco,l_ss,qss,jspin)
182 : endif
183 :
184 :
185 :
186 :
187 : !<-- Open the plot_inp file for input
188 0 : OPEN (18,file='plot_inp')
189 0 : READ(18,'(i2,5x,l1)') nplot,xsf
190 : ! If xsf is specified we create an input file for xcrysden
191 0 : IF (nplot.ge.2)
192 : & CALL juDFT_error
193 : + ("plots one by one, please, this is not charge density"
194 0 : + ,calledby="wann_plot")
195 : !<-- Loop over all plots
196 0 : DO nplo=1,nplot
197 : ! the defaults
198 0 : twodim = .TRUE.;cartesian=.TRUE.;grid=(/100,100,100/)
199 0 : vec1 = (/0.,0.,0./);vec2=(/0.,0.,0./);vec3=(/0.,0.,0./)
200 0 : zero = (/0.,0.,0./);filename="default"
201 0 : READ(18,plot)
202 0 : IF (twodim.AND.ANY(grid(1:2)<1))
203 : + CALL juDFT_error("Illegal grid size in plot",calledby
204 0 : + ="wann_plot")
205 0 : IF (.NOT.twodim.AND.ANY(grid<1))
206 : + CALL juDFT_error("Illegal grid size in plot",calledby
207 0 : + ="wann_plot")
208 0 : IF (twodim) grid(3) = 1
209 : !calculate cartesian coordinates if needed
210 0 : IF (.NOT.cartesian) THEN
211 0 : vec1=matmul(amat,vec1)
212 0 : vec2=matmul(amat,vec2)
213 0 : vec3=matmul(amat,vec3)
214 0 : zero=matmul(amat,zero)
215 : ENDIF
216 : !Open the file
217 0 : IF (filename =="default") WRITE(filename,'(a,i2)') "plot",nplo
218 : c..loop by the bands
219 0 : bands:DO nbn = nbmin,nbmax
220 :
221 0 : IF (xsf) THEN
222 0 : write (name1,22) ikpt,nbn,jspin
223 : 22 format (i5.5,'.',i3.3,'.real.',i1,'.xsf')
224 0 : write (name2,23) ikpt,nbn,jspin
225 : 23 format (i5.5,'.',i3.3,'.imag.',i1,'.xsf')
226 0 : write (name3,24) ikpt,nbn,jspin
227 : 24 format (i5.5,'.',i3.3,'.absv.',i1,'.xsf')
228 0 : OPEN(55,file=name1)
229 0 : CALL xsf_WRITE_atoms(55,atoms,film,amat)
230 0 : OPEN(56,file=name2)
231 0 : CALL xsf_WRITE_atoms(56,atoms,film,amat)
232 0 : OPEN(57,file=name3)
233 0 : CALL xsf_WRITE_atoms(57,atoms,film,amat)
234 : CALL xsf_WRITE_header(55,twodim,filename,(vec1),
235 : & (vec2),(vec3),zero
236 0 : $ ,grid)
237 : CALL xsf_WRITE_header(56,twodim,filename,(vec1),
238 : & (vec2),(vec3),zero
239 0 : $ ,grid)
240 : CALL xsf_WRITE_header(57,twodim,filename,(vec1),
241 : & (vec2),(vec3),zero
242 0 : $ ,grid)
243 : ELSE
244 0 : IF (nbn.EQ.nbmin) THEN
245 0 : WRITE (vandername,201) ikpt,jspin
246 0 : IF(l_noco) WRITE(vandername,202)ikpt,jspin
247 0 : IF(l_sgwf) WRITE(vandername,202)index_kq,jspin
248 : 201 FORMAT ('UNK',i5.5,'.',i1)
249 : 202 FORMAT ('RNK',i5.5,'.',i1)
250 0 : OPEN(55,file=vandername)
251 0 : IF(.NOT.l_sgwf) THEN
252 0 : WRITE(55,7) grid(1),grid(2),grid(3),ikpt,nslibd
253 : ELSE
254 0 : WRITE(55,7) grid(1),grid(2),grid(3),index_kq,nslibd
255 : ENDIF
256 : 7 FORMAT (5i4)
257 : ENDIF
258 : ENDIF
259 :
260 0 : if(film)then
261 0 : fas=-bkpt(3)*bmat(3,3)*z1
262 0 : factor=cmplx(cos(fas),sin(fas))
263 : else
264 : factor=cmplx(1.0,0.0)
265 : endif
266 :
267 0 : DO iz = 0,grid(3)-1
268 0 : DO iy = 0,grid(2)-1
269 0 : xloop:DO ix = 0,grid(1)-1
270 : point = zero+vec1*REAL(ix)/grid(1)+vec2*REAL(iy)
271 0 : $ /grid(2)
272 0 : IF (.NOT.twodim) point = point+vec3*REAL(iz)/grid(3)
273 : !Check if the point is in MT-sphere
274 :
275 0 : ii1 = 3
276 0 : ii2 = 3
277 0 : ii3 = 3
278 0 : IF (film ) ii3 = 0
279 :
280 0 : DO i1 = -ii1,ii1
281 0 : DO i2 = -ii2,ii2
282 0 : DO i3 = -ii3,ii3
283 0 : pt = point+MATMUL(amat,(/i1,i2,i3/))
284 0 : na = 0
285 0 : DO nt = 1,ntype
286 0 : DO nq = 1,neq(nt)
287 0 : na = na + 1
288 0 : s = SQRT(dot_PRODUCT(pos(:,na)-pt,pos(:,na)-pt))
289 0 : IF (s<rmsh(jri(nt),nt)) THEN
290 : CALL wann_real(
291 : > pt,nt,na,0,1,bkpt,nbn,
292 : > n3d,nmzxyd,n2d,ntypsd,lmaxd,jmtd,
293 : > natd,ntypd,nmzd,nop,nop2,mrot,tau,invtab,
294 : > nq3,nvac,invs,z1,delz,nmz,nmzxy,nq2,
295 : > lmax,rmsh,jri,pos,ngopr,ntypsy,nvd,
296 : > omtil,amat,bmat,nlod,llod,nlo,llo,
297 : > ff,gg,flo,acof(nbn,:,:),bcof(nbn,:,:),
298 : > ccof(:,nbn,:,:),zMat,
299 : > nv,k1,k2,k3,lmd,nbasfcn,l_ss,qss,jspin,addnoco,
300 0 : < xdnout)
301 0 : xdnout=xdnout*factor
302 0 : IF (xsf) THEN
303 0 : WRITE(55,*) real(xdnout)
304 0 : WRITE(56,*) aimag(xdnout)
305 0 : WRITE(57,*) real(xdnout*conjg(xdnout))
306 : ELSE
307 0 : WRITE(55,8) real(xdnout),aimag(xdnout)
308 : ENDIF
309 0 : CYCLE xloop
310 : ENDIF
311 : ENDDO
312 : ENDDO !nt
313 : ENDDO
314 : ENDDO
315 : ENDDO !i1
316 : !Check for point in vacuum
317 0 : IF (film.AND.ABS(point(3))>=z1) THEN
318 0 : ivac=1
319 0 : if (point(3).lt. 0.0)ivac=2
320 0 : jvac=ivac
321 0 : if(nvac==1)jvac=1
322 : call wann_plot_vac(point,z1,nmzd,nv2d,n3d,nvac,
323 : > nmz,delz,bmat,bbmat,evac,bkpt,vz,jspin,k1,k2,k3,nvd,
324 : > nbasfcn,neigd,nv,omtil,nslibd,ac(:,nbn,ivac),
325 : & bc(:,nbn,ivac),
326 0 : & u(:,:,jvac),ue(:,:,jvac),xdnout)
327 0 : xdnout=xdnout*factor
328 : if(real(xdnout).gt.9.0 .or.real(xdnout).lt.-9.0
329 0 : & .or.aimag(xdnout).gt.9.0 .or. aimag(xdnout).lt.-9.0)then
330 0 : xdnout=cmplx(0.0,0.0)
331 0 : print*,"vac-problem at z=",point(3)
332 : endif
333 : c CALL wann_real(
334 : c > point,0,0,1,0,bkpt,
335 : c > n3d,nmzxyd,n2d,ntypsd,lmaxd,jmtd,
336 : c > natd,ntypd,nmzd,nop,nop2,mrot,tau,invtab,
337 : c > nq3,nvac,invs,z1,delz,nmz,nmzxy,nq2,
338 : c > lmax,rmsh,jri,pos,ngopr,ntypsy,nvd,
339 : c > omtil,amat,bmat,nlod,llod,nlo,llo,
340 : c > ff,gg,flo,acof(nbn,:,:),bcof(nbn,:,:),
341 : c > ccof(:,nbn,:,:),z(:,nbn),
342 : c > nv,k1,k2,k3,lmd,nbasfcn,
343 : c < xdnout)
344 0 : IF (xsf) THEN
345 0 : WRITE(55,*) real(xdnout)
346 0 : WRITE(56,*) aimag(xdnout)
347 0 : WRITE(57,*) real(xdnout*conjg(xdnout))
348 : ELSE
349 0 : WRITE(55,8) real(xdnout),aimag(xdnout)
350 : ENDIF
351 : CYCLE xloop
352 : END IF
353 :
354 :
355 : CALL wann_real(
356 : > point,0,0,0,2,bkpt,nbn,
357 : > n3d,nmzxyd,n2d,ntypsd,lmaxd,jmtd,
358 : > natd,ntypd,nmzd,nop,nop2,mrot,tau,invtab,
359 : > nq3,nvac,invs,z1,delz,nmz,nmzxy,nq2,
360 : > lmax,rmsh,jri,pos,ngopr,ntypsy,nvd,
361 : > omtil,amat,bmat,nlod,llod,nlo,llo,
362 : > ff,gg,flo,acof(nbn,:,:),bcof(nbn,:,:),
363 : > ccof(:,nbn,:,:),zMat,
364 : > nv,k1,k2,k3,lmd,nbasfcn,l_ss,qss,jspin,addnoco,
365 0 : < xdnout)
366 0 : xdnout=xdnout*factor
367 0 : IF (xsf) THEN
368 0 : WRITE(55,*) real(xdnout)
369 0 : WRITE(56,*) aimag(xdnout)
370 0 : WRITE(57,*) real(xdnout*conjg(xdnout))
371 : ELSE
372 0 : WRITE(55,8) real(xdnout),aimag(xdnout)
373 : ENDIF
374 : ENDDO xloop
375 : ENDDO
376 : ENDDO !z-loop
377 :
378 0 : IF (xsf) THEN
379 0 : CALL xsf_WRITE_endblock(55,twodim)
380 0 : CALL xsf_WRITE_endblock(56,twodim)
381 0 : CALL xsf_WRITE_endblock(57,twodim)
382 0 : CLOSE (55) ; CLOSE (56) ; CLOSE (57)
383 : ENDIF
384 : c..end of the loop by the bands
385 : ENDDO bands
386 : ENDDO !nplot
387 0 : CLOSE(18)
388 0 : IF (.not.xsf) CLOSE(55)
389 : 8 FORMAT (f20.12,1x,f20.12)
390 :
391 :
392 0 : call timestop("wann_plot")
393 0 : RETURN
394 0 : END SUBROUTINE wann_plot
395 : !------------------------------------------
396 :
397 : END MODULE m_wann_plot
|
