Nikolaos I.
Sigalas
*ab,
Stefanos D.
Anogiannakis
bc,
Doros N.
Theodorou
bc and
Alexey V.
Lyulin
ad
aSoft Matter and Biological Physics Group, Department of Applied Physics, Technische Universiteit Eindhoven, 5600, MB, Eindhoven, The Netherlands. E-mail: n.sigalas@tue.nl
bDutch Polymer Institute, P.O. Box 902, 5600 AX Eindhoven, The Netherlands
cComputational Materials Science and Engineering Group, School of Chemical Engineering, National Technical University of Athens (NTUA), GR-15780 Athens, Greece
dCenter for Computational Energy Research (CCER), P.O. Box 513, 5600 MB Eindhoven, The Netherlands
First published on 30th March 2022
Understanding the process–property relations of helical polymers using molecular simulations has been an attractive research field over the years. Specifically, isotactic polypropylene still remains a challenge for current computational experimentation, as it exhibits phenomena such as crystallization that emerge on large spatial and temporal scales. Coarse-graining is an efficient technique for approaching such phenomena, although previous coarse-grained models lack in preserving important atomistic and structural details. In this paper we develop a new coarse-grained model, based on the popular MARTINI force field, that is able to reproduce the helical behavior of isotactic polypropylene. To test the model, the predicted statistical and structural properties (characteristic ratio, density, entanglement molecular weight, solubility parameter in the melt) are compared with previous simulation results and available experimental data. For the development of the new coarse-grained force field, a single unperturbed chain Monte Carlo algorithm has been implemented: an efficient algorithm which samples conformations representative of a melt by simulating just a single chain.
Important insights into both structural and dynamical properties of polypropylene can be provided by modern molecular simulations, using molecular-dynamics (MD)11–13 and Monte Carlo (MC) techniques.14,15 Despite considerable progress, phenomena such as the development of polymorphism16 are still not accessible via current simulation computer power. Furthermore, the understanding of complex and slow relaxations in polypropylene melts still remains a challenge.
A promising strategy for exploring in silico the large spatial and temporal scales of polypropylene structural and dynamic behavior is coarse-graining. In a coarse-grained model, a number of heavy atoms is combined into a single coarse-grained interaction site (bead). In several coarse-grained approaches used so far—bottom-up,17–20 top-down21–23 and hybrid24,25—two main characteristics can be distinguished: (a) mapping the atomistic structure onto the coarse-grained level and (b) development of a coarse-grained force field that describes the interactions among beads. Back in 1995, the first attempt at coarse-graining polypropylene was made by Schweizer et al.,26 who tried to implement PRISM theory27 using the rotational isomeric state model (RIS) on different coarse-grained schemes. Then, Mattice and coworkers15,28 developed a coarse-grained iPP model which was placed on a high coordination lattice, with each site corresponding to an iPP monomer. Later, based on the MARTINI force field, a coarse-grained scheme specifically for polypropylene was proposed by Panizon et al.,29 with each bead comprising a methyl group, a methine group and two half-methylene groups (one bead per repeat unit). Recent studies use even coarser schemes, choosing a Kuhn monomer to represent a bead.30,31
The challenge in polypropylene coarse-graining is connected to its structure. Specifically, polypropylene is a branched vinyl polymer which exhibits a tendency to form helices. These helical conformations are drastically affected by the stereo-regularity of the branch. For instance, while isotactic polypropylene blocks can adopt a helical conformation upon cooling down below the melting point, atactic blocks remain amorphous even at low temperatures. Short helical configurations are also present in polypropylene melts. Upon coarse-graining, as the resolution becomes less fine, the methyl branches disappear and their influence on conformation must be accounted for indirectly. Until now there is no off-lattice coarse-grained model that describes polypropylene's helicity.
The purpose of the present study, which focuses on isotactic polypropylene (iPP), is twofold: first, we introduce and check the stability of the iPP helical structure based on existing coarse-grained representations. Secondly, we implement a new coarse-graining method and test the force field produced for iPP. An iPP helix can be represented at the coarse-grained level with distinct handedness and its structure can be preserved. To develop such an iPP coarse-grained model that preserves helicity, a hybrid approach is followed, in the sense that united-atom simulations are carried out for the parametrization of the bonded interactions, while the non-bonded interaction potential is tailored to reproduce experimental results. The novelty of the method lies in that a single unperturbed chain Monte Carlo algorithm is employed at the united-atom level as a reference in order to extract the bonded coarse-grained interactions. The new coarse-grained model has been developed having crystallization in mind. Simulations of crystallization will require a separate study and are not within the scope of this paper. Only some preliminary results are provided in the Discussion section to show that the model has much potential for providing insight into the crystallization of polypropylene.
The paper is organized as follows. In Section 2 we present the models and the simulation techniques that are used in the present study. In Section 3 we develop the new HELical Isotactic Polypropylene Potential (HELIPP) at the coarse-grained level, and validate it by comparing against other existing experimental and simulation data. Finally, in Section 4 we discuss the limitations of the new HELIPP and the possibility to tackle the remaining iPP modelling challenges.
Monte Carlo and molecular-dynamics simulations were performed, as explained later, for a single iPP chain (both in a perfect helical structure and coiled), as well as for iPP melts. All the systems, except for the perfect iPP helical structures, were generated using the Amorphous Builder plugin of the Materials and Process Simulation (MAPS) platform developed by Scienomics.34 The perfect iPP helical structures at the UA level were created based on diffraction data from Immirzi and Iannelli,35 in the same way as it has been recently described by Romanos and Theodorou.36
A perfect iPP helical structure, with three iPP monomers per helical turn, can be either left-handed (LH) or right-handed (RH). At the UA level a RH helix can be spotted as sequences of trans (180°) (CH2–CH–CH2–CH) and gauche− (−60°) (CH–CH2–CH–CH2) dihedral angles along a chain (Fig. 1c) resulting in an up arrangement,37 while a LH helix is spotted as a sequence of gauche+ (+60°) (CH2–CH–CH2–CH) and trans (180°) (CH–CH2–CH–CH2) dihedral angles resulting in a down arrangement.37 For the dihedral angle definition, we make use of the IUPAC convention.38 The two types of helices are enantiomers, meaning that one is the mirror image of the other. It is worth noting that high energy helices do exist and can be defined as a sequence of (trans, gauche+) and (gauche−, trans).39 In order to represent an iPP helical structure at the CG level, the UA RH and LH helices are mapped based on the mapping proposed by Panizon et al.29 In Fig. 1d, it is clearly seen that the helical structure of the iPP chain persists at the CG level as well. The dihedral angle pattern for a CG helix should also be defined: a RH CG helix is denoted by a sequence of +100° dihedral angles, and a LH CG helix is denoted by a sequence of −100° dihedral angles. At the CG level there is no distinction between high energy helices and low energy helices.
On both levels (UA and CG), each force field describes explicitly all stretching, bending and torsional interactions. In a full melt simulation, non-bonded interactions would be calculated for atoms both along the same chain (intramolecular) and between different chains (intermolecular). Based on Flory's random coil hypothesis,49 however, a chain in a melt behaves as if it is subject to local interactions only. For the purpose of sampling single unperturbed chain conformations representative of the melt, nonlocal intramolecular interactions between topologically distant segments and all intermolecular interactions can be omitted. Thus, to sample a single unperturbed chain, one needs to define precisely what constitutes local interactions. Tzounis et al.44 have introduced a way of defining local interactions by systematically varying the range of nonbonded interactions along a chain (Fig. 3). This range is set by a parameter Δnpair, which determines with how many neighboring monomers a specific monomer interacts and is estimated based on an empirical criterion.44 Specifically, the value of Δnpair which maximizes the chain characteristic ratio is taken as representative of the unperturbed state.44 By incorporating these local interactions we are able to sample conformations representative of a melt by simulating just a single chain. Δnpair was set equal to 3 for both UA and CG simulations based on previous work on polypropylene.50
![]() | ||
Fig. 3 Definition of local interactions between monomers for single chain MC simulations. Δnpair corresponds to the number of repeat units to the left and to the right, with which a specific repeat unit interacts. All the atoms belonging to the central repeat unit, marked with a rectangle, interact with all the atoms belonging to the neighboring repeat units, marked with ellipses. Δnpair = 3 was determined appropriate for polypropylene.50 |
In Fig. 4a, the evolution of the dihedral energy of the RH and LH UA helices in time is plotted based on the united-atom representation. It is clearly seen that the dihedral energies for each helix fluctuate between the same energy levels. The fluctuations do not alter the helical structure. Next, the two helices were mapped onto a CG representation and a simulation of duration 400 ps was performed. In this case, the dihedral energy for the RH and LH CG helices came out significantly different (Fig. 4b). Even though the helicity is preserved, despite fluctuations, the RH CG helix is energetically higher than the LH helix. Therefore, the MARTINI force field29 favors one helical type (the LH helix) over the other.
![]() | ||
Fig. 4 Time evolution of the dihedral energy of single RH and LH helices comprised of 18 monomers each: (a) at the united-atom level using TraPPE-UA32 at 10 K (b) at the coarse-grained level using MARTINI29 at 10 K. In each panel the black line represents the RH helix, while the red line represents the LH helix. The total simulation time was 400 ps for the UA and the CG simulation. The scaling in each graph is different. |
To illustrate this further, we performed an MD simulation starting from an iPP CG melt of 40 chains of 50 monomers per chain. The iPP melt was first energy-minimized using the steepest descent algorithm51 until the maximum force was smaller than 100 kJ mol−1 mm−1. The algorithm converged after approximately 500 steps. Subsequently, a NVT simulation followed for 50 ns and, finally, a NPT simulation was carried out for 1000 ns. The final 500 ns were used for production. The simulation temperature (500 K) is well above the melting temperature for iPP oligomers that has been calculated at 411 K.52 The simulated radius of gyration (1.48 nm) and density (1048 kg m−3) are in excellent agreement with the results of Panizon et al.29 Note that this density deviates significantly from available experimental data (710–764 kg m−3).53
In the CG iPP melt simulated with the MARTINI force field29 we have analyzed the sequences of helical segments along a chain, in order to quantify the length of helical sequences formed spontaneously in the melt. More specifically, a CG segment was characterized as helical if |ϕ + 100°| < 25° (LH) or |ϕ − 100°| < 25° (RH), where ϕ denotes the dihedral angle based on the IUPAC convention.38 A sequence of such dihedrals along the chain is called helical sequence in what follows. In Fig. 5a, the distribution of lengths of helical sequences after 250 ns of NPT MD is plotted for the MARTINI29 model. The percentages in parentheses indicate the fraction of the atoms belonging to either a LH or a RH helix. The initial configuration produced by MAPS had equally distributed helices. As the CG simulation proceeded, the population of the LH helices increased drastically above 40%, while that of RH ones decreased to about 15%. During the whole simulation, the helical sequences remained shorter than 11, the threshold beyond which ordering is initiated.54
![]() | ||
Fig. 5 MD simulation using MARTINI29 at the CG level for an iPP melt of 40 chains of 50 monomers per chain at 500 K (a) distribution of helical sequence lengths extracted after 250 ns of NPT run. The RH helix is represented in orange and the LH helix in green. The percentages in parentheses indicate the fractions of atoms belonging to either a LH or a RH helix. (b) The probability density function (PDF) P(ϕ) of CG effective dihedral angles. The blue curve presents the results from the present simulations, while the dashed green curve presents the results from Panizon et al.29 The different probabilities are averaged over all configurations. |
The overall behavior of the helices is also displayed in the probability density function (PDF) of different CG effective dihedral angles plotted in Fig. 5b. Two maxima can be seen, one at ϕ = −100°, corresponding to the LH helix, and another one at ϕ = +100°, corresponding to the RH helix. The ϕ = −100° peak is clearly higher and thus, by using MARTINI,29 the LH helices in a melt are overpopulated.
![]() | (1) |
Following Brandrup et al.,57 the iPP chain radius of gyration Rg was estimated as
![]() | (2) |
In addition, the helical sequence length distribution for UA iPP is plotted (Fig. 6a), as also done by Yamamoto.39 The iPP configuration was extracted after 5 × 108 MC steps. It is important to note that the populations of both types of helices are equally distributed. Despite having a long chain of 2000 monomers, the helical sequences do not grow to more than 12 monomers per sequence. Thus, this UA model has no preference for any of the helical types, preserves the helical behavior, and reproduces the conformational properties of realistic iPP very well.
![]() | ||
Fig. 6 Monte Carlo simulations performed for a single chain of 2000 UA monomers at 500 K using the TraPPE-UA force field: (a) Distribution of helical sequence length. The RH helix is represented with orange and the LH helix with green. The percentages in parentheses represent the fraction of atoms that belong to either a LH or a RH helix. The configuration was extracted after 5 × 108 MC steps. (b) Dihedral effective potential U(ϕ) vs. CG effective dihedral angles. The single unperturbed chain configurations were initially sampled at the UA level using TraPPE-UA and then mapped onto the CG level to obtain the potential of mean force of CG torsion angles. The red curve represents the mapped CG trajectory from TraPPE-UA. A Ryckaert–Bellemans function was fitted to the potential of mean force U(ϕ), and is shown by the black dashed curve. The MARTINI dihedral effective potential U(ϕ) is shown by the blue curve.29 |
The next step was to use TraPPE-UA results as a reference in order to parametrize the dihedral effective potential for the new CG force field which would preserve the equal distribution of both LH and RH helices. For this purpose, the UA trajectory was mapped onto the CG level. From this mapped CG trajectory the PDF P(ϕ) of dihedral angles was extracted and, subsequently, was converted to the dihedral effective potential, U(ϕ) = −RTln(cP(ϕ)). The constant c was selected arbitrarily equal to 1000°, as it does not have any effect on the simulation. Then, a Ryckaert–Bellemans function (eqn (5)) was fitted to U(ϕ), as shown in Fig. 6b. In the same figure, the MARTINI dihedral effective potential U(ϕ) is also plotted, where at ϕ = −100°, corresponding to LH helix, the dihedral energy is significantly lower than the dihedral energy at ϕ = 100°, corresponding to RH. The new CG force field comprises effective stretching, bending, torsional and non-bonded interactions, with the parameters presented in Table 1. The reader should note that any coupling between successive effective torsion angles is ignored in the CG model. In other words, the total CG effective torsional potential of a chain is modeled as a sum of terms, each depending on a single effective torsion angle.
Stretching interactions (eqn (3)) |
l
o = 0.298 nm, kb = 48![]() |
Bending interactions (eqn (4)) |
θ 0 = 119° kθ = 78 kJ mol−1 rad−2 |
Torsional interactions (eqn (5)) |
a 0 = −8.08655 kJ mol−1, a1 = −2.88516 kJ mol−1, a2 = 16.61891 kJ mol−1 |
a 3 = 0.79249 kJ mol−1, a4 = −8.79016 kJ mol−1, a5 = 1.33047 kJ mol−1 |
Non-bonded interactions (eqn (6)) |
ε = 2.625 kJ mol−1, σ = 0.53 nm |
![]() | (3) |
![]() | (4) |
![]() | (5) |
![]() | (6) |
![]() | ||
Fig. 7 The characteristic ratio Cnvs. the number Nb of real chemical backbone bonds from Monte Carlo simulations for a single chain of 2000 monomers on both UA and CG levels at 500 K. At the UA level the TraPPE-UA was used, plotted with a red dash-dot line, while at the CG level the new HELIPP potential was used, plotted with black line. Results obtained with the existing MARTINI force field at 500 K29 are also shown as a green dashed line. Experimental values55 at 460 K are shown with the blue dots. |
Next, the PDF of CG effective valence angles is plotted in Fig. 8a. Looking at the HELIPP curve, a peak is observed at 119°, while the TraPPE-UA curve has a maximum at 109°, which is the most probable one, and a second weaker shoulder at 145°. As the ReB potential used for describing bending interactions58 is not able to reproduce a double-peak curve, the two curves do not match each other. The HELIPP curve actually replaces the UA bending angle distribution with a smoothed, unimodal approximate distribution.
In Fig. 8b, the PDF of CG effective dihedral angles is also plotted. In both curves, two maxima can be distinguished. Each peak represents a helical type, both are equally probable. In the regions [−180°, −150°] and [180°, 150°], representing the trans states, the HELIPP curve exhibits a slightly higher probability. On the contrary, in the region [−30°, 30°], which refers to the cis state, the HELIPP curve indicates a lower probability. CG beads, which have a parameter σ from the LJ interactions larger on average than the same parameter of the united atoms they are constituted from (see Appendix 1), experience slightly larger repulsion when approaching each other, which leads to the slightly less probable cis state.
Furthermore, we performed MD simulations of a CG melt of 40 chains of 500 beads each using this new HELIPP force field. In order to equilibrate the system, an energy minimization was undertaken until the maximum force was smaller than 10 kJ mol−1 nm−1, followed by a NVT MD of 5 ns and a NPT MD of duration 5 ms. The simulation temperature T = 500 K was well above the melting temperature of high molar mass iPP (T = 460 K) measured by differential scanning calorimetry.9 The simulation pressure P was set to 1 bar using the Parrinello–Rahman barostat.41 The energy minimization converged after 25000 steps.
We compare the MD simulation to the Single Chain Monte Carlo simulation with the same HELIPP force field. The PDF for CG effective valence and dihedral angles extracted from the MD simulation match perfectly (not shown here) the ones derived from CG single unperturbed chain MC sampling and shown in Fig. 8. The characteristic ratio from the CG melt MD simulation was estimated equal to 5.87 ± 0.09, close to the value 5.7 from the single CG chain MC simulation. Moreover, the density was measured as 730 kg m−3, in a good agreement with experimental results (710–764 kg m−3).53
The helical sequence length was also studied. Both types of helices are now equally probable with approximately 30% of the beads belonging to RH helix and 30% to a LH helix during the whole simulation. Neither of them grows longer than 12 monomers per helical sequence. So, HELIPP is able to reproduce LH and RH helices in equal amounts and that holds even for rather long simulations on the μs scale.
Subsequently, using this new coarse-grained model, the entanglement molecular weight Me was calculated by performing a topological analysis on the trajectory produced during the last 1 μs of the simulation. For the topological analysis, the CReTA (Contour Reduction Topological Analysis) algorithm59 was used. This algorithm maps an entangled polymer system to its corresponding system of primitive paths, having the same topological constraints as the original system. Me was computed equal to 4432 ± 121 g mol−1, whereas in experiments carried out for iPP melt Me was estimated between 5100–5500 g mol−160 and 6900 g mol−1.61 So HELIPP produces slightly more entangled system in comparison to the experiments, but no comparison was made with other simulations, as no data were available.
Finally, we have calculated the cohesive energy of the iPP melt which can be expressed by the Hildebrand solubility parameter δ.62δ can be obtained as
![]() | (7) |
In order to calculate δ for different temperatures, we performed multiple NPT MD simulations for 10 ns for a system of 40 chains of 50 monomers per chain using the HELIPP and TraPPE-UA force fields. The simulation temperatures were T = 400, 450 and 500 K, close to the temperature at which the HELIPP force field has been developed. The temperature dependence of the solubility parameter is plotted in Fig. 9. It can be seen that at T = 475 K, δ = 14.2 MPa1/2 obtained from HELIPP and δ = 10 MPa1/2 from TraPPE-UA force field. Comparing with other experimental (δ = 15.2 MPa1/2)63 and simulation data (δ = 14.8 MPa1/2)56 at T = 475 K, HELIPP seems to reproduce the solubility parameter much better than TraPPE-UA. At T = 298 K, the experimental value64 varies from 16.8 to 18.8 MPa1/2, while the present estimation with HELIPP is 15.5 MPa1/2, showing a larger deviation at this temperature. So HELIPP reproduces well the solubility parameter above the melting temperature, even better than simulations at the UA level.
![]() | ||
Fig. 9 The solubility parameter plotted at different temperatures for an iPP of 40 chains of 50 monomers per chain. The solubility was calculated both at the UA level with TraPPE-UA (black squares) and at the CG level with HELIPP (green dots) using molecular-dynamics simulations. The experimental values obtained from Maier et al.63 and from van Krevelen64 are plotted with a magenta inverse triangle and a red rhombus, respectively, while atomistic simulation results obtained by Logotheti et al.56 are plotted with a blue triangle. |
The representation of the RH and LH helices at the CG level and preservation of the handedness have been proved essential for the modeling of iPP melts. The existing MARTINI CG force field for iPP29 is not able to reproduce correctly the two different helical types and overestimates the stiffness of the chains and the density. By introducing the CG iPP helices and tailoring the dihedral effective potential so that none of the helices is favored over the other, we managed to have a more accurate estimation for the iPP chain characteristic ratio, C∞ = 5.7. Moreover, by re-adjusting the Lennard-Jones σ based on a geometric approximation, as discussed in Appendix 1, the density ρ = 730 kg m−3 was also improved. This force field was applied to an iPP melt but seems promising also for iPP semi-crystalline morphology. Below we present a snapshot of a lamella formation after an iPP coarse-grained system was submitted to stretch-induced crystallization using the HELIPP. Specifically, a system consisting of 20 chains of 1000 CG monomers each was equilibrated in the NPT statistical ensemble at T = 500 K and P = 1 bar using the Nosé–Hoover thermostat and the Parrinello–Rahman barostat.41 Then, it was cooled down to T = 380 K with cooling rate 10 K ns−1 and stretched at the same temperature to 5 times its initial length in the NėxxPyyPzzT statistical ensemble with stretching rate ėxx = 107 s−1. Finally, it was annealed in the NLxxPyyPzzT statistical ensemble at T = 380 K for 400 ns (Fig. 10).
![]() | ||
Fig. 10 A lamella snapshot of a CG system of 20 chains of 1000 CG monomers each using the HELIPP after a stretch-induced crystallization simulation. |
iPP crystals consist of lamellae which are formed by RH and LH helices in equal amounts. This characteristic is possible to maintain now at the CG level, although attention should be paid to two main aspects: (1) the presented force field has not been tested thoroughly for temperatures below the melting point. (2) At these low temperatures the reproduction of the different iPP crystal phases is crucial, and also has not been tested yet.
Single unperturbed chain MC simulations have been used for the fast and effective estimation of the stiffness of polymer melts but, as far as we know, not for the development of a force field up until now. The CG stretching, bending and torsional interactions can be accurately derived from a single chain MC simulation. In the present study we tailored only dihedral and non-bonded interactions, based on the approach implemented by Panizon et al.29 For bending interactions we adopted the Restricted Bending (ReB) potential. In addition to the work presented above, we have performed CG simulations with tabulated potentials both for bending and torsional interactions which perfectly match the PDF of the CG effective valence and dihedral angles derived from TraPPE-UA single chain MC simulation (Fig. 8a). The results using tabulated potentials can be found in the ESI.† The main reasons we did not use these potentials for the production simulations were: (1) tabulated potentials can not be easily implemented in different simulation packages. (2) The maximum timestep we had to use with the tabulated potentials was 5 fs, while with ReB we use a timestep of 20 fs. For these reasons, we incorporated ReB and Ryckaert–Bellemans potential into HELIPP. As far as the non-bonded interactions are concerned, here we used a geometric approach which proved to be effective (see Appendix 1). This geometric approach is recommended to be tested also for other polymer systems. Another method that could be used for mapping non-bonded interactions is iterative Boltzmann inversion,65 but, overall, the single chain Monte Carlo method can serve as a new coarse-graining method for mapping a variety of polymer melt models to the coarse-grained level.
Finally, we demonstrate the computational significance of the coarse-graining by comparing a UA and a CG simulations in terms of time. Coarse-grained models are far more efficient computationally than models with higher resolution. To check this, we simulated 20 iPP chains of 500 monomers per chain at the UA level using TraPPE-UA and at the CG level using the HELIPP. The simulations were carried out on a cluster made of 504 nodes comprised of 128 AMD Rome 7H12 processors, each with 2.6 GHz and 2 GB per core. The estimation made on 32 cores on 1 node for the CG performance was 4000 ns day−1, while for the UA the performance was 100 ns day−1. So, the CG iPP model performs 40 times faster than the UA iPP model.
rCH = {0,0,0} |
rCH2 = { −0.129, 0, −0.084} |
rCH3 = {0, 0.122, 0.094} |
rCG = {0, 0.044, 0.001} |
In Fig. 11, the radius of the coarse-grained LJ sphere σCG/CH2 that is tangential to the LJ sphere of the methylene group and has its center at the position of the coarse-grained bead is calculated,
σ CG/CH2 = ‖rCG − rCH2‖ + σCH2 = 0.56 nm |
Similarly, σCG/CH3 = 0.49 nm and σCG/CH = 0.51 nm. All the distances at which the LJ potential becomes zero for the united atoms are given by Pütz.33
Finally, σCG is calculated as the average of:
Footnote |
† Electronic supplementary information (ESI) available: Molecular-dynamics simulations performed using tabulated coarse-grained valence and dihedral angle potentials for isotactic polypropylene melt. See DOI: 10.1039/d2sm00200k |
This journal is © The Royal Society of Chemistry 2022 |