G protein-coupled receptors (GPCRs) are membrane proteins that transduce the signals of extracellular ligands, such as hormones, neurotransmitters and metabolites, through an intracellular response via G proteins. They are abundant in human physiology and approximately 34% of the marketed drugs target a GPCR. Upon activation, the receptor undergoes conformational changes to accommodate the binding of the G protein. Our insights in the structural determinants of ligand binding and receptor activation have increased tremendously over the past decade. This has largely been a cause of the growing amount of experimentally determined structures, which provide crucial insights in ligand binding mechanisms. However, in a typical drug design project it is unlikely that such structures can be generated for all ligands. In those cases, computationally derived models of the protein-ligand complex can be generated. Rigorous free energy calculations such as the free energy perturbation (FEP) method, can subsequently provide the missing link between those structures and experimental ligand binding data, and provide further insights in the binding mechanism. In this thesis, two workflows are presented to calculate free energies of binding for ligands to wildtype (QligFEP) and mutant (QresFEP) receptors. Both methods were tested on a set of solvation free energies of sidechain mimics. QligFEP was furthermore applied on three protein-ligand binding datasets, including pair comparisons between topologically unrelated molecules (scaffold hopping). QresFEP was used to calculate protein-ligand binding affinities to mutants of the neuropeptide Y1, and to predicte the effect of receptor modifications on the thermal stability of T4 lysozyme. The remainder of this work focussed on the application of these protocols in the design, synthesis and molecular pharmacology of ligands for the family of adenosine receptor (ARs). These receptors, involved in many physiological processes such as promotion of sleep (caffeine is a well-known inhibitor), have recently been pursued as drug targets in immuno-oncology. QligFEP was used in the design of novel series of antagonists for the A3AR and A2BAR. QresFEP was used to study ligand binding to the A1AR and in a multidisciplinary approach to characterize binding to the orphan receptor GPR139. Both approaches were combined to design a series of A2AAR antagonist, and to propose a binding mode later confirmed by new crystal structures. Finally, a new application of FEP is introduced based on conformational equilibria between the active and inactive A2AAR, to elucidate the regulation mechanism of receptor activation by ligands and receptor mutations