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

diatomic molecule example

parent fe3af47b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/diatomic_molecule.ipynb 0 → 100644
+ 262
− 0
View file @ 50ff803b
%% Cell type:code id:cb509096-42c6-4a45-8dc4-a8eed3116e67 tags:
``` python
import kwant
import numpy as np
import matplotlib.pyplot as plt
from codes import utils, hf, kwant_examples
from tqdm import tqdm
from itertools import product
```
%% Cell type:code id:d31cbfea-18ea-454e-8a63-d706a85cd3fc tags:
``` python
hamiltonian_0 = np.block([
[0 * np.eye(2), np.eye(2)],
[np.eye(2), 0 * np.eye(2)]
])
hamiltonian_0 = np.expand_dims(hamiltonian_0, axis=0)
hopping_vecs = np.array([[0,]])
```
%% Cell type:code id:46e26a1c-36bd-48bd-89e5-5a7faffa3d1e tags:
``` python
hamiltonian_0.shape
```
%% Output
(1, 4, 4)
%% Cell type:code id:b39a2976-7c35-4670-83ef-12157bd3fc0e tags:
``` python
vals, vecs = np.linalg.eigh(hamiltonian_0)
plt.plot(vals, 'o')
plt.show()
```
%% Output
%% Cell type:code id:41bd9f60-8f29-4e7c-a0c4-a0bbf66445b2 tags:
``` python
def compute_vals(
H_int,
hamiltonian_0=hamiltonian_0,
filling=2,
tol=1e-5,
mixing=0.01,
order=10,
guess=None
):
# Generate guess on the same grid
if guess is None:
guess = utils.generate_guess([0], hopping_vecs, ndof=hamiltonian_0.shape[-1], scale=1)
else:
guess += np.max(guess) * utils.generate_guess([0], hopping_vecs, ndof=hamiltonian_0.shape[-1], scale=0.1)
# Find groundstate Hamiltonian on the same grid
h = hf.find_groundstate_ham(
H_int=H_int,
filling=filling,
hamiltonians_0=hamiltonian_0,
tol=tol,
guess=guess,
mixing=mixing,
order=order,
)
# Diagonalize groundstate Hamiltonian
vals, vecs = np.linalg.eigh(h)
# Extract dense-grid Fermi energy
E_F = utils.get_fermi_energy(vals, filling)
return vals - E_F
```
%% Cell type:code id:32b9e7c5-db12-44f9-930c-21e5494404b8 tags:
``` python
def compute_phase_diagram(
Us,
Vs,
tol=1e-5,
mixing=0.1,
order=5,
):
# onsite interactions
onsite_int = np.block(
[[np.ones((2, 2)), np.zeros((2, 2))], [np.zeros((2, 2)), np.ones((2, 2))]]
)
onsite_int = np.expand_dims(onsite_int, axis=0)
# Nearest-neighbor interactions
nn_int = np.block(
[[np.zeros((2, 2)), np.ones((2, 2))], [np.zeros((2, 2)), np.zeros((2, 2))]]
)
nn_int = np.expand_dims(nn_int, axis=0)
vals = []
for U in tqdm(Us):
gap_U = []
vals_U = []
for V in Vs:
H_int = U * onsite_int + V * nn_int
_vals = compute_vals(
H_int=H_int,
tol=tol,
mixing=mixing,
order=order,
)
vals_U.append(_vals)
vals.append(vals_U)
return np.asarray(vals)
```
%% Cell type:code id:6a8c08a9-7e31-420b-b6b4-709abfb26793 tags:
``` python
# Interaction strengths
Us = np.linspace(0, 5, 20, endpoint=True)
Vs = np.linspace(0, 1, 20, endpoint=True)
vals = compute_phase_diagram(Us, Vs, tol=1e-5)
```
%% Output
100%|██████████| 20/20 [01:02<00:00, 3.15s/it]
%% Cell type:code id:e17fc96c-c463-4e1f-8250-c254d761b92a tags:
``` python
import xarray as xr
ds = xr.Dataset(
data_vars=dict(
vals=(["Us", "Vs", "n"], vals[:,:,0,:]),
),
coords=dict(Us=Us, Vs=Vs, n=np.arange(vals.shape[-1])),
)
```
%% Cell type:code id:868cf368-45a0-465e-b042-6182ff8b6998 tags:
``` python
ds.vals.plot(col='n')
plt.show()
```
%% Output
%% Cell type:code id:0cb395cd-84d1-49b4-89dd-da7a2d09c8d0 tags:
``` python
ds.to_netcdf('./data/diatomic_molecule_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