This is my starting point (f being the function to maximize, g and h being the constraints, a and b being the lagrange multipliers and nk being the number of particles in level/energy k):

f(n0, ..., n6) = ln(6!/(n0!*n1!*n2!*n3!*n4!*n5!*n6!))

g(n0, ..., n6) = 0*n0 + 1*n1 + 2*n2 + 3*n3 + 4*n4 + 5*n5 + 6*n6 = 6

h(n0, ..., n6) = n0 + n1 + n2 + n3 + n4 + n5 + n6 = 6

gradient(f) = a*gradient(g) + b*gradient(h)

I know ln(gamma(x))' = digamma(x)