00001 PROGRAM WANNIER
00002 use m_rdinput
00003 use fft_3d
00004 use m_bvector, only: generate_bvectors
00005 use m_frmsf_wan, only: wrt_frmsf_wan
00006 use m_dos, only: calculate_dos
00007 use m_eigenstate, only: calculate_eigenstate
00008 use m_wrt_frmsf, only: wrt_frmsf
00009 use m_wrt_model
00010 use m_dmx, only: calculate_dmx
00011 use m_fat_band, only: calc_fat_band
00012 include "config.h"
00013
00014 call read_input
00015 NGAUSS=n_occ
00016
00017
00018
00019
00020
00021
00022
00023
00024 if(flg_specify_energywindow)then
00025 E_LOWER=E_LOWER/au
00026 E_UPPER=E_UPPER/au
00027 write(6,*)'E_LOWER in au',E_LOWER
00028 write(6,*)'E_UPPER in au',E_UPPER
00029 endif
00030
00031 E_LOWER_inner=E_LOWER_inner/au
00032 E_UPPER_inner=E_UPPER_inner/au
00033 write(6,*)'E_LOWER_inner in au',E_LOWER_inner
00034 write(6,*)'E_UPPER_inner in au',E_UPPER_inner
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044 WRITE(6,*)
00045 WRITE(6,*)'=========================================='
00046 WRITE(6,*)'=== WANNIER FUNCTION CALCULATION START ==='
00047 WRITE(6,*)'=========================================='
00048 WRITE(6,*)
00049
00050
00051 OPEN(117,FILE='./dir-wfn/dat.bandcalc')
00052 rewind(117)
00053 read(117,*) Ecut_for_psi
00054 read(117,*) FermiEnergy
00055 read(117,*) Etot
00056 write(6,*)'Ecut_for_psi=',Ecut_for_psi
00057 write(6,*)'FermiEnergy=',FermiEnergy
00058 write(6,*)'Etot=',Etot
00059
00060
00061 OPEN(105,FILE='./dir-wfn/dat.lattice')
00062 REWIND(105)
00063 READ(105,*) a1(1),a1(2),a1(3)
00064 READ(105,*) a2(1),a2(2),a2(3)
00065 READ(105,*) a3(1),a3(2),a3(3)
00066 CLOSE(105)
00067
00068 call OUTER_PRODUCT(a2(1),a3(1),b1(1))
00069 VOLUME=a1(1)*b1(1)+a1(2)*b1(2)+a1(3)*b1(3)
00070 b1(:)=b1(:)*tpi/VOLUME
00071 call OUTER_PRODUCT(a3(1),a1(1),b2(1))
00072 b2(:)=b2(:)*tpi/VOLUME
00073 call OUTER_PRODUCT(a1(1),a2(1),b3(1))
00074 b3(:)=b3(:)*tpi/VOLUME
00075
00076 write(6,*)
00077 write(6,*)'=============='
00078 write(6,*)'LATTICE VECTOR'
00079 write(6,*)'=============='
00080 write(6,*)
00081 write(6,*) a1(1), a1(2), a1(3)
00082 write(6,*) a2(1), a2(2), a2(3)
00083 write(6,*) a3(1), a3(2), a3(3)
00084 write(6,*)
00085 write(6,*)'VOLUME OF UNIT CELL=',VOLUME
00086 write(6,*)
00087 write(6,*)'========================='
00088 write(6,*)'RECIPROCAL LATTICE VECTOR'
00089 write(6,*)'========================='
00090 write(6,*)
00091 write(6,*) b1(1), b1(2), b1(3)
00092 write(6,*) b2(1), b2(2), b2(3)
00093 write(6,*) b3(1), b3(2), b3(3)
00094 write(6,*)
00095 write(6,*)'======================'
00096 write(6,*)'CHECK OF (A_VEC*B_VEC)'
00097 write(6,*)'======================'
00098 write(6,*)
00099 write(6,*)'a1*b1=',b1(1)*a1(1)+b1(2)*a1(2)+b1(3)*a1(3)
00100 write(6,*)'a2*b2=',b2(1)*a2(1)+b2(2)*a2(2)+b2(3)*a2(3)
00101 write(6,*)'a3*b3=',b3(1)*a3(1)+b3(2)*a3(2)+b3(3)*a3(3)
00102 call est_latparam(a1(1),a2(1),a3(1),a,b,c,alp,bet,gmm)
00103 write(6,'(6f15.10)') a,b,c,alp,bet,gmm
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 OPEN(108,FILE='./dir-wfn/dat.atom_position')
00129 rewind(108)
00130 read(108,*) nsi
00131 allocate(chemical_species(nsi))
00132 allocate(asi(3,nsi));asi=0.0d0
00133 do i=1,nsi
00134 read(108,*) chemical_species(i),(asi(j,i),j=1,3)
00135 enddo
00136 close(108)
00137
00138
00139 OPEN(100,FILE='./dir-wfn/dat.symmetry')
00140 rewind(100)
00141 read(100,*) nsymq
00142 read(100,*) nnp
00143 write(6,*)'nsymq=',nsymq
00144 write(6,*)'nnp=',nnp
00145 allocate(rg(3,3,nsymq));rg=0
00146 allocate(pg(3,nsymq));pg=0
00147 allocate(rginv(3,3,nsymq));rginv=0.0d0
00148 do iop=1,nsymq
00149 read(100,*)((rg(i,j,iop),i=1,3),j=1,3)
00150 read(100,*)(pg(i,iop),i=1,3)
00151 enddo
00152 close(100)
00153
00154 write(6,*)
00155 write(6,*)'========'
00156 write(6,*)'SYMMETRY'
00157 write(6,*)'========'
00158 write(6,*)
00159 do iop=1,nsymq
00160 write(6,'(a5,i4)')'sym=',iop
00161 write(6,'(a10)')'rg and pg'
00162 do i=1,3
00163 write(6,'(3i5,5x,i5)')(rg(i,j,iop),j=1,3),pg(i,iop)
00164 enddo
00165 write(6,*)
00166 enddo
00167
00168 rginv=rg
00169 do iop=1,nsymq
00170 call invmat(3,rginv(1,1,iop))
00171 enddo
00172
00173 do iop=1,nsymq
00174 write(6,'(a5,i4)')'sym=',iop
00175 write(6,'(a10)')'rginv'
00176 do i=1,3
00177 write(6,'(3f10.5)')(rginv(i,j,iop),j=1,3)
00178 enddo
00179 write(6,*)
00180 enddo
00181 write(6,*)'finish read dat.symmetry'
00182
00183
00184 OPEN(101,FILE='./dir-wfn/dat.sample-k')
00185 rewind(101)
00186 read(101,*) Nk_irr
00187 allocate(SKI(3,Nk_irr));SKI(:,:)=0.0D0
00188 do ik=1,Nk_irr
00189 read(101,*)(SKI(i,IK),i=1,3)
00190 enddo
00191 close(101)
00192
00193
00194
00195
00196
00197
00198 OPEN(132,FILE='./dir-wfn/dat.nkm')
00199 allocate(NGI(Nk_irr));NGI(:)=0
00200 rewind(132)
00201 do ik=1,Nk_irr
00202 read(132,*) NGI(ik)
00203 enddo
00204 close(132)
00205 NTG=maxval(abs(NGI(:)))
00206 write(6,'(a10,i10)')'NTG=',NTG
00207
00208
00209 OPEN(104,FILE='./dir-wfn/dat.kg')
00210 rewind(104)
00211 allocate(KGI(3,NTG,Nk_irr));KGI=0
00212 do ik=1,Nk_irr
00213 read(104,*) NKG_tmp
00214 if(NKG_tmp/=NGI(ik)) then
00215 write(6,*)'ERROR; STOP; NKG_tmp should be NGI(ik)'
00216 write(6,*)'NKG_tmp=',NKG_tmp,'NGI(ik)=',NGI(ik)
00217 stop
00218 endif
00219 do ig=1,NKG_tmp
00220 read(104,*)(KGI(i,ig,ik),i=1,3)
00221 enddo
00222 enddo
00223 close(104)
00224
00225 L1=maxval(abs(KGI(1,:,:)))+1; write(6,*)'L1=',L1
00226 L2=maxval(abs(KGI(2,:,:)))+1; write(6,*)'L2=',L2
00227 L3=maxval(abs(KGI(3,:,:)))+1; write(6,*)'L3=',L3
00228 allocate(packing(-L1:L1,-L2:L2,-L3:L3,Nk_irr))
00229 packing(:,:,:,:)=0
00230 do ik=1,Nk_irr
00231 do ig=1,NGI(ik)
00232 i1=KGI(1,ig,ik); j1=KGI(2,ig,ik); k1=KGI(3,ig,ik)
00233 packing(i1,j1,k1,ik)=ig
00234 enddo
00235 enddo
00236
00237 allocate(LKGI(NTG,Nk_irr));LKGI=0.0d0
00238 do ik=1,Nk_irr
00239 do ig=1,NGI(ik)
00240 ktmp(1)=(SKI(1,ik)+dble(KGI(1,ig,ik)))*b1(1)
00241 + +(SKI(2,ik)+dble(KGI(2,ig,ik)))*b2(1)
00242 + +(SKI(3,ik)+dble(KGI(3,ig,ik)))*b3(1)
00243
00244 ktmp(2)=(SKI(1,ik)+dble(KGI(1,ig,ik)))*b1(2)
00245 + +(SKI(2,ik)+dble(KGI(2,ig,ik)))*b2(2)
00246 + +(SKI(3,ik)+dble(KGI(3,ig,ik)))*b3(2)
00247
00248 ktmp(3)=(SKI(1,ik)+dble(KGI(1,ig,ik)))*b1(3)
00249 + +(SKI(2,ik)+dble(KGI(2,ig,ik)))*b2(3)
00250 + +(SKI(3,ik)+dble(KGI(3,ig,ik)))*b3(3)
00251 LKGI(ig,ik)=ktmp(1)**2+ktmp(2)**2+ktmp(3)**2
00252 enddo
00253 write(6,*) maxval(LKGI(:,ik))
00254 enddo
00255 write(6,*)
00256 write(6,*) maxval(LKGI(:,:))
00257 Ecut_for_psi=maxval(LKGI(:,:))+1.0d-8
00258 write(6,*)'Ecut_for_psi=',Ecut_for_psi
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269 htmp=dsqrt(dot_product(a1,a1))
00270 h1(:)=a1(:)/htmp
00271 htmp=dsqrt(dot_product(a2,a2))
00272 h2(:)=a2(:)/htmp
00273 htmp=dsqrt(dot_product(a3,a3))
00274 h3(:)=a3(:)/htmp
00275 d1=abs(dot_product(b1,h1))
00276 d2=abs(dot_product(b2,h2))
00277 d3=abs(dot_product(b3,h3))
00278 qwf=2.0d0*dsqrt(Ecut_for_psi)
00279
00280
00281
00282
00283
00284 nwx2=algn235(int(qwf/d1)+1)
00285 nwy2=algn235(int(qwf/d2)+1)
00286 nwz2=algn235(int(qwf/d3)+1)
00287
00288 nfft1=nwx2+1;nfft2=nwy2+1;nfft3=nwz2+1
00289 Nl123=nfft1*nfft2*nfft3
00290 write(6,'(a20,i10)')'nwx2:',nwx2
00291 write(6,'(a20,i10)')'nwy2:',nwy2
00292 write(6,'(a20,i10)')'nwz2:',nwz2
00293 write(6,'(a20,i10)')'nfft1:',nfft1
00294 write(6,'(a20,i10)')'nfft2:',nfft2
00295 write(6,'(a20,i10)')'nfft3:',nfft3
00296 write(6,'(a20,i10)')'NL123:',Nl123
00297
00298
00299
00300
00301 OPEN(111,FILE='./dir-wfn/dat.eigenvalue')
00302 rewind(111)
00303 read(111,*) NBAND
00304 allocate(E_EIGI(NBAND,Nk_irr));E_EIGI=0.0d0
00305 do ik=1,Nk_irr
00306 do ib=1,NBAND
00307 read(111,*) E_EIGI(ib,ik)
00308 enddo
00309 enddo
00310
00311
00312
00313
00314
00315
00316
00317 if(flg_specify_energywindow)then
00318 call calc_NB_start_end(NK_irr,NBAND,E_EIGI(1,1),
00319 + E_LOWER,E_UPPER,NB_start,NB_end)
00320
00321 endif
00322 if(flg_specify_bandindices)then
00323 write(6,'(3a)') "WARNING : ",
00324 + "flg_specify_bandindices=.true. : ",
00325 + "setting NB_start,NB_end to bandindex_min,bandindex_max"
00326 NB_start=bandindex_min
00327 NB_end=bandindex_max
00328 endif
00329
00330 NTB_glb=NBAND
00331 NTB=NB_end-NB_start+1
00332 write(6,'(a20,i10)')'NTB_glb:',NTB_glb
00333 write(6,'(a20,i10)')'NB_start:',NB_start
00334 write(6,'(a20,i10)')'NB_end:',NB_end
00335 write(6,'(a20,i10)')'NTB:',NTB
00336 deallocate(E_EIGI)
00337
00338
00339 OPEN(111,FILE='./dir-wfn/dat.eigenvalue')
00340 rewind(111)
00341 read(111,*) NBAND
00342 allocate(E_EIGI(NTB,Nk_irr));E_EIGI=0.0d0
00343 do ik=1,Nk_irr
00344 do ib=1,NBAND
00345 if((ib.ge.NB_start).and.(ib.le.NB_end)) then
00346 read(111,*) E_EIGI(ib-NB_start+1,ik)
00347 else
00348 read(111,*)
00349 endif
00350 enddo
00351 enddo
00352
00353 call est_NTK(Nk_irr,Nsymq,SKI(1,1),rg(1,1,1),NTK)
00354 write(6,'(a20,i10)')'Estimated NTK:',NTK
00355
00356
00357 allocate(SK0(3,NTK));SK0(:,:)=0.0d0
00358 allocate(numirr(NTK));numirr(:)=0
00359 allocate(numrot(NTK));numrot(:)=0
00360 allocate(trs(NTK));trs(:)=0
00361 allocate(RW(3,NTK));RW(:,:)=0
00362
00363
00364
00365
00366
00367
00368
00369 jk=0
00370 do ik=1,Nk_irr
00371
00372
00373
00374 do iop=1,Nsymq
00375
00376 ktmp(:)=0.0d0;RWtmp(:)=0
00377 ktmp(1)=rg(1,1,iop)*SKI(1,ik)+rg(1,2,iop)*SKI(2,ik)
00378 + +rg(1,3,iop)*SKI(3,ik)
00379 ktmp(2)=rg(2,1,iop)*SKI(1,ik)+rg(2,2,iop)*SKI(2,ik)
00380 + +rg(2,3,iop)*SKI(3,ik)
00381 ktmp(3)=rg(3,1,iop)*SKI(1,ik)+rg(3,2,iop)*SKI(2,ik)
00382 + +rg(3,3,iop)*SKI(3,ik)
00383 call kcheck(ktmp(1),RWtmp(1))
00384 do iik=1,jk
00385 if(abs(SK0(1,iik)-ktmp(1))<1.0d-4.and.
00386 + abs(SK0(2,iik)-ktmp(2))<1.0d-4.and.
00387 + abs(SK0(3,iik)-ktmp(3))<1.0d-4)goto 1000
00388 enddo
00389 jk=jk+1
00390
00391 SK0(:,jk)=ktmp(:)
00392 numirr(jk)=ik;numrot(jk)=iop;trs(jk)=1;RW(:,jk)=RWtmp(:)
00393
00394
00395 1000 ktmp(:)=0.0d0;RWtmp(:)=0
00396 ktmp(1)=rg(1,1,iop)*SKI(1,ik)+rg(1,2,iop)*SKI(2,ik)
00397 + +rg(1,3,iop)*SKI(3,ik)
00398 ktmp(2)=rg(2,1,iop)*SKI(1,ik)+rg(2,2,iop)*SKI(2,ik)
00399 + +rg(2,3,iop)*SKI(3,ik)
00400 ktmp(3)=rg(3,1,iop)*SKI(1,ik)+rg(3,2,iop)*SKI(2,ik)
00401 + +rg(3,3,iop)*SKI(3,ik)
00402 call kcheck_trs(ktmp(1),RWtmp(1))
00403 do iik=1,jk
00404 if(abs(SK0(1,iik)-(-ktmp(1)))<1.0d-4.and.
00405 + abs(SK0(2,iik)-(-ktmp(2)))<1.0d-4.and.
00406 + abs(SK0(3,iik)-(-ktmp(3)))<1.0d-4)goto 2000
00407 enddo
00408 jk=jk+1
00409
00410 SK0(:,jk)=-ktmp(:)
00411 numirr(jk)=ik;numrot(jk)=iop;trs(jk)=-1;RW(:,jk)=RWtmp(:)
00412 2000 enddo
00413 enddo
00414
00415 if(NTK/=jk) then
00416 write(6,*)'ERROR;STOP;NTK should be jk'
00417 write(6,*)'NTK=',NTK,'jk=',jk;STOP
00418 endif
00419
00420 write(6,*)
00421 write(6,*)'====================='
00422 write(6,*)'WRITE SAMPLE K-POINTS'
00423 write(6,*)'====================='
00424 do ik=1,NTK
00425 WRITE(6,'(I5,3F15.10)') ik,SK0(:,ik)
00426 enddo
00427
00428
00429
00430
00431
00432 call est_nkbi(NTK,SK0(1,1),nkb1,nkb2,nkb3)
00433 write(6,*)
00434 write(6,'(a24,4i10)')'nkb1,nkb2,nkb3,NTK=',nkb1,nkb2,nkb3,NTK
00435 Na1=nkb1/2;Na2=nkb2/2;Na3=nkb3/2
00436 write(6,'(a24,3i10)')'Na1,Na2,Na3=',Na1,Na2,Na3
00437
00438 call system('rm -rf dir-wan')
00439 call system('mkdir dir-wan')
00440
00441
00442 if(CALC_FERMISURFACE)then
00443 write(6,*)
00444 write(6,'(a30)')'write Fermi surface'
00445 write(6,*)
00446 allocate(E_EIG(NTB,NTK));E_EIG(:,:)=0.0d0
00447 do jk=1,NTK
00448 if(trs(jk)==1) then
00449 ik=numirr(jk)
00450 E_EIG(:,jk)=E_EIGI(:,ik)
00451 elseif(trs(jk)==-1)then
00452 ik=numirr(jk)
00453 E_EIG(:,jk)=E_EIGI(:,ik)
00454 endif
00455 enddo
00456 filename='./dir-wan/dat.frmsf'
00457 call wrt_frmsf(NTB,NTK,nkb1,nkb2,nkb3,E_EIG(1,1),SK0(1,1),
00458 + FermiEnergy,filename,b1(1),b2(1),b3(1))
00459 deallocate(E_EIG)
00460 endif
00461
00462
00463 if(CALC_GLOBAL_DOS)then
00464 write(6,*)
00465 write(6,'(a30)')'write global DOS'
00466 write(6,*)
00467 allocate(E_EIG(NTB,NTK));E_EIG(:,:)=0.0d0
00468 allocate(V_EIG(NTB,NTB,NTK));V_EIG(:,:,:)=0.0d0
00469 do jk=1,NTK
00470 if(trs(jk)==1) then
00471 ik=numirr(jk)
00472 E_EIG(:,jk)=E_EIGI(:,ik)
00473 elseif(trs(jk)==-1)then
00474 ik=numirr(jk)
00475 E_EIG(:,jk)=E_EIGI(:,ik)
00476 endif
00477 enddo
00478 filename='./dir-wan/dat.dos.global'
00479 call calculate_dos(NTB,NTK,nkb1,nkb2,nkb3,0.0d0,FermiEnergy,
00480 + filename,b1(1),b2(1),b3(1),SK0(1,1),E_EIG(1,1),V_EIG(1,1,1))
00481 deallocate(E_EIG,V_EIG)
00482 endif
00483
00484
00485 OPEN(102,FILE='./dir-wfn/dat.wfn',FORM='unformatted')
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500 rewind(102)
00501 read(102) ncomp
00502 if(ncomp/=1)then
00503 write(6,*)'This program not suport ncomp/=1; then stop'
00504 stop
00505 endif
00506 allocate(CIR(NTG,NTB,Nk_irr));CIR=0.0d0
00507 do ik=1,Nk_irr
00508 do ib=1,NBAND
00509 if((ib.ge.NB_start).and.(ib.le.NB_end)) then
00510 read(102)(CIR(IG,ib-NB_start+1,ik),IG=1,NGI(ik))
00511 else
00512 read(102)
00513 endif
00514 enddo
00515 enddo
00516 close(102)
00517
00518
00519 allocate(C0(NTG,NTB,NTK));C0(:,:,:)=0.0d0
00520 allocate(E_EIG(NTB,NTK));E_EIG(:,:)=0.0d0
00521 allocate(NB0(NTK));NB0(:)=0
00522 allocate(NG0(NTK));NG0(:)=0
00523 allocate(KG0(3,NTG,NTK));KG0(:,:,:)=0
00524 allocate(KGtmp(3,NTG));KGtmp(:,:)=0
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535 do jk=1,NTK
00536 if(trs(jk)==1) then
00537 ik=numirr(jk);iop=numrot(jk)
00538 ktmp(1)=rg(1,1,iop)*SKI(1,ik)+rg(1,2,iop)*SKI(2,ik)
00539 + +rg(1,3,iop)*SKI(3,ik)+dble(RW(1,jk))
00540 ktmp(2)=rg(2,1,iop)*SKI(1,ik)+rg(2,2,iop)*SKI(2,ik)
00541 + +rg(2,3,iop)*SKI(3,ik)+dble(RW(2,jk))
00542 ktmp(3)=rg(3,1,iop)*SKI(1,ik)+rg(3,2,iop)*SKI(2,ik)
00543 + +rg(3,3,iop)*SKI(3,ik)+dble(RW(3,jk))
00544 E_EIG(:,jk)=E_EIGI(:,ik)
00545 NB0(jk)=NBAND
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570 call make_KG0(NTG,b1(1),b2(1),b3(1),Ecut_for_psi,
00571 + ktmp(1),ktmp(2),ktmp(3),KG0(1,1,jk),NG_for_psi)
00572 if(NG_for_psi/=NGI(ik)) then
00573 write(6,*)'ERROR; STOP; NG_for_psi should be NGI(ik)'
00574 write(6,*)'NG_for_psi=',NG_for_psi,'NG0(ik)=',NGI(ik)
00575 write(6,*)'ik,jk',ik,jk;STOP
00576 endif
00577 NG0(jk)=NG_for_psi
00578 call make_C0(NTG,NTB,trs(jk),NG0(jk),KG0(1,1,jk),RW(1,jk),
00579 + rginv(1,1,iop),pg(1,iop),nnp,L1,L2,L3,packing(-L1,-L2,-L3,ik),
00580 + CIR(1,1,ik),C0(1,1,jk))
00581
00582 elseif(trs(jk)==-1)then
00583 ik=numirr(jk);iop=numrot(jk)
00584 ktmp(1)=rg(1,1,iop)*SKI(1,ik)+rg(1,2,iop)*SKI(2,ik)
00585 + +rg(1,3,iop)*SKI(3,ik)+dble(RW(1,jk))
00586 ktmp(2)=rg(2,1,iop)*SKI(1,ik)+rg(2,2,iop)*SKI(2,ik)
00587 + +rg(2,3,iop)*SKI(3,ik)+dble(RW(2,jk))
00588 ktmp(3)=rg(3,1,iop)*SKI(1,ik)+rg(3,2,iop)*SKI(2,ik)
00589 + +rg(3,3,iop)*SKI(3,ik)+dble(RW(3,jk))
00590
00591 E_EIG(:,jk)=E_EIGI(:,ik)
00592 NB0(jk)=NBAND
00593 KGtmp(:,:)=0
00594 call make_KG0(NTG,b1(1),b2(1),b3(1),Ecut_for_psi,
00595 + ktmp(1),ktmp(2),ktmp(3),KGtmp(1,1),NG_for_psi)
00596 if(NG_for_psi/=NGI(ik))then
00597 write(6,*)'ERROR; STOP; NG_for_psi should be NGI(ik)'
00598 write(6,*)'NG_for_psi=',NG_for_psi,'NGI(ik)=',NGI(ik);STOP
00599 endif
00600 NG0(jk)=NG_for_psi
00601 KG0(:,:,jk)=-KGtmp(:,:)
00602 call make_C0(NTG,NTB,trs(jk),NG0(jk),KGtmp(1,1),RW(1,jk),
00603 + rginv(1,1,iop),pg(1,iop),nnp,L1,L2,L3,packing(-L1,-L2,-L3,ik),
00604 + CIR(1,1,ik),C0(1,1,jk))
00605 endif
00606 enddo
00607
00608 WRITE(6,*)
00609 write(6,*)'==================='
00610 write(6,*)'CALCULATED K-POINTS'
00611 write(6,*)'==================='
00612 WRITE(6,*)
00613 do ik=1,NTK
00614 write(6,'(i5,3f15.10,i5)') ik,SK0(:,ik),trs(ik)
00615 enddo
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631 allocate(N_BAND(NTK));N_BAND(:)=0
00632 allocate(N_BAND_BTM(NTK));N_BAND_BTM(:)=0
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647 write(6,*)
00648 write(6,'(2a)') "###STEP : ",
00649 + "Calculate number of bands inside outer window"
00650 write(6,*)
00651 if(flg_specify_energywindow)then
00652 do ik=1,NTK
00653 SUM_INT=0
00654 do ib=1,NTB
00655 IF(E_LOWER<=E_EIG(ib,ik).AND.E_EIG(ib,ik)<=E_UPPER)THEN
00656 IF(SUM_INT==0) N_BAND_BTM(ik)=ib-1
00657 SUM_INT=SUM_INT+1
00658 write(6,'(a,2i8,f8.4,a,2f8.4)') "ik,ib=",
00659 + ik,ib,E_EIG(ib,ik),
00660 + " inside ",E_LOWER,E_UPPER
00661 else
00662 write(6,'(a,2i8,f8.4,a,2f8.4)') "ik,ib=",
00663 + ik,ib,E_EIG(ib,ik),
00664 + " outside ",E_LOWER,E_UPPER
00665 ENDIF
00666 enddo
00667 N_BAND(ik)=SUM_INT
00668 enddo
00669 endif
00670 if(flg_specify_bandindices)then
00671 do ik=1,NTK
00672 N_BAND(ik)=NTB
00673 N_BAND_BTM(ik)=0
00674 enddo
00675 endif
00676
00677 write(6,*)
00678 write(6,*)'===================================='
00679 write(6,*)'BAND INFO INSIDE OUTER ENERGY WINDOW'
00680 write(6,*)'===================================='
00681 write(6,*)
00682
00683 if(maxval(N_BAND)/=NTB)then
00684 write(6,'(a)')'maxval(N_BAND) SHOULD BE THE SAME AS NTB; STOP'
00685 endif
00686
00687 do ik=1,NTK
00688
00689 if(N_BAND(ik)<N_wannier)then
00690 write(6,'(a)')'energy window setting: WRONG'
00691 write(6,'(a)')'N_BAND(ik) SHOULD BE LARGER THAN N_wannier: STOP'
00692 write(6,'(a,3i8)')'ik, N_BAND(ik), N_wannier',
00693 + ik, N_BAND(ik), N_wannier
00694 stop
00695 endif
00696 write(6,*) ik,N_BAND_BTM(ik),N_BAND(ik)
00697 enddo
00698 write(6,*)
00699
00700
00701 OPEN(149,FILE='./dir-wan/dat.ns-nb')
00702 rewind(149)
00703 do ik=1,NTK
00704 write(149,*) N_BAND_BTM(ik)+NB_start-1,N_BAND(ik)
00705 enddo
00706
00707
00708 if(CALC_REAL_SPACE_BLOCH)then
00709 write(6,*)
00710 write(6,*)'================================'
00711 write(6,*)'=== VISUALIZE BLOCH FUNCTION ==='
00712 write(6,*)'================================'
00713 write(6,*)
00714
00715 nfft1=nwx2+1;nfft2=nwy2+1;nfft3=nwz2+1
00716 Nl123=nfft1*nfft2*nfft3
00717 write(6,'(a20,3i10)')'nwx2,nwy2,nwz2:',nwx2,nwy2,nwz2
00718 write(6,'(a20,3i10)')'nfft1,nfft2,nfft2:',nfft1,nfft2,nfft3
00719 write(6,'(a20,i10)')'Nl123:',Nl123
00720 call fft3_init(nwx2,nwy2,nwz2,nfft1,nfft2,nfft3,fs)
00721 allocate(fftwk(Nl123*2),stat=err)
00722 allocate(wfunc(Nl123*2),stat=err)
00723
00724 write(6,'(a20,3f15.10)')'calc_k:',calc_k
00725 call search_ik(NTK,SK0(1,1),calc_k(1),ik)
00726 write(6,'(a20,i10)')'ik:',ik
00727
00728 amax=1; bmax=1; cmax=1
00729 if(CALC_BLOCH_VIS_RANGE)then
00730 write(6,*)
00731 write(6,'(a40)')'CALCULATE BLOCH VISUALIZATION RANGE'
00732 write(6,*)
00733 call search_vis_range(calc_k(1),amin,amax,bmin,bmax,cmin,cmax)
00734 endif
00735
00736 write(6,'(a20,3i10)')'supercell:',amax,bmax,cmax
00737 allocate(BF_REALSPACE((amin)*nwx2:(amax)*nwx2-1,
00738 + (bmin)*nwy2:(bmax)*nwy2-1,
00739 + (cmin)*nwz2:(cmax)*nwz2-1))
00740
00741 do ib=1,N_BAND(ik)
00742 BF_REALSPACE=0.0d0
00743 call make_pw(C0(1,N_BAND_BTM(ik)+ib,ik),wfunc(1),fftwk(1),
00744 + NG0(ik),KG0(1,1,ik),NTG,nwx2,nwy2,nwz2,nfft1,nfft2,Nl123,fs,
00745 + amin,amax,bmin,bmax,cmin,cmax,SK0(1,ik),
00746 + BF_REALSPACE(amin*nwx2,bmin*nwy2,cmin*nwz2))
00747 call phase_fix_bf(nwx2,nwy2,nwz2,amin,amax,bmin,bmax,
00748 + cmin,cmax,BF_REALSPACE(amin*nwx2,bmin*nwy2,cmin*nwz2))
00749
00750
00751
00752
00753
00754
00755 write(filename,"('dat.bf-realspace-',i3.3,'.grd')")ib
00756 OPEN(116,FILE=filename)
00757 REWIND(116)
00758 na_grid=nwx2
00759 nb_grid=nwy2
00760 nc_grid=nwz2
00761 aa1=a1*dble(amax-amin)
00762 aa2=a2*dble(bmax-bmin)
00763 aa3=a3*dble(cmax-cmin)
00764 call est_latparam(aa1(1),aa2(1),aa3(1),a,b,c,alp,bet,gmm)
00765 write(116,'(a40,f10.5)')'Electron density and energy (eV)',
00766 + E_EIG(N_BAND_BTM(ik)+ib,ik)*au
00767 write(116,'(6f15.10)') a*bohr,b*bohr,c*bohr,alp,bet,gmm
00768 write(116,'(3I5)')
00769 + (amax-amin)*na_grid,(bmax-bmin)*nb_grid,(cmax-cmin)*nc_grid
00770 write(116,'(6g25.16)')(((dble(BF_REALSPACE(i1,i2,i3)),
00771 + i3=(cmin)*nc_grid,(cmax)*nc_grid-1),
00772 + i2=(bmin)*nb_grid,(bmax)*nb_grid-1),
00773 + i1=(amin)*na_grid,(amax)*na_grid-1)
00774 CLOSE(116)
00775 enddo
00776 deallocate(BF_REALSPACE,fftwk,wfunc)
00777
00778
00779 call gencif(a1(1),a2(1),a3(1),amin,amax,bmin,bmax,cmin,cmax,
00780 + nsi,chemical_species(1),asi(1,1))
00781
00782 write(6,*)
00783 write(6,*)'==================================='
00784 write(6,*)'FINIDH VISUALIZE BLOCH FUNCTION ==='
00785 write(6,*)'!!!!!! We stop calculation !!!!!!!!'
00786 write(6,*)'==================================='
00787 write(6,*)
00788 stop
00789 endif
00790
00791 WRITE(6,*)
00792 WRITE(6,*)'========================================='
00793 WRITE(6,*)'=== WE NOW CALCULATE WANNIER FUNCTION ==='
00794 WRITE(6,*)'========================================='
00795 WRITE(6,*)
00796
00797
00798
00799
00800
00801 select case(icell)
00802 case(0)
00803
00804
00805
00806 allocate(b_LAT(3,NBMAX));b_LAT(:,:)=0.0D0
00807 allocate(VEC_b(3,NBMAX));VEC_b(:,:)=0.0d0
00808 allocate(WEIGHT_b(NBMAX));WEIGHT_b(:)=0.0D0
00809 call generate_bvectors(b1(1),b2(1),b3(1),nkb1,nkb2,nkb3,
00810 + NBMAX,maxshell,verbosity,NB,b_LAT(1,1),VEC_b(1,1),WEIGHT_b(1))
00811
00812 case(1)
00813 NB=6
00814 allocate(b_LAT(3,NB));b_LAT(:,:)=0.0D0
00815 allocate(VEC_b(3,NB));VEC_b(:,:)=0.0d0
00816 allocate(WEIGHT_b(NB));WEIGHT_b(:)=0.0D0
00817 b_LAT(:,:) = 0.0D0
00818 b_LAT(1,1) = (1.0D0/DBLE(nkb1))
00819 b_LAT(2,2) = (1.0D0/DBLE(nkb1))
00820 b_LAT(3,3) = (1.0D0/DBLE(nkb1))
00821 b_LAT(1,4) = -(1.0D0/DBLE(nkb1))
00822 b_LAT(2,5) = -(1.0D0/DBLE(nkb1))
00823 b_LAT(3,6) = -(1.0D0/DBLE(nkb1))
00824 VEC_b(:,:) = 0.0D0
00825 VEC_b(:,1) = (1.0D0/DBLE(nkb1)) * b1(:)
00826 VEC_b(:,2) = (1.0D0/DBLE(nkb1)) * b2(:)
00827 VEC_b(:,3) = (1.0D0/DBLE(nkb1)) * b3(:)
00828 VEC_b(:,4) = -(1.0D0/DBLE(nkb1)) * b1(:)
00829 VEC_b(:,5) = -(1.0D0/DBLE(nkb1)) * b2(:)
00830 VEC_b(:,6) = -(1.0D0/DBLE(nkb1)) * b3(:)
00831 VEC_b2=VEC_b(1,1)**2+VEC_b(2,1)**2+VEC_b(3,1)**2
00832 WEIGHT_b(:) = 3.0d0/dble(NB)/VEC_b2
00833
00834
00835 case(2)
00836 NB=8
00837 allocate(b_LAT(3,NB));b_LAT(:,:)=0.0D0
00838 allocate(VEC_b(3,NB));VEC_b(:,:)=0.0d0
00839 allocate(WEIGHT_b(NB));WEIGHT_b(:)=0.0D0
00840 b_LAT(:,:) = 0.0D0
00841 b_LAT(1,1) = (1.0D0/DBLE(nkb1))
00842 b_LAT(2,2) = (1.0D0/DBLE(nkb1))
00843 b_LAT(3,3) = (1.0D0/DBLE(nkb1))
00844 b_LAT(1,4) = (1.0D0/DBLE(nkb1))
00845 b_LAT(2,4) = (1.0D0/DBLE(nkb1))
00846 b_LAT(3,4) = (1.0D0/DBLE(nkb1))
00847 b_LAT(1,5) = -(1.0D0/DBLE(nkb1))
00848 b_LAT(2,6) = -(1.0D0/DBLE(nkb1))
00849 b_LAT(3,7) = -(1.0D0/DBLE(nkb1))
00850 b_LAT(1,8) = -(1.0D0/DBLE(nkb1))
00851 b_LAT(2,8) = -(1.0D0/DBLE(nkb1))
00852 b_LAT(3,8) = -(1.0D0/DBLE(nkb1))
00853 VEC_b(:,:) = 0.0D0
00854 VEC_b(:,1) = (1.0D0/DBLE(nkb1)) * b1(:)
00855 VEC_b(:,2) = (1.0D0/DBLE(nkb1)) * b2(:)
00856 VEC_b(:,3) = (1.0D0/DBLE(nkb1)) * b3(:)
00857 VEC_b(:,4) = (1.0D0/DBLE(nkb1)) *(b1(:)+b2(:)+b3(:))
00858 VEC_b(:,5) = -(1.0D0/DBLE(nkb1)) * b1(:)
00859 VEC_b(:,6) = -(1.0D0/DBLE(nkb1)) * b2(:)
00860 VEC_b(:,7) = -(1.0D0/DBLE(nkb1)) * b3(:)
00861 VEC_b(:,8) = -(1.0D0/DBLE(nkb1)) *(b1(:)+b2(:)+b3(:))
00862 VEC_b2=VEC_b(1,1)**2+VEC_b(2,1)**2+VEC_b(3,1)**2
00863 WEIGHT_b(:) = 3.0d0/dble(NB)/VEC_b2
00864
00865
00866 case(3)
00867 NB=12
00868 allocate(b_LAT(3,NB));b_LAT(:,:)=0.0D0
00869 allocate(VEC_b(3,NB));VEC_b(:,:)=0.0d0
00870 allocate(WEIGHT_b(NB));WEIGHT_b(:)=0.0D0
00871 b_LAT(:,:) = 0.0D0
00872 b_LAT(1, 1) = (1.0D0/DBLE(nkb1))
00873 b_LAT(2, 2) = (1.0D0/DBLE(nkb1))
00874 b_LAT(3, 3) = (1.0D0/DBLE(nkb1))
00875 b_LAT(1, 4) = (1.0D0/DBLE(nkb1))
00876 b_LAT(2, 4) = -(1.0D0/DBLE(nkb1))
00877 b_LAT(2, 5) = (1.0D0/DBLE(nkb1))
00878 b_LAT(3, 5) = -(1.0D0/DBLE(nkb1))
00879 b_LAT(3, 6) = (1.0D0/DBLE(nkb1))
00880 b_LAT(1, 6) = -(1.0D0/DBLE(nkb1))
00881 b_LAT(1, 7) = -(1.0D0/DBLE(nkb1))
00882 b_LAT(2, 8) = -(1.0D0/DBLE(nkb1))
00883 b_LAT(3, 9) = -(1.0D0/DBLE(nkb1))
00884 b_LAT(1,10) = -(1.0D0/DBLE(nkb1))
00885 b_LAT(2,10) = (1.0D0/DBLE(nkb1))
00886 b_LAT(2,11) = -(1.0D0/DBLE(nkb1))
00887 b_LAT(3,11) = (1.0D0/DBLE(nkb1))
00888 b_LAT(3,12) = -(1.0D0/DBLE(nkb1))
00889 b_LAT(1,12) = (1.0D0/DBLE(nkb1))
00890 VEC_b(:,:) = 0.0D0
00891 VEC_b(:, 1) = (1.0D0/DBLE(nkb1)) * b1(:)
00892 VEC_b(:, 2) = (1.0D0/DBLE(nkb1)) * b2(:)
00893 VEC_b(:, 3) = (1.0D0/DBLE(nkb1)) * b3(:)
00894 VEC_b(:, 4) = (1.0D0/DBLE(nkb1)) *(b1(:)-b2(:))
00895 VEC_b(:, 5) = (1.0D0/DBLE(nkb1)) *(b2(:)-b3(:))
00896 VEC_b(:, 6) = (1.0D0/DBLE(nkb1)) *(b3(:)-b1(:))
00897 VEC_b(:, 7) = -(1.0D0/DBLE(nkb1)) * b1(:)
00898 VEC_b(:, 8) = -(1.0D0/DBLE(nkb1)) * b2(:)
00899 VEC_b(:, 9) = -(1.0D0/DBLE(nkb1)) * b3(:)
00900 VEC_b(:,10) = -(1.0D0/DBLE(nkb1)) *(b1(:)-b2(:))
00901 VEC_b(:,11) = -(1.0D0/DBLE(nkb1)) *(b2(:)-b3(:))
00902 VEC_b(:,12) = -(1.0D0/DBLE(nkb1)) *(b3(:)-b1(:))
00903 VEC_b2=VEC_b(1,1)**2+VEC_b(2,1)**2+VEC_b(3,1)**2
00904 WEIGHT_b(:) = 3.0d0/dble(NB)/VEC_b2
00905
00906
00907 case(4)
00908 NB=8
00909 allocate(b_LAT(3,NB));b_LAT(:,:)=0.0D0
00910 allocate(VEC_b(3,NB));VEC_b(:,:)=0.0d0
00911 allocate(WEIGHT_b(NB));WEIGHT_b(:)=0.0D0
00912 b_LAT(:,:) = 0.0D0
00913 b_LAT(1,1) = (1.0D0/DBLE(nkb1))
00914 b_LAT(2,2) = (1.0D0/DBLE(nkb2))
00915 b_LAT(1,3) = (1.0D0/DBLE(nkb1))
00916 b_LAT(2,3) = -(1.0D0/DBLE(nkb2))
00917
00918 b_LAT(1,4) = -(1.0D0/DBLE(nkb1))
00919 b_LAT(2,5) = -(1.0D0/DBLE(nkb2))
00920 b_LAT(1,6) = -(1.0D0/DBLE(nkb1))
00921 b_LAT(2,6) = (1.0D0/DBLE(nkb2))
00922
00923 b_LAT(3,7) = (1.0D0/DBLE(nkb3))
00924 b_LAT(3,8) = -(1.0D0/DBLE(nkb3))
00925 VEC_b(:,:) = 0.0D0
00926 VEC_b(:,1) = (1.0D0/DBLE(nkb1)) * b1(:)
00927 VEC_b(:,2) = (1.0D0/DBLE(nkb2)) * b2(:)
00928 VEC_b(:,3) = (1.0D0/DBLE(nkb1)) * b1(:)
00929 + - (1.0D0/DBLE(nkb2)) * b2(:)
00930
00931 VEC_b(:,4) = -(1.0D0/DBLE(nkb1)) * b1(:)
00932 VEC_b(:,5) = -(1.0D0/DBLE(nkb2)) * b2(:)
00933 VEC_b(:,6) = -(1.0D0/DBLE(nkb1)) * b1(:)
00934 + + (1.0D0/DBLE(nkb2)) * b2(:)
00935
00936 VEC_b(:,7) = (1.0D0/DBLE(nkb3)) * b3(:)
00937 VEC_b(:,8) = -(1.0D0/DBLE(nkb3)) * b3(:)
00938 La = DSQRT( a1(1)**2 + a1(2)**2 + a1(3)**2 )
00939 Lc = DSQRT( a3(1)**2 + a3(2)**2 + a3(3)**2 )
00940 WEIGHT_b(:) = 0.0D0
00941 WEIGHT_b(1) = (1.0D0/4.0D0)*(La*nkb1/tpi)**2
00942 WEIGHT_b(2) = (1.0D0/4.0D0)*(La*nkb1/tpi)**2
00943 WEIGHT_b(3) = (1.0D0/4.0D0)*(La*nkb1/tpi)**2
00944 WEIGHT_b(4) = (1.0D0/4.0D0)*(La*nkb1/tpi)**2
00945 WEIGHT_b(5) = (1.0D0/4.0D0)*(La*nkb1/tpi)**2
00946 WEIGHT_b(6) = (1.0D0/4.0D0)*(La*nkb1/tpi)**2
00947 WEIGHT_b(7) = (1.0D0/2.0D0)*(Lc*nkb3/tpi)**2
00948 WEIGHT_b(8) = (1.0D0/2.0D0)*(Lc*nkb3/tpi)**2
00949
00950 case(5)
00951 NB=6
00952 allocate(b_LAT(3,NB));b_LAT(:,:)=0.0D0
00953 allocate(VEC_b(3,NB));VEC_b(:,:)=0.0d0
00954 allocate(WEIGHT_b(NB));WEIGHT_b(:)=0.0D0
00955 b_LAT(:,:) = 0.0D0
00956 b_LAT(1,1) = (1.0D0/DBLE(nkb1))
00957 b_LAT(2,2) = (1.0D0/DBLE(nkb2))
00958 b_LAT(3,3) = (1.0D0/DBLE(nkb3))
00959 b_LAT(1,4) = -(1.0D0/DBLE(nkb1))
00960 b_LAT(2,5) = -(1.0D0/DBLE(nkb2))
00961 b_LAT(3,6) = -(1.0D0/DBLE(nkb3))
00962 VEC_b(:,:) = 0.0D0
00963 VEC_b(:,1) = (1.0D0/DBLE(nkb1)) * b1(:)
00964 VEC_b(:,2) = (1.0D0/DBLE(nkb2)) * b2(:)
00965 VEC_b(:,3) = (1.0D0/DBLE(nkb3)) * b3(:)
00966 VEC_b(:,4) = -(1.0D0/DBLE(nkb1)) * b1(:)
00967 VEC_b(:,5) = -(1.0D0/DBLE(nkb2)) * b2(:)
00968 VEC_b(:,6) = -(1.0D0/DBLE(nkb3)) * b3(:)
00969 La = DSQRT( a1(1)**2 + a1(2)**2 + a1(3)**2 )
00970 Lb = DSQRT( a2(1)**2 + a2(2)**2 + a2(3)**2 )
00971 Lc = DSQRT( a3(1)**2 + a3(2)**2 + a3(3)**2 )
00972 WEIGHT_b(:) = 0.0D0
00973 WEIGHT_b(1) = (1.0D0/2.0D0)*(La*nkb1/tpi)**2
00974 WEIGHT_b(2) = (1.0D0/2.0D0)*(Lb*nkb2/tpi)**2
00975 WEIGHT_b(3) = (1.0D0/2.0D0)*(Lc*nkb3/tpi)**2
00976 WEIGHT_b(4) = (1.0D0/2.0D0)*(La*nkb1/tpi)**2
00977 WEIGHT_b(5) = (1.0D0/2.0D0)*(Lb*nkb2/tpi)**2
00978 WEIGHT_b(6) = (1.0D0/2.0D0)*(Lc*nkb3/tpi)**2
00979
00980 case(6)
00981 NB=8
00982 allocate(b_LAT(3,NB));b_LAT(:,:)=0.0D0
00983 allocate(VEC_b(3,NB));VEC_b(:,:)=0.0d0
00984 allocate(WEIGHT_b(NB));WEIGHT_b(:)=0.0D0
00985 b_LAT(:,:) = 0.0D0
00986 b_LAT(1,1) = (1.0D0/DBLE(nkb1))
00987 b_LAT(3,2) = (1.0D0/DBLE(nkb3))
00988 b_LAT(1,3) = (1.0D0/DBLE(nkb1))
00989 b_LAT(3,3) = -(1.0D0/DBLE(nkb3))
00990 b_LAT(1,4) = -(1.0D0/DBLE(nkb1))
00991 b_LAT(3,5) = -(1.0D0/DBLE(nkb3))
00992 b_LAT(1,6) = -(1.0D0/DBLE(nkb1))
00993 b_LAT(3,6) = (1.0D0/DBLE(nkb3))
00994 b_LAT(2,7) = (1.0D0/DBLE(nkb2))
00995 b_LAT(2,8) = -(1.0D0/DBLE(nkb2))
00996 VEC_b(:,:) = 0.0D0
00997 VEC_b(:,1) = (1.0D0/DBLE(nkb1))*b1(:)
00998 VEC_b(:,2) = (1.0D0/DBLE(nkb3))*b3(:)
00999 VEC_b(:,3) = VEC_b(:,1)-VEC_b(:,2)
01000 VEC_b(:,4) =-VEC_b(:,1)
01001 VEC_b(:,5) =-VEC_b(:,2)
01002 VEC_b(:,6) =-VEC_b(:,3)
01003 VEC_b(:,7) = (1.0D0/DBLE(nkb2))*b2(:)
01004 VEC_b(:,8) =-VEC_b(:,7)
01005 s = VEC_b(1,1)
01006 t = VEC_b(3,1)
01007 u = VEC_b(3,2)
01008 v = VEC_b(2,7)
01009 WEIGHT_b(:) = 0.0D0
01010 WEIGHT_b(1) = (u-t)/(2.0d0*(s**2)*u)
01011 WEIGHT_b(2) = (s**2-t*(u-t))/(2.0d0*(s**2)*(u**2))
01012 WEIGHT_b(3) = t/(2.0d0*(s**2)*u)
01013 WEIGHT_b(4) = WEIGHT_b(1)
01014 WEIGHT_b(5) = WEIGHT_b(2)
01015 WEIGHT_b(6) = WEIGHT_b(3)
01016 WEIGHT_b(7) = 1.0d0/(2.0d0*(v**2))
01017 WEIGHT_b(8) = WEIGHT_b(7)
01018
01019 case(8)
01020 write(6,*)
01021 write(6,*)'==============================='
01022 write(6,*)' body centered tetragonal cell '
01023 write(6,*)'==============================='
01024 write(6,*)
01025 NB=10
01026 allocate(b_LAT(3,NB));b_LAT(:,:)=0.0D0
01027 allocate(VEC_b(3,NB));VEC_b(:,:)=0.0d0
01028 allocate(WEIGHT_b(NB));WEIGHT_b(:)=0.0D0
01029 if(nkb1/=nkb2)stop'nkb1/=nkb2'
01030 if(nkb1/=nkb3)stop'nkb1/=nkb3'
01031 if(nkb2/=nkb3)stop'nkb2/=nkb3'
01032 b_LAT(:,:)=0.0D0
01033 b_LAT(1,1)=(1.0D0/DBLE(nkb1))
01034 b_LAT(2,1)=(1.0D0/DBLE(nkb2))
01035 b_LAT(3,1)=-(1.0D0/DBLE(nkb3))
01036 b_LAT(1,2)=(1.0D0/DBLE(nkb1))
01037 b_LAT(2,3)=(1.0D0/DBLE(nkb2))
01038 b_LAT(1,4)=(1.0D0/DBLE(nkb1))
01039 b_LAT(3,4)=-(1.0D0/DBLE(nkb3))
01040 b_LAT(2,5)=(1.0D0/DBLE(nkb2))
01041 b_LAT(3,5)=-(1.0D0/DBLE(nkb3))
01042 b_LAT(:,6)=-b_LAT(:,1)
01043 b_LAT(:,7)=-b_LAT(:,2)
01044 b_LAT(:,8)=-b_LAT(:,3)
01045 b_LAT(:,9)=-b_LAT(:,4)
01046 b_LAT(:,10)=-b_LAT(:,5)
01047 VEC_b(:,:)=0.0D0
01048 VEC_b(:,1)=(1.0D0/DBLE(nkb1))*(b1(:)+b2(:)-b3(:))
01049 VEC_b(:,2)=(1.0D0/DBLE(nkb1))*b1(:)
01050 VEC_b(:,3)=(1.0D0/DBLE(nkb2))*b2(:)
01051 VEC_b(:,4)=(1.0D0/DBLE(nkb1))*b1(:)-(1.0D0/DBLE(nkb3))*b3(:)
01052 VEC_b(:,5)=(1.0D0/DBLE(nkb2))*b2(:)-(1.0D0/DBLE(nkb3))*b3(:)
01053 VEC_b(:,6)=-VEC_b(:,1)
01054 VEC_b(:,7)=-VEC_b(:,2)
01055 VEC_b(:,8)=-VEC_b(:,3)
01056 VEC_b(:,9)=-VEC_b(:,4)
01057 VEC_b(:,10)=-VEC_b(:,5)
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095 NBh=NB/2
01096 allocate(VEC_d(6));VEC_d(:)=0.0d0
01097 allocate(aa(6,6));aa(:,:)=0.0d0
01098 allocate(bb(6,6));bb(:,:)=0.0d0
01099 do ib=1,Nbh
01100 ab=0
01101 do ix=1,3
01102 do jx=ix,3
01103 ab=ab+1
01104 aa(ib,ab)=VEC_b(ix,ib)*VEC_b(jx,ib)
01105 enddo
01106 enddo
01107 enddo
01108
01109
01110 do i=1,6
01111 do j=1,6
01112 bb(j,i)=aa(i,j)
01113 enddo
01114 enddo
01115
01116
01117
01118
01119
01120 call inv_ge_lapack(6,NBh,bb(1,1))
01121
01122
01123
01124
01125
01126
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139
01140 VEC_d(1)=0.5d0
01141 VEC_d(2)=0.0d0
01142 VEC_d(3)=0.0d0
01143 VEC_d(4)=0.5d0
01144 VEC_d(5)=0.0d0
01145 VEC_d(6)=0.5d0
01146
01147 WEIGHT_b(:)=0.0D0
01148 do ib=1,Nbh
01149 ab=0
01150 SUM_REAL=0.0d0
01151 do ix=1,3
01152 do jx=ix,3
01153 ab=ab+1
01154
01155 SUM_REAL=SUM_REAL+VEC_d(ab)*bb(ib,ab)
01156 enddo
01157 enddo
01158 WEIGHT_b(ib)=SUM_REAL
01159 enddo
01160 WEIGHT_b(6)=WEIGHT_b(1)
01161 WEIGHT_b(7)=WEIGHT_b(2)
01162 WEIGHT_b(8)=WEIGHT_b(3)
01163 WEIGHT_b(9)=WEIGHT_b(4)
01164 WEIGHT_b(10)=WEIGHT_b(5)
01165
01166 case(7)
01167 NB=12
01168 allocate(b_LAT(3,NB));b_LAT(:,:)=0.0D0
01169 allocate(VEC_b(3,NB));VEC_b(:,:)=0.0d0
01170 allocate(WEIGHT_b(NB));WEIGHT_b(:)=0.0D0
01171 NBh=NB/2
01172 b_LAT(:,:) = 0.0D0
01173 b_LAT(1,1) = (1.0D0/DBLE(nkb1))
01174 b_LAT(2,2) = (1.0D0/DBLE(nkb2))
01175 b_LAT(3,3) = (1.0D0/DBLE(nkb3))
01176 b_LAT(1,4) = (1.0D0/DBLE(nkb1))
01177 b_LAT(2,4) = -(1.0D0/DBLE(nkb2))
01178 b_LAT(2,5) = (1.0D0/DBLE(nkb2))
01179 b_LAT(3,5) = -(1.0D0/DBLE(nkb3))
01180 b_LAT(3,6) = (1.0D0/DBLE(nkb3))
01181 b_LAT(1,6) = -(1.0D0/DBLE(nkb1))
01182 b_LAT(:,7) = -b_LAT(:,1)
01183 b_LAT(:,8) = -b_LAT(:,2)
01184 b_LAT(:,9) = -b_LAT(:,3)
01185 b_LAT(:,10)= -b_LAT(:,4)
01186 b_LAT(:,11)= -b_LAT(:,5)
01187 b_LAT(:,12)= -b_LAT(:,6)
01188 VEC_b(:,:) = 0.0D0
01189 VEC_b(:,1) = (1.0D0/DBLE(nkb1))*b1(:)
01190 VEC_b(:,2) = (1.0D0/DBLE(nkb2))*b2(:)
01191 VEC_b(:,3) = (1.0D0/DBLE(nkb3))*b3(:)
01192 VEC_b(:,4) = VEC_b(:,1)-VEC_b(:,2)
01193 VEC_b(:,5) = VEC_b(:,2)-VEC_b(:,3)
01194 VEC_b(:,6) = VEC_b(:,3)-VEC_b(:,1)
01195 VEC_b(:,7) =-VEC_b(:,1)
01196 VEC_b(:,8) =-VEC_b(:,2)
01197 VEC_b(:,9) =-VEC_b(:,3)
01198 VEC_b(:,10)=-VEC_b(:,4)
01199 VEC_b(:,11)=-VEC_b(:,5)
01200 VEC_b(:,12)=-VEC_b(:,6)
01201 allocate(VEC_d(NBh));VEC_d(:)=0.0d0
01202 allocate(bb(NBh,NBh));bb(:,:)=0.0d0
01203 do ib=1,Nbh
01204 ab=0
01205 do ix=1,3
01206 do jx=ix,3
01207 ab=ab+1
01208 bb(ib,ab)=VEC_b(ix,ib)*VEC_b(jx,ib)
01209 enddo
01210 enddo
01211 enddo
01212 call invmat(Nbh,bb(1,1))
01213 VEC_d(1)=0.5d0
01214 VEC_d(2)=0.0d0
01215 VEC_d(3)=0.0d0
01216 VEC_d(4)=0.5d0
01217 VEC_d(5)=0.0d0
01218 VEC_d(6)=0.5d0
01219 WEIGHT_b(:)=0.0D0
01220 do ib=1,Nbh
01221 ab=0
01222 SUM_REAL=0.0d0
01223 do ix=1,3
01224 do jx=ix,3
01225 ab=ab+1
01226 SUM_REAL=SUM_REAL+VEC_d(ab)*bb(ab,ib)
01227 enddo
01228 enddo
01229 WEIGHT_b(ib)=SUM_REAL
01230 enddo
01231 WEIGHT_b(7) = WEIGHT_b(1)
01232 WEIGHT_b(8) = WEIGHT_b(2)
01233 WEIGHT_b(9) = WEIGHT_b(3)
01234 WEIGHT_b(10)= WEIGHT_b(4)
01235 WEIGHT_b(11)= WEIGHT_b(5)
01236 WEIGHT_b(12)= WEIGHT_b(6)
01237 end select
01238
01239 do ibvec=1,NB
01240 write(6,*) b_LAT(:,ibvec)
01241 enddo
01242 write(6,*)
01243 SUM_REAL=0.0D0
01244 Do ibvec=1,NB
01245 SUM_REAL=SUM_REAL+WEIGHT_b(ibvec)
01246 ENDDO
01247 WEIGHT=SUM_REAL
01248
01249 write(6,*)'WEIGHT_b=',WEIGHT_b(:)
01250 write(6,*)'WEIGHT =',WEIGHT
01251
01252
01253
01254
01255
01256 write(6,*)
01257 write(6,*)'SUM RULE CHECK OF B-VEC'
01258 write(6,*)
01259 Do IC=1,3
01260 Do JC=1,3
01261 SUM_REAL=0.0D0
01262 Do ibvec=1,NB
01263 SUM_REAL
01264 + =SUM_REAL+WEIGHT_b(ibvec)*VEC_b(IC,ibvec)*VEC_b(JC,ibvec)
01265 ENDDO
01266 write(6,*) IC,JC,SUM_REAL
01267 ENDDO
01268 ENDDO
01269
01270
01271
01272
01273
01274 WRITE(6,*)
01275 WRITE(6,*)'==================='
01276 WRITE(6,*)' INTERSTATE MATRIX '
01277 WRITE(6,*)'==================='
01278 WRITE(6,*)
01279 allocate(OVERLAP(NTB,NTB,NTK,NB));OVERLAP=0.0D0
01280 allocate(KPT(NTK,NB));KPT=0
01281 nvx=nwx2+1;nvy=nwy2+1;nvz=nwz2+1
01282
01283
01284
01285 allocate(C0_BRA(0:NTG,NTB));C0_BRA(0,:)=0.0D0
01286 allocate(C0_KET(0:NTG,NTB));C0_KET(0,:)=0.0D0
01287 allocate(SHIFT_b(3));SHIFT_b=0
01288 allocate(OVERLAP_TMP(NTB,NTB));OVERLAP_TMP=0.0d0
01289
01290 do ik=1,NTK
01291
01292 C0_BRA(0,:)=0.0D0
01293 do ig=1,NG0(ik)
01294 C0_BRA(ig,:)=C0(ig,:,ik)
01295 enddo
01296
01297 do ibvec=1,NB
01298 call make_shift_vector(NTK,NB,ik,ibvec,SK0(1,1),b_LAT(1,1),
01299 + SHIFT_b(1),ik_ib)
01300 KPT(ik,ibvec)=ik_ib
01301
01302 C0_KET(0,:)=0.0D0
01303 do ig=1,NG0(ik_ib)
01304 C0_KET(ig,:)=C0(ig,:,ik_ib)
01305 enddo
01306
01307 CALL CALC_OVERLAP(NTK,NTB,NTG,nvx,nvy,nvz,ik,ik_ib,NG0(1),
01308 + KG0(1,1,1),N_BAND(1),SHIFT_b(1),C0_BRA(0,1),C0_KET(0,1),
01309 + OVERLAP_TMP(1,1),N_BAND_BTM(1))
01310 OVERLAP(:,:,ik,ibvec)=OVERLAP_TMP(:,:)
01311
01312 enddo
01313 enddo
01314
01315 deallocate(C0_BRA,C0_KET,SHIFT_b,OVERLAP_TMP)
01316
01317 write(6,*)'finish calc ISM'
01318
01319 WRITE(6,*)
01320 write(6,*)'================================================='
01321 write(6,*)' LOG{det[M(k,k+b)]} only when M is square matrix '
01322 write(6,*)'================================================='
01323 WRITE(6,*)
01324 do ik=1,NTK
01325 do ibvec=1,NB
01326 ik_ib=KPT(ik,ibvec)
01327 if(N_BAND(ik)==N_BAND(ik_ib))then
01328 nm=N_BAND(ik)
01329 allocate(mat_tmp(nm,nm));mat_tmp=0.0d0
01330 mat_tmp(1:nm,1:nm)=OVERLAP(1:nm,1:nm,ik,ibvec)
01331 call calcdet(nm,mat_tmp(1,1),DET)
01332 write(6,'(a8,i6,a8,i6,a8,2f15.10)')
01333 + 'ik=',ik,'ik_ib=',ik_ib,'log=',log(DET)
01334 deallocate(mat_tmp)
01335 endif
01336 enddo
01337 enddo
01338
01339
01340
01341
01342
01343 allocate(orbtype(nigs))
01344 allocate(ALPHA_GAUSS(nigs));ALPHA_GAUSS(:)=0.0d0
01345 allocate(TAU_GAUSS(3,nigs));TAU_GAUSS(:,:)=0.0d0
01346 allocate(NORM_GAUSS(nigs));NORM_GAUSS(:)=0.0d0
01347 allocate(loc_x(3,nigs));loc_x=0.0d0
01348 allocate(loc_y(3,nigs));loc_y=0.0d0
01349 allocate(loc_z(3,nigs));loc_z=0.0d0
01350
01351 do ig=1,nigs
01352 orbtype(ig)=vec_ini(ig)%orb
01353 ALPHA_GAUSS(ig)=vec_ini(ig)%a
01354 TAU_GAUSS(1,ig)=vec_ini(ig)%x
01355 TAU_GAUSS(2,ig)=vec_ini(ig)%y
01356 TAU_GAUSS(3,ig)=vec_ini(ig)%z
01357 loc_x(:,ig)=vec_ini(ig)%lx(:)
01358 loc_y(:,ig)=vec_ini(ig)%ly(:)
01359 loc_z(:,ig)=vec_ini(ig)%lz(:)
01360 enddo
01361
01362
01363 allocate(LGAUSS(nigs));LGAUSS(:)=0
01364 allocate(MGAUSS(nigs));MGAUSS(:)=0
01365 do ig=1,nigs
01366 if(orbtype(ig)=='s')then
01367 LGAUSS(ig)=1;MGAUSS(ig)=1
01368 endif
01369 if(orbtype(ig)=='px')then
01370 LGAUSS(ig)=2;MGAUSS(ig)=1
01371 endif
01372 if(orbtype(ig)=='py')then
01373 LGAUSS(ig)=2;MGAUSS(ig)=2
01374 endif
01375 if(orbtype(ig)=='pz')then
01376 LGAUSS(ig)=2;MGAUSS(ig)=3
01377 endif
01378 if(orbtype(ig)=='dxy')then
01379 LGAUSS(ig)=3;MGAUSS(ig)=1
01380 endif
01381 if(orbtype(ig)=='dyz')then
01382 LGAUSS(ig)=3;MGAUSS(ig)=2
01383 endif
01384 if(orbtype(ig)=='dz2')then
01385 LGAUSS(ig)=3;MGAUSS(ig)=3
01386 endif
01387 if(orbtype(ig)=='dzx')then
01388 LGAUSS(ig)=3;MGAUSS(ig)=4
01389 endif
01390 if(orbtype(ig)=='dx2')then
01391 LGAUSS(ig)=3;MGAUSS(ig)=5
01392 endif
01393 enddo
01394
01395
01396 Do ig=1,nigs
01397 iL=LGAUSS(ig)
01398 IF(iL==1)THEN
01399 NORM_GAUSS(ig)
01400 + =(8.0D0*ALPHA_GAUSS(ig)/pi)**0.25D0
01401 + *(4.0D0*ALPHA_GAUSS(ig)/1.0D0)**0.50D0
01402 ENDIF
01403 IF(iL==2)THEN
01404 NORM_GAUSS(ig)
01405 + =(8.0D0*ALPHA_GAUSS(ig)/pi)**0.25D0
01406 + *(4.0D0*ALPHA_GAUSS(ig)/1.0D0)**0.50D0
01407 + *(4.0D0*ALPHA_GAUSS(ig)/3.0D0)**0.50D0
01408 ENDIF
01409 IF(iL==3)THEN
01410 NORM_GAUSS(ig)
01411 + =(8.0D0*ALPHA_GAUSS(ig)/pi)**0.25D0
01412 + *(4.0D0*ALPHA_GAUSS(ig)/1.0D0)**0.50D0
01413 + *(4.0D0*ALPHA_GAUSS(ig)/3.0D0)**0.50D0
01414 + *(4.0D0*ALPHA_GAUSS(ig)/5.0D0)**0.50D0
01415 ENDIF
01416 ENDDO
01417
01418 write(6,*)
01419 write(6,*)'=================='
01420 write(6,*)'INITIAL GUESS DATA'
01421 write(6,*)'=================='
01422 write(6,*)
01423 do ig=1,nigs
01424 write(6,'(2i5,4f15.10)')
01425 + LGAUSS(ig),MGAUSS(ig),ALPHA_GAUSS(ig),TAU_GAUSS(:,ig)
01426 enddo
01427 write(6,*)
01428 write(6,*)'=================='
01429 write(6,*)'COORD IN CARTESIAN'
01430 write(6,*)'=================='
01431 write(6,*)
01432 do ig=1,nigs
01433 tmp(:)=TAU_GAUSS(:,ig)
01434 TAU_GAUSS(:,ig)=tmp(1)*a1(:)+tmp(2)*a2(:)+tmp(3)*a3(:)
01435 write(6,'(3f20.12)') TAU_GAUSS(:,ig)
01436 enddo
01437 write(6,*)
01438
01439
01440 write(6,*)
01441 write(6,*)'=========='
01442 write(6,*)'CALC A_MAT'
01443 write(6,*)'=========='
01444 write(6,*)
01445 allocate(A_MAT(NTB,nigs,NTK));A_MAT=0.0d0
01446
01447
01448 allocate(A_TMP(NTB,nigs));A_TMP=0.0d0
01449
01450 do ik=1,NTK
01451 call make_amat(nigs,NTB,NTG,NG0(ik),SK0(1,ik),KG0(1,1,ik),
01452 + b1(1),b2(1),b3(1),TAU_GAUSS(1,1),ALPHA_GAUSS(1),
01453 + loc_x(1,1),loc_y(1,1),loc_z(1,1),LGAUSS(1),MGAUSS(1),
01454 + N_BAND(ik),C0(1,1,ik),N_BAND_BTM(ik),NORM_GAUSS(1),VOLUME,
01455 + A_TMP(1,1))
01456 A_MAT(:,:,ik)=A_TMP(:,:)
01457 enddo
01458
01459 deallocate(A_TMP)
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530 allocate(A_TMP(NTB,nigs));A_TMP=0.0d0
01531 do ik=1,NTK
01532 A_TMP(:,:)=A_MAT(:,:,ik)
01533 do ib=1,N_BAND(ik)
01534 do ig=1,NGAUSS
01535 SUM_CMPX=0.0D0
01536 do jg=1,nigs
01537 SUM_CMPX=SUM_CMPX+A_TMP(ib,jg)*B_MAT(jg,ig)
01538 ENDDO
01539 A_MAT(ib,ig,ik)=SUM_CMPX
01540 ENDDO
01541 ENDDO
01542 ENDDO
01543 write(6,*)'finish calc A_MAT'
01544
01545 allocate(S_MAT(NGAUSS,NGAUSS));S_MAT(:,:)=0.0D0
01546 allocate(S_TMP(NGAUSS,NGAUSS));S_TMP(:,:)=0.0d0
01547 allocate(S_EIG(NGAUSS));S_EIG(:)=0.0D0
01548 allocate(X_MAT(NGAUSS,NGAUSS));X_MAT(:,:)=0.0D0
01549 allocate(S_HALF_MAT(NGAUSS,NGAUSS));S_HALF_MAT(:,:)=0.0D0
01550 allocate(C_MAT(NTB,NGAUSS,NTK));C_MAT(:,:,:)=0.0D0
01551
01552 do ik=1,NTK
01553
01554 S_MAT(:,:)=0.0d0
01555 do ig=1,NGAUSS
01556 do jg=1,NGAUSS
01557 SUM_CMPX=0.0D0
01558 do ib=1,N_BAND(ik)
01559 SUM_CMPX=SUM_CMPX+CONJG(A_MAT(ib,ig,ik))*A_MAT(ib,jg,ik)
01560 enddo
01561 S_MAT(ig,jg)=SUM_CMPX
01562 enddo
01563 enddo
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573 S_EIG(:)=0.0D0
01574 S_TMP(:,:)=S_MAT(:,:)
01575 call diagV(NGAUSS,S_TMP(1,1),S_EIG(1))
01576 X_MAT(:,:)=S_TMP(:,:)
01577 do ig=1,NGAUSS
01578 if(S_EIG(ig)<=0.0D0)then
01579 write(6,*)'WARNING diag S_MAT'
01580 write(6,*) ig,S_EIG(ig)
01581 S_EIG(ig)=1.0D-10
01582 endif
01583 enddo
01584
01585 S_HALF_MAT(:,:)=0.0D0
01586 S_EIG(:)=1.0D0/DSQRT(S_EIG(:))
01587 do ig=1,NGAUSS
01588 do jg=1,NGAUSS
01589 SUM_CMPX=0.0D0
01590 do kg=1,NGAUSS
01591 SUM_CMPX=SUM_CMPX+X_MAT(ig,kg)*S_EIG(kg)*CONJG(X_MAT(jg,kg))
01592 enddo
01593 S_HALF_MAT(ig,jg)=SUM_CMPX
01594 enddo
01595 enddo
01596
01597 do ib=1,N_BAND(ik)
01598 do jg=1,NGAUSS
01599 SUM_CMPX=0.0D0
01600 do kg=1,NGAUSS
01601 SUM_CMPX=SUM_CMPX+A_MAT(ib,kg,ik)*S_HALF_MAT(kg,jg)
01602 enddo
01603 C_MAT(ib,jg,ik)=SUM_CMPX
01604 enddo
01605 enddo
01606 enddo
01607 deallocate(S_MAT,S_TMP,S_EIG,X_MAT,S_HALF_MAT)
01608
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638 if(set_inner_window)then
01639 allocate(N_BAND_inner(NTK));N_BAND_inner=0
01640 allocate(N_BAND_BTM_inner(NTK));N_BAND_BTM_inner=0
01641 allocate(inner(NTB,NTK));inner=0
01642 write(6,*)
01643 write(6,*)'===================================='
01644 write(6,*)'BAND INFO INSIDE INNER ENERGY WINDOW'
01645 write(6,*)'===================================='
01646 write(6,*)
01647 do ik=1,NTK
01648 SUM_INT=0
01649 do ib=1,N_BAND(ik)
01650 if(E_LOWER_inner<=E_EIG(N_BAND_BTM(ik)+ib,ik).and.
01651 + E_EIG(N_BAND_BTM(ik)+ib,ik)<=E_UPPER_inner)then
01652 inner(ib,ik)=1
01653 SUM_INT=SUM_INT+1
01654 if(SUM_INT==1) N_BAND_BTM_inner(ik)=ib-1
01655 endif
01656 enddo
01657 N_BAND_inner(ik)=SUM_INT
01658
01659
01660 if(N_BAND_inner(ik)>N_wannier)then
01661 write(6,'(a)')'inner energy window setting: WRONG'
01662 write(6,'(a)')'NBAND_inner SHOULD BE SMALLER THAN N_wannier'
01663 write(6,'(a,3i8)')'ik, N_BAND_inner(ik), N_wannier',
01664 + ik, N_BAND_inner(ik), N_wannier
01665 stop
01666 endif
01667 write(6,*) ik,N_BAND_BTM_inner(ik),N_BAND_inner(ik)
01668 enddo
01669 write(6,*)
01670
01671
01672 do ik=1,NTK
01673
01674 allocate(P_MAT(NTB,NTB));P_MAT=0.0d0
01675 do ib=1,N_BAND(ik)
01676 do jb=1,N_BAND(ik)
01677 SUM_CMPX=0.0d0
01678 do jg=1,NGAUSS
01679 SUM_CMPX=SUM_CMPX+C_MAT(ib,jg,ik)*CONJG(C_MAT(jb,jg,ik))
01680 enddo
01681 P_MAT(ib,jb)=SUM_CMPX
01682 enddo
01683 enddo
01684 allocate(Qin(NTB,NTB));Qin=0.0d0
01685 do ib=1,N_BAND(ik)
01686 if(inner(ib,ik)==1)cycle
01687 Qin(ib,ib)=1.0d0
01688 enddo
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698 allocate(X_MAT(NTB,NTB));X_MAT=0.0d0
01699 do ib=1,N_BAND(ik)
01700 do jb=1,N_BAND(ik)
01701 SUM_CMPX=0.0d0
01702 do k=1,N_BAND(ik)
01703 do l=1,N_BAND(ik)
01704 SUM_CMPX=SUM_CMPX+CONJG(Qin(k,ib))*P_MAT(k,l)*Qin(l,jb)
01705 enddo
01706 enddo
01707 X_MAT(ib,jb)=SUM_CMPX
01708 enddo
01709 enddo
01710
01711 nm=N_BAND(ik)
01712 allocate(mat_tmp(nm,nm));mat_tmp(:,:)=0.0d0
01713 allocate(eig_tmp(nm));eig_tmp(:)=0.0d0
01714 mat_tmp(1:nm,1:nm)=X_MAT(1:nm,1:nm)
01715 call diagV(nm,mat_tmp(1,1),eig_tmp(1))
01716
01717 X_MAT=0.0d0
01718 do ib=1,nm
01719 X_MAT(1:nm,ib)=mat_tmp(1:nm,nm-ib+1)
01720 enddo
01721 deallocate(mat_tmp,eig_tmp)
01722
01723 do jg=1,NGAUSS-N_BAND_inner(ik)
01724 C_MAT(:,jg,ik)=X_MAT(:,jg)
01725 enddo
01726 do jg=1,N_BAND_inner(ik)
01727 C_MAT(:,NGAUSS-N_BAND_inner(ik)+jg,ik)=0.0d0
01728 enddo
01729 do ib=1,N_BAND_inner(ik)
01730 i=N_BAND_BTM_inner(ik)+ib
01731 j=NGAUSS-N_BAND_inner(ik)+ib
01732 C_MAT(i,j,ik)=1.0d0
01733 enddo
01734 deallocate(P_MAT,Qin,X_MAT)
01735 enddo
01736 else
01737 allocate(N_BAND_inner(NTK));N_BAND_inner=0
01738 allocate(N_BAND_BTM_inner(NTK));N_BAND_BTM_inner=0
01739 allocate(inner(NTB,NTK));inner=0
01740 endif
01741
01742
01743
01744
01745
01746
01747 WRITE(6,*)
01748 WRITE(6,*)'==========================================='
01749 WRITE(6,*)'=== MINIMIZATION OF SPILLAGE FUNCTIONAL ==='
01750 WRITE(6,*)'==========================================='
01751 WRITE(6,*)
01752
01753 allocate(Z_EIG(NTB));Z_EIG(:)=0.0D0
01754 allocate(Z_MAT(NTB,NTB,NTK));Z_MAT(:,:,:)=0.0D0
01755 allocate(Z_TMP(NTB,NTB));Z_TMP(:,:)=0.0D0
01756 allocate(C_TMP(NTB,NTB));C_TMP(:,:)=0.0D0
01757 allocate(M_MAT(NTB,NTB,NTK,NB));M_MAT(:,:,:,:)=0.0D0
01758
01759
01760
01761
01762
01763
01764 I_SCF=0
01765 DEL_OMEGA_I=1.0D0
01766 OMEGA_I_OLD=0.0D0
01767 do while(DEL_OMEGA_I>EPS_SPILLAGE)
01768
01769 M_MAT(:,:,:,:)=0.0D0
01770 do ik=1,NTK
01771 do ibvec=1,NB
01772 ik_ib=KPT(ik,ibvec)
01773 do ib=1,N_BAND(ik)
01774 do jg=1,NGAUSS
01775 SUM_CMPX=0.0D0
01776 do kb=1,N_BAND(ik_ib)
01777 SUM_CMPX=SUM_CMPX+OVERLAP(ib,kb,ik,ibvec)*C_MAT(kb,jg,ik_ib)
01778 enddo
01779 M_MAT(ib,jg,ik,ibvec)=SUM_CMPX
01780 enddo
01781 enddo
01782 enddo
01783 enddo
01784
01785 SUM_REAL=0.0D0
01786 OMEGA_I_NEW=0.0D0
01787 C_MAT(:,:,:)=0.0D0
01788 do ik=1,NTK
01789
01790 Z_TMP(:,:)=0.0D0
01791 do ib=1,N_BAND(ik)
01792 do jb=1,N_BAND(ik)
01793 SUM_CMPX=0.0D0
01794 do ibvec=1,NB
01795 do jg=1,NGAUSS
01796 SUM_CMPX
01797 + =SUM_CMPX
01798 + +WEIGHT_b(ibvec)
01799 + *M_MAT(ib,jg,ik,ibvec)*CONJG(M_MAT(jb,jg,ik,ibvec))
01800 enddo
01801 enddo
01802 Z_TMP(ib,jb)=SUM_CMPX
01803 enddo
01804 enddo
01805
01806 IF(I_SCF==0) THEN
01807 Z_MAT(:,:,ik)=Z_TMP(:,:)
01808 ELSE
01809 Z_MAT(:,:,ik)=DAMP*Z_TMP(:,:)+(1.0D0-DAMP)*Z_MAT(:,:,ik)
01810 ENDIF
01811
01812 nm=N_BAND(ik)-N_BAND_inner(ik)
01813 allocate(mat_tmp(nm,nm));mat_tmp(:,:)=0.0d0
01814 allocate(eig_tmp(nm));eig_tmp(:)=0.0d0
01815 i=0
01816 do ib=1,N_BAND(ik)
01817 if(inner(ib,ik)==1)cycle
01818 i=i+1
01819 j=0
01820 do jb=1,N_BAND(ik)
01821 if(inner(jb,ik)==1)cycle
01822 j=j+1
01823 mat_tmp(i,j)=Z_MAT(ib,jb,ik)
01824 enddo
01825 enddo
01826
01827 eig_tmp(:)=0.0d0
01828
01829
01830
01831 if(nm>0)then
01832 call diagV(nm,mat_tmp(1,1),eig_tmp(1))
01833 endif
01834
01835
01836 Z_EIG(:)=0.0D0
01837 C_TMP(:,:)=0.0d0
01838 do ib=1,nm
01839 Z_EIG(ib)=eig_tmp(nm-ib+1)
01840 C_TMP(1:nm,ib)=mat_tmp(1:nm,nm-ib+1)
01841 enddo
01842 deallocate(mat_tmp,eig_tmp)
01843
01844 i=0
01845 do ib=1,N_BAND(ik)
01846 if(inner(ib,ik)==0)then
01847 i=i+1
01848
01849 C_MAT(ib,1:NGAUSS-N_BAND_inner(ik),ik)
01850 + =C_TMP(i ,1:NGAUSS-N_BAND_inner(ik))
01851 else
01852 C_MAT(ib,:,ik)=0.0d0
01853 endif
01854 enddo
01855 do jg=1,N_BAND_inner(ik)
01856 C_MAT(:,NGAUSS-N_BAND_inner(ik)+jg,ik)=0.0d0
01857 enddo
01858 do ib=1,N_BAND_inner(ik)
01859 i=N_BAND_BTM_inner(ik)+ib
01860 j=NGAUSS-N_BAND_inner(ik)+ib
01861 C_MAT(i,j,ik)=1.0d0
01862 enddo
01863
01864
01865 do ig=1,NGAUSS
01866 do jb=1,N_BAND(ik)
01867 do ib=1,N_BAND(ik)
01868 SUM_REAL
01869 + =SUM_REAL
01870 + +dble(conjg(C_MAT(ib,ig,ik))*Z_MAT(ib,jb,ik)*C_MAT(jb,ig,ik))
01871 enddo
01872 enddo
01873 enddo
01874 OMEGA_I_NEW=OMEGA_I_NEW+dble(NGAUSS)*WEIGHT
01875
01876
01877
01878
01879
01880
01881
01882 enddo
01883 OMEGA_I_NEW=(OMEGA_I_NEW-SUM_REAL)/DBLE(NTK)
01884 DEL_OMEGA_I=DABS(OMEGA_I_NEW-OMEGA_I_OLD)
01885 OMEGA_I_OLD=OMEGA_I_NEW
01886 I_SCF=I_SCF+1
01887 write(6,'(A8,I5,A15,F20.10,A10,F20.10)')'I_SCF=',I_SCF,
01888 + 'DEL_OMEGA_I=',DEL_OMEGA_I,'OMEGA_I',OMEGA_I_NEW
01889 enddo
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
01900
01901
01902
01903
01904 allocate(PSEUDO_EIG(NGAUSS,NTK));PSEUDO_EIG(:,:)=0.0D0
01905 allocate(PSEUDO_MAT(NGAUSS,NGAUSS,NTK));PSEUDO_MAT(:,:,:)=0.0D0
01906 allocate(P_EIG(NGAUSS));P_EIG(:)=0.0d0
01907 allocate(P_MAT_IN(NGAUSS,NGAUSS));P_MAT_IN(:,:)=0.0D0
01908 allocate(P_MAT_OUT(NGAUSS,NGAUSS));P_MAT_OUT(:,:)=0.0D0
01909
01910 do ik=1,NTK
01911
01912 P_MAT_IN(:,:)=0.0D0
01913 do ig=1,NGAUSS
01914 do jg=1,NGAUSS
01915 SUM_CMPX=0.0D0
01916 do kb=1,N_BAND(ik)
01917 SUM_CMPX
01918 + =SUM_CMPX
01919 + +CONJG(C_MAT(kb,ig,ik))
01920 + *E_EIG(kb+N_BAND_BTM(ik),ik)
01921 + *C_MAT(kb,jg,ik)
01922 enddo
01923 P_MAT_IN(ig,jg)=SUM_CMPX
01924 enddo
01925 enddo
01926
01927 P_EIG(:)=0.0D0
01928 call diagV(NGAUSS,P_MAT_IN(1,1),P_EIG(1))
01929 P_MAT_OUT(:,:)=P_MAT_IN(:,:)
01930 PSEUDO_MAT(:,:,ik)=P_MAT_OUT(:,:)
01931 PSEUDO_EIG(:,ik)=P_EIG(:)
01932 enddo
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943 allocate(S_MAT(NGAUSS,NGAUSS));S_MAT(:,:)=0.0D0
01944 allocate(S_TMP(NGAUSS,NGAUSS));S_TMP(:,:)=0.0d0
01945 allocate(S_EIG(NGAUSS));S_EIG(:)=0.0D0
01946 allocate(X_MAT(NGAUSS,NGAUSS));X_MAT(:,:)=0.0D0
01947 allocate(S_HALF_MAT(NGAUSS,NGAUSS));S_HALF_MAT(:,:)=0.0D0
01948 Z_MAT(:,:,:)=0.0D0
01949 do ik=1,NTK
01950
01951 A_TMP(:,:)=0.0D0
01952 do ib=1,n_occ
01953 do jg=1,NGAUSS
01954 SUM_CMPX=0.0D0
01955 do kb=1,N_BAND(ik)
01956 do lg=1,NGAUSS
01957 SUM_CMPX
01958 + =SUM_CMPX
01959 + +CONJG(PSEUDO_MAT(lg,ib,ik))
01960 + *CONJG(C_MAT(kb,lg,ik))
01961 + *A_MAT(kb,jg,ik)
01962 ENDDO
01963 ENDDO
01964 A_TMP(ib,jg)=SUM_CMPX
01965 ENDDO
01966 ENDDO
01967
01968 S_MAT(:,:)=0.0D0
01969 Do ig=1,NGAUSS
01970 Do jg=1,NGAUSS
01971 SUM_CMPX=0.0D0
01972 Do kb=1,n_occ
01973 SUM_CMPX
01974 + =SUM_CMPX
01975 + +CONJG(A_TMP(kb,ig))
01976 + *A_TMP(kb,jg)
01977 ENDDO
01978 S_MAT(ig,jg)=SUM_CMPX
01979 ENDDO
01980 ENDDO
01981
01982 S_EIG(:)=0.0D0
01983 call diagV(NGAUSS,S_MAT(1,1),S_EIG(1))
01984 X_MAT(:,:)=S_MAT(:,:)
01985 Do ig=1,NGAUSS
01986 IF(S_EIG(ig)<=0.0D0) THEN
01987 write(6,*) ig,S_EIG(ig)
01988 S_EIG(ig)=1.0D-10
01989 ENDIF
01990 ENDDO
01991
01992 S_EIG(:)=1.0D0/DSQRT(S_EIG(:))
01993 S_HALF_MAT(:,:)=0.0D0
01994 Do ig=1,NGAUSS
01995 Do jg=1,NGAUSS
01996 SUM_CMPX=0.0D0
01997 Do kg=1,NGAUSS
01998 SUM_CMPX
01999 + =SUM_CMPX
02000 + +X_MAT(ig,kg)
02001 + *S_EIG(kg)
02002 + *CONJG(X_MAT(jg,kg))
02003 ENDDO
02004 S_HALF_MAT(ig,jg)=SUM_CMPX
02005 ENDDO
02006 ENDDO
02007
02008 A_MAT(:,:,ik)=0.0D0
02009 Do ib=1,n_occ
02010 Do jg=1,NGAUSS
02011 SUM_CMPX=0.0D0
02012 Do kg=1,NGAUSS
02013 SUM_CMPX
02014 + =SUM_CMPX
02015 + +A_TMP(ib,kg)
02016 + *S_HALF_MAT(kg,jg)
02017 ENDDO
02018 A_MAT(ib,jg,ik)=SUM_CMPX
02019 ENDDO
02020 ENDDO
02021
02022 Do ib=1,N_BAND(ik)
02023 Do jg=1,NGAUSS
02024 SUM_CMPX=0.0D0
02025 Do kg=1,NGAUSS
02026 Do llb=1,n_occ
02027 SUM_CMPX
02028 + =SUM_CMPX
02029 + +C_MAT(ib,kg,ik)
02030 + *PSEUDO_MAT(kg,llb,ik)
02031 + *A_MAT(llb,jg,ik)
02032 ENDDO
02033 ENDDO
02034 Z_MAT(ib,jg,ik)=SUM_CMPX
02035 ENDDO
02036 ENDDO
02037 ENDDO
02038 deallocate(S_MAT,S_TMP,S_EIG,X_MAT,S_HALF_MAT)
02039
02040
02041
02042
02043
02044
02045
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065 WRITE(6,*)
02066 WRITE(6,*)'====================================='
02067 WRITE(6,*)'INTERSTATE MATRIX; LOG{det[M(k,k+b)]}'
02068 WRITE(6,*)'====================================='
02069 WRITE(6,*)
02070
02071 M_MAT(:,:,:,:)=OVERLAP(:,:,:,:)
02072 OVERLAP(:,:,:,:)=0.0D0
02073 Do ik=1,NTK
02074 Do ibvec=1,NB
02075 ik_ib=KPT(ik,ibvec)
02076 Do i_band=1,n_occ
02077 Do j_band=1,n_occ
02078 SUM_CMPX=0.0D0
02079 Do k_band=1,N_BAND(ik)
02080 Do l_band=1,N_BAND(ik_ib)
02081 SUM_CMPX
02082 + =SUM_CMPX
02083 + +CONJG(Z_MAT(k_band,i_band,ik))
02084 + * M_MAT(k_band,l_band,ik,ibvec)
02085 + * Z_MAT(l_band,j_band,ik_ib)
02086 ENDDO
02087 ENDDO
02088 OVERLAP(i_band,j_band,ik,ibvec)=SUM_CMPX
02089 ENDDO
02090 ENDDO
02091
02092
02093 nm=n_occ
02094 allocate(mat_tmp(nm,nm));mat_tmp(:,:)=0.0d0
02095 mat_tmp(1:nm,1:nm)=OVERLAP(1:nm,1:nm,ik,ibvec)
02096 call calcdet(nm,mat_tmp(1,1),DET)
02097 deallocate(mat_tmp)
02098
02099 write(6,'(a8,i5,a8,i5,a8,2f15.10)')
02100 + 'IK=',ik,'IK_IB=',ik_ib,'LOG=',LOG(DET)
02101 ENDDO
02102 WRITE(6,*)
02103 ENDDO
02104
02105
02106
02107
02108
02109 allocate(WF_CHARGE(n_occ,NTK,NB));WF_CHARGE(:,:,:)=0.0d0
02110 allocate(TR_GRAD(NTK));TR_GRAD(:)=0.0d0
02111 allocate(TR_GRAD_OLD(NTK));TR_GRAD_OLD(:)=0.0d0
02112 allocate(GRADIENT(n_occ,n_occ,NTK));GRADIENT(:,:,:)=0.0d0
02113 allocate(DIRECTION(n_occ,n_occ,NTK));DIRECTION(:,:,:)=0.0d0
02114 allocate(DIRECTION2(n_occ,n_occ,NTK));DIRECTION2(:,:,:)=0.0d0
02115 allocate(UNITARY(n_occ,n_occ,NTK));UNITARY(:,:,:)=0.0d0
02116 allocate(U_ORG(n_occ,n_occ,NTK));U_ORG(:,:,:)=0.0d0
02117 allocate(U_OLD(n_occ,n_occ,NTK));U_OLD(:,:,:)=0.0d0
02118
02119 ILS=1
02120 CALL CALC_OMEGA(NTK,NTB,n_occ,NB,VEC_b,ILS,WEIGHT_b,OVERLAP,
02121 + WF_CHARGE,OMEGA_I,OMEGA_OD,OMEGA_D)
02122
02123 write(6,*)'OMEGA_I =',OMEGA_I
02124 write(6,*)'OMEGA_OD=',OMEGA_OD
02125 write(6,*)'OMEGA_D =',OMEGA_D
02126 write(6,*)'OMEGA_total=',OMEGA_I+OMEGA_D+OMEGA_OD
02127 WRITE(6,*)
02128 WRITE(6,*)'========================================='
02129 WRITE(6,*)'=== MINIMIZATION OF SPREAD FUNCTIONAL ==='
02130 WRITE(6,*)'========================================='
02131 WRITE(6,*)
02132
02133 M_MAT(:,:,:,:)=OVERLAP(:,:,:,:)
02134 DEL_SPREAD=1.0d0
02135 SPREAD_OLD=0.0D0
02136 I_STEP=1
02137 SPREAD=OMEGA_I+OMEGA_OD+OMEGA_D
02138 do while(DEL_SPREAD>EPS_SPREAD)
02139
02140
02141 ILS=0
02142 DELTA(1)=0.0D0
02143 OMEGA(1)=SPREAD
02144 write(6,*)'ILS=',ILS,DELTA(1),OMEGA(1)
02145
02146
02147 CALL CALC_GRADIENT(NTK,NTB,n_occ,NB,WF_CHARGE,WEIGHT_b,OVERLAP,
02148 + tci,GRADIENT,TR_GRAD)
02149 DIRECTION(:,:,:)=GRADIENT(:,:,:)
02150
02151
02152
02153
02154
02155
02156
02157
02158
02159
02160
02161
02162
02163
02164
02165 DIRECTION2(:,:,:)=0.0D0
02166 Do ik=1,NTK
02167 Do i_band=1,n_occ
02168 Do j_band=1,n_occ
02169 SUM_CMPX=0.0D0
02170 Do k_band=1,n_occ
02171 SUM_CMPX
02172 + =SUM_CMPX
02173 + +DIRECTION(i_band,k_band,ik)
02174 + *DIRECTION(k_band,j_band,ik)
02175 ENDDO
02176 DIRECTION2(i_band,j_band,ik)=SUM_CMPX
02177 ENDDO
02178 ENDDO
02179 ENDDO
02180
02181
02182
02183 ILS=1
02184 STEP_LENGTH=MAX_STEP_LENGTH/4.0D0/WEIGHT
02185 DELTA(2)=STEP_LENGTH
02186 CALL UNITARY_GEN(NTK,n_occ,ILS,DIRECTION(1,1,1),DIRECTION2(1,1,1),
02187 + STEP_LENGTH,UNITARY(1,1,1),U_ORG(1,1,1),U_OLD(1,1,1))
02188 OVERLAP(:,:,:,:)=M_MAT(:,:,:,:)
02189 CALL UPDATE_OVERLAP(NTK,NTB,n_occ,NB,KPT,UNITARY,OVERLAP)
02190 CALL CALC_OMEGA(NTK,NTB,n_occ,NB,VEC_b,ILS,
02191 + WEIGHT_b,OVERLAP,WF_CHARGE,OMEGA_I,OMEGA_OD,OMEGA_D)
02192 OMEGA(2)=OMEGA_I+OMEGA_D+OMEGA_OD
02193 write(6,*)'ILS=',ILS,DELTA(2),OMEGA(2)
02194
02195 if(abs(OMEGA(1)-OMEGA(2))<1.0d-10)then
02196
02197 M_MAT(:,:,:,:)=OVERLAP(:,:,:,:)
02198
02199 CALL A_MAT_UPDATE(NTK,NTB,NGAUSS,n_occ,UNITARY,A_MAT)
02200
02201 SPREAD_OLD=OMEGA(1)
02202 SPREAD=OMEGA(2)
02203 DEL_SPREAD=ABS(SPREAD-SPREAD_OLD)
02204 SPREAD_OLD=OMEGA(2)
02205 write(6,*)
02206 write(6,'(a47)')'=============================================='
02207 write(6,'(a47)')' SPREAD HAS ALREADY CONVERGED IN INITIAL STEP '
02208 write(6,'(a47)')'=============================================='
02209 write(6,*)
02210 write(6,'(a26,i8,2f20.10)')'I_STEP SPREAD DEL_SPREAD:',
02211 + I_STEP,SPREAD,DEL_SPREAD
02212 write(6,*)
02213 I_STEP=I_STEP+1
02214 else
02215 do while(OMEGA(1)>OMEGA(2))
02216 ILS=ILS+1
02217 DELTA(3) = DELTA(2)+STEP_LENGTH
02218 CALL UNITARY_GEN(NTK,n_occ,ILS,DIRECTION(1,1,1),
02219 + DIRECTION2(1,1,1),STEP_LENGTH,UNITARY(1,1,1),
02220 + U_ORG(1,1,1),U_OLD(1,1,1))
02221 OVERLAP(:,:,:,:)=M_MAT(:,:,:,:)
02222 CALL UPDATE_OVERLAP(NTK,NTB,n_occ,NB,
02223 + KPT,UNITARY,OVERLAP)
02224 CALL CALC_OMEGA(NTK,NTB,n_occ,NB,VEC_b,ILS,
02225 + WEIGHT_b,OVERLAP,WF_CHARGE,
02226 + OMEGA_I,OMEGA_OD,OMEGA_D)
02227 OMEGA(3)=OMEGA_I+OMEGA_D+OMEGA_OD
02228 write(6,*) 'ILS=',ILS,DELTA(3),OMEGA(3)
02229 IF(OMEGA(2)<OMEGA(3)) EXIT
02230 DELTA(1)=DELTA(2)
02231 OMEGA(1)=OMEGA(2)
02232 DELTA(2)=DELTA(3)
02233 OMEGA(2)=OMEGA(3)
02234 enddo
02235
02236 DELTA(4)=0.5D0*(OMEGA(1)*(DELTA(2)**2-DELTA(3)**2)
02237 + +OMEGA(2)*(DELTA(3)**2-DELTA(1)**2)
02238 + +OMEGA(3)*(DELTA(1)**2-DELTA(2)**2))
02239 + /(OMEGA(1)*(DELTA(2)-DELTA(3))
02240 + +OMEGA(2)*(DELTA(3)-DELTA(1))
02241 + +OMEGA(3)*(DELTA(1)-DELTA(2)))
02242 ILS=1
02243 STEP_LENGTH=DELTA(4)
02244 CALL UNITARY_GEN(NTK,n_occ,ILS,DIRECTION(1,1,1),
02245 + DIRECTION2(1,1,1),STEP_LENGTH,UNITARY(1,1,1),
02246 + U_ORG(1,1,1),U_OLD(1,1,1))
02247 OVERLAP(:,:,:,:) = M_MAT(:,:,:,:)
02248 CALL UPDATE_OVERLAP(NTK,NTB,n_occ,NB,KPT,UNITARY,OVERLAP)
02249 CALL CALC_OMEGA(NTK,NTB,n_occ,NB,VEC_b,ILS,
02250 + WEIGHT_b,OVERLAP,WF_CHARGE,OMEGA_I,OMEGA_OD,OMEGA_D)
02251 OMEGA(4) = OMEGA_I+OMEGA_D+OMEGA_OD
02252 write(6,*)'DELTA_OPT=',DELTA(4),'OMEGA_OPT=',OMEGA(4)
02253 write(6,*)'OMEGA_I =',OMEGA_I
02254 write(6,*)'OMEGA_OD=',OMEGA_OD
02255 write(6,*)'OMEGA_D =',OMEGA_D
02256
02257 M_MAT(:,:,:,:)=OVERLAP(:,:,:,:)
02258
02259 CALL A_MAT_UPDATE(NTK,NTB,NGAUSS,n_occ,UNITARY,A_MAT)
02260
02261 SPREAD=OMEGA(4)
02262 DEL_SPREAD=ABS(SPREAD-SPREAD_OLD)
02263 SPREAD_OLD=OMEGA(4)
02264 write(6,'(a26,i8,2f20.10)')'I_STEP SPREAD DEL_SPREAD:',
02265 + I_STEP,SPREAD,DEL_SPREAD
02266 write(6,*)
02267 I_STEP=I_STEP+1
02268 endif
02269
02270 ENDDO
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291 allocate(wannier_center(3,n_occ));wannier_center(:,:)=0.0d0
02292 wannier_center(:,:)=0.0D0
02293 Do i_band=1,n_occ
02294 tmp(:)=0.0D0
02295 Do ik=1,NTK
02296 Do ibvec=1,NB
02297 tmp(:)
02298 + =tmp(:)
02299 + +WEIGHT_b(ibvec)
02300 + *VEC_b(:,ibvec)
02301 + *AIMAG(LOG(OVERLAP(i_band,i_band,ik,ibvec)))
02302 ENDDO
02303 ENDDO
02304 wannier_center(:,i_band)=-tmp(:)/DBLE(NTK)
02305 ENDDO
02306
02307
02308 OPEN(134,FILE='./dir-wan/dat.wan-center')
02309 rewind(134)
02310 write(134,'(a)')'#Wannier center'
02311 write(134,'(a)')'#1:x, 2:y, 3:z (in xyz coord)'
02312 do i_band=1,n_occ
02313 write(134,"(3F20.10)")(wannier_center(i,i_band),i=1,3)
02314 enddo
02315
02316
02317 OPEN(113,FILE='./dir-wan/dat.wan',FORM='unformatted')
02318 allocate(C0_TMP_1(NTG,NTB))
02319 allocate(C0_TMP_2(NTG,NTB))
02320 REWIND(113)
02321 write(113) n_occ
02322 do ik=1,NTK
02323 A_TMP(:,:)=0.0D0
02324 do ib=1,N_BAND(ik)
02325 do jb=1,n_occ
02326 SUM_CMPX=0.0D0
02327 do k=1,n_occ
02328 do l=1,n_occ
02329 SUM_CMPX
02330 + =SUM_CMPX
02331 + +C_MAT(ib,k,ik)*PSEUDO_MAT(k,l,ik)*A_MAT(l,jb,ik)
02332 enddo
02333 enddo
02334 A_TMP(ib,jb)=SUM_CMPX
02335 enddo
02336 enddo
02337 C0_TMP_1(:,:)=C0(:,:,ik)
02338 C0_TMP_2(:,:)=0.0D0
02339
02340
02341
02342
02343
02344
02345
02346 do ig=1,NG0(ik)
02347 do jb=1,n_occ
02348 SUM_CMPX=0.0D0
02349 do kb=1,N_BAND(ik)
02350 SUM_CMPX
02351 + =SUM_CMPX
02352 + +C0_TMP_1(ig,kb+N_BAND_BTM(ik))*A_TMP(kb,jb)
02353 enddo
02354 C0_TMP_2(ig,jb)=SUM_CMPX
02355 enddo
02356 enddo
02357 write(113)((C0_TMP_2(ig,ib),ig=1,NG0(ik)),ib=1,n_occ)
02358
02359
02360
02361
02362
02363
02364 enddo
02365
02366
02367 ierr=CHDIR("./dir-wan")
02368 call system('pwd')
02369 REWIND(113)
02370 read(113) NWF
02371 allocate(C0WN(NTG,NWF,NTK));C0WN(:,:,:)=0.0D0
02372 do ik=1,NTK
02373 read(113)((C0WN(ig,jb,ik),ig=1,NG0(ik)),jb=1,NWF)
02374
02375
02376
02377
02378
02379
02380 enddo
02381 write(6,*)'FINISH REDING C0_WN'
02382 ierr=CHDIR("..")
02383 call system('pwd')
02384
02385
02386 OPEN(150,FILE='./dir-wan/dat.umat')
02387 rewind(150)
02388 Mb=maxval(N_BAND)
02389 write(6,*)"Mb=",Mb
02390 write(150,*) n_occ
02391 allocate(UNT(Mb,n_occ,NTK));UNT(:,:,:)=0.0d0
02392 do ik=1,NTK
02393 do ib=1,N_BAND(ik)
02394 do jb=1,n_occ
02395 SUM_CMPX=0.0d0
02396 do ig=1,NG0(ik)
02397 SUM_CMPX
02398 + =SUM_CMPX
02399 + +conjg(C0(ig,ib+N_BAND_BTM(ik),ik))*C0WN(ig,jb,ik)
02400 enddo
02401 UNT(ib,jb,ik)=SUM_CMPX
02402 enddo
02403 enddo
02404 enddo
02405 deallocate(C0WN)
02406 do ik=1,NTK
02407 do ib=1,N_BAND(ik)
02408 write(150,*)(UNT(ib,jb,ik),jb=1,n_occ)
02409 enddo
02410 enddo
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426
02427
02428
02429
02430
02431
02432
02433
02434
02435
02436 WRITE(6,*)
02437 WRITE(6,*)'============================='
02438 WRITE(6,*)'UNITALITY CHECK OF UMAT:<a|b>'
02439 WRITE(6,*)'============================='
02440 WRITE(6,*)
02441 do ik=1,NTK
02442 nm=N_BAND(ik)
02443 allocate(mat_tmp(nm,nm));mat_tmp=0.0d0
02444 do ib=1,nm
02445 do jb=1,nm
02446 SUM_CMPX=0.0D0
02447 do iw=1,n_occ
02448 SUM_CMPX=SUM_CMPX+CONJG(UNT(ib,iw,ik))*UNT(jb,iw,ik)
02449 enddo
02450 mat_tmp(ib,jb)=SUM_CMPX
02451 enddo
02452 enddo
02453
02454 SUM_CMPX=0.0d0
02455 do ib=1,nm
02456 SUM_CMPX=SUM_CMPX+mat_tmp(ib,ib)
02457 enddo
02458 write(6,'(a8,x,i8,2f15.10)')'ik=',ik,SUM_CMPX
02459 write(6,'(300f10.5)')(mat_tmp(ib,ib),ib=1,nm)
02460 write(6,*)
02461
02462 if(N_BAND_inner(ik)/=0)then
02463 do ib=1,nm
02464 write(6,'(300f10.5)')(mat_tmp(ib,jb),jb=1,nm)
02465 enddo
02466 write(6,*)
02467 endif
02468
02469 deallocate(mat_tmp)
02470 enddo
02471
02472
02473
02474
02475
02476 allocate(H_MAT_K(n_occ,n_occ,NTK));H_MAT_K(:,:,:)=0.0D0
02477 allocate(E_TMP(n_occ));E_TMP(:)=0.0d0
02478 allocate(H_TMP_IN(n_occ,n_occ));H_TMP_IN(:,:)=0.0d0
02479 allocate(H_MAT_R(n_occ,n_occ,-Na1:Na1,-Na2:Na2,-Na3:Na3))
02480 allocate(WEIGHT_R(-Na1:Na1,-Na2:Na2,-Na3:Na3))
02481
02482
02483 NSK_BAND_DISP=Ndiv*(N_sym_points-1)+1
02484
02485
02486
02487 allocate(H_TMP_OUT(n_occ,n_occ,NSK_BAND_DISP));H_TMP_OUT(:,:,:)=0.0d0
02488
02489 allocate(SK_BAND_DISP(3,NSK_BAND_DISP));SK_BAND_DISP(:,:)=0.0d0
02490 allocate(E_BAND_DISP(n_occ,NSK_BAND_DISP));E_BAND_DISP(:,:)=0.0d0
02491 call makekpts(Ndiv,N_sym_points,NSK_BAND_DISP,SK_sym_pts(1,1),
02492 + SK_BAND_DISP(1,1))
02493 write(6,*)
02494 write(6,*)'====================================='
02495 write(6,*)' CALCULATED kpts FOR BAND DISPERSION '
02496 write(6,*)'====================================='
02497 write(6,*)
02498 write(6,'(a20,i15)')'NSK_BAND_DISP=',NSK_BAND_DISP
02499 write(6,*)
02500 do ik=1,NSK_BAND_DISP
02501 write(6,*) SK_BAND_DISP(:,ik)
02502 enddo
02503
02504
02505 H_MAT_K(:,:,:)=0.0D0
02506 Do ik=1,NTK
02507 Do i_band=1,n_occ
02508 Do j_band=1,n_occ
02509 SUM_CMPX=0.0D0
02510 Do k_band=1,n_occ
02511 SUM_CMPX
02512 + =SUM_CMPX
02513 + +CONJG(A_MAT(k_band,i_band,ik))
02514 + *PSEUDO_EIG(k_band,ik)
02515 + *A_MAT(k_band,j_band,ik)
02516 ENDDO
02517 H_MAT_K(i_band,j_band,ik)=SUM_CMPX
02518 ENDDO
02519 ENDDO
02520 ENDDO
02521
02522
02523 H_MAT_R(:,:,:,:,:)=0.0D0
02524 Do ia1=-Na1,Na1
02525 Do ia2=-Na2,Na2
02526 Do ia3=-Na3,Na3
02527 Do i_band=1,n_occ
02528 Do j_band=1,n_occ
02529 SUM_CMPX=0.0D0
02530 Do ik=1,NTK
02531 PHASE=tpi*(SK0(1,ik)*DBLE(ia1)
02532 + +SK0(2,ik)*DBLE(ia2)
02533 + +SK0(3,ik)*DBLE(ia3))
02534 PHASE_FACTOR=EXP(-ci*PHASE)
02535 SUM_CMPX
02536 + =SUM_CMPX
02537 + +H_MAT_K(i_band,j_band,ik)
02538 + *PHASE_FACTOR
02539 ENDDO
02540 H_MAT_R(i_band,j_band,ia1,ia2,ia3)=SUM_CMPX/DBLE(NTK)
02541 ENDDO
02542 ENDDO
02543 ENDDO
02544 ENDDO
02545 ENDDO
02546
02547
02548
02549
02550
02551 OPEN(122,FILE='./dir-wan/dat.h_mat_r')
02552 REWIND(122)
02553 write(122,'(a)')'#transfer'
02554 write(122,'(a)')'#1:R1, 2:R2, 3:R3 (lattice vector)'
02555 write(122,'(a)')'#1:i, 2:j, 3:Re(t_ij) [eV], 4:Im(t_ij) [eV]'
02556 do ia1=-Na1,Na1
02557 do ia2=-Na2,Na2
02558 do ia3=-Na3,Na3
02559 write(122,*) ia1,ia2,ia3
02560 do ib=1,n_occ
02561 do jb=1,n_occ
02562 write(122,'(i5,i5,2f20.10)')
02563 + ib,jb,H_MAT_R(ib,jb,ia1,ia2,ia3)*au
02564 enddo
02565 enddo
02566 write(122,*)
02567 enddo
02568 enddo
02569 enddo
02570 CLOSE(122)
02571
02572 write(6,*)
02573 write(6,*)'============='
02574 write(6,*)'H_MAT_R in eV'
02575 write(6,*)'============='
02576 write(6,*)
02577 do ia1=-Na1,Na1
02578 do ia2=-Na2,Na2
02579 do ia3=-Na3,Na3
02580 write(6,*) ia1,ia2,ia3
02581 do ib=1,n_occ
02582 write(6,'(300f15.8)')
02583 + (DBLE(H_MAT_R(ib,jb,ia1,ia2,ia3)*au),jb=1,n_occ)
02584 enddo
02585 write(6,*)
02586 enddo
02587 enddo
02588 enddo
02589
02590
02591 write(6,'(a40)')'WEIGHT CALCULATED'
02592 WEIGHT_R=1.0d0; SUM_REAL=0.0d0
02593 do ia1=-Na1,Na1
02594 do ia2=-Na2,Na2
02595 do ia3=-Na3,Na3
02596 if(abs(ia1)==Na1.and.mod(nkb1,2)==0.and.Na1/=0)
02597 + WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0
02598 if(abs(ia2)==Na2.and.mod(nkb2,2)==0.and.Na2/=0)
02599 + WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0
02600 if(abs(ia3)==Na3.and.mod(nkb3,2)==0.and.Na3/=0)
02601 + WEIGHT_R(ia1,ia2,ia3)=WEIGHT_R(ia1,ia2,ia3)*0.5d0
02602 SUM_REAL=SUM_REAL+WEIGHT_R(ia1,ia2,ia3)
02603 enddo
02604 enddo
02605 enddo
02606 write(6,'(a40,f15.8,i8)')'SUM_WEIGHT,NTK',SUM_REAL,NTK
02607 if(abs(SUM_REAL-dble(NTK))>1.0d-6)then
02608 write(6,'(a40)')'SUM_WEIGHT/=NTK'
02609 stop
02610 endif
02611
02612
02613 write(6,*)
02614 write(6,'(a40)')'+++ write transfers in dir-model +++'
02615 write(6,*)
02616 call system('rm -rf dir-model')
02617 call system('mkdir dir-model')
02618
02619 call wrt_model_hr(Na1,Na2,Na3,n_occ,H_MAT_R(1,1,-Na1,-Na2,-Na3),
02620 + WEIGHT_R(-Na1,-Na2,-Na3))
02621 call wrt_model_wcenter(n_occ,a1(1),a2(1),a3(1),
02622 + wannier_center(1,1))
02623 call wrt_model_SK_BAND_DISP(Ndiv,N_sym_points,NSK_BAND_DISP,
02624 + SK_BAND_DISP(1,1))
02625 call wrt_model_SK0(NTK,SK0(1,1))
02626 call wrt_model_ef(FermiEnergy)
02627
02628
02629
02630 if(elec_num==0.0d0)then
02631 call calculate_dmx(NTB,NTK,NWF,nkb1,nkb2,nkb3,Na1,Na2,Na3,
02632 + FermiEnergy,a1(1),a2(1),a3(1),b1(1),b2(1),b3(1),SK0(1,1),
02633 + E_EIG(1,1),UNT(1,1,1))
02634 endif
02635 deallocate(E_EIG)
02636
02637
02638 do ik=1,NSK_BAND_DISP
02639 H_TMP_IN(:,:)=0.0D0
02640 do ib=1,n_occ
02641 do jb=1,n_occ
02642 SUM_CMPX=0.0D0
02643 do ia1=-Na1,Na1
02644 do ia2=-Na2,Na2
02645 do ia3=-Na3,Na3
02646
02647
02648
02649 call search_Rmin(ia1,ia2,ia3,nkb1,nkb2,nkb3,
02650 + a1(1),a2(1),a3(1),ia1min,ia2min,ia3min)
02651 PHASE=tpi*(SK_BAND_DISP(1,ik)*DBLE(ia1min)
02652 + +SK_BAND_DISP(2,ik)*DBLE(ia2min)
02653 + +SK_BAND_DISP(3,ik)*DBLE(ia3min))
02654
02655 PHASE_FACTOR=EXP(ci*PHASE)*WEIGHT_R(ia1,ia2,ia3)
02656 SUM_CMPX=SUM_CMPX+H_MAT_R(ib,jb,ia1,ia2,ia3)*PHASE_FACTOR
02657 enddo
02658 enddo
02659 enddo
02660 H_TMP_IN(ib,jb)=SUM_CMPX
02661 enddo
02662 enddo
02663
02664 E_TMP(:)=0.0D0
02665 call diagV(n_occ,H_TMP_IN(1,1),E_TMP(1))
02666
02667
02668
02669 H_TMP_OUT(:,:,ik)=H_TMP_IN(:,:)
02670
02671 E_BAND_DISP(:,ik)=E_TMP(:)
02672 ENDDO
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682 allocate(DIST_K(NSK_BAND_DISP));DIST_K(:)=0.0d0
02683 allocate(dist(0:N_sym_points-1));dist(:)=0.0d0
02684
02685 dist(0)=0.0d0
02686 do ix=1,N_sym_points-1
02687 ks=(ix-1)*Ndiv+1
02688 ke=(ix)*Ndiv+1
02689 do ik=ks,ke
02690 DIST_B(:)
02691 + =(SK_BAND_DISP(1,ik)-SK_BAND_DISP(1,ks))*b1(:)
02692 + +(SK_BAND_DISP(2,ik)-SK_BAND_DISP(2,ks))*b2(:)
02693 + +(SK_BAND_DISP(3,ik)-SK_BAND_DISP(3,ks))*b3(:)
02694 DIST_KSPACE
02695 + =DSQRT(DIST_B(1)**2+DIST_B(2)**2+DIST_B(3)**2)
02696 DIST_K(ik)=DIST_KSPACE+dist(ix-1)
02697 enddo
02698 dist(ix)=dist(ix-1)+DIST_KSPACE
02699 enddo
02700
02701
02702
02703 OPEN(114,FILE='./dir-wan/dat.iband')
02704 write(114,'(a)')'#Wannier interpolaed band'
02705 write(114,'(a)')'#1:k, 2:Energy [eV]'
02706 REVERSE=.TRUE.
02707 do ib=1,n_occ
02708 if(REVERSE) then
02709 do ik=1,NSK_BAND_DISP
02710 write(114,*) DIST_K(ik)/DIST_K(NSK_BAND_DISP),
02711 + E_BAND_DISP(ib,ik)*au
02712 enddo
02713 REVERSE=.FALSE.
02714 else
02715 do ik=NSK_BAND_DISP,1,-1
02716 write(114,*) DIST_K(ik)/DIST_K(NSK_BAND_DISP),
02717 + E_BAND_DISP(ib,ik)*au
02718 enddo
02719 REVERSE=.TRUE.
02720 endif
02721 enddo
02722
02723
02724
02725
02726
02727
02728
02729
02730
02731 do ik=1,NSK_BAND_DISP
02732 H_TMP_IN=0.0d0
02733 do ib=1,n_occ
02734 do jb=1,n_occ
02735 H_TMP_IN(jb,ib)=H_TMP_OUT(ib,jb,ik)
02736 enddo
02737 enddo
02738 H_TMP_OUT(:,:,ik)=H_TMP_IN(:,:)
02739 enddo
02740
02741 call calc_fat_band(NWF,NSK_BAND_DISP,DIST_K(1),E_BAND_DISP(1,1),
02742 + H_TMP_OUT(1,1,1))
02743
02744
02745 IF(ALL(dense(1:3)/=0)) THEN
02746 CALL wrt_frmsf_wan(Na1,Na2,Na3,nkb1,nkb2,nkb3,
02747 + a1,a2,a3,b1,b2,b3,FermiEnergy,
02748 + WEIGHT_R,H_MAT_R)
02749 ENDIF
02750
02751
02752
02753
02754
02755 write(6,*)
02756 write(6,'(a40)')'==============================='
02757 write(6,'(a40)')' CALCULATE WANNIER DOS AND DMX '
02758 write(6,'(a40)')'==============================='
02759 write(6,*)
02760 allocate(E_EIG(NWF,NTK)); E_EIG(:,:)=0.0d0
02761 allocate(V_EIG(NWF,NWF,NTK)); V_EIG(:,:,:)=0.0d0
02762 call calculate_eigenstate(NWF,Na1,Na2,Na3,
02763 + H_MAT_R(1,1,-Na1,-Na2,-Na3),NTK,nkb1,nkb2,nkb3,
02764 + a1(1),a2(1),a3(1),1,NTK,SK0(1,1),E_EIG(1,1),V_EIG(1,1,1))
02765 filename='./dir-wan/dat.dos.wan'
02766
02767
02768
02769 if(elec_num==0.0d0)then
02770 call calculate_dos(NWF,NTK,nkb1,nkb2,nkb3,elec_num,FermiEnergy,
02771 + filename,b1(1),b2(1),b3(1),SK0(1,1),E_EIG(1,1),V_EIG(1,1,1))
02772 endif
02773
02774
02775
02776 if(elec_num/=0.0d0)then
02777 call calculate_dos(NWF,NTK,nkb1,nkb2,nkb3,elec_num,FermiEnergy,
02778 + filename,b1(1),b2(1),b3(1),SK0(1,1),E_EIG(1,1),V_EIG(1,1,1))
02779 call calculate_dmx(NWF,NTK,NWF,nkb1,nkb2,nkb3,Na1,Na2,Na3,
02780 + FermiEnergy,a1(1),a2(1),a3(1),b1(1),b2(1),b3(1),SK0(1,1),
02781 + E_EIG(1,1),V_EIG(1,1,1))
02782 endif
02783 deallocate(E_EIG,V_EIG)
02784
02785
02786
02787
02788
02789 if(CALC_REAL_SPACE_WANNIER)then
02790
02791 write(6,*)
02792 write(6,'(a40)')'============================='
02793 write(6,'(a40)')' CALCULATE REALSPACE WANNIER '
02794 write(6,'(a40)')'============================='
02795 write(6,*)
02796 nfft1=nwx2+1;nfft2=nwy2+1;nfft3=nwz2+1
02797 Nl123=nfft1*nfft2*nfft3
02798 write(6,'(a40,i5,i5,i5)')'nwx2,nwy2,nwz2',nwx2,nwy2,nwz2
02799 write(6,'(a40,i5,i5,i5)')'nfft1,nfft2,nfft2',nfft1,nfft2,nfft3
02800 write(6,'(a40,i10)')'Nl123',Nl123
02801 call fft3_init(nwx2,nwy2,nwz2,nfft1,nfft2,nfft3,fs)
02802
02803 mem_size=dble((xmax-xmin)*nwx2+(ymax-ymin)*nwy2+(zmax-zmin)*nwz2)
02804 mem_size=mem_size*16.0d0/1024.0d0/1024.0d0/1024.0d0
02805 write(6,'(a40,f20.15)')'mem size realspace wannier (GB)',mem_size
02806 if(mem_size>1.0d0)then
02807 write(6,'(a40,f20.15)')'mem size > 1 GB; skipp'
02808 go to 9999
02809 endif
02810
02811 allocate(WANNIER_REALSPACE((xmin)*nwx2:(xmax)*nwx2-1,
02812 + (ymin)*nwy2:(ymax)*nwy2-1,
02813 + (zmin)*nwz2:(zmax)*nwz2-1))
02814
02815
02816 ierr=CHDIR("./dir-wan")
02817 call system('pwd')
02818 REWIND(113)
02819 read(113) NWF
02820 allocate(C0WN(NTG,NWF,NTK));C0WN(:,:,:)=0.0D0
02821 do ik=1,NTK
02822 read(113)((C0WN(ig,jb,ik),ig=1,NG0(ik)),jb=1,NWF)
02823
02824
02825
02826
02827
02828
02829 enddo
02830
02831
02832 do ix=1,N_write_wannier
02833 iw=wrt_list(ix)
02834 WANNIER_REALSPACE=0.0d0
02835
02836 allocate(C0_I(NTG));C0_I=0.0D0
02837 allocate(fftwk(Nl123*2),stat=err)
02838 allocate(wfunc(Nl123*2),stat=err)
02839 allocate(pw((xmin)*nwx2:(xmax)*nwx2-1,
02840 + (ymin)*nwy2:(ymax)*nwy2-1,
02841 + (zmin)*nwz2:(zmax)*nwz2-1));pw=0.0d0
02842
02843 do ik=1,NTK
02844 C0_I(:)=C0WN(:,iw,ik)
02845 call make_pw(C0_I(1),wfunc(1),fftwk(1),NG0(ik),KG0(1,1,ik),NTG,
02846 + nwx2,nwy2,nwz2,nfft1,nfft2,Nl123,fs,
02847 + xmin,xmax,ymin,ymax,zmin,zmax,SK0(1,ik),
02848 + pw(xmin*nwx2,ymin*nwy2,zmin*nwz2))
02849 enddo
02850
02851
02852 WANNIER_REALSPACE=WANNIER_REALSPACE+pw
02853
02854 deallocate(C0_I,fftwk,wfunc,pw)
02855
02856
02857
02858 WANNIER_REALSPACE=WANNIER_REALSPACE/DBLE(NTK)/DSQRT(VOLUME)
02859
02860
02861
02862 write(filename,"('dat.wan-realspace-',i3.3,'.grd')") iw
02863 OPEN(116,FILE=filename)
02864 REWIND(116)
02865 na_grid=nwx2
02866 nb_grid=nwy2
02867 nc_grid=nwz2
02868 aa1=a1*dble(xmax-xmin)
02869 aa2=a2*dble(ymax-ymin)
02870 aa3=a3*dble(zmax-zmin)
02871 call est_latparam(aa1(1),aa2(1),aa3(1),a,b,c,alp,bet,gmm)
02872 write(116,*)'Electron density'
02873 write(116,'(6f15.10)')
02874 + a*bohr,b*bohr,c*bohr,alp,bet,gmm
02875 write(116,'(3I5)')
02876 + (xmax-xmin)*na_grid,(ymax-ymin)*nb_grid,(zmax-zmin)*nc_grid
02877 write(116,"(6g25.16)") (((dble(WANNIER_REALSPACE(i1,i2,i3)),
02878 + i3=(zmin)*nc_grid,(zmax)*nc_grid-1),
02879 + i2=(ymin)*nb_grid,(ymax)*nb_grid-1),
02880 + i1=(xmin)*na_grid,(xmax)*na_grid-1)
02881 CLOSE(116)
02882
02883
02884
02885
02886
02887
02888
02889
02890
02891
02892
02893 enddo
02894
02895
02896 call gencif(a1(1),a2(1),a3(1),xmin,xmax,ymin,ymax,zmin,zmax,
02897 + nsi,chemical_species(1),asi(1,1))
02898
02899 ierr=CHDIR("..")
02900 call system('pwd')
02901 9999 endif
02902
02903 STOP
02904 END
