The calculation could be massively speedup by using the partial symmetrization technique.
This is implemented in the AtomicSymmetries.jl package, but would require that Julia completely handles the symmetries.
In this way, symmetrization could be carried out by using replicas only on the violated symmetries by the perturbation.
Even better, this feature of replica symmetries could be disabled at all, and then just applied a-posteriori by averaging equivalent response functions (instead of applying replica, we run in parallel equivalent responses and mutually correct the application after each average is taken).
This would improve by a factor N the scaling cost of the TDSCHA code.