00001 subroutine kcheck(ktmp,RWtmp)
00002 implicit none
00003 real(8),intent(inout)::ktmp(3)
00004 integer,intent(out)::RWtmp(3)
00005 real(8),parameter::dlt_BZ=1.0d-6
00006
00007 if(ktmp(1)>1.50d0+dlt_BZ)then
00008 ktmp(1)=ktmp(1)-2.0d0
00009 RWtmp(1)=-2
00010 endif
00011 if(ktmp(1)>0.50d0+dlt_BZ)then
00012 ktmp(1)=ktmp(1)-1.0d0
00013 RWtmp(1)=-1
00014 endif
00015 if(ktmp(1)<=-1.50d0+dlt_BZ)then
00016 ktmp(1)=ktmp(1)+2.0d0
00017 RWtmp(1)=2
00018 endif
00019 if(ktmp(1)<=-0.50d0+dlt_BZ)then
00020 ktmp(1)=ktmp(1)+1.0d0
00021 RWtmp(1)=1
00022 endif
00023
00024 if(ktmp(2)>1.50d0+dlt_BZ)then
00025 ktmp(2)=ktmp(2)-2.0d0
00026 RWtmp(2)=-2
00027 endif
00028 if(ktmp(2)>0.50d0+dlt_BZ)then
00029 ktmp(2)=ktmp(2)-1.0d0
00030 RWtmp(2)=-1
00031 endif
00032 if(ktmp(2)<=-1.50d0+dlt_BZ)then
00033 ktmp(2)=ktmp(2)+2.0d0
00034 RWtmp(2)=2
00035 endif
00036 if(ktmp(2)<=-0.50d0+dlt_BZ)then
00037 ktmp(2)=ktmp(2)+1.0d0
00038 RWtmp(2)=1
00039 endif
00040
00041 if(ktmp(3)>1.50d0+dlt_BZ)then
00042 ktmp(3)=ktmp(3)-2.0d0
00043 RWtmp(3)=-2
00044 endif
00045 if(ktmp(3)>0.50d0+dlt_BZ)then
00046 ktmp(3)=ktmp(3)-1.0d0
00047 RWtmp(3)=-1
00048 endif
00049 if(ktmp(3)<=-1.50d0+dlt_BZ)then
00050 ktmp(3)=ktmp(3)+2.0d0
00051 RWtmp(3)=2
00052 endif
00053 if(ktmp(3)<=-0.50d0+dlt_BZ)then
00054 ktmp(3)=ktmp(3)+1.0d0
00055 RWtmp(3)=1
00056 endif
00057 return
00058 end
00059
00060 subroutine kcheck_trs(ktmp,RWtmp)
00061 implicit none
00062 real(8),intent(inout)::ktmp(3)
00063 integer,intent(out)::RWtmp(3)
00064 real(8),parameter::dlt_BZ=-1.0d-6
00065
00066 if(ktmp(1)>=1.50d0+dlt_BZ)then
00067 ktmp(1)=ktmp(1)-2.0d0
00068 RWtmp(1)=-2
00069 endif
00070 if(ktmp(1)>=0.50d0+dlt_BZ)then
00071 ktmp(1)=ktmp(1)-1.0d0
00072 RWtmp(1)=-1
00073 endif
00074 if(ktmp(1)<-1.50d0+dlt_BZ)then
00075 ktmp(1)=ktmp(1)+2.0d0
00076 RWtmp(1)=2
00077 endif
00078 if(ktmp(1)<-0.50d0+dlt_BZ)then
00079 ktmp(1)=ktmp(1)+1.0d0
00080 RWtmp(1)=1
00081 endif
00082
00083 if(ktmp(2)>=1.50d0+dlt_BZ)then
00084 ktmp(2)=ktmp(2)-2.0d0
00085 RWtmp(2)=-2
00086 endif
00087 if(ktmp(2)>=0.50d0+dlt_BZ)then
00088 ktmp(2)=ktmp(2)-1.0d0
00089 RWtmp(2)=-1
00090 endif
00091 if(ktmp(2)<-1.50d0+dlt_BZ)then
00092 ktmp(2)=ktmp(2)+2.0d0
00093 RWtmp(2)=2
00094 endif
00095 if(ktmp(2)<-0.50d0+dlt_BZ)then
00096 ktmp(2)=ktmp(2)+1.0d0
00097 RWtmp(2)=1
00098 endif
00099
00100 if(ktmp(3)>=1.50d0+dlt_BZ)then
00101 ktmp(3)=ktmp(3)-2.0d0
00102 RWtmp(3)=-2
00103 endif
00104 if(ktmp(3)>=0.50d0+dlt_BZ)then
00105 ktmp(3)=ktmp(3)-1.0d0
00106 RWtmp(3)=-1
00107 endif
00108 if(ktmp(3)<-1.50d0+dlt_BZ)then
00109 ktmp(3)=ktmp(3)+2.0d0
00110 RWtmp(3)=2
00111 endif
00112 if(ktmp(3)<-0.50d0+dlt_BZ)then
00113 ktmp(3)=ktmp(3)+1.0d0
00114 RWtmp(3)=1
00115 endif
00116 return
00117 end
00118
00119 subroutine make_KG0(NTG,b1,b2,b3,Gcut,q1,q2,q3,KG0,NG)
00120 implicit none
00121 integer,intent(in)::NTG
00122 real(8),intent(in)::b1(3),b2(3),b3(3)
00123 real(8),intent(in)::Gcut,q1,q2,q3
00124 integer,intent(out)::KG0(3,NTG),NG
00125 integer::igL,igL1,igL2,igL3
00126 real(8)::qgL(3),qgL2
00127 integer,parameter::NGL1=150
00128 integer,parameter::NGL2=150
00129 integer,parameter::NGL3=150
00130
00131 igL=0
00132 do igL1=-NGL1,NGL1
00133 do igL2=-NGL2,NGL2
00134 do igL3=-NGL3,NGL3
00135 qgL(:)=(q1+dble(igL1))*b1(:)+(q2+dble(igL2))*b2(:)+(q3+dble(igL3))*b3(:)
00136 qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00137 if(qgL2<=Gcut)then
00138 igL=igL+1
00139 KG0(1,igL)=igL1
00140 KG0(2,igL)=igL2
00141 KG0(3,igL)=igL3
00142 endif
00143 enddo
00144 enddo
00145 enddo
00146 NG=igL
00147
00148 RETURN
00149 END
00150
00151 subroutine est_nkbi(N,SK,nkb1,nkb2,nkb3)
00152 implicit none
00153 integer::N,nkb1,nkb2,nkb3,NTK
00154 real(8)::SK(3,N)
00155 integer::i
00156 real(8)::x
00157
00158 x=1.0d0
00159 do i=1,N
00160 if(abs(SK(1,i))<1.0d-7)cycle
00161 if(abs(SK(1,i))<x)then
00162 x=abs(SK(1,i))
00163 endif
00164 enddo
00165 nkb1=nint(1.0d0/x)
00166
00167 x=1.0d0
00168 do i=1,N
00169 if(abs(SK(2,i))<1.0d-7)cycle
00170 if(abs(SK(2,i))<x)then
00171 x=abs(SK(2,i))
00172 endif
00173 enddo
00174 nkb2=nint(1.0d0/x)
00175
00176 x=1.0d0
00177 do i=1,N
00178 if(abs(SK(3,i))<1.0d-7)cycle
00179 if(abs(SK(3,i))<x)then
00180 x=abs(SK(3,i))
00181 endif
00182 enddo
00183 nkb3=nint(1.0d0/x)
00184
00185 NTK=nkb1*nkb2*nkb3
00186
00187
00188
00189 return
00190 end
00191
00192 subroutine est_NTK(Nk_irr,Nsymq,SKI,rg,NTK)
00193
00194
00195
00196 implicit none
00197 integer::Nk_irr,Nsymq,N
00198 real(8)::SKI(3,Nk_irr)
00199 integer::rg(3,3,Nsymq)
00200 real(8),allocatable::SK0(:,:)
00201 real(8)::ktmp(3)
00202 integer::RWtmp(3)
00203 integer::jk,ik,iop,iik
00204 integer::NTK
00205
00206 N=Nk_irr*Nsymq*2
00207
00208
00209
00210
00211
00212 allocate(SK0(3,N));SK0(:,:)=0.0d0
00213 jk=0
00214 do ik=1,Nk_irr
00215 do iop=1,Nsymq
00216 ktmp(:)=0.0d0; RWtmp(:)=0
00217 ktmp(1)=dble(rg(1,1,iop))*SKI(1,ik)+dble(rg(1,2,iop))*SKI(2,ik)+dble(rg(1,3,iop))*SKI(3,ik)
00218 ktmp(2)=dble(rg(2,1,iop))*SKI(1,ik)+dble(rg(2,2,iop))*SKI(2,ik)+dble(rg(2,3,iop))*SKI(3,ik)
00219 ktmp(3)=dble(rg(3,1,iop))*SKI(1,ik)+dble(rg(3,2,iop))*SKI(2,ik)+dble(rg(3,3,iop))*SKI(3,ik)
00220 call kcheck(ktmp(1),RWtmp(1))
00221 do iik=1,jk
00222 if(abs(SK0(1,iik)-ktmp(1))<1.0d-4.and.abs(SK0(2,iik)-ktmp(2))<1.0d-4.and.abs(SK0(3,iik)-ktmp(3))<1.0d-4) goto 1000
00223 enddo
00224 jk=jk+1
00225 SK0(:,jk)=ktmp(:)
00226 1000 ktmp(:)=0.0d0; RWtmp(:)=0
00227 ktmp(1)=dble(rg(1,1,iop))*SKI(1,ik)+dble(rg(1,2,iop))*SKI(2,ik)+dble(rg(1,3,iop))*SKI(3,ik)
00228 ktmp(2)=dble(rg(2,1,iop))*SKI(1,ik)+dble(rg(2,2,iop))*SKI(2,ik)+dble(rg(2,3,iop))*SKI(3,ik)
00229 ktmp(3)=dble(rg(3,1,iop))*SKI(1,ik)+dble(rg(3,2,iop))*SKI(2,ik)+dble(rg(3,3,iop))*SKI(3,ik)
00230 call kcheck_trs(ktmp(1),RWtmp(1))
00231 do iik=1,jk
00232 if(abs(SK0(1,iik)-(-ktmp(1)))<1.0d-4.and.abs(SK0(2,iik)-(-ktmp(2)))<1.0d-4.and.abs(SK0(3,iik)-(-ktmp(3)))<1.0d-4) goto 2000
00233 enddo
00234 jk=jk+1
00235 SK0(:,jk)=-ktmp(:)
00236 2000 enddo
00237 enddo
00238
00239 NTK=jk
00240 if(NTK>N)then
00241 write(6,*)'Estimated NTK is too large; stop'
00242 write(6,*)'NTK, N=',NTK, N
00243 stop
00244 endif
00245
00246
00247
00248 return
00249 end
00250
00251 subroutine est_latparam(a1,a2,a3,a,b,c,alp,bet,gmm)
00252 implicit none
00253 real(8)::a1(3),a2(3),a3(3)
00254 real(8)::a,b,c,alp,bet,gmm
00255 real(8),parameter::pi=DACOS(-1.0d0)
00256
00257 a=0.0d0
00258 b=0.0d0
00259 c=0.0d0
00260 a=dsqrt(a1(1)**2+a1(2)**2+a1(3)**2)
00261 b=dsqrt(a2(1)**2+a2(2)**2+a2(3)**2)
00262 c=dsqrt(a3(1)**2+a3(2)**2+a3(3)**2)
00263 alp=(a2(1)*a3(1)+a2(2)*a3(2)+a2(3)*a3(3))/b/c
00264 bet=(a3(1)*a1(1)+a3(2)*a1(2)+a3(3)*a1(3))/c/a
00265 gmm=(a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3))/a/b
00266 alp=dacos(alp)*180.0d0/pi
00267 bet=dacos(bet)*180.0d0/pi
00268 gmm=dacos(gmm)*180.0d0/pi
00269
00270
00271
00272 return
00273 end
00274
00275 integer function algn235(inr)
00276 implicit none
00277 integer,intent(in)::inr
00278 integer:: nr,m2,m3,m5,info
00279 nr=inr
00280 call fctck(nr,m2,m3,m5,info)
00281 do while (info .eq. 1)
00282 nr = nr + 1
00283 call fctck(nr,m2,m3,m5,info)
00284 end do
00285 algn235 = nr
00286 return
00287 end function algn235
00288
00289 subroutine fctck(n,m2,m3,m5,info)
00290 implicit none
00291 integer,intent(in):: n
00292 integer,intent(out):: m2, m3, m5, info
00293 integer:: i
00294 i=n
00295 m2 = 0
00296 m3 = 0
00297 m5 = 0
00298 info = 0
00299 do while (i .ne. 1)
00300 if (mod(i,2) .eq. 0) then
00301 m2 = m2 + 1
00302 i = i / 2
00303 else if (mod(i,3) .eq. 0) then
00304 m3 = m3 + 1
00305 i = i / 3
00306 else if (mod(i,5) .eq. 0) then
00307 m5 = m5 + 1
00308 i = i / 5
00309 else
00310 info = 1
00311 exit
00312 end if
00313 end do
00314 return
00315 end subroutine fctck
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372 subroutine make_C0(NTG,NTB,itrs,NG,KGtmp,RWtmp,rginvtmp,pgtmp,nnp,L1,L2,L3,packtmp,C0irr,C0_K)
00373 implicit none
00374 integer::NTG,NTB,itrs,NG,L1,L2,L3,nnp
00375 integer::KGtmp(3,NTG),RWtmp(3),pgtmp(3)
00376 integer::packtmp(-L1:L1,-L2:L2,-L3:L3)
00377 integer::ig,jg,i1,i2,i3,j1,j2,j3,k1,k2,k3,ib
00378 real(8)::rginvtmp(3,3),phase
00379 complex(8)::C0irr(NTG,NTB)
00380 complex(8)::pf
00381 complex(8),intent(out)::C0_K(NTG,NTB)
00382 real(8),parameter::tpi=2.0d0*dacos(-1.0d0)
00383 complex(8),parameter::ci=(0.0D0,1.0D0)
00384
00385 C0_K(:,:)=0.0d0
00386 write(6,*)'nnp=',nnp
00387 select case(itrs)
00388 case(1)
00389 do ig=1,NG
00390 i1=KGtmp(1,ig); j1=KGtmp(2,ig); k1=KGtmp(3,ig)
00391 i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3)
00392 i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2+int(rginvtmp(1,3))*k2
00393 j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2+int(rginvtmp(2,3))*k2
00394 k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2+int(rginvtmp(3,3))*k2
00395 jg=packtmp(i3,j3,k3)
00396 phase=tpi*(dble(i1)*dble(pgtmp(1))+dble(j1)*dble(pgtmp(2))+dble(k1)*dble(pgtmp(3)))
00397 pf=exp(-ci*phase/dble(nnp))
00398 write(6,*)'pf=',pf
00399
00400 do ib=1,NTB
00401 write(6,*) ig,jg
00402 C0_K(ig,ib)=C0irr(jg,ib)*pf
00403 enddo
00404 enddo
00405 case(-1)
00406 do ig=1,NG
00407 i1=KGtmp(1,ig); j1=KGtmp(2,ig); k1=KGtmp(3,ig)
00408 i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3)
00409 i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2+int(rginvtmp(1,3))*k2
00410 j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2+int(rginvtmp(2,3))*k2
00411 k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2+int(rginvtmp(3,3))*k2
00412 jg=packtmp(i3,j3,k3)
00413 phase=tpi*(dble(i1)*dble(pgtmp(1))+dble(j1)*dble(pgtmp(2))+dble(k1)*dble(pgtmp(3)))
00414 pf=exp(-ci*phase/dble(nnp))
00415
00416 do ib=1,NTB
00417 write(6,*) ig,jg
00418 C0_K(ig,ib)=C0irr(jg,ib)*pf
00419 enddo
00420 enddo
00421 C0_K(:,:)=conjg(C0_K(:,:))
00422 end select
00423
00424 return
00425 end
00426
00427 subroutine make_C0_for_given_band(NTG,itrs,NG,KGtmp,RWtmp,rginvtmp,pgtmp,nnp,L1,L2,L3,packtmp,OCCtmp,C0_K)
00428 implicit none
00429 integer::NTG,L1,L2,L3,nnp
00430 integer::itrs
00431 integer::NG
00432 integer::KGtmp(3,NTG)
00433 integer::RWtmp(3)
00434 real(8)::rginvtmp(3,3)
00435 integer::pgtmp(3)
00436 integer::packtmp(-L1:L1,-L2:L2,-L3:L3)
00437
00438
00439
00440
00441 complex(4)::OCCtmp(NTG)
00442
00443 integer::ig,jg,i1,i2,i3,j1,j2,j3,k1,k2,k3
00444 real(8),parameter::pi=dacos(-1.0d0)
00445 real(8),parameter::tpi=2.0d0*pi
00446 complex(8),parameter::ci=(0.0D0,1.0D0)
00447 real(8)::phase
00448 complex(8)::pf
00449 complex(8)::C0_K(NTG)
00450
00451 C0_K(:)=0.0d0
00452 select case(itrs)
00453 case(1)
00454 do ig=1,NG
00455 i1=KGtmp(1,ig); j1=KGtmp(2,ig); k1=KGtmp(3,ig)
00456 i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3)
00457 i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2+int(rginvtmp(1,3))*k2
00458 j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2+int(rginvtmp(2,3))*k2
00459 k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2+int(rginvtmp(3,3))*k2
00460 jg=packtmp(i3,j3,k3)
00461 phase=tpi*(dble(i1)*dble(pgtmp(1))+dble(j1)*dble(pgtmp(2))+dble(k1)*dble(pgtmp(3)))
00462 pf=exp(-ci*phase/dble(nnp))
00463 C0_K(ig)=OCCtmp(jg)*pf
00464 enddo
00465 case(-1)
00466 do ig=1,NG
00467 i1=-KGtmp(1,ig); j1=-KGtmp(2,ig); k1=-KGtmp(3,ig)
00468 i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3)
00469 i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2+int(rginvtmp(1,3))*k2
00470 j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2+int(rginvtmp(2,3))*k2
00471 k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2+int(rginvtmp(3,3))*k2
00472 jg=packtmp(i3,j3,k3)
00473 phase=tpi*(dble(i1)*dble(pgtmp(1))+dble(j1)*dble(pgtmp(2))+dble(k1)*dble(pgtmp(3)))
00474 pf=exp(-ci*phase/dble(nnp))
00475 C0_K(ig)=OCCtmp(jg)*pf
00476 enddo
00477 C0_K(:)=conjg(C0_K(:))
00478 end select
00479 return
00480 end