.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/example_00_ed_calculator.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_example_00_ed_calculator.py: Exact diagonalization ===================================== Here we show how to find the eigenvalues and eigenvectors of a many-body Hamiltonian of fermions with Coulomb interactions. We then determine their spin and orbital angular momentum and how this changes when we switch on spin-orbit coupling. .. GENERATED FROM PYTHON SOURCE LINES 12-13 Import the necessary modules. .. GENERATED FROM PYTHON SOURCE LINES 13-18 .. code-block:: Python import numpy as np from math import comb import matplotlib.pyplot as plt import edrixs .. GENERATED FROM PYTHON SOURCE LINES 19-24 Parameters ------------------------------------------------------------------------------ Define the orbital angular momentum number :math:`l=1` (i.e. a `p` shell), the number of spin-orbitals, the occupancy and the Slater integrals. :math:`F^{k}` with :math:`k=0,2`: .. GENERATED FROM PYTHON SOURCE LINES 24-29 .. code-block:: Python l = 1 norb = 6 noccu = 2 F0, F2 = 4.0, 1.0 .. GENERATED FROM PYTHON SOURCE LINES 30-37 Coulomb interactions ------------------------------------------------------------------------------ The Coulomb interactions in EDRIXS are described by a tensor. Understanding this in full is complicated and requires careful consideration of the symmetry of the interactions. See :ref:`sphx_glr_auto_examples_example_09_Coulomb.py` for more discussion if desired. EDRIXS can construct the matrix via .. GENERATED FROM PYTHON SOURCE LINES 37-39 .. code-block:: Python umat = edrixs.get_umat_slater('p', F0, F2) .. GENERATED FROM PYTHON SOURCE LINES 40-50 Create basis ------------------------------------------------------------------------------ Now we build the binary form of the Fock basis :math:`|F\rangle` (we consider it preferable to use the standard :math:`F` and trust the reader to avoid confusing it with the interaction parameters.) The Fock basis is the simplest legitimate form for the basis and it consists of a series of 1s and 0s where 1 means occupied and 0 means empty. These are in order up, down, up, down, up, down. For computational efficiency, these are encoded as integers, with orbital zero in the least significant bit. We print occupations in orbital order. .. GENERATED FROM PYTHON SOURCE LINES 50-55 .. code-block:: Python basis = edrixs.get_fock_basis_int(norb, noccu) print("Integer\tOccupations (orbital 0 first)") for state in basis.basis_int: occupations = "".join(str((state >> orbital) & 1) for orbital in range(norb)) print(f"{state:2d}\t{occupations}") .. rst-class:: sphx-glr-script-out .. code-block:: none Integer Occupations (orbital 0 first) 3 110000 5 101000 9 100100 17 100010 33 100001 6 011000 10 010100 18 010010 34 010001 12 001100 20 001010 36 001001 24 000110 40 000101 48 000011 .. GENERATED FROM PYTHON SOURCE LINES 56-59 We expect the number of these states to be given by the mathematical combination of two electrons distributed among six states (three spin-orbitals with two spins per orbital). .. GENERATED FROM PYTHON SOURCE LINES 59-62 .. code-block:: Python message = ("We predict C(norb={}, noccu={})={:.0f} states and we got {:d}, " "which is reassuring!") print(message.format(norb, noccu, comb(norb, noccu), len(basis))) .. rst-class:: sphx-glr-script-out .. code-block:: none We predict C(norb=6, noccu=2)=15 states and we got 15, which is reassuring! .. GENERATED FROM PYTHON SOURCE LINES 63-65 Note that in more complicated problems with both valence and core electrons, the edrixs convention is to list the valence electrons first. .. GENERATED FROM PYTHON SOURCE LINES 67-89 Transform interactions into Fock basis ------------------------------------------------------------------------------ edrixs works by initially creating a Hamiltonian matrix :math:`\hat{H}` in the single particle basis and then transforming into our chosen Fock basis. The Hamiltonian contains a two fermion term built from the single particle matrix :math:`t_{ij}` (passed as :code:`emat`) .. math:: \hat{H}_{2} = \langle F_l|\sum_{ij}t_{ij}\hat{f}_{i}^{\dagger} \hat{f}_{j}|F_r\rangle and a four fermion term built from the interaction tensor :math:`U_{ijkl}` (passed as :code:`umat`) .. math:: \hat{H}_{4} = \langle F_l|\sum_{ijkl}U_{ijkl}\hat{f}_{i}^{\dagger} \hat{f}_{j}^{\dagger} \hat{f}_{k}\hat{f}_{l}|F_r\rangle For now we just build the :code:`umat` piece. :code:`backend='dense'` selects full matrices including the zeros. :code:`build_op` returns a callable; ed calls it to construct the matrix when needed. .. GENERATED FROM PYTHON SOURCE LINES 89-92 .. code-block:: Python backend = 'dense' H = edrixs.build_op(None, umat, basis, backend=backend) .. GENERATED FROM PYTHON SOURCE LINES 93-98 Diagonalize the matrix ------------------------------------------------------------------------------ Diagonalization returns eigenvalues :code:`e` and eigenvectors :code:`v` where eigenvalue :code:`e[i]` corresponds to eigenvector :code:`v[:,i]`. .. GENERATED FROM PYTHON SOURCE LINES 98-103 .. code-block:: Python e, v = edrixs.ed(H, num_evals=len(basis)) print("{} eigenvalues and {} eigenvectors {} elements long.".format(len(e), v.shape[1], v.shape[0])) .. rst-class:: sphx-glr-script-out .. code-block:: none 15 eigenvalues and 15 eigenvectors 15 elements long. .. GENERATED FROM PYTHON SOURCE LINES 104-111 Computing expectation values ------------------------------------------------------------------------------ To interpret the results, it is informative to compute the expectation values related to the spin :math:`\mathbf{S}`, orbital :math:`\mathbf{L}`, and total :math:`\mathbf{J}` angular momentum. We first load the relevant matrices for these quantities for a `p` atomic shell. We need to specify that we would like to include spin when loading the orbital operator. .. GENERATED FROM PYTHON SOURCE LINES 111-115 .. code-block:: Python orb_mom = edrixs.get_orb_momentum(l, ispin=True) spin_mom = edrixs.get_spin_momentum(l) tot_mom = orb_mom + spin_mom .. GENERATED FROM PYTHON SOURCE LINES 116-120 We again transform these matrices to our Fock basis to build the operators. Each is a one-body operator, so it is passed as the :code:`emat` argument with :code:`umat=None`. They have three Cartesian components, which we build one at a time. .. GENERATED FROM PYTHON SOURCE LINES 120-123 .. code-block:: Python opL, opS, opJ = [[edrixs.build_op(component, None, basis, backend=backend)() for component in op] for op in [orb_mom, spin_mom, tot_mom]] .. GENERATED FROM PYTHON SOURCE LINES 124-133 Recall that quantum mechanics forbids us from knowing all three Cartesian components of angular momentum at once, so we want to compute the squares of these operators i.e. .. math:: \mathbf{S}^2 = S^2_x + S^2_y + S^2_z\\ \mathbf{L}^2 = L^2_x + L^2_y + L^2_z\\ \mathbf{J}^2 = J^2_x + J^2_y + J^2_z .. GENERATED FROM PYTHON SOURCE LINES 133-137 .. code-block:: Python L2 = np.dot(opL[0], opL[0]) + np.dot(opL[1], opL[1]) + np.dot(opL[2], opL[2]) S2 = np.dot(opS[0], opS[0]) + np.dot(opS[1], opS[1]) + np.dot(opS[2], opS[2]) J2 = np.dot(opJ[0], opJ[0]) + np.dot(opJ[1], opJ[1]) + np.dot(opJ[2], opJ[2]) .. GENERATED FROM PYTHON SOURCE LINES 138-141 Remember that the eigenvalues of :math:`\mathbf{S}^2` are in the form :math:`S(S+1)` etc. and that they can be obtained by calculating the projection of the operators onto our eigenvectors. .. GENERATED FROM PYTHON SOURCE LINES 141-144 .. code-block:: Python L2_val = edrixs.cb_op(L2, v).diagonal().real S2_val = edrixs.cb_op(S2, v).diagonal().real J2_val = edrixs.cb_op(J2, v).diagonal().real .. GENERATED FROM PYTHON SOURCE LINES 145-147 We can determine the degeneracy of the eigenvalues numerically and print out the values as follows .. GENERATED FROM PYTHON SOURCE LINES 147-155 .. code-block:: Python e = np.round(e, decimals=6) degeneracy = [sum(eval == e) for eval in e] header = "{:<3s}\t{:>8s}\t{:>8s}\t{:>8s}\t{:>8s}" print(header.format("# ", "E ", "S(S+1)", "L(L+1)", "Degen.")) for i, eigenvalue in enumerate(e): values_list = [i, eigenvalue, S2_val[i], L2_val[i], degeneracy[i]] print("{:<3d}\t{:8.3f}\t{:8.3f}\t{:8.3f}\t{:>3d}".format(*values_list)) .. rst-class:: sphx-glr-script-out .. code-block:: none # E S(S+1) L(L+1) Degen. 0 3.800 2.000 2.000 9 1 3.800 2.000 2.000 9 2 3.800 2.000 2.000 9 3 3.800 2.000 2.000 9 4 3.800 2.000 2.000 9 5 3.800 2.000 2.000 9 6 3.800 2.000 2.000 9 7 3.800 2.000 2.000 9 8 3.800 2.000 2.000 9 9 4.040 0.000 6.000 5 10 4.040 0.000 6.000 5 11 4.040 0.000 6.000 5 12 4.040 0.000 6.000 5 13 4.040 0.000 6.000 5 14 4.400 0.000 -0.000 1 .. GENERATED FROM PYTHON SOURCE LINES 156-165 We see :math:`S=0` and :math:`S=1` states coming from the two combinations of the spin 1/2 particles. :math:`L` can take values of 0, 1, 2. Remember that spin states have degeneracy of :math:`2S+1` and the same is true for orbital states. We must multiply these :math:`S` and :math:`L` degeneracies to get the total degeneracy. Since these particles are fermions, the overall state must be antisymmetric, which dictates the allowed combinations of :math:`S` and :math:`L`. .. GENERATED FROM PYTHON SOURCE LINES 167-170 Energy level diagram ------------------------------------------------------------------------------ Let us show our findings graphically .. GENERATED FROM PYTHON SOURCE LINES 170-190 .. code-block:: Python fig, ax = plt.subplots() for i, eigenvalue in enumerate(np.unique(e)): art = ax.plot([0, 1], [eigenvalue, eigenvalue], '-', color='C{}'.format(i)) ind = np.where(eigenvalue == e)[0][0] L = (-1 + np.sqrt(1 + 4*L2_val[ind]))/2 S = (-1 + np.sqrt(1 + 4*S2_val[ind]))/2 message = "L={:.0f}, S={:.0f} ({:.0f})" ax.text(1, eigenvalue, message.format(L, S, degeneracy[ind]), horizontalalignment='right', verticalalignment='bottom', color='C{}'.format(i)) ax.set_ylabel('Energy') for loc in ['right', 'top', 'bottom']: ax.spines[loc].set_visible(False) ax.yaxis.set_ticks_position('left') ax.set_xticks([]) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_example_00_ed_calculator_001.png :alt: example 00 ed calculator :srcset: /auto_examples/images/sphx_glr_example_00_ed_calculator_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 191-195 We see Hund's rules in action! Rule 1 says that the highest spin :math:`S=1` state has the lowest energy. Of the two :math:`S=0` states, the one with the larger orbital angular momentum (here :math:`L=2`) is lower in energy, following rule 2. .. GENERATED FROM PYTHON SOURCE LINES 197-204 Spin orbit coupling ------------------------------------------------------------------------------ For fun, we can see how this changes when we add spin orbit coupling (SOC). This is a two-fermion operator that we create, transform into the Fock basis and add to the prior Hamiltonian. To make things easy, let us make the SOC small so that the LS coupling approximation is valid and we can still track the states. .. GENERATED FROM PYTHON SOURCE LINES 204-207 .. code-block:: Python soc = edrixs.atom_hsoc('p', 0.1) H2 = edrixs.build_op(soc, umat, basis, backend=backend) .. GENERATED FROM PYTHON SOURCE LINES 208-209 Then, we redo the diagonalization and print the results. .. GENERATED FROM PYTHON SOURCE LINES 209-223 .. code-block:: Python e2, v2 = edrixs.ed(H2, num_evals=len(basis)) e2 = np.round(e2, decimals=6) degeneracy2 = [sum(eval == e2) for eval in e2] print() message = "With SOC\n {:<3s}\t{:>8s}\t{:>8s}\t{:>8s}\t{:>8s}\t{:>8s}" print(message.format("#", "E", "S(S+1)", "L(L+1)", "J(J+1)", "degen.")) J2_val_soc = edrixs.cb_op(J2, v2).diagonal().real L2_val_soc = edrixs.cb_op(L2, v2).diagonal().real S2_val_soc = edrixs.cb_op(S2, v2).diagonal().real for i, eigenvalue in enumerate(e2): values_list = [i, eigenvalue, S2_val_soc[i], L2_val_soc[i], J2_val_soc[i], degeneracy2[i]] print("{:<3d}\t{:8.3f}\t{:8.3f}\t{:8.3f}\t{:8.3f}\t{:8.3f}".format(*values_list)) .. rst-class:: sphx-glr-script-out .. code-block:: none With SOC # E S(S+1) L(L+1) J(J+1) degen. 0 3.673 1.927 1.927 -0.000 1.000 1 3.750 2.000 2.000 2.000 3.000 2 3.750 2.000 2.000 2.000 3.000 3 3.750 2.000 2.000 2.000 3.000 4 3.827 1.802 2.396 6.000 5.000 5 3.827 1.802 2.396 6.000 5.000 6 3.827 1.802 2.396 6.000 5.000 7 3.827 1.802 2.396 6.000 5.000 8 3.827 1.802 2.396 6.000 5.000 9 4.063 0.198 5.604 6.000 5.000 10 4.063 0.198 5.604 6.000 5.000 11 4.063 0.198 5.604 6.000 5.000 12 4.063 0.198 5.604 6.000 5.000 13 4.063 0.198 5.604 6.000 5.000 14 4.427 0.073 0.073 -0.000 1.000 .. GENERATED FROM PYTHON SOURCE LINES 224-225 and we make an equivalent energy level diagram. .. GENERATED FROM PYTHON SOURCE LINES 225-245 .. code-block:: Python fig, ax = plt.subplots() for i, eigenvalue in enumerate(np.unique(e2)): art = ax.plot([0, 1], [eigenvalue, eigenvalue], '-', color='C{}'.format(i)) ind = np.where(eigenvalue == e2)[0][0] J = (-1 + np.sqrt(1+4*J2_val_soc[ind]))/2 message = "J={:.0f} ({:.0f})" ax.text(1, eigenvalue, message.format(J, degeneracy2[ind]), horizontalalignment='right', verticalalignment='bottom', color='C{}'.format(i)) ax.set_ylabel('Energy') for loc in ['right', 'top', 'bottom']: ax.spines[loc].set_visible(False) ax.yaxis.set_ticks_position('left') ax.set_xticks([]) plt.show() .. image-sg:: /auto_examples/images/sphx_glr_example_00_ed_calculator_002.png :alt: example 00 ed calculator :srcset: /auto_examples/images/sphx_glr_example_00_ed_calculator_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 246-250 It is clear that we have split the :math:`S=1` state, which branches into three states from :math:`J=|L-S|, |L-S|+1, ..., |L+S|`. Since the shell is less than half full, Hund's third rule dictates that the smaller :math:`J` states have the lower energies. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.797 seconds) .. _sphx_glr_download_auto_examples_example_00_ed_calculator.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_00_ed_calculator.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_00_ed_calculator.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_00_ed_calculator.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_