Density functional theory (DFT) is not only an accurate but also a widely used theory for describing the quantum-mechanical electronic structure of matter. In this approach, the intractable problem of interacting electrons is simplified to a tractable problem of non-interacting electrons moving in an effective potential. Even with this simplification, DFT remains extremely computationally expensive. In particular, DFT scales cubically with respect to the number of atoms, which restricts the size of systems that can be studied. Orbital free density functional theory (OF-DFT) represents a simplification of DFT applicable to metallic systems that behave like a free-electron gas. Current implementations of OF-DFT employ the plane-wave basis, the global nature of the basis prevents the efficient use of modern high-performance computer archi- tectures. We present a real-space formulation and higher-order finite-difference implementation of periodic Orbital-free Density Functional Theory (OF-DFT). Specifically, utilizing a local reformulation of the electrostatic and kernel terms, we develop a gener- alized framework suitable for performing OF-DFT simulations with different variants of the electronic kinetic energy. In particular, we develop a self-consistent field (SCF) type fixed-point method for calculations involving linear-response kinetic energy functionals. In doing so, we make the calculation of the electronic ground-state and forces on the nuclei amenable to computations that altogether scale linearly with the number of atoms. We develop a parallel implementation of our method using Portable, Extensible Toolkit for scientific computations (PETSc) suite of data structures and routines. The communication between processors is handled via the Message Passing Interface(MPI). We implement this formulation using the finite-difference discretization, us- ing which we demonstrate that higher-order finite-differences can achieve relatively large convergence rates with respect to mesh-size in both the energies and forces. Additionally, we establish that the fixed-point iteration converges rapidly, and that it can be further accelerated using extrapolation techniques like Anderson mixing. We verify the accuracy of our results by comparing the energies and forces with plane-wave methods for selected examples, one of which is the vacancy formation energy in Aluminum. Overall, we demonstrate that the proposed formulation and implementation is an attractive choice for performing OF-DFT calculations.