  We propose a new method to compute magnetic surfaces that are parametrized in
Boozer coordinates for vacuum magnetic fields. We also propose a measure for
quasi-symmetry on the computed surfaces and use it to design coils that
generate a magnetic field that is quasi-symmetric on those surfaces. The
rotational transform of the field and complexity measures for the coils are
also controlled in the design problem.
  Using an adjoint approach, we are able to obtain analytic derivatives for
this optimization problem, yielding an efficient gradient-based algorithm.
Starting from an initial coil set that presents nested magnetic surfaces for a
large fraction of the volume, our method converges rapidly to coil systems
generating fields with excellent quasi-symmetry and low particle losses. In
particular for low complexity coils, we are able to significantly improve the
performance compared to coils obtained from the standard two-stage approach,
e.g.~reduce losses of fusion-produced alpha particles born at half-radius from
$17.7\%$ to $6.6\%$. We also demonstrate 16-coil configurations with alpha loss
< $1\%$ and neoclassical transport magnitude $\epsilon_{\mathrm{eff}}^{3/2}$
less than approximately $5\times 10^{-9}.$
