Skip to content

Trapped Ions

oqd_trical.mechanical.trappedions

TrappedIons

Bases: Base

Object representing a system of trapped ions.

Parameters:

Name Type Description Default
N int

Number of ions.

 required
ps Potential

Potentials on the system.

 ()

Other Parameters:

Name Type Description
dim int

Dimension of system.
l float

Length scale of system.
m float

Mass of ions
q float

Charge of ions
Source code in oqd-trical/src/oqd_trical/mechanical/trappedions.py 
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
class TrappedIons(Base):
    """
    Object representing a system of trapped ions.

    Args:
        N (int): Number of ions.
        ps (Potential): Potentials on the system.

    Keyword Args:
        dim (int): Dimension of system.
        l (float): Length scale of system.
        m (float): Mass of ions
        q (float): Charge of ions

    """

    def __init__(self, N, *ps, **kwargs):
        super(TrappedIons, self).__init__()

        params = {
            "dim": ps[0].dim if "dim" in ps[0].__dict__.keys() else 3,
            "l": 1e-6,
            "m": 171 * cst.m_u,
            "q": cst.e,
        }
        params.update(kwargs)
        self.__dict__.update(params)

        self.N = N
        self.cp = CoulombPotential(N, dim=self.dim, q=self.q)

        for p in ps:
            p.update_params(N=N)
        self.ps = np.array(ps)

        self.fp = self.cp + self.ps.sum()
        pass

    def equilibrium_position(self, opt=dflt_opt, **kwargs):
        """
        Function that calculates the equilibrium position of the ions.

        Args:
            opt (Callable): Generator of the appropriate optimization function for finding the equilibrium position.

        Returns:
            x_ep (np.ndarray[float]): Equilibrium position of the ions.
        """
        ndcp = self.cp.nondimensionalize(self.l)
        ndps = np.array([p.nondimensionalize(self.l) for p in self.ps])
        ndfp = ndcp + ndps.sum()

        def _ndfp(x):
            return ndfp(x.reshape(self.dim, self.N).transpose())

        self.x_ep = (
            opt(self, **kwargs)(_ndfp).reshape(self.dim, self.N).transpose() * self.l
        )
        return self.x_ep

    def normal_modes(self, block_sort=False):
        """
        Function that calculates the normal modes of the system.

        Args:
            block_sort (bool): Indicator to apply block sorting.

        Returns:
            w (np.ndarray[float]): Normal mode frequencies of the system.
            b (np.ndarray[float]): Normal mode eigenvectors of the system.
        """
        ndcp = self.cp.nondimensionalize(self.l)
        ndps = np.array([p.nondimensionalize(self.l) for p in self.ps])

        def hess_phi(x):
            return np.array([ndp.hessian()(x) for ndp in np.append(ndps, ndcp)]).sum(0)

        if "x_ep" not in self.__dict__.keys():
            self.equilibrium_position()

        hess_phi_x_ep = hess_phi(self.x_ep / self.l)

        if isinstance(self.m, float):
            A = hess_phi_x_ep / self.m
            w, b = np.linalg.eigh(A)
        else:
            A = np.einsum(
                "ij,i,j->ij",
                hess_phi_x_ep,
                1 / np.tile(np.sqrt(self.m), 3),
                1 / np.tile(np.sqrt(self.m), 3),
            )
            w, b = np.linalg.eigh(A)
            b = np.einsum("im,i->im", b, 1 / np.tile(np.sqrt(self.m), 3))
            b = b / np.linalg.norm(b, axis=0)

        w = np.sqrt(w * cst.k_e * cst.e**2 / self.l**3)

        _b = np.round(np.copy(b), 3).transpose()
        s = (np.sign(_b).sum(1) >= 0) * 2 - 1
        b = np.einsum("n,in->in", s, b)

        if block_sort:
            n = np.round(
                np.array(
                    [
                        np.linalg.norm(
                            b[i * self.N : (i + 1) * self.N].transpose(), axis=-1
                        )
                        for i in range(self.dim)
                    ]
                )
            ).astype(int)
            idcs = np.lexsort(np.concatenate(((-w).reshape(1, -1), n)))
        else:
            idcs = np.argsort(np.flip(w, axis=0))

        self.w, self.b, self.A = (
            w[idcs],
            b[:, idcs],
            A * cst.k_e * cst.e**2 / self.l**3,
        )
        return self.w, self.b

    def principal_axes(self, tol=1e-3):
        """
        Function that calculates the principle axes of the system.

        Args:
            tol (float): Tolerance for evaluating orthogonality of principal axes.

        Returns:
            x_pa (np.ndarray[float]): Principle axes of the system.
            w_pa (np.ndarray[float]): Normal mode frequencies of the system.
            b_pa (np.ndarray[float]): Normal mode eigenvectors of the system in the principal axes coordinate system.
        """
        if np.isin(np.array(["w", "b"]), np.array(self.__dict__.keys())).sum() != 2:
            self.normal_modes()

        xs = self.b.reshape(self.dim, -1).transpose()
        x_pa = orthonormal_subset(xs, tol=tol)

        assert len(x_pa) == self.dim, "Principle axes do not exist"

        _x_pa = np.round(np.copy(x_pa), 3)
        s = (np.sign(_x_pa).sum(1) >= 0) * 2 - 1
        x_pa = np.einsum("n,ni->ni", s, x_pa)

        _x_pa = np.round(np.copy(x_pa), 3)
        idcs = np.lexsort(
            np.concatenate((_x_pa.transpose(), np.abs(_x_pa).transpose()))
        )

        x_pa = x_pa[idcs]

        b_pa = np.einsum("ij,jn->in", x_pa, self.b.reshape(self.dim, -1)).reshape(
            self.dim * self.N, -1
        )
        w_pa = self.w

        _b = np.round(np.copy(b_pa), 3).transpose()
        s = (np.sign(_b).sum(1) >= 0) * 2 - 1
        b_pa = np.einsum("n,in->in", s, b_pa)

        n = np.round(
            np.array(
                [
                    np.linalg.norm(
                        b_pa[i * self.N : (i + 1) * self.N].transpose(), axis=-1
                    )
                    for i in range(self.dim)
                ]
            )
        ).astype(int)
        idcs = np.lexsort(np.concatenate(((-w_pa).reshape(1, -1), n)))

        self.x_pa, self.w_pa, self.b_pa = x_pa, w_pa[idcs], b_pa[:, idcs]
        return self.x_pa, self.w_pa, self.b_pa

    def update_params(self, **kwargs):
        """
        Updates parameters, i.e. params attribute, of a TrappedIons object.

        Args:
            dim (float): Dimension of the system.
            l (float): Length scale of the system.
            m (float): Mass of an ion.
            q (float: Charge of an ion.
        """
        self.params.update(kwargs)
        self.__dict__.update(self.params)
        pass

    def mode_ion_coupling(self):
        if (
            np.isin(
                np.array(["w_pa", "b_pa", "x_pa"]), np.array(self.__dict__.keys())
            ).sum()
            != 3
        ):
            self.principal_axes()

        mic = np.zeros((3 * self.N, 3 * self.N, 3 * self.N))
        idcs = np.triu_indices(3 * self.N, k=1)
        mic[:, idcs[0], idcs[1]] = mic[:, idcs[1], idcs[0]] = (
            self.b_pa[np.array(idcs), :].prod(axis=0).transpose()
        )
        self.mic = mic
        return self.mic

    pass

equilibrium_position(opt=dflt_opt, **kwargs)

Function that calculates the equilibrium position of the ions.

Parameters:

Name Type Description Default
opt Callable

Generator of the appropriate optimization function for finding the equilibrium position.

 dflt_opt

Returns:

Name Type Description
x_ep ndarray[float]

Equilibrium position of the ions.
Source code in oqd-trical/src/oqd_trical/mechanical/trappedions.py 
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
def equilibrium_position(self, opt=dflt_opt, **kwargs):
    """
    Function that calculates the equilibrium position of the ions.

    Args:
        opt (Callable): Generator of the appropriate optimization function for finding the equilibrium position.

    Returns:
        x_ep (np.ndarray[float]): Equilibrium position of the ions.
    """
    ndcp = self.cp.nondimensionalize(self.l)
    ndps = np.array([p.nondimensionalize(self.l) for p in self.ps])
    ndfp = ndcp + ndps.sum()

    def _ndfp(x):
        return ndfp(x.reshape(self.dim, self.N).transpose())

    self.x_ep = (
        opt(self, **kwargs)(_ndfp).reshape(self.dim, self.N).transpose() * self.l
    )
    return self.x_ep

normal_modes(block_sort=False)

Function that calculates the normal modes of the system.

Parameters:

Name Type Description Default
block_sort bool

Indicator to apply block sorting.

 False

Returns:

Name Type Description
w ndarray[float]

Normal mode frequencies of the system.
b ndarray[float]

Normal mode eigenvectors of the system.
Source code in oqd-trical/src/oqd_trical/mechanical/trappedions.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
def normal_modes(self, block_sort=False):
    """
    Function that calculates the normal modes of the system.

    Args:
        block_sort (bool): Indicator to apply block sorting.

    Returns:
        w (np.ndarray[float]): Normal mode frequencies of the system.
        b (np.ndarray[float]): Normal mode eigenvectors of the system.
    """
    ndcp = self.cp.nondimensionalize(self.l)
    ndps = np.array([p.nondimensionalize(self.l) for p in self.ps])

    def hess_phi(x):
        return np.array([ndp.hessian()(x) for ndp in np.append(ndps, ndcp)]).sum(0)

    if "x_ep" not in self.__dict__.keys():
        self.equilibrium_position()

    hess_phi_x_ep = hess_phi(self.x_ep / self.l)

    if isinstance(self.m, float):
        A = hess_phi_x_ep / self.m
        w, b = np.linalg.eigh(A)
    else:
        A = np.einsum(
            "ij,i,j->ij",
            hess_phi_x_ep,
            1 / np.tile(np.sqrt(self.m), 3),
            1 / np.tile(np.sqrt(self.m), 3),
        )
        w, b = np.linalg.eigh(A)
        b = np.einsum("im,i->im", b, 1 / np.tile(np.sqrt(self.m), 3))
        b = b / np.linalg.norm(b, axis=0)

    w = np.sqrt(w * cst.k_e * cst.e**2 / self.l**3)

    _b = np.round(np.copy(b), 3).transpose()
    s = (np.sign(_b).sum(1) >= 0) * 2 - 1
    b = np.einsum("n,in->in", s, b)

    if block_sort:
        n = np.round(
            np.array(
                [
                    np.linalg.norm(
                        b[i * self.N : (i + 1) * self.N].transpose(), axis=-1
                    )
                    for i in range(self.dim)
                ]
            )
        ).astype(int)
        idcs = np.lexsort(np.concatenate(((-w).reshape(1, -1), n)))
    else:
        idcs = np.argsort(np.flip(w, axis=0))

    self.w, self.b, self.A = (
        w[idcs],
        b[:, idcs],
        A * cst.k_e * cst.e**2 / self.l**3,
    )
    return self.w, self.b

principal_axes(tol=0.001)

Function that calculates the principle axes of the system.

Parameters:

Name Type Description Default
tol float

Tolerance for evaluating orthogonality of principal axes.

 0.001

Returns:

Name Type Description
x_pa ndarray[float]

Principle axes of the system.
w_pa ndarray[float]

Normal mode frequencies of the system.
b_pa ndarray[float]

Normal mode eigenvectors of the system in the principal axes coordinate system.
Source code in oqd-trical/src/oqd_trical/mechanical/trappedions.py 
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
def principal_axes(self, tol=1e-3):
    """
    Function that calculates the principle axes of the system.

    Args:
        tol (float): Tolerance for evaluating orthogonality of principal axes.

    Returns:
        x_pa (np.ndarray[float]): Principle axes of the system.
        w_pa (np.ndarray[float]): Normal mode frequencies of the system.
        b_pa (np.ndarray[float]): Normal mode eigenvectors of the system in the principal axes coordinate system.
    """
    if np.isin(np.array(["w", "b"]), np.array(self.__dict__.keys())).sum() != 2:
        self.normal_modes()

    xs = self.b.reshape(self.dim, -1).transpose()
    x_pa = orthonormal_subset(xs, tol=tol)

    assert len(x_pa) == self.dim, "Principle axes do not exist"

    _x_pa = np.round(np.copy(x_pa), 3)
    s = (np.sign(_x_pa).sum(1) >= 0) * 2 - 1
    x_pa = np.einsum("n,ni->ni", s, x_pa)

    _x_pa = np.round(np.copy(x_pa), 3)
    idcs = np.lexsort(
        np.concatenate((_x_pa.transpose(), np.abs(_x_pa).transpose()))
    )

    x_pa = x_pa[idcs]

    b_pa = np.einsum("ij,jn->in", x_pa, self.b.reshape(self.dim, -1)).reshape(
        self.dim * self.N, -1
    )
    w_pa = self.w

    _b = np.round(np.copy(b_pa), 3).transpose()
    s = (np.sign(_b).sum(1) >= 0) * 2 - 1
    b_pa = np.einsum("n,in->in", s, b_pa)

    n = np.round(
        np.array(
            [
                np.linalg.norm(
                    b_pa[i * self.N : (i + 1) * self.N].transpose(), axis=-1
                )
                for i in range(self.dim)
            ]
        )
    ).astype(int)
    idcs = np.lexsort(np.concatenate(((-w_pa).reshape(1, -1), n)))

    self.x_pa, self.w_pa, self.b_pa = x_pa, w_pa[idcs], b_pa[:, idcs]
    return self.x_pa, self.w_pa, self.b_pa

update_params(**kwargs)

Updates parameters, i.e. params attribute, of a TrappedIons object.

Parameters:

Name Type Description Default
dim float

Dimension of the system.

 required
l float

Length scale of the system.

 required
m float

Mass of an ion.

 required
q float

Charge of an ion.

 required
Source code in oqd-trical/src/oqd_trical/mechanical/trappedions.py 
207
208
209
210
211
212
213
214
215
216
217
218
219
def update_params(self, **kwargs):
    """
    Updates parameters, i.e. params attribute, of a TrappedIons object.

    Args:
        dim (float): Dimension of the system.
        l (float): Length scale of the system.
        m (float): Mass of an ion.
        q (float: Charge of an ion.
    """
    self.params.update(kwargs)
    self.__dict__.update(self.params)
    pass