Variational inequality (VI) generalizes many mathematical programming problems and has a wide variety of applications. One class of VI solution methods is to reformulate a VI into a normal map nonsmooth equation system, which is then solved using nonsmooth equation-solving techniques. In this article, we propose a first practical approach for furnishing B-subdifferential elements of the normal map, which in turn enables solving the normal map equation system using variants of the B-subdifferential-based nonsmooth Newton method. It is shown that our new method requires less stringent conditions to achieve local convergence than some other established methods, and thus guarantees convergence in certain cases where other methods may fail. We compute a B-subdifferential element using the LD-derivative, which is a recently established generalized derivative concept. In our new approach, an LD-derivative is computed by solving a sequence of strictly convex quadratic programs, which can be terminated early under certain conditions. Numerical examples are provided to illustrate the convergence properties of our new method, based on a proof-of-concept implementation in Python.