From 25d0d6ca27714583c2cfc7b48b74667b39a5b60d Mon Sep 17 00:00:00 2001 From: xh125 Date: Sat, 9 Oct 2021 16:10:33 +0800 Subject: [PATCH] change dEa2_dQ2 --- src_complex/cc_fssh.f90 | 10 ++++---- src_complex/dynamics.f90 | 47 ++-------------------------------- src_complex/fssh.f90 | 8 +++--- src_complex/lvcsh.f90 | 40 ++++++++++++----------------- src_complex/sc_fssh.f90 | 8 +++--- src_complex/surfacecom.f90 | 2 +- src_complex/surfacehopping.f90 | 26 +++++++++++-------- 7 files changed, 48 insertions(+), 93 deletions(-) diff --git a/src_complex/cc_fssh.f90 b/src_complex/cc_fssh.f90 index cb7b77d..a8955e4 100644 --- a/src_complex/cc_fssh.f90 +++ b/src_complex/cc_fssh.f90 @@ -87,7 +87,7 @@ subroutine nonadiabatic_transition_ccfssh(lfeedback,nfre,nfre_sh,nmodes,nq,& integer , intent(out) :: surface_type real(kind=dp),intent(in) :: EE(nfre) complex(kind=dpc),intent(in) :: P(nfre,nfre),P0(nfre,nfre) - complex(kind=dpc),intent(in) :: DD(nfre_sh,nfre_sh,nmodes,nq) + complex(kind=dpc),intent(in) :: DD(2,nfre_sh,nmodes,nq) real(kind=dp),intent(out):: S_bi(nfre_sh) real(kind=dp),intent(in) :: GG(nfre_sh) complex(kind=dpc),intent(inout) :: VV(nmodes,nq) @@ -138,8 +138,8 @@ subroutine nonadiabatic_transition_ccfssh(lfeedback,nfre,nfre_sh,nmodes,nq,& sumdd = 0.0 do iq=1,nq do imode=1,nmodes - sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(isurface,ifre,imode,iq)) ! A - sumdd = sumdd + ABS(DD(isurface,ifre,imode,iq))**2 ! B + sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(1,ifre,imode,iq)) ! A + sumdd = sumdd + ABS(DD(1,ifre,imode,iq))**2 ! B enddo enddo detaE = EE(isurface_a)-EE(ifre) ! C @@ -152,7 +152,7 @@ subroutine nonadiabatic_transition_ccfssh(lfeedback,nfre,nfre_sh,nmodes,nq,& flagd = sumvd/sumdd*(-1.0+dsqrt(flagd)) do iq=1,nq do imode=1,nmodes - if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(isurface,ifre,imode,iq) + if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(1,ifre,imode,iq) enddo enddo isurface = isurface_k @@ -176,7 +176,7 @@ subroutine nonadiabatic_transition_ccfssh(lfeedback,nfre,nfre_sh,nmodes,nq,& flagd = sumvd/sumdd*(-1.0+dsqrt(flagd)) do iq=1,nq do imode=1,nmodes - if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(isurface,ifre,imode,iq) + if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(1,ifre,imode,iq) enddo enddo isurface = isurface_k diff --git a/src_complex/dynamics.f90 b/src_complex/dynamics.f90 index 132b832..558c9a1 100644 --- a/src_complex/dynamics.f90 +++ b/src_complex/dynamics.f90 @@ -27,40 +27,6 @@ subroutine set_gamma(nmodes,nq,gamma,ld_fric,wf,ld_gamma) end subroutine - subroutine get_dEa_dQ(nmodes,nq,nband,nk,P_nk,gmnvkq,isurface,dEa_dQ) - use elph2,only : iminusq - implicit none - integer,intent(in) :: nmodes,nq,nband,nk - integer,intent(in) :: isurface - complex(kind=dpc),intent(in) :: P_nk(nband,nk,nband*nk) - complex(kind=dpc),intent(in) :: gmnvkq(nband,nband,nmodes,nk,nq) - complex(kind=dpc),intent(out) :: dEa_dQ(nmodes,nq) - - integer :: iq,imode,ik,ikq,iband1,iband2 - complex(kind=dpc) :: epc - complex(kind=dpc) :: cmp_tmp1,cmp_tmp2 - logical :: ltmp = .true. - - ! F_qv = -dEa_dQ_-qv (2.60) - dEa_dQ = czero - do iq=1,nq - do imode=1,nmodes - do ik=1,nk - ikq = kqmap(ik,iq) - do iband1=1,nband ! m - do iband2=1,nband ! n - epc = gmnvkq(iband1,iband2,imode,ik,iq) - if(epc /= czero) then - dEa_dQ(imode,iminusq(iq)) = dEa_dQ(imode,iminusq(iq)) + & - epc*CONJG(P_nk(iband1,ikq,isurface))*P_nk(iband2,ik,isurface) - endif - enddo - enddo - enddo - enddo - enddo - - end subroutine get_dEa_dQ !=======================================================================! != rk4 method to obtain coordinate and velocitie after a time interval =! @@ -135,14 +101,6 @@ subroutine derivs_nuclei(nmodes,nq,dEa_dQ,wf,ld_gamma,xx,vv,dx,dv) dv = -wf**2*xx - dEa_dQ - ld_gamma*vv !(2.60) PPT-94 - ! F_qv = -dEa_dQ_-qv - !do iq=1,nq - ! do imode=1,nmodes - ! dx(imode,iq) = vv(imode,iq) - ! dv(imode,iq) = -wf(imode,iq)**2*xx(imode,iq)-dEa_dQ(imode,iq)-ld_gamma(imode,iq)*vv(imode,iq) - ! enddo - !enddo - endsubroutine derivs_nuclei @@ -219,12 +177,11 @@ subroutine get_dEa2_dQ2(nmodes,nq,nfre,nfre_sh,isurface,EE,dd,dEa2_dQ2) integer,intent(in) :: nmodes,nq,nfre,nfre_sh integer,intent(in) :: isurface real(kind=dp),intent(in) :: EE(nfre) - complex(kind=dpc),intent(in) :: dd(nfre_sh,nmodes,nq) + complex(kind=dpc),intent(in) :: dd(2,nfre_sh,nmodes,nq) real(kind=dp),intent(out) :: dEa2_dQ2(nmodes,nq) integer :: iq,imode,ifre - ! dEa2_dQ2 = dEa2/d(Q_qv*)d(Q_qv) ! ref : (2.60) PPT-95 dEa2_dQ2 = 0.0 do iq=1,nq @@ -232,7 +189,7 @@ subroutine get_dEa2_dQ2(nmodes,nq,nfre,nfre_sh,isurface,EE,dd,dEa2_dQ2) do ifre=1,nfre_sh if(ifre /= isurface) then dEa2_dQ2(imode,iq) = dEa2_dQ2(imode,iq) + & - (EE(isurface)-EE(ifre))*(DD(ifre,imode,iq)*CONJG(DD(ifre,imode,iq)))*2.0 + (EE(isurface)-EE(ifre))*REAL((DD(1,ifre,imode,iq)*DD(2,ifre,imode,iq)))*2.0 endif enddo enddo diff --git a/src_complex/fssh.f90 b/src_complex/fssh.f90 index a8cdab5..ef7919b 100644 --- a/src_complex/fssh.f90 +++ b/src_complex/fssh.f90 @@ -13,7 +13,7 @@ subroutine nonadiabatic_transition_fssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurfac integer , intent(inout) :: isurface real(kind=dp),intent(in) :: EE(nfre) complex(kind=dpc),intent(in) :: P(nfre,nfre) - complex(kind=dpc),intent(in) :: DD(nfre_sh,nmodes,nq) + complex(kind=dpc),intent(in) :: DD(2,nfre_sh,nmodes,nq) real(kind=dp),intent(in) :: GG(nfre_sh) complex(kind=dp),intent(inout) :: VV(nmodes,nq) @@ -37,8 +37,8 @@ subroutine nonadiabatic_transition_fssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurfac sumdd = 0.0 do iq=1,nq do imode=1,nmodes - sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(ifre,imode,iq)) ! A - sumdd = sumdd + ABS(DD(ifre,imode,iq))**2 ! B + sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(1,ifre,imode,iq)) ! A + sumdd = sumdd + ABS(DD(1,ifre,imode,iq))**2 ! B enddo enddo detaE = EE(isurface)-EE(ifre) !C @@ -51,7 +51,7 @@ subroutine nonadiabatic_transition_fssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurfac flagd = (sumvd/sumdd)*(-1.0+dsqrt(flagd)) ! effective scratting time do iq=1,nq do imode=1,nmodes - if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(ifre,imode,iq) + if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(1,ifre,imode,iq) enddo enddo isurface = ifre diff --git a/src_complex/lvcsh.f90 b/src_complex/lvcsh.f90 index 3fa3108..240da99 100644 --- a/src_complex/lvcsh.f90 +++ b/src_complex/lvcsh.f90 @@ -80,7 +80,7 @@ program lvcsh use cell_base,only : bg use date_and_times,only : get_date_and_time use io ,only : stdout,io_time,time1,time2,time_,msg,io_error - use dynamics,only : set_gamma,get_dEa_dQ,get_dEa2_dQ2,rk4_nuclei,& + use dynamics,only : set_gamma,get_dEa2_dQ2,rk4_nuclei,& rk4_electron_diabatic,ADD_BATH_EFFECT,pre_md use memory_report,only : MB,GB,complex_size,real_size,int_size,ram,& print_memory @@ -183,6 +183,7 @@ program lvcsh c_e_nk = reshape(c_e,(/ neband,nktotf /)) call calculate_nonadiabatic_coupling(nmodes,nqtotf,neband,nktotf,E_e,P_e_nk,gmnvkq_e,nefre_sh,iesurface,d_e) + dEa_dQ_e = d_e(1,iesurface,:,:) E0_e = E_e;P0_e=P_e;P0_e_nk=P_e_nk;d0_e=d_e;w0_e=w_e endif @@ -199,6 +200,7 @@ program lvcsh c_h_nk = reshape(c_h,(/ nhband,nktotf /)) call calculate_nonadiabatic_coupling(nmodes,nqtotf,nhband,nktotf,E_h,P_h_nk,gmnvkq_h,nhfre_sh,ihsurface,d_h) + dEa_dQ_h = d_h(1,ihsurface,:,:) E0_h = E_h;P0_h=P_h;P0_h_nk=P_h_nk;d0_h=d_h;w0_h=w_h endif @@ -222,22 +224,15 @@ program lvcsh !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! dEa_dQ = 0.0 !dEa_dQ in time t0 - if(l_dEa_dQ) then - if(lelecsh) then - call get_dEa_dQ(nmodes,nqtotf,neband,nktotf,P0_e_nk,gmnvkq_e,iesurface,dEa_dQ_e) - dEa_dQ = dEa_dQ + dEa_dQ_e - endif - if(lholesh) then - call get_dEa_dQ(nmodes,nqtotf,nhband,nktotf,P0_h_nk,gmnvkq_h,ihsurface,dEa_dQ_h) - dEa_dQ = dEa_dQ + dEa_dQ_h - endif - endif + if(lelecsh) dEa_dQ = dEa_dQ + dEa_dQ_e + if(lholesh) dEa_dQ = dEa_dQ + dEa_dQ_h !==============================! != update phQ,phP =! !==============================! !use rk4 to calculate the dynamical of phonon normal modes !update phQ,phP to time t0+dt + if(.not. l_dEa_dQ) dEa_dQ = 0.0 call rk4_nuclei(nmodes,nqtotf,dEa_dQ,ld_gamma,wf,phQ,phP,dt) @@ -266,7 +261,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,nqtotf,neband,nktotf,E_e,p_e,gmnvkq_e,nefre_sh,iesurface,d_e) - + dEa_dQ_e = d_e(1,iesurface,:,:) ! use p_e in time t0+dt, to convert c_e(t0+dt) to w_e(t0+dt) call convert_diabatic_adiabatic(nefre,p_e,c_e,w_e) @@ -326,7 +321,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,nqtotf,nhband,nktotf,E_h,p_h,gmnvkq_h,nhfre_sh,ihsurface,d_h) - + dEa_dQ_h = d_h(1,ihsurface,:,:) ! use p_h in time t0+dt, to convert c_h(t0+dt) to w_h(t0+dt) call convert_diabatic_adiabatic(nhfre,p_h,c_h,w_h) @@ -368,20 +363,19 @@ program lvcsh !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! dEa2_dQ2 = 0.0 !dEa2_dQ2 in time t0 - if(l_dEa2_dQ2) then - if(lelecsh) then - call get_dEa2_dQ2(nmodes,nqtotf,nefre,nefre_sh,iesurface,E0_e,d0_e,dEa2_dQ2_e) - dEa2_dQ2 = dEa2_dQ2 + dEa2_dQ2_e - endif - if(lholesh) then - call get_dEa2_dQ2(nmodes,nqtotf,nhfre,nhfre_sh,ihsurface,E0_h,d0_h,dEa2_dQ2_h) - dEa2_dQ2 = dEa2_dQ2 + dEa2_dQ2_h - endif + if(lelecsh) then + call get_dEa2_dQ2(nmodes,nqtotf,nefre,nefre_sh,iesurface,E0_e,d0_e,dEa2_dQ2_e) + dEa2_dQ2 = dEa2_dQ2 + dEa2_dQ2_e + endif + if(lholesh) then + call get_dEa2_dQ2(nmodes,nqtotf,nhfre,nhfre_sh,ihsurface,E0_h,d0_h,dEa2_dQ2_h) + dEa2_dQ2 = dEa2_dQ2 + dEa2_dQ2_h endif !===================! != add bath effect =! - !===================! + !===================! + if(.not. l_dEa2_dQ2) dEa2_dQ2 = 0.0 call add_bath_effect(nmodes,nqtotf,wf,ld_gamma,temp,dEa2_dQ2,dt,l_ph_quantum,phQ,phP) diff --git a/src_complex/sc_fssh.f90 b/src_complex/sc_fssh.f90 index 5bc0ce7..4d49fe3 100644 --- a/src_complex/sc_fssh.f90 +++ b/src_complex/sc_fssh.f90 @@ -70,7 +70,7 @@ subroutine nonadiabatic_transition_scfssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurf integer , intent(inout) :: isurface real(kind=dp),intent(in) :: EE(nfre) complex(kind=dpc),intent(in) :: P(nfre,nfre) - complex(kind=dpc),intent(in) :: DD(nfre_sh,nfre_sh,nmodes,nq) + complex(kind=dpc),intent(in) :: DD(2,nfre_sh,nmodes,nq) real(kind=dp),intent(in) :: GG(nfre_sh) complex(kind=dpc),intent(inout) :: VV(nmodes,nq) @@ -94,8 +94,8 @@ subroutine nonadiabatic_transition_scfssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurf sumdd = 0.0 do iq=1,nq do imode=1,nmodes - sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(isurface,ifre,imode,iq)) ! A - sumdd = sumdd + ABS(DD(isurface,ifre,imode,iq))**2 ! B + sumvd = sumvd + REAL(CONJG(VV(imode,iq))*DD(1,ifre,imode,iq)) ! A + sumdd = sumdd + ABS(DD(1,ifre,imode,iq))**2 ! B enddo enddo detaE = EE(isurface)-EE(ifre) ! C @@ -107,7 +107,7 @@ subroutine nonadiabatic_transition_scfssh(lfeedback,nfre,nfre_sh,nq,nmodes,isurf flagd = sumvd/sumdd*(-1.0+dsqrt(flagd)) do iq=1,nq do imode=1,nmodes - if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(isurface,ifre,imode,iq) + if(lfeedback) VV(imode,iq) = VV(imode,iq) + flagd*dd(1,ifre,imode,iq) enddo enddo isurface = ifre diff --git a/src_complex/surfacecom.f90 b/src_complex/surfacecom.f90 index 7235724..5e23b6c 100644 --- a/src_complex/surfacecom.f90 +++ b/src_complex/surfacecom.f90 @@ -79,7 +79,7 @@ module surfacecom real(kind=dp) :: E_ph_CA_sum,E_ph_QA_sum ! averager crystal energy and temperature T,in classical and quantum . - complex(kind=dpc),allocatable :: d_e(:,:,:),d_h(:,:,:),d0_e(:,:,:),d0_h(:,:,:) + complex(kind=dpc),allocatable :: d_e(:,:,:,:),d_h(:,:,:,:),d0_e(:,:,:,:),d0_h(:,:,:,:) real(kind=dp),allocatable :: g_e(:) ,g_h(:) real(kind=dp),allocatable :: g1_e(:),g1_h(:) diff --git a/src_complex/surfacehopping.f90 b/src_complex/surfacehopping.f90 index fe9e6a1..ea69498 100644 --- a/src_complex/surfacehopping.f90 +++ b/src_complex/surfacehopping.f90 @@ -91,9 +91,9 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq) if(ierr /=0) call errore('surfacehopping','Error allocating P_e_nk',1) allocate(P0_e_nk(neband,nktotf,nefre),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating P0_e_nk',1) - allocate(d_e(nefre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ajk + allocate(d_e(2,nefre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ajk,d_jak if(ierr /=0) call errore('surfacehopping','Error allocating d_e',1) - ram=real_size*nmodes*nq*nefre_sh + ram=real_size*nmodes*nq*nefre_sh*2 call print_memory("d_e",ram) allocate(dEa_dQ_e(nmodes,nq)) allocate(dEa2_dQ2_e(nmodes,nq)) @@ -149,7 +149,7 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq) if(ierr /=0) call errore('surfacehopping','Error allocating E0_e',1) allocate(P0_e(nefre,nefre),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating P0_e',1) - allocate(d0_e(nefre_sh,nmodes,nq),stat=ierr,errmsg=msg) + allocate(d0_e(2,nefre_sh,nmodes,nq),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating d0_e',1) if(methodsh == "CC-FSSH") then @@ -168,7 +168,7 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq) if(ierr /=0) call errore('surfacehopping','Error allocating P_h_nk',1) allocate(P0_h_nk(nhband,nktotf,nhfre),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating P0_h_nk',1) - allocate(d_h(nhfre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ajk + allocate(d_h(2,nhfre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ajk if(ierr /=0) call errore('surfacehopping','Error allocating d_h',1) ram=real_size*nmodes*nq*nhfre_sh call print_memory("d_h",ram) @@ -225,7 +225,7 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq) allocate(P0_h(nhfre,nhfre),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating P0_h',1) - allocate(d0_h(nhfre_sh,nmodes,nq),stat=ierr,errmsg=msg) + allocate(d0_h(2,nhfre_sh,nmodes,nq),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating d0_h',1) if(methodsh == "CC-FSSH") then @@ -356,7 +356,7 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfr real(kind=dp),intent(in) :: ee(nband*nk) complex(kind=dpc),intent(in) :: p_nk(nband,nk,nband*nk) complex(kind=dpc),intent(in) :: gmnvkq(nband,nband,nmodes,nk,nq) - complex(kind=dpc),intent(out):: dd(nfre_sh,nmodes,nq) + complex(kind=dpc),intent(out):: dd(2,nfre_sh,nmodes,nq) integer :: nfre,ifre,iq,imode integer :: ik,ikq,m,n @@ -374,8 +374,11 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfr do n=1,nband ! n epc = gmnvkq(m,n,imode,ik,iq) if(epc /= czero) then - do ifre=1,nfre - dd(ifre,imode,iq) = dd(ifre,imode,iq)+epc*CONJG(p_nk(m,ikq,isurface))*p_nk(n,ik,ifre) + do ifre=1,nfre + ! d_ajk + dd(1,ifre,imode,iq) = dd(1,ifre,imode,iq)+epc*CONJG(p_nk(m,ikq,isurface))*p_nk(n,ik,ifre) + ! d_jak + dd(2,ifre,imode,iq) = dd(2,ifre,imode,iq)+epc*CONJG(p_nk(m,ikq,ifre))*p_nk(n,ik,isurface) enddo endif enddo @@ -386,7 +389,8 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfr do ifre=1,nfre if(ifre/= isurface) then - dd(ifre,:,:) = dd(ifre,:,:)/(ee(ifre)-ee(isurface)) + dd(1,ifre,:,:) = dd(1,ifre,:,:)/(ee(ifre)-ee(isurface)) + dd(2,ifre,:,:) = dd(2,ifre,:,:)/(ee(isurface)-ee(ifre)) endif enddo @@ -416,7 +420,7 @@ subroutine calculate_hopping_probability(isurface,nfre,nfre_sh,nmodes,nq,WW,VV,d integer,intent(in) :: isurface,nfre,nfre_sh,nmodes,nq complex(kind=dpc),intent(in) :: WW(nfre) complex(kind=dpc),intent(in) :: VV(nmodes,nq) - complex(kind=dpc),intent(in) :: dd(nfre_sh,nmodes,nq) + complex(kind=dpc),intent(in) :: dd(2,nfre_sh,nmodes,nq) real(kind=dp),intent(in) :: tt real(kind=dp),intent(out) :: gg(nfre_sh) real(kind=dp),intent(out) :: gg1(nfre_sh) @@ -433,7 +437,7 @@ subroutine calculate_hopping_probability(isurface,nfre,nfre_sh,nmodes,nq,WW,VV,d sumvd = czero do iq=1,nq do imode=1,nmodes - sumvd = sumvd+VV(imode,iq)*dd(ifre,imode,iq) + sumvd = sumvd+VV(imode,iq)*dd(1,ifre,imode,iq) enddo enddo ! in adiabatic representation：the switching probabilities from the active surface isurface to another surface iefre