Because extending the dimensionality beyond two coordinates is difficult from both the implementation and parameterization perspective, Brooks and co-workers introduced a multisite representation, where a separate λ i,k -coordinate is assigned to each physical state k of a titratable group i. For a residue with multiple "chemically"-coupled titratable sites each λ i,k -coordinate has the same state at λ i,k = 0, while at λ i,k = 1 the group is in one of the n i possible protonation states (i.e., state k) of residue i. The state at λ i,k = 0 (for all coordinates k of residue i) is the same for all λ-coordinates but does not correspond to a physical protonation state of the residue, and neither do states for which the sum of the λ i,k -s is not equal to 1 (Figure ). To restrict sampling to the (hyper-)plane connecting the physical states, the sum of λ i,k -s is constrained ( k λ i,k = 1). Since the λ-dynamics implementation in a fork of GROMACS relies on linear λ-coordinates, rather than the auxiliary circular coordinates that would fulfill the constraints by construction, we apply a constraint on the sum of λ-s. To efficiently apply this constraint we use an analytical expression to solve a generalized version of charge constraint introduced by Donnini et al. (see Appendix).