Here is a partly analytic solution done in Maple 2017.3.
First, the result of

```
f := -(1-x)*(1-y)*z*ln((1-x)*(1-y)*z/((1-x)*(1-y)*z+(1-x)*y*(1-z)+x*(1-y)*(1-z)))-
(1-x)*y*(1-z)*ln((1-x)*y*(1-z)/((1-x)*(1-y)*z+(1-x)*y*(1-z)+x*(1-y)*(1-z)))-
x*(1-y)*(1-z)*ln(x*(1-y)*(1-z)/((1-x)*(1-y)*z+(1-x)*y*(1-z)+x*(1-y)*(1-z)))-
(1-x)*y*z*ln((1-x)*y*z/((1-x)*y*z+x*(1-y)*z+x*y*(1-z)))-
x*(1-y)*z*ln(x*(1-y)*z/((1-x)*y*z+x*(1-y)*z+x*y*(1-z)))-
x*y*(1-z)*ln(x*y*(1-z)/((1-x)*y*z+x*(1-y)*z+x*y*(1-z))):
DirectSearch:-GlobalOptima(f, {x >= 0, y >= 0, z >= 0, x <= 1, y <= 1, z <= 1}, maximize);
```

$$[.823959216501083, [x = .500000007051357, y = .500000002907430, z = .499999998207894], 448]
$$
means the global maximum of $f$ on the cube $[0,1]^3$ is reached at the point
$x = .500000007051357, y = .500000002907430, z = .499999998207894$ with the absolute error
$10^{-6}$ (see the help to DirectSearch).
Second, the result of

```
DirectSearch:-SolveEquations([diff(f, x) = 0, diff(f, y) = 0, diff(f, z) = 0],
{x >= 0, y >= 0, z >= 0, x <= 1, y <= 1, z <= 1}, AllSolutions, solutions = 5);
```

$$ \left[ \begin {array}{cccc} { 1.23259516440783095\times 10^{-32}}&
\left[ \begin {array}{c} -{ 1.11022302462515654\times 10^{-16}}
\\ 0.0\\ 0.0\end {array}
\right] &[x= 0.500000000000000000,y= 0.499999999999999944,z=
0.499999999999999944]&154\end {array} \right]
$$
means there is only one critical point of $f$ inside the cube $[0,1]^3$ and this point is close to
$x=\frac 1 2,y= \frac 1 2,z=\frac 1 2$. It is easy to verify that the
point $x=\frac 1 2,y= \frac 1 2,z=\frac 1 2$ is a critical point of $f$.
Third, the results of

```
Student[MultivariateCalculus]:- SecondDerivativeTest(f, [x,y,z] = [1/2, 1/2, 1/2],output='hessian');
```

$$ \left[ \begin {array}{ccc} -8/3&-\ln \left( 3 \right) +4/3&-\ln
\left( 3 \right) +4/3\\ -\ln \left( 3 \right) +4/3
&-8/3&-\ln \left( 3 \right) +4/3\\-\ln \left( 3
\right) +4/3&-\ln \left( 3 \right) +4/3&-8/3\end {array} \right]
$$
and

```
Student[MultivariateCalculus]:- SecondDerivativeTest(f, [x, y, z] = [1/2, 1/2, 1/2]);
```

${\it LocalMin}=[],{\it LocalMax}=[[1/2,1/2,1/2]],{\it Saddle}=[]$

mean that the critical point $x=\frac 1 2,y= \frac 1 2,z=\frac 1 2$ is the maximum point.
Summing up the above, we conclude, that the global maximum is reached at.

7more comments