Summary: A first-order-accurate-in-time, finite difference scheme is proposed and analyzed for the Ericksen-Leslie system, which describes the evolution of nematic liquid crystals. For the penalty function to approximate the constraint \(|\mathbf{d}|=1\), a convex-concave decomposition for the corresponding energy functional is applied. In addition, appropriate semi-implicit treatments are adopted for the convection terms, for both the velocity vector and orientation vector, as well as the coupled elastic stress terms. In turn, all the semi-implicit terms can be represented as a linear operator of a vector potential, and its combination with the convex splitting discretization for the penalty function leads to a unique solvability analysis for the proposed numerical scheme. Furthermore, a careful estimate reveals an unconditional energy stability of the numerical system, composed of the kinematic energy and internal elastic energies. More importantly, we provide an optimal rate convergence analysis and error estimate for the numerical scheme. In addition, a nonlinear iteration solver is outlined, and the numerical accuracy test results are presented, which confirm the optimal rate convergence estimate.