In this work, we construct a novel numerical method for solving the multi-marginal optimal transport problems with Coulomb cost. This type of optimal transport problems arises in quantum physics and plays an important role in understanding the strongly correlated quantum systems. With a Monge-like ansatz, the orginal high-dimensional problems are transferred into mathematical programmings with generalized complementarity constraints, and thus the curse of dimensionality is surmounted. However, the latter ones are themselves hard to deal with from both theoretical and practical perspective. Moreover in the presence of nonconvexity, brute-force searching for global solutions becomes prohibitive as the problem size grows large. To this end, we propose a global optimization approach for solving the nonconvex optimization problems, by exploiting an efficient proximal block coordinate descent local solver and an initialization subroutine based on hierarchical grid refinements. We provide numerical simulations on some typical physical systems to show the efficiency of our approach. The results match well with both theoretical predictions and physical intuitions, and give the first visualization of optimal transport maps for some two dimensional systems.