$$\newcommand{\ket}[1]{\left|{#1}\right\rangle}$$
$$\newcommand{\bra}[1]{\left\langle{#1}\right|}$$


# Explanation of First-Generation QM Model

# Disclaimer

Ultimately, I would like to understand, and be able to explain, the quantum mechanics behind the entire process of simulating an NMR spectrum. For now, here is a "recipe" of the steps to arrive at the spin Hamiltonian, and how its eigensolution can be used to calculate frequencies and intensities.

Two sources in particular enabled this:
1. Materials by Ilya Kuprov at SpinDynamics.org, particularly [Module I, Lecture 5](http://spindynamics.org/Spin-Dynamics---Part-I---Lecture-05.php) and the Matlab code of [Module II, Lecture 05](http://spindynamics.org/Spin-Dynamics---Part-II---Lecture-05.php) and [06](http://spindynamics.org/Spin-Dynamics---Part-II---Lecture-06.php).
2. [Materials](http://www.users.csbsju.edu/~frioux/workinprogress.html#Spectroscopy) by Frank Rioux at St. John's University and College of St. Benedict. In particular, [*ABC Proton NMR Using Tensor Algebra*](http://www.users.csbsju.edu/~frioux/nmr/ABC-NMR-Tensor.pdf) was very helpful.


```python
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
```


<style>.container { width:100% !important; }</style>



```python
import numpy as np
from math import sqrt
from scipy.linalg import eigh
from scipy.sparse import kron, csc_matrix, csr_matrix, lil_matrix, bmat

import bokeh.io
import bokeh.plotting

```

## Constructing the Hamiltonian From Scratch

Start with the Pauli matrices:
\begin{align}
\sigma_x = \begin{pmatrix}0& \frac{1}{2}\\ \frac{1}{2}&0\end{pmatrix}, 
\sigma_y = \begin{pmatrix}0& -\frac{i}{2}\\ \frac{i}{2}&0\end{pmatrix}, 
\sigma_z = \begin{pmatrix}\frac{1}{2}& 0\\ 0&-\frac{1}{2}\end{pmatrix}
\end{align}

plus the identity matrix $I = \begin{pmatrix}1&0\\0&1\end{pmatrix}$


```python
sigma_x = np.matrix([[0, 1 / 2], [1 / 2, 0]])
sigma_y = np.matrix([[0, -1j / 2], [1j / 2, 0]])
sigma_z = np.matrix([[1 / 2, 0], [0, -1 / 2]])
unit = np.matrix([[1, 0], [0, 1]])
```

The required inputs are a list of frequencies $\nu_i$ and a matrix of $J_{ij}$ coupling constants:


```python
freqlist = [10.0, 20.0]
couplings = np.array([[0, 5], [5, 0]])

```

Let's break the original hamiltonian function down into smaller functions and explain their roles:


```python
def hamiltonian(freqlist, couplings):
    """Calculate the Hamiltonian for a spin system (isotropic liquid).
    
    Arguments:
        freqlist: a list of n chemical shifts (in Hz)
        couplings: an n x n array of J coupling constants (in Hz)
    Return:
        H: numpy.ndarray spin Hamiltonian
    """
    Lx, Ly, Lz = create_krons(freqlist)
    Lproduct = cartesian_products(Lx, Ly, Lz)
    H_zeeman = hamiltonian_diagonal(freqlist, Lz)
    H_J = hamiltonian_off_diagonal(couplings, Lproduct)
    H = H_zeeman + H_J
    return H
```

### Step 1: Each spin gets its own $L_x$, $L_y$ and $L_z$ operators.

These are formed from Kronecker products between $\sigma_{x/y/z}$ and $I$ operators.

Each individual product, for n spins, uses 1$\sigma_{x/y/z}$ and (n - 1) $I$ operators. They all differ in where in the sequence the $\sigma_{x/y/z}$ operator is placed.

For 3 spins, and using $L_z$ for example:

\begin{align}
L_{z_1} &= \sigma_z \otimes I \otimes I\\
L_{z_2} &= I \otimes \sigma_z \otimes I\\
L_{z_3} &= I \otimes I \otimes \sigma_z
\end{align}


In the Python code, these individual $L_{x/y/z_n}$ operators get stored in a \[0, n\] array, Ln. 

*{I'm not sure why a 2-D array of one row was used. It's possible that this could be simplified to a regular 1-D array, but it's also possible that there was an issue with some other operation, such as creating the Cartesian products. Need to check this.}*



```python
def create_krons(freqlist):
    nspins = len(freqlist)
    # The following empty arrays will be used to store the
    # Cartesian spin operators.
    Lx = np.empty((1, nspins), dtype='object')
    Ly = np.empty((1, nspins), dtype='object')
    Lz = np.empty((1, nspins), dtype='object')

    for n in range(nspins):
        Lx[0, n] = 1
        Ly[0, n] = 1
        Lz[0, n] = 1
        for k in range(nspins):
            if k == n:                                  # Diagonal element
                Lx[0, n] = np.kron(Lx[0, n], sigma_x)
                Ly[0, n] = np.kron(Ly[0, n], sigma_y)
                Lz[0, n] = np.kron(Lz[0, n], sigma_z)
            else:                                       # Off-diagonal element
                Lx[0, n] = np.kron(Lx[0, n], unit)
                Ly[0, n] = np.kron(Ly[0, n], unit)
                Lz[0, n] = np.kron(Lz[0, n], unit)
                
    return Lx, Ly, Lz
```

### Step 2: Create the sums of cartesian products of $L$ operators.

Eventually, the off-diagonal components of the Hamiltonian  $H$ require calculating Cartesian products of the $L$ operators. In an attempt to hopefully "vectorize" these for faster computation, all of these products were calculated at once. Then, when a particular result is required (e.g. $L_{x_1}L_{x_2}+L_{y_1}L_{y_2}+L_{z_1}L_{z_2}$) it can be plucked like a bonbon out of a tray when desired.

First, the $L_x, L_y, $ and $L_z$ operators were compiled into a 3 x n array of operators:
```python
Lcol = np.vstack((Lx, Ly, Lz)).real
```
created:
\begin{align}
L_{col} = \begin{pmatrix}
L_{x_1}& L_{x_2}&\dots & L_{x_n}\\ 
L_{y_1}& L_{y_2}&\dots & L_{y_n}\\
L_{z_1}& L_{z_2}&\dots & L_{z_n}
\end{pmatrix}
\end{align}

Its transpose created the n x 3 array of operators:

\begin{align}
L_{row} = \begin{pmatrix}
L_{x_1}& L_{y_1}&L_{z_1}\\ 
L_{x_2}& L_{y_2}&L_{z_2}\\
\vdots&\vdots&\vdots\\
L_{x_n}& L_{y_n}&L_{z_n}
\end{pmatrix}
\end{align}

The product of these two arrays gives an array of the Cartesian products:

\begin{align}
L_{product}&= L_{row} \cdot L_{col} \\
&=\Tiny\begin{pmatrix}
L_{x_1}L_{x_1}+L_{y_1}L_{y_1}+L_{z_1}L_{z_1}&L_{x_1}L_{x_2}+L_{y_1}L_{y_2}+L_{z_1}L_{z_2}&\dots&L_{x_1}L_{x_n}+L_{y_1}L_{y_n}+L_{z_1}L_{z_n}\\
L_{x_2}L_{x_1}+L_{y_2}L_{y_1}+L_{z_2}L_{z_1}&L_{x_2}L_{x_2}+L_{y_2}L_{y_2}+L_{z_2}L_{z_2}&\dots&L_{x_2}L_{x_n}+L_{y_2}L_{y_n}+L_{z_2}L_{z_n}\\
\vdots& &\ddots& \\
L_{x_n}L_{x_1}+L_{y_n}L_{y_1}+L_{z_n}L_{z_1}&L_{x_n}L_{x_2}+L_{y_n}L_{y_2}+L_{z_n}L_{z_2}&\dots&L_{x_n}L_{x_n}+L_{y_n}L_{y_n}+L_{z_n}L_{z_n}\\
\end{pmatrix}
\end{align}






```python
def cartesian_products(Lx, Ly, Lz):
    Lcol = np.vstack((Lx, Ly, Lz)).real
    Lrow = Lcol.T  # As opposed to sparse version of code, this works!
    Lproduct = np.dot(Lrow, Lcol)
    return Lproduct
```

### Step 3: Add the Zeeman (on-diagonal) terms to the Hamiltonian.

\begin{align}
H_{Zeeman} = \sum_{i=1}^n \nu_i L_{z_i}
\end{align}


```python
def hamiltonian_diagonal(freqlist, Lz):
    nspins = len(freqlist)
    # Hamiltonian operator
    H = np.zeros((2**nspins, 2**nspins))

    # Add Zeeman interactions:
    for n in range(nspins):
        H = H + freqlist[n] * Lz[0, n]
    
    return H


```

### Step 4: Add the J-coupling (off-diagonal) terms to the Hamiltonian.

\begin{align}
H_J &= \sum_{i=1}^n \sum_{j=1}^n \frac{J_{ij}}{2} (L_{x_i}L_{x_j}+L_{y_i}L_{y_j}+L_{z_i}L_{z_j})\\
H &= H_{Zeeman} + H_J
\end{align}

In an attempt to vectorize the calculation for better performance, this was accomplished by doing in-place multiplication of the matrix of $J$ coupling constants and the Cartesian products $L_{product}$.


```python
def hamiltonian_off_diagonal(couplings, Lproduct):
    
    nspins = len(couplings[0])
    # Hamiltonian operator
    H = np.zeros((2**nspins, 2**nspins))
    
    # Scalar couplings

    # Testing with MATLAB discovered J must be /2.
    # Believe it is related to the fact that in the SpinDynamics.org simulation
    # freqs are *2pi, but Js by pi only.
    scalars = 0.5 * couplings
    scalars = np.multiply(scalars, Lproduct)
    for n in range(nspins):
        for k in range(nspins):
            H += scalars[n, k].real

    return H
```

## Extracting Signal Frequencies and Intensities From the Hamiltonian

To simulate a "modern" NMR experiment, a 90° pulse and FID acquisition is simulated, followed by Fourier transform. This is the approach used in Kuprov's Matlab code, and should be the required approach for any experiment requiring a more elaborate pulse sequence. 

For a simple NMR spectrum, we can adipt a "continuous wave spectrometer" approach. We can find the resonance frequencies and their relative intensities directly from the spin Hamiltonian. The time-independent Schrodinger equation $H\Psi = E\Psi$ is solved for eigenvectors and corresponding eigenvalues. 

For each $\psi_i$, the eigenvectors are the coefficients $c_n$ for each pure spin state. For a two-spin system, for example,

$\psi_i = c_1\ket{\alpha\alpha} + c_2\ket{\alpha\beta} + c_3\ket{\beta\alpha} + c_4\ket{\beta\beta}$.

and the corresponding eigenvector would be \begin{bmatrix}c_1\\c_2\\c_3\\c_4\end{bmatrix}

For a one-spin system, the two states for "spin-up" ($\ \ket{\uparrow}$ or $\ket{\alpha}$) and for "spin-down" ($\ \ket{\downarrow}$ or $\ket{\beta}$) are represented by vectors $\begin{bmatrix}1\\0\end{bmatrix}$ and $\begin{bmatrix}0\\1\end{bmatrix}$, respectively. For "pure" multiple-spin states, their vectors are obtained by taking tensor products of these vectors. For example:

\begin{align}
\ket{\alpha\alpha} &=
\begin{bmatrix}1\\0\end{bmatrix}\otimes\begin{bmatrix}1\\0\end{bmatrix} = \begin{bmatrix}1\\0\\0\\0\end{bmatrix}
\\
\ket{\alpha\beta} &= \begin{bmatrix}1\\0\end{bmatrix}\otimes\begin{bmatrix}0\\1\end{bmatrix} = \begin{bmatrix}0\\1\\0\\0\end{bmatrix}
\\
\ket{\beta\alpha} &= \begin{bmatrix}0\\1\end{bmatrix}\otimes\begin{bmatrix}1\\0\end{bmatrix} = \begin{bmatrix}0\\0\\1\\0\end{bmatrix}
\\
\ket{\beta\beta} &= \begin{bmatrix}0\\1\end{bmatrix}\otimes\begin{bmatrix}0\\1\end{bmatrix} = \begin{bmatrix}0\\0\\0\\1\end{bmatrix}
\end{align}

A (coincidental?) consequence of this is that the index for $H$, expressed in binary form as a series of 0s and 1s, is the eigenvector for the associated pure spin state (cf. Rioux's *ABC Proton NMR Using Tensor Algebra*). 

Since allowed transitions change the total spin of a system by $\pm$ 1, this is analogous to transitions only being allowed between spin states whose binary indices only differ at one bit. In computing terms, if the Hamming weight of the two indices differ by exactly 1, the transition is allowed. 


```python
def popcount(n=0):
    """
    Computes the popcount (binary Hamming weight) of integer n
    input:
        :param n: an integer
    returns:
        popcount of integer (binary Hamming weight)

    """
    return bin(n).count('1')


def is_allowed(m=0, n=0):
    """
    determines if a transition between two spin states is allowed or forbidden.
    The transition is allowed if one and only one spin (i.e. bit) changes
    input: integers whose binary codes for a spin state
        :param n:
        :param m:
    output: 1 = allowed, 0 = forbidden

    """
    return popcount(m ^ n) == 1
```

Knowing this, we can create a transition probability matrix $T$, where $T_{ij} = 1$ if a transition between states $i$ and $j$ are allowed, and $0$ if not.


```python
def transition_matrix(n):
    """
    Creates a matrix of allowed transitions.
    The integers 0-n, in their binary form, code for a spin state (alpha/beta).
    The (i,j) cells in the matrix indicate whether a transition from spin state
    i to spin state j is allowed or forbidden.
    See the is_allowed function for more information.

    input:
        :param n: size of the n,n matrix (i.e. number of possible spin states)

    :returns: a transition matrix that can be used to compute the intensity of
    allowed transitions.
    """
    # function was optimized by only calculating upper triangle and then adding
    # the lower.
    T = lil_matrix((n, n))  # sparse matrix created
    for i in range(n - 1):
        for j in range(i + 1, n):
            if is_allowed(i, j):
                T[i, j] = 1
    T = T + T.T
    return T
```

The eigenvector solutions for the Hamiltonian include two pure states ("all-up/$\alpha$" and "all-down/$\beta$", plus mixed states. We can construct a matrix $V_{col}$ where each column of the matrix is an eigenvector solution, in their indexed order:

\begin{align}
V_{col} = \begin{pmatrix}
\ket{\psi_1} &\ket{\psi_2} &\dots &\ket{\psi_n}
\end{pmatrix}
=\begin{pmatrix}
\begin{bmatrix}c_1\\c_2\\\vdots\\c_n\end{bmatrix}_1
&\begin{bmatrix}c_1\\c_2\\\vdots\\c_n\end{bmatrix}_2
&\dots
&\begin{bmatrix}c_1\\c_2\\\vdots\\c_n\end{bmatrix}_n
\end{pmatrix}
\end{align}

and where its transpose $V_{row} = V_{col}^T$ has an eigenvector for each row:

\begin{align}
V_{row}=\begin{pmatrix}
\bra{\psi_1} \\
\bra{\psi_2} \\
\vdots\\
\bra{\psi_n} \\
\end{pmatrix}
=\begin{pmatrix}
\begin{bmatrix}c_1&c_2&\dots&c_n\end{bmatrix}_1\\
\begin{bmatrix}c_1&c_2&\dots&c_n\end{bmatrix}_2\\
\vdots\\
\begin{bmatrix}c_1&c_2&\dots&c_n\end{bmatrix}_n\\
\end{pmatrix}
\end{align}


The intensity matrix $I$ can be obtained by taking $V_{row}\cdot T \cdot V_{col}$ and squaring it element-wise, so that $I_{ij}$ is the relative probability of a transition between the $\psi_i$ and $\psi_j$ states. The difference in energy between the two states gives the frequency in Hz.


```python
def simsignals(H, nspins):
    """
    Solves the spin Hamiltonian H and returns a list of (frequency, intensity)
    tuples. Nuclei must be spin-1/2.
    Inputs:
        :param H: a Hamiltonian array
        :param nspins: number of nuclei
    :return spectrum: a list of (frequency, intensity) tuples.
    """
    # This routine was optimized for speed by vectorizing the intensity
    # calculations, replacing a nested-for signal-by-signal calculation.
    # Considering that hamiltonian was dramatically faster when refactored to
    # use arrays instead of sparse matrices, consider an array refactor to this
    # function as well.

    # The eigensolution calculation apparently must be done on a dense matrix,
    # because eig functions on sparse matrices can't return all answers?!
    # Using eigh so that answers have only real components and no residual small
    # unreal components b/c of rounding errors
    E, V = np.linalg.eigh(H)    # V will be eigenvectors, v will be frequencies
    print(E)
    # Eigh still leaves residual 0j terms, so:
    V = np.asmatrix(V.real)
    print(V)

    # Calculate signal intensities
    Vcol = csc_matrix(V)
    print('Vcol:')
    print(Vcol.todense())
    Vrow = csr_matrix(Vcol.T)
    print('Vrow:')
    print(Vrow.todense())
    m = 2 ** nspins
    T = transition_matrix(m)
    print('T:')
    print(T.todense())
    print('T•Vcol:')
    mid_t = T * Vcol
    print(mid_t.todense())
    I = Vrow * T * Vcol
    print('I:')
    print(I.todense())
    print('I squared:')
    I = np.square(I.todense())
    print(I)

    spectrum = []
    for i in range(m - 1):
        for j in range(i + 1, m):
            if I[i, j] > 0.01:  # consider making this minimum intensity
                                # cutoff a function arg, for flexibility
                v = abs(E[i] - E[j])
                spectrum.append((v, I[i, j]))

    return spectrum
```


```python
def hamiltonian_with_prints(freqlist, couplings):
    # The following empty arrays will be used to store the
    # Cartesian spin operators.
    Lx = np.empty((1, nspins), dtype='object')
    Ly = np.empty((1, nspins), dtype='object')
    Lz = np.empty((1, nspins), dtype='object')

    for n in range(nspins):
        Lx[0, n] = 1
        Ly[0, n] = 1
        Lz[0, n] = 1
        for k in range(nspins):
            if k == n:                                  # Diagonal element
                Lx[0, n] = np.kron(Lx[0, n], sigma_x)
                Ly[0, n] = np.kron(Ly[0, n], sigma_y)
                Lz[0, n] = np.kron(Lz[0, n], sigma_z)
            else:                                       # Off-diagonal element
                Lx[0, n] = np.kron(Lx[0, n], unit)
                Ly[0, n] = np.kron(Ly[0, n], unit)
                Lz[0, n] = np.kron(Lz[0, n], unit)

    Lcol = np.vstack((Lx, Ly, Lz)).real
    Lrow = Lcol.T  # As opposed to sparse version of code, this works!
    Lproduct = np.dot(Lrow, Lcol)
    print(Lcol)
    print('-'*10)
    print(Lrow)
    print('-'*10)
    print(Lproduct)

    # Hamiltonian operator
    H = np.zeros((2**nspins, 2**nspins))

    # Add Zeeman interactions:
    for n in range(nspins):
        H = H + freqlist[n] * Lz[0, n]

    # Scalar couplings

    # Testing with MATLAB discovered J must be /2.
    # Believe it is related to the fact that in the SpinDynamics.org simulation
    # freqs are *2pi, but Js by pi only.
    scalars = 0.5 * couplings
    scalars = np.multiply(scalars, Lproduct)
    for n in range(nspins):
        for k in range(nspins):
            H += scalars[n, k].real
            
    print('Lz: ', Lz)

    return H
    
```


```python
H = hamiltonian(freqlist, couplings)
H
```




    matrix([[ 16.25,   0.  ,   0.  ,   0.  ],
            [  0.  ,  -6.25,   2.5 ,   0.  ],
            [  0.  ,   2.5 ,   3.75,   0.  ],
            [  0.  ,   0.  ,   0.  , -13.75]])




```python
np.set_printoptions(threshold=np.nan)
nspins = len(freqlist)
simsignals(H, nspins)
```

    [-13.75        -6.84016994   4.34016994  16.25      ]
    [[ 0.          0.         -0.          1.        ]
     [ 0.         -0.97324899 -0.22975292  0.        ]
     [ 0.          0.22975292 -0.97324899  0.        ]
     [ 1.          0.         -0.          0.        ]]
    Vcol:
    [[ 0.          0.          0.          1.        ]
     [ 0.         -0.97324899 -0.22975292  0.        ]
     [ 0.          0.22975292 -0.97324899  0.        ]
     [ 1.          0.          0.          0.        ]]
    Vrow:
    [[ 0.          0.          0.          1.        ]
     [ 0.         -0.97324899  0.22975292  0.        ]
     [ 0.         -0.22975292 -0.97324899  0.        ]
     [ 1.          0.          0.          0.        ]]
    T:
    [[ 0.  1.  1.  0.]
     [ 1.  0.  0.  1.]
     [ 1.  0.  0.  1.]
     [ 0.  1.  1.  0.]]
    T•Vcol:
    [[ 0.         -0.74349607 -1.20300191  0.        ]
     [ 1.          0.          0.          1.        ]
     [ 1.          0.          0.          1.        ]
     [ 0.         -0.74349607 -1.20300191  0.        ]]
    I:
    [[ 0.         -0.74349607 -1.20300191  0.        ]
     [-0.74349607  0.          0.         -0.74349607]
     [-1.20300191  0.          0.         -1.20300191]
     [ 0.         -0.74349607 -1.20300191  0.        ]]
    I squared:
    [[ 0.         0.5527864  1.4472136  0.       ]
     [ 0.5527864  0.         0.         0.5527864]
     [ 1.4472136  0.         0.         1.4472136]
     [ 0.         0.5527864  1.4472136  0.       ]]





    [(6.9098300562505255, 0.55278640450004191),
     (18.090169943749473, 1.4472135954999577),
     (23.090169943749473, 0.55278640450004191),
     (11.909830056250525, 1.4472135954999577)]



What happens if we only have one nucleus?


```python
freqlist = [10.0]
couplings = np.array([[0]])
nspins = len(freqlist)
```

two uncoupled spins:


```python
freqlist = [10.0, 20.0]
couplings = np.array([[0, 0], [0, 0]])
nspins = len(freqlist)
```


```python
def spin2():
    v = np.array([150-7.5, 150+7.5])
    J = lil_matrix((2, 2))
    J[0, 1] = 12
    J = J + J.T
    return v, J
```


```python
v, J = spin2()
J = J.todense()
print(v)
print(J)
H = hamiltonian_with_prints(v, J)
H
```

    [ 142.5  157.5]
    [[  0.  12.]
     [ 12.   0.]]
    [[ matrix([[ 0. ,  0. ,  0.5,  0. ],
            [ 0. ,  0. ,  0. ,  0.5],
            [ 0.5,  0. ,  0. ,  0. ],
            [ 0. ,  0.5,  0. ,  0. ]])
      matrix([[ 0. ,  0.5,  0. ,  0. ],
            [ 0.5,  0. ,  0. ,  0. ],
            [ 0. ,  0. ,  0. ,  0.5],
            [ 0. ,  0. ,  0.5,  0. ]])]
     [ matrix([[ 0.+0.j ,  0.+0.j ,  0.-0.5j,  0.+0.j ],
            [ 0.+0.j ,  0.+0.j ,  0.+0.j ,  0.-0.5j],
            [ 0.+0.5j,  0.+0.j ,  0.+0.j ,  0.+0.j ],
            [ 0.+0.j ,  0.+0.5j,  0.+0.j ,  0.+0.j ]])
      matrix([[ 0.+0.j ,  0.-0.5j,  0.+0.j ,  0.-0.j ],
            [ 0.+0.5j,  0.+0.j ,  0.+0.j ,  0.+0.j ],
            [ 0.+0.j ,  0.-0.j ,  0.+0.j ,  0.-0.5j],
            [ 0.+0.j ,  0.+0.j ,  0.+0.5j,  0.+0.j ]])]
     [ matrix([[ 0.5,  0. ,  0. ,  0. ],
            [ 0. ,  0.5,  0. ,  0. ],
            [ 0. ,  0. , -0.5, -0. ],
            [ 0. ,  0. , -0. , -0.5]])
      matrix([[ 0.5,  0. ,  0. ,  0. ],
            [ 0. , -0.5,  0. , -0. ],
            [ 0. ,  0. ,  0.5,  0. ],
            [ 0. , -0. ,  0. , -0.5]])]]
    ----------
    [[ matrix([[ 0. ,  0. ,  0.5,  0. ],
            [ 0. ,  0. ,  0. ,  0.5],
            [ 0.5,  0. ,  0. ,  0. ],
            [ 0. ,  0.5,  0. ,  0. ]])
      matrix([[ 0.+0.j ,  0.+0.j ,  0.-0.5j,  0.+0.j ],
            [ 0.+0.j ,  0.+0.j ,  0.+0.j ,  0.-0.5j],
            [ 0.+0.5j,  0.+0.j ,  0.+0.j ,  0.+0.j ],
            [ 0.+0.j ,  0.+0.5j,  0.+0.j ,  0.+0.j ]])
      matrix([[ 0.5,  0. ,  0. ,  0. ],
            [ 0. ,  0.5,  0. ,  0. ],
            [ 0. ,  0. , -0.5, -0. ],
            [ 0. ,  0. , -0. , -0.5]])]
     [ matrix([[ 0. ,  0.5,  0. ,  0. ],
            [ 0.5,  0. ,  0. ,  0. ],
            [ 0. ,  0. ,  0. ,  0.5],
            [ 0. ,  0. ,  0.5,  0. ]])
      matrix([[ 0.+0.j ,  0.-0.5j,  0.+0.j ,  0.-0.j ],
            [ 0.+0.5j,  0.+0.j ,  0.+0.j ,  0.+0.j ],
            [ 0.+0.j ,  0.-0.j ,  0.+0.j ,  0.-0.5j],
            [ 0.+0.j ,  0.+0.j ,  0.+0.5j,  0.+0.j ]])
      matrix([[ 0.5,  0. ,  0. ,  0. ],
            [ 0. , -0.5,  0. , -0. ],
            [ 0. ,  0. ,  0.5,  0. ],
            [ 0. , -0. ,  0. , -0.5]])]]
    ----------
    [[ matrix([[ 0.75+0.j,  0.00+0.j,  0.00+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.75+0.j,  0.00+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.00+0.j,  0.75+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.00+0.j,  0.00+0.j,  0.75+0.j]])
      matrix([[ 0.25+0.j,  0.00+0.j,  0.00+0.j,  0.00+0.j],
            [ 0.00+0.j, -0.25+0.j,  0.50+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.50+0.j, -0.25+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.00+0.j,  0.00+0.j,  0.25+0.j]])]
     [ matrix([[ 0.25+0.j,  0.00+0.j,  0.00+0.j,  0.00+0.j],
            [ 0.00+0.j, -0.25+0.j,  0.50+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.50+0.j, -0.25+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.00+0.j,  0.00+0.j,  0.25+0.j]])
      matrix([[ 0.75+0.j,  0.00+0.j,  0.00+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.75+0.j,  0.00+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.00+0.j,  0.75+0.j,  0.00+0.j],
            [ 0.00+0.j,  0.00+0.j,  0.00+0.j,  0.75+0.j]])]]
    Lz:  [[ matrix([[ 0.5,  0. ,  0. ,  0. ],
            [ 0. ,  0.5,  0. ,  0. ],
            [ 0. ,  0. , -0.5, -0. ],
            [ 0. ,  0. , -0. , -0.5]])
      matrix([[ 0.5,  0. ,  0. ,  0. ],
            [ 0. , -0.5,  0. , -0. ],
            [ 0. ,  0. ,  0.5,  0. ],
            [ 0. , -0. ,  0. , -0.5]])]]





    matrix([[ 153. ,    0. ,    0. ,    0. ],
            [   0. ,  -10.5,    6. ,    0. ],
            [   0. ,    6. ,    4.5,    0. ],
            [   0. ,    0. ,    0. , -147. ]])




```python

```

01286bb4d67e115a4ad93d3042a9cc9ef5aaa16a
ipynb
Jupyter Notebook
jupyter/.ipynb_checkpoints/qm_explanation-checkpoint.ipynb
sametz/nmrtools
dd36bd21fca1169e371bfdaff82ab721235e1d7c
MIT
jupyter/.ipynb_checkpoints/qm_explanation-checkpoint.ipynb
sametz/nmrtools
dd36bd21fca1169e371bfdaff82ab721235e1d7c
MIT
jupyter/.ipynb_checkpoints/qm_explanation-checkpoint.ipynb
sametz/nmrtools
dd36bd21fca1169e371bfdaff82ab721235e1d7c
MIT
Qwen/Qwen-72B
1. YES
2. YES