As molecular dynamics simulations continue to provide important insights into biomolecular structure and function, a growing demand for increasing the time span of the simulations is emerging. Our focus here is developing a new algorithm, LIN (Langevin/implicit-Euler/normal mode), that combines normal-mode and implicit-integration techniques, for large time step biomolecular applications. In the normal-mode phase of LIN, we solve an approximate linearized Langevin formulation to resolve the rapidly varying components of the motion. In the implicit phase, we resolve the remaining components of the motion by numerical integration with the implicit-Euler scheme. Developments of the normal-mode phase of LIN are discussed in this paper. Specifically, we solve two crucial issues of the method. The first involves how to choose and how often to update the Hessian approximation for the linearized Langevin equation. This approximation must be computationally feasible and physically reasonable to capture the motion in the higher end of the vibrational spectrum. Three such general Hessian approximations are discussed. The related issue-the frequency of the Hessian update-is analyzed by projecting the motion onto the different vibrational modes. This analysis demonstrates that a one-picosecond interval is reasonable for updating the Hessian in the model system examined here. In this connection, we illustrate that the high-frequency motions are highly localized white the low-frequency motions are delocalized. We also show rigorously that the mode amplitudes are inversely proportional to the frequency (consistent with the equipartition theorem), with 90% of the displacement fluctuations coming from a very small group of low-frequency modes. Anharmonic effects essentially influence the low-frequency modes. The second issue involves how to solve tl;e linearized Langevin equation at large timesteps correctly, where the usual discretized formulation of the random force is invalid. This is accomplished by using analytic expressions for the distributions associated with positions and velocities of the individual oscillators as a function of frequency, obtained as the solution of the corresponding Fokker-Planck equation. We apply LIN with these developments to the nucleic acid component deoxycytidine with timesteps ranging from 100 to 1000 fs. We demonstrate that LIN is stable in these simulations, with energies fluctuating about the same values-and possessing overall similar dynamical features-in comparison to 1 fs explicit simulations, though the fluctuations are significantly larger at larger timesteps. Moreover, continuous dynamics is maintained, and pathway information can be obtained. Computational performance is competitive only at very large time steps : a gain factor of 3-4 is obtained for runs with 1000 fs time steps. Larger gains may be achieved for biomolecules, where sparsity and parallelization can be exploited significantly.