A (P)HMC algorithm for N_f = 2 + 1 + 1 flavours of twisted mass fermions

We present a detailed design of a (P)HMC simulation algorithm for N_f = 2 + 1 + 1 maximally twisted Wilson quark flavours. The algorithm retains even/odd and mass-shift preconditionings combined with multiple Molecular Dynamics time scales for both the light mass degenerate, u and d, quarks and the heavy mass non-degenerate, s and c, quarks. Various non-standard aspects of the algorithm are discussed, among which those connected to the use of a polynomial approximation for the inverse (square root of the squared) Dirac matrix in the s and c quark sector.