scqubits - pypi Package Compare versions

Comparing version 3.1.0 to 3.1.1

optional-requirements.txt

 h5py>=2.10
 pytest
 ipywidgets
-pathos
+pathos>=0.3.0
 typing_extensions

PKG-INFO

 Metadata-Version: 2.1
 Name: scqubits
-Version: 3.1.0
+Version: 3.1.1
 Summary: scqubits: superconducting qubits in Python
@@ -5,0 +5,0 @@ Home-page: https://scqubits.readthedocs.io

scqubits.egg-info/PKG-INFO

 Metadata-Version: 2.1
 Name: scqubits
-Version: 3.1.0
+Version: 3.1.1
 Summary: scqubits: superconducting qubits in Python
@@ -5,0 +5,0 @@ Home-page: https://scqubits.readthedocs.io

scqubits.egg-info/SOURCES.txt

@@ -79,3 +79,2 @@ LICENSE
 scqubits/tests/test_zeropi.py
-scqubits/tests/data/.fullzeropi_1.hdf5
 scqubits/tests/data/circuit_DFC.yaml
@@ -82,0 +81,0 @@ scqubits/tests/data/circuit_fluxonium.yaml

scqubits/core/central_dispatch.py

@@ -17,3 +17,2 @@ # central_dispatch.py
 import weakref
-
 from types import MethodType
@@ -20,0 +19,0 @@ from weakref import WeakKeyDictionary

scqubits/core/circuit_utils.py

@@ -15,3 +15,3 @@ # circuit-utils.py
 
-from typing import TYPE_CHECKING, Any, Callable, Dict, List, Union
+from typing import TYPE_CHECKING, Any, Callable, List, Union
 
@@ -245,3 +245,3 @@ import numpy as np
     matrix = sparse.dia_matrix(
-        (np.ones(dim_theta), [1]),
+        (np.ones(dim_theta), [-1]),
         shape=(dim_theta, dim_theta),
@@ -258,3 +258,3 @@ ).tocsc()
     matrix = sparse.dia_matrix(
-        (np.ones(dim_theta), [-1]),
+        (np.ones(dim_theta), [1]),
         shape=(dim_theta, dim_theta),
@@ -261,0 +261,0 @@ ).tocsc()

scqubits/core/cos2phi_qubit.py

@@ -26,3 +26,3 @@ # cos2phi_qubit.py
 from scipy import sparse
-from scipy.sparse import csc_matrix, dia_matrix
+from scipy.sparse import coo_matrix, csc_matrix, dia_matrix
 
@@ -38,3 +38,2 @@ import scqubits.core.constants as constants
 import scqubits.io_utils.fileio_serializers as serializers
-import scqubits.settings as settings
 import scqubits.utils.plotting as plot
@@ -47,22 +46,32 @@ import scqubits.utils.spectrum_utils as utils
 
-# -Cosine two phi qubit noise class
+# - Cosine-2-phi qubit noise class ------------------------------------------------------------------------------------
 class NoisyCos2PhiQubit(NoisySystem, ABC):
     @abstractmethod
-    def phi_1_operator(self) -> csc_matrix:
+    def phi_1_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
         pass
 
     @abstractmethod
-    def phi_2_operator(self) -> csc_matrix:
+    def phi_2_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
         pass
 
     @abstractmethod
-    def n_1_operator(self) -> csc_matrix:
+    def n_1_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
         pass
 
     @abstractmethod
-    def n_2_operator(self) -> csc_matrix:
+    def n_2_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
         pass
 
     @abstractmethod
-    def n_zeta_operator(self) -> csc_matrix:
+    def n_zeta_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
         pass
@@ -607,7 +616,28 @@
 
-    def phi_operator(self) -> csc_matrix:
-        r"""Returns :math:`\phi` operator"""
-        return self._kron3(
+    def phi_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Returns the :math:`\phi` operator in the native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns :math:`\phi` operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns :math:`\phi` operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns :math:`\phi` operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            :math:`\phi` operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\phi` operator has dimensions of truncated_dim
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen, :math:`\phi`
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
+        native = self._kron3(
             self._phi_operator(), self._identity_zeta(), self._identity_theta()
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -626,7 +656,29 @@ def _n_phi_operator(self) -> csc_matrix:
 
-    def n_phi_operator(self) -> csc_matrix:
-        r"""Returns :math:`n_\phi` operator"""
-        return self._kron3(
+    # changed expected from csc_matrix to Union[ndarray, coo_matrix]
+    def n_phi_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, coo_matrix]:
+        r"""
+        Returns the :math:`n_\phi` operator in the harmonic oscillator or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns :math:`n_\phi` operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns :math:`n_\phi` operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns :math:`n_\phi` operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            :math:`n_\phi` operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`n_\phi` operator has dimensions of truncated_dim
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen, :math:`n_\phi`
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray :math:`\zeta`.
+        """
+        native = self._kron3(
             self._n_phi_operator(), self._identity_zeta(), self._identity_theta()
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -645,7 +697,28 @@ def _zeta_operator(self) -> csc_matrix:
 
-    def zeta_operator(self) -> csc_matrix:
-        r"""Returns :math:`\zeta` operator"""
-        return self._kron3(
+    def zeta_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Returns the :math:`\zeta`  operator in the harmonic oscillator or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns :math:`\zeta` operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns :math:`\zeta`  operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns :math:`\zeta`  operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            :math:`\zeta`  operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\zeta`  operator has dimensions of truncated_dim
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen, :math:`\zeta`
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
+        native = self._kron3(
             self._identity_phi(), self._zeta_operator(), self._identity_theta()
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -664,7 +737,28 @@ def _n_zeta_operator(self) -> csc_matrix:
 
-    def n_zeta_operator(self) -> csc_matrix:
-        r"""Returns :math:`n_\zeta` operator"""
-        return self._kron3(
+    def n_zeta_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Returns the :math:`n_\zeta`  operator in the harmonic oscillator or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns :math:`n_\zeta` operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns :math:`n_\zeta`  operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns :math:`n_\zeta`  operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            :math:`n_\zeta`  operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`n_\zeta`  operator has dimensions of truncated_dim
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen, :math:`n_\zeta`
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
+        native = self._kron3(
             self._identity_phi(), self._n_zeta_operator(), self._identity_theta()
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -707,7 +801,28 @@ def _exp_i_phi_operator(self) -> csc_matrix:
 
-    def n_theta_operator(self) -> csc_matrix:
-        r"""Returns :math:`n_\theta` operator"""
-        return self._kron3(
+    def n_theta_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Returns the :math:`n_\theta`  operator in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns :math:`n_\theta` operator in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns :math:`n_\theta`  operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns :math:`n_\theta`  operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            :math:`n_\theta`  operator in chosen basis. If charge basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`n_\theta`  operator has dimensions of truncated_dim
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen, :math:`n_\theta`
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
+        native = self._kron3(
             self._identity_phi(), self._identity_zeta(), self._n_theta_operator()
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -737,3 +852,3 @@ def _cos_theta_operator(self) -> csc_matrix:
             * sparse.dia_matrix(
-                (np.ones(self._dim_theta()), [1]),
+                (np.ones(self._dim_theta()), [-1]),
                 shape=(self._dim_theta(), self._dim_theta()),
@@ -745,3 +860,3 @@ ).tocsc()
             * sparse.dia_matrix(
-                (np.ones(self._dim_theta()), [-1]),
+                (np.ones(self._dim_theta()), [1]),
                 shape=(self._dim_theta(), self._dim_theta()),
@@ -785,7 +900,24 @@ ).tocsc()
 
-    def hamiltonian(self) -> csc_matrix:
+    def hamiltonian(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> csc_matrix:
+        """Returns Hamiltonian in basis obtained by employing harmonic basis for
+        :math:`\\phi, \\zeta` and charge basis for :math:`\\theta` or in the eigenenerg basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns Hamiltonian in basis obtained by employing harmonic basis for
+            :math:`\\phi, \\zeta` and charge basis for :math:`\\theta`.
+            If `True`, the energy eigenspectrum is computed, returns Hamiltonian in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns Hamiltonian in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Hamiltonian in chosen basis as csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, the Hamiltonian has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, Hamiltonian has dimensions of m x m,
+            for m given eigenvectors.
         """
-        Returns Hamiltonian in basis obtained by employing harmonic basis for
-        :math:`\\phi, \\zeta` and charge basis for :math:`\\theta`.
-        """
         phi_osc_mat = self._kron3(
@@ -858,3 +990,3 @@ op.number_sparse(self._dim_phi(), self.phi_plasma()),
 
-        return (
+        hamiltonian_mat = (
             phi_osc_mat
@@ -868,2 +1000,6 @@ + zeta_osc_mat
         )
+        native = hamiltonian_mat.tocsc()
+        return self.process_hamiltonian(
+            native_hamiltonian=native, energy_esys=energy_esys
+        )
 
@@ -1093,23 +1229,129 @@ def _evals_calc(self, evals_count) -> ndarray:
 
-    def phi_1_operator(self) -> csc_matrix:
-        """Returns operator representing the phase across inductor 1"""
-        return self.zeta_operator() - self.phi_operator()
+    def phi_1_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        """
+        Returns operator representing the phase across inductor 1 in harmonic oscillator or eigenenergy basis.
 
-    def phi_2_operator(self) -> csc_matrix:
-        """Returns operator representing the phase across inductor 2"""
-        return -self.zeta_operator() - self.phi_operator()
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the harmonic oscillator basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
 
-    def n_1_operator(self) -> csc_matrix:
-        """Returns operator representing the charge difference across junction 1"""
-        return 0.5 * self.n_phi_operator() + 0.5 * (
+        Returns
+        -------
+            Operator in chosen basis. If harmonic oscillator basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
+        native = self.zeta_operator() - self.phi_operator()
+        return self.process_op(native_op=native, energy_esys=energy_esys)
+
+    def phi_2_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        """
+        Returns operator representing the phase across inductor 2 in harmonic oscillator or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the harmonic oscillator basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If harmonic oscillator basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
+        native = -self.zeta_operator() - self.phi_operator()
+        return self.process_op(native_op=native, energy_esys=energy_esys)
+
+    def n_1_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        """
+        Returns operator representing the charge difference across junction 1 in native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the harmonic oscillator basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
+        native = 0.5 * self.n_phi_operator() + 0.5 * (
             self.n_theta_operator() - self.n_zeta_operator()
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def n_2_operator(self) -> csc_matrix:
-        """Returns operator representing the charge difference across junction 2"""
-        return 0.5 * self.n_phi_operator() - 0.5 * (
+    def n_2_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        """
+        Returns operator representing the charge difference across junction 2 in native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the harmonic oscillator basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
+        native = 0.5 * self.n_phi_operator() - 0.5 * (
             self.n_theta_operator() - self.n_zeta_operator()
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def d_hamiltonian_d_flux(self) -> csc_matrix:
+    def d_hamiltonian_d_flux(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        flux in the native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
         phi_flux_term = self._sin_phi_operator() * np.cos(
@@ -1139,5 +1381,28 @@ self.flux * np.pi
         )
-        return junction_mat + dis_junction_mat
+        native = junction_mat + dis_junction_mat
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def d_hamiltonian_d_EJ(self) -> csc_matrix:
+    def d_hamiltonian_d_EJ(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        EJ in the native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
         phi_flux_term = self._cos_phi_operator() * np.cos(
@@ -1160,6 +1425,29 @@ self.flux * np.pi
         )
-        return junction_mat + dis_junction_mat
+        native = junction_mat + dis_junction_mat
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def d_hamiltonian_d_ng(self) -> csc_matrix:
-        return (
+    def d_hamiltonian_d_ng(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        ng in the native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
+        """
+        native = (
             4 * self.dCJ * self._disordered_ecj() * self.n_phi_operator()
@@ -1170,1 +1458,2 @@ - 4
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)

scqubits/core/discretization.py

@@ -27,2 +27,3 @@ # discretization.py
 
+
 FIRST_STENCIL_COEFFS: Dict[int, List[float]] = {
@@ -29,0 +30,0 @@ 3: [-1 / 2, 0.0, 1 / 2],

scqubits/core/flux_qubit.py

@@ -16,3 +16,3 @@ # flux_qubit.py
 from abc import ABC, abstractmethod
-from typing import Any, Dict, List, Optional, Tuple
+from typing import Any, Dict, List, Optional, Tuple, Union
 
@@ -41,11 +41,17 @@ import numpy as np
     @abstractmethod
-    def d_hamiltonian_d_EJ1(self) -> ndarray:
+    def d_hamiltonian_d_EJ1(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
         pass
 
     @abstractmethod
-    def d_hamiltonian_d_EJ2(self) -> ndarray:
+    def d_hamiltonian_d_EJ2(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
         pass
 
     @abstractmethod
-    def d_hamiltonian_d_EJ3(self) -> ndarray:
+    def d_hamiltonian_d_EJ3(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
         pass
@@ -408,3 +414,3 @@
         evals = sp.linalg.eigh(
-            hamiltonian_mat, eigvals=(0, evals_count - 1), eigvals_only=True
+            hamiltonian_mat, subset_by_index=(0, evals_count - 1), eigvals_only=True
         )
@@ -416,3 +422,3 @@ return np.sort(evals)
         evals, evecs = sp.linalg.eigh(
-            hamiltonian_mat, eigvals=(0, evals_count - 1), eigvals_only=False
+            hamiltonian_mat, subset_by_index=(0, evals_count - 1), eigvals_only=False
         )
@@ -507,25 +513,106 @@ evals, evecs = spec_utils.order_eigensystem(evals, evecs)
 
-    def hamiltonian(self) -> ndarray:
-        """Return Hamiltonian in basis obtained by employing charge basis for both
-        degrees of freedom"""
-        return self.kineticmat() + self.potentialmat()
+    def hamiltonian(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Return Hamiltonian in the basis obtained by employing charge basis for both
+        degrees of freedom or in the eigenenergy basis.
 
-    def d_hamiltonian_d_EJ1(self) -> ndarray:
-        """Returns operator representing a derivative of the Hamiltonian with
-        respect to EJ1."""
-        return -0.5 * np.kron(
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns Hamiltonian in the basis obtained by employing charge basis for both degrees of freedom.
+            If `True`, the energy eigenspectrum is computed, returns Hamiltonian in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns Hamiltonian in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Hamiltonian in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, the Hamiltonian has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, Hamiltonian has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        hamiltonian_mat = self.kineticmat() + self.potentialmat()
+        return self.process_hamiltonian(
+            native_hamiltonian=hamiltonian_mat, energy_esys=energy_esys
+        )
+
+    def d_hamiltonian_d_EJ1(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        EJ1 in the native Hamiltonian basis or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native Hamiltonian basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        native = -0.5 * np.kron(
             self._exp_i_phi_operator() + self._exp_i_phi_operator().T, self._identity()
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def d_hamiltonian_d_EJ2(self) -> ndarray:
-        """Returns operator representing a derivative of the Hamiltonian with
-        respect to EJ2."""
-        return -0.5 * np.kron(
+    def d_hamiltonian_d_EJ2(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        EJ2 in the native Hamiltonian basis or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native Hamiltonian basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        native = -0.5 * np.kron(
             self._identity(), self._exp_i_phi_operator() + self._exp_i_phi_operator().T
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def d_hamiltonian_d_EJ3(self) -> ndarray:
-        """Returns operator representing a derivative of the Hamiltonian with
-        respect to EJ3."""
-        return (
+    def d_hamiltonian_d_EJ3(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        EJ3 in the native Hamiltonian basis or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native Hamiltonian basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        native = (
             -0.5
@@ -543,2 +630,3 @@ * (
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -552,3 +640,3 @@ def _n_operator(self) -> ndarray:
         off_diag_elements = np.ones(dim - 1, dtype=np.complex_)
-        e_iphi_matrix = np.diag(off_diag_elements, k=1)
+        e_iphi_matrix = np.diag(off_diag_elements, k=-1)
         return e_iphi_matrix
@@ -560,41 +648,199 @@
 
-    def n_1_operator(self) -> ndarray:
-        r"""Return charge number operator conjugate to :math:`\phi_1`"""
-        return np.kron(self._n_operator(), self._identity())
+    def n_1_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        r"""
+        Returns the charge number operator conjugate to :math:`\phi_1` in the charge? or eigenenergy basis.
 
-    def n_2_operator(self) -> ndarray:
-        r"""Return charge number operator conjugate to :math:`\phi_2`"""
-        return np.kron(self._identity(), self._n_operator())
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns the charge number operator conjugate to :math:`\phi_1` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns the charge number operator conjugate to :math:`\phi_1` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns the charge number operator conjugate to :math:`\phi_1` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
 
-    def exp_i_phi_1_operator(self) -> ndarray:
-        r"""Return operator :math:`e^{i\phi_1}` in the charge basis."""
-        return np.kron(self._exp_i_phi_operator(), self._identity())
+        Returns
+        -------
+            Charge number operator conjugate to :math:`\phi_1` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m, for m given eigenvectors.
+        """
+        native = np.kron(self._n_operator(), self._identity())
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def exp_i_phi_2_operator(self) -> ndarray:
-        r"""Return operator :math:`e^{i\phi_2}` in the charge basis."""
-        return np.kron(self._identity(), self._exp_i_phi_operator())
+    def n_2_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        r"""
+        Returns the charge number operator conjugate to :math:`\phi_2` in the charge? or eigenenergy basis.
 
-    def cos_phi_1_operator(self) -> ndarray:
-        """Return operator :math:`\\cos \\phi_1` in the charge basis"""
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns the charge number operator conjugate to :math:`\phi_2` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns the charge number operator conjugate to :math:`\phi_2` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns the charge number operator conjugate to :math:`\phi_2` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Charge number operator conjugate to :math:`\phi_2` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m, for m given eigenvectors.
+        """
+        native = np.kron(self._identity(), self._n_operator())
+        return self.process_op(native_op=native, energy_esys=energy_esys)
+
+    def exp_i_phi_1_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        r"""
+        Returns operator :math:`e^{i\phi_1}` in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator :math:`e^{i\phi_1}` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator :math:`e^{i\phi_1}` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator :math:`e^{i\phi_1}` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator :math:`e^{i\phi_1}` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`e^{i\phi_1}` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`e^{i\phi_1}` has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        native = np.kron(self._exp_i_phi_operator(), self._identity())
+        return self.process_op(native_op=native, energy_esys=energy_esys)
+
+    def exp_i_phi_2_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        r"""
+        Returns operator :math:`e^{i\phi_2}` in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator :math:`e^{i\phi_2}` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator :math:`e^{i\phi_2}` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator :math:`e^{i\phi_2}` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator :math:`e^{i\phi_2}` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`e^{i\phi_2}` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`e^{i\phi_2}` has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        native = np.kron(self._identity(), self._exp_i_phi_operator())
+        return self.process_op(native_op=native, energy_esys=energy_esys)
+
+    def cos_phi_1_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator :math:`\\cos \\phi_1` in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator :math:`\\cos \\phi_1` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator :math:`\\cos \\phi_1` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator :math:`\\cos \\phi_1` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator :math:`\\cos \\phi_1` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\\cos \\phi_1` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`\\cos \\phi_1` has dimensions of m x m,
+            for m given eigenvectors.
+        """
         cos_op = 0.5 * self.exp_i_phi_1_operator()
         cos_op += cos_op.T
-        return cos_op
+        native = cos_op
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def cos_phi_2_operator(self) -> ndarray:
-        """Return operator :math:`\\cos \\phi_2` in the charge basis"""
+    def cos_phi_2_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator :math:`\\cos \\phi_2` in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator :math:`\\cos \\phi_2` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator :math:`\\cos \\phi_2` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator :math:`\\cos \\phi_2` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator :math:`\\cos \\phi_2` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\\cos \\phi_2` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`\\cos \\phi_2` has dimensions of m x m,
+            for m given eigenvectors.
+        """
         cos_op = 0.5 * self.exp_i_phi_2_operator()
         cos_op += cos_op.T
-        return cos_op
+        native = cos_op
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def sin_phi_1_operator(self) -> ndarray:
-        """Return operator :math:`\\sin \\phi_1` in the charge basis"""
+    def sin_phi_1_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator :math:`\\sin \\phi_1` in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator :math:`\\sin \\phi_1` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator :math:`\\sin \\phi_1` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator :math:`\\sin \\phi_1` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator :math:`\\sin \\phi_1` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\\sin \\phi_1` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`\\sin \\phi_1` has dimensions of m x m,
+            for m given eigenvectors.
+        """
         sin_op = -1j * 0.5 * self.exp_i_phi_1_operator()
         sin_op += sin_op.conj().T
-        return sin_op
+        native = sin_op
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def sin_phi_2_operator(self) -> ndarray:
-        """Return operator :math:`\\sin \\phi_2` in the charge basis"""
+    def sin_phi_2_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator :math:`\\sin \\phi_2` in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator :math:`\\sin \\phi_2` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator :math:`\\sin \\phi_1` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator :math:`\\sin \\phi_1` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator :math:`\\sin \\phi_2` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\\sin \\phi_2` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`\\sin \\phi_2` has dimensions of m x m,
+            for m given eigenvectors.
+        """
         sin_op = -1j * 0.5 * self.exp_i_phi_2_operator()
         sin_op += sin_op.conj().T
-        return sin_op
+        native = sin_op
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -601,0 +847,0 @@ def plot_potential(

scqubits/core/fluxonium.py

@@ -142,10 +142,25 @@ # fluxonium.py
 
-    def phi_operator(self) -> ndarray:
+    def phi_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
         """
+        Returns the phi operator in the LC harmonic oscillator or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns phi operator in the LC harmonic oscillator basis.
+            If `True`, the energy eigenspectrum is computed, returns phi operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns phi operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            Returns the phi operator in the LC harmonic oscillator basis
+            Phi operator in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, phi operator has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, phi operator has dimensions of m x m,
+            for m given eigenvectors.
         """
         dimension = self.hilbertdim()
-        return (
+        native = (
             (op.creation(dimension) + op.annihilation(dimension))
@@ -156,11 +171,27 @@ * self.phi_osc()
 
-    def n_operator(self) -> ndarray:
+        return self.process_op(native_op=native, energy_esys=energy_esys)
+
+    def n_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
         """
+        Returns the :math:`n = - i d/d\\phi` operator in the LC harmonic oscillator or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns the :math:`n = - i d/d\\phi` operator in the LC harmonic oscillator basis.
+            If `True`, the energy eigenspectrum is computed, returns the :math:`n = - i d/d\\phi` operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns the :math:`n = - i d/d\\phi` operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            Returns the :math:`n = - i d/d\\phi` operator in the LC harmonic
-            oscillator basis
+            Operator :math:`n = - i d/d\\phi` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`n = - i d/d\\phi` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`n = - i d/d\\phi` has dimensions of
+            m x m, for m given eigenvectors.
         """
         dimension = self.hilbertdim()
-        return (
+        native = (
             1j
@@ -170,39 +201,113 @@ * (op.creation(dimension) - op.annihilation(dimension))
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def exp_i_phi_operator(self, alpha: float = 1.0, beta: float = 0.0) -> ndarray:
+    def exp_i_phi_operator(
+        self,
+        alpha: float = 1.0,
+        beta: float = 0.0,
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> ndarray:
         """
+        Returns the :math:`e^{i (\\alpha \\phi + \\beta) }` operator, with :math:`\\alpha` and :math:`\\beta` being
+        numbers, in the LC harmonic oscillator or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns the :math:`e^{i (\\alpha \\phi + \\beta) }` operator in the LC harmonic
+            oscillator basis. If `True`, the energy eigenspectrum is computed, returns the
+            :math:`e^{i (\\alpha \\phi + \\beta) }` operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns the :math:`e^{i (\\alpha \\phi + \\beta) }` operator in the energy eigenbasis, and does not have to
+            recalculate eigenspectrum.
+
         Returns
         -------
-            Returns the :math:`e^{i (\\alpha \\phi + \\beta) }` operator in the
-            LC harmonic oscillator basis,
-            with :math:`\\alpha` and :math:`\\beta` being numbers
+            Operator :math:`e^{i (\\alpha \\phi + \\beta) }` in chosen basis as ndarray. If the eigenenergy basis is
+            chosen, unless energy_esys is specified, :math:`e^{i (\\alpha \\phi + \\beta) }` has dimensions of
+            `truncated_dim`x `truncated_dim`. Otherwise, if eigenenergy basis is chosen,
+            :math:`e^{i (\\alpha \\phi + \\beta) }` has dimensions of m x m, for m given eigenvectors.
         """
         exponent = 1j * (alpha * self.phi_operator())
-        return sp.linalg.expm(exponent) * cmath.exp(1j * beta)
+        native = sp.linalg.expm(exponent) * cmath.exp(1j * beta)
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def cos_phi_operator(self, alpha: float = 1.0, beta: float = 0.0) -> ndarray:
+    def cos_phi_operator(
+        self,
+        alpha: float = 1.0,
+        beta: float = 0.0,
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> ndarray:
         """
+        Returns the :math:`\\cos (\\alpha \\phi + \\beta)` operator with :math:`\\alpha` and :math:`\\beta` being
+        numbers, in the LC harmonic oscillator or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns the :math:`\\cos (\\alpha \\phi + \\beta)` operator in the LC harmonic oscillator basis.
+            If `True`, the energy eigenspectrum is computed, returns the :math:`\\cos (\\alpha \\phi + \\beta)` operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns the :math:`\\cos (\\alpha \\phi + \\beta)` operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            Returns the :math:`\\cos (\\alpha \\phi + \\beta)` operator in the LC
-            harmonic oscillator basis,
-            with :math:`\\alpha` and :math:`\\beta` being numbers
+            Operator :math:`\\cos (\\alpha \\phi + \\beta)` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\\cos (\\alpha \\phi + \\beta)` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`\\cos (\\alpha \\phi + \\beta)` has dimensions of m x m, for m given eigenvectors.
         """
         argument = alpha * self.phi_operator() + beta * np.eye(self.hilbertdim())
-        return sp.linalg.cosm(argument)
+        native = sp.linalg.cosm(argument)
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def sin_phi_operator(self, alpha: float = 1.0, beta: float = 0.0) -> ndarray:
+    def sin_phi_operator(
+        self,
+        alpha: float = 1.0,
+        beta: float = 0.0,
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> ndarray:
         """
+        Returns the :math:`\\sin (\\alpha \\phi + \\beta)` operator with :math:`\\alpha` and :math:`\\beta` being
+        numbers, in the LC harmonic oscillator or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns the :math:`\\sin (\\alpha \\phi + \\beta)` operator in the LC harmonic oscillator basis.
+            If `True`, the energy eigenspectrum is computed, returns the :math:`\\sin (\\alpha \\phi + \\beta)` operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns the :math:`\\sin (\\alpha \\phi + \\beta)` operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            Returns the :math:`\\sin (\\alpha \\phi + \\beta)` operator in the
-            LC harmonic oscillator basis
-            with :math:`\\alpha` and :math:`\\beta` being numbers
+            Operator :math:`\\sin (\\alpha \\phi + \\beta)` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\\sin (\\alpha \\phi + \\beta)` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`\\sin (\\alpha \\phi + \\beta)` has dimensions of m x m, for m given eigenvectors.
         """
         argument = alpha * self.phi_operator() + beta * np.eye(self.hilbertdim())
-        return sp.linalg.sinm(argument)
+        native = sp.linalg.sinm(argument)
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def hamiltonian(self) -> ndarray:  # follow Zhu et al., PRB 87, 024510 (2013)
-        """Construct Hamiltonian matrix in harmonic-oscillator basis, following Zhu
-        et al., PRB 87, 024510 (2013)"""
+    def hamiltonian(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:  # follow Zhu et al., PRB 87, 024510 (2013)
+        """
+        Constructs Hamiltonian matrix in harmonic-oscillator, following Zhu
+        et al., PRB 87, 024510 (2013), or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns Hamiltonian in the harmonic-oscillator basis.
+            If `True`, the energy eigenspectrum is computed, returns Hamiltonian in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns Hamiltonian in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Hamiltonian in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, the Hamiltonian has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, Hamiltonian has dimensions of m x m,
+            for m given eigenvectors.
+        """
         dimension = self.hilbertdim()
@@ -215,20 +320,64 @@ diag_elements = [(i + 0.5) * self.plasma_energy() for i in range(dimension)]
         hamiltonian_mat = lc_osc_matrix - self.EJ * cos_matrix
-        return hamiltonian_mat
+        return self.process_hamiltonian(
+            native_hamiltonian=hamiltonian_mat, energy_esys=energy_esys
+        )
 
-    def d_hamiltonian_d_EJ(self) -> ndarray:
-        """Returns operator representing a derivative of the Hamiltonian with respect
-        to `EJ`.
+    def d_hamiltonian_d_EJ(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        EJ in the harmonic-oscillator or eigenenergy basis. The flux is grouped as in the Hamiltonian.
 
-        The flux is grouped as in the Hamiltonian.
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m,
+            for m given eigenvectors.
         """
-        return -self.cos_phi_operator(1, 2 * np.pi * self.flux)
+        native = -self.cos_phi_operator(1, 2 * np.pi * self.flux)
 
-    def d_hamiltonian_d_flux(self) -> ndarray:
-        """Returns operator representing a derivative of the Hamiltonian with respect
-        to `flux`.
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-        Flux is grouped as in the Hamiltonian.
+    # def d_hamiltonian_d_flux(self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False) -> ndarray:
+    #     """Returns operator representing a derivative of the Hamiltonian with respect
+    #     to flux.
+    #
+    #     The flux is grouped as in the Hamiltonian."""
+    #     return -2 * np.pi * self.EJ * self.sin_phi_operator(1, 2 * np.pi * self.flux, energy_esys=energy_esys)
+
+    def d_hamiltonian_d_flux(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
         """
-        return -2 * np.pi * self.EJ * self.sin_phi_operator(1, 2 * np.pi * self.flux)
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        flux in the harmonic-oscillator or eigenenergy basis. The flux is grouped as in the Hamiltonian.
 
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        native = -2 * np.pi * self.EJ * self.sin_phi_operator(1, 2 * np.pi * self.flux)
+        return self.process_op(native_op=native, energy_esys=energy_esys)
+
     def hilbertdim(self) -> int:
@@ -235,0 +384,0 @@ """

scqubits/core/generic_qubit.py

@@ -27,3 +27,3 @@ # generic_qubit.py
 
-# —generic qubit (two-level system)——————————————————————————————————————————————
+# -generic qubit (two-level system)----------------------------------------------
 
@@ -30,0 +30,0 @@

scqubits/core/hilbert_space.py

@@ -318,3 +318,3 @@ # hilbert_space.py
     """Class holding information about the full Hilbert space, usually composed of
-    multiple subsys_list. The class provides methods to turn subsystem operators into
+    multiple subsystems. The class provides methods to turn subsystem operators into
     operators acting on the full Hilbert space, and establishes the interface to
@@ -335,2 +335,3 @@ qutip. Returned operators are of the `qutip.Qobj` type. The class also provides
 
+    _lookup_exists = False
     osc_subsys_list = descriptors.ReadOnlyProperty(OscillatorList)
@@ -353,3 +354,3 @@ qbt_subsys_list = descriptors.ReadOnlyProperty(QubitList)
             )
-        self._subsystems: Tuple[QuantumSys, ...] = tuple(subsystem_list)
+        self._subsystems: List[QuantumSys] = subsystem_list
         self._subsys_by_id_str = {
@@ -367,3 +368,2 @@ obj._id_str: self[index] for index, obj in enumerate(self)
 
-        self._lookup: Optional[spec_lookup.SpectrumLookupAdapter] = None
         self._osc_subsys_list = [
@@ -447,8 +447,3 @@ subsys for subsys in self if isinstance(subsys, osc.Oscillator)
     @property
-    def lookup(self):
-        """[Legacy] supporting old lookup interface."""
-        return self._lookup
-
-    @property
-    def hilbertspace(self) -> "HilbertSpace":
+    def hilbertspace(self) -> HilbertSpace:
         """[Legacy] Auxiliary reference to self for compatibility with
@@ -460,2 +455,5 @@ SpectrumLookupMixin
     @property
+    @utils.DeprecationMessage(
+        "`subsys_list` is deprecated. Use `subsystem_list` instead."
+    )
     def subsys_list(self) -> List[QuantumSys]:
@@ -471,3 +469,3 @@ return list(self._subsystems)
     @classmethod
-    def deserialize(cls, io_data: "IOData") -> "HilbertSpace":
+    def deserialize(cls, io_data: "IOData") -> HilbertSpace:
         """
@@ -480,5 +478,4 @@ Take the given IOData and return an instance of the described class,
         data = alldata_dict.pop("_data", {})
-        new_hilbertspace: "HilbertSpace" = cls(**alldata_dict)
+        new_hilbertspace: HilbertSpace = cls(**alldata_dict)
         new_hilbertspace._data = data
-        new_hilbertspace._lookup = spec_lookup.SpectrumLookupAdapter(new_hilbertspace)
         return new_hilbertspace
@@ -512,3 +509,3 @@
     @classmethod
-    def create(cls) -> "HilbertSpace":
+    def create(cls) -> HilbertSpace:
         hilbertspace = cls([])
@@ -524,12 +521,12 @@ scqubits.ui.hspace_widget.create_hilbertspace_widget(hilbertspace.__init__)
             self.broadcast("HILBERTSPACE_UPDATE")
-            if self._lookup:
-                self._lookup._out_of_sync = True
+            if self.lookup_exists():
+                self._out_of_sync = True
         elif event == "INTERACTIONTERM_UPDATE" and sender in self.interaction_list:
             self.broadcast("HILBERTSPACE_UPDATE")
-            if self._lookup:
-                self._lookup._out_of_sync = True
+            if self.lookup_exists():
+                self._out_of_sync = True
         elif event == "INTERACTIONLIST_UPDATE" and sender is self:
             self.broadcast("HILBERTSPACE_UPDATE")
-            if self._lookup:
-                self._lookup._out_of_sync = True
+            if self.lookup_exists():
+                self._out_of_sync = True
 
@@ -546,3 +543,3 @@ ###################################################################################
     @property
-    def subsystem_list(self) -> Tuple[QuantumSys, ...]:
+    def subsystem_list(self) -> List[QuantumSys]:
         return self._subsystems
@@ -562,3 +559,3 @@
     def subsystem_count(self) -> int:
-        """Returns number of subsys_list composing the joint Hilbert space"""
+        """Returns number of subsystems composing the joint Hilbert space"""
         return len(self._subsystems)
@@ -570,3 +567,3 @@
     def generate_lookup(self, update_subsystem_indices: List[int] = None) -> None:
-
+        self._lookup_exists = True
         bare_esys_dict = self.generate_bare_esys(
@@ -590,5 +587,7 @@ update_subsystem_indices=update_subsystem_indices
         )
-        self._lookup = spec_lookup.SpectrumLookupAdapter(self)
 
-    def generate_bare_esys(self, update_subsystem_indices: List[int] = None) -> None:
+    def lookup_exists(self) -> bool:
+        return self._lookup_exists
+
+    def generate_bare_esys(self, update_subsystem_indices: List[int] = None) -> dict:
         # update all the subsystems when update_subsystem_indices is set to None
@@ -638,3 +637,3 @@ if update_subsystem_indices is None:
         """Calculates eigenvalues of the full Hamiltonian using
-        `qutip.Qob.eigenenergies()`.
+        `qutip.Qobj.eigenenergies()`.
 
@@ -658,3 +657,3 @@ Parameters
         """Calculates eigenvalues and eigenvectors of the full Hamiltonian using
-        `qutip.Qob.eigenstates()`.
+        `qutip.Qobj.eigenstates()`.
 
@@ -732,3 +731,3 @@ Parameters
         -------
-            composite Hamiltonian composed of bare Hamiltonians of subsys_list
+            composite Hamiltonian composed of bare Hamiltonians of subsystems
             independent of the external parameter
@@ -771,3 +770,3 @@ """
                 operator_list.append(
-                    term.hamiltonian(self.subsys_list, bare_esys=bare_esys)
+                    term.hamiltonian(self.subsystem_list, bare_esys=bare_esys)
                 )
@@ -798,3 +797,3 @@ else:
         diag_qt_op = qt.Qobj(inpt=np.diagflat(evals[0:evals_count]))  # type:ignore
-        return spec_utils.identity_wrap(diag_qt_op, subsystem, self.subsys_list)
+        return spec_utils.identity_wrap(diag_qt_op, subsystem, self.subsystem_list)
 
@@ -808,3 +807,3 @@ ###################################################################################
         the `Qobj` operator for the full Hilbert space (perform wrapping in
-        identities for other subsys_list).
+        identities for other subsystems).
 
@@ -822,3 +821,3 @@ Parameters
         diag_matrix[index, index] = diag_elements
-        return spec_utils.identity_wrap(diag_matrix, subsystem, self.subsys_list)
+        return spec_utils.identity_wrap(diag_matrix, subsystem, self.subsystem_list)
 
@@ -837,3 +836,3 @@ def hubbard_operator(self, j: int, k: int, subsystem: QuantumSys) -> Qobj:
         operator = qt.states.basis(dim, j) * qt.states.basis(dim, k).dag()
-        return spec_utils.identity_wrap(operator, subsystem, self.subsys_list)
+        return spec_utils.identity_wrap(operator, subsystem, self.subsystem_list)
 
@@ -850,3 +849,3 @@ def annihilate(self, subsystem: QuantumSys) -> Qobj:
         operator = qt.destroy(dim)
-        return spec_utils.identity_wrap(operator, subsystem, self.subsys_list)
+        return spec_utils.identity_wrap(operator, subsystem, self.subsystem_list)
 
@@ -950,2 +949,59 @@ ###################################################################################
 
+    def standardize_eigenvector_phases(self) -> None:
+        """
+        Standardize the phases of the (dressed) eigenvectors.
+        """
+        for idx, evec in enumerate(self._data["evecs"][0]):
+            phase = spec_utils.extract_phase(evec.data.toarray())
+            self._data["evecs"][0][idx] = evec * np.exp(-1j * phase)
+
+    def op_in_dressed_eigenbasis(self, **kwargs) -> Qobj:
+        """
+        Express a subsystem operator in the dressed eigenbasis of the full system
+        (as opposed to both the "native basis" or "bare eigenbasis" of the subsystem).
+        `op_in_dressed_eigenbasis(...)` offers two different interfaces:
+
+        1. subsystem operators may be expressed as Callables
+
+            signature::
+
+                .op_in_dressed_eigenbasis(op=<Callable>)
+
+        2. subsystem operators may be passed as arrays, along with the
+           corresponding subsystem. In this case the user must additionally
+           specify if the operator is in the native, subsystem-internal
+           basis or the subsystem bare eigenbasis::
+
+                .op_in_dressed_eigenbasis(op=(<ndarray>, <subsys>),
+                                          op_in_bare_eigenbasis=<Bool>)
+        """
+        op_callable_or_tuple = kwargs.pop("op")
+        if isinstance(op_callable_or_tuple, Callable):
+            subsys_index, op = self._parse_op(op_callable_or_tuple)
+            return self._op_in_dressed_eigenbasis(
+                op, subsys_index, op_in_bare_eigenbasis=False
+            )
+        else:
+            op, subsys = op_callable_or_tuple
+            op_in_bare_eigenbasis = kwargs.pop("op_in_bare_eigenbasis", False)
+            subsys_index = self.get_subsys_index(subsys)
+            return self._op_in_dressed_eigenbasis(
+                op, subsys_index, op_in_bare_eigenbasis
+            )
+
+    def _op_in_dressed_eigenbasis(
+        self, op: ndarray, subsys_index: int, op_in_bare_eigenbasis: bool = False
+    ) -> Qobj:
+        bare_evecs = self._data["bare_evecs"][subsys_index][0]
+        id_wrapped_op = spec_utils.identity_wrap(
+            op,
+            self.subsystem_list[subsys_index],
+            self.subsystem_list,
+            op_in_eigenbasis=op_in_bare_eigenbasis,
+            evecs=bare_evecs,
+        )
+        dressed_evecs = self._data["evecs"][0]
+        dressed_op = id_wrapped_op.transform(dressed_evecs)
+        return dressed_op
+
     ###################################################################################
@@ -1015,4 +1071,4 @@ # HilbertSpace: add interaction and parsing arguments to .add_interaction
             )
-        if self._lookup is not None:
-            self._lookup._out_of_sync = True
+        if self.lookup_exists():
+            self._out_of_sync = True
 
@@ -1019,0 +1075,0 @@ self.interaction_list.append(interaction)

scqubits/core/namedslots_array.py

@@ -22,2 +22,3 @@ # namedslots_array.py
 
+from matplotlib import rc_context
 from matplotlib.axes import Axes
@@ -610,2 +611,3 @@ from matplotlib.figure import Figure
 
+    @rc_context(settings.matplotlib_settings)
     def plot(self, **kwargs) -> Tuple[Figure, Axes]:
@@ -612,0 +614,0 @@ if len(self._parameters) != 1:

scqubits/core/noise.py

@@ -20,5 +20,5 @@ # noise.py
 import matplotlib.pyplot as plt
+import matplotlib as mpl
 import numpy as np
 import scipy as sp
-import scipy.constants
 
@@ -36,3 +36,5 @@ from matplotlib.axes import Axes
 from scqubits.core.storage import SpectrumData
+from scqubits.settings import matplotlib_settings
 
+
 # flag that lets us show a warning about the default t1 behavior
@@ -123,2 +125,3 @@ # (i.e., total=True setting) only once. Using the standard warnings
 
+    @mpl.rc_context(matplotlib_settings)
     def plot_coherence_vs_paramvals(
@@ -253,3 +256,2 @@ self,
         for channel_idx, noise_channel in enumerate(noise_channels):  # type:ignore
-
             # case 1: noise_channel is a string representing the noise method
@@ -338,2 +340,3 @@ if isinstance(noise_channel, str):
 
+    @mpl.rc_context(matplotlib_settings)
     def plot_t1_effective_vs_paramvals(
@@ -493,2 +496,3 @@ self,
 
+    @mpl.rc_context(matplotlib_settings)
     def plot_t2_effective_vs_paramvals(
@@ -575,3 +579,2 @@ self,
         if spectrum_data is None:
-
             # We have to figure out the largest energy level involved in the
@@ -674,6 +677,4 @@ # calculations, to know how many levels we need from the diagonalization.
         for n, noise_channel in enumerate(noise_channels):
-
             # noise_channel is a string representing the noise method
             if isinstance(noise_channel, str):
-
                 noise_channel_method = noise_channel
@@ -700,3 +701,2 @@
             elif isinstance(noise_channel, tuple):
-
                 noise_channel_method = noise_channel[0]
@@ -703,0 +703,0 @@

scqubits/core/oscillator.py

@@ -73,3 +73,3 @@ # oscillator.py
 
-# —Oscillator class———————————————————————————————————————————————————————————————————
+# -Oscillator class-------------------------------------------------------------------
 
@@ -190,3 +190,3 @@
 
-# —KerrOscillator class———————————————————————————————————————————————————————————————————
+# -KerrOscillator class-------------------------------------------------------------------
 
@@ -223,3 +223,2 @@
     ) -> None:
-
         self.K: float = K  # type:ignore
@@ -226,0 +225,0 @@

scqubits/core/param_sweep.py

@@ -118,2 +118,3 @@ # param_sweep.py
 
+    _lookup_exists = False
     _parameters = descriptors.WatchedProperty(Parameters, "PARAMETERSWEEP_UPDATE")
@@ -608,3 +609,3 @@ _evals_count = descriptors.WatchedProperty(int, "PARAMETERSWEEP_UPDATE")
 
-        self._preslicing_reset()
+        self.reset_preslicing()
 
@@ -733,3 +734,3 @@ if not as_specdata:
 
-        self._preslicing_reset()
+        self.reset_preslicing()
 
@@ -892,3 +893,3 @@ if coloring == "plain":
         if set to False (default), bare product states and dressed eigenstates are
-        identified if `\|<psi_bare\|psi_dressed>\|^2 > 0.5`; if True,
+        identified if `|<psi_bare|psi_dressed>|^2 > 0.5`; if True,
         then identification will always take place based on which bare product state
@@ -962,3 +963,3 @@ has the maximum overlap
         self._out_of_sync = False
-        self._preslicing_reset()
+        self.reset_preslicing()
 
@@ -1002,2 +1003,3 @@ dispatch.CENTRAL_DISPATCH.register("PARAMETERSWEEP_UPDATE", self)
 
+        self._lookup_exists = True
         if self._deepcopy:
@@ -1334,2 +1336,3 @@ stored_hilbertspace = copy.deepcopy(self.hilbertspace)
     ) -> None:
+        self._lookup_exists = True
         self._parameters = Parameters(paramvals_by_name)
@@ -1336,0 +1339,0 @@ self._hilbertspace = hilbertspace

scqubits/core/qubit_base.py

@@ -33,2 +33,3 @@ # qubit_base.py
 import matplotlib.pyplot as plt
+import matplotlib as mpl
 import numpy as np
@@ -40,2 +41,3 @@ import scipy as sp
 from numpy import ndarray
+from scipy.sparse import csc_matrix, dia_matrix
 
@@ -52,3 +54,3 @@ import scqubits.core.constants as constants
 from scqubits.core.storage import DataStore, SpectrumData
-from scqubits.settings import IN_IPYTHON
+from scqubits.settings import IN_IPYTHON, matplotlib_settings
 from scqubits.utils.cpu_switch import get_map_method
@@ -78,3 +80,3 @@ from scqubits.utils.misc import InfoBar, process_which
 
-# —Generic quantum system container and Qubit base class——————————————————————————————
+# -Generic quantum system container and Qubit base class------------------------------
 
@@ -236,3 +238,3 @@
 
-# —QubitBaseClass———————————————————————————————————————————————————————————————————————————————————————————————————————
+# -QubitBaseClass-------------------------------------------------------------------------------------------------------
 
@@ -258,3 +260,3 @@
         evals = sp.linalg.eigh(
-            hamiltonian_mat, eigvals_only=True, eigvals=(0, evals_count - 1)
+            hamiltonian_mat, eigvals_only=True, subset_by_index=(0, evals_count - 1)
         )
@@ -266,3 +268,3 @@ return np.sort(evals)
         evals, evecs = sp.linalg.eigh(
-            hamiltonian_mat, eigvals_only=False, eigvals=(0, evals_count - 1)
+            hamiltonian_mat, eigvals_only=False, subset_by_index=(0, evals_count - 1)
         )
@@ -375,2 +377,79 @@ evals, evecs = order_eigensystem(evals, evecs)
 
+    def process_op(
+        self,
+        native_op: Union[ndarray, csc_matrix],
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> Union[ndarray, csc_matrix]:
+        """Processes the operator `native_op`: either hand back `native_op` unchanged, or transform it into the
+        energy eigenbasis. (Native basis refers to the basis used internally by each qubit, e.g., charge basis in the
+        case of `Transmon`.
+
+        Parameters
+        ----------
+        native_op:
+            operator in native basis
+        energy_esys:
+            If `False` (default), returns operator in the native basis
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis
+            if energy_esys is the energy eigenspectrum, in the form of a tuple containing two ndarrays
+            (eigenvalues and energy eigenvectors), returns operator in the energy eigenbasis,
+            and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            `native_op` either unchanged or transformed into eigenenergy basis
+        """
+        if isinstance(energy_esys, bool):
+            if not energy_esys:
+                return native_op
+            esys = self.eigensys(evals_count=self.truncated_dim)
+        else:
+            esys = energy_esys
+        evectors = esys[1][:, : self.truncated_dim]
+        return get_matrixelement_table(native_op, evectors)
+
+    def process_hamiltonian(
+        self,
+        native_hamiltonian: Union[ndarray, csc_matrix],
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> Union[ndarray, csc_matrix]:
+        """Return qubit Hamiltonian in chosen basis: either return unchanged (i.e., in native basis) or transform
+        into eigenenergy basis
+
+        Parameters
+        ----------
+        native_hamiltonian:
+            Hamiltonian in native basis
+        energy_esys:
+            If `False` (default), returns Hamiltonian in the native basis
+            If `True`, the energy eigenspectrum is computed, returns Hamiltonian in the energy eigenbasis
+            if energy_esys is the energy eigenspectrum, in the form of a tuple containing two ndarrays
+            (eigenvalues and energy eigenvectors), returns Hamiltonian in the energy eigenbasis,
+            and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Hamiltonian, either unchanged in native basis, or transformed into eigenenergy basis
+        """
+        if isinstance(energy_esys, bool):
+            if not energy_esys:
+                return native_hamiltonian
+            esys = self.eigensys(evals_count=self.truncated_dim)
+        else:
+            esys = energy_esys
+        evals = esys[0][: self.truncated_dim]
+        if isinstance(native_hamiltonian, ndarray):
+            return np.diag(evals)
+        return dia_matrix(evals).tocsc()
+
+    def anharmonicity(self) -> float:
+        """Returns the qubit's anharmonicity, (E_2 - E_1) - (E_1 - E_0)."""
+        energies = self.eigenvals(evals_count=3)
+        return energies[2] - 2 * energies[1] + energies[0]
+
+    def E01(self) -> float:
+        """Returns the qubit's fundamental energy splitting, E_1 - E_0."""
+        energies = self.eigenvals(evals_count=2)
+        return energies[1] - energies[0]
+
     @overload
@@ -749,8 +828,11 @@ def matrixelement_table(
 
+        paramval_before = getattr(self, param_name)
         assert spectrumdata.state_table is not None
         for index, paramval in enumerate(param_vals):
             evecs = spectrumdata.state_table[index]
+            setattr(self, param_name, paramval)
             matelem_table[index] = self.matrixelement_table(
                 operator, evecs=evecs, evals_count=evals_count
             )
+        setattr(self, param_name, paramval_before)
 
@@ -760,2 +842,3 @@ spectrumdata.matrixelem_table = matelem_table
 
+    @mpl.rc_context(matplotlib_settings)
     def plot_evals_vs_paramvals(
@@ -802,2 +885,3 @@ self,
 
+    @mpl.rc_context(matplotlib_settings)
     def plot_dispersion_vs_paramvals(
@@ -878,2 +962,3 @@ self,
 
+    @mpl.rc_context(matplotlib_settings)
     def plot_matrixelements(
@@ -932,2 +1017,3 @@ self,
 
+    @mpl.rc_context(matplotlib_settings)
     def plot_matelem_vs_paramvals(
@@ -1009,3 +1095,3 @@ self,
 
-# —QubitBaseClass1d——————————————————————————————————————————————————————————————————
+# -QubitBaseClass1d------------------------------------------------------------------
 
@@ -1055,2 +1141,3 @@
 
+    @mpl.rc_context(matplotlib_settings)
     def plot_wavefunction(
@@ -1057,0 +1144,0 @@ self,

scqubits/core/spec_lookup.py

@@ -14,2 +14,3 @@ # spec_lookup.py
 import itertools
+import numbers
 
@@ -38,4 +39,2 @@ from typing import TYPE_CHECKING, Any, Dict, List, Optional, Tuple, Union
     from scqubits.core.param_sweep import Parameters
-    from scqubits.core.qubit_base import QuantumSystem
-    from scqubits.io_utils.fileio_qutip import QutipEigenstates
     from scqubits.utils.typedefs import QuantumSys
@@ -61,188 +60,2 @@
 
-class SpectrumLookupAdapter:
-    """Bridges earlier functionality provided by the `SpectrumLookup` class to
-    universal use of the `SpectrumLookupMixin` class (originally developed for
-    `ParameterSweep` and then extended to `HilbertSpace`)."""
-
-    def __init__(self, hilbertspace):
-        self.hilbertspace = hilbertspace
-
-    @property
-    def _out_of_sync(self):
-        return self.hilbertspace._out_of_sync
-
-    @_out_of_sync.setter
-    def _out_of_sync(self, value: bool):
-        self.hilbertspace._out_of_sync = value
-
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.run is deprecated. "
-        "Use <HilbertSpace>.generate_lookup instead."
-    )
-    def run(self) -> None:
-        self.hilbertspace.generate_lookup()
-
-    @utils.check_sync_status
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.dressed_index is "
-        "deprecated. Use <HilbertSpace>.dressed_index instead."
-    )
-    def dressed_index(self, bare_labels: Tuple[int, ...]) -> Union[int, None]:
-        """
-        For given bare product state return the corresponding dressed-state index.
-
-        Parameters
-        ----------
-        bare_labels:
-            bare_labels = (index, index2, ...)
-
-        Returns
-        -------
-            dressed state index closest to the specified bare state
-        """
-        return self.hilbertspace.dressed_index(bare_labels)
-
-    @utils.check_sync_status
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.bare_index is "
-        "deprecated. Use <HilbertSpace>.bare_index instead."
-    )
-    def bare_index(self, dressed_index: int) -> Union[Tuple[int, ...], None]:
-        """
-        For given dressed index, look up the corresponding bare index.
-
-        Returns
-        -------
-            Bare state specification in tuple form. Example: (1,0,3) means subsystem 1
-            is in bare state 1, subsystem 2 in bare state 0, and subsystem 3 in bare
-            state 3.
-        """
-        return self.hilbertspace.bare_index(dressed_index)
-
-    @utils.check_sync_status
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.dressed_eigenstates is "
-        "deprecated. Use <HilbertSpace>['evecs'] instead."
-    )
-    def dressed_eigenstates(self) -> List["QutipEigenstates"]:
-        """
-        Return the list of dressed eigenvectors
-
-        Returns
-        -------
-            dressed eigenvectors for the external parameter fixed to the value
-            indicated by the provided index
-        """
-        return self.hilbertspace["evecs"]
-
-    @utils.check_sync_status
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.dressed_eigenenergies is "
-        "deprecated. Use <HilbertSpace>['evals'] instead."
-    )
-    def dressed_eigenenergies(self) -> ndarray:
-        """
-        Return the array of dressed eigenenergies
-
-
-        Returns
-        -------
-            dressed eigenenergies for the external parameter fixed to the value
-            indicated by the provided index
-        """
-        return self.hilbertspace["evals"]
-
-    @utils.check_sync_status
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.energy_bare_index is "
-        "deprecated. Use "
-        "<HilbertSpace>.energy_by_bare_index instead."
-    )
-    def energy_bare_index(self, bare_tuple: Tuple[int, ...]) -> Union[float, None]:
-        """
-        Look up dressed energy most closely corresponding to the given bare-state labels
-
-        Parameters
-        ----------
-        bare_tuple:
-            bare state indices
-
-        Returns
-        -------
-            dressed energy, if lookup successful
-        """
-        return self.hilbertspace.energy_by_bare_index(bare_tuple)
-
-    @utils.check_sync_status
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.energy_dressed_index is "
-        "deprecated. Use "
-        "<HilbertSpace>.energy_by_dressed_index instead."
-    )
-    def energy_dressed_index(self, dressed_index: int) -> float:
-        """
-        Look up the dressed eigenenergy belonging to the given dressed index.
-
-        Parameters
-        ----------
-        dressed_index:
-            index of dressed state of interest
-
-        Returns
-        -------
-            dressed energy
-        """
-        return self.hilbertspace.energy_by_dressed_index(dressed_index)
-
-    @utils.check_sync_status
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.bare_eigenstates is deprecated. "
-        "Use <HilbertSpace>.bare_eigenstates instead."
-    )
-    def bare_eigenstates(self, subsys: "QuantumSystem") -> ndarray:
-        """
-        Return ndarray of bare eigenstates for given subsystem.
-        Eigenstates are expressed in the basis internal to the subsystem.
-        """
-        return self.hilbertspace.bare_eigenstates(subsys)
-
-    @utils.check_sync_status
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.bare_eigenenergies is deprecated. "
-        "Use <HilbertSpace>.bare_eigenvals instead."
-    )
-    def bare_eigenenergies(self, subsys: "QuantumSystem") -> ndarray:
-        """
-        Return list of bare eigenenergies for given subsystem.
-
-        Parameters
-        ----------
-        subsys:
-            Hilbert space subsystem for which bare eigendata is to be looked up
-
-        Returns
-        -------
-            bare eigenenergies for the specified subsystem
-        """
-        return self.hilbertspace.bare_eigenvals(subsys)
-
-    @utils.DeprecationMessage(
-        "<HilbertSpace>.lookup.bare_productstate is deprecated. "
-        "Use <HilbertSpace>.bare_productstate instead."
-    )
-    def bare_productstate(self, bare_index: Tuple[int, ...]) -> Qobj:
-        """
-        Return the bare product state specified by `bare_index`.
-
-        Parameters
-        ----------
-        bare_index:
-
-        Returns
-        -------
-            ket in full Hilbert space
-        """
-        return self.hilbertspace.bare_productstate(bare_index)
-
-
 class SpectrumLookupMixin(MixinCompatible):
@@ -255,2 +68,17 @@ """
 
+    _inside_hilbertspace = False
+
+    def __init_subclass__(cls):
+        super().__init_subclass__()
+        if cls.__name__ == "HilbertSpace":
+            cls._inside_hilbertspace = True
+        else:
+            cls._inside_hilbertspace = False
+
+    def reset_preslicing(self):
+        if self._inside_hilbertspace:
+            self._current_param_indices = 0
+        else:
+            self._current_param_indices = slice(None, None, None)
+
     @property
@@ -336,2 +164,3 @@ def _bare_product_states_labels(self) -> List[Tuple[int, ...]]:
 
+    @utils.check_lookup_exists
     @utils.check_sync_status
@@ -364,2 +193,3 @@ def dressed_index(
 
+    @utils.check_lookup_exists
     @utils.check_sync_status
@@ -439,2 +269,3 @@ def bare_index(
 
+    @utils.check_lookup_exists
     @utils.check_sync_status
@@ -446,3 +277,3 @@ def energy_by_bare_index(
         param_indices: Optional[NpIndices] = None,
-    ) -> NamedSlotsNdarray:  # the return value may also be np.nan
+    ) -> Union[float, NamedSlotsNdarray]:  # the return value may also be np.nan
         """
@@ -469,3 +300,3 @@ Look up dressed energy most closely corresponding to the given bare-state labels
             return np.nan  # type:ignore
-        if isinstance(dressed_index, int):
+        if isinstance(dressed_index, numbers.Number):
             energy = self["evals"][param_indices + (dressed_index,)]
@@ -477,3 +308,3 @@ if subtract_ground:
         dressed_index = np.asarray(dressed_index)
-        energies = np.empty_like(dressed_index)
+        energies = np.empty_like(dressed_index, dtype=np.float_)
         it = np.nditer(dressed_index, flags=["multi_index", "refs_ok"])
@@ -494,2 +325,3 @@ sliced_energies = self["evals"][param_indices]
 
+    @utils.check_lookup_exists
     @utils.check_sync_status
@@ -501,3 +333,3 @@ def energy_by_dressed_index(
         param_indices: Optional[Tuple[int, ...]] = None,
-    ) -> float:
+    ) -> Union[float, NamedSlotsNdarray]:
         """
@@ -526,2 +358,3 @@ Look up the dressed eigenenergy belonging to the given dressed index,
 
+    @utils.check_lookup_exists
     @utils.check_sync_status
@@ -536,9 +369,10 @@ def bare_eigenstates(
         Eigenstates are expressed in the basis internal to the subsystems. Usually to be
-        with pre-slicing.
+        used with pre-slicing when part of `ParameterSweep`.
         """
         param_indices_tuple = self.set_npindextuple(param_indices)
         subsys_index = self.hilbertspace.get_subsys_index(subsys)
-        self._current_param_indices = slice(None, None, None)
+        self.reset_preslicing()
         return self["bare_evecs"][subsys_index][param_indices_tuple]
 
+    @utils.check_lookup_exists
     @utils.check_sync_status
@@ -568,3 +402,3 @@ def bare_eigenvals(
         subsys_index = self.hilbertspace.get_subsys_index(subsys)
-        self._current_param_indices = slice(None, None, None)
+        self.reset_preslicing()
         return self["bare_evals"][subsys_index][param_indices_tuple]
@@ -571,0 +405,0 @@

scqubits/core/storage.py

@@ -29,3 +29,3 @@ # storage.py
 
-# —WaveFunction class———————————————————————————————————————————————————————————————————
+# -WaveFunction class-------------------------------------------------------------------
 
@@ -78,3 +78,3 @@
 
-# —WaveFunctionOnGrid class—————————————————————————————————————————————————————————————
+# -WaveFunctionOnGrid class-------------------------------------------------------------
 
@@ -104,3 +104,3 @@
 
-# —BaseData class———————————————————————————————————————————————————————————————————————
+# -BaseData class-----------------------------------------------------------------------
 
@@ -166,3 +166,3 @@
 
-# —SpectrumData class———————————————————————————————————————————————————————————————————
+# -SpectrumData class-------------------------------------------------------------------
 
@@ -169,0 +169,0 @@

scqubits/core/transmon.py

@@ -42,3 +42,3 @@ # transmon.py
 
-# —Cooper pair box / transmon——————————————————————————————————————————————
+# -Cooper pair box / transmon----------------------------------------------
 
@@ -153,28 +153,162 @@
 
-    def n_operator(self) -> ndarray:
-        """Returns charge operator `n` in the charge basis"""
+    @staticmethod
+    def find_EJ_EC(
+        E01: float, anharmonicity: float, ng=0, ncut=30
+    ) -> Tuple[float, float]:
+        """
+        Finds the EJ and EC values given a qubit splitting `E01` and `anharmonicity`.
+
+        Parameters
+        ----------
+            E01:
+                qubit transition energy
+            anharmonicity:
+                absolute qubit anharmonicity, (E2-E1) - (E1-E0)
+            ng:
+                offset charge (default: 0)
+            ncut:
+                charge number cutoff (default: 30)
+
+        Returns
+        -------
+            A tuple of the EJ and EC values representing the best fit.
+        """
+        tmon = Transmon(EJ=10.0, EC=0.1, ng=ng, ncut=ncut)
+        start_EJ_EC = np.array([tmon.EJ, tmon.EC])
+
+        def cost_func(EJ_EC: Tuple[float, float]) -> float:
+            EJ, EC = EJ_EC
+            tmon.EJ = EJ
+            tmon.EC = EC
+            energies = tmon.eigenvals(evals_count=3)
+            computed_E01 = energies[1] - energies[0]
+            computed_anharmonicity = energies[2] - energies[1] - computed_E01
+            cost = (E01 - computed_E01) ** 2
+            cost += (anharmonicity - computed_anharmonicity) ** 2
+            return cost
+
+        return sp.optimize.minimize(cost_func, start_EJ_EC).x
+
+    def n_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns charge operator n in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns charge operator n in the charge basis.
+            If `True`, energy eigenspectrum is computed, returns charge operator n in the energy eigenbasis.
+            If `energy_esys = esys`, where `esys` is a tuple containing two ndarrays (eigenvalues and energy
+            eigenvectors), returns charge operator n in the energy eigenbasis, and does not have to recalculate the
+            eigenspectrum.
+
+        Returns
+        -------
+            Charge operator n in chosen basis as ndarray.
+            For `energy_esys=True`, n has dimensions of `truncated_dim` x `truncated_dim`.
+            If an actual eigensystem is handed to `energy_sys`, then `n` has dimensions of m x m,
+            where m is the number of given eigenvectors.
+        """
         diag_elements = np.arange(-self.ncut, self.ncut + 1, 1)
-        return np.diag(diag_elements)
+        native = np.diag(diag_elements)
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def exp_i_phi_operator(self) -> ndarray:
-        """Returns operator :math:`e^{i\\varphi}` in the charge basis"""
+    def exp_i_phi_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator :math:`e^{i\\varphi}` in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator :math:`e^{i\\varphi}` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator :math:`e^{i\\varphi}` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator :math:`e^{i\\varphi}` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator :math:`e^{i\\varphi}` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`e^{i\\varphi}` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`e^{i\\varphi}` has dimensions of m x m,
+            for m given eigenvectors.
+        """
         dimension = self.hilbertdim()
         entries = np.repeat(1.0, dimension - 1)
         exp_op = np.diag(entries, -1)
-        return exp_op
+        return self.process_op(native_op=exp_op, energy_esys=energy_esys)
 
-    def cos_phi_operator(self) -> ndarray:
-        """Returns operator :math:`\\cos \\varphi` in the charge basis"""
+    def cos_phi_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator :math:`\\cos \\varphi` in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator :math:`\\cos \\varphi` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator :math:`\\cos \\varphi` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator :math:`\\cos \\varphi` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator :math:`\\cos \\varphi` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\\cos \\varphi` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`\\cos \\varphi` has dimensions of m x m,
+            for m given eigenvectors.
+        """
         cos_op = 0.5 * self.exp_i_phi_operator()
         cos_op += cos_op.T
-        return cos_op
+        return self.process_op(native_op=cos_op, energy_esys=energy_esys)
 
-    def sin_phi_operator(self) -> ndarray:
-        """Returns operator :math:`\\sin \\varphi` in the charge basis"""
+    def sin_phi_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator :math:`\\sin \\varphi` in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator :math:`\\sin \\varphi` in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator :math:`\\sin \\varphi` in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator :math:`\\sin \\varphi` in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator :math:`\\sin \\varphi` in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless energy_esys is specified, :math:`\\sin \\varphi` has dimensions of truncated_dim
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, :math:`\\sin \\varphi` has dimensions of m x m,
+            for m given eigenvectors.
+        """
         sin_op = -1j * 0.5 * self.exp_i_phi_operator()
         sin_op += sin_op.conjugate().T
-        return sin_op
+        return self.process_op(native_op=sin_op, energy_esys=energy_esys)
 
-    def hamiltonian(self) -> ndarray:
-        """Returns Hamiltonian in charge basis"""
+    def hamiltonian(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns Hamiltonian in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns Hamiltonian in the charge basis.
+            If `True`, the energy eigenspectrum is computed; returns Hamiltonian in the energy eigenbasis.
+            If `energy_esys = esys`, where `esys` is a tuple containing two ndarrays (eigenvalues and energy
+            eigenvectors); then return the Hamiltonian in the energy eigenbasis, do not recalculate eigenspectrum.
+
+        Returns
+        -------
+            Hamiltonian in chosen basis as ndarray. For `energy_esys=False`, the Hamiltonian has dimensions of
+            `truncated_dim` x `truncated_dim`. For `energy_sys=esys`, the Hamiltonian has dimensions of m x m,
+            for m given eigenvectors.
+        """
         dimension = self.hilbertdim()
@@ -190,14 +324,56 @@ hamiltonian_mat = np.diag(
         hamiltonian_mat[ind + 1, ind] = -self.EJ / 2.0
-        return hamiltonian_mat
+        return self.process_hamiltonian(
+            native_hamiltonian=hamiltonian_mat, energy_esys=energy_esys
+        )
 
-    def d_hamiltonian_d_ng(self) -> ndarray:
-        """Returns operator representing a derivative of the Hamiltonian with respect to
-        charge offset `ng`."""
-        return -8 * self.EC * self.n_operator()
+    def d_hamiltonian_d_ng(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        charge offset `ng` in the charge or eigenenergy basis.
 
-    def d_hamiltonian_d_EJ(self) -> ndarray:
-        """Returns operator representing a derivative of the Hamiltonian with respect
-        to EJ."""
-        return -self.cos_phi_operator()
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where `esys` is a tuple containing two ndarrays (eigenvalues and energy
+            eigenvectors), returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
 
+        Returns
+        -------
+            Operator in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        native = -8 * self.EC * self.n_operator(energy_esys=energy_esys)
+        return self.process_op(native_op=native, energy_esys=energy_esys)
+
+    def d_hamiltonian_d_EJ(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        EJ in the charge or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where `esys` is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        native = -self.cos_phi_operator()
+        return self.process_op(native_op=native, energy_esys=energy_esys)
+
     def hilbertdim(self) -> int:
@@ -409,3 +585,3 @@ """Returns Hilbert space dimension"""
 
-# — Flux-tunable Cooper pair box / transmon———————————————————————————————————————————
+# - Flux-tunable Cooper pair box / transmon-------------------------------------------
 
@@ -509,6 +685,31 @@
 
-    def d_hamiltonian_d_flux(self) -> ndarray:
-        """Returns operator representing a derivative of the Hamiltonian with respect
-        to `flux`."""
-        return (
+    def d_hamiltonian_d_flux(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> ndarray:
+        """
+        Returns operator representing a derivative of the Hamiltonian with respect to
+        `flux` in the charge or eigenenergy basis.
+
+        Here, the derivative is taken with respect to flux before the qubit's :math:`\phi` degree of 
+        freedom in the Hamiltonian is shifted by a flux-dependent quantity :math:`\varphi_{0}` 
+        (see Eq. 2.17 and surrounding text in PRA 76, 042319 (2007)). Then only after the flux 
+        derivative is taken, both the Hamiltonian as well as its flux derivative are assumed to 
+        be shifted by :math:`\varphi_{0}`.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the charge basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis as ndarray. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, operator has dimensions of m x m,
+            for m given eigenvectors.
+        """
+        native = (
             np.pi
@@ -524,3 +725,10 @@ * self.EJmax
             * self.cos_phi_operator()
+            - np.pi * self.EJmax * self.d 
+            / np.sqrt(
+                np.cos(np.pi * self.flux) ** 2
+                + self.d**2 * np.sin(np.pi * self.flux) ** 2
+            )
+            * self.sin_phi_operator()
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -527,0 +735,0 @@ def _compute_dispersion(

scqubits/core/zeropi_full.py

@@ -297,6 +297,8 @@ # zeropi_full.py
     def hamiltonian(
-        self, return_parts: bool = False
-    ) -> Union[csc_matrix, Tuple[csc_matrix, ndarray, ndarray, float]]:
-        """Returns Hamiltonian in basis obtained by discretizing phi, employing charge basis for theta, and Fock
-        basis for zeta.
+        self,
+        return_parts: bool = False,
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> Union[csc_matrix, Tuple[csc_matrix, ndarray, ndarray, ndarray]]:
+        r"""Returns Hamiltonian in basis obtained by discretizing :math:`\phi`, employing
+        charge basis for :math:`\theta`, and Fock basis for :math:`\zeta`, or in the eigenenergy basis.
 
@@ -306,3 +308,17 @@ Parameters
         return_parts:
-            If set to true, `hamiltonian` returns [hamiltonian, evals, evecs, g_coupling_matrix]
+            If set to true, `hamiltonian` returns
+            [hamiltonian, evals, evecs, g_coupling_matrix]
+        energy_esys:
+            If `False` (default), returns Hamiltonian in native Hamiltonian basis
+            If `True`, the energy eigenspectrum is computed, returns Hamiltonian in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns Hamiltonian in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Hamiltonian in chosen basis as csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, the Hamiltonian has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, Hamiltonian has dimensions of m x m,
+            for m given eigenvectors.
+
         """
@@ -338,44 +354,122 @@ zeropi_dim = self.zeropi_cutoff
         ) + sparse.kron(zeropi_coupling.conjugate().T, op.creation_sparse(zeta_dim))
-
+        hmtocsc = hamiltonian_mat.tocsc()
         if return_parts:
-            return hamiltonian_mat.tocsc(), zeropi_evals, zeropi_evecs, gmat
+            return (
+                self.process_hamiltonian(
+                    native_hamiltonian=hmtocsc, energy_esys=energy_esys
+                ),
+                zeropi_evals,
+                zeropi_evecs,
+                gmat,
+            )
+        return self.process_hamiltonian(
+            native_hamiltonian=hmtocsc, energy_esys=energy_esys
+        )
 
-        return hamiltonian_mat.tocsc()
+    def d_hamiltonian_d_flux(
+        self,
+        zeropi_evecs: ndarray = None,
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Calculates a derivative of the Hamiltonian w.r.t flux, at the current value of flux,
+        as stored in the object. The returned operator is in the product basis or eigenenergy basis.
 
-    def d_hamiltonian_d_flux(self, zeropi_evecs: ndarray = None) -> csc_matrix:
-        r"""Calculates a derivative of the Hamiltonian w.r.t flux, at the current value of flux,
-        as stored in the object. The returned operator is in the product basis
+        Helper method _zeropi_operator_in_product_basis is employed which converts
+        a zeropi operator into one in the product basis. If user already has zeropi eigenvectors, user can input
+        as zeropi_evecs=ndarray.
 
-        The flux is assumed to be given in the units of the ratio \Phi_{ext}/\Phi_0.
-        So if \frac{\partial H}{ \partial \Phi_{\rm ext}}, is needed, the expr returned
-        by this function, needs to be multiplied by 1/\Phi_0.
+        Parameters
+        ----------
+        zeropi_evecs:
+            If None (default), helper method _zeropi_operator_in_product_basis calculates zeropi eigenvectors and uses
+            them to convert operator to the product basis, if product basis is chosen.
+            If zeropi_evecs = zeropievecs, where zeropievecs is an ndarray, and product basis is chosen,
+            helper method _zeropi_operator_in_product_basis uses zeropievecs to convert the operator
+            to the product basis.
+        energy_esys:
+            If `False` (default), returns operator in the product basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
 
         Returns
         -------
-            matrix representing the derivative of the Hamiltonian
+            Operator in chosen basis. If product basis is chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return self._zeropi_operator_in_product_basis(
+        native = self._zeropi_operator_in_product_basis(
             self._zeropi.d_hamiltonian_d_flux(), zeropi_evecs=zeropi_evecs
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def d_hamiltonian_d_EJ(self, zeropi_evecs: ndarray = None) -> csc_matrix:
-        r"""Calculates a derivative of the Hamiltonian w.r.t EJ.
+    def d_hamiltonian_d_EJ(
+        self,
+        zeropi_evecs: ndarray = None,
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Calculates a derivative of the Hamiltonian w.r.t EJ. The returned operator is in the
+        product basis or eigenenergy basis.
 
+        Helper method _zeropi_operator_in_product_basis is employed which converts
+        a zeropi operator into one in the product basis. If user already has zeropi eigenvectors, user can input
+        as zeropi_evecs=ndarray.
+
+        Parameters
+        ----------
+        zeropi_evecs:
+            If None (default), helper method _zeropi_operator_in_product_basis calculates zeropi eigenvectors and uses
+            them to convert operator to the product basis, if product basis is chosen.
+            If zeropi_evecs = zeropievecs, where zeropievecs is an ndarray, and product basis is chosen,
+            helper method _zeropi_operator_in_product_basis uses zeropievecs to convert the operator
+            to the product basis.
+        energy_esys:
+            If `False` (default), returns operator in the product basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            matrix representing the derivative of the Hamiltonian
+            Operator in chosen basis. If product basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return self._zeropi_operator_in_product_basis(
+        native = self._zeropi_operator_in_product_basis(
             self._zeropi.d_hamiltonian_d_EJ(), zeropi_evecs=zeropi_evecs
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def d_hamiltonian_d_ng(self) -> csc_matrix:
-        r"""Calculates a derivative of the Hamiltonian w.r.t ng.
-        as stored in the object.
+    def d_hamiltonian_d_ng(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Calculates a derivative of the Hamiltonian w.r.t ng.
+        Returns matrix representing a derivative of the Hamiltonian in the native Hamiltonian basis
+        or eigenenergy basis.
 
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            matrix representing the derivative of the Hamiltonian
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return -8 * self.EC * self.n_theta_operator()
+        native = -8 * self.EC * self.n_theta_operator()
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -412,29 +506,118 @@ def _zeropi_operator_in_product_basis(
 
-    def i_d_dphi_operator(self, zeropi_evecs: ndarray = None) -> csc_matrix:
+    def i_d_dphi_operator(
+        self,
+        zeropi_evecs: ndarray = None,
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> Union[ndarray, csc_matrix]:
         r"""
-        Operator :math:`i d/d\phi`.
+        Returns :math:`i d/d\phi` operator in the product or eigenenergy basis.
+
+        Helper method _zeropi_operator_in_product_basis is employed which converts
+        a zeropi operator into one in the product basis. If user already has zeropi eigenvectors, user can input
+        as zeropi_evecs=ndarray.
+
+        Parameters
+        ----------
+        zeropi_evecs:
+            If None (default), helper method _zeropi_operator_in_product_basis calculates zeropi eigenvectors and uses
+            them to convert operator to the product basis, if product basis is chosen.
+            If zeropi_evecs = zeropievecs, where zeropievecs is an ndarray, and product basis is chosen,
+            helper method _zeropi_operator_in_product_basis uses zeropievecs to convert the operator
+            to the product basis.
+        energy_esys:
+            If `False` (default), returns operator in the product basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If product basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return self._zeropi_operator_in_product_basis(
+        native = self._zeropi_operator_in_product_basis(
             self._zeropi.i_d_dphi_operator(), zeropi_evecs=zeropi_evecs
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def n_theta_operator(self, zeropi_evecs: ndarray = None) -> csc_matrix:
+    def n_theta_operator(
+        self,
+        zeropi_evecs: ndarray = None,
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> Union[ndarray, csc_matrix]:
         r"""
-        Operator :math:`n_\theta`.
+        Returns :math:`n_\theta`  operator in the product or eigenenergy basis.
+
+        Helper method _zeropi_operator_in_product_basis is employed which converts
+        a zeropi operator into one in the product basis. If user already has zeropi eigenvectors, user can input
+        as zeropi_evecs=ndarray.
+
+        Parameters
+        ----------
+        zeropi_evecs:
+            If None (default), helper method _zeropi_operator_in_product_basis calculates zeropi eigenvectors and uses
+            them to convert operator to the product basis, if product basis is chosen.
+            If zeropi_evecs = zeropievecs, where zeropievecs is an ndarray, and product basis is chosen,
+            helper method _zeropi_operator_in_product_basis uses zeropievecs to convert the operator
+            to the product basis.
+        energy_esys:
+            If `False` (default), returns operator in the product basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If product basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return self._zeropi_operator_in_product_basis(
+        native = self._zeropi_operator_in_product_basis(
             self._zeropi.n_theta_operator(), zeropi_evecs=zeropi_evecs
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def phi_operator(self, zeropi_evecs: ndarray = None) -> csc_matrix:
+    def phi_operator(
+        self,
+        zeropi_evecs: ndarray = None,
+        energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False,
+    ) -> Union[ndarray, csc_matrix]:
         r"""
-        Operator :math:`\phi`.
+        Returns :math:`\phi`  operator in the product or eigenenergy basis.
 
+        Helper method _zeropi_operator_in_product_basis is employed which converts
+        a zeropi operator into one in the product basis. If user already has zeropi eigenvectors, user can input
+        as zeropi_evecs=ndarray.
+
+        Parameters
+        ----------
+        zeropi_evecs:
+            If None (default), helper method _zeropi_operator_in_product_basis calculates zeropi eigenvectors and uses
+            them to convert operator to the product basis, if product basis is chosen.
+            If zeropi_evecs = zeropievecs, where zeropievecs is an ndarray, and product basis is chosen,
+            helper method _zeropi_operator_in_product_basis uses zeropievecs to convert the operator
+            to the product basis.
+        energy_esys:
+            If `False` (default), returns operator in the product basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            scipy.sparse.csc_matrix
+            Operator in chosen basis. If product basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return self._zeropi_operator_in_product_basis(
+        native = self._zeropi_operator_in_product_basis(
             self._zeropi.phi_operator(), zeropi_evecs=zeropi_evecs
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -503,5 +686,7 @@ def hilbertdim(self) -> int:
     ) -> ndarray:
-        """Returns a matrix of coupling strengths g_{ll'} [cmp. Dempster et al., text above Eq. (17)], using the states
-        from 'zeropi_states'. If `zeropi_states==None`, then a set of `self.zeropi` eigenstates is calculated. Only in
-        that case is `which` used for the eigenstate number (and hence the coupling matrix size).
+        """Returns a matrix of coupling strengths g_{ll'} [cmp. Dempster et al., text
+        above Eq. (17)], using the states from 'zeropi_states'. If
+        `zeropi_states==None`, then a set of `self.zeropi` eigenstates is calculated.
+        Only in that case is `which` used for the eigenstate number (and hence the
+        coupling matrix size).
         """
@@ -508,0 +693,0 @@ if evals_count is None:

scqubits/core/zeropi.py

@@ -370,14 +370,34 @@ # zeropi.py
 
-    def hamiltonian(self) -> csc_matrix:
-        """Calculates Hamiltonian in basis obtained by discretizing phi and employing
-        charge basis for theta.
+    def hamiltonian(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> csc_matrix:
+        r"""
+        Calculates Hamiltonian in basis obtained by discretizing :math:`\phi` and employing
+        charge basis for :math:`\theta` or in the eigenenergy basis. Returns matrix representing
+        the potential energy operator.
 
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns Hamiltonian in basis obtained by discretizing :math:`\phi` and employing
+            charge basis for :math:`\theta`.
+            If `True`, the energy eigenspectrum is computed, returns Hamiltonian in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns Hamiltonian in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            matrix representing the potential energy operator
+            Hamiltonian in chosen basis as csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, the Hamiltonian has dimensions of `truncated_dim`
+            x `truncated_dim`. Otherwise, if eigenenergy basis is chosen, Hamiltonian has dimensions of m x m,
+            for m given eigenvectors.
         """
-        return self.sparse_kinetic_mat() + self.sparse_potential_mat()
+        hamiltonian_mat = self.sparse_kinetic_mat() + self.sparse_potential_mat()
+        return self.process_hamiltonian(
+            native_hamiltonian=hamiltonian_mat, energy_esys=energy_esys
+        )
 
     def sparse_d_potential_d_flux_mat(self) -> csc_matrix:
-        r"""Calculates a of the potential energy w.r.t flux, at the current value of
+        r"""
+        Calculates a derivative of the potential energy w.r.t flux, at the current value of
         flux, as stored in the object.
@@ -405,15 +425,30 @@
 
-    def d_hamiltonian_d_flux(self) -> csc_matrix:
-        r"""Calculates a derivative of the Hamiltonian w.r.t flux, at the current value
-        of flux, as stored in the object.
+    def d_hamiltonian_d_flux(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Calculates a derivative of the Hamiltonian w.r.t flux, at the current value
+        of flux, as stored in the object. The flux is assumed to be given in the units
+        of the ratio :math:`\Phi_{ext}/\Phi_0`.
+        Returns matrix representing a derivative of the Hamiltonian in the native Hamiltonian basis
+        or eigenenergy basis.
 
-        The flux is assumed to be given in the units of the ratio \Phi_{ext}/\Phi_0.
-        So if \frac{\partial H}{ \partial \Phi_{\rm ext}}, is needed, the expr returned
-        by this function, needs to be multiplied by 1/\Phi_0.
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
 
         Returns
         -------
-            matrix representing the derivative of the Hamiltonian
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return self.sparse_d_potential_d_flux_mat()
+        native = self.sparse_d_potential_d_flux_mat()
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -433,20 +468,55 @@ def sparse_d_potential_d_EJ_mat(self) -> csc_matrix:
 
-    def d_hamiltonian_d_EJ(self) -> csc_matrix:
-        r"""Calculates a derivative of the Hamiltonian w.r.t EJ.
+    def d_hamiltonian_d_EJ(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Calculates a derivative of the Hamiltonian w.r.t EJ.
+        Returns matrix representing a derivative of the Hamiltonian in the native Hamiltonian basis
+        or eigenenergy basis.
 
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            matrix representing the derivative of the Hamiltonian
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return self.sparse_d_potential_d_EJ_mat()
+        native = self.sparse_d_potential_d_EJ_mat()
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def d_hamiltonian_d_ng(self) -> csc_matrix:
-        r"""Calculates a derivative of the Hamiltonian w.r.t ng.
-        as stored in the object.
+    def d_hamiltonian_d_ng(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
+        r"""
+        Calculates a derivative of the Hamiltonian w.r.t ng as stored in the object.
+        Returns matrix representing a derivative of the Hamiltonian in the native Hamiltonian basis
+        or eigenenergy basis.
 
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
         Returns
         -------
-            matrix representing the derivative of the Hamiltonian
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return -8 * self.EC * self.n_theta_operator()
+        native = -8 * self.EC * self.n_theta_operator()
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -488,11 +558,48 @@ def _identity_phi(self) -> csc_matrix:
 
-    def phi_operator(self) -> csc_matrix:
+    def phi_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
         r"""
-        Operator :math:`\phi`.
+        Returns :math:`\phi` operator in the native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return sparse.kron(self._phi_operator(), self._identity_theta(), format="csc")
+        native = sparse.kron(self._phi_operator(), self._identity_theta(), format="csc")
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
-    def n_theta_operator(self) -> csc_matrix:
+    def n_theta_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
         r"""
-        Operator :math:`n_\theta`.
+        Returns :math:`n_\theta` operator in the native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
@@ -504,3 +611,4 @@ dim_theta = 2 * self.ncut + 1
         ).tocsc()
-        return sparse.kron(self._identity_phi(), n_theta_matrix, format="csc")
+        native = sparse.kron(self._identity_phi(), n_theta_matrix, format="csc")
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -549,9 +657,28 @@ def _sin_phi_operator(self, x: float = 0) -> csc_matrix:
 
-    def cos_theta_operator(self) -> csc_matrix:
+    def cos_theta_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
         r"""
-        Operator :math:`\cos(\theta)`.
+        Returns :math:`\cos(\theta)` operator in the native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return sparse.kron(
+        native = sparse.kron(
             self._identity_phi(), self._cos_theta_operator(), format="csc"
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -568,6 +695,6 @@ def _sin_theta_operator(self) -> csc_matrix:
                 sparse.dia_matrix(
-                    ([1.0] * dim_theta, [1]), shape=(dim_theta, dim_theta)
+                    ([1.0] * dim_theta, [-1]), shape=(dim_theta, dim_theta)
                 )
                 - sparse.dia_matrix(
-                    ([1.0] * dim_theta, [-1]), shape=(dim_theta, dim_theta)
+                    ([1.0] * dim_theta, [1]), shape=(dim_theta, dim_theta)
                 )
@@ -578,9 +705,28 @@ ).tocsc()
 
-    def sin_theta_operator(self) -> csc_matrix:
+    def sin_theta_operator(
+        self, energy_esys: Union[bool, Tuple[ndarray, ndarray]] = False
+    ) -> Union[ndarray, csc_matrix]:
         r"""
-        Operator :math:`\sin(\theta)`.
+        Returns :math:`\sin(\theta)` operator in the native or eigenenergy basis.
+
+        Parameters
+        ----------
+        energy_esys:
+            If `False` (default), returns operator in the native basis.
+            If `True`, the energy eigenspectrum is computed, returns operator in the energy eigenbasis.
+            If `energy_esys = esys`, where esys is a tuple containing two ndarrays (eigenvalues and energy eigenvectors),
+            returns operator in the energy eigenbasis, and does not have to recalculate eigenspectrum.
+
+        Returns
+        -------
+            Operator in chosen basis. If native basis chosen, operator
+            returned as a csc_matrix. If the eigenenergy basis is chosen,
+            unless `energy_esys` is specified, operator has dimensions of `truncated_dim`
+            x truncated_dim, and is returned as an ndarray. Otherwise, if eigenenergy basis is chosen,
+            operator has dimensions of m x m, for m given eigenvectors, and is returned as an ndarray.
         """
-        return sparse.kron(
+        native = sparse.kron(
             self._identity_phi(), self._sin_theta_operator(), format="csc"
         )
+        return self.process_op(native_op=native, energy_esys=energy_esys)
 
@@ -587,0 +733,0 @@ def plot_potential(

scqubits/explorer/explorer_panels.py

@@ -15,4 +15,6 @@ # explorer_panels.py
 
+# import matplotlib as mlp
 import numpy as np
 
+from matplotlib import rc_context
 from matplotlib.axes import Axes
@@ -27,2 +29,3 @@ from matplotlib.figure import Figure
 from scqubits.core.oscillator import Oscillator
+from scqubits.settings import matplotlib_settings
 from scqubits.utils.misc import tuple_to_short_str
@@ -35,2 +38,3 @@
 
+@rc_context(matplotlib_settings)
 def display_bare_spectrum(
@@ -61,2 +65,3 @@ sweep: "ParameterSweep",
 
+@rc_context(matplotlib_settings)
 def display_anharmonicity(
@@ -82,2 +87,3 @@ sweep: "ParameterSweep",
 
+@rc_context(matplotlib_settings)
 def display_matrixelements(
@@ -111,2 +117,3 @@ sweep: "ParameterSweep",
 
+@rc_context(matplotlib_settings)
 def display_matrixelement_sweep(
@@ -158,2 +165,3 @@ sweep: "ParameterSweep",
 
+@rc_context(matplotlib_settings)
 def display_bare_wavefunctions(
@@ -180,2 +188,3 @@ sweep: "ParameterSweep",
 
+@rc_context(matplotlib_settings)
 def display_transitions(
@@ -206,2 +215,3 @@ sweep: "ParameterSweep",
 
+@rc_context(matplotlib_settings)
 def display_cross_kerr(
@@ -219,3 +229,3 @@ sweep: "ParameterSweep",
         title = f"ac Stark: {subsys1.id_str} + {subsys2.id_str}"
-        ylabel = f"ac Stark shift $\chi^{{{subsys1_index},{subsys2_index}}}$"
+        ylabel = rf"ac Stark shift $\chi^{{{subsys1_index},{subsys2_index}}}$"
         levels_list = [1]
@@ -263,2 +273,3 @@ kerr_data = sweep["chi"][subsys1_index, subsys2_index]
 
+@rc_context(matplotlib_settings)
 def display_self_kerr(
@@ -265,0 +276,0 @@ sweep: "ParameterSweep",

scqubits/settings.py

@@ -96,28 +96,37 @@ """
 # Matplotlib options -------------------------------------------------------------------
-# set custom matplotlib color cycle
-mpl.rcParams["axes.prop_cycle"] = cycler(
-    color=[
-        "#016E82",
-        "#333795",
-        "#2E5EAC",
-        "#4498D3",
-        "#CD85B9",
-        "#45C3D1",
-        "#AA1D3F",
-        "#F47752",
-        "#19B35A",
-        "#EDE83B",
-        "#ABD379",
-        "#F9E6BE",
-    ]
-)
+# set matplotlib defaults for use in @mpl.rc_context
+off_black = "0.2"
+matplotlib_settings = {
+    "axes.prop_cycle": cycler(
+        color=[
+            "#016E82",
+            "#333795",
+            "#2E5EAC",
+            "#4498D3",
+            "#CD85B9",
+            "#45C3D1",
+            "#AA1D3F",
+            "#F47752",
+            "#19B35A",
+            "#EDE83B",
+            "#ABD379",
+            "#F9E6BE",
+        ]
+    ),
+    "font.family": "IBM Plex Sans, Roboto, Arial, Helvetica, DejaVu Sans",
+    "font.size": 11,
+    "font.weight": 500,
+    "axes.labelsize": 11,
+    "axes.titlesize": 11,
+    "xtick.labelsize": 10,
+    "ytick.labelsize": 10,
+    "xtick.labelcolor": off_black,
+    "ytick.labelcolor": off_black,
+    "xtick.color": off_black,
+    "ytick.color": off_black,
+    "axes.labelcolor": off_black,
+    "axes.edgecolor": off_black,
+    "axes.titlecolor": off_black
+}
 
-# set matplotlib defaults
-mpl.rcParams["font.family"] = "sans-serif"
-mpl.rcParams["font.sans-serif"] = "Roboto, Arial, Helvetica, DejaVu Sans"
-mpl.rcParams["font.size"] = 11
-mpl.rcParams["axes.labelsize"] = 11
-mpl.rcParams["axes.titlesize"] = 11
-mpl.rcParams["xtick.labelsize"] = 10
-mpl.rcParams["ytick.labelsize"] = 10
 
@@ -124,0 +133,0 @@ # toggle top and right axes on and off

scqubits/tests/test_hilbertspace.py

@@ -18,2 +18,3 @@ # test_hilbertspace.py
 import scqubits as scq
+import qutip as qt
 
@@ -387,1 +388,88 @@ from scqubits.core.hilbert_space import HilbertSpace
         hilbertspace_new = scq.HilbertSpace.create()
+
+    def test_HilbertSpace_op_in_dressed_basis(self):
+        E_osc_a = 4.0
+        E_osc_b = 3.2
+        g = 0.01
+        Delta = E_osc_a - E_osc_b
+        truncated_dim = 4
+        theta = 0.5 * np.arctan(2 * g / Delta)
+        osc_a = scq.Oscillator(E_osc=E_osc_a, truncated_dim=truncated_dim)
+        osc_b = scq.Oscillator(E_osc=E_osc_b, truncated_dim=truncated_dim)
+        hilbert_space = scq.HilbertSpace([osc_a, osc_b])
+        hilbert_space.add_interaction(
+            g=g,
+            op1=osc_a.creation_operator,
+            op2=osc_b.annihilation_operator,
+            add_hc=True,
+        )
+        hilbert_space.generate_lookup()
+        hilbert_space.standardize_eigenvector_phases()
+        # analytic answer for the dressed operator based on a Bogoliubov transformation
+        a_id_wrap = scq.utils.spectrum_utils.identity_wrap(
+            osc_a.annihilation_operator(), osc_a, hilbert_space.subsystem_list
+        )
+        b_id_wrap = scq.utils.spectrum_utils.identity_wrap(
+            osc_b.annihilation_operator(), osc_b, hilbert_space.subsystem_list
+        )
+        analytic_op = np.cos(theta) * a_id_wrap - np.sin(theta) * b_id_wrap
+        # need to order this operator according to the dressed indices for later
+        # comparison with operators expressed in the dressed basis
+        ordered_bare_indices = [
+            hilbert_space.bare_index(idx) for idx in range(truncated_dim**2)
+        ]
+        ordered_basis_states = [
+            qt.tensor(qt.basis(truncated_dim, idx_a), qt.basis(truncated_dim, idx_b))
+            for (idx_a, idx_b) in ordered_bare_indices
+        ]
+        # consider only matrix elements unaffected by the truncation level
+        analytic_op_ordered = qt.Qobj(
+            analytic_op.transform(ordered_basis_states)[0:10, 0:10]
+        )
+        op1 = qt.Qobj(
+            hilbert_space.op_in_dressed_eigenbasis(op=osc_a.annihilation_operator)[
+                0:10, 0:10
+            ]
+        )
+        op2 = qt.Qobj(
+            hilbert_space.op_in_dressed_eigenbasis(
+                op=(osc_a.annihilation_operator(), osc_a)
+            )[0:10, 0:10]
+        )
+        op3 = qt.Qobj(
+            hilbert_space.op_in_dressed_eigenbasis(
+                op=(osc_a.annihilation_operator(), osc_a), op_in_bare_eigenbasis=True
+            )[0:10, 0:10]
+        )
+        op4 = qt.Qobj(
+            hilbert_space.op_in_dressed_eigenbasis(
+                op=(osc_a.annihilation_operator(), osc_a), op_in_bare_eigenbasis=False
+            )[0:10, 0:10]
+        )
+        assert analytic_op_ordered == op1
+        assert analytic_op_ordered == op2
+        assert analytic_op_ordered == op3
+        assert analytic_op_ordered == op4
+
+    def test_HilbertSpace_op_in_dressed_basis_native_vs_bare_basis(self):
+        E_osc = 4.0
+        g = 0.01
+        truncated_dim = 4
+        tmon = scq.Transmon(
+            EJ=10.0, EC=0.2, ng=0.0, ncut=15, truncated_dim=truncated_dim
+        )
+        osc = scq.Oscillator(E_osc=E_osc, truncated_dim=truncated_dim)
+        hilbert_space = scq.HilbertSpace([tmon, osc])
+        hilbert_space.add_interaction(
+            g=g,
+            op1=tmon.n_operator,
+            op2=osc.annihilation_operator,
+            add_hc=True,
+        )
+        hilbert_space.generate_lookup()
+        op1 = hilbert_space.op_in_dressed_eigenbasis(op=tmon.n_operator)
+        n_op_bare_eigenbasis_v2 = tmon.n_operator(energy_esys=True)
+        op2 = hilbert_space.op_in_dressed_eigenbasis(
+            op=(n_op_bare_eigenbasis_v2, tmon), op_in_bare_eigenbasis=True
+        )
+        assert op1 == op2

scqubits/ui/explorer_widget.py

@@ -39,2 +39,3 @@ # explorer_widget.py
 from scqubits.utils import misc as utils
+from scqubits.settings import matplotlib_settings
 
@@ -94,3 +95,2 @@ if TYPE_CHECKING:
         width=str(pixels) + "px",
-        align="top",
         border="1px solid lightgrey",
@@ -169,2 +169,3 @@ padding="10px 10px 10px 10px",
 
+    @matplotlib.rc_context(matplotlib_settings)
     def build_figure_and_axes_table(self) -> Tuple[Figure, np.ndarray]:
@@ -430,2 +431,3 @@ # the %inline and %widget backends somehow scale differently; try to compensate
 
+    @matplotlib.rc_context(matplotlib_settings)
     def build_ui_figure_display(self):
@@ -445,2 +447,3 @@ if _HAS_WIDGET_BACKEND:
 
+    @matplotlib.rc_context(matplotlib_settings)
     def display_panel(
@@ -662,2 +665,3 @@ self,
 
+    @matplotlib.rc_context(matplotlib_settings)
     def update_layout_and_plots(self: "Explorer", change):
@@ -696,2 +700,3 @@ panels = self.get_panels_list()
 
+    @matplotlib.rc_context(matplotlib_settings)
     def update_plots(self: "Explorer", change):
@@ -822,6 +827,3 @@ if not hasattr(self, "fig"):
             self.ui["transitions"]["initial_bare_dressed_toggle"] = ToggleButtons(
-                options=["bare", "dressed"],
-                value="bare",
-                description="",
-                disable=False,
+                options=["bare", "dressed"], value="bare", description=""
             )
@@ -828,0 +830,0 @@ self.ui["transitions"][

scqubits/ui/qubit_widget.py

@@ -68,3 +68,2 @@ # qubit_widget.py
             description="",
-            disabled=False,
             layout=ipywidgets.Layout(width="150px"),
@@ -71,0 +70,0 @@ )

scqubits/utils/misc.py

@@ -165,2 +165,16 @@ # misc.py
 
+def check_lookup_exists(func: Callable) -> Callable:
+    @functools.wraps(func)
+    def wrapper(self, *args, **kwargs):
+        if not self._lookup_exists:
+            raise Exception(
+                "Lookup data not found. For HilbertSpace: data must be generated with .generate_lookup(). "
+                "For ParameterSweep: data should be automatically generated unless disabled manually. In "
+                "the latter case, apply .run()."
+            )
+        return func(self, *args, **kwargs)
+
+    return wrapper
+
+
 class DeprecationMessage:
@@ -167,0 +181,0 @@ """Decorator class, producing an adjustable warning and info upon usage of the

scqubits/utils/plot_defaults.py

@@ -61,3 +61,3 @@ # plot_defaults.py
 
-    # If energy spacing is used for scaling, fill no more than this  fraction of the
+    # If energy spacing is used for scaling, fill no more than this fraction of the
     # spacing.
@@ -79,4 +79,4 @@ FILLING_FRACTION = 0.9
     e_range = e_max - e_min
-    y_min = np.min(potential_vals)
-    y_max = e_max + 0.3 * e_range
+    y_min = np.min(potential_vals)  # lowest value of potential energy
+    y_max = e_max + 0.3 * e_range  # maximum eigenenergy plus padding
     y_range = y_max - y_min
@@ -83,0 +83,0 @@

scqubits/utils/plot_utils.py

@@ -27,2 +27,3 @@ # plot_utils.py
 from scqubits.utils import plot_defaults as defaults
+from scqubits.settings import matplotlib_settings
 
@@ -50,2 +51,3 @@ if TYPE_CHECKING:
 
[e-mail]_context(matplotlib_settings)
 def _extract_kwargs_options(
@@ -83,2 +85,3 @@ kwargs: Dict[str, Any], plot_type: str, direct_plot_options: Dict[str, Any] = None
 
[e-mail]_context(matplotlib_settings)
 def _process_options(
@@ -129,2 +132,3 @@ figure: Figure, axes: Axes, opts: Dict[str, Any] = None, **kwargs
 
[e-mail]_context(matplotlib_settings)
 def _process_special_option(figure: Figure, axes: Axes, key: str, value: Any) -> None:
@@ -148,2 +152,3 @@ """Processes a single 'special' option, i.e., one internal to scqubits and not to be
 
[e-mail]_context(matplotlib_settings)
 def despine_axes(axes: Axes) -> None:
@@ -159,2 +164,3 @@ # Hide the right and top spines
 
[e-mail]_context(matplotlib_settings)
 def scale_wavefunctions(
@@ -175,2 +181,3 @@ wavefunc_list: List["WaveFunction"],
 
[e-mail]_context(matplotlib_settings)
 def plot_wavefunction_to_axes(
@@ -189,2 +196,3 @@ axes: Axes, wavefunction: "WaveFunction", energy_offset: float, **kwargs
 
[e-mail]_context(matplotlib_settings)
 def plot_potential_to_axes(
@@ -210,2 +218,3 @@ axes: Axes,
 
[e-mail]_context(matplotlib_settings)
 def add_numbers_to_axes(
@@ -228,2 +237,3 @@ axes: Axes, matrix: ndarray, modefunc: Callable, fontsize: int = 8
 
[e-mail]_context(matplotlib_settings)
 def color_normalize(vals, mode: str) -> Tuple[float, float, mpl.colors.Normalize]:
@@ -230,0 +240,0 @@ minval = min(vals)

scqubits/utils/plotting.py

@@ -36,2 +36,3 @@ # plotting.py
 )
+from scqubits.settings import matplotlib_settings
 
@@ -49,2 +50,3 @@ if TYPE_CHECKING:
 
[e-mail]_context(matplotlib_settings)
 def wavefunction1d(
@@ -94,2 +96,3 @@ wavefuncs: Union["WaveFunction", "List[WaveFunction]"],
 
[e-mail]_context(matplotlib_settings)
 def wavefunction1d_nopotential(
@@ -133,2 +136,3 @@ wavefuncs: Union["WaveFunction", "List[WaveFunction]"],
 
[e-mail]_context(matplotlib_settings)
 def wavefunction1d_discrete(wavefunc: "WaveFunction", **kwargs) -> Tuple[Figure, Axes]:
@@ -163,2 +167,3 @@ """
 
[e-mail]_context(matplotlib_settings)
 def wavefunction2d(
@@ -217,2 +222,3 @@ wavefunc: "WaveFunctionOnGrid", zero_calibrate: bool = False, **kwargs
 
[e-mail]_context(matplotlib_settings)
 def contours(
@@ -270,2 +276,3 @@ x_vals: Union[List[float], np.ndarray],
 
[e-mail]_context(matplotlib_settings)
 def matrix(
@@ -309,2 +316,3 @@ data_matrix: np.ndarray, mode: str = "abs", show_numbers: bool = False, **kwargs
 
[e-mail]_context(matplotlib_settings)
 def matrix_skyscraper(
@@ -368,2 +376,6 @@ matrix: np.ndarray, mode: str = "abs", **kwargs
 
+    axes.tick_params(axis='x', pad=-5)
+    axes.tick_params(axis='y', pad=-5)
+    axes.tick_params(axis='z', pad=-2)
+
     _process_options(fig, axes, opts=defaults.matrix(), **kwargs)
@@ -374,2 +386,3 @@
 
[e-mail]_context(matplotlib_settings)
 def matrix2d(
@@ -425,3 +438,3 @@ matrix: np.ndarray,
         box_height_inches = fig_height / matrix.shape[0]
-        font_size = min(box_width_inches, box_height_inches) * 12
+        font_size = min(box_width_inches, box_height_inches) * 11
         add_numbers_to_axes(axes, matrix, modefunction, fontsize=font_size)
@@ -444,2 +457,3 @@
 
[e-mail]_context(matplotlib_settings)
 def data_vs_paramvals(
@@ -503,2 +517,3 @@ xdata: np.ndarray,
 
[e-mail]_context(matplotlib_settings)
 def evals_vs_paramvals(
@@ -550,2 +565,3 @@ specdata: "SpectrumData",
 
[e-mail]_context(matplotlib_settings)
 def matelem_vs_paramvals(
@@ -591,3 +607,3 @@ specdata: "SpectrumData",
 
-    for (row, col) in index_pairs:
+    for row, col in index_pairs:
         y_vals = modefunction(specdata.matrixelem_table[:, row, col])
@@ -604,4 +620,4 @@ axes.plot(
     else:
-        axes.legend(loc="center left", bbox_to_anchor=(1, 0.5))
+        axes.legend(loc="center left", bbox_to_anchor=(1, 0.5), frameon=False)
     _process_options(fig, axes, opts=defaults.matelem_vs_paramvals(specdata), **kwargs)
     return fig, axes

scqubits/utils/spectrum_utils.py

@@ -121,4 +121,4 @@ # spectrum_utils.py
 
-    Summing up to the midpoint only is to address the  danger that the sum is
-    actually zero, which may is the case for odd wavefunctions taken over an interval
+    Summing up to the midpoint only is to address the danger that the sum is
+    actually zero, which may be the case for odd wavefunctions taken over an interval
     centered at zero.
@@ -125,0 +125,0 @@ """

scqubits/version.py

 # THIS FILE IS GENERATED FROM scqubits SETUP.PY
-short_version = '3.1.0'
-version = '3.1.0'
+short_version = '3.1.1'
+version = '3.1.1'
 release = True

setup.py

@@ -56,3 +56,3 @@ """scqubits: superconducting qubits in Python
 MINOR = 1
-MICRO = 0
+MICRO = 1
 ISRELEASED = True
@@ -59,0 +59,0 @@

Alert delta unavailable

Currently unable to show alert delta for PyPI packages.

Improved metrics

Lines of code

25235

5.19%

Worsened metrics

Total package byte prevSize

8149033

-5.39%

Number of package files

121

-0.82%