James van Meter provided me with an optimised version of McLachlan. He states:
You are not taking full advantage of the chi=exp(-2phi) variable. There are several terms you divide by chi or chi^2 in expressions with overall factors exp(-4phi). I rewrote the BSSN equations to make these cancellations before coding. So where you have an expression of the form chi^2(A+B/chi^2), I have chi^2A+B. This gives a slight but noticeable advantage in both accuracy and performance.
I added Hamiltonian-constraint-damping terms due to Duez et al. These terms don't seem to be well-known but they are effective.
I added a Gamma-constraint-damping term due to Yo et al.
I enforce det(g)=1.
I have not tested this new version yet, but James suggests to include it into the Einstein Toolkit.