.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/example_01_crystal_field.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_01_crystal_field.py: Crystal fields ===================================== This example explains how to implement crystal fields in edrixs. .. GENERATED FROM PYTHON SOURCE LINES 8-9 We need to import these modules. .. GENERATED FROM PYTHON SOURCE LINES 9-16 .. code-block:: Python import edrixs import numpy as np import scipy from edrixs.fock_basis import FockBasisSpec np.set_printoptions(precision=2, suppress=True, linewidth=90) .. GENERATED FROM PYTHON SOURCE LINES 17-23 Crystal field matrices ------------------------------------------------------------------------------ Let us start by considering the common case of a :math:`d` atomic shell in a cubic crystal field. This is controlled by parameter :math:`10D_q` and is described in terms of a matrix which we will assign to :code:`cfmat`. edrixs can make this matrix via .. GENERATED FROM PYTHON SOURCE LINES 23-26 .. code-block:: Python ten_dq = 10 cfmat = edrixs.angular_momentum.cf_cubic_d(ten_dq) .. GENERATED FROM PYTHON SOURCE LINES 27-34 Note this matrix is in a complex harmonic basis :math:`Y^m_l` where :math:`m` goes from :math:`-l, -l+1, ..., l-1, l`. There is an up spin and a down spin for each :math:`Y^m_l`. This matrix is not diagonal in the complex harmonic basis, but it would be diagonal in the real harmonic basis :math:`d_{3z^2-r^2}, d_{xz}, d_{yz}, d_{x^2-y^2}, d_{xy}`. Let us diagonalize this matrix as a check and print out the energies and their degeneracies. .. GENERATED FROM PYTHON SOURCE LINES 34-45 .. code-block:: Python evals_i, evecs_i = np.linalg.eigh(cfmat) evals_i = evals_i.round(decimals=6) unique_e = np.unique(evals_i) degeneracies = [sum(evalue == evals_i) for evalue in unique_e] print("E \tDegeneracy") for evalue, degenvalue in zip(unique_e, degeneracies): print("{:.1f}\t{:.0f}".format(evalue, degenvalue)) print("{} distinct energies".format(len(unique_e))) .. rst-class:: sphx-glr-script-out .. code-block:: none E Degeneracy -4.0 6 6.0 4 2 distinct energies .. GENERATED FROM PYTHON SOURCE LINES 46-48 This makes sense! We see two different energies split by :math:`10D_q=10`. Let us look at the six columns corresponding to the lower energy eigenvalues. .. GENERATED FROM PYTHON SOURCE LINES 48-51 .. code-block:: Python print(evecs_i[:, :6].real) .. rst-class:: sphx-glr-script-out .. code-block:: none [[ 0. 0. 0. 0. -0.71 -0. ] [ 0. 0. 0. 0. 0. -0.71] [ 0. 0. -1. 0. 0. 0. ] [ 0. 0. 0. -1. 0. 0. ] [ 0. 0. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0. 0. ] [ 1. 0. 0. 0. 0. 0. ] [ 0. 1. 0. 0. 0. 0. ] [ 0. 0. 0. 0. 0.71 0. ] [ 0. 0. 0. 0. -0. 0.71]] .. GENERATED FROM PYTHON SOURCE LINES 52-55 These are the set of so-called :math:`t_{2g}` orbitals, composed of :math:`Y^2_2, Y^{-2}_2, Y^{1}_2, Y^{-1}_2`. The rest of the eigenvectors (the last four) are .. GENERATED FROM PYTHON SOURCE LINES 55-57 .. code-block:: Python print(evecs_i[:, 6:].real) .. rst-class:: sphx-glr-script-out .. code-block:: none [[-0.71 0. 0. 0. ] [ 0. -0.71 0. 0. ] [ 0. 0. 0. 0. ] [ 0. 0. 0. 0. ] [ 0. 0. 1. 0. ] [ 0. 0. 0. 1. ] [ 0. 0. 0. 0. ] [ 0. 0. 0. 0. ] [-0.71 0. 0. 0. ] [ 0. -0.71 0. 0. ]] .. GENERATED FROM PYTHON SOURCE LINES 58-70 These are the set of so-called :math:`e_{g}` orbitals, composed of :math:`Y^2_2, Y^{-2}_2, Y^{0}_2`. We can use edrixs to prove that :code:`cfmat` would be diagonal in the real harmonic basis. An operator :math:`\hat{O}` can be transformed into an operator in another basis :math:`\hat{O}^{\prime}` using a unitary transformation matrix :math:`T` as .. math:: \hat{O}^{\prime} = (T)^{\dagger} \hat{O} (T). This is computed as follows .. GENERATED FROM PYTHON SOURCE LINES 70-73 .. code-block:: Python cfmat_rhb = edrixs.cb_op(cfmat, edrixs.tmat_c2r('d', ispin=True)) print(cfmat_rhb.real) .. rst-class:: sphx-glr-script-out .. code-block:: none [[ 6. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 6. 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. -4. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. -4. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. -4. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. -4. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 6. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 6. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. -4. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0. 0. -4.]] .. GENERATED FROM PYTHON SOURCE LINES 74-79 where :code:`edrixs.tmat_c2r('d', ispin=True)` is the transformation matrix. We needed to tell edrixs that we are working with a :math:`d`-shell and that it should include spin. We could also have transformed :code:`evecs_i` to see how these eigenvectors are composed of the real harmonic basis. We will see an example of this later. .. GENERATED FROM PYTHON SOURCE LINES 81-85 Crystal field on an atom ------------------------------------------------------------------------------ To simulate the solid state, we need to combine the crystal field with Coulomb interactions. Let us choose an atomic model for Ni. .. GENERATED FROM PYTHON SOURCE LINES 85-89 .. code-block:: Python norb = 10 noccu = 8 atomic_slater = edrixs.get_atom_data('Ni', '3d', noccu, edge='L3')['slater_i'] .. GENERATED FROM PYTHON SOURCE LINES 90-93 Let us implement a tetragonal crystal field, for which we need to pass :code:`d1` the splitting of :math:`d_{yz}/d_{zx}` and :math:`d_{xy}` and :code:`d3` the splitting of :math:`d_{3z^2-r^2}` and :math:`d_{x^2-y^2}`. .. GENERATED FROM PYTHON SOURCE LINES 93-95 .. code-block:: Python ten_dq, d1, d3 = 2.5, 0.9, .2 .. GENERATED FROM PYTHON SOURCE LINES 96-107 Since this is a purely atomic calculation with no core shell, we assemble the single-particle ingredients directly rather than going through :func:`~edrixs.models.model_1v1c`. :class:`~edrixs.fock_basis.FockBasisSpec` describes the :math:`d^8` Fock space (8 electrons in 10 spin-orbitals), :func:`~edrixs.angular_momentum.cf_tetragonal_d` builds the one-body crystal-field matrix and :func:`~edrixs.coulomb_utensor.get_umat_slater` the Coulomb tensor, both in the default complex-harmonic basis. :func:`~edrixs.solvers.build_op` then constructs the many-body Hamiltonian and :func:`~edrixs.solvers.ed` diagonalizes it. We put this procedure into a function, with the option to scale the Coulomb interactions. .. GENERATED FROM PYTHON SOURCE LINES 107-120 .. code-block:: Python def diagonalize(scaleU=1): basis_i = FockBasisSpec.from_args(10, 8) slater_i = [term[1] * scaleU for term in atomic_slater[:3]] cfmat = edrixs.angular_momentum.cf_tetragonal_d(ten_dq=ten_dq, d1=d1, d3=d3) umat_i = edrixs.get_umat_slater('d', *slater_i) hmat_i = edrixs.build_op(cfmat, umat_i, basis_i, backend='dense') evals_i, evecs_i = edrixs.ed(hmat_i, num_evals=len(basis_i)) evals_i = evals_i - np.min(evals_i) # define ground state as zero energy return evals_i, evecs_i, basis_i .. GENERATED FROM PYTHON SOURCE LINES 121-126 Let us look what happens when we run the function with the Coulomb interactions switched off and check the degeneracy of the output. Look at this python `string formatting tutorial `_ if the code is confusing. .. GENERATED FROM PYTHON SOURCE LINES 126-136 .. code-block:: Python evals_i, evecs_i, basis_i = diagonalize(scaleU=0) evals_i = evals_i.round(decimals=6) unique_e = np.unique(evals_i) degeneracies = [sum(evalue == evals_i) for evalue in unique_e] print("E \tDegeneracy") for evalue, degenvalue in zip(unique_e, degeneracies): print("{:.1f}\t{:.0f}".format(evalue, degenvalue)) print("{} distinct energies".format(len(unique_e))) .. rst-class:: sphx-glr-script-out .. code-block:: none E Degeneracy 0.0 1 0.2 4 0.4 1 2.5 4 2.7 4 3.4 8 3.6 8 5.0 1 5.9 8 6.8 6 10 distinct energies .. GENERATED FROM PYTHON SOURCE LINES 137-151 We see 10 distinct energies, which is the number of ways one can arrange two holes among 4 energy levels -- which makes sense as the tetragonal field involves four levels :math:`zx/zy, xy, 3z^2-r^2, x^2-y^2`. To see what is going on in more detail, we can also calculate the expectation values of the occupancy number of the orbitals :math:`3z^2-r^2, zx, zy, x^2-y^2, xy`. To create the operator, first write the matrix in the real harmonics basis :math:`|3z^2-r^2,\uparrow>`, :math:`|3z^2-r^2,\downarrow>`, :math:`|zx,\uparrow>`, :math:`|zx,\downarrow>`, etc. In this basis, they take a simple form: only the diagonal terms have element 1. We therefore make a 3D empty array and assign the diagonal as 1. Check out the `NumPy indexing notes `_ if needed. .. GENERATED FROM PYTHON SOURCE LINES 151-155 .. code-block:: Python nd_real_harmoic_basis = np.zeros((norb, norb, norb), dtype=complex) indx = np.arange(norb) nd_real_harmoic_basis[indx, indx, indx] = 1 .. GENERATED FROM PYTHON SOURCE LINES 156-158 Recalling the necessity to put everything in the same basis, we transform into the complex harmonic basis and then transform into our Fock basis .. GENERATED FROM PYTHON SOURCE LINES 158-173 .. code-block:: Python nd_complex_harmoic_basis = edrixs.cb_op(nd_real_harmoic_basis, edrixs.tmat_r2c('d', True)) def transform_occupancy_operators(evecs_i): op = [edrixs.cb_op( edrixs.build_op(nd_vec, None, basis_i, backend='dense')(), evecs_i) for nd_vec in nd_complex_harmoic_basis] return op nd_op = transform_occupancy_operators(evecs_i) .. GENERATED FROM PYTHON SOURCE LINES 174-178 We apply the operator and print out as follows. Check the `numpy docs `_ if the details of how the spin pairs have been added up is not immediately transparent. .. GENERATED FROM PYTHON SOURCE LINES 178-188 .. code-block:: Python nd_expt = np.array([nd_vec.diagonal().real for nd_vec in nd_op]) message = "{:>3s}" + "\t{:>6s}"*5 print(message.format(*"E 3z2-r2 zx zy x2-y2 xy".split(" "))) message = "{:>3.1f}" + "\t{:>6.1f}"*5 for evalue, row in zip(evals_i, nd_expt.T): spin_pairs = row.reshape(-1, 2).sum(1) print(message.format(evalue, *spin_pairs)) .. rst-class:: sphx-glr-script-out .. code-block:: none E 3z2-r2 zx zy x2-y2 xy 0.0 2.0 2.0 2.0 0.0 2.0 0.2 1.0 2.0 2.0 1.0 2.0 0.2 1.0 2.0 2.0 1.0 2.0 0.2 1.0 2.0 2.0 1.0 2.0 0.2 1.0 2.0 2.0 1.0 2.0 0.4 0.0 2.0 2.0 2.0 2.0 2.5 2.0 2.0 2.0 1.0 1.0 2.5 2.0 2.0 2.0 1.0 1.0 2.5 2.0 2.0 2.0 1.0 1.0 2.5 2.0 2.0 2.0 1.0 1.0 2.7 1.0 2.0 2.0 2.0 1.0 2.7 1.0 2.0 2.0 2.0 1.0 2.7 1.0 2.0 2.0 2.0 1.0 2.7 1.0 2.0 2.0 2.0 1.0 3.4 2.0 1.5 1.5 1.0 2.0 3.4 2.0 1.5 1.5 1.0 2.0 3.4 2.0 1.5 1.5 1.0 2.0 3.4 2.0 1.5 1.5 1.0 2.0 3.4 2.0 1.5 1.5 1.0 2.0 3.4 2.0 1.5 1.5 1.0 2.0 3.4 2.0 1.5 1.5 1.0 2.0 3.4 2.0 1.5 1.5 1.0 2.0 3.6 1.0 1.8 1.2 2.0 2.0 3.6 1.0 1.5 1.5 2.0 2.0 3.6 1.0 1.7 1.3 2.0 2.0 3.6 1.0 1.7 1.3 2.0 2.0 3.6 1.0 1.2 1.8 2.0 2.0 3.6 1.0 1.5 1.5 2.0 2.0 3.6 1.0 1.5 1.5 2.0 2.0 3.6 1.0 1.1 1.9 2.0 2.0 5.0 2.0 2.0 2.0 2.0 0.0 5.9 2.0 1.5 1.5 2.0 1.0 5.9 2.0 1.5 1.5 2.0 1.0 5.9 2.0 1.5 1.5 2.0 1.0 5.9 2.0 1.5 1.5 2.0 1.0 5.9 2.0 1.5 1.5 2.0 1.0 5.9 2.0 1.5 1.5 2.0 1.0 5.9 2.0 1.5 1.5 2.0 1.0 5.9 2.0 1.5 1.5 2.0 1.0 6.8 2.0 0.8 1.2 2.0 2.0 6.8 2.0 1.0 1.0 2.0 2.0 6.8 2.0 1.0 1.0 2.0 2.0 6.8 2.0 1.0 1.0 2.0 2.0 6.8 2.0 1.2 0.8 2.0 2.0 6.8 2.0 1.0 1.0 2.0 2.0 .. GENERATED FROM PYTHON SOURCE LINES 189-193 The lowest energy state involves putting both holes in the :math:`x^2-y^2` orbital, which makes sense. Now, let us redo the procedure including Coulomb repulsion, which imposes an energy cost to putting multiple electrons in the same orbital. .. GENERATED FROM PYTHON SOURCE LINES 193-207 .. code-block:: Python evals_i, evecs_i, _ = diagonalize(scaleU=1) nd_op = transform_occupancy_operators(evecs_i) nd_expt = np.array([nd_vec.diagonal().real for nd_vec in nd_op]) message = "{:>3s}" + "\t{:>6s}"*5 print(message.format(*"E 3z2-r2 zx zy x2-y2 xy".split(" "))) message = "{:>3.1f}" + "\t{:>6.1f}"*5 for evalue, row in zip(evals_i, nd_expt.T): spin_pairs = row.reshape(-1, 2).sum(1) print(message.format(evalue, *spin_pairs)) .. rst-class:: sphx-glr-script-out .. code-block:: none E 3z2-r2 zx zy x2-y2 xy 0.0 1.0 2.0 2.0 1.0 2.0 0.0 1.0 2.0 2.0 1.0 2.0 0.0 1.0 2.0 2.0 1.0 2.0 2.4 1.2 2.0 2.0 0.8 2.0 2.5 1.0 2.0 2.0 1.0 2.0 2.5 1.0 2.0 2.0 2.0 1.0 2.5 1.0 2.0 2.0 2.0 1.0 2.5 1.0 2.0 2.0 2.0 1.0 3.2 1.8 1.5 1.5 1.2 2.0 3.2 1.8 1.3 1.7 1.2 2.0 3.2 1.8 1.7 1.3 1.2 2.0 3.2 1.8 1.4 1.6 1.2 2.0 3.2 1.8 1.5 1.5 1.2 2.0 3.2 1.8 1.5 1.5 1.2 2.0 4.0 2.0 1.9 1.9 1.1 1.1 4.0 2.0 1.9 1.9 1.1 1.1 4.0 2.0 1.9 1.9 1.1 1.1 4.3 0.9 2.0 2.0 1.2 1.9 4.7 1.4 1.6 1.4 1.9 1.7 4.7 1.4 1.4 1.6 1.9 1.7 4.7 1.4 1.5 1.5 1.9 1.7 4.7 1.4 1.6 1.4 1.9 1.7 4.7 1.4 1.4 1.6 1.9 1.7 4.7 1.4 1.6 1.4 1.9 1.7 4.9 1.0 2.0 2.0 2.0 1.0 5.5 2.0 2.0 2.0 1.0 1.0 5.6 1.8 1.5 1.5 1.2 2.0 5.6 1.8 1.5 1.5 1.2 2.0 6.5 1.2 1.5 1.5 1.8 2.0 6.5 1.2 1.5 1.5 1.8 2.0 6.8 1.8 1.5 1.5 1.9 1.3 6.8 1.8 1.5 1.5 1.9 1.3 6.8 1.8 1.5 1.5 1.9 1.3 6.8 1.8 1.5 1.5 1.9 1.3 6.8 1.8 1.5 1.5 1.9 1.3 6.8 1.8 1.4 1.6 1.9 1.3 7.4 2.0 1.1 1.1 1.9 1.9 7.4 2.0 1.1 1.1 1.9 1.9 7.4 2.0 1.1 1.1 1.9 1.9 8.0 2.0 1.8 1.8 2.0 0.4 8.5 2.0 1.5 1.5 2.0 1.0 8.5 2.0 1.5 1.5 2.0 1.0 9.3 2.0 1.0 1.0 2.0 2.0 9.4 2.0 1.0 1.0 2.0 2.0 12.8 1.9 1.2 1.2 1.9 1.7 .. GENERATED FROM PYTHON SOURCE LINES 208-213 Now the lowest energy state involves splitting the holes between the two highest energy :math:`x^2-y^2` and :math:`3z^2-r^2` orbitals. i.e. we have gone from low-spin to high spin. Working out the balance between these two states is often one of the first things one wants to determine upon the discovery of an interesting new material [1]_. .. GENERATED FROM PYTHON SOURCE LINES 215-218 .. rubric:: Footnotes .. [1] M. Rossi et al., `arXiv:2011.00595 (2021) `_. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.041 seconds) .. _sphx_glr_download_auto_examples_example_01_crystal_field.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: example_01_crystal_field.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: example_01_crystal_field.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: example_01_crystal_field.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_