# -*- coding: utf-8 -*-
################################################################################
#
# solid_dmft - A versatile python wrapper to perform DFT+DMFT calculations
# utilizing the TRIQS software library
#
# Copyright (C) 2018-2020, ETH Zurich
# Copyright (C) 2021, The Simons Foundation
# authors: A. Hampel, M. Merkel, and S. Beck
#
# solid_dmft 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.
#
# solid_dmft 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
# solid_dmft (in the file COPYING.txt in this directory). If not, see
# <http://www.gnu.org/licenses/>.
#
################################################################################
# system
import numpy as np
# triqs
import triqs.utility.mpi as mpi
from h5 import HDFArchive
from triqs.gfs import MeshImFreq, MeshImTime, Gf, make_gf_from_fourier
from triqs.gfs.descriptors import Fourier
from triqs.gfs.tools import inverse
from . import legendre_filter
[docs]
def _mixed_blocks(bgf, deg_shells):
"""
groups of block names whose first entry is stored in the Broyden vector;
the other entries of a group are copies of the first. Blocks that are in
no degenerate shell form a group of their own.
"""
groups = [list(deg_block) for deg_block in deg_shells] if deg_shells else []
in_shell = {bl for group in groups for bl in group}
groups += [[bl] for bl in bgf.indices if bl not in in_shell]
return groups
[docs]
def _bgf_to_vec(bgf, deg_shells=None):
"""
flattens a BlockGf to a 1d numpy array
if deg_shells is given, only the first block of each degenerate shell is stored
"""
return np.hstack([bgf[group[0]].data.flatten() for group in _mixed_blocks(bgf, deg_shells)])
def _vec_to_bgf(vec, bgf, deg_shells=None):
G = bgf.copy()
start = 0
for group in _mixed_blocks(G, deg_shells):
bl = group[0]
end = start + G[bl].data.size
G[bl].data[...] = vec[start:end].reshape(G[bl].data.shape)
start = end
for other_block in group[1:]:
G[other_block] << G[bl]
return G
[docs]
def _broyden_update(broyler, general_params, deg_shells, G0_freq, G0_freq_previous):
"""
calculates the broyden update of G0 using the algorithm presented in:
doi.org/10.1103/PhysRevB.80.125125 (Rok Zitko) &
doi.org/10.1103/PhysRevB.38.12807 (D.D. Johnson)
optimize function to find roots: F(G0)=dmft_step(G0) - G0 = 0
where dmft_step gives back the normal defined G0=inverse(Gloc^-1 - Sigma)
Parameters
----------
broyler: dict
dict containing the broyden variables; an empty history starts a new
Broyden sequence with a linear mixing step
general_paramters: general params dict
Returns
-------
"""
alpha = general_params['g0_mix']
# hard code w_0
w_0 = 0.01
G0_broyden_update = G0_freq.copy()
# build here leg compression
if type(G0_freq.mesh) is MeshImFreq and 'n_l' in general_params:
G0 = legendre_filter.apply(make_gf_from_fourier(G0_freq), order=general_params['n_l'] )
G0_prev = legendre_filter.apply(make_gf_from_fourier(G0_freq_previous), order=general_params['n_l'] )
else:
G0 = G0_freq.copy()
G0_prev = G0_freq_previous.copy()
# empty history: store the initial G0 first
if not broyler['V']:
# V holds the G0 as 1d numpy vectors
broyler['V'].append(_bgf_to_vec(G0_prev, deg_shells))
# F = dmft_step_G0 - V[-1]
broyler['F'].append(_bgf_to_vec(G0, deg_shells)-broyler['V'][-1])
# for the second iteration the broyden vec corresponds to simple linear mixing with alpha
broyden_vec = np.zeros(broyler['V'][0].shape[0],dtype='complex')
# for all other iteration calc V_mp1 = V[-1] + alpha * F[-1] - broyden_vec
# broyden vec is sum(n=1 to it-1) gamma_m_n U^n
# gamma_m_n = sum(k=1 to it-1) c_k^m beta_k_n^m
# gamma_m_n is a number for each iteration and U^n is of dim G0 for each iteration
else:
# get size of broyden update vector
dim = broyler['V'][0].shape[0]
broyden_vec = np.zeros(dim,dtype='complex')
# update broyden input params
broyler['F'].append(_bgf_to_vec(G0, deg_shells)-broyler['V'][-1])
norm_F = np.linalg.norm(broyler['F'][-1]-broyler['F'][-2])
print('broyden norm: {:.3f}'.format(norm_F))
# dF is normed differences between last two functions
broyler['dF'].append( (broyler['F'][-1]-broyler['F'][-2])/norm_F )
# dV the normed difference between the two last input G0
broyler['dV'].append( (broyler['V'][-1] - broyler['V'][-2])/norm_F )
m = len(broyler['dF'])
# now build broyden correction vector
#only keep last m iterations
mem_iter = general_params['broy_max_it']-1
start = 0
mat_dim = m
if m > mem_iter and not general_params['broy_max_it'] == -1:
for n in range(0, m-mem_iter):
start += 1
mat_dim -= 1
# first create beta_m
A = np.zeros((mat_dim,mat_dim),dtype='complex')
gamma = np.zeros(mat_dim,dtype='complex')
for k_i in range(start,m):
for n_i in range(start,m):
A[k_i-start,n_i-start] = np.dot(broyler['dF'][n_i].conjugate().transpose(),broyler['dF'][k_i])
if n_i == k_i:
A[k_i-start,n_i-start] += w_0**2
beta_m = np.linalg.inv(A)
for n_i in range(start,m):
U_n = alpha * broyler['dF'][n_i] + broyler['dV'][n_i]
# construct gamma
for k_i in range(start,m):
c_k = np.dot(broyler['dF'][k_i].conjugate().transpose(),broyler['F'][-1])
gamma[n_i-start] += c_k*beta_m[k_i-start,n_i-start]
broyden_vec += gamma[n_i-start] * U_n
# update input G0
V_mp1 = broyler['V'][-1] + alpha * broyler['F'][-1] - broyden_vec
broyler['V'].append(V_mp1)
# convert Gl back to G0_freq
if type(G0_freq.mesh) is MeshImFreq and 'n_l' in general_params:
G0_l_update = _vec_to_bgf(broyler['V'][-1], G0, deg_shells)
for block, g in G0_l_update:
g.enforce_discontinuity(np.identity(g.target_shape[0]))
G0_broyden_update[block].set_from_legendre(g)
else:
G0_broyden_update << _vec_to_bgf(broyler['V'][-1], G0, deg_shells)
return G0_broyden_update, broyler
def mix_g0(solver, general_params, icrsh, archive, G0_freq_previous, deg_shell):
if general_params['g0_mix_type'] == 'linear' and general_params['g0_mix'] < 1.0:
mpi.report('\nlinear mixing G0 with previous iteration by factor {:.3f}\n'.format(general_params['g0_mix']))
solver.G0_freq << (general_params['g0_mix'] * solver.G0_freq
+ (1-general_params['g0_mix']) * G0_freq_previous)
elif general_params['g0_mix_type'] == 'broyden':
mpi.report('\nbroyden mixing G0 with previous iteration with alpha {:.3f}'.format(general_params['g0_mix']))
mpi.report('\n############\n!!!! WARNING !!!! broyden mixing is still in early testing stage ! Use with caution.\n############\n')
# TODO implement broyden mixing for Sigma
if mpi.is_master_node():
with HDFArchive(archive, 'a') as ar:
broyler = ar['DMFT_results']['broyler']
# calculate the next G0 via broyden scheme
G0_broyden_update, broyler[icrsh] = _broyden_update(broyler[icrsh], general_params,
deg_shell, solver.G0_freq,
G0_freq_previous)
# store broyden update to h5 archive
ar['DMFT_results']['broyler'] = broyler
solver.G0_freq << G0_broyden_update
solver.G0_freq << mpi.bcast(solver.G0_freq)
return solver
def mix_sigma(general_params, n_inequiv_shells, solvers, Sigma_freq_previous):
if general_params['sigma_mix'] < 1.0 and mpi.is_master_node():
print('mixing sigma with previous iteration by factor {:.3f}\n'.format(general_params['sigma_mix']))
for icrsh in range(n_inequiv_shells):
solvers[icrsh].Sigma_freq << (general_params['sigma_mix'] * solvers[icrsh].Sigma_freq
+ (1-general_params['sigma_mix']) * Sigma_freq_previous[icrsh])
for icrsh in range(n_inequiv_shells):
solvers[icrsh].Sigma_freq << mpi.bcast(solvers[icrsh].Sigma_freq)
return solvers