<p>Under compressive loading, rocks commonly exhibit tensile–shear mixed-mode fracture behavior. Accurate prediction of crack propagation is essential for controlling rock mass instability and related failures. This paper presents a modified coupled peri-dynamics least squares minimization and finite element method (PDLSM-FEM) for simulating mixed-mode crack propagation in rocks. The maximum principal stress criterion and the Mohr–Coulomb criterion are integrated into the PDLSM-FEM framework to distinguish between tensile and shear failure modes. Meanwhile, potential bond-breaking regions are identified based on the directions of maximum principal stress and maximum shear stress, and a maximum allowable number of bond failures per load step (<i>K</i>_top) are introduced to suppress unstable crack growth, thereby establishing a stable implicit quasi-static solution scheme. The proposed method is applied to several benchmark problems and a series of uniaxial compression tests on fissured rock specimens. Numerical results are compared qualitatively and quantitatively with experimental data. The results show that the method accurately captures the initiation, propagation, and coalescence of cracks in rock, without requiring volume correction or handling boundary effects. It also distinguishes tensile and shear cracks automatically. These findings provide technical support for predicting complex fracture behaviors in rock materials.</p>