Loss function for multivariate regression where relationship between outputs matters
MSE uses the euclidean distance function. You could simply use the Pearson correlation coefficient as similarity measure, so define a new distance function that decreases when the coefficient increases ($exp(-\rho)$ for instance).

Another possibility is to use the Wasserstein distance (or Wasserstein metric, or earth mover's distance). See the definitions on Wikipedia, especially this one which is not too complicated.
This metric is designed to compute a distance between distributions (continuous) or histograms (binned / discrete). It does not relate to correlation, but rather to the shape of the vector when plotted as an histogram. This is not what you asked for but it might be what you're actually looking for.