where Both in L BKNHC A and L BKNHC B there appears an operator with the form where ( k, i ) = ( χ, ξ ) for L$_{A}$ and ( k, i ) = ( η, ζ ) for L$_{B}$ . Again following the derivation in Appendix A, we find The function ( τ p$_{k}$ 2 m$_{k}$ ) − $^{1}$sinh [ τ p$_{k}$ 2 m$_{k}$ ] is treated through an eighth order expansion [37]. The propagators with α = A, B, C can now be introduced. One possible reversible measure-preserving inte- gration algorithm for the BKNHC chain thermostat is then In pseudo-code form, we have the resulting integration algorithm: • q → q + τ 4 p m η → η + τ 4 p$_{η}$ m$_{η}$ p$_{ζ}$ → p$_{ζ}$e − τ 4 p η m η + τ $_{4}$F$_{p}$$_{ζ}$ e − τ 4 p η 2 m η ( τ 4 p$_{η}$ 2 m$_{η}$ ) − $^{1}$sinh [ τ 4 p$_{η}$ 2 m$_{η}$ ]      : U BKNHC B ( τ 4 ) 21 U ( τ ) BKNHC = U BKNHC B ( τ 4 ) U BKNHC C ( τ 2 ) U BKNHC B ( τ 4 ) × U BKNHC A τ ) × U BKNHC B ( τ 4 ) U BKNHC C ( τ 2 ) U BKNHC B ( τ 4 ) . (57) U BKNHC α ( τ ) = exp [ τ ˜ L BKNHC α ] (56) e $^{τL$_{i}$}$p$_{i}$ = p$_{i}$e − τ p k m k + τF$_{p}$$_{i}$e − τ p k 2 m k ( τ p$_{k}$ 2 m$_{k}$ ) − 1 sinh [ τ p$_{k}$ 2 m$_{k}$ ] . (55) L$_{i}$ = ( − p$_{k}$ m$_{k}$ p$_{i}$ + F$_{p}$$_{i}$ ) ∂ ∂p$_{i}$ , (54) F ( q ) = − ∂V ∂q (53a) F$_{p}$$_{ξ}$ = q ∂V ∂q − k$_{B}$T (53b) F$_{p}$$_{ζ}$ = p 2 m − k$_{B}$T (53c) F$_{p}$$_{η}$ = p 2 ζ m$_{ζ}$ − k$_{B}$T (53d) F$_{p}$$_{χ}$ = p 2 ξ m$_{ξ}$ − k$_{B}$T . (53e)