Source code for triqs.operators.util.observables

# Copyright (c) 2015-2017 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
# Copyright (c) 2015-2017 Centre national de la recherche scientifique (CNRS)
# Copyright (c) 2020-2022 Simons Foundation
# Copyright (c) 2015-2017 Igor Krivenko
#
# This program 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.
#
# This program 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 may obtain a copy of the License at
#     https:#www.gnu.org/licenses/gpl-3.0.txt
#
# Authors: Michel Ferrero, Igor Krivenko, Nils Wentzell

"""Operators for commonly used observables."""

import numpy as np
from triqs.operators.operators import Operator, n, c_dag, c
from triqs.operators.util.op_struct import get_mkind
from triqs.operators.util.U_matrix import spherical_to_cubic
from itertools import product

pauli_matrix = {'x' : np.array([[0,1],[1,0]]),
                'y' : np.array([[0,-1j],[1j,0]]),
                'z' : np.array([[1,0],[0,-1]]),
                '+' : np.array([[0,2],[0,0]]),
                '-' : np.array([[0,0],[2,0]])}

# Helper function for backward compat and improved error messages
[docs] def check_backward_compat(fname,n_orb,orb_names): if orb_names is not None: raise RuntimeError("Argument orb_names is no longer supported. Please provide n_orb instead.") if isinstance(n_orb, list): import warnings warnings.warn("Positional argument n_orb of function {} should provide the number of orbitals, not a list".format(fname)) n_orb = len(n_orb) return n_orb
[docs] def N_op(spin_names, n_orb, off_diag = None, map_operator_structure = None, orb_names = None): r""" Create an operator of the total number of particles. .. math:: \hat N = \sum_{i\sigma} a_{i\sigma}^\dagger a_{i\sigma}. Parameters ---------- spin_names : list of strings Names of the spins, e.g. ['up','down']. n_orb : int Number of orbitals. off_diag : boolean Do we have (orbital) off-diagonal elements? If yes, the operators and blocks are denoted by ('spin', 'orbital'), otherwise by ('spin_orbital',0). map_operator_structure : dict Mapping of names of GF blocks names from one convention to another, e.g. {('up', 0): ('up_0', 0), ('down', 0): ('down_0',0)}. If provided, the operators and blocks are denoted by the mapping of ``('spin', 'orbital')``. Returns ------- N : Operator The total number of particles. """ n_orb = check_backward_compat("N_op",n_orb,orb_names) mkind = get_mkind(off_diag,map_operator_structure) N = Operator() for sn, o in product(spin_names,range(n_orb)): N += n(*mkind(sn,o)) return N
[docs] def S_op(component, spin_names, n_orb, off_diag = None, map_operator_structure = None, orb_names = None): r""" Create a component of the spin vector operator. .. math:: \hat S_{x,y,z} = \frac{1}{2}\sum_{i\sigma\sigma'} a^\dagger_{i\sigma} \mathbf{\tau}^{x,y,z}_{\sigma\sigma'} a_{i\sigma'}, \quad\hat S_\pm = \hat S_x \pm i \hat S_y. Parameters ---------- component : string Component to be created, one of 'x', 'y', 'z', '+', or '-'. spin_names : list of strings Names of the spins, e.g. ['up','down']. n_orb : int Number of orbitals. off_diag : boolean Do we have (orbital) off-diagonal elements? If yes, the operators and blocks are denoted by ('spin', 'orbital'), otherwise by ('spin_orbital',0). map_operator_structure : dict Mapping of names of GF blocks names from one convention to another, e.g. {('up', 0): ('up_0', 0), ('down', 0): ('down_0',0)}. If provided, the operators and blocks are denoted by the mapping of ``('spin', 'orbital')``. Returns ------- S : Operator The component of the spin vector operator. """ n_orb = check_backward_compat("S_op",n_orb,orb_names) mkind = get_mkind(off_diag,map_operator_structure) pm = pauli_matrix[component] S = Operator() spin_range = list(range(len(spin_names))) for n1, n2 in product(spin_range,spin_range): for o in range(n_orb): S += 0.5 * c_dag(*mkind(spin_names[n1],o)) * pm[n1,n2] * c(*mkind(spin_names[n2],o)) return S
[docs] def S2_op(spin_names, n_orb, off_diag = None, map_operator_structure = None, orb_names = None): r""" Create the square of the total spin operator. .. math:: \hat S^2 = \hat S_x^2 + \hat S_y^2 + \hat S_z^2. Parameters ---------- spin_names : list of strings Names of the spins, e.g. ['up','down']. n_orb : int Number of orbitals. off_diag : boolean Do we have (orbital) off-diagonal elements? If yes, the operators and blocks are denoted by ('spin', 'orbital'), otherwise by ('spin_orbital',0). map_operator_structure : dict Mapping of names of GF blocks names from one convention to another, e.g. {('up', 0): ('up_0', 0), ('down', 0): ('down_0',0)}. If provided, the operators and blocks are denoted by the mapping of ``('spin', 'orbital')``. Returns ------- S2 : Operator The square of the total spin operator. """ n_orb = check_backward_compat("S2_op",n_orb,orb_names) Sz, Sp, Sm = [S_op(k,spin_names,n_orb,off_diag,map_operator_structure) for k in ('z','+','-')] return Sz*Sz + 0.5*(Sp*Sm + Sm*Sp)
[docs] def L_op(component, spin_names, n_orb, off_diag = None, map_operator_structure = None, basis='spherical', T=None, orb_names = None): r""" Create a component of the orbital momentum vector operator. .. math:: \hat L_{z,+,-} &= \sum_{ii'\sigma} a^\dagger_{i\sigma} L^{z,+,-}_{ii'} a_{i'\sigma},\\ \hat L_x &= \frac{1}{2}(\hat L_+ + \hat L_-),\ \hat L_y = \frac{1}{2i}(\hat L_+ - \hat L_-),\\ L^z_{ii'} &= i\delta_{i,i'}, \ L^+_{ii'} = \delta_{i,i'+1}\sqrt{l(l+1)-i'(i'+1)}, \ L^+_{ii'} = \delta_{i,i'-1}\sqrt{l(l+1)-i'(i'-1)}. Parameters ---------- component : string Component to be created, one of 'x', 'y', 'z', '+', or '-'. spin_names : list of strings Names of the spins, e.g. ['up','down']. n_orb : int Number of orbitals. off_diag : boolean Do we have (orbital) off-diagonal elements? If yes, the operators and blocks are denoted by ('spin', 'orbital'), otherwise by ('spin_orbital',0). map_operator_structure : dict Mapping of names of GF blocks names from one convention to another, e.g. {('up', 0): ('up_0', 0), ('down', 0): ('down_0',0)}. If provided, the operators and blocks are denoted by the mapping of ``('spin', 'orbital')``. basis : string, optional The basis in which the interaction matrix should be computed. Takes the values - 'spherical': spherical harmonics, - 'cubic': cubic harmonics (valid only for the integer orbital momenta, i.e. for odd values of n_orb), - 'other': other basis type as given by the transformation matrix T. T : real/complex numpy array, optional Transformation matrix for basis change. Must be provided if basis='other'. Returns ------- L : Operator The component of the orbital momentum vector operator. """ n_orb = check_backward_compat("L_op",n_orb,orb_names) l = (n_orb-1)/2.0 L_melem_dict = {'z' : lambda m,mp: m if np.isclose(m,mp) else 0, '+' : lambda m,mp: np.sqrt(l*(l+1)-mp*(mp+1)) if np.isclose(m,mp+1) else 0, '-' : lambda m,mp: np.sqrt(l*(l+1)-mp*(mp-1)) if np.isclose(m,mp-1) else 0, 'x' : lambda m,mp: 0.5*(L_melem_dict['+'](m,mp) + L_melem_dict['-'](m,mp)), 'y' : lambda m,mp: -0.5j*(L_melem_dict['+'](m,mp) - L_melem_dict['-'](m,mp))} L_melem = L_melem_dict[component] orb_range = list(range(int(2*l+1))) L_matrix = np.array([[L_melem(o1-l,o2-l) for o2 in orb_range] for o1 in orb_range]) # Transform from spherical basis if needed if basis == "cubic": if not np.isclose(np.mod(l,1),0): raise ValueError("L_op: cubic basis is only defined for the integer orbital momenta.") T = spherical_to_cubic(int(l)) if basis == "other" and T is None: raise ValueError("L_op: provide T for other bases.") if T is not None: L_matrix = np.einsum("ij,jk,kl",np.conj(T),L_matrix,np.transpose(T)) mkind = get_mkind(off_diag,map_operator_structure) L = Operator() for sn in spin_names: for o1, o2 in product(orb_range,orb_range): L += c_dag(*mkind(sn,o1)) * L_matrix[o1,o2] * c(*mkind(sn,o2)) return L
[docs] def L2_op(spin_names, n_orb, off_diag = None, map_operator_structure = None, basis='spherical', T=None, orb_names = None): r""" Create the square of the orbital momentum operator. .. math:: \hat L^2 = \hat L_x^2 + \hat L_y^2 + \hat L_z^2. Parameters ---------- spin_names : list of strings Names of the spins, e.g. ['up','down']. n_orb : int Number of orbitals. off_diag : boolean Do we have (orbital) off-diagonal elements? If yes, the operators and blocks are denoted by ('spin', 'orbital'), otherwise by ('spin_orbital',0). map_operator_structure : dict Mapping of names of GF blocks names from one convention to another, e.g. {('up', 0): ('up_0', 0), ('down', 0): ('down_0',0)}. If provided, the operators and blocks are denoted by the mapping of ``('spin', 'orbital')``. basis : string, optional The basis in which the interaction matrix should be computed. Takes the values - 'spherical': spherical harmonics, - 'cubic': cubic harmonics (valid only for the integer orbital momenta, i.e. for odd values of n_orb), - 'other': other basis type as given by the transformation matrix T. T : real/complex numpy array, optional Transformation matrix for basis change. Must be provided if basis='other'. Returns ------- L2 : Operator The square of the orbital momentum operator. """ n_orb = check_backward_compat("L2_op",n_orb,orb_names) Lz, Lp, Lm = [L_op(k,spin_names,n_orb,off_diag, map_operator_structure, basis, T) for k in ('z','+','-')] return Lz*Lz + 0.5*(Lp*Lm + Lm*Lp)
[docs] def LS_op(spin_names, n_orb, off_diag = None, map_operator_structure = None, basis='spherical', T=None, orb_names = None): r""" Create a spin-orbital coupling operator. .. math:: \mathbf{\hat L\cdot\hat S} = \hat L_x \hat S_x + \hat L_y \hat S_y + \hat L_z \hat S_z. Parameters ---------- spin_names : list of strings Names of the spins, e.g. ['up','down']. n_orb : int Number of orbitals. off_diag : boolean Do we have (orbital) off-diagonal elements? If yes, the operators and blocks are denoted by ('spin', 'orbital'), otherwise by ('spin_orbital',0). map_operator_structure : dict Mapping of names of GF blocks names from one convention to another, e.g. {('up', 0): ('up_0', 0), ('down', 0): ('down_0',0)}. If provided, the operators and blocks are denoted by the mapping of ``('spin', 'orbital')``. basis : string, optional The basis in which the interaction matrix should be computed. Takes the values - 'spherical': spherical harmonics, - 'cubic': cubic harmonics (valid only for the integer orbital momenta, i.e. for odd values of n_orb), - 'other': other basis type as given by the transformation matrix T. T : real/complex numpy array, optional Transformation matrix for basis change. Must be provided if basis='other'. Returns ------- LS : Operator The spin-orbital coupling operator. """ n_orb = check_backward_compat("LS_op",n_orb,orb_names) Sz, Sp, Sm = [S_op(k,spin_names,n_orb,off_diag,map_operator_structure) for k in ('z','+','-')] Lz, Lp, Lm = [L_op(k,spin_names,n_orb,off_diag,map_operator_structure, basis, T) for k in ('z','+','-')] return Lz*Sz + 0.5*(Lp*Sm + Lm*Sp)