We focus on the equation of state in neutron stars with the color-spin molecular dynamics. This is based on the constituent quark model, and reproduces the quark confinement/deconfinement phenomena by solving the degrees of freedom of color and spin for each quark. The equation of state for dense matter is studied, and it replicates the saturation properties: symmetric energy, $L$-parameter, and incompressibility around nuclear density. The resultant mass-radius relations are consistent with the observations such as the gravitational wave observations, NICER etc. Our numerical results suggest that de-confined quark matter appears in the core of neutron stars via crossover. Although the current constraints from the observations are not enough to conclude whether quark matter appears at high-density region, our method would help to understand high-density material properties inside neutron stars in the future.