From 5cde05f05360c80137ff073c72873e58dcd47db0 Mon Sep 17 00:00:00 2001
From: Antonio Manesco <am@antoniomanesco.org>
Date: Fri, 16 Feb 2024 12:44:55 +0100
Subject: [PATCH] add constant correction
---
codes/hf.py | 11 +++++------
1 file changed, 5 insertions(+), 6 deletions(-)
diff --git a/codes/hf.py b/codes/hf.py
index 860fd87..18e9735 100644
--- a/codes/hf.py
+++ b/codes/hf.py
@@ -98,16 +98,15 @@ def compute_mf(rho, H_int):
nk = rho.shape[0]
dim = len(rho.shape) - 2
-
+
if dim > 0:
H0_int = H_int[*([0]*dim)]
- rho0 = rho[*([0]*dim)]
local_density = np.diag(np.average(rho, axis=tuple([i for i in range(dim)])))
exchange_mf = convolution(rho, H_int) * nk ** (-dim)
direct_mf = np.diag(np.einsum("i,ij->j", local_density, H0_int))
- dc_energy_direct = np.einsum("i, j, ij->", local_density, local_density, H0_int)
- dc_energy_cross = np.einsum("ij, ji, ij->", H0_int, rho0, rho0)
- dc_energy = 2 * dc_energy_direct - dc_energy_cross
+ dc_direct = local_density.T @ H0_int @ local_density
+ dc_exchange = np.einsum('kij, kji', exchange_mf, rho) * nk ** (-dim)
+ dc_energy = 0.5*(-dc_exchange + dc_direct) * np.eye(direct_mf.shape[-1])
else:
local_density = np.diag(rho)
exchange_mf = rho * H_int
@@ -115,7 +114,7 @@ def compute_mf(rho, H_int):
dc_energy_direct = np.diag(np.einsum("ij, i, j->", H_int, local_density, local_density))
dc_energy_cross = np.diag(np.einsum("ij, ij, ji->", H_int, rho, rho))
dc_energy = 2 * dc_energy_direct - dc_energy_cross
- return direct_mf - exchange_mf - dc_energy
+ return direct_mf - exchange_mf# - dc_energy * np.eye(direct_mf.shape[-1])
def total_energy(h, rho):
"""
--
GitLab