diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 6c9fa27..517deb0 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -18,7 +18,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - julia-version: ['1'] + julia-version: ['1', 'nightly'] julia-arch: [x64] os: [ubuntu-latest, windows-latest, macOS-latest] exclude: diff --git a/Project.toml b/Project.toml index d92df86..194d003 100644 --- a/Project.toml +++ b/Project.toml @@ -7,14 +7,11 @@ version = "0.5.0" FunctionWrappers = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" -Optim = "429524aa-4258-5aef-a3af-852621145aeb" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PolynomialRoots = "3a141323-8675-5d76-9d11-e1df1406c778" QuantumOptics = "6e0679c1-51ea-5a7c-ac74-d61b76210b0c" QuantumOpticsBase = "4f57444f-1401-5e15-980d-4471b28d5678" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" WignerSymbols = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" diff --git a/README.md b/README.md index e1ee94b..cb71593 100644 --- a/README.md +++ b/README.md @@ -38,8 +38,7 @@ julia> using Pkg julia> Pkg.add("IonSim") julia> Pkg.add("IJulia") ``` -The last line adds Jupyter compatibility, so you'll have the option to start a Julia kernel from Jupyter notebook. This is the suggested way -to interact with IonSim. +The last line adds Jupyter compatibility, so you'll have the option to start a Julia kernel from Jupyter notebook. ### Updating IonSim diff --git a/configs/ions/be9.yaml b/configs/ions/be9.yaml new file mode 100644 index 0000000..4edca12 --- /dev/null +++ b/configs/ions/be9.yaml @@ -0,0 +1,36 @@ +name: "Be9" +shortname: "⁹Be" +mass: 1.496508080073e-26 +charge: 1 +nuclearspin: 0.5 +full_level_structure: + S1/2f=1: {n: 2, l: 0, j: 0.5, f: 1, E: 0.78126104631e+9} + S1/2f=2: {n: 2, l: 0, j: 0.5, f: 2, E: -0.468756627786e+9} + P1/2f=1: {n: 2, l: 1, j: 0.5, f: 1, E: 957.2010729076436e+12} + P1/2f=2: {n: 2, l: 1, j: 0.5, f: 2, E: 957.2008357076436e+12} + P3/2f=0: {n: 2, l: 1, j: 1.5, f: 0, E: 957.3965669467407e+12} + P3/2f=1: {n: 2, l: 1, j: 1.5, f: 1, E: 957.3965659267407e+12} + P3/2f=2: {n: 2, l: 1, j: 1.5, f: 2, E: 957.3965638867406e+12} + P3/2f=3: {n: 2, l: 1, j: 1.5, f: 3, E: 957.3965608267406e+12} +full_transitions: + - {from: S1/2f=1, to: P1/2f=1, multipole: E1, einsteinA: 19.4000000e+9} + - {from: S1/2f=1, to: P1/2f=2, multipole: E1, einsteinA: 19.4000000e+9} + - {from: S1/2f=2, to: P1/2f=1, multipole: E1, einsteinA: 19.4000000e+9} + - {from: S1/2f=2, to: P1/2f=2, multipole: E1, einsteinA: 19.4000000e+9} + - {from: S1/2f=1, to: P3/2f=0, multipole: E1, einsteinA: 19.4000000e+9} + - {from: S1/2f=1, to: P3/2f=1, multipole: E1, einsteinA: 19.4000000e+9} + - {from: S1/2f=1, to: P3/2f=2, multipole: E1, einsteinA: 19.4000000e+9} + - {from: S1/2f=2, to: P1/2f=1, multipole: E1, einsteinA: 19.4000000e+9} + - {from: S1/2f=2, to: P1/2f=2, multipole: E1, einsteinA: 19.4000000e+9} + - {from: S1/2f=2, to: P1/2f=3, multipole: E1, einsteinA: 19.4000000e+9} + +# Optional fields +default_sublevel_selection: +- ["S1/2f=1", "all"] +- ["S1/2f=2", "all"] +- ["P1/2f=1", "all"] +- ["P1/2f=2", "all"] +- ["P3/2f=0", "all"] +- ["P3/2f=1", "all"] +- ["P3/2f=2", "all"] +- ["P3/2f=3", "all"] diff --git a/configs/ions/ca40.yaml b/configs/ions/ca40.yaml new file mode 100644 index 0000000..32099ea --- /dev/null +++ b/configs/ions/ca40.yaml @@ -0,0 +1,28 @@ +name: "Ca40" +shortname: "⁴⁰Ca" +mass: 6.635943757345042e-26 +charge: 1 +nuclearspin: 0 +full_level_structure: + S1/2: {n: 4, l: 0, j: 0.5, f: 0.5 , E: 0.0} + D3/2: {n: 3, l: 2, j: 1.5, f: 1.5 , E: 4.09335071228e+14} + D5/2: {n: 3, l: 2, j: 2.5, f: 2.5 , E: 4.1115503183857306e+14} + P1/2: {n: 4, l: 1, j: 0.5, f: 0.5 , E: 7.554e+14} + P3/2: {n: 4, l: 1, j: 1.5, f: 1.5 , E: 7.621e+14} +full_transitions: + - {from: "S1/2", to: "D5/2", multipole: "E2", einsteinA: 8.562e-1} + - {from: "S1/2", to: "P1/2", multipole: "E1", einsteinA: 1.299e+8} + - {from: "D3/2", to: "P1/2", multipole: "E1", einsteinA: 1.060e+7} + - {from: "S1/2", to: "D3/2", multipole: "E2", einsteinA: 9.259e-1} + - {from: "S1/2", to: "P3/2", multipole: "E1", einsteinA: 1.351e+8} + - {from: "D3/2", to: "P3/2", multipole: "E1", einsteinA: 1.110e+6} + - {from: "D5/2", to: "P3/2", multipole: "E1", einsteinA: 9.901e+6} + +# Optional fields +default_sublevel_selection: +- ["S1/2", "all"] +- ["D5/2", "all"] + +gfactors: + S1/2: 2.00225664 + D5/2: 1.2003340 \ No newline at end of file diff --git a/configs/ions/mg25.yaml b/configs/ions/mg25.yaml new file mode 100644 index 0000000..5316337 --- /dev/null +++ b/configs/ions/mg25.yaml @@ -0,0 +1,37 @@ +name: "Mg25" +shortname: "²⁵Mg" +mass: 4.1489954812e-26 +charge: 1 +nuclearspin: 2.5 +full_level_structure: + S1/2f=2: {n: 3, l: 0, j: 0.5, f: 2, E: 1.043445158e+9} + S1/2f=3: {n: 3, l: 0, j: 0.5, f: 3, E: -0.745317970e+9} + P1/2f=2: {n: 3, l: 1, j: 0.5, f: 2, E: 1069.3408690519877e+12} + P1/2f=3: {n: 3, l: 1, j: 0.5, f: 3, E: 1069.3405639519879e+12} + P3/2f=1: {n: 3, l: 1, j: 1.5, f: 1, E: 1072.0853411194748e+12} + P3/2f=2: {n: 3, l: 1, j: 1.5, f: 2, E: 1072.0853033394746e+12} + P3/2f=3: {n: 3, l: 1, j: 1.5, f: 3, E: 1072.0852466694748e+12} + P3/2f=4: {n: 3, l: 1, j: 1.5, f: 4, E: 1072.0851711094747e+12} + +full_transitions: + - {from: "S1/2f=2", to: "P1/2f=2", multipole: "E1", einsteinA: 41.3e+6} + - {from: "S1/2f=2", to: "P1/2f=3", multipole: "E1", einsteinA: 41.3e+6} + - {from: "S1/2f=3", to: "P1/2f=2", multipole: "E1", einsteinA: 41.3e+6} + - {from: "S1/2f=3", to: "P1/2f=3", multipole: "E1", einsteinA: 41.3e+6} + - {from: "S1/2f=2", to: "P3/2f=1", multipole: "E1", einsteinA: 41.3e+6} + - {from: "S1/2f=2", to: "P3/2f=2", multipole: "E1", einsteinA: 41.3e+6} + - {from: "S1/2f=2", to: "P3/2f=3", multipole: "E1", einsteinA: 41.3e+6} + - {from: "S1/2f=3", to: "P1/2f=2", multipole: "E1", einsteinA: 41.3e+6} + - {from: "S1/2f=3", to: "P1/2f=3", multipole: "E1", einsteinA: 41.3e+6} + - {from: "S1/2f=3", to: "P1/2f=4", multipole: "E1", einsteinA: 41.3e+6} + +# Optional fields +default_sublevel_selection: + - ["S1/2f=2", "all"] + - ["S1/2f=3", "all"] + - ["P1/2f=2", "all"] + - ["P1/2f=3", "all"] + - ["P3/2f=1", "all"] + - ["P3/2f=2", "all"] + - ["P3/2f=3", "all"] + - ["P3/2f=4", "all"] \ No newline at end of file diff --git a/configs/ions/yb171.yaml b/configs/ions/yb171.yaml new file mode 100644 index 0000000..f6f1741 --- /dev/null +++ b/configs/ions/yb171.yaml @@ -0,0 +1,80 @@ +name: "Yb171" +shortname: "¹⁷¹Yb" +mass: 28.8384644689030595108e-26 +charge: 1 +nuclearspin: 0.5 +full_level_structure: + "S1/2f=0": {n: 6, l: 0, j: 0.5, f: 0, E: -9.48210909315e+9} + "S1/2f=1": {n: 6, l: 0, j: 0.5, f: 1, E: 3.16070303105e+9} + "D3/2f=0": {n: 5, l: 2, j: 1.5, f: 1, E: 6.88342460964640e+14} + "D3/2f=1": {n: 5, l: 2, j: 1.5, f: 2, E: 6.88350470564640e+14} + "D5/2f=2": {n: 5, l: 2, j: 2.5, f: 2, E: 7.29477895985202e+14} + "D5/2f=3": {n: 5, l: 2, j: 2.5, f: 3, E: 7.29474121985202e+14} + "P1/2f=0": {n: 6, l: 1, j: 0.5, f: 0, E: 8.112913749003559e+14} + "P1/2f=1": {n: 6, l: 1, j: 0.5, f: 1, E: 8.112934798003559e+14} + "[3/2]1/2f=0": {n: 6, l: null, j: 0.5, f: 0, E: 1.008917341058788e+15} + "[3/2]1/2f=1": {n: 6, l: null, j: 0.5, f: 1, E: 1.008917341058788e+15} + "[3/2]3/2f=1": {n: 6, l: null, j: 1.5, f: 1, E: 8.621425511314839e+14} + "[3/2]3/2f=2": {n: 6, l: null, j: 1.5, f: 2, E: 8.621425511314839e+14} + "[5/2]5/2f=2": {n: 6, l: null, j: 2.5, f: 2, E: 9.70461163716380e+14} + "[5/2]5/2f=3": {n: 6, l: null, j: 2.5, f: 3, E: 9.70461163716380e+14} + "F7/2f=3": {n: 6, l: 3, j: 3.5, f: 3, E: 6.42115934728750e+14} + "F7/2f=4": {n: 6, l: 3, j: 3.5, f: 4, E: 6.42119554728750e+14} + +full_transitions: + # Laser lines + # 411nm + - {from: "S1/2f=0", to: "D5/2f=2", multipole: "E2", einsteinA: 22} + - {from: "S1/2f=1", to: "D5/2f=2", multipole: "E2", einsteinA: 22} + - {from: "S1/2f=1", to: "D5/2f=3", multipole: "E2", einsteinA: 22} + # 369nm + - {from: "S1/2f=0", to: "P1/2f=1", multipole: "E1", einsteinA: 1.155e+8} + - {from: "S1/2f=1", to: "P1/2f=0", multipole: "E1", einsteinA: 1.155e+8} + - {from: "S1/2f=1", to: "P1/2f=1", multipole: "E1", einsteinA: 1.155e+8} + # 935nm + - {from: "D3/2f=1", to: "[3/2]1/2f=0", multipole: "E1", einsteinA: 1.2e+5} + - {from: "D3/2f=1", to: "[3/2]1/2f=1", multipole: "E1", einsteinA: 1.2e+5} + - {from: "D3/2f=2", to: "[3/2]1/2f=1", multipole: "E1", einsteinA: 1.2e+5} + # 760nm + - {from: "F7/2f=3", to: "[3/2]3/2f=1", multipole: "E2", einsteinA: 5.0e+4} + - {from: "F7/2f=3", to: "[3/2]3/2f=2", multipole: "E2", einsteinA: 5.0e+4} + - {from: "F7/2f=4", to: "[3/2]3/2f=2", multipole: "E2", einsteinA: 5.0e+4} + # 638nm + # 976nm + # 861nm + # Decay lines + # P1/2 -> D3/2 + - {from: "P1/2f=0", to: "D3/2f=1", multipole: "E1", einsteinA: 5.7e+5} + - {from: "P1/2f=1", to: "D3/2f=2", multipole: "E1", einsteinA: 5.7e+5} + - {from: "P1/2f=1", to: "D3/2f=2", multipole: "E1", einsteinA: 5.7e+5} + # [3/2]1/2 -> S1/2 + - {from: "[3/2]1/2f=0", to: "S1/2f=1", multipole: "E1", einsteinA: 8.05e+7} + - {from: "[3/2]1/2f=1", to: "S1/2f=0", multipole: "E1", einsteinA: 8.05e+7} + - {from: "[3/2]1/2f=1", to: "S1/2f=1", multipole: "E1", einsteinA: 8.05e+7} + # [3/2]3/2 -> S1/2 + - {from: "[3/2]1/2f=0", to: "S1/2f=1", multipole: "E1", einsteinA: 5.125e+7} + - {from: "[3/2]1/2f=1", to: "S1/2f=0", multipole: "E1", einsteinA: 5.125e+7} + - {from: "[3/2]1/2f=1", to: "S1/2f=1", multipole: "E1", einsteinA: 5.125e+7} + +# Optional fields +default_sublevel_selection: + - ["S1/2f=0", "all"] + - ["S1/2f=1", "all"] + - ["P1/2f=0", "all"] + - ["P1/2f=1", "all"] + +gfactors: + "S1/2f=0": 1.998 + "S1/2f=1": 1.998 + "D5/2f=2": 1.202 + "D5/2f=3": 1.202 + "F7/2f=3": 1.145 + "F7/2f=4": 1.145 + "[3/2]1/2f=0": 1.32 + "[3/2]1/2f=1": 1.32 + "[3/2]3/2f=1": 1.44 + "[3/2]3/2f=2": 1.44 + +nonlinear_zeeman: +- {level: "S1/2f=0", sublevel: 0, coeffs: [0.0, 0.0, -155.305]} +- {level: "S1/2f=1", sublevel: 0, coeffs: [0.0, 0.0, 155.305]} diff --git a/src/IonSim.jl b/src/IonSim.jl index 80116e0..f713f27 100644 --- a/src/IonSim.jl +++ b/src/IonSim.jl @@ -34,6 +34,7 @@ abstract type IonSimBasis <: Basis end export IonSimBasis include("constants.jl") +include("ionproperties.jl") include("ions.jl") include("vibrationalmodes.jl") include("lasers.jl") @@ -42,7 +43,6 @@ include("chambers.jl") include("operators.jl") include("hamiltonians.jl") include("timeevolution.jl") -include("species/_include_species.jl") module analytical include("analyticfunctions.jl") diff --git a/src/hamiltonians.jl b/src/hamiltonians.jl index ae5132c..c7b2bb0 100644 --- a/src/hamiltonians.jl +++ b/src/hamiltonians.jl @@ -14,7 +14,7 @@ Constructs the Hamiltonian for `chamber` as a function of time. Return type is a `h(t::Real, ψ)` that, itself, returns a `QuantumOptics.SparseOperator`. **args** -* `timescale`: e.g. a value of 1e-6 will take time to be in μs +* `timescale`: e.g. a value of 1e-6 will take time to be in ``\\mu s`` * `lamb_dicke_order`: Only consider terms that change the phonon number by up to this value. If this is an `Int`, then the cutoff is applied to all modes. If this is a `Vector{Int}`, then `lamb_dicke_order[i]` is applied to the iᵗʰ mode, according to the order in diff --git a/src/ionproperties.jl b/src/ionproperties.jl new file mode 100644 index 0000000..9c776a1 --- /dev/null +++ b/src/ionproperties.jl @@ -0,0 +1,165 @@ +# This file contains functionality for tracking and manipulating ion properties. +module Properties # Open to better names for this module. + +import YAML +using LinearAlgebra: dot +using OrderedCollections + +export IonProperties, loadfromconfig + +function roundnearesthalf(value::Number)::Rational + halves::Integer = 2 * value + return halves // 2 +end + + +""" +IonProperties type. +**Required keywords** +* `mass`: Mass of ion in kg +* `charge`: Charge of ion in units of elementary charge +* `full_level_structure`: OrderedDict describing properties of each energy level + * `key::String`: Name of energy level. Spectroscopic notation encouraged, e.g. `"S1/2,f=1"` + * `value::NamedTuple(:n, :l, :j, :f, :E)`: Quantum numbers `n`, `l`, `j`, `f`, and energy `E` (in Hz) +* `full_transitions`: Dict of all allowed transitions between energy levels + * `key::Tuple{String,String}` Pair of levels, ordered (lower, upper) in energy + * `value::NamedTuple(:multipole, :einsteinA)`: Leading-order multipole of the transition (e.g. `"E1"`, `"E2"`) and Einstein A coefficient (between fine structure levels only; hyperfine factors are calculated when needed) + + **Optional keywords** + * `default_sublevel_selection`: Default value of `selected_sublevels` argument in Ion constructor + * `gfactors`: `Dict(level::String => g::Real)` Custom Landé g-factors, if contributions from higher-than-first-order perturbations are desired + * `nonlinear_zeeman`: `Dict` describing nonlinear contributions to Zeeman shift of certain sublevels + * `key::Tuple{String,Real}`: sublevel name + * `value::Function(B::Real)`: Nonlinear term(s) of Zeeman shift. Full Zeeman shift will be calculated as the sum of the usual linear term and this function +""" +Base.@kwdef struct IonProperties # Apparently Parameters.jl is preferred for this kind of thing. + shortname::String + mass::Number + charge::Number + nuclearspin::Number + full_level_structure::OrderedDict{ + String, + NamedTuple{(:n, :l, :j, :f, :E), T} where T <: Tuple + } + full_transitions::Dict{ + Tuple{String, String}, + NamedTuple{(:multipole, :einsteinA), Tuple{String, Float64}} + } + default_sublevel_selection::Union{Vector{Tuple{String, String}}, Missing} = missing + gfactors::Union{Dict{String, Number}, Missing} = missing + nonlinear_zeeman::Union{Dict{Tuple{String, Real}, Function}, Missing} = missing +end + + +""" +Helper function to convert the ``full_level_structure``, as defined in the +config, to the format expected by ``IonProperties`` +""" +function process_full_level_structure(from_config::OrderedDict)::OrderedDict + to_return = OrderedDict{String, NamedTuple{(:n, :l, :j, :f, :E)}}() + for (level_name, level_data) in from_config + to_return[level_name] = ( + n=level_data["n"], + l=level_data["l"], + j=roundnearesthalf(level_data["j"]), + f=roundnearesthalf(level_data["f"]), + E=level_data["E"] + ) + end + + return to_return +end + + +""" +Helper function to convert the ``full_transitions``, as defined in the config, +to the format expected by ``IonProperties``. +""" +function process_full_transitions(from_config::Vector)::Dict + full_transitions = + Dict{Tuple{String, String}, @NamedTuple{multipole::String, einsteinA::Real}}() + + for transition_data in from_config + key = (transition_data["from"], transition_data["to"]) + full_transitions[key] = + (multipole=transition_data["multipole"], einsteinA=transition_data["einsteinA"]) + end + + return full_transitions +end + + +""" +Helper function to convert the nonlinear zeeman entries defined in the config +to the format expected by ``IonProperties`` +""" +function process_nonlinear_zeeman(from_config::Vector)::Dict + nonlinear_zeemans = Dict{Tuple{String, Integer}, Function}() + for nonlinear_shift in from_config + # Because of the way Julia refreshes bindings with each iteration, there + # is no need to create closures. Wonder if this is actually efficient. + key = (nonlinear_shift["level"], nonlinear_shift["sublevel"]) + coeffs = nonlinear_shift["coeffs"] + nonlinear_zeemans[key] = B -> dot([B^pwr for pwr in 0:(length(coeffs)-1)], coeffs) + end + + return nonlinear_zeemans +end +process_nonlinear_zeeman(::Missing) = missing + + +"""Instantiate from config file""" +function loadfromconfig(config_filepath::String) + config = YAML.load_file(config_filepath, dicttype=OrderedDict{String, Any}) + return _createfromconfigobject(config) +end + + +function _createfromconfigobject(config::OrderedDict) + full_level_structure = process_full_level_structure(config["full_level_structure"]) + full_transitions = process_full_transitions(config["full_transitions"]) + + # Optional arguments. + if haskey(config, "default_sublevel_selection") + default_sublevel_selection = + [Tuple(selection) for selection in config["default_sublevel_selection"]] + else + default_sublevel_selection = missing + end + gfactors = get(config, "gfactors", missing) + nonlinear_zeeman = process_nonlinear_zeeman(get(config, "nonlinear_zeeman", missing)) + + return IonProperties( + shortname=config["shortname"], + mass=config["mass"], + charge=config["charge"], + nuclearspin=roundnearesthalf(config["nuclearspin"]), + full_level_structure=full_level_structure, + full_transitions=full_transitions, + default_sublevel_selection=default_sublevel_selection, + gfactors=gfactors, + nonlinear_zeeman=nonlinear_zeeman + ) + +end + + +DEFAULT_PROPERTIES = Dict() + + +# Define one-shot function to just populate the DEFAULT_PROPERTIES dictionary. +function _populate_properties() + config_dir = joinpath(dirname(@__DIR__), "configs", "ions") + for ion_fname::String in readdir(config_dir) + full_path = joinpath(config_dir, ion_fname) + config_yaml = YAML.load_file(full_path, dicttype=OrderedDict{String, Any}) + ion_name::String = config_yaml["name"] + DEFAULT_PROPERTIES[ion_name] = _createfromconfigobject(config_yaml) + end +end + + +_populate_properties() + + +end # Properties module diff --git a/src/ions.jl b/src/ions.jl index 9191cb7..63be749 100644 --- a/src/ions.jl +++ b/src/ions.jl @@ -1,9 +1,11 @@ +import YAML # >.> is this non-standard? using WignerSymbols: wigner3j, wigner6j -using LinearAlgebra: cross +using OrderedCollections +using LinearAlgebra: cross, dot using IonSim.PhysicalConstants +using IonSim.Properties: IonProperties export Ion, - IonProperties, speciesproperties, sublevels, sublevelaliases, @@ -33,19 +35,85 @@ export Ion, einsteinA, transitionmultipole, lifetime, - matrixelement, - IonInstance + matrixelement ############################################################################################# # Ion - the physical parameters defining an ion's structure ############################################################################################# - """ - Ion +Ion(selected_sublevels::Vector{Tuple}[, manualshift::Dict]) + +`selected_sublevels` specifies which energy sublevels will be present in the Hilbert space of this Ion instance, as a subset of all possible sublevels. + +Each element of `selected_sublevels` is a 2-element Tuple `(level::String, sublevels)`, with the first element being the name of a level and the second specifying which sublevels should be included. +Allowed sublevels are those whose magnetic quantum number `m` is in the set {`-f`, `-f+1`, `-f+2`, ... `f-1`, `f`}, where `f` is the total angular momentum quantum number of `level`. +For each `level` specified there are three allowed options to specify the set of `sublevels` to include: +* `sublevels::Real`: Includes only one `m = sublevels` +* `sublevels::Vector{Real}`: Includes all sublevels whose magnetic quantum number `m` is in `sublevels` +* `sublevels = "all"`: Includes all allowed sublevels +If an element of `selected_sublevels` utilizes the first option, specifying a single `m`, one may optionally may this a 3-element tuple instead: `(level::String, m::Real, alias::String)`, assinging this particular sublevel the alias `alias`. -The physical parameters defining an isolated ion's internal structure. +Omission of a level in `selected_sublevels` will exclude all sublevels. + +**Fields** +* `speciesproperties::IonProperties`: Contains constants specifying parameters specific to species +* `sublevels`::Vector{Tuple{String,Real}}: List of all sublevels present in the Hilbert space +* `sublevelaliases::Dict{String,Tuple}`: Dict specifying aliases assigned to sublevels, in the format `alias => sublevel` +* `shape`::Vector{Int}: Dimension of the Hilbert space +* `manualshift::OrderedDict`: A dictionary with keys denoting the selected levels and values, a real number for describing a shift of the level's energy. This is just a convenient way to add manual shifts to the simulation, such as Stark shifts off of energy levels not present in the Hilbert space, without additional resources +* `ionnumber`: When the ion is added to an `IonTrap`, this value keeps track of its order +* `ionposition`: When the ion is added to an `IonTrap`, this value keeps track of its physical position in meters """ -abstract type Ion <: IonSimBasis end +mutable struct Ion <: IonSimBasis + # fields + speciesproperties::IonProperties + sublevels::Vector{Tuple{String, Real}} + sublevelaliases::Dict{String, Tuple} + shape::Vector{Int} + manualshift::OrderedDict{Tuple, Real} + ionnumber::Union{Int, Missing} + ionposition::Union{Real, Missing} + + # constructors (overrides default) + function Ion(properties, selected_sublevels=nothing, manualshift=Dict()) + (sublevels, aliases) = _construct_sublevels(selected_sublevels, properties) + shape = [length(sublevels)] + manualshift_full = _construct_manualshift(manualshift, sublevels) + return new( + properties, + sublevels, + aliases, + shape, + manualshift_full, + missing, + missing + ) + end + function Ion( + speciesproperties, + sublevels, + sublevelaliases, + shape, + manualshift, + ionnumber, + ionposition + ) + sublevels = deepcopy(sublevels) + sublevelaliases = deepcopy(sublevelaliases) + shape = copy(shape) + manualshift = deepcopy(manualshift) + return new( + speciesproperties, + sublevels, + sublevelaliases, + shape, + manualshift, + ionnumber, + ionposition + ) + end +end + ############################################################################################# # Object fields @@ -232,14 +300,12 @@ end """ levels(I::Ion) - Returns array of all energy levels of `I`. """ levels(I::Ion) = unique([sublevel[1] for sublevel in sublevels(I)]) """ sublevelalias(I::Ion, sublevel::Tuple{String,Real}) - Returns the alias assined to `sublevel` of `I`. If no alias is assigned, returns `nothing`. """ function sublevelalias(I::Ion, sublevel::Tuple{String, Real}) @@ -279,17 +345,15 @@ function level(I::Ion, alias::String) return sl[1] end -#= -quantumnumbers is written to be able to accept either a level or sublevel in the second -argument. Since both levels and aliases are strings, multidispatach can't tell the difference, -so the second method distinguishes these cases with an if statement. -=# +# quantumnumbers is written to be able to accept either a level or sublevel in the second argument +# Since both levels and aliases are strings, multidispatach can't tell the difference, so the second method distinguishes these cases with an if statement. """ quantumnumbers(I::Ion, level::String) quantumnumbers(I::Ion, sublevel) - Returns the quantum numbers of an energy level or sublevel of `I` as a `NamedTuple`. + If second argument is a level, returns `(:n, :i, :s, :l, :j, :f)` + If second argument is a sublevel, returns `(:n, :i, :s, :l, :j, :f, :m)` """ function quantumnumbers(I::Ion, sublevel::Tuple{String, Real}) @@ -308,7 +372,7 @@ function quantumnumbers(I::Ion, sublevel::Tuple{String, Real}) return NamedTuple{names}(values) end function quantumnumbers(I::Ion, level_or_alias::String) - # If the second argument is a String, it's either a level name or the alias of a sublevel + # If the second argument is a String, it could be either a level name or the alias of a sublevel if level_or_alias in levels(I) # Second argument is a level name. Leave out the m quantum number levelstruct = speciesproperties(I).full_level_structure[level_or_alias] @@ -330,7 +394,6 @@ end """ landegj(l::Real, j::Real, s::Real=1//2) - Landé g-factor of fine structure energy level **args** @@ -343,7 +406,6 @@ landegj(l::Real, j::Real, s::Real=1 // 2) = """ landegf(l::Real, j::Real, f::Real, i::Real, s::Real=1//2) - Landé g-factor of hyperfine energy level **args** @@ -358,7 +420,6 @@ landegf(l::Real, j::Real, f::Real, i::Real, s::Real=1 // 2) = landegf(qnums::NamedTuple) = landegf(qnums.l, qnums.j, qnums.f, qnums.i, qnums.s) """ landegf(I::Ion, level::String) - `landegf` for the quantum numbers of `level` in `I`. """ function landegf(I::Ion, level::String) @@ -383,8 +444,9 @@ manualshift(I::Ion, alias::String) = manualshift(I, sublevel(I, alias)) """ zeemanshift(I::Ion, sublevel, B::Real) Returns the Zeeman shift at a magnetic field of `B` of `sublevel` of `I`. -If `sublevel` has a custom g-factor defined, then this is used. Otherwise, `landegf` is used -to compute the Landé g-factor. + +If `sublevel` has a custom g-factor defined, then this is used. Otherwise, `landegf` is used to compute the Landé g-factor. + Zeeman shift calculated as ``ΔE = (μ_B/ħ) ⋅ g_f ⋅ B ⋅ m / 2π`` """ function zeemanshift(I::Ion, sublevel::Tuple{String, Real}, B::Real) @@ -405,11 +467,8 @@ zeemanshift(B::Real, qnums::NamedTuple) = zeemanshift(B, qnums.l, qnums.j, qnums.f, qnums.m, qnums.i, qnums.s) zeemanshift(I::Ion, alias::String, B::Real) = zeemanshift(I, sublevel(I, alias), B) -#= -This function is written to be able to accept either a level or sublevel in the second -argument. Since both levels and aliases are strings, multidispatach can't tell the difference, -so the second method distinguishes these cases with an if statement. -=# +# This function is written to be able to accept either a level or sublevel in the second argument +# Since both levels and aliases are strings, multidispatach can't tell the difference, so the second method distinguishes these cases with an if statement. """ energy(I::Ion, sublevel; B=0, ignore_manualshift=false) Returns energy of `sublevel` of `I`. A Zeeman shift may be included by setting the value of the magnetic field `B`. The manual shift may be omitted by setting `ignore_manualshift=true`. @@ -423,7 +482,6 @@ function energy(I::Ion, sublevel::Tuple{String, Real}; B=0, ignore_manualshift=f end """ energy(I::Ion, level::String) - Returns the energy of `level` of `I`. """ function energy(I::Ion, level_or_alias::String; B=0, ignore_manualshift=false) @@ -469,7 +527,6 @@ end """ leveltransitions(I::Ion) - Returns all allowed transitions between levels of `I` as a vector of `Tuple{String,String}`. """ function leveltransitions(I::Ion) @@ -485,9 +542,7 @@ end """ subleveltransitions(I::Ion) - -Returns all allowed transitions between sublevels of `I` as a vector of `Tuple{S,S}` where -`S=Tuple{String,Real}`. +Returns all allowed transitions between sublevels of `I` as a vector of `Tuple{S,S}` where `S=Tuple{String,Real}`. """ function subleveltransitions(I::Ion) list = [] @@ -501,9 +556,7 @@ function subleveltransitions(I::Ion) m1 = sl1[2] m2 = sl2[2] if abs(m2 - m1) <= parse(Int, multipole[2]) - #= Only add to list of sublevel transitions if Δm is not larger than the - transition multipole allows (1 for E1, 2 for E2, etc) - =# + # Only add to list of sublevel transitions if Δm is not larger than the transition multipole allows (1 for E1, 2 for E2, etc) push!(list, (sl1, sl2)) end end @@ -534,14 +587,12 @@ end """ lifetime(I::Ion, level::String) +Computes lifetime of `level` by summing the transition rates out of `level`. -Computes lifetime of `level` by summing the transition rates out of `level`. The sum is taken - over all levels that the species may have, rather than the levels present in the instance `I`. +The sum is taken over all levels that the species may have, rather than the levels present in the instance `I`. """ function lifetime(I::Ion, level::String) - @assert level in keys(speciesproperties(I).full_level_structure) ( - "Ion species $(typeof(I)) does not contain level $level" - ) + @assert level in keys(speciesproperties(I).full_level_structure) "Ion species $(typeof(I)) does not contain level $level" totaltransitionrate = 0.0 for (transition, info) in speciesproperties(I).full_transitions if transition[2] == level @@ -599,8 +650,7 @@ function matrixelement( a = cross(Bhat_array, [0, 0, 1]) / norm(cross(Bhat_array, [0, 0, 1])) theta = acos(Bhat_array[3]) amatrix = [0 -a[3] a[2]; a[3] 0 -a[1]; -a[2] a[1] 0] - # Rotation matrix in axis-angle representation (axis=a, angle=theta) - R = eye3 + sin(theta) * amatrix + (1 - cos(theta)) * amatrix^2 + R = eye3 + sin(theta) * amatrix + (1 - cos(theta)) * amatrix^2 # Rotation matrix in axis-angle representation (axis=a, angle=theta) end ϵhat_rotated = R * ϵhat_array khat_rotated = R * khat_array @@ -636,10 +686,7 @@ function matrixelement( return units_factor * hyperfine_factor * geometric_factor / 2π end else - @error ( - "calculation of atomic transition matrix element for transition type $multipole " * - "not currently supported" - ) + @error "calculation of atomic transition matrix element for transition type $multipole not currently supported" end end function matrixelement( @@ -693,10 +740,7 @@ function _construct_sublevels(selected_sublevels, properties) if !ismissing(properties.default_sublevel_selection) selected_sublevels = properties.default_sublevel_selection else - @error ( - "no level structure specified in constructor, and no default level structure " * - "specified for this ion species" - ) + @error "no level structure specified in constructor, and no default level structure specified for this ion species" end elseif selected_sublevels == "all" selected_sublevels = [(sublevel, "all") for sublevel in keys(full_level_structure)] @@ -732,9 +776,7 @@ function _construct_sublevels(selected_sublevels, properties) end for m in selectedms m = Rational(m) - @assert m in m_allowed ( - "Zeeman sublevel m = $m not allowed for state $level with f = $f" - ) + @assert m in m_allowed "Zeeman sublevel m = $m not allowed for state $level with f = $f" @assert (level, m) ∉ sublevels "repeated instance of sublevel $m in state $level" push!(sublevels, (level, m)) end @@ -768,71 +810,6 @@ function Base.:(==)(b1::T, b2::T) where {T <: Ion} ) end -""" -IonProperties type. - -**Required keywords** -* `mass`: Mass of ion in kg -* `charge`: Charge of ion in units of elementary charge -* `full_level_structure`: OrderedDict describing properties of each energy level - * `key::String`: Name of energy level. Spectroscopic notation encouraged, e.g. `"S1/2,f=1"` - * `value::NamedTuple(:n, :l, :j, :f, :E)`: Quantum numbers `n`, `l`, `j`, `f`, and energy - `E` (in Hz) -* `full_transitions`: Dict of all allowed transitions between energy levels - * `key::Tuple{String,String}` Pair of levels, ordered (lower, upper) in energy - * `value::NamedTuple(:multipole, :einsteinA)`: Leading-order multipole of the transition - (e.g. `"E1"`, `"E2"`) and Einstein A coefficient (between fine structure levels only; - hyperfine factors are calculated when needed) - - **Optional keywords** - * `default_sublevel_selection`: Default value of `selected_sublevels` argument in Ion - constructor - * `gfactors`: `Dict(level::String => g::Real)` Custom Landé g-factors, if contributions from - higher-than-first-order perturbations are desired - * `nonlinear_zeeman`: `Dict` describing nonlinear contributions to Zeeman shift of certain - sublevels - * `key::Tuple{String,Real}`: sublevel name - * `value::Function(B::Real)`: Nonlinear term(s) of Zeeman shift. Full Zeeman shift will be - calculated as the sum of the usual linear term and this function -""" -struct IonProperties - shortname::String - mass::Number - charge::Number - nuclearspin::Number - full_level_structure::OrderedDict{ - String, - NamedTuple{(:n, :l, :j, :f, :E), T} where T <: Tuple - } - full_transitions::Dict{ - Tuple{String, String}, - NamedTuple{(:multipole, :einsteinA), Tuple{String, Float64}} - } - default_sublevel_selection::Union{Vector{Tuple{String, String}}, Missing} - gfactors::Union{Dict{String, Number}, Missing} - nonlinear_zeeman::Union{Dict{Tuple{String, Real}, Function}, Missing} -end -IonProperties(; - shortname, - mass, - charge, - nuclearspin, - full_level_structure, - full_transitions, - default_sublevel_selection=missing, - gfactors=missing, - nonlinear_zeeman=missing -) = IonProperties( - shortname, - mass, - charge, - nuclearspin, - full_level_structure, - full_transitions, - default_sublevel_selection, - gfactors, - nonlinear_zeeman -) function Base.:(==)(p1::IonProperties, p2::IonProperties) for fn in fieldnames(IonProperties) @@ -849,87 +826,8 @@ function Base.:(==)(p1::IonProperties, p2::IonProperties) return true end -""" -IonInstance(selected_sublevels::Vector{Tuple}[, manualshift::Dict]) -Ion instance of some species - -`selected_sublevels` specifies which energy sublevels will be present in the Hilbert space of -this Ion instance, as a subset of all possible sublevels. - -Each element of `selected_sublevels` is a 2-element Tuple `(level::String, sublevels)`, with the first element being the name of a level and the second specifying which sublevels should be included. -Allowed sublevels are those whose magnetic quantum number `m` is in the set {`-f`, `-f+1`, `-f+2`, ... `f-1`, `f`}, where `f` is the total angular momentum quantum number of `level`. -For each `level` specified there are three allowed options to specify the set of `sublevels` to include: -* `sublevels::Real`: Includes only one `m = sublevels` -* `sublevels::Vector{Real}`: Includes all sublevels whose magnetic quantum number `m` is in `sublevels` -* `sublevels = "all"`: Includes all allowed sublevels -If an element of `selected_sublevels` utilizes the first option, specifying a single `m`, one may optionally may this a 3-element tuple instead: `(level::String, m::Real, alias::String)`, assinging this particular sublevel the alias `alias`. - -Omission of a level in `selected_sublevels` will exclude all sublevels. - -**Fields** -* `speciesproperties::IonProperties`: Contains constants specifying parameters specific to species -* `sublevels`::Vector{Tuple{String,Real}}: List of all sublevels present in the Hilbert space -* `sublevelaliases::Dict{String,Tuple}`: Dict specifying aliases assigned to sublevels, in the format `alias => sublevel` -* `shape`::Vector{Int}: Dimension of the Hilbert space -* `manualshift::OrderedDict`: A dictionary with keys denoting the selected levels and values, a real number for describing a shift of the level's energy. This is just a convenient way to add manual shifts to the simulation, such as Stark shifts off of energy levels not present in the Hilbert space, without additional resources -* `ionnumber`: When the ion is added to an `IonTrap`, this value keeps track of its order -* `ionposition`: When the ion is added to an `IonTrap`, this value keeps track of its physical position in meters -""" -mutable struct IonInstance{Species <: Any} <: Ion - # fields - speciesproperties::IonProperties - sublevels::Vector{Tuple{String, Real}} - sublevelaliases::Dict{String, Tuple} - shape::Vector{Int} - manualshift::OrderedDict{Tuple, Real} - ionnumber::Union{Int, Missing} - ionposition::Union{Real, Missing} - - # constructors (overrides default) - function IonInstance{Species}( - properties, - selected_sublevels=nothing, - manualshift=Dict() - ) where {Species <: Any} - (sublevels, aliases) = _construct_sublevels(selected_sublevels, properties) - shape = [length(sublevels)] - manualshift_full = _construct_manualshift(manualshift, sublevels) - return new{Species}( - properties, - sublevels, - aliases, - shape, - manualshift_full, - missing, - missing - ) - end - function IonInstance{Species}( - speciesproperties, - sublevels, - sublevelaliases, - shape, - manualshift, - ionnumber, - ionposition - ) where {Species <: Any} - sublevels = deepcopy(sublevels) - sublevelaliases = deepcopy(sublevelaliases) - shape = copy(shape) - manualshift = deepcopy(manualshift) - return new{Species}( - speciesproperties, - sublevels, - sublevelaliases, - shape, - manualshift, - ionnumber, - ionposition - ) - end -end -function Base.print(ion::IonInstance) +function Base.print(ion::Ion) println(ion.speciesproperties.shortname) for sublevel in sublevels(ion) if sublevel in values(sublevelaliases(ion)) @@ -941,12 +839,4 @@ function Base.print(ion::IonInstance) end end -Base.show(io::IO, I::IonInstance) = print(io, I.speciesproperties.shortname) # suppress long output - -function Base.:(==)(C1::Ion, C2::Ion) - if (mass(C1) == mass(C2)) && (charg(C1) == charg(C2)) - return true - else - return false - end -end +Base.show(io::IO, I::Ion) = println(io, I.speciesproperties.shortname) # suppress long output diff --git a/src/iontraps.jl b/src/iontraps.jl index 95a65a5..e8ef279 100644 --- a/src/iontraps.jl +++ b/src/iontraps.jl @@ -1,15 +1,13 @@ using LinearAlgebra: eigen using NLsolve: nlsolve -using Statistics: mean -using Optim using .PhysicalConstants: e, ϵ₀ -using Plots -import YAML export IonTrap, ions, linear_equilibrium_positions, + Anm, LinearChain, + characteristic_length_scale, comfrequencies, selectedmodes, full_normal_mode_description, @@ -17,16 +15,12 @@ export IonTrap, xmodes, ymodes, zmodes, - modecutoff!, - length, - ionpositions, - LinearChain_fromyaml, - visualize + modecutoff! """ IonTrap -A physical configuration of ions. Stores a collection of ions and information about their -relative center of mass motion. +Physical configuration of ions. Stores a collection of ions and information about the +interactions of their center of mass motion. """ abstract type IonTrap end @@ -38,338 +32,63 @@ ions(I::IonTrap) = I.ions ############################################################################################# """ - LinearChain(; - ions, comfrequencies::NamedTuple{(:x,:y,:z)}, - selectedmodes::NamedTuple{(:x,:y,:z), N::Int=10 - ) - -Generates and stores all of the information necessary to describe a collection of ions trapped -in a 3D harmonic potential and forming a linear coulomb crystal. Relevant to calculating the -normal mode structure of the linear chain is the charge, mass and ordering of the ions. - -**user-defined fields** -* `ions`: An iterable list of ions that compose the linear Coulomb crystal -* `comfrequencies::NamedTuple{(:x,:y,:z)`: Describes the COM frequencies - `(x=ν_x, y=ν_y, z=ν_z)`. The z-axis is taken to be along the crystal's axis of - symmetry. Note that these only correspond to true COM modes if all ions have the same - mass, otherwise this is a misnomer and corresponds to the smallest (largest) eigenfrequency - mode in the axial (radial) directions. Note: we assume that `ν_x > ν_y > ν_z`. -* `selectedmodes::NamedTuple{(:x,:y,:z)}`: e.g. `selectedmodes=(x=[1], y=[1, 2:3], z=[:])`. - Specifies the axis and a list of integers which correspond to the ``i^{th}`` farthest - mode away from the COM (see note above about meaning of "COM") for that axis. For example, - `selectedmodes=(z=[2])` would specify the axial stretch mode. These are the modes that will - be modeled in the chain.Note: `selectedmodes=(x=[],y=[],z=[1])`, `selectedmodes=(y=[],z=[1])` - and `selectedmodes=(;z=[1])` are all acceptable and equivalent. -* `N::Int=10`: Optionally specify the Hilbert space dimension for each normal mode. -**derived fields** -* `ionpositions::Vector{<:Real}`: The relative positions of the ions in meters. -""" -struct LinearChain <: IonTrap - ions::Tuple{Vararg{Ion}} - comfrequencies::NamedTuple{(:x, :y, :z)} - selectedmodes::NamedTuple{(:x, :y, :z), Tuple{Vararg{Vector{VibrationalMode}, 3}}} - ionpositions::Vector{<:Real} - - function LinearChain(; - ions, - comfrequencies::NamedTuple, - selectedmodes::NamedTuple, - N::Int=10 - ) - try - ions = collect(ions) - catch - AssertionError("`ions` must be an iterable collection of `<:Ion`.") |> throw - end - num_ions = length(ions) - selectedmodes = _construct_vibrational_modes(selectedmodes, num_ions) - - # shouldn't need this after - warn = nothing - for i in 1:num_ions, j in (i+1):num_ions - if ions[j] ≡ ions[i] - ions[j] = copy(ions[i]) - if isnothing(warn) - warn = "Some ions point to the same thing. Making copies." - @warn warn - end - end - end - - normalmodes = full_normal_mode_description(ions, comfrequencies) - vm = ( - x=Vector{VibrationalMode}(undef, 0), - y=Vector{VibrationalMode}(undef, 0), - z=Vector{VibrationalMode}(undef, 0) - ) - axes = [x̂, ŷ, ẑ] - for (i, modes) in enumerate(selectedmodes), mode in modes - push!(vm[i], VibrationalMode(normalmodes[i][mode]..., axis=axes[i], N=N)) - end - l = linear_equilibrium_positions(num_ions) - kz = searchfor_trappingpotential_parameters(ions, comfrequencies)[1] - l0 = (e^2 / (4π * ϵ₀ * kz))^(1 / 3) - ionpositions = l * l0 - _setionproperties(ions, ionpositions) - return new(ions |> Tuple, comfrequencies, vm, ionpositions) - end - - function LinearChain(ions, vm, ionpositions) - _setionproperties(ions, ionpositions) - new(ions |> Tuple, (x=NaN, y=NaN, z=NaN), vm, ionpositions) - end -end - -""" - LinearChain_fromyaml(ions::Vector{<:Ion}, yaml::String; N::Int=10) - -Load normal mode structure from the specified `yaml` file. It's up to the user to enforce -physicality. -```` -x: - - frequency: 1e6 - mode: [0.1, 0.5, 0.3, 0.8] - - frequency: 2e6 - mode: [0.3, 0.6, 0.5, 3] -y: - - frequency: 8e6 - mode: [1, 1, 1, 1] -ionpositions: [-1, -0.5, -0.25, 5] -```` -It is up to the user to enforce the physicality of this structure. + linear_equilibrium_positions(N::Int) +Returns the scaled equilibrium positions of `N` ions in a harmonic potential, assuming that +all ions have the same mass. +[ref](https://doi.org/10.1007/s003400050373) """ -function LinearChain_fromyaml(; ions::Vector{<:Ion}, yaml::String, N::Int=10) - num_ions = length(ions) - axes = ["x", "y", "z"] - axes_dict = Dict("x" => x̂, "y" => ŷ, "z" => ẑ) - # comfrequencies = (x=NaN, y=NaN, z=NaN) - yaml = YAML.load_file(yaml) - k = keys(yaml) - @assert "ionpositions" in k ("yml file must have a key 'ionpositions'") - ionpositions = yaml["ionpositions"] - @assert length(ionpositions) == num_ions ( - "length(ionpositions) must equal length(ions)" - ) - vm = Dict("x" => [], "y" => [], "z" => []) - for axis in axes - if axis in k - axismodes = yaml[axis] - for mode in axismodes - if !("frequency" in keys(mode)) - @warn "no frequency specified for a mode in $axis" - continue - elseif !("mode" in keys(mode)) - @warn "no mode specified for frequency $(mode["frequency"])" - continue - elseif !(length(mode["mode"]) == num_ions) - @warn "length(mode) ≠ length(ions) for $(mode["frequency"]), $(mode["mode"])" - continue - end - push!(vm[axis], (mode["frequency"], mode["mode"])) - end - end - end - modes = (x=vm["x"], y=vm["y"], z=vm["z"]) - vm = ( - x=Vector{VibrationalMode}(undef, 0), - y=Vector{VibrationalMode}(undef, 0), - z=Vector{VibrationalMode}(undef, 0) - ) - for (i, axis) in enumerate(axes), mode in modes[Symbol(axis)] - push!(vm[i], VibrationalMode(mode..., axis=axes_dict[axis], N=N)) - end - println(vm) - LinearChain(ions, vm, ionpositions) - # LinearChain( - # ions=ions, comfrequencies=comfrequencies, selectedmodes=modes, - # ionpositions=ionpositions - # ) -end - -function _setionproperties(ions, ionpositions) - for (i, ion) in enumerate(ions) - Core.setproperty!(ion, :ionnumber, i) - Core.setproperty!(ion, :ionposition, ionpositions[i]) - end -end -############################################################################################# -# LinearChain - helper functions -############################################################################################# - -_sparsify!(x, eps) = @. x[abs(x) 0)) ( - "Can't have negative and positive charges in the same chain." - ) +function linear_equilibrium_positions(N::Int) function f!(F, x, N) for i in 1:N F[i] = ( - Q[i] * x[i] - sum([Q[i] * Q[j] / (x[i] - x[j])^2 for j in 1:(i-1)]) + sum([Q[i] * Q[j] / (x[i] - x[j])^2 for j in (i+1):N]) + x[i] - sum([1 / (x[i] - x[j])^2 for j in 1:(i-1)]) + sum([1 / (x[i] - x[j])^2 for j in (i+1):N]) ) end end - # equilibrium positions follow an empirical relationship used for picking starting point - # (see ref) + function j!(J, x, N) + for i in 1:N, j in 1:N + if i ≡ j + J[i, j] = ( + 1 + 2 * sum([1 / (x[i] - x[j])^3 for j in 1:(i-1)]) - 2 * sum([1 / (x[i] - x[j])^3 for j in (i+1):N]) + ) + else + J[i, j] = ( + -2 * sum([1 / (x[i] - x[j])^3 for j in 1:(i-1)]) + 2 * sum([1 / (x[i] - x[j])^3 for j in (i+1):N]) + ) + end + end + end + # see eq.8 in the ref to see where (2.018/N^0.559) comes from if isodd(N) initial_x = [(2.018 / N^0.559) * i for i in (-N÷2):(N÷2)] else initial_x = [(2.018 / N^0.559) * i for i in filter(x -> x ≠ 0, collect((-N÷2):(N÷2)))] end - sol = nlsolve((F, x) -> f!(F, x, N), initial_x, method=:newton) - # clean up values - sol = round.(sol.zero, digits=6) - for i in eachindex(sol) # remove annoying signed zeros - if sol[i] == 0 - sol[i] = 0 - end - end - return sol + sol = nlsolve((F, x) -> f!(F, x, N), (J, x) -> j!(J, x, N), initial_x, method=:newton) + return sol.zero end -linear_equilibrium_positions(ions::Vector{<:Ion}) = linear_equilibrium_positions( - length(ions), - [ion.speciesproperties.charge for ion in ions] -) - -linear_equilibrium_positions(N::Int) = linear_equilibrium_positions(N, ones(Int, N)) - -#= -This is the Jacobian for the force of the ions in the chain due to all of the other ions. -Diagonalize it to find the normal modes of motion. -=# -function diagonalize_Kij( - M_actual::Vector{<:Real}, - Q::Vector{<:Int}, - axis::NamedTuple{(:x, :y, :z)}, - kz::Real, - P::Real, - kr::Real; - optimize=false -) - N = length(M_actual) - maxM = maximum(M_actual) - M = M_actual / maxM # scale M_actual, so masses are order unity - # P *= maxM # scale P by Mtot since it gets divided by M[i] - - #= - kz = ka - kr/2 - kxⱼ = P/Mⱼ - ka/2 + kr = -kz/2 + 3kr/4 - kyⱼ = P/Mⱼ - ka/2 - kr/2 = -kz/2 - 3kr/4 - - Where P is defined s.t. the pseudopotential frequency for jth ion is (√(Qⱼ * P)/Mⱼ) * kz. - And ka is defined s.t. the axial DC confining potential frequency is √(Qⱼ * ka / Mⱼ) * kz. - (Similar for kr, but this is an additional DC quadrupole that is confining in radial direction). - =# - if axis == x̂ - a = 1 - kDC = -kz / 2 + 3kr / 4 # kx - elseif axis == ŷ - a = 1 - kDC = -kz / 2 - 3kr / 4 # ky - elseif axis == ẑ - a = -2 - P = 0 - kDC = kz # kz - end - - l = linear_equilibrium_positions(N, Q) - K = Array{Real}(undef, N, N) - for i in 1:N, j in 1:N - if i ≡ j - K[i, j] = - Q[i] * (P / M[i] + kDC) / kz - - a * sum([Q[i] * Q[p] / abs(l[j] - l[p])^3 for p in 1:N if p != j]) - else - K[i, j] = a * Q[i] * Q[j] / abs(l[j] - l[i])^3 - end - # We switch to mass-scaled coordinates as is standard for linear spring-mass systems. - K[i, j] /= √(M[i] * M[j]) - end - - bvalues, bvectors = eigen(K) - bcleanedvalues = [] - # When arg for sqrt is negative, we still want to return something real so the optimizer will work. - # But we need to make sure it returns a *definitely incorrect* but not super divergent value. - # This method seems to work, but seems a bit sketchy. - for value in bvalues - if kz > 0 && value > 0 - push!(bcleanedvalues, sqrt.(value * kz / maxM)) - # push!(bcleanedvalues, sqrt.(value * kz)) - else - # we make it negative so it is *definitely incorrect* - push!(bcleanedvalues, -sqrt.(abs(value * kz / maxM))) - # push!(bcleanedvalues, -sqrt.(abs(value * kz))) - end - end - - # When optimizing, we are just trying to match the computed eigenfrequencies with the target - # eigenfrequencies, so no need to compute the eigenvectors - if optimize - if axis == ẑ - return minimum(bcleanedvalues) - else - return maximum(bcleanedvalues) - end - end - - for i in 1:N - v = view(bvectors, :, i) - _sparsify!(v, 1e-6) - for j in 1:N - v ./= M[j] - end - # normalization takes care of any remaining global scalings introduced for convenience - normalize!(v) - v .*= sign(v[1]) - end - - sortedresults = - sort([(bcleanedvalues[i], bvectors[:, i]) for i in 1:length(bcleanedvalues)]) - if axis == ẑ - # axial modes are sorted from lowest to highest eigenfrequency - return sortedresults - else - # this ordering is reversed for radial modes - return reverse(sortedresults) - end -end +""" + characteristic_length_scale(M::Real, ν::Real) +Returns the characteristic length scale for a linear chain of identical ions of mass `M` +and with axial trap frequency 2π × ν. +""" +characteristic_length_scale(M::Real, ν::Real) = (e^2 / (4π * ϵ₀ * M * (2π * ν)^2))^(1 / 3) -diagonalize_Kij( - ions::Vector{<:Ion}, - axis::NamedTuple{(:x, :y, :z)}, - kz::Real, - P::Real, - kr::Real; - optimize=false -) = diagonalize_Kij( - [mass(ion) for ion in ions], - [ion.speciesproperties.charge for ion in ions], - axis, - kz, - P, - kr, - optimize=optimize -) +""" + Anm(N::Real, com::NamedTuple{(:x,:y,:z)}, axis::NamedTuple{(:x,:y,:z)}) +Computes the normal modes and corresponding trap frequencies along a particular `axis` for a +collection of `N` ions in a linear Coloumb crystal and returns an array of tuples with first +element the frequency of the normal mode and 2nd element the corresponding eigenvector. -#= -This performs the same function as diagonalize_Kij, but for a homogeneous chain where -all of the ions have the same mass. In that case, one can write the Jacobian Kij in -terms of the characteristic_frequencies, which, in this case are the uniqueCOM frequencies. -Then we can just diagonalize Kij directly instead of having to optimize. -=# -function diagonalize_Kij_for_homogeneous_chain( - N::Int, - com::NamedTuple{(:x, :y, :z)}, - axis::NamedTuple{(:x, :y, :z)} -) +`com` should be a `NamedTuple` of COM frequences for the different axes: +`(x<:Real, y<:Real, z<:Real)`, where the ``z``-axis is taken to be parallel to the axis of +the crystal. +""" +function Anm(N::Int, com::NamedTuple{(:x, :y, :z)}, axis::NamedTuple{(:x, :y, :z)}) + @assert axis in [x̂, ŷ, ẑ] "axis can be x̂, ŷ, ẑ" axis == ẑ ? a = 2 : a = -1 l = linear_equilibrium_positions(N) axis_dict = Dict([(x̂, :x), (ŷ, :y), (ẑ, :z)]) @@ -392,9 +111,7 @@ function diagonalize_Kij_for_homogeneous_chain( a1 = [sqrt(i) * com.z for i in a.values] a2 = a.vectors for i in 1:size(a2, 2) - v = view(a2, :, i) - _sparsify!(v, 1e-6) - v .*= sign(v[1]) + _sparsify!(view(a2, :, i), 1e-5) end if axis == ẑ return [(a1[i], a2[:, i]) for i in 1:length(a1)] @@ -404,187 +121,81 @@ function diagonalize_Kij_for_homogeneous_chain( end end +_sparsify!(x, eps) = @. x[abs(x) :x, ŷ => :y) - computed_eigenfrequencies = [] - for i in eachindex(axes) - axis = axes[i] - computed_eigenfrequency = diagonalize_Kij(M, Q, axis, kz, P, kr, optimize=true) - push!(computed_eigenfrequencies, computed_eigenfrequency) - end - difference = 0 - for i in eachindex(computed_eigenfrequencies) - difference += - abs(computed_eigenfrequencies[i] - target_eigenfrequencies[axes_dict[axes[i]]]) - end - return difference -end +""" + LinearChain(; + ions::Vector{Ion}, + comfrequencies::NamedTuple{(:x,:y,:z)}, + selectedmodes::NamedTuple{(:x,:y,:z),Tuple{Vararg{Vector{VibrationalMode},3}}} + ) -sequentially_diagonalize( - ions::Vector{<:Ion}, - target_eigenfrequencies::NamedTuple{(:x, :y, :z)}, - kz::Real, - P::Real, - kr::Real, - just_y=false -) = sequentially_diagonalize( - [mass(ion) for ion in ions], - [ion.speciesproperties.charge for ion in ions], - target_eigenfrequencies, - kz, - P, - kr, - just_y=just_y -) +Contains all of the information necessary to describe a collection of ions trapped in a 3D +harmonic potential and forming a linear coulomb crystal. -function searchfor_trappingpotential_parameters( - M::Vector{<:Real}, - Q::Vector{<:Int}, - target_eigenfrequencies::NamedTuple{(:x, :y, :z)} -) - @assert length(M) == length(Q) "length(M) ≠ length(Q)" - maxM = maximum(M) - value = diagonalize_Kij(M, Q, ẑ, maxM, 0, 0; optimize=true) - kz = (target_eigenfrequencies.z / value)^2 * maxM - lowerbounds = [-1, -1] - upperbounds = [Inf, Inf] - difference = 0 - if target_eigenfrequencies.x != target_eigenfrequencies.y - initial = [10kz, kz / 2] - res = optimize( - x -> - sequentially_diagonalize(M, Q, target_eigenfrequencies, kz, x[1], x[2]), - lowerbounds, - upperbounds, - initial, - Fminbox(NelderMead()), - Optim.Options(g_abstol=1e-12, g_reltol=1e-6), +**user-defined fields** +* `ions::Vector{Ion}`: a list of ions that compose the linear Coulomb crystal +* `comfrequencies::NamedTuple{(:x,:y,:z),Tuple{Vararg{Vector{VibrationalMode},3}}}`: + Describes the COM frequencies `(x=ν_x, y=ν_y, z=ν_z)`. The ``z``-axis is taken to be + parallel to the crystal's symmetry axis. +* `selectedmodes::NamedTuple{(:x,:y,:z)}`: e.g. `selectedmodes=(x=[1], y=[2], z=[1,2])`. + Specifies the axis and a list of integers which correspond to the ``i^{th}`` farthest + mode away from the COM for that axis. For example, `selectedmodes=(z=[2])` would + specify the axial stretch mode. These are the modes that will be modeled in the chain. + Note: `selectedmodes=(x=[],y=[],z=[1])`, `selectedmodes=(y=[],z=[1])` + and `selectedmodes=(;z=[1])` are all acceptable and equivalent. +""" +struct LinearChain <: IonTrap # Note: this is not a mutable struct + ions::Vector{<:Ion} + comfrequencies::NamedTuple{(:x, :y, :z)} + selectedmodes::NamedTuple{(:x, :y, :z), Tuple{Vararg{Vector{VibrationalMode}, 3}}} + function LinearChain(; + ions, + comfrequencies::NamedTuple, + selectedmodes::NamedTuple, + N=10 + ) + selectedmodes = _construct_vibrational_modes(selectedmodes) + warn = nothing + for i in 1:length(ions), j in (i+1):length(ions) + @assert typeof(ions[i]) == typeof(ions[j]) "multispecies chains not yet supported; all ions in chain must be of same species" + if ions[j] ≡ ions[i] + ions[j] = copy(ions[i]) + if isnothing(warn) + warn = "Some ions point to the same thing. Making copies." + @warn warn + end + end + end + nions = length(ions) + A = ( + x=Anm(nions, comfrequencies, x̂), + y=Anm(nions, comfrequencies, ŷ), + z=Anm(nions, comfrequencies, ẑ) ) - xfreqs = [i[1] for i in diagonalize_Kij(M, Q, x̂, kz, res.minimizer...)] - yfreqs = [i[1] for i in diagonalize_Kij(M, Q, ŷ, kz, res.minimizer...)] - difference += abs(target_eigenfrequencies.x - xfreqs[1]) - difference += abs(target_eigenfrequencies.y - yfreqs[1]) - res = [kz, res.minimizer...] - else - initial = [10kz] - res = optimize( - x -> sequentially_diagonalize(M, Q, target_eigenfrequencies, kz, x[1], 0), - lowerbounds[1:1], - upperbounds[1:1], - initial[1:1], - Fminbox(NelderMead()), - Optim.Options(g_abstol=1e-6, g_reltol=1e-6), + vm = ( + x=Vector{VibrationalMode}(undef, 0), + y=Vector{VibrationalMode}(undef, 0), + z=Vector{VibrationalMode}(undef, 0) ) - xfreqs = [i[1] for i in diagonalize_Kij(M, Q, x̂, kz, res.minimizer[1], 0)] - yfreqs = xfreqs - difference += abs(target_eigenfrequencies.x - xfreqs[1]) - res = [kz, res.minimizer..., 0] + r = [x̂, ŷ, ẑ] + for (i, modes) in enumerate(selectedmodes), mode in modes + push!(vm[i], VibrationalMode(A[i][mode]..., axis=r[i], N=N)) + end + l = linear_equilibrium_positions(length(ions)) + l0 = characteristic_length_scale(mass(ions[1]), comfrequencies.z) # Needs to be changed when allowing for multi-species chains. Current workaround takes the mass of only the first ion to define the characteristic length scale. + for (i, ion) in enumerate(ions) + Core.setproperty!(ion, :ionnumber, i) + Core.setproperty!(ion, :ionposition, l[i] * l0) + end + return new(ions, comfrequencies, vm) end - - only_positive_eigenvalues = all(vcat(xfreqs, yfreqs) .> 0) - only_nonnegative_trapping_parameters = all(res .>= 0) - error_string = ( - "Solver failed to find a normal mode structure compatible with `characteristicfrequencies`." - ) - # make sure optimal values are actually close to target_eigenfrequencies - @assert difference < 1e-3 error_string * " [1]" - # make sure there are no negative eigenfrequencies - @assert only_positive_eigenvalues error_string * " [2]" - # enfore convention that all trapping potentials are nonnegative - @assert only_nonnegative_trapping_parameters error_string * " [3]" - # make sure kz > kr, otherwise there is not confinment in z-direction - @assert res[1] > res[3] error_string * " [4]" - - return res # This is kz, P, kr end -searchfor_trappingpotential_parameters( - ions::Vector{<:Ion}, - target_eigenfrequencies::NamedTuple{(:x, :y, :z)} -) = searchfor_trappingpotential_parameters( - [mass(ion) for ion in ions], - [ion.speciesproperties.charge for ion in ions], - target_eigenfrequencies -) - - ############################################################################################# # Object fields ############################################################################################# -""" - full_normal_mode_description(chain::LinearChain) - -For each axis, this contains an array of tuples where the first element is a vibrational -frequency [Hz] and the second element is a vector describing the participation of each ion -at that vibrational frequency (i.e. the normal mode eigenvector corresponding to that -eigenfrequency). -""" -full_normal_mode_description(lc::LinearChain) = - full_normal_mode_description(ions(lc), comfrequencies(lc)) - -""" - full_normal_mode_description(ions, comfreqs::NamedTuple{(:x, :y, :z)}) - -Same thing but with an iterable list of `ions` and `NamedTuple` of COM frequencies -explicitly given. -""" -function full_normal_mode_description(ions, comfreqs::NamedTuple{(:x, :y, :z)}) - try - ions = collect(ions) - catch - AssertionError("`ions` must be an iterable collection of `<:Ion`.") |> throw - end - @assert !isnan(comfreqs.x) "This function doesn't work for user defined mode structure." - trappingparams = searchfor_trappingpotential_parameters(ions, comfreqs) - normalmodes = ( - x=diagonalize_Kij(ions, x̂, trappingparams...), - y=diagonalize_Kij(ions, ŷ, trappingparams...), - z=diagonalize_Kij(ions, ẑ, trappingparams...) - ) - return normalmodes -end - -""" - function full_normal_mode_description( - M::Vector{<:Real}, Q::Vector{<:Int}, comfreqs::NamedTuple{(:x, :y, :z)} - ) -Same thing but explicitly provide the masses `M` and charges `Q` of the ions. -""" -function full_normal_mode_description( - M::Vector{<:Real}, - Q::Vector{<:Int}, - comfreqs::NamedTuple{(:x, :y, :z)} -) - @assert length(M) == length(Q) "length(M) ≠ length(Q)" - trappingparams = searchfor_trappingpotential_parameters(M, Q, comfreqs) - normalmodes = ( - x=diagonalize_Kij(M, Q, x̂, trappingparams...), - y=diagonalize_Kij(M, Q, ŷ, trappingparams...), - z=diagonalize_Kij(M, Q, ẑ, trappingparams...) - ) - return normalmodes -end - """ comfrequencies(chain::LinearChain) Returns `chain.comfrequencies` @@ -597,6 +208,23 @@ Returns `chain.selectedmodes` """ selectedmodes(chain::LinearChain) = chain.selectedmodes +""" + full_normal_mode_description(chain::LinearChain)::NamedTuple{(:x,:y,:z)} +For each axis, this contains an array of tuples where the first element is a vibrational +frequency [Hz] and the second element is a vector describing the corresponding normalized +normal mode structure. +""" +function full_normal_mode_description(chain::LinearChain) + nions = length(ions(chain)) + com_freqs = comfrequencies(chain) + A = ( + x=Anm(nions, com_freqs, x̂), + y=Anm(nions, com_freqs, ŷ), + z=Anm(nions, com_freqs, ẑ) + ) + return A +end + """ modes(lc::LinearChain) Returns an array of all of the selected `VibrationalModes` in the `LinearChain`. @@ -608,29 +236,23 @@ end """ xmodes(lc::LinearChain) -Returns an array of all of the selected `VibrationalModes` in the x-direction in the -`LinearChain`. +Returns an array of all of the selected `VibrationalModes` in the x-direction in the `LinearChain`. """ xmodes(lc::LinearChain) = selectedmodes(lc).x - """ ymodes(lc::LinearChain) -Returns an array of all of the selected `VibrationalModes` in the y-direction in the -`LinearChain`. +Returns an array of all of the selected `VibrationalModes` in the y-direction in the `LinearChain`. """ ymodes(lc::LinearChain) = selectedmodes(lc).y - """ zmodes(lc::LinearChain) -Returns an array of all of the selected `VibrationalModes` in the z-direction in the -`LinearChain`. +Returns an array of all of the selected `VibrationalModes` in the z-direction in the `LinearChain`. """ zmodes(lc::LinearChain) = selectedmodes(lc).z """ modecutoff!(lc::LinearChain, N::Int) -Sets the upper bound of the Hilbert space of all `VibrationalMode`s in `lc` to be the Fock -state `N`. +Sets the upper bound of the Hilbert space of all `VibrationalMode`s in `lc` to be the Fock state `N`. """ function modecutoff!(lc::LinearChain, N::Int) for mode in modes(lc) @@ -661,278 +283,18 @@ function Base.show(io::IO, lc::LinearChain) # suppress long output end # Takes e.g. (y=[1]) to (x=[], y=[1], z=[]) -function _construct_vibrational_modes(x, num_ions) +function _construct_vibrational_modes(x) k = collect(keys(x)) xyz = [:x, :y, :z] - @assert isnothing(findfirst(x -> x ∉ xyz, k)) ( - "keys of `selectedmodes` must be `:x`, `:y` or `:z`" - ) + @assert isnothing(findfirst(x -> x ∉ xyz, k)) "keys of `selectedmodes` must be `:x`, `:y` or `:z`" indxs = findall(x -> x ∉ k, xyz) values = [] for i in 1:3 if i in indxs push!(values, Int[]) else - el = x[xyz[i]] - if !isempty(el) - # allows slices like (y=[1, 2:5], ) - elscreened = [] - for k in eachindex(el) - if typeof(el[k]) == Colon - @assert length(el) == 1 "Improper indexing." - elscreened = 1:num_ions - else - push!(elscreened, el[k]) - end - end - flattenedx = reduce( - vcat, - [typeof(i) <: UnitRange ? collect(i) : i for i in elscreened] - ) - else - flattenedx = [] - end - push!(values, flattenedx) + push!(values, x[xyz[i]]) end end return (; zip(xyz, values)...) end - -Base.length(lc::LinearChain) = length(lc.ions) - -""" - ionpositions(chain::LinearChain) -Returns the positions of the ions in `chain` in meters. -""" -ionpositions(chain) = chain.ionpositions - - -############################################################################################# -# Visualization -############################################################################################# - -""" - visualize(lc::LinearChain, axis::String, modes::Vector{<:Int}, format="bars") - -Visualize the normal mode structure of a `lc` as a Plots.jl plot. `axis` can either be one -of "x", "y", "z" or a NamedTuple{(:x,:y,:z)}. `modes` is a vector of indices -selecting which modes should be included in the plot. Slicing is supported. `format` can be -either "bars" or "circles." - -Note that the indexing refers to the full normal mode description and not the subset -`selectedmodes(lc)`. -""" -function visualize(lc::LinearChain, axis, modes; format="bars", legend=true) - ion_list = ions(lc) - axes_dict = Dict("x" => :x, "y" => :y, "z" => :z, x̂ => :x, ŷ => :y, ẑ => :z) - direction = axes_dict[axis] - radialplane = direction in [:x, :y] ? true : false - fnm = full_normal_mode_description(ion_list, comfrequencies(lc))[direction] - num_ions = length(ion_list) - unique_ion_labels = unique([ion.speciesproperties.shortname for ion in ion_list]) - color_list = [] - for ion in ion_list - push!(color_list, findfirst(ion.speciesproperties.shortname .== unique_ion_labels)) - end - - indextypes = map(typeof, modes) - if Colon in indextypes - indices = [1:num_ions...] - elseif UnitRange{Int} in indextypes - indices = reduce(vcat, [typeof(i) <: UnitRange ? collect(i) : i for i in modes]) - else - indices = modes - end - modes = fnm[indices] - subplots = [] - legend_plot = plot() - for label in unique_ion_labels - plot!( - [], - [], - color_palette=palette(:tab10), - label=" " * label, - lw=10, - legend=:top, - legendcolumns=1, - grid=false, - showaxis=false, - legendfontsize=15, - foreground_color_legend=nothing - ) - end - for mode in modes - frequency = round(mode[1] / 1e6, digits=3) - v = mode[2] - if format == "bars" - push!( - subplots, - bar( - v, - xaxis=false, - title="$frequency MHz", - legend=false, - color=color_list, - color_palette=palette(:tab10), - label=false, - yticks=false, - showaxis=false - ) - ) - hline!([0], label=false, lc=:black) - elseif format == "circles" - v = map(x -> x == 0 ? NaN : x, v) - offset = 0.6 - v = map(x -> x + sign(x) * offset, v) - xpos = linear_equilibrium_positions(num_ions) - GR.setarrowsize(0.75) - if !radialplane - push!( - subplots, - quiver( - xpos, - zeros(num_ions), - quiver=(v / 3, zeros(num_ions)), - lw=2, - color=:black - ) - ) - else - push!( - subplots, - quiver( - xpos, - zeros(num_ions), - quiver=(zeros(num_ions), v), - lw=2, - color=:black - ) - ) - end - scatter!( - xpos, - zeros(length(xpos)), - markersize=15, - bg=:white, - showaxis=false, - grid=false, - markerstrokewidth=0, - ylim=radialplane ? [minimum([v..., 0] .- 0.5), maximum([v..., 0] .+ 0.5)] : - [-0.5, 0.5], - xlim=[xpos[1] - 0.5, xpos[end] + 0.5], - size=radialplane ? (500, 500) : (1200, 300), - title="$frequency MHz", - legend=false, - color_palette=palette(:tab10), - markercolor=color_list - ) - end - end - nrows = 2 - divbythree = length(modes) ÷ nrows - modthree = length(modes) % nrows - if divbythree == 0 - l = (modthree + 1, 1) - elseif modthree == 0 - if legend - l = @layout [grid(1, 1){0.1w} grid(nrows, divbythree)] - else - l = @layout [grid(nrows, divbythree)] - end - else - if legend - l = @layout [grid(1, 1){0.1w} grid(nrows, divbythree) grid(modthree, 1)] - else - l = @layout [grid(nrows, divbythree) grid(modthree, 1)] - end - end - if legend - p = plot(legend_plot, subplots..., layout=l) - else - p = plot(subplots..., layout=l) - end - display(p) -end - -""" - visualize(vm::VibrationalMode; format="bars") - -Visualize the normal mode structure of a `vm` as a Plots.jl plot. `format` can be either -`bar` or `ions`. `format` can be either "bars" or "circles." -""" -function visualize(vm::VibrationalMode; format="bars") - axes_dict = Dict(x̂ => :x, ŷ => :y, ẑ => :z) - direction = axes_dict[vm.axis] - radialplane = direction in [:x, :y] ? true : false - frequency = round(vm.ν / 1e6, digits=3) - v = vm.modestructure - if format == "bars" - bar( - v, - xticks=(1:4, ones(length(v))), - xaxis=false, - ylabel="Relative participation\nof ion in mode", - title="$frequency MHz", - legend=false - ) - hline!([0], label=false, lc=:black) - elseif format == "circles" - v = map(x -> x == 0 ? NaN : x, v) - N = length(v) - xpos = linear_equilibrium_positions(N) - GR.setarrowsize(0.75) - scatter( - xpos, - zeros(length(xpos)), - markersize=15, - bg=:white, - showaxis=false, - grid=false, - markerstrokecolor=:white, - markercolor=:red3, - ylim=radialplane ? [minimum([v..., 0] .- 0.5), maximum([v..., 0] .+ 0.5)] : - [-0.5, 0.5], - xlim=[xpos[1] - 0.5, xpos[end] + 0.5], - size=(1200, 300), - title="$frequency MHz", - legend=false, - ) - offset_vector = copy(v) - for i in eachindex(offset_vector) - if offset_vector[i] > 0 - offset_vector[i] += 0.1 - else - offset_vector[i] -= 0.1 - end - end - if !radialplane - quiver!( - xpos, - zeros(N), - quiver=(offset_vector / 2, zeros(N)), - lw=2, - color=:black - ) - else - quiver!(xpos, zeros(N), quiver=(zeros(N), offset_vector), lw=3, color=:black) - end - end - display(current()) -end - -""" - visualize( - vm::Tuple{Float64, Vector{Float64}}, axis::NamedTuple{(:x, :y, :z)}; format="bars" - ) - -Same thing but input a normal mode description as a tuple with first element the eigenfrequency -and second the eigenvector. -""" -function visualize( - vm::Tuple{Float64, Vector{Float64}}, - axis::NamedTuple{(:x, :y, :z)}; - format="bars" -) - VM = VibrationalMode(vm..., axis=axis, N=1) - visualize(VM; format=format) -end diff --git a/src/species/_include_species.jl b/src/species/_include_species.jl deleted file mode 100644 index 51c6510..0000000 --- a/src/species/_include_species.jl +++ /dev/null @@ -1,4 +0,0 @@ -include("be9.jl") -include("ca40.jl") -include("mg25.jl") -include("yb171.jl") diff --git a/src/species/be9.jl b/src/species/be9.jl deleted file mode 100644 index 13e1cc9..0000000 --- a/src/species/be9.jl +++ /dev/null @@ -1,48 +0,0 @@ -export Be9 - -const properties_be9 = IonProperties( - shortname="⁹Be", - mass=1.496508080073e-26, - charge=1, - nuclearspin=3 // 2, - full_level_structure=OrderedDict( - "S1/2f=1" => (n=2, l=0, j=1 // 2, f=1, E=0.78126104631e9), - "S1/2f=2" => (n=2, l=0, j=1 // 2, f=2, E=-0.468756627786e9), - "P1/2f=1" => (n=2, l=1, j=1 // 2, f=1, E=957.2010729076436e12), - "P1/2f=2" => (n=2, l=1, j=1 // 2, f=2, E=957.2008357076436e12), - "P3/2f=0" => (n=2, l=1, j=3 // 2, f=0, E=957.3965669467407e12), - "P3/2f=1" => (n=2, l=1, j=3 // 2, f=1, E=957.3965659267407e12), - "P3/2f=2" => (n=2, l=1, j=3 // 2, f=2, E=957.3965638867406e12), - "P3/2f=3" => (n=2, l=1, j=3 // 2, f=3, E=957.3965608267406e12), - ), - full_transitions=Dict( - ("S1/2f=1", "P1/2f=1") => (multipole="E1", einsteinA=19.4e6), - ("S1/2f=1", "P1/2f=2") => (multipole="E1", einsteinA=19.4e6), - ("S1/2f=2", "P1/2f=1") => (multipole="E1", einsteinA=19.4e6), - ("S1/2f=2", "P1/2f=2") => (multipole="E1", einsteinA=19.4e6), - ("S1/2f=1", "P3/2f=0") => (multipole="E1", einsteinA=19.4e6), - ("S1/2f=1", "P3/2f=1") => (multipole="E1", einsteinA=19.4e6), - ("S1/2f=1", "P3/2f=2") => (multipole="E1", einsteinA=19.4e6), - ("S1/2f=2", "P1/2f=1") => (multipole="E1", einsteinA=19.4e6), - ("S1/2f=2", "P1/2f=2") => (multipole="E1", einsteinA=19.4e6), - ("S1/2f=2", "P1/2f=3") => (multipole="E1", einsteinA=19.4e6), - ), - - # Optional fields - default_sublevel_selection=[ - ("S1/2f=1", "all"), - ("S1/2f=2", "all"), - ("P1/2f=1", "all"), - ("P1/2f=2", "all"), - ("P3/2f=0", "all"), - ("P3/2f=1", "all"), - ("P3/2f=2", "all"), - ("P3/2f=3", "all"), - ], -) - -# boilerplate code -IonInstance{:Be9}(selected_sublevels=nothing, manualshift=Dict()) = - IonInstance{:Be9}(properties_be9, selected_sublevels, manualshift) - -Be9 = IonInstance{:Be9} diff --git a/src/species/ca40.jl b/src/species/ca40.jl deleted file mode 100644 index c9a2c68..0000000 --- a/src/species/ca40.jl +++ /dev/null @@ -1,37 +0,0 @@ -export Ca40 - -const properties_ca40 = IonProperties(; - shortname="⁴⁰Ca", - mass=6.635943757345042e-26, - charge=1, - nuclearspin=0, - full_level_structure=OrderedDict( - "S1/2" => (n=4, l=0, j=1 // 2, f=1 // 2, E=0), - "D3/2" => (n=3, l=2, j=3 // 2, f=3 // 2, E=4.09335071228e14), - "D5/2" => (n=3, l=2, j=5 // 2, f=5 // 2, E=4.1115503183857306e14), - "P1/2" => (n=4, l=1, j=1 // 2, f=1 // 2, E=7.554e14), - "P3/2" => (n=4, l=1, j=3 // 2, f=3 // 2, E=7.621e14), - ), - full_transitions=Dict( - ("S1/2", "D5/2") => (multipole="E2", einsteinA=8.562e-1), - ("S1/2", "P1/2") => (multipole="E1", einsteinA=1.299e8), - ("D3/2", "P1/2") => (multipole="E1", einsteinA=1.060e7), - ("S1/2", "D3/2") => (multipole="E2", einsteinA=9.259e-1), - ("S1/2", "P3/2") => (multipole="E1", einsteinA=1.351e8), - ("D3/2", "P3/2") => (multipole="E1", einsteinA=1.110e6), - ("D5/2", "P3/2") => (multipole="E1", einsteinA=9.901e6), - ), - - # Optional fields - default_sublevel_selection=[("S1/2", "all"), ("D5/2", "all"),], - gfactors=Dict("S1/2" => 2.00225664, "D5/2" => 1.2003340) - - #nonlinear_zeeman = Dict(("S1/2", -1//2) => B->1.3e-4*B^2, - # ("D5/2", -5//2) => B->4.5e-4*B^2) # Syntax example, not numerically accurate -) - -# boilerplate code -IonInstance{:Ca40}(selected_sublevels=nothing, manualshift=Dict()) = - IonInstance{:Ca40}(properties_ca40, selected_sublevels, manualshift) - -Ca40 = IonInstance{:Ca40} diff --git a/src/species/mg25.jl b/src/species/mg25.jl deleted file mode 100644 index b27ce16..0000000 --- a/src/species/mg25.jl +++ /dev/null @@ -1,47 +0,0 @@ -export Mg25 - -const properties_mg25 = IonProperties( - shortname="²⁵Mg", - mass=4.1489954812e-26, - charge=1, - nuclearspin=5 // 2, - full_level_structure=OrderedDict( - "S1/2f=2" => (n=3, l=0, j=1 // 2, f=2, E=1.043445158e9), - "S1/2f=3" => (n=3, l=0, j=1 // 2, f=3, E=-0.745317970e9), - "P1/2f=2" => (n=3, l=1, j=1 // 2, f=2, E=1069.3408690519877e12), - "P1/2f=3" => (n=3, l=1, j=1 // 2, f=3, E=1069.3405639519879e12), - "P3/2f=1" => (n=3, l=1, j=3 // 2, f=1, E=1072.0853411194748e12), - "P3/2f=2" => (n=3, l=1, j=3 // 2, f=2, E=1072.0853033394746e12), - "P3/2f=3" => (n=3, l=1, j=3 // 2, f=3, E=1072.0852466694748e12), - "P3/2f=4" => (n=3, l=1, j=3 // 2, f=4, E=1072.0851711094747e12), - ), - full_transitions=Dict( - ("S1/2f=2", "P1/2f=2") => (multipole="E1", einsteinA=41.3e6), - ("S1/2f=2", "P1/2f=3") => (multipole="E1", einsteinA=41.3e6), - ("S1/2f=3", "P1/2f=2") => (multipole="E1", einsteinA=41.3e6), - ("S1/2f=3", "P1/2f=3") => (multipole="E1", einsteinA=41.3e6), - ("S1/2f=2", "P3/2f=1") => (multipole="E1", einsteinA=41.3e6), - ("S1/2f=2", "P3/2f=2") => (multipole="E1", einsteinA=41.3e6), - ("S1/2f=2", "P3/2f=3") => (multipole="E1", einsteinA=41.3e6), - ("S1/2f=3", "P1/2f=2") => (multipole="E1", einsteinA=41.3e6), - ("S1/2f=3", "P1/2f=3") => (multipole="E1", einsteinA=41.3e6), - ("S1/2f=3", "P1/2f=4") => (multipole="E1", einsteinA=41.3e6), - ), - # Optional fields - default_sublevel_selection=[ - ("S1/2f=2", "all"), - ("S1/2f=3", "all"), - ("P1/2f=2", "all"), - ("P1/2f=3", "all"), - ("P3/2f=1", "all"), - ("P3/2f=2", "all"), - ("P3/2f=3", "all"), - ("P3/2f=4", "all"), - ], -) - -# boilerplate code -IonInstance{:Mg25}(selected_sublevels=nothing, manualshift=Dict()) = - IonInstance{:Mg25}(properties_mg25, selected_sublevels, manualshift) - -Mg25 = IonInstance{:Mg25} diff --git a/src/species/yb171.jl b/src/species/yb171.jl deleted file mode 100644 index 273ced9..0000000 --- a/src/species/yb171.jl +++ /dev/null @@ -1,91 +0,0 @@ -export Yb171 - -const properties_yb171 = IonProperties( - shortname="¹⁷¹Yb", - mass=28.8384644689030595108e-26, - charge=1, - nuclearspin=1 // 2, - full_level_structure=OrderedDict( - # should this be negative, or should I increase the energy of every state by 9GHz? - "S1/2f=0" => (n=6, l=0, j=1 // 2, f=0, E=-9.48210909315e9), - "S1/2f=1" => (n=6, l=0, j=1 // 2, f=1, E=3.16070303105e9), - "D3/2f=0" => (n=5, l=2, j=3 // 2, f=1, E=6.88342460964640e14), - "D3/2f=1" => (n=5, l=2, j=3 // 2, f=2, E=6.88350470564640e14), - "D5/2f=2" => (n=5, l=2, j=5 // 2, f=2, E=7.29477895985202e14), - "D5/2f=3" => (n=5, l=2, j=5 // 2, f=3, E=7.29474121985202e14), - "P1/2f=0" => (n=6, l=1, j=1 // 2, f=0, E=8.112913749003559e14), - "P1/2f=1" => (n=6, l=1, j=1 // 2, f=1, E=8.112934798003559e14), - "[3/2]1/2f=0" => (n=6, l=nothing, j=1 // 2, f=0, E=1.008917341058788e15), - "[3/2]1/2f=1" => (n=6, l=nothing, j=1 // 2, f=1, E=1.008917341058788e15), - "[3/2]3/2f=1" => (n=6, l=nothing, j=3 // 2, f=1, E=8.621425511314839e14), - "[3/2]3/2f=2" => (n=6, l=nothing, j=3 // 2, f=2, E=8.621425511314839e14), - "[5/2]5/2f=2" => (n=6, l=nothing, j=5 // 2, f=2, E=9.70461163716380e14), - "[5/2]5/2f=3" => (n=6, l=nothing, j=5 // 2, f=3, E=9.70461163716380e14), - "F7/2f=3" => (n=6, l=3, j=7 // 2, f=3, E=6.42115934728750e14), - "F7/2f=4" => (n=6, l=3, j=7 // 2, f=4, E=6.42119554728750e14), - ), - full_transitions=Dict( - # Laser lines - # 411nm - ("S1/2f=0", "D5/2f=2") => (multipole="E2", einsteinA=22), - ("S1/2f=1", "D5/2f=2") => (multipole="E2", einsteinA=22), - ("S1/2f=1", "D5/2f=3") => (multipole="E2", einsteinA=22), - # 369nm - ("S1/2f=0", "P1/2f=1") => (multipole="E1", einsteinA=1.155e8), - ("S1/2f=1", "P1/2f=0") => (multipole="E1", einsteinA=1.155e8), - ("S1/2f=1", "P1/2f=1") => (multipole="E1", einsteinA=1.155e8), - # 935nm - ("D3/2f=1", "[3/2]1/2f=0") => (multipole="E1", einsteinA=1.2e5), - ("D3/2f=1", "[3/2]1/2f=1") => (multipole="E1", einsteinA=1.2e5), - ("D3/2f=2", "[3/2]1/2f=1") => (multipole="E1", einsteinA=1.2e5), - # 760nm - ("F7/2f=3", "[3/2]3/2f=1") => (multipole="E2", einsteinA=5e4), - ("F7/2f=3", "[3/2]3/2f=2") => (multipole="E2", einsteinA=5e4), - ("F7/2f=4", "[3/2]3/2f=2") => (multipole="E2", einsteinA=5e4), - # 638nm - # 976nm - # 861nm - # Decay lines - # P1/2 -> D3/2 - ("P1/2f=0", "D3/2f=1") => (multipole="E1", einsteinA=5.7e5), - ("P1/2f=1", "D3/2f=2") => (multipole="E1", einsteinA=5.7e5), - ("P1/2f=1", "D3/2f=2") => (multipole="E1", einsteinA=5.7e5), - # [3/2]1/2 -> S1/2 - ("[3/2]1/2f=0", "S1/2f=1") => (multipole="E1", einsteinA=8.05e7), - ("[3/2]1/2f=1", "S1/2f=0") => (multipole="E1", einsteinA=8.05e7), - ("[3/2]1/2f=1", "S1/2f=1") => (multipole="E1", einsteinA=8.05e7), - # [3/2]3/2 -> S1/2 - ("[3/2]1/2f=0", "S1/2f=1") => (multipole="E1", einsteinA=5.125e7), - ("[3/2]1/2f=1", "S1/2f=0") => (multipole="E1", einsteinA=5.125e7), - ("[3/2]1/2f=1", "S1/2f=1") => (multipole="E1", einsteinA=5.125e7), - ), - # Optional fields - default_sublevel_selection=[ - ("S1/2f=0", "all"), - ("S1/2f=1", "all"), - ("P1/2f=0", "all"), - ("P1/2f=1", "all") - ], - gfactors=Dict( - "S1/2f=0" => 1.998, - "S1/2f=1" => 1.998, - "D5/2f=2" => 1.202, - "D5/2f=3" => 1.202, - "F7/2f=3" => 1.145, - "F7/2f=4" => 1.145, - "[3/2]1/2f=0" => 1.32, - "[3/2]1/2f=1" => 1.32, - "[3/2]3/2f=1" => 1.44, - "[3/2]3/2f=2" => 1.44, - ), - nonlinear_zeeman=Dict( - ("S1/2f=0", 0) => B -> -155.305 * B^2, - ("S1/2f=1", 0) => B -> 155.305 * B^2 - ) -) - -# boilerplate code -IonInstance{:Yb171}(selected_sublevels=nothing, manualshift=Dict()) = - IonInstance{:Yb171}(properties_yb171, selected_sublevels, manualshift) - -Yb171 = IonInstance{:Yb171} diff --git a/test/runtests.jl b/test/runtests.jl index f6710bb..a83c5b0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,8 @@ +# Might be a bad idea, but it's convenient! If there's a more julia-esque idiom +# for test fixtures, that's much preferred. +using IonSim.Properties: loadfromconfig +CA40_PROPERTIES = loadfromconfig("../configs/ions/ca40.yaml") + include("test_analytic_functions.jl") include("test_constants.jl") include("test_operators.jl") @@ -6,5 +11,6 @@ include("test_vibrational_modes.jl") include("test_ions.jl") include("test_ion_traps.jl") include("test_chambers.jl") -# include("test_hamiltonians.jl") +include("test_hamiltonians.jl") include("test_dynamics.jl") +include("test_ionproperties.jl") diff --git a/test/test_chambers.jl b/test/test_chambers.jl index 1bfb5e3..620743c 100644 --- a/test/test_chambers.jl +++ b/test/test_chambers.jl @@ -2,11 +2,12 @@ using QuantumOptics: NLevelBasis, CompositeBasis, FockBasis using Test, IonSim using IonSim.PhysicalConstants: ħ, c using Suppressor +using IonSim.Properties: loadfromconfig @suppress_err begin # set up system - C = Ca40() + C = Ion(loadfromconfig("../configs/ions/ca40.yaml"), nothing) λ = transitionwavelength(C, ("S1/2", "D5/2")) L1 = Laser(λ=λ) L2 = Laser(λ=λ) diff --git a/test/test_dynamics.jl b/test/test_dynamics.jl index 023d459..8e1bbfc 100644 --- a/test/test_dynamics.jl +++ b/test/test_dynamics.jl @@ -1,12 +1,13 @@ using QuantumOptics using Test, IonSim using IonSim.analytical +using IonSim.Properties: loadfromconfig using Suppressor @suppress_err begin @testset "Dynamics -- Rabi" begin # carrier - C = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) + C = Ion(CA40_PROPERTIES, [("S1/2", -1 / 2), ("D5/2", -1 / 2)]) L = Laser() chain = LinearChain( ions=[C], @@ -127,7 +128,7 @@ using Suppressor # could use some better tests # use ϕ(t) to set laser frequency to resonance - C = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) + C = Ion(CA40_PROPERTIES, [("S1/2", -1 / 2), ("D5/2", -1 / 2)]) L = Laser() chain = LinearChain( ions=[C], @@ -199,7 +200,7 @@ using Suppressor @testset "molmer-sorensen" begin # slow gate - C = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) + C = Ion(CA40_PROPERTIES, [("S1/2", -1 / 2), ("D5/2", -1 / 2)]) L1 = Laser() L2 = Laser() chain = LinearChain( diff --git a/test/test_hamiltonians.jl b/test/test_hamiltonians.jl index 946ef1b..6e65d35 100644 --- a/test/test_hamiltonians.jl +++ b/test/test_hamiltonians.jl @@ -109,7 +109,10 @@ end @testset "hamiltonians -- parameter arrays" begin # setup system - C = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) + C = Ion( + loadfromconfig("../configs/ions/ca40.yaml"), + [("S1/2", -1 / 2), ("D5/2", -1 / 2)] + ) C1 = copy(C) L1 = Laser() @@ -220,7 +223,7 @@ end # tests for _setup_global_B_hamiltonian # setup system - C = Ca40() + C = Ion(CA40_PROPERTIES, nothing) L = Laser() L.λ = transitionwavelength(C, ("S1/2", "D5/2")) chain = LinearChain( @@ -321,7 +324,7 @@ end # _setup_base_hamiltonian # Let's make sure that _setup_base_hamiltonian is recording the appropriate indices. - C = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) + C = Ion(CA40_PROPERTIES, [("S1/2", -1 / 2), ("D5/2", -1 / 2)]) L = Laser() L.λ = transitionwavelength(C, ("S1/2", "D5/2")) chain = LinearChain( @@ -357,7 +360,7 @@ end @test length(unique([cidxs; c])) - 1 == (N^2 - N) == length(cidxs) ## test for two modes - C1 = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) + C1 = Ion(CA40_PROPERTIES, [("S1/2", -1 / 2), ("D5/2", -1 / 2)]) L1 = Laser() L1.λ = transitionwavelength(C1, ("S1/2", "D5/2")) chain1 = LinearChain( @@ -408,8 +411,8 @@ end # * Hamiltonian built with hamiltonians.jl # define ion, laser, chain, trap - C_a = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) - C_b = Ca40([("S1/2", -1 / 2), ("D5/2", -1 / 2)]) + C_a = Ion(CA40_PROPERTIES, [("S1/2", -1 / 2), ("D5/2", -1 / 2)]) + C_b = Ion(CA40_PROPERTIES, [("S1/2", -1 / 2), ("D5/2", -1 / 2)]) L = Laser() L.pointing = [(1, 1.0), (2, 1.0)] chain = LinearChain( @@ -466,8 +469,8 @@ end H = hamiltonian(T, timescale=1e-6, lamb_dicke_order=101) H1 = hamiltonian(T, timescale=1e-6, lamb_dicke_order=101, time_dependent_eta=true) # test similarity at a random time input - @test norm(qoH(tp).data - H(tp, 0).data) < 1e-2 - @test norm(qoH(tp).data - H1(tp, 0).data) < 1e-2 + @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4 + @test norm(qoH(tp).data - H1(tp, 0).data) < 1e-4 # Case 1b: analytic solution (as opposed to truncated solution) mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="analytic") @@ -486,7 +489,7 @@ end displacement="analytic", time_dependent_eta=true ) - @test norm(qoH(tp).data - H(tp, 0).data) < 1e-2 + @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4 @test H1(tp, 0).data ≈ H(tp, 0).data # Case 2a: full hamiltonian (w/o conj_repeated_indices) @@ -500,8 +503,8 @@ end ) mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="truncated") mode_op2(t; η) = displace(mode2, im * η * exp(im * 2π * √3 * t), method="truncated") - @test norm(qoH(tp).data - H(tp, 0).data) < 1e-2 - @test norm(qoH(tp).data - H1(tp, 0).data) < 1e-2 + @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4 + @test norm(qoH(tp).data - H1(tp, 0).data) < 1e-4 # Case 2b: Like 2a, but with analytic solution mode_op1(t; η) = displace(mode1, im * η * exp(im * 2π * t), method="analytic") @@ -522,7 +525,7 @@ end time_dependent_eta=true, rwa_cutoff=1e10 ) - @test norm(qoH(tp).data - H(tp, 0).data) < 1e-2 + @test norm(qoH(tp).data - H(tp, 0).data) < 1e-4 @test H1(tp, 0).data ≈ H(tp, 0).data # Case 3: test Hamiltonian with zero Lamb-Dicke value, i.e. no vibrational modes diff --git a/test/test_ion_traps.jl b/test/test_ion_traps.jl index ab006d9..b868e3c 100644 --- a/test/test_ion_traps.jl +++ b/test/test_ion_traps.jl @@ -1,17 +1,13 @@ using Test, IonSim +# using IonSim.PhysicalConstants using Suppressor -using IonSim: - linear_equilibrium_positions, - searchfor_trappingpotential_parameters, - diagonalize_Kij, - diagonalize_Kij_for_homogeneous_chain @suppress_err begin @testset "ion_traps -- LinearChain" begin - C = Ca40() + C = Ion(CA40_PROPERTIES, nothing) lc = LinearChain( ions=[C, C, C, C], - comfrequencies=(x=5e6, y=5e6, z=1e6), + comfrequencies=(x=5, y=5, z=1), selectedmodes=(y=[1], z=[4]) ) @test ions(lc) == lc.ions @@ -24,13 +20,11 @@ using IonSim: @test [ionnumber(I) for I in lc.ions] == [1, 2, 3, 4] # cmode -> pre-evaluated equilibrium positions for four ion chain - cmode = [ - -2.1766240336602492e-5, - -6.883427582857696e-6, - 6.883427582857696e-6, - 2.1766240336602492e-5 - ] - @test ionpositions(lc) ≈ cmode rtol = 1e-6 + cmode = [-0.06392393573, -0.0202155287427, 0.0202155287427, 0.0639239357] + if ionposition(lc.ions[1]) > 0 + cmode .*= -1 + end + @test [ionposition(I) for I in lc.ions] ≈ cmode rtol = 1e-6 # should get warning if same ion is input multiple times to ion kwarg warning = "Some ions point to the same thing. Making copies." @@ -50,8 +44,6 @@ using IonSim: # Make sure there are no errors with print/show @suppress print(lc) show(lc) - - @test length(lc) == length(lc.ions) end @testset "ion_traps -- general" begin @@ -61,191 +53,23 @@ using IonSim: pos = [posL; -1 .* reverse(posL)] @test any(isapprox.(linear_equilibrium_positions(10), pos, rtol=1e-4)) + # and test calculation of characterstic length scale for linear chain, equal mass + C = Ion(CA40_PROPERTIES, nothing) + @test characteristic_length_scale(mass(C), 1e6) ≈ 4.449042804354206e-6 + + # and do the same for Anm, which computes the normal modes + @test_throws AssertionError Anm(2, (x=0.5, y=0.5, z=1), (x=1, y=0, z=0)) + cst = [-0.2132, 0.6742, -0.6742, 0.2132] + freq, mode = Anm(4, (x=2, y=2, z=1), (x=0, y=0, z=1))[end] + @test freq ≈ √9.308 rtol = 1e-4 + if mode[1] > 0 + cst = -1 .* cst + end + @test any(isapprox.(mode, cst, rtol=1e-4)) + # test '_sparsify', which should drop small values from an array x = [1e-6, 1e-5] IonSim._sparsify!(x, 2e-6) @test any(x .== [0, 1e-5]) - - # ions field must be an iterable - C = Ca40() - @test_throws AssertionError LinearChain( - ions=C, - comfrequencies=(x=4, y=4, z=1), - selectedmodes=(x=[], y=[], z=[1, 2]) - ) - @test_throws AssertionError full_normal_mode_description(C, (x=1, y=1, z=1)) - end - - @testset "ion_configurations -- mixed species" begin - #################################################################################### - # mixed-species linear chain tests - #################################################################################### - c = Ca40() - b = Be9() - y = Yb171() - - # homogeneous chain - - ion_list = [c for _ in 1:rand(2:7)] - randval = rand() * 1e6 - target_eigenfrequencies = (x=randval, y=randval, z=randval / 10) - - # compare optimization routine to james formula for homogeneous case https://doi.org/10.1007/s003400050373 - res = searchfor_trappingpotential_parameters(ion_list, target_eigenfrequencies) - nmstructure_x_direction = diagonalize_Kij(ion_list, x̂, res...) - nmstructure_y_direction = diagonalize_Kij(ion_list, ŷ, res...) - nmstructure_z_direction = diagonalize_Kij(ion_list, ẑ, res...) - - for (i, e) in enumerate(nmstructure_x_direction) - @test e[1] ≈ diagonalize_Kij_for_homogeneous_chain( - length(ion_list), - target_eigenfrequencies, - x̂ - )[i][1] - @test e[2] ≈ diagonalize_Kij_for_homogeneous_chain( - length(ion_list), - target_eigenfrequencies, - x̂ - )[i][2] - end - - for (i, e) in enumerate(nmstructure_y_direction) - @test e[1] ≈ diagonalize_Kij_for_homogeneous_chain( - length(ion_list), - target_eigenfrequencies, - ŷ - )[i][1] - @test e[2] ≈ diagonalize_Kij_for_homogeneous_chain( - length(ion_list), - target_eigenfrequencies, - ŷ - )[i][2] - end - - for (i, e) in enumerate(nmstructure_z_direction) - @test e[1] ≈ diagonalize_Kij_for_homogeneous_chain( - length(ion_list), - target_eigenfrequencies, - ẑ - )[i][1] - @test e[2] ≈ diagonalize_Kij_for_homogeneous_chain( - length(ion_list), - target_eigenfrequencies, - ẑ - )[i][2] - end - - # test for assertion error if x > y > z not obeyed - @test_throws AssertionError searchfor_trappingpotential_parameters( - ion_list, - ( - x=target_eigenfrequencies.z, - y=target_eigenfrequencies.x, - z=target_eigenfrequencies.y - ) - ) - - # heterogeneous chain - - #= - Compare to exact formulas for axial eigenvalues and eigenvectors of a two-ion mixed-species crystal: - Eur. Phys. J. D 13, 261–269 (2001) https://doi.org/10.1007/s100530170275 - =# - - function two_ion_eigenvalues(k, masses) - m = minimum(masses) - M = maximum(masses) - μ = M / m - C1 = 1 + (1 / μ) - C2 = sqrt(1 + 1 / μ^2 - 1 / μ) - return sqrt.((k / m) * [C1 - C2, C1 + C2]) - end - - function two_ion_eigenvectors(masses) - m = minimum(masses) - M = maximum(masses) - μ = M / m - C = √(1 + μ^2 - μ) - q1 = normalize([√μ, 1 - μ + C]) - q2 = normalize([√μ, 1 - μ - C]) - if masses[2] > masses[1] - reverse!(q1) - reverse!(q2) - end - q1 *= sign(q1[1]) - q2 *= sign(q2[1]) - return [q1, q2] - end - - for ion_list in [[c, b], [b, c], [c, y], [y, c], [b, y], [y, b]] - target_eigenfrequencies = (x=20e6, y=20e6, z=1e6) # have to make these huge for [y, b] case - res = searchfor_trappingpotential_parameters(ion_list, target_eigenfrequencies) - computed_values = diagonalize_Kij(ion_list, ẑ, res...) - computed_eigenvalues = [i[1] for i in computed_values] - computed_eigenvectors = [i[2] for i in computed_values] - - M = [mass(ion) for ion in ion_list] - analytic_eigenvectors = two_ion_eigenvectors(M) - analytic_eigenvalues = two_ion_eigenvalues(res[1], M) - - # for some reason this isn't working here, even though it works fine in a notebook - # also returns the same values for every itration of the loop - # @test analytic_eigenvalues ≈ computed_eigenvalues rtol=1e-3 - @test analytic_eigenvectors ≈ computed_eigenvectors rtol = 1e-3 - end - - #= - Compare to eigenvalues from PHYSICAL REVIEW A103, 012610 (2021) - =# - - M = [mass(y), mass(y), mass(y), mass(y), 2.28997e-25] - Q = ones(Int, length(M)) - target_eigenfrequencies = (x=3.039e6, y=3.039e6, z=0.158e6) - res = searchfor_trappingpotential_parameters(M, Q, target_eigenfrequencies) - computed_values = diagonalize_Kij(M, Q, ẑ, res...) - axial_eigenfrequencies = [round(i[1] / 1e6, digits=3) for i in computed_values] - - axial_values_from_PRA = [0.158, 0.28, 0.388, 0.48, 0.57] - @test axial_eigenfrequencies ≈ axial_values_from_PRA - - # Note: we don't test radial modes against this PRA since there is an unresolved - # quantitative (but not qualitative difference). Normal modes are the same, - # eigenfrequencies are gapped between modes that mostly vibrate in one or the other - # species and the general scaling seems to agree. but the gap width is larger for - # our calculation and I'm not sure why - end - - @testset "load from yaml" begin - c = Ca40() - chain = LinearChain_fromyaml( - ions=[c, c, c, c], - yaml="test_load_LinearChain_fromyaml.yaml", - ) - @test isnan(chain.comfrequencies.x) - @test chain.ionpositions == [1, 2, 3, 4] - @test chain.ions == ions(chain) - @test_throws AssertionError full_normal_mode_description(chain) - lc = LinearChain_fromyaml( - ions=[c, c, c, c], - yaml="test_load_LinearChain_fromyaml_broken.yaml", - ) - @test length(lc.selectedmodes.x) == 0 - end - - @testset "visualize-iontraps" begin - # for now just run all of the plots and make sure nothing errors - c = Ca40() - b = Be9() - y = Yb171() - chain = LinearChain( - ions=[c, b, b, c, y, y], - comfrequencies=(x=20e6, y=20e6, z=0.1e6), - selectedmodes=(; x=[1, 3], z=[1, 3:4]) - ) - visualize(chain, ẑ, [:], format="circles") - visualize(chain, x̂, [1]) - visualize(chain, ẑ, [1:4], format="circles") - visualize(xmodes(chain)[1], format="circles") - visualize(xmodes(chain)[1], format="bar") end end # end suppress diff --git a/test/test_ionproperties.jl b/test/test_ionproperties.jl new file mode 100644 index 0000000..4d6da5f --- /dev/null +++ b/test/test_ionproperties.jl @@ -0,0 +1,25 @@ +using Test, IonSim +using IonSim.Properties: process_nonlinear_zeeman, loadfromconfig +using Suppressor + +@suppress_err begin + @testset "properties -- Verify nonlinear zeeman shift mapping" begin + nonlinear_zeeman_config = [ + Dict("level" => "S1/2=0", "sublevel" => 0, "coeffs" => [1.0, 2.0, 3.0]), + Dict("level" => "S1/2=1", "sublevel" => 0, "coeffs" => [10.0, 20.0, 30.0]) + ] + + processed_functions = process_nonlinear_zeeman(nonlinear_zeeman_config) + + # Really, want to verify that the functions mapped to different + # level/sublevel pairs are actually distinct from each other. + @test isapprox(processed_functions["S1/2=0", 0](2.0), 17.0, rtol=1e-5) + @test isapprox(processed_functions["S1/2=1", 0](2.0), 170.0, rtol=1e-5) + end + + @testset "properties -- Verify that all provided ion configs parse" begin + for ion_fname in readdir("../configs/ions") + loadfromconfig(joinpath("../configs/ions", ion_fname)) + end + end +end diff --git a/test/test_ions.jl b/test/test_ions.jl index 244b504..c7de367 100644 --- a/test/test_ions.jl +++ b/test/test_ions.jl @@ -1,11 +1,12 @@ using QuantumOptics: NLevelBasis, nlevelstate using Test, IonSim +using IonSim.Properties: IonProperties using IonSim.PhysicalConstants using Suppressor @suppress_err begin @testset "ions -- general" begin - C = Ca40() + C = Ion(CA40_PROPERTIES, nothing) # test for required fields @test typeof(speciesproperties(C)) <: IonProperties @test typeof(sublevels(C)) == Vector{Tuple{String, Real}} @@ -17,7 +18,7 @@ using Suppressor @test typeof(ionposition(C)) == Missing # test == - C1 = Ca40() + C1 = Ion(CA40_PROPERTIES, nothing) @test C1 == C # test aliases @@ -66,16 +67,8 @@ using Suppressor @test transitionwavelength(C, ("S", "D")) ≈ 7.29147e-7 end - @testset "ions -- species" begin - # attempt to instantiate all Ion subtypes (use default sublevel selection) - for s in [Be9, Ca40, Mg25, Yb171] - ion = s() - @test typeof(ion) <: Ion - end - end - @testset "ions -- Ca40" begin - C = Ca40() + C = Ion(CA40_PROPERTIES, nothing) # set some aliases for convenience sublevelalias!(C, ("S1/2", -1 / 2), "S") @@ -109,7 +102,7 @@ using Suppressor @test_throws AssertionError C[""] # # test indexing - C1 = Ca40([("S1/2", -1 / 2), ("D5/2", -5 / 2)]) + C1 = Ion(CA40_PROPERTIES, [("S1/2", -1 / 2), ("D5/2", -5 / 2)]) sublevelalias!(C1, ("S1/2", -1 / 2), "S") sublevelalias!(C1, ("D5/2", -5 / 2), "D") @test C1[("S1/2", -1 / 2)].data == ComplexF64[1; 0] diff --git a/test/test_load_LinearChain_fromyaml.yaml b/test/test_load_LinearChain_fromyaml.yaml deleted file mode 100644 index 2bd55a9..0000000 --- a/test/test_load_LinearChain_fromyaml.yaml +++ /dev/null @@ -1,9 +0,0 @@ -x: - - frequency: 1e6 - mode: [0.1, 0.5, 0.3, 0.8] - - frequency: 2e6 - mode: [0.3, 0.6, 0.5, 3] -y: - - frequency: 8e6 - mode: [1, 1, 1, 1] -ionpositions: [1,2,3,4] \ No newline at end of file diff --git a/test/test_load_LinearChain_fromyaml_broken.yaml b/test/test_load_LinearChain_fromyaml_broken.yaml deleted file mode 100644 index aef45f3..0000000 --- a/test/test_load_LinearChain_fromyaml_broken.yaml +++ /dev/null @@ -1,8 +0,0 @@ -x: - - frequency: 1e6 - mode: [0.1, 0.5, 0.3, 0.8, 1] - - frequency: 2e6 -y: - - frequency: 8e6 - mode: [1, 1, 1, 1] -ionpositions: [1,2,3,4] \ No newline at end of file diff --git a/test/test_operators.jl b/test/test_operators.jl index 69cc1cd..671ecd5 100644 --- a/test/test_operators.jl +++ b/test/test_operators.jl @@ -6,7 +6,8 @@ using Suppressor @suppress_err begin # setup system - C = Ca40([("S1/2", -1 // 2), ("D5/2", -1 // 2)]) + print(Base.Filesystem.pwd()) + C = Ion(CA40_PROPERTIES, [("S1/2", -1 // 2), ("D5/2", -1 // 2)]) chain = LinearChain( ions=[C, C], comfrequencies=(x=2, y=2, z=1),