Generalized Shadow Hybrid Monte Carlo (GSHMC) is a method for molecular simulations that rigorously alternates Monte Carlo sampling from a canonical ensemble with integration of trajectories using Molecular Dynamics (MD). While conventional hybrid Monte Carlo methods completely re-sample particle’s velocities between MD trajectories, our method suggests a partial velocity update procedure which keeps a part of the dynamic information throughout the simulation. We use shadow (modified) Hamiltonians, the asymptotic expansions in powers of the discretization parameter corresponding to timestep, which are conserved by symplectic integrators to higher accuracy than true Hamiltonians. We present the implementation of this method into the highly efficient MD code GROMACS and demonstrate its performance and accuracy on computationally expensive systems like proteins in comparison with the molecular dynamics techniques already available in GROMACS. We take advantage of the state-of-the-art algorithms adopted in the code, leading to an optimal implementation of the method. Our implementation introduces virtually no overhead and can accurately recreate complex biological processes, including rare event dynamics, saving much computational time compared with the conventional simulation methods.