The thermodynamics of association of fullerene [C₆₀] and water-soluble fullerene derivatives, i.e., fullerols [C ₆₀ (OH)_n, n = 0, 2, 4, 8 and 12] in aqueous solutions have been studied using molecular dynamics simulations. The potentials of mean force (PMFs) bring out the tendency of aggregation of these nanostructuresin water. The extent of hydroxylation seems to have a minor effect on the depth of the contact minima (the first minimum in the PMFs). The positions of the subsequent minima and maxima in the PMFs change with the size of the solute molecules. Higher stability of the contact state of highly hydroxylated fullerols is due to the van der Waals interactions whereas intermolecular solute-solvent hydrogen bonding nearly flattens the PMFs beyond the 2nd minima for higher fullerols. The solvent contributions to the PMFs for all the solute particles studied here are positive. Entropic and enthalpic contributions to the association of solute molecules are calculated in the isothermal-isobaric (NPT) ensemble. We find that the contact pair formation is governed by entropy with the enthalpic contributions being highly unfavorable, whereas the solvent assisted and solvent separated configurations show entropy-enthalpy compensation.