mirror of
https://github.com/triqs/dft_tools
synced 2024-12-21 11:53:41 +01:00
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
This commit is contained in:
Priyanka Seth
parent
6fadfecf33
commit
f803c13285
@ -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 <http://www.gnu.org/licenses/>.
|
||||
#
|
||||
################################################################################
|
||||
|
||||
# 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
|
||||
|
||||
@ -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']
|
||||
|
||||
127
python/hamiltonians.py
Normal file
127
python/hamiltonians.py
Normal file
@ -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
|
||||
@ -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 <http://www.gnu.org/licenses/>.
|
||||
#
|
||||
################################################################################
|
||||
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.")
|
||||
|
||||
Loading…
Reference in New Issue
Block a user