NF = 4
op = NewOperator("Number",NF,{1},{1},{0.1+I}) + NewOperator("U",NF,{0},{1},{5}) + 3
rho = {{0.7,0.3+I,0,0},{0.3-I,0.4,0,0},{0,0,0,0},{0,0,0,0}}

print("Full Operator:")
print(op)
print("\nDensity:")
print(rho)
print("\nMeanFieldOperator:")
print( MeanFieldOperator(op, rho) )
print("\nMeanFieldOperator with electron self-interaction:")
print( MeanFieldOperator(op, rho, {{"AddDFTSelfInteraction",true}}) )