Simulating and Predicting Adsorption of Organic Pollutants onto Black Phosphorus Nanomaterials

Abstract

Layered black phosphorus (BP) has exhibited exciting application prospects in diverse fields. Adsorption of organics onto BP may influence environmental behavior and toxicities of both organic pollutants and BP nanomaterials. However, contributions of various intermolecular interactions to the adsorption remain unclear, and values of adsorption parameters such as adsorption energies (E_ad) and adsorption equilibrium constants (K) are lacking. Herein, molecular dynamic (MD) and density functional theory (DFT) was adopted to calculate E_ad and K values. The calculated E_ad and K values for organics adsorbed onto graphene were compared with experimental ones, so as to confirm the reliability of the calculation methods. Polyparameter linear free energy relationship (pp-LFER) models on E_ad and logK were developed to estimate contributions of different intermolecular interactions to the adsorption. The adsorption in the gaseous phase was found to be more favorable than in the aqueous phase, as the adsorbates need to overcome cohesive energies of water molecules onto BP. The affinity of the aromatics to BP was comparable to that of graphene. The pp-LFER models performed well for predicting the E_ad and K values, with external explained variance ranging from 0.90 to 0.97, and can serve as effective tools to rank adsorption capacities of organics onto BP.

1. Introduction

Layered black phosphorus (BP) has become a nanomaterial star since 2014 in diverse fields, including optoelectronics, energy storage, sensor and biomedical applications, etc. [1,2,3]. Because of its superior semiconductor properties [4], BP has been igniting an upsurge in laboratory investigations. Recently, with developments in production techniques, the production costs of BP have been largely reduced to <1 $g⁻¹ [5]. In the foreseeable future, BP or derivative products are promising for massive industrial productions, and commercial or engineering applications.

BP possesses a large theoretical specific surface area of 2400 m² g⁻¹ close to that of graphene [6], and exhibits an extremely high drug-loading ability [7]. Accordingly, BP, once released into the environment, may adsorb a variety of organic and inorganic pollutants. The interactions between BP and organic/inorganic pollutants may have impacts on environmental behavior, bioavailability, toxicology and ecological risks of both BP and organic/inorganic pollutants [8,9,10], as can be inferred from previous studies concerning adsorption of organic chemicals onto carbon nanomaterials such as carbon nanotubes and graphene [11,12,13]. It is known that BP can be oxidized by oxygen that attacks the lone-pair electrons of phosphorus [14]. The adsorption of organic chemicals may protect BP against oxidation by interacting with the lone-pair electrons and influence the environmental lifetime of BP. Therefore, it is important to investigate the adsorption of organic chemicals onto BP for understanding and assessing their environmental risks, and it is also of interest to understand the distinctions in adsorption mechanisms (i.e., intermolecular interactions) between BP and the widely investigated graphene [15,16].

Adsorption energies (E_ad) and equilibrium constants (K) are essential parameters characterizing the adsorption of chemicals onto nanomaterials. Due to the large and ever-increasing number of environmental organic chemicals, it is unrealistic and impossible to exhaustively determine the values of these parameters experimentally. Quantitative structure-activity relationship (QSAR) models may serve as alternatives to the empirical values of the parameters [16,17,18,19,20]. To the best of the authors’ knowledge, there have been no available QSAR models that can be employed to predict the adsorption parameters of organics onto BP.

Compared with the comprehensive experimental adsorption data onto carbon nanomaterials (e.g., carbon nanotubes, graphene and graphene oxides) [15,21,22,23,24,25,26,27], adsorption of only nine organic compounds (benzene, dioxane, cyclohexane, acetonitrile, nitromethane, acetone, tetrachloromethane, methylene blue and congo red) onto BP has been experimentally investigated to obtain E_ad or K so far [6,28]. Thus, the development of QSARs for predicting adsorption onto BP is seriously impeded by a realistic problem of insufficient experimental data on E_ad or K. This can be solved by computations based on density functional theory (DFT) or molecular dynamics (MD) [24,25,29,30].

DFT can accurately calculate E_ad of a static adsorption configuration and has been successfully applied to probe adsorption of different organic chemicals onto carbon nanomaterials [24,25,31]. Lazar et al. [24] found that the DFT calculated E_ad values for adsorption of seven organic molecules onto graphene were in excellent agreement with the experimental ones. MD is capable of providing dynamic evolution of the interactions between nanomaterials and organic chemicals in aqueous environments, which is helpful for explaining the adsorption mechanisms. With MD calculation, Tang et al. [32] found aromatic compounds exhibited a strong preference for the edges and wrinkles of graphene oxides. Besides, MD integrated with an adaptive biasing force (ABF) method has been proved to be reliable and efficient in calculating K values and has been successfully used to evaluate adsorption affinities of different organic compounds onto carbon nanomaterials [29,30].

Theoretically, adsorption of organic compounds onto BP is governed concurrently by multiple intermolecular interactions, similar to the adsorption of various organics onto carbon nanomaterials, such as electrostatic, van der Waals, hydrophobic interactions and hydrogen bonds [12,15,16]. Polyparameter linear free energy relationships (pp-LFERs), which form a theoretical basis of QSARs and employ a set of Abraham molecular descriptors (E, S, A, B, V and L) [33], have been successfully employed to elucidate relative contributions of the different interactions to the overall adsorption onto carbon nanomaterials [16,18,34]. However, contributions of the different interactions to the overall adsorption onto BP have not been estimated yet.

In this study, the DFT and MD methods were evaluated by comparing the calculated E_ad and K values with empirical ones for the adsorption of organics onto graphene. E_ad and K values of a series of organic compounds (aliphatic hydrocarbons, benzene and its derivatives, polycyclic aromatic hydrocarbons) adsorbed onto BP were calculated by the DFT and MD methods, respectively. The objectives of this study were three-fold: (1) to develop pp-LFER models on E_ad and K from the DFT or MD results so as to fill the data gap of BP adsorption, and (2) to elucidate adsorption mechanisms of environmental organic chemicals onto BP by quantifying the relative contributions of the various interactions to the overall adsorption and (3) to unveil the differences in adsorption between BP and graphene for organic compounds by a comprehensive comparison of the pp-LFER models.

2. Computational Methods

2.1. DFT Calculation

As listed in Table S1 in supporting information (SI), 41 aliphatic and aromatic compounds with diverse functional groups were selected as model adsorbates. As these compounds were often selected in previous adsorption studies on carbon nanomaterials [16,30], the comparison on the different adsorbing nanomaterials can be implemented.

A 6 × 5 × 1 BP supercell containing 120 phosphorus atoms was built as an adsorbent model. Periodic boundary conditions were applied in the adsorption system. A box of 19.88 × 21.87 × 22.16 ų was adopted to eliminate interactions of the adsorbates with BP in neighboring periodic structures. Similar phosphorene models have been implemented to study their interactions with 2,3,7,8-tetrachlorodibenzo-p-dioxin and amino acids [35,36]. The structural details about the graphene model were given in the Supplementary Materials.

All the DFT calculation was carried out with the Dmol³ program [37,38] in the Materials Studio package. A Perdew-Burke-Ernzerhof generalized gradient approximation (GGA-PBE) [39] method was employed to calculate exchange-correlation energies, and van der Waals interactions were treated by an empirical correction method proposed by Grimme (DFT-D2) [40]. These methods have also been well implemented to investigate the adsorption of graphene and boron nitride [16,41]. A double-numeric quality with a polarization function (DNP) basis set (comparable to the 6-31G** basis set) [42,43] was used. DFT semi-core pseudopotentials treatments were employed to deal with the relativistic effects of bromine atoms [44]. A 4 × 4 × 1 k-grid mesh was adopted for the Brillouin zone integration [45] that was required for calculating the energy by DFT. Conductor-like screening model (COSMO) [46], an implicit solvation model, was utilized to mimic aqueous environments by employing a dielectric constant of 78.54 F∙m⁻¹.

E_ad was calculated as:

E_ad = E_complex − E_adsorbent − E_adsorbate

(1)

where E stands for the total energy of a configuration, including kinetic energy, static potential energy, Coulomb energy and exchange-correlation energy [47,48]; E_complex, E_adsorbent, and E_adsorbate are total energies of the adsorption complexes, adsorbents (i.e., BP and graphene) and adsorbates, respectively. The more negative is the E_ad, the stronger is the interactions between the adsorbents and adsorbates.

2.2. MD Simulation

The MD calculations were carried out with NAMD 2.12 [49]. Force field parameters for the adsorbates were generated according to the CHARMM General Force Field (CGenFF) [50] parameters with the ParamChem Web interface [51,52]. The generated parameter files contain “penalty scores” that indicate whether the parameters can be directly used. When the penalty value is >10, the parameters need to be further optimized [52]. The Force Field Toolkit (ffTK) [53] used for developing CHARMM-compatible parameters was employed to further refine the parameters for compounds with penalty scores of >10.

Force field parameters of BP were derived from Blankschtein et al. [54], which accurately reproduced experimental crystal structures and mechanical properties. The force field parameters and structural information for graphene were detailed in the Supplementary Materials. A TIP3P water model [55] adopted in the CHARMM force field was employed for simulating aqueous environments.

All the simulation was run with a 2 fs time-step in the isothermal-isobaric ensemble at 1.01325 × 10⁵ Pa and 300 K by the Langevin piston and Langevin thermostat method [56]. The particle-mesh Ewald algorithm [57] with 1.2 Å grids was employed to deal with electrostatic interactions. A cutoff for truncating non-bonded interactions was set to 9 Å according to previous adsorption simulations of compounds onto graphene and carbon nanotubes [26,30]. The system was subject to energy minimizations of 10⁴ steps and equilibrations of 0.5 ns. An MD integrated with the ABF method [58,59], as implemented in the Colvars module [60] of NAMD 2.12, was employed to compute the potential of mean force (PMF) to obtain the free energies G(Z) along a transition coordinate (Z). As indicated by Figure S1a, Z was defined as the vertical distance between the center of mass for the adsorbates and of the first layer of BP [26]. VMD 1.9.3 was adopted for the analysis and visualization [61].

Before the MD simulation, one adsorbate molecule was arranged on the top of a three-layer BP adsorbent in a periodic box of 3.6 × 3.6 × 6.5 nm³ (Figure S1b). According to previous adsorption simulations about the carbon nanomaterials [26,30], Z was estimated to range from 3 to 15 Å in this study. A single window with the Z interval was adopted as the sampling space. The force was sampled in bins of 0.05 Å width along Z. G(Z) was normalized to zero at a distance of 14 ≤ Z ≤ 15 Å, where the energies reached a plateau as there were no interactions between BP and the adsorbates (Figure S2). The convergence of G(Z) in the free energy calculation was checked as detailed in the Supplementary Materials.

The difference in the adsorption free energy (ΔG_MD) that represents the energy required to detach an adsorbate vertically from the surface, was calculated as:

ΔG_MD = G(Z_min) − G(Z_far)

(2)

where G(Z_min) and G(Z_far) stands for the minimum G value, and the value for the adsorbate far from the nanomaterial surface (Figure S2), respectively.

K (mL/g) was calculated according to the following equation that has been applied to compute K of organic compounds onto carbon nanomaterials [29,30]:

$$K=S_{\mathrm{A}}{\int }^{ }{\mathit{e}}^{-\mathit{G}\left(\mathit{Z}\right)/\mathit{k}\mathit{T}}\mathit{d}\mathit{Z}$$

(3)

where k is Boltzmann constant (k = 1.381 × 10⁻²³ J K⁻¹), T is temperature (T = 300 K), and S_A represents specific surface areas of adsorbents. The theoretical specific surface areas (S_A(BP) = 2400 m² g⁻¹; S_A(graphene) = 2630 m² g⁻¹) [6,62] were adopted.

2.3. pp-LFER Modeling

The values of the Abraham descriptors were collected from the UFZ-LSER database [63]. The model adsorbates with the calculated logK and |E_ad| values were randomly split into a training and a validation set (Table S1) with a ratio of 4:1. Multiple linear regressions (MLR) were utilized for establishing the pp-LFER models. The goodness of fit and robustness of the model were evaluated by adjusted determination coefficient (R_adj²), root-mean-square error (RMSE), leave-one-out cross-validated Q²_LOO and k-fold cross-validation Q²_kfold (5 fold, 5000 repetitions) [64]. External explained variance (Q²_ext) [64] was used to estimate the model predictive power. Application domains for the models were characterized with Williams plots [65].

3. Results and Discussion

3.1. Reliability of Computational Methods

As experimental data (E_ad and K) on adsorption of BP are of deficiency, reliabilities of the DFT and MD methods were evaluated by taking the widely investigated graphene as an example. As can be seen from Figure 1a, the DFT calculated E_ad values linearly correlate (r = 0.98) with the experimental E_ad(exp) values [24,66] for five compounds adsorbed onto graphene (Table S2), and the slope is close to 1. The absolute errors between the experimental and calculated E_ad values range from 0.1 to 0.8 kcal/mol (Table S2), which are within experimental errors (1 kcal/mol) for E_ad [66]. The E_ad values obtained from ab initio molecular dynamics for 6 organic molecules onto graphene were also in excellent agreement with experimental values with small errors (0.1~0.9 kcal/mol) [24]. Thus, it can be concluded that the DFT method in this study can give a viable estimation on the E_ad values.

It is known that macroscopically and empirically determined logK is theoretically related with the thermodynamic quantities ΔG, enthalpy (ΔH), entropy (ΔS), temperature (T) and universal gas constant (R) through the equation [67]:

logK = −ΔG/2.303RT + logC = −(ΔH − TΔS)/2.303RT + logC

(4)

where the slope term −1/2.303RT can be calculated as −0.73 mol kcal⁻¹ at T = 300 K. The term logC in Equation (4) is a constant of investigated systems consisting of water with a specific volume (V, mL) and adsorbent with a specific mass (m, g). As detailed in the Supplementary Materials, C = V/m (mL/g). From another point of view, the terms S_A and Z in Equation (3) jointly lead to the constant term of Equation (4).

The ΔG_MD and logK_cal values (Table S3) for adsorption of 15 aromatic compounds onto graphene were calculated by the MD integrated with the ABF method. logK_cal and experimentally determined logK_exp [20] for the aromatics correlate with ΔG_MD significantly, respectively:

logK_cal = −0.70ΔG_MD − 0.8, n = 15, r = 0.99, p < 0.01

(5)

logK_exp = −0.68ΔG_MD + 0.3, n = 15, r = 0.81, p < 0.01

(6)

It can be seen that the slopes in Equations (5) and (6) are close to the theoretical value (−0.73) of Equation (4), indicating Equations (5) and (6) conform to the theoretical relationship of Equation (4). As indicated by Figure 1b, logK_exp significantly correlates (r = 0.84) with the logK_cal too. Comer et al. [30] also reported that logK_cal significantly correlated with corresponding experimental values for adsorption of aromatics with a variety of functional groups onto carbon nanotubes (r = 0.90). Therefore, it can be inferred that the MD simulation coupled with ABF in this study can give reliable estimations for logK or ΔG_MD.

Although the current study indicates the experimental logK_exp and E_ad values correlate well with the MD or DFT calculated ones, it deserves mentioning that possible biases due to wrong experimental determinations could not be addressed by the models. Although the DFT or MD methods can estimate the adsorption parameters (E_ad and logK), it can be impractical to calculate the adsorption parameters for all organics owing to the enormous computational cost. Therefore, it is necessary to develop high-throughput prediction models for the adsorption parameters.

3.2. G(Z) and E_ad of the Adsorbates onto BP

The penalty scores of force field parameters for methyl-2-methylbenzoate, 2,4-dinitrotoluene, malonic acid and 1,2-dinitrobenzene are >10. To improve the accuracy of G(Z), the force field parameters are further refined by the ffTK [53], and the refined parameters are listed in Tables S4–S10.

Figure 2 (taking phenanthrene as an example) shows that variation in G(Z) along the adsorption distances exhibits a typical V-shape. It indicates that the adsorption distance of 4.4 Å between the BP surface and phenanthrene is a crest, blow and above which repulsion and attraction become dominant, respectively. G(Z_min) of the other adsorbates also occur at about Z = 4.4 Å. Azhagiya Singam et al. [26] investigated adsorption thermodynamics of a diverse set of aromatic compounds onto graphene by MD, and found that Z_min ranged from 3.5 to 3.8 Å.

ΔG_MD for all the adsorbates ranges from −3.8 to −18.8 kcal/mol in the gaseous phase and from −2.2 to −12.7 kcal/mol in the aqueous phase (Table S11), indicating BP can capture the aromatic and aliphatic compounds spontaneously. The DFT calculated E_ad values range from −7.5 to −28.7 kcal/mol and −4.7 to −24.7 kcal/mol in gaseous and aqueous phases, respectively (Table S1). The MD calculated logK values range from 2.1 to 12.5, and 1.0 to 8.1, in gaseous and aqueous phases, respectively. logK correlates with ΔG_MD or E_ad as follows:

logK_aqueous = −0.68ΔG_MD(aqueous) − 0.7, n = 41, r = 0.99, p < 0.01

(7)

logK_gaseous = −0.70ΔG_MD(gaseous) − 0.7, n = 41, r = 0.99, p < 0.01

(8)

logK_aqueous = −0.32E_ad(aqueous) − 1.0, n = 41, r = 0.86, p < 0.01

(9)

logK_gaseous = −0.48E_ad(gaseous) − 2.1, n = 41, r = 0.93, p < 0.01

(10)

It is under the expectation that the calculated logK_aqueous or logK_gaseous correlated with ΔG_MD(aqueous) or ΔG_MD(gaseous) significantly, since the K values were derived from the MD calculated G(Z) with Equation (3). Similar to Equations (5) and (6), the slopes of Equations (7) and (8) are close to the theoretical value −0.73 mol kcal⁻¹, indicating Equations (7) and (8) also conform to the theoretical relationship of Equation (4).

Besides, since E_ad(aqueous) or E_ad(gaseous) are essentially ΔH in Equation (4), and they do not carry information on ΔS in the adsorption, it is not unexpected that the statistical significance of Equations (9) and (10) as indicated by the r values, is inferior to that of Equations (7) and (8).

It deserves mentioning that the intercepts of Equations (7) and (8) are close to that of Equation (5). According to Equation (3), the intercepts are relevant with logS_A. Since the difference (0.04) in theoretical logS_A values between pristine BP and graphene is negligible, it is rational that the intercepts of Equations (5), (7) and (8) are similar. As detailed in the Supplementary Materials, the theoretical intercept of Equation (4) calculated from the MD simulations is −1.0 (Table S12). The intercepts in Equations (5), (7) and (8) are generally close to the theoretical value.

Nevertheless, the intercepts of Equations (5) and (6) are quite disparate. The disparity may lie in the difference in S_A between the MD simulation and the experiments. Compared with the theoretical S_A value of pristine graphene in the MD simulation, the S_A values in the experiment could not be uniform, let alone the logK_exp values obtained from different experimental studies [20]. For the same batch of adsorption experiments, the S_A values of graphene are within a certain distribution range. Thus, it can be inferred that the nonuniform S_A values in the experiments lead to the intercept of Equation (6).

The ΔG_MD values in the gaseous phase are 1.6~6.1 kcal/mol lower than those in the aqueous phase (Table S11) and the DFT calculated E_ad values in the gaseous phase are 1.0~6.9 kcal/mol lower than those in the aqueous phase (Table S1), indicating water plays a negative role in the adsorption. The trajectory of the MD simulations shows that water molecules form a stable hydrogen bond network on the BP surface (Figure S3) and reduce the accessibility of the adsorbates to the BP surface. Previous MD simulations on adsorption of biomolecules onto graphene oxide and TiO₂ nanomaterials also indicated the competing behavior of water [68,69].

The MD calculation was also performed to examine the binding strength of water onto BP in the gaseous phase. ΔG_MD of water is −1.7 kcal/mol, which is higher than those (−2.2~−12.7 kcal/mol) of the 41 adsorbates. Therefore, it is thermodynamically favorable for the adsorbate to replace water molecules during adsorption onto the BP surface in the aqueous phase (Figure S4), although the process requires overcoming extra energy barriers. Previous DFT calculations also indicated that water can impede the adsorption of aromatics onto graphene oxide by interacting with the hydroxy and epoxy groups of graphene oxides [32].

3.3. Adsorption Configurations

Dihedral angles between the plane formed by the aromatic rings of the adsorbates and the BP surface, vertical distances between the adsorbates and the BP surface, and corresponding interaction energies calculated from MD, are shown in Figure 3, taking phenanthrene as an example. The results for the other 40 adsorbates are shown in Table S13. The dominant adsorption configurations were identified by analyzing the interaction energy of each configuration in the trajectory, as all the accessible orientations of the adsorbates have been sampled during the G(Z) calculation. Taking phenanthrene as an example, all the configurations sampled for the G(Z) calculation are shown in Figure S5.

The most energetically favorable adsorption configuration is the adsorbate roughly parallelling to the BP in a distance of around 4.1 Å. The DFT results also show that the adsorbates are roughly parallel to the BP surface at a distance of about 4.2 Å (Table S14). The parallel adsorption configurations were also observed for adsorption of 2,3,7,8-tetrachlorodibenzo-p-dioxin, nucleobases and amino acids onto phosphorene [35,36,70].

For adsorption onto graphene, carbon nanotubes and boron nitride nanosheets, DFT or MD simulations also showed that the aromatic rings of adsorbates were nearly parallel to the nanomaterial surface [25,30,41]. Generally, the adsorption of the compounds onto these nanomaterials is mainly physically, which is dominated by van der Waals forces. It can be seen that the parallel configuration can maximize the contact areas between the nanomaterials and the adsorbates, and thus enhance the van der Waals interactions.

Besides, BP was found to deform to different extents after the adsorption. Root mean square displacements (RMSD, detailed in Table S15) between the initial and the optimized adsorption geometries of the BP were adopted to describe the deviation/deformation. RMSD of BP adsorbing acetaldehyde is 0.01 nm, and the deformation is the slightest. RMSD of BP adsorbing pyrene is 0.21 nm, and the deformation is the largest. Acetaldehyde and pyrene correspond to the highest and the lowest values of E_ad and ΔG_MD among the 41 compounds, respectively (Tables S1 and S11). Previous studies also revealed that the wrinkles of graphene and graphene oxides can promote their adsorption for organic pollutants [15,32,71].

Pearson coefficients between RMSD and E_ad or ΔG_MD in gaseous and aqueous phase range from −0.48 to −0.58, indicating that E_ad or ΔG_MD decreases with the increase of the deformation. Zhao et al. [6] also found that BP had wrinkles after adsorbing methylene blue or congo red, and the deformation can increase adsorption sites or adsorption capacities.

3.4. pp-LFER Models

The constructed pp-LFER models are as follows:

Gaseous phase: logK = 0.2 − 0.0076E − 1.1S − 1.1A + 1.4B + 1.5L
n_tra = 33, R²_adj = 0.96, RMSE_tra = 0.35, Q²_LOO = 0.95, Q²_{kfold (k = 5, 5000)} = 0.95,
n_ext = 8, RMSE_ext = 0.35, Q²_ext = 0.97

(11)

Aqueous phase: logK = −1.7 + 0.65E + 0.75S + 0.048A − 0.095B + 4.0V
n_tra = 33, R²_adj = 0.87, RMSE_tra = 0.46, Q²_LOO = 0.83, Q²_{kfold (k = 5, 5000)} = 0.82,
n_ext = 8, RMSE_ext = 0.43, Q²_ext = 0.90

(12)

Gaseous phase: |E_ad| = 6.0 + 1.4E − 7.4S + 1.2A + 6.6B + 3.4L
n_tra = 33, R²_adj = 0.97, RMSE_tra = 0.62, Q²_LOO = 0.96, Q²_{kfold (k = 5, 5000)} = 0.96,
n_ext = 8, RMSE_ext = 0.61, Q²_ext = 0.97

(13)

Aqueous phase: |E_ad| = 1.5 + 3.5E − 3.8S + 0.89A + 1.7B + 13V
n_tra = 33, R²_adj = 0.96, RMSE_tra = 0.66, Q²_LOO = 0.95, Q²_{kfold (k = 5, 5000)} = 0.94,
n_ext = 8, RMSE_ext = 0.73, Q²_ext = 0.96

(14)

The |E_ad| or logK values predicted by the pp-LFER models agree well with the DFT or MD calculation values (Figure S6). Values of the statistical parameters (R²_adj, Q²_LOO, Q²_kfold and Q²_ext) indicate that these models possess high goodness-of-fit, robustness and predictive abilities [72]. As shown by the Williams plots (Figure S7), there are no outliers or high leverage points in the models. Thus, the pp-LFER models can offer an efficient way for high-throughput estimating the adsorption parameters (E_ad and logK) of other aliphatic and aromatic compounds structurally similar to the training set compounds.

Table 1. Relative contributions of intermolecular interactions to the overall adsorption of organic chemicals onto BP in an aqueous phase.

Term	|Ead|	logK	a logK		b logK
			TΔS	|Ead|	ΔGH2O	ΔGMD (Gaseous)
eE	23%	17%	−47%	62%	−5%	21%
sS	−26%	20%	93%	−70%	−8%	30%
aA	1%	0	−3%	3%	2%	−2%
bB	3%	−1%	−11%	9%	−1%	−1%
vV	88%	111%	−126%	239%	−49%	162%
c Sum	89%	147%	−94%	243%	−61%	210%

^a ΔG_MD(aqueous) was decomposed with Equation (4) into TΔS and E_ad using the results of MD and DFT, TΔS and E_ad were further decomposed into various intermolecular interactions by the pp-LFERs. ^b ΔG_MD(aqueous) was decomposed into ΔG_H2O (ΔG_H2O = ΔG_MD(aqueous) – ΔG_MD(gaseous)) and ΔG_MD(gaseous). ΔG_H2O and ΔG_MD(gaseous) were also decomposed into various intermolecular interactions by the pp-LFERs. ^c The sum consists of the five terms with the explicit interactions and does not include the contributions of the intercept term in the pp-LFERs. The contribution is average fractions of the terms in the sum of the six terms in the pp-LFERs (detailed in the Supplementary Materials).

and Table S16 list the relative contributions of intermolecular interactions to the overall adsorption. It can be seen that the sS term contributes negatively (−26%) to |E_ad|, and positively (20%) to logK. The contributions of the other terms to |E_ad| are not equal to those to logK either. It can be concluded that the relative contributions of a specific interaction to |E_ad| and logK are quite disparate. As can be seen from Equation (4), the disparity can be due to the absence of the entropic term (TΔS) in E_ad of the adsorption.

When the contribution of each pp-LFER term to logK is decomposed into the enthalpic (|E_ad|) and entropic (TΔS) components, the disparity can be explained explicitly. For example, the contribution of the sS term to logK is 20%, which is approximately equal to the sum of the enthalpic (−70%) and entropic (93%) contributions. Thus, the contribution of the sS term to logK is opposite to that to |E_ad|, since the negative enthalpic contribution is exceeded by the positive entropic contribution.

The total enthalpic contribution of the five terms is 243%, which is far greater than that of entropy (−94%). The result indicates that enthalpy dominates the adsorption, which echoes the significant correlations between logK and E_ad as indicated by Equations (9) and (10). Shen et al. [73] also found that the enthalpy dominates the adsorption of nitroaromatics on multiwall carbon nanotubes by experimental studies.

logK can also be decomposed into the contribution of water (ΔG_H2O) and of the interaction (ΔG_MD(gaseous)) between BP and the compounds in the gaseous phase. The total contribution of water is −61%, in which vV term (−49%) is the most significant factor. The vV term of ΔG_H2O may represent energy required to form a water cavity onto the BP surface to accommodate the adsorbates [33,67]. Thus, it can be concluded that the negative contribution of water to the adsorption mainly arises from the cavity effect.

Previous pp-LFER models on graphene, graphene oxides, carbon nanotubes were also adopted to analyze the contributions of the different interactions [17,19,20,29,74]. However, the previous studies did not distinguish the contributions of enthalpy, entropy and water to logK. In the current study, the contributions of entropy, enthalpy and water to the interactions were further estimated by the pp-LFERs integrated with the MD and DFT results.

3.5. Comparisons with Graphene

Table S17 lists the ΔG_MD and logK values for adsorption of 30 aromatic compounds onto BP and graphene in gaseous and aqueous phases. The logK values (Figure 4) for BP range from 1.9 to 12.5, which are close to those (1.6~12.4) for graphene. The ΔG_MD differences between BP and graphene range from −0.9 to 0.3 kcal/mol in gaseous and aqueous phases. The results indicate that the affinity of the aromatics to BP is comparable to that to graphene. Lazar et al. [28] also found that experimental ΔH of benzene adsorption onto BP was equal to that onto graphene.

The contributions of the different interactions to the overall adsorption for the 30 aromatics onto BP and graphene are shown in Figure S8. The contribution of each term to BP is close to that of graphene, explaining why the adsorption affinity of BP to aromatic hydrocarbons was comparable to that of graphene. The models in this study indicated that no matter for BP or graphene, the dispersion interactions (lL and vV) dominate the adsorption in the gaseous phase and the aqueous phase. Previous pp-LFER models for the adsorption of organics onto different carbon nanomaterials including graphene, carbon nanotubes and graphene oxides (

Table 2. Comparison of pp-LFER models developed in the current and previous studies.

No.	Nanomaterial	Phase	Ntrain	R2train	Prediction Model
1	SWCNT [34]	aqueous	30	0.87	logK = −1.3 + 0.40E + 0.36S + 0.93A − 3.9B + 2.8V
2	MWCNT [17]	aqueous	29	0.83	logK = −4.3 + 0.61S + 0.050A − 0.48B + 4.5V
3	MWCNT [75]	aqueous	28	0.93	logK = −1.3 + 0.043E + 1.7S − 0.37A − 2.7B + 4.1V
4	Graphene [20]	aqueous	29	0.89	logK = −1.4 + 0.11E + 1.4S + 0.42A − 3.8B + 2.2V
5	Graphene [29]	aqueous	35	0.88	logK = −1.8εα − 1.2εβ + 1.3q+ − 1.5q− + 1.0V − 1.6π + 42
6	Graphene oxide [20]	aqueous	36	0.84	logK = −1.4 + 0.29E + 0.28S − 0.19A − 2.6B + 2.6V
7	Graphene oxide [74]	aqueous	36	0.92	logK = −1.7 + 0.93E + 0.060S − 0.38A − 1.9B + 2.2V
8	BP (this study)	aqueous	33	0.87	logK = −1.7 + 0.65E + 0.75S + 0.048A – 0.095B + 4.0V
9	BP (this study)	gaseous	33	0.96	logK = 0.2 − 0.0076E − 1.1S − 1.1A + 1.4B + 1.5L
10	Graphene (this study)	aqueous	30	0.86	logK = −1.6 + 0.42E + 1.0S + 0.26A − 0.78B + 4.0V
11	Graphene (this study)	gaseous	30	0.97	logK = 0.052 − 0.043E − 0.62S − 0.78A − 0.36B + 1.5L

Multi-walled carbon nanotubes (MWCNT); single-walled carbon nanotubes (SWCNT); black phosphorus (BP); N_train: the number of compounds used for developing the model; R²_train: determination coefficient of training sets.

) also indicated that the dispersion interactions play crucial roles.

Although the number and category of the compounds in the reported models differ from those in this study, the pp-LFER models for predicting K values of BP and graphene (0.86~0.97) developed in the current study outperformed or was comparable to the previous ones (0.83~0.93) in terms of R²_train. Thus, the pp-LFER models in the current study not only are valuable for understanding the contributions of the different intermolecular interactions but also reliable in predicting logK values for estimating adsorption capacities.

4. Conclusions

In this study, DFT and MD calculations were combined with pp-LFER models to investigate the adsorption of aromatic chemicals onto BP and unveil differences in adsorption mechanisms between BP and the widely investigated graphene. The results indicate that the DFT and MD simulation can give reliable E_ad, logK and ΔG_MD values. The dispersion interactions dominate the adsorption of the aromatics in the gaseous phase and aqueous phase. As extra energy barriers need to be overcome to form a water cavity onto the BP surface to accommodate the adsorbates, the adsorption in the gaseous phase is more favorable than in the aqueous phase. The affinity of the aromatics to BP is comparable to that of graphene. The pp-LFER models can serve as effective tools to rank adsorption capacities and to estimate adsorption parameters of organics onto BP.