Photopolymerization is widely used for creating materials for biomedical applications, lithography, and fast 3D fabrication. The resin composition and curing protocol during polymerization define a high-dimensional parameter space that dictates the reaction kinetics, network structures, and physical properties of polymerized materials. But a quantitative map from the input parameter space to the transient and final properties does not exist. Here, a computational method is presented to simulate network growth as a stochastic process using reaction-diffusion master equations. The evolving diffusivity of reacting species during polymerization is modeled using effective medium method, thereby making the rates of reaction and fluctuation events sensitive to the spatial configuration of a growing network. The numerical studies show that this method can correlate fabrication parameters to polymer properties like crosslink heterogeneity, relaxation time, phonon density of states, and bulk modulus, in silico. Hence, stochastic simulation of polymerization is broadly applicable to accelerate the design of polymeric materials with targeted physical properties for specific applications.Photopolymerization of networks affords a rapid, solvent-free approach to produce polymers and composites, which is used for fabricating microfluidic sensors, [1,2] biomedical implants, [3,4] and recently for fast 3D printing. [5] Two sets of parameters control the kinetics and properties of photopolymerization: 1) chemical composition of the resin (i.e., rate constants, interaction energies, of monomers and initiators), and 2) curing conditions (e.g., light source, curing time, temperature). The prevailing question for polymer design is: how to select polymerization parameters for the production of materials with a specific set of material properties? A quantitative answer to this question is the key for simulation-based design of polymers. [6] But the changing diffusivity during polymerization and the complex relationship between fabrication parameters and the emergent material properties (e.g., phonon density of states) have been significant impediments so far.Traditional simulations performed for polymerization kinetics use ordinary differential equations, called reaction rate equations (RREs), [6][7][8][9][10] which only determine the degree of double-bond conversion (DC) during the reaction; no information regarding spatial organization of the polymer network is obtained. Decrease in diffusivity is incorporated in RREs by treating the rate constants as empirical functions of DC, [7] which cannot account for differences in network structures developed at the same DC. So, results from RREs are insufficient for computing mechanical and thermal properties of the cured polymer. Previous efforts at stochastic simulation of polymerization omitted either the fluctuation of the network or the diminishing diffusivity. [11][12][13][14][15] Neglecting these details fail to simulate autoacceleration, post-cure polymerization, and the network structu...