Skip to content

Commit

Permalink
add resort poential energy surface subroutine
Browse files Browse the repository at this point in the history
  • Loading branch information
@xh125
xh125 committed Jun 25, 2021
1 parent c6cfa9a commit aec5669
Show file tree
Hide file tree
Showing 4 changed files with 92 additions and 20 deletions.
65 changes: 61 additions & 4 deletions src/hamiltonian.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,9 @@ module hamiltonian
real(kind=dp),allocatable :: H0_e_nk(:,:,:,:),H_e_nk(:,:,:,:)
real(kind=dp),allocatable :: H0_h(:,:),H_h(:,:)
real(kind=dp),allocatable :: H0_h_nk(:,:,:,:),H_h_nk(:,:,:,:)

real(kind=dp),allocatable :: H_e_eq(:,:),E_e_eq(:),P_e_eq(:,:) !for position of equilibrium
real(kind=dp),allocatable :: H_h_eq(:,:),E_h_eq(:),P_h_eq(:,:) !for position of equilibrium

real(kind=dp),allocatable :: gmnvkq_e(:,:,:,:,:)
real(kind=dp),allocatable :: gmnvkq_h(:,:,:,:,:)
Expand All @@ -25,7 +28,6 @@ module hamiltonian
E0_e(:),P0_e(:,:),P0_e_nk(:,:,:)
real(kind=dp),allocatable :: E_h(:),P_h(:,:),P_h_nk(:,:,:),&
E0_h(:),P0_h(:,:),P0_h_nk(:,:,:)
real(kind=dp),allocatable :: E_h_(:)
!哈密顿量的本征值与本征矢
integer :: ierr
contains
Expand All @@ -52,7 +54,16 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min
allocate(H_e_nk(neband,nk,neband,nk),stat=ierr)
if(ierr /=0) call errore('hamiltonian','Error allocating H_e_nk',1)
H_e_nk = 0.0

allocate(E_e_eq(nefre),stat=ierr)
if(ierr /=0) call errore('hamiltonian','Error allocating E_e_eq',1)
E_e_eq = 0.0
allocate(P_e_eq(nefre,nefre),stat=ierr)
if(ierr /=0) call errore('hamiltonian','Error allocating P_e_eq',1)
P_e_eq = 0.0
allocate(H_e_eq(nefre,nefre),stat=ierr)
if(ierr /=0) call errore('hamiltonian','Error allocating H_e_eq',1)
H_e_eq = 0.0

endif

if(lholesh) then
Expand All @@ -72,7 +83,17 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min
H_h = 0.0
allocate(H_h_nk(nhband,nk,nhband,nk),stat=ierr)
if(ierr /=0) call errore('hamiltonian','Error allocating H_h_nk',1)
H_h_nk = 0.0
H_h_nk = 0.0
allocate(E_h_eq(nhfre),stat=ierr)
if(ierr /=0) call errore('Energy of hamiltonian for position of equilibrium ','Error allocating E_h_eq',1)
E_h_eq = 0.0
allocate(P_h_eq(nhfre,nhfre),stat=ierr)
if(ierr /=0) call errore('hamiltonian','Error allocating P_h_eq',1)
P_h_eq = 0.0
allocate(H_h_eq(nhfre,nhfre),stat=ierr)
if(ierr /=0) call errore('hamiltonian for position of equilibrium ','Error allocating E_h_eq',1)
H_h_eq = 0.0

endif

end subroutine allocate_hamiltonian
Expand Down Expand Up @@ -129,7 +150,7 @@ subroutine set_H_nk(nband,nk,nmodes,nq,ph_Q,wf,gmnvkq_eh,H0_nk,H_nk)
do iband=1,nband !|iband,ik>
do jband=1,nband !|jband,ikq>
H_nk(jband,ikq,iband,ik) = H_nk(jband,ikq,iband,ik)+&
ph_Q(nu,iq)*sqrt(2.0*wf(nu,iq)/nq)*gmnvkq_eh(iband,jband,nu,ik,iq)
gmnvkq_eh(iband,jband,nu,ik,iq)*ph_Q(nu,iq)*sqrt(2.0*wf(nu,iq)/nq)
enddo
enddo
enddo
Expand Down Expand Up @@ -158,4 +179,40 @@ subroutine calculate_eigen_energy_state(nfre,H,ee,pp)

end subroutine calculate_eigen_energy_state

subroutine resort_eigen_energy_stat(nfre,ee,pp,ee_eq,pp_eq)
implicit none
integer,intent(in) :: nfre
real(kind=dp),intent(inout) :: ee(nfre),pp(nfre,nfre)
real(kind=dp),intent(in) :: ee_eq(nfre),pp_eq(nfre,nfre)

real(kind=dp),allocatable :: p_tmp(:),pdotp(:)
real(kind=dp) :: e_tmp,flad

integer :: ifre,jfre,cfre(1)

if(.not. allocated(p_tmp)) allocate(p_tmp(nfre))
if(.not. allocated(pdotp)) allocate(pdotp(nfre))


do ifre =1 ,nfre
flad = ABS(SUM(pp(:,ifre)*pp_eq(:,ifre)))
if(flad >= 0.5) then
cycle
else
pdotp = 0.0
do jfre=ifre,nfre
pdotp(jfre) = ABS(SUM(pp(:,jfre)*pp_eq(:,ifre)))
enddo
cfre = Maxloc(pdotp)
e_tmp= ee(ifre)
p_tmp= pp(:,ifre)
ee(ifre) = ee(cfre(1))
pp(:,ifre) = pp(:,cfre(1))
ee(cfre(1))= e_tmp
pp(:,cfre(1)) = p_tmp
endif
enddo

end subroutine

end module hamiltonian
17 changes: 13 additions & 4 deletions src/lvcsh.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,11 +15,11 @@ program lvcsh
!=====================================================================================================!
!= system interaction with environment by system-bath interactions =!
!=====================================================================================================!
!= Last updata 2021-05.18 Version:0.1.3 =!
!= Last updata 2021-06.25 Version:0.1.5 =!
!= Developed by XieHua at department of physic, USTC;[email protected] =!
!=====================================================================================================!
!! Author: HuaXie
!! Version: v0.1.1
!! Version: v0.1.5
!! License: GNU
!!=================================================================================================
use mkl_service
Expand All @@ -38,9 +38,9 @@ program lvcsh
calculation,verbosity,nnode,ncore,naver_sum
use hamiltonian,only : nefre,neband,H_e,H_e_nk,E_e,P_e,P_e_nk,P0_e_nk,gmnvkq_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,H0_h_nk,E0_h,P0_h,&
E_h_,&
H_e_eq,E_e_eq,P_e_eq,H_h_eq,E_h_eq,P_h_eq,&
allocate_hamiltonian,set_H_nk,set_H0_nk,&
calculate_eigen_energy_state
calculate_eigen_energy_state,resort_eigen_energy_stat
use randoms,only : init_random_seed
use lasercom,only : fwhm,w_laser
use getwcvk,only : get_Wcvk
Expand Down Expand Up @@ -112,10 +112,15 @@ program lvcsh
call allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min,ihband_max)
if(lelecsh) then
call set_H0_nk(nktotf,neband,H0_e_nk,ieband_min,gmnvkq_e)
H_e_eq = reshape(H0_e_nk,(/ nefre,nefre /))
call calculate_eigen_energy_state(nefre,H_e_eq,E_e_eq,P_e_eq)
endif
if(lholesh) then
call set_H0_nk(nktotf,nhband,H0_h_nk,ihband_min,gmnvkq_h)
H0_h_nk = -1.0 * H0_h_nk
H_h_eq = reshape(H0_h_nk,(/ nhfre,nhfre /))
call calculate_eigen_energy_state(nhfre,H_h_eq,E_h_eq,P_h_eq)

gmnvkq_h = -1.0 * gmnvkq_h
endif
!get H0_e_nk(neband,nktotf,neband,nktotf),gmnvkq_e(neband,neband,nktotf,nmodes,nqtotf)
Expand Down Expand Up @@ -254,6 +259,7 @@ program lvcsh
call set_H_nk(neband,nktotf,nmodes,nqtotf,phQ,wf,gmnvkq_e,H0_e_nk,H_e_nk)
H_e = reshape(H_e_nk,(/ nefre,nefre /))
call calculate_eigen_energy_state(nefre,H_e,E_e,P_e)
call resort_eigen_energy_stat(nefre,E_e,P_e,E_e_eq,P_e_eq)
P_e_nk = reshape(P_e,(/ neband,nktotf,nefre /))
call convert_diabatic_adiabatic(nefre,P_e,c_e,w_e)
call init_surface(nefre,w_e,iesurface)
Expand All @@ -267,6 +273,7 @@ program lvcsh
call set_H_nk(nhband,nktotf,nmodes,nqtotf,phQ,wf,gmnvkq_h,H0_h_nk,H_h_nk)
H_h = reshape(H_h_nk,(/ nhfre,nhfre /))
call calculate_eigen_energy_state(nhfre,H_h,E_h,P_h)
call resort_eigen_energy_stat(nhfre,E_h,P_h,E_h_eq,P_h_eq)
P_h_nk = reshape(P_h,(/ nhband,nktotf,nhfre /))
call convert_diabatic_adiabatic(nhfre,P_h,c_h,w_h)
call init_surface(nhfre,w_h,ihsurface)
Expand Down Expand Up @@ -454,6 +461,7 @@ program lvcsh
H_e = reshape(H_e_nk,(/ nefre,nefre /))
! update E_e,P_e to time t0+dt
call calculate_eigen_energy_state(nefre,H_e,E_e,P_e)
call resort_eigen_energy_stat(nefre,E_e,P_e,E_e_eq,P_e_eq)
P_e_nk = reshape(P_e,(/ neband,nktotf,nefre /))

! Calculate non-adiabatic coupling vectors with the Hellmann-Feynman theorem.
Expand Down Expand Up @@ -511,6 +519,7 @@ program lvcsh
H_h = reshape(H_h_nk,(/ nhfre,nhfre /))
! update E_h,P_h in time t0+dt
call calculate_eigen_energy_state(nhfre,H_h,E_h,P_h)
call resort_eigen_energy_stat(nhfre,E_h,P_h,E_h_eq,P_h_eq)
P_h_nk = reshape(P_h,(/ nhband,nktotf,nhfre /))

! Calculate non-adiabatic coupling vectors with the Hellmann-Feynman theorem.
Expand Down
15 changes: 8 additions & 7 deletions src/readepw.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1057,20 +1057,21 @@ subroutine readepwout(fepwout)
enddo


do iq=1,nqtotf
do nu=1,nmodes
if(wf(nu,iq)<= 0.0) then
gmnvkq(:,:,nu,:,iq) = 0.0
endif
enddo
enddo
!do iq=1,nqtotf
! do nu=1,nmodes
! if(wf(nu,iq)<= 0.0) then
! gmnvkq(:,:,nu,:,iq) = 0.0
! endif
! enddo
!enddo

do iq=1,nqtotf
do nu=1,nmodes
if(wf(nu,iq)<0.0) then
write(stdout,"(A,I5,1X,A,3(F12.6,1X),A8,I5,A1,F12.6,A3)") &
"Carefully!!! the energy of phonon in iq=",iq,"(coord.:",(xqf(ipol,iq),ipol=1,3),") modes=",nu,"=",wf(nu,iq)*ryd2mev,"meV"
write(stdout,"(A)") "Setting asr_type='simple' in epw.in could solve gamma point problem."
write(stdout,"(A)") "Take a callfully check for the ph.x calculation."
wf(nu,iq)=0.0
endif
enddo
Expand Down
15 changes: 10 additions & 5 deletions src/surfacehopping.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ module surfacehopping
use elph2,only : wf,nktotf,nbndfst,ibndmin,ibndmax
use hamiltonian,only : nphfre,neband,nhband,nefre,nhfre,&
E_e,P_e,P_e_nk,E0_e,P0_e,P0_e_nk,&
E_h,E_h_,P_h,P_h_nk,E0_h,P0_h,P0_h_nk
E_h,P_h,P_h_nk,E0_h,P0_h,P0_h_nk
use parameters, only : nsnap,naver,ncore,nnode
use surfacecom, only : iesurface,ihsurface,esurface_type,hsurface_type,&
phQ,phP,phQ0,phP0,phK,phU,SUM_phU,SUM_phK,SUM_phK0,SUM_phE,&
Expand Down Expand Up @@ -73,6 +73,8 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq)
if(ierr /=0) call errore('surfacehopping','Error allocating P_e',1)
allocate(P_e_nk(neband,nktotf,nefre),stat=ierr)
if(ierr /=0) call errore('surfacehopping','Error allocating P_e_nk',1)
allocate(P0_e_nk(neband,nktotf,nefre),stat=ierr)
if(ierr /=0) call errore('surfacehopping','Error allocating P0_e_nk',1)
allocate(d_e(nefre,nefre,nmodes,nq),stat=ierr) !d_ijk
if(ierr /=0) call errore('surfacehopping','Error allocating d_e',1)
allocate(dEa_dQ_e(nmodes,nq))
Expand Down Expand Up @@ -142,6 +144,8 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq)
if(ierr /=0) call errore('surfacehopping','Error allocating P_h',1)
allocate(P_h_nk(nhband,nktotf,nhfre),stat=ierr)
if(ierr /=0) call errore('surfacehopping','Error allocating P_h_nk',1)
allocate(P0_h_nk(nhband,nktotf,nhfre),stat=ierr)
if(ierr /=0) call errore('surfacehopping','Error allocating P0_h_nk',1)
allocate(d_h(nhfre,nhfre,nmodes,nq),stat=ierr) !d_ijk
if(ierr /=0) call errore('surfacehopping','Error allocating d_h',1)
allocate(dEa_dQ_h(nmodes,nq))
Expand Down Expand Up @@ -190,8 +194,8 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq)

allocate(E0_h(1:nhfre),stat=ierr)
if(ierr /=0) call errore('surfacehopping','Error allocating E0_h',1)
allocate(E_h_(1:nhfre),stat=ierr)
if(ierr /=0) call errore('surfacehopping','Error allocating E_h_',1)
!allocate(E_h_(1:nhfre),stat=ierr)
!if(ierr /=0) call errore('surfacehopping','Error allocating E_h_',1)

allocate(P0_h(nhfre,nhfre),stat=ierr)
if(ierr /=0) call errore('surfacehopping','Error allocating P0_h',1)
Expand Down Expand Up @@ -328,19 +332,20 @@ subroutine add_decoherence(C,Ekin,dt,nfre,isurface,E,ww)

integer :: ifre
real(kind=dp) :: tau
real(kind=dp) :: flad
real(kind=dp) :: flad,factor

tau = 0.0
flad = 0.0
do ifre=1,nfre
if(ifre /= isurface) then
tau = (1.0 + C/Ekin)/abs(E(ifre)-E(isurface))
factor = exp(-1.0*dt/tau)
ww(ifre) = ww(ifre)*exp(-1.0*dt/tau)
flad = flad + ww(ifre)*CONJG(ww(ifre))
endif
enddo

ww(isurface) = ww(isurface)*sqrt(1.0 - flad)/abs(ww(isurface))
ww(isurface) = ww(isurface)*sqrt(1.0 - flad)/sqrt(ww(isurface)*CONJG(ww(isurface)))

end subroutine add_decoherence

Expand Down

0 comments on commit aec5669

Please sign in to comment.