Skip to content
Snippets Groups Projects
Commit b893b3c1 authored by Antonio Manesco's avatar Antonio Manesco
Browse files

clean up and fix k-space sampling to include K and K'

parent eecaf3ba
No related branches found
No related tags found
1 merge request!3create solvers and interface modules
Hide whitespace changes
Inline Side-by-side
examples/graphene_extended_hubbard.ipynb
+ 50
− 65
View file @ b893b3c1
%% Cell type:code id:cb509096-42c6-4a45-8dc4-a8eed3116e67 tags: %% Cell type:code id:cb509096-42c6-4a45-8dc4-a8eed3116e67 tags:
``` python ``` python
import kwant import kwant
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from codes import utils, hf, kwant_examples, model from codes import utils, hf, kwant_examples, model, interface, solvers
from tqdm import tqdm from tqdm import tqdm
``` ```
%% Cell type:markdown id:75721173-f9af-4b9b-a694-7a0bdf65c307 tags: %% Cell type:markdown id:75721173-f9af-4b9b-a694-7a0bdf65c307 tags:
Now we show the interface with `kwant`. We start by using `kwant` to build two tight-binding systems with translational symmetry: Now we show the interface with `kwant`. We start by using `kwant` to build two tight-binding systems with translational symmetry:
* graphene; * graphene;
* a dummy `kwant.Builder` that encodes the interaction matrix. * a dummy `kwant.Builder` that encodes the interaction matrix.
See [`kwant_examples`](./codes/kwant_examples.py) to verify how these two steps are done. See [`kwant_examples`](./codes/kwant_examples.py) to verify how these two steps are done.
%% Cell type:code id: tags:
``` python
def flat_to_matrix(flat, shape):
matrix = np.zeros(shape, dtype=complex)
matrix[:, :, *np.triu_indices(shape[-1])] = flat.reshape(shape[0], shape[1], -1)
lower_triangle_wo_diag = matrix.transpose(0, 1, 3, 2)[
:, :, *np.tril_indices(shape[-1], k=-1)
]
matrix[:, :, *np.tril_indices(shape[-1], k=-1)] = lower_triangle_wo_diag.reshape(
shape[0], shape[1], -1
).conj()
return matrix
def matrix_to_flat(matrix):
matrix[
:,
:,
*np.triu_indices(matrix.shape[-1]),
].flatten()
def upper_triangle_cost(delta_mf_flatten, model):
mfk_shape = model.mf_k.shape
# From flat to non-flat:
mf_delta_shaped = flat_to_matrix(delta_mf_flatten, mfk_shape)
# Doing the update
mf_shaped = model.mf_k + mf_delta_shaped
_, model.mf_k = hf.updated_matrices(mf_k=mf_shaped, model=model)
delta_mf = mf_shaped - model.mf_k
# Now flattening again
return matrix_to_flat(delta_mf)
```
%% Cell type:code id:9ebc3cce-0338-4616-8021-9fecee76f178 tags: %% Cell type:code id:9ebc3cce-0338-4616-8021-9fecee76f178 tags:
``` python ``` python
# Create translationally-invariant `kwant.Builder` # Create translationally-invariant `kwant.Builder`
graphene_builder, int_builder = kwant_examples.graphene_extended_hubbard() graphene_builder, int_builder = kwant_examples.graphene_extended_hubbard()
``` ```
%% Cell type:markdown id:8f004476-fc3b-4c50-808c-4a636ce17c03 tags: %% Cell type:markdown id:8f004476-fc3b-4c50-808c-4a636ce17c03 tags:
We then use `utils.builder2tb_model` to parse the `kwant.Builder` to a `tb_model` that we will use in the self-consistent calculations. We then use `utils.builder2tb_model` to parse the `kwant.Builder` to a `tb_model` that we will use in the self-consistent calculations.
%% Cell type:code id:c2dd6c3c-d9bb-4833-b1a0-d28c5d5a935c tags: %% Cell type:code id:c2dd6c3c-d9bb-4833-b1a0-d28c5d5a935c tags:
``` python ``` python
tb_model = utils.builder2tb_model(graphene_builder) tb_model = utils.builder2tb_model(graphene_builder)
``` ```
%% Cell type:markdown id:075df6f6-9311-4122-8b1e-2d709058e574 tags: %% Cell type:markdown id:075df6f6-9311-4122-8b1e-2d709058e574 tags:
Note that the self-consistent loop is performed on a coarse k-point grid, and thus not necessarily appropriate to compute observables. We thus use `utils.kgrid_hamiltonian` to evaluate the Hamiltonian on a denser k-point grid and compute the gap. Note that the self-consistent loop is performed on a coarse k-point grid, and thus not necessarily appropriate to compute observables. We thus use `utils.kgrid_hamiltonian` to evaluate the Hamiltonian on a denser k-point grid and compute the gap.
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
params = dict(U=1, V=1) params = dict(U=1, V=1)
int_model = utils.builder2tb_model(int_builder, params) int_model = utils.builder2tb_model(int_builder, params)
mf_model = model.Model(tb_model, int_model=int_model) mf_model = model.Model(tb_model, int_model=int_model)
``` ```
%% Cell type:code id:41bd9f60-8f29-4e7c-a0c4-a0bbf66445b2 tags: %% Cell type:code id:41bd9f60-8f29-4e7c-a0c4-a0bbf66445b2 tags:
``` python ``` python
from scipy import optimize
def compute_gap( def compute_gap(
model, model,
nk, nk,
nk_dense, nk_dense,
filling=2, filling=2,
guess=None, guess=None,
): ):
model.guess = guess model.guess = guess
# Find groundstate Hamiltonian on the same grid # Find groundstate Hamiltonian on the same grid
mf_model = hf.find_groundstate_ham( mf_model, guess = interface.find_groundstate_ham(
model, model,
filling=filling, filling=filling,
nk=nk, nk=nk,
cutoff_Vk=1, return_mf=True,
optimizer_kwargs={'M':0} optimizer_kwargs={}
) )
mf_k = utils.kgrid_hamiltonian( # only used for diagonalization to get gap mf_k = utils.kgrid_hamiltonian( # only used for diagonalization to get gap
nk=nk_dense, nk=nk_dense,
hk=utils.model2hk(tb_model=mf_model), hk=utils.model2hk(tb_model=mf_model),
dim=2 dim=2
) )
# Diagonalize groundstate Hamiltonian # Diagonalize groundstate Hamiltonian
vals, _ = np.linalg.eigh(mf_k) vals, _ = np.linalg.eigh(mf_k)
# Extract dense-grid Fermi energy # Extract dense-grid Fermi energy
E_F = utils.get_fermi_energy(vals, filling) E_F = utils.get_fermi_energy(vals, filling)
gap = utils.calc_gap(vals, E_F) gap = utils.calc_gap(vals, E_F)
# the guess was kind of unclear # the guess was kind of unclear
return gap, None#, mf_model return gap, guess
``` ```
%% Cell type:markdown id:718bc267-0899-4d45-8592-deabd6849a75 tags: %% Cell type:markdown id:718bc267-0899-4d45-8592-deabd6849a75 tags:
Finally, we also parse `int_builder` with the wanted interaction strength. Note that we pass a `params` dictionary to evaluate the Hamiltonian with `kwant`. Finally, we also parse `int_builder` with the wanted interaction strength. Note that we pass a `params` dictionary to evaluate the Hamiltonian with `kwant`.
%% Cell type:code id:32b9e7c5-db12-44f9-930c-21e5494404b8 tags: %% Cell type:code id:32b9e7c5-db12-44f9-930c-21e5494404b8 tags:
``` python ``` python
def compute_phase_diagram(Us, Vs, nk, nk_dense): def compute_phase_diagram(Us, Vs, nk, nk_dense):
gap = [] gap = []
for U in tqdm(Us): for U in tqdm(Us):
guess = None guess = None
gap_U = [] gap_U = []
for V in Vs:#tqdm(Vs): for V in Vs:
params = dict(U=U, V=V) params = dict(U=U, V=V)
int_model = utils.builder2tb_model(int_builder, params) int_model = utils.builder2tb_model(int_builder, params)
mf_model = model.Model(tb_model, int_model=int_model) mf_model = model.Model(tb_model, int_model=int_model)
_gap, guess = compute_gap( _gap, guess = compute_gap(
model=mf_model, model=mf_model,
nk=nk, nk=nk,
nk_dense=nk_dense, nk_dense=nk_dense,
guess=guess, guess=guess,
) )
gap_U.append(_gap) gap_U.append(_gap)
gap.append(gap_U) gap.append(gap_U)
return np.asarray(gap, dtype=float) return np.asarray(gap, dtype=float)
``` ```
%% Cell type:markdown id:f1eba14e-e006-4162-885f-3302a92a21eb tags: %% Cell type:markdown id:f1eba14e-e006-4162-885f-3302a92a21eb tags:
**Warning:** this phase diagram calculation takes about one hour. **Warning:** this phase diagram calculation takes about one hour.
%% Cell type:code id:6a8c08a9-7e31-420b-b6b4-709abfb26793 tags: %% Cell type:code id:6a8c08a9-7e31-420b-b6b4-709abfb26793 tags:
``` python ``` python
# Generate dense-grid k-points # Generate dense-grid k-points
# Interaction strengths # Interaction strengths
nk = 20 nk = 9
nk_dense = 30 nk_dense = 30
Us = np.linspace(0, 3, 10, endpoint=True) Us = np.linspace(0, 3, 10, endpoint=True)
Vs = np.linspace(0, 1.5, 10, endpoint=True) Vs = np.linspace(0, 1.5, 10, endpoint=True)
gap = compute_phase_diagram(Us, Vs, nk=nk, nk_dense=nk_dense) gap = compute_phase_diagram(Us, Vs, nk=nk, nk_dense=nk_dense)
``` ```
%% Output %% Output
0%| | 0/10 [00:00<?, ?it/s] 100%|██████████| 10/10 [00:15<00:00, 1.59s/it]
100%|██████████| 10/10 [16:50<00:00, 101.02s/it]
%% Cell type:code id:e17fc96c-c463-4e1f-8250-c254d761b92a tags: %% Cell type:code id:e17fc96c-c463-4e1f-8250-c254d761b92a tags:
``` python ``` python
import xarray as xr import xarray as xr
gap_da = xr.DataArray(data=gap, coords=dict(Us=Us, Vs=Vs)) gap_da = xr.DataArray(data=gap, coords=dict(Us=Us, Vs=Vs))
``` ```
%% Cell type:markdown id:288f3d8e-8432-4697-be78-5156026f9fac tags: %% Cell type:markdown id:288f3d8e-8432-4697-be78-5156026f9fac tags:
We note that the gap openings coincide with the phase transitions from gapless to charge density wave or antiferromagnetic groundstates as predicted in [arXiv:1204.4531](https://arxiv.org/abs/1204.4531). We note that the gap openings coincide with the phase transitions from gapless to charge density wave or antiferromagnetic groundstates as predicted in [arXiv:1204.4531](https://arxiv.org/abs/1204.4531).
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
np.log10(gap_da).plot(x="Us", y="Vs", vmin=-3, vmax=1) np.log10(gap_da).plot(x="Us", y="Vs", vmin=-3, vmax=1)
``` ```
%% Output %% Output
<matplotlib.collections.QuadMesh at 0x13158f150> <matplotlib.collections.QuadMesh at 0x130c85390>
%% Cell type:code id:4f71f6f2 tags:
``` python
np.log10(gap_da).plot(x="Us", y="Vs", vmin=-3, vmax=1)
```
%% Output
%% Cell type:code id:0cb395cd-84d1-49b4-89dd-da7a2d09c8d0 tags: %% Cell type:code id:0cb395cd-84d1-49b4-89dd-da7a2d09c8d0 tags:
``` python ``` python
gap_da.to_netcdf("./data/graphene_example.nc") gap_da.to_netcdf("./data/graphene_example.nc")
``` ```
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment