Question: Creation and annihilation operators in Sum with indefinite upper bond

I have several questions concerning manipulations with creation and annihilation operators. I would be grateful if  besides answering my questions you would give me some suggestions for improving the code.

Background of the physic problem
I have a crystal lattice with several sites per unit cell. Every site can holds one, two or none of electrons. Electron at every site can be at one of two states: spin up and spin down. Electrons are Fermi particles so anticommutative basis is used. Index 'i' (1<=i<=f) corresponds to i's lattice cell, index 'j' (1<=j<=3) corresponds to j's site within i'th lattice cell, index 'sigma' corresponds to spin of the electron (if sigma=1 then spin is down, if sigma=2 then spin is up).

Code

>restart;
>with Physics :
>Setup mathematicalnotation = true :
>Physics:-Setup anticommutativeprefix = Psi :

>#Below the annihilation operators are defined. They depend on indicies [i,j,sigma].
>am:=proc(i, j, sigma)
>local k: : posint :
>k:=6*(i-1)+2*(j-1)+sigma;
>Annihilation(Psi, k, notation = explicit);
>end proc:

>#Below the creation operators are defined. They depend onindicies [i, j, sigma].
>ap:=proc(i, j, sigma)
>local k :: posint :
>k:=6*(i-1)+2*(j-1)+sigma;
>Creation(Psi, k, notation = explicit);
>end proc:

>#Here the occupation number at state i, j, sigma is programmed
>nn:=proc( i, j, sigma)
>ap(i, j, sigma).am(i, j, sigma)
>end proc:

>#Now I want to define operator of total number of particles.
># I have tried the:
>#N:=sum( sum( sum( nn( i, j, sigma), sigma = 1 ..2), j = 1 ..3), i = 1 ..f) ;
>#or
>#N:=Sum( Sum( Sum( nn( i, j, sigma), sigma = 1 ..2), j = 1 ..3), i = 1 ..f) ;
>#but both give an error.

>#The only working way is
>f:=3;   # 'f' have to be defined explicitly
>N:=add( add( add( nn( i, j, sigma), sigma = 1 ..2), j = 1 ..3), i = 1 ..f) ;

Is there a way to define operator 'N' with indefinite 'f'?
Is there a way to have operators nn in output of N rather then a^+a?
Is there a tool for simplification of the operators?

I'd like to have something like
>simplify(nn(1,2,1).nn(1,2,1))=nn(1.2.1)

Thank you