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

remove line label

parent fa38b895
No related branches found
No related tags found
1 merge request!3create solvers and interface modules
Hide whitespace changes
Inline Side-by-side
examples/tEnergyTest.ipynb
+ 1
− 1
View file @ e3f1ee69
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
from codes import utils, model, interface, solvers, hf
```
%% Cell type:code id: tags:
``` python
U0 = 1
nk = 100
filling = 2
hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
tb_model = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
def model_U(U):
int_model = {
(0,): U * np.kron(np.ones((2, 2)), np.eye(2)),
}
return model.Model(tb_model=tb_model, int_model=int_model)
model0 = model_U(U0)
model0.vectors = [*model0.int_model.keys()]
model0.random_guess(model0.vectors)
model0.kgrid_evaluation(nk=100)
```
%% Cell type:code id: tags:
``` python
hamiltonians_0 = utils.kgrid_hamiltonian(
nk=nk, hk=utils.model2hk(tb_model=tb_model), dim=1
)
def groundstate(U):
tb_mf_k = interface.find_groundstate_ham(
model_U(U),
filling=filling,
nk=nk,
solver=solvers.rspace_solver,
cost_function=solvers.real_space_cost,
return_kspace=True,
)
vals, vecs = np.linalg.eigh(tb_mf_k)
EF = utils.get_fermi_energy(vals, filling)
densityMatrix = hf.density_matrix(vals, vecs, EF)
return tb_mf_k - EF * np.eye(hamiltonians_0.shape[-1]), densityMatrix
tb_mf0, densityMatrix0 = groundstate(U0)
mf0 = tb_mf0 - hamiltonians_0
@np.vectorize
def mfRescaled(alpha, mf0=mf0):
hamiltonian = hamiltonians_0 + mf0 * alpha
vals, vecs = np.linalg.eigh(hamiltonian)
EF = utils.get_fermi_energy(vals, filling)
densityMatrix = hf.density_matrix(vals, vecs, EF)
return hf.total_energy(tb_mf0, densityMatrix)
```
%% Cell type:code id: tags:
``` python
alphas = np.linspace(0.5, 1.5, 21)
plt.plot(alphas, mfRescaled(alphas), label="MF groundstate from U=1")
plt.plot(alphas, mfRescaled(alphas))
plt.axvline(x=1, c="k", ls="--")
plt.ylabel("Total Energy")
plt.xlabel(r"$\alpha$")
plt.legend()
```
%% Output
<matplotlib.legend.Legend at 0x123524050>
%% Cell type:code id: tags:
``` python
```
......
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