By applying the least squares solution for the middle wave analysis inside the Riemann fan, we construct a four-state Harten-Lax-van Leer (HLL) Riemann solver for numerical simulation of magneto-hydrodynamics (MHD). First, we revisit the two-state HLL scheme and obtain the two outer intermediate states across the two out-bounding fast waves through Rankine-Hugoniot (R-H) conditions; Second, the two inner intermediate states are calculated using a geometric interpretation of the R-H conditions across the middle contact wave. This newly constructed four-state HLL solver contains different wave structures from those of the HLLD Riemann solver; namely, the two Alfvén waves are replaced as the two combination waves originated from the merging of Alfvén and slow waves inside the Riemann fan. As we tested, this solver resolves the MHD discontinuities well, and has better capture ability than the HLLD solver for the slow waves, although it appears more diffusive than the latter in the situations where the slow waves are not solely generated. Overall, the new solver has the similar accuracy as the HLLD solver, thus it is suitable for the calculation of numerical fluxes for the Godunov-type numerical simulation of MHD equations where the slow waves are expected to be resolved.