Adsorption of wastewater pollutants on amorphous TiO₂: an atomistic simulation study

Abstract

The increasing presence of polluting chemicals in man-made wastewater poses significant environmental and health risks. Advanced oxidation processes, particularly those involving photocatalytic materials like titanium dioxide (TiO₂), offer a promising solution for degrading these pollutants. This study employs force field molecular dynamics simulations to investigate the interactions between pollutants, the TiO₂ surface, water and ions, aiming to elucidate their role in the adsorption process. The results reveal that the protonation state of pollutants significantly influences their contact with the TiO₂ surface, with negatively charged species showing a higher affinity for the surface's active sites, especially those containing carboxylate groups. The formation of hydrogen bond networks affects the stability of these contacts positively, while the tendency of some pollutants to aggregate hinders surface contacts. Furthermore, we observe cations (Na⁺) to alter the surface-near environment in a typical electrical double-layer manner, as well as to participate in pollutant adsorption and aggregation. These findings provide insights into the adsorption features triggering the initial pollutant degradation on amorphous TiO₂, which could enhance the design of more efficient wastewater treatment technologies.

---

1 Introduction

The increasing levels of polluting chemicals in human wastewater are a growing concern as they can be harmful to us and our environment.^1–3 These pollutants originate from various sources and applications, vary in size and structure, and can be challenging to remove during water remediation processes.^1,4,5 Consequently, they can reach groundwater and rivers, interfering with the natural environment and endangering the balance of ecosystems.^2,6

Advanced oxidation processes (AOPs) are promising technologies for removing organic pollutants. In the focus of AOP research are photocatalytic materials such as TiO₂, which are excited under illumination, especially in the UV range. Energies greater than the band gap promote electrons from the valence band to the conduction band, creating electron–hole pairs. In aqueous solutions, electrons and holes can either directly reduce or oxidize pollutant molecules directly or form hydroxyl and superoxide radicals with surrounding water and oxygen, which also attack and degrade the pollutants.^3,5,7–9

The degradation of pollutants using AOPs and TiO₂ has already been achieved in experiments. Numerous studies have investigated degradation pathways, intermediates, and influencing factors, including pollutant species, solvent nature, pH value, available oxygen, illumination, photocatalyst type, concentration, and aggregation.^3,5,7–17 Additionally, the composition of wastewater, the characteristics of treatment plants, and the concentration of biomass can also alter the pollutant removal rate.¹⁸

Different TiO₂ phases like anatase or rutile possess different photocatalytic activities¹⁹ and pollutant degradation efficiencies.²⁰ Colloidal suspensions of nanoparticles are typically preferred for degradation processes over film-like structures due to their increased surface area,²¹ although they can aggregate under certain conditions, thereby reducing the available surface area and blocking light accessibility.²² The stability of a colloidal suspension can be rationalized by the Derjaguin–Landau–Verwey–Overbeek (DLVO) theory, where electrostatic double-layer repulsion and van der Waals attraction between the nanoparticle surfaces are considered. In wastewater, the surface of the nanoparticles is terminated by molecularly adsorbed or dissociated water molecules, resulting in terminal hydroxyl groups.^23–25 The charge distribution at the TiO₂ surface highly depends on the solvent and the pH. At pH values close to the point of zero charge (PZC, between 5 and 7 for rutile in water solutions at moderate ionic strengths), TiO₂ surfaces are amphoteric, presenting both net positive and net negative charged sites in a number and ratio that depends on the pH.²⁶ Counterions present in the solvent screen the charged particle surface, thereby lowering the electrostatic double-layer repulsion. As a result, coagulation of the particles is promoted above a critical counterion concentration.^27,28

The local environment on these charged surfaces can be described by the Graham Model of the electrical double-layer.²⁹ For titanium dioxide surfaces, the inner Helmholtz layer consists of highly oriented, rigid water molecules with embedded counterions that are in direct contact with the surface.^25,26,30 The degree of water orientation and the exact dimensions of the layer depend strongly on the crystalline facet. So far, studies have examined crystalline models indicating that the rutile (110) facet has a denser inner Helmholtz layer than all other surfaces.³¹

The outer Helmholtz layer is located at the position of a second shell of hydrated ions.²⁹ Theoretical studies on titanium dioxide explore the surface's impact in terms of the orientation of the water layer relative to the normal. For all TiO₂ facets, the surface impact range varies within 5–8 Å depending on the surface charge density.^26,31–33 These length scales can be interpreted as 2–3 layers of oriented water and fall within the range of the outer Helmholtz layer. The highest degree of orientation occurs on the rutile (110) facet,³¹ as the first water layer forms through subsequent surface titanium and water oxygen contacts.³⁴

The diffuse layer thickness is determined on the one hand by the surface charge density and on the other hand directly depends on the ionic strength of the solvent. Under physiological conditions, the thickness is estimated to be about 8 Å,^26,33 completely embedded in the range of structured water layers. At lower salt concentrations, the thickness of the diffuse layer increases significantly up to a range beyond the local structural features.²⁶ Therefore, the electric double-layer affects molecular adsorption processes, including the initial approach of molecules to the surface and the formation of localized contacts.

Not only does the surface change with pH,²⁵ but the pollutants also vary. Many wastewater pollutants are small organic molecules with neutral or negatively charged functional groups at pH levels between 5 and 8. Their charge states can influence adsorption and degradation in several ways. On the one hand, the net charges on the surface can promote or hinder the adsorption of charged pollutants, affecting direct degradation pathways. On the other hand, the indirect degradation through water splitting occurs via the formation of hydroxyl or superoxide radicals under either basic or acidic conditions respectively.⁹

This makes predicting how pH affects the degradation rate far from trivial. Many studies report a decrease in degradation efficiency at pH values above the pollutant's pK_a, i.e. when the pollutants are in an anionic form.^35–37 While this is often explained by the lower adsorption tendency of negatively charged pollutants on an overall (slightly) negatively charged surface, we note that the amphoteric character of TiO₂ and the presence of dissolved cations may still allow for direct adsorption at locally positive surface terminal groups. Furthermore, other forces besides electrostatics, particularly a match between the hydrophilic/hydrophobic patterns of the surface and the pollutant at the molecular level, are known to influence the binding affinity of organic molecules to the surface.^38,39 In fact, some studies indicate a higher degradation efficiency for pollutants at pH values above their pK_a.¹⁶ Under these conditions, the charge delocalization over ring systems of aromatic pollutants can also attract radicals,^2,40 favoring an indirect degradation pathway.

Understanding model systems, such as those studied in laboratories for decades, becomes more complex when considering remediation with AOPs due to the wide variety in wastewater's chemical composition and potential degradation pathways. This leads to many uncertainties that need clarification of the underlying molecular mechanisms. Clarifying these mechanisms would help optimize efficiency, identify reaction pathways and degradation intermediates, and support practical applications. Typically, assessments of specific use cases involve extensive experimental testing and numerical simulations to understand the mechanisms at play.

In this work, we investigate the complex interplay of degradation process participants, including pollutants, the TiO₂ surface, water, and ions, using force field molecular dynamics (FF MD) simulations. The purpose is to clarify the influence of parameters such as surface and pollutant structure, ions, protonation state, and solvent, and to establish guidelines for selecting effective parameter combinations for other applications. We focus on the prerequisite for subsequent chemical degradation: the formation of stable (i.e., lasting long enough to enable redox reactions) direct or indirect contacts between pollutants and surface active sites. Eleven common pollutant substances from various applications are selected for this study (Table 1), representing a subset of frequently encountered wastewater pollutant classes. These are detailed in the following section.

Table 1 Modelled wastewater pollutants. In this paper the pollutants are mentioned mainly by their abbreviation code. To give an overview of the compound characteristics, the typical application field, their pK_a, employed structure, appearance in different water bodies (c_nat) and guideline values of limit concentrations (c_lim) in aqueous environments are summarized here. Beneath calculated dipole moments, a dipole moment from literature is given for comparison. SW – surface water, GW – ground water, WI – wastewater treatment plant (WWTP) influent, WE – WWTP effluent, CPE – chemical plant effluent, ND – non-detectable

Pollutant	Abbreviation	pKa	Structure	Dipole (D)	c nat (µg L^−1)	c lim (µg L^−1)
a Dipole moment compared at different protonation states, details in SI A.6. b Uslu et al.:^65 PNEC – lowest predicted no effect concentrations for invertebrates, algae, fish, plants and microorganisms are reported. c Schulman et al.:^44 AWQC – ambient water quality criteria for human based on health based limits for non-cancer effects, risk-specific doses for linear carcinogens, drinking intake, fish intake and bioaccumulation factor. d Pais and Nascimento:^66 guideline value for human based on tolerable daily intake, body weight, fraction in drinking water, drinking intake and NOAEL (no-observed adverse effect level). e Brotzmann et al.:^67 NOEC – no observed effect concentration in 120 h zebrafish embryo study (acute fish embryo toxicity test). f Wang et al.:^68 EU drinking water limit. g Feltracco et al.:^53 max. annual mean surface water concentration set by the Italian Legislative for each pesticide. h WHO:^69 WHO guidelines for drinking water quality.
Pharmaceutics
Acetaminophen	ACE^0	9.4^41		2.2	2.3 [SW]^41	9200^b
				2.1–2.4 ^42	62.0 [WE]^43
Acetylsalicylic acid	ASA^−1	3.5^44		7.1	1.5 [WE]^45	480^c
				—^a	0.4 [SW]^46
Carbamazepine	CBZ^0	2.3^47		3.6	4.6 [WE]^43	500^b
		13.9^47		3.6 ^47	1.7 [SW]^48
					0.2 [GW]^18
Ibuprofen	IBU^−1	4.9^41		13.5	14.6 [WE]^43	2800^d
				—^a	5.0 [SW]^41	3840^b
Valproic acid	VLP^−1	4.5^49		6.6	0.1 [WI]^50	2080^e
				—^a	ND [WE]^50
Agriculture
Glyphosate	GGG^−2	2.0^51		10.3	0.2 [SW]^53	0.1^f^,^g
		2.6^51		10.9 ^52	441.5 [SW]^54
		5.6^51
		10.6^51
Tetracycline	TTR^−1	3.3^55		13.1	3.6 [WE]^43	90^b
		7.7^55		—^a	0.05 [SW]^56
		9.7^55
	TTR^−2			19.3
				—^a
Every-day life
Caffeine	CAF^0	0.7^57		3.8	51.7 [WE]^43	87 000^b
		14.0^57		3.6–3.7 ^58	1.7 [SW]^41
					0.7 [GW]^18
Salicylic acid	SLC^−1	2.9^59		9.3	10.1 [WE]^43	190^c
		13.6^60		—^a	0.3 [SW]^41	134 000^b
Industry
Benzene	BEN^0	—		0.0	179.0 [CPE]^62	10^h
				0.0 ^61	80.0 [GW]^62
Dichloromethane	DCL^0	—		1.7	743.0 [SW]^64	20^h
				1.6 ^63	3600.0 [GW]^64

---

2 Investigated pollutant substances

We choose substances from diverse application fields: pharmaceuticals, agriculture, every-day life products, and industrial waste (Table 1).

Pharmaceuticals are frequently detected in our waste, ground, and surface waters, including acetaminophen (ACE),^8,43,46 known as paracetamol, and acetylsalicylic acid (ASA),⁴⁶ known as aspirin, as well as ibuprofen (IBU)^18,56 and carbamazepine (CBZ).^18,46,56 The analgesics ACE, ASA, and IBU are widely used to treat pain, fever, and inflammation^5,8,70 and are readily available worldwide. Carbamazepine and valproic acid (VLP) are anticonvulsants used to treat epileptic seizures.^13,67,70 Although VLP is found only in low amounts in wastewater influents, it is commonly prescribed and has a high toxicity potential.⁵⁰ Therefore, its behavior in AOPs is of interest.¹⁷ The antibiotic tetracycline (TTR) is used to treat bacterial infections in human medicine as well as in veterinary medicine for poultry, husbandry, and fishery.^71–73 Since it is one of the most widely used drugs, high concentrations of this substance can be found in wastewater treatment plants (WWTPs) serving pharmaceutical factories.⁷³ Its high excretion rate (up to 80%⁷²), aqueous solubility and environmental stability⁷¹ are causing growing concern.

A polluting substance widely used in agriculture is the broad-spectrum herbicide glyphosate (GGG).^74,75 By inhibiting the plant enzyme EPSP, it damages plants that compete with crops. The impact of GGG on animal and human cells is still under investigation, but its potential to disrupt the endocrine system in human cells was demonstrated by Gasnier et al.⁷⁶ in 2009. GGG can disperse through various means, including wind, rain, and water erosion. It has low mobility in soil,⁷⁵ but can be found in surface and groundwater samples at concentrations exceeding the EU drinking water limit.⁶⁸ The initial degradation of GGG can be relatively quick; however, reaction intermediates like AMPA (aminomethylphosphonic acid) tend to be more persistent than GGG itself and are highly toxic.⁷⁵

We also include in our study every-day life products that can enter the environment through excretion and wash water.⁷⁷ We consider the stimulant caffeine (CAF), released from tea and coffee consumption,^58,77,78 and salicylic acid (SLC), an additive in cosmetics and food,¹¹ as well as a precursor and metabolite of ASA.^11,44 Both can be found in significant concentrations in wastewater near urban areas.¹¹ Benzene (BEN) and dichloromethane (DCL) are used as organic solvents and processing materials in the chemical industry.^{15,62,64,78,79} They are released into the environment through production emissions, effluents, and the use of products containing these chemicals.^69,80 While only suggestive evidence of carcinogenicity is found for DCL, BEN is confirmed to be carcinogenic to humans.⁶⁹

This selection of wastewater pollutants will ensure a comprehensive screening for common impurities in water bodies. TiO₂ has already demonstrated the potential to degrade most of these pollutants photocatalytically under controlled laboratory conditions.^{3,5,7,9–17,81}

3 Methods

3.1 Pollutant models

All pollutant molecules (Table 1) considered here are organic molecules with a mass less than 444.4 g mol⁻¹ (ref. 82) (TTR⁰). Thus, their chemical groups resemble those found in biological macromolecules, which have been studied with FF MD for decades. The force field parameters are typically taken from AMBER99SB and GAFF,^83,84 with adjustments made to the partial charges. The partial charges are calculated using the DFT software ORCA v4.2.1.⁸⁵ First, the molecular geometry is optimized in vacuum using the B3LYP hybrid functional and the cc-pVTZ basis set. In the next step, the grid-based CHELPG charges are computed from the electrostatic potential (ESP) with the same functional and basis set. These ESP charges are then refined to ensure equal charges for equivalent atoms, based on the restrained electrostatic potential (RESP) approach. A table listing all RESP charges used is provided in the SI A.1.

The geometry optimizations and CHELPG calculations of GGG⁻² and TTR⁻² were performed using implicit CPCM water (conductor-like polarizable continuum model) to prevent spurious proton transfers to the charged group that can occur during optimization in vacuum. For TTR⁻², an additional distance restraint was included for the bond between O and H3 (atom naming scheme: SI A.1.6) in the geometry optimization. This approach stabilized the structures of TTR⁻² and GGG⁻² reported by Zuo et al.⁸⁶ and Sprankle et al.,⁸⁷ respectively (cf. SI A.2). Since TTR presents a pK_a value of approximately 7.7,⁵⁵ both TTR⁻¹ and TTR⁻² species are present in solution at pH 7.4, so both protonation states are included in our study.

3.2 TiO₂ model

An adjusted amorphous TiO₂ model is used, originating from Derr et al.⁸⁸ It consists of a rutile crystal cut at a high Miller index, relaxed, and annealed at 500 K, which allows for the recombination of surface atoms in a thin amorphous surface region, as commonly present in TiO₂ nanoparticles used in photocatalytic degradation experiments. The intramolecular FF parameters for TiO₂ are taken from Schneider and Colombi Ciacchi.⁸⁹ The intermolecular parameters for interactions between TiO₂, water, and biomolecules were derived in Schneider and Colombi Ciacchi⁹⁰ and are used here with the net-charge corrections from Derr et al.⁸⁸ This parameter set has been applied successfully in previous studies by Laube et al.⁹¹ and Balzaretti et al.⁹²

The surface in contact with TIP3P⁹³ water is terminated according to the following rules. (1) Ti atoms, that are less than 5-fold coordinated by O atoms after the annealing step, are terminated with one hydroxyl group (TNX/OTX/HTX residues). An equal number of H atoms is placed on surface O atoms, forming bridging hydroxyl groups (ONX/HBX residues). In this way, TIP3P water can adsorb in molecular form on the less than 6-fold coordinated Ti atoms during subsequent MD simulations, consistent with the findings of Schneider et al.^89,90 (2) A small amount of dissociative adsorption at these sites is accounted for by terminating some of the 5-fold coordinated Ti atoms with additional hydroxyl groups (HTI) and an equal number of surface O atoms with H atoms (HTT). The number of these sites is adjusted based on a pH-dependent density of surface charged groups. They are taken from Köppen and Langel²⁶ who used the modified multisite complexation model (MUSIC) for crystalline rutile surface models. The partial charges and FF parameters of these OH and H terminal groups are the same as those of TIP3P water. (3) Net positive and negative sites (formed by the depletion of hydroxyl groups and protons, respectively, as explained below) are then distributed on the surface to model its amphoteric character according to the MUSIC model. The total density of charged sites is set at 2.4 nm⁻² at the PZC of 6.0,²⁶ plus an excess of negative charges amounting to −0.055 C m⁻² at the simulated pH of 7.4.⁸⁸ To model the formation of integer +1 and −1 charges, respectively, upon removal of HTI hydroxyl groups and HTT protons (which carry partial charges equivalent to TIP3P water), the partial charges of the first and second neighboring Ti or O atoms of these groups are adjusted as shown in Fig. 1, resulting in positive and negative TI2 and OTT residues, respectively. The total amount and surface density of all these surface-terminating residues is reported in Table 2. All partial charges, FF parameters, and a snapshot of the TiO₂ surface model are reported in SI A.3.

Fig. 1 (A) Overview of the simulation cell with the TiO₂ surface (gray), pollutant molecules (here CAF⁰), Na⁺ (orange), and water (transparent background within blue boundaries). (B) Active sites at the TiO₂ surface with atom names. Hydroxylated undercoordinated titanium surface atoms form HTI groups with OH⁻ from water dissociation. In the corresponding dehydroxilated form (TI2), the charges on the first and the second neighbor (Ti_t, O_b, and O_s) are adjusted (Δq_Ti, and Δq_O in units of the electron charge, e). (C) The water's H⁺ can protonate bridging oxygens and create an HTT group. Its deprotonated group (OTT) has changed charges.

Table 2 Number and density of TiO₂ sites for the symmetric amorphous slab model. The numbers of active site types (HTI/TI2 and HTT/OTT) are adjusted to a pH value of 7.4 according to Köppen and Langel.²⁶ Since the surface is amorphous, less than five-fold coordinated surface titanium atoms are saturated with TIP3P hydroxyls (TNX), resembling HTI groups in structure. The remaining proton is attached to an undercoordinated surface oxygen (ONX), resembling HTT. TNX and ONX are called inactive sites since they are not adjusted according to the pH. The distribution of sites results in a surface charge density of −0.055 C m⁻² (ref. 88)

Site	Number	Site density (nm^−2)
HTI	120	0.768
TI2	254	1.626
HTT	66	0.422
OTT	308	1.971
TNX	234	1.498
ONX	234	1.498

---

In all simulations, an inversion-symmetric surface slab model is used to prevent artificial macroscopic dipole moments in the simulation cell. To achieve this, the pristine surface model of a single TiO₂/water interface is duplicated, rotated by 180° around the x-axis, and shifted along the z-axis perpendicular to the interface as needed to create a continuous slab, shown in Fig. 1A. The final slab measures 8.16 nm × 9.53 nm × 4.02 nm along the x-, y- and z-axes, respectively. Since the titanium atoms Ti_sat and Ti_t have similar charge ranges and locations, they are grouped under Ti_t in the analysis to simplify the presentation of results. Similarly, O_sat and O_s are combined under O_s.

3.3 Combined models and simulation setup

Our simulations incorporate four components: pollutants, the surface model, ions, and water (Fig. 1A). These are assembled in the input files using tLeap from AmberTools v12 and v18,⁹⁴ then converted from AMBER to GROMACS input format with the acpype.py script.⁹⁵ All FF MD simulations are performed using GROMACS v2022.6.^96,97 In each setup, 25 molecules of each pollutant species are placed above the TiO₂ surface to ensure thorough sampling of surface/molecule interactions and adsorption modes. This results in a concentration of ∼50 mM. At the beginning of each simulation, the molecules are arranged in rows of 5 × 5, with their centers of geometry separated by 12 Å (15 Å for TTR) and at ∼10 Å from the surface. To prevent the results from depending on the initial orientation, each molecule in this arrangement is rotated by about 15° (360°/N_molecules) around the x-axis relative to the previous molecule before starting the MD simulations.

The models are solvated with TIP3P water in a periodically repeated orthorhombic box measuring 8.16 nm × 9.53 nm × 15.81 nm, with the surface slab aligned in the xy plane and sufficient water bulk phase between the slab and its periodic images in the z direction. The systems are neutralized with sodium ions, resulting in a number of ions that depend on the pollutant (POL) protonation state: POL⁰ models include 54 (neutralizing the surface charge only), POL⁻¹ 79 (neutralizing surface charge and 25 single negatively charged pollutants), and POL⁻² 104 sodium ions (neutralizing surface charge and 25 double negatively charged pollutants), corresponding to concentrations of 0.10 M, 0.14 M and 0.19 M, respectively. A simulation without pollutants is run as a reference system. In the absence of POL, the setup includes the TiO₂ surface slab, water molecules, 108 sodium ions (0.20 M), and 54 chloride ions (0.10 M). In this case, a symmetric ion distribution on both sides of the surface slab is ensured by bias wall constraints placed in the middle of the water region, which act selectively on the ions with a maximum force constant of 1.500 kJ mol⁻¹. These constraints are implemented using the PLUMED package v2.8.3 and v2.9.0⁹⁸ (see SI A.4).

A series of conscientious equilibration steps is carried out prior to the production MD runs. An initial energy minimization is followed by equilibration in the NVT and then the NpT ensemble with position restraints applied to the pollutants (1000 kJ mol⁻¹ nm⁻²) and the surface core atoms (10 000 kJ mol⁻¹ nm⁻²). The surface terminal groups are not constrained and thus free to move. The pollutants' position restraints are then removed, and another series of energy minimization, NVT, and NpT equilibration runs is performed. Finally, the system is annealed through a short MD run that increases the temperature from 300 K to 500 K over 300 ps, maintains it at 500 K for 400 ps, and then decreases it back to 300 K over 300 ps. The equilibration protocol is detailed in the SI A.5. Production runs are performed in the NVT ensemble at 300 K and 1 bar, each lasting 100 ns. All simulations are run under periodic boundary conditions applied in the x, y, and z directions, with van der Waals and real-space electrostatic cut-offs set at 1.2 nm. The energy is minimized using the steepest descent algorithm. We employ the LINCS algorithm to constrain hydrogen bonds, the Verlet cut-off scheme, and the particle mesh Ewald (PME) method for long-range electrostatics. We use the V-rescale thermostat in the NVT ensemble and an additional semi-isotropic Parrinello–Rahman barostat in the NpT ensemble. To prevent spurious deformation effects in the xy plane of the TiO₂ slab, the compressibilities are set to 4.5 × 10⁻⁹ bar⁻¹ along the x and y directions and to 4.5 × 10⁻⁵ bar⁻¹ along the z direction.

We visualize our systems using VMD v1.9.3⁹⁹ and analyze the trajectories with GROMACS and Python3, mainly using the packages MDAnalysis v2.3.0,^100,101 NumPy v1.23.3, pandas v1.5.1 and Matplotlib v3.6.2. Molecular electrostatic potentials on solvent-excluding surfaces are calculated with Jmol v16.1.47.¹⁰²

4 Interactions

4.1 TiO₂ in aqueous solution

First, the interactions between water, the TiO₂ surface, and ions are studied in the absence of pollutant molecules. Due to the amorphous character of the TiO₂ surface, its density profile exhibits an irregular and somewhat rough decay along the z axis compared to crystalline surfaces (Fig. 2A). In the analysis protocols, we arbitrarily set the surface plane (0 Å distance) where the normalized TiO₂ surface density equals 0.5 relative to its average bulk value. Water molecules occupy the gaps between the terminal groups on the surface and penetrate deeper into the structure than Na⁺ and Cl⁻ ions. The density profile of water relative to the nearest surface atoms (Fig. 2B) displays two maxima for the O (at 2.0 Å and 2.7 Å) and H (at 1.8 Å and 2.7 Å) atoms, respectively. The closer proximity of H atoms to the surface compared to O atoms indicates that a substantial amount of water forms hydrogen bonds with the surface. The variation in water density relative to the bulk extends up to ∼6 Å away from the surface. These results qualitatively align with the early simulations by Köppen and Langel²⁶ of the crystalline rutile (100) surface, although in that case water layering induced by the surface extended up to 9 Å. The higher range of oriented water is anticipated due to the crystalline structure of the surface. Our findings also agree with Bischoff et al.,³⁰ who observed that the influence of charged amorphous TiO₂ surfaces reaches into the water phase up to three molecular layers at salt concentrations above 100 µM. Consequently, we defined our criterion for direct contact as a relative distance of 2.5 Å or less from the surface, corresponding to entities within the first water layer and inner Helmholtz layer. The criterion for indirect contact is a relative distance greater than 2.5 Å up to 6 Å to include molecules within the second and third water layers.

Fig. 2 Density profiles in simulations without pollutants. (A) Water, ion, and TiO₂ density profiles along the z-direction with z = 0 set where the TiO₂ density reaches a normalized density of 0.5. The TiO₂ density profile reflects the roughness of the amorphous surface. Water and Na⁺ ions approach and can penetrate the surface. Cl⁻ ions are more separated. TiO₂ and water densities are normalized with respect to their bulk density, ion densities to their maximum value. (B) Density profiles for ions and water's oxygen and hydrogen atoms as a function of the distance to the closest TiO₂ atom. The double peaks in the OW and HW profiles emphasize the water layers near the surface. All profiles are normalized w.r.t. their maximum value. (C) Contact frequencies of the ions with the atoms of the surface as defined in Fig. 1. Note that the contact frequencies with Ti_t and Ti_b are zero, while Cl⁻ are rarely found near hydrogen atoms (H_t, H_s), but Na⁺ quite frequently near oxygen atoms (O_t, O_s, O_b). In all the panels, quantities with index [0] or [POL] refer to the simulations without or with pollutants, respectively.

In the absolute density profiles (Fig. 2A) Na⁺ ions are abundant near the surface and in its grooves, peaking first at approximately −0.1 Å and then forming multiple broader peaks up to ∼2.8 Å from the surface. On the other hand, the relative density (Fig. 2B) shows a strong Na⁺ peak at 2.3 Å, followed by a faint second peak around 4.3 Å. They form close contacts with O atoms on the surface and of the terminal hydroxyl groups (O_s, O_b, O_t species, Fig. 2C). This suggests that Na⁺ ions are located in the closest layer to the surface and the uneven surface topology causes the different maximum heights in the density profile over the cell height.

Cl⁻ ions gather near the surface, but at a greater distance, with an initial density peak at ∼3.2 Å and a broader distribution (Fig. 2A). A similar pattern appears when examining the profiles relative to the nearest surface atoms (Fig. 2B), where Cl⁻ ions show two peaks at 2.8 Å and 3.8 Å. The ion adsorption fractions (Fig. 2C) show that chloride ions have almost no direct contact with the surface, instead residing in the outer Helmholtz layer. In the few cases where they approach the surface, they attach to H_s and H_t atoms that are part of the hydroxyl groups (HTT and HTI residues, both 0.2%). Considering the position of the density minima of both water and ions at about 3.3 Å as the edge of the outer Helmholtz plane, we see that the vast majority of Na⁺ ions are located in the inner Helmholtz plane, while the first peak of Cl⁻ ions is inside the outer Helmholtz plane, whereas the second peak of Cl⁻ ions accumulates nearby as part of the diffusion layer. An enrichment of ions forming inner-sphere complexes in direct contact with the TiO₂ surface is known from previous works,^{26,30,103,104} as well as the Cl⁻ enrichment in the diffuse layer.^26,104 Notably, ions in direct surface contact interacting with the H atoms of terminal and bridging hydroxyl groups have also been observed in simulations of neutral rutile (110) surfaces.¹⁰⁵

In summary, the surface presents a diverse array of terminal groups capable of engaging in ion contacts, particularly with positively charged species, resulting in a highly structured electrical double-layer. Furthermore, negatively charged species can closely approach the inner Helmholtz layer while remaining fully solvated, despite the overall negative charge of the surface at pH levels above the PZC. We therefore anticipate various forces to influence the adsorption of pollutant molecules, which will need to compete with adsorbed ions and water to either penetrate the inner Helmholtz layer and make direct surface contact or interact more loosely with the surface as outer Helmholtz complexes. Before addressing these issues, we will investigate the POL molecules in bulk water with salt ions in the next section.

4.2 Pollutants in aqueous solution

The selected pollutants are small organic molecules with 5 to 56 atoms and protonation states ranging from 0 to −2. When solvated in water, they display neutral and charged regions that interact with water molecules and ions. These regions can be visualized by mapping their molecular electrostatic potential (MEP) onto their solvent-accessible surface area (SASA, Fig. 3A). As expected, the neutral pollutants (ACE⁰, BEN⁰, CAF⁰, and CBZ⁰) show a surface potential close to zero. In contrast, the deprotonated pollutants (ASA⁻¹, IBU⁻¹, SLC⁻¹, TTR⁻¹, VLP⁻¹, GGG⁻² and TTR⁻²) have negative surface potentials not only around the charged functional groups, but also distributed across the whole molecule, reflecting the presence of delocalized electrons in the aromatic molecules. The only positive surface potential occurs in DCL⁰ around its hydrogen atoms, compensating the negative partial charged localized on the Cl atoms.

Fig. 3 (A) Pollutants' molecular electrostatic potential (MEP) mapped onto the solvent accessible surface area (SASA). The color ranges from −0.81 eV (red) to +0.81 eV (blue). The SASA area is displayed beneath the molecule in nm². (B) Sodium attachment sites at the pollutant functional groups. Since only positive Na⁺ ions are present in the systems with pollutants, they get closest to the highly negative functional groups (i.e. carboxylate groups, phosphonate, ketone, hydroxyl, carboxamid). (C) RDF of sodium ions at the functional groups of GGG⁻² and TTR⁻².

The calculated dipole moments for the neutral species range from 0.0 D (BEN⁰) to 3.8 D (CAF⁰) (Table 1). The dipole moments of the negatively charged pollutants are larger, ranging from 6.6 D (VLP⁻¹) to 19.3 D (TTR⁻²). These calculated dipole moments fall within the ranges reported in the literature (cf.Table 1 and SI A.6.1).

Fig. 3B shows the direct contact fractions for Na⁺ ions with the pollutants' functional groups. As described in Section 4.1, direct contacts are defined when ions approach molecules at distances less than 2.5 Å. As expected, Na⁺ primarily interacts with the functional groups of negatively charged pollutants that have stronger localized charges. For example, GGG⁻² has a maximum negative MEP of −0.808 eV, mainly located at its phosphonate and carboxylate groups, while TTR⁻² has a maximum of −0.723 eV, mainly located at its ketone, hydroxyl and carboxamide groups. All these groups bind to Na⁺ ions with high frequency, particularly in the case of well-exposed anionic groups (such as phosphonate groups in GGG⁻²) and to a lesser extent for anionic groups surrounded by positive or neutral partial charges (such as ketone groups in TTR⁻²).

Consistent with the observation above, the RDF profiles of Na⁺ around oxygen containing anionic groups (Fig. 3C) of GGG⁻² and TTR⁻² reveal very strong peaks at phosphonate, at distances of around 2.32 Å for the first maximum, followed by two additional but smaller maxima at 3.68 Å and 4.08 Å. Similar trends with smaller maxima are observed for the carboxylate groups (first maximum at 2.36 Å), the ketone groups (first maximum at 2.34 Å), the carboxamide groups (first maximum at 2.32 Å) and hydroxyl groups (first maximum at 2.4 Å, second maximum at 3.26 Å). The first neighbor signals in the RDFs of all anionic groups surrounding sodium ions reveal a distance close to the typical ion-water oxygen distance of 2.2 Ź⁰⁶ for fully hydrated ions. Therefore, oxygen atoms of the functional groups replace at least one water molecule in the sodium hydration shell to form permanent direct contacts. For further investigation, the ions must be considered, as they can form stable neutralized complexes with anionic pollutants and thus directly influence the adsorption properties.

4.3 Pollutants on TiO₂

4.3.1 Contact frequencies. To study the parameters of pollutants interacting with the TiO₂ surface, the contact behavior with surface atoms is analyzed. We have identified two nominal categories of pollutants based on their observed behavior in the MD simulations: (i) forming surface contact and (ii) clustering. The categorization depends on whether the pollutant molecules can form direct contact interactions with the surface or with each other, respectively. The contact criteria are consistent with those set in Section 4.1.

The surface contact frequency for all pollutants is shown in Fig. 4, measured as the percentage of time during which direct or indirect contacts occur over the 100 ns trajectories, averaged across all 25 pollutant molecules in each simulation. They are derived from more detailed contact maps over time, which can be found in the SI Section A.7.

Fig. 4 Pollutant surface contact percentage over the whole trajectory length (100 ns). (A) Neutrally charged pollutants show the lowest direct contact behavior (< 3.3%) with CBZ⁰ showing nearly no contacts at all. Indirect contact levels range from 1.1–7.0%. (B) Single negatively charged species show a slightly higher contact rate (3.6–16.4%) with also higher rates of indirectly adsorbed molecules (6.9–10.6%) than neutrally charged pollutants. (C) The highest probability for direct interaction is observed for twice negatively charged pollutants (GGG⁻² with 27.8%, TTR⁻² with 47.8%), whereas the indirect contact rates are 37.2% for GGG⁻² and 9.6% for TTR⁻². (D) Surface atom participation in direct contact events. Especially hydrogen atoms (H_s, H_t) generate a high contact probability. (E) Adsorption probabilities of the pollutants functional groups associated with direct contacts. Negatively charged groups like the phosphonate, carboxylate, carboxamide, hydroxyl and ketone show high contact rates. (F) Adsorption probabilities of functional groups at surface atoms. Next to the hydrogen, bulk and terminal oxygens, only terminal titanium atoms are involved in minimal distances with carboxylate and carboxamide groups. The functional groups are namely amino (–NH₂), benzol (–C₆H_x), carbohydrate (–C_xH_x), carbonylic acid (–OCO), carboxamide (–CONH₂), carboxylate (–COO⁻), chloride (–Cl), hydroxyl (–OH), keton (–CO), methyl (–CH₃), nitrogen (–N), peptide (–CONH), phosphonate (–PO₃²⁻) and purin (–C₅N₄H).

Neutral pollutants (Fig. 4A) show the lowest direct contact frequencies, ranging from 0.2% for CBZ⁰ and 0.3% for CAF⁰ to 3.3% for ACE⁰. For these molecules, the indirect contact frequency is also very low for CBZ⁰ (1.1%), but reaches the range of 4.2% to 7.0% for CAF⁰, DCL⁰, BEN⁰ and ACE⁰.

A slightly higher frequency of indirect contact is observed for all monovalent negatively charged pollutants, ranging from 6.9% (ASA⁻¹) to 10.6% (TTR⁻¹). Notably, the indirect contact frequency is much higher than the direct contact frequencies for neutral molecules, and slightly higher for singly negatively charged molecules, with the exception of TTR⁻¹, which shows a larger direct contact frequency (16.4%) than indirect. TTR⁻¹ has the highest direct contact frequency of all monovalent anionic pollutants. Strikingly different and much higher are the direct contact frequencies of both pollutants with a net charge of −2, GGG⁻² and TTR⁻², with 27.8% and 47.8%, respectively. The indirect contact frequencies of the divalent anionic pollutants are also noteworthy: while GGG⁻² exhibits even more molecules in indirect contact (37.2%) than in direct, only 9.6% of TTR⁻² are in indirect contact with the surface.

In experimental studies, degradation is often investigated, rather than adsorption. Even if so, most studies do not report their solution pH or the amount of salts used to tune the pH, which both alter the interface environment. Ma et al.¹⁰⁷ reported removal efficiencies of 42.8% to 58.9% from 30 to 150 minutes of equilibration time for TTR on TiO₂, which is comparable to the direct contact frequency of 47.8% observed here for TTR⁻². On the contrary Yang et al.¹⁰⁸ only found 13.3% of adsorbed TTR. Neither studies report the pH used in their experiments, even though this can greatly influence the amount of adsorption, as shown by Chen and Liu:⁹ they determined adsorption amounts of glyphosate onto TiO₂ ranging from 3.5 × 10⁻⁶ mol g⁻¹ at a pH of 10 to 9.8 × 10⁻⁶ mol g⁻¹ at a pH of 2. At pH values of 6 to 8, the adsorption amount is around 6 × 10⁻⁶ mol g⁻¹, resulting in around 24% of adsorbed molecules, which is in line with our findings. Bianchi et al.¹⁰ used the same setup to compare adsorption and degradation of ACE and ASA, and found negligible amounts of adsorption of ACE⁰ and ASA⁻¹ (1.3% and 2.6%, respectively) on anatase TiO₂ powder in dark conditions in deionized water at pH 6.6. Although the percentages in our simulations are higher, the same trends are observed: Both species adsorb onto the surface in small amounts, with ASA⁻¹ to a greater extent than ACE⁰.

4.3.2 Binding motives. The frequency of direct contact is further broken down into all-atom contact maps of surface atoms and adsorbed pollutant species (Fig. 4D), of pollutant functional groups (Fig. 4E) and of functional groups in contact with surface atoms (Fig. 4F). The atom-resolved distance analysis shows a contact signal as soon as any atom of the addressed functional group is in contact with any atom of the surface.

At the surface, H atoms of terminal (H_s) and bridging (H_t) hydroxyl groups form many contacts, especially with GGG⁻² and TTR⁻² (about 8%–16%), but also with monovalent anionic pollutants (between 2% and 5%) and to a lesser extent with neutral pollutants. Even though H_s and H_t carry the same partial charges, the terminal OH groups (H_t) more frequently act as contact points than the bridging hydroxyl groups (H_s). This can be attributed to the slightly different occurrences of the two groups in the TiO₂ model (H_s : H_t = 300 : 354), but possibly also because of the better solvent exposure of the H_t atoms. Instead, the negatively charged O atoms of bridging hydroxyl groups (O_s) are rarely involved in contacts except with GGG⁻², TTR⁻¹ and TTR⁻². In general, the more accessible O atoms of the TiO₂ surface (O_b) and terminal O atoms (O_t) interact extensively with several pollutants, in particular ACE⁰, GGG⁻², TTR⁻¹ and TTR⁻². Ti_t surface atoms are only involved in contacts with GGG⁻² and TTR⁻². In the case of TTR⁻², contacts are formed via the carboxamide group. In the case of GGG⁻², contacts to Ti_t atoms are formed only via the carboxylate group (Fig. 4F). This is remarkable in view of the recent results reported in Balzaretti et al.,⁹² where we found that GGG⁻² was able to form bonds with partially covalent character towards the surface Ti atoms of rutile (110) via both the carboxylate and the phosphonate groups. The here-observed absence of phosphonate-mediated contacts can be attributed either to the difference surface topography (amorphous rather than crystalline) or to molecular competition, as discussed in Section 4.3.4.

In general, the involvement of phosphonate (present in GGG⁻²) and carboxylate groups (present in five pollutants) occurs with frequencies as high as 10%. All functional groups of the pollutants involved in contacts and their corresponding frequencies are reported in Fig. 4E (see SI Section A.1 for a detailed description of the pollutants functional groups). In TTR⁻¹ and TTR⁻², carboxamide, ketone, and hydroxyl groups also show pronounced involvement, with frequency peaks at 16%, 12% and 17%, respectively. Interestingly, also (to a lesser extent) ACE⁰ exhibits the most contact with the surface through its hydroxyl group. The high contact frequency for the amino group of GGG⁻² points towards findings from a previous study,⁹² where we investigated the GGG⁻² adsorption orientation on a pristine rutile (110) surface, finding it mostly in bidentated conformations. Through the negatively charged ends with a positively charged amino group in the middle, GGG⁻² was able to follow local surface charge patterns, resulting in a stable, laid-down position. The high contact frequency of GGGs hydrocarbon could be attributed to the lying position induced by the other functional groups. Surprisingly, the hydrocarbon of VLP⁻¹ shows a nearly threefold higher contact frequency than its carboxylate group. Smaller, but still significant frequencies (below 2.8%) can be observed for benzol rings and methyl groups.

Atom-resolved density profiles of six selected pollutants (ASA⁻¹, GGG⁻², IBU⁻¹, SLC⁻¹, TTR⁻², VLP⁻¹) along the distance to the TiO₂ surface provide direct insight into the embedding of functional groups within the well-ordered near-surface water layers and preferential adsorption orientation of pollutants (Fig. 5). As a mean of comparison, the densities of water, Na_[0]⁺ and Cl_[0]⁻ of the reference simulation are also plotted (gray lines). The densities are calculated for each atom of the functional groups of the respective pollutant and averaged over the simulation time. The so-obtained profiles account for the different sizes of the functional groups, which is relevant for adsorption on amorphous materials surfaces. All densities are normalized to their maximum value (water normalized to its bulk value).

Fig. 5 Ion and surface densities over z axis as well as pollutant atom densities close to the surface (upper half), divided into functional groups. The densities of the reference calculation (from Fig. 2A) are shown in the background. Heights over surface are calculated with an offset to match the reference surface density height. Shown are the values for ASA⁻¹ (A), GGG⁻² (B), IBU⁻¹ (C), SLC⁻¹ (D), TTR⁻² (E) and VLP⁻¹ (F).

In general, we observe that the molecules remain mainly within the outer Helmholtz and diffusive layer, whereas the inner Helmholtz layer (up to ∼2.5 Å) is mostly comprised of water and Na⁺ ions filling the molecular-scale crevices of the amorphous surface. The only exceptions are carboxylate groups, which can penetrate well into the inner Helmholtz layer. Especially GGG⁻², and, to a lesser extent, ASA⁻¹, IBU⁻¹, SLC⁻¹ and VLP⁻¹, provide density signals in the first water layer above the TiO₂ surface.

Regarding GGG⁻², the densities of phosphonate and carboxylate span a broad range: from 1.8 Å to 9.0 Å for phosphonate and from 0.5 Å to 7.0 Å for carboxylate. The deeper protrusion of the carboxylate into the amorphous surface, combined with a slightly higher contact frequency, indicates that its smaller size allows it to pierce more easily through the inner Helmholtz layer. The amino and hydrocarbon atoms are located at the center of the overall density profile (from 3.5 Å to 7 Å). These profiles suggest a tilted adsorption mode, with the carboxylate group well within the inner Helmholtz layer, whereas the majority of the phosphonate group is found in the diffusive layer. This is confirmed by the relative distances of the GGG⁻² functional groups (cf. SI A.8), where the carboxylate group (first peak at 1.9 Å) moves closer to the surface than the phosphonate. A representative snapshot in Fig. 6B reveals cooperative adsorption of multiple Na⁺ ions and multiple GGG⁻² molecules forming mutual salt bridges among themselves and with the O atoms of the TiO₂ surface.

Fig. 6 Typical adsorption geometries of pollutants at the TiO₂ surface. Shown are snapshots for ASA⁻¹ (A), GGG⁻² (B), IBU⁻¹ (C), SLC⁻¹ (D), TTR⁻² (E) and VLP⁻¹ (F). It is evident that the adsorption takes place mainly between the deprotonated functional group and one hydrogen atom of a TiO₂ surface site or over a connecting ion. Distances indicate direct contacts between the pollutant, ions and the surface and are shown in Å.

Comparing the functional group contact frequencies of VLP⁻¹ with its height density distribution, it appears that the carboxylate, rather than the hydrocarbon, dominates the contact. An explanation could be that the carboxylate group drives the attraction, but the uneven surface morphology causes the molecule to hit the crests with its bulky hydrocarbon before reaching the equilibrium distances between the carboxylate and the surface.

The other functional groups (hydrocarbon, benzol, methyl, hydroxyl, carbonylic acid) in ASA⁻¹, IBU⁻¹, SLC⁻¹ and VLP⁻¹ follow with growing surface distance in accordance with their molecule structure. They indicate preferential adsorption orientations dictated by the formation of contacts between the carboxylate groups and the H atoms of the terminal hydroxyl groups of the surface (see also Fig. 4D and E). These preferential orientations and contacts are indeed visible in representative simulation snapshots of the adsorbed pollutants (Fig. 6A, C, D and F). Notable is the direct involvement of Na⁺ in mediating the contacts, as shown for SLC⁻¹ or TTR⁻² in Fig. 6D and E.

The density profiles of TTR⁻² are indicative of its distinctive adsorption mode. Here, the hydroxyl, ketone, and carboxamide groups are closest to the surface at similar distances (between 1.5 Å to 5 Å), while the nitrogen atom is mostly found at distances around 7.5 Å. As the former groups are mainly found on one side of the molecule, this indicates that TTR⁻² primarily adsorbs to the surface with this side. This can be confirmed when looking at representative snapshots of the MD trajectory (Fig. 6E), where the closest contacts are formed via ketone, hydroxyl, and carboxamide groups, again with direct involvement of Na⁺ ions.

Paying closer attention to the continuous contact time τ_con (Fig. 7A–C), a clear trend similarity with the contact frequency (cf.Fig. 4) is evident. Especially the molecules with highest contact frequencies within their respective protonation state, namely ACE⁰, TTR⁻¹ and TTR⁻², also show the longest continuous retention times at the surface (48 ps, 116 ps and 217 ps, respectively). While CBZ⁰ has the lowest overall contact frequency with the surface, its retention time is surprisingly high, especially compared to other pollutants. The second species that opposes the trend is GGG⁻², where τ_con is short (81 ps) compared to TTR⁻². To understand the mismatch in contact frequency, binding motifs and retention time, we calculated the number of contacts for the two molecules with the highest contact frequency, namely GGG⁻² and TTR⁻² (Fig. 7D–G). For each pollutant, one representative molecule is chosen in contact with the surface and another fully solvated, respectively. The contacts are provided per molecule for water and other pollutant molecules, but per atom for surface and ions. Each grey window highlights 2 ns for different states, respectively. With a total of 5000 frames in that time span, the average contact numbers are listed in the SI A.9.

Fig. 7 Continuous contact time τ_con (A)–(C) and number of water, ion, surface and other pollutant contacts (D)–(F) of two GGG⁻² and TTR⁻² molecules. (A)–(C) The averaged continuous direct and indirect contact time per pollutants in ps reveals similarity to the overall contact frequency. (D)–(H) Running average (window: 200 ps) of number of contacts for two individual molecules of GGG⁻² (D)-(E) and TTR⁻² (F)-(G) over 50 ns. Grey areas show 2 ns windows that are evaluated as different states of the molecules. Insets show exemplary snapshots at these times.

GGG⁻² shows a loss of two water molecules in solution when an ion approaches it (drop of 19 to 17 water molecules). When it contacts the surface, it loses two to three more water molecules, depending on the presence of ions and adsorption geometry. The liberation of water molecules upon contact with TiO₂ could indicate an entropy gain, favoring the adsorbed state of GGG⁻².

TTR⁻² has around 24 to 25 water contacts in solution, also losing two of them when forming a contact with an ion (Fig. 7F). Up to three ions can be found with this individual TTR⁻² molecule, which is in contact with another TTR⁻² at the same time, and the number of water contacts decreases to 19, probably due to the molecular rearrangement of water, ions and the other TTR⁻². At the surface, the TTR⁻² has few contacts to Na⁺, but water contacts between 23 and 25 still. Less water is liberated upon adsorption of TTR⁻², thereby indicating that entropic effects may be less important here. Even more so, we can see water inside the grooves of the surface beneath TTR⁻², which leads us to the conclusion that water could stabilize the adsorbed configuration onto TiO₂. The number of H-bridges per frame (cf. SI A.10) between the pollutant molecules and TiO₂ confirms a big network of water around TTR⁻² molecules.

In summary, hydrogen atoms of the surface terminal groups serve as the primary anchoring points for negatively charged groups of the pollutants. The formation of ionic bridges to the adsorbing molecules via surface-bound Na⁺ ions is often observed. Direct pollutant-surface contacts with the surface Ti atoms are exclusively formed by carboxylate groups. This aligns with other studies reporting that carboxylate preferentially adsorbs on crystalline titania surfaces rather than water,^92,109,110 as a prerequisite for surface-induced molecular degradation processes.

4.3.3 Pollutant clustering. As the clustering of the pollutants is a competing factor for the direct surface adsorption, we now proceed to assess the clustering propensity of the pollutants. This can be quantified by analyzing the occurrence of intermolecular contacts during the MD simulations. Non-clustering pollutants (Fig. 8A) are hardly ever observed to undergo (temporary) contacts with more than three other molecules at a time and have a cluster-size distribution with a value of 1 over 50% of the time, indicating a higher probability of being isolated than in aggregates. Considered cluster contacts are relative distances smaller than 4 Å to account for nucleic acid ring stacking distances (3.4 Ź¹¹). Only in GGG⁻² distances up to 6 Å are included, since the addition of ion bridges partly extends the distances between individual GGG⁻² molecules beyond 4 Å, even though they still belong to the same cluster. Additionally, the formation of clusters will be discussed through the evolution of the combined solvent-accessible surface area (SASA) of pollutants and TiO₂ (Fig. 9A and B, Fig. SI A.11.1).

Fig. 8 Intermolecular pollutant competition based on cluster size analysis. Both diagrams show the proportion of cluster sizes in discrete numbers of molecules per cluster over the whole trajectory (100 ns), sorted by their tendency to cluster. (A) With more than 50% of simulation time as a single molecule, the pollutants are categorized as non-clustering. ACE⁰, ASA⁻¹, BEN⁰, DCL⁰, IBU⁻¹, SLC⁻¹, TTR⁻² and VLP⁻¹ are non-clustering, all of them depict a steep decrease in cluster size fractions with rising cluster size. Only TTR⁻² is found in dimers more often than the other molecules in this category. (B) CAF⁰, CBZ⁰, GGG⁻² and TTR⁻¹ are classified as clustering pollutants. While GGG⁻² is mostly found alone or clusters of 9 molecules, TTR⁻¹ can be observed in clusters up to 9 molecules in significant numbers. CAF⁰ and CBZ⁰ mostly exhibit large clusters with 20 or more molecules. (C)–(F) Snapshots of clustering pollutants. (C) TTR⁻¹ forms amorphous clusters that can create direct contacts with the surface. (D) CAF⁰ shows ladder like structures with purine ring stacking. (E) CBZ⁰ also shows big amorphous clusters like TTR⁻¹ but resides in solution. (F) GGG⁻² forms cluster that are mainly connected over Na⁺.

Fig. 9 Solvent accessible surface area (SASA, (A) and (B)), time resolved number of contacts (C) and (D) and contact matrix (E) and (F) of GGG⁻² and TTR⁻², respectively. (A) and (B) The combined pollutant and TiO₂ SASA of the initial geometry before equilibration (single point data at −1 ns) and over the production run time. In both cases, the SASA is decreasing. (C) and (D) The sum of the number of pollutant contacts with water, Na⁺, TiO₂ and itself over the production run for all 25 pollutant molecules in the cell. For both pollutants, a jump in surface contacts is visible around 20 ns. GGG⁻² also has a steep increase in contacts with itself and a decrease in water contacts, which indicates clustering. (E) and (F) Matrix of contact number drifts over the entire simulation time of 100 ns. The arrows are scaled w.r.t. the amount of contacts.

We found that ACE⁰, ASA⁻¹, BEN⁰, DCL⁰, IBU⁻¹, SLC⁻¹, TTR⁻² and VLP⁻¹ preferentially stay in solution or at the surface as single molecules (Fig. 8). With the exception of TTR⁻², their corresponding SASA reveals constant mean values throughout the entire simulation (Fig. SI A.11.1). However, the fluctuations for non-adsorbing species such as BEN⁰ and DCL⁰ are significantly smaller than those of adsorbing species such as ACE⁰, ASA⁻¹, IBU⁻¹, SLC⁻¹ and VLP⁻¹. The SASA of TTR⁻² shows an overall decrease, which can be attributed to surface adsorption.

The clustering pollutants (CAF⁰, CBZ⁰, GGG⁻² and TTR⁻¹) show clear tendencies to form larger molecular aggregates (Fig. 8B). Here, the SASAs show significant decreases compared to their initial configurations. A clear distinction is observed in the development of SASAs of CAF⁰ and CBZ⁰ compared to GGG⁻². CBZ⁰ tends to form very large aggregates (Fig. 8E) and is found in clusters of 20 or more molecules during about 75% of the simulation time. CAF⁰ forms ladder-like structures by ring-stacking (Fig. 8D). During the dynamics, the molecules frequently stack and break apart, covering sizes from 10 and 14 molecules, but they are mostly found as a long aggregate of nearly all molecules in the simulation cell (66.4%). The driving force for the formation of the molecular stack, which is observed already during the MD equilibration stage, appears to be hydrophobic collapse.

GGG⁻², although staying as a single molecule for 26.1% of the simulation time, forms ionic clusters in which the molecules interact with each other through Na⁺-bridges (Fig. 8F). This is observed also when the molecules are in contact with the surface, forming adsorbed aggregates. Unlike the hydrophobic cluster formation, the formation of ionic clusters takes longer and a continuous reduction of the SASA is observed up to 20 ns into the production run (Fig. 9A). A steady-state level of the SASA for the remaining simulation time points to stabilized ionic clusters of 9 molecules, which cover nearly a third (29.3%) of the simulation time (Fig. 8B). Larger clusters are rarely observed.

TTR⁻¹ is mostly found alone or in small clusters of 2 to 8 molecules, with a probability ranging from 8.4% to 17.9%. The trajectory reveals that TTR⁻¹ clusters are amorphous and can approach the surface to form direct contacts (Fig. 8C), without disrupting the aggregate. However, a significant reduction in the SASA in the magnitude of CAF⁰'s and CBZ⁰'s reduction is observed for TTR⁻¹ at the beginning of the production run with respect to their initial configurations. During the simulation, the significant oscillations could be correlated to the reaggregation and readsorption of molecular clusters.

Our simulations thus reveal multiple driving forces for the formation of pollutant clusters, of both hydrophilic (salt bridges) and hydrophobic (ring stacking) nature. In the first case, observed for GGG⁻², Na⁺ ions drive the association of small clusters and stay in their core. Since GGG⁻² is negatively charged on both ends, it is arranged in a negatively charged shell that can approach the surface and make direct contacts through the free anionic group. In the last few decades, several studies have investigated the influence of water hardness, pH, salt formulation, and additives on the effectiveness of glyphosate^112–115 in plants. Water hardness, which is due to increased Ca²⁺, Mg²⁺ and other ion concentrations, was in fact found to lower the effectiveness. This was attributed to the tendency of deprotonated glyphosate to form salt crystals that impede uptake in plants. Our simulations, despite involving only monovalent cations, are in line with this observation, indicating a strong tendency to form salt clusters that, due to their size, remain in solution but could still limit the uptake of glyphosate by plants to some extend.

The striking stacking behavior of CAF⁰ has already been observed in other molecular simulation studies^58,116–118 and can be compared to the self-aggregation of purine bases in DNA. Regarding CBZ⁰, the high content of neutral MEP area indicates that it is primarily hydrophobic, which thus drives a hydrophobic collapse¹¹⁹ to minimize contact with the polar solvent. We assume that the aggregation of TTR⁻¹ also results from hydrophobic effects, as described by Nabavi et al.¹²⁰ The larger charge of TTR⁻², instead, prevents the formation of clusters larger than three molecules, so that adsorption to the surface as single molecules prevails. More detailed investigations of the clustering behaviors of the pollutants will be presented in a follow-up study.

4.3.4 Adsorption and clustering in competition. Next, we discuss the most frequently adsorbed pollutants GGG⁻², TTR⁻¹ and TTR⁻². In both protonation states, TTR shows longer surface-retention times, while GGG⁻² has a higher tendency to form ionic clusters, but still exhibits high contact frequencies. The SASA and number of contacts between the pollutants, water, ions and the surface can provide further information on the driving forces. Here, we present these results for the double negatively charged pollutant species (GGG⁻², TTR⁻²). In Fig. 9A and B, the development of the SASA over the simulation time is shown in direct comparison with the number of contacts formed with the pollutant molecules (Fig. 9C and D).

A clear decrease in SASA is visible for both pollutants. The GGG⁻² system shows a decrease and convergence within 20 ns, while the TTR⁻² system shows a more gradual decrease throughout the entire simulation. The reduction of the SASA for GGG⁻² is accompanied by a sharp increase in pollutant self-interaction from 20 up to about 40 contacts for all molecules (1.6 per molecule) and approximately 20 TiO₂ surface contacts (0.8 per molecule). Unlike these observed plateaus, the number of ionic contacts increases steadily and the number of contacts with water decreases over the entire simulation.

Fig. 9E shows contact number drifts comparing the first (t = 0 ns to 1 ns) to the last nanosecond (t = 99 ns to 100 ns) of the production run. In GGG⁻², a nearly equal amount of new GGG⁻²–surface, ion–surface and self-contacts are formed. This is balanced by a great loss of GGG⁻²–water contacts and a significant amount of surface–water, that is replaced by GGG⁻² and Na⁺. As previous results showed (Fig. 3B and C and 6B), the formation of GGG⁻²–Na⁺ cluster with the cation as a bridge between two phosphonate groups is evident and explains the absence of the phosphonate signal to surface metal ions in Fig. 4F. This deviates from our results in Balzaretti et al.,⁹² where a bidentate adsorption configuration involving the carboxylate and phosphonate group was found. We report in Balzaretti et al.⁹² an electronic enrichment in the Ti_surf–O_GGG bond, that causes electronic depletion of the neighboring C–C and P–C bonds and forms the prerequisite for surface induced degradation of GGG⁻². The preferred degradation pathway, either the AMPA or the sarcosine route, has been discussed to depend on the presence of light conditions⁹² and on the pH and concentration of radical species.⁷⁵ We thus conclude that the presence of ions due to the salt concentration in wastewater directly influences the degradation pathways and that the obstruction of the phosphonate could explain a probability shift towards the AMPA degradation pathway in experimental studies.¹²¹

TTR⁻² shows stable numbers of contacts for ions and molecules, but an increase of up to 50 molecule–surface contacts (Fig. 9D, 2 contacts per molecule). This points towards a strong mono-molecular adsorption of TTR⁻² (Fig. 6E) compared to the ionic clustered adsorption of GGG⁻² (Fig. 6B). This can also be seen in the contact matrix of TTR⁻² (Fig. 9F): in contrast to GGG⁻², the number of TTR⁻²–water contacts shows a reverse, positive trend. This speaks for a high involvement of water in the surface adsorption of TTR⁻². In conclusion, the adsorption of GGG⁻²–ion clusters leads to a significant release of water, hence an increase in entropy, compared to the mono-molecular TTR⁻², and promotes the GGG⁻² adsorption even though it is prone to aggregate.

5 Discussion

Our simulations reveal that the net charge of the pollutants (sum of protonation states of individual functional groups) plays a key role in determining their adsorption propensities and orientations. First, although the PZC of the TiO₂ surface is lower than the simulated pH and thus the bare surface has a net negative surface potential, the strong adsorption of water molecules and Na⁺ ions promotes the formation of an electrical double-layer^26,30,122 with an overall positive zeta potential, as shown by the accumulation of Cl⁻ ions in the diffuse layer. It is well known that water and cations in the inner Helmholtz layer, because of their reduced diffusivity and permittivity, effectively screen surface charges and thin the electrostatic double-layer near the surface.^27,123,124 For this reason, direct surface contacts can be made with pollutants carrying a net negative charge, which settle, like Cl⁻ ions, at the outer Helmholtz layer.

Second, most direct surface contacts are formed through net or partially negatively charged functional groups, with some (such as carboxylate groups) penetrating into the inner Helmholtz layer, thus stabilizing the surface–molecule interactions. This aligns with a number of theoretical and experimental works by Köppen and Langel,³⁴ Balzaretti et al.,⁹² Balajka et al.,¹⁰⁹ Wen and Selloni,¹¹⁰ and Raman and Selloni,¹²⁵ which found formic and acetic acid as well as carboxylate groups of bigger organic molecules inside the interfacial water layer of solvated crystalline anatase and rutile surfaces. In this layer, the pollutants are either hydrogen-bonded to terminal hydroxyl groups, which are always present due to the amphoteric nature of TiO₂ surfaces,¹²⁶ form salt bridges mediated by the adsorbed Na⁺ ions or attach directly to surface metal ions. Since these contacts are essential for the surface-induced molecule degradation, it is noteworthy that salt concentration and surface topography directly influence the preferred adsorption configuration compared to crystalline, clean systems. While the carboxylate group remains directly bound to the surface, the phosphonate group tends to form ionic clusters at the expense of surface adsorption. However, the only pollutant with a region of the molecular surface at positive electrostatic potential, DCL⁰, forms hardly any contact with the surface (Fig. 4D). Notably, no correlation was observed between the molecular dipole strength and direct contact probability in our simulations.

In studies by Che Abdul Rahim et al.,¹²⁷ Chen et al.,¹²⁸ and Farner Budarz et al.,²⁸ a decrease in pollutant degradation efficiency was observed with increasing salt concentration, which is attributed to the anion content in the solution. Chen et al.¹²⁸ and Farner Budarz et al.²⁸ found only a slight influence of Cl⁻ on degradation rates compared to other anions. While an increase in anion concentration could lead to competition for surface adsorption by anionic pollutants, we cannot draw definitive conclusions on this from our simulations because the high negative charge density (carried by the anionic pollutants) in the simulation cell prevents adding more explicit Cl⁻ without an unphysically high ionic strength. We also note that increasing the salt concentration beyond a certain limit would promote the agglomeration of TiO₂ particles according to colloidal theories.¹²³ Therefore, the salt concentration in experimental works should be chosen as a trade-off between electrical double-layer screening ability and particle suspension stability.

Our findings suggest that the capacity to form hydrogen bonds (number of H-acceptors) and the ability to aggregate in clusters significantly influence the adsorption capability. This can be observed in the behavior of TTR⁻¹, which shows a higher level of adsorption than smaller molecules that also carry a single negative charge.

Also the direct comparison of GGG⁻² and TTR⁻² as well as TTR⁻¹ and TTR⁻² indicates the greater significance of H-acceptor sites and the clustering behavior, which is, in other words, molecular competition: TTR⁻² adsorbs for a longer overall simulation time than the smaller GGG⁻², which can be attributed to its higher number of H-acceptors, and the limited molecular competition, but its counterpart TTR⁻¹ probably adsorbs less due to its higher tendency to aggregate and the single negative protonation state. However, TTR⁻¹ and GGG⁻² clusters are still able to form surface contacts, which is not the case for other pollutants such as CAF⁰ and CBZ⁰. In neutral pollutants, the clustering driving force caused by the hydrophobic collapse of molecules in water^{58,116–119,129,130} appears to dominate and prevent the formation of direct surface contacts entirely. Nevertheless, a significant number of indirect contacts are also formed by these and other neutral pollutants (see Fig. 4). In these cases, degradation could still occur through indirect mechanisms mediated by hydroxyl radicals formed near the surface. In fact, BEN⁰, CAF⁰, CBZ⁰, and DCL⁰ have been found experimentally to be degraded by TiO₂.^11,13,15,81

The simulations further reveal that individual form factors of the pollutants, such as size, bulkiness or the local distribution of functional groups within the molecule, don't significantly influence adsorption. Instead, the presence of specific negatively charged functional groups in the pollutant molecule is identified as the main factor driving adsorption. However, the sheet-like, two-dimensional carboxylate group of GGG⁻² reaches deeper into the water layers than the phosphonate group or carboxylate groups of other pollutants. Although prone to clustering, GGG⁻² can penetrate the water layer at the surface. It achieves a high contact rate with many, but shorter surface contacts by repeatedly overcoming the water barrier through its small size and strong negative net charge.

Overall, we observe contact patterns similar to those in crystalline surface studies,^34,125,131 but scattered across the uneven amorphous surface and with a narrower width of oriented water, resulting in a thinner electrochemical double-layer. This validates the transfer of knowledge from crystalline to amorphous surfaces, providing a more complex picture with these differences in mind.

We finally observe that MD studies using amorphous models are more challenging than those with crystalline surfaces. However, comparing with experimental studies and applications, which mainly rely on inexpensive amorphous powders, is more accurate than comparisons based on crystal interfaces.

6 Conclusions

The present study highlights the complexity of interactions governing the essential first step of pollutant degradation on TiO₂ surfaces, which is adsorption with sufficiently long residence times. Our model system remains simple compared to real wastewater, which contains many more components than the four studied here, specifically TiO₂, water, cations, and a pollutant species. Nevertheless, our MD simulations already identify a first set of rules of thumb that could potentially enhance the degradation efficiency of TiO₂-based methods.

We anticipate that thorough adsorption of polar and anionic pollutants will occur at neutral pH, especially those containing carboxylate groups. Additionally, a high number of hydrogen-bond acceptor sites on the pollutant is beneficial, as it can facilitate the anchoring to the surface through a favorable hydrogen-bond network of water molecules and thus lengthening the retention time.

Also, the ability to aggregate (molecule–molecule competition) significantly influences the propensity for adsorption. Moreover, the form factor of negative functional groups plays a subtle yet important role, as they are more capable of piercing through the inner Helmholtz layer to establish direct contact with the surface.

Most of these effects are heavily influenced by the exact protonation state of the molecule and surface, which can be adjusted by changing the pH of the solution. Future studies will thus be devoted to understanding the influence of pH on the adsorption frequencies and modes of pollutants.

Author contributions

Maria von Einem: simulation protocol development, parametrization of molecule charges, FF and DFT simulation performance, analysis design and programming, writing – original draft. Filippo Balzaretti: analysis design. Manuela Romero: simulation protocol development, collection of representative pollutant species. Wilke Dononelli: interpretation of DFT data, discussion of results. Lucio Colombi Ciacchi: supervision, interpretation of results, writing – review & editing. Giancarlo Franzese: supervision, interpretation of results, writing – review & editing. Susan Köppen-Hannemann: conceptualization of the study, supervision, interpretation of results, writing – review & editing.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Supplementary information (SI) contains a list of tables with partial charges of all pollutants, a detailed description of the simulated system cells, additional descriptions of the analysis calculation, contact maps of all individual pollutant molecules on titanium dioxide and additional results. See DOI: https://doi.org/10.1039/d5ma01085c.

The data for this article, including all input and simulation setups as tpr files and trajectories with a reduced sampling frequency by a factor of ten are available at NOMAD at DOI: https://doi.org/10.17172/NOMAD/2025.09.20-2. In addition, a complete set of data of the ESP partial charge calculations are available at NOMAD, see DOI: https://doi.org/10.17172/NOMAD/2025.09.20-1.

Acknowledgements

We are very grateful for the funding of this project through the DFG research training group 2247 ‘QM3-quantum mechanical materials modeling’. We acknowledge financial support through the APF project ‘Materials on Demand’ within the ‘Humans on Mars’ Initiative funded by the Federal State of Bremen and the University of Bremen. WD gratefully acknowledges funding by the Central Research Development Fund of the University of Bremen. GF acknowledges funding by MCIN/AEI/10.13039/501100011033/ERDF/EU [grant numbers PID2021-124297NB-C31; PID2024-157478NB-C31]. MVE and GF acknowledge funding through the Project HPC-EUROPA3 (INFRAIA-2016-1-730897) of the EC Research Innovation Action under the H2020 Programme. This provided support for close research collaborations between the University of Bremen and the University of Barcelona, as well as the provision of computational resources and technical support by the Barcelona Supercomputing Center (BSC). MVE, FB and SK acknowledge computational resources that have been allocated on the supercomputers of the National High Performance Computing Alliance (NHR, formerly Northern German Supercomputing Alliance, HLRN).

References

1. F. A. Caliman and M. Gavrilescu, Clean: Soil, Air, Water, 2009, 37, 277–303.
2. X. Zhang, F. Wu, X. Wu, P. Chen and N. Deng, J. Hazard. Mater., 2008, 157, 300–307.
3. A.-M. Abdel-Wahab, A.-S. Al-Shirbini, O. Mohamed and O. Nasr, J. Photochem. Photobiol., A, 2017, 347, 186–198.
4. J. Borràs-Ferrís, R. Sánchez-Tovar, E. Blasco-Tamarit, M. J. Muñoz-Portero, R. M. Fernández-Domene and J. García-Antón, Nanomaterials, 2019, 9, 583.
5. D. Mukherjee, A. K. Ray and S. Barghi, Processes, 2016, 4, 13.
6. N. Bolong, A. F. Ismail, M. R. Salim and T. Matsuura, Desalination, 2009, 239, 229–246.
7. E. Moctezuma, E. Leyva, C. A. Aguilar, R. A. Luna and C. Montalvo, J. Hazard. Mater., 2012, 243, 130–138.
8. M. D. G. de Luna, J. C.-T. Lin, M. J. N. Gotostos and M.-C. Lu, Sustainable Environ. Res., 2016, 26, 161–167.
9. S. Chen and Y. Liu, Chemosphere, 2007, 67, 1010–1017.
10. C. L. Bianchi, B. Sacchi, C. Pirola, F. Demartin, G. Cerrato, S. Morandi and V. Capucci, Environ. Sci. Pollut. Res., 2017, 24, 12646–12654.
11. M. K. Arfanis, P. Adamou, N. G. Moustakas, T. M. Triantis, A. G. Kontos and P. Falaras, Chem. Eng. J., 2017, 310, 525–536.
12. A. Achilleos, E. Hapeshi, N. P. Xekoukoulotakis, D. Mantzavinos and D. Fatta-Kassinos, Sep. Sci. Technol., 2010, 45, 1564–1570.
13. A. Carabin, P. Drogui and D. Robert, J. Taiwan Inst. Chem. Eng., 2015, 54, 109–117.
14. J.-K. Im, H.-S. Son, Y.-M. Kang and K.-D. Zoh, Water Environ. Res., 2012, 84, 554–561.
15. J. Borisch, S. Pilkenton, M. L. Miller, D. Raftery and J. S. Francisco, J. Phys. Chem. B, 2004, 108, 5640–5646.
16. X.-D. Zhu, Y.-J. Wang, R.-J. Sun and D.-M. Zhou, Chemosphere, 2013, 92, 925–932.
17. D. Haranaka-Funai, F. Didier, J. Giménez, P. Marco, S. Esplugas and A. Machulek-Junior, Chem. Eng. J., 2017, 327, 656–665.
18. S. Castiglioni, E. Davoli, F. Riva, M. Palmiotto, P. Camporini, A. Manenti and E. Zuccato, Water Res., 2018, 131, 287–298.
19. M. Xu, Y. Gao, E. M. Moreno, M. Kunst, M. Muhler, Y. Wang, H. Idriss and C. Wöll, Phys. Rev. Lett., 2011, 106, 138302.
20. A. Gautam, A. Kshirsagar, R. Biswas, S. Banerjee and P. K. Khanna, RSC Adv., 2016, 6, 2746–2759.
21. D. Daskalova, G. Aguila Flores, U. Plachetka, M. Möller, J. Wolters, T. Wintgens and M. C. Lemme, ACS Appl. Nano Mater., 2023, 6, 15204–15212.
22. D. Jassby, J. Farner Budarz and M. Wiesner, Environ. Sci. Technol., 2012, 46, 6934–6941.
23. U. Diebold, Surf. Sci. Rep., 2003, 48, 53–229.
24. L. Huang, K. E. Gubbins, L. Li and X. Lu, Langmuir, 2014, 30, 14832–14840.
25. C. Zhang, M. F. Calegari Andrade, Z. K. Goldsmith, A. S. Raman, Y. Li, P. M. Piaggi, X. Wu, R. Car and A. Selloni, Nat. Commun., 2024, 15, 10270.
26. S. Köppen and W. Langel, Surf. Sci., 2006, 600, 2040–2050.
27. H.-J. Butt, K. Graf and M. Kappl, Physics and Chemistry of Interfaces, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, 2013.
28. J. Farner Budarz, A. Turolla, A. F. Piasecki, J.-Y. Bottero, M. Antonelli and M. R. Wiesner, Langmuir, 2017, 33, 2770–2779.
29. W. R. Fawcett, Liquids, Solutions, and Interfaces: From Classical Macroscopic Descriptions to Modern Microscopic Details, Oxford University Press, New York, 2004.
30. M. Bischoff, D. Biriukov, M. Predota, S. Roke and A. Marchioro, J. Phys. Chem. C, 2020, 124, 10961–10974.
31. S. M. Fatemi and S. J. Fatemi, J. Mol. Graphics Modell., 2022, 112, 108123.
32. A. YazdanYar, U. Aschauer and P. Bowen, Colloids Surf., B, 2018, 161, 563–577.
33. W. Friedrichs and W. Langel, Biointerphases, 2014, 9, 031006.
34. S. Köppen and W. Langel, Phys. Chem. Chem. Phys., 2008, 10, 1907–1915.
35. L. Yang, L. E. Yu and M. B. Ray, Water Res., 2008, 42, 3480–3488.
36. N. Jallouli, K. Elghniji, H. Trabelsi and M. Ksibi, Arabian J. Chem., 2017, 10, S3640–S3645.
37. M. O. Miranda, W. E. Cabral Cavalcanti, F. F. Barbosa, J. Antonio de Sousa, F. Ivan da Silva, S. B. C. Pergher and T. P. Braga, RSC Adv., 2021, 11, 27720–27733.
38. J. Schneider and L. Colombi Ciacchi, J. Am. Chem. Soc., 2012, 134, 2407–2413.
39. F. Anwar, A. M. Varghese, S. Kuppireddy, A. Gotzias, M. Khaleel, K. Wang and G. N. Karanikolos, Chem. Eng. J., 2025, 524, 169411.
40. P. Wang, P.-S. Yap and T.-T. Lim, Appl. Catal., A, 2011, 399, 252–261.
41. B. Petrie, R. Barden and B. Kasprzyk-Hordern, Water Res., 2015, 72, 3–27.
42. R. Alajlani and A. Alssadi, Open J. Phys. Chem., 2019, 09, 216–220.
43. N. H. Tran, M. Reinhard and K. Y.-H. Gin, Water Res., 2018, 133, 182–207.
44. L. J. Schulman, E. V. Sargent, B. D. Naumann, E. C. Faria, D. G. Dolan and J. P. Wargo, Hum. Ecol. Risk Assess., 2002, 8, 657–680.
45. T. A. Ternes, Water Res., 1998, 32, 3245–3260.
46. E. P. Munzhelele, R. Mudzielwana, W. B. Ayinde and W. M. Gitari, Water, 2024, 16, 796.
47. X. Wei, X. Bao, J. Wu, C. Li, Y. Shi, J. Chen, B. Lv and B. Zhu, RSC Adv., 2018, 8, 10396–10408.
48. L. Dsikowitzky, J. Schwarzbauer and R. Littke, Chemosphere, 2004, 57, 1289–1300.
49. J. Ulander and A. D. J. Haymet, Biophys. J., 2003, 85, 3475–3484.
50. J. T. Yu, E. J. Bouwer and M. Coelhan, Agric. Water Manage., 2006, 86, 72–80.
51. K.-C. Wang, S.-M. Chen, J.-F. Hsu, S.-G. Cheng and C.-K. Lee, J. Chromatogr. B:Anal. Technol. Biomed. Life Sci., 2008, 876, 211–218.
52. O. V. de Oliveira, J. D. dos Santos, J. C. F. Silva, L. T. Costa, M. F. Ferreira Junior and E. D. F. Franca, Orbital:Electron. J. Chem., 2017, 9, 175–180.
53. M. Feltracco, E. Barbaro, E. Morabito, R. Zangrando, R. Piazza, C. Barbante and A. Gambaro, Environ. Sci. Pollut. Res., 2022, 29, 16383–16391.
54. M. A. Álvarez Bayona, A. Maturana Córdoba, R. J. Gallardo Amaya and A. Muñoz Acevedo, Front. Environ. Sci., 2022, 10, 941836.
55. L. J. Leeson, J. E. Krueger and R. A. Nash, Tetrahedron Lett., 1963, 4, 1155–1160.
56. L. E. Lesser, A. Mora, C. Moreau, J. Mahlknecht, A. Hernández-Antonio, A. I. Ramírez and H. Barrios-Piña, Chemosphere, 2018, 198, 510–521.
57. A. Vandeponseele, M. Draye, C. Piot and G. Chatel, Clean Technol., 2021, 3, 335–350.
58. L. Tavagnacco, S. Di Fonzo, F. D’Amico, C. Masciovecchio, J. W. Brady and A. Cesàro, Phys. Chem. Chem. Phys., 2016, 18, 13478–13486.
59. L. Geiser, Y. Henchoz, A. Galland, P.-A. Carrupt and J.-L. Veuthey, J. Sep. Sci., 2005, 28, 2374–2380.
60. R. Porto, G. de Tommaso and E. Furia, Ann. Chim., 2005, 95, 551–558.
61. P. Pulay, G. Fogarasi and J. E. Boggs, J. Chem. Phys., 1981, 74, 3999–4014.
62. WHO, Benzene in Drinking-water. Background document for development of WHO Guidelines for Drinking-water Quality, 2003, https://cdn.who.int/media/docs/default-source/chemical-safety/benzene-in-drinking-water.pdf?sfvrsn=f4e4aca7_1ua=1, (accessed December 2024).
63. CRC handbook of chemistry and physics, ed. D. R. Lide, CRC Press, Boca Raton, 90th edn, 2009.
64. WHO, Dichloromethane in drinking-water. Background document for development of WHO Guidelines for Drinking-water Quality, 2003, https://cdn.who.int/media/docs/default-source/wash-documents/wash-chemicals/dichloromethane.pdf?sfvrsn=877218cb_4, (accessed December 2024).
65. M. O. Uslu, S. Jasim, A. Arvai, J. Bewtra and N. Biswas, Ozone: Sci. Eng., 2013, 35, 249–262.
66. M. C. N. Pais and E. D. S. Nascimento, Braz. J. Pharm. Sci., 2018, 54, 17320.
67. K. Brotzmann, A. Wolterbeek, D. Kroese and T. Braunbeck, Arch. Toxicol., 2021, 95, 641–657.
68. S. Wang, B. Seiwert, M. Kästner, A. Miltner, A. Schäffer, T. Reemtsma, Q. Yang and K. M. Nowak, Water Res., 2016, 99, 91–100.
69. WHO, Guidelines for drinking-water quality: incorporating the first and second addenda, World Health Organization, Geneva, 4th edn, 2022.
70. I. Georgaki, E. Vasilaki and N. Katsarakis, Am. J. Anal. Chem., 2014, 5, 518–534.
71. N. Qi, P. Wang, C. Wang and Y. Ao, J. Hazard. Mater., 2018, 341, 187–197.
72. F. Saadati, N. Keramati and M. M. Ghazi, Crit. Rev. Environ. Sci. Technol., 2016, 46, 757–782.
73. G. H. Safari, M. Hoseini, M. Seyedsalehi, H. Kamani, J. Jaafari and A. H. Mahvi, Int. J. Environ. Sci. Technol., 2015, 12, 603–616.
74. N. Suciu, E. Russo, M. Calliera, G. P. Luciani, M. Trevisan and E. Capri, Sci. Total Environ., 2023, 866, 161171.
75. M. Sadatsharifi, D. W. Ingersoll and M. Purgel, Environ. Sci.: Processes Impacts, 2021, 23, 1018–1028.
76. C. Gasnier, C. Dumont, N. Benachour, E. Clair, M.-C. Chagnon and G.-E. Séralini, Toxicology, 2009, 262, 184–191.
77. B. Savun-Hekimoglu, Z. Eren and N. H. Ince, Sustainability, 2020, 12, 10314.
78. I. Skarmoutsos, I. D. Petsalakis and J. Samios, Ind. Eng. Chem. Res., 2021, 60, 11834–11847.
79. C. Gui, G. Li, M. Song and Z. Lei, Sep. Purif. Technol., 2023, 311, 123281.
80. B. Huang, C. Lei, C. Wei and G. Zeng, Environ. Int., 2014, 71, 118–138.
81. H. Einaga, S. Futamura and T. Ibusuki, Phys. Chem. Chem. Phys., 1999, 1, 4903–4908.
82. National Center for Biotechnology Information, PubChem. Compound Summary for CID 54675776, https://pubchem.ncbi.nlm.nih.gov/compound/Tetracycline, (accessed December 2024).
83. V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg and C. Simmerling, Proteins: Struct., Funct., Bioinf., 2006, 65, 712–725.
84. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
85. F. Neese, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2012, 2, 73–78.
86. W. Zuo, N. Li, B. Chen, C. Zhang, Q. Li and M. Yan, Sci. Total Environ., 2020, 726, 138432.
87. P. Sprankle, W. F. Meggitt and D. Penner, Weed Sci., 1975, 23, 229–234.
88. L. Derr, N. Hildebrand, S. Köppen, S. Kunze, L. Treccani, R. Dringen, K. Rezwan and L. Colombi Ciacchi, Biointerphases, 2016, 11, 011007.
89. J. Schneider and L. Colombi Ciacchi, Surf. Sci., 2010, 604, 1105–1115.
90. J. Schneider and L. Colombi Ciacchi, J. Chem. Theory Comput., 2011, 7, 473–484.
91. J. Laube, S. Salameh, M. Kappl, L. Mädler and L. Colombi Ciacchi, Langmuir, 2015, 31, 11288–11295.
92. F. Balzaretti, M. von Einem, L. Gerhards, E. Macke, T. Neudecker, L. Colombi Ciacchi, W. Dononelli and S. Köppen-Hannemann, J. Phys. Chem. C, 2025, 129, 13910–13922.
93. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
94. D. Case, I. Ben-Shalom, S. Brozell, D. Cerutti, T. Cheatham, III, V. Cruzeiro, R. D. T. A. Darden, D. Ghoreishi, M. Gilson, H. Gohlke, A. Goetz, D. Greene, R. Harris, N. Homeyer, S. I. Y. Huang, A. Kovalenko, T. Kurtzman, T. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, D. Mermelstein, K. Merz, Y. Miao, G. Monard, C. Nguyen, H. Nguyen, I. Omelyan, A. Onufriev, F. Pan, R. Qi, D. Roe, A. Roitberg, C. Sagui, S. Schott-Verdugo, J. Shen, C. Simmerling, J. Smith, R. Salomon-Ferrer, J. Swails, R. Walker, J. Wang, H. Wei, R. Wolf, X. Wu, L. Xiao, D. York and P. Kollman, AMBER 2018, University of California, San Francisco, 2018.
95. A. W. Sousa da Silva and W. F. Vranken, BMC Res. Notes, 2012, 5, 367.
96. D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
97. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
98. The PLUMED consortium, Nat. Methods, 2019, 16, 670–673.
99. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
100. R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, J. Domański, D. L. Dotson, S. Buchoux, I. M. Kenney and O. Beckstein, Proceedings of the 15th Python in Science Conference, presented in part at SciPy 2016, Austin, Texas, July, 2016.
101. N. Michaud-Agrawal, E. J. Denning, T. B. Woolf and O. Beckstein, J. Comput. Chem., 2011, 32, 2319–2327.
102. Jmol: an open-source Java viewer for chemical structures in 3D, https://www.jmol.org.
103. Z. Zhang, P. Fenter, L. Cheng, N. C. Sturchio, M. J. Bedzyk, M. Predota, A. Bandura, J. D. Kubicki, S. N. Lvov, P. T. Cummings, A. A. Chialvo, M. K. Ridley, P. Bénézeth, L. Anovitz, D. A. Palmer, M. L. Machesky and D. J. Wesolowski, Langmuir, 2004, 20, 4954–4969.
104. K. Bourikas, C. Kordulis and A. Lycourghiotis, Chem. Rev., 2014, 114, 9754–9823.
105. M. Predota, Z. Zhang, P. Fenter, D. J. Wesolowski and P. T. Cummings, J. Phys. Chem. B, 2004, 108, 12061–12072.
106. K. Kiyohara and Y. Kawai, J. Chem. Phys., 2019, 151, 104704.
107. S. Ma, Y. Qin, K. Sun, J. Ahmed, W. Tian and Z. Ma, Materials, 2024, 17, 2930.
108. Y. Yang, L. Chen, M. Li and X. Qiu, Ind. Eng. Chem. Res., 2025, 64, 6956–6967.
109. J. Balajka, M. A. Hines, W. J. I. DeBenedetti, M. Komora, J. Pavelec, M. Schmid and U. Diebold, Science, 2018, 361, 786–789.
110. B. Wen and A. Selloni, J. Phys. Chem. Lett., 2021, 12, 6840–6846.
111. R. E. Dickerson, Nucleic Acids Res., 1989, 17, 1797–1803.
112. J. D. Nalewaja, B. Devilliers and R. Matysiak, Weed Res., 1996, 36, 241–247.
113. I. Travlos, N. Cheimona and D. Bilalis, Agronomy, 2017, 7, 60.
114. J. P. A. R. da Cunha, R. P. Palma, A. C. de Oliveira, M. G. Marques and C. B. D. Alvarenga, Eng. Agric., 2020, 40, 555–560.
115. J. Schönherr and L. Schreiber, J. Agric. Food Chem., 2004, 52, 6546–6551.
116. L. Tavagnacco, U. Schnupf, P. E. Mason, M.-L. Saboungi, A. Cesàro and J. W. Brady, J. Phys. Chem. B, 2011, 115, 10957–10966.
117. L. Tavagnacco, Y. Gerelli, A. Cesàro and J. W. Brady, J. Phys. Chem. B, 2016, 120, 9987–9996.
118. H. Ramézani, I. Ellien, Z. El Oufir, N. Mathieu, S. Delpeux and S. K. Bhatia, Colloids Surf., A, 2023, 673, 131645.
119. A. Gadelmeier, M. Macht and D. Zahn, J. Pharm. Sci., 2022, 111, 2898–2906.
120. M. Nabavi, M. R. Housaindokht, M. R. Bozorgmehr and A. Sadeghi, J. Iran. Chem. Soc., 2020, 17, 2997–3007.
121. S. H. Bai and S. M. Ogbourne, Environ. Sci. Pollut. Res., 2016, 23, 18988–19001.
122. L. Agosta, E. G. Brandt and A. P. Lyubartsev, J. Chem. Phys., 2017, 147, 024704.
123. R. A. French, A. R. Jacobson, B. Kim, S. L. Isley, R. L. Penn and P. C. Baveye, Environ. Sci. Technol., 2009, 43, 1354–1359.
124. B. Chu, D. Biriukov, M. Bischoff, M. Predota, S. Roke and A. Marchioro, Faraday Discuss., 2023, 246, 407–425.
125. A. S. Raman and A. Selloni, Angew. Chem., Int. Ed., 2025, 64, e202507721.
126. K. Lau, B. Giera, S. Barcikowski and S. Reichenberger, Nanoscale, 2024, 16, 2552–2564.
127. A. N. Che Abdul Rahim, S. Yamada, H. Bonkohara, S. Mestre, T. Imai, Y.-T. Hung and I. Kumakiri, Int. J. Environ. Res. Public Health, 2022, 19, 15736.
128. H. Y. Chen, O. Zahraa and M. Bouchy, J. Photochem. Photobiol., A, 1997, 108, 37–44.
129. C. A. Hunter, J. F. McCabe and A. Spitaleri, CrystEngComm, 2012, 14, 7115–7117.
130. I. Rosbottom, T. N. H. Cheng and J. Y. Y. Heng, Chem. Eng. Technol., 2020, 43, 1152–1159.
131. V. Carravetta, S. Monti and W. Zhang, Theor. Chem. Acc., 2009, 123, 299–309.

---