We analyze domino tilings of the two-periodic Aztec diamond by means of matrix valued orthogonal polynomials that we obtain from a reformulation of the Aztec diamond as a non-intersecting path model with periodic transition matrices. In a more general framework we express the correlation kernel for the underlying determinantal point process as a double contour integral that contains the reproducing kernel of matrix valued orthogonal polynomials. We use the Riemann-Hilbert problem to simplify this formula for the case of the two-periodic Aztec diamond. In the large size limit we recover the three phases of the model known as solid, liquid and gas. We describe the fine asymptotics for the gas phase and at the cusp points of the liquid-gas boundary, thereby complementing and extending results of Chhita and Johansson.