Does it matter which ensemble you use for a molecular dynamics simulation? Shouldn't they in principle give the same results? Why, historically, was it hard (and desirable) to develop simulations in ensembles other than the microcanonical (NVE) ensemble?
It's very hard to answer this definitively unless you are a bit more specific.
However, to speak somewhat generally, it is the case that NVE and NVT ensembles become equivalent in certain thermodynamic limits (namely an infinite number of particles). In practice, however, one chooses the ensemble based on the free energy you are interested in sampling, or the experiment you are interested in comparing to.
If you are modelling a gas-phase reaction, you probably want the NVE ensemble. Unless there is a buffer gas, then you probably want the NVT. If the process takes place in the liquid phase, you probably want NPT. In my experience, though, it is very rare that NVT and NPT simulations are not both carried out if you know you care about the NPT simulation. They don't tend to give all that different of results for the most part.
It also depends a lot on what you are calculating. For instance, if you want to simulate the infrared spectrum of a liquid, then you will probably equilibrate the system in the NVT ensemble at your desired temperature and then carry out an NVE simulation beginning from that equilibrated state. This is because thermostats are designed to de-correlate velocities which will destroy your spectrum calculated from correlation functions.
If you are explicitly calculating free energies, then you just choose your ensemble for the free energy you want to compare to. NVE gives you back the internal energy. NVT gives you back the helmholtz free energy. NPT gives you back the Gibbs free energy.
Why, historically, was it hard (and desirable) to develop simulations in ensembles other than the microcanonical (NVE) ensemble
You would have to provide some reference for this, because it seems somewhat doubtful to me that it was ever difficult to develop NVT simulations. As far as I can tell, some of the first MD simulations knew they needed to do something with the velocities in order to keep the temperature constant throughout a simulation. As time went one, more sophisticated thermostats were developed, and it was shown more rigorously that give the right fluctuations to actually sample the canonical ensemble. I don't think this idea was ever particularly difficult to understand though.
As to why it was desriable to sample other ensembles, I think that is because it was obvious to people that their MD simulations were not at the appropriate thermodynamic limits as to make ensembles equivalent. Even in modern simulations with many particles and periodic boundary conditions, one will usually get different results depending on the ensemble you choose to use. So, it is desirable to develop different ensembles because MD is not easily capable of reaching the same limits which exist in nature.
It definitely can matter. The best thing to do is to consider the experiment you are comparing to, even if that experiment hasn't actually been done yet. In many cases, that means you want to run dynamics in the NPT ensemble (since bench experiments tend to be at fixed pressure and temperature).
Here's an example where we know we won't get the same results between NVE and NVT: Consider the calculation of the rate (we can do this in 1D), but assume that the barrier height is just below the total energy for the NVE ensemble. The rate is 0: it can never cross the barrier.
However, in NVT with an average energy equal to the energy of your NVE ensemble, the rate will not be zero! Sometimes the thermal energy is enough to push you over the barrier.
I think I've already covered why it is desirable (we want to compare with experiment!) I'm not sure I'd describe as being all that hard to do, but:
Thermostats and barostats came after NVE because we've had Newton's equations since long before modern statistical mechanics. Newton's ideas were necessary to inspire people like Boltzmann, not to mention the early computational chemists.
Dynamics with a thermostat or barostat aren't actually "real." They're modeling what should be many different trajectories (or a much larger system) by using only one trajectory. This is how you get problems like the lack of ergodicity in Nosé-Hoover on the harmonic oscillator.
Should you worry that you're using an NVT integrator to model a system that is really NVE? Typically only if your system is small/gas phase. Larger condensed phase system have enough chaos that it doesn't matter.
Ensembles are essentially artificial constructs. In the thermodynamic limit (for an infinite system size) and as long as we avoid the neighborhood of phase transitions it is generally believed that there is an equivalence between ensembles.
A consequence of the equivalence of ensembles is that the basic thermodynamic properties of a model system may be calculated as averages in any convenient ensemble.
For more details, "Computer Simulation of Liquids" by Allen and Tildesley.