diff --git a/src_complex/dynamics.f90 b/src_complex/dynamics.f90 index 5721850..6a77f10 100644 --- a/src_complex/dynamics.f90 +++ b/src_complex/dynamics.f90 @@ -1,8 +1,8 @@ module dynamics use kinds,only : dp,dpc use constants,only : cone,czero - use parameters,only: temp - use epwcom,only : kqmap_sub + use parameters,only : temp + use epwcom,only : kqmap implicit none @@ -27,11 +27,12 @@ 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_sub + subroutine get_dEa_dQ(nmodes,nq,nband,nk,P_nk,gmnvkq,lit_gmnvkq,isurface,dEa_dQ) + use elph2,only : iminusq implicit none integer,intent(in) :: nmodes,nq,nband,nk integer,intent(in) :: isurface + real(kind=dp),intent(in) :: lit_gmnvkq 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) @@ -46,20 +47,18 @@ subroutine get_dEa_dQ(nmodes,nq,nband,nk,P_nk,gmnvkq,isurface,dEa_dQ) do iq=1,nq do imode=1,nmodes do ik=1,nk - ikq = kqmap_sub(ik,iq) - if(ikq /= 0) then + ikq = kqmap(ik,iq) do iband1=1,nband do iband2=1,nband epc = gmnvkq(iband1,iband2,imode,ik,iq) if(epc /= czero) then - dEa_dQ(imode,iminusq_sub(iq)) = dEa_dQ(imode,iminusq_sub(iq)) + & + dEa_dQ(imode,iminusq(iq)) = dEa_dQ(imode,iminusq(iq)) + & epc*CONJG(P_nk(iband1,ik,isurface))*P_nk(iband2,ikq,isurface) !dEa_dQ(imode,iq) = dEa_dQ(imode,iq) + & !epc*CONJG(P_nk(iband1,ik,isurface))*P_nk(iband2,ikq,isurface) endif enddo enddo - endif enddo enddo enddo @@ -137,7 +136,7 @@ subroutine rk4_nuclei(nmodes,nq,dEa_dQ,ld_gamma,wf,xx,vv,tt) ! ref : 1 J. Qiu, X. Bai, and L. Wang, The Journal of Physical Chemistry Letters 9 (2018) 4319. ! ref : "." subroutine derivs_nuclei(nmodes,nq,dEa_dQ,wf,ld_gamma,xx,vv,dx,dv) - use elph2,only : iminusq_sub + use elph2,only : iminusq implicit none integer,intent(in) :: nmodes,nq complex(kind=dpc),intent(in) :: dEa_dQ(nmodes,nq) @@ -231,7 +230,7 @@ subroutine rk4_electron_diabatic(nfre,HH,cc,cc0,dc1,dc2,dc3,dc4,tt) subroutine get_dEa2_dQ2(nmodes,nq,nfre,nfre_sh,isurface,EE,dd,dEa2_dQ2) - use elph2 , only : iminusq_sub + use elph2 , only : iminusq implicit none integer,intent(in) :: nmodes,nq,nfre,nfre_sh integer,intent(in) :: isurface @@ -252,8 +251,8 @@ subroutine get_dEa2_dQ2(nmodes,nq,nfre,nfre_sh,isurface,EE,dd,dEa2_dQ2) !(EE(isurface)-EE(ifre))*(DD(ifre,isurface,imode,iq)*DD(isurface,ifre,imode,iq)+& ! (DD(isurface,ifre,imode,iq)*DD(ifre,isurface,imode,iq))) dEa2_dQ2(imode,iq) = dEa2_dQ2(imode,iq) + & - (EE(isurface)-EE(ifre))*(DD(ifre,isurface,imode,iq)*DD(isurface,ifre,imode,iminusq_sub(iq))+& - (DD(isurface,ifre,imode,iq)*DD(ifre,isurface,imode,iminusq_sub(iq)))) + (EE(isurface)-EE(ifre))*(DD(ifre,isurface,imode,iq)*DD(isurface,ifre,imode,iminusq(iq))+& + (DD(isurface,ifre,imode,iq)*DD(ifre,isurface,imode,iminusq(iq)))) endif enddo enddo @@ -271,7 +270,7 @@ end subroutine get_dEa2_dQ2 SUBROUTINE ADD_BATH_EFFECT(nmodes,nq,wf,ld_gamma,temp,dEa2_dQ2,TT,l_ph_quantum,XX,VV) use kinds,only : dp,dpc use parameters,only : lit_ephonon,eps_acustic - use elph2,only : iminusq_sub + use elph2,only : iminusq use randoms,only : GAUSSIAN_RANDOM_NUMBER_FAST use constants,only : KB=>K_B_Ryd,sqrt3,sqrt5,sqrt7,ryd2mev,cone,ci,tpi implicit none @@ -332,16 +331,16 @@ SUBROUTINE ADD_BATH_EFFECT(nmodes,nq,wf,ld_gamma,temp,dEa2_dQ2,TT,l_ph_quantum,X ENDSUBROUTINE - subroutine pre_md(nmodes,nq,wf,ld_gamma,temp,phQ,phP,l_ph_quantum,pre_dt) + subroutine pre_md(nmodes,nqtotf,wf,ld_gamma,temp,phQ,phP,l_ph_quantum,pre_dt) use surfacecom,only : dEa_dQ,dEa2_dQ2,pre_nstep use constants,only : ry_to_fs,ryd2eV use io,only : stdout implicit none - integer ,intent(in) :: nmodes,nq + integer ,intent(in) :: nmodes,nqtotf real(kind=dp),intent(in) :: pre_dt,temp logical,intent(in) :: l_ph_quantum - real(kind=dp),intent(in) :: wf(nmodes,nq),ld_gamma(nmodes,nq) - complex(kind=dpc),intent(inout) :: phQ(nmodes,nq),phP(nmodes,nq) + real(kind=dp),intent(in) :: wf(nmodes,nqtotf),ld_gamma(nmodes,nqtotf) + complex(kind=dpc),intent(inout) :: phQ(nmodes,nqtotf),phP(nmodes,nqtotf) integer :: istep @@ -351,8 +350,8 @@ subroutine pre_md(nmodes,nq,wf,ld_gamma,temp,phQ,phP,l_ph_quantum,pre_dt) dEa_dQ = 0.0 dEa2_dQ2 = 0.0 do istep=1,pre_nstep - call rk4_nuclei(nmodes,nq,dEa_dQ,ld_gamma,wf,phQ,phP,pre_dt) - call add_bath_effect(nmodes,nq,wf,ld_gamma,temp,dEa2_dQ2,pre_dt,l_ph_quantum,phQ,phP) + call rk4_nuclei(nmodes,nqtotf,dEa_dQ,ld_gamma,wf,phQ,phP,pre_dt) + call add_bath_effect(nmodes,nqtotf,wf,ld_gamma,temp,dEa2_dQ2,pre_dt,l_ph_quantum,phQ,phP) enddo time = pre_nstep*pre_dt*ry_to_fs diff --git a/src_complex/elph2.f90 b/src_complex/elph2.f90 index eb2e84c..974d898 100644 --- a/src_complex/elph2.f90 +++ b/src_complex/elph2.f90 @@ -39,7 +39,6 @@ MODULE elph2 igk(:), &! Index for k+G vector igkq(:), &! Index for k+q+G vector iminusq(:), &! Index for -q vector - iminusq_sub(:), &! iminusk(:), &! Index for -k vector igk_k_all(:, :), &! Global index (in case of parallel) ngk_all(:), &! Global number of plane wave for each global k-point @@ -74,11 +73,9 @@ MODULE elph2 xqf(:, :), &! fine q point grid wqf(:), &! weights on the fine q grid etf(:, :), &! interpolated eigenvalues (nbnd, nkqf) - etf_sub(:,:), &! etf_k(:, :), &! Saved interpolated KS eigenenergies for later used in q-parallelization (nbnd, nkqf) etf_ks(:, :), &! interpolated eigenvalues (nbnd, nkqf) KS eigenvalues in the case of eig_read wf(:, :), &! interpolated eigenfrequencies - wf_sub(:,:), &! gmnvkq(:,:,:,:,:), &! !! g vectex gamma_all(:, :, :, :), &! Gamma gamma_nest(:, :), &! Nesting function in the case of q-parallelization diff --git a/src_complex/epwcom.f90 b/src_complex/epwcom.f90 index 55e1186..e673880 100644 --- a/src_complex/epwcom.f90 +++ b/src_complex/epwcom.f90 @@ -420,9 +420,7 @@ module klist_epw use kinds, only :dp implicit none integer, allocatable :: kqmap(:,:) - !! map of k+q grid into k grid - integer, allocatable :: kqmap_sub(:,:) - !! + !!!! map of k+q grid into k grid integer, allocatable :: isk_all(:) !! Spin index of each k-point (used in LSDA calculations only) INTEGER, ALLOCATABLE :: isk_loc(:) diff --git a/src_complex/getwcvk.f90 b/src_complex/getwcvk.f90 index 3c61fd0..4d50ba9 100644 --- a/src_complex/getwcvk.f90 +++ b/src_complex/getwcvk.f90 @@ -6,15 +6,14 @@ module getwcvk contains !ref : 1 S. Butscher et al., Physical Review B 72 (2005) !ref : 2 <固体物理> (9-29)(9-31) - subroutine get_Wcvk(ihband_min,ieband_max,nk_sub,fwhm,w_center) + subroutine get_Wcvk(ihband_min,ieband_max,fwhm,w_center) !得到光激发下垂直跃迁的跃迁几率 - use elph2,only : vmef,nkf,etf_sub !vmef(3,nbndsub,nbndsub,nkf) - use readepw,only : icbm - use surfacecom,only : indexk + use elph2,only : vmef,nkf !vmef(3,nbndsub,nbndsub,nkf) + use readepw,only : etf,icbm use io,only : stdout use constants,only : ryd2eV,ry_to_fs implicit none - integer , intent(in) :: ihband_min,ieband_max,nk_sub + integer , intent(in) :: ihband_min,ieband_max real(kind=dp),intent(in) :: fwhm real(kind=dp),intent(in) :: w_center @@ -30,21 +29,21 @@ subroutine get_Wcvk(ihband_min,ieband_max,nk_sub,fwhm,w_center) integer :: ivbm ivbm = icbm-1 - - allocate(W_cvk(icbm:ieband_max,ihband_min:ivbm,nk_sub),stat=ierr) + + allocate(W_cvk(icbm:ieband_max,ihband_min:ivbm,nkf),stat=ierr) if(ierr /=0) call errore('getmcvk','Error allocating W_cvk',1) !ref : 1 S. Butscher et al., Physical Review B 72 (2005) !ref : 1 S. Fernandez-Alberti et al., The Journal of Chemical Physics 137 (2012) fwhm_2T2 = fwhm**2.0/4.0*log(2.0) W_cvk = 0.0 - do ik=1,nk_sub + do ik=1,nkf do ibnd=ihband_min,ivbm do jbnd=icbm,ieband_max Evmef = 0.0 - E_mnk = etf_sub(jbnd,ik)-etf_sub(ibnd,ik) + E_mnk = etf(jbnd,ik)-etf(ibnd,ik) do ipol=1,3 - Evmef =Evmef+ (efield_cart(ipol)*(vmef(ipol,jbnd,ibnd,indexk(ik)))) + Evmef =Evmef+ (efield_cart(ipol)*(vmef(ipol,jbnd,ibnd,ik))) enddo fcw = f_w(E_mnk,w_center,fwhm_2T2) W_cvk(jbnd,ibnd,ik) = REAL(Evmef*CONJG(Evmef))*fcw @@ -52,7 +51,6 @@ subroutine get_Wcvk(ihband_min,ieband_max,nk_sub,fwhm,w_center) enddo enddo - !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! !% Write laser information %! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! diff --git a/src_complex/hamiltonian.f90 b/src_complex/hamiltonian.f90 index 90520be..dd00f13 100644 --- a/src_complex/hamiltonian.f90 +++ b/src_complex/hamiltonian.f90 @@ -1,14 +1,14 @@ module hamiltonian use kinds,only : dp,dpc - use io, only : msg,io_error,stdout + use io, only : msg,io_error,stdout use types use parameters - use elph2, only : nbndfst,epmatq,wf,wf_sub,nqtotf + use elph2, only : nbndfst,nk=>nktotf,epmatq,nq => nqtotf,wf use modes, only : nmodes use readepw,only : E_nk use constants,only : ryd2mev,cone,czero use surfacecom,only: lelecsh,lholesh,ieband_min,ieband_max,& - ihband_min,ihband_max,nq_sub + ihband_min,ihband_max use memory_report,only : MB,GB,complex_size, real_size,int_size,ram,print_memory implicit none @@ -24,36 +24,32 @@ 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 - real(kind=dp),allocatable :: E_e(:),E0_e(:),E_h(:),E0_h(:) + real(kind=dp),allocatable :: E_e(:),E0_e(:),E_h(:),E0_h(:) complex(kind=dp),allocatable :: P_e(:,:),P_e_nk(:,:,:),P0_e(:,:),P0_e_nk(:,:,:),& P_h(:,:),P_h_nk(:,:,:),P0_h(:,:),P0_h_nk(:,:,:) !哈密顿量的本征值与本征矢 integer :: ierr contains - subroutine allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihband_min,ihband_max) + subroutine allocate_hamiltonian(lelecsh,lholesh,ieband_min,ieband_max,ihband_min,ihband_max) implicit none logical,intent(in) :: lelecsh,lholesh - integer,intent(in) :: nk_sub,ieband_min,ieband_max,ihband_min,ihband_max + integer,intent(in) :: ieband_min,ieband_max,ihband_min,ihband_max if(lelecsh) then neband = ieband_max - ieband_min + 1 - nefre = neband * nk_sub - allocate(Enk_e(neband,nk_sub),stat=ierr,errmsg=msg) + nefre = neband * nk + allocate(Enk_e(neband,nk),stat=ierr,errmsg=msg) if(ierr /= 0) call io_error(msg) - allocate(gmnvkq_e(neband,neband,nmodes,nk_sub,nq_sub),stat=ierr,errmsg=msg) + allocate(gmnvkq_e(neband,neband,nmodes,nk,nq),stat=ierr,errmsg=msg) if(ierr /=0) then call errore('hamiltonian','Error allocating gmnvkq_e',1) call io_error(msg) endif gmnvkq_e = 0.0 - ram = complex_size*neband*neband*nmodes*nk_sub*nq_sub + ram = complex_size*neband*neband*nmodes*nk*nq call print_memory("gmnvkq_e",ram) allocate(H0_e(nefre,nefre),stat=ierr,errmsg=msg) if(ierr /=0) then @@ -63,7 +59,7 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihb H0_e = 0.0 ram = complex_size*nefre*nefre call print_memory("H0_e",ram) - allocate(H0_e_nk(neband,nk_sub,neband,nk_sub),stat=ierr,errmsg=msg) + allocate(H0_e_nk(neband,nk,neband,nk),stat=ierr,errmsg=msg) if(ierr /=0) then call errore('hamiltonian','Error allocating H0_e_nk',1) call io_error(msg) @@ -77,7 +73,7 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihb endif H_e = 0.0 call print_memory("H_e",ram) - allocate(H_e_nk(neband,nk_sub,neband,nk_sub),stat=ierr,errmsg=msg) + allocate(H_e_nk(neband,nk,neband,nk),stat=ierr,errmsg=msg) if(ierr /=0) then call errore('hamiltonian','Error allocating H_e_nk',1) call io_error(msg) @@ -89,16 +85,16 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihb if(lholesh) then nhband = ihband_max - ihband_min + 1 - nhfre = nhband * nk_sub - allocate(Enk_h(nhband,nk_sub),stat=ierr,errmsg=msg) + nhfre = nhband * nk + allocate(Enk_h(nhband,nk),stat=ierr,errmsg=msg) if(ierr /= 0) call io_error(msg) - allocate(gmnvkq_h(nhband,nhband,nmodes,nk_sub,nq_sub),stat=ierr,errmsg=msg) + allocate(gmnvkq_h(nhband,nhband,nmodes,nk,nq),stat=ierr,errmsg=msg) if(ierr /=0) then call errore('hamiltonian','Error allocating gmnvkq_h',1) call io_error(msg) endif gmnvkq_h = 0.0 - ram = complex_size*nhband*nhband*nmodes*nk_sub*nq_sub + ram = complex_size*nhband*nhband*nmodes*nk*nq call print_memory("gmnvkq_h",ram) allocate(H0_h(nhfre,nhfre),stat=ierr,errmsg=msg) if(ierr /=0) then @@ -108,7 +104,7 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihb H0_h = 0.0 ram = complex_size*nhfre*nhfre call print_memory("H0_h",ram) - allocate(H0_h_nk(nhband,nk_sub,nhband,nk_sub),stat=ierr,errmsg=msg) + allocate(H0_h_nk(nhband,nk,nhband,nk),stat=ierr,errmsg=msg) if(ierr /=0) then call errore('hamiltonian','Error allocating H0_h_nk',1) call io_error(msg) @@ -123,7 +119,7 @@ subroutine allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihb H_h = 0.0 ram = complex_size*nhfre*nhfre call print_memory("H_h",ram) - allocate(H_h_nk(nhband,nk_sub,nhband,nk_sub),stat=ierr,errmsg=msg) + allocate(H_h_nk(nhband,nk,nhband,nk),stat=ierr,errmsg=msg) if(ierr /=0) then call errore('hamiltonian','Error allocating H_h_nk',1) call io_error(msg) @@ -135,100 +131,40 @@ 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,lit_g,gnzfree) - use elph2, only : etf_sub,gmnvkq - use epwcom,only : kqmap_sub + subroutine set_H0_nk(nk,ibandmin,ibandmax,Enk,H0_nk,gmnvkq_eh) + use elph2, only : etf,gmnvkq 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 , intent(in) :: nk + integer , intent(in) :: ibandmin,ibandmax + real(kind=dp), intent(out) :: Enk(ibandmax-ibandmin+1,nk) + complex(kind=dp), intent(out) :: H0_nk(ibandmax-ibandmin+1,nk,ibandmax-ibandmin+1,nk) + complex(kind=dp), intent(out) :: gmnvkq_eh(ibandmax-ibandmin+1,ibandmax-ibandmin+1,nmodes,nk,nq) + integer :: nband - integer :: ik,iband,jband - integer :: m,n,nu,iq,ikq + integer :: ik,iband + integer :: m,n,nu,iq nband = ibandmax-ibandmin+1 H0_nk = 0.0 - do ik=1,nk_sub + do ik=1,nk do iband=1,nband - Enk(iband,ik) = etf_sub(iband+ibandmin-1,ik) - H0_nk(iband,ik,iband,ik)=etf_sub(iband+ibandmin-1,ik)*cone + Enk(iband,ik) = etf(iband+ibandmin-1,ik) + H0_nk(iband,ik,iband,ik)=etf(iband+ibandmin-1,ik)*cone enddo enddo - - do iq=1,nq_sub - do ik=1,nk_sub - gmnvkq_eh(:,:,:,ik,iq) = epmatq(ibandmin:ibandmax,ibandmin:ibandmax,:,indexk(ik),indexq(iq)) - 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 - - + + gmnvkq_eh = epmatq(ibandmin:ibandmax,ibandmin:ibandmax,:,:,:) + 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, 1981), ! (3.20) (6.4) !ref: 1 F. Giustino, Reviews of Modern Physics 89 (2017) 015003. ! (1) subroutine set_H_nk(nband,nk,nmodes,nq,ph_Q,gmnvkq_eh,H0_nk,H_nk) - use epwcom,only : kqmap_sub + use epwcom,only : kqmap implicit none integer , intent(in) :: nband,nk,nmodes,nq complex(kind=dpc),intent(in) :: ph_Q(nmodes,nq) @@ -243,17 +179,15 @@ subroutine set_H_nk(nband,nk,nmodes,nq,ph_Q,gmnvkq_eh,H0_nk,H_nk) do iq=1,nq do ik=1,nk - ikq=kqmap_sub(ik,iq) - if(ikq/=0) then - do nu=1,nmodes - do iband=1,nband !|iband,ik> - do jband=1,nband !|jband,ikq> - H_nk(jband,ikq,iband,ik) = H_nk(jband,ikq,iband,ik)+& - gmnvkq_eh(iband,jband,nu,ik,iq)*ph_Q(nu,iq) - enddo + ikq=kqmap(ik,iq) + do nu=1,nmodes + do iband=1,nband !|iband,ik> + do jband=1,nband !|jband,ikq> + H_nk(jband,ikq,iband,ik) = H_nk(jband,ikq,iband,ik)+& + gmnvkq_eh(iband,jband,nu,ik,iq)*ph_Q(nu,iq) enddo enddo - endif + enddo enddo enddo diff --git a/src_complex/initialsh.f90 b/src_complex/initialsh.f90 index 347b51e..ebbec97 100644 --- a/src_complex/initialsh.f90 +++ b/src_complex/initialsh.f90 @@ -80,13 +80,12 @@ end subroutine set_subband ! ref: 1 S. Fernandez-Alberti et al., The Journal of Chemical Physics 137 (2012) ! ref: 2. HuangKun <固体物理> (9-29) (9-31) - subroutine init_eh_KSstat(lelecsh,lholesh,llaser,nkf,init_ik,init_eband,init_hband,e_en,h_en) + subroutine init_eh_KSstat(lelecsh,lholesh,llaser,init_ik,init_eband,init_hband,e_en,h_en) use parameters,only : init_ikx,init_iky,init_ikz,init_kx,init_ky,init_kz use epwcom,only : nkf1,nkf2,nkf3 - use readepw,only : icbm,xkf - use surfacecom,only : ieband_min,ieband_max,ihband_min,ihband_max,& - indexk - use elph2,only : etf_sub !vmef(3,nbndsub,nbndsub,nkf) + use readepw,only : etf,icbm,xkf + use surfacecom,only : ieband_min,ieband_max,ihband_min,ihband_max,c_e_nk,c_h_nk + use elph2,only : vmef,ibndmin,ibndmax,nbndfst,nkf !vmef(3,nbndsub,nbndsub,nkf) use getwcvk,only: W_cvk !(W_cvk(icbm:ieband_max,ihband_min:ivbm,nkf) use constants,only :ryd2eV use io,only : stdout @@ -94,8 +93,7 @@ subroutine init_eh_KSstat(lelecsh,lholesh,llaser,nkf,init_ik,init_eband,init_hba logical , intent(in) :: lelecsh logical , intent(in) :: lholesh - logical ,intent(in) :: llaser - integer , intent(in) :: nkf + logical,intent(in) :: llaser integer,intent(out) :: init_ik integer,intent(inout) :: init_eband,init_hband real(kind=dp),intent(out) :: e_en,h_en @@ -129,13 +127,7 @@ subroutine init_eh_KSstat(lelecsh,lholesh,llaser,nkf,init_ik,init_eband,init_hba init_ikx = get_ik(init_kx,nkf1) init_iky = get_ik(init_ky,nkf2) init_ikz = get_ik(init_kz,nkf3) - init_ik = (init_ikx - 1) * nkf2 * nkf3 + (init_iky - 1) * nkf3 + init_ikz - sub:do ik=1,nkf - if(init_ik == indexk(ik)) then - init_ik = ik - exit sub - endif - enddo sub + init_ik = (init_ikx - 1) * nkf2 * nkf3 + (init_iky - 1) * nkf3 + init_ikz if(init_eband < icbm .or. init_eband > ieband_max) then write(stdout,"(/5X,A29,I5)") "Wrong! The init_eband set as:",init_eband write(stdout,"(5X,A29,I5,A16,I5)") "The init_eband need to set : ",icbm,"<= init_eband <=",ieband_max @@ -146,15 +138,15 @@ subroutine init_eh_KSstat(lelecsh,lholesh,llaser,nkf,init_ik,init_eband,init_hba endif endif - obsorbEn = (etf_sub(init_eband,init_ik)-etf_sub(init_hband,init_ik))*ryd2eV - e_en = etf_sub(init_eband,init_ik) - h_en = -1.0*etf_sub(init_hband,init_ik) + obsorbEn = (etf(init_eband,init_ik)-etf(init_hband,init_ik))*ryd2eV + e_en = etf(init_eband,init_ik) + h_en = -1.0*etf(init_hband,init_ik) write(stdout,"(/5X,A)") "Initial eh_KSstate: " - write(stdout,"(5X,A3,I5,1X,A10,3(F12.6,1X),A2)") "ik=",init_ik, "coord.: ( ",(xkf(ipol,indexk(init_ik)),ipol=1,3)," )" + write(stdout,"(5X,A3,I5,1X,A10,3(F12.6,1X),A2)") "ik=",init_ik, "coord.: ( ",(xkf(ipol,init_ik),ipol=1,3)," )" write(stdout,"(5X,A11,I5,A21,F12.5,A3)") "init_hband=",init_hband," Initial hole Energy:",& - etf_sub(init_hband,init_ik)*ryd2eV," eV" + etf(init_hband,init_ik)*ryd2eV," eV" write(stdout,"(5X,A11,I5,A21,F12.5,A3)") "init_eband=",init_eband," Initial elec Energy:",& - etf_sub(init_eband,init_ik)*ryd2eV," eV" + etf(init_eband,init_ik)*ryd2eV," eV" write(stdout,"(5X,A18,F12.5,A3)")"elec-hole energy= ",obsorbEn," eV" end subroutine init_eh_KSstat @@ -243,7 +235,7 @@ subroutine init_normalmode_coordinate_velocity(nmodes,nq,w,T,l_ph_quantum,ph_Q,p use randoms,only : gaussian_random_number !use parameters,only : lit_ephonon use surfacecom,only : nqv,E_ph_CA_sum,E_ph_QA_sum,eps_acustic - use elph2,only : iminusq_sub + use elph2,only : iminusq use io,only : stdout use constants,only : ryd2eV implicit none @@ -263,6 +255,7 @@ subroutine init_normalmode_coordinate_velocity(nmodes,nq,w,T,l_ph_quantum,ph_Q,p real(kind=dp) :: E_ph_class,E_ph_quantum integer :: iq,imode + allocate(nqv(nmodes,nq)) nqv = 0.0 do iq=1,nq @@ -278,7 +271,7 @@ subroutine init_normalmode_coordinate_velocity(nmodes,nq,w,T,l_ph_quantum,ph_Q,p E_ph_QA_sum = 0.0 do iq=1,nq - if(iminusq_sub(iq)>=iq) then + if(iminusq(iq)>=iq) then ! ph_Q(v,-q)=ph_Q(v,q)* ! ph_P(v,-q)=ph_P(v,q)* ! Ph_P(v,q) = d(ph_Q(v,q))/dt @@ -287,7 +280,7 @@ subroutine init_normalmode_coordinate_velocity(nmodes,nq,w,T,l_ph_quantum,ph_Q,p call random_number(theta1) call random_number(theta2) - if(iq==iminusq_sub(iq)) then + if(iq==iminusq(iq)) then theta1 = 0.0 theta2 = 0.0 endif @@ -303,10 +296,10 @@ subroutine init_normalmode_coordinate_velocity(nmodes,nq,w,T,l_ph_quantum,ph_Q,p else E_ph_class = K_B_Ryd*T ! IN class E_ph_CA_sum = E_ph_CA_sum + E_ph_class - if(iminusq_sub(iq)/=iq) E_ph_CA_sum = E_ph_CA_sum + E_ph_class + if(iminusq(iq)/=iq) E_ph_CA_sum = E_ph_CA_sum + E_ph_class E_ph_quantum = nqv(imode,iq)*womiga ! In Quantum E_ph_QA_sum = E_ph_QA_sum + E_ph_quantum - if(iminusq_sub(iq)/=iq) E_ph_QA_sum = E_ph_QA_sum + E_ph_quantum + if(iminusq(iq)/=iq) E_ph_QA_sum = E_ph_QA_sum + E_ph_quantum if(l_ph_quantum) then ph_Q(imode,iq) = gaussian_random_number(0.0d0,dsqrt(E_ph_quantum)/womiga)*cplx_tmp1 ph_P(imode,iq) = gaussian_random_number(0.0d0,dsqrt(E_ph_quantum))*cplx_tmp2 @@ -318,8 +311,8 @@ subroutine init_normalmode_coordinate_velocity(nmodes,nq,w,T,l_ph_quantum,ph_Q,p enddo else ! ph_Q(v,-q)=ph_Q(v,q)* ! ph_P(v,-q)=ph_P(v,q)* - ph_Q(:,iq) = CONJG(ph_Q(:,iminusq_sub(iq))) - ph_P(:,iq) = CONJG(ph_P(:,iminusq_sub(iq))) + ph_Q(:,iq) = CONJG(ph_Q(:,iminusq(iq))) + ph_P(:,iq) = CONJG(ph_P(:,iminusq(iq))) endif enddo @@ -379,193 +372,5 @@ subroutine init_surface(nfre,nfre_sh,ww,isurface) end subroutine init_surface - !==========================================================! - != get the subspace for non-adiabatic molecular dynamica =! - !==========================================================! - subroutine set_subspace(lelecsh,lholesh,llaser) - use parameters,only : init_ikx,init_iky,init_ikz,init_kx,init_ky,init_kz - use surfacecom,only : ieband_min,ieband_max,ihband_min,ihband_max - use elph2,only : nktotf,etf - implicit none - logical, intent(in) :: lelecsh,lholesh - logical, intent(in) :: llaser - - - - end subroutine - - !==========================================================! - != get the subkspace for electron diaganization =! - !==========================================================! - subroutine set_subkspace(lelecsh,lholesh,llaser,w_laser,nk_sub) - use parameters,only : init_ikx,init_iky,init_ikz,init_kx,init_ky,init_kz,& - init_eband,init_hband - use surfacecom,only : ieband_min,ieband_max,ihband_min,ihband_max,indexk - use epwcom,only : nkf1,nkf2,nkf3 - use readepw ,only : nvbmax,ncbmin,evbmax,ecbmin - use elph2,only : nktotf,etf,etf_sub,wf,wf_sub,ibndmin,ibndmax - implicit none - logical, intent(in) :: lelecsh,lholesh - logical, intent(in) :: llaser - real(kind=dp),intent(in) :: w_laser - integer, intent(out):: nk_sub - - integer :: ik,cband,vband,init_ik - real(kind=dp) :: initE_e,initE_h,maxwf - integer,allocatable :: indexk_(:) - - - allocate(indexk_(nktotf)) - - nk_sub = 0 - indexk_ = 0 - if(llaser) then - do ik=1,nktotf - if((etf(ncbmin,ik)-evbmax)<=w_laser .or. (ecbmin-etf(nvbmax,ik))<=w_laser) then - nk_sub = nk_sub + 1 - indexk_(nk_sub) = ik - endif - enddo - else - init_ikx = get_ik(init_kx,nkf1) - init_iky = get_ik(init_ky,nkf2) - init_ikz = get_ik(init_kz,nkf3) - init_ik = (init_ikx - 1) * nkf2 * nkf3 + (init_iky - 1) * nkf3 + init_ikz - maxwf = maxval(wf) - if(lelecsh .and. .not. lholesh) then - initE_e = etf(init_eband,init_ik) - do ik=1,nktotf - if(etf(ncbmin,ik)<=initE_e+maxwf) then - nk_sub = nk_sub + 1 - indexk_(nk_sub) = ik - endif - enddo - - elseif(.not. lelecsh .and. lholesh) then - initE_h = etf(init_hband,init_ik) - do ik=1,nktotf - if(etf(nvbmax,ik)>=initE_h-maxwf) then - nk_sub = nk_sub + 1 - indexk_(nk_sub) = ik - endif - enddo - elseif(lelecsh .and. lholesh) then - initE_e = etf(init_eband,init_ik) - initE_h = etf(init_hband,init_ik) - do ik=1,nktotf - if(etf(ncbmin,ik)<=initE_e+maxwf .or. etf(nvbmax,ik)>=initE_h-maxwf) then - nk_sub = nk_sub + 1 - indexk_(nk_sub) = ik - endif - enddo - endif - endif - - allocate(indexk(nk_sub)) - indexk = 0 - indexk = indexk_(1:nk_sub) - deallocate(indexk_) - - allocate(etf_sub(ibndmin:ibndmax,1:nk_sub)) - etf_sub = 0.0d0 - - do ik=1,nk_sub - etf_sub(ibndmin:ibndmax,ik) = etf(ibndmin:ibndmax,indexk(ik)) - enddo - - deallocate(etf) - - end subroutine set_subkspace - - subroutine set_subqspace(nk_sub,indexk,nq,nq_sub) - use epwcom,only : kqmap,kqmap_sub - use modes,only : nmodes - use elph2,only : wf,wf_sub,iminusq,iminusq_sub - use surfacecom,only : indexq - implicit none - integer,intent(in) :: nk_sub - integer,intent(in) :: indexk(nk_sub) - integer,intent(in) :: nq - integer,intent(out):: nq_sub - - integer :: iq,iq_,ik,ik_,ikq - - integer,allocatable :: indexq_(:) - - allocate(indexq_(nq)) - indexq_ = 0 - - - nq_sub = 0 - do iq=1,nq - KSUB:do ik=1,nk_sub - ikq=kqmap(indexk(ik),iq) - do ik_=1,nk_sub - if(ikq==indexk(ik_)) then - nq_sub = nq_sub + 1 - indexq_(nq_sub) = iq - exit KSUB - endif - enddo - enddo KSUB - enddo - - allocate(indexq(nq_sub)) - indexq = 0 - indexq = indexq_(1:nq_sub) - allocate(wf_sub(nmodes,nq_sub)) - wf_sub = 0.0 - - do iq=1,nq_sub - wf_sub(:,iq) = wf(:,indexq(iq)) - enddo - - deallocate(wf) - - allocate(iminusq_sub(nq_sub)) - iminusq_sub = 0 - do iq=1,nq_sub - sub :do iq_=1,nq_sub - if(iminusq(indexq(iq)) == indexq(iq_)) then - iminusq_sub(iq) = iq_ - exit sub - endif - enddo sub - enddo - deallocate(indexq_) - - allocate(kqmap_sub(nk_sub,nq_sub)) - kqmap_sub = 0 - do iq=1,nq_sub - do ik=1,nk_sub - ikq=kqmap(indexk(ik),indexq(iq)) - sub1 : do ik_=1,nk_sub - if(ikq==indexk(ik_)) then - kqmap_sub(ik,iq) = ik_ - exit sub1 - endif - enddo sub1 - enddo - enddo - end subroutine set_subqspace - - !==========================================================! - != get the subkspace of gmnvkq to build the hamiltoian =! - !==========================================================! - subroutine set_subgmnvkq() - - end subroutine - - !==========================================================! - != get the phonon subkspace for phonon dynamica =! - !==========================================================! - subroutine set_subphonon() - - end subroutine - - - - - end module initialsh \ No newline at end of file diff --git a/src_complex/int_to_char.f90 b/src_complex/int_to_char.f90 deleted file mode 100644 index a36bfd4..0000000 --- a/src_complex/int_to_char.f90 +++ /dev/null @@ -1,61 +0,0 @@ -! -! Copyright (C) 2009 Quantum ESPRESSO groups -! This file is distributed under the terms of the -! GNU General Public License. See the file `License' -! in the root directory of the present distribution, -! or http://www.gnu.org/copyleft/gpl.txt . -! - !----------------------------------------------------------------------- - FUNCTION int_to_char( i ) - !----------------------------------------------------------------------- - !! Converts an integer number of up to 6 figures into a left-justifed - !! character variable. - ! - IMPLICIT NONE - ! - INTEGER, INTENT(IN) :: i - CHARACTER (LEN=6) :: int_to_char - CHARACTER :: c - INTEGER :: n, j, nc - LOGICAL :: neg - ! - nc = 6 - ! - IF( i < 0 ) then - nc = nc - 1 - n = -i - neg = .true. - ELSE - n = i - neg = .false. - END IF - ! - j = 1 - DO WHILE( j <= nc ) - int_to_char(j:j) = CHAR( MOD( n, 10 ) + ICHAR( '0' ) ) - n = n / 10 - IF( n == 0 ) EXIT - j = j + 1 - END DO - ! - IF( j <= nc ) THEN - DO n = 1, j/2 - c = int_to_char( n : n ) - int_to_char( n : n ) = int_to_char( j-n+1 : j-n+1 ) - int_to_char( j-n+1 : j-n+1 ) = c - END DO - IF( j < nc ) int_to_char(j+1:nc) = ' ' - ELSE - int_to_char(:) = '*' - END IF - ! - IF( neg ) THEN - DO n = nc+1, 2, -1 - int_to_char(n:n) = int_to_char(n-1:n-1) - END DO - int_to_char(1:1) = '-' - END IF - ! - RETURN - ! - END FUNCTION int_to_char diff --git a/src_complex/lvcsh.f90 b/src_complex/lvcsh.f90 index af3c099..52c7431 100644 --- a/src_complex/lvcsh.f90 +++ b/src_complex/lvcsh.f90 @@ -39,25 +39,21 @@ 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,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 - + 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 use sortting,only : resort_eigen_energy_stat use randoms,only : init_random_seed use lasercom,only : fwhm,w_laser use getwcvk,only : get_Wcvk use initialsh,only : set_subband,init_normalmode_coordinate_velocity, & - init_eh_KSstat,init_stat_adiabatic,init_surface, & - set_subkspace,set_subqspace + init_eh_KSstat,init_stat_adiabatic,init_surface use write_sh_information,only : write_initial_information use surfacecom,only : methodsh,lfeedback,lit_gmnvkq,naver,nsnap,nstep,dt,& pre_nstep,pre_dt,l_ph_quantum,gamma,temp,iaver, & isnap,istep,ldecoherence,Cdecoherence, & - lelecsh,lholesh,ieband_min,ieband_max,nk_sub, & - indexk,nq_sub,indexq,& + lelecsh,lholesh,ieband_min,ieband_max, & ihband_min,ihband_max,nefre_sh,nhfre_sh,iesurface, & ihsurface,iesurface_j,ihsurface_j,c_e,c_e_nk,d_e, & d0_e,c_h,c_h_nk,d_h,d0_h,dEa_dQ, & @@ -76,11 +72,10 @@ 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 + use elph2,only : wf,xkf,nqtotf,nktotf,nbndfst,gmnvkq,etf use modes,only : nmodes use cell_base,only : bg use date_and_times,only : get_date_and_time @@ -120,42 +115,32 @@ program lvcsh call readepwout(epwoutname) call set_subband(lelecsh,lholesh,ieband_min,ieband_max,ihband_min,ihband_max) !get ieband_min,ieband_max,ihband_min,ihband_max - call set_subkspace(lelecsh,lholesh,llaser,w_laser,nk_sub) - ! get nk_sub and indexk(nk_sub), etf_sub(:,nk_sub) - call set_subqspace(nk_sub,indexk,nqtotf,nq_sub) - ! get nq_sub and indexq(nq_sub) wf_sub(nmodes,nq_sub) iminusq_sub(nq_sub) kqmap_sub(nk_sub,nq_sub) - - - call allocate_hamiltonian(lelecsh,lholesh,nk_sub,ieband_min,ieband_max,ihband_min,ihband_max) + call allocate_hamiltonian(lelecsh,lholesh,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,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) + call set_H0_nk(nktotf,ieband_min,ieband_max,Enk_e,H0_e_nk,gmnvkq_e) + !set H0_e_nk and gmnvkq_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,lit_gmnvkq,gnzfree_h) + call set_H0_nk(nktotf,ihband_min,ihband_max,Enk_h,H0_h_nk,gmnvkq_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) - allocate(gijk_h(gnzfree_h)) + !set H0_h_nk and gmnvkq_h if(nhfre_sh == 0 .or. nhfre_sh > nhfre) nhfre_sh = nhfre endif deallocate(gmnvkq) - deallocate(epmatq) - call allocatesh(methodsh,lelecsh,lholesh,nk_sub,nmodes,nq_sub) + call allocatesh(methodsh,lelecsh,lholesh,nmodes,nqtotf) if(trim(calculation) == "lvcsh") then ! set the friction coefficient of Langevin dynamica of all phonon modes. - call set_gamma(nmodes,nq_sub,gamma,ld_fric,wf_sub,ld_gamma) + call set_gamma(nmodes,nqtotf,gamma,ld_fric,wf,ld_gamma) - if(llaser) call get_Wcvk(ihband_min,ieband_max,nk_sub,fwhm,w_laser) + if(llaser) call get_Wcvk(ihband_min,ieband_max,fwhm,w_laser) ! ref : https://journals.aps.org/prb/pdf/10.1103/PhysRevB.72.045314 !get W_cvk(icband,ivband,ik) @@ -175,46 +160,45 @@ program lvcsh !==================! !!Get the initial normal mode coordinate phQ and versity phP - call init_normalmode_coordinate_velocity(nmodes,nq_sub,wf_sub,temp,l_ph_quantum,phQ,phP) + call init_normalmode_coordinate_velocity(nmodes,nqtotf,wf,temp,l_ph_quantum,phQ,phP) !应该先跑平衡后，再做电子空穴动力学计算 - call pre_md(nmodes,nq_sub,wf_sub,ld_gamma,temp,phQ,phP,l_ph_quantum,pre_dt) + call pre_md(nmodes,nqtotf,wf,ld_gamma,temp,phQ,phP,l_ph_quantum,pre_dt) !!得到初始电子和空穴的初始的KS状态 init_ik,init_eband,init_hband(in the diabatic states) - call init_eh_KSstat(lelecsh,lholesh,llaser,nk_sub,init_ik,init_eband,init_hband,init_e_en,init_h_en) + call init_eh_KSstat(lelecsh,lholesh,llaser,init_ik,init_eband,init_hband,init_e_en,init_h_en) if(lelecsh) then - call set_H_nk(neband,nk_sub,nmodes,nq_sub,phQ,gmnvkq_e,H0_e_nk,H_e_nk) + call set_H_nk(neband,nktotf,nmodes,nqtotf,phQ,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) - P_e_nk = reshape(P_e,(/ neband,nk_sub,nefre /)) + P_e_nk = reshape(P_e,(/ neband,nktotf,nefre /)) call init_stat_adiabatic(nefre,E_e,nefre_sh,init_e_en,iesurface) w_e = czero w_e(iesurface) = cone call convert_adiabatic_diabatic(nefre,P_e,w_e,c_e) - c_e_nk = reshape(c_e,(/ neband,nk_sub /)) + c_e_nk = reshape(c_e,(/ neband,nktotf /)) - 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) + call calculate_nonadiabatic_coupling(nmodes,nqtotf,neband,nktotf,E_e,P_e_nk,gmnvkq_e,lit_gmnvkq,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 if(lholesh) then - call set_H_nk(nhband,nk_sub,nmodes,nq_sub,phQ,gmnvkq_h,H0_h_nk,H_h_nk) + call set_H_nk(nhband,nktotf,nmodes,nqtotf,phQ,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) - P_h_nk = reshape(P_h,(/ nhband,nk_sub,nhfre /)) + P_h_nk = reshape(P_h,(/ nhband,nktotf,nhfre /)) call init_stat_adiabatic(nhfre,E_h,nhfre_sh,init_h_en,ihsurface) w_h = czero w_h(ihsurface) = cone call convert_adiabatic_diabatic(nefre,P_h,w_h,c_h) - c_h_nk = reshape(c_h,(/ nhband,nk_sub /)) + c_h_nk = reshape(c_h,(/ nhband,nktotf /)) - call calculate_nonadiabatic_coupling(nmodes,nq_sub,nhband,nk_sub,E_h,P_h_nk,gmnvkq_h,nhfre_sh,d_h) + call calculate_nonadiabatic_coupling(nmodes,nqtotf,nhband,nktotf,E_h,P_h_nk,gmnvkq_h,lit_gmnvkq,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 @@ -222,7 +206,7 @@ program lvcsh phQ0=phQ; phP0=phP - call write_initial_information(iaver,nmodes,nq_sub,wf_sub,phQ,phP) + call write_initial_information(iaver,nmodes,nqtotf,wf,phQ,phP) !=======================! != loop over snapshots =! @@ -240,11 +224,11 @@ program lvcsh !dEa_dQ in time t0 if(l_dEa_dQ) then if(lelecsh) then - call get_dEa_dQ(nmodes,nq_sub,neband,nk_sub,P0_e_nk,gmnvkq_e,iesurface,dEa_dQ_e) + call get_dEa_dQ(nmodes,nqtotf,neband,nktotf,P0_e_nk,gmnvkq_e,lit_gmnvkq,iesurface,dEa_dQ_e) dEa_dQ = dEa_dQ + dEa_dQ_e endif if(lholesh) then - call get_dEa_dQ(nmodes,nq_sub,nhband,nk_sub,P0_h_nk,gmnvkq_h,ihsurface,dEa_dQ_h) + call get_dEa_dQ(nmodes,nqtotf,nhband,nktotf,P0_h_nk,gmnvkq_h,lit_gmnvkq,ihsurface,dEa_dQ_h) dEa_dQ = dEa_dQ + dEa_dQ_h endif endif @@ -254,7 +238,7 @@ program lvcsh !==============================! !use rk4 to calculate the dynamical of phonon normal modes !update phQ,phP to time t0+dt - call rk4_nuclei(nmodes,nq_sub,dEa_dQ,ld_gamma,wf_sub,phQ,phP,dt) + call rk4_nuclei(nmodes,nqtotf,dEa_dQ,ld_gamma,wf,phQ,phP,dt) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! @@ -264,24 +248,24 @@ program lvcsh !electron wave function is propagated in diabatic representation !hamiltonian in time t0 - call set_H_nk(neband,nk_sub,nmodes,nq_sub,phQ0,gmnvkq_e,H0_e_nk,H_e_nk) + call set_H_nk(neband,nktotf,nmodes,nqtotf,phQ0,gmnvkq_e,H0_e_nk,H_e_nk) H_e = reshape(H_e_nk,(/ nefre,nefre /)) !update c_e to time t0+dt call rk4_electron_diabatic(nefre,H_e,c_e,cc0_e,dc1_e,dc2_e,dc3_e,dc4_e,dt) ! update H_nk in time t0+dt - call set_H_nk(neband,nk_sub,nmodes,nq_sub,phQ,gmnvkq_e,H0_e_nk,H_e_nk) + call set_H_nk(neband,nktotf,nmodes,nqtotf,phQ,gmnvkq_e,H0_e_nk,H_e_nk) 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) if(lsortpes) call resort_eigen_energy_stat(nefre,E_e,P_e,E0_e,P0_e,mix_thr) - P_e_nk = reshape(P_e,(/ neband,nk_sub,nefre /)) + P_e_nk = reshape(P_e,(/ neband,nktotf,nefre /)) ! 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,nefre_sh,d_e) + call calculate_nonadiabatic_coupling(nmodes,nqtotf,neband,nktotf,E_e,p_e,gmnvkq_e,lit_gmnvkq,nefre_sh,d_e) ! use p_e in time t0+dt, to convert c_e(t0+dt) to w_e(t0+dt) @@ -295,7 +279,7 @@ program lvcsh ! use FSSH calculation hopping probability in adiabatic representation,get g_e,g1_e - call calculate_hopping_probability(iesurface,nefre,nefre_sh,nmodes,nq_sub,w0_e,phP0,d0_e,dt,g_e,g1_e) + call calculate_hopping_probability(iesurface,nefre,nefre_sh,nmodes,nqtotf,w0_e,phP0,d0_e,dt,g_e,g1_e) !dealwith trilvial crossing,by fixed ge if(methodsh == "SC-FSSH") then @@ -310,11 +294,11 @@ program lvcsh !% change potential energy surface %! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! if(methodsh == "FSSH") then - call nonadiabatic_transition_fssh(lfeedback,nefre,nefre_sh,nq_sub,nmodes,iesurface,E0_e,P0_e,d0_e,g_e,phP) + call nonadiabatic_transition_fssh(lfeedback,nefre,nefre_sh,nqtotf,nmodes,iesurface,E0_e,P0_e,d0_e,g_e,phP) elseif( methodsh == "SC-FSSH") then - call nonadiabatic_transition_scfssh(lfeedback,nefre,nefre_sh,nq_sub,nmodes,iesurface,E0_e,P0_e,d0_e,g_e,phP) + call nonadiabatic_transition_scfssh(lfeedback,nefre,nefre_sh,nqtotf,nmodes,iesurface,E0_e,P0_e,d0_e,g_e,phP) elseif(methodsh == "CC-FSSH") then - call nonadiabatic_transition_ccfssh(lfeedback,nefre,nefre_sh,nq_sub,nmodes,iesurface,iesurface_j,& + call nonadiabatic_transition_ccfssh(lfeedback,nefre,nefre_sh,nqtotf,nmodes,iesurface,iesurface_j,& &esurface_type,E0_e,P0_e,P_e,d0_e,S_bi_e,g_e,phP) endif @@ -324,24 +308,24 @@ program lvcsh !hole wave function is propagated in diabatic representation !hamiltonian in time t0 - call set_H_nk(nhband,nk_sub,nmodes,nq_sub,phQ0,gmnvkq_h,H0_h_nk,H_h_nk) + call set_H_nk(nhband,nktotf,nmodes,nqtotf,phQ0,gmnvkq_h,H0_h_nk,H_h_nk) H_h = reshape(H_h_nk,(/ nhfre,nhfre /)) !update c_h to time t0+dt call rk4_electron_diabatic(nhfre,H_h,c_h,cc0_h,dc1_h,dc2_h,dc3_h,dc4_h,dt) ! update H_nk in time t0+dt - call set_H_nk(nhband,nk_sub,nmodes,nq_sub,phQ,gmnvkq_h,H0_h_nk,H_h_nk) + call set_H_nk(nhband,nktotf,nmodes,nqtotf,phQ,gmnvkq_h,H0_h_nk,H_h_nk) 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) if(lsortpes) call resort_eigen_energy_stat(nhfre,E_h,P_h,E0_h,P0_h,mix_thr) - P_h_nk = reshape(P_h,(/ nhband,nk_sub,nhfre /)) + P_h_nk = reshape(P_h,(/ nhband,nktotf,nhfre /)) ! 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,nhfre_sh,d_h) + call calculate_nonadiabatic_coupling(nmodes,nqtotf,nhband,nktotf,E_h,p_h,gmnvkq_h,lit_gmnvkq,nhfre_sh,d_h) ! use p_h in time t0+dt, to convert c_h(t0+dt) to w_h(t0+dt) @@ -353,7 +337,7 @@ program lvcsh ! use FSSH calculation hopping probability in adiabatic representation - call calculate_hopping_probability(ihsurface,nhfre,nhfre_sh,nmodes,nq_sub,w0_h,phP0,d0_h,dt,g_h,g1_h) + call calculate_hopping_probability(ihsurface,nhfre,nhfre_sh,nmodes,nqtotf,w0_h,phP0,d0_h,dt,g_h,g1_h) !dealwith trilvial crossing,fixed ge if(methodsh == "SC-FSSH") then @@ -368,11 +352,11 @@ program lvcsh !% change potential energy surface %! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! if(methodsh == "FSSH") then - call nonadiabatic_transition_fssh(lfeedback,nhfre,nhfre_sh,nq_sub,nmodes,ihsurface,E0_h,P0_h,d0_h,g_h,phP) + call nonadiabatic_transition_fssh(lfeedback,nhfre,nhfre_sh,nqtotf,nmodes,ihsurface,E0_h,P0_h,d0_h,g_h,phP) elseif( methodsh == "SC-FSSH") then - call nonadiabatic_transition_scfssh(lfeedback,nhfre,nhfre_sh,nq_sub,nmodes,ihsurface,E0_h,P0_h,d0_h,g_h,phP) + call nonadiabatic_transition_scfssh(lfeedback,nhfre,nhfre_sh,nqtotf,nmodes,ihsurface,E0_h,P0_h,d0_h,g_h,phP) elseif(methodsh == "CC-FSSH") then - call nonadiabatic_transition_ccfssh(lfeedback,nhfre,nhfre_sh,nq_sub,nmodes,ihsurface,ihsurface_j,& + call nonadiabatic_transition_ccfssh(lfeedback,nhfre,nhfre_sh,nqtotf,nmodes,ihsurface,ihsurface_j,& &hsurface_type,E0_h,P0_h,P_h,d0_h,S_bi_h,g_h,phP) endif @@ -386,11 +370,11 @@ program lvcsh !dEa2_dQ2 in time t0 if(l_dEa2_dQ2) then if(lelecsh) then - call get_dEa2_dQ2(nmodes,nq_sub,nefre,nefre_sh,iesurface,E0_e,d0_e,dEa2_dQ2_e) + 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,nq_sub,nhfre,nhfre_sh,ihsurface,E0_h,d0_h,dEa2_dQ2_h) + 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 endif @@ -398,7 +382,7 @@ program lvcsh !===================! != add bath effect =! !===================! - call add_bath_effect(nmodes,nq_sub,wf_sub,ld_gamma,temp,dEa2_dQ2,dt,l_ph_quantum,phQ,phP) + call add_bath_effect(nmodes,nqtotf,wf,ld_gamma,temp,dEa2_dQ2,dt,l_ph_quantum,phQ,phP) !============================! @@ -417,7 +401,7 @@ program lvcsh !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! !% Write non-adiabatic dynamica energy information every step %! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%! - phU = 0.5*(wf_sub**2)*(phQ*CONJG(phQ)) + phU = 0.5*(wf**2)*(phQ*CONJG(phQ)) phK = 0.5*phP*CONJG(phP) SUM_phK = SUM(phK) SUM_phU = SUM(phU) @@ -460,7 +444,7 @@ program lvcsh !=====================! != store information =! !=====================! - do iq=1,nq_sub + do iq=1,nqtotf do imode=1,nmodes phQsit(imode,iq,isnap) = phQsit(imode,iq,isnap)+phQ(imode,iq) phPsit(imode,iq,isnap) = phPsit(imode,iq,isnap)+phP(imode,iq) @@ -526,17 +510,17 @@ program lvcsh !====================! != save information =! !====================! - call save_phQ(nmodes,nq_sub,nsnap,phQsit) - call save_phP(nmodes,nq_sub,nsnap,phPsit) - call save_phK(nmodes,nq_sub,nsnap,phKsit) - call save_phU(nmodes,nq_sub,nsnap,phUsit) + call save_phQ(nmodes,nqtotf,nsnap,phQsit) + call save_phP(nmodes,nqtotf,nsnap,phPsit) + call save_phK(nmodes,nqtotf,nsnap,phKsit) + call save_phU(nmodes,nqtotf,nsnap,phUsit) if(lelecsh) then call save_pes(nefre,nsnap,naver,pes_one_e,pes_e,pes_e_file) call save_csit(nefre,nsnap,naver,csit_e,csit_e_file) call save_wsit(nefre,nsnap,naver,wsit_e,wsit_e_file) call save_psit(nefre,nsnap,naver,psit_e,psit_e_file) - call plot_band_occupatin_withtime(neband,nk_sub,Enk_e,nsnap,psit_e,csit_e,savedsnap,band_e_file) + call plot_band_occupatin_withtime(neband,nktotf,Enk_e,xkf,nsnap,psit_e,csit_e,savedsnap,band_e_file) endif if(lholesh) then @@ -544,7 +528,7 @@ program lvcsh call save_csit(nhfre,nsnap,naver,csit_h,csit_h_file) call save_wsit(nhfre,nsnap,naver,wsit_h,wsit_h_file) call save_psit(nhfre,nsnap,naver,psit_h,psit_h_file) - call plot_band_occupatin_withtime(nhband,nk_sub,Enk_h,nsnap,psit_h,csit_h,savedsnap,band_h_file) + call plot_band_occupatin_withtime(nhband,nktotf,Enk_h,xkf,nsnap,psit_h,csit_h,savedsnap,band_h_file) endif @@ -561,10 +545,10 @@ program lvcsh do inode=1,nnode write(stdout,"(A,I4)") "Read information in node:",inode do icore=1,ncore - call read_phQ(inode,icore,nmodes,nq_sub,nsnap,phQsit) - call read_phP(inode,icore,nmodes,nq_sub,nsnap,phPsit) - call read_phK(inode,icore,nmodes,nq_sub,nsnap,phKsit) - call read_phU(inode,icore,nmodes,nq_sub,nsnap,phUsit) + call read_phQ(inode,icore,nmodes,nqtotf,nsnap,phQsit) + call read_phP(inode,icore,nmodes,nqtotf,nsnap,phPsit) + call read_phK(inode,icore,nmodes,nqtotf,nsnap,phKsit) + call read_phU(inode,icore,nmodes,nqtotf,nsnap,phUsit) if(lelecsh) then call read_pes( inode,icore,nefre,nsnap,naver,pes_one_e,pes_e,pes_e_file) @@ -606,34 +590,34 @@ program lvcsh !!plot !!!! - call plot_phQ(nmodes,nq_sub,nsnap,phQsit) - call plot_phP(nmodes,nq_sub,nsnap,phPsit) - call plot_phK(nmodes,nq_sub,nsnap,phKsit) - call plot_phU(nmodes,nq_sub,nsnap,phUsit) - call plot_ph_temp(nmodes,nq_sub,nsnap,phKsit,phUsit) + call plot_phQ(nmodes,nqtotf,nsnap,phQsit) + call plot_phP(nmodes,nqtotf,nsnap,phPsit) + call plot_phK(nmodes,nqtotf,nsnap,phKsit) + call plot_phU(nmodes,nqtotf,nsnap,phUsit) + call plot_ph_temp(nmodes,nqtotf,nsnap,phKsit,phUsit) deallocate(phPsit,phQsit,phKsit,phUsit) if(lelecsh) then write(stdout,"(/,A)") "Plotting electron non-adiabatic dynamica Information." - call plot_pes(nefre,nefre_sh,nsnap,pes_one_e,pes_e,wsit_e,pes_e_file) + call plot_pes(nefre,nsnap,pes_one_e,pes_e,wsit_e,pes_e_file) call plot_csit(nefre,nsnap,naver,csit_e,csit_e_file) - call plot_wsit(nefre,nefre_sh,nsnap,naver,wsit_e,wsit_e_file) + call plot_wsit(nefre,nsnap,naver,wsit_e,wsit_e_file) call plot_psit(nefre,nsnap,naver,psit_e,psit_e_file) - call plot_band_occupatin_withtime(neband,nk_sub,Enk_e,nsnap,psit_e,csit_e,savedsnap,band_e_file) - call plot_eh_temp(neband,nk_sub,Enk_e,nsnap,psit_e,csit_e,savedsnap,e_temp_file) + call plot_band_occupatin_withtime(neband,nktotf,Enk_e,xkf,nsnap,psit_e,csit_e,savedsnap,band_e_file) + call plot_eh_temp(neband,nktotf,Enk_e,xkf,nsnap,psit_e,csit_e,savedsnap,e_temp_file) deallocate(pes_e,pes_one_e,csit_e,wsit_e,psit_e) endif if(lholesh) then write(stdout,"(/,A)") "Plotting hole non-adiabatic dynamica Information." - call plot_pes(nhfre,nhfre_sh,nsnap,pes_one_h,pes_h,wsit_h,pes_h_file) + call plot_pes(nhfre,nsnap,pes_one_h,pes_h,wsit_h,pes_h_file) call plot_csit(nhfre,nsnap,naver,csit_h,csit_h_file) - call plot_wsit(nhfre,nhfre_sh,nsnap,naver,wsit_h,wsit_h_file) + call plot_wsit(nhfre,nsnap,naver,wsit_h,wsit_h_file) call plot_psit(nhfre,nsnap,naver,psit_h,psit_h_file) - call plot_band_occupatin_withtime(nhband,nk_sub,Enk_h,nsnap,psit_h,csit_h,savedsnap,band_h_file) - call plot_eh_temp(nhband,nk_sub,Enk_h,nsnap,psit_h,csit_h,savedsnap,h_temp_file) + call plot_band_occupatin_withtime(nhband,nktotf,Enk_h,xkf,nsnap,psit_h,csit_h,savedsnap,band_h_file) + call plot_eh_temp(nhband,nktotf,Enk_h,xkf,nsnap,psit_h,csit_h,savedsnap,h_temp_file) deallocate(pes_h,pes_one_h,csit_h,wsit_h,psit_h) endif diff --git a/src_complex/readepw.f90 b/src_complex/readepw.f90 index 895f9ca..e138a66 100644 --- a/src_complex/readepw.f90 +++ b/src_complex/readepw.f90 @@ -126,7 +126,7 @@ subroutine readepwout(fepwout) !logical :: eig_read,epbread,epbwrite,efermi_read LOGICAL :: already_skipped !! Skipping band during the Wannierization - integer :: nbndskip = 0 + integer :: nbndskip logical :: wannierize integer :: itmp,count_piv_spin @@ -981,11 +981,14 @@ subroutine readepwout(fepwout) ef = ef /ryd2ev WRITE(stdout,'(/5x,a,f10.6,a)') 'Fermi energy coarse grid = ', ef * ryd2ev, ' eV' + !if(nelec == 0.0 .and. ieband_max==0 .and. ihband_min==0) then + ! write(stdout,"(5X,A,F8.4,A)") "WARNING! The nelec =",nelec,"and ieband_max=0 ihband_min=0" + ! write(stdout,"(5X,A)") "Need to set nelec right in LVCSH.in .OR. set lreadscfout= .true." + !endif !IF (efermi_read) THEN read(unitepwout,"(/5X,A)") ctmp read(unitepwout,"(/5X,A)") ctmp - write(stdout,"(/5X,A)") ctmp if(ctmp(1:47)=="Fermi energy is read from the input file: Ef = ") then efermi_read = .true. backspace(unitepwout) @@ -1088,8 +1091,8 @@ subroutine readepwout(fepwout) ENDIF endif - !WRITE(stdout, '(/5x," icbm(conductor band mim) = ",i6)' ) icbm - !WRITE(stdout, '(5x," ivbm(valence band max) = ",i6)' ) icbm-1 + WRITE(stdout, '(/5x," icbm(conductor band mim) = ",i6)' ) icbm + WRITE(stdout, '(5x," ivbm(valence band max) = ",i6)' ) icbm-1 ! @@ -1162,7 +1165,6 @@ subroutine readepwout(fepwout) call errore('readepw','Error allocating gmnvkq',1) call io_error(msg) endif - gmnvkq = 0.0 ram = real_size*(ibndmax-ibndmin+1)**2*nmodes*nktotf*nqtotf call print_memory("gmnvkq",ram) @@ -1171,7 +1173,6 @@ subroutine readepwout(fepwout) call errore('readepw','Error allocating epmatq',1) call io_error(msg) endif - epmatq = czero ram = complex_size*(ibndmax-ibndmin+1)**2*nmodes*nktotf*nqtotf call print_memory("epmatq",ram) @@ -1319,24 +1320,17 @@ subroutine readepwout(fepwout) ecbmin = Minval(etf(ncbmin,:)) WRITE(stdout,'(/14x,a,i5,2x,a,f9.3,a)') 'Valence band max = ', nvbmax, 'evbmax = ', evbmax , ' eV' WRITE(stdout,'(14x,a,i5,2x,a,f9.3,a/)') 'Conductor band min = ', ncbmin, 'ecbmin = ', ecbmin , ' eV' - icbm = ncbmin - !if(icbm /= ncbmin) write(stdout,"(5X,A)") "Warning! The nelec need to be set right." + if(icbm /= ncbmin) write(stdout,"(5X,A)") "Warning! The nelec need to be set right." etf = etf/ryd2eV evbmax = evbmax/ryd2eV ecbmin = ecbmin/ryd2eV - etf = etf - evbmax - ecbmin = ecbmin - evbmax - evbmax = evbmax - evbmax - wf = wf/ryd2mev gmnvkq = gmnvkq/ryd2mev epmatq = epmatq/ryd2mev eps_acustic = eps_acustic/ryd2mev - - do iq=1,nqtotf do nu = 1,nmodes gmnvkq(:,:,nu,:,iq)=sqrt(2.0*wf(nu,iq)/nqtotf)*gmnvkq(:,:,nu,:,iq) diff --git a/src_complex/saveinf.f90 b/src_complex/saveinf.f90 index c86f594..eb8d3ff 100644 --- a/src_complex/saveinf.f90 +++ b/src_complex/saveinf.f90 @@ -1,10 +1,10 @@ module saveinf use kinds,only : dp,dpc - use surfacecom,only : dt,nstep,temp,nqv,l_ph_quantum,eps_acustic,indexk,indexq + use surfacecom,only : dt,nstep,temp,nqv,l_ph_quantum,eps_acustic use io,only : io_file_unit,open_file,close_file,stdout use constants,only : ry_to_fs,maxlen,RYTOEV,ryd2meV,K_B_Ryd - use elph2,only : wf_sub - use parameters,only : outdir + use elph2,only : wf,xqf + use parameters,only : outdir implicit none integer :: iaver,isnap,ifre,iq,imode,inode,icore integer :: i @@ -23,10 +23,10 @@ module saveinf character(len=maxlen) :: phK_file = "phKsit.dat" character(len=maxlen) :: phT_file = "phLOtemp.dat" character(len=maxlen) :: band_e_file = "band_e" - character(len=maxlen) :: band_h_file = "band_h" + character(len=maxlen) :: band_h_file = "band_h" character(len=maxlen) :: e_temp_file = "temp_e" character(len=maxlen) :: h_temp_file = "temp_h" - + character(len=12) :: ctmp1,ctmp2 character(len=12),allocatable :: cphmode(:,:) character(len=12),allocatable :: cefree(:),cnk(:,:) @@ -34,13 +34,13 @@ module saveinf contains - !Write average phQ infortmation to the file: phQsit.dat + !Write average phQ infortmation to the file: phQsit.dat subroutine save_phQ(nmodes,nq,nsnap,phQsit) integer,intent(in) :: nmodes,nq,nsnap complex(kind=dpc),intent(in) :: phQsit(nmodes,nq,0:nsnap) integer :: phq_unit - character(len=maxlen) :: phQ_file_ + character(len=maxlen) :: phQ_file_ phQ_file_ = trim(outdir)//trim(adjustl(phQ_file)) phq_unit = io_file_unit() call open_file(phQ_file_,phq_unit) @@ -54,60 +54,60 @@ subroutine save_phQ(nmodes,nq,nsnap,phQsit) enddo enddo - write(phq_unit,"(A)") "Average of Normal mode coordinate for one core sample. ((phQ(imode,iq),imode=1,nmodes),iq=1,nq)" - write(phq_unit,"(A12,*(2(1X,A12)))") "time",(("Re[Q"//trim(adjustl(cphmode(imode,iq)))//"]",& + write(phq_unit,"(A)") "Average of Normal mode coordinate for one core sample. ((phQ(imode,iq),imode=1,nmodes),iq=1,nq)" + write(phq_unit,"(A12,*(2(1X,A12)))") "time",(("Re[Q"//trim(adjustl(cphmode(imode,iq)))//"]",& "Im[Q"//trim(adjustl(cphmode(imode,iq)))//"]",imode=1,nmodes),iq=1,nq) write(phq_unit,"(A12,*(2(1X,A12)))") "fs ",((" a.u. "," a.u. ",imode=1,nmodes),iq=1,nq) - write(phq_unit,"(A12,*(2(1X,F12.5)))") "Omega(meV)",(((wf_sub(imode,iq)*ryd2meV,wf_sub(imode,iq)*ryd2meV),imode=1,nmodes),iq=1,nq) + write(phq_unit,"(A12,*(2(1X,F12.5)))") "Omega(meV)",(((wf(imode,iq)*ryd2meV,wf(imode,iq)*ryd2meV),imode=1,nmodes),iq=1,nq) do isnap=0,nsnap write(phq_unit,"(F12.5,*(2(1X,E12.5)))") dt*nstep*isnap*ry_to_fs,((phQsit(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo call close_file(phQ_file_,phq_unit) - - write(stdout,"(/A,A)") "Save average phQ infortmation to the file: ",trim(phQ_file_) - + + write(stdout,"(/A,A)") "Save average phQ infortmation to the file: ",trim(phQ_file_) + end subroutine save_phQ - + subroutine read_phQ(inode,icore,nmodes,nq,nsnap,phQsit) integer,intent(in) :: inode,icore,nmodes,nq,nsnap complex(kind=dpc),intent(inout) :: phQsit(nmodes,nq,0:nsnap) integer :: phq_unit - character(len=maxlen) :: phQ_file_ - character(len=maxlen) :: ctmp1,ctmp2 - complex(kind=dpc),allocatable :: phQsit_(:,:,:) - - if(.not. allocated(phQsit_)) allocate(phQsit_(nmodes,nq,0:nsnap)) - - write(ctmp1,*) inode - write(ctmp2,*) icore - phQ_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(phQ_file)) + character(len=maxlen) :: phQ_file_ + character(len=maxlen) :: ctmp1,ctmp2 + complex(kind=dpc),allocatable :: phQsit_(:,:,:) + + if(.not. allocated(phQsit_)) allocate(phQsit_(nmodes,nq,0:nsnap)) + + write(ctmp1,*) inode + write(ctmp2,*) icore + phQ_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(phQ_file)) phq_unit = io_file_unit() call open_file(phQ_file_,phq_unit) - - do i=1,4 + + do i=1,4 read(phq_unit,*) - enddo + enddo do isnap=0,nsnap read(phq_unit,"(13X,*(2(1X,E12.5)))") ((phQsit_(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo - phQsit =phQsit + phQsit_ - + phQsit =phQsit + phQsit_ + call close_file(phQ_file_,phq_unit) - end subroutine read_phQ - + end subroutine read_phQ + subroutine plot_phQ(nmodes,nq,nsnap,phQsit) integer,intent(in) :: nmodes,nq,nsnap complex(kind=dpc),intent(in) :: phQsit(nmodes,nq,0:nsnap) integer :: phq_unit - character(len=maxlen) :: phQ_file_ - phQ_file_ = trim(outdir)//trim(adjustl(phQ_file)) + character(len=maxlen) :: phQ_file_ + phQ_file_ = trim(outdir)//trim(adjustl(phQ_file)) phq_unit = io_file_unit() call open_file(phQ_file_,phq_unit) @@ -120,32 +120,32 @@ subroutine plot_phQ(nmodes,nq,nsnap,phQsit) enddo enddo - write(phq_unit,"(A)") "Average of Normal mode coordinate for one core sample. ((phQ(imode,iq),imode=1,nmodes),iq=1,nq)" - write(phq_unit,"(A12,*(2(1X,A12)))") "time",(("Re[Q"//trim(adjustl(cphmode(imode,iq)))//"]",& + write(phq_unit,"(A)") "Average of Normal mode coordinate for one core sample. ((phQ(imode,iq),imode=1,nmodes),iq=1,nq)" + write(phq_unit,"(A12,*(2(1X,A12)))") "time",(("Re[Q"//trim(adjustl(cphmode(imode,iq)))//"]",& "Im[Q"//trim(adjustl(cphmode(imode,iq)))//"]",imode=1,nmodes),iq=1,nq) write(phq_unit,"(A12,*(2(1X,A12)))") "fs ",((" a.u. "," a.u. ",imode=1,nmodes),iq=1,nq) - write(phq_unit,"(A12,*(2(1X,F12.5)))") "Omega(meV)",(((wf_sub(imode,iq)*ryd2meV,wf_sub(imode,iq)*ryd2meV),imode=1,nmodes),iq=1,nq) + write(phq_unit,"(A12,*(2(1X,F12.5)))") "Omega(meV)",(((wf(imode,iq)*ryd2meV,wf(imode,iq)*ryd2meV),imode=1,nmodes),iq=1,nq) do isnap=0,nsnap write(phq_unit,"(F12.5,*(2(1X,E12.5)))") dt*nstep*isnap*ry_to_fs,((phQsit(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo call close_file(phQ_file_,phq_unit) - - write(stdout,"(A,A)") "Write average phQ infortmation to the file: ",trim(phQ_file_) + + write(stdout,"(A,A)") "Write average phQ infortmation to the file: ",trim(phQ_file_) end subroutine plot_phQ - + - !Write average phP infortmation to the file: phPsit.dat.gnu + !Write average phP infortmation to the file: phPsit.dat.gnu subroutine save_phP(nmodes,nq,nsnap,phPsit) integer,intent(in) :: nmodes,nq,nsnap complex(kind=dp),intent(in) :: phPsit(nmodes,nq,0:nsnap) integer :: php_unit - character(len=maxlen) :: phP_file_ + character(len=maxlen) :: phP_file_ php_unit = io_file_unit() - phP_file_ = trim(outdir)//trim(adjustl(phP_file)) + phP_file_ = trim(outdir)//trim(adjustl(phP_file)) call open_file(phP_file_,php_unit) if(.not. allocated(cphmode)) allocate(cphmode(nmodes,nq)) @@ -158,20 +158,20 @@ subroutine save_phP(nmodes,nq,nsnap,phPsit) enddo - write(php_unit,"(A)") "Average of Normal mode verlocity for one core sample.((phP(imode,iq),imode=1,nmodes),iq=1,nq)" - write(php_unit,"(A12,*(2(1X,A12)))") "time",(("Re[P"//trim(adjustl(cphmode(imode,iq)))//"]",& + write(php_unit,"(A)") "Average of Normal mode verlocity for one core sample.((phP(imode,iq),imode=1,nmodes),iq=1,nq)" + write(php_unit,"(A12,*(2(1X,A12)))") "time",(("Re[P"//trim(adjustl(cphmode(imode,iq)))//"]",& "Im[P"//trim(adjustl(cphmode(imode,iq)))//"]",imode=1,nmodes),iq=1,nq) - !write(php_unit,"(A12,*(2(1X,A12)))") "time",(("Re[phP]","Im[phP]",imode=1,nmodes),iq=1,nq) + !write(php_unit,"(A12,*(2(1X,A12)))") "time",(("Re[phP]","Im[phP]",imode=1,nmodes),iq=1,nq) write(php_unit,"(A12,*(2(1X,A12)))") "fs ",((" a.u. "," a.u. ",imode=1,nmodes),iq=1,nq) - write(php_unit,"(A12,*(2(1X,F12.5)))") "Omega(meV)",(((wf_sub(imode,iq)*ryd2meV,wf_sub(imode,iq)*ryd2meV),imode=1,nmodes),iq=1,nq) + write(php_unit,"(A12,*(2(1X,F12.5)))") "Omega(meV)",(((wf(imode,iq)*ryd2meV,wf(imode,iq)*ryd2meV),imode=1,nmodes),iq=1,nq) do isnap=0,nsnap write(php_unit,"(F12.5,*(2(1X,E12.5)))") dt*nstep*isnap*ry_to_fs,((phPsit(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo call close_file(phP_file_,php_unit) - write(stdout,"(A,A)") "Save average phP infortmation to the file: ",trim(phP_file_) - + write(stdout,"(A,A)") "Save average phP infortmation to the file: ",trim(phP_file_) + end subroutine save_phP subroutine read_phP(inode,icore,nmodes,nq,nsnap,phPsit) @@ -179,19 +179,19 @@ subroutine read_phP(inode,icore,nmodes,nq,nsnap,phPsit) complex(kind=dpc),intent(inout) :: phPsit(nmodes,nq,0:nsnap) integer :: php_unit - character(len=maxlen) :: phP_file_ - character(len=maxlen) :: ctmp1,ctmp2 - complex(kind=dpc),allocatable :: phPsit_(:,:,:) - - if(.not. allocated(phPsit_)) allocate(phPsit_(nmodes,nq,0:nsnap)) - - write(ctmp1,*) inode - write(ctmp2,*) icore - phP_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(phP_file)) + character(len=maxlen) :: phP_file_ + character(len=maxlen) :: ctmp1,ctmp2 + complex(kind=dpc),allocatable :: phPsit_(:,:,:) + + if(.not. allocated(phPsit_)) allocate(phPsit_(nmodes,nq,0:nsnap)) + + write(ctmp1,*) inode + write(ctmp2,*) icore + phP_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(phP_file)) php_unit = io_file_unit() call open_file(phP_file_,php_unit) - + do i=1,4 read(php_unit,*) enddo @@ -199,19 +199,19 @@ subroutine read_phP(inode,icore,nmodes,nq,nsnap,phPsit) read(php_unit,"(13X,*(2(1X,E12.5)))") ((phPsit_(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo - phPsit =phPsit + phPsit_ - + phPsit =phPsit + phPsit_ + call close_file(phP_file_,php_unit) - end subroutine read_phP + end subroutine read_phP subroutine plot_phP(nmodes,nq,nsnap,phPsit) integer,intent(in) :: nmodes,nq,nsnap complex(kind=dpc),intent(in) :: phPsit(nmodes,nq,0:nsnap) integer :: php_unit - character(len=maxlen) :: phP_file_ - phP_file_ = trim(outdir)//trim(adjustl(phP_file)) + character(len=maxlen) :: phP_file_ + phP_file_ = trim(outdir)//trim(adjustl(phP_file)) php_unit = io_file_unit() call open_file(phP_file_,php_unit) @@ -224,33 +224,33 @@ subroutine plot_phP(nmodes,nq,nsnap,phPsit) enddo enddo - write(php_unit,"(A)") "Average of Normal mode verlocity for one core sample.((phP(imode,iq),imode=1,nmodes),iq=1,nq)" - write(php_unit,"(A12,*(2(1X,A12)))") "time",(("Re[P"//trim(adjustl(cphmode(imode,iq)))//"]",& - "Im[P"//trim(adjustl(cphmode(imode,iq)))//"]",imode=1,nmodes),iq=1,nq) + write(php_unit,"(A)") "Average of Normal mode verlocity for one core sample.((phP(imode,iq),imode=1,nmodes),iq=1,nq)" + write(php_unit,"(A12,*(2(1X,A12)))") "time",(("Re[P"//trim(adjustl(cphmode(imode,iq)))//"]",& + "Im[P"//trim(adjustl(cphmode(imode,iq)))//"]",imode=1,nmodes),iq=1,nq) !write(php_unit,"(A12,*(2(1X,A12)))") "time",(("Re[phP]","Im[phP]",imode=1,nmodes),iq=1,nq) write(php_unit,"(A12,*(2(1X,A12)))") "fs ",((" a.u. "," a.u. ",imode=1,nmodes),iq=1,nq) - write(php_unit,"(A12,*(2(1X,F12.5)))") "Omega(meV)",(((wf_sub(imode,iq)*ryd2meV,wf_sub(imode,iq)*ryd2meV),imode=1,nmodes),iq=1,nq) + write(php_unit,"(A12,*(2(1X,F12.5)))") "Omega(meV)",(((wf(imode,iq)*ryd2meV,wf(imode,iq)*ryd2meV),imode=1,nmodes),iq=1,nq) do isnap=0,nsnap write(php_unit,"(F12.5,*(2(1X,E12.5)))") dt*nstep*isnap*ry_to_fs,((phPsit(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo call close_file(phP_file_,php_unit) - write(stdout,"(A,A)") "Write average phP infortmation to the file: ",trim(phP_file_) - - end subroutine plot_phP + write(stdout,"(A,A)") "Write average phP infortmation to the file: ",trim(phP_file_) + + end subroutine plot_phP + - - !Write average phK infortmation to the file: phKsit.dat.gnu + !Write average phK infortmation to the file: phKsit.dat.gnu subroutine save_phK(nmodes,nq,nsnap,phKsit) integer,intent(in) :: nmodes,nq,nsnap real(kind=dp),intent(in) :: phKsit(nmodes,nq,0:nsnap) integer :: phK_unit - character(len=maxlen) :: phK_file_ + character(len=maxlen) :: phK_file_ phK_unit = io_file_unit() - phK_file_ = trim(outdir)//trim(adjustl(phK_file)) + phK_file_ = trim(outdir)//trim(adjustl(phK_file)) call open_file(phK_file_,phK_unit) if(.not. allocated(cphmode)) allocate(cphmode(nmodes,nq)) @@ -262,20 +262,20 @@ subroutine save_phK(nmodes,nq,nsnap,phKsit) enddo enddo - write(phK_unit,"(A)") "Average of Normal mode kinetic energy for one core sample.SUM_phK,((phK(imode,iq),imode=1,nmodes),iq=1,nq)" - write(phK_unit,"(*(1X,A12))") "time","SUM_phK",(("K"//trim(adjustl(cphmode(imode,iq))),imode=1,nmodes),iq=1,nq) + write(phK_unit,"(A)") "Average of Normal mode kinetic energy for one core sample.SUM_phK,((phK(imode,iq),imode=1,nmodes),iq=1,nq)" + write(phK_unit,"(*(1X,A12))") "time","SUM_phK",(("K"//trim(adjustl(cphmode(imode,iq))),imode=1,nmodes),iq=1,nq) !write(phK_unit,"(*(1X,A12))") "time ","SUM_phK",(("phK(mode,q)",imode=1,nmodes),iq=1,nq) - write(phK_unit,"(*(1X,A12))") "fs "," meV ",((" meV ",imode=1,nmodes),iq=1,nq) - write(phK_unit,"(2(1X,A12),*(1X,F12.5))") "Omega(meV)","SUM_phK",((wf_sub(imode,iq)*ryd2meV,imode=1,nmodes),iq=1,nq) + write(phK_unit,"(*(1X,A12))") "fs "," meV ",((" meV ",imode=1,nmodes),iq=1,nq) + write(phK_unit,"(2(1X,A12),*(1X,F12.5))") "Omega(meV)","SUM_phK",((wf(imode,iq)*ryd2meV,imode=1,nmodes),iq=1,nq) do isnap=0,nsnap write(phK_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,SUM(phKsit(:,:,isnap))*ryd2meV,& - ((phKsit(imode,iq,isnap)*ryd2meV,imode=1,nmodes),iq=1,nq) + ((phKsit(imode,iq,isnap)*ryd2meV,imode=1,nmodes),iq=1,nq) enddo call close_file(phK_file_,phK_unit) - write(stdout,"(A,A)") "Save average phK infortmation to the file: ",trim(phK_file_) - + write(stdout,"(A,A)") "Save average phK infortmation to the file: ",trim(phK_file_) + end subroutine save_phK subroutine read_phK(inode,icore,nmodes,nq,nsnap,phKsit) @@ -283,33 +283,33 @@ subroutine read_phK(inode,icore,nmodes,nq,nsnap,phKsit) real(kind=dp),intent(inout) :: phKsit(nmodes,nq,0:nsnap) integer :: phk_unit - character(len=maxlen) :: phK_file_ - character(len=maxlen) :: ctmp1,ctmp2 - real(kind=dp),allocatable :: phKsit_(:,:,:) - - if(.not. allocated(phKsit_)) allocate(phKsit_(nmodes,nq,0:nsnap)) - phKsit_ = 0.0 + character(len=maxlen) :: phK_file_ + character(len=maxlen) :: ctmp1,ctmp2 + real(kind=dp),allocatable :: phKsit_(:,:,:) + + if(.not. allocated(phKsit_)) allocate(phKsit_(nmodes,nq,0:nsnap)) + phKsit_ = 0.0 - write(ctmp1,*) inode - write(ctmp2,*) icore - phK_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(phK_file)) + write(ctmp1,*) inode + write(ctmp2,*) icore + phK_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(phK_file)) phk_unit = io_file_unit() call open_file(phK_file_,phk_unit) - + do i=1,4 read(phK_unit,*) - enddo + enddo do isnap=0,nsnap read(phk_unit,"(26X,*(1X,E12.5))") ((phKsit_(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo - phKsit =phKsit + phKsit_ - + phKsit =phKsit + phKsit_ + call close_file(phK_file_,phk_unit) - end subroutine read_phK - + end subroutine read_phK + subroutine plot_phK(nmodes,nq,nsnap,phKsit) use dynamics,only : bolziman implicit none @@ -321,8 +321,8 @@ subroutine plot_phK(nmodes,nq,nsnap,phKsit) real(kind=dp) :: womiga integer :: phK_unit - character(len=maxlen) :: phK_file_ - + character(len=maxlen) :: phK_file_ + phK_file_ = trim(outdir)//trim(adjustl(phK_file)) phK_unit = io_file_unit() call open_file(phK_file_,phK_unit) @@ -338,13 +338,13 @@ subroutine plot_phK(nmodes,nq,nsnap,phKsit) allocate(phKsit_ave(1:nmodes,1:nq,0:nsnap)) phKsit_ave = 0.0 - - if(.not. allocated(nqv)) allocate(nqv(nmodes,nq)) + + allocate(nqv(nmodes,nq)) nqv = 0.0 do iq=1,nq do imode=1,nmodes - womiga = wf_sub(imode,iq) + womiga = wf(imode,iq) if(womiga >= eps_acustic) then nqv(imode,iq) = bolziman(womiga,temp)+0.5 endif @@ -353,7 +353,7 @@ subroutine plot_phK(nmodes,nq,nsnap,phKsit) if(l_ph_quantum) then do isnap =0,nsnap - phKsit_ave(:,:,isnap) = phKsit(:,:,isnap) - 0.5*nqv*wf_sub*ryd2meV + phKsit_ave(:,:,isnap) = phKsit(:,:,isnap) - 0.5*nqv*wf*ryd2meV enddo else phKsit_ave = phKsit - 0.5*temp*K_B_Ryd*ryd2meV @@ -361,35 +361,35 @@ subroutine plot_phK(nmodes,nq,nsnap,phKsit) if(l_ph_quantum) then write(phK_unit,"(A)") "Average of Normal mode kinetic energy - 0.5*nqv*hbar*wqv for all trajecotry.SUM_phK,((phK(imode,iq),imode=1,nmodes),iq=1,nq)" - else + else write(phK_unit,"(A)") "Average of Normal mode kinetic energy - 0.5*KB*T for all trajecotry.SUM_phK,((phK(imode,iq),imode=1,nmodes),iq=1,nq)" endif write(phK_unit,"(*(1X,A12))") "time ","phK(mode,q)",(("K"//trim(adjustl(cphmode(imode,iq))),imode=1,nmodes),iq=1,nq) - write(phK_unit,"(*(1X,A12))") "fs "," meV ",((" meV ",imode=1,nmodes),iq=1,nq) - write(phK_unit,"(2(1X,A12),*(1X,F12.5))") "Omega(meV)","SUM_phK",((wf_sub(imode,iq)*ryd2meV,imode=1,nmodes),iq=1,nq) + write(phK_unit,"(*(1X,A12))") "fs "," meV ",((" meV ",imode=1,nmodes),iq=1,nq) + write(phK_unit,"(2(1X,A12),*(1X,F12.5))") "Omega(meV)","SUM_phK",((wf(imode,iq)*ryd2meV,imode=1,nmodes),iq=1,nq) do isnap=0,nsnap write(phK_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,SUM(phKsit(:,:,isnap)),& - ((phKsit(imode,iq,isnap),imode=1,nmodes),iq=1,nq) + ((phKsit(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo call close_file(phK_file_,phK_unit) deallocate(phKsit_ave) - write(stdout,"(A,A)") "Write average phK infortmation to the file: ",trim(phK_file_) - + write(stdout,"(A,A)") "Write average phK infortmation to the file: ",trim(phK_file_) + end subroutine plot_phK - !Write average phU infortmation to the file: phUsit.dat.gnu + !Write average phU infortmation to the file: phUsit.dat.gnu subroutine save_phU(nmodes,nq,nsnap,phUsit) integer,intent(in) :: nmodes,nq,nsnap real(kind=dp),intent(in) :: phUsit(nmodes,nq,0:nsnap) integer :: phU_unit - character(len=maxlen) :: phU_file_ + character(len=maxlen) :: phU_file_ phU_unit = io_file_unit() - phU_file_ = trim(outdir)//trim(adjustl(phU_file)) + phU_file_ = trim(outdir)//trim(adjustl(phU_file)) call open_file(phU_file_,phU_unit) if(.not. allocated(cphmode)) allocate(cphmode(nmodes,nq)) @@ -400,20 +400,20 @@ subroutine save_phU(nmodes,nq,nsnap,phUsit) cphmode(imode,iq) = "("//trim(adjustl(ctmp1))//","//trim(adjustl(ctmp2))//")" enddo enddo - - write(phU_unit,"(A)") "Average of Normal mode potential energy for one core sample.SUM_phU,((phU(imode,iq),imode=1,nmodes),iq=1,nq)" - write(phU_unit,"(*(1X,A12))") "time ","SUM_phU",(("U"//trim(adjustl(cphmode(imode,iq))),imode=1,nmodes),iq=1,nq) - write(phU_unit,"(*(1X,A12))") "fs "," meV ",((" meV ",imode=1,nmodes),iq=1,nq) - write(phU_unit,"(2(1X,A12),*(1X,F12.5))") "Omega(meV)","SUM_phU",((wf_sub(imode,iq)*ryd2meV,imode=1,nmodes),iq=1,nq) + + write(phU_unit,"(A)") "Average of Normal mode potential energy for one core sample.SUM_phU,((phU(imode,iq),imode=1,nmodes),iq=1,nq)" + write(phU_unit,"(*(1X,A12))") "time ","SUM_phU",(("U"//trim(adjustl(cphmode(imode,iq))),imode=1,nmodes),iq=1,nq) + write(phU_unit,"(*(1X,A12))") "fs "," meV ",((" meV ",imode=1,nmodes),iq=1,nq) + write(phU_unit,"(2(1X,A12),*(1X,F12.5))") "Omega(meV)","SUM_phU",((wf(imode,iq)*ryd2meV,imode=1,nmodes),iq=1,nq) do isnap=0,nsnap write(phU_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,SUM(phUsit(:,:,isnap))*ryd2meV,& - ((phUsit(imode,iq,isnap)*ryd2meV,imode=1,nmodes),iq=1,nq) + ((phUsit(imode,iq,isnap)*ryd2meV,imode=1,nmodes),iq=1,nq) enddo call close_file(phU_file_,phU_unit) - write(stdout,"(A,A)") "Save average phU infortmation to the file: ",trim(phU_file_) - + write(stdout,"(A,A)") "Save average phU infortmation to the file: ",trim(phU_file_) + end subroutine save_phU subroutine read_phU(inode,icore,nmodes,nq,nsnap,phUsit) @@ -421,39 +421,39 @@ subroutine read_phU(inode,icore,nmodes,nq,nsnap,phUsit) real(kind=dp),intent(inout) :: phUsit(nmodes,nq,0:nsnap) integer :: phu_unit - character(len=maxlen) :: phU_file_ - character(len=maxlen) :: ctmp1,ctmp2 - real(kind=dp),allocatable :: phUsit_(:,:,:) - - if(.not. allocated(phUsit_)) allocate(phUsit_(nmodes,nq,0:nsnap)) - phUsit_ = 0.0 - - write(ctmp1,*) inode - write(ctmp2,*) icore - phU_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(phU_file)) + character(len=maxlen) :: phU_file_ + character(len=maxlen) :: ctmp1,ctmp2 + real(kind=dp),allocatable :: phUsit_(:,:,:) + + if(.not. allocated(phUsit_)) allocate(phUsit_(nmodes,nq,0:nsnap)) + phUsit_ = 0.0 + + write(ctmp1,*) inode + write(ctmp2,*) icore + phU_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(phU_file)) phu_unit = io_file_unit() call open_file(phU_file_,phu_unit) - + do i=1,4 read(phu_unit,*) - enddo + enddo do isnap=0,nsnap read(phu_unit,"(26X,*(1X,E12.5))") ((phUsit_(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo - phUsit =phUsit + phUsit_ - + phUsit =phUsit + phUsit_ + call close_file(phU_file_,phu_unit) - end subroutine read_phU + end subroutine read_phU subroutine plot_phU(nmodes,nq,nsnap,phUsit) integer,intent(in) :: nmodes,nq,nsnap real(kind=dp),intent(in) :: phUsit(nmodes,nq,0:nsnap) integer :: phU_unit - character(len=maxlen) :: phU_file_ - real(kind=dp),allocatable :: phUsit_ave(:,:,:) + character(len=maxlen) :: phU_file_ + real(kind=dp),allocatable :: phUsit_ave(:,:,:) @@ -475,7 +475,7 @@ subroutine plot_phU(nmodes,nq,nsnap,phUsit) if(l_ph_quantum) then do isnap =0,nsnap - phUsit_ave(:,:,isnap) = phUsit(:,:,isnap) - 0.5*nqv*wf_sub*ryd2meV + phUsit_ave(:,:,isnap) = phUsit(:,:,isnap) - 0.5*nqv*wf*ryd2meV enddo else phUsit_ave = phUsit_ave - 0.5*temp*K_B_Ryd*ryd2meV @@ -483,29 +483,28 @@ subroutine plot_phU(nmodes,nq,nsnap,phUsit) if(l_ph_quantum) then write(phU_unit,"(A)") "Average of Normal mode potential energy - 0.5*nqv*hbar*wqv for all trajecotry.SUM_phU,((phU(imode,iq),imode=1,nmodes),iq=1,nq)" - else + else write(phU_unit,"(A)") "Average of Normal mode potential energy - 0.5*KB*T for all trajecotry.SUM_phU,((phU(imode,iq),imode=1,nmodes),iq=1,nq)" endif write(phU_unit,"(*(1X,A12))") "time ","phU(mode,q)",(("U"//trim(adjustl(cphmode(imode,iq))),imode=1,nmodes),iq=1,nq) - write(phU_unit,"(*(1X,A12))") "fs "," meV ",((" meV ",imode=1,nmodes),iq=1,nq) - write(phU_unit,"(2(1X,A12),*(1X,F12.5))") "Omega(meV)","SUM_phU",((wf_sub(imode,iq)*ryd2meV,imode=1,nmodes),iq=1,nq) + write(phU_unit,"(*(1X,A12))") "fs "," meV ",((" meV ",imode=1,nmodes),iq=1,nq) + write(phU_unit,"(2(1X,A12),*(1X,F12.5))") "Omega(meV)","SUM_phU",((wf(imode,iq)*ryd2meV,imode=1,nmodes),iq=1,nq) do isnap=0,nsnap write(phU_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,SUM(phUsit_ave(:,:,isnap)),& - ((phUsit_ave(imode,iq,isnap),imode=1,nmodes),iq=1,nq) + ((phUsit_ave(imode,iq,isnap),imode=1,nmodes),iq=1,nq) enddo call close_file(phU_file_,phU_unit) deallocate(phUsit_ave) - write(stdout,"(A,A)") "Write average phU infortmation to the file: ",trim(phU_file_) - + write(stdout,"(A,A)") "Write average phU infortmation to the file: ",trim(phU_file_) + end subroutine plot_phU subroutine plot_ph_temp(nmodes,nq,nsnap,phKsit,phUsit) - use elph2,only : xqf implicit none integer , intent(in) :: nmodes,nq,nsnap real(kind=dp),intent(in) :: phKsit(nmodes,nq,0:nsnap),phUsit(nmodes,nq,0:nsnap) @@ -533,8 +532,8 @@ subroutine plot_ph_temp(nmodes,nq,nsnap,phKsit,phUsit) if(l_ph_quantum) then do isnap=0,nsnap do iq=1,nq - if(wf_sub(imode,iq)>0.0) then - phEsit(iq,isnap) = (phKsit(imode,iq,isnap)+phUsit(imode,iq,isnap))/(wf_sub(imode,iq)*ryd2meV) + if(wf(imode,iq)>0.0) then + phEsit(iq,isnap) = (phKsit(imode,iq,isnap)+phUsit(imode,iq,isnap))/(wf(imode,iq)*ryd2meV) else phEsit(iq,isnap) = 0.0 endif @@ -546,7 +545,7 @@ subroutine plot_ph_temp(nmodes,nq,nsnap,phKsit,phUsit) if(phEsit(iq,isnap)<=0.5) then phTemp(iq,isnap) = 0.0 else - phTemp(iq,isnap) = (wf_sub(nmodes,iq)/K_B_Ryd)/log((phEsit(iq,isnap)+0.5)/(phEsit(iq,isnap)-0.5)) + phTemp(iq,isnap) = (wf(nmodes,iq)/K_B_Ryd)/log((phEsit(iq,isnap)+0.5)/(phEsit(iq,isnap)-0.5)) endif enddo enddo @@ -565,9 +564,9 @@ subroutine plot_ph_temp(nmodes,nq,nsnap,phKsit,phUsit) write(phT_unit,"(6(1X,A12))") " fs ","iq","b1","b2","b3"," K " do isnap=0,nsnap do iq=1,nq - qx=xqf(1,indexq(iq)) - qy=xqf(2,indexq(iq)) - qz=xqf(3,indexq(iq)) + qx=xqf(1,iq) + qy=xqf(2,iq) + qz=xqf(3,iq) if(qx>0.5) qx=qx-1.0 if(qy>0.5) qy=qy-1.0 if(qz>0.5) qz=qz-1.0 @@ -589,17 +588,17 @@ end subroutine plot_ph_temp - !The electron(hole) population on different adiabatic wave function. + !The electron(hole) population on different adiabatic wave function. subroutine save_wsit(nfre,nsnap,naver,wsit,wsit_file) implicit none integer,intent(in) :: nfre,nsnap,naver real(kind=dp),intent(in) :: wsit(nfre,0:nsnap) character(len=*),intent(in) :: wsit_file - character(len=maxlen) :: wsit_file_ - integer :: wsit_unit + character(len=maxlen) :: wsit_file_ + integer :: wsit_unit wsit_unit = io_file_unit() - wsit_file_ = trim(outdir)//trim(adjustl(wsit_file)) + wsit_file_ = trim(outdir)//trim(adjustl(wsit_file)) call open_file(wsit_file_,wsit_unit) if(.not. allocated(cefree)) allocate(cefree(nfre)) @@ -608,18 +607,18 @@ subroutine save_wsit(nfre,nsnap,naver,wsit,wsit_file) cefree(ifre) = "("//trim(adjustl(ctmp1))//")" enddo - write(wsit_unit,"(A)") "The average electron(hole) population on different adiabatic wave function for one core sample. " - write(wsit_unit,"(*(1X,A12))") "time ",("wsit"//trim(adjustl(cefree(ifre))),ifre=1,nfre) + write(wsit_unit,"(A)") "The average electron(hole) population on different adiabatic wave function for one core sample. " + write(wsit_unit,"(*(1X,A12))") "time ",("wsit"//trim(adjustl(cefree(ifre))),ifre=1,nfre) write(wsit_unit,"((1X,A12),*(1X,I12))") "fs ",(ifre,ifre=1,nfre) - do isnap=0,nsnap + do isnap=0,nsnap write(wsit_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(wsit(ifre,isnap),ifre=1,nfre) - enddo + enddo call close_file(wsit_file_,wsit_unit) - write(stdout,"(A,A)")"Save the average electron(hole) wavefction population on & - &different adiabatic wave function to file:",trim(wsit_file_) + write(stdout,"(A,A)")"Save the average electron(hole) wavefction population on & + &different adiabatic wave function to file:",trim(wsit_file_) deallocate(cefree) @@ -629,43 +628,43 @@ subroutine read_wsit(inode,icore,nfre,nsnap,naver,wsit,wsit_file) implicit none integer,intent(in) :: nfre,nsnap,naver,inode,icore real(kind=dp),intent(inout) :: wsit(nfre,0:nsnap) - character(len=*),intent(in) :: wsit_file - + character(len=*),intent(in) :: wsit_file + integer :: wsit_unit - character(len=maxlen) :: wsit_file_ - character(len=maxlen) :: ctmp1,ctmp2 - - real(kind=dp),allocatable :: wsit_(:,:) - if(.not. allocated(wsit_)) allocate(wsit_(nfre,0:nsnap)) - - write(ctmp1,*) inode - write(ctmp2,*) icore - wsit_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(wsit_file)) + character(len=maxlen) :: wsit_file_ + character(len=maxlen) :: ctmp1,ctmp2 + + real(kind=dp),allocatable :: wsit_(:,:) + if(.not. allocated(wsit_)) allocate(wsit_(nfre,0:nsnap)) + + write(ctmp1,*) inode + write(ctmp2,*) icore + wsit_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(wsit_file)) wsit_unit = io_file_unit() call open_file(wsit_file_,wsit_unit) - + read(wsit_unit,*) read(wsit_unit,*) - read(wsit_unit,*) - do isnap=0,nsnap - read(wsit_unit,"(13X,*(1X,E12.5))") (wsit_(ifre,isnap),ifre=1,nfre) + read(wsit_unit,*) + do isnap=0,nsnap + read(wsit_unit,"(13X,*(1X,E12.5))") (wsit_(ifre,isnap),ifre=1,nfre) enddo - wsit = wsit + wsit_ - + wsit = wsit + wsit_ + call close_file(wsit_file,wsit_unit) - end subroutine read_wsit + end subroutine read_wsit - subroutine plot_wsit(nfre,nfre_sh,nsnap,naver,wsit,wsit_file) + subroutine plot_wsit(nfre,nsnap,naver,wsit,wsit_file) implicit none - integer,intent(in) :: nfre,nfre_sh,nsnap,naver + integer,intent(in) :: nfre,nsnap,naver real(kind=dp),intent(in) :: wsit(nfre,0:nsnap) character(len=*),intent(in) :: wsit_file integer :: wsit_unit - character(len=maxlen) :: wsit_file_ - wsit_file_ = trim(outdir)//trim(adjustl(wsit_file)) + character(len=maxlen) :: wsit_file_ + wsit_file_ = trim(outdir)//trim(adjustl(wsit_file)) wsit_unit = io_file_unit() call open_file(wsit_file_,wsit_unit) @@ -675,33 +674,33 @@ subroutine plot_wsit(nfre,nfre_sh,nsnap,naver,wsit,wsit_file) cefree(ifre) = "("//trim(adjustl(ctmp1))//")" enddo - write(wsit_unit,"(A)") "The average electron(hole) population on different adiabatic wave function for all trajecotry. " - write(wsit_unit,"(*(1X,A12))") "time ",("wsit"//trim(adjustl(cefree(ifre))),ifre=1,nfre_sh) - write(wsit_unit,"((1X,A12),*(1X,A12))") "fs ",("|w|2",ifre=1,nfre_sh) - do isnap=0,nsnap - write(wsit_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(wsit(ifre,isnap),ifre=1,nfre_sh) + write(wsit_unit,"(A)") "The average electron(hole) population on different adiabatic wave function for all trajecotry. " + write(wsit_unit,"(*(1X,A12))") "time ",("wsit"//trim(adjustl(cefree(ifre))),ifre=1,nfre) + write(wsit_unit,"((1X,A12),*(1X,A12))") "fs ",("|w|2",ifre=1,nfre) + do isnap=0,nsnap + write(wsit_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(wsit(ifre,isnap),ifre=1,nfre) enddo call close_file(wsit_file_,wsit_unit) - write(stdout,"(A,A)")"Write the average electron(hole) wavefction population on & - &different adiabatic wave function to file:",trim(wsit_file_) + write(stdout,"(A,A)")"Write the average electron(hole) wavefction population on & + &different adiabatic wave function to file:",trim(wsit_file_) deallocate(cefree) end subroutine plot_wsit - ! save the first trajecotry active PES and first trajecotry PES : pes_e.dat_f.gnu - ! save the average active PES and PES for all trajecotry. : pes_e.dat_average.gnu + ! save the first trajecotry active PES and first trajecotry PES : pes_e.dat_f.gnu + ! save the average active PES and PES for all trajecotry. : pes_e.dat_average.gnu subroutine save_pes(nfre,nsnap,naver,pes_one,pes,pes_file) implicit none integer , intent(in) :: nfre,nsnap,naver real(kind=dp),intent(in) :: pes_one(0:nfre,0:nsnap),pes(0:nfre,0:nsnap) character(len=*),intent(in) :: pes_file - character(len=maxlen) :: pes_file_ - integer :: pes_unit - pes_file_= trim(outdir)//trim(adjustl(pes_file)) + character(len=maxlen) :: pes_file_ + integer :: pes_unit + pes_file_= trim(outdir)//trim(adjustl(pes_file)) pes_unit = io_file_unit() call open_file(pes_file_,pes_unit) @@ -712,30 +711,30 @@ subroutine save_pes(nfre,nsnap,naver,pes_one,pes,pes_file) enddo write(pes_unit,"(A,I8,A)") "Averager of naver=",naver, " pes(ifre,isnap),ifre=0,nfre)" - write(pes_unit,"(*(1X,A12))") "time ","active_pes",("pes"//trim(adjustl(cefree(ifre))),ifre=1,nfre) + write(pes_unit,"(*(1X,A12))") "time ","active_pes",("pes"//trim(adjustl(cefree(ifre))),ifre=1,nfre) write(pes_unit,"(*(1X,A12))") "fs",("eV",ifre=0,nfre) - do isnap=0,nsnap + do isnap=0,nsnap write(pes_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(pes(ifre,isnap)*RYTOEV,ifre=0,nfre) enddo - write(stdout,"(A,A)") "Save average of active PES and PES for all trajecotry to the file:",trim(pes_file_) + write(stdout,"(A,A)") "Save average of active PES and PES for all trajecotry to the file:",trim(pes_file_) - call close_file(pes_file_,pes_unit) + call close_file(pes_file_,pes_unit) - pes_file_= trim(outdir)//trim(adjustl(pes_file))//".first.dat" + pes_file_= trim(outdir)//trim(adjustl(pes_file))//".first.dat" pes_unit = io_file_unit() call open_file(pes_file_,pes_unit) write(pes_unit,"(A,I8,A)") "One trajecotry PES of iaver=",1, " pes_one(ifre,isnap),ifre=0,nfre)" - write(pes_unit,"(*(1X,A12))") "time ","active_pes",("pesone(ifre)",ifre=1,nfre) + write(pes_unit,"(*(1X,A12))") "time ","active_pes",("pesone(ifre)",ifre=1,nfre) write(pes_unit,"(*(1X,A12))") "fs",("eV",ifre=0,nfre) - do isnap=0,nsnap + do isnap=0,nsnap write(pes_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(pes_one(ifre,isnap)*RYTOEV,ifre=0,nfre) enddo - call close_file(pes_file_,pes_unit) - write(stdout,"(A,A)") "Save the first trajecotry active PES and PES to the file:",trim(pes_file_) + call close_file(pes_file_,pes_unit) + write(stdout,"(A,A)") "Save the first trajecotry active PES and PES to the file:",trim(pes_file_) deallocate(cefree) @@ -749,56 +748,56 @@ subroutine read_pes(inode,icore,nfre,nsnap,naver,pes_one,pes,pes_file) integer :: pes_unit - character(len=maxlen) :: pes_file_ - character(len=maxlen) :: ctmp1,ctmp2 - - real(kind=dp),allocatable :: pes_(:,:) - if(.not. allocated(pes_)) allocate(pes_(0:nfre,0:nsnap)) - - write(ctmp1,*) inode - write(ctmp2,*) icore - pes_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(pes_file)) + character(len=maxlen) :: pes_file_ + character(len=maxlen) :: ctmp1,ctmp2 + + real(kind=dp),allocatable :: pes_(:,:) + if(.not. allocated(pes_)) allocate(pes_(0:nfre,0:nsnap)) + + write(ctmp1,*) inode + write(ctmp2,*) icore + pes_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(pes_file)) pes_unit = io_file_unit() call open_file(pes_file_,pes_unit) - - read(pes_unit,*) - read(pes_unit,*) + read(pes_unit,*) + read(pes_unit,*) + read(pes_unit,*) do isnap=0,nsnap - read(pes_unit,"(13X,*(1X,E12.5))") (pes_(ifre,isnap),ifre=0,nfre) + read(pes_unit,"(13X,*(1X,E12.5))") (pes_(ifre,isnap),ifre=0,nfre) enddo - - pes = pes + pes_ - - call close_file(pes_file_,pes_unit) - - pes_file_ = trim(adjustl(pes_file_))//".first.dat" + + pes = pes + pes_ + + call close_file(pes_file_,pes_unit) + + pes_file_ = trim(adjustl(pes_file_))//".first.dat" pes_unit = io_file_unit() call open_file(pes_file_,pes_unit) - - read(pes_unit,*) - read(pes_unit,*) + read(pes_unit,*) + read(pes_unit,*) + read(pes_unit,*) do isnap=0,nsnap - read(pes_unit,"(13X,*(1X,E12.5))") (pes_(ifre,isnap),ifre=0,nfre) + read(pes_unit,"(13X,*(1X,E12.5))") (pes_(ifre,isnap),ifre=0,nfre) enddo - - pes_one = pes_ - - call close_file(pes_file_,pes_unit) - + + pes_one = pes_ + + call close_file(pes_file_,pes_unit) + end subroutine read_pes - - subroutine plot_pes(nfre,nfre_sh,nsnap,pes_one,pes,wsit,pes_file) + + subroutine plot_pes(nfre,nsnap,pes_one,pes,wsit,pes_file) implicit none - integer , intent(in) :: nfre,nfre_sh,nsnap + integer , intent(in) :: nfre,nsnap real(kind=dp),intent(in) :: pes_one(0:nfre,0:nsnap),pes(0:nfre,0:nsnap),wsit(1:nfre,0:nsnap) character(len=*),intent(in) :: pes_file - - character(len=maxlen) :: pes_file_ + + character(len=maxlen) :: pes_file_ integer :: pes_unit - pes_file_ = trim(outdir)//trim(adjustl(pes_file))//"_f.dat" + pes_file_ = trim(outdir)//trim(adjustl(pes_file))//"_f.dat" if(.not. allocated(cefree)) allocate(cefree(nfre)) do ifre=1,nfre @@ -808,59 +807,59 @@ subroutine plot_pes(nfre,nfre_sh,nsnap,pes_one,pes,wsit,pes_file) pes_unit = io_file_unit() call open_file(pes_file_,pes_unit) - write(pes_unit,"(A)") "USing the first trajecotry as an example to Plotting potential Energy Surface(PES),and the active PES." - write(pes_unit,"(*(1X,A12))") "time ","active_fpes",("pes"//trim(adjustl(cefree(ifre))),ifre=1,nfre_sh) - write(pes_unit,"(*(1X,A12))") "fs "," eV " ,(" eV ",ifre=1,nfre_sh) + write(pes_unit,"(A)") "USing the first trajecotry as an example to Plotting potential Energy Surface(PES),and the active PES." + write(pes_unit,"(*(1X,A12))") "time ","active_fpes",("pes"//trim(adjustl(cefree(ifre))),ifre=1,nfre) + write(pes_unit,"(*(1X,A12))") "fs "," eV " ,(" eV ",ifre=1,nfre) do isnap=0,nsnap - write(pes_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(pes_one(ifre,isnap),ifre=0,nfre_sh) + write(pes_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(pes_one(ifre,isnap),ifre=0,nfre) enddo - + call close_file(pes_file_,pes_unit) - - write(stdout,"(A,A)") "Write the first trajecotry active PES and PES to the file:",trim(pes_file_) + + write(stdout,"(A,A)") "Write the first trajecotry active PES and PES to the file:",trim(pes_file_) - pes_file_ = trim(outdir)//trim(adjustl(pes_file))//"_average.dat" + pes_file_ = trim(outdir)//trim(adjustl(pes_file))//"_average.dat" pes_unit = io_file_unit() call open_file(pes_file_,pes_unit) - write(pes_unit,"(A)") "Plotting averager of potential Energy Surface(PES),and the averager active PES and their weight." + write(pes_unit,"(A)") "Plotting averager of potential Energy Surface(PES),and the averager active PES and their weight." write(pes_unit,"(*(1X,A12))") "time ","ave_PES","wsit ","avePES(a)","avePES(1)","avePES(nfre)" - write(pes_unit,"(*(1X,A12))") " fs "," eV ","w(ifre)**2"," eV "," eV "," eV " - - do ifre =1 , nfre_sh - !write(pes_unit,"(A,I12,1X,A12)") "#isurface=",ifre,"wsit" - do isnap=0,nsnap - if(ifre ==1 ) then - write(pes_unit,"(6(1X,E12.5))") dt*nstep*isnap*ry_to_fs,pes(ifre,isnap),wsit(ifre,isnap),& - pes(0,isnap),pes(1,isnap),pes(nfre,isnap) - else - write(pes_unit,"(3(1X,E12.5))") dt*nstep*isnap*ry_to_fs,pes(ifre,isnap),wsit(ifre,isnap) - endif - enddo - write(pes_unit,*) + write(pes_unit,"(*(1X,A12))") " fs "," eV ","w(ifre)**2"," eV "," eV "," eV " + + do ifre =1 , nfre + !write(pes_unit,"(A,I12,1X,A12)") "#isurface=",ifre,"wsit" + do isnap=0,nsnap + if(ifre ==1 ) then + write(pes_unit,"(6(1X,E12.5))") dt*nstep*isnap*ry_to_fs,pes(ifre,isnap),wsit(ifre,isnap),& + pes(0,isnap),pes(1,isnap),pes(nfre,isnap) + else + write(pes_unit,"(3(1X,E12.5))") dt*nstep*isnap*ry_to_fs,pes(ifre,isnap),wsit(ifre,isnap) + endif + enddo + write(pes_unit,*) enddo - + call close_file(pes_file_,pes_unit) - - write(stdout,"(A,A)") "Write average of active PES and PES for all trajecotry to the file:",trim(pes_file_) - + + write(stdout,"(A,A)") "Write average of active PES and PES for all trajecotry to the file:",trim(pes_file_) + deallocate(cefree) end subroutine plot_pes - + - !The electron(hole) population on different diabatic wave function. + !The electron(hole) population on different diabatic wave function. subroutine save_csit(nfre,nsnap,naver,csit,csit_file) implicit none integer,intent(in) :: nfre,nsnap,naver real(kind=dp),intent(in) :: csit(nfre,0:nsnap) character(len=*),intent(in) :: csit_file - - character(len=maxlen) :: csit_file_ + + character(len=maxlen) :: csit_file_ integer :: csit_unit csit_unit = io_file_unit() - csit_file_ = trim(outdir)//trim(adjustl(csit_file)) + csit_file_ = trim(outdir)//trim(adjustl(csit_file)) call open_file(csit_file_,csit_unit) if(.not. allocated(cefree)) allocate(cefree(nfre)) @@ -868,21 +867,21 @@ subroutine save_csit(nfre,nsnap,naver,csit,csit_file) write(ctmp1,*) ifre cefree(ifre) = "("//trim(adjustl(ctmp1))//")" enddo - - write(csit_unit,"(A)") "The average electron(hole) population on different diabatic wave function for one core sample. " - write(csit_unit,"(*(1X,A12))") "time",("csit"//trim(adjustl(cefree(ifre))),ifre=1,nfre) + + write(csit_unit,"(A)") "The average electron(hole) population on different diabatic wave function for one core sample. " + write(csit_unit,"(*(1X,A12))") "time",("csit"//trim(adjustl(cefree(ifre))),ifre=1,nfre) write(csit_unit,"((1X,A12),*(1X,I12))") "fs ",(ifre,ifre=1,nfre) - do isnap=0,nsnap + do isnap=0,nsnap write(csit_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(csit(ifre,isnap),ifre=1,nfre) enddo call close_file(csit_file_,csit_unit) - write(stdout,"(A,A)")"Save the average electron(hole) wavefction population on & - &different diabatic wave function to file:",trim(csit_file_) + write(stdout,"(A,A)")"Save the average electron(hole) wavefction population on & + &different diabatic wave function to file:",trim(csit_file_) deallocate(cefree) end subroutine save_csit - + subroutine read_csit(inode,icore,nfre,nsnap,naver,csit,csit_file) implicit none integer,intent(in) :: nfre,nsnap,naver,inode,icore @@ -890,30 +889,30 @@ subroutine read_csit(inode,icore,nfre,nsnap,naver,csit,csit_file) character(len=*),intent(in) :: csit_file integer :: csit_unit - character(len=maxlen) :: csit_file_ - character(len=maxlen) :: ctmp1,ctmp2 - - real(kind=dp),allocatable :: csit_(:,:) - if(.not. allocated(csit_)) allocate(csit_(nfre,0:nsnap)) - - write(ctmp1,*) inode - write(ctmp2,*) icore - csit_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(csit_file)) + character(len=maxlen) :: csit_file_ + character(len=maxlen) :: ctmp1,ctmp2 + + real(kind=dp),allocatable :: csit_(:,:) + if(.not. allocated(csit_)) allocate(csit_(nfre,0:nsnap)) + + write(ctmp1,*) inode + write(ctmp2,*) icore + csit_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(csit_file)) csit_unit = io_file_unit() call open_file(csit_file_,csit_unit) read(csit_unit,*) read(csit_unit,*) - read(csit_unit,*) - do isnap=0,nsnap - read(csit_unit,"(13X,*(1X,E12.5))") (csit_(ifre,isnap),ifre=1,nfre) + read(csit_unit,*) + do isnap=0,nsnap + read(csit_unit,"(13X,*(1X,E12.5))") (csit_(ifre,isnap),ifre=1,nfre) enddo - csit = csit + csit_ - + csit = csit + csit_ + call close_file(csit_file,csit_unit) - end subroutine read_csit + end subroutine read_csit subroutine plot_csit(nfre,nsnap,naver,csit,csit_file) implicit none @@ -922,8 +921,8 @@ subroutine plot_csit(nfre,nsnap,naver,csit,csit_file) character(len=*),intent(in) :: csit_file integer :: csit_unit - character(len=maxlen) :: csit_file_ - csit_file_ = trim(outdir)//trim(adjustl(csit_file)) + character(len=maxlen) :: csit_file_ + csit_file_ = trim(outdir)//trim(adjustl(csit_file)) csit_unit = io_file_unit() call open_file(csit_file_,csit_unit) @@ -933,34 +932,34 @@ subroutine plot_csit(nfre,nsnap,naver,csit,csit_file) cefree(ifre) = "("//trim(adjustl(ctmp1))//")" enddo - write(csit_unit,"(A)") "The electron(hole) population on different diabatic wave function. " - write(csit_unit,"(*(1X,A12))") "time",("csit"//trim(adjustl(cefree(ifre))),ifre=1,nfre) + write(csit_unit,"(A)") "The electron(hole) population on different diabatic wave function. " + write(csit_unit,"(*(1X,A12))") "time",("csit"//trim(adjustl(cefree(ifre))),ifre=1,nfre) write(csit_unit,"((1X,A12),*(1X,A12))") "fs ",("|c|2",ifre=1,nfre) - do isnap=0,nsnap + do isnap=0,nsnap write(csit_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(csit(ifre,isnap),ifre=1,nfre) enddo call close_file(csit_file_,csit_unit) - - write(stdout,"(A,A)")"Write the average electron(hole) wavefction population on & - &different diabatic wave function to file:",trim(csit_file_) - + + write(stdout,"(A,A)")"Write the average electron(hole) wavefction population on & + &different diabatic wave function to file:",trim(csit_file_) + deallocate(cefree) end subroutine plot_csit - !The active adiabatic PES project to diabatic states. psit(ifre)=|P(ifre,isurface)|**2 + !The active adiabatic PES project to diabatic states. psit(ifre)=|P(ifre,isurface)|**2 subroutine save_psit(nfre,nsnap,naver,psit,psit_file) implicit none integer,intent(in) :: nfre,nsnap,naver real(kind=dp),intent(in) :: psit(nfre,0:nsnap) character(len=*),intent(in) :: psit_file - - character(len=maxlen) :: psit_file_ + + character(len=maxlen) :: psit_file_ integer :: psit_unit psit_unit = io_file_unit() - psit_file_ = trim(outdir)//trim(adjustl(psit_file)) + psit_file_ = trim(outdir)//trim(adjustl(psit_file)) call open_file(psit_file_,psit_unit) if(.not. allocated(cefree)) allocate(cefree(nfre)) @@ -969,21 +968,21 @@ subroutine save_psit(nfre,nsnap,naver,psit,psit_file) cefree(ifre) = "("//trim(adjustl(ctmp1))//")" enddo - write(psit_unit,"(A)") "The average active PES project on different diabatic wave function for one core sample. psit(ifre)=P(ifre,isurface)**2 " + write(psit_unit,"(A)") "The average active PES project on different diabatic wave function for one core sample. psit(ifre)=P(ifre,isurface)**2 " write(psit_unit,"(*(1X,A12))") "time",("psit"//trim(adjustl(cefree(ifre))),ifre=1,nfre) - write(psit_unit,"((1X,A12),*(1X,I12))") "fs ",(ifre,ifre=1,nfre) + write(psit_unit,"((1X,A12),*(1X,I12))") "fs ",(ifre,ifre=1,nfre) do isnap=0,nsnap write(psit_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(psit(ifre,isnap),ifre=1,nfre) enddo call close_file(psit_file_,psit_unit) - write(stdout,"(A,A)")"Save the average electron(hole) active PES project to & - &different diabatic wave function to file:",trim(psit_file_) + write(stdout,"(A,A)")"Save the average electron(hole) active PES project to & + &different diabatic wave function to file:",trim(psit_file_) deallocate(cefree) end subroutine save_psit - + subroutine read_psit(inode,icore,nfre,nsnap,naver,psit,psit_file) implicit none integer,intent(in) :: nfre,nsnap,naver,inode,icore @@ -991,31 +990,31 @@ subroutine read_psit(inode,icore,nfre,nsnap,naver,psit,psit_file) character(len=*),intent(in) :: psit_file integer :: psit_unit - character(len=maxlen) :: psit_file_ - character(len=maxlen) :: ctmp1,ctmp2 - - real(kind=dp),allocatable :: psit_(:,:) - if(.not. allocated(psit_)) allocate(psit_(nfre,0:nsnap)) - - write(ctmp1,*) inode - write(ctmp2,*) icore - psit_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(psit_file)) - + character(len=maxlen) :: psit_file_ + character(len=maxlen) :: ctmp1,ctmp2 + + real(kind=dp),allocatable :: psit_(:,:) + if(.not. allocated(psit_)) allocate(psit_(nfre,0:nsnap)) + + write(ctmp1,*) inode + write(ctmp2,*) icore + psit_file_ = "./node"//trim(adjustl(ctmp1))//"/sample"//trim(adjustl(ctmp2))//"/"//trim(outdir)//trim(adjustl(psit_file)) + psit_unit = io_file_unit() call open_file(psit_file_,psit_unit) read(psit_unit,*) - read(psit_unit,*) - read(psit_unit,*) - do isnap=0,nsnap - read(psit_unit,"(13X,*(1X,E12.5))") (psit_(ifre,isnap),ifre=1,nfre) + read(psit_unit,*) + read(psit_unit,*) + do isnap=0,nsnap + read(psit_unit,"(13X,*(1X,E12.5))") (psit_(ifre,isnap),ifre=1,nfre) enddo - psit = psit + psit_ - + psit = psit + psit_ + call close_file(psit_file,psit_unit) - end subroutine read_psit + end subroutine read_psit subroutine plot_psit(nfre,nsnap,naver,psit,psit_file) implicit none @@ -1024,8 +1023,8 @@ subroutine plot_psit(nfre,nsnap,naver,psit,psit_file) character(len=*),intent(in) :: psit_file integer :: psit_unit - character(len=maxlen) :: psit_file_ - psit_file_ = trim(outdir)//trim(adjustl(psit_file)) + character(len=maxlen) :: psit_file_ + psit_file_ = trim(outdir)//trim(adjustl(psit_file)) psit_unit = io_file_unit() call open_file(psit_file_,psit_unit) @@ -1034,36 +1033,36 @@ subroutine plot_psit(nfre,nsnap,naver,psit,psit_file) write(ctmp1,*) ifre cefree(ifre) = "("//trim(adjustl(ctmp1))//")" enddo - - write(psit_unit,"(A)") "The average active PES project on different diabatic wave function for all trajecotry. psit(ifre)=P(ifre,isurface)**2 " + + write(psit_unit,"(A)") "The average active PES project on different diabatic wave function for all trajecotry. psit(ifre)=P(ifre,isurface)**2 " write(psit_unit,"(*(1X,A12))") "time",("psit"//trim(adjustl(cefree(ifre))),ifre=1,nfre) - write(psit_unit,"((1X,A12),*(1X,A12))") "fs ",("|P|2",ifre=1,nfre) + write(psit_unit,"((1X,A12),*(1X,A12))") "fs ",("|P|2",ifre=1,nfre) do isnap=0,nsnap write(psit_unit,"(*(1X,E12.5))") dt*nstep*isnap*ry_to_fs,(psit(ifre,isnap),ifre=1,nfre) enddo call close_file(psit_file_,psit_unit) - write(stdout,"(A,A)")"Write the average electron(hole) active PES project to & - &different diabatic wave function to file:",trim(psit_file_) - + write(stdout,"(A,A)")"Write the average electron(hole) active PES project to & + &different diabatic wave function to file:",trim(psit_file_) + deallocate(cefree) - end subroutine plot_psit - - - subroutine plot_band_occupatin_withtime(nband,nk,Enk,nsnap,psit,csit,dsnap,band_file) - use elph2,only : xkf - implicit none - integer , intent(in) :: nband,nk,dsnap - real(kind=dp),intent(in) :: Enk(nband,nk) - integer , intent(in) :: nsnap - real(kind=dp),intent(in) :: psit(1:nband*nk,0:nsnap),csit(1:nband*nk,0:nsnap) - character(len=*),intent(in) :: band_file - - character(len=maxlen) :: band_file_ - integer :: band_unit - integer :: ipol,iband,ik,ifre - real(kind=dp) :: kx,ky,kz + end subroutine plot_psit + + + subroutine plot_band_occupatin_withtime(nband,nk,Enk,xk,nsnap,psit,csit,dsnap,band_file) + implicit none + integer , intent(in) :: nband,nk,dsnap + real(kind=dp),intent(in) :: Enk(nband,nk) + real(kind=dp),intent(in) :: xk(3,nk) + integer , intent(in) :: nsnap + real(kind=dp),intent(in) :: psit(1:nband*nk,0:nsnap),csit(1:nband*nk,0:nsnap) + character(len=*),intent(in) :: band_file + + character(len=maxlen) :: band_file_ + integer :: band_unit + integer :: ipol,iband,ik,ifre + real(kind=dp) :: kx,ky,kz do iband=1,nband write(ctmp1,*) iband @@ -1073,44 +1072,44 @@ subroutine plot_band_occupatin_withtime(nband,nk,Enk,nsnap,psit,csit,dsnap,band_ write(band_unit,"(A32,I12,A)") "Carrier occupation on the iband=",iband," structure at different time" write(band_unit,"(*(1X,A12))") "iband","ik","kx","ky","kz","E_nk",("psit","csit",isnap=0,nsnap,dsnap) write(band_unit,"(6(1X,A12),*(1X,F12.2))") "index_band"," index_k","b1","b2", "b3", " eV" ,& - (dt*nstep*isnap*ry_to_fs,dt*nstep*isnap*ry_to_fs,isnap=0,nsnap,dsnap) - do ik=1,nk - ifre = (ik-1)*nband+iband - kx = xkf(1,indexk(ik)) - ky = xkf(2,indexk(ik)) - kz = xkf(3,indexk(ik)) + (dt*nstep*isnap*ry_to_fs,dt*nstep*isnap*ry_to_fs,isnap=0,nsnap,dsnap) + do ik=1,nk + ifre = (ik-1)*nband+iband + kx = xk(1,ik) + ky = xk(2,ik) + kz = xk(3,ik) if(kx > 0.5) kx=kx-1.0 if(ky > 0.5) ky=ky-1.0 if(kz > 0.5) kz=kz-1.0 - write(band_unit,"(2(1X,I12),*(1X,F12.5))") iband,ik,kx,ky,kz,Enk(iband,ik)*RYTOEV,(psit(ifre,isnap),csit(ifre,isnap),isnap=0,nsnap,dsnap) + write(band_unit,"(2(1X,I12),*(1X,F12.5))") iband,ik,kx,ky,kz,Enk(iband,ik)*RYTOEV,(psit(ifre,isnap),csit(ifre,isnap),isnap=0,nsnap,dsnap) if(kx == 0.5 .or. ky==0.5 .or. kz==0.5) then if(kx==0.5) kx=kx-1.0 if(ky==0.5) ky=ky-1.0 if(kz==0.5) kz=kz-1.0 write(band_unit,"(2(1X,I12),*(1X,F12.5))") iband,ik,kx,ky,kz,Enk(iband,ik)*RYTOEV,(psit(ifre,isnap),csit(ifre,isnap),isnap=0,nsnap,dsnap) endif - enddo - + enddo + call close_file(band_file_,band_unit) - enddo + enddo write(stdout,"(A,A)") "Write Carrier occupation on the band structure at different time to file:",trim(band_file_) - - end subroutine plot_band_occupatin_withtime + + end subroutine plot_band_occupatin_withtime - - subroutine plot_eh_temp(nband,nk,Enk,nsnap,psit,csit,dsnap,temp_file) - use elph2,only : xkf + + subroutine plot_eh_temp(nband,nk,Enk,xk,nsnap,psit,csit,dsnap,temp_file) implicit none - integer , intent(in) :: nband,nk,dsnap - real(kind=dp),intent(in) :: Enk(nband,nk) - integer , intent(in) :: nsnap - real(kind=dp),intent(in) :: psit(1:nband*nk,0:nsnap),csit(1:nband*nk,0:nsnap) - character(len=*),intent(in) :: temp_file - - character(len=maxlen) :: temp_file_ - integer :: temp_unit - integer :: ipol,iband,ik,ifre - real(kind=dp) :: kx,ky,kz + integer , intent(in) :: nband,nk,dsnap + real(kind=dp),intent(in) :: Enk(nband,nk) + real(kind=dp),intent(in) :: xk(3,nk) + integer , intent(in) :: nsnap + real(kind=dp),intent(in) :: psit(1:nband*nk,0:nsnap),csit(1:nband*nk,0:nsnap) + character(len=*),intent(in) :: temp_file + + character(len=maxlen) :: temp_file_ + integer :: temp_unit + integer :: ipol,iband,ik,ifre + real(kind=dp) :: kx,ky,kz do iband=1,nband write(ctmp1,*) iband @@ -1121,13 +1120,13 @@ subroutine plot_eh_temp(nband,nk,Enk,nsnap,psit,csit,dsnap,temp_file) write(temp_unit,"(A34,I12,A)") "Carrier occupation on the iband = ",iband," structure at different time." write(temp_unit,"(9(A12,1X))") " time "," iband "," ik "," kx "," ky "," kz "," E_nk ","psit "," csit " write(temp_unit,"(9(A12,1X))") " fs "," iband "," ik "," b1 "," b2 ", " b3 "," eV "," psit "," csit " - - do isnap =0,nsnap + + do isnap =0,nsnap do ik=1,nk ifre = (ik-1)*nband+iband - kx = xkf(1,indexk(ik)) - ky = xkf(2,indexk(ik)) - kz = xkf(3,indexk(ik)) + kx = xk(1,ik) + ky = xk(2,ik) + kz = xk(3,ik) if(kx > 0.5) kx=kx-1.0 if(ky > 0.5) ky=ky-1.0 if(kz > 0.5) kz=kz-1.0 @@ -1144,7 +1143,7 @@ subroutine plot_eh_temp(nband,nk,Enk,nsnap,psit,csit,dsnap,temp_file) write(temp_unit,*) enddo call close_file(temp_file_,temp_unit) - enddo + enddo write(stdout,"(A,A)") "Write Carrier occupation on the band structure at different time to file:",trim(temp_file_) end subroutine plot_eh_temp diff --git a/src_complex/surfacecom.f90 b/src_complex/surfacecom.f90 index 2cb2781..5e23b6c 100644 --- a/src_complex/surfacecom.f90 +++ b/src_complex/surfacecom.f90 @@ -51,12 +51,6 @@ module surfacecom integer :: ieband_min,ieband_max,ihband_min,ihband_max integer :: nefre_sh,nhfre_sh - ! get a subspace of k-points to reduce the calcution resume. - integer ::nk_sub - integer,allocatable :: indexk(:) - integer :: nq_sub - integer,allocatable :: indexq(:) - !method of surface hopping character(len=8) :: MethodSH ! FSSH, SC-FSSH, CC-FSSH diff --git a/src_complex/surfacehopping.f90 b/src_complex/surfacehopping.f90 index bf3c96a..2cc03f7 100644 --- a/src_complex/surfacehopping.f90 +++ b/src_complex/surfacehopping.f90 @@ -1,15 +1,15 @@ module surfacehopping use kinds, only : dp,dpc use constants,only : cone,czero - use epwcom,only : nkf1,nkf2,nkf3,nqf1,nqf2,nqf3,kqmap,kqmap_sub - use elph2,only : wf,nbndfst,ibndmin,ibndmax + use epwcom,only : nkf1,nkf2,nkf3,nqf1,nqf2,nqf3,kqmap + 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,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,& - phQsit,phPsit,phKsit,phUsit,ld_gamma,nqv,& + phQsit,phPsit,phKsit,phUsit,ld_gamma,& dEa_dQ,dEa_dQ_e,dEa_dQ_h,dEa2_dQ2,dEa2_dQ2_e,dEa2_dQ2_h,& d_e,g_e,g1_e,c_e,c_e_nk,w_e,w0_e,d0_e,nefre_sh,& d_h,g_h,g1_h,c_h,c_h_nk,w_h,w0_h,d0_h,nhfre_sh,& @@ -26,19 +26,17 @@ module surfacehopping contains - subroutine allocatesh(methodsh,lelecsh,lholesh,nk_sub,nmodes,nq) + subroutine allocatesh(methodsh,lelecsh,lholesh,nmodes,nq) use io,only : msg,io_error implicit none character(len=*),intent(in) :: methodsh logical,intent(in):: lelecsh,lholesh - integer,intent(in):: nmodes,nq,nk_sub + integer,intent(in):: nmodes,nq integer :: ierr allocate(ld_gamma(nmodes,nq)) ld_gamma = 0.0 - allocate(nqv(nmodes,nq)) - nqv = 0.0 allocate(phQ(nmodes,nq),stat=ierr,errmsg=msg) ! x(1:nphfre) if(ierr /=0) then call errore('surfacehopping','Error allocating phQ',1) @@ -87,9 +85,9 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nk_sub,nmodes,nq) if(ierr /=0) call errore('surfacehopping','Error allocating E_e',1) allocate(P_e(nefre,nefre),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating P_e',1) - allocate(P_e_nk(neband,nk_sub,nefre),stat=ierr,errmsg=msg) + allocate(P_e_nk(neband,nktotf,nefre),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating P_e_nk',1) - allocate(P0_e_nk(neband,nk_sub,nefre),stat=ierr,errmsg=msg) + 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,nefre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ijk if(ierr /=0) call errore('surfacehopping','Error allocating d_e',1) @@ -98,9 +96,9 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nk_sub,nmodes,nq) allocate(dEa_dQ_e(nmodes,nq)) allocate(dEa2_dQ2_e(nmodes,nq)) - allocate(c_e_nk(neband,nk_sub),stat=ierr,errmsg=msg) + allocate(c_e_nk(neband,nktotf),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating c_e_nk',1) - allocate(c_e(neband*nk_sub),stat=ierr,errmsg=msg) + allocate(c_e(neband*nktotf),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating c_e',1) allocate(pes_one_e(0:nefre,0:nsnap),pes_e(0:nefre,0:nsnap)) @@ -164,9 +162,9 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nk_sub,nmodes,nq) if(ierr /=0) call errore('surfacehopping','Error allocating E_h',1) allocate(P_h(nhfre,nhfre),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating P_h',1) - allocate(P_h_nk(nhband,nk_sub,nhfre),stat=ierr,errmsg=msg) + allocate(P_h_nk(nhband,nktotf,nhfre),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating P_h_nk',1) - allocate(P0_h_nk(nhband,nk_sub,nhfre),stat=ierr,errmsg=msg) + 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,nhfre_sh,nmodes,nq),stat=ierr,errmsg=msg) !d_ijk if(ierr /=0) call errore('surfacehopping','Error allocating d_h',1) @@ -175,9 +173,9 @@ subroutine allocatesh(methodsh,lelecsh,lholesh,nk_sub,nmodes,nq) allocate(dEa_dQ_h(nmodes,nq)) allocate(dEa2_dQ2_h(nmodes,nq)) - allocate(c_h_nk(nhband,nk_sub),stat=ierr,errmsg=msg) + allocate(c_h_nk(nhband,nktotf),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating c_h_nk',1) - allocate(c_h(nhband*nk_sub),stat=ierr,errmsg=msg) + allocate(c_h(nhband*nktotf),stat=ierr,errmsg=msg) if(ierr /=0) call errore('surfacehopping','Error allocating c_h',1) allocate(pes_one_h(0:nhfre,0:nsnap),pes_h(0:nhfre,0:nsnap)) @@ -262,10 +260,10 @@ subroutine convert_diabatic_adiabatic(nfre,pp,cc,ww) sum_ww2 = REAL(SUM(ww*CONJG(ww))) !ww=0.0d0 - !do ik=1,nk_sub + !do ik=1,nktotf ! do iband=1,nbndfst ! iefre = (ik-1)*nbndfst + iband - ! do jk=1,nk_sub + ! do jk=1,nktotf ! do jband=1,nbndfst ! ww(iefre) = ww(iefre)+pp_nk(jband,jk,iefre)*c_nk(jband,jk) ! enddo @@ -349,18 +347,19 @@ 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,nfre_sh,dd) + subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,lit_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,pnk12 + complex(kind=dpc) :: epc !nfre = nband*nk nfre = nfre_sh @@ -368,8 +367,7 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfr dd=czero do iq=1,nq do ik =1 ,nk - ikq = kqmap_sub(ik,iq) - if(ikq /= 0) then + ikq = kqmap(ik,iq) do imode=1,nmodes do iband1=1,nband do iband2=1,nband @@ -378,17 +376,16 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfr 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 enddo enddo enddo - endif enddo enddo - + do ifre=1,nfre do jfre=1,nfre if(jfre/= ifre) then @@ -401,51 +398,6 @@ subroutine calculate_nonadiabatic_coupling(nmodes,nq,nband,nk,ee,p_nk,gmnvkq,nfr 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 diff --git a/src_complex/types.f90 b/src_complex/types.f90 index 56a0b1d..7c7f0bf 100644 --- a/src_complex/types.f90 +++ b/src_complex/types.f90 @@ -1,12 +1,15 @@ module types - use kinds,only :dp,dpc + use kinds,only :dp implicit none !declare type of gmnvkq_n0 with gmnvkq>0.0 type :: gmnvkq_n0 - integer :: iqv - integer :: ink - integer :: imkq - complex(kind=dpc) :: g + integer :: m + integer :: n + integer :: v + integer :: ik + integer :: iq + integer :: ikq + real(kind=dp) :: g end type gmnvkq_n0 end module types \ No newline at end of file diff --git a/src_complex/write_sh_information.f90 b/src_complex/write_sh_information.f90 index 0b5cc9f..78f0f74 100644 --- a/src_complex/write_sh_information.f90 +++ b/src_complex/write_sh_information.f90 @@ -63,7 +63,7 @@ subroutine write_initial_information(iaver,nmodes,nqtotf,wf,phQ,phP) else write(stdout,"(/,A)") " time(fs) rt(s) hsur& & E_h(eV) T_ph(eV) U_ph(eV) E_ph(eV) E_tot(eV)" - write(stdout,"(F9.2,F9.2,I5,5(1X,F11.4))") 0.00,0.00,& + write(stdout,"(F11.2,F11.2,I5,5(1X,F11.4))") 0.00,0.00,& ihsurface,-e_h(ihsurface)*ryd2eV,& SUM_phK*ryd2eV,SUM_phU*ryd2eV,SUM_phE*ryd2eV,(E_h(ihsurface)+SUM_phE)*ryd2eV endif @@ -71,7 +71,7 @@ subroutine write_initial_information(iaver,nmodes,nqtotf,wf,phQ,phP) if(lelecsh) then write(stdout,"(/,A)") " time(fs) rt(s) esur& & E_e(eV) T_ph(eV) U_ph(eV) E_ph(eV) E_tot(eV)" - write(stdout,"(F9.2,F9.2,I5,5(1X,F11.4))") 0.00,0.00,& + write(stdout,"(F11.2,F11.2,I5,5(1X,F11.4))") 0.00,0.00,& iesurface,e_e(iesurface)*ryd2eV,& SUM_phK*ryd2eV,SUM_phU*ryd2eV,SUM_phE*ryd2eV,(E_e(iesurface)+SUM_phE)*ryd2eV else