diff --git a/examples/tutorials/APS2026/Diamond_NV_PyCCE_APS2026.ipynb b/examples/tutorials/APS2026/Diamond_NV_PyCCE_APS2026.ipynb
new file mode 100644
index 0000000..e37f33e
--- /dev/null
+++ b/examples/tutorials/APS2026/Diamond_NV_PyCCE_APS2026.ipynb
@@ -0,0 +1 @@
+{"cells":[{"cell_type":"markdown","metadata":{"pycharm":{"name":"#%% md\n"},"id":"pBT83BxJ69ko"},"source":["# **NV Center in Diamond**\n","\n","In this tutorial we will go over the main steps of running CCE calculations for the NV center in diamond with the **PyCCE** module. Those include:\n","\n","* Generating the spin bath using the `pycce.BathCell` instance.\n","* Setting up properties of the `pycce.Simulator` instance.\n","* Running the calculations with the `Simulator.compute` function.\n","\n","We will compute the Hahn-echo coherence function (with decoupling $\\pi$-pulse) using conventional CCE."]},{"cell_type":"markdown","source":["## **(0) Import packages**"],"metadata":{"id":"oGPlaQHFVg-R"}},{"cell_type":"code","source":["# Version control\n","!pip install scipy==1.13.0\n","!pip install numpy==1.26.0"],"metadata":{"id":"5Srgo2Ot88sB"},"execution_count":null,"outputs":[]},{"cell_type":"code","execution_count":null,"metadata":{"pycharm":{"name":"#%%\n"},"id":"mT7tFLXV69kp"},"outputs":[],"source":["!pip install pycce\n","!pip install ase\n","\n","import numpy as np\n","import matplotlib.pyplot as plt\n","import pandas as pd\n","\n","import sys\n","import pycce as pc\n","import ase\n","\n","from mpl_toolkits import mplot3d\n","\n","seed = 8805\n","np.random.seed(seed)\n","np.set_printoptions(suppress=True, precision=5)"]},{"cell_type":"markdown","source":["# **=======================================================================**"],"metadata":{"id":"dFiewWMWH3BV"}},{"cell_type":"markdown","metadata":{"pycharm":{"name":"#%% md\n"},"id":"lG1jgJm569kq"},"source":["## **(1) Generate the nuclear spin bath**\n","\n","There are 4 steps in generating the nuclear spin bath (i.e. the <sup>13</sup>C nuclear spins around the NV center).\n"]},{"cell_type":"markdown","source":["### **(1.1) Create BathCell**"],"metadata":{"id":"vQSQg5fHPp_U"}},{"cell_type":"markdown","source":["We first build the `BathCell` object in PyCCE, which is a precursor to generating the nuclear spin bath.\n","\n","We can generate the `BathCell` object from an `ase.Atoms` object by using the classmethod `BathCell.from_ase`."],"metadata":{"id":"zpX_VoxcPia6"}},{"cell_type":"code","execution_count":null,"metadata":{"id":"my63osab69kq"},"outputs":[],"source":["from ase.build import bulk\n","\n","# Create diamond unit cell from ase\n","diamond_ase = bulk('C', 'diamond', cubic=True)  # ASE Atoms object\n","\n","# Generate PyCCE BathCell\n","diamond_bc = pc.read_ase(diamond_ase)           # PyCCE BathCell object"]},{"cell_type":"markdown","metadata":{"id":"xD-2Fsfn69kr"},"source":["The following attributes are created with this initiallization:\n","\n","* `.cell` is an ndarray containing lattice vector information. Each **column** is a lattice vector in cartesian coordinates.\n","* `.atoms` is a dictionary with keys corresponding to the atom name, and each item is a list of the coordinates in cell coordinates."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"JKatsNXI69kr"},"outputs":[],"source":["print('Cell\\n', diamond_bc.cell)\n","print('\\nAtoms\\n', diamond_bc.atoms)"]},{"cell_type":"markdown","source":["### **(1.2) Specify isotopes in BathCell**"],"metadata":{"id":"yW_9oYmHPv39"}},{"cell_type":"markdown","metadata":{"id":"eIqO0_kP69ks"},"source":["By default, PyCCE uses the **EasySpin** database (https://easyspin.org/easyspin/documentation/isotopetable.html) containing the concentrations of all common stable isotopes of all elements, however the user can proide custom concentrations.\n","\n","Use the function `BathCell.add_isotopes` to add one (or several) isotopes of the element. Each isotope is initiallized with a ``tuple`` containing the name of the isotope and its concentration.\n","\n","Structure of the dictionary-like object:\n","```\n","{element_1: {isotope_1: concentration, isotope_2: concentration},\n"," element_2: {isotope_3: concentration, ...}}\n","```\n"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"0YU2wYk169ks"},"outputs":[],"source":["# Add types of isotopes\n","diamond_bc.add_isotopes(('13C', 0.011))"]},{"cell_type":"markdown","metadata":{"id":"U6wK851h69ks"},"source":["Isotopes may also be directly added to `BathCell.isotopes`. For example, below we are adding an isotope without a nuclear spin:"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"lv4xW3RG69ks"},"outputs":[],"source":["diamond_bc.isotopes['C']['14C'] = 0.001"]},{"cell_type":"markdown","source":["### **(1.3) Define the quantization axis of the spin defect**"],"metadata":{"id":"YPZPKUdnP2Mv"}},{"cell_type":"markdown","metadata":{"id":"iW0ViLdQ69ks"},"source":["In the `Simulator` object, everything is set in $S_z$ basis. When the quantization axis of the defect does not align with the (0, 0, 1) direction of the crystal axis, the user needs to manually define this axis.\n","\n","If one wants to specify the complete rotation of Cartesian axes, one can provide a rotation matrix to rotate the Cartesian reference frame with respect to the cell coordinates by calling the `BathCell.rotate` method."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"mh8D87RF69ks"},"outputs":[],"source":["# set z direction of the defect\n","diamond_bc.zdir = [1, 1, 1]"]},{"cell_type":"markdown","source":["### **(1.4) Generate the spin bath**"],"metadata":{"id":"CqYxT34D6cFf"}},{"cell_type":"markdown","metadata":{"id":"fJ2LS9F469ks"},"source":["To generate the spin bath, use the `BathCell.gen_supercell` method.**This is when the NV center is created!**\n","\n","`sc_size` is the linear size of the supercell in &Aring;.\n","\n","`remove` takes a ``tuple`` or ``list of tuples`` as an argument. First element of each tuple is the name of the **atom**, and the second element is the fractional coordinate. If such atoms are found in the supercell, they are removed.\n","\n","`add` takes a ``tuple`` or ``list of tuples`` as an argument. First element of each tuple is the name of the **isotope**, and the second element is the fractional coordinate. Each of the specified isotopes will be added in the final supercell at the specified locations."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"hmtu6WAn69kt"},"outputs":[],"source":["# Add the defect: remove and add atoms at positions (in cell coordinates)\n","sc_size = 200  # linear size of supercell (A)\n","diamond_bath = diamond_bc.gen_supercell(sc_size,\n","                                        remove=[('C', [0., 0, 0]),\n","                                                ('C', [0.25, 0.25, 0.25])],\n","                                        add=('14N', [0.5, 0.25, 0.25]),\n","                                        seed=seed)"]},{"cell_type":"markdown","metadata":{"id":"BlCOqLbp69kt"},"source":["Note, that because the `14C` isotope doesn't have a spin, **PyCCE** does not find it in common isotopes, and raises a warning. We have to provide `SpinType` for it separately, or define the properties as follows:"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"A2MWI_T569kt"},"outputs":[],"source":["diamond_bath['14C'].gyro = 0\n","diamond_bath['14C'].spin = 0"]},{"cell_type":"markdown","source":["Note that the `BathCell.gen_supercell` method creates a `BathArray` object. This will be discussed in (1.6)."],"metadata":{"id":"2zVqbx_yLZt2"}},{"cell_type":"markdown","source":["### **(1.5) Visualize the nuclear spin bath**"],"metadata":{"id":"UZJDWjd_C41t"}},{"cell_type":"markdown","source":["We can use `matplotlib` to visualize the spatial distribution of the spin bath."],"metadata":{"id":"DUzzHGOPDG5u"}},{"cell_type":"code","source":["# add 3D axis\n","fig = plt.figure(figsize=(6,6))\n","ax = fig.add_subplot(projection='3d')\n","\n","# plot the positions of the bath\n","ax.scatter3D(diamond_bath.x, diamond_bath.y, diamond_bath.z, c='k')\n","\n","ax.set(xlabel='x (A)', ylabel='y (A)', zlabel='z (A)');"],"metadata":{"id":"MDBOC9jc8b7i"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["It looks like a dirty blob; **that's okay**! The `BathCell.gen_supercell` routine simply creates all nuclear spins within the generated supercell (`sc_size`) and with the concentration of each isotope; **the clusters are not created yet**.\n","\n","Later, we'll define bath parameters `r_bath`, `r_dipole`, and `order`, which will generate the spin bath with clusters."],"metadata":{"id":"fa0ZzCeDERwL"}},{"cell_type":"markdown","source":["**Check your understanding**: How many <sup>13</sup>C nuclear spins should we expect to be generated by the `gen_supercell` routine?"],"metadata":{"id":"8JJwgmjyIG3e"}},{"cell_type":"markdown","source":["### **(1.6) The BathArray structure**"],"metadata":{"id":"veF-8il7PaCs"}},{"cell_type":"markdown","metadata":{"id":"mW-XXq8c69kt"},"source":["The bath spins are stored in the `BathArray` object - a subclass of `np.ndarray` with fixed datastructure:\n","\n","* `N` field `dtype('<U16')` contains the names of bath spins.\n","* `xyz` field `dtype('<f8', (3,))` contains the positions of bath spins (in &Aring;).\n","* `A` field `dtype('<f8', (3, 3))` contains the hyperfine coupling of bath spins (in kHz).\n","* `Q` field `dtype('<f8', (3, 3))` contains the quadrupole tensor of bath spins (in kHz) (Relevant for spin &ge; 1).\n","\n","All of the fields are accesible as attributes of `BathArray`.\n","Additionally, the subarrays of the specific spins are accessible with their name as indicated above.\n","\n","Upon generation of the array from the cell, the `Q` and `A` fields are empty. The Hyperfine couplings `A` will be automatically computed by the `Simulator` object, however the quadrupole couplings `Q` must be set by the user.\n","\n","The additional attributes allow one to access `SpinType` properties:\n","\n","+ `name` returns the spin name or array of spin names;\n","+ `spin` returns the value of the spin or array of ones;\n","+ `gyro` returns gyromagnetic ratios of the spins;\n","+ `q` returns quadrupole constants of the spins;\n","+ `h` returns a dictionary with user-defined additions to the Hamiltonian.\n","+ `detuning` returns detunings of the spins (See definition below).\n","\n","For example, below we print out the attributes of the first two spins (not all) in the `BathArray`."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"SErPI_RU69kt"},"outputs":[],"source":["print('Names\\n', diamond_bath[:2].N)\n","print('\\nCoordinates\\n', diamond_bath[:2].xyz)\n","print('\\nHyperfine tensors\\n', diamond_bath[:2].A)\n","print('\\nQuadrupole tensors\\n', diamond_bath[:2].Q)"]},{"cell_type":"markdown","metadata":{"id":"KhWMR90r69kt"},"source":["The properties of spin types (gyromagnetic ratio, quadrupole moment, etc.) are stored in the `BathArray.types` attribute, which is an instance of `SpinDict` containing `SpinType` classes. For most known isotopes, `SpinType` can be found in the `pycce.common_isotopes` dictionary, and is set by default (including electrons, which can be set with `N = e`). The user can add additional `SpinType` objects, by calling `BathArray.add_type` method or setting elements of `SpinDict` directly. For details of the first approach see documentation of `SpinDict.add_type` method.\n","\n","The direct setting of types is rather simple. The user can set elements of `SpinDict` with ``tuple``, containing:\n","\n","* **(spin, gyromagnetic ratio, quadrupole moment *(optional)*, detuning *(optional)*, )**\n","\n","  **OR**\n","  \n","* **(isotope, spin, gyromagnetic ratio, quadrupole moment *(optional)*, detuning *(optional)*, )**\n","\n","where:\n","\n","- **isotope** (*str*) is the name of the given spin (same one as in `N` field of `BathArray`) to define new `SpinType` object.\n","  The key of `SpinDict` **has** to be the correct name of the spin (\"isotope\" field in the tuple).\n","- **spin** (*float*) is the total spin of the given bath spin.\n","- **gyromagnetic ratio** (*float*) is the gyromagnetic ratio of the given bath spin.\n","- **quadrupole moment** (*float*) is the quadrupole moment of the given bath spin. Relevant only when electric field gradient are used\n","  to generate quadrupole couplings for spins, stored in the ``BathArray``, with ``BathArray.from_efg`` method.\n","- **detuning** (float) is an additional energy splitting for model spins, included as an extra $+\\omega \\hat S_z$ term in the Hamiltonian,\n","  where $\\omega$ is the detuning.\n","\n","Units of gyromagnetic ratio are rad / ms / G, quadrupole moments are given in barn, detunings are given in kHz."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"WT8kpnOh69kt"},"outputs":[],"source":["# Several ways to set SpinDict elements\n","diamond_bath.types['14C'] = 0, 0, 0\n","diamond_bath.types['Y'] = ('Y', 0, 0, 0)\n","diamond_bath.types['A'] = pc.SpinType('A', 0, 0, 0)\n","\n","print(diamond_bath.types)"]},{"cell_type":"markdown","source":["### **(1.7) *Try it yourself* : V<sub>B</sub><sup>-</sup> in hBN**"],"metadata":{"id":"-2JNxxFEDQkL"}},{"cell_type":"markdown","source":["Let's try generating a nuclear spin bath representative of the environment around the V<sub>B</sub><sup>-</sup> center in hBN.\n","\n","Below, we've provided the `ase.Atoms` object for the hBN lattice."],"metadata":{"id":"ArnUJrIBFmh3"}},{"cell_type":"code","source":["from ase import Atoms\n","\n","a = 2.504\n","c = 6.661\n","\n","cell = [\n","    [a, 0, 0],\n","    [-a/2, np.sqrt(3)*a/2, 0],\n","    [0, 0, c]\n","]\n","\n","symbols = ['B', 'N', 'B', 'N']\n","\n","scaled_positions = [\n","    [1/3, 2/3, 1/4],\n","    [2/3, 1/3, 1/4],\n","    [2/3, 1/3, 3/4],\n","    [1/3, 2/3, 3/4]\n","]\n","\n","hbn_bulk_ase = Atoms(\n","    symbols=symbols,\n","    scaled_positions=scaled_positions,\n","    cell=cell,\n","    pbc=True\n",")\n","\n","print(hbn_bulk_ase)\n"],"metadata":{"id":"g3qxvr3RDtgt"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Begin here\n"],"metadata":{"id":"iJnMxfWnMMc2"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# **=======================================================================**"],"metadata":{"id":"zJ9fjxYpLLDz"}},{"cell_type":"markdown","source":["## **(2) Run a CCE simulation**"],"metadata":{"id":"Fsu0BBDT4JB2"}},{"cell_type":"markdown","source":["Here we describe how to run CCE simulations."],"metadata":{"id":"TnwQJnr-VwOD"}},{"cell_type":"markdown","source":["### **(2.1) The `Simulator` class**"],"metadata":{"id":"0aa3pdW2V5xr"}},{"cell_type":"markdown","metadata":{"pycharm":{"name":"#%% md\n"},"id":"0809UK_z69kt"},"source":["The `Simulator` object enables CCE simulations.\n","\n","Main parameters to consider:\n","\n","* `spin` — Either an instance of the `CenterArray` object, or float - total spin of the central spin (assuming one central spin).\n","\n","* `position` — Array of Cartesian coordinates (in &Aring;) of the central spin(s).\n","\n","* `bath` — spin bath in any specified format. Can be either:\n","\n","  - Instance of `BathArray` class\n","  - `ndarray` with ``dtype([('N', np.str_, 16), ('xyz', np.float64, (3,))])`` containing names of bath spins (same ones as stored in self.ntype) and positions of the spins in angstroms\n","  - The name of the `.xyz` text file containing 4 columns: name of the bath spin and xyz coordinates in &Aring;\n","\n","* `r_bath` — cutoff radius (in &Aring;) around the central spin for the bath.\n","\n","* `r_dipole` — cutoff radius (in &Aring;) for the pairwise distance to consider two nuclear spins to be connected.\n","\n","* `order` — maximum size of nuclear spin clusters.\n","\n","* `magnetic_field` — applied magnetic field (in Gauss). Can also be provided during the simulation run.\n","\n","* `pulses` — number of pulses in Carr-Purcell-Meiboom-Gill (CPMG) sequence or the pulse sequence itself.\n","\n","For the full description see the documentation of the `Simulator` object (https://pycce.readthedocs.io/en/stable/simulator.html)."]},{"cell_type":"markdown","source":["### **(2.2) Visualizing the nuclear spin clusters:**"],"metadata":{"id":"bw6L1NvA406L"}},{"cell_type":"markdown","metadata":{"id":"0THSKlsn69kt"},"source":["Let's set up a \"mock\" instance of `Simulator` to visualize the smaller part of the bath around the central spin."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"u1M5pA_w69kt"},"outputs":[],"source":["# Setting the runner engine\n","mock = pc.Simulator(spin=1,\n","                    position=[0,0,0],\n","                    bath=diamond_bath,\n","                    r_bath=20,\n","                    r_dipole=6,\n","                    order=3\n","                    )"]},{"cell_type":"markdown","metadata":{"id":"FTNl2jRY69kt"},"source":["During the initiallization, depending on the provided keyword arguments several methods may be called:\n","\n","* `Simulator.read_bath` is called if keyword `bath` is provided. It may take several additional arguments:\n","\n","    * `r_bath` - cutoff distance from the qubit for the bath.\n","    * `skiprows` - if `bath` is provided as `.xyz` file, this argument tells how many rows to skip when reading the file.\n","    * `external_bath` - `BathArray` instance, which contains bath spins with pre defined hyperfines to be used.\n","    * `hyperfine` - defines the way to compute hyperfine couplings. If it is not given and `bath` doesn't contain any predefined hyperfines\n","      (`bath['A'].any() == False`) the point dipole approximation is used.\n","      Otherwise it can be an instance of ``pc.Cube`` object, or callable with signature ``func(coord, gyro, central_gyro)``,\n","      where `coord` is an array of the bath spin coordinate, gyro is the gyromagnetic ratio of bath spin,\n","      central_gyro is the gyromagnetic ratio of the central bath spin.\n","    * `types` - instance of `SpinDict` or input to create one.\n","    * `error_range` - maximum allowed distance between positions in `bath` and `external_bath` for two spins to be considered the same.\n","    * `ext_r_bath` - cutoff distance from the qubit for the `external_bath`.\n","      Useful if `external_bath` has very assymetric shape and user wants to keep the precision level of the hyperfine at different distances consistent.\n","    * `imap` - instance of the `pc.InteractionMap` class, which contain tensor of bath spin interactions.\n","      If not provided, interactions between bath spins are assumed to be the same as one of point dipoles.\n","   \n","  Generates `BathArray` object with hyperfine tensors to be used in the calculation.\n","\n","* `Simulator.generate_clusters` is called if `order` and `r_dipole` are provided. It produces a `dict` object, which contains the indexes of the bath spins in the clusters.\n","\n","We implemented the following procedure to determine the clusters:\n","  \n","  * Each bath spin $i$ forms a **cluster of one**.\n","  * Bath spins $i$ and $j$ form a **cluster of two** if there is an edge between them (distance $d_{ij} \\le$ `r_dipole`).\n","  * Bath spins $i$, $j$, and $k$ form a **cluster of three** if enough edges connect them (e.g., there are two edges $ij$ and $jk$).\n","  \n","In general, we assume that spins $\\{1..n\\}$ form clusters if they form a connected graph. Only clusters up to the size indicated by the `order` parameter (equal to the **CCE order**) are included.\n"]},{"cell_type":"markdown","metadata":{"id":"mxl7F50v69ku"},"source":["We use `matplotlib` to visualize the spatial distribution of the spin bath. The grey lines show connected pairs of nuclear spins, red dashed lines show clusters of three."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"RxGZiWpG69ku"},"outputs":[],"source":["# add 3D axis\n","fig = plt.figure(figsize=(6,6))\n","ax = fig.add_subplot(projection='3d')\n","\n","# We want to visualize the smaller bath\n","data = mock.bath\n","\n","# First plot the positions of the bath\n","colors = np.abs(data.A[:,2,2]/data.A[:,2,2].max())\n","ax.scatter3D(data.x, data.y, data.z, c=colors, cmap='viridis');\n","\n","# Plot all pairs of nuclear spins, which are contained\n","# in the calc.clusters dictionary under they key 2\n","for c in mock.clusters[2]:\n","    ax.plot3D(data.x[c], data.y[c], data.z[c], color='grey')\n","\n","# Plot all triplets of nuclear spins, which are contained\n","# in the calc.clusters dictionary under they key 3\n","for c in mock.clusters[3]:\n","    ax.plot3D(data.x[c], data.y[c], data.z[c], color='red', ls='--', lw=0.5)\n","\n","ax.set(xlabel='x (A)', ylabel='y (A)', zlabel='z (A)');"]},{"cell_type":"markdown","source":["Adjust the `r_dipole`, `r_bath`, or `order` parameters and visualize the spin bath."],"metadata":{"id":"a-0Lr3o25RIh"}},{"cell_type":"markdown","source":["### **(2.2.1) *Try it yourself* : V<sub>B</sub><sup>-</sup> in hBN**"],"metadata":{"id":"y7W4dQ4xLXC5"}},{"cell_type":"markdown","source":["Generate the spin bath clusters around the V<sub>B</sub><sup>-</sup> center in hBN. Complete (1.6) if you haven't already.\n","\n","You should see more nuclear spins and more clusters in this case (**Why?**)."],"metadata":{"id":"A-LLozSjLiLF"}},{"cell_type":"code","source":["# Begin here\n"],"metadata":{"id":"LqrboT89MKTA"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### **(2.3) Setting up a CCE simulation**"],"metadata":{"id":"XOg9XM9X46N9"}},{"cell_type":"markdown","metadata":{"id":"j6-uMsL069ku"},"source":["Now let's set up the `Simulator` object for a CCE simulation."]},{"cell_type":"code","execution_count":null,"metadata":{"pycharm":{"name":"#%%\n"},"id":"F4YdXG0b69ku"},"outputs":[],"source":["# Parameters of CCE calculations engine\n","\n","# Order of CCE aproximation\n","order = 2  # singlets and pairs of nuclear spins\n","\n","# Bath cutoff radius\n","r_bath = 40  # in A\n","\n","# Cluster cutoff radius\n","r_dipole = 8  # in A\n"]},{"cell_type":"markdown","metadata":{"id":"Xkk2h1Ho69ku"},"source":["We will use the ``CenterArray`` object to store the properties of the central spin."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"Bh4hQ9Nu69ku"},"outputs":[],"source":["# position of central spin\n","position = [0, 0, 0]\n","\n","# Qubit levels (in Sz basis)\n","alpha = [0, 0, 1]; beta = [0, 1, 0]\n","\n","# ZFS Parametters of NV center in diamond\n","D = 2.88 * 1e6 # in kHz\n","E = 0 # in kHz\n","nv = pc.CenterArray(spin=1, position=position, D=D, E=E, alpha=alpha, beta=beta)"]},{"cell_type":"markdown","source":["**Important**: The position of the central spin should match the position of the NV center you created in (1.4)!"],"metadata":{"id":"Rf5gGhj-Jr7G"}},{"cell_type":"markdown","metadata":{"id":"j4sx6wNX69ku"},"source":["The code already knows the properties of the most common nuclear spins and of the elecron spin (accessible under the name `'e'`), however the user can provide their own by calling `BathArray.add_type` method. The way to initiallize `SpinType` objects is the same as in  `SpinDict` above."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"1n96Dk8l69ku"},"outputs":[],"source":["# The code already knows most exsiting isotopes.\n","#              Bath spin types\n","#              name    spin    gyro       quadrupole (for s>1/2)\n","spin_types = [('14N',  1,      1.9338,    20.44),\n","              ('13C',  1 / 2,  6.72828),\n","              ('29Si', 1 / 2, -5.3188),]\n","diamond_bath.add_type(*spin_types)"]},{"cell_type":"markdown","metadata":{"id":"MkXWeTJZ69ku"},"source":["Now we can set up the `Simulator` object."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"rhGZThbR69ku"},"outputs":[],"source":["# Setting the runner engine\n","calc = pc.Simulator(spin=nv, bath=diamond_bath,\n","                    r_bath=r_bath, r_dipole=r_dipole, order=order)"]},{"cell_type":"markdown","metadata":{"id":"mgu8YFws69ku"},"source":["Taking advantage of subclassing `np.ndarray` we can change the hyperfine and quadrupole tensors of any nuclear spin. Here, we update the quadrupole tensor of the Nitrogen nuclear spin."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"Jg480auE69ku"},"outputs":[],"source":["# Set model quadrupole tensor at N atom\n","quad = np.asarray([[-2.5, 0, 0],\n","                   [0, -2.5, 0],\n","                   [0, 0,  5.0]]) * 1e3 * 2 * np.pi\n","\n","calc.bath['Q'][calc.bath['N'] == '14N'] = quad\n"]},{"cell_type":"markdown","metadata":{"id":"0xmJXhRR69ku"},"source":["Note, that we need to apply the boolean mask **second** because of how structured arrays work."]},{"cell_type":"markdown","source":["### **(2.4) Compute the coherence function (conventional CCE)**"],"metadata":{"id":"zEC8_tQpUw_f"}},{"cell_type":"markdown","metadata":{"pycharm":{"name":"#%% md\n"},"id":"T5ZNDit_69ku"},"source":["The general interface to compute any property with PyCCE is implemented through the `Simulator.compute` method. It takes two keyword arguments to determine which quantity to compute and how:\n","\n","* `method` can take 'cce' or 'gcce' values, and determines which method to use - conventional or generalized CCE.\n","* `quantity` can take 'coherence' or 'noise' values, and determines which quantity to compute - coherence function\n","  or autocorrelation function of the noise.\n","\n","Each of the methods can be performed with Monte Carlo bath state sampling (if `nbstates` keyword is non zero) and with interlaced averaging (If `interlaced` keyword is set to `True`).\n","\n","In the first example we use the conventional CCE method without Monte Carlo bath state sampling. In the conventional CCE method the Hamiltonian is projected on the qubit levels, and the coherence is computed from the overlap of the bath evolution, entangled with two different qubit states.\n","\n","The conventional CCE requires one argument:\n","\n","* `timespace` — time points at which the coherence function is computed.\n","\n","Additionally, one can provide the following arguments now, instead of when initializing ``Simulator`` object:\n","\n","* `pulses` — number of pulses in CPMG sequence (0 - free induction decay, 1 - Hahn echo, etc., default 0) or explicit sequence of pulses as `Sequence` class instance.\n","* `magnetic_field` — magnetic field along z-axis or vector of the magnetic field. Default (0, 0, 0)."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"BzN42yP569kx"},"outputs":[],"source":["# Time points\n","time_space = np.linspace(0, 2, 1001)  # in ms\n","\n","# Number of pulses in CPMG seq (0 = FID, 1 = HE)\n","n = 1\n","\n","# Magnetic field (Bx By Bz)\n","b = np.array([0, 0, 15])  # in G\n","\n","# Compute coherence function\n","l_conv = calc.compute(time_space, pulses=n, magnetic_field=b,\n","                      method='cce', quantity='coherence', as_delay=False)"]},{"cell_type":"markdown","source":["Let's plot the coherence function."],"metadata":{"id":"0CvPeSKlSdXb"}},{"cell_type":"code","source":["fig, ax = plt.subplots(figsize=(6,6))\n","\n","ax.plot(time_space, l_conv)\n","ax.set(xlabel='Time (ms)', ylabel='Coherence')"],"metadata":{"id":"uPmSL2IiR35b"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["There are two features of interest here:\n","\n","1.   The wiggles. These oscillations are caused by the precession of <sup>13</sup>C nuclear spins in a magnetic field.\n","2.   The overall decay. The decay is caused by flip-flop dynamics of <sup>13</sup>C pairs.\n","\n","For more information, check out [this paper](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.85.115303) on decoherence dynamics in different magnetic field regimes.\n"],"metadata":{"id":"Cm2DvKMWNEEa"}},{"cell_type":"markdown","source":["### **(2.5)** Fit $T_2$ and $n$"],"metadata":{"id":"0EGJQqfOMSBe"}},{"cell_type":"markdown","source":["We can also fit the coherence function to a stretched exponential function,\n","\n","$$\n","\\mathcal{L}(t) = \\exp\\left(- \\left(\\frac{t}{T_2}\\right)^n \\right)\n","$$\n","\n","to obtain the spin-echo coherence time ($T_2$) and the stretch exponential ($n$).\n","\n","**Note**: Because we computed the coherence function by applying a single pulse `n=1`, we fit $T_2$ (the Hahn-echo coherence time)."],"metadata":{"id":"NS-yduvVSgOx"}},{"cell_type":"code","source":["from scipy.optimize import curve_fit\n","\n","# Stretched exponential fitting function\n","def Lfunc(time, T2, n):\n","  return np.exp(-(time/T2)**n)\n","\n","# Fit T2 and n\n","param, conv = curve_fit(Lfunc, time_space, l_conv)\n","\n","print(f'T2: {round(param[0], 2)} ms')\n","print(f'n: {round(param[1], 2)}')"],"metadata":{"id":"TB-ePz2xSpJV"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["We can check the goodness of our fit."],"metadata":{"id":"rWm-gos3UX4a"}},{"cell_type":"code","source":["fig, ax = plt.subplots(figsize=(6,6))\n","\n","ax.plot(time_space, l_conv, label=\"CCE\")\n","ax.plot(time_space, Lfunc(time_space, *param), label=\"Fit\")\n","\n","ax.set(xlabel='Time (ms)', ylabel='Coherence')\n","ax.legend()"],"metadata":{"id":"NVyd-sKSUcmz"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["### **(2.6) Check for convergence**"],"metadata":{"id":"DlLxvBSEVAQe"}},{"cell_type":"markdown","metadata":{"id":"JrxdSdHS69ky"},"source":["It is important to check the convergence of the CCE simulation with respect to the `order`, `r_bath`, and `r_dipole` parameters of the `Simulator` object.\n","\n","First, define all of the parameters."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"Zt-7lioA69ky"},"outputs":[],"source":["parameters = dict(\n","    order          = 2,           # CCE order\n","    r_bath         = 40,          # Size of the bath in A\n","    r_dipole       = 8,           # Cutoff of pairwise clusters in A\n","    position       = [0, 0, 0],   # Position of central Spin\n","    alpha          = [0, 0, 1],\n","    beta           = [0, 1, 0],\n","    pulses         = 1,           # N pulses in CPMG sequence\n","    magnetic_field = [0, 0, 500]\n",")\n","\n","time_space = np.linspace(0, 2, 201)  # Time points in ms"]},{"cell_type":"markdown","metadata":{"id":"AWf_maap69ky"},"source":["We can define a little helper function to streamline the process. Note that resetting the parameters automatically recomputes the properties of the bath."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"Rub7vgQS69ky"},"outputs":[],"source":["def runner(variable, values):\n","    invalue = parameters[variable]\n","    calc = pc.Simulator(spin=1, bath=diamond_bath, **parameters)\n","    lcoh = []\n","    for v in values:\n","        setattr(calc, variable, v)\n","        l = calc.compute(time_space, method='cce',\n","                         quantity='coherence')\n","        lcoh.append(l.real)\n","    parameters[variable] = invalue\n","    df_coh = pd.DataFrame(lcoh, columns=time_space, index=values).T\n","    return df_coh"]},{"cell_type":"markdown","metadata":{"id":"E3GbKxQr69ky"},"source":["Now we can compute the coherence function at different values of the parameters:"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"t0YccUdi69kz"},"outputs":[],"source":["orders = runner('order', [1, 2, 3, 4])\n","rbs = runner('r_bath', [20, 30, 40, 50, 60])\n","rds = runner('r_dipole', [4, 6, 8, 10])"]},{"cell_type":"markdown","metadata":{"id":"mzOGHvSs69kz"},"source":["We can visualize the convergence of the coherence function with respect to different parameters:"]},{"cell_type":"code","execution_count":null,"metadata":{"id":"redGQkv269kz"},"outputs":[],"source":["fig, axes = plt.subplots(1, 3, figsize=(12, 3))\n","orders.plot(ax=axes[0], title='order')\n","rbs.plot(ax=axes[1], title='r_bath')\n","rds.plot(ax=axes[2], title='r_dipole')\n","for ax in axes:\n","    ax.set(xlabel='Time (ms)', ylabel='Coherence')\n","fig.tight_layout()"]}],"metadata":{"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.8.13"},"pycharm":{"stem_cell":{"cell_type":"raw","metadata":{"collapsed":false},"source":[]}},"colab":{"provenance":[{"file_id":"https://github.com/MICCoMpy/PyCCE/blob/master/examples/tutorials/diamond_nv.ipynb","timestamp":1770045952649}],"collapsed_sections":["lG1jgJm569kq","vQSQg5fHPp_U","yW_9oYmHPv39","YPZPKUdnP2Mv","CqYxT34D6cFf","UZJDWjd_C41t","veF-8il7PaCs","-2JNxxFEDQkL","y7W4dQ4xLXC5"]}},"nbformat":4,"nbformat_minor":0}
\ No newline at end of file
diff --git a/plugins/aiida-pycce b/plugins/aiida-pycce
new file mode 160000
index 0000000..ec1874f
--- /dev/null
+++ b/plugins/aiida-pycce
@@ -0,0 +1 @@
+Subproject commit ec1874fe15d4b3f4c5e03e19d99558e947dd2786