  The equation of State for dense matter is studied with color molecular
dynamics, in which hadron matter and quark matter are automatically
distinguished only from quark color state. The quark-quark interactions are
optimized to be consistent with saturation properties: symmetric energy,
$L-$parameter, and incompressibility around nuclear density. In the
calculations, the degrees of freedom of colors are solved at each numerical
step, although the flavors are fixed as up or down quarks. The resultant
mass-radius relations also satisfy the observational constraints such as the
gravitational wave observations, NICER, and ``{\it the two-solar mass
observations}" of neutron stars. In this model with the allowed parameter
range, deconfined 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.
