Skip to content
Snippets Groups Projects
Commit 91b40732 authored by Kostas Vilkelis's avatar Kostas Vilkelis :flamingo:
Browse files

merge from main

parent bcece409
No related branches found
No related tags found
1 merge request!14Documentation
# pymf authors
## Current pymf maintainers
- Kostas Vilkelis
- R. Johanna Zijderveld
- Anton R. Akhmerov
- Antonio L.R. Manesco
## Other contributors
- Isidora Araya Day
- José L. Lado
## Funding
BSD 2-Clause License
Copyright (c) 2024, pymf authors
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions and the following disclaimer in the documentation
and/or other materials provided with the distribution.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
# Project Name
# `pymf`
## Research Goal
## What is `pymf`?
## Research Plan
`pymf` is a Python package that performs self-consistent Hartree-Fock calculations on tight-binding models.
It aims to find the groundstate of a Hamiltonian with density-density interactions
## Working on this project
Configure the project by running
$$
\hat{H} = \hat{H_0} + \hat{V} = \sum_{ij} h_{ij} c^\dagger_{i} c_{j} + \frac{1}{2} \sum_{ij} v_{ij} \hat{n}_i \hat{n}_j,
$$
./setup
and computes the mean-field correction $\hat{V}_{\text{MF}}$ which approximates the interaction term:
$$
\hat{V} \approx \hat{V}_{\text{MF}} = \sum_{ij} \tilde{v}_{ij} c^\dagger_{i} c_{j}.
$$
For more details, refer to the [theory overview](docs/source/documentation/mf_notes.md) and [algorithm description](docs/source/documentation/algorithm.md).
## How to use `pymf`?
The calculation of a mean-field Hamiltonian is a simple 3-step process:
1. **Define**
To specify the interacting problem, use a `Model` object which collects:
- Non-interacting Hamiltonian as a tight-binding dictionary.
- Interaction Hamiltonian as a tight-binding dictionary.
- Particle filling number in the unit cell.
2. **Guess**
Construct a starting guess for the mean-field correction.
3. **Solve**
Solve for the mean-field correction using the `solver` function and add it to the non-interacting part to obtain the total mean-field Hamiltonian.
```python
import pymf
#Define
h_0 = {(0,) : onsite, (1,) : hopping, (-1,) : hopping.T.conj()}
h_int = {(0,) : onsite_interaction}
model = pymf.Model(h_0, h_int, filling=2)
#Guess
guess = pymf.generate_guess(guess_hopping_keys, ndof)
#Solve
mf_correction = pymf.solver(model, guess)
h_mf = pymf.add_tb(h_0, mf_correction)
```
For more details and examples on how to use the package, we refer to the [tutorials](docs/source/tutorial/hubbard_1d.md).
## Why `pymf`?
Here is why you should use `pymf`:
* Simple
The workflow is straightforward.
Interface with `Kwant` allows easy creation of complicated tight-binding systems and interactions.
* Extensible
`pymf`'s code is structured to be easy to understand, modify and extend.
* Optimized numerical workflow
Introduces minimal overhead to the calculation of the mean-field Hamiltonian.
## What `pymf` doesn't do (yet)
Here are some features that are not yet implemented but are planned for future releases:
- **Superconductive order parameters**. Mean-field Hamiltonians do not include pairing terms.
- **General interactions**. We allow only density-density interactions (e.g. Coulomb) which can be described by a second-order tensor.
- **Temperature effects**. Density matrix calculations are done at zero temperature.
## Installation
```
pip install pymf
```
## Citing `pymf`
If you have used `pymf` for work that has lead to a scientific publication, please cite us as:
```bibtex
@misc{pymf,
author = {Vilkelis, Kostas and Zijderveld, R. Johanna and Akhmerov, Anton R. and Manesco, Antonio L.R.},
doi = {10.5281/zenodo.11149850},
month = {5},
title = {pymf},
year = {2024}
}
```
html[data-theme="light"] {
--pst-color-primary: rgb(4, 87, 13);
--pst-color-secondary: rgb(66, 96, 150);
--pst-color-border-tippy: rgb(0, 0, 0);
--pst-color-inline-code-links: rgb(4, 87, 13);
--pst-color-link-hover: var(--pst-color-secondary)
}
html[data-theme="dark"] {
--pst-color-primary: rgb(176, 217, 162);
--pst-color-secondary: rgb(125, 166, 244);
--pst-color-border-tippy: white;
--pst-color-inline-code-links: var(--pst-color-primary);
--pst-color-link-hover: var(--pst-color-secondary)
}
.tippy-box {
background-color:var(--pst-color-surface);
color:var(--pst-color-text-base);
border: 1px solid var(--pst-color-border-tippy);
}
......@@ -88,6 +88,7 @@ autoclass_content = "both"
# a list of builtin themes.
#
html_theme = "sphinx_book_theme"
html_title = "PYMF"
html_theme_options = {
"repository_url": "https://gitlab.kwant-project.org/qt/kwant-scf",
......@@ -114,5 +115,5 @@ html_theme_options = {
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
# html_static_path = ["_static"]
# html_css_files = ["local.css"]
html_static_path = ["_static"]
html_css_files = ["css/custom.css"]
---
jupytext:
text_representation:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.14.4
kernelspec:
display_name: Python 3 (ipykernel)
language: python
name: python3
---
# Algorithm overview
## Self-consistency loop
......
---
jupytext:
text_representation:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.14.4
kernelspec:
display_name: Python 3 (ipykernel)
language: python
name: python3
---
# Theory overview
## Interacting problems
......
......@@ -11,13 +11,13 @@ kernelspec:
name: python3
---
# pymf
```{toctree}
:hidden:
:maxdepth: 1
:caption: Tutorials
tutorial/hubbard_1d.md
tutorial/graphene_example.md
```
```{toctree}
......@@ -25,24 +25,11 @@ kernelspec:
:maxdepth: 1
:caption: Documentation
mf_notes.md
algorithm.md
documentation/mf_notes.md
documentation/algorithm.md
documentation/pymf.md
```
## What is pymf?
## Why pymf?
## How does pymf work?
## What does pymf not do yet?
## Installation
```bash
pip install .
```{include} ../../README.md
:relative-docs: docs/source/
```
## Citing
## Contributing
---
jupytext:
text_representation:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.14.4
kernelspec:
display_name: Python 3 (ipykernel)
language: python
name: python3
---
# Interacting graphene
In the previous tutorial, we showed how to use `pymf` to solve a simple 1D Hubbard model with onsite interactions.
In this tutorial, we will apply `pymf` to more complex system: graphene with onsite $U$ and nearest-neighbour $V$ interactions.
The system is more complicated in every aspect: the lattice structure, dimension of the problem, complexity of the interactions.
And yet, the workflow is the same as in the previous tutorial and remains simple and straightforward.
## Building the system with `kwant`
### Non-interacting part
As in the previous tutorial, we could construct a tight-binding dictionary of graphene by hand, but instead it is much easier to use [`kwant`](https://kwant-project.org/) to build the system.
For a more detailed explanation on `kwant` see the [tutorial](https://kwant-project.org/doc/1/tutorial/graphene).
```{code-cell} ipython3
import numpy as np
import matplotlib.pyplot as plt
import kwant
import pymf
from pymf.kwant_helper import utils
s0 = np.identity(2)
sx = np.array([[0, 1], [1, 0]])
sy = np.array([[0, -1j], [1j, 0]])
sz = np.diag([1, -1])
# Create graphene lattice
graphene = kwant.lattice.general([(1, 0), (1 / 2, np.sqrt(3) / 2)],
[(0, 0), (0, 1 / np.sqrt(3))], norbs=2)
a, b = graphene.sublattices
# Create bulk system
bulk_graphene = kwant.Builder(kwant.TranslationalSymmetry(*graphene.prim_vecs))
# Set onsite energy to zero
bulk_graphene[a.shape((lambda pos: True), (0, 0))] = 0 * s0
bulk_graphene[b.shape((lambda pos: True), (0, 0))] = 0 * s0
# Add hoppings between sublattices
bulk_graphene[graphene.neighbors(1)] = s0
```
The `bulk_graphene` object is a `kwant.Builder` object that represents the non-interacting graphene system.
To convert it to a tight-binding dictionary, we use the {autolink}`~pymf.kwant_helper.utils.builder_to_tb` function:
```{code-cell} ipython3
h_0 = utils.builder_to_tb(bulk_graphene)
```
### Interacting part
We utilize `kwant` to build the interaction tight-binding dictionary as well.
To define the interactions, we need to specify two functions:
* `onsite_int(site)`: returns the onsite interaction matrix.
* `nn_int(site1, site2)`: returns the interaction matrix between `site1` and `site2`.
We feed these functions to the {autolink}`~pymf.kwant_helper.utils.build_interacting_syst` function, which constructs the `kwant.Builder` object encoding the interactions.
All we need to do is to convert this object to a tight-binding dictionary using the {autolink}`~pymf.kwant_helper.utils.builder_to_tb` function.
```{code-cell} ipython3
def onsite_int(site, U):
return U * sx
def nn_int(site1, site2, V):
return V * np.ones((2, 2))
builder_int = utils.build_interacting_syst(
builder=bulk_graphene,
lattice=graphene,
func_onsite=onsite_int,
func_hop=nn_int,
max_neighbor=1
)
params = dict(U=0.2, V=1.2)
h_int = utils.builder_to_tb(builder_int, params)
```
Because `nn_int` function returns the same interaction matrix for all site pairs, we set `max_neighbor=1` to ensure that the interaction only extends to nearest-neighbours and is zero for longer distances.
## Computing expectation values
As before, we construct {autolink}`~pymf.model.Model` object to represent the full system to be solved via the mean-field approximation.
We then generate a random guess for the mean-field solution and solve the system:
```{code-cell} ipython3
filling = 2
model = pymf.Model(h_0, h_int, filling=2)
int_keys = frozenset(h_int)
ndof = len(list(h_0.values())[0])
guess = pymf.generate_guess(int_keys, ndof)
mf_sol = pymf.solver(model, guess, nk=18)
h_full = pymf.add_tb(h_0, mf_sol)
```
To investigate the effects of interaction on systems with more than one degree of freedom, it is more useful to consider the expectation values of various operators which serve as order parameters.
For example, we can compute the charge density wave (CDW) order parameter which is defined as the difference in the charge density between the two sublattices.
To calculate operator expectation values, we first need to construct the density matrix via the {autolink}`~pymf.mf.construct_density_matrix` function.
We then feed it into {autolink}`~pymf.observables.expectation_value` function together with the operator we want to measure.
In this case, we compute the CDW order parameter by measuring the expectation value of the $\sigma_z$ operator acting on the graphene sublattice degree of freedom.
```{code-cell} ipython3
cdw_operator = {(0, 0): np.kron(sz, np.eye(2))}
rho, _ = pymf.construct_density_matrix(h_full, filling=filling, nk=40)
rho_0, _ = pymf.construct_density_matrix(h_0, filling=filling, nk=40)
cdw_order_parameter = pymf.expectation_value(rho, cdw_operator)
cdw_order_parameter_0 = pymf.expectation_value(rho_0, cdw_operator)
print(f"CDW order parameter for interacting system: {np.round(np.abs(cdw_order_parameter), 2)}")
print(f"CDW order parameter for non-interacting system: {np.round(np.abs(cdw_order_parameter_0), 2)}")
```
We see that the CDW order parameter is non-zero only for the interacting system, indicating the presence of a CDW phase.
## Graphene phase diagram
In the remaining part of this tutorial, we will utilize all the tools we have developed so far to create a phase diagram for the graphene system.
To identify phase changes, it is convenient to track the gap of the system as a function of $U$ and $V$.
To that end, we first create a function that calculates the gap of the system given the tight-binding dictionary and the Fermi energy.
```{code-cell} ipython3
def compute_gap(h, fermi_energy=0, nk=100):
kham = pymf.tb_to_kgrid(h, nk)
vals = np.linalg.eigvalsh(kham)
emax = np.max(vals[vals <= fermi_energy])
emin = np.min(vals[vals > fermi_energy])
return np.abs(emin - emax)
```
And proceed to compute the gap and the mean-field correction for a range of $U$ and $V$ values:
```{code-cell} ipython3
Us = np.linspace(0, 4, 10)
Vs = np.linspace(0, 1.5, 10)
gaps = []
mf_sols = []
for U in Us:
for V in Vs:
params = dict(U=U, V=V)
h_int = utils.builder_to_tb(builder_int, params)
model = pymf.Model(h_0, h_int, filling=filling)
guess = pymf.generate_guess(int_keys, ndof)
mf_sol = pymf.solver(model, guess, nk=18)
mf_sols.append(mf_sol)
gap = compute_gap(pymf.add_tb(h_0, mf_sol), fermi_energy=0, nk=100)
gaps.append(gap)
gaps = np.asarray(gaps, dtype=float).reshape((len(Us), len(Vs)))
mf_sols = np.asarray(mf_sols).reshape((len(Us), len(Vs)))
plt.imshow(gaps.T, extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto')
plt.colorbar()
plt.xlabel('V')
plt.ylabel('U')
plt.title('Gap')
plt.show()
```
This phase diagram has gap openings at the same places as shown in the [literature](https://arxiv.org/abs/1204.4531).
We can now use the stored results in `mf_sols` to fully map out the phase diagram with order parameters.
On top of the charge density wave (CDW), we also expect a spin density wave (SDW) in a different region of the phase diagram.
We construct the SDW order parameter with the same steps as before, but now we need to sum over the expectation values of the three Pauli matrices to account for the $SU(2)$ spin-rotation symmetry.
```{code-cell} ipython3
s_list = [sx, sy, sz]
cdw_list = []
sdw_list = []
for mf_sol in mf_sols.flatten():
rho, _ = pymf.construct_density_matrix(pymf.add_tb(h_0, mf_sol), filling=filling, nk=40)
# Compute CDW order parameter
cdw_list.append(np.abs(pymf.expectation_value(rho, cdw_operator))**2)
# Compute SDW order parameter
sdw_value = 0
for s_i in s_list:
sdw_operator_i = {(0, 0) : np.kron(sz, s_i)}
sdw_value += np.abs(pymf.expectation_value(rho, sdw_operator_i))**2
sdw_list.append(sdw_value)
cdw_list = np.asarray(cdw_list).reshape(mf_sols.shape)
sdw_list = np.asarray(sdw_list).reshape(mf_sols.shape)
```
Finally, we can combine the gap, CDW and SDW order parameters into one plot.
We naively do this by plotting the difference between CDW and SDW order parameters and indicate the gap with the transparency.
```{code-cell} ipython3
import matplotlib.ticker as mticker
normalized_gap = gaps/np.max(gaps)
plt.imshow((cdw_list - sdw_list).T, extent=(Us[0], Us[-1], Vs[0], Vs[-1]), origin='lower', aspect='auto', cmap="coolwarm", alpha=normalized_gap.T, vmin=-2.6, vmax=2.6)
plt.colorbar(ticks=[-2.6, 0, 2.6], format=mticker.FixedFormatter(['SDW', '0', 'CDW']), label='Order parameter', extend='both')
plt.xlabel('V')
plt.ylabel('U')
plt.show()
```
---
jupytext:
text_representation:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.14.4
kernelspec:
display_name: Python 3 (ipykernel)
language: python
name: python3
---
# 1D Hubbard model
## Background physics
To show the basic functionality of the package, we consider a simple interacting electronic system: a 1D chain of sites that allow nearest-neighbor tunneling with strength $t$ and on-site repulsion $U$ between two electrons if they are on the same site.
Such a model is known as the 1D [Hubbard model](https://en.wikipedia.org/wiki/Hubbard_model) and is useful for understanding the onset of insulating phases in interacting metals.
To begin, we first consider the second quantized form of the non-interacting Hamiltonian.
Because we expect the interacting ground state to be antiferromagnetic, we build a two-atom cell and name the two sublattices $A$ and $B$.
These sublattices are identical to each other in the non-interacting case $U=0$.
The non-interacting Hamiltonian reads:
$$
\hat{H_0} = - t \sum_\sigma \sum_i \left(c_{i, B, \sigma}^{\dagger}c_{i, A, \sigma} + c_{i, A, \sigma}^{\dagger}c_{i+1, B, \sigma} + \textrm{h.c}\right).
$$
where $\textrm{h.c}$ is the hermitian conjugate, $\sigma$ denotes spin ($\uparrow$ or $\downarrow$) and $c_{i, A, \sigma}^{\dagger}$ creates an electron with spin $\sigma$ in unit cell $i$ of sublattice $A$.
Next up, is the interacting part of the Hamiltonian:
$$
\hat{V} = U \sum_i \left(n_{i, A, \uparrow} n_{i, A, \downarrow} + n_{i, B, \uparrow} n_{i, B, \downarrow}\right).
$$
where $n_{i, A, \sigma} = c_{i, A, \sigma}^{\dagger}c_{i, A, \sigma}$ is the number operator for sublattice $A$ and spin $\sigma$.
The total Hamiltonian is then $\hat{H} = \hat{H_0} + \hat{V}$.
With the model defined, we can now proceed to input the Hamiltonian into the package and solve it using the mean-field approximation.
## Problem definition
### Non-interacting Hamiltonian
First, let's get the basic imports out of the way.
```{code-cell} ipython3
import numpy as np
import matplotlib.pyplot as plt
import pymf
```
Now let us translate the non-interacting Hamiltonian $\hat{H_0}$ defined above into the basic input format for the package: a **tight-binding dictionary**.
The tight-binding dictionary is a python dictionary where the keys are tuples of integers representing the hopping vectors and the values are the hopping matrices.
For example, a key `(0,)` represents the onsite term in one dimension and a key `(1,)` represents the hopping a single unit cell to the right.
In two dimensions a key `(0,0)` would represent the onsite term and `(1,0)` would represent hopping to the right in the direction of the first reciprocal lattice vector.
In the case of our 1D Hubbard model, we only have an onsite term and hopping a single unit cell to the left and right.
Thus our non-interacting Hamiltonian becomes:
```{code-cell} ipython3
hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
h_0 = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
```
Here `hopp` is the hopping matrix which we define as a kronecker product between sublattice and spin degrees of freedom: `np.array([[0, 1], [0, 0]])` corresponds to the hopping between sublattices and `np.eye(2)` leaves the spin degrees of freedom unchanged.
In the corresponding tight-binding dictionary `h_0`, the key `(0,)` contains hopping within the unit cell and the keys `(1,)` and `(-1,)` correspond to the hopping between the unit cells to the right and left respectively.
To verify the validity of `h_0`, we evaluate it in the reciprocal space using the {autolink}`~pymf.tb.transforms.tb_to_kgrid`, then diagonalize it and plot the band structure:
```{code-cell} ipython3
nk = 50 # number of k-points
ks = np.linspace(0, 2*np.pi, nk, endpoint=False)
hamiltonians_0 = pymf.tb_to_kgrid(h_0, nk)
vals, vecs = np.linalg.eigh(hamiltonians_0)
plt.plot(ks, vals, c="k")
plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
plt.xlim(0, 2 * np.pi)
plt.ylabel("$E - E_F$")
plt.xlabel("$k / a$")
plt.show()
```
which seems metallic as expected.
### Interaction Hamiltonian
We now proceed to define the interaction Hamiltonian $\hat{V}$.
To achieve this, we utilize the same tight-binding dictionary format as before.
Because the interaction Hamiltonian is on-site, it must be defined only for the key `(0,)` and only for electrons on the same sublattice with opposite spins.
Based on the kronecker product structure we defined earlier, the interaction Hamiltonian is:
```{code-cell} ipython3
U=2
s_x = np.array([[0, 1], [1, 0]])
h_int = {(0,): U * np.kron(np.eye(2), s_x),}
```
Here `s_x` is the Pauli matrix acting on the spin degrees of freedom, which ensures that the interaction is only between electrons with opposite spins whereas the `np.eye(2)` ensures that the interaction is only between electrons on the same sublattice.
### Putting it all together
To combine the non-interacting and interaction Hamiltonians, we use the {autolink}`~pymf.model.Model` class.
In addition to the Hamiltonians, we also need to specify the filling of the system --- the number of electrons per unit cell.
```{code-cell} ipython3
filling = 2
full_model = pymf.Model(h_0, h_int, filling)
```
The object `full_model` now contains all the information needed to solve the mean-field problem.
## Solving the mean-field problem
To find a mean-field solution, we first require a starting guess.
In cases where the non-interacting Hamiltonian is highly degenerate, there exists several possible mean-field solutions, many of which are local and not global minima of the energy landscape.
Therefore, the choice of the initial guess can significantly affect the final solution depending on the energy landscape.
Here the problem is simple enough that we can generate a random guess for the mean-field solution through the {autolink}`~pymf.tb.utils.generate_guess` function.
It creates a random Hermitian tight-binding dictionary based on the hopping keys provided and the number of degrees of freedom within the unit cell.
Because the mean-field solution cannot contain hoppings longer than the interaction itself, we use `h_int` keys as an input to {autolink}`~pymf.tb.utils.generate_guess`.
Finally, to solve the model, we use the {autolink}`~pymf.solvers.solver` function which by default employes a root-finding [algorithm](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.anderson.html) to find a self-consistent mean-field solution.
```{code-cell} ipython3
filling = 2
full_model = pymf.Model(h_0, h_int, filling)
guess = pymf.generate_guess(frozenset(h_int), ndof=4)
mf_sol = pymf.solver(full_model, guess, nk=nk)
```
The {autolink}`~pymf.solvers.solver` function returns only the mean-field correction to the non-interacting Hamiltonian in the same tight-binding dictionary format.
To get the full Hamiltonian, we add the mean-field correction to the non-interacting Hamiltonian and plot the band structure just as before:
```{code-cell} ipython3
h_mf = pymf.add_tb(h_0, mf_sol)
hamiltonians = pymf.tb_to_kgrid(h_mf, nk)
vals, vecs = np.linalg.eigh(hamiltonians)
plt.plot(ks, vals, c="k")
plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
plt.xlim(0, 2 * np.pi)
plt.ylabel("$E - E_F$")
plt.xlabel("$k / a$")
plt.show()
```
the band structure now shows a gap at the Fermi level, indicating that the system is in an insulating phase!
We can go further and compute the gap for a wider range of $U$ values:
```{code-cell} ipython3
def compute_sol(U, h_0, nk, filling=2):
h_int = {
(0,): U * np.kron(np.eye(2), np.ones((2, 2))),
}
guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
full_model = pymf.Model(h_0, h_int, filling)
mf_sol = pymf.solver(full_model, guess, nk=nk)
return pymf.add_tb(h_0, mf_sol)
def compute_gap_and_vals(full_sol, nk_dense, fermi_energy=0):
h_kgrid = pymf.tb_to_kgrid(full_sol, nk_dense)
vals = np.linalg.eigvalsh(h_kgrid)
emax = np.max(vals[vals <= fermi_energy])
emin = np.min(vals[vals > fermi_energy])
return np.abs(emin - emax), vals
def compute_phase_diagram(
Us,
nk,
nk_dense,
):
gaps = []
vals = []
for U in Us:
full_sol = compute_sol(U, h_0, nk)
gap, _vals = compute_gap_and_vals(full_sol, nk_dense)
gaps.append(gap)
vals.append(_vals)
return np.asarray(gaps, dtype=float), np.asarray(vals)
Us = np.linspace(0, 4, 40, endpoint=True)
gap, vals = compute_phase_diagram(Us=Us, nk=20, nk_dense=100)
plt.plot(Us, gap, c="k")
plt.xlabel("$U / t$")
plt.ylabel("$\Delta{E}/t$")
plt.show()
```
We see that at around $U=1$ the gap opens up and the system transitions from a metal to an insulator. In order to more accurately determine the size of the gap, we chose to use a denser k-grid for the diagonalization of the mean-field solution.
Source diff could not be displayed: it is too large. Options to address this: view the blob.
../codes/
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment