Both cubic and hexagonal phases of ZnS are modeled using interatomic potential based simulations and density functional theory. A new set of potential parameters is derived, showing improved behavior compared with the previous ones. Results obtained with this new potential model show very good agreement with those obtained with density functional theory calculations and with experimental results when available. To calculate crystal morphologies for both phases, we perform an extensive study of the surface energies. In the cubic phase we take into account all the nonpolar surfaces with Miller indexes 0, 1, 2, 3, and 4, and all the polar surfaces with indexes 0, 1, and 2. The nonpolar (110) surface is the most stable surface in this phase and entirely dominates the crystal morphology, which is a dodecahedron showing only the (110) surface and its equivalents. In the hexagonal phase we find that it is necessary to take into account polar surfaces to obtain the crystal morphology, which has a highly anisotropic, cylindrical-like shape, with nonpolar surfaces on the sides and polar surfaces closing the cylinder.