Abstract: We present the current status of molecular dynamics, in which hadron and quark matter are described only within quark models. In our previous calculations, we solved the time evolution of internal degrees of freedom for color. In this article, we show the results simultaneously solved the internal degrees of freedom for spin as well. Unlike conventional methods, our approach also considers the flavor conversion between u,d-quarks and electons during the calculations, which naturally achieves the beta equilibrium under the charge neutrality condition. Based on this latest formulation, the equation of state for dense matter is studied, replicating the saturation properties: symmetric energy, L−parameter, incompressibility around nuclear density. The resultant mass-radius relations are also consistent with the observations such as the gravitational wave observations, and NICER. Our numerical results suggest that deconfined quark matter appears in the core of neutron stars via crossover. The key difference from previous molecular dynamics simulations is the spontaneous emergence of clustering of quarks or baryons on a local scale. These properties would never be observed in equations of state based on the other approaches. While our formulation does not yet fully include strangeness, we expect that our method has a possibility to show the direct evidence of the quark clustering models such as “the strangeon model” or “the sexaquark model”.