  In this work we develop an implementation of the Wang\--Landau algorithm
[Phys. Rev. Lett. \textbf{86}, 2050-2053 (2001)]. This algorithm allows us to
find the density of states (DOS), a function that, for a given system,
describes the proportion of states that have a certain energy. The
implementation uses the Python language for the algorithm itself, and it can
take advantage of any library, such as the powerful LAMMPS library, for the
computation of energy. Therefore, the resulting implementation is simple and
flexible without sacrificing efficiency. This implementation also considers
recent developments in the parallelization of the code for faster computation.
We establish the soundness and effectiveness of our implementation by studying
well-known systems such as the Ising model, the Lennard\--Jones and EAM solids.
We have found that our implementation can find the DOS with very good precision
in a reasonable amount of time. Therefore, we are equipped with a very powerful
and flexible implementation that can be easily used in order to study more
realistic models of matter.
