  We present a semi-analytic binary black hole (BBH) metric approximation that
models the entire evolution of the system from inspiral to merger. The metric
is constructed as a boosted Kerr-Schild superposition following post-Newtonian
(PN) trajectories at the fourth PN order in the inspiral phase. During merger,
we interpolate the binary metric in time to a single black hole remnant with
properties obtained from numerical relativity (NR) fittings. Different from
previous approaches, the new metric can model binary black holes with arbitrary
spin direction, mass ratio, and eccentricity at any stage of their evolution in
a fast and computationally efficient way. We analyze the properties of our new
metric and we compare it with a full numerical relativity evolution. We show
that Hamiltonian constraints are well-behaved even at merger and that the mass
and spin of the black holes deviate in average only $\sim 5 \%$ compared to the
full numerical evolution. We perform a General Relativistic
magneto-hydrodynamical (GRMHD) simulation of uniform gas evolving on top of our
approximate metric. We compare it with a full numerical relativity evolution of
the fluid and Einstein's equations. We show that the properties of the gas such
as the accretion rate are remarkably similar between the two approaches,
exhibiting only $\sim 10 \%$ differences in average. The approximate metric is
five times more efficient among other computational advantages. The numerical
implementation of the metric is now open-source and optimized for numerical
work. We have also implemented this spacetime in the widely used GRMHD codes
Athena++ and EinsteinToolkit.
