This paper describes a method of simulating the operational shock (op-shock) response of hard disk drives (HDDs) and its application to improving the op-shock robustness of HDDs. This flexible multibody dynamics simulation model is based on component-mode synthesis of three components; a rotating part including rotating disks and the hub/shaft of the spindle motor, a stationary part including base plate, top cover, and other major stationary parts, and an actuator part representing the moving parts of the head actuator mounted on head sliders. These components are connected to each other by fluid dynamic bearings, pivot bearings, and air bearings, and their nonlinear characteristics are considered in the op-shock response simulation. Linear and nonlinear drive-level simulation models were built, and their accuracies were experimentally verified in terms of their predicted linear frequency response, nonlinear time historical response, and head lift-off boundary. Moreover, a parametric study was performed to improve the op-shock robustness of a 2.5-inch HDD for the head lift-off boundaries of 3920 m/s2 (400 G) and 4900 m/s2 (500 G). The study indicated that the HDD model sensitivities to changes in the stiffness of the base plate and in the thickness of the disk were substantial but they saturated as the parameters became larger. The results indicated that the head lift-off boundary of 3920 m/s2 can be reached by making small modifications to the parts design, but the 4900 m/s2 boundary can’t be reached without the whole drive system being redesigned.