00001 SUBROUTINE CALC_GRADIENT(NTK,NTB,n_occ,NB,WF_CHARGE,WEIGHT_b,
00002 + OVERLAP,tci,GRADIENT,TR_GRAD)
00003 implicit none
00004 integer::NTK,NTB,n_occ,NB
00005 integer::ik,ibvec,i_band,j_band
00006 real(8)::WF_CHARGE(n_occ,NTK,NB),WEIGHT_b(NB)
00007 real(8)::TR_GRAD(NTK)
00008 real(8)::SUM_REAL
00009 complex(8)::OVERLAP(NTB,NTB,NTK,NB)
00010 complex(8)::GRADIENT(n_occ,n_occ,NTK)
00011 complex(8)::WF_R(n_occ,n_occ,NB)
00012 complex(8)::WF_T(n_occ,n_occ,NB)
00013 complex(8)::SUM_CMPX,tci
00014
00015 GRADIENT(:,:,:)=0.0D0
00016 Do ik=1,NTK
00017 Do i_band=1,n_occ
00018 Do j_band=1,n_occ
00019 Do ibvec=1,NB
00020 WF_R(i_band,j_band,ibvec)
00021 + =OVERLAP(i_band,j_band,ik,ibvec)
00022 + *CONJG(OVERLAP(j_band,j_band,ik,ibvec))
00023 WF_T(i_band,j_band,ibvec)
00024 + =OVERLAP(i_band,j_band,ik,ibvec)
00025 + /OVERLAP(j_band,j_band,ik,ibvec)
00026 + *WF_CHARGE(j_band,ik,ibvec)
00027 ENDDO
00028 ENDDO
00029 ENDDO
00030 Do i_band=1,n_occ
00031 Do j_band=1,n_occ
00032 SUM_CMPX=0.0D0
00033 Do ibvec=1,NB
00034 SUM_CMPX
00035 + =SUM_CMPX
00036 + +4.0D0
00037 + *WEIGHT_b(ibvec)
00038 + *(((WF_R(i_band,j_band,ibvec)
00039 + -CONJG(WF_R(j_band,i_band,ibvec)))/2.0D0)
00040 + -((WF_T(i_band,j_band,ibvec)
00041 + +CONJG(WF_T(j_band,i_band,ibvec)))/tci))
00042 ENDDO
00043 GRADIENT(i_band,j_band,ik)=SUM_CMPX/CMPLX(NTK)
00044 ENDDO
00045 ENDDO
00046 ENDDO
00047
00048
00049
00050
00051
00052
00053
00054 TR_GRAD(:)=0.0D0
00055 SUM_REAL=0.0D0
00056 Do ik=1,NTK
00057 SUM_CMPX=0.0D0
00058 Do i_band=1,n_occ
00059 Do j_band=1,n_occ
00060 SUM_CMPX
00061 + =SUM_CMPX
00062 + +GRADIENT(i_band,j_band,ik)
00063 + *GRADIENT(j_band,i_band,ik)
00064 ENDDO
00065 ENDDO
00066 TR_GRAD(ik)=DBLE(SUM_CMPX)
00067 SUM_REAL=SUM_REAL+TR_GRAD(ik)
00068 ENDDO
00069 write(6,*) 'Tr(G^2)=',SUM_REAL
00070 RETURN
00071 END
00072
00073 SUBROUTINE CALC_OMEGA(NTK,NTB,n_occ,NB,VEC_b,ILS,
00074 + WEIGHT_b,OVERLAP,WF_CHARGE,OMEGA_I,OMEGA_OD,OMEGA_D)
00075 implicit none
00076 integer::NTK,NTB,n_occ,NB,ILS
00077 integer::ik,ibvec,i_band,j_band
00078 real(8)::VEC_b(3,NB),WEIGHT_b(NB)
00079 real(8)::WF_CHARGE(n_occ,NTK,NB)
00080 real(8)::OMEGA_I,OMEGA_OD,OMEGA_D
00081 real(8)::WF_POSITION(3,n_occ)
00082 real(8)::VEC_SUM_REAL(3),SUM_REAL,SUM_REAL1,SUM_REAL2
00083 complex(8)::LOG_OVERLAP
00084 complex(8)::OVERLAP(NTB,NTB,NTK,NB)
00085
00086
00087
00088 WF_POSITION(:,:) = 0.0D0
00089 Do i_band = 1 , n_occ
00090 VEC_SUM_REAL(:) = 0.0D0
00091 Do ik = 1 , NTK
00092 Do ibvec = 1 , NB
00093 LOG_OVERLAP=LOG(OVERLAP(i_band,i_band,ik,ibvec))
00094 VEC_SUM_REAL(:)
00095 + = VEC_SUM_REAL(:)
00096 + + WEIGHT_b(ibvec)
00097 + * VEC_b (:,ibvec)
00098 + * AIMAG( LOG_OVERLAP )
00099 ENDDO
00100 ENDDO
00101 WF_POSITION(:,i_band)=-VEC_SUM_REAL(:)/DBLE(NTK)
00102 ENDDO
00103
00104
00105
00106 WF_CHARGE(:,:,:) = 0.0D0
00107 Do ik = 1 , NTK
00108 Do ibvec = 1 , NB
00109 Do i_band = 1 , n_occ
00110 SUM_REAL = VEC_b(1,ibvec)*WF_POSITION(1,i_band)
00111 + + VEC_b(2,ibvec)*WF_POSITION(2,i_band)
00112 + + VEC_b(3,ibvec)*WF_POSITION(3,i_band)
00113 LOG_OVERLAP = LOG ( OVERLAP(i_band,i_band,ik,ibvec) )
00114 WF_CHARGE(i_band,ik,ibvec)
00115 + = AIMAG( LOG_OVERLAP )
00116 + + SUM_REAL
00117 ENDDO
00118 ENDDO
00119 ENDDO
00120
00121
00122
00123 SUM_REAL1 = 0.0D0
00124 Do ik = 1 , NTK
00125 Do ibvec = 1 , NB
00126 SUM_REAL2 = 0.0D0
00127 Do i_band = 1 , n_occ
00128 Do j_band = 1 , n_occ
00129 SUM_REAL2
00130 + = SUM_REAL2
00131 + + ( ABS( OVERLAP(i_band,j_band,ik,ibvec) ) )**2
00132 ENDDO
00133 ENDDO
00134 SUM_REAL1
00135 + = SUM_REAL1
00136 + + WEIGHT_b(ibvec)
00137 + * ( DBLE(n_occ) - SUM_REAL2 )
00138 ENDDO
00139 ENDDO
00140 OMEGA_I = SUM_REAL1 / DBLE(NTK)
00141
00142
00143
00144 SUM_REAL1 = 0.0D0
00145 Do ik = 1 , NTK
00146 Do ibvec = 1 , NB
00147 SUM_REAL2 = 0.0D0
00148 Do i_band = 1 , n_occ
00149 Do j_band = 1 , n_occ
00150 IF(i_band /= j_band) THEN
00151 SUM_REAL2
00152 + = SUM_REAL2
00153 + + ( ABS( OVERLAP(i_band,j_band,ik,ibvec) ) )**2
00154 ENDIF
00155 ENDDO
00156 ENDDO
00157 SUM_REAL1
00158 + = SUM_REAL1
00159 + + WEIGHT_b(ibvec)
00160 + * SUM_REAL2
00161 ENDDO
00162 ENDDO
00163 OMEGA_OD = SUM_REAL1 / DBLE(NTK)
00164
00165
00166
00167 SUM_REAL1 = 0.0D0
00168 Do ik = 1 , NTK
00169 Do ibvec = 1 , NB
00170 SUM_REAL2 = 0.0D0
00171 Do i_band = 1 , n_occ
00172 SUM_REAL2
00173 + = SUM_REAL2
00174 + + ( WF_CHARGE(i_band,ik,ibvec) )**2
00175 ENDDO
00176 SUM_REAL1
00177 + = SUM_REAL1
00178 + + WEIGHT_b(ibvec)
00179 + * SUM_REAL2
00180 ENDDO
00181 ENDDO
00182 OMEGA_D = SUM_REAL1 / DBLE(NTK)
00183
00184
00185
00186 IF(ILS == 1) THEN
00187 Do i_band = 1 , n_occ
00188 SUM_REAL1 = 0.0D0
00189 SUM_REAL2 = 0.0D0
00190 Do ik = 1 , NTK
00191 Do ibvec = 1 , NB
00192 SUM_REAL1
00193 + =SUM_REAL1
00194 + +WEIGHT_b(ibvec)
00195 + *(1.0D0-(ABS(OVERLAP(i_band,i_band,ik,ibvec)))**2)
00196 SUM_REAL2
00197 + =SUM_REAL2
00198 + +WEIGHT_b(ibvec)
00199 + *(WF_CHARGE(i_band,ik,ibvec))**2
00200 ENDDO
00201 ENDDO
00202 SUM_REAL1=SUM_REAL1/DBLE(NTK)
00203 SUM_REAL2=SUM_REAL2/DBLE(NTK)
00204 write(6,*) i_band,SUM_REAL1,SUM_REAL2
00205 ENDDO
00206 ENDIF
00207
00208 RETURN
00209 END
00210
00211 SUBROUTINE UNITARY_GEN(NTK,n_occ,ILS,
00212 + DIRECTION,DIRECTION2,STEP_LENGTH,UNITARY,U_ORG,U_OLD)
00213 implicit none
00214 integer::NTK,n_occ,ILS
00215 integer::ik,i_band,j_band,k_band
00216 real(8)::STEP_LENGTH
00217 real(8)::D_EIG(n_occ),D_INV(n_occ),D_COS(n_occ),D_SIN(n_occ)
00218 complex(8),INTENT(IN)::DIRECTION (n_occ,n_occ,NTK)
00219 complex(8),INTENT(IN)::DIRECTION2(n_occ,n_occ,NTK)
00220 complex(8)::UNITARY(n_occ,n_occ,NTK)
00221 complex(8)::U_ORG(n_occ,n_occ,NTK)
00222 complex(8)::U_OLD(n_occ,n_occ,NTK)
00223 complex(8)::D_MAT(n_occ,n_occ)
00224 complex(8)::R_MAT(n_occ,n_occ)
00225 complex(8)::RC_MAT(n_occ,n_occ)
00226 complex(8)::RS_MAT(n_occ,n_occ)
00227 complex(8)::G_MAT(n_occ,n_occ)
00228 complex(8)::SUM_CMPX
00229
00230 UNITARY(:,:,:)=0.0D0
00231 IF(ILS==1) THEN
00232 Do ik=1,NTK
00233
00234 D_MAT(:,:)=DIRECTION2(:,:,ik)
00235 D_EIG(:)=0.0D0
00236 call diagV(n_occ,D_MAT(1,1),D_EIG(1))
00237 R_MAT(:,:)=D_MAT(:,:)
00238
00239
00240 D_INV(:)=0.0D0
00241 D_COS(:)=0.0D0
00242 D_SIN(:)=0.0D0
00243 Do i_band=1,n_occ
00244 IF(D_EIG(i_band)==0.0D0) THEN
00245 D_EIG(i_band)=1.0D-20
00246 ENDIF
00247 D_EIG(i_band)=DSQRT(ABS(D_EIG(i_band)))
00248 D_INV(i_band)=1.0D0/D_EIG(i_band)
00249 D_COS(i_band)=DCOS(D_EIG(i_band)*STEP_LENGTH)
00250 D_SIN(i_band)=DSIN(D_EIG(i_band)*STEP_LENGTH)
00251 ENDDO
00252
00253
00254 RC_MAT(:,:)=0.0D0
00255 RS_MAT(:,:)=0.0D0
00256 Do i_band=1,n_occ
00257 Do j_band=1,n_occ
00258 RC_MAT(i_band,j_band)
00259 + =R_MAT(i_band,j_band)
00260 + *D_COS(j_band)
00261 RS_MAT(i_band,j_band)
00262 + =R_MAT(i_band,j_band)
00263 + *D_SIN(j_band)
00264 + *D_INV(j_band)
00265 ENDDO
00266 ENDDO
00267
00268
00269 G_MAT(:,:)=DIRECTION(:,:,ik)
00270 Do i_band=1,n_occ
00271 Do j_band=1,n_occ
00272 SUM_CMPX=0.0D0
00273 Do k_band=1,n_occ
00274 SUM_CMPX
00275 + =SUM_CMPX
00276 + +G_MAT(i_band,k_band)
00277 + *RS_MAT(k_band,j_band)
00278 ENDDO
00279 RC_MAT(i_band,j_band)
00280 + =RC_MAT(i_band,j_band)+SUM_CMPX
00281 ENDDO
00282 ENDDO
00283
00284
00285 Do i_band=1,n_occ
00286 Do j_band=1,n_occ
00287 SUM_CMPX=0.0D0
00288 Do k_band=1,n_occ
00289 SUM_CMPX
00290 + =SUM_CMPX
00291 + +RC_MAT(i_band,k_band)
00292 + *CONJG(R_MAT(j_band,k_band))
00293 ENDDO
00294 UNITARY(i_band,j_band,ik)=SUM_CMPX
00295 ENDDO
00296 ENDDO
00297 ENDDO
00298 U_ORG(:,:,:)=UNITARY(:,:,:)
00299 U_OLD(:,:,:)=UNITARY(:,:,:)
00300 ELSE
00301 Do ik=1,NTK
00302 DO i_band=1,n_occ
00303 DO j_band=1,n_occ
00304 SUM_CMPX=0.0D0
00305 Do k_band=1,n_occ
00306 SUM_CMPX
00307 + =SUM_CMPX
00308 + +U_OLD(i_band,k_band,ik)
00309 + *U_ORG(k_band,j_band,ik)
00310 ENDDO
00311 UNITARY(i_band,j_band,ik)=SUM_CMPX
00312 ENDDO
00313 ENDDO
00314 ENDDO
00315 U_OLD(:,:,:)=UNITARY(:,:,:)
00316 ENDIF
00317 RETURN
00318 END
00319
00320 SUBROUTINE UPDATE_OVERLAP(NTK,NTB,n_occ,NB,KPT,UNITARY,OVERLAP)
00321 implicit none
00322 integer::NTK,NTB,n_occ,NB
00323 integer::KPT(NTK,NB)
00324 integer::ik,ik_ib,ibvec,i_band,j_band,k_band
00325 complex(8)::UNITARY(n_occ,n_occ,NTK)
00326 complex(8)::OVERLAP(NTB,NTB,NTK,NB)
00327 complex(8)::U_MAT(n_occ,n_occ)
00328 complex(8)::SUM_CMPX
00329
00330 Do ik=1,NTK
00331 Do ibvec=1,NB
00332 ik_ib=KPT(ik,ibvec)
00333 Do i_band=1,n_occ
00334 Do j_band=1,n_occ
00335 SUM_CMPX=0.0D0
00336 Do k_band=1,n_occ
00337 SUM_CMPX
00338 + =SUM_CMPX
00339 + +OVERLAP(i_band,k_band,ik,ibvec)
00340 + *UNITARY(k_band,j_band,ik_ib)
00341 ENDDO
00342 U_MAT(i_band,j_band)=SUM_CMPX
00343 ENDDO
00344 ENDDO
00345 Do i_band=1,n_occ
00346 Do j_band=1,n_occ
00347 SUM_CMPX=0.0D0
00348 Do k_band=1,n_occ
00349 SUM_CMPX
00350 + =SUM_CMPX
00351 + +CONJG(UNITARY(k_band,i_band,ik))
00352 + *U_MAT(k_band,j_band)
00353 ENDDO
00354 OVERLAP(i_band,j_band,ik,ibvec)=SUM_CMPX
00355 ENDDO
00356 ENDDO
00357 ENDDO
00358 ENDDO
00359 RETURN
00360 END
00361
00362 SUBROUTINE A_MAT_UPDATE(NTK,NTB,NGAUSS,n_occ,UNITARY,A_MAT)
00363 implicit none
00364 integer::NTK,NTB,NGAUSS,n_occ
00365 integer::ik,i_band,j_band,k_band
00366 complex(8)::UNITARY(n_occ,n_occ,NTK)
00367 complex(8)::A_MAT(NTB,NGAUSS,NTK)
00368 complex(8)::A_TMP(NTB,NGAUSS)
00369 complex(8)::SUM_CMPX
00370
00371 Do ik=1,NTK
00372 A_TMP(:,:)=0.0D0
00373 A_TMP(:,:)=A_MAT(:,:,ik)
00374 Do i_band=1,n_occ
00375 Do j_band=1,n_occ
00376 SUM_CMPX=0.0D0
00377 Do k_band=1,n_occ
00378 SUM_CMPX
00379 + =SUM_CMPX
00380 + +A_TMP(i_band,k_band)
00381 + *UNITARY(k_band,j_band,ik)
00382 ENDDO
00383 A_MAT(i_band,j_band,ik)=SUM_CMPX
00384 ENDDO
00385 ENDDO
00386 ENDDO
00387 RETURN
00388 END
00389
00390 SUBROUTINE CALC_OVERLAP(NTK,NTB,NTG,nvx,nvy,nvz,ik,ik_ib,NG0,
00391 + KG0,N_BAND,SHIFT_b,C0_BRA,C0_KET,OVERLAP_TMP,N_BAND_BOTTOM)
00392 implicit none
00393 integer::NTK,NTB,NTG,nvx,nvy,nvz
00394 integer::NG0(NTK),KG0(3,NTG,NTK),N_BAND(NTK)
00395 integer::SHIFT_b(3),N_BAND_BOTTOM(NTK)
00396 integer::PULLBACK(NTG)
00397 integer::MAP(-nvx:nvx,-nvy:nvy,-nvz:nvz)
00398 integer::ik,ik_ib,ig,ig1,ig2,i_band,j_band
00399 complex(8)::C0_BRA(0:NTG,NTB)
00400 complex(8)::C0_KET(0:NTG,NTB)
00401 complex(8)::OVERLAP_TMP(NTB,NTB)
00402 complex(8)::ELEMENT_OVERLAP
00403
00404 MAP(:,:,:)=0
00405 do ig=1,NG0(ik)
00406 MAP(KG0(1,ig,ik)+SHIFT_b(1),
00407 + KG0(2,ig,ik)+SHIFT_b(2),
00408 + KG0(3,ig,ik)+SHIFT_b(3))=ig
00409 enddo
00410
00411 PULLBACK(:)=0
00412 do ig=1,NG0(ik_ib)
00413 PULLBACK(ig)
00414 + = MAP(KG0(1,ig,ik_ib)
00415 + ,KG0(2,ig,ik_ib)
00416 + ,KG0(3,ig,ik_ib))
00417 enddo
00418
00419 OVERLAP_TMP(:,:)=0.0D0
00420 do i_band=1,N_BAND(ik)
00421 do j_band=1,N_BAND(ik_ib)
00422 ELEMENT_OVERLAP=0.0D0
00423 do ig2=1,NG0(ik_ib)
00424 ig1=PULLBACK(ig2)
00425 ELEMENT_OVERLAP
00426 + = ELEMENT_OVERLAP
00427 + + CONJG(C0_BRA(ig1,i_band+N_BAND_BOTTOM(ik)))
00428 + * C0_KET(ig2,j_band+N_BAND_BOTTOM(ik_ib))
00429 enddo
00430 OVERLAP_TMP(i_band,j_band)=ELEMENT_OVERLAP
00431 enddo
00432 enddo
00433 RETURN
00434 END
00435
00436 subroutine kcheck(ktmp,RWtmp)
00437 implicit none
00438 real(8),intent(inout)::ktmp(3)
00439 integer,intent(out)::RWtmp(3)
00440 real(8),parameter::dlt_BZ=1.0d-6
00441
00442
00443
00444 if(ktmp(1)>1.50d0+dlt_BZ)then
00445 ktmp(1)=ktmp(1)-2.0d0;RWtmp(1)=-2
00446 endif
00447 if(ktmp(1)>0.50d0+dlt_BZ)then
00448 ktmp(1)=ktmp(1)-1.0d0;RWtmp(1)=-1
00449 endif
00450 if(ktmp(1)<=-1.50d0+dlt_BZ)then
00451 ktmp(1)=ktmp(1)+2.0d0;RWtmp(1)=2
00452 endif
00453 if(ktmp(1)<=-0.50d0+dlt_BZ)then
00454 ktmp(1)=ktmp(1)+1.0d0;RWtmp(1)=1
00455 endif
00456
00457 if(ktmp(2)>1.50d0+dlt_BZ)then
00458 ktmp(2)=ktmp(2)-2.0d0;RWtmp(2)=-2
00459 endif
00460 if(ktmp(2)>0.50d0+dlt_BZ)then
00461 ktmp(2)=ktmp(2)-1.0d0;RWtmp(2)=-1
00462 endif
00463 if(ktmp(2)<=-1.50d0+dlt_BZ)then
00464 ktmp(2)=ktmp(2)+2.0d0;RWtmp(2)=2
00465 endif
00466 if(ktmp(2)<=-0.50d0+dlt_BZ)then
00467 ktmp(2)=ktmp(2)+1.0d0;RWtmp(2)=1
00468 endif
00469
00470 if(ktmp(3)>1.50d0+dlt_BZ)then
00471 ktmp(3)=ktmp(3)-2.0d0;RWtmp(3)=-2
00472 endif
00473 if(ktmp(3)>0.50d0+dlt_BZ)then
00474 ktmp(3)=ktmp(3)-1.0d0;RWtmp(3)=-1
00475 endif
00476 if(ktmp(3)<=-1.50d0+dlt_BZ)then
00477 ktmp(3)=ktmp(3)+2.0d0;RWtmp(3)=2
00478 endif
00479 if(ktmp(3)<=-0.50d0+dlt_BZ)then
00480 ktmp(3)=ktmp(3)+1.0d0;RWtmp(3)=1
00481 endif
00482
00483
00484
00485 return
00486 end
00487
00488
00489 subroutine kcheck_trs(ktmp,RWtmp)
00490 implicit none
00491 real(8),intent(inout)::ktmp(3)
00492 integer,intent(out)::RWtmp(3)
00493 real(8),parameter::dlt_BZ=-1.0d-6
00494
00495
00496
00497 if(ktmp(1)>=1.50d0+dlt_BZ)then
00498 ktmp(1)=ktmp(1)-2.0d0;RWtmp(1)=-2
00499 endif
00500 if(ktmp(1)>=0.50d0+dlt_BZ)then
00501 ktmp(1)=ktmp(1)-1.0d0;RWtmp(1)=-1
00502 endif
00503 if(ktmp(1)<-1.50d0+dlt_BZ)then
00504 ktmp(1)=ktmp(1)+2.0d0;RWtmp(1)=2
00505 endif
00506 if(ktmp(1)<-0.50d0+dlt_BZ)then
00507 ktmp(1)=ktmp(1)+1.0d0;RWtmp(1)=1
00508 endif
00509
00510 if(ktmp(2)>=1.50d0+dlt_BZ)then
00511 ktmp(2)=ktmp(2)-2.0d0;RWtmp(2)=-2
00512 endif
00513 if(ktmp(2)>=0.50d0+dlt_BZ)then
00514 ktmp(2)=ktmp(2)-1.0d0;RWtmp(2)=-1
00515 endif
00516 if(ktmp(2)<-1.50d0+dlt_BZ)then
00517 ktmp(2)=ktmp(2)+2.0d0;RWtmp(2)=2
00518 endif
00519 if(ktmp(2)<-0.50d0+dlt_BZ)then
00520 ktmp(2)=ktmp(2)+1.0d0;RWtmp(2)=1
00521 endif
00522
00523 if(ktmp(3)>=1.50d0+dlt_BZ)then
00524 ktmp(3)=ktmp(3)-2.0d0;RWtmp(3)=-2
00525 endif
00526 if(ktmp(3)>=0.50d0+dlt_BZ)then
00527 ktmp(3)=ktmp(3)-1.0d0;RWtmp(3)=-1
00528 endif
00529 if(ktmp(3)<-1.50d0+dlt_BZ)then
00530 ktmp(3)=ktmp(3)+2.0d0;RWtmp(3)=2
00531 endif
00532 if(ktmp(3)<-0.50d0+dlt_BZ)then
00533 ktmp(3)=ktmp(3)+1.0d0;RWtmp(3)=1
00534 endif
00535
00536
00537
00538 return
00539 end
00540
00541 subroutine kcheck_old(ktmp,RWtmp)
00542 implicit none
00543 real(8),intent(inout)::ktmp(3)
00544 integer,intent(out)::RWtmp(3)
00545 if(ktmp(1)>=0.5d0)then
00546 ktmp(1)=ktmp(1)-1.0d0;RWtmp(1)=-1
00547 endif
00548 if(ktmp(1)<-0.5d0)then
00549 ktmp(1)=ktmp(1)+1.0d0;RWtmp(1)=1
00550 endif
00551 if(ktmp(2)>=0.5d0)then
00552 ktmp(2)=ktmp(2)-1.0d0;RWtmp(2)=-1
00553 endif
00554 if(ktmp(2)<-0.5d0)then
00555 ktmp(2)=ktmp(2)+1.0d0;RWtmp(2)=1
00556 endif
00557 if(ktmp(3)>=0.5d0)then
00558 ktmp(3)=ktmp(3)-1.0d0;RWtmp(3)=-1
00559 endif
00560 if(ktmp(3)<-0.5d0)then
00561 ktmp(3)=ktmp(3)+1.0d0;RWtmp(3)=1
00562 endif
00563 return
00564 end
00565
00566 subroutine make_KG0(NTG,b1,b2,b3,Gcut,q1,q2,q3,KG0,NG)
00567 implicit none
00568 integer,intent(in)::NTG
00569 real(8),intent(in)::b1(3),b2(3),b3(3)
00570 real(8),intent(in)::Gcut,q1,q2,q3
00571 integer,intent(out)::KG0(3,NTG),NG
00572 integer::igL,igL1,igL2,igL3
00573 real(8)::qgL(3),qgL2
00574 integer,parameter::NGL1=150
00575 integer,parameter::NGL2=150
00576 integer,parameter::NGL3=150
00577
00578 igL=0
00579 do igL1=-NGL1,NGL1
00580 do igL2=-NGL2,NGL2
00581 do igL3=-NGL3,NGL3
00582 qgL(:)=(q1+dble(igL1))*b1(:)+(q2+dble(igL2))*b2(:)
00583 + +(q3+dble(igL3))*b3(:)
00584 qgL2=qgL(1)**2+qgL(2)**2+qgL(3)**2
00585 if(qgL2<=Gcut) then
00586 igL=igL+1
00587 KG0(1,igL)=igL1
00588 KG0(2,igL)=igL2
00589 KG0(3,igL)=igL3
00590 endif
00591 enddo
00592 enddo
00593 enddo
00594 NG=igL
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604 RETURN
00605 END
00606
00607 subroutine make_C0(NTG,NTB,itrs,NG,KGtmp,RWtmp,rginvtmp,pgtmp,
00608 + nnp,L1,L2,L3,packtmp,C0irr,C0_K)
00609 implicit none
00610
00611 integer::NTG,NTB,itrs,NG,L1,L2,L3,nnp
00612 integer::KGtmp(3,NTG),RWtmp(3),pgtmp(3)
00613 integer::packtmp(-L1:L1,-L2:L2,-L3:L3)
00614 integer::ig,jg,i1,i2,i3,j1,j2,j3,k1,k2,k3
00615 real(8)::rginvtmp(3,3),phase
00616 complex(8)::C0irr(NTG,NTB),C0_K(NTG,NTB)
00617 complex(8)::pf
00618 real(8),parameter::pi=dacos(-1.0d0)
00619 real(8),parameter::tpi=2.0d0*pi
00620 complex(8),parameter::ci=(0.0D0,1.0D0)
00621
00622 C0_K(:,:)=0.0d0
00623 select case(itrs)
00624 case(1)
00625 do ig=1,NG
00626 i1=KGtmp(1,ig); j1=KGtmp(2,ig); k1=KGtmp(3,ig)
00627 i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3)
00628 i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2
00629 + +int(rginvtmp(1,3))*k2
00630 j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2
00631 + +int(rginvtmp(2,3))*k2
00632 k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2
00633 + +int(rginvtmp(3,3))*k2
00634 jg=packtmp(i3,j3,k3)
00635 phase=tpi*(dble(i1)*dble(pgtmp(1))
00636 + +dble(j1)*dble(pgtmp(2))
00637 + +dble(k1)*dble(pgtmp(3)))
00638 pf=exp(-ci*phase/dble(nnp))
00639 C0_K(ig,:)=C0irr(jg,:)*pf
00640 enddo
00641 case(-1)
00642 do ig=1,NG
00643 i1=KGtmp(1,ig); j1=KGtmp(2,ig); k1=KGtmp(3,ig)
00644 i2=i1+RWtmp(1); j2=j1+RWtmp(2); k2=k1+RWtmp(3)
00645 i3=int(rginvtmp(1,1))*i2+int(rginvtmp(1,2))*j2
00646 + +int(rginvtmp(1,3))*k2
00647 j3=int(rginvtmp(2,1))*i2+int(rginvtmp(2,2))*j2
00648 + +int(rginvtmp(2,3))*k2
00649 k3=int(rginvtmp(3,1))*i2+int(rginvtmp(3,2))*j2
00650 + +int(rginvtmp(3,3))*k2
00651 jg=packtmp(i3,j3,k3)
00652 phase=tpi*(dble(i1)*dble(pgtmp(1))
00653 + +dble(j1)*dble(pgtmp(2))
00654 + +dble(k1)*dble(pgtmp(3)))
00655 pf=exp(-ci*phase/dble(nnp))
00656 C0_K(ig,:)=C0irr(jg,:)*pf
00657 enddo
00658 C0_K(:,:)=conjg(C0_K(:,:))
00659 end select
00660
00661
00662
00663
00664
00665
00666
00667 return
00668 end
00669
00670 SUBROUTINE OUTER_PRODUCT(vec_x,vec_y,vec_z)
00671 implicit none
00672 real(8)::vec_x(3),vec_y(3),vec_z(3)
00673 vec_z(1)=vec_x(2)*vec_y(3)-vec_x(3)*vec_y(2)
00674 vec_z(2)=vec_x(3)*vec_y(1)-vec_x(1)*vec_y(3)
00675 vec_z(3)=vec_x(1)*vec_y(2)-vec_x(2)*vec_y(1)
00676 RETURN
00677 END
00678
00679 subroutine BZcheck(ktmp,kflg)
00680 implicit none
00681 real(8),intent(in)::ktmp(3)
00682 real(8)::kx,ky,kz
00683 integer,intent(out)::kflg
00684 kx=-ktmp(1)
00685 ky=-ktmp(2)
00686 kz=-ktmp(3)
00687 kflg=0
00688 if(kx<=-0.5d0) kflg=1
00689 if(kx>0.5d0) kflg=1
00690 if(ky<=-0.5d0) kflg=1
00691 if(ky>0.5d0) kflg=1
00692 if(kz<=-0.5d0) kflg=1
00693 if(kz>0.5d0) kflg=1
00694 return
00695 end
00696
00697
00698 subroutine calc_NB_start_end(NK,NB,EIG,EL,EU,NB_start,NB_end)
00699 implicit none
00700 integer::NK,NB
00701 real(8)::EIG(NB,NK)
00702 real(8)::EL,EU
00703 integer::NB_start,NB_end
00704 integer::Ein(NB,NK)
00705 integer::ik,ib
00706
00707 Ein=0
00708 do ik=1,NK
00709 do ib=1,NB
00710 if(EL<=EIG(ib,ik).and.EIG(ib,ik)<=EU)then
00711 Ein(ib,ik)=ib
00712 endif
00713 enddo
00714 enddo
00715 NB_end=maxval(Ein)
00716 do ik=1,NK
00717 do ib=1,NB
00718 if(Ein(ib,ik)==0) Ein(ib,ik)=10000
00719 enddo
00720 enddo
00721 NB_start=minval(Ein)
00722 return
00723 end
00724
00725
00726 subroutine make_shift_vector(NTK,NB,ik,ibvec,SK0,b_LAT,
00727 +SHIFT_b,ik_ib)
00728 implicit none
00729 integer::NTK,NB,ik,ibvec
00730 real(8)::SK0(3,NTK)
00731 real(8)::b_LAT(3,NB)
00732 integer::SHIFT_b(3)
00733 integer::ik_ib
00734 real(8)::SL(3)
00735 integer::jk
00736 real(8),parameter::dlt_BZ=1.0d-6
00737
00738 SL(:)=SK0(:,ik)+b_LAT(:,ibvec)
00739
00740
00741 SHIFT_b(:)=0
00742 IF(SL(1)<=-0.5D0+dlt_BZ)THEN
00743 SL(1)=SL(1)+1.0D0
00744 SHIFT_b(1)=-1
00745 ENDIF
00746 IF(SL(2)<=-0.5D0+dlt_BZ)THEN
00747 SL(2)=SL(2)+1.0D0
00748 SHIFT_b(2)=-1
00749 ENDIF
00750 IF(SL(3)<=-0.5D0+dlt_BZ)THEN
00751 SL(3)=SL(3)+1.0D0
00752 SHIFT_b(3)=-1
00753 ENDIF
00754 IF(SL(1)>0.5D0+dlt_BZ)THEN
00755 SL(1)=SL(1)-1.0D0
00756 SHIFT_b(1)=1
00757 ENDIF
00758 IF(SL(2)>0.5D0+dlt_BZ)THEN
00759 SL(2)=SL(2)-1.0D0
00760 SHIFT_b(2)=1
00761 ENDIF
00762 IF(SL(3)>0.5D0+dlt_BZ)THEN
00763 SL(3)=SL(3)-1.0D0
00764 SHIFT_b(3)=1
00765 ENDIF
00766
00767
00768 do jk=1,NTK
00769 IF(ABS(SK0(1,jk)-SL(1))<=1.d-4.AND.
00770 + ABS(SK0(2,jk)-SL(2))<=1.d-4.AND.
00771 + ABS(SK0(3,jk)-SL(3))<=1.d-4)THEN
00772 ik_ib=jk
00773 ENDIF
00774 enddo
00775 return
00776 end
00777
00778
00779
00780
00781
00782
00783
00784 subroutine make_amat(nigs,NTB,NTG,NG0,SK0,KG0,b1,b2,b3,
00785 + TAU_GAUSS,ALPHA_GAUSS,loc_x,loc_y,loc_z,LGAUSS,MGAUSS,
00786 + N_BAND,C0,N_BAND_BTM,NORM_GAUSS,VOLUME,A_TMP)
00787 implicit none
00788 integer::nigs,NTB,NTG
00789 integer::NG0
00790 real(8)::SK0(3)
00791 integer::KG0(3,NTG)
00792 real(8)::b1(3),b2(3),b3(3)
00793 real(8)::VOLUME
00794 real(8)::TAU_GAUSS(3,nigs)
00795 real(8)::ALPHA_GAUSS(nigs)
00796 real(8)::NORM_GAUSS(nigs)
00797 real(8)::loc_x(3,nigs),loc_y(3,nigs),loc_z(3,nigs)
00798 integer::LGAUSS(nigs),MGAUSS(nigs)
00799 integer::N_BAND,N_BAND_BTM
00800 complex(8)::C0(NTG,NTB)
00801 complex(8)::A_TMP(NTB,nigs)
00802
00803 real(8),parameter::pi=dacos(-1.0d0)
00804 complex(8),parameter::ci=(0.0D0,1.0D0)
00805 integer::j_gauss,ig,i_band,jL,jM
00806 real(8)::PHASE,ZETA1,ZETA2,YLM
00807 real(8)::VEC_KG(3,NTG)
00808 real(8)::KG2(NTG)
00809 complex(8)::PF1(NTG)
00810 real(8)::PF2(NTG)
00811 real(8)::ktmp(3),tmp(3)
00812 complex(8)::SUM_CMPX
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823 A_TMP=0.0d0
00824 do j_gauss=1,nigs
00825
00826 VEC_KG(:,:)=0.0D0
00827 PF1(:)=0.0D0
00828 PF2(:)=0.0D0
00829 KG2(:)=0.0D0
00830 do ig=1,NG0
00831 VEC_KG(:,ig)=(SK0(1)+DBLE(KG0(1,ig)))*b1(:)
00832 + +(SK0(2)+DBLE(KG0(2,ig)))*b2(:)
00833 + +(SK0(3)+DBLE(KG0(3,ig)))*b3(:)
00834 KG2(ig)=VEC_KG(1,ig)**2+VEC_KG(2,ig)**2+VEC_KG(3,ig)**2
00835 PHASE=VEC_KG(1,ig)*TAU_GAUSS(1,j_gauss)
00836 + +VEC_KG(2,ig)*TAU_GAUSS(2,j_gauss)
00837 + +VEC_KG(3,ig)*TAU_GAUSS(3,j_gauss)
00838 PF1(ig)=EXP(-ci*PHASE )
00839 PF2(ig)=EXP(-KG2(ig)/4.0D0/ALPHA_GAUSS(j_gauss))
00840 enddo
00841
00842
00843
00844
00845
00846
00847 do ig=1,NG0
00848 ktmp(:)=VEC_KG(:,ig)
00849 tmp(1)=ktmp(1)*loc_x(1,j_gauss)
00850 + +ktmp(2)*loc_x(2,j_gauss)
00851 + +ktmp(3)*loc_x(3,j_gauss)
00852 tmp(2)=ktmp(1)*loc_y(1,j_gauss)
00853 + +ktmp(2)*loc_y(2,j_gauss)
00854 + +ktmp(3)*loc_y(3,j_gauss)
00855 tmp(3)=ktmp(1)*loc_z(1,j_gauss)
00856 + +ktmp(2)*loc_z(2,j_gauss)
00857 + +ktmp(3)*loc_z(3,j_gauss)
00858 VEC_KG(:,ig)=tmp(:)
00859 enddo
00860
00861
00862
00863
00864 jL=LGAUSS(j_gauss)
00865 jM=MGAUSS(j_gauss)
00866
00867 do i_band=1,N_BAND
00868 SUM_CMPX=0.0D0
00869
00870 IF(jL==1)THEN
00871 do ig=1,NG0
00872 ZETA1=(pi/ALPHA_GAUSS(j_gauss))**1.5D0
00873 YLM=DSQRT(1.0D0/4.0D0/pi)
00874 SUM_CMPX
00875 + =SUM_CMPX
00876 + +CONJG(C0(ig,i_band+N_BAND_BTM))
00877 + *PF1(ig)
00878 + *YLM
00879 + *NORM_GAUSS(j_gauss)
00880 + *ZETA1*PF2(ig)
00881 enddo
00882 ENDIF
00883
00884 IF(jL==2.AND.jM==2)THEN
00885 do ig=1,NG0
00886 ZETA1=(pi/ALPHA_GAUSS(j_gauss))**1.5D0
00887 ZETA2=1.0D0/2.0D0/ALPHA_GAUSS(j_gauss)
00888 YLM=DSQRT(3.0D0/4.0D0/pi)*VEC_KG(2,ig)
00889 SUM_CMPX
00890 + =SUM_CMPX
00891 + +CONJG(C0(ig,i_band+N_BAND_BTM))
00892 + *PF1(ig)
00893 + *(-ci)*YLM
00894 + *NORM_GAUSS(j_gauss)
00895 + *ZETA1*ZETA2*PF2(ig)
00896 enddo
00897 ENDIF
00898
00899 IF(jL==2.AND.jM==3)THEN
00900 do ig=1,NG0
00901 ZETA1=(pi/ALPHA_GAUSS(j_gauss))**1.5D0
00902 ZETA2=1.0D0/2.0D0/ALPHA_GAUSS(j_gauss)
00903 YLM=DSQRT(3.0D0/4.0D0/pi)*VEC_KG(3,ig)
00904 SUM_CMPX
00905 + =SUM_CMPX
00906 + +CONJG(C0(ig,i_band+N_BAND_BTM))
00907 + *PF1(ig)
00908 + *(-ci)*YLM
00909 + *NORM_GAUSS(j_gauss)
00910 + *ZETA1*ZETA2*PF2(ig)
00911 enddo
00912 ENDIF
00913
00914 IF(jL==2.AND.jM==1)THEN
00915 do ig=1,NG0
00916 ZETA1=(pi/ALPHA_GAUSS(j_gauss))**1.5D0
00917 ZETA2=1.0D0/2.0D0/ALPHA_GAUSS(j_gauss)
00918 YLM=DSQRT(3.d0/4.d0/pi)*VEC_KG(1,ig)
00919 SUM_CMPX
00920 + =SUM_CMPX
00921 + +CONJG(C0(ig,i_band+N_BAND_BTM))
00922 + *PF1(ig)
00923 + *(-ci)*YLM
00924 + *NORM_GAUSS(j_gauss)
00925 + *ZETA1*ZETA2*PF2(ig)
00926 enddo
00927 ENDIF
00928
00929 IF(jL==3.AND.jM==1)THEN
00930 do ig=1,NG0
00931 ZETA1=(pi/ALPHA_GAUSS(j_gauss))**1.5D0
00932 ZETA2=(1.0D0/2.0D0/ALPHA_GAUSS(j_gauss))**2.0D0
00933 YLM=0.5d0*DSQRT(15.d0/pi)*VEC_KG(1,ig)*VEC_KG(2,ig)
00934 SUM_CMPX
00935 + =SUM_CMPX
00936 + +CONJG(C0(ig,i_band+N_BAND_BTM))
00937 + *PF1(ig)
00938 + *(-YLM)
00939 + *NORM_GAUSS(j_gauss)
00940 + *ZETA1*ZETA2*PF2(ig)
00941 enddo
00942 ENDIF
00943
00944 IF(jL==3.AND.jM==2)THEN
00945 do ig=1,NG0
00946 ZETA1=(pi/ALPHA_GAUSS(j_gauss))**1.5D0
00947 ZETA2=(1.0D0/2.0D0/ALPHA_GAUSS(j_gauss))**2.0D0
00948 YLM=(0.5d0)*DSQRT(15.d0/pi)*VEC_KG(2,ig)*VEC_KG(3,ig)
00949 SUM_CMPX
00950 + =SUM_CMPX
00951 + +CONJG(C0(ig,i_band+N_BAND_BTM))
00952 + *PF1(ig)
00953 + *(-YLM)
00954 + *NORM_GAUSS(j_gauss)
00955 + *ZETA1*ZETA2*PF2(ig)
00956 enddo
00957 ENDIF
00958
00959 IF(jL==3.AND.jM==3)THEN
00960 do ig=1,NG0
00961 ZETA1=(pi/ALPHA_GAUSS(j_gauss))**1.5D0
00962 ZETA2=(1.0D0/2.0D0/ ALPHA_GAUSS(j_gauss))**2.0D0
00963 YLM=(0.25d0)*DSQRT(5.d0/pi)*(3.d0*VEC_KG(3,ig)**2-KG2(ig))
00964 SUM_CMPX
00965 + =SUM_CMPX
00966 + +CONJG(C0(ig,i_band+N_BAND_BTM))
00967 + *PF1(ig)
00968 + *(-YLM)
00969 + *NORM_GAUSS(j_gauss)
00970 + *ZETA1*ZETA2*PF2(ig)
00971 enddo
00972 ENDIF
00973
00974 IF(jL==3.AND.jM==4)THEN
00975 do ig=1,NG0
00976 ZETA1=(pi/ALPHA_GAUSS(j_gauss))**1.5D0
00977 ZETA2=(1.0D0/2.0D0/ALPHA_GAUSS(j_gauss))**2.0D0
00978 YLM=(0.5d0)*DSQRT(15.d0/pi)*VEC_KG(3,ig)*VEC_KG(1,ig)
00979 SUM_CMPX
00980 + =SUM_CMPX
00981 + +CONJG(C0(ig,i_band+N_BAND_BTM))
00982 + *PF1(ig)
00983 + *(-YLM)
00984 + *NORM_GAUSS(j_gauss)
00985 + *ZETA1*ZETA2*PF2(ig)
00986 enddo
00987 ENDIF
00988
00989 IF(jL==3.AND.jM==5)THEN
00990 do ig=1,NG0
00991 ZETA1=(pi/ALPHA_GAUSS(j_gauss))**1.5D0
00992 ZETA2=(1.0D0/2.0D0/ALPHA_GAUSS(j_gauss))**2.0D0
00993 YLM=(0.25d0)*DSQRT(15.d0/pi)
00994 + *(VEC_KG(1,ig)**2-VEC_KG(2,ig)**2)
00995 SUM_CMPX
00996 + =SUM_CMPX
00997 + +CONJG(C0(ig,i_band+N_BAND_BTM))
00998 + *PF1(ig)
00999 + *(-YLM)
01000 + *NORM_GAUSS(j_gauss)
01001 + *ZETA1*ZETA2*PF2(ig)
01002 enddo
01003 ENDIF
01004
01005 A_TMP(i_band,j_gauss)=SUM_CMPX/DSQRT(VOLUME)
01006 enddo
01007 enddo
01008
01009
01010
01011
01012
01013
01014 return
01015 end
01016
01017 subroutine make_pw(C0_I,wfunc,fftwk,NG0,KG0,NTG,nwx2,nwy2,nwz2,
01018 + nfft1,nfft2,Nl123,fs,xmin,xmax,ymin,ymax,zmin,zmax,SK0,pw)
01019 use fft_3d
01020 implicit none
01021 type(fft3_struct)::fs
01022 integer::NTG,nwx2,nwy2,nwz2,nfft1,nfft2,Nl123
01023 integer::xmin,xmax,ymin,ymax,zmin,zmax
01024 integer::KG0(3,NTG)
01025 integer::NG0
01026 real(8)::fftwk(Nl123*2)
01027 real(8)::wfunc(Nl123*2)
01028 real(8)::SK0(3)
01029 complex(8)::C0_I(NTG)
01030 complex(8)::pw((xmin)*nwx2:(xmax)*nwx2-1,
01031 + (ymin)*nwy2:(ymax)*nwy2-1,
01032 + (zmin)*nwz2:(zmax)*nwz2-1)
01033 integer::ig,igb1,igb2,igb3,ind,ir,ir1,ir2,ir3
01034 integer::ir1_SHIFT,ir2_SHIFT,ir3_SHIFT,jr1,jr2,jr3
01035 real(8)::PHASE
01036 complex(8)::pf
01037 real(8),parameter::tpi=2.0d0*dacos(-1.0d0)
01038 complex(8),parameter::ci=(0.0D0,1.0D0)
01039 complex(8)::u_1D(Nl123)
01040 complex(8)::u_3D(0:nwx2-1,0:nwy2-1,0:nwz2-1)
01041
01042 u_3D=0.0D0
01043 u_1D=0.0D0
01044 wfunc=0.0d0
01045 fftwk=0.0d0
01046 do ig=1,NG0
01047 igb1=KG0(1,ig)
01048 igb2=KG0(2,ig)
01049 igb3=KG0(3,ig)
01050 igb1=MOD(nwx2+igb1,nwx2)+1
01051 igb2=MOD(nwy2+igb2,nwy2)+1
01052 igb3=MOD(nwz2+igb3,nwz2)+1
01053 ind=igb1+(igb2-1)*nfft1+(igb3-1)*nfft1*nfft2
01054 wfunc(ind)=dble(C0_I(ig))
01055 wfunc(ind+Nl123)=dimag(C0_I(ig))
01056 enddo
01057 call fft3_bw(fs,wfunc(1),fftwk(1))
01058 do ir=1,Nl123
01059 u_1D(ir)=cmplx(wfunc(ir),wfunc(ir+Nl123))
01060 enddo
01061
01062 do ir3=1,nwz2
01063 do ir2=1,nwy2
01064 do ir1=1,nwx2
01065 ir=ir1+(ir2-1)*nfft1+(ir3-1)*nfft1*nfft2
01066 u_3D(ir1-1,ir2-1,ir3-1)=u_1D(ir)
01067 enddo
01068 enddo
01069 enddo
01070
01071 do ir1=(xmin)*nwx2,(xmax)*nwx2-1
01072 do ir2=(ymin)*nwy2,(ymax)*nwy2-1
01073 do ir3=(zmin)*nwz2,(zmax)*nwz2-1
01074 PHASE=tpi*(DBLE(ir1)*SK0(1)/DBLE(nwx2)
01075 + +DBLE(ir2)*SK0(2)/DBLE(nwy2)
01076 + +DBLE(ir3)*SK0(3)/DBLE(nwz2))
01077 pf=EXP(ci*PHASE)
01078
01079 ir1_SHIFT=ir1
01080 ir2_SHIFT=ir2
01081 ir3_SHIFT=ir3
01082 DO WHILE(ir1_SHIFT<0)
01083 ir1_SHIFT=ir1_SHIFT+nwx2
01084 ENDDO
01085 DO WHILE(ir2_SHIFT<0)
01086 ir2_SHIFT=ir2_SHIFT+nwy2
01087 ENDDO
01088 DO WHILE(ir3_SHIFT<0)
01089 ir3_SHIFT=ir3_SHIFT+nwz2
01090 ENDDO
01091 jr1=MOD(ir1_SHIFT,nwx2)
01092 jr2=MOD(ir2_SHIFT,nwy2)
01093 jr3=MOD(ir3_SHIFT,nwz2)
01094
01095 pw(ir1,ir2,ir3)=pw(ir1,ir2,ir3)+u_3D(jr1,jr2,jr3)*pf
01096 enddo
01097 enddo
01098 enddo
01099
01100 return
01101 end
01102
01103 subroutine make_bf(C0_I,wfunc,fftwk,NG0,KG0,NTG,nwx2,nwy2,nwz2,
01104 + nfft1,nfft2,Nl123,fs,pw)
01105 use fft_3d
01106 implicit none
01107 type(fft3_struct)::fs
01108 integer::NTG,nwx2,nwy2,nwz2,nfft1,nfft2,Nl123
01109 integer::KG0(3,NTG)
01110 integer::NG0
01111 real(8)::fftwk(Nl123*2)
01112 real(8)::wfunc(Nl123*2)
01113 complex(8)::C0_I(NTG)
01114 integer::ig,igb1,igb2,igb3,ind,ir,ir1,ir2,ir3
01115 complex(8)::u_1D(Nl123)
01116 real(8)::pw(0:nwx2-1,0:nwy2-1,0:nwz2-1)
01117
01118 u_1D=0.0d0
01119 wfunc=0.0d0
01120 fftwk=0.0d0
01121 do ig=1,NG0
01122 igb1=KG0(1,ig)
01123 igb2=KG0(2,ig)
01124 igb3=KG0(3,ig)
01125 igb1=MOD(nwx2+igb1,nwx2)+1
01126 igb2=MOD(nwy2+igb2,nwy2)+1
01127 igb3=MOD(nwz2+igb3,nwz2)+1
01128 ind=igb1+(igb2-1)*nfft1+(igb3-1)*nfft1*nfft2
01129 wfunc(ind)=dble(C0_I(ig))
01130 wfunc(ind+Nl123)=dimag(C0_I(ig))
01131 enddo
01132 call fft3_bw(fs,wfunc(1),fftwk(1))
01133 do ir=1,Nl123
01134 u_1D(ir)=cmplx(wfunc(ir),wfunc(ir+Nl123))
01135 enddo
01136
01137 do ir3=1,nwz2
01138 do ir2=1,nwy2
01139 do ir1=1,nwx2
01140 ir=ir1+(ir2-1)*nfft1+(ir3-1)*nfft1*nfft2
01141 pw(ir1-1,ir2-1,ir3-1)=abs(u_1D(ir))**2
01142 enddo
01143 enddo
01144 enddo
01145
01146 return
01147 end
01148
01149 subroutine search_ik(NTK,SK0,ktmp,ik)
01150 implicit none
01151 integer::NTK
01152 real(8)::SK0(3,NTK)
01153 real(8)::ktmp(3)
01154 integer::iik,ik
01155 real(8),parameter::dlt_BZ=1.0d-5
01156
01157 ik=0
01158 do iik=1,NTK
01159
01160 if( abs(SK0(1,iik)-ktmp(1))<dlt_BZ.and.
01161 + abs(SK0(2,iik)-ktmp(2))<dlt_BZ.and.
01162 + abs(SK0(3,iik)-ktmp(3))<dlt_BZ )then
01163 ik=iik
01164 endif
01165 enddo
01166 if(ik==0)then
01167 write(6,*)'Failed search ik:'
01168 write(6,*)'There exits no k pts to be calculated:'
01169 write(6,*)'Then, stop'
01170 stop
01171 endif
01172
01173 return
01174 end
01175
01176
01177
01178 subroutine search_vis_range(ktmp,amin,amax,bmin,bmax,cmin,cmax)
01179 implicit none
01180 real(8),intent(in)::ktmp(3)
01181 integer,intent(out)::amin,amax,bmin,bmax,cmin,cmax
01182 real(8)::denom
01183 real(8),parameter::dlt_BZ=1.0d-5
01184 amin=0; amax=0; bmin=0; bmax=0; cmin=0; cmax=0
01185
01186 if(dabs(ktmp(1))<dlt_BZ)then
01187 denom=1.0d0
01188 else
01189 denom=ktmp(1)
01190 endif
01191 amax=nint(1.0d0/denom)
01192
01193 if(dabs(ktmp(2))<dlt_BZ)then
01194 denom=1.0d0
01195 else
01196 denom=ktmp(2)
01197 endif
01198 bmax=nint(1.0d0/denom)
01199
01200 if(dabs(ktmp(3))<dlt_BZ)then
01201 denom=1.0d0
01202 else
01203 denom=ktmp(3)
01204 endif
01205 cmax=nint(1.0d0/denom)
01206
01207 return
01208 end
01209
01210
01211
01212 subroutine phase_fix_bf(nwx2,nwy2,nwz2,xmin,xmax,ymin,ymax,
01213 + zmin,zmax,pw)
01214 implicit none
01215 integer,intent(in)::nwx2,nwy2,nwz2,xmin,xmax,ymin,ymax,zmin,zmax
01216 complex(8),intent(inout)::pw((xmin)*nwx2:(xmax)*nwx2-1,
01217 + (ymin)*nwy2:(ymax)*nwy2-1,
01218 + (zmin)*nwz2:(zmax)*nwz2-1)
01219 integer::ir1,ir2,ir3
01220 real(8)::x,y,r,rmax
01221 complex(8)::glb_phase
01222 complex(8),parameter::ci=(0.0d0,1.0d0)
01223
01224 rmax=1.0d-6
01225 do ir1=(xmin)*nwx2,(xmax)*nwx2-1
01226 do ir2=(ymin)*nwy2,(ymax)*nwy2-1
01227 do ir3=(zmin)*nwz2,(zmax)*nwz2-1
01228 x=dble(pw(ir1,ir2,ir3))
01229 y=dimag(pw(ir1,ir2,ir3))
01230 r=dsqrt(x**2+y**2)
01231 if(r>rmax)then
01232 glb_phase=(x-ci*y)/r
01233 rmax=r
01234 endif
01235 enddo
01236 enddo
01237 enddo
01238
01239 pw=pw*glb_phase
01240
01241 return
01242 end