Designing porous molecularly imprinted polymers via simulation of pre-polymerisation mixtures: a case study with trinitrotoluene

Abstract

Selective adsorption of hazardous micropollutants from water remains a critical challenge in sustainable materials design. Herein, we demonstrate a combined computational–experimental approach to rationally engineer molecularly imprinted polymers for targeted porosity, using 2,4,6-trinitrotoluene as a model template. By simulating pre-polymerisation mixtures of monomers, crosslinkers, and solvent using molecular dynamics, we capture key template–monomer interactions and predict the resulting porosity of the final polymer network. Surface area and free volume predictions from simulations show excellent agreement with experimental nitrogen sorption data across varying solvent compositions. Our findings highlight a fundamental trade-off between imprinting efficiency (favoured in acetonitrile-rich environments) and porous structure (promoted by dimethyl sulfoxide). We validate that pre-polymerisation simulations alone can accurately guide formulations toward high-performance materials, opening new pathways for computationally-driven design of porous polymeric adsorbents.

---

Design, System, Application

This study introduces a simulation-informed strategy to guide the design of porous molecularly imprinted polymers (MIPs) based on the structural organisation of pre-polymerisation mixtures. By modelling the interactions between monomers, crosslinkers, templates, and solvents using molecular dynamics simulations, we predict the porosity of the final polymer without needing to simulate the polymerisation reaction itself. The system under investigation targets 2,4,6-trinitrotoluene (TNT), a hazardous water pollutant, serving as a template to demonstrate the interplay between chemical recognition (imprinting) and structural functionality (porosity). Our approach captures how solvent composition modulates both binding fidelity and polymer morphology, revealing a fundamental trade-off between imprinting efficiency and surface area. The ability to anticipate and control porosity computationally, prior to synthesis, has wide-reaching applications in materials for environmental remediation, sensing, and separations. This framework is generalisable to other crosslinked porous systems and significantly reduces the trial-and-error typically associated with polymer formulation, offering a practical design route for advanced functional materials.

1 Introduction

The design of porous functional polymers for selective recognition and separation has become a central challenge in environmental remediation, sensing, and chemical engineering. Waterborne micropollutants pose a significant threat to both ecosystems and human health due to their persistence, bioaccumulation, and resistance to conventional treatment methods.¹ With chemical contamination contributing to over 800 000 annual deaths worldwide and the United Nations Sustainable Development Goals underscoring clean water access (Goals 6 and 14), there is a pressing need for advanced materials that enable more targeted treatment strategies.²

Molecularly imprinted polymers (MIPs) represent a powerful class of synthetic materials capable of selectively recognising target compounds through templated cavities formed during polymerisation.³ Foundational reviews and mechanism studies in molecular imprinting provide the broader context for this work, covering core design principles, stability/reusability in application, and surface-imprinting mechanisms.^4–6 These engineered cavities mimic the structural and functional complementarity found in biological systems and are particularly suited for applications involving trace-level analyte detection or capture. Recent “greenificated” MIP strategies also emphasise sustainability across the materials' life cycle, underscoring the value of predictive design to minimise experimental iteration.⁷ However, successful MIP synthesis is highly sensitive to the pre-polymerisation formulation—including the choice of functional monomer, crosslinker, and, critically, the porogenic solvent.^8–10 The latter not only influences the template–monomer interactions that determine imprinting fidelity, but also governs the final polymer's morphology and porosity, which are key to material performance. The overall synthetic strategy for generating the MIP materials is summarised in Fig. 1.

Fig. 1 Generalised schematic of the molecular imprinting process involving (a) functional monomers, (b) crosslinker and (c) template. The workflow illustrates (1) template–monomer complexation in solution, (2) polymerisation in the presence of porogen, (3) extraction of template and solvent to yield a porous MIP structure with accessible recognition sites and (4) rebinding.

Designing MIPs with both high binding specificity and appropriate surface area remains challenging due to the complex interplay of chemical and physical interactions in the pre-polymerisation mixture. The high dimensionality of the formulation space and the lack of predictive models has meant that MIP development is often reliant on empirical, trial-and-error methods.^11,12 In this context, molecular modelling has emerged as a valuable tool for rational MIP design, enabling insights into monomer–template interactions, solvent effects, and network formation tendencies.^9–14 While quantum chemical methods can estimate binding energies for isolated monomer–template complexes,^15,16 they typically neglect solvent influence and are limited in scope. By contrast, classical molecular dynamics (MD) simulations offer a route to explore the full pre-polymerisation environment—capturing competition, clustering, and emergent behaviour from multi-component mixtures.^10,17

Previous work has demonstrated the utility of MD simulations in screening MIP formulations, particularly in the context of pharmaceutical targets.^10,13 Recent studies have also model pre-polymerisation compatibility in full template–monomer–crosslinker–solvent systems for MIP design¹⁴ however, these do not attempt to predict final polymer porosity directly from pre-polymerisation snapshots or to map porogen-controlled surface area trends. This is the specific gap we address. Here, we extend our simulation-guided approach to explore how solvent composition affects the porosity of the final polymer network. Rather than simulating polymerisation explicitly—a computationally demanding task requiring reactive force fields—we propose that static snapshots of the pre-polymerisation mixture can be used to predict surface area trends. This is particularly relevant for thermosetting systems where porosity is strongly linked to the packing and arrangement of precursor molecules before curing.^8,10

As a case study, we examine a MIP system targeting 2,4,6-trinitrotoluene (TNT), a nitroaromatic compound of environmental concern due to its widespread industrial and military use and associated toxicity.¹⁸ In fact, the UN estimates that 10 million hectares of land worldwide are contaminated with explosives like TNT,¹⁹ with concentrations, detected in military sites, exceeding 100 μg L⁻¹, which is far above safe limits for drinking water (typically <2 μg L⁻¹).¹⁸ Moreover, TNT has been widely studied as a template for MIP-based sensors,^3,9 including fluorescent platforms, reinforcing its relevance as a model system,²⁰ but prior modelling work has not addressed how porogen choice influences polymer morphology. Building on our earlier TNT–monomer interaction studies⁹ and recent advances in porosity prediction,¹⁰ we combine MD simulations and experimental synthesis to investigate how mixed DMSO–acetonitrile (ACN) porogens affect both template–monomer interactions and the resulting material structure. Our aim is to demonstrate that pre-polymerisation simulations can provide predictive insights into the porous architecture of MIPs, thereby guiding materials design prior to synthesis.

This dual focus on chemical binding and physical morphology aligns with emerging strategies in simulation-informed polymer engineering. We anticipate that the same pre-polymerisation-snapshot approach can inform other crosslinked porous materials (e.g., HCPs, COFs, PIMs) where solvent-driven packing dictates morphology;^21–23 related simulation-guided studies in polymer porosity and transport support this transferability.^24,25

2 Computational methodology

2.1 Simulation framework for porosity-directed MIP design

To investigate how solvent composition influences the structural organisation of pre-polymerisation mixtures—and ultimately the porous morphology of the final polymer—we employed classical molecular dynamics (MD) simulations using the GROMACS package.^26–28 This approach allows the explicit treatment of all molecular species involved in MIP formulation and captures collective effects, such as hydrogen bonding networks and molecular clustering, which govern both imprinting and porosity development.¹⁰

Simulation boxes were constructed to model the initial formulation of a molecularly imprinted polymer targeting TNT. Each system contained 10 TNT molecules (template), 60 methacrylic acid (MAA, functional monomer), 250 ethylene glycol dimethacrylate (EGDMA, crosslinker), and 600 solvent molecules. The solvent composition was systematically varied across five ratios of dimethyl sulfoxide (DMSO) and acetonitrile (ACN): 0 : 100, 25 : 75, 50 : 50, 75 : 25, and 100 : 0 (v/v). These formulations mirror those used experimentally, enabling direct comparison between predicted and measured surface area trends.

Molecular interactions were described using the OPLS-AA force field with atom types, bonded parameters, and electrostatic charges obtained from the automated topology builder (ATB).²⁹ This combination has been validated in our previous MIP studies¹⁰ and provides a transferable framework for crosslinked monomer systems.

2.2 Equilibration and sampling protocol

All systems underwent a consistent equilibration protocol, previously optimised for MIP design simulations.¹⁰ After initial energy minimisation to remove steric clashes, systems were equilibrated under isothermal–isobaric (NPT) conditions at 298 K and 1 bar for 2 ns using the Berendsen thermostat and barostat. To enhance conformational sampling and capture diverse packing states, a temperature annealing step was applied by gradually heating the system to 1000 K over 20 ns. Production runs were then performed for 10 ns in the canonical (NVT) ensemble at 298 K, using the Nosé–Hoover thermostat and a 2 fs time step. All bonds involving hydrogen atoms were constrained using the LINCS algorithm. A nonbonded interaction cutoff of 1.0 nm was employed.

To ensure statistical robustness, 18 independent replicas were generated for each solvent composition, enabling ensemble averaging of structural and energetic metrics. Importantly, these simulations capture the distribution and organisation of functional monomers and crosslinkers under different solvent environments—providing a structural basis for predicting polymer morphology.

2.3 Hydrogen bonding and solvent–monomer organisation

To evaluate the chemical interactions governing imprinting fidelity, hydrogen bonding analysis was performed between the carboxylic acid of MAA (donor) and the nitro groups of TNT (acceptors). Hydrogen bonds were identified using geometric criteria (donor–acceptor distance ≤0.35 nm, H–D–A angle ≥150°) and quantified as the percentage of simulation time a given pair remained bonded. Trajectories were centred and unwrapped prior to analysis. For each solvent composition, we analysed the final 10 ns of the NVT production trajectory, sampling every 1 ps. Percent “time bonded” values in Table 1 are averages over 18 independent replicas; uncertainties are reported as ±0.1s.d. across replicas. Analysis was performed with gmx hbond and validated against an in-house script implementing identical geometric criteria. Additionally, radial distribution functions (RDFs) were computed between MAA hydroxyl hydrogens and solvent acceptor atoms (DMSO oxygen and ACN nitrile nitrogen), providing insight into competitive binding and solvent structuring around the monomer.

Table 1 Hydrogen bonding between the TNT template and MAA functional monomer, given as the percentage of simulation time (out of 10 ns) that a hydrogen bond is present. Values for the ortho and para nitro groups of TNT are averaged (TNT has two ortho and one para nitro group). Values are mean ± 0.1 s.d. across 18 replicas; geometric criteria and sampling are detailed in 2.3 Hydrogen bonding and solvent–monomer organisation

Porogen	H-bonding time (%)
(DMSO : ACN)	Ortho	Para
0 : 100	93.0	95.4
25 : 75	90.8	93.7
50 : 50	88.8	92.8
75 : 25	85.3	91.2
100 : 0	85.7	90.4

---

These analyses collectively characterise the extent to which solvent composition affects template–monomer complexation and monomer clustering—both of which are known to influence the imprinting process.^8–10

2.4 Porosity prediction from pre-polymerisation snapshots

To predict the porous characteristics of the final polymer, we extracted 20 structural snapshots from the last 5 ns of each production trajectory. These configurations were stripped of all solvent and template molecules, leaving only the monomer and crosslinker matrix. This static representation approximates the network topology immediately prior to polymerisation, assuming that the spatial arrangement of monomer units is preserved upon curing—an approach validated in prior work.^8,10

Surface area was evaluated using two complementary computational approaches:

2.4.1 MeshSA. Grand canonical Monte Carlo (GCMC) simulations using nitrogen as a probe (mimicking BET measurements) were performed with DL MONTE.^30,31 Probe-accessible volume was mapped onto a 3D grid and visualised using VMD³² to compute accessible surface area.

2.4.2 FreeSASA. Solvent-accessible surface area was calculated via the Shrake–Rupley algorithm using a 0.14 nm probe and 200 test points per atom.^33,34 This geometric approach quantifies surface exposure without requiring adsorption modelling.

Fig. 2 illustrates how solvent-accessible surface area (denoted with the dashed-blue region) was computed using spherical probe methods.

Fig. 2 Illustration of solvent-accessible surface area (SASA) estimation using 1.4 Å probe radius, green sphere (FreeSASA) rolled over the MIP atoms, red spheres. Blue mesh regions denote regions accessible to solvent prior to crosslinking.

Fractional free volume (FFV) was also determined using the GROMACS freevolume utility, which computes the unoccupied volume fraction based on van der Waals radii.

By correlating these simulated structural descriptors with experimentally measured BET surface areas (see section 4), we demonstrate that pre-polymerisation MD simulations can serve as a reliable, predictive tool for tuning porosity in MIP design. This framework eliminates the need to explicitly simulate polymerisation or assess adsorption performance—instead, it provides a rapid and computationally efficient route to formulate porous, functional polymers.

3 Experimental methodology

3.1 Materials and equipment

All reagents were used as received without further purification. The template molecule 2,4,6-trinitrotoluene (TNT), functional monomer methacrylic acid (MAA), crosslinker ethylene glycol dimethacrylate (EGDMA), and initiator azobisisobutyronitrile (AIBN) were procured from commercial suppliers. Analytical-grade acetonitrile (ACN) and dimethyl sulfoxide (DMSO) were used as porogenic solvents. Ultrapure water (18.2 MΩ cm) was used for cleaning and preparation.

Nitrogen sorption measurements were carried out using a Quantachrome Autosorb iQ instrument to assess surface area and porosity. UV-visible spectra were collected on an Agilent Cary 100 spectrophotometer to confirm template presence during synthesis. Scanning electron microscopy (SEM) imaging was performed using a Hitachi SU3900 microscope. Prior to imaging, polymer powders were sputter-coated with gold using an Edwards 150B coater.

3.2 Polymer synthesis

Polymer synthesis was carried out following the same monomer–template ratios and solvent compositions used in the simulations, enabling direct comparison between predicted and experimental morphology.

For each formulation, TNT (220 μL of a stock solution corresponding to 10 mmol) and MAA (110 μL, 60 mmol) were first mixed in a glass vial (1,6 molar ratio). EGDMA (1.10 g, 250 mmol) and AIBN (36 mg, 2.2 mmol) were added, and the mixture was dissolved in 8 mL of solvent comprising one of five ACN : DMSO volume ratios: 0 : 100, 25 : 75, 50 : 50, 75 : 25, or 100 : 0. Solutions were purged with nitrogen gas for 5 minutes, sealed, and polymerised at 60 °C for 24 hours in an oil bath.

The resulting polymers were dried at 90 °C for 24 hours, then ground and further processed with a ball mill to obtain fine, homogeneous powders suitable for surface area analysis and SEM imaging. Non-imprinted polymers (NIPs) were synthesised identically, excluding the TNT template.

3.3 Morphological characterisation

To validate the predictions from pre-polymerisation simulations, physical characterisation focused exclusively on morphological parameters: surface area and particle structure.

Nitrogen adsorption–desorption isotherms at 77 K were recorded for all MIP and NIP samples after drying. BET surface areas were calculated from the adsorption branch using the standard Brunauer–Emmett–Teller method, assuming monolayer coverage of N₂ on accessible surface sites.

SEM was employed to visualise the surface texture and morphology of the polymer particles. Samples synthesised at different DMSO : ACN ratios were imaged at magnifications of 300×, 2000×, and 10 000×. The surface topography was analysed to qualitatively assess the degree of porosity and surface roughness, which correlates with the computationally estimated free volume and solvent-accessible surface areas.§

4 Resulting modelling-informed design of porous MIPs

This section presents a simulation-informed strategy to rationally design MIPs with tailored porosity by understanding the impact of solvent composition on monomer organisation and polymer morphology. Rather than assessing adsorption performance or binding specificity, the focus here is on predicting and experimentally validating trends in polymer porosity and surface structure—essential parameters for the development of functional porous materials.

4.1 Template–monomer interactions in pre-polymerisation mixtures

We first examined how solvent composition influences hydrogen bonding between TNT and MAA, the functional monomer. In all simulated systems, TNT's three nitro groups serve as primary hydrogen bond acceptors, forming interactions with the hydroxyl hydrogen of MAA. Due to steric hindrance near the methyl group, the para-nitro group exhibited a slightly higher interaction frequency than the ortho groups (∼4% difference), consistent with prior findings on TNT–MAA complexation.⁹ This preference is illustrated schematically in Fig. 3, where percentages indicate the fraction of simulation time that each nitro group is engaged in hydrogen bonding with MAA (for a representative solvent mixture of 75 : 25 DMSO : ACN). TNT shows a slightly higher tendency to interact via the para-nitro group compared to the ortho-nitro groups, due to steric hindrance affecting the ortho positions. Quantitative hydrogen-bonding results are given in Table 1.

Fig. 3 Chemical structure of the TNT molecule, highlighting the nitro groups in the ortho and para positions (relative to the methyl group).

In addition to nitro–MAA interactions, TNT can, in principle, interact with MAA via its aromatic π system. A previous theoretical study of TNT complexation³⁵ showed that the electron-deficient aromatic ring of TNT can act as a hydrogen-bond acceptor (through its π cloud) or as a π–π stacking partner. In our simulations, we observed occasional π-hydrogen interactions, where the MAA hydroxyl forms a hydrogen bond perpendicular to the face of TNT's aromatic ring. However, these interactions were minor in frequency compared to the dominant nitro-based hydrogen bonds and are not the focus of this study (though they may be of interest for future work on MIP–TNT recognition mechanisms).

As the DMSO content increased, hydrogen bonding between TNT and MAA decreased (Table 1). In pure ACN, template–monomer interactions were sustained for ∼94% of the simulation time, while in pure DMSO this dropped to ∼88%. This trend reflects competition from the porogen itself: DMSO, being a stronger hydrogen bond acceptor than ACN, more effectively interacts with MAA's hydroxyl hydrogen, displacing TNT from potential binding interactions (Fig. 4). This solvent-driven modulation of imprinting interactions confirms previous observations that porogen choice can critically impact MIP formation pathways.^8–10

Fig. 4 Radial distribution functions (RDFs) between the primary hydrogen-bond accepting atom of DMSO (oxygen, blue curve) or ACN (nitrile nitrogen, orange curve) and the carboxylic acid hydrogen of MAA.

Additionally, MAA clustering was observed to increase in DMSO-rich environments. The aggregation of monomer units leads to decreased accessibility for both TNT and solvent molecules. Coordination analysis and RDFs revealed that in high-DMSO formulations, monomers form larger, more internally hydrogen-bonded clusters, consistent with solvophobic effects and poor solvation of MAA in DMSO.¹⁰ This self-association limits effective template complexation, suggesting that DMSO may reduce imprinting efficiency despite improving polymer morphology. To quantify this effect, we computed RDFs between the solvent acceptor atoms and the MAA hydroxyl hydrogen, data shown in Fig. 4 are for the single-solvent systems (100% DMSO and 100% ACN). The peaks near 0.17 nm indicate the presence of hydrogen-bonded solvent–MAA interactions. The larger peak for DMSO reflects stronger competition by DMSO for bonding with MAA relative to ACN. Both pure ACN and pure DMSO systems show a peak in the RDF at around 0.17–0.18 nm, corresponding to hydrogen-bonding distance. However, the peak for DMSO is higher than that for ACN, indicating that DMSO forms stronger/more frequent hydrogen bonds with MAA compared to ACN. This is consistent with DMSO's stronger competitive binding, which in turn explains the reduction in TNT–MAA hydrogen bonding as DMSO content increases.

Fig. 5 shows RDFs between MAA's hydroxyl hydrogen and the acceptor atoms of both ACN and DMSO for each solvent mixture. In both cases, the peak height decreases when going from lower DMSO content (25%) to higher DMSO content (75%), indicating fewer solvent–monomer hydrogen bonds at higher DMSO fractions. As the fraction of ACN decreases (and DMSO increases), the first peak of ACN around MAA diminishes (due to lower ACN concentration and increased competition from DMSO). One might expect the DMSO–MAA RDF peak to grow; accordingly, however, we observe that the DMSO–MAA RDF peak also decreases with higher DMSO content. This somewhat counterintuitive result indicates that beyond a certain concentration, additional DMSO does not lead to more MAA–DMSO interactions, because many of the MAA molecules are now engaged in MAA–MAA interactions (clustered) and thus less accessible even to DMSO. In other words, at high DMSO content, the monomers form self-associated clusters wherein internal MAA molecules are shielded from both TNT and solvent. This highlights DMSO's disruptive role: while it competes strongly with TNT for monomer binding, it simultaneously induces monomer clustering that limits overall hydrogen bonding in the system. The implications of this behaviour on MIP performance (in terms of binding vs. porosity) will be discussed next.

Fig. 5 RDFs of the solvent around MAA for two representatives mixed porogens. (Top) RDF between ACN nitrile nitrogen and MAA hydroxyl hydrogen in solvent mixtures of 75 : 25 (dashed orange line) and 25 : 75 (solid orange line) DMSO : ACN. (Bottom) RDF between DMSO oxygen and MAA hydroxyl hydrogen in solvent mixtures of 75 : 25 (dashed blue line) and 25 : 75 (solid blue line) DMSO : ACN.

4.2 Solvent-driven control of porosity and surface area

The morphological effects of porogen composition were quantified both experimentally and computationally. Nitrogen sorption measurements showed that polymers prepared in pure ACN were essentially nonporous, with BET surface areas below 85 m² g⁻¹. However, the introduction of just 25% DMSO increased surface area nearly fivefold and continued increases in DMSO fraction led to maximum BET values of ∼385 m² g⁻¹ in the 100% DMSO system (Table 2).

Table 2 Experimental and simulated surface areas for the polymers prepared with different porogen compositions. Experimental values are BET surface areas (MIP and NIP averaged) in m² g⁻¹. Computational values are given as relative surface areas from two methods (FreeSASA and MeshSA), normalised to the 100% DMSO case (set as 1.00). For comparison, the experimental surface areas are also expressed as relative values normalised to the 100% DMSO result

Porogen	Exp. SA		Comp. SA
(DMSO : ACN)	(m^2 g^−1)	(rel.)	FreeSASA (rel.)	MeshSA (rel.)
0 : 100	70.98	0.184	0.881	0.950
25 : 75	352.60	0.916	0.957	0.968
50 : 50	371.46	0.965	0.963	0.984
75 : 25	379.77	0.986	0.997	0.989
100 : 0	385.11	1.000	1.000	1.000

---

These trends were accurately predicted by pre-polymerisation molecular simulations. By analysing configurations from MD trajectories with solvent and template removed, two computational approaches—MeshSA (GCMC probing) and FreeSASA (geometric surface analysis)—independently replicated the experimental increase in surface area with rising DMSO content. For porous formulations (≥25% DMSO), simulated relative surface areas were within 4–6% of experimental values (Table 2). This strong agreement confirms the validity of using non-reactive, pre-polymerisation simulations to predict polymer morphology.

One exception was the 0% DMSO case, where simulations overestimated porosity. This discrepancy is likely due to the lack of polymer network collapse in the simulation: after solvent removal, the static polymer configuration is assumed to remain rigid, whereas in reality, polymers formed in poor porogens like ACN collapse into denser structures during curing and drying. Capturing such densification would require reactive MD or post-curing simulations, beyond the current scope.

Further supporting these findings, SEM imaging showed distinct morphological changes with porogen variation (Fig. 6). Polymers synthesized in ACN formed smooth, featureless surfaces, while DMSO-derived samples displayed granular textures with fine surface roughness consistent with microporosity. To further illustrate this solvent-dependent morphological trend, additional SEM micrographs for all intermediate formulations (MIP 2–MIP 4) are provided in Fig. S3 at three magnifications (300×, 2000×, and 10 000×). These images bracket the transition from the smooth morphology of pure ACN-derived MIP 1 to the highly textured surface of pure DMSO-derived MIP 5, enabling visual correlation with the corresponding BET surface area and simulated SASA/FFV values. These visual observations correlate directly with both experimental BET data and the predicted expansion of accessible void space in the simulations (Fig. 7).

Fig. 6 SEM images of polymer particles synthesised in (top row) pure ACN and (bottom row) pure DMSO, at three different magnifications (300×, 2000×, and 10 000×). The ACN-derived polymer forms larger, smoother particles, whereas the DMSO-derived polymer forms smaller, more irregular granules. At high magnification, the surface of the DMSO-derived polymer (b, bottom) shows significant roughness and fine porosity, in contrast to the relatively smooth surface of the ACN polymer (a, bottom).

Fig. 7 Visualization of the computed free volume within the MIP structures (solvent removed) for each solvent composition: 0 : 100, 25 : 75, 50 : 50, 75 : 25, and 100 : 0 DMSO : ACN. These images (from the MeshSA analysis) show regions of accessible void space (coloured surfaces) within the polymer matrix.

Fractional free volume (FFV) calculations from MD trajectories also revealed a linear increase with DMSO content (Table 3), from 0.194 (0% DMSO) to 0.248 (100% DMSO). This trend aligns with the molar volume difference between DMSO (∼71 cm³ mol⁻¹) and ACN (∼53 cm³ mol⁻¹), reinforcing the idea that solvent selection impacts monomer packing density and available void volume in the final material.

Table 3 Fractional free volume (FFV) of the polymer matrices (after removing all solvent molecules), as determined from simulation. FFV is the volume of voids divided by the total volume of the simulation cell

Porogen	FFV
(DMSO : ACN)	(dimensionless)
0 : 100	0.194
25 : 75	0.208
50 : 50	0.222
75 : 25	0.235
100 : 0	0.248

---

Across DMSO fractions, both SASA (FreeSASA/MeshSA) and FFV increase monotonically, mirroring the BET trend (Tables 2 and 3). Together with the hydrogen-bond analysis (Table 1), these data support a solvent-driven trade-off: ACN-rich mixtures favour sustained template–monomer bonding, whereas DMSO-rich mixtures favour enlarged accessible voids.

4.3 Design trade-off: binding interactions vs. porosity

Taken together, these results highlight a critical trade-off in MIP design: acetonitrile-rich mixtures promote stronger template–monomer interactions, suggesting better imprinting fidelity, while DMSO-rich mixtures induce greater porosity and surface area, likely enhancing accessibility and mass transport. An optimal MIP formulation for any application requiring both recognition and adsorption capacity must therefore balance these competing effects.

While this study did not assess rebinding performance or target removal efficiency, the structural insights gained provide a predictive framework for selecting formulations tailored to specific applications. For instance, systems requiring high surface area but tolerating lower imprinting specificity may favour DMSO-rich porogens, whereas selective sensing materials may benefit from ACN-rich environments despite lower porosity.

4.4 Generalisability of the approach

This modelling-experimental strategy for porosity prediction requires only equilibrium MD simulations of pre-polymerisation mixtures and geometric or probe-based surface analysis. No polymerisation modelling, crosslink formation, or adsorption simulations are needed, making this approach computationally tractable and broadly applicable. We anticipate that this method could be extended to other crosslinked systems such as hyper crosslinked polymers, covalent organic frameworks (COFs), and polymers of intrinsic microporosity (PIMs), where solvent choice and monomer arrangement similarly govern final morphology.^15–17

Conclusions

This study presents a computationally informed strategy for the rational design of porous molecularly imprinted polymers (MIPs), demonstrating how pre-polymerisation molecular dynamics (MD) simulations can be used to predict the morphological characteristics of crosslinked polymers prior to synthesis. By systematically varying the composition of DMSO and ACN as porogenic solvents, we explored the interplay between monomer–template interactions and polymer porosity—two key design parameters in MIP formulation.

Simulations revealed that acetonitrile-rich environments enhance specific hydrogen bonding between the TNT template and methacrylic acid (MAA) monomer, whereas DMSO-rich mixtures promote monomer clustering and reduce imprinting efficiency. However, the inclusion of DMSO significantly increased polymer surface area and free volume, as confirmed by both BET measurements and SEM analysis. These experimental trends were quantitatively reproduced by surface area calculations and free volume estimations from solvent-stripped MD configurations, validating the use of MD as a predictive tool for morphology design.

Our approach does not model polymerisation kinetics, curing-induced densification/collapse, or template rebinding; static, solvent-stripped configurations can therefore overestimate porosity in poor-porogen cases (e.g., 0% DMSO), where network densification during curing/drying is expected. Extending the workflow to other chemistries will require force-field validation and descriptor checks tailored to those monomers/porogens; nonetheless, when solvent-driven packing dominates morphology, pre-polymerisation snapshots provide actionable guidance for porosity-directed formulation.

Importantly, our approach requires only simulations of the pre-polymerisation mixture, avoiding the need for polymerisation modelling, reactive force fields, or computationally intensive curing simulations. As such, it provides a tractable and versatile screening method for tailoring the structural properties of porous polymers based on monomer, crosslinker, and solvent composition.

Beyond MIPs, this framework is applicable to a wider range of porous materials, including hyper crosslinked polymers (HCPs), covalent organic frameworks (COFs), and polymers of intrinsic microporosity (PIMs), where solvent-induced morphological control is similarly critical.^21–23 By bridging molecular-scale simulations with experimental synthesis, this work underscores the value of integrative, modelling-guided design in accelerating the development of advanced functional polymeric materials.

Author contributions

Jasmine C. Lightfoot: methodology, formal analysis, investigation, writing – original draft, visualization. William Battell: methodology, formal analysis, investigation, writing – original draft, visualization. Bernardo Castro-Dominguez: supervision, writing – review & editing, funding acquisition. Carmelo Herdes: conceptualization, methodology, project administration, supervision, writing – review & editing, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information: Data supporting this article have been included as part of the SI, which contains individual and comparative nitrogen adsorption isotherms for all MIPs and NIPs from where BET was processed and SEM for all MIPs at three magnification levels. Custom scripts for surface area analysis using FreeSASA and MeshSA are available from the corresponding author upon reasonable request. See DOI: https://doi.org/10.1039/D5ME00102A.

Acknowledgements

The authors gratefully acknowledge the UK Engineering and Physical Sciences Research Council (EPSRC) for funding this work. W. B. acknowledges support from a PhD studentship under grant number EP/R513155/1. This research was also supported in part by the EPSRC grant EP/V051083/1 (Manufacturing in Hospital: BioMed 4.0). Computational resources were provided by the University of Bath's High-Performance Computing (HPC) facility. The authors thank the technical staff at the University of Bath for assistance with nitrogen sorption analysis and scanning electron microscopy.

Notes and references

1. L. Lin, H. Yang and X. Xu, Front. Environ. Sci., 2022, 10, 880246.
2. United Nations, The Sustainable Development Goals Report 2022, United Nations, New York, 2022.
3. J. J. BelBruno, Chem. Rev., 2019, 119, 94–119.
4. L. Chen, X. Wang, W. Lu, X. Wu and J. Li, Chem. Soc. Rev., 2016, 45, 2137–2211.
5. G. Szekely, J. Kupai, M. Razali, S. Buyuktiryaki and R. Kecili, Polym. Chem., 2017, 8, 666–673.
6. J. Gauczinski, Z. Liu, X. Zhang and M. Schönhoff, Langmuir, 2010, 26, 10122–10128.
7. A. Ostovan, M. Arabi, Y. Wang, J. Li, B. Li, X. Wang and L. Chen, Adv. Mater., 2022, 34, 2203154.
8. C. Herdes and L. Sarkisov, Langmuir, 2009, 25, 5352–5359.
9. L. Bird and C. Herdes, Mol. Syst. Des. Eng., 2018, 3, 89–95.
10. W. Battell, G. Donval, B. Castro-Dominguez and C. Herdes, Mol. Syst. Des. Eng., 2022, 7, 196–204.
11. I. A. Nicholls, K. Golker, G. D. Olsson, S. Suriyanarayanan and J. G. Wiklander, Polymer, 2021, 13, 2841.
12. S. Shoravi, G. D. Olsson, B. C. G. Karlsson and I. A. Nicholls, Int. J. Mol. Sci., 2014, 15, 10622–10634.
13. E. Daniels, Y. L. Mustafa, C. Herdes and H. S. Leese, ACS Appl. Bio Mater., 2021, 4, 7243–7253.
14. A. Ostovan, M. Arabi, Y. Wang, J. Li, B. Li, X. Wang and L. Chen, Adv. Mater., 2022, 34, 2203154.
15. F. Ahmadi, E. Yawari and M. Nikbakht, J. Chromatogr. A, 2014, 1338, 9–16.
16. K. M. Muzyka, J. Nano- Electron. Phys., 2015, 7, 01017.
17. M. Perez, R. Concu, M. Ornelas, M. N. D. S. Cordeiro, M. Azenha and A. F. Silva, Spectrochim. Acta, Part A, 2016, 153, 661–668.
18. M. Naderi, M. Ghanei, M. Shohrati, A. Saburi, M. Babaei and B. Najafian, Cutaneous Ocul. Toxicol., 2013, 32, 31–34.
19. United Nations Convention to Combat Desertification, Reaping the Rewards: Financing Land Degradation Neutrality, United Nations, Bonn, 2015.
20. S. Xu, H. Lu, J. Li, X. Song, A. Wang, L. Chen and S. Han, ACS Appl. Mater. Interfaces, 2013, 5, 8146–8154.
21. K. Cho, C. W. Kang, S. H. Ryu, J. Y. Jang and S. U. Son, J. Mater. Chem. A, 2022, 10, 6950–6964.
22. M. Khakbaz, A. Ghaemi and G. M. M. Sadeghi, Polym. Chem., 2021, 12, 6962–6997.
23. N. B. McKeown and P. M. Budd, Chem. Soc. Rev., 2006, 35, 675–683.
24. M. Piccini, J. Lightfoot, B. Castro-Dominguez and A. Buchard, ACS Appl. Polym. Mater., 2021, 3, 5870–5881.
25. J. C. Lightfoot, A. Buchard, B. Castro-Dominguez and S. C. Parker, Macromolecules, 2022, 55, 498–510.
26. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
27. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
28. S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess and E. Lindahl, Bioinformatics, 2013, 29, 845–854.
29. A. K. Malde, L. Zuo, M. Breeze, M. Stroet, D. Poger, P. C. Nair, C. Oostenbrink and A. E. Mark, J. Chem. Theory Comput., 2011, 7, 4026–4037.
30. J. A. Purton, J. C. Crabtree and S. C. Parker, Mol. Simul., 2013, 39, 1240–1252.
31. A. V. Brukhno, J. Grant, T. L. Underwood, K. Stratford, S. C. Parker, J. A. Purton and N. B. Wilding, Mol. Simul., 2021, 47, 131–151.
32. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
33. S. Mitternacht, F1000Research, 2016, 5, 189.
34. S. A. Ali, M. I. Hassan, A. Islam and F. Ahmad, Curr. Protein Pept. Sci., 2014, 15, 456–476.
35. M.-Z. Ao, Z. Tao, H.-X. Liu, D.-Y. Wu and X. Wang, Comput. Theor. Chem., 2015, 1064, 25–34.

---

Footnotes

† J. C. L. and W. B. contributed equally to this work.

‡ Current address: The Hartree Centre, STFC Daresbury Laboratory, Warrington, WA4 4AD, UK.

§ TNT was extracted from the MIPs following synthesis using Soxhlet extraction with refluxing methanol. While this ensured template removal for surface area consistency, no rebinding experiments were conducted, and no adsorption performance metrics were evaluated in this study.

---