18 files changed, 339 insertions, 150 deletions
diff --git a/Makefile b/Makefile
index 79fee56..f6233fe 100644
--- a/Makefile
+++ b/Makefile
@@ -18,7 +18,7 @@
 #   if static=1 then link libkondo statically but other libraries dynamically
 STATIC=1
 
-VERSION=1.2.1
+VERSION=1.3
 
 # products of the compilation
 PROJECT_BINS= meankondo numkondo meantools kondo_preprocess meantools-convert
diff --git a/doc/meankondo-doc.html b/doc/meankondo-doc.html
index 27ffbbb..c4a7e2a 100644
--- a/doc/meankondo-doc.html
+++ b/doc/meankondo-doc.html
@@ -69,10 +69,10 @@
   </head>
 
   <body>
-    <h1 style="margin-bottom:50pt;">meankondo <span style="margin-left:10pt;font-size:18pt">v1.2</span></h1>
+    <h1 style="margin-bottom:50pt;">meankondo <span style="margin-left:10pt;font-size:18pt">v1.3</span></h1>
 
     <p>
-      This is the official documentation for <b>meankondo</b>, version 1.2. The aim of this document is not to give a technical description of how to use the various programs bundled with <b>meankondo</b>, nor is it to explain where hierarchical models come from and what their meaning is, but rather a conceptual overview of how <b>meankondo</b> approaches the computation of flow equations, and how its programs can be made to interact with one another to compute various quantities. For a more technical description, see the man pages included with the <b>meankondo</b> source code. For a more theoretical discussion of Fermionic hierarchical models, see <a href="http://ian.jauslin.org/publications/15bgj">[G.Benfatto, G.Gallavotti, I.Jauslin, 2015]</a>.
+      This is the official documentation for <b>meankondo</b>, version 1.3. The aim of this document is not to give a technical description of how to use the various programs bundled with <b>meankondo</b>, nor is it to explain where hierarchical models come from and what their meaning is, but rather a conceptual overview of how <b>meankondo</b> approaches the computation of flow equations, and how its programs can be made to interact with one another to compute various quantities. For a more technical description, see the man pages included with the <b>meankondo</b> source code. For a more theoretical discussion of Fermionic hierarchical models, see <a href="http://ian.jauslin.org/publications/15bgj">[G.Benfatto, G.Gallavotti, I.Jauslin, 2015]</a>.
     </p>
 
     <h2 style="margin-top:50pt;">Table of contents</h2>
@@ -150,7 +150,7 @@
 	<li><b>external</b>: which are organized in pairs, and are denoted by \((\Psi_i^+,\Psi_i^-)\) for \(i\in\{1,\cdots,E\}\).
 	<li><b>super-external</b>: which denoted by \(H_i\) for \(i\in\{1,\cdots,X\}\) (the only difference with external fields is that super-external fields are not in pairs, which is a seemingly innocuous difference; but super-external fields are meant to be used for different purposes as external fields (see <a href="#flow_equation_definition">Definition</a> below)).
       </ul>
-      The fields are used as a basis for a complex algebra, so that we can take products and linear combinations of fields (in other words, the concept of <i>polynomials over the fields</i> is well defined). Some of the fields (<i>Fermions</i>) anti-commute with each other (two fields \(a\) and \(b\) are said to anti-commute if \(ab\equiv-ba\)), and the rest (<i>Bosons</i>) commute. Which fields are Fermions and which are Bosons is specified in the <code>#!fields</code> entry in the configuration file. <b>(Warning: As of version 1.2, all internal fields must be Fermions.)</b>
+      The fields are used as a basis for a complex algebra, so that we can take products and linear combinations of fields (in other words, the concept of <i>polynomials over the fields</i> is well defined). Some of the fields (<i>Fermions</i>) anti-commute with each other (two fields \(a\) and \(b\) are said to anti-commute if \(ab\equiv-ba\)), and the rest (<i>Bosons</i>) commute. Which fields are Fermions and which are Bosons is specified in the <code>#!fields</code> entry in the configuration file. <b>(Warning: As of version 1.3, all internal fields must be Fermions.)</b>
     </p>
     <p>
       In the configuration file of the <b>meankondo</b> program, the fields are specified in the <code>#!fields</code> entry.
diff --git a/man/kondo_preprocess.1 b/man/kondo_preprocess.1
index 5ec1399..e6e5c59 100644
--- a/man/kondo_preprocess.1
+++ b/man/kondo_preprocess.1
@@ -1,5 +1,5 @@
 .Dd $Mdocdate: April 14 2015 $
-.Dt kondo_preprocess 1.2.1
+.Dt kondo_preprocess 1.3
 .Os
 .Sh NAME
 .Nm kondo_preprocess
@@ -94,6 +94,30 @@ defines external fields for A and B, denoted by a and b. They can be used as fie
 .D1 <axb.h>
 .Pp
 .It
+Scalar products of A's and B's may also be specified using the '<#.#>' syntax:
+.D1 <An.An>
+.D1 <Bn.Bn>
+.D1 <An.Bn>
+.D1 <An.h>
+.D1 <Bn.h>
+.Pp
+The difference between '[f #.#]' and '<#.#>' is that the former corresponds to a '#!symbols' entry whereas the latter is replaced by its corresponding polynomial when
+.Nm
+reads it (see
+.Sx meankondo Ns (1)).
+.Pp
+.It
+A vector 't=(t1,t2,t3)' of Pauli matrices (satisfying the Pauli commutation relations [ti,tj]=\\delta_{i,j}1+\\epsilon_{i,j,k}tk) is introduced as a non-commuting object. It can be used in scalar producs:
+.D1 <An.t>
+.D1 <Bn.t>
+.D1 <t.h>
+.D1 <a.t>
+.D1 <b.t>
+.Pp
+Note that the '<#,#>' must be used since these scalar products do not commute whereas '#!symbols' entries must commute (see
+.Sx meankondo Ns (1)).
+.Pp
+.It
 Furthermore, in order to simplify writing products of polynomials over each box index, if the polynomial contains a '%', then
 .Nm
 multiplies the polynomial by itself as many times as there are boxes (2^dimension times), replacing '%' with the appropriate box index. For example, if dimension=1
@@ -118,6 +142,9 @@ in which the polynomial can use the fields
 .D1 <a.h>
 .D1 <b.h>
 .D1 <axb.h>
+.D1 <t.h>
+.D1 <a.t>
+.D1 <b.t>
 defined above.
 .Pp
 Example:
diff --git a/man/meankondo.1 b/man/meankondo.1
index f7a0290..0a299d7 100644
--- a/man/meankondo.1
+++ b/man/meankondo.1
@@ -1,5 +1,5 @@
 .Dd $Mdocdate: April 13 2015 $
-.Dt meankondo 1.2.1
+.Dt meankondo 1.3
 .Os
 .Sh NAME
 .Nm meankondo
@@ -67,7 +67,7 @@ recognizes the following entries (unless explicitly mentioned, the entries below
 .It Sy #!fields
 A list of the fields of the model.
 .Pp
-The fields entry contains 4 lines which start with 'i:', 'x:', 'h:' and 'f:'. Each of these is followed by a ',' separated list of field indices, which are positive integers.
+The fields entry contains 5 lines which start with 'i:', 'x:', 'h:', 'f:' and 'a:'. Each of these is followed by a ',' separated list of field indices, which are positive integers.
 .Bl -bullet
 .It
 The indices following 'i' correspond to internal fields, which are integrated out using the Wick rule and the propagator provided in the '#!propagator' entry. Each internal field is associated a conjugate field, whose index is the opposite of the field's index (e.g. 'i:101' defines a field whose index is -101)
@@ -77,6 +77,10 @@ The indices following 'x' correspond to external fields that are associated conj
 The indices following 'h' correspond to external fields that are not associated a conjugate field. External indices may not appear as internal indices.
 .It
 The 'f' line specifies which of the internal and external indices are Fermions, i.e. which fields anti-commute. The fields appearing in the 'f' line should also either appear in the 'i' or 'x' line. WARNING: for the moment, only cases in which all of the internal fields are Fermions are supported.
+.It
+The 'a' line specifies a list of external fields listed in the 'h' entry that do not commute with each other. Specifying fields in this entry will prevent
+.Nm
+from sorting them. These fields may not be in the 'i', 'x' or 'f' entries. This entry can be used to treat cases in which the coefficients of the input polynomial are operators that do not commute. Their commutation relations may be specified in the '#!identities' entries (see below).
 .El
 .Pp
 .Em Line breaks are not ignored in this entry.
@@ -84,8 +88,9 @@ The 'f' line specifies which of the internal and external indices are Fermions,
 Example:
 .D1 i:101,102,201,202
 .D1 x:100,200
-.D1 h:301,302,303
+.D1 h:301,302,303,401,402,403
 .D1 f:100,101,102
+.D1 a:401,402,403
 .It Sy #!propagator
 The propagator of the model.
 .Pp
@@ -98,31 +103,41 @@ just as easily as propagators with symbolic entries. Such an entry means that
 .Pp
 Example:
 .D1 101;102: 1 , 102;101: -1 , 201;202: s{-1} + (-1)[l10] , 202;201: (-1)s{-1} + [l10]
-.It Sy #!symbols
-Symbolic variables used as shortcuts for more complicated expressions (optional entry).
 .Pp
-In order to simplify long expressions, symbolic variables can be defined in this entry. Each variable is assigned an index, which is a positive integer that must be different from any of the internel and external indices defined in the '#!fields' entry.
+.It Sy #!identities
+Identities satisfied by some of the fields (optional entry).
 .Pp
-The symbols entry is a ',' separated list, whose elements are of the form
-.D1 index= polynomial
-where index is the index of the variable and polynomial is the expression it stands for (see the POLYNOMIALS section below for information on how to format polynomials). Note that polynomial can contain other symbolic variables. There is no safeguard against self-referencing definitions that may cause infinite loops.
+In some cases, some of the quantities involved in a model will satisfy an identity (e.g. a vector may be of unit-norm, or non-commuting objects may satisfy non-trivial commutation relations), which should be simplified out from the flow equation.
+.Pp
+The identities entry is a ',' separated list, whose elements are of the form
+.D1 monomial=polynomial
+where monomial represents the left side of the identity and is a sequence of field indices of the form '[f index1][f index2]...' and polynomial represents the right side of the identity (see the POLYNOMIALS section below for information on how to format polynomials).
 .Pp
 Example:
-.D1 1001= (-1)[f-100][f100] + (-1)[f-101][f101] , 2001=[f-100][f100] + [f-201][f201]
+.D1 [f301][f301]=(1)+(-1)[f302][f302]+(-1)[f303][f303],
+.D1 [f401][f401]=(1),
+.D1 [f401][f402]=(s{-1})[f403],
+.D1 [f401][f403]=((-1)s{-1})[f402]
 .Pp
 This entry is optional.
 .Pp
-.It Sy #!identities
-Identities satisfied by some of the fields (optional entry).
+.It Sy #!symbols
+Symbolic variables used as shortcuts for more complicated expressions (optional entry).
 .Pp
-In some cases, some of the quantities involved in a model will satisfy an identity (e.g. a vector may be of unit-norm), which should simplified out from the flow equation.
+In order to simplify long expressions, symbolic variables can be defined in this entry. Each variable is assigned an index, which is a positive integer that must be different from any of the internal and external indices defined in the '#!fields' entry.
 .Pp
-The identities entry is a ',' separated list, whose elements are of the form
-.D1 monomial=polynomial
-where monomial represents the left side of the identity and is a sequence of field indices of the form '[f index1][f index2]...' and polynomial represents the right side of the identity (see the POLYNOMIALS section below for information on how to format polynomials).
+Seemingly similar functionality can be achieved using an '#!identity' entry (see above), though symbols are handled differently from identities. Indeed, while identities are simplified out of the polynomials as soon as they occur, symbols are only resolved when
+.Nm
+computes the mean of the input polynomial. Using symbols can thereby be a lot faster than using identities. However, as is mentioned below, symbols must commute with each other and all other fields, whereas identities can be made to be fermionic or non-commuting.
+.Pp
+The symbols entry is a ',' separated list, whose elements are of the form
+.D1 index= polynomial
+where index is the index of the variable and polynomial is the expression it stands for (see the POLYNOMIALS section below for information on how to format polynomials). Note that polynomial can contain other symbolic variables. There is no safeguard against self-referencing definitions that may cause infinite loops.
+.Pp
+WARNING: Symbols are assumed to commute with each other and all other Fermions. They should therefore not represent quantities that do not commute (e.g. odd monomials of fermions or non-commuting objects specified in the 'a:' entry in the '#!fields' entry).
 .Pp
 Example:
-.D1 [f301][f301]=(1)+(-1)[f302][f302]+(-1)[f303][f303]
+.D1 1001= (-1)[f-100][f100] + (-1)[f-101][f101] , 2001=[f-100][f100] + [f-201][f201]
 .Pp
 This entry is optional.
 .Pp
@@ -157,16 +172,16 @@ computes the mean of a monomial containing elements of different groups, it fact
 .Nm
 does not repeatedly try to pair independent fields.
 .Pp
+WARNING:
+.Nm
+assumes that the symbols and fields in each group are independent but does not check that they are. If symbols or fields that are not independent are put in different groups, or if some are in a group while others are not in any group, then the resulting flow equation may be wrong.
+.Pp
 The groups entry is a list of collections of fields or symbols of the following form
 .D1 (index1,index2,...)
 .Pp
 Example:
 .D1 (1001,1002) (2001,2002)
 .Pp
-.Em Warning:
-.Nm
-does not check that the fields in different groups are truly independent, so cases in which fields in different group have a non-vanishing propagator entry may give unexpected results.
-.Pp
 This entry is optional.
 .El
 .Pp
diff --git a/man/meantools-convert.1 b/man/meantools-convert.1
index fe021bb..4f8b191 100644
--- a/man/meantools-convert.1
+++ b/man/meantools-convert.1
@@ -1,5 +1,5 @@
 .Dd $Mdocdate: June 12 2015 $
-.Dt meantools-convert 1.2.1
+.Dt meantools-convert 1.3
 .Os
 .Sh NAME
 .Nm meantools-convert
diff --git a/man/meantools.1 b/man/meantools.1
index 30eb3ff..376a4c1 100644
--- a/man/meantools.1
+++ b/man/meantools.1
@@ -1,5 +1,5 @@
 .Dd $Mdocdate: April 14 2015 $
-.Dt meantools 1.2.1
+.Dt meantools 1.3
 .Os
 .Sh NAME
 .Nm meantools
diff --git a/man/numkondo.1 b/man/numkondo.1
index d1ae065..9506cf4 100644
--- a/man/numkondo.1
+++ b/man/numkondo.1
@@ -1,5 +1,5 @@
 .Dd $Mdocdate: April 14 2015 $
-.Dt numkondo 1.2.1
+.Dt numkondo 1.3
 .Os
 .Sh NAME
 .Nm numkondo
diff --git a/scripts/meantools-convert b/scripts/meantools-convert
index 603749e..ef87227 100755
--- a/scripts/meantools-convert
+++ b/scripts/meantools-convert
@@ -105,7 +105,7 @@ def latex_engine(argv,text):
             oneline=0
             i=i+1
 
-    return(convert_latex(text,lsym,Lsym,Csym,oneline,columns))
+    return(convert_latex(text,lsym,Lsym,Csym,oneline))
 
 # convert to C format
 def convert_C(text, lsym, Lsym, Csym, oneline):
@@ -160,7 +160,7 @@ def convert_C(text, lsym, Lsym, Csym, oneline):
     return(text+';')
 
 # convert to LaTeX format
-def convert_latex(text, lsym, Lsym, Csym, oneline, columns):
+def convert_latex(text, lsym, Lsym, Csym, oneline):
     # remove newlines
     if (oneline==0):
         text=text.replace('\n','\\\\\n')
diff --git a/src/definitions.cpp b/src/definitions.cpp
index 982f7a7..1884488 100644
--- a/src/definitions.cpp
+++ b/src/definitions.cpp
@@ -17,7 +17,7 @@ limitations under the License.
 #ifndef DEFINITIONS_GCC
 #define DEFINITIONS_GCC
 
-#define VERSION "1.2.1"
+#define VERSION "1.3"
 
 // number of entries in a configuration file
 #define ARG_COUNT 10
diff --git a/src/fields.c b/src/fields.c
index 1b221f2..bfd1f39 100644
--- a/src/fields.c
+++ b/src/fields.c
@@ -32,6 +32,7 @@ int init_Fields_Table(Fields_Table* fields){
   init_Identities(&((*fields).ids), FIELDS_SIZE);
   init_Symbols(&((*fields).symbols), FIELDS_SIZE);
   init_Int_Array(&((*fields).fermions),FIELDS_SIZE);
+  init_Int_Array(&((*fields).noncommuting),FIELDS_SIZE);
   return(0);
 }
 int free_Fields_Table(Fields_Table fields){
@@ -41,6 +42,7 @@ int free_Fields_Table(Fields_Table fields){
   free_Identities(fields.ids);
   free_Symbols(fields.symbols);
   free_Int_Array(fields.fermions);
+  free_Int_Array(fields.noncommuting);
   return(0);
 }
 
@@ -73,6 +75,16 @@ int is_fermion(int index, Fields_Table fields){
   }
 }
 
+// check whether a field is non-commuting
+int is_noncommuting(int index, Fields_Table fields){
+  if(int_array_find(abs(index), fields.noncommuting)>=0){
+    return(1);
+  }
+  else{
+    return(0);
+  }
+}
+
 
 // ------------------ Identities --------------------
 
@@ -180,13 +192,16 @@ int identities_concat(Identities input, Identities* output){
 
 // resolve the identities
 // requires both the monomials in polynomial and the ids in fields to be sorted
+// IMPORTANT: the sorting must be such that noncommuting fields must come before the other fields
 int resolve_ids(Polynomial* polynomial, Fields_Table fields){
   int i,j,k,l;
   int sign;
   int fermion_count;
+  int first_field;
   int at_least_one;
   int security;
-  Int_Array monomial;
+  Int_Array pre_monomial;
+  Int_Array post_monomial;
   Number num;
   Number tmp_num;
 
@@ -207,29 +222,38 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields){
 
       // loop over ids
       for(j=0;j<fields.ids.length;j++){
+
 	// check whether the monomial matches the id
-	if(int_array_is_subarray_ordered(fields.ids.lhs[j],(*polynomial).monomials[i])==1){
-	  init_Int_Array(&monomial, (*polynomial).monomials[i].length);
-	  
-	  // remove lhs from monomial
-	  // sign from moving the fields out of the monomial
+	first_field=int_array_is_subarray_noncommuting(fields.ids.lhs[j],(*polynomial).monomials[i],fields);
+	if(first_field>=0){
+	  init_Int_Array(&pre_monomial, (*polynomial).monomials[i].length);
+	  init_Int_Array(&post_monomial, (*polynomial).monomials[i].length);
+
+	  // add whatever is before the first field to pre
+	  for(k=0;k<first_field;k++){
+	    int_array_append((*polynomial).monomials[i].values[k],&pre_monomial);
+	  }
+
+	  // find others and move them together
+	  // sign from moving the fields
 	  sign=1;
-	  // number of Fermions to remove from the monomial
+	  // number of Fermions to jump over
 	  fermion_count=0;
-	  for(k=0,l=0;k<(*polynomial).monomials[i].length;k++){
+	  for(l=1,k=first_field+1;k<(*polynomial).monomials[i].length;k++){
 	    // check whether the field is identical to the "current" one in the id
 	    // if l is too large, then keep the field
 	    if(l>=fields.ids.lhs[j].length || (*polynomial).monomials[i].values[k]!=fields.ids.lhs[j].values[l]){
-	      int_array_append((*polynomial).monomials[i].values[k],&monomial);
-	      // sign correction
-	      if(fermion_count % 2 ==1 && is_fermion((*polynomial).monomials[i].values[k], fields)){
-		sign*=-1;
+	      // add to post
+	      int_array_append((*polynomial).monomials[i].values[k],&post_monomial);
+	      // count Fermions to jump
+	      if(is_fermion((*polynomial).monomials[i].values[k],fields)){
+		fermion_count++;
 	      }
 	    }
 	    else{
-	      // increment fermion_count
-	      if(is_fermion(fields.ids.lhs[j].values[l],fields)){
-		fermion_count++;
+	      // sign correction
+	      if(is_fermion(fields.ids.lhs[j].values[l],fields) && fermion_count % 2 == 1){
+		sign*=-1;
 	      }
 	      // increment "current" field in the id
 	      l++;
@@ -240,30 +264,33 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields){
 	  // add extra monomials (if there are more than 1)
 	  for(k=1;k<fields.ids.rhs[j].length;k++){
 	    number_prod(num, fields.ids.rhs[j].nums[k], &tmp_num);
-	    polynomial_append(monomial, (*polynomial).factors[i], tmp_num, polynomial);
+	    polynomial_append(pre_monomial, (*polynomial).factors[i], tmp_num, polynomial);
 	    free_Number(tmp_num);
 	    int_array_concat(fields.ids.rhs[j].monomials[k],(*polynomial).monomials+(*polynomial).length-1);
+	    int_array_concat(post_monomial,(*polynomial).monomials+(*polynomial).length-1);
 	    // re-sort monomial
 	    sign=1;
-	    monomial_sort((*polynomial).monomials[(*polynomial).length-1],0,(*polynomial).monomials[(*polynomial).length-1].length-1,fields,&sign);
+	    monomial_sort((*polynomial).monomials[(*polynomial).length-1],fields,&sign);
 	    number_Qprod_chain(quot(sign,1),(*polynomial).nums+(*polynomial).length-1);
 	  } 
 	  // correct i-th monomial
 	  free_Number((*polynomial).nums[i]);
 	  (*polynomial).nums[i]=number_prod_ret(num,fields.ids.rhs[j].nums[0]);
 	  free_Int_Array((*polynomial).monomials[i]);
-	  (*polynomial).monomials[i]=monomial;
+	  (*polynomial).monomials[i]=pre_monomial;
 	  int_array_concat(fields.ids.rhs[j].monomials[0],(*polynomial).monomials+i);
+	  int_array_concat(post_monomial,(*polynomial).monomials+i);
+	  free_Int_Array(post_monomial);
 	  // re-sort monomial
 	  sign=1;
-	  monomial_sort((*polynomial).monomials[i],0,(*polynomial).monomials[i].length-1,fields,&sign);
+	  monomial_sort((*polynomial).monomials[i],fields,&sign);
 	  number_Qprod_chain(quot(sign,1),(*polynomial).nums+i);
 
 	  // free num
 	  free_Number(num);
 
-	  // repeat the step (in order to perform all of the replacements if several are necessary)
-	  j--;
+	  // repeat the replacement (in order to perform all of the replacements if several are necessary)
+	  j=0;
 	  if(at_least_one==0){
 	    at_least_one=1;
 	  }
@@ -275,6 +302,62 @@ int resolve_ids(Polynomial* polynomial, Fields_Table fields){
   return(0);
 }
 
+// check whether an array is a sub-array of another
+// requires noncommuting elements to be next to each other
+// other elements may be separated, but the order must be respected
+// returns the first index of the sub-array
+// IMPORTANT: the noncommuting elements must precede all others in input and in test_array
+int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fields_Table fields){
+  int i,j;
+  int matches=0;
+  int post_nc=0;
+  int match_nc;
+  int first=-1;
+
+  // bound noncommuting elements
+  while(is_noncommuting(input.values[post_nc], fields)==1){
+    post_nc++;
+  }
+
+  for(i=0,match_nc=0;i<test_array.length;i++){
+    if(test_array.values[i]==input.values[0]){
+      match_nc=1;
+    }
+    for(j=1;j<post_nc;j++){
+      if(test_array.values[i+j]!=input.values[j]){
+	match_nc=0;
+      }
+    }
+    if(match_nc==1){
+      first=i;
+      break;
+    }
+  }
+
+  if(first<0){
+    return(-1);
+  }
+
+  if(post_nc>0){
+    matches=post_nc;
+  }
+  else{
+    matches=1;
+  }
+
+  for(i=first+1;i<test_array.length && matches<input.length;i++){
+    if(input.values[matches]==test_array.values[i]){
+      matches++;
+    }
+  }
+  if(matches==input.length){
+    return(first);
+  }
+  else{
+    return(-1);
+  }
+}
+
 
 // ------------------ Symbols --------------------
 
diff --git a/src/fields.h b/src/fields.h
index 3795d92..b6f05c2 100644
--- a/src/fields.h
+++ b/src/fields.h
@@ -29,6 +29,8 @@ int free_Fields_Table(Fields_Table fields);
 int field_type(int index, Fields_Table fields);
 // check whether a field anticommutes
 int is_fermion(int index, Fields_Table fields);
+// check whether a field is non-commuting
+int is_noncommuting(int index, Fields_Table fields);
 
 
 // init
@@ -51,6 +53,8 @@ int identities_concat(Identities input, Identities* output);
 
 // resolve the identities
 int resolve_ids(Polynomial* polynomial, Fields_Table fields);
+// check whether an array is a sub-array of another, support for noncommuting elements
+int int_array_is_subarray_noncommuting(Int_Array input, Int_Array test_array, Fields_Table fields);
 
 
 // init
diff --git a/src/kondo.c b/src/kondo.c
index 39d74bb..79b7c56 100644
--- a/src/kondo.c
+++ b/src/kondo.c
@@ -29,16 +29,17 @@ limitations under the License.
 #include "definitions.cpp"
 #include "rational.h"
 
-// dimension of A, B and h (must be <10)
+// dimension of A, B, h and t (must be <10)
 #define KONDO_DIM 3
 // number of spin components
 #define KONDO_SPIN 2
 
-// offsets for indices of A, B and h
+// offsets for indices of A, B, h and t
 // order matters for symbols table
 #define KONDO_A_OFFSET 1
 #define KONDO_B_OFFSET 2
 #define KONDO_H_OFFSET 3
+#define KONDO_T_OFFSET 4
 
 // parsing modes (from parse_file.c)
 #define PP_NULL_MODE 0
@@ -193,13 +194,19 @@ int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fiel
 
   //  parameters
   char_array_append_str("h:",str_fields);
+  // h
   for(i=0;i<KONDO_DIM;i++){
-    char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_H_OFFSET));
+    char_array_snprintf(str_fields, "%d,", 10*(i+10*KONDO_H_OFFSET));
+  }
+  // t
+  for(i=0;i<KONDO_DIM;i++){
+    char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_T_OFFSET));
     if(i<KONDO_DIM-1){
-      char_array_append(',',str_fields);
+      char_array_append(',', str_fields);
     }
   }
-  char_array_append('\n',str_fields);
+  char_array_append('\n', str_fields);
+
 
   // declare Fermions
   char_array_append_str("f:",str_fields);
@@ -226,6 +233,16 @@ int kondo_fields_table(int box_count, Char_Array* str_fields, Fields_Table* fiel
   }
   char_array_append('\n',str_fields);
 
+  // declare noncommuting
+  char_array_append_str("a:",str_fields);
+  for(i=0;i<KONDO_DIM;i++){
+    char_array_snprintf(str_fields, "%d", 10*(i+10*KONDO_T_OFFSET));
+    if(i<KONDO_DIM-1){
+      char_array_append(',',str_fields);
+    }
+  }
+  char_array_append('\n', str_fields);
+
   // parse fields table
   parse_input_fields(*str_fields, fields);
 
@@ -245,7 +262,7 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
 
   // loop over box index
   for(i=1;i<=box_count;i++){
-    // loop over letters
+    // loop over letters (A and B)
     for(j=0;j<2;j++){
       // loop over space dimension
       for(k=0;k<KONDO_DIM;k++){
@@ -255,7 +272,7 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
 	init_Char_Array(&tmp_str, 6);
 	char_array_snprintf(&tmp_str, "%c%d%d", letters[j], k, i);
 	// compute corresponding polynomial
-	kondo_resolve_ABh(tmp_str.str, &poly, *fields);
+	kondo_resolve_ABht(tmp_str.str, &poly, *fields);
 	free_Char_Array(tmp_str);
 	// write to output
 	polynomial_sprint(poly, str_symbols);
@@ -320,45 +337,6 @@ int kondo_symbols(Char_Array* str_symbols, int box_count, Fields_Table* fields){
   return(0);
 }
 
-// generate Kondo symbols (older method: one symbol for each scalar product)
-int kondo_symbols_scalarprod(Char_Array* str_symbols, int box_count, Fields_Table* fields){
-  int i,j,k;
-  Char_Array tmp_str;
-  Polynomial poly;
-  char letters[3]={'A','B','h'};
-
-  init_Char_Array(str_symbols, STR_SIZE);
-  char_array_snprintf(str_symbols, "#!symbols\n");
-
-  // loop over box index
-  for(i=1;i<=box_count;i++){
-    // loop over letters
-    for(j=0;j<3;j++){
-      for(k=0;k<3;k++){
-	// write index
-	char_array_snprintf(str_symbols, "%d=", 100*(10*(KONDO_A_OFFSET+j)+KONDO_A_OFFSET+k)+i);
-	// write the name of the scalar product
-	init_Char_Array(&tmp_str, 6);
-	char_array_snprintf(&tmp_str, "%c%d.%c%d", letters[j], i, letters[k], i);
-	// compute corresponding polynomial
-	kondo_resolve_scalar_prod(tmp_str.str, &poly, *fields);
-	free_Char_Array(tmp_str);
-	// write to output
-	polynomial_sprint(poly, str_symbols);
-	free_Polynomial(poly);
-	// add ,
-	if(i<box_count || j<2 || k<2){
-	  char_array_snprintf(str_symbols,",\n");
-	}
-      }
-    }
-  }
-
-  parse_input_symbols(*str_symbols, fields);
-
-  return(0);
-}
-
 
 // generate Kondo groups (groups of independent variables)
 int kondo_groups(Char_Array* str_groups, int box_count){
@@ -395,12 +373,21 @@ int kondo_groups(Char_Array* str_groups, int box_count){
 
 // generate Kondo identities
 // h_3^2=1-h_1^2-h_2^2
+// and Pauli matrices
 int kondo_identities(Char_Array* str_identities){
   int i;
 
   init_Char_Array(str_identities,STR_SIZE);
   char_array_snprintf(str_identities, "#!identities\n");
 
+  // Pauli matrices
+  for(i=0;i<KONDO_DIM;i++){
+    char_array_snprintf(str_identities,"[f%d][f%d]=(1),\n",10*(10*KONDO_T_OFFSET+i),10*(10*KONDO_T_OFFSET+i));
+    char_array_snprintf(str_identities,"[f%d][f%d]=(s{-1})[f%d],\n",10*(10*KONDO_T_OFFSET+i),10*(10*KONDO_T_OFFSET+(i+1)%3),10*(10*KONDO_T_OFFSET+(i+2)%3));
+    char_array_snprintf(str_identities,"[f%d][f%d]=((-1)s{-1})[f%d],\n",10*(10*KONDO_T_OFFSET+(i+2)%3),10*(10*KONDO_T_OFFSET+(i+1)%3),10*(10*KONDO_T_OFFSET+i));
+  }
+
+  // h
   char_array_snprintf(str_identities, "[f%d][f%d]=(1)",10*(KONDO_DIM-1+10*KONDO_H_OFFSET),10*(KONDO_DIM-1+10*KONDO_H_OFFSET));
   for(i=0;i<KONDO_DIM-1;i++){
     char_array_snprintf(str_identities, "+(-1)[f%d][f%d]",10*(i+10*KONDO_H_OFFSET),10*(i+10*KONDO_H_OFFSET));
@@ -809,11 +796,10 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_
 	  if(tmp_poly.length>0){
 	    for(i=0;i<tmp_poly.length;i++){
 	      if(mode==PP_FIELD_SCALAR_MODE){
-		if(offset1!=KONDO_H_OFFSET || offset2!=KONDO_H_OFFSET){
-		  int_array_append(1000*(10*offset1+offset2)+index, tmp_poly.monomials+i);
-		}
+		int_array_append(1000*(10*offset1+offset2)+index, tmp_poly.monomials+i);
 	      }
 	      else{
+		// vector product
 		int_array_append(100*(100*KONDO_A_OFFSET+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+index, tmp_poly.monomials+i);
 	      }
 	    }
@@ -822,28 +808,15 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_
 	  else{
 	    init_Int_Array(&tmp_monomial, MONOMIAL_SIZE);
 	    if(mode==PP_FIELD_SCALAR_MODE){
-	      if(offset1!=KONDO_H_OFFSET || offset2!=KONDO_H_OFFSET){
-		int_array_append(1000*(10*offset1+offset2)+index, &tmp_monomial);
-	      }
+	      int_array_append(1000*(10*offset1+offset2)+index, &tmp_monomial);
 	    }
 	    else{
+	      // vector product
 	      int_array_append(100*(100*KONDO_A_OFFSET+10*KONDO_B_OFFSET+KONDO_H_OFFSET)+index, &tmp_monomial);
 	    }
 	    init_Int_Array(&dummy_factor, 1);
 	    polynomial_append_noinit(tmp_monomial, dummy_factor, number_one(), &tmp_poly);
 	  }
-	  /* // older method in which a scalar product was expanded in A, B and h
-	  // resolve scalar product
-	  kondo_resolve_scalar_prod_symbols(buffer, &scalar_prod_poly);
-	  // add to tmp_poly
-	  if(tmp_poly.length==0){
-	    polynomial_concat(scalar_prod_poly,&tmp_poly);
-	  }
-	  else{
-	    polynomial_prod_chain(scalar_prod_poly,&tmp_poly,fields);
-	  }
-	  free_Polynomial(scalar_prod_poly);
-	  */
 	}
 	// switch back to null mode
 	mode=PP_NULL_MODE;
@@ -941,10 +914,10 @@ int parse_kondo_polynomial(Char_Array kondo_polynomial_str, Polynomial* polynomi
 }
 
 
-// read Aij, Bij, hi where i is a space dimension and j is a box index
-int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){
+// read Aij, Bij, hi, ti where i is a space dimension and j is a box index
+int kondo_resolve_ABht(char* str, Polynomial* output, Fields_Table fields){
   char* ptr;
-  // offset (A,B or H)
+  // offset (A,B, H or T)
   int offset=-1;
   // dimension
   int dim=-1;
@@ -978,6 +951,9 @@ int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){
     case 'h':
       offset=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset=KONDO_T_OFFSET;
+      break;
     default:
       // set index if dim was already set
       if(dim>=0){
@@ -1001,8 +977,8 @@ int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){
     }
   }
 
-  // h's
-  if(offset==KONDO_H_OFFSET){
+  // h's and t's
+  if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
     // external field
     init_Int_Array(&monomial,1);
     init_Int_Array(&factor,1);
@@ -1062,7 +1038,7 @@ int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields){
 // read a Kondo scalar product (generalized to vector products as well)
 int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields){
   char* ptr;
-  // offset of each term (A,B or H)
+  // offset of each term (A,B,H or T)
   int offset=-1;
   // index of each term (0,...,box_count)
   int index=0;
@@ -1091,6 +1067,9 @@ int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields
     case 'h':
       offset=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset=KONDO_T_OFFSET;
+      break;
 
     // scalar product
     case '.':
@@ -1195,8 +1174,8 @@ int kondo_polynomial_vector(int offset, int index, Polynomial (*polys)[3], Field
     init_Polynomial((*polys)+i,POLY_SIZE);
   }
 
-  // h's
-  if(offset==KONDO_H_OFFSET){
+  // h's and t's
+  if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
     // construct every component field
     for(i=0;i<KONDO_DIM;i++){
       // external field
@@ -1261,7 +1240,7 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){
   char* ptr;
   // first or second term
   int term=0;
-  // offset of each term (A,B or H)
+  // offset of each term (A,B,H or T)
   int offset[2];
   // index of each term (0,...,box_count)
   int index[2]={0,0};
@@ -1289,6 +1268,9 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){
     case 'h':
       offset[term]=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset[term]=KONDO_T_OFFSET;
+      break;
     // switch term
     case '.':
       term=1-term;
@@ -1304,13 +1286,13 @@ int kondo_resolve_scalar_prod_symbols(char* str, Polynomial* output){
     init_Int_Array(&monomial,2);
     init_Int_Array(&factor, 1);
 
-    if(offset[0]==KONDO_H_OFFSET){
+    if(offset[0]==KONDO_H_OFFSET || offset[0]==KONDO_T_OFFSET){
       int_array_append(10*(10*offset[0]+i)+index[0], &monomial);
     }
     else{
       int_array_append(100*(10*offset[0]+i)+index[0], &monomial);
     }
-    if(offset[1]==KONDO_H_OFFSET){
+    if(offset[1]==KONDO_H_OFFSET || offset[1]==KONDO_T_OFFSET){
       int_array_append(10*(10*offset[1]+i)+index[1], &monomial);
     }
     else{
@@ -1344,6 +1326,9 @@ int get_offset_index(char* str, int* offset, int* index){
     case 'h':
       *offset=KONDO_H_OFFSET;
       break;
+    case 't':
+      *offset=KONDO_T_OFFSET;
+      break;
     default:
       // char to int
       *index=*ptr-'0';
@@ -1378,6 +1363,9 @@ int get_offsets_index(char* str, int* offset1, int* offset2, int* index){
     case 'h':
       offset[term]=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset[term]=KONDO_T_OFFSET;
+      break;
     // switch term
     case '.':
       term=1-term;
@@ -1391,6 +1379,11 @@ int get_offsets_index(char* str, int* offset1, int* offset2, int* index){
   *offset1=offset[0];
   *offset2=offset[1];
 
+  // if no A's or B's, then index=0
+  if((offset[0]==KONDO_H_OFFSET || offset[0]==KONDO_T_OFFSET) && (offset[1]==KONDO_H_OFFSET || offset[1]==KONDO_T_OFFSET)){
+    *index=0;
+  }
+
   return(0);
 }
 
@@ -1429,6 +1422,9 @@ int get_symbol_index(char* str){
     case 'h':
       offset=KONDO_H_OFFSET;
       break;
+    case 't':
+      offset=KONDO_T_OFFSET;
+      break;
     default:
       // set index if dim was already set
       if(dim>=0){
@@ -1443,8 +1439,8 @@ int get_symbol_index(char* str){
   if(offset==-1){
     return(-1);
   }
-  // no symbol for h
-  if(offset==KONDO_H_OFFSET){
+  // no symbol for h or t
+  if(offset==KONDO_H_OFFSET || offset==KONDO_T_OFFSET){
     return(10*(10*offset+dim));
   }
   else{
diff --git a/src/kondo.h b/src/kondo.h
index 23756ad..c534145 100644
--- a/src/kondo.h
+++ b/src/kondo.h
@@ -55,7 +55,7 @@ int parse_kondo_polynomial_str(char* str_polynomial, Polynomial* output, Fields_
 int parse_kondo_polynomial(Char_Array kondo_polynomial_str, Polynomial* polynomial, Fields_Table fields);
 
 // read Aij, Bij, hi where i is a space dimension and j is a box index
-int kondo_resolve_ABh(char* str, Polynomial* output, Fields_Table fields);
+int kondo_resolve_ABht(char* str, Polynomial* output, Fields_Table fields);
 // read a Kondo scalar product
 int kondo_resolve_scalar_prod(char* str, Polynomial* output, Fields_Table fields);
 // compute a scalar product of polynomial vectors
diff --git a/src/mean.c b/src/mean.c
index aad7a5a..39ece56 100644
--- a/src/mean.c
+++ b/src/mean.c
@@ -42,7 +42,7 @@ int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Ma
   *out=polynomial_one();
 
   // sort first
-  monomial_sort(monomial, 0, monomial.length-1, fields, &sign);
+  monomial_sort(monomial, fields, &sign);
   polynomial_multiply_Qscalar(*out, quot(sign,1));
   // get internals
   // (*out).monomials is the first element of out but it only has 1 element
@@ -417,7 +417,7 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po
     if(check_monomial_match(tmp_monomial, fields)==1){
       // sort monomial
       sign=1;
-      monomial_sort(tmp_monomial, 0, tmp_monomial.length-1, fields, &sign);
+      monomial_sort(tmp_monomial, fields, &sign);
       number_Qprod_chain(quot(sign,1), &tmp_num);
 
       // mean
@@ -628,7 +628,7 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol
   init_Polynomial(output, MONOMIAL_SIZE);
 
   // check whether there are symbols
-  // requires the symbols to be at the end of the monomial
+  // IMPORTANT: the symbols must be at the end of the monomial
   if(monomial.length==0 || field_type(monomial.values[monomial.length-1], fields)!=FIELD_SYMBOL){
     // mean
     mean(monomial, &num_mean, fields, propagator);
@@ -639,7 +639,7 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol
     // sort into groups
     if(groups.length>0){
       sign=1;
-      monomial_sort_groups(monomial, 0, monomial.length-1, fields, groups, &sign);
+      monomial_sort_groups(monomial, fields, groups, &sign);
     }
     // construct groups and take mean
     init_Int_Array(&tmp_monomial, MONOMIAL_SIZE);
diff --git a/src/parse_file.c b/src/parse_file.c
index 6054372..19b91c6 100644
--- a/src/parse_file.c
+++ b/src/parse_file.c
@@ -44,16 +44,17 @@ limitations under the License.
 #define PP_EXTERNAL_MODE 8
 #define PP_INTERNAL_MODE 9
 #define PP_FERMIONS_MODE 10
+#define PP_NONCOMMUTING_MODE 11
 // indices
-#define PP_INDEX_MODE 11
+#define PP_INDEX_MODE 12
 // factors or monomials
-#define PP_BRACKET_MODE 12
+#define PP_BRACKET_MODE 13
 // labels
-#define PP_LABEL_MODE 13
+#define PP_LABEL_MODE 14
 // polynomial
-#define PP_POLYNOMIAL_MODE 14
+#define PP_POLYNOMIAL_MODE 15
 // group
-#define PP_GROUP_MODE 15
+#define PP_GROUP_MODE 16
 
 
 // parse fields list
@@ -101,6 +102,11 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
 	  mode=PP_FERMIONS_MODE;
 	}
 	break;
+      case 'a':
+	if(mode==PP_NULL_MODE){
+	  mode=PP_NONCOMMUTING_MODE;
+	}
+	break;
 
       // reset buffer
       case ':':
@@ -123,6 +129,9 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
 	else if(mode==PP_FERMIONS_MODE){
 	  int_array_append(i,&((*fields).fermions));
 	}
+	else if(mode==PP_NONCOMMUTING_MODE){
+	  int_array_append(i,&((*fields).noncommuting));
+	}
 	buffer_ptr=buffer;
 	*buffer_ptr='\0';
 	break;
@@ -142,6 +151,9 @@ int parse_input_fields(Char_Array str_fields, Fields_Table* fields){
 	else if(mode==PP_FERMIONS_MODE){
 	  int_array_append(i,&((*fields).fermions));
 	}
+	else if(mode==PP_NONCOMMUTING_MODE){
+	  int_array_append(i,&((*fields).noncommuting));
+	}
 	mode=PP_NULL_MODE;
 	break;
 
@@ -417,7 +429,7 @@ int parse_input_identities(Char_Array str_identities, Fields_Table* fields){
   (*fields).ids.length=0;
   for(i=0;i<tmp;i++){
     sign=1;
-    monomial_sort((*fields).ids.lhs[i], 0, (*fields).ids.lhs[i].length-1, *fields, &sign);
+    monomial_sort((*fields).ids.lhs[i], *fields, &sign);
     polynomial_simplify((*fields).ids.rhs+i, *fields);
     polynomial_multiply_Qscalar((*fields).ids.rhs[i],quot(sign,1));
   }
diff --git a/src/polynomial.c b/src/polynomial.c
index 639728a..d082b60 100644
--- a/src/polynomial.c
+++ b/src/polynomial.c
@@ -688,7 +688,7 @@ int polynomial_simplify(Polynomial* polynomial, Fields_Table fields){
   // sort monomials and factors
   for(i=0;i<(*polynomial).length;i++){
     sign=1;
-    monomial_sort((*polynomial).monomials[i],0,(*polynomial).monomials[i].length-1,fields,&sign);
+    monomial_sort((*polynomial).monomials[i],fields,&sign);
     number_Qprod_chain(quot(sign,1),(*polynomial).nums+i);
     int_array_sort((*polynomial).factors[i],0,(*polynomial).factors[i].length-1);
   }
@@ -787,7 +787,30 @@ int exchange_polynomial_terms(int i, int j, Polynomial* polynomial){
 }
 
 // sort a monomial (with sign coming from exchanging two Fermions)
-int monomial_sort(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign){
+// if the monomial contains noncommuting elements, put them at the beginning of the monomial
+int monomial_sort(Int_Array monomial, Fields_Table fields, int* sign){
+  int i,j;
+  int tmp;
+  // first index after noncommuting indices
+  int post_nc=0;
+
+  for(i=0;i<monomial.length;i++){
+    if(is_noncommuting(monomial.values[i], fields)){
+      tmp=monomial.values[i];
+      for(j=i;j>post_nc;j--){
+	monomial.values[j]=monomial.values[j-1];
+      }
+      monomial.values[post_nc]=tmp;
+      post_nc++;
+    }
+  }
+
+  monomial_sort_nonc(monomial, post_nc, monomial.length-1, fields, sign);
+
+  return(0);
+}
+// without noncommuting terms
+int monomial_sort_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign){
   int i;
   int index;
   // the pivot: middle of the monomial
@@ -812,8 +835,8 @@ int monomial_sort(Int_Array monomial, int begin, int end, Fields_Table fields, i
     exchange_monomial_terms(monomial, index, end, fields, sign);
 
     // recurse
-    monomial_sort(monomial, begin, index-1, fields, sign);
-    monomial_sort(monomial, index+1, end, fields, sign);
+    monomial_sort_nonc(monomial, begin, index-1, fields, sign);
+    monomial_sort_nonc(monomial, index+1, end, fields, sign);
   }
   return(0);
 }
@@ -872,7 +895,30 @@ int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table
 
 
 // sort a monomial by putting each group together
-int monomial_sort_groups(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign){
+// if the monomial contains noncommuting elements, put them at the beginning of the monomial
+int monomial_sort_groups(Int_Array monomial, Fields_Table fields, Groups groups, int* sign){
+  int i,j;
+  int tmp;
+  // first index after noncommuting indices
+  int post_nc=0;
+
+  for(i=0;i<monomial.length;i++){
+    if(is_noncommuting(monomial.values[i], fields)){
+      tmp=monomial.values[i];
+      for(j=post_nc;j<i;j++){
+	monomial.values[j+1]=monomial.values[j];
+      }
+      monomial.values[post_nc]=tmp;
+      post_nc++;
+    }
+  }
+
+  monomial_sort_groups_nonc(monomial, post_nc, monomial.length-1, fields, groups, sign);
+
+  return(0);
+}
+// without noncommuting terms
+int monomial_sort_groups_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign){
   int i;
   int index;
   // the pivot: middle of the monomial
@@ -897,8 +943,8 @@ int monomial_sort_groups(Int_Array monomial, int begin, int end, Fields_Table fi
     exchange_monomial_terms(monomial, index, end, fields, sign);
 
     // recurse
-    monomial_sort(monomial, begin, index-1, fields, sign);
-    monomial_sort(monomial, index+1, end, fields, sign);
+    monomial_sort_groups_nonc(monomial, begin, index-1, fields, groups, sign);
+    monomial_sort_groups_nonc(monomial, index+1, end, fields, groups, sign);
   }
   return(0);
 }
diff --git a/src/polynomial.h b/src/polynomial.h
index ec50227..6aff76d 100644
--- a/src/polynomial.h
+++ b/src/polynomial.h
@@ -99,14 +99,18 @@ int polynomial_sort(Polynomial* polynomial, int begin, int end);
 int exchange_polynomial_terms(int i, int j, Polynomial* polynomial);
 
 // sort a monomial (with sign coming from exchanging two Fermions)
-int monomial_sort(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign);
+int monomial_sort(Int_Array monomial, Fields_Table fields, int* sign);
+// without noncommuting terms
+int monomial_sort_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, int* sign);
 // order fields: parameter, external, internal
 int compare_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields);
 // exchange two fields (with sign)
 int exchange_monomial_terms(Int_Array monomial, int pos1, int pos2, Fields_Table fields, int* sign);
 
 // sort a monomial by putting each group together
-int monomial_sort_groups(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign);
+int monomial_sort_groups(Int_Array monomial, Fields_Table fields, Groups groups, int* sign);
+// without noncommuting terms
+int monomial_sort_groups_nonc(Int_Array monomial, int begin, int end, Fields_Table fields, Groups groups, int* sign);
 // order fields: group, then parameter, external, internal
 int compare_monomial_terms_groups(Int_Array monomial, int pos1, int pos2, Fields_Table fields, Groups groups);
 
diff --git a/src/types.h b/src/types.h
index 0452ebc..6b029d5 100644
--- a/src/types.h
+++ b/src/types.h
@@ -154,6 +154,8 @@ typedef struct Fields_Table{
   Symbols symbols;
   // list of anti-commuting variables (fields or symbols)
   Int_Array fermions;
+  // list of non-commuting variables (fields or symbols)
+  Int_Array noncommuting;
 } Fields_Table;
 
 // index labels