Accurate predictions of free energies of binding of ligands to proteins are challenging partly because of the nonadditivity of protein ligand interactions; i.e., the free energy of binding is the sum of numerous enthalpic and entropic contributions that cannot be separated into functional group contributions. In principle, molecular simulations methodologies that compute free energies of binding do capture nonadditivity of protein ligand interactions, but efficient protocols are necessary to compute well-converged free energies of binding that clearly resolve nonadditive effects. To this end, an efficient GPU-accelerated implementation of alchemical free energy calculations has been developed and applied to two congeneric series of ligands of the enzyme thrombin. The results show that accurate binding affinities are computed across the two congeneric series and positive coupling between nonpolar R-1 substituents and a X = NH3+ substituent is reproduced, albeit with a weaker trend than experimentally observed. By contrast, a docking methodology completely fails to capture nonadditive effects. Further analysis shows that the nonadditive effects are partly due to variations in the strength of a hydrogen-bond between the X = NH3+ ligands family and thrombin residue Gly216. However, other partially compensating interactions occur across the entire binding site, and no single interaction dictates the magnitude of the nonadditive effects for all the analyzed protein ligand complexes.