Aromatic compounds have fascinated scientists and have been studied for decades due to their unique properties of extended π-conjugation and π – π stacking. This extended conjugation leads to greater stability of aromatic systems and abnormally high heat of hydrogenation compared to other open-chained conjugated systems. On the other hand, π-stacking is directly or indirectly responsible for numerous applications ranging from material sciences to drug synthesis. The importance of this π-stacking is also found in biological systems [1], [2] where the nucleobases of DNA/RNA’s double helices make π – π stacking to stabilize the structure. It is also shown that hydrogen bonding ability depends upon the stacking orientation of aromatic nucleobases [3]. Another example of π-stacking can be found in globular proteins containing phenylalanine, tyrosine, and tryptophan amino acids, all of which have aromatic rings as side chains. Studies have shown that around 60% of aromatic side chains in the proteins show aromatic pairing, 80% of which form three or more interacting aromatic side chain networks [4]. In recent studies, the π-stacking between HOFs (Hydrogen-bonded Organic Frameworks) has been found to stabilize the materials not only in common polar solvents but also in acidic and basic mediums [5]. Numerous potential applications of such aromatic systems have led to many studies over half a century both on theoretical level [6], [7], [8], [9], [10] and experimentally [11], [12], [13], [14]. Despite making tremendous strides in this field, most theoretical studies have focused on the gaseous phase. Studies of such systems in condensed phase to understand the unique π-bonding properties are limited [15].

Fluorine, oxygen, and nitrogen are the three most electronegative atoms responsible for the vast majority of H-bonding in nature. For aromatic systems, the three most basic molecular structures we can think of with those electronegative atoms are aniline with NH2 attached, phenol with OH attached, and fluorobenzene with –F attached. Large number of aromatic compounds found both in nature and industries are derivatives of these three primitive molecules. Studying their liquid structure and intermolecular interactions is essential to understand their properties at the molecular level. The molecules also differ vastly from each other in their physical and chemical properties. Although studies have already been done for aniline [16] and phenol [17], [18], [19], [20], [21] in the gaseous state, their properties in the liquid phase are still not studied thoroughly. The studies of such arene systems in the liquid phase are thus essential to figure out how the π-interaction between two molecules is affected by adding different electronegative atoms and whether these molecules form hydrogen bonds with other aromatic hydrogens. The interaction between electronegative atoms and the π-electron cloud is also a subject of interest. Molecular dynamics (MD) studies are an integral tool to explore their properties and interacting behaviors in the liquid state. Force field based MD allows the simulation of large systems. In some cases, studies using ab initio MD, where the forces are calculated on the fly, may provide additional insights. Since these aromatic systems have large π-electron clouds and strong electronegative atoms, electron cloud diffusion and bond polarization will happen to a great extent. Hence, ab initio MD is important for such liquid phase studies, where the calculation of bonding between atoms, charge transfer, and polarization are more accurate than the force-field based MD.