# -*- coding: utf-8 -*-
################################################################################
#
# TPRF: Two-Particle Response Function (TPRF) Toolbox for TRIQS
#
# Copyright (C) 2018 by The Simons Foundation
# Author: H. U.R. Strand
#
# TPRF is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TPRF is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TPRF. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
""" TRIQS: Hartree-Fock Solver
for systems with general dispersion and local interaction
Author: Hugo U. R. Strand, hugo.strand@gmail.com (2018)
"""
# ----------------------------------------------------------------------
import sys
import itertools
import numpy as np
# ----------------------------------------------------------------------
from scipy.optimize import fsolve
from scipy.optimize import brentq
# ----------------------------------------------------------------------
from triqs.gf.block_gf import fix_gf_struct_type
import triqs.utility.mpi as mpi
# ----------------------------------------------------------------------
from triqs_tprf.rpa_tensor import get_rpa_tensor
from triqs_tprf.rpa_tensor import fundamental_operators_from_gf_struct
from triqs_tprf.OperatorUtils import is_operator_composed_of_only_fundamental_operators
# ----------------------------------------------------------------------
[docs]
class HartreeFockSolver(object):
""" Hartree-Fock solver for local interactions.
Parameters
----------
e_k : TRIQS Greens function on a Brillouin zone mesh
Single-particle dispersion.
beta : float
inverse temperature
H_int : TRIQS Operator instance
Local interaction Hamiltonian
gf_struct : list of pairs of block index and its size
gf_struct fixing orbital order between e_k and H_int
mu0 : float
chemical potential
mu_min, mu_max : float, optional
range for chemical potential search.
"""
# ------------------------------------------------------------------
def __init__(self, e_k, beta, H_int=None, gf_struct=None,
mu0=0., mu_max=None, mu_min=None):
if mpi.is_master_node():
print(self.logo())
gf_struct = fix_gf_struct_type(gf_struct)
self.mu = mu0
self.beta = beta
self.mu_max = mu_max
self.mu_min = mu_min
self.e_k = e_k.copy()
self.e_k_MF = e_k.copy()
self.n_k = len(self.e_k.mesh)
self.target_shape = self.e_k.target_shape
self.shape_ab = self.e_k.target_shape
self.shape_abcd = list(self.shape_ab) * 2
self.norb = self.target_shape[0]
self.triu_idxs = np.triu_indices(self.norb, k=1)
if mpi.is_master_node():
print('beta =', self.beta)
print(f'mu = {self.mu} -- (min/max = {self.mu_min}/{self.mu_max})')
print('bands =', self.norb)
print('n_k =', len(self.e_k.mesh))
print('H_int =', H_int)
print()
if gf_struct is None:
assert( H_int is None ), \
'Error: gf_struct = None, but H_int is not None'
if H_int is not None:
assert( gf_struct is not None ), \
'Error: H_int = None, but gf_struct is not None'
fundamental_operators = fundamental_operators_from_gf_struct(gf_struct)
if mpi.is_master_node():
print('gf_struct =', gf_struct)
print('fundamental_operators =', fundamental_operators)
print()
assert( is_operator_composed_of_only_fundamental_operators(
H_int, fundamental_operators) ), \
'Error: H_int is incompatible with gf_struct and its fundamental_operators'
self.U_abcd = get_rpa_tensor(H_int, fundamental_operators)
else:
self.U_abcd = np.zeros(self.shape_abcd)
# ------------------------------------------------------------------
def update_mean_field(self, rho_ab):
self.M = np.einsum('abcd,cd->ba', -self.U_abcd, rho_ab)
return self.M
# ------------------------------------------------------------------
def update_mean_field_dispersion(self):
self.e_k_MF.data[:] = self.e_k.data + self.M[None, ...]
# ------------------------------------------------------------------
def update_chemical_potential(self, N_target, mu0=None):
if mu0 is None:
mu0 = self.mu
e = np.linalg.eigvalsh(self.e_k_MF.data)
fermi = lambda e : 1./(np.exp(self.beta * e) + 1)
def target_function(mu):
n = np.sum(fermi(e - mu)) / self.n_k
return n - N_target
mu_min = self.mu_min if self.mu_min is not None else e.min()
mu_max = self.mu_max if self.mu_max is not None else e.max()
mu = brentq(target_function, mu_min, mu_max)
self.mu = mu
# ------------------------------------------------------------------
def update_momentum_density_matrix(self):
e, V = np.linalg.eigh(self.e_k_MF.data)
e -= self.mu
fermi = lambda e : 1./(np.exp(self.beta * e) + 1)
self.rho_kab = np.einsum('kab,kb,kcb->kac', V, fermi(e), np.conj(V))
return self.rho_kab
# ------------------------------------------------------------------
def update_density_matrix(self):
rho_kab = self.update_momentum_density_matrix()
self.rho_ab = np.sum(rho_kab, axis=0) / self.n_k
self.N_tot = np.sum(np.diag(self.rho_ab))
return self.rho_ab
# ------------------------------------------------------------------
def update_non_int_free_energy(self):
e = np.linalg.eigvalsh(self.e_k_MF.data)
e -= self.mu
self.Omega0 = -1./self.beta * np.sum( np.log(1. + np.exp(-self.beta*e)) )
self.Omega0 /= len(self.e_k_MF.mesh)
return self.Omega0
# ------------------------------------------------------------------
def update_kinetic_energy(self):
rho_kab = self.update_momentum_density_matrix()
self.E_kin = np.einsum('kab,kba->', self.e_k.data, rho_kab)
self.E_kin /= len(self.e_k.mesh)
return self.E_kin
# ------------------------------------------------------------------
def update_interaction_energy(self):
#self.E_int = 0.5 * np.einsum('aa,ab,bb->', self.rho, self.U_ab, self.rho)
self.E_int = 0.5 * np.einsum(
'ab,abcd,cd->', self.rho_ab, self.U_abcd, self.rho_ab)
return self.E_int
# ------------------------------------------------------------------
def update_total_energy(self):
E_kin = self.update_kinetic_energy()
E_int = self.update_interaction_energy()
self.E_tot = E_kin + E_int
Omega0 = self.update_non_int_free_energy()
self.Omega = Omega0 + E_int
return self.E_tot
# ------------------------------------------------------------------
def density_matrix_step(self, rho_vec, N_target=None):
rho_ab = self.vec2mat(rho_vec)
self.update_mean_field(rho_ab)
self.update_mean_field_dispersion()
if N_target is not None:
self.update_chemical_potential(N_target, mu0=self.mu)
rho_ab = self.update_density_matrix()
rho_vec = self.mat2vec(rho_ab)
return rho_vec
# ------------------------------------------------------------------
def solve_setup(self, N_target, M0=None, mu0=None):
self.N_target = N_target
if mu0 is None:
self.update_chemical_potential(self.N_target, mu0=0.)
else:
self.mu = mu0
if M0 is None:
M0 = np.zeros(self.target_shape)
M0 = np.array(M0, dtype=complex)
self.M = np.copy(M0)
self.update_mean_field_dispersion()
self.update_chemical_potential(self.N_target)
rho0_ab = self.update_density_matrix()
rho0_vec = self.mat2vec(rho0_ab)
return rho0_vec
# ------------------------------------------------------------------
[docs]
def solve_iter(self, N_target, M0=None, mu0=None,
nitermax=100, mixing=0.5, tol=1e-9):
""" Solve the HF-equations using forward recursion at fixed density.
Parameters
----------
N_target : float
Total density.
M0 : ndarray (2D), optional
Initial mean-field (0 if None).
mu0 : float, optional
Initial chemical potential.
nitermax : int, optional
Maximal number of self consistent iterations.
mixing : float, optional
Linear mixing parameter.
tol : float, optional
Convergence in relative change of the density matrix.
Returns
-------
rho : ndarray (2D)
Local density matrix.
"""
print('MF: Forward iteration')
print('nitermax =', nitermax)
print('mixing =', mixing)
print('tol =', tol)
print()
assert( mixing >= 0. )
assert( mixing <= 1. )
self.mixing = mixing
self.nitermax = nitermax
rho_vec = self.solve_setup(N_target, M0, mu0)
rho_iter = []
for idx in range(self.nitermax):
rho_iter.append(rho_vec)
rho_vec_old = np.copy(rho_vec)
rho_vec_new = self.density_matrix_step(rho_vec, N_target)
norm = np.linalg.norm(rho_vec_old)
drho = np.linalg.norm(rho_vec_old - rho_vec_new) / norm
print('MF: iter, drho = %3i, %2.2E' % (idx, drho))
if drho < tol:
print('MF: Converged drho = %3.3E\n' % drho)
break
rho_vec = (1. - mixing) * rho_vec_old + mixing * rho_vec_new
self.update_total_energy()
print(self.__str__())
rho_iter = np.array(rho_iter)
return rho_iter
# ------------------------------------------------------------------
[docs]
def solve_newton(self, N_target, M0=None, mu0=None):
""" Solve the HF-equations using a quasi Newton method at fixed density.
Parameters
----------
N_tot : float
Total density.
M0 : ndarray (2D), optional
Initial mean-field (0 if None).
mu0 : float, optional
Initial chemical potential.
"""
print('MF: Newton solver')
print()
rho0_vec = self.solve_setup(N_target, M0, mu0)
def target_function(rho_vec):
rho_vec_new = self.density_matrix_step(rho_vec, N_target)
return rho_vec_new - rho_vec
rho_vec = fsolve(target_function, rho0_vec)
self.update_total_energy()
print(self.__str__())
self.density_matrix_step(rho_vec)
# ------------------------------------------------------------------
[docs]
def solve_newton_mu(self, mu, M0=None):
""" Solve the HF-equations using a quasi Newton method at fixed chemical potential.
Parameters
----------
mu0 : float
Initial chemical potential.
M0 : ndarray (2D), optional
Initial mean-field (0 if None).
"""
print('MF: Newton solver')
print()
self.mu = mu
if M0 is None:
M0 = np.zeros(self.target_shape)
M0 = np.array(M0, dtype=complex)
self.M = np.copy(M0)
self.update_mean_field_dispersion()
rho0_ab = self.update_density_matrix()
rho0_vec = self.mat2vec(rho0_ab)
def target_function(rho_vec):
rho_vec_new = self.density_matrix_step(rho_vec)
return rho_vec_new - rho_vec
rho_vec = fsolve(target_function, rho0_vec)
self.update_total_energy()
print(self.__str__())
self.density_matrix_step(rho_vec)
# ------------------------------------------------------------------
def total_density(self):
return self.N_tot
# ------------------------------------------------------------------
def density_matrix(self):
return self.rho_ab
# ------------------------------------------------------------------
def mean_field_matrix(self):
return self.M
# ------------------------------------------------------------------
def chemical_potential(self):
return self.mu
# ------------------------------------------------------------------
def expectation_value(self, op_ab):
return np.einsum('ab,ba->', op_ab, self.rho_ab)
# ------------------------------------------------------------------
def __str__(self):
txt = 'MF: Solver state\n' + \
'beta = ' + str(self.beta) + '\n' + \
'N_tot = ' + str(self.N_tot) + '\n' + \
'E_tot = ' + str(self.E_tot) + '\n' + \
'E_int = ' + str(self.E_int) + '\n' + \
'E_kin = ' + str(self.E_kin) + '\n' + \
'Omega0 = ' + str(self.Omega0) + '\n' + \
'Omega = ' + str(self.Omega) + '\n' + \
'rho_ab =\n' + str(self.rho_ab) + '\n' + \
'mu = ' + str(self.mu) + '\n' + \
'M =\n' + str(self.M) + '\n' + \
'M - mu =\n' + str(self.M - self.mu * np.eye(self.norb)) + '\n'
return txt
# ------------------------------------------------------------------
[docs]
def mat2vec(self, mat):
r""" Converts a unitary matrix to a vector representation
with the order
1. the real diagonal entries
2. the real part of the upper triangular entries
3. the imaginary part of the upper triangular entries
"""
assert( len(mat.shape) == 2 )
diag = np.diag(mat).real
up_tri = mat[self.triu_idxs]
vec = np.concatenate((diag, up_tri.real, up_tri.imag))
return vec
# ------------------------------------------------------------------
[docs]
def vec2mat(self, vec):
""" Converts from vector representation to a unitary matrix
see mat2vec(...) for details. """
assert( len(vec.shape) == 1 )
mat = np.zeros(self.shape_ab, dtype=complex)
diag, up_tri = vec[:self.norb], vec[self.norb:]
re, im = np.split(up_tri, 2)
mat[self.triu_idxs] = re + 1.j * im
mat += np.conj(mat).T # reconstruct lower triangle
mat += np.diag(diag) # add diagonal
return mat
# ------------------------------------------------------------------
def logo(self):
if 'UTF' in sys.stdout.encoding:
logo = """
╔╦╗╦═╗╦╔═╗ ╔═╗ ┬ ┬┌─┐
║ ╠╦╝║║═╬╗╚═╗ ├─┤├┤
╩ ╩╚═╩╚═╝╚╚═╝ ┴ ┴└
TRIQS: Hartree-Fock solver
"""
else:
logo = """
_____ ___ ___ ___ ___ _ _ ___
|_ _| _ \_ _/ _ \/ __| | || | __|
| | | /| | (_) \__ \ | __ | _|
|_| |_|_\___\__\_\___/ |_||_|_|
TRIQS: Hartree-Fock solver
"""
return logo
# ----------------------------------------------------------------------
[docs]
class HartreeSolver(HartreeFockSolver):
""" Hartree solver for local interactions.
Parameters
----------
e_k : TRIQS Greens function on a Brillouin zone mesh
Single-particle dispersion.
beta : float
inverse temperature
H_int : TRIQS Operator instance
Local interaction Hamiltonian
gf_struct : list of lists of block index and list of internal orbital indices
gf_struct fixing orbital order between e_k and H_int
mu0 : float
chemical potential
mu_min, mu_max : float, optional
range for chemical potential search.
"""
# ------------------------------------------------------------------
[docs]
def mat2vec(self, mat):
assert( len(mat.shape) == 2 )
return np.diag(mat).real
# ------------------------------------------------------------------
[docs]
def vec2mat(self, vec):
assert( len(vec.shape) == 1 )
return np.diag(vec.real)
# ------------------------------------------------------------------
def logo(self):
if 'UTF' in sys.stdout.encoding:
logo = """
╔╦╗╦═╗╦╔═╗ ╔═╗ ┬ ┬
║ ╠╦╝║║═╬╗╚═╗ ├─┤
╩ ╩╚═╩╚═╝╚╚═╝ ┴ ┴
TRIQS: Hartree solver
"""
else:
logo = """
_____ ___ ___ ___ ___ _ _
|_ _| _ \_ _/ _ \/ __| | || |
| | | /| | (_) \__ \ | __ |
|_| |_|_\___\__\_\___/ |_||_|
TRIQS: Hartree solver
"""
return logo
# ----------------------------------------------------------------------