wannier_sub.F

Go to the documentation of this file.
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)) !modified 2007 6 26 
00023           WF_T(i_band,j_band,ibvec)          
00024      +   =OVERLAP(i_band,j_band,ik,ibvec)
00025      +   /OVERLAP(j_band,j_band,ik,ibvec) !modified 2007 6 26 
00026      +   *WF_CHARGE(j_band,ik,ibvec) !modified 2007 6 26 
00027          ENDDO !ibvec            
00028         ENDDO !j_band           
00029        ENDDO !i_band           
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 !ibvec            
00043          GRADIENT(i_band,j_band,ik)=SUM_CMPX/CMPLX(NTK) 
00044         ENDDO !j_band           
00045        ENDDO !i_band           
00046       ENDDO !ik       
00047 !-- block diagonal constraint --
00048 !      G(:,:,:)=GRADIENT(:,:,:) 
00049 !      GRADIENT=0.0d0 
00050 !      GRADIENT(1:2,1:2,:)=G(1:2,1:2,:) 
00051 !      GRADIENT(3:4,3:4,:)=G(3:4,3:4,:) 
00052 !      GRADIENT(5:6,5:6,:)=G(5:6,5:6,:) 
00053 !-- Tr(G^2) --                 
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!j_band   
00065        ENDDO!i_band   
00066        TR_GRAD(ik)=DBLE(SUM_CMPX)          
00067        SUM_REAL=SUM_REAL+TR_GRAD(ik)          
00068       ENDDO !ik         
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 !(1) WANNIER POSITION
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!ibvec            
00100       ENDDO!ik       
00101       WF_POSITION(:,i_band)=-VEC_SUM_REAL(:)/DBLE(NTK)   
00102       ENDDO!i_band           
00103 
00104 !(2) WANNIER CHARGE
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!i_band           
00118       ENDDO!ibvec            
00119       ENDDO!ik       
00120 
00121 !(3-a) OMEGA_I 
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!j_band           
00133       ENDDO!i_band           
00134          SUM_REAL1 
00135      + = SUM_REAL1 
00136      + + WEIGHT_b(ibvec) 
00137      + * ( DBLE(n_occ) - SUM_REAL2 )       
00138       ENDDO!ibvec            
00139       ENDDO!ik       
00140       OMEGA_I = SUM_REAL1 / DBLE(NTK)        
00141 
00142 !(3-b) OMEGA_OD 
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!j_band           
00156       ENDDO!i_band           
00157          SUM_REAL1 
00158      + = SUM_REAL1 
00159      + + WEIGHT_b(ibvec) 
00160      + * SUM_REAL2                  
00161       ENDDO!ibvec            
00162       ENDDO!ik       
00163       OMEGA_OD = SUM_REAL1 / DBLE(NTK)        
00164 
00165 !(3-c) OMEGA_D 
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!i_band           
00176          SUM_REAL1 
00177      + = SUM_REAL1 
00178      + + WEIGHT_b(ibvec) 
00179      + * SUM_REAL2                  
00180       ENDDO!ibvec            
00181       ENDDO!ik       
00182       OMEGA_D = SUM_REAL1 / DBLE(NTK)        
00183 
00184 !(4) OMEGA_I,OD AND OMEGA_D OF EACH BAND 
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!ibvec            
00201         ENDDO!ik       
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!i_band           
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 !STEP(1) diag D_MAT
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 !STEP(2) CONSTRUCT DIAGONAL MATRICES
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 !STEP(3) CONSTRUCT RC & RS MATRIX 
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 !STEP(4) CONSTRUCT RC MATRIX 
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 !STEP(5) CONSTRUCT U MATRIX 
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!ik === UNITARY MATRIX GENERATION ROOP ===     
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!j_band 
00313         ENDDO!i_band 
00314        ENDDO!ik 
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 !ibvec                   
00358       ENDDO !ik          
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 !j_band
00385        ENDDO !i_band
00386       ENDDO !ik          
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 !MAP FUNCTION
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 !PULLBACK 
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 !INTERSTATE MATRIX 
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!IG2 
00430         OVERLAP_TMP(i_band,j_band)=ELEMENT_OVERLAP      
00431        enddo!j_band     
00432       enddo!i_band     
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 !      write(6,'(a8,3f15.10)') 'before',ktmp(1),ktmp(2),ktmp(3)
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 !      write(6,'(a8,3f15.10)') 'after',ktmp(1),ktmp(2),ktmp(3)
00484 !--
00485       return
00486       end 
00487 !--
00488 !20170316 
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 !      write(6,'(a8,3f15.10)') 'before',ktmp(1),ktmp(2),ktmp(3)
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 !      write(6,'(a8,3f15.10)') 'after',ktmp(1),ktmp(2),ktmp(3)
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!100!300
00575       integer,parameter::NGL2=150!100!300 
00576       integer,parameter::NGL3=150!100!300
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 !write(6,*) "maxabs KG0_1=",maxval(abs(KG0(1,:))) 
00597 !write(6,*) "maxabs KG0_2=",maxval(abs(KG0(2,:))) 
00598 !write(6,*) "maxabs KG0_3=",maxval(abs(KG0(3,:))) 
00599 !--
00600 !      do ig=1,NG 
00601 !       write(6,'(4i8)') ig,KG0(:,ig)
00602 !      enddo 
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)!=== not time-reversal ===*      
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!ig 
00641       case(-1)!=== time-reversal ===*      
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!ig 
00658       C0_K(:,:)=conjg(C0_K(:,:))  
00659       end select 
00660 !--
00661 !      SUM_CMPX=0.0d0 
00662 !      do ig=1,NG 
00663 !       SUM_CMPX=SUM_CMPX+CONJG(C0_K(ig))*C0_K(ig) 
00664 !      enddo 
00665 !      write(6,'(a,x,2f15.10)') 'sumC0_K',SUM_CMPX 
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 !set NB_start,NB_end 
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 !SHIFT VECTOR
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 !     write(6,'(a8,3f15.10)')'SLbfr=',SL(1),SL(2),SL(3) 
00740 !     dlt_BZ=1.d-6
00741       SHIFT_b(:)=0           
00742       IF(SL(1)<=-0.5D0+dlt_BZ)THEN!20170210 
00743        SL(1)=SL(1)+1.0D0 
00744        SHIFT_b(1)=-1                
00745       ENDIF           
00746       IF(SL(2)<=-0.5D0+dlt_BZ)THEN!20170210
00747        SL(2)=SL(2)+1.0D0 
00748        SHIFT_b(2)=-1                
00749       ENDIF           
00750       IF(SL(3)<=-0.5D0+dlt_BZ)THEN!20170210 
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 !     write(6,'(a8,3f15.10)')'SLaft=',SL(1),SL(2),SL(3) 
00767 !IK_IB
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!jk               
00775       return
00776       end 
00777 !--
00778 !       call make_amat(nigs,NTB,NTG,NG0(ik),SK0(1,ik),KG0(1,1,ik),
00779 !     + b1(1),b2(1),b3(1),TAU_GAUSS(1,1),ALPHA_GAUSS(1),
00780 !     + loc_x(1,1),loc_y(1,1),loc_z(1,1),LGAUSS(1),MGAUSS(1),
00781 !     + N_BAND(ik),C0(1,1,ik),N_BAND_BTM(ik),NORM_GAUSS(1),VOLUME,
00782 !     + A_TMP(1,1)) 
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       !write(6,*) C0 
00815       !write(6,*) VOLUME  
00816       !write(6,*) ALPHA_GAUSS 
00817       !write(6,*) TAU_GAUSS 
00818       !write(6,*) NORM_GAUSS 
00819       !write(6,*) loc_x 
00820       !write(6,*) loc_y 
00821       !write(6,*) loc_z 
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!ig  
00841        !
00842        !write(6,*)'PF1'
00843        !write(6,*)PF1
00844        !write(6,*)'PF2'
00845        !write(6,*)PF2
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!ig 
00860        !
00861        !write(6,*)'VEC_KG'
00862        !write(6,*)VEC_KG 
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 !<Psi|s>
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!ig          
00882         ENDIF 
00883 !<Psi|py>
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!ig          
00897         ENDIF 
00898 !<Psi|pz>
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!ig          
00912         ENDIF 
00913 !<Psi|px>
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!ig          
00927         ENDIF 
00928 !<Psi|dxy>
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!ig          
00942         ENDIF 
00943 !<Psi|dyz>
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!ig          
00957         ENDIF 
00958 !<Psi|dz2>
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!ig          
00972         ENDIF 
00973 !<Psi|dzx>
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!ig          
00987         ENDIF 
00988 !<Psi|dx2>
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!ig          
01003         ENDIF 
01004         !--
01005         A_TMP(i_band,j_gauss)=SUM_CMPX/DSQRT(VOLUME)             
01006        enddo!i_band         
01007       enddo!j_gauss 
01008       ! 
01009       !do i_band=1,N_BAND
01010       ! write(6,'(100f12.5)')(A_TMP(i_band,j_gauss),j_gauss=1,nigs)
01011       !enddo 
01012       !write(6,*) 
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!ig 
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!ir 
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          !(ir1,ir2,ir3)->(jr1,jr2,jr3)
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          !CALC wannier(ir1,ir2,ir3)
01095          pw(ir1,ir2,ir3)=pw(ir1,ir2,ir3)+u_3D(jr1,jr2,jr3)*pf 
01096         enddo!ir3         
01097        enddo!ir2 
01098       enddo!ir1 
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!ig 
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!ir 
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        !write(6,*)SK0(:,iik) 
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!iik
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       !20180817
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       !20180820
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!ir3         
01236        enddo!ir2 
01237       enddo!ir1 
01238       !write(6,*)'rmax=',rmax 
01239       pw=pw*glb_phase 
01240       !
01241       return
01242       end

Generated on 17 Nov 2020 for respack by  doxygen 1.6.1