Here, we present an efficient algorithmic technique and implementation to solve this hard problem. Specifically, we present a sequential Monte Carlo method based on a recurrent network search algorithm to simultaneously identify a reaction network, its constituent reactant sets, and the final products that satisfy the arbitrary target properties. An important practical requirement considered here is the maintenance of diversity in the designed molecules and synthetic routes. There are always errors in the properties and synthetic routes predicted by statistical models. Therefore, the optimal solution in a model is not optimal in practice. Our method aims to exhaustively identify a wide variety of promising candidates and to present various scenarios to domain experts. The final decision is left to the expert. A Python library called Seq-Stack-Reaction has been made available on GitHub , which allows us to plug-in any forward model into the Bayesian design workflow (Figure ).