Skip to content
Snippets Groups Projects

Examples

Merged Examples
examples into main
Merged Kostas Vilkelis requested to merge examples into main
Compare
and
Show latest version
5 files
+ 0
− 1373
Compare changes
  • Side-by-side
  • Inline
Files
5
examples/1d_hubbard_totalenergy.ipynb deleted 100644 → 0
+ 0
− 331
%% Cell type:code id:cb509096-42c6-4a45-8dc4-a8eed3116e67 tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import pymf
from tqdm import tqdm
```
%% Cell type:markdown id:396d935c-146e-438c-878b-04ed70aa6d63 tags:
To simulate infinite systems, we provide the corresponding tight-binding model.
We exemplify this construction by computing the ground state of an infinite spinful chain with onsite interactions.
Because the ground state is an antiferromagnet, so we must build a two-atom cell. We name the two sublattices, $A$ and $B$. The Hamiltonian in is:
$$
H_0 = \sum_i c_{i, B}^{\dagger}c_{i, A} + c_{i, A}^{\dagger}c_{i+1, B} + h.c.
$$
We write down the spinful by simply taking $H_0(k) \otimes \mathbb{1}$.
%% Cell type:code id:5529408c-fb7f-4732-9a17-97b0718dab29 tags:
``` python
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()}
```
%% Cell type:markdown id:050f5435-6699-44bb-b31c-8ef3fa2264d4 tags:
To build the tight-binding model, we need to generate a Hamiltonian on a k-point and the corresponding hopping vectors to generate a guess. We then verify the spectrum and see that the bands indeed consistent of two bands due to the Brillouin zone folding.
%% Cell type:code id:b39a2976-7c35-4670-83ef-12157bd3fc0e tags:
``` python
# Set number of k-points
nk = 100
ks = np.linspace(0, 2*np.pi, nk, endpoint=False)
hamiltonians_0 = pymf.tb_to_khamvector(h_0, nk, ks=ks)
# Perform diagonalization
vals, vecs = np.linalg.eigh(hamiltonians_0)
# Plot data
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()
```
%% Output
%% Cell type:markdown id:6ec53b08-053b-4aad-87a6-525dd7f61687 tags:
Here, in the workflow to find the ground state, we use a helper function to build the initial guess. because we don't need a dense k-point grid in the self-consistent loop, we compute the spectrum later on a denser k-point grid.
%% Cell type:markdown id:dc59e440-1289-4735-9ae8-b04d0d13c94a tags:
Finally, we compute the eigen0alues for a set of Ualues of $U$. For this case, since the interaction is onsite only, the interaction matrix is simply
$$
H_{int} =
\left(\begin{array}{cccc}
U & U & 0 & 0\\
U & U & 0 & 0\\
0 & 0 & U & U\\
0 & 0 & U & U
\end{array}\right)~.
$$
%% Cell type:code id:32b9e7c5-db12-44f9-930c-21e5494404b8 tags:
``` python
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, optimizer_kwargs={"M":0})
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_khamvector(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 tqdm(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)
```
%% Cell type:code id:6a8c08a9-7e31-420b-b6b4-709abfb26793 tags:
``` python
# Interaction strengths
Us = np.linspace(0.5, 10, 20, endpoint=True)
nk, nk_dense = 40, 100
gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)
```
%% Output
75%|███████▌ | 15/20 [00:02<00:00, 5.07it/s]/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=9.07651e-17): result may not be accurate.
gamma = solve(self.a, df_f)
/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=8.54968e-17): result may not be accurate.
gamma = solve(self.a, df_f)
/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=8.17567e-17): result may not be accurate.
gamma = solve(self.a, df_f)
/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=8.88621e-17): result may not be accurate.
gamma = solve(self.a, df_f)
/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=9.35564e-17): result may not be accurate.
gamma = solve(self.a, df_f)
/Users/rzijderveld/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py:1072: LinAlgWarning: Ill-conditioned matrix (rcond=9.73285e-17): result may not be accurate.
gamma = solve(self.a, df_f)
85%|████████▌ | 17/20 [00:06<00:01, 2.47it/s]
---------------------------------------------------------------------------
NoConvergence Traceback (most recent call last)
/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/1169016709.py in <module>
2 Us = np.linspace(0.5, 10, 20, endpoint=True)
3 nk, nk_dense = 40, 100
----> 4 gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)
/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/3682270308.py in compute_phase_diagram(Us, nk, nk_dense)
26 vals = []
27 for U in tqdm(Us):
---> 28 full_sol = compute_sol(U, h_0, nk)
29 gap, _vals = compute_gap_and_vals(full_sol, nk_dense)
30 gaps.append(gap)
/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/3682270308.py in compute_sol(U, h_0, nk, filling)
5 guess = pymf.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
6 full_model = pymf.Model(h_0, h_int, filling)
----> 7 mf_sol = pymf.solver(full_model, guess, nk=nk)
8 return pymf.add_tb(h_0, mf_sol)
9
~/kwant-scf/pymf/solvers.py in solver(Model, mf_guess, nk, optimizer, optimizer_kwargs)
60 f = partial(cost, Model=Model, nk=nk)
61 result = rparams_to_tb(
---> 62 optimizer(f, mf_params, **optimizer_kwargs), list(Model.h_int), shape
63 )
64 fermi = calculate_fermi_energy(add_tb(Model.h_0, result), Model.filling, nk=nk)
~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py in anderson(F, xin, iter, alpha, w0, M, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, **kw)
~/opt/anaconda3/lib/python3.9/site-packages/scipy/optimize/_nonlin.py in nonlin_solve(F, x0, jacobian, iter, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, full_output, raise_exception)
239 else:
240 if raise_exception:
--> 241 raise NoConvergence(_array_like(x, x0))
242 else:
243 status = 2
NoConvergence: [ 1.86715628e+00 4.47287535e-02 -1.23356829e-04 1.20430729e-03
4.47287535e-02 2.63543104e+00 2.63234161e-03 3.71702031e-03
-1.23356829e-04 2.63234161e-03 1.86002862e+00 4.51677226e-02
1.20430729e-03 3.71702031e-03 4.51677226e-02 2.64245378e+00
0.00000000e+00 1.61192759e-01 -1.51104882e-03 1.47113079e-03
-1.61192759e-01 0.00000000e+00 5.03808388e-03 1.99828384e-03
1.51104882e-03 -5.03808388e-03 0.00000000e+00 1.64032886e-01
-1.47113079e-03 -1.99828384e-03 -1.64032886e-01 0.00000000e+00]
%% Cell type:code id:e17fc96c-c463-4e1f-8250-c254d761b92a tags:
``` python
import xarray as xr
ds = xr.Dataset(
data_vars=dict(vals=(["Us", "ks", "n"], vals), gap=(["Us"], gap)),
coords=dict(
Us=Us,
ks=np.linspace(0, 2 * np.pi, nk_dense),
n=np.arange(vals.shape[-1])
),
)
```
%% Cell type:markdown id:5a87dcc1-208b-4602-abad-a870037ec95f tags:
We observe that as the interaction strength increases, a gap opens due to the antiferromagnetic ordering.
%% Cell type:code id:50f02d06 tags:
``` python
plt.scatter(Us, gap)
```
%% Output
<matplotlib.collections.PathCollection at 0x7fe0a154a130>
%% Cell type:code id:868cf368-45a0-465e-b042-6182ff8b6998 tags:
``` python
ds.vals.plot.scatter(x="ks", hue="Us", ec=None, s=5)
plt.axhline(0, ls="--", 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()
```
%% Output
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
/var/folders/yf/2jcxwgld3l77h6y62fb5ty8ssr_5lq/T/ipykernel_78867/3492086229.py in <module>
----> 1 ds.vals.plot.scatter(x="ks", hue="Us", ec=None, s=5)
2 plt.axhline(0, ls="--", c="k")
3 plt.xticks([0, np.pi, 2 * np.pi], ["$0$", "$\pi$", "$2\pi$"])
4 plt.xlim(0, 2 * np.pi)
5 plt.ylabel("$E - E_F$")
AttributeError: '_PlotMethods' object has no attribute 'scatter'
%% Cell type:markdown id:0761ed33-c1bb-4f12-be65-cb68629f58b9 tags:
The Hartree-Fock dispersion should follow (see [these notes](https://www.cond-mat.de/events/correl11/manuscript/Lechermann.pdf))
$$
\epsilon_{HF}^{\sigma}(\mathbf{k}) = \epsilon(\mathbf{k}) + U \left(\frac{n}{2} + \sigma m\right)
$$
where $m=(\langle n_{i\uparrow} \rangle - \langle n_{i\downarrow} \rangle) / 2$ is the magnetization per atom and $n = \sum_i \langle n_i \rangle$ is the total number of atoms per cell. Thus, for the antiferromagnetic groundstate, $m=1/2$ and $n=2$. The gap thus should be $\Delta=U$. And we can confirm it indeed follows the expected trend.
%% Cell type:code id:ac2eb725-f3bd-4d5b-a509-85d0d0071958 tags:
``` python
ds.gap.plot()
plt.plot(ds.Us, ds.Us)
plt.show()
```
%% Output
%% Cell type:markdown id:06e0d356-558e-40e3-8287-d7d2e0bee8cd tags:
We can also fit
%% Cell type:code id:5499ea62-cf1b-4a13-8191-ebb73ea38704 tags:
``` python
ds.gap.polyfit(dim="Us", deg=1).polyfit_coefficients[0].data
```
%% Output
array(0.9997852)
%% Cell type:code id:0cb395cd-84d1-49b4-89dd-da7a2d09c8d0 tags:
``` python
ds.to_netcdf("./data/1d_hubbard_example.nc")
```
%% Cell type:code id:ce428241 tags:
``` python
```{code-cell} ipython3
def compute_phase_diagram(
Us,
nk,
nk_dense,
filling=2,
):
gap = []
vals = []
for U in tqdm(Us):
# onsite interactions
guess = utils.generate_guess(frozenset(h_int), len(list(h_0.values())[0]))
full_model = Model(h_0, h_int, filling)
mf_sol = solver(full_model, guess, nk=nk)
hkfunc = transforms.tb_to_kfunc(add_tb(h_0, mf_sol))
ks_dense = np.linspace(0, 2 * np.pi, nk_dense, endpoint=False)
hkarray = np.array([hkfunc(kx) for kx in ks_dense])
_vals = np.linalg.eigvalsh(hkarray)
_gap = (utils.compute_gap(add_tb(h_0, mf_sol), fermi_energy=0, nk=nk_dense))
gap.append(_gap)
vals.append(_vals)
return np.asarray(gap, dtype=float), np.asarray(vals)
import xarray as xr
ds = xr.Dataset(
data_vars=dict(vals=(["Us", "ks", "n"], vals), gap=(["Us"], gap)),
coords=dict(
Us=Us,
ks=np.linspace(0, 2 * np.pi, nk_dense),
n=np.arange(vals.shape[-1])
),
)
# Interaction strengths
Us = np.linspace(0.5, 10, 20, endpoint=True)
nk, nk_dense = 40, 100
gap, vals = compute_phase_diagram(Us=Us, nk=nk, nk_dense=nk_dense)
ds.vals.plot.scatter(x="ks", hue="Us", ec=None, s=5)
plt.axhline(0, ls="--", 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 Hartree-Fock dispersion should follow (see [these notes](https://www.cond-mat.de/events/correl11/manuscript/Lechermann.pdf))
$$
\epsilon_{HF}^{\sigma}(\mathbf{k}) = \epsilon(\mathbf{k}) + U \left(\frac{n}{2} + \sigma m\right)
$$
where $m=(\langle n_{i\uparrow} \rangle - \langle n_{i\downarrow} \rangle) / 2$ is the magnetization per atom and $n = \sum_i \langle n_i \rangle$ is the total number of atoms per cell. Thus, for the antiferromagnetic groundstate, $m=1/2$ and $n=2$. The gap thus should be $\Delta=U$. And we can confirm it indeed follows the expected trend.
```