I investigated exactly this question 20 years ago when designing a desktop GIS. We needed to find point-to-point distances interactively; our target was to do the computations in less than 1/2 second for thousands of points. Testing (on a 25 MHz 486 PC!) showed that we could compute all the distances, exactly as you describe (with the simple obvious algorithm), so quickly that it made no sense to create a more sophisticated solution, such as a quadtree structure.
For computing distances to a single "probe" point your options include (a) projecting all points using an equidistant projection centered at the probe point or (b) adopting a spherical earth model and using the Haversine formula. The first is appropriate if you need the accuracy of an ellipsoidal model. In either case the calculations are reasonably fast, probably taking less than 1000 ticks: you could query around a million points a second with a single processor.
Fast enough for you? If not, the brute-force method parallelizes easily and scales directly with the number of processors: just divide the points among the processors and then do a final comparison of the closest one found by each processor.
If you need to go faster, you can use various approximations to screen points. For example, if you are between -88 and +88 degree latitude and the nearest point found so far is 200 km away, then any point whose latitude differs from the probe point's latitude by more than 2 degrees cannot possibly be closer (because anywhere on earth, one degree of latitude exceeds about 110 km). In many cases this kind of pre-screening might enable you to process hundreds of millions of points a second.