We establish a semi-classical formula for the sum of eigenvalues of a magnetic Schrodinger operator in a three-dimensional domain with compact smooth boundary and Neumann boundary conditions. The eigenvalues we consider have eigenfunctions localized near the boundary of the domain, hence they correspond to surface states. Using relevant coordinates that straighten out the boundary, the leading order term of the energy is described in terms of the eigenvalues of model operators in the half-axis and the half-plane.