From 5b992e50d58c5988da308d050645b8c43bf1e677 Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Mon, 21 Aug 2023 11:44:51 +0100 Subject: [PATCH] Add vectorised version of `quantise_diatomic` `quantise_diatomic` is now able to process a single configuration or a vector of multiple configurations. This improves performance for multiple configurations, as only one potential-specific binding curve is generated to analyse all configurations. --- .../QuantisedDiatomic/QuantisedDiatomic.jl | 130 ++++++++++++++++-- 1 file changed, 122 insertions(+), 8 deletions(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index f2330ab46..b9c435221 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -27,6 +27,7 @@ using UnicodePlots: lineplot, lineplot!, DotCanvas, histogram using Optim: Optim using Roots: Roots using TimerOutputs: TimerOutputs, @timeit +using ProgressMeter using Interpolations: interpolate, BSpline, Cubic, Line, OnGrid, scale, hessian, knots using NQCDynamics: Simulation, Calculators, DynamicsUtils, masses @@ -230,12 +231,122 @@ function find_integral_bounds(total_energy::Real, V::EffectivePotential) return r₁, r₂ end +function binding_curve_from_structure(sim::Simulation, v::Matrix, r::Matrix; + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], + args... +) + # Separate slab from molecule if necessary. + r, slab = separate_slab_and_molecule(atom_indices, r) + v, slab_v = separate_slab_and_molecule(atom_indices, v) + environment = EvaluationEnvironment(atom_indices, size(sim), slab, austrip(height), surface_normal) + + binding_curve = calculate_binding_curve(bond_lengths, sim.calculator.model, environment) + @debug "Finished generating binding curve for the given potential." + + return binding_curve +end + +""" + quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], + max_translation=1, + show_timer=false, reset_timer=false + ) + +Quantise the vibrational and rotational degrees of freedom for the specified +positions and velocities. + +If the potential can be evaluated for the diatomic only, independent of position, +supplying a `Simulation` for just the diatomic will speed up evaluation. + +When evaluating the potential, the molecule is moved to `height` in direction `normal_vector`. +If the potential is independent of centre of mass position, this has no effect. +Otherwise, be sure to modify these parameters to give the intended behaviour. + +If a `Simulation` with a `PeriodicCell` is supplied, periodic copies of the diatomic atoms +will be used if positions are close to cell boundaries. +Set `max_translation` to the radius of surrounding unit cells to search. +(e.g. 1 if positions are already wrapped around cell boundaries) +""" +function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], + args... +) + binding_curve=binding_curve_from_structure(sim, v, r; bond_lengths=bond_lengths, height=height, surface_normal=surface_normal, atom_indices=atom_indices) + return quantise_diatomic(sim, v, r, binding_curve; height=height, surface_normal=surface_normal, atom_indices=atom_indices, args...) +end + +""" + quantise_diatomic(sim::Simulation, v::Vector{Matrix}, r::Vector{Matrix}; + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], + show_timer=false, reset_timer=false + ) + +Quantise the vibrational and rotational degrees of freedom of multiple atomic configurations +given as a vector of velocity matrices and a vector of position matrices. + +If the potential can be evaluated for the diatomic only, independent of position, +supplying a `Simulation` for just the diatomic will speed up evaluation. + +When evaluating the potential, the molecule is moved to `height` in direction `normal_vector`. +If the potential is independent of centre of mass position, this has no effect. +Otherwise, be sure to modify these parameters to give the intended behaviour. + +If a `Simulation` with a `PeriodicCell` is supplied, periodic copies of the diatomic atoms +will be used if positions are close to cell boundaries. +Set `max_translation` to the radius of surrounding unit cells to search. +(e.g. 1 if positions are already wrapped around cell boundaries) + +Specify `show_timer=true` for performance timings of the EBK quantisation process and +`reset_timer=true` to see timings for each individual quantisation. +""" +function quantise_diatomic(sim::Simulation, v::Vector{<: Matrix{<: Any}}, r::Vector{<: Matrix{<: Any}}; + bond_lengths=0.5:0.01:5.0, + height=10.0, + surface_normal=[0, 0, 1.0], + atom_indices=[1, 2], + args... +) + # Generate binding curve for all structures + binding_curve=binding_curve_from_structure(sim, v[1], r[1]; bond_lengths=bond_lengths, height=height, surface_normal=surface_normal, atom_indices=atom_indices) + ν, J=[],[] + + # Quantise all configurations in the given vectors. + @showprogress for index in eachindex(v) + try + result=quantise_diatomic(sim, v[index], r[index], binding_curve; height=height, surface_normal=surface_normal, atom_indices=atom_indices, args...) + push!(ν, result[1]) + push!(J, result[2]) + catch e + @warn "Quantisation was unsuccessful for configuration at $(index).\nThe final results vectors will show ν=-1 and J=-1 for this configuration.\n The error was: $(e)" + push!(ν, -1) + push!(J, -1) + end + end + return ν, J +end + """ - quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; - height=10, normal_vector=[0, 0, 1]) + quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix, binding_curve::BindingCurve; + show_timer=false, reset_timer=false, + height=10, normal_vector=[0, 0, 1], atom_indices=[1,2], max_translation=1) + ) Quantise the vibrational and rotational degrees of freedom for the specified -positions and velocities. +positions and velocities using the `BindingCurve` specified. +A binding curve will be automatically generated if you do not supply one. If the potential can be evaluated for the diatomic only, independent of position, supplying a `Simulation` for just the diatomic will speed up evaluation. @@ -243,11 +354,16 @@ supplying a `Simulation` for just the diatomic will speed up evaluation. When evaluating the potential, the molecule is moved to `height` in direction `normal_vector`. If the potential is independent of centre of mass position, this has no effect. Otherwise, be sure to modify these parameters to give the intended behaviour. + +If a `Simulation` with a `PeriodicCell` is supplied, periodic copies of the diatomic atoms +will be used if positions are close to cell boundaries. +Set `max_translation` to the radius of surrounding unit cells to search. +(e.g. 1 if positions are already wrapped around cell boundaries) """ -function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; - height=10.0, surface_normal=[0, 0, 1.0], atom_indices=[1, 2], +function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix, binding_curve::BindingCurve; show_timer=false, reset_timer=false, - bond_lengths=0.5:0.01:5.0, max_translation=1 + max_translation=1, height=10.0, + surface_normal=[0, 0, 1.0], atom_indices=[1, 2], ) reset_timer && TimerOutputs.reset_timer!(TIMER) @@ -282,8 +398,6 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; μ = reduced_mass(masses(sim)[atom_indices]) - binding_curve = calculate_binding_curve(bond_lengths, sim.calculator.model, environment) - @debug "k=$k\np=$p\nE=$E\nL=$L\nJ=$J\n" V = EffectivePotential(μ, J, binding_curve)