From 1deb3b242360ff5a31b0a1164393ef214cd61f72 Mon Sep 17 00:00:00 2001
From: Antonio Manesco <am@antoniomanesco.org>
Date: Tue, 19 Dec 2023 22:48:03 +0100
Subject: [PATCH] add real-space solver
---
codes/hf.py | 83 +++++++++++++++++++++++++----------------------------
1 file changed, 39 insertions(+), 44 deletions(-)
diff --git a/codes/hf.py b/codes/hf.py
index 254f1fb..4e4fb62 100644
--- a/codes/hf.py
+++ b/codes/hf.py
@@ -2,7 +2,7 @@ from scipy.ndimage import convolve
import numpy as np
import codes.utils as utils
from functools import partial
-from scipy.optimize import anderson, minimize
+from scipy import optimize
def density_matrix(vals, vecs, E_F):
"""
@@ -158,7 +158,7 @@ def updated_matrices(mf_k, model):
rho = density_matrix(vals=vals, vecs=vecs, E_F=E_F)
return rho, compute_mf(rho=rho, H_int=model.H_int) - E_F * np.eye(model.hamiltonians_0.shape[-1])
-def default_cost(mf, model):
+def rspace_solver(model, optimizer, optimizer_kwargs):
"""
Default cost function.
@@ -177,29 +177,36 @@ def default_cost(mf, model):
* Non-hermitian part of the mean-field correction.
* The commutator between the mean-field Hamiltonian and density matrix.
"""
- mf_dict = {}
- for i, key in enumerate(model.guess.keys()):
- mf_dict[key] = mf[i]
- mf = utils.kgrid_hamiltonian(
- nk=model.nk,
- hk=utils.model2hk(mf_dict),
- dim=model.dim,
- hermitian=False
+ model.kgrid_evaluation(nk=model.nk)
+ mf = np.array([*model.guess.values()])
+ shape = mf.shape
+
+ def cost_function(mf):
+ mf = utils.flat_to_matrix(utils.real_to_complex(mf), shape)
+ mf_dict = {}
+ for i, key in enumerate(model.guess.keys()):
+ mf_dict[key] = mf[i]
+ mf = utils.kgrid_hamiltonian(
+ nk=model.nk,
+ hk=utils.model2hk(mf_dict),
+ dim=model.dim,
+ hermitian=False
+ )
+ model.rho, model.mf_k = updated_matrices(mf_k=mf, model=model)
+ model.energy = total_energy(h=model.hamiltonians_0 + model.mf_k, rho=model.rho)
+ delta_mf = model.mf_k - mf
+ delta_mf = utils.hk2tb_model(delta_mf, model.vectors, model.ks)
+ delta_mf = np.array([*delta_mf.values()])
+ return utils.matrix_to_flat(utils.complex_to_real(delta_mf))
+
+ _ = optimizer(
+ cost_function,
+ utils.matrix_to_flat(utils.complex_to_real(mf)),
+ **optimizer_kwargs
)
- model.rho, model.mf_k = updated_matrices(mf_k=mf, model=model)
- model.energy = total_energy(h=model.hamiltonians_0 + model.mf_k, rho=model.rho)
- h = model.hamiltonians_0 + model.mf_k
- commutator = (h@model.rho - model.rho@h)
- n_herm = (mf - np.moveaxis(mf, -1, -2).conj())/2
- delta_mf = model.mf_k - mf
- n_herm = [*utils.hk2tb_model(n_herm, model.vectors, model.ks).values()]
- commutator = [*utils.hk2tb_model(commutator, model.vectors, model.ks).values()]
- delta_mf = [*utils.hk2tb_model(delta_mf, model.vectors, model.ks).values()]
- quantities = np.array([commutator, delta_mf, n_herm])
- idx_max = np.argmax(np.linalg.norm(quantities.reshape(3, -1), axis=-1))
- return quantities[idx_max]
-
-def kspace_solver(model, nk, optimizer, optimizer_kwargs):
+
+
+def kspace_solver(model, optimizer, optimizer_kwargs):
"""
Default cost function.
@@ -218,7 +225,7 @@ def kspace_solver(model, nk, optimizer, optimizer_kwargs):
* Non-hermitian part of the mean-field correction.
* The commutator between the mean-field Hamiltonian and density matrix.
"""
- model.kgrid_evaluation(nk=nk)
+ model.kgrid_evaluation(nk=model.nk)
def cost_function(mf):
mf = utils.flat_to_matrix(utils.real_to_complex(mf), model.mf_k.shape)
model.rho, model.mf_k = updated_matrices(mf_k=mf, model=model)
@@ -226,25 +233,12 @@ def kspace_solver(model, nk, optimizer, optimizer_kwargs):
delta_mf = model.mf_k - mf
return utils.matrix_to_flat(utils.complex_to_real(delta_mf))
- mf = optimizer(
+ _ = optimizer(
cost_function,
utils.matrix_to_flat(utils.complex_to_real(model.mf_k)),
**optimizer_kwargs
)
- mf = utils.flat_to_matrix(utils.real_to_complex(mf), model.mf_k.shape)
- h = model.hamiltonians_0 + model.mf_k
- commutator = (h@model.rho - model.rho@h)
- while np.invert(np.isclose(commutator, 0, atol=1e-15)).all():
- model.random_guess(model.vectors)
- model.kgrid_evaluation(nk=nk)
- mf = optimizer(
- cost_function,
- utils.matrix_to_flat(utils.complex_to_real(model.mf_k)),
- **optimizer_kwargs
- )
- mf = utils.flat_to_matrix(utils.real_to_complex(mf), model.mf_k.shape)
- h = model.hamiltonians_0 + mf
- commutator = (h@model.rho - model.rho@h)
+
def find_groundstate_ham(
@@ -253,8 +247,8 @@ def find_groundstate_ham(
filling,
nk=10,
solver=kspace_solver,
- optimizer=anderson,
- optimizer_kwargs={'M':0},
+ optimizer=optimize.anderson,
+ optimizer_kwargs={'verbose': False},
):
"""
Self-consistent loop to find groundstate Hamiltonian.
@@ -280,9 +274,10 @@ def find_groundstate_ham(
model.nk=nk
model.filling=filling
vectors = utils.generate_vectors(cutoff_Vk, model.dim)
- model.vectors=[*vectors, *model.tb_model.keys()]
+ model.vectors=[*model.int_model.keys()]
if model.guess is None:
model.random_guess(model.vectors)
- solver(model, nk, optimizer, optimizer_kwargs)
+ solver(model, optimizer, optimizer_kwargs)
+ model.vectors=[*model.int_model, *model.tb_model.keys()]
assert np.allclose((model.mf_k - np.moveaxis(model.mf_k, -1, -2).conj())/2, 0, atol=1e-15)
return utils.hk2tb_model(model.hamiltonians_0 + model.mf_k, model.vectors, model.ks)
--
GitLab