add many body operators in python (old version)

- The old code, which was in the solvers,
is used in several apps (dft, ctqmc_hyb).
It is really a lib, even though it is now superseded
by Igor and Andrey's library...
(to be deprecated at some point).

pytriqs/operators/CMakeLists.txt

@@ -0,0 +1,9 @@
+
+SET(PYTHON_SOURCES 
+  ${CMAKE_CURRENT_SOURCE_DIR}/__init__.py
+  ${CMAKE_CURRENT_SOURCE_DIR}/operators.py
+)
+
+install (FILES ${PYTHON_SOURCES} DESTINATION ${TRIQS_PYTHON_LIB_DEST}/operators)
+
+

pytriqs/operators/__init__.py

@@ -0,0 +1,24 @@
+
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by M. Ferrero, O. Parcollet
+#
+# 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/>.
+#
+################################################################################
+
+from operators import Operator, C, Cdag, N, state, all_pure_states, C_list, number_of_C, commutator, anti_commutator, operator_id, extend_function_on_fundamentals
+

pytriqs/operators/operators.py

@@ -0,0 +1,438 @@
+
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by M. Ferrero, O. Parcollet
+#
+# 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/>.
+#
+################################################################################
+
+import string
+from copy import deepcopy
+"""
+ State : represented by a dict : n-> lambda_n
+ n = binary decomposition is the occupation number as usual
+"""
+        
+class OpC:
+    """ C or B (?) operator """
+    def __init__(self,dagger,n,*args):
+        self.name,self.dag = args,dagger
+        self.number,self.two_num = n,2**n
+        self.char = {}
+        self.isFermion = True
+
+    WhichOrder=0
+
+    def num(self) :
+        if self.dag and OpC.WhichOrder != 0 :
+            return -self.number
+        else :
+            return self.number
+
+    def dagger(self) :
+        if self.dag : return apply(C,self.name).L[0][2][0]
+        else : return apply(Cdag,self.name).L[0][2][0]
+        
+    def __repr__(self):
+        s = reduce (lambda x,y : str(x) + ',%s'%y,self.name)
+        if self.dag:
+            return('C^(%s)'%s)
+        else :
+            return('C(%s)'%s)
+        
+    def __call__(self,state):
+        """ State : cf header of the file"""
+        res,fil = {}, self.two_num
+        def signe(p):
+            r=0
+            for i in xrange(self.number):
+                r += p & 1
+                p >>= 1
+            return 1 - 2* (r%2)
+        if self.dag :
+            for n,l in state.items() :
+                if n & fil ==0 :
+                    res[n+fil] = signe(n)*l
+        else:
+            for n,l in state.items() :
+                if n & fil !=0 :
+                    res[n-fil] = signe(n)*l
+        return res
+
+################################################################
+
+class Operator:
+    """ # L is a list of [ float, [liste of name of var],[list of OPbase to mult]] """
+    def __init__ (self,*args):
+        """Construct from xxxx"""
+        if args:
+            self.L=[ [ 1 , [] ,[ apply(OpC,args) ] ]  ]
+        else:
+            self.L=[ ]
+        #self.const = 0
+
+    def copy(self) :
+        res = Operator()
+        res.L = [ [c,vl[:],mono[:]] for c,vl,mono in self.L]
+        return res
+    
+    def remove_quadratic(self):
+        res = Operator()
+        res.L = [ [c,vl[:],mono[:]] for c,vl,mono in self.L if len(mono) >2 ]
+        return res
+
+
+    def normal_order(self):
+        """Returns an new operator put in normal order : Cdag first, C at the right.
+           Cdag are in increasing order, C in decreasing order."""
+        OpC.WhichOrder = 1
+        #res = self
+        self.__redu()
+        OpC.WhichOrder = 0
+        return self 
+    
+    def __redu(self) :
+        """ Internal use"""
+        def tri1(IN):
+            oplist = IN[2]
+            for i in range(1,len(oplist)):
+                n,c= oplist[i].num(),1
+                d = n - oplist[i-1].num()
+                if (d<0) :
+                    j,c = i-1,1
+                    while j>=0 and n<oplist[j].num() : c *= -1; j -= 1
+                    oplist.insert(j+1,oplist.pop(i))
+                    return  [ [c*IN[0], IN[1][:], oplist] ]
+                if (d==0) :
+                    if oplist[i].dag ==  oplist[i-1].dag : return []
+                    if oplist[i].dag :
+                        oplist.insert(i-1,oplist.pop(i))
+                        res = [ [ -IN[0], IN[1][:], oplist] ]
+                        l2 =oplist[0:i-1] + oplist[i+1:]
+                        res +=  [ [IN[0],IN[1][:], l2 ] ]
+                        return res
+        Fini =[]
+        while self.L :
+            Work = []
+            for l in self.L :
+                #print "TRI : input",l
+                r = tri1(l)
+                #print "TRI RETURNS",r
+                if r==None : Fini.append(l)
+                else :  Work+=r 
+            self.L = Work
+            #print "------------------------ LOOP",Work
+        self.L = Fini
+                   
+	for n,l  in enumerate(self.L):
+            for b2 in self.L[n+1:] :
+		if b2[1]==l[1] and b2[2]==l[2] :
+		    l[0] +=b2[0]
+                    b2[0]=0
+        self.L = filter( lambda x : (abs(x[0])>=1.e-14),self.L)
+
+    def dagger(self) :
+        res = Operator()
+        for coef,vl,mono in self.L :
+            res.L.append([coef.conjugate(),vl[:],[ o.dagger() for o in reversed(mono)]])
+        res.__redu()
+        return res
+
+    def is_zero(self):
+        return self.L==[] or max([abs(x[0]) for x in self.L])<1.e-13
+
+    def __add__(self,b):
+        assert type(b)==type(self)
+        res = Operator()
+        # we do not want to copy by reference
+        for l in self.L:
+            res.L.append([ l[0],list(l[1]),list(l[2]) ])
+        for l in b.L:
+            res.L.append([ l[0],list(l[1]),list(l[2]) ])
+        res.__redu()
+        return(res)
+
+    def __sub__(self,b):
+        assert type(b)==type(self)
+        res = Operator()
+        for l in self.L:
+            res.L.append([ l[0],list(l[1]),list(l[2]) ])
+        for l in b.L:
+            res.L.append([ -l[0],list(l[1]),list(l[2]) ])
+        res.__redu()
+        return(res)
+    
+    def __mul__(self,b):
+        assert type(b)==type(self)
+        res = Operator()
+        for m1 in self.L :
+            for m2 in b.L :
+                res.L.append( [m1[0]*m2[0],sorted(m1[1]+m2[1]), list(m1[2]+m2[2])] )
+        res.__redu()
+        return(res)
+
+    def __idiv__(self,b):
+        assert type(b) in [type(1), type(1.0), type(1j)]
+        for l in self.L :
+            l[0] /=b
+        #self.__redu()
+        return self
+    
+    def __div__(self,b):
+        res = self.copy()
+        res /= b
+        return(res)
+
+    def __rmul__(self,b):
+        if type(b)==type(self) :
+            return(b*self)
+        elif type(b) in [type(1),type(1j),type(1.0) ]:
+            res = Operator()
+            res.L = [ [ b*l[0],list(l[1]),list(l[2]) ] for l in self.L]
+            res.__redu()
+            return(res)
+    
+    def __repr__(self):
+        s=''
+        for bl in self.L:
+            if type(bl[0]) is not complex:
+                x,signe,s1= bl[0],1,''
+                if (x<0) : signe,x = -1,abs(x)
+                if (abs(x)<1.e-10) : signe=0
+                if (abs(x-1)>1.e-10) : s1=repr(x)+' '
+                for var in bl[1]:
+                    s1+=var+' '
+                for op in bl[2]:
+                    s1+=repr(op)
+                if s1=='' : s1='1'
+                if (signe==1) : s+='+ ' + s1 + ' '
+                if (signe==-1) : s+='- ' + s1 + ' '
+            else:
+                s += '+ (' + repr(bl[0]) + ') '
+                for op in bl[2]:
+                    s+=repr(op)
+                s += ' '
+
+        if (s=='') : return('0')
+        if (s[0]=="+") : return(s[1:])
+        return(s)
+
+    def __call__(self,state) :
+        res = {}
+        for coef,var,oplist in self.L :
+            r = state
+            # compute the action of a monomial
+            for op in reversed(oplist) :
+                r = op(r)
+            # add r to res with coef
+            for n,l in r.items() :
+                if n in res : res[n] += coef*l
+                else : res[n] = coef*l
+        return res
+
+    def make_coef_real_and_check(self) :
+        """ Transforms the coefficients """
+        def f(x) : 
+            if type(x) != type(1j) : return x
+	    assert abs(x.imag)<1.e-12
+            return x.real
+        self.L = [ (f(coef),var,oplist) for coef,var,oplist in self.L]
+
+    def number_as_C(self) :
+        """
+        If the operator is alpha C , return the number of C
+        """
+        if len(self.L)!=1 or len(self.L[0][2])!=1 :
+            raise "Internal error : len of C is not right"
+        return  self.L[0][2][0].number
+ 
+    def is_fundamental(self) :
+        """
+        If the operator is alpha C , return true
+        """
+        return len(self.L)==1 and len(self.L[0][2])==1
+
+    def as_fundamental(self) :
+        """
+        If the operator is alpha C , return the C
+        """
+        if not self.is_fundamental():
+            raise "Internal error : len of C is not right"
+        return  self.L[0][2][0]
+
+    def symmetry(self) :
+        """
+        If the operator is alpha C , return the symmetry of C
+        """
+        if len(self.L)!=1 or len(self.L[0][2])!=1 :
+            raise "Internal error : len of C is not right"
+        return  self.L[0][2][0].char
+
+ 
+_Varlist={}
+def var(x):
+    global _Varlist
+    _Varlist[x]=1
+    return(Op(1,[x],()))
+
+
+# list of all the C operators
+_C_list ={}
+
+def C_list():
+    for cdag,c,n,b in _C_list.values() :
+        yield n,cdag,c
+
+def C_list_names():
+    return _C_list.keys()
+
+def number_of_C():
+    return len(_C_list)
+
+def operator_id():
+    res = Operator()
+    res.L.append([1,[],[] ])
+    return res
+
+def C(*args):
+    """
+    C operator.
+    Usage : C(any argument) defines or retrieves a the operator::
+
+      print C('down',1) # gives the C_{down,1} operator
+
+    """
+    global _C_list
+    t = tuple(args)
+    if t not in _C_list :
+        _C_list[t] = apply(Operator,[True,len(_C_list)] + list(args)),apply(Operator,[False,len(_C_list)]+ list(args)),len(_C_list),2**len(_C_list)
+    return _C_list [t][1]
+        
+def Cdag(*args):
+    """ C^dagger operator. Usage : Cdag(any argument) defines or retrieves a the operator"""
+    global _C_list
+    t = tuple(args)
+    if t not in _C_list :
+        _C_list[t] = apply(Operator,[True,len(_C_list)]+ list(args)),apply(Operator,[False,len(_C_list)]+ list(args)),len(_C_list),2**len(_C_list)
+    return _C_list [t][0]
+
+def N(*args):
+    """ Short cut for Cdag(args)*C(args)"""
+    return(apply(Cdag,args)*apply(C,args))
+
+def state(*args):
+    """
+Input : a list of indices for C operators. e.g. State( (1,up), (1,down) )
+Returns : the corresponding indices, whose binary decomposition is given by 000100..1000
+where the 1 are at the position described in input"""
+    r = 0
+    try : 
+        for x in args :
+            r +=  _C_list[tuple(x)][3]
+    except :
+        raise "One C operator unknown !"
+    return r
+    
+def all_pure_states():
+    """ Generates all the pure states"""
+    for i in xrange(2**len(_C_list)):
+        yield {i : 1}
+        
+def commutator(A,B) :
+    C = A*B - B*A
+    return C
+
+def anti_commutator(A,B) :
+    C = A*B + B*A
+    return C
+
+def extend_function_on_fundamentals(func) :
+    """
+    Given a function func that map any C operator to another one,
+    it extends the function on the algebra of all operators.
+    Returns : a function acting an operators.
+    """
+    def f(Op) : 
+        RES = Operator()
+        for l,X,M in Op.L : # term is of type (lambda, [...], [MOMEOME])
+            r = [ func(c).as_fundamental() for c in M]
+            RES.L.append ( [l, X, r] )
+        RES._Operator__redu()
+        return RES
+    return f
+
+
+def complete_op_list_with_fundamentals (op_dict) : 
+    """
+      Given a dict of operators, it  : 
+       - finds all the fundamental operators contained in these operators
+       - add them to the dictionnary with a unique name
+       - renumbers them, without holes, with n for C, -n for the corresponding Cdagger
+       - prepares a dict : name_of_the_symmetry -> X
+          where X is a dict : FundamentalOperator.number_for_C -> Character (COMPLEX)
+       - returns (OpDict_Complete, Nfermions,Nbosons, SymCharacters,NameOpFundamentalList)
+
+    """      
+    # first find all fundamental operators
+    Opfund,NameOpFundamentalList = [],[] 
+    for name, Op in op_dict.items() :
+        for l,X,M in Op.L : # term is of type (lambda, [...], [MOMEOME])
+            Opfund += M
+    L = sorted( [ (Op.number, Op) for Op in set(Opfund) if not(Op.dag) ])
+    NF = sum ( [ 1 if c.isFermion else 0 for n,c in L]) # number of fermions
+    NB = sum ( [ 1 if not(c.isFermion) else 0 for n,c in L]) # number of bosons
+    SymCharacters={}
+    for (n,(i,OP)) in enumerate(L):
+        OP.number_for_C = n+1
+        OP.dagger().number_for_C =  - (n+1)
+        for sym_name,char in OP.dagger().char.items() : 
+            if sym_name in SymCharacters : 
+                SymCharacters[sym_name][n] = char
+            else :
+                SymCharacters[sym_name] = {n : char}
+        OP_C = C(*OP.name)
+        if len([ k for (k,v) in op_dict.items() if (v-OP_C).is_zero()]) == 0: 
+            k = 'C'    + string.join([ "%s"%x for x in OP.name],'')
+            op_dict[k] = OP_C
+            NameOpFundamentalList.append(k)
+        OP_C = Cdag(*OP.name)
+        if len([ k for (k,v) in op_dict.items() if (v-OP_C).is_zero()]) == 0: 
+            k = 'Cdag' + string.join([ "%s"%x for x in OP.name],'')
+            op_dict[k] = OP_C
+            NameOpFundamentalList.append(k)
+
+    return op_dict,NF,NB,SymCharacters,NameOpFundamentalList
+
+
+def transcribe_op_list_for_C (op_dict) : 
+    """
+      Given a dict of operators, it  : 
+       - transcribes all operators into a dict D : 
+           name - >  [  [ COMPLEX, [ int] ]     ]
+       - returns D
+
+    """      
+    D = {}
+    for name, Op in op_dict.items() :
+        opval = []
+        for l,X,M in Op.L : # term is of type (lambda, [...], [MOMEOME])
+            RES = [ c.number_for_C for c in M]
+            RES.reverse()
+            opval.append ( [l, RES] )
+        D[name] = opval
+    return D
+