Figure 2A shows the volume change trend of the asphalt binder model during a stepwise cooling process. Considering a constant volume expansion coefficient for the model in the glassy and rubbery states (i.e., the volume changes linearly with temperature), the intersection temperature of the two fitted lines from the corresponding low-temperature and high-temperature regions (LL fit and HL fit) represents T g . Figure 2B illustrates the trend of the coefficient of determination R 2 for the linear fits as the selected data points vary. Here, R 2 = 1 − ∑ i = 1 n y i − f x i 2 ∑ i = 1 n y i − y ¯ 2 . y i and y ¯ represent the data to be fitted and its mean value, f x i represents the fitted data. According to Figure 2B, the coefficient of determination R 2 is greater than 0.998 in the temperature ranges of 100–200 K and 340–400 K. However, it exhibits a relatively large decrease in value for temperatures above 180 K or below 340 K. Therefore, the linear fitting range should be based on the aforementioned temperature ranges. Based on the fitting results from 100–160 K and 340–400 K, the intersection temperature of the asphalt binder model is approximately 271.10 K, consistent with the experimental data from the Strategic Highway Research Program (SHRP) at around 273.15 K (Khabaz and Khare, 2015).

From Figure 2A, it is evident that the value of T g depends on the selected data points. In light of this, we further employ a data-driven fitting (DDF) analysis method. Based on a 10-fold cross-validation, the computed data is divided into ten parts, and the volume-temperature relationship is fitted using piecewise cubic Hermite interpolation polynomials. Since the data points in Figure 2A represent the averaged values of the model’s relaxation under the NPT ensemble, this can lead to local perturbations and overfitting when the number of data points is limited. To mitigate this overfitting issue, stability during relaxation can be addressed by randomly selecting 2-3 temperature-volume data points for DDF fitting. According to Figure 2A, the fitting model established by the DDF method aligns well with the original data. By differentiating the DDF fitting model, the rate of change of the model’s volume with temperature is obtained according to α = dV / dT . From Figure 2C, it is evident that the values in the temperature ranges of 100–160 K and 360–400 K tend towards constants, indicating a linear change in volume with temperature, consistent with the coefficient of determination for linear fitting shown in Figure 2B. These results demonstrate that the DDF model can effectively analyze and determine the temperature range for linear fitting.

Due to the different motion characteristics of molecules during the cooling and heating processes, molecular systems typically exhibit hysteresis effects in their thermal and volumetric properties during the glassy state transition (Minisini and Soldera, 2023). Hence, it is necessary to explore the influence of the simulation process on T g . For such purpose, we reheat the model to 400 K under the NPT ensemble and then re-initiated the stepwise cooling simulation. Figure 2D compares the temperature-volume curves obtained from the three stepwise cooling simulations. It is evident from the figure that the third reheating process results in a noticeable upward shift in the system’s volume-temperature curve, while the volume-temperature curves from the first and second cycles are closely aligned. These results indicate that irreversible microstructural changes occur in the model during repeated heating and cooling processes, leading to an overall increase in the model’s volume. Through derivative analysis of the DDF fitting model, the value of α in the temperature ranges of 100–160 K and 360–400 K remain constant, indicating a linear change in volume with temperature. Based on linear fitting, the value of T g exhibits a 6% fluctuation. Overall, the simulation process has an insignificant impact on the glassy state transition temperature.

The glass transition region of the system can be identified from the volume-temperature curve. Based on the results from Figure 2A, the deviation from linear behavior in the volume-temperature relationship occurs between 180 K and 340 K, corresponding to the beginning and end of the glassy state transition, respectively. This can be expressed as T g L and T g U , and the glass transition temperature threshold is defined as ∆ T g = T g U − T g L . In this study, the behavior of the system’s heat capacity with respect to temperature was further analyzed. Experimentally, it is common to define the glass transition region by scanning heat capacity. A change in the heat capacity-temperature curve indicates the beginning or end of the glassy state transition (Apostolidis et al., 2021; Soenen et al., 2013; Nikookar et al., 2022; Tayfun et al., 2015). To investigate whether this method is applicable to molecular dynamics simulation models, we conducted data sampling based on the aforementioned DDF model and the heat capacity is calculated according to C V = dU / dT , where U represents the internal potential energy and T represents the temperature. The calculation yields the results shown in Figure 3, it can be observed that there are three temperature domains during the expansion process of the simulated system, which are consistent with experimental results. However, the glass transition region is wider in the simulation, measuring 180 K, compared to experimental findings. This discrepancy arises due to the cooling rate in the simulation being approximately 10¹² times faster than the experimental rate. The glass transition is a gradual process, and during the rapid cooling process in the simulation, molecules may have frozen, but this may not be immediately reflected in the system’s heat capacity. Therefore, volume-temperature curve has been commonly employed to calculate the glassy state transition region in MD simulations and calculation of the glassy state transition region based on the volume-temperature curve. In conclusion, MD simulations can effectively represent the microscopic-scale glassy state transition phenomenon, albeit occurring over a broader temperature range.

The above defines the transition range of the glassy state from a molecular thermodynamic perspective (volume and energy), and similarly, in molecular dynamics, it is possible to capture the differences between the asphalt glassy state and rubbery state. Figure 4 depicts the mean square displacement (MSD) plot from MD simulations under the NPT ensemble for asphalt components in the rubbery state (400 K) and glassy state (100 K) with a 0.5 ns timescale. It is evident that the motion of asphalt at the two temperatures is entirely different. At 400 K, the components of asphalt are nearly in a flowing state, especially the saturated components, which undergo extensive motion within the system. Conversely, at 100 K, the asphalt components are only engaged in localized vibrations. This demonstrates the transition of asphalt components from a free-flowing state to a “frozen state.”

Since the glass transition phenomenon is directly related to molecular motion, and molecular motion is temperature-dependent, it is necessary to examine the effect of cooling rates on the glassy state transition. For macroscopic-scale experiments, the formation of the glassy state requires a relatively rapid cooling rate. Asphalt molecules, unable to rearrange slowly during cooling, consequently form the glassy state. Typically, the standard cooling rate for experiments is 20 K/min. In contrast, for different simulation systems, MD simulations currently employ cooling rates of 20 K or 25 K per cycle. Due to the disparity in time scales, the cooling rate in MD simulations differs by approximately 10¹² orders of magnitude from experimental cooling rates. Considering that the simulated cooling time for each cycle in this study is 100 ps, the corresponding cooling rate is 0.2 K/ps. To further examine the potential impact of cooling rates, five different cooling rates were selected, including 0.05, 0.10, 0.15, 0.20, and 0.25 K/ps.

The volume-temperature relationships of asphalt at different cooling rates are compared in Figure 5A, and the nonlinear trends indicate the presence of the glass transition phenomenon. As expected, the cooling rate has a more pronounced effect on the volume at lower temperatures. Specifically, at cooling rates of 0.20 and 0.25 K/ps, the volumes almost overlap within the range of 100–400 K. However, when cooling rate below 0.15 K/ps, a significant decrease in volume is observed in the lower temperature region (under 250 K). This result suggests that when cooling rate below 0.15 K/ps, asphalt molecules undergo more thorough microstructural adjustments, leading to an overall decrease in volume. To compare the effects of different cooling rates, we statistically analyzed the fitted data within the low-temperature ranges of 100–140 K and 100–180 K (corresponding to the high-temperature ranges of 360–400 K and 320–400 K). As shown in Figure 5B, at around 270 K, the fittings under cooling rates of 0.20 and 0.25 K/ps maintain good consistency, while due to more thorough microstructural relaxation at lower cooling rates, the simulated system exhibits a slightly lower value. Based on five simulation calculations, the average T g is 266.33 K with a standard deviation of 6.65 K.

To investigate the influence of cooling rate on the structure of asphalt, we selected the same model and examined the aggregation of asphalt at 100 K under different cooling rates (0.05 and 0.25 K/ps). We define free volume as the sum of the regions within which each atom can move freely. According to the coordination number, we define the polyhedron for each atom, and the volume of this polyhedron represents the free movement region of the atom. Figure 6A shows the trend of free volume variation with temperature at different cooling rates, when the cooling rate was 0.20 K/ps, the decrease in free volume is slightly larger, confirming that slow cooling leads to more complete molecular rearrangement, resulting in more uniform system molecules. Figure 6B displays the x-y direction view of the model, from which it is observed that under slow cooling conditions, the asphalt has sufficient time to rearrange, forming larger aggregations and local voids, resulting in a relatively dense system. On the other hand, it illustrates that under rapid cooling conditions, the asphalt does not have enough time to rearrange, leading to a more dispersed system and consequently a relatively larger volume at low temperatures. This observation aligns with the previous volume calculation results.

The transition of an amorphous solid (glass transition) is directly related to the motion of different molecules, and the model size affects the number of molecules, their aggregation, and kinetic behavior. To further investigate T g , we examined the cubic models with six different sizes, including 50 Å, 60 Å, 80 Å, 100 Å, 110 Å, and 130 Å, considering two cooling rates of 0.20 and 0.25 K/ps. In the process of modeling, we first set different box dimensions, then set the same density, to obtain models of different sizes. The proportions of different molecules within each model were kept consistent. To reduce calculation errors, we performed simulations for three initial models for each size with a relaxation time of 2 ns, 2.05 ns and 2.10 ns, respectively. The glass transition temperature was averaged from these three calculations. According to Figure 7, the glass transition temperature converges in general when the size increases. Convergence begins at a size of 100 Å, with a convergent value of approximately 266 K, close to the experimental value (Yang et al., 2020). This indicates that larger models can better reflect the volume variation with temperature, suggesting the use of models of 100 Å or larger for studying the performance of asphalt (Li et al., 2022).

Based on the simulation results, we further investigated the influence of model size on the glass transition region. Figure 8A illustrates the volume-temperature curves for models of 50 Å under two cooling rates. Due to the relatively small model size, the sensitivity of volume to temperature increases (microscopic pressure and temperature control exhibit relatively large fluctuations), resulting in a relatively large standard deviation for each volume data point. As the model size increases to 50 Å, the fluctuation in volume during the relaxation phase notably decreases, as shown in Figure 8B. Apart from volume fluctuations, we observed a good consistency in the volume of the simulated systems under both cooling rates. The graph clearly demonstrates the close relationship between the glass transition region and the model size. Specifically, when the model size is 50 Å, the glass transition region is small, but it significantly widens when the size increases to 100 Å. This phenomenon is possibly due to the larger molecular weight in larger models, leading to more complex entanglement between molecules and longer relaxation times, thus resulting in a wider glass transition region.