A method is developed for performing classical explicit-solvent molecular dynamics (MD) simulations at constant pH, where the protonation state of each ionizable (titratable) group in a simulated compound is allowed to fluctuate in time, depending on the instantaneous system configuration and the imposed pH. In this method, each ionizable group is treated as a mixed state, i.e., the interaction-function parameters for the group are a linear combination of those of the protonated state and those of the deprotonated state. Free protons are not handled explicitly. Instead, the extent of deprotonation of each group is relaxed towards its equilibrium value by weak coupling to a "proton bath." The method relies on precalibrated empirical functions, one for each type of ionizable group present in the simulated compound, which are obtained through multiple MD simulations of monofunctional model compounds. In this study, the method is described in detail and its application illustrated by a series of constant-pH MD simulations of small monofunctional amines. In particular, we investigate the influence of the relaxation time used in the weak-coupling scheme, the choice of appropriate model compounds for the calibration of the required empirical functions, and corrections for finite-size effects linked with the small size of the simulation box.