Potentials of mean force (PMF) between ionizable amino acid side chains (Arg, Lys, His, Glu) in the head-group area of a palmitoyl oleoyl phosphatidylcholine lipid bilayer were obtained from all-atom molecular dynamics simulations and the adaptive biasing force method. Simulations in bulk water were also performed for comparison. Side chains were constrained in collinear, stacking, and orthogonal (T-shaped) orientations. The most structured and attractive PMFs were observed for hydrogen-bonded side chains. Contact minima occurred at a distance of 2.6-3.1 angstrom between selected atoms or centers of mass with the most attractive interaction (-9.6 kcal/mol) observed between Arg(+) and Glu(-) Hydrogen bonds play a significant role in stabilizing these interactions. Interactions between like charged side chains can also be very attractive if the charges are screened by surrounding molecules or groups (e.g., the PMF value at the contact minimum for Arg(+)center dot center dot center dot Arg(+) is -7.6 kcal/mol). Like charged side chains can have contact minima as close as 3.6 angstrom. The PMFs depend strongly on the relative orientation of the side chains. In agreement with experimental studies and other simulations, we found the stacking arrangement of like charged side chains to be the most favorable orientation. Interaction energies and Lennard-Jones energies between side chains, headgroups, and water molecules were analyzed in order to rationalize the observed PMFs and their dependence on orientation. In general, the results cannot be explained by simple dielectric arguments.