"""
Coherent cluster approximation (CPA) to substitutional disorder.
For a high-level interface use `solve_root` to solve the CPA problem for
arbitrary frequencies `z`.
Solutions for fixed occupation `occ` can be obtained on the imaginary axis only
using `solve_fxdocc_root`.
Fixed occupation on the real axis is currently not support. We recommend
obtaining the chemical potential `mu` for the given `occ` using
`solve_fxdocc_root` on the imaginary axis, and then run `solve_root` with the
given `mu` on the real axis. In fact, we expect this to be more stable than
fixing the charge on the real axis directly.
"""
import logging
from functools import partial
from typing import Callable, NamedTuple
from typing import Optional as Opt
import numpy as np
from scipy import optimize
from gftool.density import chemical_potential, density_iw
LOGGER = logging.getLogger(__name__)
def _join(*args):
"""
Join arguments to 1D array for use in `optimize`.
Parameters
----------
args : array_like
Arrays to join together.
Returns
-------
joined : np.ndarray
1D array containing all flattened elements of `args`.
shapes : list of tuple of int
Shapes of all `args`. This is necessary to perform the inverse
operation `_split`.
"""
args = [np.asanyarray(ar) for ar in args]
joined = np.concatenate([ar.reshape(-1) for ar in args])
shapes = [ar.shape for ar in args]
return joined, shapes
def _split(joined, shapes):
"""Inverse operation to `_join` separating arrays."""
sizes = [np.prod(sh) for sh in shapes[:-1]]
split = np.split(joined, indices_or_sections=np.cumsum(sizes))
return [array.reshape(sh) for array, sh in zip(split, shapes)]
[docs]
def gf_cmpt_z(z, self_cpa_z, e_onsite, hilbert_trafo: Callable[[complex], complex]):
"""
Green's function for the components embedded in `self_cpa_z`.
Parameters
----------
z, self_cpa_z : (...) complex np.ndarray
Frequency points and corresponding CPA self-energy.
e_onsite : (..., N_cmpt) float of complex np.ndarray
On-site energy of the components. This can also include a local
frequency dependent self-energy of the component sites.
hilbert_trafo : Callable[[complex], complex]
Hilbert transformation of the lattice to calculate the coherent Green's
function.
Returns
-------
(..., N_cmpt) complex np.ndarray
The Green's function of the components embedded in `self_cpa_z`.
"""
gf_coher_z = hilbert_trafo(z - self_cpa_z)[..., np.newaxis]
return gf_coher_z / (1 - (e_onsite - self_cpa_z[..., np.newaxis])*gf_coher_z)
[docs]
def self_root_eq(self_cpa_z, z, e_onsite, concentration,
hilbert_trafo: Callable[[complex], complex]):
"""
Root equation r(Σ)=0 for CPA.
The root equation writes `r(Σ, z) = T(z) / (1 + T(z)*hilbert_trafo(z-Σ))`.
Parameters
----------
self_cpa_z : (...) complex np.ndarray
CPA self-energy.
z : (...) complex array_like
Frequency points.
e_onsite : (..., N_cmpt) float or complex np.ndarray
On-site energy of the components. This can also include a local
frequency dependent self-energy of the component sites.
concentration : (..., N_cmpt) float array_like
Concentration of the different components used for the average.
hilbert_trafo : Callable[[complex], complex]
Hilbert transformation of the lattice to calculate the coherent Green's
function.
Returns
-------
(...) complex np.ndarray
The result of r(Σ), if it is `0` and hence a root, `self_cpa_z` is the
correct CPA self-energy.
"""
gf_coher_z = hilbert_trafo(z - self_cpa_z)
energy_diff = e_onsite - self_cpa_z[..., np.newaxis]
T = np.sum(concentration * energy_diff / (1 - energy_diff*gf_coher_z[..., np.newaxis]), axis=-1)
return T/(1+T*gf_coher_z)
[docs]
def self_fxdpnt_eq(self_cpa_z, z, e_onsite, concentration,
hilbert_trafo: Callable[[complex], complex]):
"""
Fixed-point equation f(Σ)=Σ for CPA.
The fixed-point equation writes `f(Σ, z) = Σ + T(z) / (1 + T(z)*hilbert_trafo(z-Σ))`.
Parameters
----------
self_cpa_z : (...) complex np.ndarray
CPA self-energy.
z : (...) complex array_like
Frequency points.
e_onsite : (..., N_cmpt) float complex np.ndarray
On-site energy of the components. This can also include a local
frequency dependent self-energy of the component sites.
concentration : (..., N_cmpt) float array_like
Concentration of the different components used for the average.
hilbert_trafo : Callable[[complex], complex]
Hilbert transformation of the lattice to calculate the coherent Green's
function.
Returns
-------
(..., N_z) complex np.ndarray
The new self-energy f(Σ), if it is Σ again and hence a fixed-point,
`self_cpa_z_new` is the correct CPA self-energy.
"""
return self_cpa_z + self_root_eq(self_cpa_z, z, e_onsite=e_onsite, concentration=concentration,
hilbert_trafo=hilbert_trafo)
[docs]
def restrict_self_root_eq(self_cpa_z, *args, **kwds):
"""Wrap `self_root_eq` to restrict the solutions to `self_cpa_z.imag > 0`."""
unphysical = self_cpa_z.imag > 0
if np.all(~unphysical): # no need for restrictions
return self_root_eq(self_cpa_z, *args, **kwds)
distance = self_cpa_z.imag[unphysical].copy() # distance to physical solution
# print('>', max(distance))
self_cpa_z.imag[unphysical] = 0
root = np.asanyarray(self_root_eq(self_cpa_z, *args, **kwds))
root[unphysical] *= (1 + distance) # linearly enlarge residues
# kill unphysical roots
root.real[unphysical] += 1e-3 * distance * np.where(root.real[unphysical] >= 0, 1, -1)
root.imag[unphysical] += 1e-3 * distance * np.where(root.imag[unphysical] >= 0, 1, -1)
return root
[docs]
def solve_root(z, e_onsite, concentration, hilbert_trafo: Callable[[complex], complex],
self_cpa_z0=None, restricted=True, **root_kwds):
"""
Determine the CPA self-energy by solving the root problem.
Note that the result should be checked, whether the obtained solution is
physical.
Parameters
----------
z : (...) complex array_like
Frequency points.
e_onsite : (..., N_cmpt) float or complex np.ndarray
On-site energy of the components. This can also include a local
frequency dependent self-energy of the component sites.
concentration : (..., N_cmpt) float array_like
Concentration of the different components used for the average.
hilbert_trafo : Callable[[complex], complex]
Hilbert transformation of the lattice to calculate the coherent Green's
function.
self_cpa_z0 : (...) complex np.ndarray, optional
Starting guess for CPA self-energy.
restricted : bool, optional
Whether `self_cpa_z` is restricted to `self_cpa_z.imag <= 0`. (default: True)
Note, that even if `restricted=True`, the imaginary part can get negative
within tolerance. This should be removed by hand if necessary.
**root_kwds
Additional arguments passed to `scipy.optimize.root`.
`method` can be used to choose a solver.
`options=dict(fatol=tol)` can be specified to set the desired tolerance
`tol`.
Returns
-------
(...) complex np.ndarray
The CPA self-energy as the root of `self_root_eq`.
Raises
------
RuntimeError
If unable to find a solution.
Notes
-----
For `restricted=True` root-serach, we made good experience with the methods
`'anderson'`, `'krylov'` and `'df-sane'`.
For `restricted=False`, we made made good experience with the method `'broyden2'`.
Examples
--------
>>> from functools import partial
>>> parameter = dict(
... e_onsite=[-0.3, 0.3],
... concentration=[0.3, 0.7],
... hilbert_trafo=partial(gt.bethe_gf_z, half_bandwidth=1),
... )
>>> ww = np.linspace(-1.5, 1.5, num=5000) + 1e-10j
>>> self_cpa_ww = gt.cpa.solve_root(ww, **parameter)
>>> del parameter['concentration']
>>> gf_cmpt_ww = gt.cpa.gf_cmpt_z(ww, self_cpa_ww, **parameter)
>>> import matplotlib.pyplot as plt
>>> __ = plt.plot(ww.real, -1./np.pi*gf_cmpt_ww[..., 0].imag)
>>> __ = plt.plot(ww.real, -1./np.pi*gf_cmpt_ww[..., 1].imag)
>>> plt.show()
"""
concentration = np.array(concentration)
if self_cpa_z0 is None: # static average + 0j to make it complex array
self_cpa_z0 = np.sum(e_onsite * concentration, axis=-1) + 0j
self_cpa_z0, __ = np.broadcast_arrays(self_cpa_z0, z)
else: # make sure that `self_cpa_z0` has right shape of output for root
output = np.broadcast(z, e_onsite[..., 0], concentration[..., 0], self_cpa_z0)
self_cpa_z0 = np.broadcast_to(self_cpa_z0, shape=output.shape)
root_eq = partial(restrict_self_root_eq if restricted else self_root_eq,
z=z, e_onsite=e_onsite, concentration=concentration,
hilbert_trafo=hilbert_trafo)
root_kwds.setdefault("method", "anderson" if restricted else "broyden2")
sol = optimize.root(root_eq, x0=self_cpa_z0, **root_kwds)
if not sol.success:
raise RuntimeError(sol.message)
return sol.x
[docs]
class RootFxdocc(NamedTuple):
"""
CPA solution for the self-energy root-equation for fixed occupation.
Parameters
----------
self_cpa : np.ndarray or complex
The CPA self-energy.
mu : float
Chemical potential.
"""
self_cpa: np.ndarray #: The CPA self-energy.
mu: float #: Chemical potential.
[docs]
def solve_fxdocc_root(iws, e_onsite, concentration, hilbert_trafo: Callable[[complex], complex],
beta: float, occ: Opt[float] = None, self_cpa_iw0=None, mu0: float = 0,
weights=1, n_fit=0, restricted=True, **root_kwds) -> RootFxdocc:
"""
Determine the CPA self-energy by solving the root problem for fixed `occ`.
Parameters
----------
iws : (N_iw) complex array_like
Positive fermionic Matsubara frequencies.
e_onsite : (N_cmpt) float or (..., N_iw, N_cmpt) complex np.ndarray
On-site energy of the components. This can also include a local
frequency dependent self-energy of the component sites.
If multiple non-frequency dependent on-site energies should be
considered simultaneously, pass an on-site energy with `N_z=1`:
`e_onsite[..., np.newaxis, :]`.
concentration : (..., N_cmpt) float array_like
Concentration of the different components used for the average.
hilbert_trafo : Callable[[complex], complex]
Hilbert transformation of the lattice to calculate the coherent Green's
function.
beta : float
Inverse temperature.
occ : float
Total occupation.
self_cpa_iw0, mu0 : (..., N_iw) complex np.ndarray and float, optional
Starting guess for CPA self-energy and chemical potential.
`self_cpa_iw0` implicitly contains the chemical potential `mu0`,
thus they should match.
Returns
-------
root.self_cpa : (..., N_iw) complex np.ndarray
The CPA self-energy as the root of `self_root_eq`.
root.mu : float
Chemical potential for the given occupation `occ`.
Other Parameters
----------------
weights : (N_iw) float np.ndarray, optional
Passed to `gftool.density_iw`.
Residues of the frequencies with respect to the residues of the
Matsubara frequencies `1/beta`. (default: 1.)
For Padé frequencies this needs to be provided.
n_fit : int, optional
Passed to `gftool.density_iw`.
Number of additionally fitted moments. If Padé frequencies
are used, this is typically not necessary (default: 0).
restricted : bool, optional
Whether `self_cpa_z` is restricted to `self_cpa_z.imag <= 0`. (default: True)
Note, that even if `restricted=True`, the imaginary part can get negative
within tolerance. This should be removed by hand if necessary.
**root_kwds
Additional arguments passed to `scipy.optimize.root`.
`method` can be used to choose a solver.
`options=dict(fatol=tol)` can be specified to set the desired tolerance
`tol`.
Raises
------
RuntimeError
If unable to find a solution.
See Also
--------
solve_root
Examples
--------
>>> from functools import partial
>>> beta = 30
>>> e_onsite = [-0.3, 0.3]
>>> conc = [0.3, 0.7]
>>> hilbert = partial(gt.bethe_gf_z, half_bandwidth=1)
>>> occ = 0.5,
>>> iws = gt.matsubara_frequencies(range(1024), beta=30)
>>> self_cpa_iw, mu = gt.cpa.solve_fxdocc_root(iws, e_onsite, conc,
... hilbert, occ=occ, beta=beta)
>>> import matplotlib.pyplot as plt
>>> __ = plt.plot(iws.imag, self_cpa_iw.imag, '+--')
>>> __ = plt.axhline(np.average(e_onsite, weights=conc) - mu)
>>> __ = plt.plot(iws.imag, self_cpa_iw.real, 'x--')
>>> plt.show()
check occupation
>>> gf_coher_iw = hilbert(iws - self_cpa_iw)
>>> gt.density_iw(iws, gf_coher_iw, beta=beta, moments=[1, self_cpa_iw[-1].real])
np.float64(0.499999...)
check CPA
>>> self_compare = gt.cpa.solve_root(iws, np.array(e_onsite)-mu, conc,
... hilbert_trafo=hilbert)
>>> np.allclose(self_cpa_iw, self_compare, atol=1e-5)
True
"""
concentration = np.asarray(concentration)[..., np.newaxis, :]
e_onsite = np.asarray(e_onsite)
if self_cpa_iw0 is None: # static average + 0j to make it complex array
self_cpa_iw0 = np.sum(e_onsite * concentration, axis=-1) - mu0 + 0j
self_cpa_iw0, __ = np.broadcast_arrays(self_cpa_iw0, iws)
self_cpa_nomu = self_cpa_iw0 + mu0 # strip contribution of mu
# TODO: use on-site energy to estimate m2+mu, which only has to be adjusted by mu
m1 = np.ones_like(self_cpa_iw0[..., -1].real)
def _occ_diff(x):
gf_coher_iw = hilbert_trafo(iws - x)
m2 = x[..., -1].real # for large iws, real part should static part
occ_root = density_iw(iws, gf_iw=gf_coher_iw, beta=beta, weights=weights,
moments=np.stack([m1, m2], axis=-1), n_fit=n_fit).sum()
return occ_root - occ
mu = chemical_potential(lambda mu: _occ_diff(self_cpa_nomu - mu), mu0=mu0)
LOGGER.debug("VCA chemical potential: %s", mu)
# one iteration gives the ATA: average t-matrix approximation
self_cpa_nomu = self_fxdpnt_eq(self_cpa_nomu - mu, iws, e_onsite - mu,
concentration, hilbert_trafo) + mu
mu = chemical_potential(lambda mu: _occ_diff(self_cpa_nomu - mu), mu0=mu)
LOGGER.debug("ATA chemical potential: %s", mu)
x0, shapes = _join([mu], self_cpa_nomu.real, self_cpa_nomu.imag)
self_root_eq_ = partial(restrict_self_root_eq if restricted else self_root_eq,
z=iws, concentration=concentration, hilbert_trafo=hilbert_trafo)
def root_eq(mu_selfcpa):
mu, self_cpa_re, self_cpa_im = _split(mu_selfcpa, shapes)
self_cpa = self_cpa_re + 1j*self_cpa_im - mu # add contribution of mu
self_root = self_root_eq_(self_cpa, e_onsite=e_onsite - mu)
occ_root = _occ_diff(self_cpa)
return _join([self_root.size*occ_root],
self_root.real, self_root.imag)[0]
root_kwds.setdefault("method", "krylov")
LOGGER.debug('Search CPA self-energy root')
if 'callback' not in root_kwds and LOGGER.isEnabledFor(logging.DEBUG):
# setup LOGGER if no 'callback' is provided
root_kwds['callback'] = lambda _, f: LOGGER.debug(
'Residue: mu=%+6g cpa=%6g', f[0], np.linalg.norm(f[1:])
)
sol = optimize.root(root_eq, x0=x0, **root_kwds)
LOGGER.debug("CPA self-energy root found after %s iterations.", sol.nit)
if not sol.success:
raise RuntimeError(sol.message)
mu, self_cpa_re, self_cpa_im = _split(sol.x, shapes)
self_cpa = self_cpa_re - mu + 1j*self_cpa_im # add contribution of mu
LOGGER.debug("CPA chemical potential: %s", mu.item())
return RootFxdocc(self_cpa, mu=mu.item())