We calculate the ghost two-point function in Coulomb gauge QCD with a simple
model vacuum gluon wavefunction using Monte Carlo integration. This approach
extends the previous analytic studies of the ghost propagator with this ansatz,
where a ladder-rainbow expansion was unavoidable for calculating the path
integral over gluon field configurations. The new approach allows us to study
the possible critical behavior of the coupling constant, as well as the Coulomb
potential derived from the ghost dressing function. We demonstrate that IR
enhancement of the ghost correlator or Coulomb form factor fails to
quantitatively reproduce confinement using Gaussian vacuum wavefunctional