For your first request, I modify my code and add a Print to show only the last two items.
v=Select[Tuples[{0,1},16],
Total[Take[#,{1,4}]]==...==1&];
Print[Drop[v,22]];
which displays
{{1,0,0,0,0,1,0,0,0,0,0,1,0,0,1,0},{1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1}}
which corresponds to
{{z0,z5,z11,z14},{z0,z5,z10,z15}}
Is that what you were asking for?
For your second item, do I correctly understand that you are removing the assignments giving values to D00, D11, D22 and D33? If I do that then do I make any changes to the assignments to TC0, TC1, TC2 and TC3 which use D00, D11, D22 and D33? Do I understand that I then use the 8 constraints that follow to create the 24 alternatives as before? If so then I believe the current Select matches your equations you have given for this.
I do not believe I understand your third item? Are there now 32 binary variables instead of 16?
Note: If it might make my code any easier to understand there are two items.
Mathematica starts list indexes at one, not zero. You named your variables starting at zero. So z0 corresponds to #[[1]], z1 corresponds to #[[2]], etc.
Next, you can replace
{P00,P11,P22,P33,D00,D11,D22,D33}=
{z0 P0+z1 P1+z2 P2+z3 P3, z4 P0+z5 P1+z6 P2+z7 P3, z8 P0+z9 P1+z10 P2+z11 P3,
z12 P0+z13 P1+z14 P2+z15 P3, z0 D0+z1 D1+z2 D2+z3 D3, z4 D0+z5 D1+z6 D2+z7 D3,
z8 D0+z9 D1+z10 D2+z11 D3, z12 D0+z13 D1+z14 D2+z15 D3}/.r[[i]];
with
P00=z0 P0+z1 P1+z2 P2+z3 P3/.r[[i]];
P11=z4 P0+z5 P1+z6 P2+z7 P3/.r[[i]];
P22=z8 P0+z9 P1+z10 P2+z11 P3/.r[[i]];
P33=z12 P0+z13 P1+z14 P2+z15 P3/.r[[i]];
D00=z0 D0+z1 D1+z2 D2+z3 D3/.r[[i]];
D11=z4 D0+z5 D1+z6 D2+z7 D3/.r[[i]];
D22=z8 D0+z9 D1+z10 D2+z11 D3/.r[[i]];
D33=z12 D0+z13 D1+z14 D2+z15 D3/.r[[i]];
if that is any easier for you to understand.
