Low-density parity-check (LDPC) codes are among the most prominent error-correction schemes. They find application to fortify various modern storage, communication, and computing systems. Protograph-based (PB) LDPC codes offer many degrees of freedom in the code design and enable fast encoding and decoding. In particular, spatially-coupled (SC) and multi-dimensional (MD) circulant-based codes are PB-LDPC codes with excellent performance. Efficient finite-length (FL) algorithms are required in order to effectively exploit the available degrees of freedom offered by SC partitioning, lifting, and MD relocations. In this paper, we propose a novel Markov chain Monte Carlo (MCMC or MC$^2$) method to perform this FL optimization, addressing the removal of short cycles. While iterating, we draw samples from a defined distribution where the probability decreases as the number of short cycles from the previous iteration increases. We analyze our MC$^2$ method theoretically as we prove the invariance of the Markov chain where each state represents a possible partitioning or lifting arrangement. Via our simulations, we then fit the distribution of the number of cycles resulting from a given arrangement on a Gaussian distribution. We derive estimates for cycle counts that are close to the actual counts. Furthermore, we derive the order of the expected number of iterations required by our approach to reach a local minimum as well as the size of the Markov chain recurrent class. Our approach is compatible with code design techniques based on gradient-descent. Numerical results show that our MC$^2$ method generates SC codes with remarkably less number of short cycles compared with the current state-of-the-art. Moreover, to reach the same number of cycles, our method requires orders of magnitude less overall time compared with the available literature methods.