Skip to content
Snippets Groups Projects
Commit 4ed37a68 authored by Kostas Vilkelis's avatar Kostas Vilkelis :flamingo:
Browse files

push total energy check

parent 17a7647d
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 0 → 100644
+ 155
− 0
View file @ 4ed37a68
%% 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
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()}
# define interaction
def intModel(U):
model = {
(0,): U * np.kron(np.ones((2, 2)), np.eye(2)),
}
return model
```
%% Cell type:code id: tags:
``` python
U0 = 1
nk = 100
filling = 2
hamiltonians_0 = utils.kgrid_hamiltonian(
nk=nk, hk=utils.model2hk(tb_model=tb_model), dim=1
)
def groundstate(U):
tb_mf = model.Model(tb_model=tb_model, int_model=intModel(U))
tb_mf_k = interface.find_groundstate_ham(
tb_mf,
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, densityMatrix
@np.vectorize
def groundstateEnergy(U):
_, densityMatrix = groundstate(U)
Vk = utils.model2hk(tb_model=intModel(U))
H_int = utils.kgrid_hamiltonian(nk=nk, hk=Vk, dim=1)
mf = hf.compute_mf(densityMatrix, H_int)
return hf.total_energy(mf + hamiltonians_0, densityMatrix)
tb_mf0, densityMatrix0 = groundstate(U0)
print(groundstateEnergy(U0))
@np.vectorize
def otherEnergy(U):
# Total Energy with density matrix from the groundstate at U0
Vk = utils.model2hk(tb_model=intModel(U))
H_int = utils.kgrid_hamiltonian(nk=nk, hk=Vk, dim=1)
mf = hf.compute_mf(densityMatrix0, H_int)
return hf.total_energy(mf + hamiltonians_0, densityMatrix0)
```
%% Output
-2.2252238082028337
%% Cell type:code id: tags:
``` python
U_list = np.linspace(U0 - 0.2, U0 + 0.2, 10)
plt.plot(U_list, otherEnergy(U_list), label='MF groundstate from U=1')
plt.plot(U_list, groundstateEnergy(U_list), label='MF groundstate')
plt.axvline(x=1, c='k', ls='--')
plt.ylabel('Total Energy')
plt.xlabel('U')
plt.legend()
```
%% Output
<matplotlib.legend.Legend at 0x7ff797ad0450>
%% 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