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

cleanup notebook

parent 6d3bf614
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
− 57
View file @ d7e96d28
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from codes import utils, model, interface, solvers, hf from codes import utils, model, interface, solvers, hf
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
U0 = 1 U0 = 1
nk = 1001 nk = 1001
filling = 2 filling = 2
hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(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()} tb_model = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
def model_U(U): def model_U(U):
int_model = { int_model = {
(0,): U * np.kron(np.eye(2), np.ones((2, 2))), (0,): U * np.kron(np.eye(2), np.ones((2, 2))),
} }
return model.Model(tb_model=tb_model, int_model=int_model) return model.Model(tb_model=tb_model, int_model=int_model)
model0 = model_U(U0) model0 = model_U(U0)
model0.vectors = [*model0.int_model.keys()] model0.vectors = [*model0.int_model.keys()]
model0.random_guess(model0.vectors) model0.random_guess(model0.vectors)
model0.kgrid_evaluation(nk=100) model0.kgrid_evaluation(nk=100)
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
hamiltonians_0 = utils.kgrid_hamiltonian( hamiltonians_0 = utils.kgrid_hamiltonian(
nk=nk, hk=utils.model2hk(tb_model=tb_model), dim=1 nk=nk, hk=utils.model2hk(tb_model=tb_model), dim=1
) )
def groundstate(U): def groundstate(U):
_model=model_U(U) _model=model_U(U)
tb_mf_k = interface.find_groundstate_ham( tb_mf_k = interface.find_groundstate_ham(
_model, _model,
filling=filling, filling=filling,
nk=nk, nk=nk,
solver=solvers.rspace_solver, solver=solvers.rspace_solver,
cost_function=solvers.real_space_cost, cost_function=solvers.real_space_cost,
return_kspace=True, return_kspace=True,
optimizer_kwargs={}, optimizer_kwargs={},
) )
vals, vecs = np.linalg.eigh(tb_mf_k) vals, vecs = np.linalg.eigh(tb_mf_k)
EF = utils.get_fermi_energy(vals, filling) EF = utils.get_fermi_energy(vals, filling)
densityMatrix = hf.density_matrix(vals, vecs, EF) densityMatrix = hf.density_matrix(vals, vecs, EF)
return tb_mf_k, densityMatrix, _model.H_int return tb_mf_k, densityMatrix, _model.H_int
tb_mf0, densityMatrix0, H_int0 = groundstate(U0) tb_mf0, densityMatrix0, H_int0 = groundstate(U0)
mf0 = tb_mf0 - hamiltonians_0 mf0 = tb_mf0 - hamiltonians_0
@np.vectorize @np.vectorize
def mfRescaled(alpha, mf0=mf0, H_int0=H_int0): def mfRescaled(alpha, mf0=mf0, H_int0=H_int0):
hamiltonian = hamiltonians_0 + mf0 * alpha hamiltonian = hamiltonians_0 + mf0 * alpha
vals, vecs = np.linalg.eigh(hamiltonian) vals, vecs = np.linalg.eigh(hamiltonian)
EF = utils.get_fermi_energy(vals, filling) EF = utils.get_fermi_energy(vals, filling)
densityMatrix = hf.density_matrix(vals, vecs, EF) densityMatrix = hf.density_matrix(vals, vecs, EF)
return hf.total_energy( return hf.total_energy(
hamiltonians_0 + np.sign(alpha) * mf0, hamiltonians_0 + np.sign(alpha) * mf0,
densityMatrix, densityMatrix,
) )
``` ```
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
alphas = np.linspace(-5, 5, 201) alphas = np.linspace(-5, 5, 201)
plt.plot(alphas, mfRescaled(alphas), 'o') plt.plot(alphas, mfRescaled(alphas), 'o')
plt.plot(-alphas, mfRescaled(alphas), 'o') plt.plot(-alphas, mfRescaled(alphas), 'o')
plt.axvline(x=1, c="k", ls="--") plt.axvline(x=1, c="k", ls="--")
plt.axvline(x=-1, c="k", ls="--") plt.axvline(x=-1, c="k", ls="--")
plt.ylabel("Total Energy") plt.ylabel("Total Energy")
plt.xlabel(r"$\alpha$") plt.xlabel(r"$\alpha$")
# plt.ylim(-4.6, -4.5) # plt.ylim(-4.6, -4.5)
plt.show() plt.show()
``` ```
%% Output %% Output
%% Cell type:code id: tags:
``` python
alphas = np.linspace(-5, 5, 201)
plt.plot(alphas, mfRescaled(alphas), 'o')
plt.plot(-alphas, mfRescaled(alphas), 'o')
plt.axvline(x=1, c="k", ls="--")
plt.axvline(x=-1, c="k", ls="--")
plt.ylabel("Total Energy")
plt.xlabel(r"$\alpha$")
# plt.ylim(-4.6, -4.5)
plt.show()
```
%% Output
%% Cell type:code id: tags:
``` python
alphas = np.linspace(-5, 5, 201)
plt.plot(alphas, mfRescaled(alphas), 'o')
plt.plot(-alphas, mfRescaled(alphas), 'o')
plt.axvline(x=1, c="k", ls="--")
plt.axvline(x=-1, c="k", ls="--")
plt.ylabel("Total Energy")
plt.xlabel(r"$\alpha$")
# plt.ylim(-4.6, -4.5)
plt.show()
```
%% Output
%% Cell type:code id: tags: %% Cell type:code id: tags:
``` python ``` python
plt.plot(alphas[:-1], np.diff(mfRescaled(alphas))) plt.plot(alphas[:-1], np.diff(mfRescaled(alphas)))
# plt.plot(-alphas[:-1], -np.diff(mfRescaled(alphas))) # plt.plot(-alphas[:-1], -np.diff(mfRescaled(alphas)))
plt.axhline(0, ls='--', c='k') plt.axhline(0, ls='--', c='k')
plt.axvline(x=1, c="k", ls="--") plt.axvline(x=1, c="k", ls="--")
plt.axvline(x=-1, c="k", ls="--") plt.axvline(x=-1, c="k", ls="--")
# plt.ylim(-4, -2) # plt.ylim(-4, -2)
plt.show() plt.show()
``` ```
%% Output %% Output
......
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