  We describe a stochastic series expansion (SSE) quantum Monte Carlo method
for a two-dimensional S=1/2 XY-model (or, equivalently, hard-core bosons at
half-filling) which in addition to the standard pair interaction J includes a
four-particle term K that flips spins on a square plaquette. The model has
three ordered ground state phases; for K/J<8 it has long-range xy spin order
(superfluid bosons), for K/J>15 it has staggered spin order in the z direction
(charge-density-wave), and between these phases it is in a state with columnar
order in the bond and plaquette energy densities. We discuss an implementation
of directed-loop updates for the SSE simulations of this model and also
introduce a "multi-branch" cluster update which significantly reduces the
autocorrelation times for large K/J. In addition to the pure J-K model, which
in the z basis has only off-diagonal terms, we also discuss modifications of
the algorithm needed when various diagonal interactions are included.
