Line data Source code
1 : c*************************c
2 : c routine to set up the c
3 : c composite matrices as c
4 : c input for wannier90 c
5 : c*************************c
6 : module m_wann_gwf_commat
7 : USE m_juDFT
8 : implicit none
9 : contains
10 :
11 0 : subroutine wann_gwf_commat(nkpts,nntot,bpt,nqpts,nntot_q,
12 0 : > bpt_q,gb,gb_q,latt_const_q,
13 : > l_unformatted,l_amn,l_mmn,l_dim,
14 0 : > nparampts,param_vec)
15 : use m_wann_gwf_tools
16 : use m_wann_gwf_auxovlp
17 : use m_constants, only : pimach
18 : use m_intgr, only : intgz0
19 :
20 : implicit none
21 :
22 : ! input parameters
23 : logical,intent(in) :: l_unformatted,l_amn,l_mmn,l_dim(3)
24 : integer,intent(in) :: nkpts,nntot,nqpts,nntot_q,nparampts
25 : integer,intent(in) :: bpt(nntot,nkpts),bpt_q(nntot_q,nqpts)
26 : integer,intent(in) :: gb(3,nntot,nkpts),gb_q(3,nntot_q,nqpts)
27 : real, intent(in) :: latt_const_q,param_vec(3,nparampts)
28 :
29 : ! allocatable arrays
30 0 : integer,allocatable :: gb_kq(:,:,:)
31 0 : integer,allocatable :: bpt_kq(:,:)
32 0 : integer,allocatable :: g1(:,:,:),g2(:,:)
33 0 : integer,allocatable :: b1(:,:),b2(:)
34 0 : real, allocatable :: xxr(:),xxi(:)
35 0 : complex,allocatable :: tmp_a(:,:,:),tmp_m(:,:,:,:)
36 0 : complex,allocatable :: mmnk(:,:,:,:),mmnq(:,:,:,:)
37 0 : complex,allocatable :: mk(:,:,:),mq(:,:,:)
38 0 : complex,allocatable :: mmn(:,:,:,:)
39 0 : complex,allocatable :: amn_arti2(:)
40 0 : character(len=20),allocatable :: feig(:),famn(:)
41 0 : character(len=20),allocatable :: fmmn(:),fmmn2(:)
42 :
43 : ! further variables
44 : integer :: nbnd,nwfs,kb,qb,nn_kq,kqb,arr_len,shift(3)
45 : integer :: i,j,ikqpt,ikqpt_b,nk,nq,g(3)
46 : integer :: nbnd_arti,nqpts_arti,nntot_arti
47 : integer :: nkqpts,nntot_kq,nn,ikqpt_help
48 : integer :: dump,nn_arti,nkp,k,b,q,nwf,n,is,ie,nx
49 : real :: mmn_r,mmn_i,tau,tpi,a_arti,b_arti,eig
50 : real :: dx,x0,sig,center,amnr,amni,norm,bf
51 : complex :: mmn_arti,amn_arti,amn,psi,wftrial
52 : logical :: l_exist,l_miss,l_fe,l_fa,l_fm,l_fm2,l_proj
53 : character(len=20) :: fq
54 : character(len=32) :: fmt,fmt2
55 :
56 0 : call timestart("wann_gwf_commat")
57 0 : write(*,*)'create HDWFs input for w90...'
58 :
59 : ! arr_len=3; shifty=0; shiftz=0
60 : ! if(l_dim(1))then
61 : ! arr_len=arr_len+1
62 : ! shifty=shifty+1
63 : ! shiftz=shiftz+1
64 : ! endif
65 : ! if(l_dim(2))then
66 : ! arr_len=arr_len+1
67 : ! shiftz=shiftz+1
68 : ! endif
69 : ! if(l_dim(3))arr_len=arr_len+1
70 :
71 0 : call get_dimension(l_dim,arr_len)
72 0 : call get_shift(l_dim,shift)
73 :
74 0 : write(*,*)'dimension:',arr_len
75 0 : write(*,*)'x?,y?,z? :',l_dim(1:3)
76 0 : if(arr_len.le.3) call juDFT_error("dimension<4",
77 0 : > calledby='wann_gwf_commat')
78 :
79 0 : nkqpts=nkpts*nqpts
80 0 : tpi = 2.0*pimach()
81 0 : a_arti = latt_const_q
82 0 : b_arti = 0.98*latt_const_q
83 : ! b_arti = 0.5*latt_const_q
84 : ! b_arti = 0.7*latt_const_q
85 : ! b_arti = 0.3*latt_const_q
86 : ! b_arti = 0.05*latt_const_q
87 0 : l_fe = .true.
88 0 : l_fa = .true.
89 0 : l_fm = .true.
90 0 : l_fm2 = .true.
91 0 : l_exist= .true.
92 :
93 : ! compute auxiliary overlaps and projections
94 0 : tau = 1.0*tpi/real(nqpts)/a_arti
95 0 : mmn_arti = 2.0*tpi*tpi*sin(tau*b_arti/2.0)
96 0 : mmn_arti = mmn_arti/(tpi*tpi-tau*tau*b_arti*b_arti)
97 0 : mmn_arti = mmn_arti/(tau*b_arti)
98 : c mmn_arti = 0.5
99 0 : write(*,*)'ov(0) =',mmn_arti
100 : c CALL wann_gwf_auxovlp(param_vec(:,1),param_vec(:,2),
101 : c > latt_const_q,mmn_arti)
102 :
103 0 : amn_arti = cmplx(1.0,0.0)
104 :
105 0 : open(777,file='ma_aux')
106 0 : write(777,'(2(f24.18))')real(mmn_arti),aimag(mmn_arti)
107 0 : write(777,'(2(f24.18))')real(amn_arti),aimag(amn_arti)
108 0 : write(777,*)a_arti
109 0 : write(777,*)b_arti
110 0 : close(777)
111 :
112 0 : bf = b_arti/a_arti
113 0 : sig = 1.00!1.2!0.05
114 0 : center = 0.0
115 : !norm = sqrt(sqrt(2./pimach())/sig)*sqrt(2./bf)
116 0 : nx=301*nqpts+1
117 0 : dx = 1./real(nx-1)*nqpts
118 0 : allocate(xxr(nx),xxi(nx))
119 0 : allocate(amn_arti2(nqpts))
120 0 : open(555,file='debug_wfaux')
121 0 : do q=1,nqpts
122 0 : do i=1,nx
123 : !x0 = -bf/2. + (i-1)*dx*bf
124 : !xx(i) = exp(-x0*x0/sig/sig)*cos(pimach()/bf*x0)
125 0 : x0 = -nqpts/2. + (i-1)*dx
126 0 : psi=wfpsi(x0,param_vec(3,q),bf)
127 0 : wftrial = trialorb(x0,center,sig)
128 0 : xxr(i) = real(conjg(wftrial)*psi)
129 0 : xxi(i) = aimag(conjg(wftrial)*psi)
130 0 : write(555,'(6(f10.6,1x))')param_vec(3,q),x0,
131 0 : > real(psi),aimag(psi),real(wftrial),aimag(wftrial)
132 : enddo
133 : !call intgz0(xx,dx,nx,amnr,.false.)
134 : !amn_arti2(q) = amnr*norm
135 0 : call intgz0(xxr,dx,nx,amnr,.false.)
136 0 : call intgz0(xxi,dx,nx,amni,.false.)
137 0 : amn_arti2(q) = cmplx(amnr,amni)
138 : enddo
139 0 : close(555)
140 :
141 0 : open(777,file='debug_amn')
142 0 : do q=1,nqpts
143 0 : write(777,*)q,amn_arti2(q)
144 : enddo
145 0 : close(777)
146 0 : deallocate(xxr,xxi)
147 :
148 : ! are all necessary files present?
149 0 : allocate(feig(nqpts),famn(nqpts),fmmn(nqpts),fmmn2(nqpts))
150 0 : l_miss=.false.
151 0 : do q=1,nqpts
152 0 : write(fq,'(i4.4)')q
153 0 : feig(q) = 'WF1_'//trim(fq)//'.eig' ! energies
154 0 : famn(q) = 'WF1_'//trim(fq)//'.amn' ! projections
155 0 : fmmn(q) = 'WF1_'//trim(fq)//'.mmn' ! k-overlaps
156 0 : fmmn2(q)= 'param_'//trim(fq)//'.mmn' ! q-overlaps
157 :
158 0 : inquire(file=feig(q) ,exist=l_fe )
159 0 : if(l_amn) inquire(file=famn(q) ,exist=l_fa )
160 0 : if(l_mmn) inquire(file=fmmn(q) ,exist=l_fm )
161 0 : if(l_mmn) inquire(file=fmmn2(q),exist=l_fm2)
162 :
163 0 : if(.not.l_fe ) write(*,*)'missing: ',feig(q)
164 0 : if(.not.l_fa ) write(*,*)'missing: ',famn(q)
165 0 : if(.not.l_fm ) write(*,*)'missing: ',fmmn(q)
166 0 : if(.not.l_fm2) write(*,*)'missing: ',fmmn2(q)
167 0 : if(.not.(l_fe.and.l_fa.and.l_fm.and.l_fm2))l_miss=.true.
168 : enddo
169 0 : if(l_mmn) inquire(file='bkqpts',exist=l_exist)
170 0 : if(.not.l_exist) write(*,*)'missing: bkqpts'
171 0 : inquire(file='proj',exist=l_proj)
172 0 : if(.not.l_proj) write(*,*)'missing: proj'
173 0 : if(l_miss.or.(.not.l_exist).or.(.not.l_proj))
174 0 : > call juDFT_error("missing file(s) for HDWFs")
175 :
176 : ! get number of bands and wfs from proj
177 0 : open(405,file='proj',status='old')
178 0 : read(405,*)nwfs,nbnd
179 0 : close(405)
180 0 : write(*,*)'nbnd=',nbnd
181 0 : write(*,*)'nwfs=',nwfs
182 :
183 :
184 : c*****************************c
185 : c COMPOSITE .EIG FILE c
186 : c*****************************c
187 0 : open(305,file='WF1_gwf.eig')
188 0 : do q=1,nqpts
189 0 : open(405,file=feig(q))
190 0 : do k=1,nkpts
191 0 : ikqpt=get_index_kq(k,q,nkpts)
192 0 : do i=1,nbnd
193 0 : read(405,*)b,nk,eig
194 0 : write(305,'(2i12,f19.13)')b,ikqpt,eig
195 : enddo
196 : enddo
197 0 : close(405)!,status='delete')
198 : enddo
199 0 : close(305)
200 :
201 :
202 : c*****************************c
203 : c COMPOSITE .AMN FILE c
204 : c*****************************c
205 0 : if(.not.l_amn) goto 100 ! skip amn part
206 :
207 0 : if(l_unformatted) then
208 0 : allocate(tmp_a(nbnd,nwfs,nkqpts))
209 0 : do q=1,nqpts
210 0 : is = get_index_kq(1,q,nkpts)
211 0 : ie = get_index_kq(nkpts,q,nkpts)
212 0 : open(405,file=famn(q),form='unformatted')
213 0 : read(405)b,k,nwf
214 0 : read(405)tmp_a(:,:,is:ie)
215 0 : do i=1,nwfs
216 0 : tmp_a(:,i,is:ie) = tmp_a(:,i,is:ie)*amn_arti2(q)
217 : enddo
218 0 : close(405)!,status='delete')
219 : enddo
220 :
221 0 : open(306,file='WF1_gwf.amn',form='unformatted')
222 0 : write(306)nbnd,nkqpts,nwfs
223 0 : write(306)tmp_a!*amn_arti
224 0 : close(306)
225 0 : deallocate(tmp_a)
226 : else
227 0 : open(306,file='WF1_gwf.amn')
228 0 : write(306,*)'Projections for HDWFs'
229 0 : write(306,'(i5,i7,i5)')nbnd,nkqpts,nwfs
230 :
231 0 : do q=1,nqpts
232 0 : open(405,file=famn(q))
233 0 : read(405,*)!title
234 0 : read(405,*)b,k,nwf
235 0 : do k=1,nkpts
236 0 : ikqpt=get_index_kq(k,q,nkpts)
237 0 : do nwf=1,nwfs
238 0 : do i=1,nbnd
239 0 : read(405,*)b,n,nk,mmn_r,mmn_i
240 0 : amn = cmplx(mmn_r,mmn_i)
241 0 : amn = amn*amn_arti2(q) !*amn_arti
242 : write(306,'(i5,i5,i7,3x,2f18.12)')
243 0 : > b,n,ikqpt,real(amn),aimag(amn)
244 : enddo
245 : enddo
246 : enddo
247 0 : close(405,status='delete')
248 : enddo
249 0 : close(306)
250 : endif!l_unformatted
251 :
252 : 100 continue
253 :
254 :
255 : c*****************************c
256 : c COMPOSITE .MMN FILE c
257 : c*****************************c
258 0 : if(.not.l_mmn) goto 200 ! skip mmn part
259 :
260 0 : write(fmt,'(a,i1,a)')'(2i6,3x,',arr_len,'i4)'
261 0 : write(fmt2,'(a,i1,a)')'(i7,i7,3x,',arr_len,'i4)'
262 :
263 0 : open (202,file='bkqpts',form='formatted',status='old')
264 0 : rewind (202)
265 0 : read (202,'(i4)') nntot_kq
266 0 : write (*,*) 'nntot_kq=',nntot_kq
267 0 : allocate ( gb_kq(arr_len,nntot_kq,nkqpts),
268 0 : & bpt_kq(nntot_kq,nkqpts))
269 0 : do ikqpt=1,nkqpts
270 0 : do nn=1,nntot_kq
271 : read (202,fmt)
272 0 : & ikqpt_help,bpt_kq(nn,ikqpt),(gb_kq(i,nn,ikqpt),i=1,arr_len)
273 : enddo
274 : enddo
275 0 : close (202)
276 :
277 0 : if(l_unformatted) then
278 0 : allocate(mmnk(nbnd,nbnd,nntot,nkpts))
279 0 : allocate(mmnq(nbnd,nbnd,nntot_q,nkpts))
280 0 : allocate(g1(3,nntot,nkpts),g2(3,nntot_q))
281 0 : allocate(b1(nntot,nkpts),b2(nntot_q))
282 0 : allocate(mmn(nbnd,nbnd,nntot_kq,nkqpts))
283 :
284 0 : do q=1,nqpts
285 0 : is = get_index_kq(1,q,nkpts)
286 0 : ie = get_index_kq(nkpts,q,nkpts)
287 :
288 0 : open(405,file=fmmn(q),form='unformatted')
289 0 : read(405)b,k,nn
290 0 : read(405)b1,g1
291 0 : read(405)mmnk
292 0 : close(405,status='delete')
293 :
294 0 : open(405,file=fmmn2(q),form='unformatted')
295 0 : read(405)b,k,nn
296 0 : read(405)b2,g2
297 0 : read(405)mmnq
298 0 : close(405,status='delete')
299 0 : mmnq = mmnq*conjg(mmn_arti)
300 :
301 0 : do k=1,nkpts
302 0 : ikqpt=get_index_kq(k,q,nkpts)
303 0 : do nn=1,nntot_kq
304 0 : kqb = bpt_kq(nn,ikqpt)
305 0 : kb=get_index_k(kqb,nkpts)
306 0 : qb=get_index_q(kqb,nkpts)
307 0 : if(q.eq.qb) then ! k-overlap
308 : nk=get_index_nn_k(bpt(:,k),nntot,kb,
309 : > gb_kq(:,nn,ikqpt),
310 0 : > gb(1:3,1:nntot,k),arr_len)
311 0 : mmn(:,:,nn,ikqpt)=mmnk(:,:,nk,k)
312 : !write(*,*)'k neig',k,kb,nk
313 0 : elseif(k.eq.kb) then ! q-overlap
314 : nq=get_index_nn_q(bpt_q(:,q),nntot_q,qb,
315 : > gb_kq(:,nn,ikqpt),
316 : > gb_q(1:3,1:nntot_q,q),
317 0 : > arr_len,shift,l_dim)
318 0 : mmn(:,:,nn,ikqpt)=mmnq(:,:,nq,k)
319 : !write(*,*)'q neig',q,qb,nq
320 : else ! otherwise problem
321 : call juDFT_error("problem in overlap mmn gwf",
322 0 : > calledby="wann_gwf_commat")
323 : endif
324 : enddo!nn
325 : enddo!k
326 : enddo!q
327 0 : deallocate(b1,b2,g1,g2,mmnk,mmnq)
328 :
329 0 : open(307,file='WF1_gwf.mmn',form='unformatted')
330 0 : write(307)nbnd,nkqpts,nntot_kq
331 0 : write(307)bpt_kq,gb_kq
332 0 : write(307)mmn
333 0 : close(307)
334 0 : deallocate(mmn)
335 : else
336 0 : allocate(mk(nbnd,nbnd,nntot))
337 0 : allocate(mq(nbnd,nbnd,nntot_q))
338 0 : open(307,file='WF1_gwf.mmn')
339 0 : write(307,*)'Overlaps for HDWFs'
340 0 : write(307,'(i5,i7,i5)')nbnd,nkqpts,nntot_kq
341 :
342 0 : do q=1,nqpts
343 0 : open(405,file=fmmn(q))
344 0 : read(405,*)!title
345 0 : read(405,*)b,k,nn
346 0 : open(406,file=fmmn2(q))
347 0 : read(406,*)!title
348 0 : read(406,*)b,k,nn
349 0 : do k=1,nkpts
350 0 : ikqpt=get_index_kq(k,q,nkpts)
351 :
352 0 : do nn=1,nntot ! read k-overlaps
353 0 : kb = bpt(nn,k)
354 0 : read(405,*)nkp,b,g(1:3)
355 0 : do i=1,nbnd
356 0 : do j=1,nbnd
357 0 : read(405,*)mmn_r,mmn_i
358 0 : mk(j,i,nn)=cmplx(mmn_r,mmn_i)
359 : enddo
360 : enddo
361 : enddo
362 :
363 0 : do nn=1,nntot_q ! read q-overlaps
364 0 : qb = bpt_q(nn,q)
365 0 : read(406,*)nkp,b,g(1:3)
366 0 : do i=1,nbnd
367 0 : do j=1,nbnd
368 0 : read(406,*)mmn_r,mmn_i
369 0 : mq(j,i,nn)=cmplx(mmn_r,mmn_i)
370 : enddo
371 : enddo
372 : enddo
373 0 : mq = mq*conjg(mmn_arti)
374 :
375 0 : do nn=1,nntot_kq ! write composite overlaps
376 0 : kqb = bpt_kq(nn,ikqpt)
377 0 : kb=get_index_k(kqb,nkpts)
378 0 : qb=get_index_q(kqb,nkpts)
379 0 : write(307,fmt2)ikqpt,kqb,gb_kq(1:arr_len,nn,ikqpt)
380 0 : if(q.eq.qb) then ! k-overlap
381 : nk=get_index_nn_k(bpt(:,k),nntot,kb,
382 : > gb_kq(:,nn,ikqpt),
383 0 : > gb(1:3,1:nntot,k),arr_len)
384 0 : do i=1,nbnd
385 0 : do j=1,nbnd
386 : write(307,'(2f24.18)')
387 0 : > real(mk(j,i,nk)),aimag(mk(j,i,nk))
388 : enddo
389 : enddo
390 0 : elseif(k.eq.kb) then ! q-overlap
391 : nq=get_index_nn_q(bpt_q(:,q),nntot_q,qb,
392 : > gb_kq(:,nn,ikqpt),
393 : > gb_q(1:3,1:nntot_q,q),
394 0 : > arr_len,shift,l_dim)
395 0 : do i=1,nbnd
396 0 : do j=1,nbnd
397 : write(307,'(2f24.18)')
398 0 : > real(mq(j,i,nq)),aimag(mq(j,i,nq))
399 : enddo
400 : enddo
401 : else ! otherwise problem
402 : call juDFT_error("problem in overlap mmn gwf",
403 0 : > calledby="wann_gwf_commat")
404 : endif
405 : enddo!nn
406 : enddo!k
407 0 : close(406,status='delete')
408 0 : close(405,status='delete')
409 : enddo!q
410 0 : close(307)
411 0 : deallocate(mk,mq)
412 : endif!l_unformatted
413 :
414 0 : deallocate(gb_kq,bpt_kq)
415 : 200 continue
416 0 : deallocate(feig,famn,fmmn,fmmn2)
417 0 : deallocate(amn_arti2)
418 :
419 0 : call timestop("wann_gwf_commat")
420 : contains
421 :
422 :
423 0 : complex function wfu(x0,qpt,bb)
424 : implicit none
425 : real,intent(in) :: x0,qpt,bb
426 : integer :: i0
427 : real :: qx,xshift
428 : complex :: pha
429 :
430 0 : wfu = cmplx(0.,0.)
431 0 : i0=floor(x0+0.5)
432 0 : xshift=x0-real(i0)
433 0 : if(abs(xshift) .ge. bb/2.) return
434 :
435 0 : qx = 2.*pimach()*qpt*xshift
436 0 : pha = cmplx( cos(qx), -sin(qx) )
437 0 : wfu = pha*sqrt(2./bb)*cos(pimach()*xshift/bb)
438 0 : end function wfu
439 :
440 0 : complex function wfpsi(x0,qpt,bb)
441 : implicit none
442 : real,intent(in) :: x0,qpt,bb
443 : complex :: pha
444 : real :: qx
445 :
446 0 : qx = 2.*pimach()*qpt*x0
447 0 : pha = cmplx( cos(qx), sin(qx) )
448 0 : wfpsi = pha*wfu(x0,qpt,bb)
449 0 : end function wfpsi
450 :
451 :
452 0 : real function trialorb(x0,center,sigma)
453 : implicit none
454 : real,intent(in) :: x0,center,sigma
455 : real :: xp
456 :
457 0 : xp = (x0-center)/sigma
458 0 : trialorb = exp(-xp*xp)*sqrt(sqrt(2./pimach())/sigma)
459 0 : end function trialorb
460 :
461 : end subroutine wann_gwf_commat
462 : end module m_wann_gwf_commat
|
