Conventional deterministic design sensitivity analysis (DSA) is carried out by taking the variation of molecular dynamics (MD) trajectory with respect to design variables. The deterministic DSA requires an additional time integration for design sensitivity computation, which usually results in numerically unstability or divergence as the time integration marches in chaotic systems. Thus, we propose a stochastic design sensitivity formulation for MD based on the linear response theory, which observes the probability variations with small changes of Hamiltonian systems near thermodynamic equilibrium. In addition, considering the common workflow where an equilibration process is performed in an ensemble different from the production ensemble, we also present an efficiency strategy that consistently accounts for the equilibration effects on the production ensemble, and thereby further reduce the computational costs of sensitivity evaluation. The computed design sensitivity provides reasonable values compared with averaged finite difference results even though long-time simulations. We obtained the optimal design for the minimal thermal conductivity of graphene by locating isotopes and antidots in the graphene using a gradient-based design optimization method, which requires a stochastic design sensitivity. Thermal conductivity is determined via the Green-Kubo expression. The obtained optimal structure is an aperiodic isotope superlattice structure with the antidots located at interfaces. In general, the acoustic phonons and the low-frequency optical phonons in dispersion curves have a significant impact on the thermal conductivity. Also, the slope of the phonon curves account for the group velocity. We notice that the slopes of the acoustic phonons and the low-frequency optical phonons get decreased as the design optimization progresses. Also, investigating the effects of the stochastic design optimization on the density of states (DOS), the DOS of the optimal model are smaller than that of the pristine one in low-frequency ranges. REFERENCES [1] J. Kim, S. Cho, “Determination of optimal potential parameters for the self-assembly of various lattice structures,” Nanocomposites, 9(1), 18-29, 2022. [2] S. Cho, M. Oh, W. Park, H. Kim, “Stochastic Design Sensitivity and Optimization for Maximum Bulk Modulus using Molecular Dynamics Simulation,” Materials Today Conference 2025, Spain, 2025. 6. 23 ~ 26.