diff --git a/src/sh_integ.F90 b/src/sh_integ.F90 index c932ae49..da4141f1 100644 --- a/src/sh_integ.F90 +++ b/src/sh_integ.F90 @@ -10,6 +10,7 @@ module mod_sh_integ public :: sh_set_initialwf, sh_set_energy_shift public :: sh_select_integrator, sh_integrate_wf public :: sh_decoherence_correction, check_popsum, sh_TFS_transmat + public :: shwit_switch ! Restart functionality public :: sh_write_phase_bin, sh_read_phase_bin @@ -445,6 +446,20 @@ subroutine sh_set_initialwf(initial_state) cel_re(initial_state) = 1.0D0 end subroutine sh_set_initialwf + ! A horrible hack for SHwIT, we just switch + ! switch the coefficients between S1-S0 + ! TODO: Do something less gross + subroutine shwit_switch() + real(DP) :: tmp + tmp = cel_re(1) + cel_re(1) = cel_re(2) + cel_re(2) = tmp + + tmp = cel_im(1) + cel_im(1) = cel_im(2) + cel_im(2) = tmp + end subroutine shwit_switch + subroutine sh_set_energy_shift(potential_energy) real(DP), intent(in) :: potential_energy diff --git a/src/surfacehop.F90 b/src/surfacehop.F90 index 921c9aeb..8c7476be 100644 --- a/src/surfacehop.F90 +++ b/src/surfacehop.F90 @@ -77,7 +77,12 @@ module mod_sh ! Stop simulation if total energy drift since the beginning exceeds energydriftthr real(DP) :: energydriftthr = 1.0D0 - ! Special case for adiabatic TDDFT, terminate when close to S1-S0 crossing + ! Surface Hopping with Induced Transitions, + ! typically used with methods that can't describe S1S0 intersections, + ! such as TDDFT or ADC2. + ! In this case, we force a hop to S0 when the energy difference + ! between S1 and S0 drops below user defined threshold. + logical :: SHwIT = .false. real(DP) :: dE_S0S1_thr = 0.0D0 !eV ! NA Couplings real(DP), allocatable :: nacx(:, :, :), nacy(:, :, :), nacz(:, :, :) @@ -100,7 +105,7 @@ module mod_sh namelist /sh/ istate_init, nstate, substep, deltae, integ, couplings, nohop, phase, decoh_alpha, popthr, ignore_state, & nac_accu1, nac_accu2, popsumthr, energydifthr, energydriftthr, velocity_rescaling, revmom, & - dE_S0S1_thr, correct_decoherence + shwit, dE_S0S1_thr, correct_decoherence save contains @@ -252,7 +257,7 @@ subroutine check_sh_parameters() write (stdout, '(A)') 'Using approximate Baeck-An couplings.' case ('none') inac = 2 - write (stdout, '(A)') 'Ignoring nonadaiabatic couplings.' + write (stdout, '(A)') 'Ignoring nonadiabatic couplings.' case default write (stderr, '(A)') 'Parameter "couplings" must be "analytic", "baeck-an" or "none".' error = .true. @@ -279,7 +284,7 @@ subroutine check_sh_parameters() error = .true. end if - if (inac == 2 .and. nohop == 0) then + if (inac == 2 .and. nohop == 0 .and. .not. shwit) then write (stdout, '(A)') 'WARNING: For simulations without couplings(="none") hopping probability cannot be determined.' nohop = 1 end if @@ -821,6 +826,12 @@ subroutine surfacehop(x, y, z, vx, vy, vz, vx_old, vy_old, vz_old, dt, eclas) !itp loop end do + if (SHwIT) then + call shwit_check(en_array(1), en_array(2), dE_S0S1_thr, & + vx, vy, vz, eclas) + end if + + ! set tocalc array for the next step call set_tocalc() @@ -1211,15 +1222,6 @@ subroutine check_energy(vx_old, vy_old, vz_old, vx, vy, vz) real(DP), intent(in) :: vx_old(:, :), vy_old(:, :), vz_old(:, :) real(DP) :: ekin, ekin_old, entot, entot_old - ! Special case for running AIMD with TDDFT: - ! End the simulation when S1-S0 energy difference drops - ! below a certain small threshold. - if (dE_S0S1_thr > 0.0D0 .and. nstate >= 2) then - call check_S0S1_gap(en_S0=en_array(1), & - en_S1=en_array(2), & - threshold_ev=dE_S0S1_thr) - end if - ekin = ekin_v(vx, vy, vz) ekin_old = ekin_v(vx_old, vy_old, vz_old) @@ -1242,26 +1244,29 @@ subroutine check_energy(vx_old, vy_old, vz_old, vx, vy, vz) end subroutine check_energy - subroutine check_S0S1_gap(en_S0, en_S1, threshold_ev) - use mod_general, only: STOP_SIMULATION + subroutine shwit_check(en_S0, en_S1, threshold_ev, vx, vy, vz, eclas) use mod_const, only: AUtoEV real(DP), intent(in) :: en_S0, en_S1 real(DP), intent(in) :: threshold_ev + real(DP), intent(inout) :: vx(:, :), vy(:, :), vz(:, :), eclas real(DP) :: dE_S0S1 + ! We only hop from S1 to S0 + if (istate /= 2) return + dE_S0S1 = (en_S1 - en_S0) * AUtoEV if (dE_S0S1 < threshold_ev) then - write (stdout, *) 'S1 - S0 gap dropped below threshold!' + write (stdout, *) 'SHwIT: S1 - S0 gap dropped below threshold!' write (stdout, *) dE_S0S1, ' < ', threshold_ev - ! TODO: This is an unfortunate hack, - ! but we can't directly stop here since we need the last step - ! to be written to the output files. - write (stdout, *) 'Simulation will stop at the end of this step.' + write (stdout, *) "Let's jump to S0 and continue!" ! This global flag is checked in the main loop in abin.F90 - STOP_SIMULATION = .true. + ! STOP_SIMULATION = .true. + call try_hop_simple_rescale(vx, vy, vz, 2, 1, eclas) + ! Horrible hack to exchange c_el coefficients between S1 - S0 + call shwit_switch() end if - end subroutine check_S0S1_gap + end subroutine shwit_check subroutine check_energydrift(vx, vy, vz) use mod_const, only: AUtoEV diff --git a/tests/SH_S0S1/input.in b/tests/SH_S0S1/input.in index 4b514e2b..6cd5954b 100644 --- a/tests/SH_S0S1/input.in +++ b/tests/SH_S0S1/input.in @@ -19,17 +19,18 @@ inose=0, / &sh +shwit=.true. dE_S0S1_thr=9.0 -istate_init=3, +istate_init=2, nstate=3, substep=100, deltae=100., integ='butcher', -couplings='analytic', ! non-adiabatic coupling terms 'analytic', 'baeck-an', 'none' +couplings='none', ! non-adiabatic coupling terms 'analytic', 'baeck-an', 'none' nohop=0, decoh_alpha=0.1 popthr=0.01 -energydifthr=0.50 -energydriftthr=0.50 +energydifthr=1.00 +energydriftthr=1.00 correct_decoherence=.false. /