Polychlorinated Biphenyls Interactions with Water—Characterization Based on the Analysis of Non-Covalent Interactions and Energy Partitioning

: Although polychlorinated biphenyls (PCBs) have been long banned from use, both they and dioxins are still considered persistent organic pollutants. The reason is twofold: their relative inertness (especially to oxidation) and their ability to accumulate in fat tissue. The current study sheds light on the interactions of PCBs with water, chlorine, and chlorine dioxide. Necessary insight is gained from Atoms in Molecules (AIM) and Non-Covalent Interactions (NCI) index and analyses for the 1:1 complexes of PCBs with water and chlorine (Cl 2 ) molecules. Further, Symmetry-Adapted Perturbation Theory (SAPT) calculations reveal the strength and nature of the intermolecular interactions, and the presence of halogen bonding is demonstrated in AIM, NCI, and SAPT studies. The stability of water, chlorine, and chlorine dioxide (ClO 2 ) complexes with PCBs is discussed using the supramolecular MP2 approach. Finally, analysis of microsolvation shells of PCBs showed the origins of the hydrophobicity and environmental persistence of these chemicals. Our results are applicable to the sustainability of water treatment strategies providing a description of forces and interactions at the molecular level.


Introduction
Persistent organic pollutants (POPs), according to the text of the United Nations Stockholm Convention of 22nd May 2001, are chemicals that possess toxic properties, resist degradation, bioaccumulate, and are transported, through air, water, and migratory species, across international boundaries and are deposited far from their place of release, where they accumulate in terrestrial and aquatic ecosystems [1]. Although the list of POPs is open to additions and amendments, the most widespread classification of the considered chemicals is known as "the dirty dozen": aldrin, chlordane, DDT, dieldrin, endrin, heptachlor, hexachlorobenzene, mirex, toxaphene, polychlorinated biphenyls, polychlorinated dibenzo-p-dioxins, and polychlorinated dibenzofurans [2]. It is easy to observe that most of these chemicals contain chlorine atoms in their molecular structure. This fact adds to their generally hydrophobic nature, leading to chemical stability, inertness, and bioaccumulation in fatty tissues. The presence of -Cl atom in the structure of the discussed compounds increases the spectrum of possible interactions with the aquatic environment. In addition to intermolecular hydrogen bonds (HBs) and electrostatic interactions, there can also be halogen bonds. The van der Waals forces present in interactions of chlorine-containing compounds with other species are also enhanced by large electron clouds of chlorine atoms. This influences the intermolecular interactions of POPs with the environment, including water. In turn, these interactions are, in the long-term, responsible for the fate of POPs in the environment, including transport and bioaccumulation. The above reasons prompted us to investigate the diversity of non-covalent interactions [3] in a selected group of POPs-polychlorinated biphenyls. Among the oxidizing gases for water treatment, chlorine, ozone, and chlorine dioxide, the latter is the least common, despite some interesting properties. Chlorine dioxide, ClO 2 , is a gaseous oxidant used in water treatment as an alternative to Cl 2 [13]. It possesses some advantages over Cl 2 , e.g., it is effective not only below neutral pH, but also in moderately alkaline waters, and its action on water impurities generates lesser amounts of (possibly carcinogenic) halomethanes [14]. On the other hand, ClO 2 can undergo disproportionation into ClO − 2 and ClO − 3 , which limits its use in various jurisdictions [15]. PCBs are resistant to oxidation with either Cl 2 or ClO 2 ; however, there is a possibility of the formation of halogen bonds between these species. Halogen bonding with PCB chlorine atoms acting as donors has been implicated in the modulation of enzymatic activity of iodothyronine deiodinases, with possible result of thyroid activity disruption [8]. As already noted in the first paragraph of the Introduction, investigation of halogen-bonded (as competitive to hydrogen-bonded) systems involving PCBs would shed light on the environmental fate of these persistent organic pollutants in a variety of environments, including contaminated water and wastewaters, and thereby contribute to improving the sustainable management of water resources.
Modeling of PCB distribution, transport, and environmental fate has been carried out from diverse points of view. From the macroscale, they focused on the distribution of PCBs in water, sediments, and organisms [16,17], including phytoremediation techniques [18]. Models of interactions between PCBs and soil, as well as attempts to determine the efficiency of PCB removal using steam, use water/organic phase partition coefficients [19,20]. Such partition coefficients and Henry law constants can be subjected to the quantitative structure-activity relationship investigations [21], and theoretical chemistry approaches have been used to assess these parameters linking molecular structure with macroscopic properties [22,23]. A quantum-chemical study of atmospheric oxidation of PCB-47 has shown the potential of theoretical methods in elucidating the fate of PCBs in the gaseous phase [24]. On the other hand, the role of weak hydrophobic (dispersion) interactions between a PCB molecule and a fluorescent sensor molecule, benzopyrene, was successfully applied in a new approach to PCB detection in a water environment [25]. Moreover, PCBs, regarded as relatively inert compounds, nevertheless have been shown to undergo metabolism in human cell lines leading to dechlorinated species [26]. This overview of the literature shows that the molecular-level description of the PCBs intermolecular interactions has a profound impact on their environmental fate, but quantum-chemical investigations have not yet deeply addressed this subject. In the current study, aimed to tackle this issue, PCBs interactions with water, chlorine, and chlorine dioxide are discussed based on quantum chemistry approaches. The Atoms in Molecules (AIM) theory [27] and Non-Covalent Interactions (NCI) index [28,29] were employed for non-covalent interactions characteristics and electron density distribution and estimation of the hydrogen/halogen bond energy. Symmetry-Adapted Perturbation Theory (SAPT) [30] was applied for the energy decomposition to give insight into energy components and the strength of the interactions.
Our studies will provide an overview of the nature and energy scale of interactions between PCBs and water, as well as chlorine and chlorine dioxide, which will enable the rational design of water purification substances. The scientific value of this study is contained in the detailed description of the differences in the interaction of PCBs with water and chlorine-containing species. The molecular-level explanation of the propensity of persistent organic pollutants to accumulate in a hydrophobic environment is a step toward a broader understanding of the problem of water pollutants. The current society is searching for new forms of water treatment in view of green chemistry. Therefore, our study indicates perspectives on the use of computational chemistry to this effect. The design of catalysts and scavengers capable of handling relatively inert POP molecules will benefit from the employment of computational strategies.

Structure Description and Energy Partitioning
Initial structures of selected PCBs, water, chlorine, and chlorine dioxide were drawn with the Avogadro 1.2.0 program [31] and optimized using the default MMFF94 force field. Further quantum-chemical calculations on the PCBs and the small molecules were carried out at the RI-MP2/def2-TZVP level of theory [32][33][34] using the ORCA 5.0.3 (Max-Planck Institut für Kohlenforschung, Mühlheim an der Ruhr, Germany) software [35,36]. The geometry optimization was followed by vibrational spectrum calculation with the numerical finite-difference method at the RI-MP2/def2-TZVP level of theory as well to confirm the nature of the located stationary points; in all cases, no imaginary wavenumbers were found.
Microsolvated models were generated using the classical molecular force-field (GAFF) [37] and atomic charges from the AM1-BCC model [38] based on the Amber-Tools2021 and Amber20 software (University of California, San Francisco, USA) [39]. Each PCB molecule was then solvated using the TIP3P water model [40] with 8 Å solute-cell boundary buffer distance, resulting in five cubic simulation boxes of 33 Å edge. Each solvated system underwent the following procedure: 1000 steps of steepest descent minimization, 20 ps of NVT thermalization at 300 K, 80 ps of NpT equilibration at 300 K and 1 atm, and finally, the NVT production run of 5 ns. The NpT equilibration phase led to the final density of the simulation cells varying from 0.987 to 0.991 g cm −3 . The time step of 1 fs and non-bonded cutoff of 10 Å were used in the classical MD runs. The classical MD trajectories were used solely to generate the microsolvation models and their results were not further analyzed. Post-processing of these trajectories was performed with the use of the VMD 1.9.3 program (University of Illinois at Urbana-Champaign, Urbana, IL, USA) [41]. The final frames of each of the classical MD production runs were then used to extract 1:1 PCB-water complexes, two for each PCB: one with a water molecule along the line provided by the C-Cl bond (in later sections denoted by the number 1 in the abbreviated designation of the complexes), and one with a water molecule above the aromatic ring (denoted with the number 2). A simple atom substitution (chlorine in place of the oxygen atom in water) was then used to construct PCB-Cl 2 and PCB-ClO 2 complexes manually. These three series of structures were subjected to structure optimization at the RI-MP2 level, as described in the first paragraph of this section. After confirming the correctness of the nature of the stationary points (no imaginary wavenumbers in the vibrational spectrum), wavefunction files (WFX format) of non-radical singlet ground state complexes were generated for further analysis using RI-MP2 relaxed electron densities. The RI-MP2/def2-TZVP supramolecular binding energies were corrected for the Basis Set Superposition Error using the counterpoise scheme [42]. The nature of interactions in the complexes 1:1 PCB-water and PCB-Cl 2 was investigated using the Symmetry-Adapted Perturbation Theory (SAPT) method [30]. The structures optimized previously at the RI-MP2/def2-TZVP level of theory, as described above, were taken and subjected to the density-fitting accelerated SAPT [43] interaction energy calculation at two levels: SAPT2 level with the same def2-TZVP basis set (reported in the Supplementary Materials only) and SAPT2+ level with Dunning aug-cc-pVDZ basis set (reported in the main text). The SAPT energy terms are defined in detail in the original SAPT review article [30] as well as in the paper providing implementation details in the PSI4 package (see Ref. [43], Equation (8) for SAPT2 and Equation (9) for SAPT2+); in particular, SAPT2+ contains corrections to the dispersion term resulting from the firstand second-order treatment of intramonomer electron correlation. The energy partitioning used in the SAPT framework consists of the following specific terms of well-defined physical meaning: • electrostatic energy defined as Coulombic force between the molecular skeletons and electron densities of the non-interacting, unperturbed monomers, without polarization effects, • Pauli repulsion exchange energy, which is a purely quantum-mechanical effect related to the orbital overlap, calculated using the unperturbed wavefunctions of the monomers, • induction term gathering the effects of mutual polarization of the monomers on the electron density distribution within them, including corrections to the electrostatic and exchange terms, • dispersion energy related to the instantaneous multipoles-fluctuations of the electron density-between the monomers.
These four terms, when combined, provide the fifth value of the total interaction energy at the SAPT2 level. The interaction energy is calculated by definition using the structures of the monomers taken from the dimer so that the effects of structural relaxation are not present. The SAPT calculations were prepared with PSI4 1.3.2 program [44]. The last considered approach to the structural investigations was to generate larger microsolvation models by including the water shell of molecules closer than 3.5 Å to the solvated PCB molecule in five snapshots taken (at regular intervals of 1 ns starting at 0.5 ns of the production run) from classical molecular dynamics for each of the five PCBs, obtaining 25 microsolvated models. An example of such a shell is shown in Figure 1. This postprocessing step and preparation of the corresponding figure were carried out with the VMD 1.9.3 program [41]. The PCB-water interaction energy for the generated microsolvation models was then calculated at the RI-MP2/def2-TZVP level. These quantum-chemical calculations were performed with the ORCA 5.0.3 package [35,36]. An example of the microsolvation shell of PCB169 as a solute-snapshot from the classical MD run. Thicker lines denote the PCB and water molecules closer to the solute than 3.5 Å. Color coding for atoms: gray-carbon, white-hydrogen, red-oxygen, green -chlorine.

Non-Covalent Interactions (NCI) Description
Atoms in Molecules (AIM) theory [27] was applied in the study to investigate electron density distribution and molecular features derived from it. The electron density and its Laplacian were computed at Bond Critical Points (BCPs) involved in the non-covalent interactions. Additionally, the potential electronic energy density was calculated to give an insight into the strength of the interaction. Finally, to estimate interaction energy the Espinosa equation [45]: E int = − 1 2 V CP (r) was used, where V CP (r) is the potential energy density at the relevant BCP. The AIM analysis was performed with the assistance of the MultiWFN program [46,47]. Intermolecular, non-covalent interactions are manifested by bringing the atoms in closer than expected contact, which is usually assessed on the basis of the conventional van der Waals atomic radii [48]. This approach is well established, especially in crystal diffraction studies, where the accurate determination of the hydrogen atoms' positions can often be challenging. Computational approaches provide a unique opportunity to appeal to the properties of electron density itself. This forms the basis for the Atoms in Molecules (AIM) theory [27], and this framework is employed in the current study for the 1:1 PCB complexes with water and chlorine.
As a complementary method to the AIM, the well-known Reduced Density Gradient (RDG) analysis was performed [28,29] using the MultiWFN 3.8.dev program (Beijing Kein Research Center for Natural Sciences, Beijing, China) [46,47]. It is a visualization index showing interatomic interactions based on the electron density and its reduced gradient. It is possible to estimate the strength of the non-covalent interactions in three-dimensional space. Non-Covalent Interactions Index (NCI) was employed to detect all possible secondary bonds in the studied complexes. The main building block of this method is the RDG, s(r) [49], a dimensionless quantity which was introduced in the context of generalized gradient approximation (GGA) to the exchange energy in DFT and serves as a correction to the electron uniform gas model [50]. The properties of RDG have been thoroughly tested and, on the basis of small values of it and electronic density, it is possible to easily separate non-covalent interactions occurring in the studied system from covalent ones. A further distinction of the NCIs is performed using the sign of λ 2 -non-bonded and bonded interactions are characterized by positive and negative values of λ 2 , respectively. The pictorial representation is built upon this idea and while repulsive interactions are presented using red color, the attractive ones can be discerned by its characteristic green to blue (when stronger) RDG isosurfaces.
The proposed theoretical approaches offer the advantage of deeper insight into the nature of the intermolecular interactions of the PCBs. A classical supramolecular approach to calculating the intermolecular energies provides only a general information on the stability of a given structure, while the SAPT methodology allows for obtaining physical partitioning of the interaction energy into terms that can be interpreted from the physicochemical point of view. The AIM theory is a conceptual framework to analyze the strength and nature of the bonding interactions present in the system using an observable quantityelectron density. The NCI index adds the possibility of three-dimensional visualization of regions involved in diverse types of interactions, expanding the AIM capability so that not only the existing, but also possible interactions can be investigated. The advantages of these methods are also their accuracy and relatively low computational complexity coupled with the possibility of using freely available software tools.

Non-Covalent Bonding Analysis of the Studied 1:1 PCB Complexes
The AIM theoretical framework detects interatomic interactions by locating Bond Critical Points (BCPs) and associated bond paths linking the interacting atoms through the BCPs. Properties of bonding are revealed in the values of electron density and related parameters at the BCP. Pictorial representation of the electron density topology, including the BCPs as small green spheres, is shown in Figure 2 for the PCB-water complexes, and the parameters of the relevant intermolecular BCPs are provided in Table 1. For PCB-chlorine complexes, similarly arranged Figure 3 and Table 2 are also presented and discussed. The Figures and Tables use a shorthand labeling of the PCB complexes: letters from a to e denote the PCB according to Scheme 1, the digit 1 or 2 differentiates between two variants of the complex.
The complexes of PCB15 and PCB11 (each ring has only one chlorine substituent) with water exhibit a preference of the water molecule to interact rather with -H than -Cl functions, even if the geometry optimization started from the O-H· · ·Cl-C arrangementsee the structures a1 and b1 in Figure 2. The more chlorine-rich PCBs are more likely to form water-chlorine atom interactions, visible in the structures c1, d1, and e1. Systems with classical, medium, or strong hydrogen bonds (HBs) prefer arrangements in which the role of HBs is maximized; here, we observe that the competitive and alternative bond paths are formed, e.g., in the a2 structure: both O-H· · ·C and C-H· · ·O contacts are visible. Further, the presence of at least two closely positioned -Cl atoms in PCB77 and PCB169 leads to interesting bifurcated O-H· · ·Cl-C contacts, where two chlorine atoms are engaged by the same hydrogen atom of the water molecule-see the structures d1 and e1 in Figure 2. The presence of BCPs and bond paths indicates, however, merely the existence of interaction, while more details can be inferred from the numerical values of the relevant BCP properties. Noteworthy is the fact that averaged sums of estimated energies for all complexes are equal to 1.222 kcal/mol for PCB-water and 1.539 kcal/mol for PCB-chlorine ones. As seen in Table 1, the electron density ρ at the BCPs is very low and does not exceed 0.01 e · a −3 0 . In each case, however, the value of ρ is larger than the lowest boundary, 0.002 e · a −3 0 , assumed in the criteria of Koch and Popelier [51]. The values of electron density Laplacian are positive, according to expectation, but in many cases are lower than the threshold of 0.024 e · a −5 0 given in [51]. The simple but effective approach of Espinosa [45] allowed us to estimate the contributions of individual interactions to the total interaction energy between the PCB and water molecule. The E1 values given in Table 2 are between 0.617 and 1.937 kcal mol −1 , placing the intermolecular HBs at the lower (weaker) range of energies.
The PCB-chlorine complexes also require the PCB to be chlorine-rich to exhibit significant Cl· · ·Cl contacts, as shown in Figure 3 (see especially the structure e1). On the other hand, there is a specific class of interactions, labeled in Table 2 as H-C· · ·Cl and indicated in Figure 3 by the bond path between a carbon atom of PCB and one of the Cl atoms of the Cl 2 molecule. In these cases, the electron density and its Laplacian at the corresponding BCP are consistently higher than for the weak C-H· · ·Cl hydrogen bonds-e.g., for the d2 structure, the corresponding values of ρ are 0.014 and 0.005 e · a −3 0 . This also means that such interactions of electron clouds of the Cl atom and the aromatic ring contribute ca. 2.5 kcal mol −1 to the interaction energy, according to the Espinosa model [45]. Returning to the structure e1, exhibiting bifurcated Cl· · ·Cl halogen bonds, these interactions are similar in strength and in topological parameters to the weak HBs noted for the PCB-water complexes. The contributions of the two halogen bonds, calculated with the Espinosa equation, combine to the value of 1.738 kcal mol −1 , in a surprisingly good agreement with the SAPT2+ interaction energy of −1.732 kcal mol −1 reported in the next section. To sum up, one can conclude that the overall picture of AIM investigations for 1:1 models reveals that the Cl 2 molecules compare favorably to the H 2 O when it comes to interacting with various PCBs.
Let us start the discussion regarding non-covalent bonds present in the studied systems in the context of the Non-Covalent Interaction index (NCI) method. It should be noted that this method lacks the quantitative treatment of the system-we cannot easily estimate the interaction energies of the studied interactions; thus, only qualitative discussion can be provided. For H 2 O-PCB complexes (see Figure 4), the more diverse, although weaker in strength when compared to Cl 2 -PCB systems, the pattern of non-covalent interactions can be seen. Aside from weak hydrogen (C-H· · ·O-H or H-C· · ·H-O) and halogen bonds that are present in a2, b1, b2, c2, d2, e1, e2, and c1, O-H· · ·Cl interactions in c2, d2 systems can also be observed. Regarding Cl 2 -PCB complexes, one can note that the vast majority of interactions presented in Figure 5 are of weak van der Waals type. In virtually every presented complex, the donor-acceptor interactions between Cl 2 and the aromatic rings were present-on the other hand, hydrogen bonds of type R-C-H· · ·Cl, were also detected in structures b1, c1, c2, d1, and d2. It is visible that the substitution of the aromatic rings with more than one chlorine substituent leads to the beefing up of the π-Cl interactions-most probably via the change of the quadrupole moment of the aromatic motif from negative to positive value. Interesting, and unnoticed by AIM analysis, was the presence of repulsive intramolecular C-H· · ·H-C between benzene-like frameworks and Cl· · ·Cl interactions. Subfigure labels (a-e) identify the PCB compound according to Scheme 1, while numbers 1 and 2 differentiate between two initial structures subjected to optimization. Color coding for atoms: gray-carbon, white-hydrogen, red-oxygen, green-chlorine. Bond critical points are labeled as small green spheres. Intermolecular bond paths are drawn as grey lines. Table 1. AIM-derived properties at BCPs of the investigated 1:1 complexes of chosen PCBs and water molecules. E1 is a bond energy computed using the Espinosa equation. The dimensions of given quantities are as follows: electron density, ρ BCP , is given in e · a −3 0 atomic units and the Laplacian of electron density, ∇ 2 ρ BCP , is in e · a −5 0 units. V CP denotes potential energy density.

Interaction Energy Analysis in 1:1 PCB Complexes with Water, Chlorine, and Chlorine Dioxide
Bonding and interaction energies of the singlet spin state 1:1 complexes of PCBs with water and Cl 2 molecules are reported in Table 3. Bonding energies obtained at the RI-MP2/def2-TZVP level are calculated with respect to the optimized monomers, while the SAPT2+/aug-cc-pVDZ interaction energies do not include structural relaxation. Additional results from the SAPT2/def2-TZVP level are provided in Table S1 of the Supplementary Material, and they follow the same trends as the values in Table 3. The most general observation related to these energy values concerns the comparison of H 2 O and Cl 2 as the interaction donors/acceptors. It is seen that, in general, complexes of the studied PCBs with chlorine are stronger than the corresponding PCB-water systems. SAPT2+, besides similar treatment of electrostatics, induction and exchange terms (at the secondorder perturbational level with respect to intramonomer electron correlation) to SAPT2, treats dispersion at the level similar to MP4 and, thus, provides a more accurate estimation of that contribution. The importance of the more thorough level is especially visible when one considers results for chlorine complexes, where dispersion contributions are definitely the most important part of the interaction energy. Providing a more robust discussion of the presented complexes, let us start with the comparison of obtained sums of interaction energies with the binding energy provided by the supramolecular method. It can be seen that for the PCB-water series of the complexes, the RI-MP2 based approach systematically underestimates binding energy when compared to SAPT2+/aDZ level, while for the PCB-Cl 2 series, both under-and overestimation occurs. Moreover, when one considers E elst +E exch as the antisymmetrized reference state electrostatic energy, it will be evident that its sum is certainly not attractive but even repulsive. Thus, as one may observe in Table 3, the most important contributions to the interaction energy stem from the dispersion and the induction energy for both investigated series of complexes. It is especially noticeable for Cl 2 complexes, where, due to significantly larger van der Waals radii of chlorine atoms (and, thus, larger intermonomer charge density overlaps), the shortrange exchange interactions dominate above all other terms. The same argument can be put forward when the dispersion contributions will be taken into consideration-generally, their contribution to the interaction energy of Cl 2 complexes is almost twice larger in magnitude than its counterpart of H 2 O complexes. For c1, d1, and e1 water complexes, one can observe significantly smaller values of exchange energy-in these cases, water molecules are located in the plane with the associated PCBs and form hydrogen (c1, d1) and halogen bonds (e1), respectively. Similar findings were obtained for e1, Cl 2 -PCB complex, where also the halogen bond is formed and, due to lower overlap of atomic densities, the exchange and dispersion contributions are significantly lower in magnitude.
A final remark is necessary for the PCB-ClO 2 complexes, which are spin doublets. For this reason, only their binding energies, evaluated at the RI-MP2/def2-TZVP level, are reported. Severe difficulties in localizing the minima on the potential energy surface have emerged, and only five structures were finally optimized. Their binding energies, grouped at the bottom of Table 3, are negative (unfavorable interactions) and on the order of 25-30 kcal mol −1 . Such values indicate that the ClO 2 molecules are not forming stable structures with the PCBs. The spin density, according to Mulliken and Löwdin population analysis schemes, stays totally on the ClO 2 molecule. One cannot, therefore, hope for effective contact of PCBs with the oxidizing agent, chlorine dioxide. Table 3. Binding or interaction energies in kcal mol −1 between PCBs and water, chlorine or chlorine dioxide in the studied 1:1 complexes optimized at the MP2/def2-TZVP level of theory. Results of supramolecular MP2/def2-TZVP (binding energy with counterpoise correction) and perturbational SAPT2+/aug-cc-pVDZ (interaction energy) approaches. The contributions summing up to the total SAPT2+ energy are: Elst = electrostatics; Exch = Exchange; Ind = Induction; Disp = Dispersion.

Interactions between PCBs and Their Solvation Shells
High persistence of PCBs in the environment and their ability to bioaccumulate are attributed, among other factors, to their hydrophobicity and generally weak interactions with water. The latter fact was investigated in detail with the combination of classical MD, serving as the source of a representative ensemble of PCB solvation shell structures, and RI-MP2/def2-TZVP interaction energy calculations. As already noted in the Computational Methodology, the extent of the microsolvation shell was chosen to encompass water molecules with the nearest distance to the PCB less or equal to 3.5 Å. This ensured that the shell covered the first solvation sphere of the PCB, without additional water molecules not in direct interaction with the solute. For each of the five PCBs, five classical MD snapshots were analyzed-their structures are presented in Figure S1 of the Supplementary Material. As seen in Table 4, the number of water molecules surrounding a given PCB was between 21 and 35, with a median of 29. No obvious correlation between the number of chlorine atoms in the PCB molecule and the size of the microsolvation shell is visible in the gathered data. While the total interaction energies at a given snapshot are relatively large, between −7.8 and −18.3 kcal mol −1 (with one outlier, an exceptionally weak interaction of −0.52 kcal mol −1 ), very low values of the interaction energy "per molecule" are seen, i.e., E int /n, where n is the number of water molecules in the microsolvation shell. Such an average provides an estimate of the contribution of each single water molecule when the shell is built around the solute. The contributions gathered in Table 4 are indeed small, ca. 0.5 kcal mol −1 stability gain for each attached water molecule, in comparison to the electronic energies (no zero point and temperature effects included) of classical water dimer structure (−5.0 kcal mol −1 [52]) and a set of other possible diverse (H 2 O) 2 structures (from −5.0 to −2.4 kcal mol −1 [53]). The reported interaction energy values explain the hydrophobicity of PCBs-preparation of a cavity for the solute molecule, together with the associated breakup of the water-water hydrogen bonding network, is energetically more costly than the stability gain due to the PCB-water interactions. The value of interaction energy per one water molecule is also consistently several times smaller than most of the values reported in the SAPT study where optimized structures were employed (see Table 3). This indicates that interaction "hot spots" are scarce on the boundary of the PCB molecules-thus, most of the microsolvation shell water molecules are passive towards the solute. The net result is high hydrophobicity and inertness of the PCB molecules in a water environment, a molecular-scale origin of the macroscopic behavior of PCBs as persistent organic pollutants.

Conclusions
We examined five PCBs interactions with the water molecule, chlorine, and chlorine dioxide. Subsequently, the models with solvation shells were constructed to reproduce the water environment. SAPT method has allowed the energy decomposition giving insight into the energy components. It was found that the dispersion is the most significant stabilizing factor in the studied complexes. The use of AIM and NCI methods enabled us to look at the strength and type of interactions on the basis of electron density. A variety of weak non-covalent interactions have been detected and their molecular features have been determined at BCPs. The strength of the interactions was estimated using Espinosa equations as well. The applied methodologies consistently show that the PCB-Cl 2 complexes are slightly, but unmistakably stronger than the PCB-H 2 O systems; the average binding energy of the former is equal to 1.539 kcal/mol and for the latter equals 1.222 kcal/mol.
The theoretical findings of the current study provide a quantitative description of the PCBs interactions with various compounds, and highlight the molecular scale origins of the persistence of PCBs in the natural environment. The results of the research can be used in the design of new compounds and methods for water purification. Such methods must tackle two problems: • PCBs do not interact with water and accumulate in non-polar environments, • the relatively inert C-Cl bonds must be activated.
It will hopefully be possible in near future to carry out a de novo design or directed evolution aimed at artificial enzymes recognizing aromatic molecules with halogen functions, activating the C-Cl bonds, and transforming the PCBs into hydrocarbons and chloride ions. Such enzymes will probably employ halogen bonding as one of the PCB recognition techniques in its active site. A currently available alternative is a two-phase system, in which the PCBs are washed out from the aqueous phase into the organic hydrophobic medium, and catalysts carry out the PCB degradation. The design of such catalysts should take into account the fact that PCBs interact via weak non-covalent interactions.
Supplementary Materials: The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/su141912529/s1, Table S1: Binding or interaction energies in kcal mol −1 between PCBs and water or chlorine in the studied 1:1 complexes. Results of supramolecular MP2 (binding energy) and perturbational SAPT2 (interaction energy) approaches with def2-TZVP basis set; Figure S1: Microsolvation shells-snapshots from the classical MD runs of the studied PCBs in water, used in the solvation shell interaction energy study at the RI-MP2/def2-TZVP level; Archive file: atomic coordinates of the investigated stable complexes optimized at the RI-MP2/def2-TZVP level.

Data Availability Statement:
The data presented in this study are available in the article itself.

Acknowledgments:
The authors gratefully acknowledge generous grants of computer time provided by the Academic Computing Centre Cyfronet-Kraków (Prometheus supercomputer, part of the PL-Grid infrastructure) as well as by the Poznań Supercomputing and Networking Center (PSNC). Additionally, A.J. and J.J.P. acknowledge the ARCHER2 UK National Supercomputing Service (https: //www.archer2.ac.uk, accessed on 5 March 2022) for generous CPU time and facilities within the framework of the DECI-17 access programme.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: