We employ the molecular theory of solvation also known as the three-dimensional reference interaction site model with the Kovalenko–Hirata closure relation (3D-RISM-KH) to study adsorption of several heterocyclic aromatic hydrocarbons (acridine, benzothiophene, carbazole, dibenzothiophene, indole, and phenanthridine), which are intended to represent bitumen fragments on the surfaces of a single-sheet kaolinite nanoparticle in cyclohexane and toluene solvents. In addition to adsorption, we also examine solvent-mediated effective interactions between two kaolinite nanoparticles in organic and aqueous solutions. The proposed adsorption model serves as a simple prototype of suspensions formed during nonaqueous extraction of bitumen from oil sands of the Athabasca Basin. Using the proposed computational approach, excess adsorption isotherms of bitumen fragments on two different faces of kaolinite were calculated and compared in cyclohexane and toluene solvents at several temperatures. Almost all the studied molecules show strong preference for adsorption on the octahedral aluminum hydroxide surface of kaolinite. While adsorption of the molecules on the tetrahedral silicon oxide surface is weaker, it is still significant compared to the octahedral surface and for some of the studied molecules adsorption appears to be stronger than for the octahedral surface. Due to the better surface and bulk solvation properties of toluene, the adsorption of the bitumen fragments in toluene is weaker than in cyclohexane. Potentials of mean force between kaolinite nanoplatelets in cyclohexane at several temperatures were calculated using the 3D-RISM-KH molecular theory of solvation. Evaluation of effective interactions by conventional molecular simulation techniques generally requires costly free energy integration and therefore is usually avoided by adopting simple and often unrealistic semiempirical effective potentials. On the other hand, the 3D-RISM-KH molecular theory of solvation offers a fully atomistic description with all-atom force fields and notable performance advantages over molecular simulations. Obtained in this way potentials of mean force in organic solvents exhibit complex oscillating behavior and generally possess deeper main minima and smaller aggregation barriers than the corresponding potentials in aqueous solution.