Skip to content

Commit

Permalink
add nk
Browse files Browse the repository at this point in the history
  • Loading branch information
@xh125
xh125 committed Sep 28, 2021
1 parent b2f9b5f commit 39bc85e
Show file tree
Hide file tree
Showing 4 changed files with 137 additions and 30 deletions.
69 changes: 64 additions & 5 deletions src_complex/hamiltonian.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ module hamiltonian
use io, only : msg,io_error,stdout
use types
use parameters
use elph2, only : nbndfst,epmatq,wf
use elph2, only : nbndfst,epmatq,wf,wf_sub,nqtotf
use modes, only : nmodes
use readepw,only : E_nk
use constants,only : ryd2mev,cone,czero
Expand All @@ -24,6 +24,10 @@ module hamiltonian
complex(kind=dp),allocatable :: gmnvkq_e(:,:,:,:,:)
complex(kind=dp),allocatable :: gmnvkq_h(:,:,:,:,:)

type(gmnvkq_n0),allocatable :: gijk_e(:),gijk_h(:)

integer :: gnzfree_e,gnzfree_h

integer :: ngfre_e,ngfre_h

!平衡位置哈密顿量，电声耦合项，总的H和临时H
Expand Down Expand Up @@ -131,18 +135,21 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihb

end subroutine allocate_hamiltonian

subroutine set_H0_nk(nk_sub,indexk,nq_sub,indexq,ibandmin,ibandmax,Enk,H0_nk,gmnvkq_eh)
subroutine set_H0_nk(nk_sub,indexk,nq_sub,indexq,ibandmin,ibandmax,Enk,H0_nk,gmnvkq_eh,lit_g,gnzfree)
use elph2, only : etf_sub,gmnvkq
use epwcom,only : kqmap_sub
implicit none
integer , intent(in) :: nk_sub,nq_sub
integer , intent(in) :: indexk(nk_sub),indexq(nq_sub),ibandmin,ibandmax
real(kind=dp), intent(out) :: Enk(ibandmax-ibandmin+1,nk_sub)
complex(kind=dp), intent(out) :: H0_nk(ibandmax-ibandmin+1,nk_sub,ibandmax-ibandmin+1,nk_sub)
complex(kind=dp), intent(out) :: gmnvkq_eh(ibandmax-ibandmin+1,ibandmax-ibandmin+1,nmodes,nk_sub,nq_sub)

integer , intent(out) :: gnzfree
real(kind=dp) , intent(in) :: lit_g

integer :: nband
integer :: ik,iband
integer :: m,n,nu,iq
integer :: ik,iband,jband
integer :: m,n,nu,iq,ikq

nband = ibandmax-ibandmin+1

Expand All @@ -160,9 +167,61 @@ subroutine set_H0_nk(nk_sub,indexk,nq_sub,indexq,ibandmin,ibandmax,Enk,H0_nk,gmn
enddo
enddo

gnzfree = 0
do iq=1,nq_sub
do ik=1,nk_sub
ikq = kqmap_sub(ik,iq)
if(ikq /= 0) then
do nu=1,nmodes
do iband=1,nband
do jband=1,nband
if(ABS(gmnvkq_eh(jband,iband,nu,ik,iq)) > sqrt(2.0*wf_sub(nu,iq)/nqtotf)*lit_g) gnzfree = gnzfree +1
enddo
enddo
enddo
endif
enddo
enddo


end subroutine set_H0_nk

subroutine set_gijk(nband,nmodes,nk_sub,nq_sub,gmnvkq,lit_g,gnzfree,gijk)
use epwcom,only : kqmap_sub
implicit none
integer, intent(in) :: nband,nmodes,nk_sub,nq_sub,gnzfree
complex(kind=dpc),intent(in) :: gmnvkq(nband,nband,nmodes,nk_sub,nq_sub)
real(kind=dp) , intent(in) :: lit_g
type(gmnvkq_n0),intent(out) :: gijk(gnzfree)

integer :: iq,ik,ikq,nu,iband,jband,ink,iqv,imkq
integer :: ig

ig = 0
do iq=1,nq_sub
do ik=1,nk_sub
ikq = kqmap_sub(ik,iq)
if(ikq /= 0) then
do nu=1,nmodes
do iband=1,nband
do jband=1,nband
if(ABS(gmnvkq(jband,iband,nu,ik,iq)) > sqrt(2.0*wf_sub(nu,iq)/nqtotf)*lit_g) then
ig = ig +1
gijk(ig)% iqv = (iq-1)*nmodes + nu
gijk(ig)% ink = (ik-1)*nband + iband
gijk(ig)% imkq= (ikq-1)*nband + jband
gijk(ig)% g = gmnvkq(jband,iband,nu,ik,iq)
endif
enddo
enddo
enddo
endif
enddo
enddo



end subroutine

!ref: 1 G. GRIMvall, <The electron-phonon interaction in metals by Goran Grimvall (z-lib.org).pdf> 1981),
! (3.20) (6.4)
Expand Down
31 changes: 19 additions & 12 deletions src_complex/lvcsh.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,10 +39,12 @@ program lvcsh
init_hband,init_e_en,init_h_en,mix_thr,lsortpes, &
calculation,verbosity,naver_sum,savedsnap,ncore
use hamiltonian,only : nefre,neband,H_e,H_e_nk,E_e,P_e,P_e_nk,P0_e_nk, &
gmnvkq_e,Enk_e,H0_e_nk,E0_e,P0_e,nhfre,nhband,H_h, &
H_h_nk,E_h,P_h,P_h_nk,P0_h_nk,gmnvkq_h,Enk_h, &
H0_h_nk,E0_h,P0_h,allocate_hamiltonian,set_H_nk, &
set_H0_nk,calculate_eigen_energy_state
gmnvkq_e,gnzfree_e,gijk_e,Enk_e,H0_e_nk,E0_e,P0_e, &
nhfre,nhband,H_h,H_h_nk,E_h,P_h,P_h_nk,P0_h_nk, &
gmnvkq_h,gnzfree_h,gijk_h,Enk_h,H0_h_nk,E0_h,P0_h, &
allocate_hamiltonian,set_H_nk, set_gijk, &
set_H0_nk,calculate_eigen_energy_state

use sortting,only : resort_eigen_energy_stat
use randoms,only : init_random_seed
use lasercom,only : fwhm,w_laser
Expand Down Expand Up @@ -74,7 +76,8 @@ program lvcsh
dc3_e,dc4_e,&
w_h,w0_h,g_h,g1_h,hsurface_type,cc0_h,dc1_h,dc2_h, &
dc3_h,dc4_h,allocatesh,convert_diabatic_adiabatic, &
calculate_nonadiabatic_coupling, &
calculate_nonadiabatic_coupling, &
calculate_nonadiabatic_coupling_,&
calculate_hopping_probability, &
convert_adiabatic_diabatic,add_decoherence
use elph2,only : wf_sub,xkf,nqtotf,nktotf,nbndfst,gmnvkq,etf,epmatq
Expand Down Expand Up @@ -126,16 +129,19 @@ program lvcsh
call allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihband_min,ihband_max)

if(lelecsh) then
call set_H0_nk(nk_sub,indexk,nq_sub,indexq,ieband_min,ieband_max,Enk_e,H0_e_nk,gmnvkq_e)
call set_H0_nk(nk_sub,indexk,nq_sub,indexq,ieband_min,ieband_max,Enk_e,H0_e_nk,gmnvkq_e,lit_gmnvkq,gnzfree_e)
!set Enk_e(nband,nk_sub) H0_e_nk and gmnvkq_e
allocate(gijk_e(gnzfree_e))
call set_gijk(neband,nmodes,nk_sub,nq_sub,gmnvkq_e,lit_gmnvkq,gnzfree_e,gijk_e)
if(nefre_sh == 0 .or. nefre_sh > nefre) nefre_sh = nefre
endif

if(lholesh) then
call set_H0_nk(nk_sub,indexk,nq_sub,indexq,ihband_min,ihband_max,Enk_h,H0_h_nk,gmnvkq_h)
call set_H0_nk(nk_sub,indexk,nq_sub,indexq,ihband_min,ihband_max,Enk_h,H0_h_nk,gmnvkq_h,lit_gmnvkq,gnzfree_h)
H0_h_nk = -1.0 * H0_h_nk
gmnvkq_h = -1.0 * gmnvkq_h
!set Enk_h(nband,nk_sub) H0_h_nk and gmnvkq_h(:::,nk_sub,nq_sub)
!set Enk_h(nband,nk_sub) H0_h_nk and gmnvkq_h(:::,nk_sub,nq_sub)
allocate(gijk_h(gnzfree_h))
if(nhfre_sh == 0 .or. nhfre_sh > nhfre) nhfre_sh = nhfre
endif

Expand Down Expand Up @@ -191,7 +197,8 @@ program lvcsh
call convert_adiabatic_diabatic(nefre,P_e,w_e,c_e)
c_e_nk = reshape(c_e,(/ neband,nk_sub /))

call calculate_nonadiabatic_coupling(nmodes,nq_sub,neband,nk_sub,E_e,P_e_nk,gmnvkq_e,lit_gmnvkq,nefre_sh,d_e)
call calculate_nonadiabatic_coupling(nmodes,nq_sub,neband,nk_sub,E_e,P_e_nk,gmnvkq_e,nefre_sh,d_e)
!call calculate_nonadiabatic_coupling_(nmodes,nq_sub,neband,nk_sub,E_e,P_e_nk,gnzfree_e,gijk_e,nefre_sh,d_e)
E0_e = E_e;P0_e=P_e;P0_e_nk=P_e_nk;d0_e=d_e;w0_e=w_e
endif

Expand All @@ -207,7 +214,7 @@ program lvcsh
call convert_adiabatic_diabatic(nefre,P_h,w_h,c_h)
c_h_nk = reshape(c_h,(/ nhband,nk_sub /))

call calculate_nonadiabatic_coupling(nmodes,nq_sub,nhband,nk_sub,E_h,P_h_nk,gmnvkq_h,lit_gmnvkq,nhfre_sh,d_h)
call calculate_nonadiabatic_coupling(nmodes,nq_sub,nhband,nk_sub,E_h,P_h_nk,gmnvkq_h,nhfre_sh,d_h)
E0_h = E_h;P0_h=P_h;P0_h_nk=P_h_nk;d0_h=d_h;w0_h=w_h
endif

Expand Down Expand Up @@ -274,7 +281,7 @@ program lvcsh

! Calculate non-adiabatic coupling vectors with the Hellmann-Feynman theorem.
! update d_e in time t0+dt
call calculate_nonadiabatic_coupling(nmodes,nq_sub,neband,nk_sub,E_e,p_e,gmnvkq_e,lit_gmnvkq,nefre_sh,d_e)
call calculate_nonadiabatic_coupling(nmodes,nq_sub,neband,nk_sub,E_e,p_e,gmnvkq_e,nefre_sh,d_e)


! use p_e in time t0+dt, to convert c_e(t0+dt) to w_e(t0+dt)
Expand Down Expand Up @@ -334,7 +341,7 @@ program lvcsh

! Calculate non-adiabatic coupling vectors with the Hellmann-Feynman theorem.
! update d_h in time t0+dt
call calculate_nonadiabatic_coupling(nmodes,nq_sub,nhband,nk_sub,E_h,p_h,gmnvkq_h,lit_gmnvkq,nhfre_sh,d_h)
call calculate_nonadiabatic_coupling(nmodes,nq_sub,nhband,nk_sub,E_h,p_h,gmnvkq_h,nhfre_sh,d_h)


! use p_h in time t0+dt, to convert c_h(t0+dt) to w_h(t0+dt)
Expand Down
54 changes: 49 additions & 5 deletions src_complex/surfacehopping.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -349,19 +349,18 @@ end subroutine add_decoherence
!==================================================================!
! ref : PPT-91
! The most time-consuming part of the program
subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,lit_gmnvkq,nfre_sh,dd)
subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfre_sh,dd)
use kinds,only : dp
implicit none
integer, intent(in) :: nmodes,nq,nband,nk,nfre_sh
real(kind=dp),intent(in) :: ee(nband*nk)
complex(kind=dpc),intent(in) :: p_nk(nband,nk,nband*nk)
real(kind=dp),intent(in) :: lit_gmnvkq
complex(kind=dpc),intent(in) :: gmnvkq(nband,nband,nmodes,nk,nq)
complex(kind=dpc),intent(out):: dd(nfre_sh,nfre_sh,nmodes,nq)

integer :: nfre,ifre,jfre,iq,imode ,igfre
integer :: ik,ikq,iband1,iband2
complex(kind=dpc) :: epc
complex(kind=dpc) :: epc,pnk12

!nfre = nband*nk
nfre = nfre_sh
Expand All @@ -379,7 +378,7 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,lit
do ifre=1,nfre
do jfre=1,nfre
dd(ifre,jfre,imode,iq) = dd(ifre,jfre,imode,iq)+&
&epc*CONJG(p_nk(iband1,ik,ifre))*p_nk(iband2,ikq,jfre)
epc*CONJG(p_nk(iband1,ik,ifre))*p_nk(iband2,ikq,jfre)
enddo
enddo
endif
Expand All @@ -389,7 +388,7 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,lit
endif
enddo
enddo

do ifre=1,nfre
do jfre=1,nfre
if(jfre/= ifre) then
Expand All @@ -402,6 +401,51 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,lit

end subroutine calculate_nonadiabatic_coupling

subroutine calculate_nonadiabatic_coupling_(nmodes,nq,nband,nk,ee,p_nk,gnzfree,gijk,nfre_sh,dd)
use types
use kinds,only : dp
implicit none
integer, intent(in) :: nmodes,nq,nband,nk,nfre_sh,gnzfree
real(kind=dp),intent(in) :: ee(nband*nk)
complex(kind=dpc),intent(in) :: p_nk(nband*nk,nband*nk)
type(gmnvkq_n0),intent(in) :: gijk(gnzfree)
complex(kind=dpc),intent(out):: dd(nfre_sh,nfre_sh,nmodes*nq)

integer :: nfre,ifre,jfre,iq,imode ,ig
integer :: ik,ikq,iband1,iband2 ,iqv,ink,imkq
complex(kind=dpc) :: epc,pnk12,g

nfre = nfre_sh

dd=czero
do ig=1,gnzfree
iqv = gijk(ig) % iqv
ink = gijk(ig) % ink
imkq= gijk(ig) % imkq
epc = gijk(ig) % g
do ifre=1,nfre
do jfre=1,nfre
dd(ifre,jfre,iqv) = dd(ifre,jfre,iqv)+&
epc*CONJG(p_nk(ink,ifre))*p_nk(imkq,jfre)
enddo
enddo

enddo



do ifre=1,nfre
do jfre=1,nfre
if(jfre/= ifre) then
dd(ifre,jfre,:) = dd(ifre,jfre,:)/(ee(jfre)-ee(ifre))
else
!dd(ifre,jfre,:,:) = czero
endif
enddo
enddo

end subroutine calculate_nonadiabatic_coupling_


function bolziman(womiga,temp)
use kinds ,only : dp
Expand Down
13 changes: 5 additions & 8 deletions src_complex/types.f90
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,12 @@
module types
use kinds,only :dp
use kinds,only :dp,dpc
implicit none
!declare type of gmnvkq_n0 with gmnvkq>0.0
type :: gmnvkq_n0
integer :: m
integer :: n
integer :: v
integer :: ik
integer :: iq
integer :: ikq
real(kind=dp) :: g
integer :: iqv
integer :: ink
integer :: imkq
complex(kind=dpc) :: g
end type gmnvkq_n0

end module types

0 comments on commit 39bc85e

Please sign in to comment.