Thanks for answering, but actually the I::ImaginaryI does not work.
Moreover, given a combination of tensors XXX,YYY and 2 parameter a,b, it does not recognise that
(XXX) a (YYY) b 
is the same as 
(XXX) a b (YYY) 
and
a b (XXX) (YYY)
etc...
check the final expression of 
 {a,b,c,d,k,l,m,n,p,q,r,s,u,v,e,f,g,h,i,z}::Indices(position=free);
    {a,b,c,d,k,l,m,n,p,q,r,s,u,v,e,f,g,h,i,z}::Integer(0..2). 
    \eta_{a b}::Metric(signature=1). 
    \delta_{m n}::KroneckerDelta. 
    \eta_{a? b?}::Symmetric. 
    \eta_{a? b?}::InverseMetric. 
    \eta_{m}^{n}::KroneckerDelta. 
    \eta^{m}_{n}::KroneckerDelta. 
    x::Coordinate. 
    \Gamma{#}::GammaMatrix(metric=\eta).
    x::Coordinate; 
    \nabla{#}::Derivative; 
    F_{a b}::AntiSymmetric. 
    F_{a b}::Depends(x). 
    F_{a b}::Depends(\nabla{#}). 
    F1_{a b}::AntiSymmetric. 
    F1_{a b}::Depends(x). 
    F1_{a b}::Depends(\nabla{#}). 
    H_{a b}::AntiSymmetric. 
    H_{a b}::Depends(x). 
    H_{a b}::Depends(\nabla{#}).
ex:= -3 (ricci / 12 + 2 A**2 H_{m n} H^{m n})(2 F_{a b} F^{a b} 
+ 16 A I N R**(-1) H_{p q} F^{p q}  
- 32 A**2 N**2 R**(-2) H_{p q} H^{p q}  ) 
+ 2 (ricci / 12 + 2 A**2 H_{m n} H^{m n})  F_{a b} F^{a b}   
+ ( I / 2 F^{c d} - 2 A N R**(-1) H^{c d}) ( -8 I F^{a b} (1/4 R_{a b c d} 
+ 8  A**2 H_{a c} H_{b d}) 
- 2 I  \nabla^{a}{\nabla_{a}{F_{c d}}} 
+ 8 A N R**(-1) \nabla^{a}{\nabla_{a}{ H_{c d} }} )
- 64 A**2 ( ( I / 2 F^{c d} - 2 A N R**(-1) H^{c d}) ( I / 2 F_{c d}
 - 2 A N R**(-1) H_{c d}) H_{a b} H^{a b} 
 + ( I / 2 F^{c d} H_{c d} - 2 A N R**(-1) H_{c d} H^{c d})  ( I / 2 F^{a b} H_{a b} 
- 2 A N R**(-1) H_{a b} H^{a b}) )
+ 4/15  F^{b}_{n} F^{n d} R^{a}_{b a d} - 2/15 F_{a b} F_{c d} R^{a b c d});
    distribute(_); 
    expand_power(_); 
    substitute(_, $A -> I*R/8$ ); 
    substitute(_, $I*I -> -1$ ); 
    substitute(_, $I*I -> -1$ ); 
    collect_factors(_);
    substitute(_, $F_{a b} -> Q F1_{a b} - N H_{a b}$ ); 
    distribute(_); 
    unwrap(_); 
    collect_factors(_);
    substitute(_, $N**2 -> W$ ); 
    substitute(_, $N -> 0$ ); 
    substitute(_, $W -> N**2$ ); 
    collect_factors(_); 
    rename_dummies(_); 
    canonicalise(_);
    # 3d property substitute(_, $H^{a b} H^{c d} H_{a c} H_{b d} -> 1/2 H^{a b} H_{a b} H^{c d} H_{c d}  $ ); 
    collect_factors(_); 
    rename_dummies(_); 
    canonicalise(_);
    # use the trick with commutators to rewrite nablas as ricci/Riemann (use F1,H EOM's and discrard total derivatives) 
    substitute(_, $ F1^{c d} \nabla^{a}{\nabla_{a}{F1_{c d}}} -> -  ( - 2 R^{a}_{b a c} F1^{b m} F1^{c}_{m} + R_{a b c d} F1^{a b} F1^{c d}  )  $ ); 
    substitute(_, $ H^{c d} \nabla^{a}{\nabla_{a}{H_{c d}}} -> -  ( - 2 R^{a}_{b a c} H^{b m} H^{c}_{m} + R_{a b c d} H^{a b} H^{c d}  )  $ ); 
    distribute(_); 
    collect_factors(_); 
    rename_dummies(_); 
    canonicalise(_); 
    collect_factors(_);
    substitute(_, $F1^{a b} F1_{a b} -> F1^2$ ); 
    substitute(_, $H^{a b} H_{a b} -> H^2$ );
    collect_factors(_); 
    rename_dummies(_); 
    canonicalise(_); 
    substitute(_, $H^{a b} H_{a b} -> H^2$ );
    canonicalise(_);