electronic structure code with fused gradient kernel for multiple gradient roots

when working with NAMD code, a common request pattern for gradient is that the gradients for all active roots are requested(for gradient diff velocity rescaling)

one can imagine a fused gradient kernel where quantities like integral derivative and orbital response reusable for all roots are computed on the fly and contracted with the cphf RHS or effective density of all roots to build partial sum

is any ES code already allowing this without a modification of source code?