From 8c7bbe980096f6dea8d44c28f22ba54a9aae8adf Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Tue, 11 Jul 2023 13:15:01 +0100 Subject: [PATCH 01/14] Correct for periodicity in quantise_diatomic --- .../QuantisedDiatomic/QuantisedDiatomic.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index 5197f13a3..394210e70 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -265,11 +265,20 @@ Otherwise, be sure to modify these parameters to give the intended behaviour. function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; height=10.0, surface_normal=[0, 0, 1.0], atom_indices=[1, 2], show_timer=false, reset_timer=false, - bond_lengths=0.5:0.01:5.0 + bond_lengths=0.5:0.01:5.0, max_translation=1 ) reset_timer && TimerOutputs.reset_timer!(TIMER) + if typeof(sim.cell)==PeriodicCell + # If the simulation used a `PeriodicCell`, translate `atom_indices` so they are at their minimum distance. (This is necessary if atoms were translated back into the original unit cell) + translations=[[i,j,k] for i in -max_translation:max_translation for j in -max_translation:max_translation for k in -max_translation:max_translation] + which_translation=argmin([norm(abs.(r[:,atom_indices[1]]-r[:,atom_indices[2]]+sim.cell.vectors*operation)) for operation in translations]) + # Translate one atom for minimal distance. + r[:,atom_indices[2]].=r[:,atom_indices[2]]+sim.cell.vectors*translations[which_translation] + @info "Atom 2 was moved as the simulation uses periodic boundary conditions" + end + r, slab = separate_slab_and_molecule(atom_indices, r) environment = EvaluationEnvironment(atom_indices, size(sim), slab, austrip(height), surface_normal) From 8bceb71ce1b2050cfee2dd213d9bb7b2129af100 Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Tue, 11 Jul 2023 14:57:46 +0100 Subject: [PATCH 02/14] Split velocity array into slab & molecule as well --- src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index 394210e70..9a593c583 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -280,6 +280,7 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; end 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) r_com = subtract_centre_of_mass(r, masses(sim)[atom_indices]) From a182247c6b5c86af0cab7ff0d5a92c3e6fc5e163 Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Wed, 12 Jul 2023 08:50:23 +0100 Subject: [PATCH 03/14] Checking for PeriodicCell now type invariant --- src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index 9a593c583..1268d1dfd 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -270,7 +270,7 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; reset_timer && TimerOutputs.reset_timer!(TIMER) - if typeof(sim.cell)==PeriodicCell + if isa(sim.cell,PeriodicCell) # If the simulation used a `PeriodicCell`, translate `atom_indices` so they are at their minimum distance. (This is necessary if atoms were translated back into the original unit cell) translations=[[i,j,k] for i in -max_translation:max_translation for j in -max_translation:max_translation for k in -max_translation:max_translation] which_translation=argmin([norm(abs.(r[:,atom_indices[1]]-r[:,atom_indices[2]]+sim.cell.vectors*operation)) for operation in translations]) From ba6bc2e9a25215123b98a99c8ec99f08fcfa4c3b Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Thu, 13 Jul 2023 09:18:21 +0100 Subject: [PATCH 04/14] Partial fix for #292 --- .../QuantisedDiatomic/QuantisedDiatomic.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index 1268d1dfd..ddbc3435c 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -265,7 +265,7 @@ Otherwise, be sure to modify these parameters to give the intended behaviour. function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; height=10.0, surface_normal=[0, 0, 1.0], atom_indices=[1, 2], show_timer=false, reset_timer=false, - bond_lengths=0.5:0.01:5.0, max_translation=1 + bond_lengths=0.5:0.01:5.0, max_translation=1, debug=false ) reset_timer && TimerOutputs.reset_timer!(TIMER) @@ -275,8 +275,10 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; translations=[[i,j,k] for i in -max_translation:max_translation for j in -max_translation:max_translation for k in -max_translation:max_translation] which_translation=argmin([norm(abs.(r[:,atom_indices[1]]-r[:,atom_indices[2]]+sim.cell.vectors*operation)) for operation in translations]) # Translate one atom for minimal distance. - r[:,atom_indices[2]].=r[:,atom_indices[2]]+sim.cell.vectors*translations[which_translation] - @info "Atom 2 was moved as the simulation uses periodic boundary conditions" + if translations[which_translation]!=[0,0,0] + r[:,atom_indices[2]].=r[:,atom_indices[2]]+sim.cell.vectors*translations[which_translation] + @info "Using a periodic copy of atom "*string(atom_indices[end])*" to bring it closer to atom "*string(atom_indices[begin]) + end end r, slab = separate_slab_and_molecule(atom_indices, r) @@ -287,7 +289,7 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; v_com = subtract_centre_of_mass(v, masses(sim)[atom_indices]) p_com = v_com .* masses(sim)[atom_indices]' - k = DynamicsUtils.classical_kinetic_energy(sim, v_com) + k = DynamicsUtils.classical_kinetic_energy(masses(sim)[atom_indices], v_com) p = calculate_diatomic_energy(bond_length(r_com), sim.calculator.model, environment) E = k + p @@ -298,6 +300,10 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; binding_curve = calculate_binding_curve(bond_lengths, sim.calculator.model, environment) + if debug + println(k,"\n", p,"\n", E,"\n", L,"\n", J) + end + V = EffectivePotential(μ, J, binding_curve) r₁, r₂ = @timeit TIMER "Finding bounds" find_integral_bounds(E, V) From 49e927df1f6c0222502401e91587d7d225488690 Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Tue, 15 Aug 2023 10:56:59 +0100 Subject: [PATCH 05/14] Add debug info --- .../QuantisedDiatomic/QuantisedDiatomic.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index ddbc3435c..aabebc734 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -265,7 +265,7 @@ Otherwise, be sure to modify these parameters to give the intended behaviour. function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; height=10.0, surface_normal=[0, 0, 1.0], atom_indices=[1, 2], show_timer=false, reset_timer=false, - bond_lengths=0.5:0.01:5.0, max_translation=1, debug=false + bond_lengths=0.5:0.01:5.0, max_translation=1 ) reset_timer && TimerOutputs.reset_timer!(TIMER) @@ -277,7 +277,7 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; # Translate one atom for minimal distance. if translations[which_translation]!=[0,0,0] r[:,atom_indices[2]].=r[:,atom_indices[2]]+sim.cell.vectors*translations[which_translation] - @info "Using a periodic copy of atom "*string(atom_indices[end])*" to bring it closer to atom "*string(atom_indices[begin]) + @debug "Using a periodic copy of atom "*string(atom_indices[end])*" to bring it closer to atom "*string(atom_indices[begin]) end end @@ -300,10 +300,8 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; binding_curve = calculate_binding_curve(bond_lengths, sim.calculator.model, environment) - if debug - println(k,"\n", p,"\n", E,"\n", L,"\n", J) - end - + @debug "k=$k\np=$p\nE=$E\nL=$L\nJ=$J\n" + V = EffectivePotential(μ, J, binding_curve) r₁, r₂ = @timeit TIMER "Finding bounds" find_integral_bounds(E, V) From e4cc66fee8b44f8c4ab1086993c48c905010fea2 Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Tue, 15 Aug 2023 11:02:48 +0100 Subject: [PATCH 06/14] More debug info --- src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index aabebc734..29faaef8c 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -317,6 +317,7 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; TimerOutputs.complement!(TIMER) show_timer && show(TIMER) + @debug "Found ν=$ν" return round(Int, ν), J end From 45eaa64a982fdbc433af1fab70144086687cb2da Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Wed, 16 Aug 2023 09:24:23 +0100 Subject: [PATCH 07/14] Possibly translating in the wrong direction --- src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index 29faaef8c..a853468ca 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -273,7 +273,7 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; if isa(sim.cell,PeriodicCell) # If the simulation used a `PeriodicCell`, translate `atom_indices` so they are at their minimum distance. (This is necessary if atoms were translated back into the original unit cell) translations=[[i,j,k] for i in -max_translation:max_translation for j in -max_translation:max_translation for k in -max_translation:max_translation] - which_translation=argmin([norm(abs.(r[:,atom_indices[1]]-r[:,atom_indices[2]]+sim.cell.vectors*operation)) for operation in translations]) + which_translation=argmin([norm(abs.(r[:,atom_indices[2]]-r[:,atom_indices[1]]+sim.cell.vectors*operation)) for operation in translations]) # Translate one atom for minimal distance. if translations[which_translation]!=[0,0,0] r[:,atom_indices[2]].=r[:,atom_indices[2]]+sim.cell.vectors*translations[which_translation] From 16f804c3fd528d097dd0b40d5dc1264bdda5852b Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Wed, 16 Aug 2023 09:58:50 +0100 Subject: [PATCH 08/14] More debug outputs --- .../QuantisedDiatomic/QuantisedDiatomic.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index a853468ca..e59d62216 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -256,7 +256,10 @@ end height=10, normal_vector=[0, 0, 1]) Quantise the vibrational and rotational degrees of freedom for the specified -positions and velocities +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. @@ -285,6 +288,8 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix; v, slab_v = separate_slab_and_molecule(atom_indices, v) environment = EvaluationEnvironment(atom_indices, size(sim), slab, austrip(height), surface_normal) + @debug "After PBC check and separation from slab, diatomic positions are:\n$(r)" + r_com = subtract_centre_of_mass(r, masses(sim)[atom_indices]) v_com = subtract_centre_of_mass(v, masses(sim)[atom_indices]) p_com = v_com .* masses(sim)[atom_indices]' From 11c0e18a6574292e2c2f2c664d11cd9a575a2e67 Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Fri, 18 Aug 2023 15:38:26 +0100 Subject: [PATCH 09/14] Fix unrelated typo --- src/DynamicsOutputs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DynamicsOutputs.jl b/src/DynamicsOutputs.jl index 8fec42162..c147374dd 100644 --- a/src/DynamicsOutputs.jl +++ b/src/DynamicsOutputs.jl @@ -64,7 +64,7 @@ export OutputPosition Output the position of the ring polymer centroid at each timestep during the trajectory. """ -OutputCentroidPosition(sol, i) = [get_centroid(get_position(u)) for u in sol.u] +OutputCentroidPosition(sol, i) = [get_centroid(get_positions(u)) for u in sol.u] export OutputCentroidPosition """ From e57ca70bd5cd57d1120c0d2fc73d0f23fd6b3a19 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Fri, 18 Aug 2023 16:45:24 +0100 Subject: [PATCH 10/14] Improve accuracy of equilibrium bond length --- .../QuantisedDiatomic/binding_curve.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/InitialConditions/QuantisedDiatomic/binding_curve.jl b/src/InitialConditions/QuantisedDiatomic/binding_curve.jl index 8b82f1bc1..48aa1a07d 100644 --- a/src/InitialConditions/QuantisedDiatomic/binding_curve.jl +++ b/src/InitialConditions/QuantisedDiatomic/binding_curve.jl @@ -1,3 +1,4 @@ +using Optim: Optim struct BindingCurve{T,B,F} bond_lengths::B @@ -11,9 +12,8 @@ function calculate_binding_curve( bond_lengths::AbstractVector, model::AdiabaticModel, environment::EvaluationEnvironment ) potential = calculate_diatomic_energy.(bond_lengths, model, environment) # Calculate binding curve - potential_minimum, index = findmin(potential) - equilibrium_bond_length = bond_lengths[index] fit = fit_binding_curve(bond_lengths, potential) + equilibrium_bond_length, potential_minimum = find_minimum(fit) return BindingCurve(bond_lengths, potential, equilibrium_bond_length, potential_minimum, fit) end @@ -23,6 +23,12 @@ function fit_binding_curve(bond_lengths, binding_curve) return sitp end +function find_minimum(fit) + grid = knots(fit) + results = Optim.optimize(fit, minimum(grid), maximum(grid), Optim.Brent()) + return Optim.minimizer(results), Optim.minimum(results) +end + function calculate_force_constant(binding_curve) return hessian(binding_curve.fit, binding_curve.equilibrium_bond_length)[1] end From cc1414f5a6bbd319d2dfb846436664838e87ce59 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Fri, 18 Aug 2023 17:07:43 +0100 Subject: [PATCH 11/14] Fix missing import --- src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index e59d62216..d482c16e4 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -27,7 +27,7 @@ using UnicodePlots: lineplot, lineplot!, DotCanvas, histogram using Optim: Optim using Roots: Roots using TimerOutputs: TimerOutputs, @timeit -using Interpolations: interpolate, BSpline, Cubic, Line, OnGrid, scale, hessian +using Interpolations: interpolate, BSpline, Cubic, Line, OnGrid, scale, hessian, knots using NQCDynamics: Simulation, Calculators, DynamicsUtils, masses using NQCBase: Atoms, PeriodicCell, InfiniteCell From d996f6df5f53cae628ce105de603d4ab7830251f Mon Sep 17 00:00:00 2001 From: Alexander Spears Date: Mon, 21 Aug 2023 11:44:51 +0100 Subject: [PATCH 12/14] 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 d482c16e4..2751596c3 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 @@ -251,12 +252,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. @@ -264,11 +375,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) @@ -303,8 +419,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) From c08efa9a04451a9643329f3632c7a392fd12a806 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Mon, 21 Aug 2023 17:20:47 +0100 Subject: [PATCH 13/14] Delay rounding of J values until return --- .../QuantisedDiatomic/QuantisedDiatomic.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index 2751596c3..8b075835a 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -44,9 +44,9 @@ include("energy_evaluation.jl") include("binding_curve.jl") include("random_configuration.jl") -struct EffectivePotential{T,B,F} +struct EffectivePotential{JType,T,B,F} μ::T - J::Int + J::JType binding_curve::BindingCurve{T,B,F} end @@ -415,7 +415,7 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix, binding_curve: E = k + p L = total_angular_momentum(r_com, p_com) - J = round(Int, (sqrt(1 + 4 * L^2) - 1) / 2) # L^2 = J(J+1)ħ^2 + J = (sqrt(1+4*L^2) - 1) / 2 # L^2 = J(J+1)ħ^2 μ = reduced_mass(masses(sim)[atom_indices]) @@ -437,7 +437,7 @@ function quantise_diatomic(sim::Simulation, v::Matrix, r::Matrix, binding_curve: show_timer && show(TIMER) @debug "Found ν=$ν" - return round(Int, ν), J + return round(Int, ν), round(Int, J) end function quantise_1D_vibration(model::AdiabaticModel, μ::Real, r::Real, v::Real; From 9786fdd1953e0235d85bd3477f8ce844f07e1085 Mon Sep 17 00:00:00 2001 From: James Gardner Date: Fri, 25 Aug 2023 10:20:26 +0100 Subject: [PATCH 14/14] Fix type definition --- src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl index 8b075835a..403f0039e 100644 --- a/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl +++ b/src/InitialConditions/QuantisedDiatomic/QuantisedDiatomic.jl @@ -58,9 +58,9 @@ function (effective_potential::EffectivePotential)(r) return rotational + potential end -struct RadialMomentum{T,B,F} +struct RadialMomentum{JType,T,B,F} total_energy::T - V::EffectivePotential{T,B,F} + V::EffectivePotential{JType,T,B,F} end function (radial_momentum::RadialMomentum)(r)