Let A be an expansive dilation on ℝn and φ : ℝn × [0, ∞) → [0, ∞) an anisotropic Musielak–Orlicz function. Let HAφ(ℝn) be the anisotropic Hardy space of Musielak–Orlicz type defined via the grand maximal function. In this article, the authors establish its molecular characterization via the atomic characterization of HAφ(ℝn). The molecules introduced in this article have the vanishing moments up to order s and the range of s in the isotropic case (namely, A := 2In×n) coincides with the range of well-known classical molecules and, moreover, even for the isotropic Hardy space Hp(ℝn) with p ∈ (0, 1] (in this case, A := 2In×n, φ(x, t) := tp for all x ∈ ℝn and t ∈ [0, ∞)), this molecular characterization is also new. As an application, the authors obtain the boundedness of anisotropic Calderón–Zygmund operators from HAφ(ℝn) to Lφ(ℝn) or from HAφ(ℝn) to itself.