In all of the simulations, the same simulation parameters recommended for the CHARMM force field were used. The Coulomb interactions were computed with the PME method with a real-space cut-off of 1.2 nm and a grid spacing of 0.14 nm. The van der Waals interactions were modeled with the Lennard-Jones potential, smoothly switched to zero in the range from 1.0 to 1.2 nm. A constant temperature was maintained at 300 K using the v-rescale thermostat with a time constant of 0.5 ps -1 , and a constant pressure was maintained at 1 bar using the c-rescale barostat with a relaxation time of 5.0 ps. For GLIC a semiisotropic barostat was used. In all simulations, the leapfrog integrator and a time step of 2 fs was used. The bond lengths to hydrogen atoms were constrained with the LINCS algorithm, and the internal degrees of freedom of the water molecules were constrained with the SETTLE algorithm. Next, details of the constant pH MD simulations constructed and performed by applying phbuilder are given. For the protein systems, the aspartic acid (Asp), glutamic acid (Glu), and histidine (His) residues were considered as titratable groups. The reference pK a values were 3.65 and 4.25 for Asp and Glu, respectively, and 6.53 and 6.92 for His HSD (singly protonated at the N δ ) and His HSE (singly protonated at the N ϵ ), respectively. For the λ-particles, the mass was set to 5 atomic units, and the temperature was maintained at 300 K by using a separate v-rescale thermostat with a time constant of 2.0 ps -1 . The barrier height of the biasing potential barrier was set to 7.5 kJ/mol.