In Section \[sec:MC\] we confirm our analytic predictions via extensive classical MC simulations of the generalized RK wave functions. For the study of the line of critical points at low doping we employ the canonical generalized geometric cluster algorithm (GGCA), whereas the first-order transition is studied via grand-canonical Metropolis-type simulations. The GGCA algorithm enables us to study relatively large systems, up to $400 \times 400$, for a range of dopings, $0 \leq x \leq 0.06$, and to investigate the finite-size scaling behavior. The accessible range of system sizes should be compared to what can be reached for full quantum models, away from the RK condition, where available methods, such as exact diagonalization and Green function Monte Carlo, allow the study of only very small systems with few holes.[@syljuasen05; @Poilblanc06]

In the undoped case we confirm the existence of a Kosterlitz-Thouless transition from a line of critical points to an ordered columnar state, as found in the work of Alet and coworkers.[@alet05; @alet06] We study the scaling behavior of the columnar and orientational order parameters and of their susceptibilities. We also use a mapping of the orientational order parameter of the interacting classical dimer model to the staggered polarization operator obtained 