From d0a1dcd731b5cdea4d7dd75807f2ff386940994d Mon Sep 17 00:00:00 2001 From: Reinhard Maurer Date: Sat, 9 Dec 2023 23:49:12 +0000 Subject: [PATCH 1/3] added new rescaling method ":vinversion" to abstract SurfaceHopping method. during a frustrated hop, the velocity along the NAC is inverted --- .../SurfaceHoppingMethods/surface_hopping.jl | 52 +++++++++++++++---- 1 file changed, 41 insertions(+), 11 deletions(-) diff --git a/src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl b/src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl index ec0f4239e..f6a574005 100644 --- a/src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl +++ b/src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl @@ -50,7 +50,7 @@ Rescale the velocity in the direction of the nonadiabatic coupling. [HammesSchiffer1994](@cite) """ function rescale_velocity!(sim::AbstractSimulation{<:SurfaceHopping}, u)::Bool - sim.method.rescaling === :off && return true + sim.method.rescaling === :off && return true #no rescaling so always accept hop new_state, old_state = unpack_states(sim) velocity = DynamicsUtils.get_hopping_velocity(sim, DynamicsUtils.get_velocities(u)) @@ -63,17 +63,29 @@ function rescale_velocity!(sim::AbstractSimulation{<:SurfaceHopping}, u)::Bool c = calculate_potential_energy_change(eigs, new_state, old_state) discriminant = b^2 - 4a * c - discriminant < 0 && return false # Frustrated hop with insufficient kinetic energy - - root = sqrt(discriminant) - if b < 0 - γ = (b + root) / 2a - else - γ = (b - root) / 2a + if discriminant < 0 # frustrated hop + if sim.method.rescaling == :standard + # Frustrated hop with insufficient kinetic energy, no velocity inversion + # Follow 1990 Tully recipe and do nothing + return false + elseif sim.method.rescaling == :vinversion + # Frustrated hop with insufficient kinetic energy + # perform inversion of the velocity component along the nonadiabatic coupling vector + frustrated_hop_invert_velocity!(sim, DynamicsUtils.get_velocities(u), d) + else + throw(error("This mode of rescaling is not implemented")) + end + else #sufficient energy for hopping + root = sqrt(discriminant) + if b < 0 + γ = (b + root) / 2a + else + γ = (b - root) / 2a + end + perform_rescaling!(sim, DynamicsUtils.get_velocities(u), γ, d) + + return true end - perform_rescaling!(sim, DynamicsUtils.get_velocities(u), γ, d) - - return true end """ @@ -123,6 +135,24 @@ function perform_rescaling!( return nothing end +""" + frustrated_hop_invert_velocity!( + sim::AbstractSimulation{<:SurfaceHopping}, velocity, d + ) + +Measures the component of velocity along the nonadiabatic coupling vector and inverts that component. +""" +function frustrated_hop_invert_velocity!( + sim::AbstractSimulation{<:SurfaceHopping}, velocity, d +) + dn = d/norm(d) + γ = dot(velocity,dn) + for I in CartesianIndices(dn) + velocity[I] -= 2γ * dn[I] + end + return nothing +end + function calculate_potential_energy_change(eigs::AbstractVector, new_state::Integer, current_state::Integer) return eigs[new_state] - eigs[current_state] end From 3bf00ff66972ddb7b9defd9bddd4f92f12c55693 Mon Sep 17 00:00:00 2001 From: Reinhard Maurer Date: Mon, 11 Dec 2023 16:32:18 +0000 Subject: [PATCH 2/3] Update src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl Co-authored-by: James Gardner <38594562+jamesgardner1421@users.noreply.github.com> --- src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl b/src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl index f6a574005..544c485b6 100644 --- a/src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl +++ b/src/DynamicsMethods/SurfaceHoppingMethods/surface_hopping.jl @@ -145,7 +145,7 @@ Measures the component of velocity along the nonadiabatic coupling vector and in function frustrated_hop_invert_velocity!( sim::AbstractSimulation{<:SurfaceHopping}, velocity, d ) - dn = d/norm(d) + dn = LinearAlgebra.normalize(d) γ = dot(velocity,dn) for I in CartesianIndices(dn) velocity[I] -= 2γ * dn[I] From 4098956d70066a33c2c8cad0d7b116df4b325ce4 Mon Sep 17 00:00:00 2001 From: Reinhard Maurer Date: Tue, 19 Dec 2023 11:36:27 +0000 Subject: [PATCH 3/3] added version of frustrated velocity inversion for Ring polymers --- src/DynamicsMethods/SurfaceHoppingMethods/rpsh.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/DynamicsMethods/SurfaceHoppingMethods/rpsh.jl b/src/DynamicsMethods/SurfaceHoppingMethods/rpsh.jl index 4e4315200..936bfb6e0 100644 --- a/src/DynamicsMethods/SurfaceHoppingMethods/rpsh.jl +++ b/src/DynamicsMethods/SurfaceHoppingMethods/rpsh.jl @@ -48,6 +48,19 @@ function perform_rescaling!( return nothing end +function frustrated_hop_invert_velocity!( + sim::RingPolymerSimulation{<:SurfaceHopping}, velocity, d +) + #invert center of mass velocity of ring polymer + dn = LinearAlgebra.normalize(d) + vcom = sum(velocity, dims=3) / size(velocity,3) + γ = dot(vcom,dn) + for I in CartesianIndices(dn) + velocity[I,:] -= 2γ * dn[I] + end + return nothing +end + function DynamicsUtils.classical_potential_energy(sim::RingPolymerSimulation{<:FSSH}, u) all_eigs = Calculators.get_eigen(sim.calculator, DynamicsUtils.get_positions(u)) potential = sum(eigs.values[u.state] for eigs in all_eigs)