  Boltzmann distributions over lattices are pervasive in Computational Physics.
Sampling them becomes increasingly difficult with increasing lattice-size,
especially near critical regions, e.g., phase transitions in statistical
systems and continuum limits in lattice field theory. Machine learning-based
methods, such as normalizing flows, have demonstrated the ability to alleviate
these issues for small lattices. However, scaling to large lattices is a major
challenge. We present a novel approach called Parallelizable Block
Metropolis-within-Gibbs (PBMG) for generating samples for any lattice model. It
factorizes the joint distribution of the lattice into local parametric kernels,
thereby allowing efficient sampling of very large lattices. We validate our
approach on the XY model and the Scalar {\phi}^4 theory. PBMG achieves high
acceptance rates and less correlated samples, and the observable statistics
estimated from the samples match the ground truth. Moreover, PBMG significantly
speeds up inference for large lattices as compared to HMC and plain MCMC
algorithms.
