Source code for dmft_tools.greens_functions_mixer

# -*- 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.gf import MeshImFreq, MeshImTime, Gf, make_gf_from_fourier
from triqs.gf.descriptors import Fourier
from triqs.gf.tools import inverse

from . import legendre_filter


[docs] def _bgf_to_vec(bgf, deg_shells=None): """ flattens a BlockGf to a 1d numpy array if deg_shells is given, only non-deg shells are given back """ if deg_shells: deg_vecs = [] for deg_block in deg_shells: deg_vecs.append(bgf[deg_block[0]]) vec = np.hstack(gf.data.flatten() for gf in deg_vecs) else: vec = np.hstack(bgf[bl].data.flatten() for bl in bgf.indices) return vec
def _vec_to_bgf(vec, bgf, deg_shells=None): G = bgf.copy() if deg_shells: start = 0 for deg_block in deg_shells: # only the first element was stored in numpy vec bl = deg_block[0] points = len([w.value for w in G[bl].mesh]) * np.prod(list(G[bl].target_shape)) end = start+points G[bl].data[:,:,:] = vec[start:end].reshape(G[bl].data.shape) start = end for other_blocks in deg_block[1:]: G[other_blocks] << G[bl] else: start = 0 for i, bl in enumerate(G.indices): points = len([w.value for w in G[bl].mesh]) * np.prod(list(G[bl].target_shape)) end = start+points G[bl].data[:,:,:] = vec[start:end].reshape(G[bl].data.shape) start = end return G
[docs] def _broyden_update(it, 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 ---------- it : int iteration number broyler: dict dict containing the broyden variables 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() # if this is the first time broyden update is called (it==2) # store init G0 first if it == 2: # 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') # define m as it-1-1 (last minus for idx starting at 0) m = it-2 # 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 ) # 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, it, 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(it, 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