Prove $\frac{x+y+z}{3}+\frac{3}{\frac1x+\frac1y+\frac1z}\geq5\sqrt[3]{\frac{xyz}{16}}$

Since we have a homogeneous inequality, we may assume $xyz=1$ without loss of generality, then find the infimum of $AM(x,y,z)+HM(x,y,z)$ over the set $$D=\left\{(x,y,z)\in\mathbb{R}_{+}^3: xyz=1\right\}.$$ So we want: $$ \inf_{\mathbb{R}^+\times\mathbb{R}^+}\left(\frac{x+y+\frac{1}{xy}}{3}+\frac{3}{xy+\frac{1}{x}+\frac{1}{y}}\right)=\inf_{\mathbb{R}^+\times\mathbb{R}^+}\left(\frac{x^2 y+xy^2+1}{3xy}+\frac{3xy}{x^2y^2+x+y}\right).$$ Let $f(x,y)$ be the last function. It is a symmetric function, hence it is enough to study its behaviour over the set $0<y\leq x$. So, let $y=kx$ with $k\in(0,1]$. We have: $$ g_k(x)=f(x,kx) = \frac{1}{3kx^2}+\frac{(k+1)x}{3}+\frac{3kx}{1+k+k^2 x^3} $$ and disregarding the last (positive) term, by the AM-GM inequality: $$ g_k(x)\geq \frac{1}{3kx^2}+\frac{(k+1)x}{6}+\frac{(k+1)x}{6} \geq 3\sqrt[3]{\frac{(k+1)^2}{108k}}$$ where the last term is a decreasing function over $(0,1]$. We just need a little improvement of the last inequality to prove the original claim. So, let us study $f(x,y)$ under another perspective. Assuming $x+y=k$ we have: $$ f(x,y) = \frac{k+\frac{1}{xy}}{3}+\frac{3}{\frac{k}{xy}+xy}=\frac{k+\frac{1}{x(k-x)}}{3}+\frac{3}{\frac{k}{x(k-x)}+x(k-x)} $$ and by putting $x=\lambda k, y=(1-\lambda)k$ we are left with: $$ \frac{k+\frac{1}{k^2\lambda(1-\lambda)}}{3}+\frac{3}{\frac{1}{k\lambda(1-\lambda)}+k^2\lambda(1-\lambda)}=\color{red}{\frac{k+\frac{1}{k^2 L}}{3}+\frac{3}{\frac{1}{kL}+k^2 L}=h(L,k)}$$ where $L$ ranges over $\left(0,\frac{1}{4}\right)$ and $k$ ranges over $\mathbb{R}_+$. The stationary points of $h(L,k)$ for a fixed $k$ occur at $$ L \in\left\{\pm\frac{1}{\sqrt{2} k^{3/2}},\pm\frac{1}{\sqrt{5} k^{3/2}}\right\} $$ and $$ h\left(\frac{1}{\sqrt{2} k^{3/2}},k\right) = \frac{4}{3}\sqrt{\frac{2}{k}}+\frac{k}{3},\qquad h\left(\frac{1}{\sqrt{5} k^{3/2}},k\right) = \frac{5}{6}\sqrt{\frac{5}{k}}+\frac{k}{3},$$ so the AM-GM inequality proves the claim.

In facts, the original inequality is optimal and our approach shows that equality is attained by $$ \color{red}{(x,y,z) = \lambda\cdot (1,4,4)} $$ and cyclic shifts.