The Atacama Large Millimeter/Submillimeter Array has revolutionized the field of dust polarization in protoplanetary disks across multiple wavelengths. Previous observations and empirical modeling have suggested multiple mechanisms of dust polarization toward HL Tau, including grain alignment and dust scattering. However, a detailed modeling of dust polarization based on grain alignment physics is not yet available. Here, using an updated POLARIS code, we perform numerical modeling of dust polarization arising from both grain alignment by the magnetically enhanced radiative torque mechanism and self-scattering to reproduce the HL Tau polarization observed at three wavelengths of 0.87, 1.3, and 3.1 mm. Our modeling results show that the observed multiwavelength polarization could be reproduced only when large grains contain embedded iron inclusions and those with slow internal relaxation must have wrong internal alignment (i.e., the grain’s major axis parallel to its angular momentum). The abundance of iron embedded inside grains in the form of clusters is constrained to be ≳16%, and the number of iron atoms per cluster is N
cl ∼ 9 × 102. Maximum grain sizes probed at wavelengths of λ = 0.87, 1.3, and 3.1 mm are constrained at ∼60, 80, and 90 μm, respectively.