We propose a multivariate k-centers functional clustering algorithm for the multivariate functional data. We assume that clusters can be defined via functional principal components subspace projection for each variable. A newly observed subject with multivariate functions is classified into a best-predicted cluster by minimizing a weighted distance measure, which is a weighted sum of discrepancies in observed functions and their corresponding projections onto the subspaces for all variables, among all the clusters. The weight of the proposed algorithm is flexible and can be chosen by the objective of clustering. The proposed method can take the means and modes of variation differentials among groups of each variable into account simultaneously. Numerical performance of the proposed method is demonstrated by simulation studies, with an application to a data example.