Let n ≥ 2 and Ω be a bounded Lipschitz domain of R n . Assume that L R is a second-order divergence form elliptic operator having real-valued, bounded, symmetric, and measurable coefficients on L 2 (Ω) with the Robin boundary condition. In this article, via first obtaining the Hölder estimate of the heat kernels of L R , the authors establish a new atomic characterization of the Hardy space H L R p (Ω) associated with L R . Using this, the authors further show that, for any given p ∈ (n n + δ 0 , 1 ] , H z p (Ω) + L ∞ (Ω) = H L N p (Ω) = H L R p (Ω) ⫋ H L D p (Ω) = H r p (Ω) , where H L D p (Ω) and H L N p (Ω) denote the Hardy spaces on Ω associated with the corresponding elliptic operators respectively having the Dirichlet and the Neumann boundary conditions, H z p (Ω) and H r p (Ω) respectively denote the "supported type" and the "restricted type" Hardy spaces on Ω , and δ 0 ∈ (0 , 1 ] is the critical index depending on the operators L D , L N , and L R . The authors then obtain the boundedness of the Riesz transform ∇ L R - 1 / 2 on the Lebesgue space L p (Ω) when p ∈ (1 , ∞) [if p > 2 , some extra assumptions are needed] and its boundedness from H L R p (Ω) to L p (Ω) when p ∈ (0 , 1 ] or to H r p (Ω) when p ∈ (n n + 1 , 1 ] . As applications, the authors further obtain the global regularity estimates, in L p (Ω) when p ∈ (0 , p 0) and in H r p (Ω) when p ∈ (n n + 1 , 1 ] , for the inhomogeneous Robin problem of L R on Ω , where p 0 ∈ (2 , ∞) is a constant depending only on n, Ω , and the operator L R . The main novelties of these results are that the range (0 , p 0) of p for the global regularity estimates in the scale of L p (Ω) is sharp and that, in some sense, the space X : = H L R 1 (Ω) is also optimal to guarantee both the boundedness of ∇ L R - 1 / 2 from X to L 1 (Ω) or to H r 1 (Ω) and the global regularity estimate ‖ ∇ u ‖ L n n - 1 (Ω ; R n) ≤ C ‖ f ‖ X for inhomogeneous Robin problems with C being a positive constant independent of both u and f. [ABSTRACT FROM AUTHOR]