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.Comment: Contribution to the proceedings of Lattice 2005 (Algorithms and
Machines), 6 pages, 1 figur