  We present BIFROST, an extended version of the GPU-accelerated hierarchical
fourth-order forward symplectic integrator code FROST. BIFROST (BInaries in
FROST) can efficiently evolve collisional stellar systems with arbitrary binary
fractions up to $f_\mathrm{bin}=100\%$ by using secular and regularised
integration for binaries, triples, multiple systems or small clusters around
black holes within the fourth-order forward integrator framework.
Post-Newtonian (PN) terms up to order PN3.5 are included in the equations of
motion of compact subsystems with optional three-body and spin-dependent terms.
PN1.0 terms for interactions with black holes are computed everywhere in the
simulation domain. The code has several merger criteria (gravitational-wave
inspirals, tidal disruption events and stellar and compact object collisions)
with the addition of relativistic recoil kicks for compact object mergers. We
show that for systems with $N$ particles the scaling of the code remains good
up to $N_\mathrm{GPU} \sim 40\times N / 10^6$ GPUs and that the increasing
binary fractions up to 100 per cent hardly increase the code running time (less
than a factor $\sim 1.5$). We also validate the numerical accuracy of BIFROST
by presenting a number of star clusters simulations the most extreme ones
including a core collapse and a merger of two intermediate mass black holes
with a relativistic recoil kick.
