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-r
ainbow 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.
