{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"nbsphinx": "hidden"
},
"outputs": [],
"source": [
"import matplotlib\n",
"matplotlib.use(\"Pdf\")\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"# set matplotlib parameters\n",
"params = {'backend': 'ps',\n",
" 'axes.labelsize': 13,\n",
" 'font.size': 13,\n",
" 'legend.fontsize': 13,\n",
" 'xtick.labelsize': 13,\n",
" 'ytick.labelsize': 13,\n",
" 'text.usetex': False,\n",
" 'text.latex.preamble': \"\\\\usepackage{mathpazo}, \\\\usepackage{amsmath}\"}\n",
"plt.rcParams.update(params)\n",
"\n",
"import warnings \n",
"warnings.filterwarnings(\"ignore\") #ignore some matplotlib warnings\n",
"\n",
"# numpy\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Multiplet analysis & particle number histograms\n",
"\n",
"Utilizing the cthyb functionality to measure the [many body density matrix](https://triqs.github.io/cthyb/unstable/guide/static_observables_notebook.html), it is possible to analyze the multiplet structure of the impurity states, e.g. occupation probabilities of eigenstates, find the associated particle number, and obtaining the spin state of those eigenstates. In principle all observables can be examined that are good quantum numbers of the system. More information of the measurement of the density matrix in cthyb can be found here: [cthyb measure static observables](https://triqs.github.io/cthyb/unstable/guide/static_observables_notebook.html).\n",
"\n",
"This tutorial will demonstrate, based on the above mentioned [tutorial](https://triqs.github.io/cthyb/unstable/guide/static_observables_notebook.html), how to analyze the multiplet structure from a measured density matrix and the `h_loc_diag` object of the solver. For more information on how to measure those properties please take a look first [here](https://triqs.github.io/cthyb/unstable/guide/static_observables_notebook.html). This tutorial will not give the reader a full introduction to the physical background of the discussed matter, it is more designed for someone who is already familiar with the matter itself, seeking for a guide on how to perform such analysis with triqs.\n",
"\n",
"Let us consider a single-orbital Anderson impurity problem. At first, we import all necessary modules, define some input parameters and construct a cthyb Solver object in the usual way:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting run with 1 MPI rank(s) at : 2020-09-18 07:57:58.642240\n"
]
}
],
"source": [
"from h5 import HDFArchive\n",
"from triqs.gf.descriptors import Fourier\n",
"from triqs.gf import Gf, MeshImFreq, iOmega_n, inverse, GfImTime, BlockGf, Wilson\n",
"from triqs.operators import c, c_dag, n\n",
"\n",
"# redirect cthyb c++ output to notebook\n",
"from triqs.utility.redirect import start_redirect\n",
"start_redirect()\n",
"\n",
"from triqs_cthyb.solver import Solver"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Parameters\n",
"D = 1.0 # Half-bandwidth of the bath\n",
"V = 0.4 # Hybridisation amplitude\n",
"U = 4.0 # Coulomb interaction\n",
"e_f = -U/2 # Local energy level\n",
"beta = 1 # Inverse temperature\n",
"\n",
"# Construct the impurity solver with the inverse temperature\n",
"# and the structure of the Green's functions\n",
"S = Solver(beta = beta, gf_struct = [ ('up',[0]), ('down',[0]) ])\n",
"\n",
"# Initialize the non-interacting Green's function S.G0_iw\n",
"S.G0_iw['up'] << inverse(iOmega_n - e_f - V**2 * Wilson(D))\n",
"S.G0_iw['down'] << inverse(iOmega_n - e_f - V**2 * Wilson(D))\n",
"\n",
"# setting up a local interaction Hamiltonian\n",
"h_int = U * n('up',0) * n('down',0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we instruct the ``solve`` function to accumulate the density matrix by passing the solver parameters ``measure_density_matrix = True`` and ``use_norm_as_weight = True``. See also [cthyb solver parameters](https://triqs.github.io/cthyb/latest/reference/solve_parameters.html) for more information.\n",
"\n",
"Now we run the solver:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"╔╦╗╦═╗╦╔═╗ ╔═╗ ┌─┐┌┬┐┬ ┬┬ ┬┌┐ \n",
" ║ ╠╦╝║║═╬╗╚═╗ │ │ ├─┤└┬┘├┴┐\n",
" ╩ ╩╚═╩╚═╝╚╚═╝ └─┘ ┴ ┴ ┴ ┴ └─┘\n",
"\n",
"The local Hamiltonian of the problem:\n",
"-2*c_dag('down',0)*c('down',0) + -2*c_dag('up',0)*c('up',0) + 4*c_dag('down',0)*c_dag('up',0)*c('up',0)*c('down',0)\n",
"Using autopartition algorithm to partition the local Hilbert space\n",
"Found 4 subspaces.\n",
"\n",
"Warming up ...\n",
"07:57:58 12% ETA 00:00:00 cycle 1279 of 10000\n",
"\n",
"\n",
"\n",
"Accumulating ...\n",
"07:57:59 0% ETA 00:00:30 cycle 1620 of 500000\n",
"07:58:01 6% ETA 00:00:31 cycle 31153 of 500000\n",
"07:58:04 14% ETA 00:00:27 cycle 72956 of 500000\n",
"07:58:07 24% ETA 00:00:23 cycle 124185 of 500000\n",
"07:58:11 37% ETA 00:00:19 cycle 189226 of 500000\n",
"07:58:16 54% ETA 00:00:14 cycle 270958 of 500000\n",
"07:58:22 73% ETA 00:00:08 cycle 365660 of 500000\n",
"07:58:30 98% ETA 00:00:00 cycle 490624 of 500000\n",
"\n",
"\n",
"[Rank 0] Collect results: Waiting for all mpi-threads to finish accumulating...\n",
"[Rank 0] Timings for all measures:\n",
"Measure | seconds \n",
"Average sign | 0.0188487 \n",
"Density Matrix for local static observable | 0.171204 \n",
"G_tau measure | 0.0382508 \n",
"Total measure time | 0.228303 \n",
"[Rank 0] Acceptance rate for all moves:\n",
"Move set Insert two operators: 0.0767329\n",
" Move Insert Delta_up: 0.0767082\n",
" Move Insert Delta_down: 0.0767575\n",
"Move set Remove two operators: 0.0767529\n",
" Move Remove Delta_up: 0.076729\n",
" Move Remove Delta_down: 0.0767768\n",
"Move set Insert four operators: 0.0018322\n",
" Move Insert Delta_up_up: 0.000226071\n",
" Move Insert Delta_up_down: 0.00346066\n",
" Move Insert Delta_down_up: 0.00341485\n",
" Move Insert Delta_down_down: 0.000228812\n",
"Move set Remove four operators: 0.00182204\n",
" Move Remove Delta_up_up: 0.00023108\n",
" Move Remove Delta_up_down: 0.00338986\n",
" Move Remove Delta_down_up: 0.00344724\n",
" Move Remove Delta_down_down: 0.000219039\n",
"Move Shift one operator: 0.130805\n",
"[Rank 0] Warmup lasted: 0.654128 seconds [00:00:00]\n",
"[Rank 0] Simulation lasted: 31.2514 seconds [00:00:31]\n",
"[Rank 0] Number of measures: 500000\n",
"Total number of measures: 500000\n",
"Average sign: 1\n"
]
}
],
"source": [
"# Run the solver\n",
"S.solve(h_int = h_int, # interaction Hamiltonian\n",
" n_cycles = 500000, # Number of QMC cycles\n",
" length_cycle = 200, # Length of one cycle\n",
" n_warmup_cycles = 10000, # Warmup cycles\n",
" measure_density_matrix = True, # Measure the reduced density matrix\n",
" use_norm_as_weight = True) # Required to measure the density matrix"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Analysis of results\n",
"\n",
"Next we do the actual analysis of the multiplet structure. First we need to setup the operators, that we want to measure:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from triqs.operators.util import make_operator_real\n",
"from triqs.operators.util.observables import S_op, S2_op\n",
"from triqs.atom_diag import quantum_number_eigenvalues"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, we are going to take a look at $S_z, \\ S^2,$ and $N$. The two central objects we need for the analysis are:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# accumulated many body density matrix\n",
"rho = S.density_matrix\n",
"# the atom diag object containing information about the H_loc entering the solver\n",
"h_loc_diag = S.h_loc_diagonalization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`rho` is the many body density matrix, which stores the occupation probability for each eigenstate of the impurity Hamiltonian, and `h_loc_diag` contains information about the impurity Hamiltonian, e.g. eigenstates, eigenvalues, Fock states ... \n",
"\n",
"To evaluate the expectation values of $S_z, \\ S^2,$ and $N$ values we first create the operators:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# first get fundamental operators of H_loc from atom_diag object\n",
"occ_operators = [n(*op) for op in h_loc_diag.fops]\n",
"\n",
"# construct total occupation operator from list\n",
"N_op = sum(occ_operators)\n",
"\n",
"# create S2 and Sz operator\n",
"spin_names = ['up', 'down']\n",
"orb_names = [0]\n",
"\n",
"# The spin operators are created from orbital and spin names\n",
"# the flag off_diag determines if blocks of Gf are down_0 / up_0 ,\n",
"# or down / up \n",
"S2 = S2_op(spin_names, orb_names, off_diag=True)\n",
"Sz=S_op('z', spin_names, orb_names, off_diag=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we evaluate the operators in the impurity eigenstates by using the [triqs.atom_diag](https://triqs.github.io/triqs/latest/documentation/manual/triqs/atom_diag/contents.html) function [quantum_number_eigenvalues](https://triqs.github.io/triqs/latest/documentation/cpp_api/triqs/atom_diag/quantum_number_eigenvalues.html). The object will be sorted by the subspaces contained in h_loc_diag:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"# get eigenvalues of S2 and Sz\n",
"S2_states = quantum_number_eigenvalues(S2, h_loc_diag)\n",
"Sz_states = quantum_number_eigenvalues(Sz, h_loc_diag)\n",
"\n",
"# get particle numbers from h_loc_diag\n",
"particle_numbers = quantum_number_eigenvalues(N_op, h_loc_diag)\n",
"N_max = int(max(map(max, particle_numbers)))"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[array([1.]), array([1.]), array([0.]), array([2.])]\n"
]
}
],
"source": [
"print(particle_numbers)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One can see, that `particle_numbers` contains now 4 numpy arrays, representing the 4 subspaces of the impurity problem. Those are actually identified by different particle numbers $0,1,$ and $2$. Each array has the size of number of eigenstates in that particular subspace.\n",
"\n",
"Now we construct the eigenstates for each eigenvalue of `h_loc_diag` in their corresponding subspace, which can be also labeled by the quantum numbers that we just extracted. Finally, `rho` will contain the actual probability that this eigenstate is occupied:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"subspace No.: 0\n",
"EV No: 0 , N= 1 , E=0.00000 , prob:0.43381 , S(S+1): 0.75 , ms: 0.50 ,\n",
" eigenstate: +1.0000|01>\n",
"----\n",
"subspace No.: 1\n",
"EV No: 0 , N= 1 , E=0.00000 , prob:0.43393 , S(S+1): 0.75 , ms: -0.50 ,\n",
" eigenstate: +1.0000|10>\n",
"----\n",
"subspace No.: 2\n",
"EV No: 0 , N= 0 , E=2.00000 , prob:0.06625 , S(S+1): 0.00 , ms: 0.00 ,\n",
" eigenstate: +1.0000|00>\n",
"----\n",
"subspace No.: 3\n",
"EV No: 0 , N= 2 , E=2.00000 , prob:0.06601 , S(S+1): 0.00 , ms: 0.00 ,\n",
" eigenstate: +1.0000|11>\n",
"----\n",
"multiplet occupation by particle number: \n",
" [0.06624796 0.86774326 0.06600878]\n",
"sum of multiplet occ: 1.00000\n"
]
}
],
"source": [
"multiplet_occ = np.zeros(N_max + 1)\n",
"\n",
"# we loop through all subspaces, which are just indexed as integers\n",
"for sub in range(0,h_loc_diag.n_subspaces): \n",
"\n",
" print(\"subspace No.:\", sub)\n",
"\n",
" # first get Fock space spanning the subspace\n",
" fs_states = []\n",
" for ind, fs in enumerate(h_loc_diag.fock_states[sub]):\n",
" # get state in binary representation\n",
" state = bin(int(fs))[2:].rjust(N_max, '0')\n",
" fs_states.append(\"|\"+state+\">\")\n",
"\n",
" # now we loop over all eigenstates within the subspace\n",
" for ind in range(h_loc_diag.get_subspace_dim(sub)):\n",
"\n",
" # get particle number\n",
" # carefully here to not cast to int as the particle number\n",
" # can be something like 1.999999996 and would get then 1! \n",
" particle_number = int(round(particle_numbers[sub][ind]))\n",
" # here you can see that the object obtained from quantum_number_eigenvalues\n",
" # is structured in the same way as h_loc_diag\n",
" \n",
" # get spin operator values in the same way\n",
" ms=Sz_states[sub][ind]\n",
" s_square=S2_states[sub][ind]\n",
" \n",
" # energy of state\n",
" eng=h_loc_diag.energies[sub][ind]\n",
"\n",
" # probability that this state is occupied is obtained from rho\n",
" prob = rho[sub][ind,ind]\n",
" \n",
" # add to multiplet histogram sorted by particle number\n",
" multiplet_occ[particle_number] += prob\n",
" \n",
" # construct eigenvector in Fock state basis by rotation \n",
" ev_state = ''\n",
" for i, elem in enumerate(h_loc_diag.unitary_matrices[sub][:,ind]):\n",
" ev_state += ' {:+1.4f}'.format(elem)+fs_states[i]\n",
"\n",
" # print info obtained\n",
" print(\"EV No:\", ind,\n",
" \", N=\",particle_number,\n",
" \", E=\"+\"{:1.5f}\".format(eng) ,\n",
" \", prob:\"+\"{:1.5f}\".format(prob),\n",
" \", S(S+1): \"+\"{:1.2f}\".format(s_square),\n",
" \", ms: \"+\"{:1.2f}\".format(ms),\n",
" \",\\n eigenstate:\", ev_state)\n",
" \n",
" print('----')\n",
"\n",
"print('multiplet occupation by particle number: \\n', multiplet_occ)\n",
"print('sum of multiplet occ: '+\"{:1.5f}\".format(np.sum(multiplet_occ)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It can be observed that the main contribution is coming from eigenstates with $N=1$ with $m_s=\\pm 0.5$ as expected, and that the Fock states are the eigenvectors of the impurity system. The order of the orbitals in the Fock states is given by the order of operators in `h_loc_diag.fops` starting from the right site, with the first half representing one spin-channel and the second half representing the other spin-channel. Note, that in general each eigenstate can be a linear combination of Fock states. This can be important in sub-spaces with dimension $>1$, when manually using the information stored in `rho`.\n",
"\n",
"The occupations can now be ploted has histogram."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### plot particle number histogram"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, (ax1) = plt.subplots(1,1,figsize=(6,4))\n",
"\n",
"ax1.bar(range(0,len(multiplet_occ)) ,multiplet_occ, width=0.2)\n",
"\n",
"\n",
"ax1.set_xlabel(r'N')\n",
"ax1.set_ylabel(r'prob amplitude')\n",
"\n",
"ax1.set_xticks(list(range(0,len(multiplet_occ)))) \n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It can be observed, that due to the high temperature of $\\beta=1$ the two excited states are slightly populated, but the $N=1$ states have the highest occupation probability."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## More complex example\n",
"\n",
"The previous example was quite simple. We also provide a function doing the analysis with a few extra checks, storing the result as [Pandas DataFrame](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html): "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"pd.set_option('display.width', 130)\n",
"\n",
"from triqs_cthyb.multiplet_tools import multiplet_analysis"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" Computes operator expectation values from measured\n",
" density matrix and h_loc_diag object from cthyb solver.\n",
"\n",
" Measures N, Sz, S(S+1), constructs eigenstates in Fock state\n",
" basis, assigns the according probabilities from rho for all\n",
" impurity eigenstates, and stores the data in a Panda DataFrame.\n",
" For more information check the guide in the documentation of cthyb\n",
" regarding `Multiplet analysis & particle number histograms`.\n",
"\n",
" Warning: function assumes that N, Sz, and S(S+1) are good\n",
" quantum numbers, labeling the eigenstates of the system.\n",
"\n",
" Parameters\n",
" ----------\n",
" rho : numpy array\n",
" measured density matrix from cthyb or equivalent solver, structured in subspaces\n",
" h_loc_diag: triqs atom_diag object\n",
" contains information about the local Hamiltonian (Hloc_0 + H_int)\n",
" orb_names : list of int\n",
" list of orbital indices\n",
" spin_names : list of string\n",
" list of strings containing the spin channel names\n",
" off_diag: boolean\n",
" determines whether blocks of Gf are named up_0 (false) or just up (true)\n",
"\n",
" Returns:\n",
" --------\n",
" res : Panda DataFrame\n",
" containing all results structured\n",
" \n"
]
}
],
"source": [
"print(multiplet_analysis.__doc__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can use this function to analyze a stored result from a realistic five orbital calculation with an occupation of ~8.5 electrons:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"with HDFArchive('density_matrix_example.h5','r') as ar:\n",
" rho_5orb = ar['rho']\n",
" h_loc_diag_5orb = ar['h_loc_diag']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calling the function does the whole analysis at once:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"res_5orb = multiplet_analysis(rho = rho_5orb, \n",
" h_loc_diag = h_loc_diag_5orb, \n",
" # five orbitals\n",
" orb_names = [0, 1, 2, 3, 4],\n",
" # two spin channels\n",
" spin_names = ['up','down'],\n",
" # and no off-diag elements have been used in the solver\n",
" off_diag = False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"print the 15 eigenstates with highest occupation probability"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Sub# EV# N energy prob S2 m_s |m_s| state\n",
"3 1 1 9 0.144618 0.185854 0.75 0.5 0.5 +1.0000|1111111101> +0.0000|1111111110>\n",
"1 0 1 9 0.144618 0.185387 0.75 -0.5 0.5 +1.0000|1110111111> +0.0000|1111011111>\n",
"12 4 0 8 0.048852 0.063419 2.00 1.0 1.0 +1.0000|1111111100>\n",
"4 2 0 8 0.048852 0.063355 2.00 -0.0 0.0 -0.0000|0111101111> -0.0000|1011110111> +0.00...\n",
"11 3 0 8 0.048852 0.063293 2.00 -1.0 1.0 +1.0000|1110011111>\n",
"55 20 0 10 6.036870 0.043536 0.00 0.0 0.0 +1.0000|1111111111>\n",
"0 0 0 9 0.000000 0.042266 0.75 -0.5 0.5 +0.0000|1110111111> +1.0000|1111011111>\n",
"2 1 0 9 0.000000 0.042165 0.75 0.5 0.5 +0.0000|1111111101> +1.0000|1111111110>\n",
"5 2 1 8 1.190955 0.034021 0.00 -0.0 0.0 -0.2346|0111101111> -0.2346|1011110111> +0.00...\n",
"6 2 2 8 1.206474 0.030471 0.00 -0.0 0.0 +0.0686|0111101111> -0.0686|1011110111> -0.19...\n",
"7 2 3 8 1.832629 0.023962 0.00 -0.0 0.0 +0.4135|0111101111> -0.4135|1011110111> +0.04...\n",
"38 14 0 8 0.382981 0.008635 2.00 1.0 1.0 -0.8096|1111101101> -0.5841|1111101110> -0.05...\n",
"35 13 0 8 0.382981 0.008633 2.00 1.0 1.0 +0.0581|1111101011> +0.8096|1111110101> -0.58...\n",
"23 11 0 8 0.382981 0.008632 2.00 -0.0 0.0 -0.0411|0111111011> -0.5725|1011111101> +0.41...\n",
"29 12 0 8 0.382981 0.008625 2.00 0.0 0.0 +0.5725|0111111101> +0.4130|0111111110> +0.04...\n"
]
}
],
"source": [
"print(res_5orb.sort_values('prob', ascending=False)[:15])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This shows a nice representation of the top occupied eigenstates and their quantum numbers.\n",
"\n",
"A side note: as earlier mentioned, one can see that the eigenstates are linear combinations of Fock states here. This becomes important when for example analyzing the occupation of a single orbital in that given subspace. Let's assume we want to know what the probability is of the first orbital to be occupied in the $N=8$ particle sector. To obtain the result one needs to evaluate $Tr(\\rho n_0)$ in the given subspace. $n_0$ is diagonal in the Fock basis, but it is not necessarily diagonal in the eigenbasis of the system. Therefore, one needs to first rotate $\\rho$ into the Fock basis: $U \\cdot \\rho \\cdot U^\\dagger$, with $U$ beeing `h_loc_diag.unitary_matrices[sub]`, and then one can obtain the occupation probability of this orbital from the diagonal elements of $\\rho$ in the Fock basis. \n",
"\n",
"\n",
"The results obtained from the multiplet analysis can now be easily accessed for further post-processing. For example the occupation probabilities per particle number are obtained by running the command:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"N\n",
"0 0.000000e+00\n",
"1 0.000000e+00\n",
"2 0.000000e+00\n",
"3 0.000000e+00\n",
"4 9.937759e-08\n",
"5 2.319210e-05\n",
"6 1.790289e-03\n",
"7 5.183764e-02\n",
"8 4.129959e-01\n",
"9 4.898165e-01\n",
"10 4.353634e-02\n",
"Name: prob, dtype: float64\n"
]
}
],
"source": [
"print(res_5orb.groupby(['N']).sum()['prob'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting the data is also achieved in similar ways. Here, we plot the particle number histogram, as well as the particle number histogram resolved by spin state:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, (ax1) = plt.subplots(1,1,figsize=(8,5))\n",
"\n",
"mult_occ_5orb = res_5orb.groupby(['N']).sum()['prob']\n",
"mult_occ_5orb.plot.bar(y='prob' ,ax=ax1, rot=0)\n",
"\n",
"ax1.set_xlabel(r'N')\n",
"ax1.set_ylabel(r'prob amplitude')\n",
"plt.show()\n",
"\n",
"# split into ms occupations\n",
"fig, (ax1) = plt.subplots(1,1,figsize=(8,5))\n",
"\n",
"spin_occ_five = res_5orb.groupby(['N', '|m_s|']).sum()\n",
"pivot_df = spin_occ_five.pivot_table(index='N', columns='|m_s|', values='prob')\n",
"pivot_df.plot.bar(stacked = True, rot=0, ax = ax1)\n",
"\n",
"ax1.set_ylabel(r'prob amplitude') \n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.2"
}
},
"nbformat": 4,
"nbformat_minor": 4
}