  Most recently, the path integral molecular dynamics has been successfully
used to consider the thermodynamics of single-component identical bosons and
fermions. In this work, the path integral molecular dynamics is developed to
simulate the thermodynamics, Green's function and momentum distribution of
two-component bosons in three dimensions. As an example of our general method,
we consider the thermodynamics of up to sixteen bosons in a three-dimensional
harmonic trap. For noninteracting spinor bosons, our simulation shows a bump in
the heat capacity. As the repulsive interaction strength increases, however, we
find the gradual disappearance of the bump in the heat capacity. We believe
this simulation result can be tested by ultracold spinor bosons with optical
lattices and magnetic-field Feshbach resonance to tune the inter-particle
interaction. We also calculate Green's function and momentum distribution of
spinor bosons. Our work facilitates the exact numerical simulation of spinor
bosons, whose property is one of the major problems in ultracold Bose gases.
