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