Many integral equations for numerically predicting pair correlation functions g(2) in classical fluids have been proposed, each based on its own ``closure approximation.'' We have investigated a procedure for evaluating such closures, particularly their capacity to describe clustering under the influence of attractive interparticle forces. Our approach utilizes the Gaussian core model in a large closed system which, under sign change of the coupling constant λ, undergoes a collapse to form a single compact aggregate with known properties. In the infinite system limit this phenomenon causes inverse temperature (β) expansions, such as that for g(2), to diverge. By contrast available closures lead to convergent expansions for the Gaussian core model, but with a ``critical point'' on the negative real βλ axis. Specific calculations have been performed illustrating this behavior for the BGYK and the PY integral equations. We suggest that the proximity of the artifactual negative-axis critical point to the origin (where the collapse singularity should actually appear) provides a measure of accuracy of closure approximations.