Matrix valued continuation & Hardy optimization
[1]:
%matplotlib inline
import sys, os
import numpy as np
import scipy
from triqs.plot.mpl_interface import plt,oplot
from h5 import HDFArchive
from triqs.atom_diag import *
from triqs.gf import *
from triqs.operators import c, c_dag, n, dagger
from itertools import product
from triqs.lattice.tight_binding import TBLattice
from triqs.lattice.utils import k_space_path
from triqs_tprf.lattice import lattice_dyson_g0_wk, lattice_dyson_g0_fk
from triqs_Nevanlinna import Solver
from triqs_maxent import *
import seaborn as sns
import scienceplots
plt.style.use(['science','notebook'])
# sns.set_palette('muted')
Warning: could not identify MPI environment!
Starting serial run at: 2024-08-08 15:49:26.812916
/mnt/sw/nix/store/29h1dijh98y9ar6n8hxv78v8zz2pqfzf-python-3.11.7-view/lib/python3.11/site-packages/numpy/core/getlimits.py:549: UserWarning: The value of the smallest subnormal for <class 'numpy.float64'> type is zero.
setattr(self, word, getattr(machar, word).flat[0])
/mnt/sw/nix/store/29h1dijh98y9ar6n8hxv78v8zz2pqfzf-python-3.11.7-view/lib/python3.11/site-packages/numpy/core/getlimits.py:89: UserWarning: The value of the smallest subnormal for <class 'numpy.float64'> type is zero.
return self._float_to_str(self.smallest_subnormal)
/mnt/sw/nix/store/29h1dijh98y9ar6n8hxv78v8zz2pqfzf-python-3.11.7-view/lib/python3.11/site-packages/numpy/core/getlimits.py:549: UserWarning: The value of the smallest subnormal for <class 'numpy.float32'> type is zero.
setattr(self, word, getattr(machar, word).flat[0])
/mnt/sw/nix/store/29h1dijh98y9ar6n8hxv78v8zz2pqfzf-python-3.11.7-view/lib/python3.11/site-packages/numpy/core/getlimits.py:89: UserWarning: The value of the smallest subnormal for <class 'numpy.float32'> type is zero.
return self._float_to_str(self.smallest_subnormal)
Setup simple two orbital 2D Hubbard model on square lattice
[2]:
n_orb = 2
loc = np.zeros((n_orb,n_orb))
for i in range(n_orb):
for j in range(n_orb):
if i != 0 and i==j:
loc[i,j] = 0.2
if j > i or j < i:
loc[i,j] = -0.4
#nearest neighbor hopping
if n_orb == 2:
t = np.array([[-1.,0.1],[0.1,-0.4]])
else:
t = -1.0 * np.eye(n_orb)
tp = 0.1 * np.eye(n_orb) #next nearest neighbor hopping
hop= { (0,0) : loc,
(1,0) : t,
(-1,0) : t,
(0,1) : t,
(0,-1) : t,
(1,1) : tp,
(-1,-1): tp,
(1,-1) : tp,
(-1,1) : tp}
TB = TBLattice(units = [(1, 0, 0) , (0, 1, 0)], hoppings = hop, orbital_positions= [(0., 0., 0.)]*n_orb)
plot dispersion along high-symmetry lines
[3]:
n_pts = 101
G = np.array([ 0.00, 0.00, 0.00])
M = np.array([0.5, 0.5, 0.0])
X = np.array([0.5, 0.0, 0.0])
R = np.array([0.5, 0.5, 0.5])
paths = [(R, G), (G, X), (X, M), (M, G)]
kvecs, k, ticks = k_space_path(paths, num=n_pts, bz=TB.bz)
e_mat = TB.fourier(kvecs)
e_val = TB.dispersion(kvecs)
eps_k = {'k': k, 'K': ticks, 'eval': e_val, 'emat' : e_mat}
[4]:
fig, ax = plt.subplots(1,1, figsize=(8,4), dpi=100, squeeze=False)
ax = ax.reshape(-1)
for band in range(eps_k['eval'].shape[1]):
ax[0].plot(eps_k['k'], eps_k['eval'][:,band].real, color='C0', zorder=2.5)
ax[0].axhline(y=0,zorder=2,color='gray',alpha=0.5,ls='--')
ax[0].set_xticks(eps_k['K'])
ax[0].set_xticklabels([r'R', '$\Gamma$', 'X', 'M', r'$\Gamma$'])
ax[0].set_xlim([eps_k['K'].min(), eps_k['K'].max()])
ax[0].set_ylabel(r'$\omega$ (eV)')
plt.show()
Setup lattice Green’s function and calculate \(G_{loc}\) on Matsubara axis
[5]:
k_grid = (100,100,1)
k_mesh = TB.get_kmesh(n_k = k_grid)
e_k = TB.fourier(k_mesh)
eps_k = TB.dispersion(k_mesh)
mu = 0.
beta = 100
n_iw = 700
iw_mesh = MeshImFreq(beta=beta, S='Fermion', n_max=n_iw)
G0_iwk = lattice_dyson_g0_wk(mu=mu, e_k=e_k, mesh=iw_mesh)
calculate local Green’s function
[6]:
G_Γ_iw = G0_iwk[:, Idx(0,0,0)]
G_X_iw = G0_iwk[:, Idx(k_grid[1]-1,0,0)]
G_loc_iw = Gf(mesh=iw_mesh, target_shape=(n_orb,n_orb))
iw_arr = np.array(list(iw_mesh.values()))
for k in G0_iwk.mesh[1]:
G_loc_iw[:] += G0_iwk[:,k]
G_loc_iw = G_loc_iw / len(k_mesh)
[7]:
fig, ax = plt.subplots(n_orb, n_orb, figsize=(6*n_orb,3*n_orb), dpi=100, squeeze=False,sharex=True)
shp = G_Γ_iw.target_shape
for i in range(n_orb):
for j in range(n_orb):
ax[i,j].oplot(G_loc_iw[i,j].imag, label=f'[{i},{j}] sum')
ax[i,j].oplot(G_Γ_iw[i,j].imag, label=f'[{i},{j}]@$\Gamma$')
ax[i,j].oplot(G_X_iw[i,j].imag, label=f'[{i},{j}]@X')
ax[i,j].set_xlim(0,30)
if i == shp[0]-1:
ax[i,j].set_xlabel(r'$i\omega_n$')
else:
ax[i,j].set_xlabel(r'')
if j == 0:
ax[i,j].set_ylabel(r'Im $G (i\omega_n)$ (eV)')
else:
ax[i,j].set_ylabel(r'')
plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
The local Green’s function is off-diagonal in orbital space.
Calculate reference data directly on real frequency axis
[8]:
n_w = 1001
delta = 0.1
w_window = (-7, 7)
w_mesh = MeshReFreq(window=w_window, n_w=n_w)
G0_wk = lattice_dyson_g0_fk(mu=mu, e_k=e_k, mesh=w_mesh, delta=delta)
G_Γ_w = G0_wk[:, Idx(0,0,0)]
G_X_w = G0_wk[:, Idx(k_grid[1]-1,0,0)]
G_loc_w = Gf(mesh=w_mesh, target_shape=(n_orb,n_orb))
w_arr = np.array(list(w_mesh.values()))
for k in G0_wk.mesh[1]:
G_loc_w[:] += G0_wk[:,k]
G_loc_w = G_loc_w / len(k_mesh)
Analytic continuation with Pade, Nevanlinna, and MaxEnt
Pade
[9]:
G_loc_w_P = G_loc_w.copy()
G_loc_w_P.set_from_pade(G_loc_iw, n_points=n_iw, freq_offset=delta)
Nevanlinna & Caratheodory
[10]:
# setup Nevanlinna kernel solver
nevan_mesh = MeshReFreq(window=w_window, n_w=n_w//2)
solver = Solver(kernel='NEVANLINNA', precision=128)
solver.solve(G_loc_iw)
G_loc_w_N = solver.evaluate(nevan_mesh, eta=delta)
This is Nevanlinna analytical continuation. All off-diagonal elements will be ignored.
[11]:
# setup matrix valued Nevanlinna kernel solver
solver = Solver(kernel='CARATHEODORY', precision=128)
solver.solve(G_loc_iw)
G_loc_w_CN = solver.evaluate(nevan_mesh, eta=delta)
run MaxEnt
[12]:
tm = PoormanMaxEnt(use_hermiticity=True)
tm.set_G_iw(G_loc_iw)
tm.set_error(1.e-5)
tm.omega = HyperbolicOmegaMesh(omega_min=w_window[0], omega_max=w_window[1], n_points=100)
tm.alpha_mesh = LogAlphaMesh(alpha_min=1e-4, alpha_max=1e2, n_points=30)
result = tm.run()
k_diag = tm.maxent_diagonal.K
k_off_diag = tm.maxent_offdiagonal.K
print(result.analyzer_results[0][0]['LineFitAnalyzer']['alpha_index'])
print(result.analyzer_results[0][1]['LineFitAnalyzer']['alpha_index'])
print(result.analyzer_results[1][1]['LineFitAnalyzer']['alpha_index'])
w_mesh_arr = np.linspace(w_mesh.w_min, w_mesh.w_max, len(w_mesh))
G_loc_w_ME = np.zeros((n_orb,n_orb,len(w_mesh_arr)))
for i in range(n_orb):
for j in range(n_orb):
G_loc_w_ME[i,j,:] = np.interp(w_mesh_arr, np.array(result.omega), result.get_A_out('LineFitAnalyzer')[i,j])
Calculating diagonal elements.
Calling MaxEnt for element 0 0
2024-08-08 15:49:50.522895
MaxEnt run
TRIQS application maxent
Copyright(C) 2018 Gernot J. Kraberger
Copyright(C) 2018 Simons Foundation
Authors: Gernot J. Kraberger and Manuel Zingl
This program comes with ABSOLUTELY NO WARRANTY.
This is free software, and you are welcome to redistributeit under certain conditions; see file LICENSE.
Please cite this code and the appropriate original papers (see documentation).
Minimal chi2: 0.37113432732748763
scaling alpha by a factor 4201 (number of data points)
alpha[ 0] = 4.20100000e+05, chi2 = 1.16730199e+04, n_iter= 243
alpha[ 1] = 2.60889217e+05, chi2 = 6.48908374e+03, n_iter= 14
alpha[ 2] = 1.62016624e+05, chi2 = 3.58268230e+03, n_iter= 19
alpha[ 3] = 1.00615068e+05, chi2 = 1.97629995e+03, n_iter= 17
alpha[ 4] = 6.24836620e+04, chi2 = 1.09106752e+03, n_iter= 19
alpha[ 5] = 3.88034127e+04, chi2 = 5.99958550e+02, n_iter= 21
alpha[ 6] = 2.40975767e+04, chi2 = 3.27708715e+02, n_iter= 24
alpha[ 7] = 1.49650034e+04, chi2 = 1.78957155e+02, n_iter= 26
alpha[ 8] = 9.29352064e+03, chi2 = 9.87981823e+01, n_iter= 29
alpha[ 9] = 5.77143377e+03, chi2 = 5.55351358e+01, n_iter= 32
alpha[10] = 3.58415815e+03, chi2 = 3.18850251e+01, n_iter= 37
alpha[11] = 2.22582293e+03, chi2 = 1.87995508e+01, n_iter= 41
alpha[12] = 1.38227375e+03, chi2 = 1.14435629e+01, n_iter= 48
alpha[13] = 8.58415417e+02, chi2 = 7.17101663e+00, n_iter= 58
alpha[14] = 5.33090517e+02, chi2 = 4.59177848e+00, n_iter= 73
alpha[15] = 3.31058243e+02, chi2 = 3.01273928e+00, n_iter= 85
alpha[16] = 2.05592778e+02, chi2 = 2.06250136e+00, n_iter= 97
alpha[17] = 1.27676598e+02, chi2 = 1.50327290e+00, n_iter= 115
alpha[18] = 7.92893304e+01, chi2 = 1.17470857e+00, n_iter= 139
alpha[19] = 4.92400175e+01, chi2 = 9.78901860e-01, n_iter= 161
alpha[20] = 3.05788851e+01, chi2 = 8.61868408e-01, n_iter= 190
alpha[21] = 1.89900057e+01, chi2 = 7.93202792e-01, n_iter= 221
alpha[22] = 1.17931153e+01, chi2 = 7.53699263e-01, n_iter= 266
alpha[23] = 7.32372438e+00, chi2 = 7.30813043e-01, n_iter= 318
alpha[24] = 4.54815692e+00, chi2 = 7.16960474e-01, n_iter= 404
alpha[25] = 2.82448250e+00, chi2 = 7.07904454e-01, n_iter= 518
alpha[26] = 1.75405149e+00, chi2 = 7.00374462e-01, n_iter= 787
alpha[27] = 1.08929569e+00, chi2 = 6.92638321e-01, n_iter= 854
alpha[28] = 6.76471078e-01, chi2 = 6.86081810e-01, n_iter= 1000!
alpha[29] = 4.20100000e-01, chi2 = 6.80797843e-01, n_iter= 1000!
! ... The minimizer did not converge. Results might be wrong.
MaxEnt loop finished in 0:00:37.055558
Calling MaxEnt for element 1 1
Minimal chi2: 0.05463864042652783
scaling alpha by a factor 4201 (number of data points)
alpha[ 0] = 4.20100000e+05, chi2 = 2.91457175e+03, n_iter= 428
alpha[ 1] = 2.60889217e+05, chi2 = 1.41567503e+03, n_iter= 12
alpha[ 2] = 1.62016624e+05, chi2 = 7.18660282e+02, n_iter= 12
alpha[ 3] = 1.00615068e+05, chi2 = 3.81185054e+02, n_iter= 13
alpha[ 4] = 6.24836620e+04, chi2 = 2.12421862e+02, n_iter= 14
alpha[ 5] = 3.88034127e+04, chi2 = 1.26010197e+02, n_iter= 14
alpha[ 6] = 2.40975767e+04, chi2 = 7.99349340e+01, n_iter= 16
alpha[ 7] = 1.49650034e+04, chi2 = 5.34691132e+01, n_iter= 19
alpha[ 8] = 9.29352064e+03, chi2 = 3.66639535e+01, n_iter= 24
alpha[ 9] = 5.77143377e+03, chi2 = 2.49757946e+01, n_iter= 30
alpha[10] = 3.58415815e+03, chi2 = 1.65116350e+01, n_iter= 38
alpha[11] = 2.22582293e+03, chi2 = 1.05024813e+01, n_iter= 46
alpha[12] = 1.38227375e+03, chi2 = 6.47542138e+00, n_iter= 53
alpha[13] = 8.58415417e+02, chi2 = 3.94927055e+00, n_iter= 61
alpha[14] = 5.33090517e+02, chi2 = 2.43448245e+00, n_iter= 80
alpha[15] = 3.31058243e+02, chi2 = 1.52805650e+00, n_iter= 92
alpha[16] = 2.05592778e+02, chi2 = 9.67192877e-01, n_iter= 115
alpha[17] = 1.27676598e+02, chi2 = 6.10734199e-01, n_iter= 147
alpha[18] = 7.92893304e+01, chi2 = 3.87884255e-01, n_iter= 176
alpha[19] = 4.92400175e+01, chi2 = 2.56416965e-01, n_iter= 188
alpha[20] = 3.05788851e+01, chi2 = 1.84413321e-01, n_iter= 213
alpha[21] = 1.89900057e+01, chi2 = 1.47412141e-01, n_iter= 211
alpha[22] = 1.17931153e+01, chi2 = 1.29031391e-01, n_iter= 230
alpha[23] = 7.32372438e+00, chi2 = 1.19864347e-01, n_iter= 306
alpha[24] = 4.54815692e+00, chi2 = 1.15096619e-01, n_iter= 323
alpha[25] = 2.82448250e+00, chi2 = 1.12349399e-01, n_iter= 398
alpha[26] = 1.75405149e+00, chi2 = 1.10082715e-01, n_iter= 546
alpha[27] = 1.08929569e+00, chi2 = 1.08113392e-01, n_iter= 519
alpha[28] = 6.76471078e-01, chi2 = 1.06817567e-01, n_iter= 595
alpha[29] = 4.20100000e-01, chi2 = 1.06020333e-01, n_iter= 724
MaxEnt loop finished in 0:00:31.812850
Calculating off-diagonal elements using default model from diagonal solution
Calling MaxEnt for element 0 1
2024-08-08 15:51:00.628124
MaxEnt run
TRIQS application maxent
Copyright(C) 2018 Gernot J. Kraberger
Copyright(C) 2018 Simons Foundation
Authors: Gernot J. Kraberger and Manuel Zingl
This program comes with ABSOLUTELY NO WARRANTY.
This is free software, and you are welcome to redistributeit under certain conditions; see file LICENSE.
Please cite this code and the appropriate original papers (see documentation).
Minimal chi2: 0.2524130033415118
scaling alpha by a factor 4201 (number of data points)
alpha[ 0] = 4.20100000e+05, chi2 = 2.12746519e+02, n_iter= 5
alpha[ 1] = 2.60889217e+05, chi2 = 1.15037153e+02, n_iter= 4
alpha[ 2] = 1.62016624e+05, chi2 = 6.16385648e+01, n_iter= 4
alpha[ 3] = 1.00615068e+05, chi2 = 3.42476314e+01, n_iter= 4
alpha[ 4] = 6.24836620e+04, chi2 = 2.05439008e+01, n_iter= 4
alpha[ 5] = 3.88034127e+04, chi2 = 1.34474239e+01, n_iter= 4
alpha[ 6] = 2.40975767e+04, chi2 = 9.32054995e+00, n_iter= 5
alpha[ 7] = 1.49650034e+04, chi2 = 6.50461581e+00, n_iter= 7
alpha[ 8] = 9.29352064e+03, chi2 = 4.47131187e+00, n_iter= 8
alpha[ 9] = 5.77143377e+03, chi2 = 3.10803541e+00, n_iter= 10
alpha[10] = 3.58415815e+03, chi2 = 2.24361509e+00, n_iter= 11
alpha[11] = 2.22582293e+03, chi2 = 1.67781559e+00, n_iter= 12
alpha[12] = 1.38227375e+03, chi2 = 1.28091836e+00, n_iter= 14
alpha[13] = 8.58415417e+02, chi2 = 9.93840445e-01, n_iter= 16
alpha[14] = 5.33090517e+02, chi2 = 7.92590341e-01, n_iter= 19
alpha[15] = 3.31058243e+02, chi2 = 6.60831703e-01, n_iter= 23
alpha[16] = 2.05592778e+02, chi2 = 5.80658369e-01, n_iter= 27
alpha[17] = 1.27676598e+02, chi2 = 5.34656810e-01, n_iter= 30
alpha[18] = 7.92893304e+01, chi2 = 5.09232651e-01, n_iter= 34
alpha[19] = 4.92400175e+01, chi2 = 4.95355478e-01, n_iter= 38
alpha[20] = 3.05788851e+01, chi2 = 4.87553900e-01, n_iter= 46
alpha[21] = 1.89900057e+01, chi2 = 4.82615747e-01, n_iter= 60
alpha[22] = 1.17931153e+01, chi2 = 4.78701449e-01, n_iter= 77
alpha[23] = 7.32372438e+00, chi2 = 4.75026539e-01, n_iter= 114
alpha[24] = 4.54815692e+00, chi2 = 4.71457460e-01, n_iter= 147
alpha[25] = 2.82448250e+00, chi2 = 4.66160378e-01, n_iter= 218
alpha[26] = 1.75405149e+00, chi2 = 4.59323377e-01, n_iter= 267
alpha[27] = 1.08929569e+00, chi2 = 4.52652176e-01, n_iter= 409
alpha[28] = 6.76471078e-01, chi2 = 4.46353696e-01, n_iter= 656
alpha[29] = 4.20100000e-01, chi2 = 4.39988779e-01, n_iter= 946
MaxEnt loop finished in 0:00:19.608640
Element 1 0 not calculated, can be determined from hermiticity
17
14
21
Final plot for comparison
[13]:
fig, ax = plt.subplots(n_orb, n_orb, figsize=(8*n_orb,4*n_orb), dpi=100, squeeze=False,sharex=True)
shp = G_loc_w.target_shape
for i in range(n_orb):
for j in range(n_orb):
# plotting results
ax[i,j].oplot(-1*G_loc_w[i,j].imag, lw=5, label=f'[{i},{j}] ref')
ax[i,j].oplot(-1*G_loc_w_N[i,j].imag, label=f'[{i},{j}] Nev')
ax[i,j].oplot(-1*G_loc_w_CN[i,j].imag, label=f'[{i},{j}] Car')
ax[i,j].oplot(-1*G_loc_w_P[i,j].imag, label=f'[{i},{j}] Pade')
ax[i,j].plot(w_mesh_arr,np.pi*G_loc_w_ME[i,j], label=f'[{i},{j}] ME')
ax[i,j].set_xlim(w_window)
ax[i,j].legend()
if i == shp[0]-1:
ax[i,j].set_xlabel(r'$\omega$')
else:
ax[i,j].set_xlabel(r'')
if j == 0:
ax[i,j].set_ylabel(r'Im $G (\omega)$ (eV)')
else:
ax[i,j].set_ylabel(r'')
plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
fig, ax = plt.subplots(n_orb, n_orb, figsize=(8*n_orb,4*n_orb), dpi=100, squeeze=False,sharex=True)
for i in range(n_orb):
for j in range(n_orb):
# plotting results
ax[i,j].oplot(G_loc_w[i,j].real, lw=5, label=f'[{i},{j}] ref')
ax[i,j].oplot(G_loc_w_N[i,j].real, label=f'[{i},{j}] Nev')
ax[i,j].oplot(G_loc_w_CN[i,j].real, label=f'[{i},{j}] Car')
ax[i,j].oplot(G_loc_w_P[i,j].real, label=f'[{i},{j}] Pade')
ax[i,j].set_xlim(w_window)
if i == shp[0]-1:
ax[i,j].set_xlabel(r'$\omega$')
else:
ax[i,j].set_xlabel(r'')
if j == 0:
ax[i,j].set_ylabel(r'Re $G (\omega)$ (eV)')
else:
ax[i,j].set_ylabel(r'')
plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
Hardy optimization
Both Nevanlinna and Caratheodory continuation work well for finite systems. However, for infinite systems / continuous spectra the continuation will return a highly oscillating function. To address this issue an additional optimization of the continuation can be performed.
In the Nevanlinna package such a Hardy function optimization is implemented. As described in PRL 126 056402 (2021) the continued fraction form of the Nevanlinna factorization has a freedom in choosing \(\theta_{M+1}\), the last element in the factorization:
\[G(z) = - h^{-1}\left[\theta(z)\right] \nonumber
\theta(z) = \frac{a(z) \theta_{M+1}(z) +b(z)}{c(z)\theta_{M+1}(z) +d(z)},\]
where \(a(z)\), \(b(z)\), \(c(z)\) and \(d(z)\) are Nevanlinna factorization coefficients, and \(h^{-1}\) is the inverse Mobius transform.
To perform the optimization the user has to define target function for the optimization problem, that takes a retarded Green’s function as an input. The optimization is then performed using SciPy.
First we repeat the continuation using the normal Nevanlinna kernel of the [0,0] orbital of the Green’s function from above.
[14]:
# setup Nevanlinna kernel solver
nevan_mesh = MeshReFreq(window=w_window, n_w=n_w//2)
solver = Solver(kernel='NEVANLINNA')
solver.solve(G_loc_iw[0:1,0:1])
G_w_non_opt = solver.evaluate(nevan_mesh, delta)
This is Nevanlinna analytical continuation. All off-diagonal elements will be ignored.
Now we perform the optimization, by defining first a target function that optimizes the smoothness of the resulting Green’s function. We then perform the optimization for a few eta values:
[15]:
nevan_mesh_arr = np.array([x for x in nevan_mesh.values()])
Gw = G_loc_iw[0:1,0:1]
eta_list = [0.4, 0.2, 0.1]
def SmoothnessOptimization (Gw):
lmbda = 1e-4
diff = 0.0
# loop over all orbitals
for i in range(Gw.data.shape[1]):
# extract spectral function
Aw = -Gw.data[:,i,i].imag/np.pi
# compute estimate for second derivative
Aw2 = (Aw[:-2] - 2*Aw[1:-1] + Aw[2:])/(Gw.mesh.delta**2)
# evaluate normalization criteria weight
diff += np.abs(1 - scipy.integrate.simpson(Aw, x=nevan_mesh_arr))**2
# evaluate smoothness criteria weight
diff += lmbda * scipy.integrate.simpson(Aw2**2, x=nevan_mesh_arr[1:-1])
return diff
# perform Hardy optimization
G_w_opt_res = []
for eta in eta_list:
print(f'eta={eta}')
G_w_opt_res.append(solver.optimize(nevan_mesh, eta, SmoothnessOptimization, gtol=1e-4, maxiter=1000, verbose=True))
print('-------')
eta=0.4
Optimization terminated successfully.
The optimizer performed 0 functions evaluations.
-------
eta=0.2
Optimization terminated successfully.
The optimizer performed 12 functions evaluations.
-------
eta=0.1
Optimization terminated successfully.
The optimizer performed 115 functions evaluations.
-------
Now we plot all results in one plot:
[16]:
fig, ax = plt.subplots(1,1, figsize=(8,4), dpi=150, squeeze=False)
ax = ax.reshape(-1)
ax[0].oplot(-1*G_loc_w[0,0].imag, lw=5, label=f'ref')
ax[0].oplot(-1*G_w_non_opt.imag, label=f'Nev')
for eta, G_w_opt in zip(eta_list, G_w_opt_res):
ax[0].oplot(-1*G_w_opt.imag, label=fr'$\eta=${eta}')
ax[0].set_xlim(w_window)
ax[0].set_xlabel(r'$\omega$')
ax[0].set_ylabel(r'Im $G (\omega)$ (eV)')
plt.legend()
plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
It becomes clear that the larger eta values will introduce more smoothness. However, this goes hand in hand with a lost in features. An optimal eta has to be choosen carefully.