@ecterrab

As I mentioned above I need to collect terms with the same operator part: that is

x1*ap[3].am[4]+x2*ap[3].am[4] -> (x1+x2)*ap[3].am[4]

In the previous message I showed that simplify(z,size) doesn't solve the problem.

Besides that I need to get the coefficients near the operator part:
A*ap[i1].ap[i2]. ... .ap[in].am[j1].am[j2]. ... .am[jm] -> A

Unfortunately, I didn't manage to find a way of using Physics:-Coefficients for this purpose. I tried numerous variations of the following:
......

z:=x1*ap[1].am[1]+x2*ap[2].am[2]:
ApAm1 := [seq(ap[j], j = 1 .. 2), seq(am[j], j = 1 .. 2)]; 
Coefficients(z, ApAm1, all);

......
but unsuccessfully.

More precisely, I need to compose a system of equations A[alpha]=0 where A[alpha] are coefficients in front of operator parts:
z:=sum_{\alpha}(A[alpha]*ap[i1].ap[i2]. ... .ap[in].am[j1].am[j2]. ... .am[jm]).

I tried the following:

restart; 
with(Physics); with(Library); Setup(mathematicalnotation = true); 
Physics:-Setup(anticommutativeprefix = psi); 

for i to 4 do 
ap[i] := Creation(psi, i, notation = explicit); 
am[i] := Annihilation(psi, i, notation = explicit) 
end do; 

system_creation := proc (z)::list; 
local i, x, syst; 
syst := {}: 
for i to nops(z) do 
x := op(1, op(i, z)):
syst := syst union{`if`(type(x, Or(symbol, numeric, `+`)), x, 1) = 0} :
end do; 
syst 
end proc: 

z := D1*ap[1].am[2].ap[4].am[3]+A*ap[2].am[2].ap[1].am[1]+B*ap[1].am[1].ap[2].am[2]+B*ap[3].am[3].ap[4].am[4]+ap[1].(D2*am[2].ap[4].am[3])+3*ap[1].am[1]+R*ap[1].am[1]+ap[2].am[2]+(K[1]+L)*ap[3]/(C+5*D).am[4]; 

z2 := system_creation(z)

but the result is




instead of



I understand why it is so, but I can't get the correct result. It seems that Physics Package lack the tool for getting the coefficients of the operator expression. Besides that it whould be great if you introduce Physics:-Library:-Collect (the tool that collects terms with the same operator part).

Unfortunately, simplify/size doesn't always allow to factor out common factors.

restart; with(Physics); with(Library); Setup(mathematicalnotation = true); 
Physics:-Setup(anticommutativeprefix = psi); 
Physics:-Version(); 

for i to 2 do 
ap[i] := Creation(psi, i, notation = explicit); 
am[i] := Annihilation(psi, i, notation = explicit) 
end do; 

type(ap, Library:-PhysicsType:-ExtendedQuantumOperator); 
3*ap[1].am[1]+ap[1].am[1];
simplify(%, size);

Result:




The result is ok.

But!

z2 := 3*ap[1].am[1]+R*ap[1].am[1]+K/(C+D)*ap[2]).am[2]; 
z2 := simplify(z2, size)

Result:



That is wrong. The correct result is

Unfortunately, simplify/size doesn't always allow to factor out common factors.

restart; with(Physics); with(Library); Setup(mathematicalnotation = true); 
Physics:-Setup(anticommutativeprefix = psi); 
Physics:-Version(); 

for i to 2 do 
ap[i] := Creation(psi, i, notation = explicit); 
am[i] := Annihilation(psi, i, notation = explicit) 
end do; 

type(ap, Library:-PhysicsType:-ExtendedQuantumOperator); 
3*ap[1].am[1]+ap[1].am[1];
simplify(%, size);

Result:




The result is ok.

But!

z2 := 3*ap[1].am[1]+R*ap[1].am[1]+K/(C+D)*ap[2]).am[2]; 
z2 := simplify(z2, size)

Result:



That is wrong. The correct result is

SortProducts works better but still not properly.

restart; 
with(Physics): with(Library): Setup(mathematicalnotation = true): 
Physics:-Setup(anticommutativeprefix = psi):
Physics:-Version(); 

for i to 4 do 
ap[i] := Creation(psi, i, notation = explicit); 
am[i] := Annihilation(psi, i, notation = explicit) 
end do:

ApAm1 := [seq(ap[j], j = 1 .. 4), seq(am[j], j = 1 .. 4)]: 

z := D*ap[2].am[2].am[1].ap[1]; 
z := SortProducts(z, ApAm1, useanticommutator);

Result:







But I would like to see

SortProducts works better but still not properly.

restart; 
with(Physics): with(Library): Setup(mathematicalnotation = true): 
Physics:-Setup(anticommutativeprefix = psi):
Physics:-Version(); 

for i to 4 do 
ap[i] := Creation(psi, i, notation = explicit); 
am[i] := Annihilation(psi, i, notation = explicit) 
end do:

ApAm1 := [seq(ap[j], j = 1 .. 4), seq(am[j], j = 1 .. 4)]: 

z := D*ap[2].am[2].am[1].ap[1]; 
z := SortProducts(z, ApAm1, useanticommutator);

Result:







But I would like to see

First of all, sorry for unclear question. I'll try to formulate my problem without ambiguities.

The result is a bit embarrassing to me. Algebra rules for fermionic creation and annihilation operators are {a_i,a_j^+}=\delta_{i,j} and {a_i,a_j}=0 ({x,y} is anticommutator of x and y). Maple shows only {a_i,a_i^+}=1 but I suppose that all other rules are implied.

Now I restrict myself only with consideration of the case ApAm1. I expected to see the transformation of the term B*ap[1].am[1].ap[2].am[2] into -B*ap[1].ap[2].am[1].am[2], not into B*ap[1].ap[2].am[1].am[2] as Maple shows. The same concerns the third term: C*ap[3].am[3].ap[4].am[4]=-C*ap[3].ap[4].am[3].am[4].

It is either a bug or I missed something.

Besides that I expected Maple (Physic package) to move creation/anihilation operators step by step by one position:

A*ap[2].am[2].ap[1].am[1]=A*ap[2].(-ap[1].am[2]+{am[2].ap[1]}).am[1]=-A*ap[2].ap[1].am[2].am[1]=-A*(-ap[1].ap[2]+{ap[1],ap[2]}).am[2].am[1]=A*ap[1].ap[2].am[2].am[1]=A*ap[1].ap[2].(-am[1].am[2]+{am[1],am[2]})=-A*ap[1].ap[2].am[1].am[2]

Finally A*ap[2].am[2].ap[1].am[1]->-A*ap[1].ap[2].am[1].am[2]

If we have ap[2].am[2].am[1].ap[1] I would be happy to see

ap[2].am[2].am[1].ap[1]=ap[2].am[2].(-ap[1].am[1]+1)= ... =ap[1].ap[2].am[1].am[2]+ap[2].am[2]

First of all, sorry for unclear question. I'll try to formulate my problem without ambiguities.

The result is a bit embarrassing to me. Algebra rules for fermionic creation and annihilation operators are {a_i,a_j^+}=\delta_{i,j} and {a_i,a_j}=0 ({x,y} is anticommutator of x and y). Maple shows only {a_i,a_i^+}=1 but I suppose that all other rules are implied.

Now I restrict myself only with consideration of the case ApAm1. I expected to see the transformation of the term B*ap[1].am[1].ap[2].am[2] into -B*ap[1].ap[2].am[1].am[2], not into B*ap[1].ap[2].am[1].am[2] as Maple shows. The same concerns the third term: C*ap[3].am[3].ap[4].am[4]=-C*ap[3].ap[4].am[3].am[4].

It is either a bug or I missed something.

Besides that I expected Maple (Physic package) to move creation/anihilation operators step by step by one position:

A*ap[2].am[2].ap[1].am[1]=A*ap[2].(-ap[1].am[2]+{am[2].ap[1]}).am[1]=-A*ap[2].ap[1].am[2].am[1]=-A*(-ap[1].ap[2]+{ap[1],ap[2]}).am[2].am[1]=A*ap[1].ap[2].am[2].am[1]=A*ap[1].ap[2].(-am[1].am[2]+{am[1],am[2]})=-A*ap[1].ap[2].am[1].am[2]

Finally A*ap[2].am[2].ap[1].am[1]->-A*ap[1].ap[2].am[1].am[2]

If we have ap[2].am[2].am[1].ap[1] I would be happy to see

ap[2].am[2].am[1].ap[1]=ap[2].am[2].(-ap[1].am[1]+1)= ... =ap[1].ap[2].am[1].am[2]+ap[2].am[2]

@ecterrab 

1) Simplify works perfect
2) As concerns summation dummy parameter. It is not just a curiosity. My final program will be cumbersome and it would be useful to have a sum with 'safe' dummy parameter. 'Safe' means that I am allowed to forget, not to keep in mind the internal structure of constructions containing sum. However it causes the lost of some functionality. I decided not to restrict my code in functionality. I will avoid assigning the summation dummy.

Looking forward for update where Annihilation and Creation operators can accept non-numerical quantum-number-positions returning unevaluated.

@ecterrab 

1) Simplify works perfect
2) As concerns summation dummy parameter. It is not just a curiosity. My final program will be cumbersome and it would be useful to have a sum with 'safe' dummy parameter. 'Safe' means that I am allowed to forget, not to keep in mind the internal structure of constructions containing sum. However it causes the lost of some functionality. I decided not to restrict my code in functionality. I will avoid assigning the summation dummy.

Looking forward for update where Annihilation and Creation operators can accept non-numerical quantum-number-positions returning unevaluated.

@Alejandro Jakubi 

in 2D only this works:

help('`if`');

@Alejandro Jakubi 

in 2D only this works:

help('`if`');

The Physics package is of the latest version:



Ok. The following code works:

restart;
with(Physics):
Setup(mathematicalnotation = true):
Physics:-Setup(anticommutativeprefix = psi):
i := 3;

am[i] := Annihilation(psi, i):
ap[i] := Creation(psi, i):
n := ap[i].am[i]:
z := n.n:

Simplify(z);
Simplify(am[i]^2);
Simplify(ap[i].ap[i]);

Output:


   0
   0

But this one doesn't:

restart;
with(Physics):
Setup(mathematicalnotation = true):
Physics:-Setup(anticommutativeprefix = psi):
 
am := proc (i)
if [i]::(Not(list(posint))) then return 'procname(args)' end if;
Annihilation(psi, i, notation = explicit)
end proc:

ap := proc (i)
if [i]::(Not(list(posint))) then return 'procname(args)' end if;
Creation(psi, i, notation = explicit)
end proc:

nn := proc (i)
if [i]::(Not(list(posint))) then return 'procname(args)' end if;
ap(i).am(i)
end proc:

Setup(op = {am, ap, nn}):

i := 3:

z := nn(i).nn(i);
Simplify(z);
Simplify(am(i).am(i));


Output:

Error, (in unknown) invalid subscript selector
Error, (in unknown) invalid subscript selector


I don't know how to workaround this problem.

The Physics package is of the latest version:



Ok. The following code works:

restart;
with(Physics):
Setup(mathematicalnotation = true):
Physics:-Setup(anticommutativeprefix = psi):
i := 3;

am[i] := Annihilation(psi, i):
ap[i] := Creation(psi, i):
n := ap[i].am[i]:
z := n.n:

Simplify(z);
Simplify(am[i]^2);
Simplify(ap[i].ap[i]);

Output:


   0
   0

But this one doesn't:

restart;
with(Physics):
Setup(mathematicalnotation = true):
Physics:-Setup(anticommutativeprefix = psi):
 
am := proc (i)
if [i]::(Not(list(posint))) then return 'procname(args)' end if;
Annihilation(psi, i, notation = explicit)
end proc:

ap := proc (i)
if [i]::(Not(list(posint))) then return 'procname(args)' end if;
Creation(psi, i, notation = explicit)
end proc:

nn := proc (i)
if [i]::(Not(list(posint))) then return 'procname(args)' end if;
ap(i).am(i)
end proc:

Setup(op = {am, ap, nn}):

i := 3:

z := nn(i).nn(i);
Simplify(z);
Simplify(am(i).am(i));


Output:

Error, (in unknown) invalid subscript selector
Error, (in unknown) invalid subscript selector


I don't know how to workaround this problem.

@ecterrab 

Thank you for explanations. I understood everything.

A couple of thoughts.

Physics:-Version():



1) A way of workingaround of the problem with index 'i' within sum.

One can make the variables i, j, sigma to be local:

N := proc (f)
local i, j, sigma;
sum(add(add(nn(i, j, sigma), sigma = 1 .. 2), j = 1 .. 3), i = 1 .. f)
end proc:

instead of

N:=sum(add(add(nn(i, j, sigma), sigma = 1 .. 2), j = 1 .. 3), i = 1 .. f):

Now using of i, j, k in body of the program can be done without warnings.

However it kill the tricks you demonstrated:

......

K:=Ket(psi,1,0,1,1,0,0,0,1,1,0,0,0):
NK:=N(f).K:
summand:=op(1,NK):
nn_operands := map2(op, 1, summand);
[seq(nn_operands, i = 1 .. 2)];

gives




instead of





And some functionality lost. It is clear why it happens.

2) Your wrote:
"...However, because these am, ap and nn defined above are returning unevaluated when any of i, j, sigma are not numbers, the internal routines have no way to know these am, ap, nn programs are quantum operators..."

Look at the following situation:

am1 := proc (i)
if [i]::(Not(list(nonnegint))) then return 'procname(args)' end if;
Annihilation(psi, i, notation = explicit)
end proc:
type(am1, Library:-PhysicsType:-ExtendedQuantumOperator);

result:

false

as expected. Now:

am2 := proc (i)
Annihilation(psi, i, notation = explicit)
end proc:
type(am2, Library:-PhysicsType:-ExtendedQuantumOperator);

result:

false

again. However according to your words I expected to see 'true'.

Another example. If we previously declared


then

N := sum(add(add(nn(i, j, sigma), sigma = 1 .. 2), j = 1 .. 3), i = 1 .. f):
type(N, Library:-PhysicsType:-ExtendedQuantumOperator);

gives true

but

N:=f->sum(add(add(nn(i, j, sigma), sigma = 1 .. 2), j = 1 .. 3), i = 1 .. f):
type(N, Library:-PhysicsType:-ExtendedQuantumOperator);

surprisingly gives false. Of couse, one can easy fix it by Setup(op=N) but it would be more comfortable if Maple (Physics library) do it itself.

3) Unfortunately 'simplify' doesn't work:

restart;
with(Physics);
Setup(mathematicalnotation = true, anticommutativeprefix = psi):

am := Annihilation(psi, notation = explicit):
z := am.am:
simplify(z);

gives

@ecterrab 

Thank you for explanations. I understood everything.

A couple of thoughts.

Physics:-Version():



1) A way of workingaround of the problem with index 'i' within sum.

One can make the variables i, j, sigma to be local:

N := proc (f)
local i, j, sigma;
sum(add(add(nn(i, j, sigma), sigma = 1 .. 2), j = 1 .. 3), i = 1 .. f)
end proc:

instead of

N:=sum(add(add(nn(i, j, sigma), sigma = 1 .. 2), j = 1 .. 3), i = 1 .. f):

Now using of i, j, k in body of the program can be done without warnings.

However it kill the tricks you demonstrated:

......

K:=Ket(psi,1,0,1,1,0,0,0,1,1,0,0,0):
NK:=N(f).K:
summand:=op(1,NK):
nn_operands := map2(op, 1, summand);
[seq(nn_operands, i = 1 .. 2)];

gives




instead of





And some functionality lost. It is clear why it happens.

2) Your wrote:
"...However, because these am, ap and nn defined above are returning unevaluated when any of i, j, sigma are not numbers, the internal routines have no way to know these am, ap, nn programs are quantum operators..."

Look at the following situation:

am1 := proc (i)
if [i]::(Not(list(nonnegint))) then return 'procname(args)' end if;
Annihilation(psi, i, notation = explicit)
end proc:
type(am1, Library:-PhysicsType:-ExtendedQuantumOperator);

result:

false

as expected. Now:

am2 := proc (i)
Annihilation(psi, i, notation = explicit)
end proc:
type(am2, Library:-PhysicsType:-ExtendedQuantumOperator);

result:

false

again. However according to your words I expected to see 'true'.

Another example. If we previously declared


then

N := sum(add(add(nn(i, j, sigma), sigma = 1 .. 2), j = 1 .. 3), i = 1 .. f):
type(N, Library:-PhysicsType:-ExtendedQuantumOperator);

gives true

but

N:=f->sum(add(add(nn(i, j, sigma), sigma = 1 .. 2), j = 1 .. 3), i = 1 .. f):
type(N, Library:-PhysicsType:-ExtendedQuantumOperator);

surprisingly gives false. Of couse, one can easy fix it by Setup(op=N) but it would be more comfortable if Maple (Physics library) do it itself.

3) Unfortunately 'simplify' doesn't work:

restart;
with(Physics);
Setup(mathematicalnotation = true, anticommutativeprefix = psi):

am := Annihilation(psi, notation = explicit):
z := am.am:
simplify(z);

gives