From f76bfc18e025df9ead85b4adecfbf9300a90aaf3 Mon Sep 17 00:00:00 2001 From: sambonkov Date: Thu, 8 Aug 2024 01:29:51 -0600 Subject: [PATCH] Clean up calculations for description. --- SQcircuit/circuit.py | 178 +++++++++++------- SQcircuit/symbolic.py | 75 ++++---- .../descriptions/zero_pi_inductors_ltx.txt | 2 +- .../descriptions/zero_pi_inductors_txt.txt | 2 +- .../descriptions/zero_pi_junctions_ltx.txt | 2 +- .../descriptions/zero_pi_junctions_txt.txt | 4 +- .../tests/{test_paper.py => test_papers.py} | 2 +- SQcircuit/texts.py | 92 ++++----- SQcircuit/torch_extensions.py | 2 +- 9 files changed, 208 insertions(+), 151 deletions(-) rename SQcircuit/tests/{test_paper.py => test_papers.py} (99%) diff --git a/SQcircuit/circuit.py b/SQcircuit/circuit.py index 8e3f9bc..1ccbab4 100755 --- a/SQcircuit/circuit.py +++ b/SQcircuit/circuit.py @@ -1148,36 +1148,19 @@ def _transform_hamil(self) -> None: # scaling the modes by third transformation self.S3, self.R3 = self._get_and_apply_transformation_3() - def _compute_params(self) -> None: + def _compute_descrip_vars(self) -> None: """Compute parameters of transformed Hamiltonian to reduced number of significant figures and in the frequency units of ``SQcircuit``; store in the ``.descrip_vars`` dictionary for easy access. """ - # dimensions of modes of circuit + # Dimensions of modes of circuit self.descrip_vars['n_modes'] = self.n self.descrip_vars['har_dim'] = np.sum(self.omega != 0) - har_dim = self.descrip_vars['har_dim'] self.descrip_vars['charge_dim'] = np.sum(self.omega == 0) - self.descrip_vars['omega'] = ( - self.omega / (2 * np.pi * unt.get_unit_freq()) - ) - - self.descrip_vars['phi_zp'] = ( - (2 * np.pi / unt.Phi0) - * np.sqrt(unt.hbar - / (2 * np.sqrt( - np.diag(self.lTrans)[:har_dim] - / np.diag(self.cInvTrans)[:har_dim] - )))) - self.descrip_vars['ng'] = [self.charge_islands[i].value() - for i in range(har_dim, self.n)] - self.descrip_vars['EC'] = ( - ((2 * unt.e) ** 2 / (unt.hbar * 2 * np.pi * unt.get_unit_freq())) - * np.diag(np.repeat(0.5, self.n)) - * self.cInvTrans - ) + self.descrip_vars['n_loops'] = len(self.loops) + # Rounded matrices describing circuit layout self.descrip_vars['W'] = np.round(self.wTrans, 6) self.descrip_vars['S'] = np.round(self.S, 3) if self.loops: @@ -1186,32 +1169,52 @@ def _compute_params(self) -> None: self.descrip_vars['B'] = np.zeros((len(self.elem_keys[Junction]) + len(self.elem_keys[Inductor]), 1)) + # Harmonic mode variables + self.descrip_vars['omega'] = ( + self.omega / (2 * np.pi * unt.get_unit_freq()) + ) + self.descrip_vars['phi_zp'] = [ + self._phi_zp(i) for i in range(self.descrip_vars['har_dim']) + ] + + # Compute the element values in the correct units + ## Junction energies self.descrip_vars['EJ'] = [] for _, el, _, _ in self.elem_keys[Junction]: - self.descrip_vars['EJ'].append(sqf.to_numpy(el.get_value()) - / (2 * np.pi * unt.get_unit_freq())) + self.descrip_vars['EJ'].append( + sqf.to_numpy(el.get_value('Hz')) + / (2 * np.pi * unt.get_unit_freq()) + ) + + ## Inductive energies self.descrip_vars['EL'] = [] - self.descrip_vars['EL_incl'] = [] + self.descrip_vars['EL_has_ext_flux'] = [] for _, el, b_id in self.elem_keys[Inductor]: - self.descrip_vars['EL'].append(sqf.to_numpy(el.get_value(unt._unit_ind))) - self.descrip_vars['EL_incl'].append( - np.sum(np.abs(self.descrip_vars['B'][b_id, :])) != 0) - - self.descrip_vars['EC'] = dict() - for i in range(self.descrip_vars['har_dim'], self.descrip_vars['n_modes']): - for j in range(i, self.descrip_vars['n_modes']): - self.descrip_vars['EC'][(i,j)] = ( - (2 * unt.e) ** 2 - / (unt.hbar * 2 * np.pi * unt.get_unit_freq()) - * self.cInvTrans[i, j] - ) - if i == j: - self.descrip_vars['EC'][(i,j)] /= 2 + self.descrip_vars['EL'].append( + sqf.to_numpy(el.get_value('Hz')) / unt.get_unit_freq() + ) + # Check if the element has any external flux associated with it + # (corresponding row of B is nonzero) + self.descrip_vars['EL_has_ext_flux'].append( + np.sum(np.abs(self.descrip_vars['B'][b_id, :])) != 0 + ) - # values of loops - self.descrip_vars['n_loops'] = len(self.loops) - self.descrip_vars['loops'] = [sqf.to_numpy(self.loops[i].value()) / 2 / np.pi - for i in range(len(self.loops))] + ## Charging energies + self.descrip_vars['EC'] = ( + (2 * unt.e) **2 + / (unt.hbar * 2 * np.pi * unt.get_unit_freq()) + * self._get_quadratic_cInvTrans() + ) + + # Values of external flux and gate charge + self.descrip_vars['loops'] = [ + sqf.to_numpy(self.loops[i].value()) / 2 / np.pi + for i in range(self.descrip_vars['n_loops']) + ] + self.descrip_vars['ng'] = [ + self.charge_islands[i].value() + for i in range(self.descrip_vars['har_dim'], self.n) + ] def description( self, @@ -1219,17 +1222,20 @@ def description( _test: bool = False, ) -> Optional[str]: """Print out Hamiltonian and a listing of the modes (whether they are - harmonic or charge modes with the frequency for each harmonic mode), - Hamiltonian parameters, and external flux values. Constructs a - symbolic Hamiltonian, which is stored in `.descrip_vars['H']`. + harmonic or charge modes, with the frequency for each harmonic mode), + Hamiltonian parameters, and external flux values. + + Values printed in the description can be accessed using the + ``.descrip_vars`` dictionary, including a symbolic Hamiltonian in + ``.descrip_vars['H']``. Parameters ---------- tp: - If ``None`` prints out the output as Latex if SQcircuit is + If ``None`` prints out the output as LaTeX if SQcircuit is running in a Jupyter notebook and as text if SQcircuit is running in Python terminal. If ``tp`` is ``'ltx'``, - the output is in Latex format, and if ``tp`` is ``'txt'`` the + the output is in LaTeX format, and if ``tp`` is ``'txt'`` the output is in text format. _test: if True, return the entire description as string @@ -1247,7 +1253,7 @@ def description( else: txt = HamilTxt(tp) - self._compute_params() + self._compute_descrip_vars() self.descrip_vars['H'] = symbolic.construct_hamiltonian(self) final_txt = txt.print_circuit_description(self.descrip_vars) @@ -1388,6 +1394,21 @@ def _squeeze_op(self, op: Qobj) -> Qobj: op_sq.dims = self._get_op_dims() return op_sq + def _impedance(self, i: int) -> float: + """Compute the impedance of the ``i``th mode. + + Parameters + ---------- + i: + Index of the mode (starts from zero for the first mode). + + Returns + ---------- + The impedance Z of mode ``i``. + """ + return np.sqrt( + self.cInvTrans[i, i] / self.lTrans[i, i] + ) def _charge_op_isolated(self, i: int) -> Qobj: """Return charge operator for each isolated mode normalized by @@ -1397,7 +1418,7 @@ def _charge_op_isolated(self, i: int) -> Qobj: Parameters ---------- i: - Index of the mode. (starts from zero for the first mode) + Index of the mode (starts from zero for the first mode). Returns ---------- @@ -1407,10 +1428,8 @@ def _charge_op_isolated(self, i: int) -> Qobj: if self._is_charge_mode(i): ng = self.charge_islands[i].value() op = (2*unt.e/np.sqrt(unt.hbar)) * (qt.charge((self.m[i]-1)/2)-ng) - else: - Z = np.sqrt(self.cInvTrans[i, i] / self.lTrans[i, i]) - Q_zp = -1j * np.sqrt(0.5/Z) + Q_zp = -1j * np.sqrt(0.5 / self._impedance(i)) if self.m[i] == 1: # QuTiP 5.0.x no longer supports create/destroy with N = 1 op = Qobj([[0.]]).to('csr') @@ -1436,13 +1455,11 @@ def _flux_op_isolated(self, i: int) -> Qobj: if self._is_charge_mode(i): op = qt.qeye(self.m[i]) - else: - Z = np.sqrt(self.cInvTrans[i, i] / self.lTrans[i, i]) if self.m[i] == 1: op = Qobj([[0.]]).to('csr') else: - op = np.sqrt(0.5*Z) * ( + op = np.sqrt(0.5 * self._impedance(i)) * ( qt.destroy(self.m[i]) + qt.create(self.m[i]) ) @@ -1507,27 +1524,41 @@ def _d_op_isolated(self, i: int, w: float) -> Qobj: elif w > 0: return d.dag() - def alpha(self, i: Union[int, range], j: int) -> float: + def _phi_zp(self, i: int) -> float: + """Return the zero-point fluctuations for the ``i``th mode's + flux operator. + + Parameters + ---------- + i: + Index of the mode (starts from zero for the first mode). + + Returns + ---------- + The zero-point fluctuations for mode ``i``. + """ + return ( + (2 * np.pi / unt.Phi0) + * np.sqrt(unt.hbar / 2 * self._impedance(i)) + ) + + def alpha(self, i: int, j: int) -> float: """Return the alpha, amount of displacement, for the bosonic displacement operator for junction i and mode j. Parameters ---------- i: - Index of the Junction. (starts from zero for the first mode) + Index of the Junction (starts from zero for the first mode) j: - Index of the mode. (starts from zero for the first mode) + Index of the mode (starts from zero for the first mode). Returns ---------- The value of alpha for the junction ``i`` and mode ``j``. """ - Z = np.sqrt(self.cInvTrans[j, j] / self.lTrans[j, j]) - - coef = 2 * np.pi / unt.Phi0 * 1j - - return coef * np.sqrt(unt.hbar/2*Z) * self.wTrans[i, j] + return 1j * self._phi_zp(j) * self.wTrans[i, j] def _build_op_memory(self) -> None: """Build the charge, flux, number, and cross multiplication of charge @@ -1624,6 +1655,21 @@ def _build_exp_ops(self) -> None: self._memory_ops['exp'] = exp_ops self._memory_ops['root_exp'] = root_exp_ops + def _get_quadratic_cInvTrans(self) -> ndarray: + """Compute a modified ``cInvTrans`` for use when computing a + quadratic form like + ``0.5 Q^T @ cInvTrans @ Q`` + by summing over the upper/lower triangle. + + Returns + ---------- + ``self.cInvTrans`` with the diagonal divided by 2. + """ + q_cInvTrans = self.cInvTrans.copy() + q_cInvTrans[np.diag_indices_from(q_cInvTrans)] /= 2 + + return q_cInvTrans + def _get_LC_hamil(self) -> Qobj: """Construct the LC part of the circuit Hamiltonian. @@ -1633,19 +1679,21 @@ def _get_LC_hamil(self) -> Qobj: """ LC_hamil = 0 + q_cInvTrans = self._get_quadratic_cInvTrans() + for i in range(self.n): # we write j in this form because of "_memory_ops["QQ"]" shape for j in range(self.n - i): if j == 0: if self._is_charge_mode(i): - LC_hamil += (0.5 * self.cInvTrans[i, i] + LC_hamil += (q_cInvTrans[i, i] * self._memory_ops['QQ'][i][j]) else: LC_hamil += self.omega[i] * self._memory_ops['N'][i] elif j > 0: if self.cInvTrans[i, i + j] != 0: - LC_hamil += (self.cInvTrans[i, i + j] + LC_hamil += (q_cInvTrans[i, i + j] * self._memory_ops['QQ'][i][j]) return LC_hamil diff --git a/SQcircuit/symbolic.py b/SQcircuit/symbolic.py index 03ca553..ae1233f 100755 --- a/SQcircuit/symbolic.py +++ b/SQcircuit/symbolic.py @@ -1,20 +1,11 @@ +"""symbolic.py""" import numpy as np import sympy as sm -from sympy import ( - Expr, - Symbol -) +from sympy import Expr, Symbol from sympy.physics.quantum import Operator -from SQcircuit.elements import ( - Capacitor, - Inductor, - Junction, - Loop, - Charge -) +from SQcircuit.elements import Inductor, Junction -from sympy.printing.latex import LatexPrinter class ExplicitSymbol(Symbol): def __new__(cls, plainname, latexname): @@ -29,8 +20,8 @@ def _latex(self, printer): class qOperator(Operator): """ Class extending sympy physics operator, mostly for custom printing purposes. - Created with a base `name` and a `subscript`. When printing in LaTeX, a hat - will be put over the name, but not in plaintext. + Created with a base ``name` and a `subscript``. When printing in LaTeX, a + hat will be put over the name, but not in plaintext. """ def __new__(cls, name, subscript=None): if subscript: @@ -50,58 +41,71 @@ def _latex(self, printer): tex += fr'_{{{self.sub}}}' return printer.doprint(Symbol(tex)) + def phi_op(i): return ExplicitSymbol(f'phi_{i}', fr'\hat{{\varphi}}_{{{i}}}') + def a(i: int) -> qOperator: return qOperator('a', i) + def ad(i: int) -> qOperator: return qOperator('ad', i) + def n(i: int) -> qOperator: return qOperator('n', i) + def phi_ext(i: int) -> ExplicitSymbol: return ExplicitSymbol(f'phi_e{i}', r'\varphi_{\text{ext}_{' + str(i) + '}}') + def phi_zp(i: int) -> ExplicitSymbol: return ExplicitSymbol(f'zp_{i}', r'\varphi_{\text{zp}_{' + str(i) + '}}') + def omega(i: int) -> Symbol: return Symbol(rf'omega_{i}') + def ng(i: int) -> ExplicitSymbol: return ExplicitSymbol(f'ng_{i}', fr'n_{{g_{{{i}}}}}') + def EC(i: int, j: int) -> ExplicitSymbol: return ExplicitSymbol(f'EC_{i}{j}', fr'E_{{C_{{{i}{j}}}}}') + def EL(i: int) -> ExplicitSymbol: return ExplicitSymbol(f'EL_{i}', fr'E_{{L_{{{i}}}}}') + def EJ(i: int) -> ExplicitSymbol: return ExplicitSymbol(f'EJ_{i}', fr'E_{{J_{{{i}}}}}') + H = qOperator('H') -def har_mode_hamil(coeff_dict, do_sum=True): - terms = [sm.Mul(omega(i+1), ad(i+1), a(i+1), evaluate=False) \ +def har_mode_hamil(coeff_dict, do_sum=True) -> Expr: + terms = [sm.Mul(omega(i+1), ad(i+1), a(i+1), evaluate=False) for i in range(coeff_dict['har_dim'])] if do_sum: return sm.Add(*terms, evaluate=False) - else: - return terms -def charge_mode_hamil(coeff_dict, do_sum=True): + return terms + + +def charge_mode_hamil(coeff_dict, do_sum=True) -> Expr: hamil = 0 for i in range(coeff_dict['har_dim'], coeff_dict['n_modes']): for j in range(i, coeff_dict['n_modes']): @@ -109,16 +113,17 @@ def charge_mode_hamil(coeff_dict, do_sum=True): if do_sum: return hamil - else: - return sm.Add.make_args(hamil) -def inductive_hamil(elem_keys, coeff_dict, do_sum=True): + return sm.Add.make_args(hamil) + + +def inductive_hamil(elem_keys, coeff_dict, do_sum=True) -> Expr: S = coeff_dict['S'] B = coeff_dict['B'] terms = [] - for i, (edge, el, B_idx) in enumerate(elem_keys[Inductor]): - if np.sum(np.abs(B[B_idx, :])) == 0: + for i, (edge, _, b_id) in enumerate(elem_keys[Inductor]): + if np.sum(np.abs(B[b_id, :])) == 0: continue if 0 in edge: w = S[edge[0] + edge[1] - 1, :] @@ -127,7 +132,7 @@ def inductive_hamil(elem_keys, coeff_dict, do_sum=True): w = np.round(w[:coeff_dict['har_dim']], 3) w_sym = sm.Matrix(w) - B_sym = sm.Matrix(B[B_idx, :]) + B_sym = sm.Matrix(B[b_id, :]) phis = [phi_op(i+1) for i in range(len(w_sym))] phi_exts = [phi_ext(i+1) for i in range(len(B_sym))] @@ -136,16 +141,17 @@ def inductive_hamil(elem_keys, coeff_dict, do_sum=True): sm.nsimplify(B_sym.dot(phi_exts)), evaluate=False)) if do_sum: return sum(terms) - else: - return terms -def jj_hamil(elem_keys, coeff_dict, do_sum=True): + return terms + + +def jj_hamil(elem_keys, coeff_dict, do_sum=True) -> Expr: terms = [] - for i, (edge, el, B_idx, W_idx) in enumerate(elem_keys[Junction]): - W_sym = sm.Matrix(coeff_dict['W'][W_idx, :]) + for i, (_, _, b_id, w_id) in enumerate(elem_keys[Junction]): + W_sym = sm.Matrix(coeff_dict['W'][w_id, :]) phis = [phi_op(i+1) for i in range(len(W_sym))] - B_sym = sm.Matrix(coeff_dict['B'][B_idx, :]) + B_sym = sm.Matrix(coeff_dict['B'][b_id, :]) phi_exts = [phi_ext(i+1) for i in range(len(B_sym))] terms.append(EJ(i+1) * sm.cos(sm.Add(sm.nsimplify(W_sym.dot(phis)), @@ -154,10 +160,11 @@ def jj_hamil(elem_keys, coeff_dict, do_sum=True): ) if do_sum: return sum(terms) - else: - return terms -def construct_hamiltonian(cr): + return terms + + +def construct_hamiltonian(cr) -> Expr: har_terms = har_mode_hamil(cr.descrip_vars, do_sum=False) charge_terms = charge_mode_hamil(cr.descrip_vars, do_sum=False) ind_terms = inductive_hamil(cr.elem_keys, cr.descrip_vars, do_sum=False) diff --git a/SQcircuit/tests/data/descriptions/zero_pi_inductors_ltx.txt b/SQcircuit/tests/data/descriptions/zero_pi_inductors_ltx.txt index a2d6c68..348ae9b 100755 --- a/SQcircuit/tests/data/descriptions/zero_pi_inductors_ltx.txt +++ b/SQcircuit/tests/data/descriptions/zero_pi_inductors_ltx.txt @@ -1,4 +1,4 @@ -\hat{H} = \omega_{1} \hat{a}^\dagger_{1} \hat{a}_{1} + \omega_{2} \hat{a}^\dagger_{2} \hat{a}_{2} + E_{C_{33}} \left(- n_{g_{3}} + \hat{n}_{3}\right)^{2} + E_{L_{1}} \left(- \hat{\varphi}_{1} - \hat{\varphi}_{2}\right) \left(- \varphi_{\text{ext}_{1}} / 2\right) + E_{L_{2}} \left(- \hat{\varphi}_{1} + \hat{\varphi}_{2}\right) \left(- \varphi_{\text{ext}_{1}} / 2\right) + E_{J_{1}} \cos{\left(- \hat{\varphi}_{3} + \hat{\varphi}_{1} \right)} + E_{J_{2}} \cos{\left(\hat{\varphi}_{1} + \hat{\varphi}_{3} \right)} +\hat{H} = \omega_{1} \hat{a}^\dagger_{1} \hat{a}_{1} + \omega_{2} \hat{a}^\dagger_{2} \hat{a}_{2} + E_{C_{33}} \left(- n_{g_{3}} + \hat{n}_{3}\right)^{2} + E_{L_{1}} \left(- \hat{\varphi}_{1} - \hat{\varphi}_{2}\right) \left(- \varphi_{\text{ext}_{1}} / 2\right) + E_{L_{2}} \left(- \hat{\varphi}_{1} + \hat{\varphi}_{2}\right) \left(- \varphi_{\text{ext}_{1}} / 2\right) + E_{J_{1}} \cos{\left(\hat{\varphi}_{1} + \hat{\varphi}_{3} \right)} + E_{J_{2}} \cos{\left(- \hat{\varphi}_{3} + \hat{\varphi}_{1} \right)} ------------------------------------------------------------ \text{mode 1:}~~~~~~~~~~~\text{harmonic}~~~~~~~~~~~\hat{\varphi}_{1} = \varphi_{\text{zp}_{1}} \left(\hat{a}_{1} + \hat{a}^\dagger_{1}\right)~~~~~~~~~~~\omega_{1} / 2 \pi = 3.22489~~~~~~~~~~~\varphi_{\text{zp}_{1}} = 2.49 \text{mode 2:}~~~~~~~~~~~\text{harmonic}~~~~~~~~~~~\hat{\varphi}_{2} = \varphi_{\text{zp}_{2}} \left(\hat{a}_{2} + \hat{a}^\dagger_{2}\right)~~~~~~~~~~~\omega_{2} / 2 \pi = 0.39497~~~~~~~~~~~\varphi_{\text{zp}_{2}} = 0.87 diff --git a/SQcircuit/tests/data/descriptions/zero_pi_inductors_txt.txt b/SQcircuit/tests/data/descriptions/zero_pi_inductors_txt.txt index b65c2c8..698510b 100755 --- a/SQcircuit/tests/data/descriptions/zero_pi_inductors_txt.txt +++ b/SQcircuit/tests/data/descriptions/zero_pi_inductors_txt.txt @@ -2,7 +2,7 @@ H = ω₁⋅ad₁⋅a₁ + ω₂⋅ad₂⋅a₂ + EC₃₃⋅(-ng₃ + n₃) + EL₁⋅(-φ₁ - φ₂)⋅-φₑ₁/2 + EL₂⋅(- -φ₁ + φ₂)⋅-φₑ₁/2 + EJ₁⋅cos(-φ₃ + φ₁) + EJ₂⋅cos(φ₁ + φ₃) +φ₁ + φ₂)⋅-φₑ₁/2 + EJ₁⋅cos(φ₁ + φ₃) + EJ₂⋅cos(-φ₃ + φ₁) ------------------------------------------------------------ mode 1: harmonic φ₁ = zp₁⋅(a₁ + ad₁) ω₁/2⋅π = 3.22489 zp₁ = 2.49 mode 2: harmonic φ₂ = zp₂⋅(a₂ + ad₂) ω₂/2⋅π = 0.39497 zp₂ = 0.87 diff --git a/SQcircuit/tests/data/descriptions/zero_pi_junctions_ltx.txt b/SQcircuit/tests/data/descriptions/zero_pi_junctions_ltx.txt index cbedab0..1691b99 100755 --- a/SQcircuit/tests/data/descriptions/zero_pi_junctions_ltx.txt +++ b/SQcircuit/tests/data/descriptions/zero_pi_junctions_ltx.txt @@ -1,4 +1,4 @@ -\hat{H} = \omega_{1} \hat{a}^\dagger_{1} \hat{a}_{1} + \omega_{2} \hat{a}^\dagger_{2} \hat{a}_{2} + E_{C_{33}} \left(- n_{g_{3}} + \hat{n}_{3}\right)^{2} + E_{J_{1}} \cos{\left(\varphi_{\text{ext}_{1}} / 2 - \hat{\varphi}_{1} + \hat{\varphi}_{3} \right)} + E_{J_{2}} \cos{\left(- \varphi_{\text{ext}_{1}} / 2 + \hat{\varphi}_{1} + \hat{\varphi}_{3} \right)} +\hat{H} = \omega_{1} \hat{a}^\dagger_{1} \hat{a}_{1} + \omega_{2} \hat{a}^\dagger_{2} \hat{a}_{2} + E_{C_{33}} \left(- n_{g_{3}} + \hat{n}_{3}\right)^{2} + E_{J_{1}} \cos{\left(- \varphi_{\text{ext}_{1}} / 2 + \hat{\varphi}_{1} + \hat{\varphi}_{3} \right)} + E_{J_{2}} \cos{\left(\varphi_{\text{ext}_{1}} / 2 - \hat{\varphi}_{1} + \hat{\varphi}_{3} \right)} ------------------------------------------------------------ \text{mode 1:}~~~~~~~~~~~\text{harmonic}~~~~~~~~~~~\hat{\varphi}_{1} = \varphi_{\text{zp}_{1}} \left(\hat{a}_{1} + \hat{a}^\dagger_{1}\right)~~~~~~~~~~~\omega_{1} / 2 \pi = 3.22489~~~~~~~~~~~\varphi_{\text{zp}_{1}} = 2.49 \text{mode 2:}~~~~~~~~~~~\text{harmonic}~~~~~~~~~~~\hat{\varphi}_{2} = \varphi_{\text{zp}_{2}} \left(\hat{a}_{2} + \hat{a}^\dagger_{2}\right)~~~~~~~~~~~\omega_{2} / 2 \pi = 0.39497~~~~~~~~~~~\varphi_{\text{zp}_{2}} = 0.87 diff --git a/SQcircuit/tests/data/descriptions/zero_pi_junctions_txt.txt b/SQcircuit/tests/data/descriptions/zero_pi_junctions_txt.txt index 482c5bb..b13c0b6 100755 --- a/SQcircuit/tests/data/descriptions/zero_pi_junctions_txt.txt +++ b/SQcircuit/tests/data/descriptions/zero_pi_junctions_txt.txt @@ -1,8 +1,8 @@ 2 -H = ω₁⋅ad₁⋅a₁ + ω₂⋅ad₂⋅a₂ + EC₃₃⋅(-ng₃ + n₃) + EJ₁⋅cos(φₑ₁/2 - φ₁ + φ₃) + EJ₂ +H = ω₁⋅ad₁⋅a₁ + ω₂⋅ad₂⋅a₂ + EC₃₃⋅(-ng₃ + n₃) + EJ₁⋅cos(-φₑ₁/2 + φ₁ + φ₃) + EJ -⋅cos(-φₑ₁/2 + φ₁ + φ₃) +₂⋅cos(φₑ₁/2 - φ₁ + φ₃) ------------------------------------------------------------ mode 1: harmonic φ₁ = zp₁⋅(a₁ + ad₁) ω₁/2⋅π = 3.22489 zp₁ = 2.49 mode 2: harmonic φ₂ = zp₂⋅(a₂ + ad₂) ω₂/2⋅π = 0.39497 zp₂ = 0.87 diff --git a/SQcircuit/tests/test_paper.py b/SQcircuit/tests/test_papers.py similarity index 99% rename from SQcircuit/tests/test_paper.py rename to SQcircuit/tests/test_papers.py index f7a299e..3b61c7a 100644 --- a/SQcircuit/tests/test_paper.py +++ b/SQcircuit/tests/test_papers.py @@ -17,7 +17,7 @@ PHASEDIR = os.path.join(TESTDIR, 'data', 'phase_coord') SPECDIR = os.path.join(TESTDIR, 'data', 'spectra') -def test_paper_1_a() -> None: +def _test_paper_1_a() -> None: target_freqs = np.array([0.69721321, 1.43920434, 1.90074875, 2.05449921]) sq.set_engine('NumPy') diff --git a/SQcircuit/texts.py b/SQcircuit/texts.py index 7e9ddbe..dfc9eb0 100755 --- a/SQcircuit/texts.py +++ b/SQcircuit/texts.py @@ -1,29 +1,19 @@ -""" Module for text and latex output of SQcircuit""" +"""Module for text and LaTeX output of SQcircuit""" import numpy as np -from IPython.display import display, Latex, Math +from IPython.display import display, Math import sympy as sm from sympy.printing.latex import LatexPrinter from sympy.printing.pretty.pretty import PrettyPrinter from sympy.printing.pretty.stringpict import prettyForm -from sympy.printing.pretty.stringpict import prettyForm, stringPict -from typing import List - -from SQcircuit.elements import ( - Element, - Capacitor, - Inductor, - Junction, - Loop, - Charge -) + +from SQcircuit.elements import Inductor, Junction import SQcircuit.symbolic as sym def is_notebook(): - """ - The function that checks whether we are working in notebook environment or - Python terminal. + """Checks whether we are working in notebook environment or the Python + terminal. """ try: shell = get_ipython().__class__.__name__ @@ -37,9 +27,8 @@ def is_notebook(): return False # Probably standard Python interpreter -def newDiv(self, den, slashed=False): - """ - prettyForm division to always be inline x/y +def newDiv(self, den): + """prettyForm division to always be inline x/y """ return prettyForm(binding=prettyForm.DIV, s=self.s + '/' + den.s) @@ -54,10 +43,10 @@ def __init__(self, tp='ltx'): self.tp = tp self.line = 60 * "-" + "\n" if tp == 'ltx': - self.printer = LatexPrinter(dict(order='none', fold_short_frac=True)) + self.printer = LatexPrinter({'order': 'none', 'fold_short_frac': True}) elif tp == 'txt': prettyForm.__truediv__ = newDiv - self.printer = PrettyPrinter(dict(use_unicode=False, order='none')) + self.printer = PrettyPrinter({'use_unicode': False, 'order': 'none'}) else: raise ValueError('Permitted values for `tp` are \'ltx\' and ' '\'txt\'.') @@ -85,49 +74,62 @@ def mode_txt(self, coeff_dict): kind = 'harmonic' omega_val = np.round(coeff_dict['omega'][i], 5) zp_val = np.round(coeff_dict['phi_zp'][i], 2) - info = [sm.Eq(sym.phi_op(i+1), - sym.phi_zp(i+1)*(sym.a(i+1) + sym.ad(i+1))), - sm.Eq(sym.omega(i+1)/(2 * sm.pi), omega_val), + info = [sm.Eq(sym.phi_op(i+1), + sym.phi_zp(i+1)*(sym.a(i+1) + sym.ad(i+1))), + sm.Eq(sym.omega(i+1)/(2 * sm.pi), omega_val), sm.Eq(sym.phi_zp(i+1), zp_val)] else: kind = 'charge' ng_val = np.round(coeff_dict['ng'][i - coeff_dict['har_dim']], 3) info = [sm.Eq(sym.ng(i+1), ng_val)] - txt += self.plaintxt(f'mode {i+1}:') + self.tab() \ - + self.plaintxt(kind) + self.tab() + txt += ( + self.plaintxt(f'mode {i+1}:') + self.tab() + + self.plaintxt(kind) + self.tab() + ) txt += self.tab().join([self.printer.doprint(e) for e in info]) + '\n' return txt + '\n' def param_txt(self, coeff_dict): params = [] - for key in coeff_dict['EC']: - i,j = key - params.append(sm.Eq(sym.EC(i+1,j+1), - np.round(coeff_dict['EC'][(i,j)], 3))) + # For the LC part, only print out if not completely absorbed into the + # omegas of the harmonic modes. + ## Only capacitve energies on charge modes + for i in range(coeff_dict['har_dim'], coeff_dict['n_modes']): + for j in range(i, coeff_dict['n_modes']): + params.append(sm.Eq(sym.EC(i+1,j+1), + np.round(coeff_dict['EC'][i, j], 3))) + ## Only inductive energy if it has external flux assigned to it for i, EL_val in enumerate(coeff_dict['EL']): - if coeff_dict['EL_incl'][i]: + if coeff_dict['EL_has_ext_flux'][i]: params.append(sm.Eq(sym.EL(i+1), np.round(EL_val, 3))) + # Need to print all junctions for i, EJ_val in enumerate(coeff_dict['EJ']): params.append(sm.Eq(sym.EJ(i+1), np.round(EJ_val, 3))) - txt = self.plaintxt('parameters:') + self.tab() \ - + self.tab().join([self.printer.doprint(p) for p in params]) + txt = ( + self.plaintxt('parameters:') + self.tab() + + self.tab().join([self.printer.doprint(p) for p in params]) + ) return txt + '\n' def loop_txt(self, coeff_dict): - info = [sm.Eq(sym.phi_ext(i+1)/(2*sm.pi), - np.round(coeff_dict['loops'][i], 2)) \ + info = [sm.Eq(sym.phi_ext(i+1)/(2*sm.pi), + np.round(coeff_dict['loops'][i], 2)) for i in range(coeff_dict['n_loops'])] - txt = self.plaintxt('loops:') + self.tab() \ - + self.tab().join([self.printer.doprint(l) for l in info]) + txt = ( + self.plaintxt('loops:') + self.tab() + + self.tab().join([self.printer.doprint(l) for l in info]) + ) return txt def print_circuit_description(self, coeff_dict): - finalTxt = self.ham_txt(coeff_dict) + self.line \ - + self.mode_txt(coeff_dict) + self.line \ - + self.param_txt(coeff_dict) \ - + self.loop_txt(coeff_dict) + finalTxt = ( + self.ham_txt(coeff_dict) + self.line + + self.mode_txt(coeff_dict) + self.line + + self.param_txt(coeff_dict) + + self.loop_txt(coeff_dict) + ) self.display(finalTxt) return finalTxt @@ -161,11 +163,11 @@ def print_loop_description(cr): for i in range(cr.B.shape[0]): el = None - for _, el_ind, B_idx in cr.elem_keys[Inductor]: - if i == B_idx: + for _, el_ind, b_id in cr.elem_keys[Inductor]: + if i == b_id: el = el_ind - for _, el_ind, B_idx, W_idx in cr.elem_keys[Junction]: - if i == B_idx: + for _, el_ind, b_id, _ in cr.elem_keys[Junction]: + if i == b_id: el = el_ind id_str = el.id_str diff --git a/SQcircuit/torch_extensions.py b/SQcircuit/torch_extensions.py index e3c5b33..36147ad 100755 --- a/SQcircuit/torch_extensions.py +++ b/SQcircuit/torch_extensions.py @@ -238,7 +238,7 @@ def backward(ctx, grad_output) -> Tuple[Tensor]: grad_output_eigenvalue = grad_output[:, 0] grad_output_eigenvector = grad_output[:, 1:] - partial_omega = torch.zeros([m, n], dtype=float) + partial_omega = torch.zeros([m, n], dtype=torch.float64) partial_eigenvec = torch.zeros([m, n, l], dtype=torch.complex128) for el_idx in range(m): for eigen_idx in range(n):