A variational formulation of the non-smooth contact dynamics method is proposed to address the dynamic response of historical masonry structures modeled as systems of 3D rigid blocks and subjected to ground excitation. Upon assuming a unilateral-frictional contact law between the blocks, the equations of motions are formulated in a time-discrete impulse theorem format in the unknown block velocities and contact impulses. The variational structure of the problem to be solved at each time step is proven. On that basis, the numerical method requires at each time step to perform a collision detection that identifies antagonist contact points based on the given structural configuration, to solve a second-order conic programming problem that outputs block velocities and contact impulses, and to update the structural configuration for the solution to advance in time. As a merit of the formulation, large-scale problems can be robustly and efficiently addressed thanks to the convex setting of the time-step optimization problem. Numerical results are presented to test the computational performances of the proposed approach. Benchmark problems provide numerical evidence that the formulation is consistent with event-driven solutions based on the classical Housner impact model. The dynamic response, failure domains, and fragility functions of real-size masonry structures are then explored under ground impulse or earthquake excitation. The obtained results prove the reliability of the present computational method for the dynamic analysis and seismic assessment of historical masonry constructions of engineering interest.