The LS and HS geometry-optimized clusters consisting of the complexes surrounded by 17-20 water clusters were used to calculate all 3-body interaction energy terms that encompasses the cluster and 2 waters, using the same level of theory as the geometry optimizations. The combinatorial task of obtaining single point energies of the right cluster-combinations was automated with a python class available at , using ASEs ORCA calculator-interface. Basis set superposition errors were addressed by employing the full counterpoise-correction, such that electron-less basis functions were placed on all nuclear positions of atoms that do not explicitly enter into the single-point calculation at hand. For instance, to get the potential interaction energy of the dimer comprised of molecule 1 and 2, E(r 1 , r 2 ) of the cluster comprised of the complex and 17 waters, ghost atoms were placed on r 3 , . . . , r 18 .