From f803c13285c3fa38e7d67f30d8191ee775458917 Mon Sep 17 00:00:00 2001 From: Priyanka Seth Date: Mon, 22 Sep 2014 19:22:17 +0200 Subject: [PATCH] New files and changes * hamiltonians contains functions to create basic hamiltonians * Us are now computed by a compact python script * solver_multiband sets up the local problem for LDA+DMFT --- python/U_matrix.py | 287 +++++++++++++++------- python/__init__.py | 12 +- python/hamiltonians.py | 127 ++++++++++ python/solver_multiband.py | 476 +++---------------------------------- 4 files changed, 367 insertions(+), 535 deletions(-) create mode 100644 python/hamiltonians.py diff --git a/python/U_matrix.py b/python/U_matrix.py index a2500c34..846e7222 100644 --- a/python/U_matrix.py +++ b/python/U_matrix.py @@ -1,100 +1,213 @@ - -################################################################################ -# -# TRIQS: a Toolbox for Research in Interacting Quantum Systems -# -# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola -# -# TRIQS 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. -# -# TRIQS 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 -# TRIQS. If not, see . -# -################################################################################ - -# calculates the four index U matrix - -import numpy -from types import * from math import sqrt -import copy -from vertex import u4ind -#from pytriqs.applications.dft.vertex import u4ind - -class Umatrix: - """calculates, stores, and manipulates the four index U matrix""" - - def __init__(self, l, U_interact=0, J_hund=0): - - self.l = l - self.U_av = U_interact - self.J = J_hund - - self.N = 2*l+1 # multiplicity - - #self.Ucmplx = numpy.zeros([self.N,self.N,self.N,self.N],numpy.float_) - #self.Ucubic = numpy.zeros([self.N,self.N,self.N,self.N],numpy.float_) +from scipy.misc import factorial as fact +from itertools import product +import numpy as np - def __call__(self, T = None, rcl = None): - """calculates the four index matrix. Slater parameters can be provided in rcl, - and a transformation matrix from complex harmonics to a specified other representation (e.g. cubic). - If T is not given, use standard complex harmonics.""" +# The interaction matrix in desired basis +# U^{spherical}_{m1 m2 m3 m4} = \sum_{k=0}^{2l} F_k angular_matrix_element(l, k, m1, m2, m3, m4) +def U_matrix(l, radial_integrals=None, U_int=None, J_hund=None, basis="spherical", T=None): + """Calculate U matrix being given either radial_integrals or U_int and J_hund. + l = angular momentum of shell being treated (l=2 for d shell, l=3 for f shell) + radial_integrals = [F0,F2,F4,..] (default None) + U_int, J_hund = values to use to compute radial_integrals (default None), + basis = "spherical", "cubic", or "other", + T = transformation matrix from spherical to desired basis, if basis='other' (default None)""" - if rcl is None: rcl = self.get_rcl(self.U_av,self.J,self.l) + # Check all necessary information is present and consistent + if radial_integrals == None and (U_int == None and J_hund == None): + raise ValueError("U_matrix: provide either the radial_integrals or U_int and J_hund.") + if radial_integrals == None and (U_int != None and J_hund != None): + radial_integrals = U_J_to_radial_integrals(l, U_int, J_hund) + if radial_integrals != None and (U_int != None and J_hund != None): + if len(radial_integrals)-1 != l: + raise ValueError("U_matrix: inconsistency in l and number of radial_integrals provided.") + if (radial_integrals - U_J_to_radial_integrals(l, U_int, J_hund)).any() != 0.0: + print "Warning: U_matrix: radial_integrals provided do not match U_int and J_hund. Using radial_integrals to calculate spherical U_matrix." - if (T is None): - TM = numpy.identity(self.N,numpy.complex_) + # Full interaction matrix + # Basis of spherical harmonics Y_{-2}, Y_{-1}, Y_{0}, Y_{1}, Y_{2} + # U^{spherical}_{m1 m2 m3 m4} = \sum_{k=0}^{2l} F_k angular_matrix_element(l, k, m1, m2, m3, m4) + U_matrix = np.zeros((2*l+1,2*l+1,2*l+1,2*l+1),dtype=float) + + m_range = range(-l,l+1) + for n, F in enumerate(radial_integrals): + k = 2*n + for m1, m2, m3, m4 in product(m_range,m_range,m_range,m_range): + U_matrix[m1+l,m2+l,m3+l,m4+l] += F * angular_matrix_element(l,k,m1,m2,m3,m4) + + # Transform from spherical basis if needed + if basis == "cubic": T = spherical_to_cubic(l) + if basis == "other" and T == None: + raise ValueError("U_matrix: provide T for other bases.") + if T != None: transform_U_matrix(U_matrix, T) + + return U_matrix + +# Convert full 4-index U matrix to 2-index density-density form +def reduce_4index_to_2index(U_4index): + """Reduces the four-index matrix to two-index matrices.""" + + size = len(U_4index) # 2l+1 + U = np.zeros((size,size),dtype=float) # matrix for same spin + Uprime = np.zeros((size,size),dtype=float) # matrix for opposite spin + + m_range = range(size) + for m,mp in product(m_range,m_range): + U[m,mp] = U_4index[m,mp,m,mp].real - U_4index[m,mp,mp,m].real + Uprime[m,mp] = U_4index[m,mp,m,mp].real + + return U, Uprime + +# Construct the 2-index matrices for the density-density form +def U_matrix_kanamori(n_orb, U_int=None, J_hund=None): + """Calculate the Kanamori U and Uprime matrices.""" + + U = np.zeros((n_orb,n_orb),dtype=float) # matrix for same spin + Uprime = np.zeros((n_orb,n_orb),dtype=float) # matrix for opposite spin + + m_range = range(n_orb) + for m,mp in product(m_range,m_range): + if m == mp: + Uprime[m,mp] = U_int else: - TM = T - self.Nmat = len(TM) + U[m,mp] = U_int - 3.0*J_hund + Uprime[m,mp] = U_int - 2.0*J_hund - self.Ufull = u4ind(rcl,TM) - - - - def reduce_matrix(self): - """Reduces the four-index matrix to two-index matrices.""" + return U, Uprime - if (self.N==self.Nmat): - self.U = numpy.zeros([self.N,self.N],numpy.float_) # matrix for same spin - self.Up = numpy.zeros([self.N,self.N],numpy.float_) # matrix for opposite spin +# Get t2g or eg components +def t2g_submatrix(U): + """Return only the t2g part of the full d-manifold two- or four-index U matrix.""" + return subarray(U, len(U.shape)*[(0,1,3)]) - for m in range(self.N): - for mp in range(self.N): - self.U[m,mp] = self.Ufull[m,mp,m,mp].real - self.Ufull[m,mp,mp,m].real - self.Up[m,mp] = self.Ufull[m,mp,m,mp].real - else: - self.U = numpy.zeros([self.Nmat,self.Nmat],numpy.float_) # matrix - for m in range(self.Nmat): - for mp in range(self.Nmat): - self.U[m,mp] = self.Ufull[m,mp,m,mp].real - self.Ufull[m,mp,mp,m].real - - - +def eg_submatrix(U): + """Return only the eg part of the full d-manifold two- or four-index U matrix.""" + return subarray(U, len(U.shape)*[(2,4)]) +# Transform the interaction matrix into another basis +def transform_U_matrix(U_matrix, T): + """Transform the interaction matrix into another basis by applying matrix T.""" + return np.einsum("ij,kl,jlmo,mn,op",np.conj(T),np.conj(T),U_matrix,np.transpose(T),np.transpose(T)) - def get_rcl(self, U_int, J_hund, l): +# Rotation matrices: complex harmonics to cubic harmonics +# Complex harmonics basis: ..., Y_{-2}, Y_{-1}, Y_{0}, Y_{1}, Y_{2}, ... +def spherical_to_cubic(l): + """Returns the spherical harmonics to cubic harmonics rotation matrix.""" + size = 2*l+1 + T = np.zeros((size,size),dtype=complex) + if l == 0: + cubic_names = ("s") + elif l == 1: + cubic_names = ("x","y","z") + T[0,0] = 1.0/sqrt(2); T[0,2] = -1.0/sqrt(2) + T[1,0] = 1j/sqrt(2); T[1,2] = 1j/sqrt(2) + T[2,1] = 1.0 + elif l == 2: + cubic_names = ("xy","yz","z^2","xz","x^2-y^2") + T[0,0] = 1j/sqrt(2); T[0,4] = -1j/sqrt(2) + T[1,1] = 1j/sqrt(2); T[1,3] = 1j/sqrt(2) + T[2,2] = 1.0 + T[3,1] = 1.0/sqrt(2); T[3,3] = -1.0/sqrt(2) + T[4,0] = 1.0/sqrt(2); T[4,4] = 1.0/sqrt(2) + elif l == 3: + cubic_names = ("x(x^2-3y^2)","z(x^2-y^2)","xz^2","z^3","yz^2","xyz","y(3x^2-y^2)") + T[0,0] = 1.0/sqrt(2); T[0,6] = -1.0/sqrt(2) + T[1,1] = 1.0/sqrt(2); T[1,5] = 1.0/sqrt(2) + T[2,2] = 1.0/sqrt(2); T[2,4] = -1.0/sqrt(2) + T[3,5] = 1.0 + T[4,2] = 1j/sqrt(2); T[4,4] = 1j/sqrt(2) + T[5,1] = 1j/sqrt(2); T[5,5] = -1j/sqrt(2) + T[6,0] = 1j/sqrt(2); T[6,6] = 1j/sqrt(2) + else: raise ValueError("spherical_to_cubic: implemented only for l=0,1,2,3") - #rcl = numpy.array([0.0, 0.0, 0.0, 0.0],numpy.float_) - xx = l+1 - rcl = numpy.zeros([xx],numpy.float_) - if(l==2): - rcl[0] = U_int - rcl[1] = J_hund * 14.0 / (1.0 + 0.63) - rcl[2] = 0.630 * rcl[1] - elif(l==3): - rcl[0] = U_int - rcl[1] = 6435.0 * J_hund / (286.0 + 195.0 * 451.0 / 675.0 + 250.0 * 1001.0 / 2025.0) - rcl[2] = 451.0 * rcl[1] / 675.0 - rcl[3] = 1001.0 * rcl[1] / 2025.0 + return T - return rcl +# Names of cubic harmonics +def cubic_names(l): + if l == 0 or l == 's': + return ("s") + elif l == 1 or l == 'p': + return ("x","y","z") + elif l == 2 or l == 'd': + return ("xy","yz","z^2","xz","x^2-y^2") + elif l == 't2g': + return ("xy","yz","xz") + elif l == 'eg': + return ("z^2","x^2-y^2") + elif l == 3 or l == 'f': + return ("x(x^2-3y^2)","z(x^2-y^2)","xz^2","z^3","yz^2","xyz","y(3x^2-y^2)") + else: raise ValueError("cubic_names: implemented only for l=0,1,2,3") + +# Convert U,J -> radial integrals F_k +def U_J_to_radial_integrals(l, U_int, J_hund): + """Determines the radial integrals F_k from U_int and J_hund.""" + + F = np.zeros((l+1),dtype=float) + if l == 2: + F[0] = U_int + F[1] = J_hund * 14.0 / (1.0 + 0.63) + F[2] = 0.630 * F[1] + elif l == 3: + F[0] = U_int + F[1] = 6435.0 * J_hund / (286.0 + 195.0 * 451.0 / 675.0 + 250.0 * 1001.0 / 2025.0) + F[2] = 451.0 * F[1] / 675.0 + F[3] = 1001.0 * F[1] / 2025.0 + else: raise ValueError("U_J_to_radial_integrals: implemented only for l=2,3") + + return F + +# Angular matrix elements of particle-particle interaction +# (2l+1)^2 ((l 0) (k 0) (l 0))^2 \sum_{q=-k}^{k} (-1)^{m1+m2+q} ((l -m1) (k q) (l m3)) ((l -m2) (k -q) (l m4)) +def angular_matrix_element(l, k, m1, m2, m3, m4): + result = 0 + for q in range(-k,k+1): + result += three_j_symbol((l,-m1),(k,q),(l,m3))*three_j_symbol((l,-m2),(k,-q),(l,m4))*(-1.0 if (m1+q+m2) % 2 else 1.0) + result *= (2*l+1)**2 * (three_j_symbol((l,0),(k,0),(l,0))**2) + return result + +# Wigner 3-j symbols +# ((j1 m1) (j2 m2) (j3 m3)) +def three_j_symbol(jm1, jm2, jm3): + j1, m1 = jm1 + j2, m2 = jm2 + j3, m3 = jm3 + + if (m1+m2+m3 != 0 or + m1 < -j1 or m1 > j1 or + m2 < -j2 or m2 > j2 or + m3 < -j3 or m3 > j3 or + j3 > j1 + j2 or + j3 < abs(j1-j2)): + return .0 + + result = -1.0 if (j1-j2-m3) % 2 else 1.0 + result *= sqrt(fact(j1+j2-j3)*fact(j1-j2+j3)*fact(-j1+j2+j3)/fact(j1+j2+j3+1)) + result *= sqrt(fact(j1-m1)*fact(j1+m1)*fact(j2-m2)*fact(j2+m2)*fact(j3-m3)*fact(j3+m3)) + + t_min = max(j2-j3-m1,j1-j3+m2,0) + t_max = min(j1-m1,j2+m2,j1+j2-j3) + + t_sum = 0 + for t in range(t_min,t_max+1): + t_sum += (-1.0 if t % 2 else 1.0)/(fact(t)*fact(j3-j2+m1+t)*fact(j3-j1-m2+t)*fact(j1+j2-j3-t)*fact(j1-m1-t)*fact(j2+m2-t)) + + result *= t_sum + return result + +# Clebsch-Gordan coefficients +# < j1 m1 j2 m2 | j3 m3 > = (-1)^{j1-j2+m3} \sqrt{2j3+1} ((j1 m1) (j2 m2) (j3 -m3)) +def clebsch_gordan(jm1, jm2, jm3): + norm = sqrt(2*jm3[0]+1)*(-1 if jm1[0]-jm2[0]+jm3[1] % 2 else 1) + return norm*three_j_symbol(jm1,jm2,(jm3[0],-jm3[1])) + +# Create subarray containing columns in idxlist +# e.g. idxlist = [(0),(2,3),(0,1,2,3)] gives +# column 0 for 1st dim, +# columns 2 and 3 for 2nd dim, +# columns 0,1,2 and 3 for 3rd dim. +#def subarray(a,idxlist,n=len(a.shape)-1) : +def subarray(a,idxlist,n=None) : + if n == None: n = len(a.shape)-1 + sa = a[tuple(slice(x) for x in a.shape[:n]) + (idxlist[n],)] + return subarray(sa,idxlist, n-1) if n > 0 else sa diff --git a/python/__init__.py b/python/__init__.py index 757ecb23..3d2edeb9 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -23,9 +23,13 @@ from sumk_lda import SumkLDA from symmetry import Symmetry from sumk_lda_tools import SumkLDATools -from U_matrix import Umatrix +from U_matrix import * from converters import * +from hamiltonians import * -__all__ =['SumkLDA','Symmetry','SumkLDATools','Umatrix','Wien2kConverter','HkConverter'] - - +__all__=['SumkLDA','Symmetry','SumkLDATools','Wien2kConverter','HkConverter', +'U_J_to_radial_integrals', 'U_matrix', 'U_matrix_kanamori', +'angular_matrix_element', 'clebsch_gordan', 'cubic_names', 'eg_submatrix', +'reduce_4index_to_2index', 'spherical_to_cubic', 't2g_submatrix', +'three_j_symbol', 'transform_U_matrix', +'h_loc_slater','h_loc_kanamori','h_loc_density','get_mkind'] diff --git a/python/hamiltonians.py b/python/hamiltonians.py new file mode 100644 index 00000000..f24708a4 --- /dev/null +++ b/python/hamiltonians.py @@ -0,0 +1,127 @@ +from pytriqs.operators.operators2 import * +from itertools import product + +# Define commonly-used Hamiltonians here: Slater, Kanamori, density-density + +def h_loc_slater(spin_names,orb_names,orb_hyb,U_matrix,H_dump=None): + + if H_dump: + H_dump_file = open(H_dump,'w') + H = Operator() + mkind = get_mkind(orb_hyb) + for s1, s2 in product(spin_names,spin_names): + for a1, a2, a3, a4 in product(orb_names,orb_names,orb_names,orb_names): + U_val = U_matrix[orb_names.index(a1),orb_names.index(a2),orb_names.index(a3),orb_names.index(a4)] + if abs(U_val.imag) > 1e-10: + raise RuntimeError("Matrix elements of U are not real. Are you using a cubic basis?") + + H_term = 0.5 * U_val.real * c_dag(*mkind(s1,a1)) * c_dag(*mkind(s2,a2)) * c(*mkind(s2,a4)) * c(*mkind(s1,a3)) + H += H_term + + # Dump terms of H + if H_dump and not H_term.is_zero(): + H_dump_file.write(mkind(s1,a1)[0] + '\t') + H_dump_file.write(mkind(s2,a2)[0] + '\t') + H_dump_file.write(mkind(s2,a3)[0] + '\t') + H_dump_file.write(mkind(s1,a4)[0] + '\t') + H_dump_file.write(str(U_val.real) + '\n') + + return H + +def h_loc_kanamori(spin_names,orb_names,orb_hyb,U,Uprime,J_hund,H_dump=None): + + if H_dump: + H_dump_file = open(H_dump,'w') + H = Operator() + mkind = get_mkind(orb_hyb) + + # density terms: + for s1, s2 in product(spin_names,spin_names): + for a1, a2 in product(orb_names,orb_names): + if (s1==s2): + U_val = U[orb_names.index(a1),orb_names.index(a2)] + else: + U_val = Uprime[orb_names.index(a1),orb_names.index(a2)] + + H_term = 0.5 * U_val * n(*mkind(s1,a1)) * n(*mkind(s2,a2)) + H += H_term + + # Dump terms of H + if H_dump and not H_term.is_zero(): + H_dump_file.write("Density-density terms" + '\n') + H_dump_file.write(mkind(s1,a1)[0] + '\t') + H_dump_file.write(mkind(s2,a2)[0] + '\t') + H_dump_file.write(str(U_val) + '\n') + + # spin-flip terms: + for s1, s2 in product(spin_names,spin_names): + if (s1==s2): + continue + for a1, a2 in product(orb_names,orb_names): + if (a1==a2): + continue + H_term = -0.5 * J_hund * c_dag(*mkind(s1,a1)) * c(*mkind(s2,a1)) * c_dag(*mkind(s2,a2)) * c(*mkind(s1,a2)) + H += H_term + + # Dump terms of H + if H_dump and not H_term.is_zero(): + H_dump_file.write("Spin-flip terms" + '\n') + H_dump_file.write(mkind(s1,a1)[0] + '\t') + H_dump_file.write(mkind(s2,a2)[0] + '\t') + H_dump_file.write(mkind(s2,a3)[0] + '\t') + H_dump_file.write(mkind(s1,a4)[0] + '\t') + H_dump_file.write(str(-J_hund) + '\n') + + # pair-hopping terms: + for s1, s2 in product(spin_names,spin_names): + if (s1==s2): + continue + for a1, a2 in product(orb_names,orb_names): + if (a1==a2): + continue + H_term = 0.5 * J_hund * c_dag(*mkind(s1,a1)) * c_dag(*mkind(s2,a1)) * c(*mkind(s2,a2)) * c(*mkind(s1,a2)) + H += H_term + + # Dump terms of H + if H_dump and not H_term.is_zero(): + H_dump_file.write("Pair-hopping terms" + '\n') + H_dump_file.write(mkind(s1,a1)[0] + '\t') + H_dump_file.write(mkind(s2,a2)[0] + '\t') + H_dump_file.write(mkind(s2,a3)[0] + '\t') + H_dump_file.write(mkind(s1,a4)[0] + '\t') + H_dump_file.write(str(-J_hund) + '\n') + + return H + +def h_loc_density(spin_names,orb_names,orb_hyb,U,Uprime,H_dump=None): + + if H_dump: + H_dump_file = open(H_dump,'w') + H = Operator() + mkind = get_mkind(orb_hyb) + for s1, s2 in product(spin_names,spin_names): + for a1, a2 in product(orb_names,orb_names): + if (s1==s2): + U_val = U[orb_names.index(a1),orb_names.index(a2)] + else: + U_val = Uprime[orb_names.index(a1),orb_names.index(a2)] + + H_term = 0.5 * U_val * n(*mkind(s1,a1)) * n(*mkind(s2,a2)) + H += H_term + + # Dump terms of H + if H_dump and not H_term.is_zero(): + H_dump_file.write(mkind(s1,a1)[0] + '\t') + H_dump_file.write(mkind(s2,a2)[0] + '\t') + H_dump_file.write(str(U_val) + '\n') + + return H + +# Set function to make index for GF blocks given spin sn and orbital name on +def get_mkind(orb_hyb): + if orb_hyb: + mkind = lambda sn, on: (sn, on) + else: + mkind = lambda sn, on: (sn+'_'+on, 0) + + return mkind diff --git a/python/solver_multiband.py b/python/solver_multiband.py index e9605b97..8edafd33 100644 --- a/python/solver_multiband.py +++ b/python/solver_multiband.py @@ -1,453 +1,41 @@ +from itertools import product +from hamiltonians import * -################################################################################ -# -# TRIQS: a Toolbox for Research in Interacting Quantum Systems -# -# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola -# -# TRIQS 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. -# -# TRIQS 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 -# TRIQS. If not, see . -# -################################################################################ +class LocalProblem(): + def __init__(self,spin_names,orb_names,orb_hyb,h_loc_type,**h_loc_params): -from pytriqs.operators import * -from pytriqs.applications.impurity_solvers.cthyb_matrix import Solver -from pytriqs.applications.dft.U_matrix import * + self.spin_names = spin_names + self.orb_names = orb_names + self.orb_hyb = orb_hyb + self.h_loc_type = h_loc_type + self.h_loc_params = h_loc_params -import pytriqs.utility.mpi as mpi -from types import * -import numpy + self.gf_struct = self.get_gf_struct() + self.h_loc = self.get_h_loc() -def sum_list(L): - """ Can sum any list""" - return reduce(lambda x, y: x+y, L) if len(L)>0 else [] + # Set block structure of GF + def get_gf_struct(self): + gf_struct = {} + if self.orb_hyb: # outer blocks are spin blocks + for sn in self.spin_names: + gf_struct[sn] = [int(i) for i in self.orb_names] + else: # outer blocks are spin-orbital blocks + for sn, an in product(self.spin_names,self.orb_names): + gf_struct[sn+'_'+an] = [0] -######################################### -# -# Solver for the Multi-Band problem -# -######################################### + return gf_struct + # Pick desired Hamiltonian + def get_h_loc(self): -class SolverMultiBand(Solver): - """ - This is a general solver for a multiband local Hamiltonian. - Calling arguments: - beta = inverse temperature - n_orb = Number of local orbitals - U_interact = Average Coulomb interaction - J_hund = Hund coupling - use_spinflip = true/false - use_pairhop = true/false - use_matrix: Use the interaction matrix calculated from the Slater integrals - is use_matrix, you need also: - l: angular momentum of the orbital, l=2 is d - T: Transformation matrix for U vertex. If not present, use standard complex harmonics - - """ - - def __init__(self, beta, n_orb, gf_struct = False, map = False, omega_max = None): - - self.n_orb = n_orb - - # either get or construct gf_struct - if (gf_struct): - assert map, "give also the mapping!" - self.map = map + if self.h_loc_type == "slater": + return h_loc_slater(self.spin_names,self.orb_names,self.orb_hyb,**self.h_loc_params) # h_loc_params must include U_matrix, and optionally H_dump + elif self.h_loc_type == "kanamori": + return h_loc_kanamori(self.spin_names,self.orb_names,self.orb_hyb,**self.h_loc_params) # h_loc_params must include U, Uprime, J_hund, and optionally H_dump + elif self.h_loc_type == "density": + return h_loc_density(self.spin_names,self.orb_names,self.orb_hyb,**self.h_loc_params) # h_loc_params must include U, Uprime, and optionally H_dump + elif self.h_loc_type == "other": + return self.h_loc_params["h_loc"] # user provides h_loc with argument h_loc else: - # standard gf_struct and map - gf_struct = [ ('%s'%(ud),[n for n in range(n_orb)]) for ud in ['up','down'] ] - self.map = {'up' : ['up' for v in range(n_orb)], 'down' : ['down' for v in range(n_orb)]} - - # now initialize the solver with the stuff given above: - if (omega_max is None): - Solver.__init__(self, beta = beta, gf_struct = gf_struct) - else: - n_w = int(omega_max*beta/(2.*numpy.pi)) - Solver.__init__(self, beta = beta, gf_struct = gf_struct, n_w = n_w) - - - - def solve(self, U_interact=None, J_hund=None, use_spinflip=False, - use_matrix = True, l=2, T=None, dim_reps=None, irep=None, deg_orbs = [], sl_int = None, **params): - - self.use_spinflip = use_spinflip - self.U, self.Up, self.U4ind, self.offset = set_U_matrix(U_interact,J_hund,self.n_orb,l,use_matrix,T,sl_int,use_spinflip,dim_reps,irep) - - # define mapping of indices: - self.map_ind={} - for nm in self.map: - bl_names = self.map[nm] - block = [] - for a,al in self.gf_struct: - if a in bl_names: block.append(al) - - self.map_ind[nm] = range(self.n_orb) - i = 0 - for al in block: - cnt = 0 - for a in range(len(al)): - self.map_ind[nm][i] = cnt - i = i+1 - cnt = cnt+1 - - # set the Hamiltonian - if (use_spinflip==False): - Hamiltonian = self.__set_hamiltonian_density() - else: - if (use_matrix): - #Hamiltonian = self.__set_full_hamiltonian_slater() - Hamiltonian = self.__set_spinflip_hamiltonian_slater() - else: - Hamiltonian = self.__set_full_hamiltonian_kanamori(J_hund = J_hund) - - # set the Quantum numbers - Quantum_Numbers = self.__set_quantum_numbers(self.gf_struct) - - # Determine if there are only blocs of size 1: - self.blocssizeone = True - for ib in self.gf_struct: - if (len(ib[1])>1): self.blocssizeone = False - - nc = params.pop("n_cycles",10000) - if ((self.blocssizeone) and (self.use_spinflip==False)): - use_seg = True - else: - use_seg = False - #gm = self.set_global_moves(deg_orbs) - - Solver.solve(self,H_local = Hamiltonian, quantum_numbers = Quantum_Numbers, n_cycles = nc, use_segment_picture = use_seg, **params) - - - def set_global_moves(self, deg_orbs, factor=0.05): - # Sets some global moves given orbital degeneracies: - - strbl = '' - strind = '' - inddone = [] - - for orbs in deg_orbs: - ln = len(orbs) - orbsorted = sorted(orbs) - for ii in range(ln): - if (strbl!=''): strbl += ',' - bl1 = orbsorted[ii] - bl2 = orbsorted[(ii+1)%ln] - ind1 = [ll for ll in self.Sigma[bl1].indices ] - ind2 = [ll for ll in self.Sigma[bl2].indices ] - - strbl += "'" + bl1 + "':'" + bl2 + "'" - for kk, ind in enumerate(ind1): - if not (ind in inddone): - if (strind!=''): strind += ',' - strind += '%s:%s'%(ind1[kk],ind2[kk]) - inddone.append(ind) - - - if len(deg_orbs)>0: - str = 'Global_Moves = [ (%s, lambda (a,alpha,dag) : ({ '%factor + strbl + ' }[a], {' + strind + '}[alpha], dag) )]' - exec str - return Global_Moves - else: - return [] - - - def __set_hamiltonian_density(self): - # density-density Hamiltonian: - - spinblocs = [v for v in self.map] - #print spinblocs - Hamiltonian = N(self.map[spinblocs[0]][0],0) # initialize it - - for sp1 in spinblocs: - for sp2 in spinblocs: - for i in range(self.n_orb): - for j in range(self.n_orb): - if (sp1==sp2): - Hamiltonian += 0.5 * self.U[self.offset+i,self.offset+j] * N(self.map[sp1][i],self.map_ind[sp1][i]) * N(self.map[sp2][j],self.map_ind[sp2][j]) - else: - Hamiltonian += 0.5 * self.Up[self.offset+i,self.offset+j] * N(self.map[sp1][i],self.map_ind[sp1][i]) * N(self.map[sp2][j],self.map_ind[sp2][j]) - - Hamiltonian -= N(self.map[spinblocs[0]][0],0) # substract the initializing value - - return Hamiltonian - - - def __set_full_hamiltonian_slater(self): - - spinblocs = [v for v in self.map] - Hamiltonian = N(self.map[spinblocs[0]][0],0) # initialize it - #print "Starting..." - # use the full 4-index U-matrix: - #for sp1 in spinblocs: - # for sp2 in spinblocs: - for m1 in range(self.n_orb): - for m2 in range(self.n_orb): - for m3 in range(self.n_orb): - for m4 in range(self.n_orb): - if (abs(self.U4ind[self.offset+m1,self.offset+m2,self.offset+m3,self.offset+m4])>0.00001): - for sp1 in spinblocs: - for sp2 in spinblocs: - #print sp1,sp2,m1,m2,m3,m4 - Hamiltonian += 0.5 * self.U4ind[self.offset+m1,self.offset+m2,self.offset+m3,self.offset+m4] * \ - Cdag(self.map[sp1][m1],self.map_ind[sp1][m1]) * Cdag(self.map[sp2][m2],self.map_ind[sp2][m2]) * C(self.map[sp2][m4],self.map_ind[sp2][m4]) * C(self.map[sp1][m3],self.map_ind[sp1][m3]) - #print "end..." - Hamiltonian -= N(self.map[spinblocs[0]][0],0) # substract the initializing value - - return Hamiltonian - - - def __set_spinflip_hamiltonian_slater(self): - """Takes only spin-flip and pair-hopping terms""" - - spinblocs = [v for v in self.map] - Hamiltonian = N(self.map[spinblocs[0]][0],0) # initialize it - #print "Starting..." - # use the full 4-index U-matrix: - #for sp1 in spinblocs: - # for sp2 in spinblocs: - for m1 in range(self.n_orb): - for m2 in range(self.n_orb): - for m3 in range(self.n_orb): - for m4 in range(self.n_orb): - if ((abs(self.U4ind[self.offset+m1,self.offset+m2,self.offset+m3,self.offset+m4])>0.00001) and - ( ((m1==m2)and(m3==m4)) or ((m1==m3)and(m2==m4)) or ((m1==m4)and(m2==m3)) ) ): - for sp1 in spinblocs: - for sp2 in spinblocs: - #print sp1,sp2,m1,m2,m3,m4 - Hamiltonian += 0.5 * self.U4ind[self.offset+m1,self.offset+m2,self.offset+m3,self.offset+m4] * \ - Cdag(self.map[sp1][m1],self.map_ind[sp1][m1]) * Cdag(self.map[sp2][m2],self.map_ind[sp2][m2]) * C(self.map[sp2][m4],self.map_ind[sp2][m4]) * C(self.map[sp1][m3],self.map_ind[sp1][m3]) - #print "end..." - Hamiltonian -= N(self.map[spinblocs[0]][0],0) # substract the initializing value - - return Hamiltonian - - - - def __set_full_hamiltonian_kanamori(self,J_hund): - - spinblocs = [v for v in self.map] - assert len(spinblocs)==2,"spinflips in Kanamori representation only implemented for up/down structure!" - - Hamiltonian = N(self.map[spinblocs[0]][0],0) # initialize it - - # density terms: - for sp1 in spinblocs: - for sp2 in spinblocs: - for i in range(self.n_orb): - for j in range(self.n_orb): - if (sp1==sp2): - Hamiltonian += 0.5 * self.U[self.offset+i,self.offset+j] * N(self.map[sp1][i],self.map_ind[sp1][i]) * N(self.map[sp2][j],self.map_ind[sp2][j]) - else: - Hamiltonian += 0.5 * self.Up[self.offset+i,self.offset+j] * N(self.map[sp1][i],self.map_ind[sp1][i]) * N(self.map[sp2][j],self.map_ind[sp2][j]) - - # spinflip term: - sp1 = spinblocs[0] - sp2 = spinblocs[1] - for i in range(self.n_orb-1): - for j in range(i+1,self.n_orb): - Hamiltonian -= J_hund * ( Cdag(self.map[sp1][i],self.map_ind[sp1][i]) * C(self.map[sp2][i],self.map_ind[sp2][i]) * Cdag(self.map[sp2][j],self.map_ind[sp2][j]) * C(self.map[sp1][j],self.map_ind[sp1][j]) ) # first term - Hamiltonian -= J_hund * ( Cdag(self.map[sp2][i],self.map_ind[sp2][i]) * C(self.map[sp1][i],self.map_ind[sp1][i]) * Cdag(self.map[sp1][j],self.map_ind[sp1][j]) * C(self.map[sp2][j],self.map_ind[sp2][j]) ) # second term - - # pairhop terms: - for i in range(self.n_orb-1): - for j in range(i+1,self.n_orb): - Hamiltonian -= J_hund * ( Cdag(self.map[sp1][i],self.map_ind[sp1][i]) * Cdag(self.map[sp2][i],self.map_ind[sp2][i]) * C(self.map[sp1][j],self.map_ind[sp1][j]) * C(self.map[sp2][j],self.map_ind[sp2][j]) ) # first term - Hamiltonian -= J_hund * ( Cdag(self.map[sp2][j],self.map_ind[sp2][j]) * Cdag(self.map[sp1][j],self.map_ind[sp1][j]) * C(self.map[sp2][i],self.map_ind[sp2][i]) * C(self.map[sp1][i],self.map_ind[sp1][i]) ) # second term - - Hamiltonian -= N(self.map[spinblocs[0]][0],0) # substract the initializing value - - return Hamiltonian - - - def __set_quantum_numbers(self,gf_struct): - - QN = {} - spinblocs = [v for v in self.map] - - # Define the quantum numbers: - if (self.use_spinflip) : - Ntot = sum_list( [ N(self.map[s][i],self.map_ind[s][i]) for s in spinblocs for i in range(self.n_orb) ] ) - QN['NtotQN'] = Ntot - #QN['Ntot'] = sum_list( [ N(self.map[s][i],i) for s in spinblocs for i in range(self.n_orb) ] ) - if (len(spinblocs)==2): - # Assuming up/down structure: - Sz = sum_list( [ N(self.map[spinblocs[0]][i],self.map_ind[spinblocs[0]][i])-N(self.map[spinblocs[1]][i],self.map_ind[spinblocs[1]][i]) for i in range(self.n_orb) ] ) - QN['SzQN'] = Sz - # new quantum number: works only if there are only spin-flip and pair hopping, not any more complicated things - for i in range(self.n_orb): - QN['Sz2_%s'%i] = (N(self.map[spinblocs[0]][i],self.map_ind[spinblocs[0]][i])-N(self.map[spinblocs[1]][i],self.map_ind[spinblocs[1]][i])) * (N(self.map[spinblocs[0]][i],self.map_ind[spinblocs[0]][i])-N(self.map[spinblocs[1]][i],self.map_ind[spinblocs[1]][i])) - - else : - for ibl in range(len(gf_struct)): - QN['N%s'%gf_struct[ibl][0]] = sum_list( [ N(gf_struct[ibl][0],gf_struct[ibl][1][i]) for i in range(len(gf_struct[ibl][1])) ] ) - - return QN - - - def fit_tails(self): - """Fits the tails using the constant value for the Re Sigma calculated from F=Sigma*G. - Works only for blocks of size one.""" - - #if (len(self.gf_struct)==2*self.n_orb): - if (self.blocssizeone): - spinblocs = [v for v in self.map] - mpi.report("Fitting tails manually") - - known_coeff = numpy.zeros([1,1,2],numpy.float_) - msh = [x.imag for x in self.G[self.map[spinblocs[0]][0]].mesh ] - fit_start = msh[self.fitting_Frequency_Start] - fit_stop = msh[self.N_Frequencies_Accumulated] - - # Fit the tail of G just to get the density - for n,g in self.G: - g.fitTail([[[0,0,1]]],7,fit_start,2*fit_stop) - densmat = self.G.density() - - for sig1 in spinblocs: - for i in range(self.n_orb): - - coeff = 0.0 - - for sig2 in spinblocs: - for j in range(self.n_orb): - if (sig1==sig2): - coeff += self.U[self.offset+i,self.offset+j] * densmat[self.map[sig1][j]][0,0].real - else: - coeff += self.Up[self.offset+i,self.offset+j] * densmat[self.map[sig2][j]][0,0].real - - known_coeff[0,0,1] = coeff - self.Sigma[self.map[sig1][i]].fitTail(fixed_coef = known_coeff, order_max = 3, fit_start = fit_start, fit_stop = fit_stop) - - else: - - for n,sig in self.Sigma: - - known_coeff = numpy.zeros([sig.N1,sig.N2,1],numpy.float_) - msh = [x.imag for x in sig.mesh] - fit_start = msh[self.fitting_Frequency_Start] - fit_stop = msh[self.N_Frequencies_Accumulated] - - sig.fitTail(fixed_coef = known_coeff, order_max = 3, fit_start = fit_start, fit_stop = fit_stop) - - - - -class SolverMultiBandOld(SolverMultiBand): - """ - Old MultiBand Solver construct - """ - - def __init__(self, Beta, Norb, U_interact=None, J_Hund=None, GFStruct=False, map=False, use_spinflip=False, - useMatrix = True, l=2, T=None, dimreps=None, irep=None, deg_orbs = [], Sl_Int = None): - - SolverMultiBand.__init__(self, beta=Beta, n_orb=Norb, gf_struct=GFStruct, map=map) - self.U_interact = U_interact - self.J_Hund = J_Hund - self.use_spinflip = use_spinflip - self.useMatrix = useMatrix - self.l = l - self.T = T - self.dimreps = dimreps - self.irep = irep - self.deg_orbs = deg_orbs - self.Sl_Int = Sl_Int - self.gen_keys = copy.deepcopy(self.__dict__) - - msg = """ -********************************************************************************** - Warning: You are using the old constructor for the solver. Beware that this will - be deprecated in future versions. Please check the documentation. -********************************************************************************** -""" - mpi.report(msg) - - - def Solve(self): - - params = copy.deepcopy(self.__dict__) - for i in self.gen_keys: self.params.pop(i) - self.params.pop("gen_keys") - self.solve(self, U_interact=self.U_interact, J_hund=self.J_Hund, use_spinflip=self.use_spinflip, - use_matrix = self.useMatrix, l=self.l, T=self.T, dim_reps=self.dimreps, irep=self.irep, - deg_orbs = self.deg_orbs, sl_int = self.Sl_Int, **params) - - - - - -def set_U_matrix(U_interact,J_hund,n_orb,l,use_matrix=True,T=None,sl_int=None,use_spinflip=False,dim_reps=None,irep=None): - """ Set up the interaction vertex""" - - offset = 0 - U4ind = None - U = None - Up = None - if (use_matrix): - if not (sl_int is None): - Umat = Umatrix(l=l) - assert len(sl_int)==(l+1),"sl_int has the wrong length" - if (type(sl_int)==ListType): - Rcl = numpy.array(sl_int) - else: - Rcl = sl_int - Umat(T=T,rcl=Rcl) - else: - if ((U_interact==None)and(J_hund==None)): - mpi.report("Give U,J or Slater integrals!!!") - assert 0 - Umat = Umatrix(U_interact=U_interact, J_hund=J_hund, l=l) - Umat(T=T) - - Umat.reduce_matrix() - if (Umat.N==Umat.Nmat): - # Transformation T is of size 2l+1 - U = Umat.U - Up = Umat.Up - else: - # Transformation is of size 2(2l+1) - U = Umat.U - # now we have the reduced matrices U and Up, we need it for tail fitting anyways - - if (use_spinflip): - #Take the 4index Umatrix - # check for imaginary matrix elements: - if (abs(Umat.Ufull.imag)>0.0001).any(): - mpi.report("WARNING: complex interaction matrix!! Ignoring imaginary part for the moment!") - mpi.report("If you want to change this, look into Wien2k/solver_multiband.py") - U4ind = Umat.Ufull.real - - # this will be changed for arbitrary irep: - # use only one subgroup of orbitals? - if not (irep is None): - #print irep, dim_reps - assert not (dim_reps is None), "Dimensions of the representatives are missing!" - assert n_orb==dim_reps[irep-1],"Dimensions of dimrep and n_orb do not fit!" - for ii in range(irep-1): - offset += dim_reps[ii] - else: - if ((U_interact==None)and(J_hund==None)): - mpi.report("For Kanamori representation, give U and J!!") - assert 0 - U = numpy.zeros([n_orb,n_orb],numpy.float_) - Up = numpy.zeros([n_orb,n_orb],numpy.float_) - for i in range(n_orb): - for j in range(n_orb): - if (i==j): - Up[i,i] = U_interact - else: - Up[i,j] = U_interact - 2.0*J_hund - U[i,j] = U_interact - 3.0*J_hund - - return U, Up, U4ind, offset + raise RuntimeError("Hamiltonian type not implemented.")