A 3D limit analysis model specialized in the analysis at the collapse of massive masonry structures subjected to seismic loads is proposed. The approach belongs to the wide family of distinct element methods, because is based on a discretization of the structure by means of infinitely resistant hexahedron elements and quadrilateral interfaces where all dissipation occurs. Any failure surface -after proper linearization- can be adopted for interfaces, as for instance a Mohr-Coulomb failure criterion with tension and compression cutoffs. Having discretized the continuum in such a way that internal dissipation is allowed only between adjoining elements, the corresponding linear programming problem necessary to estimate the collapse load, can be obtained through both a kinematic and static approach. The first is more straightforward and the primal problem provides a failure mechanism, interfaces plastic multipliers and collapse multipliers. By virtue of the duality theorem, the second lower bound approach is associated with the same failure load, also providing the internal forces acting on interfaces. The model is validated on three examples of historical value, namely-three pagodas in China, all characterized by complex 3D geometries and massive walls. The results show how the procedure proposed is able to provide fast estimations of the load carrying capacity under the application of horizontal seismic actions and how the corresponding failure mechanism is rather complex, thus precluding the utilization of simplified approaches on idealized geometries.