diff --git a/configure b/configure index 56674f6058..03b99c1f1e 100755 --- a/configure +++ b/configure @@ -1,6 +1,6 @@ #! /bin/sh # Guess values for system-dependent variables and create Makefiles. -# Generated by GNU Autoconf 2.71 for Yambo 5.2.0 r.23132 h.d2d5779d5. +# Generated by GNU Autoconf 2.71 for Yambo 5.2.0 r.23133 h.161205154. # # Report bugs to . # @@ -610,8 +610,8 @@ MAKEFLAGS= # Identity of this package. PACKAGE_NAME='Yambo' PACKAGE_TARNAME='yambo' -PACKAGE_VERSION='5.2.0 r.23132 h.d2d5779d5' -PACKAGE_STRING='Yambo 5.2.0 r.23132 h.d2d5779d5' +PACKAGE_VERSION='5.2.0 r.23133 h.161205154' +PACKAGE_STRING='Yambo 5.2.0 r.23133 h.161205154' PACKAGE_BUGREPORT='yambo@yambo-code.org' PACKAGE_URL='' @@ -1600,7 +1600,7 @@ if test "$ac_init_help" = "long"; then # Omit some internal or obsolete options to make the list less imposing. # This message is too long to be a string in the A/UX 3.1 sh. cat <<_ACEOF -\`configure' configures Yambo 5.2.0 r.23132 h.d2d5779d5 to adapt to many kinds of systems. +\`configure' configures Yambo 5.2.0 r.23133 h.161205154 to adapt to many kinds of systems. Usage: $0 [OPTION]... [VAR=VALUE]... @@ -1666,7 +1666,7 @@ fi if test -n "$ac_init_help"; then case $ac_init_help in - short | recursive ) echo "Configuration of Yambo 5.2.0 r.23132 h.d2d5779d5:";; + short | recursive ) echo "Configuration of Yambo 5.2.0 r.23133 h.161205154:";; esac cat <<\_ACEOF @@ -1876,7 +1876,7 @@ fi test -n "$ac_init_help" && exit $ac_status if $ac_init_version; then cat <<\_ACEOF -Yambo configure 5.2.0 r.23132 h.d2d5779d5 +Yambo configure 5.2.0 r.23133 h.161205154 generated by GNU Autoconf 2.71 Copyright (C) 2021 Free Software Foundation, Inc. @@ -2505,7 +2505,7 @@ cat >config.log <<_ACEOF This file contains any messages produced by compilers while running configure, to aid debugging if configure makes a mistake. -It was created by Yambo $as_me 5.2.0 r.23132 h.d2d5779d5, which was +It was created by Yambo $as_me 5.2.0 r.23133 h.161205154, which was generated by GNU Autoconf 2.71. Invocation command line was $ $0$ac_configure_args_raw @@ -3263,8 +3263,8 @@ ac_compiler_gnu=$ac_cv_c_compiler_gnu SVERSION="5" SSUBVERSION="2" SPATCHLEVEL="0" -SREVISION="23132" -SHASH="d2d5779d5" +SREVISION="23133" +SHASH="161205154" @@ -16578,7 +16578,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1 # report actual input values of CONFIG_FILES etc. instead of their # values after options handling. ac_log=" -This file was extended by Yambo $as_me 5.2.0 r.23132 h.d2d5779d5, which was +This file was extended by Yambo $as_me 5.2.0 r.23133 h.161205154, which was generated by GNU Autoconf 2.71. Invocation command line was CONFIG_FILES = $CONFIG_FILES @@ -16642,7 +16642,7 @@ ac_cs_config_escaped=`printf "%s\n" "$ac_cs_config" | sed "s/^ //; s/'/'\\\\\\\\ cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1 ac_cs_config='$ac_cs_config_escaped' ac_cs_version="\\ -Yambo config.status 5.2.0 r.23132 h.d2d5779d5 +Yambo config.status 5.2.0 r.23133 h.161205154 configured by $0, generated by GNU Autoconf 2.71, with options \\"\$ac_cs_config\\" diff --git a/include/version/version.m4 b/include/version/version.m4 index 6902b4d309..0c7f183362 100644 --- a/include/version/version.m4 +++ b/include/version/version.m4 @@ -1,9 +1,9 @@ -AC_INIT(Yambo, 5.2.0 r.23132 h.d2d5779d5, yambo@yambo-code.org) +AC_INIT(Yambo, 5.2.0 r.23133 h.161205154, yambo@yambo-code.org) SVERSION="5" SSUBVERSION="2" SPATCHLEVEL="0" -SREVISION="23132" -SHASH="d2d5779d5" +SREVISION="23133" +SHASH="161205154" AC_SUBST(SVERSION) AC_SUBST(SSUBVERSION) AC_SUBST(SPATCHLEVEL) diff --git a/src/interface/INIT_activate.F b/src/interface/INIT_activate.F index ea00cfed11..5b3023f184 100644 --- a/src/interface/INIT_activate.F +++ b/src/interface/INIT_activate.F @@ -532,7 +532,7 @@ subroutine INIT_activate() ! call initactivate(1,'RTBands Integrator SlowDynamics LinearRegime GrKind TwoAlph RTehEny ScattTresh') call initactivate(1,'RADLifeTime RADmagnific PhLifeTime DephCVonly DephEThresh DephTRange') - call initactivate(1,'RTstep NEsteps NETime DipoleEtresh SPINprojected') + call initactivate(1,'RTstep RKquestTresh NEsteps NETime DipoleEtresh SPINprojected') call initactivate(1,'IOtime IOCachetime') ! ! n_ext_fields is by default 0. It is set to the input number of fields when a command line option is provided (below in init.F) diff --git a/src/interface/INIT_load.F b/src/interface/INIT_load.F index c43c0405a3..ad586a16e6 100644 --- a/src/interface/INIT_load.F +++ b/src/interface/INIT_load.F @@ -68,12 +68,12 @@ subroutine INIT_load(defs,en,q,k,X,Xw,Dip) use stderr, ONLY:intc use RT_occupations,ONLY:RT_RTA_temp,RT_RTA_chem use RT_lifetimes, ONLY:RT_life_extrap_times_INPUT - use real_time, ONLY:RT_step, Integrator_name, RAD_LifeTime, Phase_LifeTime, & + use real_time, ONLY:RT_step,Integrator_name, RAD_LifeTime, Phase_LifeTime, & & NE_tot_time,RT_bands,RT_eh_en,RT_scatt_tresh,Gr_kind, Integrator_slow_approx, & & NE_MEM_treshold,NE_time_step_update_jump_INPUT, & & two_alpha,DbGd_EE_percent,NE_time_step_update_jump, & & NE_initial_time_step_update,NE_step_update_treshold, & -& RT_MAX_step,RAD_magnification, & +& RT_MAX_step,RAD_magnification,RK_quest_treshold, & & RT_deph_deg_thresh,RT_step_manual_prefactor,NE_time_step_update_manual use fields, ONLY:Efield,n_ext_fields_max,n_ext_fields use plasma, ONLY:PLASMA_redux_percent @@ -411,15 +411,16 @@ subroutine INIT_load(defs,en,q,k,X,Xw,Dip) call it(defs,'TwoAlpha', '[RT] C_nk ~ alpha*Gamma_nk^2 ',two_alpha,Verb_level=V_real_time) call it(defs,'GrKind', '[RT] G-ret kind: Lorentzian (QP)/ Hyperbolic QP_secant (HS)',Gr_kind,Verb_level=V_real_time) call it(defs,'RADLifeTime','[RT] Radiative life-time (if negative Yambo sets it equal to Phase_LifeTime in NL)', & -& RAD_LifeTime,unit=Time_unit(1)) +& RAD_LifeTime,unit=Time_unit(1),Verb_level=V_real_time) call it(defs,'RADmagnific','[RT] Radiative life-time magnification',RAD_magnification,Verb_level=V_real_time) - call it(defs,'PhLifeTime', '[RT] Constant Dephasing Time',Phase_LifeTime,unit=Time_unit(1)) + call it(defs,'PhLifeTime', '[RT] Constant Dephasing Time',Phase_LifeTime,unit=Time_unit(1),Verb_level=V_real_time) call it(defs,'DephEThresh', '[RT] Threshold on the energy difference between two states to dephase them',& - & RT_deph_deg_thresh,unit=E_unit,Verb_level=V_real_time) +& RT_deph_deg_thresh,unit=E_unit,Verb_level=V_real_time) ! ! Dynamics ! - call it(defs,'RTstep', '[RT] Real Time step length',RT_step,unit=Time_unit(3)) + call it(defs,'RTstep', '[RT] Real Time step: >0 absolute, <0 used as seed for the Runge-Kutta optimization procedure.',& +& RT_step,unit=Time_unit(3)) call it(defs,'NETime', '[RT] Simulation Time',NE_tot_time,unit=Time_unit(2)) call it(defs,'dTupdateTimeSet','[RT] Time for manual deltaT update',NE_time_step_update_manual,& & unit=Time_unit(1),verb_level=V_real_time) @@ -428,11 +429,12 @@ subroutine INIT_load(defs,en,q,k,X,Xw,Dip) call it(defs,'dTupdateJump','[RT] Time betweem two deltaT updates', NE_time_step_update_jump_INPUT,& & unit=Time_unit(1),verb_level=V_real_time) call it(defs,'dTupdateTresh','[RT][o/o] Treshold of deltaT updates', NE_step_update_treshold ,verb_level=V_real_time) + call it(defs,'RKquestTresh','[RT][o/oo] Treshold used to identify the best dT/integrator',RK_quest_treshold,verb_level=V_real_time) call it(defs,'dT_MAX', '[RT] Maximum value for the time-dependent dT', RT_MAX_step ,verb_level=V_real_time,unit=Time_unit(1)) call it(defs,'dT_SET', '[RT] Prefactor for manual dT update', RT_step_manual_prefactor,verb_level=V_real_time) - call it(defs,'Integrator', '[RT] Time-evolution (TE) integrator: "EULER", "RK2/RK4", "HEUN"',Integrator_name) + call it(defs,'Integrator', '[RT] Time-evolution RK integrator: "EULER/EMP/RK2/RK3/H3/RK42/RK41/best"',Integrator_name) call it(defs,'SlowDynamics','[RT] Slow-TE approximations: "EXPn" (n<=6), "EXP" (n=3), "INV+ACC", "EXP+ACC", "INV+DIAG"',& -& Integrator_slow_approx,verb_level=V_real_time) +& Integrator_slow_approx,verb_level=V_real_time) call it('f',defs,'LinearRegime','[RT] Adopt specific tricks in the EULER case assuming the linear-regime',verb_level=V_real_time) ! IO_times=(/CARR_RT_IO_t%INTERVAL_time_INPUT,Gless_RESTART_RT_IO_t%INTERVAL_time_INPUT,OUTPUT_RT_IO_t%INTERVAL_time_INPUT/) diff --git a/src/modules/SET_defaults.F b/src/modules/SET_defaults.F index dcb0a4dae7..c7255fcbe8 100644 --- a/src/modules/SET_defaults.F +++ b/src/modules/SET_defaults.F @@ -604,7 +604,7 @@ subroutine SET_defaults(INSTR,IND,OD,COM_DIR) NE_time_step_update_jump_INPUT = 0._SP NE_initial_time_step_update = 0._SP NE_step_update_treshold = 0.1_SP - NE_time_step_update_manual =-1.0_SP + NE_time_step_update_manual =-1.0_SP*FS2AUT ! CACHE_OBS_steps = 0 CACHE_OBS_steps_now = 0 @@ -624,7 +624,7 @@ subroutine SET_defaults(INSTR,IND,OD,COM_DIR) RAD_LifeTime = -1._SP*FS2AUT ! Automatic set equal to the dephasing #endif Phase_LifeTime = 0._SP - RT_deph_deg_thresh = 1.E-5_SP + RT_deph_deg_thresh = 1.E-5_SP/HA2EV RT_deph_cv_only =.FALSE. ! NE_MEM_treshold=0._SP diff --git a/src/modules/mod_real_time.F b/src/modules/mod_real_time.F index 256fe131b2..37ea86ef29 100644 --- a/src/modules/mod_real_time.F +++ b/src/modules/mod_real_time.F @@ -224,6 +224,7 @@ module real_time integer, parameter :: N_RK_integrators=7 type(RK_integrator_t) :: RK_integrators(N_RK_integrators) type(RK_integrator_t) :: RK_integrator + real(SP) :: RK_quest_treshold=1._SP ! ! ... Slowly varying dynamics logical :: l_slow_TE=.FALSE. ! Slow time-evolution => INV/EXPn/LINEAR_RESP approximations @@ -239,12 +240,6 @@ module real_time ! ... INV related logical :: l_RT_DIAG=.FALSE. ! Used in INV (RT_INV_step_diago) ! - ! - ! ... Coefficients in the implicit/explicit propagation procedure - real(SP) :: a_tableau(4) - real(SP) :: b_tableau(4) - real(SP) :: c_tableau(4) - ! ! ... names (referring to RKn) and approximation (Physics based) character(schlen) :: Integrator_name character(schlen) :: Integrator_slow_approx='no' @@ -368,6 +363,7 @@ subroutine RK_integrator_alloc(I,N,short,name) I%short=short I%name=name I%N=N + if (allocated(I%a)) return YAMBO_ALLOC(I%a,(I%N,I%N)) YAMBO_ALLOC(I%b,(I%N)) YAMBO_ALLOC(I%c,(I%N)) diff --git a/src/real_time_drivers/RT_driver.F b/src/real_time_drivers/RT_driver.F index d52de41c3f..a327dac337 100644 --- a/src/real_time_drivers/RT_driver.F +++ b/src/real_time_drivers/RT_driver.F @@ -131,7 +131,7 @@ subroutine RT_driver(E,X,k,q,Dip) ! call RT_apply(RT_bands,E,k,what="G",VERBOSE=.true.) ! - if (STRING_same(Integrator_name,"best")) call Runge_Kutta_optimization(E) + if (RT_step<0._SP.or.STRING_same(Integrator_name,"best")) call Runge_Kutta_step_and_order_estimation(E) ! call section('=','Initialization') !################################# diff --git a/src/real_time_initialize/.objects b/src/real_time_initialize/.objects index 947c084e82..644213530b 100644 --- a/src/real_time_initialize/.objects +++ b/src/real_time_initialize/.objects @@ -2,7 +2,7 @@ ELPH_objs = RT_ELPH_initialize.o #endif #if defined _RT -RT_head_objs = RT_initialize.o RT_start_and_restart.o Runge_Kutta_optimization.o +RT_head_objs = RT_initialize.o RT_start_and_restart.o Runge_Kutta_step_and_order_estimation.o RT_foot_objs = RT_Field_Commensurable_Frequencies.o RT_Dephasing_Matrix.o RT_G_lesser_init.o RT_occupations_update.o #endif objs = $(RT_head_objs) RT_occupations_and_levels_init.o $(RT_foot_objs) $(ELPH_objs) diff --git a/src/real_time_initialize/Runge_Kutta_optimization.F b/src/real_time_initialize/Runge_Kutta_optimization.F deleted file mode 100644 index 4f2a7f1409..0000000000 --- a/src/real_time_initialize/Runge_Kutta_optimization.F +++ /dev/null @@ -1,221 +0,0 @@ -! -! License-Identifier: GPL -! -! Copyright (C) 2024 the Yambo Team -! -! Authors (see AUTHORS file for details): AM -! -subroutine Runge_Kutta_optimization(E) - ! - use pars, ONLY:SP,cI - use real_time, ONLY:RT_bands,RT_step,NE_tot_time,NE_i_time,RT_dyn_step,& -& RK_integrators,N_RK_integrators,N_RK_integrators - use stderr, ONLY:STRING_same - use electrons, ONLY:levels - use R_lattice, ONLY:nXkibz - use fields, ONLY:gauge_field,Efield,n_ext_fields,Efield_strength - use units, ONLY:AUT2FS,SPEED_OF_LIGHT - use vec_operate, ONLY:normalize_v - use OUTPUT, ONLY:OUTPUT_driver - use LIVE_t, ONLY:live_timing - ! -#include - ! - type(levels) :: E - ! - ! Work Space - ! - type(gauge_field) :: A - real(SP) :: Damp,delta_E_max,delta_E - real(SP),allocatable:: Time(:),Ef(:) - integer :: N_steps,i_f,ik,ib,ibp - ! - ! .. iterative dT optimization - integer :: i_iter,N_iter_steps,it,itp,i_best_int,i_best_iter,i_int,m,n - integer,parameter :: N_max_iterations=20 - complex(SP),allocatable:: Y(:,:,:) - real(SP) :: TRESH,ERR(N_RK_integrators,N_max_iterations),dT_iter - complex(SP) :: Y_trial,M_trial - ! - call section('=','Better Runge-Kutta integrator estimation') - !############################################################ - ! - ! Initial step & damping - !------------------------ - RT_step=1./AUT2FS/1.E3_SP ! 1as - TRESH=10.E-3 - RT_dyn_step=RT_step - !Damp=100.0_SP/HA2EV/1000. ! 50 meV - ! - ! Fields init - !------------- - do i_f=1,n_ext_fields - if (STRING_same(Efield(i_f)%ef_name,'none')) cycle - Efield(i_f)%versor =normalize_v(Efield(i_f)%versor) - Efield(i_f)%amplitude =Efield_strength(Efield(i_f)%intensity,'AU') - Efield(i_f)%t_initial_indx=max(nint(Efield(i_f)%t_initial/RT_step)+1,2) - Efield(i_f)%t_initial =real(Efield(i_f)%t_initial_indx-1)*RT_step - enddo - ! - ! Steps & Allocs - !---------------- - N_steps=NE_tot_time/RT_step - YAMBO_ALLOC(Time,(N_steps+1)) - YAMBO_ALLOC(Ef,(N_steps+1)) - YAMBO_ALLOC(Y,(N_RK_integrators,N_max_iterations,N_steps+1)) - ! - ! Elemental frequencies - !======================= - ! - ! .. Electronic transitions - !---------------------------- - delta_E_max=-1. - do ik=1,nXkibz - do ib=RT_bands(1),RT_bands(2) - do ibp=RT_bands(1),RT_bands(2) - if (ib==ibp) cycle - delta_E=abs(E%E(ib,ik,1)-E%E(ibp,ik,1)) - delta_E_max=max(delta_E,delta_E_max) - enddo - enddo - enddo - ! - ! .. Efield - !----------- - NE_i_time=0 - do it=1,N_steps+1 - Time(it)=RT_step*(it-1) - NE_i_time=NE_i_time+1 - call RT_Ext_fields(A,Time(it),-1) - Ef(it)=dot_product(real(A%vecpot_vel),(/1.,1.,1./))/SPEED_OF_LIGHT/sqrt(3._SP) - enddo - ! - ! RK integrators - !---------------- - call Runge_Kutta_integrators(.FALSE.) - ! - ! Integrations - !-------------- - call live_timing('RK integrations',N_RK_integrators) - ! - ERR=10._SP*TRESH - ! - do i_int=1,N_RK_integrators - ! - ITER_LOOP: do i_iter=1,N_max_iterations - ! - Y(i_int,i_iter,1)=1._SP - ! - if (i_iter==1) dT_iter=RT_step - if (i_iter> 1) dT_iter=5._SP*real(i_iter-1,SP)*RT_step - ! - it=1 - do while (it=N_steps) cycle - Y(i_int,i_iter,it+itp)=Y(i_int,i_iter,it)+M_trial*itp*RT_step - enddo - it=it+N_iter_steps - enddo - ! - if (i_iter>1) then - call EVAL_the_diff(i_int,i_iter,ERR(i_int,i_iter)) - if (ERR(i_int,i_iter)>TRESH) exit ITER_LOOP - endif - ! - enddo ITER_LOOP - ! - write (*,*) i_int,ERR(i_int,:) - call live_timing(steps=1) - ! - enddo - call live_timing( ) - ! - do i_int=N_RK_integrators,1,-1 - do i_iter=2,N_max_iterations-1 - if (ERR(i_int,i_iter)TRESH) then - i_best_int=i_int - i_best_iter=i_int - endif - enddo - enddo - ! - write (*,*) trim(RK_integrators(i_best_int)%short) - write (*,*) 5._SP*real(i_best_iter-1,SP)*RT_step - stop - ! - call OUTPUT_driver("RK_check",action="open") - do it=1,N_steps+1 - call OUTPUT_driver("RK_check",TITLES=(/"Time"/),R_VALUES=(/Time(it)/),UNIT="fs",Kind="Time") - do i_int=1,N_RK_integrators - call OUTPUT_driver("RK_check",TITLES=(/trim(RK_integrators(i_int)%short)//" dT=1as"/),R_VALUES=(/aimag(Y(i_int,1,it))/)) - !call OUTPUT_driver("RK_check",TITLES=(/trim(RK_integrators(i_int)%short)//" dT=1as"/),R_VALUES=(/aimag(Y(i_int,1,it))/)) - enddo - call OUTPUT_driver("RK_check",action="write") - enddo - call OUTPUT_driver("RK_check",action="close") - ! - stop - ! - contains - ! - real(SP) function S(T) - ! - real(SP)::T - integer ::itp - ! - do itp=1,N_steps+1 - if (T>=Time(itp).and.TDIFF_max) DIFF_max=DIFF - if (abs(aimag(Y(I,1,it)))>VAL_max) VAL_max=abs(aimag(Y(I,1,it))) - enddo - ERR=DIFF_max/VAL_max - end subroutine - ! - subroutine SOLVE_it(I,F,it,Fp) - ! - real(SP) :: T_eval - integer :: I,it - complex(SP) :: K(4),Y_RK(4),F,Fp - ! - Y_RK=0._SP - K =0._SP - ! - do n=1,RK_integrators(I)%N - ! - Y_RK(n)=F - if (n>1) then - do m=1,n-1 - Y_RK(n)=Y_RK(n)+dT_iter*RK_integrators(I)%a(n,m)*K(m) - enddo - endif - ! - T_eval=Time(it)+RK_integrators(I)%c(n)*dT_iter - K(n) =-cI*S(T_eval)*Y_RK(n) - ! - enddo - ! - Fp=F - do n=1,RK_integrators(I)%N - Fp=Fp+RK_integrators(I)%b(n)*K(n)*dT_iter - enddo - ! - end subroutine - ! -end subroutine Runge_Kutta_optimization diff --git a/src/real_time_initialize/Runge_Kutta_step_and_order_estimation.F b/src/real_time_initialize/Runge_Kutta_step_and_order_estimation.F new file mode 100644 index 0000000000..c347891872 --- /dev/null +++ b/src/real_time_initialize/Runge_Kutta_step_and_order_estimation.F @@ -0,0 +1,292 @@ +! +! License-Identifier: GPL +! +! Copyright (C) 2024 the Yambo Team +! +! Authors (see AUTHORS file for details): AM +! +subroutine Runge_Kutta_step_and_order_estimation(E) + ! + use pars, ONLY:SP,cI,lchlen + use real_time, ONLY:RT_bands,RT_step,NE_tot_time,NE_i_time,RT_dyn_step,Integrator_name,& +& RK_integrators,N_RK_integrators,N_RK_integrators,RK_quest_treshold + use stderr, ONLY:STRING_same,real2ch,set_real_printed_length,intc + use electrons, ONLY:levels + use R_lattice, ONLY:nXkibz + use fields, ONLY:gauge_field,Efield,n_ext_fields,Efield_strength + use units, ONLY:AUT2FS,SPEED_OF_LIGHT + use vec_operate, ONLY:normalize_v + use OUTPUT, ONLY:OUTPUT_driver + use com, ONLY:msg + use parallel_m, ONLY:PP_indexes,myid,PP_indexes_reset + use parallel_int, ONLY:PP_redux_wait,PARALLEL_index + use LIVE_t, ONLY:live_timing + ! +#include + ! + type(levels) ::E + ! + ! Work Space + ! + type(gauge_field) ::A + real(SP) ::Damp,delta_E_max,delta_E + real(SP),allocatable::Time(:),Ef(:) + integer ::N_steps,i_f,ik,ib,ibp + ! + ! .. iterative dT optimization + integer,parameter ::N_max_iterations=20 + integer ::i_iter,N_iter_steps,it,itp,i_int,m,n,i_int_ref(1),dT_table(N_max_iterations),i_best_int + complex(SP),allocatable::Y(:,:,:) + real(SP) ::ERR(N_RK_integrators,N_max_iterations),dT_iter(N_max_iterations),& +& dT(N_RK_integrators),TRESH + complex(SP) ::Y_trial,M_trial + character(lchlen) ::string + type(PP_indexes) ::px + ! + call section('=','Runge-Kutta order and dT estimation') + !####################################################### + ! + ! Initial step & damping + !------------------------ + if (abs(RT_step)>0) then + RT_step=abs(RT_step) + else + RT_step=1./AUT2FS/1.E3_SP ! 1as + endif + RT_dyn_step=RT_step + TRESH=RK_quest_treshold/1000._SP + !Damp=100.0_SP/HA2EV/1000. ! 50 meV + ! + ! Fields init + !------------- + do i_f=1,n_ext_fields + if (STRING_same(Efield(i_f)%ef_name,'none')) cycle + Efield(i_f)%versor =normalize_v(Efield(i_f)%versor) + Efield(i_f)%amplitude =Efield_strength(Efield(i_f)%intensity,'AU') + Efield(i_f)%t_initial_indx=max(nint(Efield(i_f)%t_initial/RT_step)+1,2) + Efield(i_f)%t_initial =real(Efield(i_f)%t_initial_indx-1)*RT_step + enddo + ! + ! Steps & Allocs + !---------------- + N_steps=NE_tot_time/RT_step + YAMBO_ALLOC(Time,(N_steps+1)) + YAMBO_ALLOC(Ef,(N_steps+1)) + YAMBO_ALLOC(Y,(N_RK_integrators,N_max_iterations,N_steps+1)) + ! + ! Elemental frequencies + !======================= + ! + ! .. Electronic transitions + !---------------------------- + delta_E_max=-1. + do ik=1,nXkibz + do ib=RT_bands(1),RT_bands(2) + do ibp=RT_bands(1),RT_bands(2) + if (ib==ibp) cycle + delta_E=abs(E%E(ib,ik,1)-E%E(ibp,ik,1)) + delta_E_max=max(delta_E,delta_E_max) + enddo + enddo + enddo + ! + ! .. Efield + !----------- + NE_i_time=0 + do it=1,N_steps+1 + Time(it)=RT_step*(it-1) + NE_i_time=NE_i_time+1 + call RT_Ext_fields(A,Time(it),-1) + Ef(it)=dot_product(real(A%vecpot_vel),(/1.,1.,1./))/SPEED_OF_LIGHT/sqrt(3._SP) + enddo + ! + ! RK integrators + !---------------- + call Runge_Kutta_integrators(.FALSE.) + ! + ! Test time-steps + !----------------- + do i_iter=1,N_max_iterations + if (i_iter==1) dT_iter(i_iter)=RT_step + if (i_iter> 1) dT_iter(i_iter)=5._SP*real(i_iter-1,SP)*RT_step + enddo + ! + ! Integrations + !============== + ERR =0._SP + i_int_ref=0 + ! + ! Parallel distribution + !----------------------- + call PP_indexes_reset(px) + call PARALLEL_index(px,(/N_RK_integrators/)) + if (px%n_of_elements(myid+1)>0) call live_timing('dT optimization',px%n_of_elements(myid+1)) + ! + do i_int=1,N_RK_integrators + ! + if (.not.px%element_1D(i_int)) cycle + ! + ERR(i_int,2:)=10._SP*TRESH + ! + if (.not.STRING_same(Integrator_name,"best")) then + if (.not.STRING_same(Integrator_name,RK_integrators(i_int)%short)) then + call live_timing(steps=1) + cycle + endif + i_int_ref=i_int + endif + ! + ITER_LOOP: do i_iter=1,N_max_iterations + ! + Y(i_int,i_iter,1)=1._SP + ! + it=1 + do while (it=N_steps) cycle + Y(i_int,i_iter,it+itp)=Y(i_int,i_iter,it)+M_trial*itp*RT_step + enddo + it=it+N_iter_steps + enddo + ! + if (i_iter>1) then + call EVAL_the_diff(i_int,i_iter,ERR(i_int,i_iter)) + if (ERR(i_int,i_iter)>TRESH) exit ITER_LOOP + endif + ! + enddo ITER_LOOP + ! + call live_timing(steps=1) + ! + enddo + ! + call live_timing(steps=1) + ! + call PP_redux_wait(ERR) + call PP_redux_wait(i_int_ref) + if (i_int_ref(1)==0) i_int_ref=7 + ! + if (all(ERR(:,2:)>TRESH)) call warning("Impossible to reduce the time-step. Using user provided.") + ! + ! dT table + !---------- + dT=-1._SP + dT_table=0 + do i_int=1,N_RK_integrators + do i_iter=2,N_max_iterations-1 + if (ERR(i_int,i_iter)TRESH) then + dT(i_int)=dT_iter(i_iter) + dT_table(i_int)=i_iter + endif + enddo + enddo + ! + ! Output... + !----------- + ! + ! ... CHECK file + string="Runge_Kutta_optimization" + call OUTPUT_driver(trim(string),action="open") + n=nint(NE_tot_time/RT_step) + do it=1,N_steps+1 + call OUTPUT_driver(trim(string),TITLES=(/"Time"/),R_VALUES=(/Time(it)/),UNIT="fs",Kind="Time") + call set_real_printed_length(f_length=5) + call OUTPUT_driver(trim(string),TITLES=(/trim(RK_integrators(i_int_ref(1))%short)//" [dT=1as]"/),& +& R_VALUES=(/aimag(Y(i_int_ref(1),1,it))/)) + call set_real_printed_length( ) + do i_int=1,N_RK_integrators + if (dT_Table(i_int)==0) cycle + if (nint(NE_tot_time/dT(i_int))*RK_integrators(i_int)%N=Time(itp).and.TDIFF_max) DIFF_max=DIFF + if (abs(aimag(Y(I,1,it)))>VAL_max) VAL_max=abs(aimag(Y(I,1,it))) + enddo + ERR=DIFF_max/VAL_max + end subroutine + ! + subroutine SOLVE_it(I,it,F,Fp,dT_) + ! + real(SP) :: T_eval,dT_ + integer :: I,it + complex(SP) :: K(4),Y_RK(4),F,Fp + ! + Y_RK=0._SP + K =0._SP + ! + do n=1,RK_integrators(I)%N + ! + Y_RK(n)=F + if (n>1) then + do m=1,n-1 + Y_RK(n)=Y_RK(n)+dT_*RK_integrators(I)%a(n,m)*K(m) + enddo + endif + ! + T_eval=Time(it)+RK_integrators(I)%c(n)*dT_ + K(n) =-cI*S(T_eval)*Y_RK(n) + ! + enddo + ! + Fp=F + do n=1,RK_integrators(I)%N + Fp=Fp+RK_integrators(I)%b(n)*K(n)*dT_ + enddo + ! + end subroutine + ! +end subroutine Runge_Kutta_step_and_order_estimation diff --git a/src/real_time_propagation/.objects b/src/real_time_propagation/.objects index 77d8ff1455..58ca99aafd 100644 --- a/src/real_time_propagation/.objects +++ b/src/real_time_propagation/.objects @@ -7,6 +7,6 @@ PHEL_objects = RT_PHEL_Scattering_step.o RUNGE_KUTTA_objects = Runge_Kutta_integrators.o Runge_Kutta_initialize.o Runge_Kutta_driver.o Runge_Kutta_kernel.o #if defined _RT objs = RT_Scattering_step.o RT_propagate_fields.o RT_Ext_fields.o RT_occupations_eval.o \ - $(RUNGE_KUTTA_objects) \ + $(RUNGE_KUTTA_objects) RT_G_symmetrization.o \ RT_MEMORY_index.o RT_time_step_update.o RT_IO_type_time_steps.o $(PHEL_objects) #endif diff --git a/src/real_time_propagation/RT_G_symmetrization.F b/src/real_time_propagation/RT_G_symmetrization.F new file mode 100644 index 0000000000..175599765a --- /dev/null +++ b/src/real_time_propagation/RT_G_symmetrization.F @@ -0,0 +1,44 @@ +! +! License-Identifier: GPL +! +! Copyright (C) 2017 The Yambo Team +! +! Authors (see AUTHORS file for details): AM CA DS +! +subroutine RT_G_symmetrization(ik,G) + ! + use pars, ONLY:SP + use real_time, ONLY:RT_bands,G_lesser_reference + use electrons, ONLY:spin_occ + ! + implicit none + ! + integer, intent(in) :: ik + complex(SP), intent(inout) :: G(RT_bands(1):RT_bands(2),RT_bands(1):RT_bands(2)) + ! + integer :: ib,ibp + real(SP) :: E_occ_tmp,H_occ_tmp + complex(SP) :: tmp + ! + do ib=RT_bands(1),RT_bands(2) + ! + G(ib,ib)=cmplx(0._SP,aimag(G(ib,ib))) + ! + E_occ_tmp= aimag(G(ib,ib))+ aimag(G_lesser_reference(ib,ib,ik)) + H_occ_tmp=-aimag(G(ib,ib))+ ( spin_occ-aimag(G_lesser_reference(ib,ib,ik)) ) + ! + if ( E_occ_tmp<0._SP .or. H_occ_tmp>spin_occ ) G(ib,ib)= G_lesser_reference(ib,ib,ik) + if ( H_occ_tmp<0._SP .or. E_occ_tmp>spin_occ ) G(ib,ib)=cmplx(0._SP,spin_occ)-G_lesser_reference(ib,ib,ik) + ! + do ibp=ib+1,RT_bands(2) + ! + tmp=( G(ib,ibp)-conjg(G(ibp,ib)) )/2._SP + ! + G(ib,ibp)= tmp + G(ibp,ib)=-conjg(tmp) + ! + enddo + ! + enddo + ! +end subroutine diff --git a/src/real_time_propagation/RT_Glesser_evolve.F b/src/real_time_propagation/RT_Glesser_evolve.F deleted file mode 100644 index 39c28900d5..0000000000 --- a/src/real_time_propagation/RT_Glesser_evolve.F +++ /dev/null @@ -1,92 +0,0 @@ -! -! License-Identifier: GPL -! -! Copyright (C) 2017 The Yambo Team -! -! Authors (see AUTHORS file for details): AM CA DS -! -subroutine RT_Glesser_evolve(En,kpt,dG_old,dG_in,dG_out,dT) - ! - ! dG_old = dG(T0) --> Time at first iteration in integrator with more than one step - ! dG_in = dG(T ) --> Effective time at the present step, i.e. for RK2 at step 2, T=T0+dt/2 - ! dG_out = dG(T0+dt) - ! - ! T0+dt = T+dt' - ! - use pars, ONLY:SP,cZERO - use electrons, ONLY:levels,spin_occ - use R_lattice, ONLY:bz_samp - use parallel_m, ONLY:PAR_IND_Xk_ibz,PAR_G_k_range - use real_time, ONLY:G_lesser_reference,RT_bands,RT_nk,RTibz,l_slow_TE - ! - implicit none - ! - type(bz_samp), intent(in) :: kpt - type(levels), intent(in) :: En - ! - real(SP), intent(in) :: dT - complex(SP), intent(in) :: dG_old(RT_bands(1):RT_bands(2),RT_bands(1):RT_bands(2),PAR_G_k_range(1):PAR_G_k_range(2)) - complex(SP), intent(in) :: dG_in(RT_bands(1):RT_bands(2),RT_bands(1):RT_bands(2),PAR_G_k_range(1):PAR_G_k_range(2)) - complex(SP), intent(out) :: dG_out(RT_bands(1):RT_bands(2),RT_bands(1):RT_bands(2),PAR_G_k_range(1):PAR_G_k_range(2)) - ! - ! Work space - ! - integer :: ik - ! - dG_out=cZERO - ! - do ik=1,RT_nk - ! - if( .not.PAR_IND_Xk_ibz%element_1D(RTibz%k_map(ik)) ) cycle - ! - if ( l_slow_TE) call RT_slow_dynamics_driver(dG_old(:,:,ik),dG_in(:,:,ik),dG_out(:,:,ik),ik,dT) - if (.not.l_slow_TE) call RT_EULER_step(dG_in(:,:,ik),dG_out(:,:,ik),ik,En%nbf,dT) - ! - ! Carriers relaxation step: - ! - ! G(T0+dt) = G(T0) - i dt * Sigma_relax(T) - ! - call RT_Scattering_step(dG_old(:,:,ik),dG_in(:,:,ik),dG_out(:,:,ik),ik,dT) - ! - ! G(T+dt) = G(t)+dG(T) - dG_out(:,:,ik)=dG_out(:,:,ik)+dG_in(:,:,ik) - ! - call RT_G_symmetrization(dG_out(:,:,ik),ik) - ! - enddo - ! - contains - ! - subroutine RT_G_symmetrization(G_inout,ik) - ! - integer, intent(in) :: ik - complex(SP), intent(inout) :: G_inout(RT_bands(1):RT_bands(2),RT_bands(1):RT_bands(2)) - ! - integer :: ib,ibp - real(SP) :: E_occ_tmp,H_occ_tmp - complex(SP) :: tmp - ! - do ib=RT_bands(1),RT_bands(2) - ! - G_inout(ib,ib)=cmplx(0._SP,aimag(G_inout(ib,ib))) - ! - E_occ_tmp= aimag(G_inout(ib,ib))+ aimag(G_lesser_reference(ib,ib,ik)) - H_occ_tmp=-aimag(G_inout(ib,ib))+ ( spin_occ-aimag(G_lesser_reference(ib,ib,ik)) ) - ! - if ( E_occ_tmp<0._SP .or. H_occ_tmp>spin_occ ) G_inout(ib,ib)= G_lesser_reference(ib,ib,ik) - if ( H_occ_tmp<0._SP .or. E_occ_tmp>spin_occ ) G_inout(ib,ib)=cmplx(0._SP,spin_occ)-G_lesser_reference(ib,ib,ik) - ! - do ibp=ib+1,RT_bands(2) - ! - tmp=( G_inout(ib,ibp)-conjg(G_inout(ibp,ib)) )/2._SP - ! - G_inout(ib,ibp)= tmp - G_inout(ibp,ib)=-conjg(tmp) - ! - enddo - ! - enddo - ! - end subroutine - ! -end subroutine RT_Glesser_evolve diff --git a/src/real_time_propagation/Runge_Kutta_driver.F b/src/real_time_propagation/Runge_Kutta_driver.F index 7cf9384ad8..8d9c2c625d 100644 --- a/src/real_time_propagation/Runge_Kutta_driver.F +++ b/src/real_time_propagation/Runge_Kutta_driver.F @@ -75,7 +75,7 @@ subroutine Runge_Kutta_driver(G_tpdt,dG_tpdt,dG_t,A_tpdt,A_t,E,k,q) ! call timing('RT integrator',OPR='stop') ! - ! dG_n => New Hamiltonian@(t+c_n dt) Eq. 11b of Overleaf Notes ( in the SB) + ! dG_n => New Hamiltonian@(t+c_n dt) Eq. 11b of Overleaf Notes (in the SB) !--------------------------------------------------------------------------- ! ! ... E_ext(t+c*dt) (needed by RT_Hamiltonian) @@ -111,6 +111,7 @@ subroutine Runge_Kutta_driver(G_tpdt,dG_tpdt,dG_t,A_tpdt,A_t,E,k,q) eval_T=NE_time+RT_dyn_step ! do ik=1,RT_nk + ! if( .not.PAR_IND_Xk_ibz%element_1D(RTibz%k_map(ik)) ) cycle ! ! The following step(s) are in the RF @@ -121,10 +122,12 @@ subroutine Runge_Kutta_driver(G_tpdt,dG_tpdt,dG_t,A_tpdt,A_t,E,k,q) enddo ! ! Final rotation and sum in the SB - !------------------------------------------- + !---------------------------------- dG_tpdt(:,:,ik)=dG_tpdt(:,:,ik)*exp(-cI*H_RF(:,:,ik,1)*eval_T) G_tpdt(:,:,ik) =G_lesser_reference(:,:,ik)+dG_tpdt(:,:,ik) ! + call RT_G_symmetrization(ik,G_tpdt) + ! enddo ! ! ... E_ext(t+dt) diff --git a/src/real_time_propagation/Runge_Kutta_initialize.F b/src/real_time_propagation/Runge_Kutta_initialize.F index 850ec45db7..2923f4c6d8 100644 --- a/src/real_time_propagation/Runge_Kutta_initialize.F +++ b/src/real_time_propagation/Runge_Kutta_initialize.F @@ -79,12 +79,14 @@ subroutine Runge_Kutta_initialize() ! Rebuild Integrator Name(s) !=========================== Integrator_name=RK_integrators(i_active_integrator)%short - Integrator_slow_approx="" - if (l_slow_TE) then - if (l_RT_EXP ) Integrator_slow_approx = trim(Integrator_slow_approx)//" EXP"//trim(intc(Integrator_exp_order)) - if (l_RT_INV ) Integrator_slow_approx = trim(Integrator_slow_approx)//" INV" - if (l_RT_ACC ) Integrator_slow_approx = trim(Integrator_slow_approx)//"+ACCURATE" - if (l_RT_DIAG ) Integrator_slow_approx = trim(Integrator_slow_approx)//"+DIAGO" + if (.not.STRING_match(Integrator_slow_approx,"no")) then + Integrator_slow_approx="" + if (l_slow_TE) then + if (l_RT_EXP ) Integrator_slow_approx = trim(Integrator_slow_approx)//" EXP"//trim(intc(Integrator_exp_order)) + if (l_RT_INV ) Integrator_slow_approx = trim(Integrator_slow_approx)//" INV" + if (l_RT_ACC ) Integrator_slow_approx = trim(Integrator_slow_approx)//"+ACCURATE" + if (l_RT_DIAG ) Integrator_slow_approx = trim(Integrator_slow_approx)//"+DIAGO" + endif endif ! ! Copy the active integrator in the RT_integrator type