PSD-95/discs-large/ZO-1 (PDZ) domains form a large family of adaptor proteins that bind to the C-terminal tails of their binding partner proteins. Via extensive molecular dynamics simulations and alchemical free energy calculations, we characterized the binding modi of phosphorylated and unphosphorylated EQVSAV peptides and of a EQVEAV phosphate mimic to the hPTP1E PDZ2 and MAGI1 PDZ1 domains. The simulations reproduced the well-known binding characteristics such as tight coordination of the peptidic carboxyl tail and pronounced hydrogen bonding between the peptide backbone and the backbone atoms of a β-sheet in PDZ. Overall, coordination by hPTP1E PDZ2 appeared tighter than by MAGI1 PDZ1. Simulations of wild-type PDZ and arginine mutants suggest that contacts with Arg79/85 in hPTP1E/MAGI1 are more important for the EQVEAV peptide than for EQVSAV. Alchemical free energy calculations and PaCS-MD simulations could well reproduce the difference in binding free energy between unphosphorylated EQVSAV and EQVEAV peptides and the absolute binding free energy of EQVSAV. However, likely due to small force field inaccuracies, the simulations erroneously favored binding of the phosphorylated peptide instead of its unphosphorylated counterpart, which is in contrast to the experiment.