Density functional theory and molecular statics (MS) were performed to parametrize and validate the force field, and molecular dynamics was performed to simulate cation-π interac-tions in the solvated environment close to room temperature. GAMESS and LAMMPS were used for DFT, and MS and MD, respectively. In both MS and MD, force fields for TBT and TriTBT were taken from our previous studies. 2.1 Gas phase: Density functional theory and molecular statics All reference binding energies for parameterization were calculated in the gas phase using DFT. The optimized geometries and binding energies of the systems were obtained using the B3LYP/6-311++G(d,p) level of theory with Grimm's dispersion corrections and with consideration of counterpoise correction for basis set superpositions errors (BSSE). The BSSE corrections were considered to ensure the accuracy of the calculated binding energies by accounting for artificial stabilization arising from basis set limitations. The B3LYP functional has been widely used for studying cation-π interactions due to its balance between accuracy and computational efficiency. Different exchange-correlation functions rescale the binding energy by 10-15%; 32 however, the qualitative shape of the binding energy profiles is not altered by the choice of the exchange-correlation functional form. Therefore, we did not modify these functions to calculate the open parameters. Previous studies have demonstrated that higher-level exchange-correlation functions, such as ωB97X-D/6-311++G(d,p), produce results that are comparable to those obtained with B3LYP/6-311++G(d,p). However, B3LYP/6-311++G(d,p) is less computationally demanding, making it an ideal choice for calculating interaction energies for numerous configurations. Therefore, we opted to use the B3LYP/6-311++G(d,p) level of theory for all our calculations. The binding energy calculation procedure is discussed in the parameterization strategy and verification section.