  We study the dynamical evolution of globular clusters containing populations
of primordial binaries, using our newly updated Monte Carlo cluster evolution
code with the inclusion of direct integration of binary scattering
interactions. We describe the modifications we have made to the code, as well
as improvements we have made to the core Monte Carlo method. We present several
test calculations to verify the validity of the new code, and perform many
comparisons with previous analytical and numerical work in the literature. We
simulate the evolution of a large grid of models, with a wide range of initial
cluster profiles, and with binary fractions ranging from 0 to 1, and compare
with observations of Galactic globular clusters. We find that our code yields
very good agreement with direct N-body simulations of clusters with primordial
binaries, but yields some results that differ significantly from other
approximate methods. Notably, the direct integration of binary interactions
reduces their energy generation rate relative to the simple recipes used in
Paper III, and yields smaller core radii. Our results for the structural
parameters of clusters during the binary-burning phase are now in the tail of
the range of parameters for observed clusters, implying that either clusters
are born significantly more or less centrally concentrated than has been
previously considered, or that there are additional physical processes beyond
two-body relaxation and binary interactions that affect the structural
characteristics of clusters.
