I described the method we use and some of the computational difficulties in some previous threads
https://news.ycombinator.com/item?id=15780514
https://news.ycombinator.com/item?id=12048170
But, if you're comfortable with an argument from authority, suffice it to say that much of the computational power in the largest supercomputing centers all across the world are dedicated to quantum many-body problems.
If you want equilibrium properties at temperature T [non equilibrium properties, like real-time dynamics have additional exponentially-bad computational intractabilities, known as the sign problem], you want to evaluate a partition function
Z = tr[ exp( - H / T ) ]
where H is the Hamiltonian that, given an eigenstate, has an the state's energy as the eigenvalue. The space H acts on grows exponentially with the number of particles in your simulation [simulations to determine crystal structures are typically in the canonical ensemble; if that doesn't mean anything to you, that's OK, but it's something you can read further on; basically it means: you try to take the number of particles to be large, rather than, for example, fix a density].We have importance-sampling Markov Chain Monte Carlo techniques for evaluating Z. For some problems it succeeds with flying colors, although demonstrating that you're in the thermodynamic [infinite number of particles and infinite volume] limit requires calculating at a variety of volumes and particle numbers, for extrapolation. If you just do random updates to your state vector, most updates will be rejected; use HMC to update the whole state at once [HMC = Hybrid Monte Carlo, or Hamiltonian Monte Carlo, if you're not a computational physicist, strangely].
This approach is not naively parallelizable, and the machines where HMC is performed requires extremely high-performance interconnects, think 10-100 times faster than what's available in a standard datacenter.
The method is inherently statistical: give me more computing time and I can shrink the error bars. But if the problem has certain technical difficulties (or is in real-time, for example, rather than equilibrium), you encounter the sign problem, where the partition function you want to evaluate now becomes
Z = tr[ exp( i H / T ) ]
and you need an exponentially large sample to resolve the intricate cancellations that arise from summing just a bunch of phases.