NF = 4

rho = {{0.8,0.3,0,0},{0.3,0.7,0,0},{0,0,0.4,0.5*I},{0,0,-0.5*I,0.2}}

print("Rho:")
print(rho)

op1 = NewOperator("Number",NF,{1},{1},{0.1+I}) + NewOperator("U",NF,{0},{1},{5}) + 3
op2 = NewOperator("Number",NF,{1},{1},{0.00001 + I}) + NewOperator("Number",NF,{2},{2},{1}) + NewOperator("U",NF,{0},{1},{0},{1},{2},{3},{2},{3},{3}) + NewOperator("U",NF,{0},{1},{5})

print("\n\n\n")
print("op1:")
print(op1)
print("PartialMeanFieldOperator(op1, rho, {0,1})")
print(PartialMeanFieldOperator(op1, rho, {0,1}))
print("\n\n\n")
print("op2:")
print(op2)
print("PartialMeanFieldOperator(op2, rho, {0,1})")
print(PartialMeanFieldOperator(op2, rho, {0,1}))
print("PartialMeanFieldOperator(op2, rho, {0,1}, {{\"AddDFTSelfInteraction\",true}})")
print(PartialMeanFieldOperator(op2, rho, {0,1}, {{"AddDFTSelfInteraction",true}}))
print("PartialMeanFieldOperator(op2, rho, {2})")
print(PartialMeanFieldOperator(op2, rho, {2}))