We consider the problem of optimizing by sampling under multiple black-box constraints in nano-material design. We leverage the posterior regularization framework and show that the constraint satisfaction problem can be formulated as sampling from a Gibbs distribution. The main challenges come from the black-box nature of the constraints obtained by solving complex and expensive PDEs. To circumvent these issues, we introduce Surrogate-based Constrained Langevin dynamics for black-box sampling. We devise two approaches for learning surrogate gradients of the black-box functions: first, by using zero-order gradients approximations; and second, by approximating the Langevin gradients with deep neural networks. We prove the convergence of both approaches when the target distribution is -concave and smooth. We also show the effectiveness of our approaches over Bayesian optimization in designing optimal nano-porous material configurations that achieve low thermal conductivity and reasonable mechanical stability.