Anticooperative self-assembly of perylene diimide dyes in water unveiled by advanced molecular dynamics simulations

Abstract

The self-assembly mechanism and the dynamics in water of π-stacked columnar aggregates of N,N-bis(2-(trimethylammonium)ethylene)perylene-3,4,9,10-tetracarboxylic acid bis-imide (PDI) have been investigated using classical molecular dynamics (MD) simulations based on an accurate quantum mechanically derived force field. Nano-scale aggregates composed of up to twelve monomers were studied through a combination of bottom-up and top-down approaches, complemented by enhanced sampling MD techniques. Our simulations provide unprecedented atomistic information on the anti-cooperative aggregation mechanism, consistent with experimental findings for PDIs with bulky imide substituents. Umbrella sampling simulations coupled with in-depth thermodynamic analysis reveal, indeed, that entropic contributions dominate the aggregation process, resulting in a non-linear growth mechanism that saturates with increasing aggregate size. Replica exchange MD runs confirm that experimental properties of these systems can be effectively described considering species up to tetramers or hexamers. Although larger aggregates may form, they represent only a minor fraction of the observed species. The stability and dynamics of PDI aggregates in water were further analyzed by means of specific geometric supramolecular descriptors designed to capture relative monomer motions. This analysis indicates that dimers, rather than monomers, act as the primary building blocks, with even-numbered aggregates behaving like strongly-bound interacting dimers. Monitoring water structure around PDI aggregates unveils a rather complex scenario, where the hydration of charged side chains is persistent across all aggregate sizes, while steric shielding within π–π stacks progressively reduces core solvation, though carbonyl oxygens maintain hydrogen-bond interactions with water. Analysis of hydrogen-bond statistics confirms a gradual decrease in solvation per monomer with aggregation, reaching a plateau at larger stack sizes. Complementary dynamical investigations demonstrated that water expulsion does not follow a simple concerted process, but instead proceeds through partially solvated intermediates reflecting the energetic cost of disrupting interfacial hydrogen bonding. Together, these findings highlight the critical entropic role of solvent reorganization in driving aggregation and may provide mechanistic insights relevant for optimizing photocatalytic activity through rational monomer design.

---

1. Introduction

Over the past few decades, the self-assembly of polycyclic aromatic molecules via non-covalent interactions to form functional supramolecular nanostructures has opened new and promising avenues in the field of materials chemistry. In this context, perylene diimide (PDI) and its derivatives are widely employed^1–7 due to their excellent chemical stability, intense optical absorption, and efficient charge transport capabilities.^8–11 PDI-based materials have found successful technological applications in biomedicine,^12–14 as well as in organic electronics,^12,15–19 including organic light-emitting diodes (OLEDs),²⁰ organic field-effect transistors (OFETs),²¹ organic photovoltaic (OPV) cells,^22,23 and photocatalysis.^24–26 The great interest in PDI aggregates and their unique characteristics stems from their stiff planar backbone and extended conjugated core, which strongly favor intermolecular π–π interactions, leading to the formation of columnar-like one-dimensional nanostructures.^3,27 The optoelectronic and packing properties of such PDI aggregates can be finely tuned by functionalization at the imide, bay, and ortho carbon atom positions.^28–30 Specifically, imide substitutions influence self-assembly behavior independently of the intrinsic optoelectronic properties of the perylene core, providing an effective strategy for controlling the morphology, the type (J, H) and the mechanisms of aggregation. By exploiting covalent linking as well as a complex interplay of hydrophilic, hydrophobic, π–π and hydrogen bond (HB) interactions and metal-coordination, indeed, a large variety of supramolecular PDI aggregates, with peculiar structural, optical, redox and charge transport properties, have been reported in the literature.^5,7,31–45 In a more general context, three principal models have been proposed to describe supramolecular growth mechanisms: (i) isodesmic, where the binding strength remains constant regardless of aggregate size (i.e., nucleation and elongation equilibria are indistinguishable); (ii) cooperative, where the binding strength increases as the aggregate size increases (i.e., elongation equilibria are more favorable than nucleation equilibria); and (iii) anti-cooperative, where the binding strength decreases as the aggregate size increases (i.e., elongation is less thermodynamically favored than nucleation).⁴⁶

Recent works have revealed that achieving fine-tuned morphological, optoelectronic, and charge transport properties requires coupling molecular design of standalone dyes with precise knowledge of the self-assembly mechanism as a function of temperature and environmental factors, such as solvent mixture and pH.^47–51 From this perspective, it is crucial to achieve a molecular-level understanding of the driving forces that govern self-assembly. In this context, the ability of atomistic simulations to precisely characterize the geometrical arrangements of individual monomeric units within a supramolecular system represents an important complement to experimental techniques. These latter often lack spatial, temporal, and energetic resolution, particularly when it comes to resolving binding energies and distinguishing the subtle energetic interplay between different molecular interaction types.^52–54 Additionally, the collection of configurational ensembles, capable of representing with a reliable statistic the detailed structure of such aggregates, also paves the way to an accurate modeling of their optoelectronic properties.⁵⁵ Computational simulations, and in particular classical molecular dynamics (MD), have been extensively applied to study molecular packing in solid state and in solution, providing atomistic insights into dynamics, disorder, and nonequilibrium phenomena at the nanoscale.^56–58 Due to the significant interest in PDI-based materials, several MD-based investigations, often combined with quantum chemical calculations and experimental measurements, have been conducted to extract and predict the structural, thermodynamic, and optoelectronic properties of their supramolecular aggregates.^{47,55,59–75} The reliability of MD simulations, however, is grounded in the accuracy of the employed interatomic potentials, i.e., the force fields (FFs). While standard FFs are usually considered accurate for large biomolecules,⁷⁶ most MD studies on supramolecular materials typically reoptimize an initial set of bonded and nonbonded parameters extracted from general-purpose transferable FFs (e.g., OPLS-AA, GAFF, GROMOS, CHARMM).^56,77–83 A robust approach to improve FF accuracy is to rely on quantum mechanical (QM) calculations to parametrize the classical FF model potential terms.^81,84–86 This strategy, while lacking in transferability, offers notable advantages in terms of predictivity, chemical specificity, and reliability under different thermodynamic conditions. These advantages are significant when one aims to rationalize specific structure–property relationships and predict supramolecular architectures based on the monomeric molecular structure. Nonetheless, such Quantum Mechanically Derived Force Fields (QMD-FFs)^81,84 have been mainly applied on homogeneous assemblies,^53,81–83 where supramolecular aggregates were built from a random arrangement of identical molecular units.

When simulating self-assembling processes in solution, two approaches can be employed: the aforementioned bottom-up strategy, where the aggregate spontaneously assembles from a random arrangement of monomers, or a top-down scheme, where pre-stacked configurations are optimized, and the packing dynamics are monitored. To achieve a robust assessment of the stability and properties of supramolecular aggregates in solution, both protocols should be followed, integrating the information retrieved from each scheme. In fact, while the top-down approach allows for the study of realistic, large-sized aggregates, it may be biased by the chosen initial conditions, partially trapping the system in metastable states, due to insufficient sampling of the conformational space, which is related to ergodic limitations. On the other side, a bottom-up strategy may not guarantee to achieve a complete self-assembly of large aggregates even after a simulation time of hundreds of nanoseconds, due to similar constraints on sampling. For this reason, enhanced sampling techniques (e.g., umbrella sampling, steered MD, and replica exchange MD) are often employed.⁸⁷

Here we report a detailed investigation of N,N-bis(2-(trimethylammonium)ethylene) perylene 3,4,9,10-tetracarboxylic acid bis-imide (hereafter termed PDI, see Fig. 1) self-assembly in water, providing insights into the formation and dynamics of their 1D columnar aggregates, and thoroughly analyzing their interaction with the solution environment. We rely on extensive MD simulations, considering 3, 4, 8, and 12 monomeric units, based on an accurate QMD-FF^84,86 previously developed by some of us and applied to study the optical properties of PDI monomers and dimers.^55,88 A bottom-up approach was adopted for the smallest systems, while results for pre-stacked (top-down scheme) 8- and 12-monomeric systems were complemented with enhanced sampling simulations. Umbrella sampling calculations and a detailed thermodynamic analysis revealed that the aggregation process is entropically driven and follows a non-linear trend, with saturation observed as the aggregate size increases. This behavior is indicative of an anti-cooperative mechanism that disfavors the formation of large supramolecular assemblies. This was confirmed by additional replica exchange MD runs, which indicate that only relatively small aggregates, composed by up to seven units, are likely to self-assemble, whereas larger structures were less likely to form.

Fig. 1 N,N-Bis(2-(trimethylammonium)ethylene) perylene 3,4,9,10-tetracarboxylic acid bis-imide (PDI).

Detailed information on the dynamics of the PDI units within the supramolecular systems, as well as on their collective dynamics, was obtained defining a comprehensive set of geometrical supramolecular descriptors and analyzing their evolution and correlations throughout MD simulations of different aggregate models. Note that similar approaches, although based on smaller sets of geometrical parameters, have been previously proposed to correlate optical and charge transport properties to structural features of π–π stacks of PDI.^66,69 The proposed structural analysis protocol utilizes geometrical parameters to precisely track monomer motion within aggregates. When complemented by electronic structure, excited state dynamics, and optical property calculations, it can provide detailed information on the structure–property relationships in supramolecular systems.

The paper is structured as follows: Section 2 outlines the computational details, Section 3 presents and discusses the results, and Section 4 summarizes our findings. We highlight an extended set of complementary data described in the SI, containing additional details about supramolecular descriptors, aggregation theromdynamics and the role of solvent structure.

2. Computational details

2.1. Aggregate models and QMD-FF

We considered aggregates models composed of three (PDI₃), four (PDI₄), eight (PDI₈) and twelve (PDI₁₂) PDI molecules in water, hereafter indicated as PDI_n@H₂O and shown in Fig. S1. A bottom-up approach was used to model self-assembly of the smaller species (PDI₃ and PDI₄), by starting MD simulations from monomers being randomly distributed in solvated boxes of 45³ ų, containing 2976 and 2942 water molecules, respectively. For the larger systems, as we could not observe complete self-assembly within our simulation times (500 ns), and we could not establish whether this observation describes the physical situation in detail or is affected by the limited simulation time, we started from π-stacked columnar aggregate models,⁵⁹ in boxes of 95³ ų and 120³ ų filled with 32296, for PDI₈, and 56 358 water molecules, for PDI₁₂. To maintain charge neutrality, 2 Cl⁻ counterions were added for each PDI; this choice assures internal coherence with MD simulations previously performed for the PDI monomer⁸⁸ and dimer (PDI₂).⁵⁵ The MD trajectories on PDI₂, previously studied in ref. 55, are re-analyzed here, following the new protocol and are used to evaluate the effect of system size on supramolecular dynamics. To observe spontaneous aggregation of the larger PDI₈ and PDI₁₂ aggregates, we instead decided to resort to advanced sampling methods (see Section 2.3.2).

All intramolecular QMD-FF parameters employed here were obtained following the JOYCE protocol^84,86 and DFT data previously obtained for the PDI monomer, and validated against experimental spectra recorded in dilute solutions.⁸⁸ For intermolecular terms, point charges were derived from the same DFT density used for the JOYCE procedure and obtained using the RESP protocol by the Antechamber suite,⁸⁹ while Lennard–Jones (LJ) parameters were transferred from OPLS.⁹⁰ The accuracy of intermolecular interactions was assessed in a previous work on the dimer,⁵⁵ through comparison of stacking energy profiles with ab initio reference data, as well as by successfully reproducing the experimental spectra of both dilute and concentrated PDI solutions. Water molecules were simulated with the TIP3P model, while Cl⁻ parameters were taken from the OPLS force field.⁹⁰

2.2. Molecular dynamics simulations

MD simulations were performed under periodic boundary conditions using the GROMACS 5.1 engine.⁹¹ The PDI_n@H₂O systems were first minimized to avoid possible bad contacts and then thermally equilibrated for 2 ns at 300 K in the NVT ensemble. We then performed equilibration steps in the NPT ensemble for 1 ns at 1 atm and 300 K with a time step of 1 fs, using the LINCS⁹² algorithm to fix bond distances involving hydrogen atoms. Subsequently, the production runs were computed differently depending on the size of the aggregate. MD production runs for PDI₂@H₂O, PDI₃@H₂O and PDI₄@@H₂O were carried out for 1 μs with a timestep of 1 fs at 300 K of temperature and 1 bar pressure. Conversely, for the larger PDI₈@H₂O and PDI₁₂@H₂O systems, 500 ns trajectories with a timestep of 1 fs were carried out. At the beginning of each production run, no specific conformational constraint was imposed on the PDI units, except for bond distances as in the equilibration steps. Temperature and pressure coupling were described using the Parrinello–Rahman and v-rescale schemes, with coupling constants of 0.1 and 1 ps. For short-range interactions, a cutoff radius of 11 Å was set up, while long-range electrostatics was accounted through the particle mesh Ewald (PME) procedure. These settings were maintained in all MD runs, including enhanced sampling simulations described in Section 2.3.

To conduct an in-depth structural analysis of the MD trajectories for the smaller systems, we extracted 50 000 conformations equally spaced in time, with a frequency interval of 10 ps, from the last 500 ns of the trajectory. Only 5000 conformations were instead considered for the octamer and dodecamer, whose trajectories were sampled considering intervals of 100 ps. Both sampling frequencies allowed us to get a detailed and accurate representation of the motion of the PDI units and of the interactions occurring within the aggregates. Finally, we also investigated the solvent structure around the columnar aggregates by computing the radial distribution function (RDF) g(r), between selected atoms of the solute and the surrounding solvent molecules.⁹³ Additionally, we plotted the three-dimensional solvent density distribution (SDF) around relevant regions of the solvated aggregates, i.e. either with respect to the PDI scaffold, or around the flexible chains. In addition to RDF and SDF analyses, we quantified the number of HBs formed between water molecules and the PDI carbonyl oxygens as well as the side-chain nitrogens, also monitoring the number of water molecules confined between PDI cores during the self-aggregation process. These descriptors capture how solute–solvent interactions evolve upon aggregation, with full methodological details provided in Section D of the SI, while the corresponding results are discussed in Sections 3.1 and 3.3.

2.3. Enhanced sampling

2.3.1. Umbrella sampling simulations. The PDI self-assembly mechanism was investigated using umbrella sampling (US) simulations.⁹⁴ The systems were initially energy-minimized and equilibrated using the NVT and NPT ensembles to ensure stability. To this aim, two to six PDI units were solvated and equilibrated at atmospheric pressure and room temperature through consecutive energy minimizations, NVT and NPT runs. The equilibrated systems have box dimensions spacing from ca. 40 × 40 × 60 ų, for the 2 PDI unit system, up to 40 × 40 × 80 ų, for the 6 PDI unit systems. The umbrella sampling was applied to each equilibrated system, using the distance between centres of mass (COM) of two cofacial PDI monomers as reaction coordinate, applying a harmonic potential with force constant of 1000 kJ mol⁻¹ nm⁻² as a constraint. To ensure proper sampling, independent umbrella sampling runs were carried out for 10 ns in 200 configurations in 25 equidistant windows, spanning the reaction coordinate range from 2.30 Å to 15 Å. All such calculations were run at constant pressure and temperature. The final potential of mean force (PMF) profiles were computed using the Weighted Histogram Analysis Method (WHAM), as implemented in the GROMACS software package⁹¹ by the wham module. The uncertainty related to the PMF calculations was evaluated using a Bayesian block bootstrapping technique,⁹⁵ performing 500 bootstrapping iterations to ensure robust statistical analysis.

2.3.2. HREMD simulations. To investigate spontaneous formation of larger stacks, we exploited advanced sampling simulations, similarly to what recently reported in literature to investigate aggregates of similar organic species.^96,97 Concretely, we resorted to Hamiltonian Replica Exchange MD (HREMD), also carried out with the GROMACS software,⁹⁸ extending its capabilities through the open–source, community–developed PLUMED library.^99–101 HREMD runs were performed only for a system constituted by 48 PDI monomer units in a water box. In HREMD simulations, non–bonded (LJ and charges) and torsional parameters of a (sub)system are down–scaled in the various replicas, hence decreasing energy barriers separating local minima, allowing for an exploration of wider regions of the free energy surface. Roughly speaking, a “selective heating” is applied to a specific portion of the system that, in our case, is the solute (PDI units). A swap between contiguous replicas is possible when their energies are close enough, and a trajectory collects all visited nuclear arrangements, which come from the different sets of parameters used in all the replicas. We used 16 replicas, with effective temperatures between 300 K and 800 K, corresponding to (the scaling factor for non–bonded and torsional parameters) between 0.375 and 1, following a geometric distribution.⁹⁸ For each replica, we ran a 465 ns simulation, totaling 7.44 μs, and obtaining an exchange probability between contiguous replicas between 10% and 40%, as recommended in literature.^102,103 Analysis of these simulations was eventually carried out identifying pairs of PDI units whose centers of mass are within 5 Å from one another, and grouping them into clusters of increasing dimensions.

2.4. Aggregation thermodynamics

To gain a deeper insight into the self-assembly process, we further analyzed the US trajectories focusing on the thermodynamics subtending the stack elongation. In this framework, the aggregation enthalpy ΔH_aggr was evaluated as the enthalpy difference between the elongated stack and a matched reference, as discussed in detail in Section D of the SI. The corresponding entropy term was then obtained from

TΔS_aggr = ΔH_aggr − ΔG,

where ΔG is the binding free energy derived from the US free-energy profiles. To rationalize the origin of the enthalpic trend, we further analysed it in terms of the stepwise process that leads to the stack formation (see eqn (S3) and (S4)). To this end, we decomposed the stepwise potential-energy change into three contributions (Section D.2, eqn (S5)–(S7)): (i) ΔE^PDI–PDI_aggr, describing intra-stack interactions between chromophores, including both stabilizing π–π contacts and less favorable side-chain contacts; (ii) ΔE^wat−wat_aggr, associated with solvent–solvent interactions, primarily reflecting the HBs formed among water molecules released into the bulk; and (iii) ΔE^PDI−wat_aggr, quantifying the change in solute–solvent interactions and thus the loss of hydration upon stacking. In qualitative terms, both ΔE^PDI−PDI_aggr and ΔE^wat−wat_aggr are stabilizing, whereas ΔE^PDI−wat_aggr is strongly positive, reflecting the energetic penalty of reduced solvation as discussed in Section 3.1 in detail. Taken together, these opposing contributions establish the enthalpic profile of column elongation, where the driving force arises from a delicate balance between intra-stack stabilization and solvent reorganization on the one hand, and the loss of PDI–water interactions on the other.

2.5. Descriptors of the supramolecular structure

An effective way to analyze the supramolecular structure and dynamics of a PDI aggregate is in terms of geometrical descriptors.^55,73,75 As described in the following, in this work we integrate the set of descriptors previously adopted by some of us⁵⁵ for PDI₂ with additional parameters to better describe the self-assembly of larger aggregates.

As shown in (a) and (b) panels of Fig. 2, a molecular reference system is first assigned to each PDI unit k. The first two axes, û_k,∥ and û_k,s, both lay on the aromatic scaffold, respectively defining (see also Fig. S2 in SI) the PDI's long (∥) and short (s) axes, while the third one, û_k,⊥, is defined as

û_k,⊥ = û_k,∥ × û_k,s (1)

and is hence perpendicular to the aromatic plane. Together with the center of mass position of each unit (), the three reference axes û_k,x (x = ∥, s, ⊥) were employed to define a collection of supramolecular descriptors, aimed to analyze the structure of large PDI aggregates of N_PDI (N_PDI = 2, 4, 8, 12), in terms of dimer contributions. Hence, as displayed in Fig. 2c and d, four positional (, , and ) and three orientational supramolecular descriptors (α_kl, γ_kl and ) were defined for each dimer pair k,l within the aggregate, as detailed in the following. Within this framework, the distance vector between two PDI units k and l is computed as

Fig. 2 Definition of the geometrical descriptors employed to describe PDI's supramolecular assembly. In top panels, the molecular reference systems (û_∥ and û_s, û_⊥, see also Fig. S2 in the SI) are shown for PDI units k and l (green): (a) side view, (b) top view. In central panels, selected supramolecular descriptors are shown for the same units. (c) Inter-unit vector (), stacking vector (), rolling (γ_kl) and tumbling angles (). (d) sliding () and shifting () vectors and spinning angle (α_kl). Bottom panels display the connection between each descriptor and the positional, (e), and orientational, (f), supramolecular displacements (see also Fig. S3–S9 for more details).

As shown in Fig. 2e, to separately follow the stacking, sliding and shifting motions of unit l with respect to k (see also Fig. S3–S5), the displacements along the three molecular reference axes are respectively defined as:

and

Similarly, as displayed in Fig. 2f and reported in more detail in Fig. S6–S8 in the SI, the rotational motion of each PDI unit with respect to another can be decomposed in terms of spinning, rolling and tumbling, and thus described by the yaw, roll and pitch angles computed from the molecular descriptors as:

α_kl = arccos (û_k,∥·û_l,∥) (6)

γ_kl = arccos (û_k,⊥·û_l,s) (7)

and

Finally, to obtain a deeper insight on the relation of the aggregate size with its stability, two additional descriptors were included in the set: the columnar vector and, for each PDI unit k within the aggregate, the dephasing angle δ_k. As shown in Fig. S9 in the SI, the former can be defined for each considered aggregate as

while the latter is computed for each unit k as

For the sake of clarity, Table 1 provides a summary of the geometrical descriptors discussed thus far and illustrated in Fig. 2 and S3–S8 in SI, along with their reference values obtained by averaging over 5000 snapshots from the PDI₂@H₂O trajectory.⁵⁵

Table 1 List of the employed structural descriptors and their average values computed for the a and b PDI units along the PDI₂@H₂O trajectory⁵⁵

Descriptor	Reference value
ρ ab (Å)	4.1 ± 0.3
R ^πab (Å)	2.9 ± 0.5
R ^∥ab (Å)	2.0 ± 0.7
R ^Sab (Å)	1.8 ± 0.7
α ab (°)	32 ± 7
γ ab (°)	90 ± 5
(°)	90 ± 3
δ b (°)	25 ± 10

---

3. Results and discussion

3.1. PDI self-assembly

We first investigate the mechanism of PDI π-stacked self-assembly in water by performing MD runs on smaller systems (PDI₂@H₂O, PDI₃@H₂O and PDI₄@H₂O), whose monomers were randomly distributed in the simulation box and thereafter solvated with ∼3 × 10³ water molecules, corresponding to concentrations (of the species PDI) of 0.036 M, 0.054 M, and 0.072 M, respectively.⁵⁵

Aggregation was monitored evaluating the time evolution of the intermolecular distance defined in eqn (2) (see Fig. 2), hereafter referred to simply as ρ. The data shown in Fig. 3 indicate for all systems the sequential formation of stable π-aggregates within the considered timescale (50 ns). The intermolecular distance between two neighboring units stabilizes in the range of 3.7 Å < ρ < 3.9 Å, consistently with the expected experimental value for π-stacked systems.^104,105 The top panel of Fig. 3a indicates that, in PDI₂@H₂O, a dimer self-assembles within 1 ns. In line with results of previous MD simulations,⁶¹ where the spontaneous formation of trimeric species was observed within 50 ns, addition of the third monomer leads to stepwise aggregation, with the dimer assembling in less than 2 ns and the third PDI unit joining the aggregate after more than 30 ns (Fig. 3a, central panel). As expected considering the increased concentration with respect to the PDI₃@H₂O system, tetramer formation in PDI₄@H₂O is observed within 26 ns (see Fig. 3a, bottom panel). The latter exhibits the same sequential assembly pathway, in which a dimer forms within 1 ns, followed by the formation of a stable trimer after approximately 7 ns. This result should hint at the high chance that self-aggregation is occurring.

Fig. 3 (a) Time evolution of ρ (Å) along the first 50 ns of MD production runs for PDI₂@H₂O (top), PDI₃@H₂O (middle), and PDI₄@H₂O (bottom). (b) Graphical representation of ρ (Å), corresponding to the intermolecular distance between geometrical centers O_k and O_l of two adjacent monomers (see eqn (2)).

As described in the Methods Section, to tackle the self-assembly of larger aggregates we resorted to enhanced sampling techniques. First, we carried out US calculations, summarized in Fig. 4a, which clearly indicate that the PDI self-assembly mechanism exhibits saturation, where the binding free energy ΔG stabilizes at an energetic plateau as the size of the aggregates increases (see also Table 2). This finding points to a reduction of the probability of the formation of very large aggregates, instead suggesting that self-assembly occurs through an anti-cooperative mechanism, as experimentally observed in cases where PDI molecules bear sterically demanding substituents or when strong H-bond interactions significantly favor dimer formation.¹⁰⁶ Indeed, the values gathered in Table S1 in the SI, confirm that the dimer is the species showing the strongest interaction, with a calculated binding free energy ΔG of −46 kJ mol⁻¹, in line with results of previous reports.^62,107 As the number of PDI units in the aggregate increases, ΔG steadily decreases (absolute terms), as expected for bulky-substituted PDIs.^108,109 For example, ΔG between a monomer aggregating with a pentamer drops to −35.5 kJ mol⁻¹ (see Table S1 and Fig. 4a). As discussed in ref. 106, when dimerization is favored over stepwise monomer addition in the formation of larger aggregates, a predominance of even-numbered aggregates is observed, as the dimer serves as the primary building block. The calculated binding free energies listed in Table S1 support this dimer-based aggregation mechanism. Specifically, PDI₄ formation is favored through binding of two dimers (−40.0 kJ mol⁻¹) compared to addition of a monomer to PDI₃ (−38.0 kJ mol⁻¹). Similarly, hexamer formation is more favorable when PDI₂ is added to PDI₄, compared to the aggregation of monomer + pentamer or PDI₃ + PDI₃.

Fig. 4 (a) Computed binding free energy (ΔG, kJ mol⁻¹) with error bars as a function of the number of assembled PDI monomers (data in Table S1). The black line represents an exponential fit of binding energies, yielding a saturation ΔG of −34.3 kJ mol⁻¹ (dashed black line). (b) Fraction of species identified in HREMD simulations (in blue) and fraction of snapshots considering only the largest aggregate (in red).

Table 2 Aggregation enthalpy and entropy contributions to the total binding free energy ΔG, computed for the self-assembled dimer, tetramer, and hexamer aggregates. All values are reported in kJ mol⁻¹. For tetramer and hexamer, both ΔH and ΔG were obtained by considering the stepwise aggregation pathway derived from US MD simulations, as detailed in Table S1. The uncertainty on TΔS_aggr was estimated by summing the uncertainties on ΔH_aggr and ΔG

n	ΔHaggr	TΔSaggr	ΔG
2	40 ± 32	86 ± 34	−46 ± 2
4	72 ± 40	110 ± 42	−38 ± 2
6	86 ± 46	122 ± 48	−36 ± 2

---

The thermodynamic analysis of the US trajectories, reported in detail in Section D of the SI and summarized in Table 2, shows that aggregation is enthalpically disfavored (ΔH_aggr > 0) but proceeds thanks to a favorable entropy term (TΔS_aggr > 0), mainly coming from the release of solvating water molecules into the bulk.¹¹⁰ For the sake of simplicity, we focused on stepwise N → N + 1 aggregation events (dimer, tetramer, hexamer), defined as the transformation of a reference system consisting of an N-mer and a solvated monomer (S_N;M) into the elongated stack S_N+1. In practice, this corresponds to the addition of one monomer to a pre-formed stack or, in the case of the dimer, to the association of two monomers. This approach provides a qualitative yet representative picture of the energetics governing column elongation, capturing the essential balance of interactions even if not encompassing all possible assembly pathways. As the stack grows, both the enthalpic penalty and the entropy gain tend towards a plateau: the former reflects the saturation of the interaction–energy balance (see Table S7), while the latter mirrors the progressively smaller number of interfacial waters released into the bulk (Fig. S47 and S48). As a consequence, ΔG becomes nearly size-independent, consistently with the free-energy saturation reported in Fig. 4 (fit to ∼–34.3 kJ mol⁻¹). This leveling-off implies that the driving force for elongation weakens with aggregate size, rationalizing the anti-cooperative nature of the self-assembly and the limited extension of the columnar stacks.

As the column length increases, each PDI unit establishes progressively fewer HBs with the surrounding water molecules, and this reduced hydration correlates with the increase of the PDI–water contribution to the aggregation energy, which rises from about 88 kJ mol⁻¹ in the dimer to nearly 238 kJ mol⁻¹ in the hexamer. Such contribution is directly connected to the progressive loss of hydration upon stacking, and is decisive in determining the positive sign of the total potential energy balance reported in eqn (S7). The effect of the reduced solvation upon aggregation therefore has a two-fold impact on the free energy of the assembly process: on the one hand, the reduction of PDI–water interactions reflects the increasingly unfavorable enthalpic contribution, on the other hand, with fewer interfacial waters available, the entropy gain from their release into bulk also decreases as the column elongates.

To confirm the US results, we complemented the above information with HREMD simulations. The starting configuration consisted in a previously equilibrated system containing 48 solvated PDI units, which could be identified as 13 monomers, one self-assembled trimer and 16 dimers. Analysis of HREMD trajectories from the 16 replica simulations reveals the formation of aggregates containing up to 10 π-stacked PDI units. Although it is possible that increasing solute concentration in the simulation box, the length of these simulations, or sampling lower λ values, even larger aggregates could be obtained, it is clear from Fig. 4b that the probability of formation of decamers is very low (below ≈0.05% of the clusters observed). Furthermore, species up to hexamers account for ≈99% of all species observed in solution, with a majority of monomers (≈37%), dimers (≈30%), trimers (≈20%) and tetramers (≈10%). Examples of the largest aggregates obtained are displayed in Fig. S10 in the SI. Since these results suggest that observed experimental properties should be mainly attributed to smaller assemblies, the search for larger aggregates was deemed as not necessary, hence avoiding a considerable computational burden. Fig. 4b also displays the largest aggregate formed in each snapshot. In this case, we see that there is essentially no situation in which only monomers or dimers are present. In ≈37% of the cases, the predominant species is the tetramer, followed by cases in which pentamers (≈32%), and hexamers (≈15%) are the largest species present, hence confirming the lower probability of observing larger aggregates such as heptamers (≈9%) or octamers (≈1%). In conclusion, the enhanced exploration of the system configurational space confirms that the system spontaneously self-assembles, yet through a dimer-based anti-cooperative mechanism, which disfavors the formation of large columnar aggregates, likely limiting the size of the stacked supramolecular structure to less than ten units. All simulations were carried out under standard conditions (1 atm, 300 K), and it remains to be seen whether the observed anti-cooperative mechanism would hold under different thermodynamic regimes. Investigating the role of temperature may offer insight into potential changes in aggregation behavior—a question that could be explored in future studies.

3.2. Structure and dynamics of aggregates in water

The conformational dynamics of a pair of consecutive PDI units within aggregates of varying sizes can serve as a further indicator of the strength of π–π stacking, which directly affects the material's optical and photophysical properties.⁵⁵ We first examine the population of two of the geometrical descriptors displayed in Fig. 2 along MD runs, namely the distance between the centroids of two adjacent monomers, ρ, and its projection along the stacking axis, R_π. We performed this analysis for all neighboring pairs within each columnar aggregate considered in this work. The complete set of data is reported in Fig. S12–S14 in the SI, while Fig. 5 highlights only the most relevant findings. Specifically, in Fig. 5, the left panels show the ρ and R^π probability distributions for the PDI₂@H₂O and PDI₄@H₂O trajectories, representing the dimer and selected PDI pairs within the tetramer, respectively. The right panels display the corresponding 2D heat maps, revealing correlations between these parameters. R^π distribution peaks at approximately 2.6 Å, and its shape remains similar in all pairs, showing a significantly more pronounced broadening with respect to ρ, with a long and pronounced tail at larger distances. Consistently with the values reported for similar π–π stacked aggregates, the most probable distance ρ is ca. 3.8 Å. Interestingly, at difference with R^π, a noticeable shape variation in ρ distributions is evident when moving from the dimer (Fig. 5a) to the tetramer pairs (Fig. 5b–d). Visual inspection is supported by analysis of the first moment (and standard deviation) of the ρ distributions (see Table S4 in the SI) which are 4.12 (0.34) Å for PDI₂@H₂O (Fig. 5a), and 3.87 (0.26), 3.95 (0.32), and 3.90 (0.28) Å, for the three dimers in PDI₄@H₂O in Fig. 5b–d respectively. In particular, the larger first moment and broader signal obtained for PDI₂@H₂O indicates that units within the isolated dimer experience greater relative mobility compared to those in larger aggregates. Although this could be expected for the central pair, it is noticeable that even when positioned at the ends of a tetrameric column, the conformational dynamics of such PDI pairs is also reduced with respect to PDI₂@H₂O, thus indicating a stabilizing effect from the extended π-stacked environment. This trend is further supported by Fig. S12–S14 in the SI and also clearly visible in the 2D heat map displayed in the right panel of Fig. 5a, which shows a pronounced asymmetrical distribution and regions of low-to-medium color intensity. On the contrary, for dimeric pairs within the tetramer, high probability (red) areas appear, again pointing to a more constrained dynamics in the larger aggregate. It is also interesting to examine the slight differences between different pairs, separately displayed in the (b) to (d) panels. In agreement with ΔG calculations discussed in the previous section, the tetramer appears to be constituted by two tightly bound dimers (a–b and c–d), which exhibit similar ρ distributions (residual differences, e.g. in first moment and standard deviation, are only due to a non perfect statistical sampling) and heat maps, indicating relatively stronger stacking compared to the internal b–c pair (whose ρ distribution has larger first moment and standard deviation). While similar considerations apply to the octamer (Fig. S13 in SI), the dodecamer shows the presence of both strongly stacked tetramers and dimers (Fig. S14 in SI).

Fig. 5 Distributions (left) and 2D heat maps (right) of ρ and R^π (see Fig. 2 for definition), computed along either the PDI₂@H₂O and PDI₄@H₂O trajectories for (a) the PDI₂@H₂O dimer; (b) the a–b PDI pair (see inset) within the tetramer; (c) the central b–c pair; (d) the c–d pair.

To better characterize the mutual arrangement of PDI units within the column and shed further light on their internal dynamics, we examined the correlation between ρ and the other geometrical supramolecular descriptors defined in the Methods Section. Results are again displayed as 2D heat maps in Fig. 6 and 7, where we selected specific descriptors to highlight the most relevant findings. The complete dataset is reported in the SI, Fig. S15–S31. To corroborate the model of “strongly interacting dimers” within even numbered aggregates, we report in Fig. 6 the data calculated for the external a–b pair monitored in the PDI₂@H₂O, PDI₄@H₂O and PDI₁₂@H₂O systems, while those for the inner b–c pair within either the tetramer or the dodecamer are displayed in Fig. 7. The first two columns of both Figures show the ρ, R^∥ and ρ, R^S distributions, respectively. In general, the distribution obtained for the sliding displacement along the PDI long axis (R^∥) always shows a maximum around 2.2 Å, while the most probable distance between two neighboring units, ρ, lays within 3.6 Å and 3.9 Å, depending on both the dimension of the aggregate and the considered pair. A similar behavior was registered for the R^S, ρ distribution, indicating that the energy paths along either the sliding or the shifting directions are equivalent, and a “circular” motion within the π-plane should be expected. Indeed, the distinctly asymmetric shape observed in all cases further suggests that while the distance ρ between the centers of mass of two neighboring units remains largely unchanged during the dynamics (confirming the stability of the aggregates), greater mobility is allowed for sliding and shifting displacements within planes perpendicular to the column axis. Interestingly, in line with the results displayed in Fig. 5, the conformational dynamics of the solvated dimer in PDI₂@H₂O explores wider regions of both the ρ, R^∥ and ρ, R^S landscapes, resulting in less defined and broadened heat maps displayed in Fig. 6a. Moreover, we note that in the dodecamer, sliding and shifting motions are larger with respect to the tetramer, again pointing to weaker aggregation as the size of the system increases. Turning to the comparison between external and inner PDI pairs, the lower populations displayed for the latter in the (a) and (b) panels of Fig. 7 confirm that, within both the tetramer and the dodecamer, b–c inner pairs exhibit a more pronounced mobility, consistently with the reduced stacking strength reported in Tables S2 and S3.

Fig. 6 2D heat maps of the distributions of the ρ distances with the geometrical descriptors defined in Fig. 2, and computed for the external (a–b) pair along the MD trajectories produced for (a) PDI₂@H₂O, (b) PDI₄@H₂O, (c) PDI₁₂@H₂O. In each panel, columns refer, from left to right, to the sliding (R^∥), shifting (R^S), spinning (α) and tumbling (γ′) motions. For the sake of clarity, a graphical representation of the displacement associated with each geometrical descriptor is displayed in the top panel.

Fig. 7 2D heat maps of the distributions of the ρ distances with the geometrical descriptors defined in Fig. 2, and computed for the inner (b–c) pair along the MD trajectories produced for (a) PDI₈@H₂O, (b) PDI₁₂@H₂O. In each panel, columns refer, from left to right, to the sliding (R^∥), shifting (R^S), spinning (α) and tumbling (γ′) motions. For the sake of clarity, a graphical representation of the displacement associated with each geometrical descriptor is displayed in the top panel.

After assessing the stability of the stacked self-assemblies through the above analyses of positional descriptors, additional insights into the reciprocal orientation of neighboring units within the stacks can be gained by examining the distribution of spinning and tumbling angles, as shown in the two rightmost columns of Fig. 6 and 7. As far as the former (α, see Fig. 2) is concerned, heat maps indicate that the long axes of the considered pair are never found completely aligned (α = 0°). They have maxima at about 30° and 150°, hence pointing toward staggered arrangements.^55,59 In line with analyses above, comparison of the distribution of the spinning angle achieved for either the outer (a–b, c–d) or the inner (b–c) pairs within the tetramer suggests that the system is well described as two interacting dimers. In PDI₈@H₂O (see Fig. S22 in SI), as well as in PDI₁₂@H₂O (Fig. S23), the inner dimer shows a more constrained dynamics, where only one minimum is thoroughly explored. This is likely a consequence of the increased steric hindrance caused by the side chains, which restricts the spinning of the inner PDI units within the large column. Moving to tumbling descriptors, it is not surprising that these movements are rather limited in π-stacked columnar aggregates. We report the ρ vs. γ′ correlation maps in the right panels of Fig. 6 and 7, while all data concerning monomers rotational dynamics within the columnar stack are gathered in Section B.3.2 in the SI. Although, as shown in Fig. 6 and 7, the tumbling motion between two neighboring PDI units is more pronounced in the solvated dimer, both in the case of the outer a–b pairs and of the inner b–c pairs, the non-negligible γ′ population around 90° suggests a residual dynamics around the tumbling (and rolling, see γ in Fig. S24–S26) axes. The limited rotational dynamics upon stacking corresponds to a significant reduction of the configurational space explored by each pair, thereby contributing to the entropic penalty that ultimately leads to the anti-cooperative assembly mechanism. The reduced effectiveness of the self-assembly process with the increasing dimensions of the aggregate is also evident when considering the dephasing angle δ, the supramolecular descriptors defined in Fig. S9. Fig. S30 and S31 report 2D heat maps for the distribution functions of the ρ and δ descriptors along the dynamics of all considered stacked systems. It should be noticed that in a perfect stack of N units: (i) the column height h_π is , where represents the pair stacking distance, as ρ and R_π are identical in a perfect stack, and (ii) δ_k is zero for each unit k, since each lays parallel to the columnar vector . However, our simulations reveal that the δ distributions reported in Fig. S30 and S31 exhibit a significant population extending up to 30°, induced by the residual rolling and tumbling motions of each PDI unit within the stack. Although the aggregates deviate from a perfect column, as evidenced by the angular variations which point to a tilted stack, the distance (ρ) between the first and last unit remains comparable to that observed in the tetramer case. This supports the idea that the PDI column follows a dimer-based dynamic behavior. Consequently, the system exhibits anti-cooperative behavior, as increasing the column size does not enhance stability, which would otherwise be expected in a cooperative assembly.

3.3. Analysis of the first solvation shell

Obtaining atomistic insights into the structure of solvent/aggregates interface is of paramount importance when considering their use in various technological applications as, for instance, in photocatalysis. Monomer functionalization has, indeed, direct impact on the optoelectronic and redox properties but directly influences self-assembling and solvation properties. Thus, characterizing solvent-accessible catalytic sites within the nanostructure provides rational guidelines for optimizing the topology of electron/hole separation and, consequently, boosting photocatalytic activity.¹¹¹ Atomistic details of aggregate/water interactions and the solvation shells surrounding π-stacked columns were obtained by calculating appropriately selected solute–solvent pair correlation functions, g_ab(r), reported in Section C in the SI, and by plotting the three-dimensional spatial density function (SDF) of the water proton, shown in Fig. 8. Since similar features were found for all considered aggregates, we limit the discussion to the tetramer, taken as representative of larger aggregates, while the complete set of data can be found in the SI. Panels (a) and (b) of Fig. 8 show clearly that a well-defined water shell surrounds both charged side chains and the solvated monomer π-core. Conversely, for the tetramer (see panels (c) to (e)), water only externally wraps the PDI π-stacked system, mainly interacting with side chains and, to a minor extent, with external PDI units. To characterize the interface between the solvent and the aggregate in more detail, we computed g_ab(r) between water hydrogen or oxygen atoms (H_w, O_w) and selected positions on the PDIs, namely the geometrical center of central ring (R_k, see Fig. S32 in the SI) of the k-th stacked unit, the Nitrogen atom at the imide (N₁), the charged ammonium Nitrogen (N₂) and the carbonyl Oxygen (O). On the one hand, by looking at Fig. S37, g_R1(4) − H_w points to a non-negligible probability to find the proton within 4 Å from the core, whereas g_R2(3) − H_w is still null in that range, confirming the scarce solvation of inner units. Nonetheless, each PDI monomer in the aggregates still maintains core/water interactions through a relatively strong HB with the carbonyl oxygen O, whose radial distribution function (see Fig. S38 and S39) shows the insurgence of peaks at around 2 Å. This interaction is also clearly visible in panel (c) of Fig. 8.

Fig. 8 Three-dimensional SDF maps of hydrogen and oxygen atoms from water solvent around solute selected moieties: (a) the monomer aromatic core; (b) the monomer's chain; (c) the tetramer aromatic cores. (d) One pendant chain of an external PDI of the tetramer; (e) one chain of an internal PDI of the tetramer as well. Isovalue δ = 12.

Prompted by the outcomes of the previously discussed thermodynamic analysis, which pointed to the key role of PDI–water interactions in the aggregation process, we further investigated the solvation patterns by quantifying the number of HBs formed between the carbonyl oxygen atoms and surrounding water molecules across different aggregate sizes. As shown in Fig. S47, the number of H-bonds per monomer gradually decreases with aggregation size and reaches a plateau around 1.35. This confirms that solvation is present throughout the stack, though slightly reduced in inner units due to steric shielding caused by π–π stacking. We then examined solvation at the polar side chains by analyzing the interactions between the positively charged nitrogens (N₂) and water molecules. As shown in Fig. S48, the average number of N₂–water contacts per monomer remains high across all aggregate sizes, with only a slight decrease compared to the carbonyl oxygens. This indicates that side-chain solvation is largely preserved even in larger stacks, albeit with a modest reduction as the column grows. This picture is consistent with the SDF plots shown in panels (d) and (e) of Fig. 8, which confirm the persistent hydration of the side chains.

Finally, the pair distribution functions calculated for the imide N (N₁) are always found to be less structured, even at short distances, reflecting a reduced solvent exposure due to their proximity to the aromatic core.

To complement this static view, we also examined how water molecules are forced out the inter-layer region during dimer formation. Fig. 9 shows the joint evolution of the interplanar distance, R^π_ab, and the number of interlayer waters. Rather than a simple monotonic and continuous expulsion flow, the system behaves as a bellows: the PDI cores first approach each other to form a quasi-parallel arrangement, then transiently re-separate, and only later collapse into the final stacked configuration. This pathway reflects the energetic cost of displacing interfacial water molecules, which must overcome their favorable hydrogen-bond interactions before being released into bulk. As a result, water expulsion proceeds through intermediate, partially solvated states rather than a single concerted event. Although these observations provide mechanistic insight, we stress that the description remains qualitative, and a statistically converged characterization would require longer trajectories and extended sampling.

Fig. 9 Time evolution of the stacking distance R_p between two PDI cores along a 1 ns MD trajectory, color–coded by the number of confined water molecules N_w. Red dashed boxes mark the time intervals from which representative snapshots (1–4) are extracted and shown as insets, illustrating the progressive expulsion of interlayer waters during π-stack formation. Water oxygens are displayed as red spheres, PDI cores as orange planes, and side chains in stick representation.

4. Conclusions and final remarks

In this study, we investigate the self-assembly dynamics and the resulting supramolecular structure of columnar stacked aggregates of N,N-bis(2-(trimethylammonium)ethylene) perylene 3,4,9,10-tetracarboxylic acid bis-imide (PDI) in water, relying on extended classical molecular dynamics simulations. These simulations are based on an accurate QMD-FF that was previously parametrized by some of us for the PDI monomer,⁸⁸ and further validated for the dimer.⁵⁵ Here, the self-assembly mechanism of trimeric and tetrameric supramolecular columnar aggregates has been investigated using a bottom-up approach. For larger aggregates, octamer and dodecamer, we used a top-down scheme, starting simulations from pre-assembled larger π-stacked columns. The presence of these larger structures was established using enhanced sampling MD techniques. In fact, while a trimer and a tetramer were found to spontaneously self-assemble in around 50 ns, HREMD simulations allowed the formation of larger aggregates, up to a decamer. In line with experimental evidence on self-assembly of PDI bearing bulky side chains,¹⁰⁶ we highlight, through enhanced sampling simulations and free-energy calculations, that PDI self-assembly follows an anti-cooperative, dimer-based mechanism, in which dimers act as the primary building blocks for larger aggregates. In-depth thermodynamic analysis reveals that PDI aggregation is enthalpically disfavored but driven by a favorable entropy gain originating from the release of solvating water molecules. In a model case of a step-wise aggregation process, we shows that both enthalpic penalties and entropic benefits progressively plateau with increasing stack size, leading to a size-independent free energy. This saturation explains the anti-cooperative nature of the process and the limited extension of columnar stacks. Consequently, this suggests that experimental properties observed for these systems could be described by taking into account up to tetramers or hexamers. Although larger aggregates might be formed, these represent a small fraction of the species observed.

Once the self-assembled PDI columns are formed, their stability and dynamics in water was analyzed in terms of the relative motions of their building blocks, defining a complete set of geometrical supramolecular descriptors, accounting for both translational and rotational motions within the π stacks. Regardless of the size of the aggregate, we found that the average distance between PDI falls within the range of 3.7–3.9 Å, in agreement with experimental data for similar π-stacked structures.^104,105 Moreover, analysis of the relative freedom of motion of adjacent dimers within different-sized aggregates further confirms that the dimer, rather than the monomer, serves as the primary “building block”. Notably, we observed that the dynamics of even-numbered aggregates resembles those of strongly interacting dimers.

Atomistic analysis of solvent–aggregate interfaces highlights the pivotal role of solvation in dictating both structural organization and aggregation energetics. While side chains remain persistently hydrated across aggregate sizes, solvation of the π-cores and carbonyl oxygens decreases with stacking, plateauing at ∼1.35 HBs per monomer. This reduction reflects steric shielding within inner units and directly correlates with the enthalpic cost of aggregation. Pair correlation and spatial density analyses confirm that hydration is localized at external units and polar groups, whereas imide nitrogens remain poorly solvated. Moreover, water expulsion during stacking follows a stepwise, bellows-like mechanism, requiring transient rearrangements to overcome the energetic penalty of disrupting hydrogen-bonded interfacial waters. Overall, these findings establish that solvent structuring and water displacement are central to the anti-cooperative nature of PDI self-assembly and provide molecular-level guidelines for tuning aggregate interfaces in photocatalytic applications.

Author contributions

Marta Cantina: Investigation, formal analysis, visualization, writing – original draft. Daniele Padula: Investigation, formal analysis, visualization, writing – review & editing. Alekos Segalina: Investigation, writing – review & editing. Samuele Giannini: Resources, investigation, writing – review & editing. Fabrizio Santoro: Conceptualization, writing – review & editing. Giacomo Prampolini: Conceptualization, software, methodology, supervision, review & editing. Mariachiara Pastore: Project administration, resources, conceptualization, supervision, review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

QMD-FF topologies, equilibrated configurations, and MD specifics are available for hydrated PDI dimers, trimers, tetramers, octamers and dodecamers are available free of charge in GROMACS format at https://doi.org/10.5281/zenodo.15166891.

Supplementary information: Section A outlines the aggregates and geometrical descriptors that characterize their structural features. Section B focuses on self-assembly and supramolecular dynamics, including the binding and interaction energies of different aggregates and their translational, rotational, and columnar dephasing dynamics. Section C presents the solute-solvent radial distribution functions (g_ab(r)) for the PDI aggregates, highlighting insights into the solvation environment and its interaction with the aggregates. Section D report a detailed thermodynamic analysis of the aggregation process. See DOI: https://doi.org/10.1039/d5nr02723c.

Acknowledgements

M. P. and M. C. acknowledge HPC resources from mesocentre EXPLOR of the University of Lorraine (project 2018CPMXX0602). D. P. acknowledges funding from Programma CN00000013 from ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU – PNRR, Missione 4 Componente 2 Investimento 1.4, SPOKE 7 Materials & Molecular Sciences (CUP B93C22000620006), funding (PSR2024 – F-NF) and computational resources kindly provided by Università di Siena (hpc@dbcf). S. G., G. P. and F. S. thank the support of ICSC-Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union-NextGenerationEU-PNRR, Missione 4 Componente 2 Investimento 1.4. All authors gratefully acknowledge Professor Javier Cerezo for the many useful discussions and the stimulating ideas.

References

1. M. Al Kobaisi, S. V. Bhosale, K. Latham, A. M. Raynor and S. V. Bhosale, Chem. Rev., 2016, 116, 11685–11796.
2. J. Wu, M. Peng, M. Mu, J. Li and M. Yin, Supramol. Mater., 2023, 2, 100031.
3. S. Chen, P. Slattum, C. Wang and L. Zang, Chem. Rev., 2015, 115, 11967–11998.
4. K. Ariga, M. Nishikawa, T. M. J. Takeya, L. K. Shrestha and J. P. Hill, Sci. Technol. Adv. Mater., 2019, 20, 51–95.
5. N. Kihal, A. Nazemi and S. Bourgault, Nanomaterials, 2022, 12, 2079–4991.
6. A. Ashcraft, K. Liu, A. Mukhopadhyay, V. Paulino, C. Liu, B. Bernard, D. Husainy, T. Phan and J.-H. Olivier, Angew. Chem., Int. Ed., 2020, 59, 7487–7493.
7. F. Würthner, Chem. Commun., 2004, 1564–1579.
8. B. Zhang, I. Lyskov, L. J. Wilson, R. P. Sabatini, A. Manian, H. Soleimaninejad, J. M. White, T. A. Smith, G. Lakhwani, D. J. Jones, K. P. Ghiggino, S. P. Russo and W. W. H. Wong, J. Mater. Chem. C, 2020, 8, 8953–8961.
9. S. Berardi, V. Cristino, M. Canton, R. Boaretto, R. Argazzi, E. Benazzi, L. Ganzer, R. Borrego Varillas, G. Cerullo, Z. Syrgiannis, F. Rigodanza, M. Prato, C. A. Bignozzi and S. Caramori, J. Phys. Chem. C, 2017, 121, 17737–17745.
10. P. M. Kazmaier and R. Hoffmann, J. Am. Chem. Soc., 1994, 116, 9684–9691.
11. H. Ben, G. Yan, Y. Wang, H. Zeng, Y. Wu, F. Lin, J. Zhao, W. Du, S. Zhang, S. Zhou, J. Pu, M. Ye, H. Ji and L. Lv, Molecules, 2024, 29, 1420–3049.
12. Z. Yang and X. Chen, Acc. Chem. Res., 2019, 52, 1245–1254.
13. H. Wang, K.-F. Xue, Y. Yang, H. Hu, J.-F. Xu and X. Zhang, J. Am. Chem. Soc., 2022, 144, 2360–2367.
14. M. Sun, K. Müllen and M. Yin, Chem. Soc. Rev., 2016, 45, 1513–1528.
15. Q. Jiang, J. Zhang, Z. Mao, Y. Yao, D. Zhao, Y. Jia, D. Hu and Y. Ma, Adv. Mater., 2022, 34, 2108103.
16. Z.-G. Zhang, B. Qi, Z. Jin, D. Chi, Z. Qi, Y. Li and J. Wang, Energy Environ. Sci., 2014, 7, 1966–1973.
17. J. Cao and S. Yang, RSC Adv., 2022, 12, 6966–6973.
18. E. Kozma and M. Catellani, Dyes Pigm., 2013, 98, 160–179.
19. C. Yan, S. Barlow, Z. Wang, H. Yan, A. K.-Y. Jen, S. R. Marder and X. Zhan, Nat. Rev. Mater., 2018, 3, 18003.
20. C. W. Tang and S. A. VanSlyke, Appl. Phys. Lett., 1987, 51, 913–915.
21. S. Chen, Z. Li, Y. Qiao and Y. Song, J. Mater. Chem. C, 2021, 9, 1126–1149.
22. G. Zhang, J. Zhao, P. C. Y. Chow, K. Jiang, J. Zhang, Z. Zhu, J. Zhang, F. Huang and H. Yan, Chem. Rev., 2018, 118, 3447–3507.
23. J. E. Anthony, Chem. Mater., 2011, 23, 583–590.
24. Y. Li, X. Zhang and D. Liu, J. Photochem. Photobiol., C, 2021, 48, 100436.
25. B. Yang, L. Lu, S. Liu, W. Cheng, H. Liu, C. Huang, X. Meng, R. D. Rodriguez and X. Jia, J. Mater. Chem. A, 2024, 12, 3807–3843.
26. X. Lin, Y. Hao, Y. Gong, P. Zhou, D. Ma, Z. Liu, Y. Sun, H. Sun, Y. Chen and S. Jia, et al. , Nat. Commun., 2024, 15, 5047.
27. F. S. Kim, G. Ren and S. A. Jenekhe, Chem. Mater., 2011, 23, 682–732.
28. M. Hecht and F. Wurthner, Acc. Chem. Res., 2020, 54, 642–653.
29. D. Bialas, E. Kirchner, M. I. Röhr and F. Würthner, J. Am. Chem. Soc., 2021, 143, 4500–4518.
30. F. Marin, A. Zappi, D. Melucci and L. Maini, Mol. Syst. Des. Eng., 2023, 8, 500–515.
31. J. Zhou, L. Xue, Y. Shi, X. Li, Q. Xue and S. Wang, Langmuir, 2012, 28, 14386–14394.
32. Z. Chen, V. Stepanenko, V. Dehm and P. Prins, Chem. – Eur. J., 2007, 13, 436–449.
33. F. Würthner, C. R. Saha-Möller, B. Fimmel, S. Ogi, P. Leowanawat and D. Schmidt, Chem. Rev., 2016, 116, 962–1052.
34. Y. Hou, Z. Zhang, S. Lu, J. Yuan, Q. Zhu, W.-P. Chen, S. Ling, X. Li, Y.-Z. Zheng, K. Zhu and M. Zhang, J. Am. Chem. Soc., 2020, 142, 18763–18768.
35. S. Kumagai, H. Ishii, G. Watanabe, C. P. Yu, S. Watanabe, J. Takeya and T. Okamoto, Acc. Chem. Res., 2022, 55, 660–672.
36. I. Biran, S. Rosenne, H. Weissman, Y. Tsarfati, L. Houben and B. Rybtchinski, Cryst. Growth Des., 2022, 22, 6647–6655.
37. M. R. Hansen, R. Graf, S. Sekharan and D. Sebastiani, J. Am. Chem. Soc., 2009, 131, 5251–5256.
38. N. J. Hestand and F. C. Spano, J. Chem. Phys., 2015, 143, 244707.
39. T. Seki, X. Lin and S. Yagai, Asian J. Org. Chem., 2013, 2, 708–724.
40. S. Parida, S. K. Patra and S. Mishra, ChemPhysChem, 2022, 23, e202200361.
41. F. Würthner, Z. Chen, V. Dehm and V. Stepanenko, Chem. Commun., 2006, 11, 1188–1190.
42. J. Wang, A. Kulago, W. R. Browne and B. L. Feringa, J. Am. Chem. Soc., 2010, 132, 4191–4196.
43. P. Singh Bisht, R. Garg, N. Nakka and A. K. Mondal, J. Phys. Chem. Lett., 2024, 15, 6605–6610.
44. Y. Liu, Z. Li, M.-W. Wang, J. Chan, G. Liu, Z. Wang and W. Jiang, J. Am. Chem. Soc., 2024, 146, 5295–5304.
45. C. Otsuka, S. Takahashi, A. Isobe, T. Saito, T. Aizawa, R. Tsuchida, S. Yamashita, K. Harano, H. Hanayama, N. Shimizu, H. Takagi, R. Haruki, L. Liu, M. J. Hollamby, T. Ohkubo and S. Yagai, J. Am. Chem. Soc., 2023, 145, 22563–22576.
46. R. S. Wilson-Kovacs, X. Fang, M. J. Hagemann, H. E. Symons and C. F. Faul, Chem. – Eur. J., 2022, 28, e202103443.
47. K. Bag, R. Halder, B. Jana and S. Malik, J. Phys. Chem. C, 2019, 123, 6241–6249.
48. X. Yang, X. Xu and H.-F. Ji, J. Phys. Chem. B, 2008, 112, 7196–7202.
49. E. R. Draper, B. J. Greeves, M. Barrow, R. Schweins, M. A. Zwijnenburg and D. J. Adams, Chem, 2017, 2, 716–731.
50. R. E. Ginesi, N. R. Murray, R. M. Dalgliesh, J. Doutch and E. R. Draper, Chem. – Eur. J., 2023, 29, e202301042.
51. K. Gayen, S. Paul, S. Hazra and A. Banerjee, Langmuir, 2021, 37, 9577–9587.
52. P. W. J. M. Frederix, I. Patmanidis and S. J. Marrink, Chem. Soc. Rev., 2018, 47, 3470–3489.
53. G. Prampolini, L. Greff da Silveira, J. G. Vilhena and P. R. Livotto, J. Phys. Chem. Lett., 2022, 13, 243–250.
54. A. S. Bondarenko, I. Patmanidis, R. Alessandri, P. C. T. Souza, T. L. C. Jansen, A. H. de Vries, S. J. Marrink and J. Knoester, Chem. Sci., 2020, 11, 11514–11524.
55. A. Segalina, D. Aranda, J. A. Green, V. Cristino, S. Caramori, G. Prampolini, M. Pastore and F. Santoro, J. Chem. Theory Comput., 2022, 18, 3718–3736.
56. V. Bhat, C. P. Callaway and C. Risko, Chem. Rev., 2023, 123, 7498–7547.
57. S. Pal and B. C. Ray, Molecular dynamics simulation of nanostructured materials: An understanding of mechanical behavior, CRC Press, 2020.
58. P. Clancy, Chem. Mater., 2011, 23, 522–543.
59. A. Segalina, X. Assfeld, A. Monari and M. Pastore, J. Phys. Chem. C, 2019, 123, 6427–6437.
60. C. Naranjo, A. Doncel-Giménez, R. Gómez, J. Aragó, E. Ortí and L. Sánchez, Chem. Sci., 2023, 14, 9900–9909.
61. G. Yu and M. R. Wilson, J. Mol. Liq., 2022, 345, 118210.
62. J. Baz and N. Hansen, J. Phys. Chem. C, 2019, 123, 8027–8036.
63. C. B. Markegard, A. Mazaheripour, J.-M. Jocson, A. M. Burke, M. N. Dickson, A. A. Gorodetsky and H. D. Nguyen, J. Phys. Chem. B, 2015, 119, 11459–11465.
64. A. Bartlett, C. B. Markegard, D. J. Dibble, A. A. Gorodetsky, S. Sharifzadeh and H. D. Nguyen, Synth. Met., 2019, 253, 146–152.
65. M. A. Mattson, T. D. Green, P. T. Lake, M. McCullagh and A. T. Krummel, J. Phys. Chem. B, 2018, 122, 4891–4900.
66. E. R. Draper, L. Wilbraham, D. J. Adams, M. Wallace, R. Schweins and M. A. Zwijnenburg, Nanoscale, 2019, 11, 15917–15928.
67. N. Oldani, V. M. Freixas, D. Ondarse-Alvarez, S. Sharifzadeh, T. Gibson, S. Tretiak and S. Fernandez-Alberti, J. Chem. Theory Comput., 2024, 20, 5820–5828.
68. A. F. de Moura and M. Trsic, J. Phys. Chem. B, 2005, 109, 4032–4041.
69. J. Idé, R. Méreau, L. Ducasse, F. Castet, Y. Olivier, N. Martinelli, J. Cornil and D. Beljonne, J. Phys. Chem. B, 2011, 115, 5593–5603.
70. V. Marcon, D. W. Breiby, W. Pisula, J. Dahl, J. Kirkpatrick, S. Patwardhan, F. Grozema and D. Andrienko, J. Am. Chem. Soc., 2009, 131, 11426–11432.
71. N. Meftahi, A. Manian, A. J. Christofferson, I. Lyskov and S. P. Russo, J. Chem. Phys., 2020, 153, 064108.
72. S. S. Panda, K. Shmilovich, N. S. Herringer, N. Marin, A. L. Ferguson and J. D. Tovar, Langmuir, 2021, 37, 8594–8606.
73. Y. Geng, H.-B. Li, S.-X. Wu and Z.-M. Su, J. Mater. Chem., 2012, 22, 20840–20851.
74. M. Hollfelder and S. Gekle, J. Phys. Chem. B, 2015, 119, 10216–10223.
75. M. A. Martínez, D. Aranda, E. Ortí, J. Aragó and L. Sánchez, Org. Chem. Front., 2023, 10, 1959–1967.
76. M. Bonomi and C. Camilloni, Biomolecular Simulations, Springer, 2019.
77. J. Wildman, P. Repiščák, M. J. Paterson and I. Galbraith, J. Chem. Theor. Comput., 2016, 12, 3813–3824.
78. M. Casalegno, A. Famulari and S. V. Meille, Macromolecules, 2022, 55, 2398–2412.
79. L. Jiang, D. M. Rogers, J. D. Hirst and H. Do, J. Chem. Theory Comput., 2020, 16, 5150–5162.
80. C. M. Wolf, K. H. Kanekal, Y. Y. Yimer, M. Tyagi, S. Omar-Diallo, V. Pakhnyuk, C. K. Luscombe, J. Pfaendtner and L. D. Pozzo, Soft Matter, 2019, 15, 5067–5083.
81. J. G. Vilhena, L. Greff da Silveira, P. R. Livotto, I. Cacelli and G. Prampolini, J. Chem. Theory Comput., 2021, 17, 4449–4464.
82. D. Padula, A. Landi and G. Prampolini, Energy Adv., 2023, 2, 1215–1224.
83. S. Giannini, J. Cerda, G. Prampolini, F. Santoro and D. Beljonne, J. Mater. Chem. C, 2024, 12, 10009–10028.
84. I. Cacelli and G. Prampolini, J. Chem. Theory Comput., 2007, 3, 1803–1817.
85. G. König and S. Riniker, Interface Focus, 2020, 10, 20190121.
86. S. Giannini, P. M. Martinez, A. Semmeq, J. P. Galvez, A. Piras, A. Landi, D. Padula, J. G. Vilhena, J. Cerezo and G. Prampolini, J. Chem. Theory Comput., 2025, 21, 3156–3175.
87. B. J. Williams-Noonan, A. Kamboukos, N. Todorova and I. Yarovsky, Chem. Phys. Rev., 2023, 4, 021304.
88. A. Segalina, J. Cerezo, G. Prampolini, F. Santoro and M. Pastore, J. Chem. Theory Comput., 2020, 16, 7061–7077.
89. J. Wang, W. Wang, P. A. Kollman and D. A. Case, J. Mol. Graphics Modell., 2006, 25, 247–260.
90. W. L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc., 1988, 110, 1657–1666.
91. 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.
92. B. Hess, B. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
93. L. Greff da Silveira, P. R. Livotto, D. Padula, J. G. Vilhena and G. Prampolini, J. Chem. Theory Comput., 2022, 18, 6905–6919.
94. M. J. Ahrens, L. E. Sinks, B. Rybtchinski, W. Liu, B. A. Jones, J. M. Giaimo, A. V. Gusev, A. J. Goshe, D. M. Tiede and M. R. Wasielewski, J. Am. Chem. Soc., 2004, 126, 8284–8294.
95. J. S. Hub, B. L. de Groot and D. van der Spoel, J. Chem. Theory Comput., 2010, 6, 3713–3720.
96. D. Giavazzi, M. F. Schumacher, L. Grisanti, M. Anzola, F. Di Maiolo, J. Zablocki, A. Lützen, M. Schiek and A. Painelli, J. Mater. Chem. C, 2023, 11, 8307–8321.
97. N. Romagnoli and D. Padula, J. Phys. Chem. C, 2024, 48, 19901–19911.
98. G. Bussi, Mol. Phys., 2013, 112, 379–384.
99. The PLUMED Consortium, Nat. Methods, 2019, 16, 670–673.
100. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.
101. M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, F. Pietrucci, R. A. Broglia and M. Parrinello, Comput. Phys. Commun., 2009, 180, 1961–1972.
102. R. Qi, G. Wei, B. Ma and R. Nussinov, in Replica Exchange Molecular Dynamics: A Practical Application Protocol with Solutions to Common Problems and a Peptide Aggregation and Self-Assembly Example, Springer, New York, 2018, pp. 101–119.
103. D. Sindhikara, Y. Meng and A. E. Roitberg, J. Chem. Phys., 2008, 128, 024103.
104. A. E. Clark, C. Qin and A. D. Q. Li, J. Am. Chem. Soc., 2007, 129, 7586–7595.
105. S. Ghosh, X.-Q. Li, V. Stepanenko and F. Würthner, Chem. – Eur. J., 2008, 14, 11343–11357.
106. J. Gershberg, F. Fennel, T. H. Rehm, S. Lochbrunner and F. Würthner, Chem. Sci., 2016, 7, 1729–1737.
107. W. E. Ford, J. Photochem., 1986, 34, 43–54.
108. C. Shao, M. Grüne, M. Stolte and F. Würthner, Chem. – Eur. J., 2012, 18, 13665–13677.
109. Z. Chen, U. Baumeister, C. Tschierske and F. Würthner, Chem. – Eur. J., 2007, 13, 450–465.
110. P. P. Syamala and F. Würthner, Chem. – Eur. J., 2020, 26, 8426–8434.
111. M. Sachs, R. S. Sprick, D. Pearce, S. A. Hillman, A. Monti, A. A. Guilbert, N. J. Brownbill, S. Dimitrov, X. Shi and F. Blanc, et al. , Nat. Commun., 2018, 9, 4968.

---