<p>Optimal exploitation of supercomputing resources for the evaluation of electrostatic forces remains a challenge in molecular dynamics simulations of very large systems. The most efficient methods are currently based on particle-mesh Ewald sums and achieve semi-logarithmic scaling in the number of particles. These methods solve the problem in reciprocal space, requiring extensive use of Fast Fourier transforms (FFTs). While highly efficient in many contexts, FFTs may encounter scalability challenges at very large processor counts due to their communication requirements. To mitigate these problems, the development and scalable coding of real-space approaches to solve the Poisson equation on a grid is an active field of research. In this work, we introduce a novel real-space approach that provides some advantages over alternatives. Our method exploits an extended Lagrangian in which the values of the field at the grid points are treated as auxiliary variables of zero inertia and the discretized Poisson equation is enforced as a dynamical constraint. The solution of the constraints leads to a linear system - different from those appearing in other real-space approaches - that can be efficiently solved via state-of-the-art iterative methods. The method inherits the numerical scaling of the adopted iterative solver, e.g. linear with a multigrid (MG) approach, but converges with fewer cycles. We analyze this approach considering realistic simulations of molten NaCl that validate its ability to reproduce structural and transport properties. Using this non-trivial benchmark, we demonstrate linear scaling and illustrate some features of our algorithm.</p>