  A versatile and numerically inexpensive method is presented allowing the
accurate calculation of phase diagrams for bosonic lattice models. By treating
clusters within the Gutzwiller theory, a surprisingly good description of
quantum fluctuations beyond the mean-field theory is achieved approaching
quantum Monte-Carlo predictions for large clusters. Applying this powerful
method to the Bose-Hubbard model, we demonstrate that it yields precise results
for the superfluid to Mott-insulator transition in square, honeycomb, and cubic
lattices. Due to the exact treatment within a cluster, the method can be
effortlessly adapted to more complicated Hamiltonians in the fast progressing
field of optical lattice experiments. This includes state- and site-dependent
superlattices, large confined atomic systems and disordered potentials, as well
as various types of extended Hubbard models. Furthermore, the approach allows
an excellent treatment of systems with arbitrary filling factors. We discuss
the perspectives that allow for the computation of large, spatially-varying
lattices, low-lying excitations, and time evolution.
