Analysing π–π-stacking interactions in lignin nanoparticles from molecular simulations – insights and lessons learned

Abstract

A special molecular association is π–π-stacking, driven by weak interactions within aromatic compounds. The π–π-stacking interactions can occur in either a sandwich-like or T-shaped manner. In this study, a method to recognise π–π-stacking from classical molecular dynamics trajectories is developed. By applying three criteria, the method is tested for simple lignin dimer, tetramer and octamer systems, with all G units and β-O4′ linkages. The criteria are geometric and based on distance between ring centroids, the angle between the planes of the two rings and the lateral displacement of the rings. In addition, a wide-angle X-ray scattering (WAXS) profile was calculated from a tetramer system, in agreement with previous experimental results. However, when the WAXS peak assigned to sandwich-shaped stacking was analysed in-depth, it was found to mainly be caused by other intramolecular structural motifs involving e.g. the α-carbon and ring carbons, rather than π–π-stacking. This finding is important for future analyses of WAXS profiles originating from lignin-based materials and shows the strength of combining X-ray scattering methods with molecular modelling.

---

Lignin is a heterogeneous, aromatic biopolymer that is found in plant cell walls where it provides structural integrity and resistance to degradation. Its chemical structure consists of interconnected aromatic rings that can interact by so-called π–π interactions. These non-covalent interactions are weak, but assumed to play a crucial role in lignin self-assembly, solubility, and reactivity. They are hypothesised to influence both biological processes, such as the growth of lignin macromolecules in the cell wall in softwood,¹ as well as industrial processes.

The arrangement of aromatic rings allowing for π–π-interaction can be face-to-face (full or partial, also called sandwich or parallel-displaced) or face-to-side (T-shaped). The face-to-face motifs are sometimes distinguished as J-type (head-to-tail fashion, leading to a redshift in UV-vis absorption spectra) and H-type (full stacking, characterized by a blueshift in UV-vis absorption spectra). By the use of UV-vis absorption spectroscopy, a predominance of J-type π–π-stacking in lignosulfonates as well as alkali lignin has been confirmed.^2,3

Based on what was seen in X-ray diffraction patterns of di- and trilignols,^4,5 X-ray scattering techniques such as wide angle X-ray scattering (WAXS) or X-ray diffraction (XRD) on lignin-based materials show features at distances that correspond well to typical π–π-stacking in what is being referred to as sandwich and T-shaped formations for benzene dimers.^6–9 In addition, WAXS studies have been performed on polymer blends including kraft lignin derivatives¹⁰ and lignosulfonates.¹¹ Besides experimental studies, π–π-stacking was scarcely analysed in computational studies of lignin using Density Functional Theory (DFT) and Molecular Dynamics (MD) simulations.^12,13 It is tempting to consider π–π-stacking a binary property, either existing or not, allowing for quantitative comparisons between molecular models and e.g. WAXS data.

Despite being a (somewhat controversial, overused and maybe misinterpreted^14,15) phenomenon arising from quantum mechanical interactions, it is useful to recognise plausible π–π-stacking from classical MD by applying three simultaneous geometric criteria: the distance between ring center of mass (COM), the angle between the planes of the two rings and the lateral displacement between the rings. However, to recognise for which range of distances, angles and displacements an arrangement is considered π–π-stacking in lignin using atomistic simulations, is a non-trivial task and in contrast to similar analysis of hydrogen bonds using MD simulations there are no standard criteria available. Moreover, although it has been established as a structural feature that can be observed by X-ray scattering techniques, a direct comparison between lignin scattering patterns obtained from molecular models and experimental lignin samples has so far not been made.

The purpose of the present study is twofold. On one hand, a methodology and analysis tool is developed to systematically investigate the presence of π–π-stacking from MD trajectories of lignin. The analysis is performed on model systems of solvated lignin nanoparticles made up of dimers, tetramers and octamers. Geometric criteria are explored and applied to detect and distinguish inter- and intramolecular π–π-stacking from MD simulations. On the other hand, by computing X-ray scattering patterns from atomistic simulations, the contribution of π–π-arrangements to the scattering peaks will be analysed.

1 Methods

1.1 Modeled lignin systems

Analysis is based on trajectories from MD simulations of clusters of lignin oligomers (dimers, tetramers and octamers) in water. A summary of simulated systems is presented in Table 1, divided into the set of simulations used for studying geometric criteria and simulated WAXS spectra respectively. In the criteria simulations, the number of lignin molecules varied between 50 and 200, which means that each simulation box contained 400 aromatic rings, independent of the size of the molecules. The WAXS simulations will be described separately in Section 1.3.1.

Table 1 Simulated systems

	Type	Number of lignins	Number of waters	Box size [nm^3]
Criteria	Dimer	200	29 160	1000
	Tetramer	100	29 085	1000
	Octamer	50	29 022	1000
WAXS	Tetramer	1000	10 670	1325

---

All lignin models used in the current work consist of G-monomers connected with β-O4′ linkages. The chemical structures of a dimer, tetramer and octamer are shown in Fig. S1.† All simulations were performed with GROMACS 2021–2023 (ref. 16) using the CHARMM force fields for lignin¹⁷ and TIP3P water model.¹⁸ The topologies were converted to GROMACS format using the TopoTools plugin¹⁹ for VMD.²⁰

1.2 Simulation protocol and analysis

1.2.1 Geometric analysis. Lignin molecules were randomly placed in a simulation box (10 × 10 × 10 nm³) with periodic boundary conditions. Thereafter, the box was filled with water, aiming at a density of 1000 kg m⁻³.

Energy minimization was performed with steepest descent algorithm. Subsequently, the NPT production run had a total length of 200 ns of which equilibration took place during the first 100 ns and the last 100 ns trajectory was used for analysis. During the equilibration phase, lignin clusters were formed. The time step of the production run was 2 fs, progressing with the leap-frog algorithm. Using the NPT ensemble, the Bussi–Donadio–Parrinello thermostat²¹ and the c-rescale barostat²² maintained 300 K (τ_T = 1 ps) and 1 bar (τ_P = 2 ps), respectively. Covalent hydrogen bonds were constrained by P-LINCS.²³

The Lennard-Jones and Coulomb interactions were truncated at 1.2 nm. The Particle-Mesh Ewald method was used to calculate the long-range electrostatic interactions.^24,25 Dispersion corrections for energy and pressure were applied.

1.3 Analysis of geometric criteria of π–π-stacking

Within a MD simulation, the coordinates of every atom in the system is readily available as a function of simulation time. The presence of π–π-stacking (both sandwich and T-shaped) was determined based on three geometric parameters extracted from these trajectories. The parameters are defined as follows (see Fig. 1):

(1) Aromatic rings COM distance (d)

(2) Angle between aromatic rings plane normals (α)

(3) Aromatic rings COM lateral displacement (l)

Fig. 1 Graphical representation of two aromatic rings and the geometric criteria of their relative position and orientation used to define π–π-stacking.

Radial distribution functions (RDFs) were calculated with the tool according to eqn (1), where 〈ρ_B(r)〉 is the particle density of type B at a distance r around particles A, and 〈ρ_B〉_local is the particle density of type B averaged over all spheres around particles A with radius r_max. In this case, particles A and B are the COM of the aromatic rings. An RDF was calculated for each possible pairwise combination and grouped into intra- (neighbouring and non-neighbouring) and intermolecular contributions.

1.3.1 Computed X-ray scattering. X-ray scattering data was computed from MD trajectories using the Debyer package.²⁶ A system of 1000 lignin tetramers and water to moisture content of 25.01% per weight was used to calculate the WAXS spectra and all associated analysis (simulation 20a in Hartwig-Nair et al.¹³), see WAXS in Table 1. This approach was chosen because systems with higher water content formed lignin clusters that induced an oscillating pattern in the WAXS signal, attributable to the shape of the clusters. The intensity was calculated every 5 ns from the last 50 ns of the simulation and then averaged. A cutoff distance of half the simulation box was used, covering a similar range of the scattering vector q as in experiments. Note that no hydrogens were included in the analysis, as their low electron density makes their contribution to the WAXS signal negligible.

In order to identify the contribution from π–π-stacking to the calculated WAXS intensities, one of every aromatic ring that was found to be part of a sandwich-shaped π–π-stacking pair was removed and the WAXS pattern was recalculated. To support the data, RDFs according to eqn (1) of all and subsets of non-hydrogen atoms were computed.

2 Results

In the following section, results related to the process of selecting geometric criteria is presented first, thereafter computed WAXS spectra.

2.1 Defining geometric criteria for π–π-stacking

Given the three geometric parameters, d, α and l, there were two questions to be investigated.

• What does the distribution of these three parameters look like for the tetramer system?

• How do the distributions compare with theoretical general criteria for π–π-stacking and lignin π–π-stacking observed in DFT simulations?

2.1.1 Distance d. The first parameter of interest was the characteristic distance between aromatic rings. The COM distance between pairs of aromatic rings in the entire tetramer system during a simulation are presented in Fig. 2 as the RDF between the center of mass of aromatic rings, according to eqn (1). The RDF is separated into intramolecular and intermolecular contributions.

Fig. 2 Radial distribution function g(r) of the COM of the aromatic rings in the tetramer simulation described in Table 1 with separated contributions of intra- and inter-molecular distances.

Intramolecularly, there are distinct peaks at approximately 3.9 and 6.8 Å, which can be attributed to neighbouring aromatic rings. In a previous study of lignin dimers, these distances were found to correspond to a folded (4.0 Å) and unfolded or stretched conformation (6.8 Å) of a dimer with β-O4′ linkage in water.²⁷ The intramolecular contributions are further separated into contributions from neighbouring and non-neighbouring rings in Fig. 3. This separation confirms the attribution of the distinct peaks at 3.9 and 6.8 Å to neighbouring rings in folded and stretched conformation, while non-neighbouring rings show a broad peak around 6 Å.

Fig. 3 Radial distribution function g(r) of the COM of the aromatic rings in the tetramer simulation described in Table 1 with separated contributions from intramolecular neighbouring and non-neighbouring rings.

The characteristic distances are not enough to directly indicate the presence of sandwich or T-shaped π–π-stacked single lignin molecules or stacking within lignin aggregates. Therefore, it is of interest to compare RDF peaks with COM distances determined from DFT calculations of lignin oligomers.¹² According to DFT calculations, the RDF peak around 4.0 Å could also mean that sandwich π–π-stacking is taking place. From DFT, T-shaped stacking is expected to appear around 5.0 Å but at this distance no clearly visible peak is seen in the RDF. To further explore the conformational preferences at distances expected for sandwich-shaped and T-shaped stacking, distance ranges were now chosen, based on DFT values in combination with the apparent peaks in the RDF. The ranges are presented in Table 2.

Table 2 Selected cutoff criteria for d

	d (Å)
Sandwich	3.2–4.6
T-shaped	4.0–6.0

---

2.1.2 Angle α and lateral displacement l. After setting the distance ranges, it was investigated whether there are conformational states with angles α and lateral displacements l that correspond to values expected for sandwich or T-shaped π–π-stacking. To this end, the angles and lateral displacements of all conformations that meet the distance criteria in the tetramer simulations are shown in Fig. 4 and 5.

Fig. 4 (Left) 2D plot of all angles α and lateral displacements l recorded for rings with the distance d defined as for sandwich-type π–π-stacking. The red square indicates cutoff-values suggested for α and l. (Right) Snapshot showing an exemplary conformation of a tetramer that meets the criteria of a sandwich-shaped π–π-stacking arrangement.

Fig. 5 (Left) 2D plot of all angles α and lateral displacements l recorded for rings with the distance d defined as for T-shaped π–π-stacking. The red square indicates cutoff-values suggested for α and l. (Right) Snapshot showing an exemplary conformation of a tetramer that meets the criteria of a T-shaped π–π-stacking arrangement.

In Fig. 4, a peak is seen around α = 30° and l = 1.2 Å, meaning that this is a favourable configuration. In a similar way, Fig. 5 shows the angle α versus l of aromatic rings at distances 4.0 to 6.0 Å, as expected for T-shaped arrangements with a maximum around α = 40° and l = 1.2 Å. This is likely, at least partially, corresponding to the peak at similar angle and displacement in Fig. 4, due to the overlap in the distance ranges. A second less clear peak is seen at a higher angle and displacement, approximately at α = 90° and l = 2.5 Å. While the angle of this peak falls within the criteria of T-shaped π–π-stacking, the lateral displacement is almost corresponding to the width of an entire ring and therefore not likely reflecting this type of interaction. The peak itself was found to primarily be due to two neighbouring rings of the tetramers interacting intramolecularly, or intermolecular arrangement within the lignin cluster.

2.1.3 Setting cutoff criteria of α and l. Now, given the populations and conformational preferences of rings within two certain distance ranges from each other, as shown in Fig. 4 and 5, the second question was to investigate how these distributions compare with geometries expected for lignin π–π-stacking.

For the case of sandwich-stacked lignin, the peak seen in Fig. 4 is located within ranges of α and l that would be accepted as sandwich-lateral displaced stacking from theory.⁶ Initial cutoff criteria for α and l were therefore chosen to cover the peak, and the ranges are indicated as a red box in Fig. 4 and the values are presented in Table 3.

Table 3 Selected cutoff criteria for α and l

	α (°)	l (Å)
Sandwich	10.0–50.0	0.0–2.0
T-shaped	70.0–90.0	0.0–2.0

---

For T-shaped stacked lignin, the lack of a clear population around the expected value of α (ideally 90°) and lateral displacements that are less than the size of a ring (approximately 2.8 Å) makes it necessary to base the choice of initial cutoff criteria on theoretical values. Any stacked conformation should appear as an occurrence at angles close to 90° and at a small lateral displacement. Cut-off values suggested in this work are again indicated with a red box in Fig. 5 and presented in Table 3. The lack of a peak indicates that, while there might be T-shaped stackings as the example in the snapshot in Fig. 5, there is no clear preference for such motifs.

2.1.4 Sensitivity to criterion selection. The sensitivity of recorded number of π–π-stackings with respect to the parameters d, α and l was tested. The analysis is reported in detail in the ESI.† By systematic variation of the three criteria for the tetramer system, it was found that the calculated number of sandwich and T-shaped π–π-stacking exhibits in general low sensitivity to the initial choice of criteria (Tables 2, 3 and Fig. S2†). Therefore, these criteria were employed for an extended investigation within the systems of dimers and octamers. It is noted that if lateral displacements between 2.0–3.0 Å were allowed, the weak second population peak in Fig. 5 would result in a significant increase in number of accepted T-shaped stackings.

2.1.5 Application for dimers, tetramers and octamers. When the criteria were defined, the number of stackings was determined for the first set of simulations of dimers, tetramers and octamers (Table 1). The absolute average number of motifs that meet the criteria in each simulation are presented in Fig. 6. The total number of aromatic rings (n) and possible pairs (n_pairs) is the same in all systems , however the fraction of intramolecular vs. intermolecular ring pairs varies due to the different molecular sizes, making it less likely to form an intermolecular pair within systems of larger molecules. The numbers of possible intramolecular pairs are 200 for dimers, 600 for tetramers and 1400 for octamers. The results show that intramolecularly, there are more sandwich π–π-stackings compared to T-shaped π–π-stackings independent of the molecular weight of the lignin molecules. However, the difference increases drastically when comparing dimers to the larger lignin molecules. T-shaped stacking is more common inter- than intramolecularly, and the absolute values are not varying much with increasing molecular size.

Fig. 6 Absolute number of recorded intra- and intermolecular sandwich and T-shaped stacking occurrences.

2.1.6 Comparison to density functional theory. Based on structure files from a previously published DFT study of lignin hexamers consisting of all G units connected with β-O 4′ linkages,¹² values for d, α and l found in these molecules are listed in Table 4. It is noted that the ranges of α and l for sandwich-shaped stacking from DFT are 16 ≤ α ≤ 31° and 0.7 ≤ l ≤ 1.4 Å, which all fall within the criteria shown in Table 3.

Table 4 Geometric parameters from density functional theory

	d (Å)	α (°)	l (Å)
Sandwich	4.0	30	1.4
	3.7	16	0.7
	4.0	31	1.4
T-shaped	4.8	85	1.3

---

2.2 Comparison to wide angle X-ray scattering

The use of X-ray scattering techniques to quantify structural features, specifically π–π-stacking interactions in lignin-based thermosets,^7–9 highlights the importance of comparing experimental X-ray scattering patterns with those calculated from MD simulations. To this end, a significantly more lignin-rich system was analysed, described in Table 1 under WAXS. The calculated WAXS pattern of the entire system (excluding hydrogen atoms) is presented in Fig. 7, together with a 1D azimuthally integrated intensity profile. A broad asymmetric peak can be seen, which decays slower toward lower q values (to the left) from the peak position compared to higher. This is similar to that seen in experiments, with the maximum peak around 1.6 Å⁻¹. The distribution indicates an ordering of distances d = 2π/q around these q-values, which has previously been assigned to distances of π–π-stacked aromatic rings, as lignin is otherwise disordered. In experimental studies,^7–9 two Gaussian functions are commonly fitted to the broad lignin specific peak, assigned to sandwich and T-shaped π–π-stacking respectively. Fig. 7 therefore also includes fitted Gaussian functions.

Fig. 7 Simulated WAXS: 2D scattering (left) and 1D integrated peak at relevant q-values including two fitted Gaussian peaks in the same way as previous experiments (right).

Recalculated WAXS spectra with removed sandwich-shaped π–π-stacking pairs as explained in Section 1.2 effectively removes any sandwich stacking contribution to the WAXS signal. It is seen in Fig. 8 that when doing this, the WAXS peak is still present with the same overall shape, meaning that the sandwich-shaped π–π-stacked rings in the studied system are not the main contributor to the peak. Therefore, the question arose: what are the structural features that give rise to the peak?

Fig. 8 Simulated WAXS based on all non-hydrogen atoms in the system (blue line) and all non-hydrogen atoms except carbons that are part of sandwich-shaped π–π-stacking rings (orange line).

Since WAXS intensities are directly related to the spatial coordinates of all atoms, RDFs can help identify structural features that contribute to scattering. The RDF of all non-hydrogen atoms (Fig. 9a) reveals distinct pair correlation within the range of the broad WAXS peak. Notably, the first RDF peak aligns closely with the WAXS intensity maximum. The chemical structure of the simulated lignin is relatively simple and regular. Besides the aromatic rings, each unit has an α-carbon as part of the linkage (see the tetramer in Fig. S1†). As a subsequent step, the corresponding RDFs for ring carbons and ring and α-carbons were investigated separately (Fig. 9b and c). The ring and α-carbon RDF shows a prominent peak at ∼3.8 Å and a secondary peak at ∼4.3 Å, not observed for the ring carbons alone.

Fig. 9 RDF and graphics indicating included atoms: (a) all non-hydrogen atoms, (b) ring carbons and (c) ring and α-carbons.

Closer inspection of intramolecular distances (Fig. 9c) show that the smaller RDF peak at ∼4.3 Å corresponds to the distance between the α-carbon and the furthest ring carbon (green line). The more intense peak at ∼3.8 Å arises from two equivalent distances between the α-carbon and the two carbons adjacent to the aforementioned ring carbon, leading to approximately double the frequency in the RDF. These longer distances are absent in the RDF for ring carbons alone (Fig. 9b) where the maximum distance within the ring is ∼2.8 Å.

Similarly, the two oxygen atoms covalently bonded to the aromatic ring also contribute to peaks at comparable distances in the RDF. These distances are slightly shorter than the carbon equivalents, due to the smaller atomic radius of oxygen, but the higher electron density of oxygen increases their impact on the WAXS signal. These structural motifs are consistent across G-type monomers, owing to the conserved bonding pattern. Additional distribution in the 4–5 Å range arises from atoms located at a greater distance from the ring, as seen in the RDF of all non-hydrogen atoms.

Since sandwich-stacked rings were found not to be the primary contributors for the WAXS peak in the simulations, an alternative explanation lies in the intramolecular atomic distances, as revealed by the RDF analysis. This is comparable to how X-ray scattering data taken on liquid water – which despite being amorphous – exhibits a distinct scattering peak at q-values corresponding to oxygen–oxygen distances between neighbouring water molecules. These distances arise from the short-range order imposed by hydrogen bonding. Narten and Levy²⁸ were the first to use X-ray scattering to extract the RDF of liquid water and identified a nearest-neighbour oxygen–oxygen distance around 2.8 Å. Here, RDF peaks from recurrent, covalently defined interatomic distances of lignin were observed in the WAXS peak range, and thus must be considered when interpreting X-ray scattering data.

This observation is in line with the results presented in a publication by Li and Sarkanen,¹⁰ where it was concluded that two components of the amorphous peak from methylated kraft and dioxane lignin may only in part be assigned to arrangements of aromatic rings that are respectively parallel and edge-on in relation to each other. That statement is based on the comparison to other aromatic compounds, and hereby strengthened by atomistic modelling of lignin.

3 Conclusions

A method applying three geometric criteria (distance, angle and lateral displacement) was developed to analyse the presence of π–π-stacking motifs within MD simulations of lignin. The method was tested on systems of dimers, tetramers and octamers in water showing that for all molecular sizes, intramolecular neighbouring rings are the most prone to form π–π-stacking. The method also showed that for molecules larger than dimers, the sandwich-shaped feature dominates. Intermolecular π–π-stacking is less common, and the T-shaped π–π-stacking feature does not seem to be a dominant structural motif in the studied systems.

Moreover, simulated WAXS patterns together with detailed analysis of radial distribution functions, suggest that the WAXS peak around a q-value of 1.5 Å⁻¹ is predominantly due to intramolecular structural features, for instance between the α-carbon and ring carbons, rather than sandwich-shaped π–π-stacking. This finding is important for future analyses of WAXS profiles originating from lignin-based materials and shows the strength of combining advanced scattering techniques with atomistic models.

Data availability

This study was carried out using publicly available scripts and data from https://github.com/KlaraHackenstrass/pi-pi_stacking.

Author contributions

Klara Hackenstrass: methodology, software (molecular dynamics simulations), validation, formal analysis (π–π-stacking), data curation, visualization, writing – review & editing. Nil Tabudlong Jonasson: methodology, software (code development for data analysis), validation, formal analysis (X-ray scattering pattern), data curation, visualization, writing – review & editing. Marie Hartwig-Nair: conceptualization, investigation (analysis with respect to experimental data), resources. Tomas Rosén: supervision regarding analysis of X-ray scattering data, writing – review & editing. Sara Florisson: conceptualisation, project administration, supervision, writing – review & editing, funding acquisition. Malin Wohlert: conceptualisation, project administration, supervision, writing – original draft, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by funding from VINNOVA project number 2021-02086. Computational resources were provided by National Academic Infrastructure for Supercomputing in Sweden (NAISS) under project NAISS 2024/5-82, whose support is gratefully acknowledged. This work was supported by the Swedish Research Council (Vetenskapsrådet) under grant number 2023-04705.

Notes and references

1. N. Terashima, M. Yoshida, J. Hafrén, K. Fukushima and U. Westermark, Holzforschung, 2012, 66, 907–915.
2. Y. Deng, X. Feng, D. Yang, C. Yi and X. Qiu, Bioresources, 2012, 7, 1145–1156.
3. Y. Deng, X. Feng, M. Zhou, Y. Qian, H. Yu and X. Qiu, Biomacromolecules, 2011, 12, 1116–1125.
4. R. Stomberg, K. Lundquist, J. Koziol, F. Müller and M. Sjöström, Acta Chem. Scand., 1987, 41b, 304–309.
5. J.-P. Roblin, H. Duran, E. Duran, L. Gorrichon and B. Donnadieu, Chem.–Eur. J., 2000, 6, 1229–1235.
6. M. O. Sinnokrot, E. F. Valeev and C. D. Sherrill, J. Am. Chem. Soc., 2002, 124, 10887–10893.
7. M. E. Jawerth, C. J. Brett, C. Terrier, P. T. Larsson, M. Lawoko, S. V. Roth, S. Lundmark and M. Johansson, ACS Appl. Polym. Mater., 2020, 2, 668–676.
8. I. Ribca, M. E. Jawerth, C. J. Brett, M. Lawoko, M. Schwartzkopf, A. Chumakov, S. V. Roth and M. Johansson, ACS Sustain. Chem. Eng., 2021, 9, 1692–1702.
9. I. Ribca, B. Sochor, M. Betker, S. V. Roth, M. Lawoko, O. Sevastyanova, M. A. Meier and M. Johansson, Eur. Polym. J., 2023, 194, 112141.
10. Y. Li and S. Sarkanen, Macromolecules, 2005, 38, 2296–2306.
11. Y.-Y. Wang, Y.-r. Chen and S. Sarkanen, Green Chem., 2015, 17, 5069–5078.
12. I. V. Pylypchuk, M. Karlsson, P. A. Lindén, M. E. Lindström, T. Elder, O. Sevastyanova and M. Lawoko, Green Chem., 2023, 25, 4415–4428.
13. M. Hartwig-Nair, A. Nasedkin, K. Hackenstrass, E. De Santis, S. Florisson and M. Wohlert, Wood Sci. Technol., 2025, 59, 21,  DOI:10.1007/s00226-024-01624-4.
14. K. Carter-Fenk and J. M. Herbert, Phys. Chem. Chem. Phys., 2020, 22, 24870–24886.
15. C. R. Martinez and B. L. Iverson, Chem. Sci., 2012, 3, 2191–2201.
16. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1–2, 19–25.
17. J. V. Vermaas, L. Petridis, J. Ralph, M. F. Crowley and G. T. Beckham, Green Chem., 2019, 21, 109–122.
18. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
19. J. V. Vermaas, D. J. Hardy, J. E. Stone, E. Tajkhorshid and A. Kohlmeyer, J. Chem. Inf. Model., 2016, 56, 1112–1116.
20. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graph., 1996, 14, 33–38.
21. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
22. M. Bernetti and G. Bussi, J. Chem. Phys., 2020, 153, 114107.
23. B. Hess, J. Chem. Theory Comput., 2008, 4, 116–122.
24. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
25. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
26. Debyer: X-ray diffraction pattern simulator, 2022, https://github.com/hampusak/debyer, accessed: 2025-04-08.
27. K. Hackenstrass, M. Hasani and M. Wohlert, Holzforschung, 2024, 78, 98–108.
28. A. Narten and H. Levy, J. Chem. Phys., 1971, 55, 2263–2269.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5fd00052a

‡ These authors contributed equally to this work.

---