This study develops a molecular dynamics (MD)-based numerical method implemented in LAMMPS to simulate the experimental stress–strain curves from tensile tests of epoxy thermosets. We focus on crosslinked structures composed of diglycidyl ether of bisphenol A (DGEBA) as the base resin and 4,4'-diaminodiphenylsulfone (4,4’-DDS) as the curing agent. The novelty of our proposed approach is that it incorporates covalent bond rupture to reproduce the failure behavior of crosslinked structures under tensile loading. First, molecular models for DGEBA and 4,4’-DDS are created using the Polymer Consistent Force Field (PCFF) and General AMBER Force Field (GAFF). A model containing approximately 30,000 atoms is constructed, and the ratio of reactive functional groups is adjusted to achieve stoichiometric equivalence. The system is equilibrated by performing a 1 ns relaxation calculation at a temperature of 300 K and a pressure of 1 atm, resulting in a liquid-like equilibrium state. The chemical reaction between the epoxy and amino groups is simulated to form a crosslinked network. After crosslinking, an additional 1 ns relaxation simulation is performed at a temperature of 300 K and a pressure of 1 atm to stabilize the system. Uniaxial tensile simulations are then performed at a strain rate of 5 × 10⁹ /s. Bond rupture is incorporated by defining a threshold at which a covalent bond breaks when its length exceeds a specified multiple of its equilibrium length. By adjusting this threshold, our proposed approach accurately captures Young’s modulus, tensile strength, and fracture behavior, including the characteristic brittle fracture of epoxy resins. It successfully reflects experimental trends in these properties for both the PCFF and GAFF force fields. This framework will be extended to efficiently perform tensile simulations for a wide range of epoxy curing materials in the future.