From time dependent incorporation of molecular building blocks to application properties for inorganic and organic three-dimensional network polymers

The three-dimensional configurational arrangement of natural and synthetic network materials determines their application range. Control of the real time incorporation of each building block, hence, all functional groups is desired so that we can regulate macroscopic properties from the molecular level onwards. Here we interconnect kinetic Monte Carlo simulations from the field of chemical kinetics and molecular dynamic simulations from the field of physics. We visualize for (in)organic network material synthesis how the initial building blocks interact timewise and spatially, accounting for variations in inter- and intramolecular chemical reactivity, diffusivity, segmental compositions, branch/network point locations, and defects. We use the kinetic and three-dimensional structural information to construct structure-property relationships based on molecular descriptors such as the molecular pore size or dangling chain distribution, differentiating between ideal and non-ideal structural elements. The generic nature is illustrated by constructing such relationships for the synthesis of organosilica, epoxy-amine and Diels-Alder

Essential for network materials is the degree of three-dimensional (3D) configurational ordering or arrangement, as it determines their macroscopic property and application range. 1-4 Highly ordered carbon atom (C) arrangements possess e.g. an exceptional behavior with diamond (sp 3 hybridization; Fig. 1a) and graphite (sp 2 hybridization) respectively known for their hardness and lubricating potential. [5][6] Another example is a synthetic zeolite with its ordered molecular pore structure allowing shape selectivity of catalyzed chemical reactions. 7,8 Fig. 1: Network materials with decreasing order of their 3D structural configuration. a perfectly ordered carbon (C) atoms in sp 3 hybridization, as in diamonds. b highly ordered silica atoms with oxygen (O) bridges in sp 3 hybridization with a limited number of structural defects, as in solar cells. Examples of such defects are unreacted functional groups (FG 1 = OEt, FG 2 = OH; Et= ethyl), leading to incomplete crosslinking points (CP's; CPi = CP with i crosslinks). c highly random configurations of (co)monomer units (mixture of sp3 and sp2 hybridization) containing atoms such as C, O and nitrogen (N), and functional groups leading to network molecules with different segment lengths and many structural defects (e.g. dangling chains: short and long branches (S/LCB's), and inverse monomer insertions), as in hydrogels; examples of molecular pores (mP's) are colored in c in light orange; the building blocks are symbolized by colored spheres that represent atoms (in a and b) and a chemical moiety such as a targeted comonomer unit A/B or side reaction defect D (in c).
Less ordered configurations also exist with a prominent role for natural polymers such as cellulose and lignin forming plant cell walls, 9 and network polymers synthesized from oil-derived or renewable (co)monomers. [10][11][12][13][14][15][16] For example, organosilica networks as employed for solar cells, separation

Conceptual framework from molecular to material and ultimately application level
Experimental methodologies are typically delivering information at scales lager than the molecular one with chromatography and mass spectrometry not automatically applicable and analytical methods such as solid-state nuclear magnetic resonance spectrometry providing only bulk information. 17,19,20 Sometimes average characteristics such as the overall comonomer content and the number of branching points per (co)monomer unit can be measured for specially modified reactants 19  Only rarely 2D connectivities of the CP's or clustering in (graph) node notation have been aimed at 21,22 but in any case a 3D network configuration has not been calculated. Most promising is our recently developed matrix-based kinetic Monte Carlo (kMC) methodology, which upon experimental validation considering average characteristics (e.g. the number average chain length and dispersity) 16,23 and inputted individual kinetic parameters from independent experimental techniques, 24 allows to access (co)monomer sequences of a representative number of individual (e.g. 10 4 -10 5 ) linear or slightly branched chains. 14,25,26 Also in the physics research field, simplifications have been made with force-field-based molecular dynamic (MD) simulations mostly predicting ideal defect-free 3D network structures at high monomer conversion often neglecting molecular kinetic constraints. A key MD assumption is the free movement of atoms or atom groups until the thermodynamically equilibrated configuration is reached, typically recognizing a limited number of chemical bond types and thus ignoring many side (e.g. intramolecular) reactions and diffusivity constraints. 10,[27][28][29][30][31][32] For the equilibrium settling, MD simulations implicitly assume a generally incorrect large real time. Free radical based organic networks are e.g. produced on a minute to hour scale with altering dynamic viscosity (η), leading to time and chain length dependent diffusivities. This prevents the final configuration from reaching thermodynamic equilibrium for the complete 3D network, causing structural defect formation which influences the material properties. 3,33 Fig. 2: Main concepts of framework to design network polymer synthesis starting at the molecular scale; matrix-based kinetic Monte Carlo (kMC) and molecular dynamic (MD) simulations are interconnected to visualize the incorporation of building block by building block so each functional group (FG) at any synthesis time t; illustration for organosilica network synthesis (Fig. 1b) commencing from tetraethylorthosilicate (TEOS). a molecular scale with chemical rate coefficients for main and side reactions (kchem values; Arrhenius parameters A and Ea; reaction distance σ; polymerization temperature Tp and pressure pp). b chemical kinetics with calculation of inter/intramolecular rates (rinter/intra; (mol L -1 ) s -1 ) from concentrations of reactants, intermediates and products (C values) accounting for interplay of inter/intramolecular reactions and diffusivities (D: diffusion coefficient; η: dynamic viscosity; kdiff: diffusion rate coefficient). c kMC data storage and update of compositions, including segments, structural defects, and connectivities of crosslinking points (CP's). Combined with MD data on bond/dihedral angles and bond lengths 3D configurations for each molecule result at any t accounting for non-idealities thus defects. d in silico derived molecular network characteristics up to topological scale at any t; examples are CP distribution (at least 3 crosslinks; CP3+CP4), molar mass distribution, and molecular pore size distribution. e characteristics from d serve as input to construct structureproperty relationships that evaluate the network material performance, with a possible differentiation between ideal and non-ideal structural network elements. Example of hydrophilicity calculation based on the presence of OH functional groups near the surface (see Supplementary Methods).
In the present work, we connect the fields of chemical kinetics and physics by developing a generic framework according to the principles in the Methods section to visualize the spatial kinetic growth of network polymers at the molecular level so that we can fully grasp the birth of their 3D nature and the further inter/intramolecular interactions and chemical modifications. Matrix-based kMC and MD simulations are uniquely combined to store at any synthesis time t the molecular information on individual (segment) compositions, functional groups, bond lengths, bond and dihedral angles, and structural defects. As shown in Fig. 2, the step by step incorporation of each building block thus functional group is accounted for while grasping at various length scales (top row) the essential fundamental phenomena (bottom row) through step-wise experimental validation. Focus is on the synthesis of organosilica networks (Fig. 1b), commencing with a sufficiently large number of tetraethylorthosilicate (TEOS) molecules as kMC ensemble (e.g. 10 5 -10 6 ).
As shown in Fig. 2a, the initial TEOS molecules can be hydrolyzed multiple times (top of box) and take part in subsequent condensation reactions, resulting in the formation of Si-O-Si crosslinks (bottom of box). In general, a differentiation is made between inter-and intramolecular crosslinking reactions altering the local network make-up. The intermolecular reactions bridge network molecules so that a large(r) network molecule is formed but also deplete the linear molecules by their inclusion as dangling chains in an already existing network molecule or by the creation of branched still non-crosslinked species. The intramolecular reactions increase the crosslinking degree within a given network molecule through cyclization. They are more likely at higher polymer network yields 34 and lead to structural defects often denoted as internal loops. 35 Every reaction step in Fig. 2a, as defined by the involved functional groups, is characterized by Arrhenius parameters (A and Ea values) so that at a given polymerization temperature (Tp) and pressure (pp) chemical rate coefficients (kchem values) describing intrinsic chemical reactivities can be obtained. Regression to experimental data using monofunctional analogues or limited to low crosslinking yields is most suited to determine kchem values, linking via kinetic modeling analysis the dominance of specific reactions to specific time regimes for certain experimental responses. The number of unreacted groups of a network moiety (or building block) can influence kchem values of the remaining functional groups, as shown in Supplementary Figure 11.
Besides the kchem values in Fig. 2a the concentration (C) variations of all molecule types are required thus from monomer to dimer and ultimately network molecules, addressing the spectrum in functional group types and combinations. As shown in Fig. 2b, inter-and intermolecular reaction rates (rinter/intra) are calculated that are corrected for diffusional limitations. For larger η values, the intermolecular reactivities can be even fully dictated by chain length and CP dependent diffusivities (Dinter values). 14 Intermolecular apparent rate coefficients (kapp instead of kchem values) are therefore introduced to grasp the interaction of intrinsic kinetics and diffusional constraints. 36 Intramolecular Si-O-Si crosslinking in turn becomes impossible if the distance between the functional groups is too large (red vs. green arrows in Fig. 3b), since the intramolecular diffusion coefficients on the functional group level (Dintra values) are then too low. Fundamental distance rules accounting for rigidity, as codictated by the number of CP's between functional groups, are therefore introduced.
During the calculation of rinter/intra values, which are function of kchem and C variations, we keep track of the connectivity history of all molecules in the kMC ensemble through a composite topology matrix T.
In this way, we can retrieve at any t extensive molecular distributed data regarding (i) the composition of remaining and formed linear/branched molecules, and (ii) the segment compositions, including structural defects, and the connectivities of each separate network molecule ( Fig. 2c; top). Combined with MD simulation input a 3D representation of these individual (network) molecules is within reach.
MD simulations allow to access the spectrum of bond and dihedral angles, and bond lengths that are thermodynamically feasible, 10 as illustrated in Fig. 2c (bottom). Hence, by following each connectivity in the kMC matrix T the local environment can be scanned and a feasible thus stable 3D structural element can be generated if no structural defect is detected. Upon doing so alternative 3D visualization rules are needed to ensure the correct representation of the non-ideality in the 3D network structure. The unique combination of kMC and MD data (Fig. 2a-c) therefore enables to visualize the 3D configuration of network molecules as a function of t covering both ideal and non-ideal building block incorporations.
Notably MD simulations deliver input data for the kMC algorithm and, hence, any further development in the physics field, e.g. the improved description of polymer-medium interactions, can be directly translated in the current multi-scale development. (ii) mass chain length distribution (CLD: monomer unit: building block), including in inset for t2 the distinction between sol (linear/loosely branched chains) and gel (precipitated molecules) with for molecules as cut-off to belong to gel the ratio of number of CP's to number of monomer units larger than 0.5; (iii) number molecular pore size distribution with the pore size related to number of monomer units in the pore; 10% error bar.
As highlighted in Fig. 2d-e, the detailed 3D molecular information from Fig. 2a-c can be utilized to construct fundamental structure-property relationships in which the dependent variables are detailed molecular network descriptors coming from explicit distributions. Examples of such distributions are the CP connectivity, dangling chain, functional group distance, and molecular pore size distribution, with a possible differentiation between the contribution of ideal and non-ideal structural network elements. Upon recording experimental macroscopic responses for instance at various t values we can thus identify the main network molecular descriptors for material design and future development so that the synthesis conditions can be tuned toward a controlled balance of main and side reactions, and of chemistry and viscosity variations. This t dependency constitutes the strength of the framework, as macroscopic property design is achieved whenever desired not restricted to current simplified cases of limiting crosslinking or ideally constructed high yield (theoretical) polymer networks. We can see in Fig. 3 that at low yield a gradual transition from a linear to a slightly crosslinked polymer is obtained that is fully converted in a network material at high yield, as evidenced by the increasing fraction of network molecules with more than 3 CP's containing at least 3 crosslinks. A high degree of structural heterogeneity is always established. At small times (t1) the CP connectivity distribution is lessdefined particularly for molecules with less crosslinking, whereas at larger times (t2 and t3) the nonideality is also established at the tail as evident upon comparison with the inset results (fintra=0) displaying different contributions of more crosslinked molecules. Such tail variation is also valid for the chain length distribution, indicative of a continuous compositional variation. The molecular pore size distribution first develops a bimodal character (t1 to t2), as fintra increases from 0 to a significant value so that more pores are formed from smaller molecules, to then at larger time (t3) loose this bimodality accompanied by a shift to the left as many smaller pores are created by a dominance of intramolecular reactions. Similar non-trivial molecular property dynamics are depicted in Supplementary Figure 20, addressing the CP distance distribution, and the average crosslinking density and number/mass chain length variation. Hence, t dependencies strongly alter the molecular build-up so that synthesis conditions are essential in network material design.

Discussion
As shown in Fig. 4, even upon focusing only on conventional macroscopic properties such as the contact angle, the storage modulus, and the swelling degree bearing in mind applications such as membrane separation, 12 coatings, 38 and drug delivery, 13 the developed framework provides added value in molecular scale interpretation and design. In Fig. 4a, emphasis is still on the organosilica case dealing with 3 high yield (70%) 3D ensembles as in silico generated for 3 initial sets of conditions that differ in the initial H2O amount but with otherwise conditions as in Fig. 3. For one such ensemble (initial water to TEOS molar ratio of 4) the 3D structure is depicted in the left part of Fig. 4a. Post-processing using density-based outlier detection, 39 as explained in the Supplementary Discussion, allows calculating the number of OH groups near the surface so that the hydrophilicity can be simulated as normalized with respect to the total number of functional groups. These simulated hydrophilicities are displayed in the right part of Fig. 4a (left axis; black symbols) and correlate with measured contact angles (right axis; orange symbols). 40 A larger simulated hydrophilicity implies a better interaction with the surface thus better wetting so that the experimental contact angle decreases. A closer inspection of Fig. 4a, however, reveals that the conventional contact angle is only sensitive in the region of lower to intermediate hydrophilicities (up to 60%), which is associated with initial water to TEOS molar ratios smaller than or equal to 4. The simulation results although show that higher initial molar ratios are experimentally useful if one wants to maximize the availability of OH groups near the surface (> 60%), as for instance crucial for molecular scale driven separation. 12,41 In Fig Table 1). Relation of initial molar ratio of water (H2O) to tetraethylorthosilicate (TEOS) for synthesis (298 K) and simulated hydrophilicity at 70% yield (density-based outlier detection: 39 Supplementary Methods section for relative contribution of OH near surface) to understand the change in experimental contact angles 40 (orange; right axis) b application for epoxy-amine curing case (network chemistry 2; initial and model parameters: Supplementary Table  2). Explanation of temporal variation of experimentally recorded variation of storage modulus (orange; right axis; this work; initial storage modulus due to the flexible cuvette) based on simulated concentration of network molecules with more than x (3-6) crosslinking points (CP's) containing at least three crosslinks (left axis; CP3 + CP4) and relative importance of intramolecular reactions (upper graph; fintra); c application for Diels-Alder based case (network chemistry 3; initial and model parameters: Supplementary Table 3). Relation of stoichiometry for initial building blocks for synthesis at 353 K and simulated number average of CP's with at least 3 crosslinks (black; left axis) and the number average molecular pore size (blue; left axis) at final yield to understand experimental swelling data (orange; right axis; this work; toluene; Supplementary Figure 32a Hence, by accessing the contribution of specific molecules either with a given number of CP's of a certain type or a certain level of intramolecular defects we can understand the traditional variation of the storage modulus in Fig. 4b form the molecular scale onwards. The current work thus opens the pathway to control the contribution of specific molecules to the overall macroscopic behavior in this case material stiffness.
In Fig. 4c (right part), we further relate molecular scale information to the swelling degree, which is a key conventional polymer network macroscopic property in the fields of hydrogels and drug delivery. 42 We focus on previously established rapidly linking ambient temperature Diels-Alder chemistry by the Barner-Kowollik group 43 Moreover, the left part of Fig. 4c shows the 3D structural incorporation of OH groups as originally present in the bifunctional monomeric building blocks (colored in purple; molar ratio of 2). Such groups are a readily accessible handle for chemical modification (e.g. drug loading) but we cannot assume that they are characterized by a controlled functional group distance distribution, which is a typical assumption if one cannot make a dedicated link to the prior synthesis on the individual functional group level. As shown in Fig. 5a, for the idealized case of fintra = 0 (open purple symbols; molar ratio of 2; additional simulation results in Supplementary Figure 24), a broad OH group distance distribution results representing a wide variety of network molecules. In practice, as only accessible with our framework, intramolecular reactions occur leading to loop formation so that the OH group distance distribution becomes more complex with a shift to the left and a switch to a bimodal behavior (full purple symbols in Fig. 5a).

Outlook
Our framework provides a detailed characterization and design of both organic and inorganic polymer network materials from the molecular to the application level. We can truly evaluate if a certain network synthesis is close to the targeted structure or not, and we can run the platform over a wide range of reaction conditions to identify the suited reaction conditions, which is purely experimental extremely cost-intensive. The current work also highlights that we need to revise previously developed theories based on average characteristics not explicitly coming from distributions and ideal network representations, as network molecular species are too heterogenous in composition and structure to be described by a single ideal or average molecule.
We can thus perform analysis at the functional group level, which is relevant for existing synthesis platforms but also emerging fields, for instance network materials from supramolecular chemistry in which it can be expected that only guest-host interactions for individual well-positioned functional groups lead to advanced sensing and molecular recognition. 45,46 Another example are the more recently developed dynamic recyclable networks in which material reshaping can be fully designed if one knows at which 3D position the exchangeable functional groups are located. 47,48 In general, based on any desired molecular distributed descriptor, we are capable to fundamentally understand which type of molecules or structural elements allow to design which material property.

Experimental framework
To illustrate the generic character of the developed framework focus is on three network chemistries: (i) organosilica network synthesis (chemistry 1), epoxy-amine curing (chemistry 2), and Diels-Alder network formation (chemistry 3). As explained below for chemistry 1 literature data were employed and for chemistry 2 and 3 experimental data were recorded in the present work.  For chemistry 2, bisphenol-F diglycidyl ether (DGEBF; Sigma-Aldrich, used as received) and ethylene diamine (EDA; Sigma-Aldrich, used as received) were used as precursors to perform the epoxy-amine curing at 298 K. The curing was analyzed with near-infrared (nIR) spectroscopy employing an in-house available Perkin Elmer Lambda 900 UV/VIS/NIR Spectrometer. Both the DGEBF and EDA precursor were thoroughly mixed with equimolar functional group concentrations before being brought into a nIRtransparent cuvette. nIR spectra of the reaction mixture were subsequently recorded at frequent time intervals from 4000 to 10000 cm -1 with a resolution of 3 nm. The concentration of primary, secondary and tertiary amine, epoxide and hydroxyl groups were analyzed as a function of synthesis time using Spectragryph software, according to the procedure previously reported. 49 In this procedure, the spectra are internally normalized by the aromatic ring peak at 4620 cm -1 and then the primary amine and epoxide For the small molecule experiment, sorbic alcohol (reactant 1; 0.549 mg, 1 eq.) and the RAFT agent (reactant 2; 2.02 mg, 1 eq.) were dissolved in toluene (0.5 mL; Thermo Fisher Scientific, after drying and purification with SP-1 Stand Alone Solvent Purification System LC Technology Solutions Inc.).
The vial was closed and purged with a constant argon flow for 5 minutes. The reaction mixture was heated up to 353 K and maintained at that temperature. Conversion was determined using 1 H NMR spectroscopy (600 MHz; chloroform-d3, Sigma-Aldrich).
Network formation (phase 3b) was subsequently performed considering a tetrafunctional and bifunctional monomeric building block related to phase 3a, respectively denoted as reactant 3 and reactant 4. Benzene-1,2,4,5-tetrayltetrakis(methylene)tetrakis((diethoxyphosphoryl)-methanedithioate) (reactant 3: tetrafunctional monomer) was synthesized as follows. In a dry 10 mL Schlenk tube with a magnetic stirrer, sodium hydride (suspension in paraffin oil; 90.0 mg, 1.74 mmol, 1 eq.; Sigma-Aldrich) was added and vigorously mixed with 3 mL heptane (Thermo Fisher Scientific). After a few minutes, the liquid phase was removed and heptane was added again. The procedure was repeated three times until a white powder of sodium hydride was obtained. The residual liquid was removed under vacuum conditions and dry THF (2 mL) was added. To that suspension a solution of diethyl phosphite (240 mg, 1.74 mmol, 1 eq.; Sigma-Aldrich) in THF (1 mL) was added dropwise. After addition, the mixture was allowed to stir for 1 h and then heated to 339 K for 2 minutes. After cooling the reaction mixture in an acetone-liquid nitrogen bath, CS2 (660 mg, 8.70 mmol, 5 eq.; Sigma-Aldrich) was added over a period For network formation, tetrafunctional monomer (reactant 3; 26.5 mg, 0.027 mmol, 1 eq.), different equivalents of bifunctional monomer (reactant 4; 6.8 mg, 1 eq.; 10.2 mg, 1.5 eq.; 11.9 mg, 1.75 eq.; 13.6 mg, 2 eq.) and toluene (0.025 mL) were added in a 1 mL crimp vial. The vial was closed and purged with a constant argon flow for 5 minutes. Afterwards, the reaction mixture was heated up to 353 K for 24 h and analysis was performed. The networks did not dissolve in toluene and swelling tests were done in that solvent for 3 h. The percentual degree of swelling was calculated according to the ratio of 100 (ws-wd) and wd, with ws the mass of the sample after 3 h swelling and wd the mass of sample prior to swelling. The related experimental data are included in Fig. 4c. Images of the synthesized networks are shown in Supplementary Figure 32a  Under stoichiometric conditions a yellow network is obtained, while going more and more to offstoichiometric conditions (excess of tetrafunctional linker) a pink color is observed due to the presence of unreacted dithioester moieties of the tetrafunctional linker.
Moreover, in phase 3c, the network synthesis was also performed with Fmoc-loaded bifunctional monomer (reactant 5; short notation bifunctional linker-Fmoc, with the main results shown in Figure 5.  Figure 33). The respective network (r = 2; 1.3 mg or r = 1.75; 1.5 mg or r = 1.5; 3.2 mg) was added into the cuvette and washed with DMF (three times). Afterwards, the cleavage solution (2 mL; 2% 1,8-diazabicyclo [5.4.0]undec-7-ene (DBU; Sigma Aldrich) in DMF) was added, the cuvette was inverted twice to ensure homogenic mixing and the UV-Vis spectrum was immediately measured (first measurement after ca. 10 s). This was repeated until the final measuring time was reached.

Modeling framework
As conceptually illustrated in Fig. 2  given by introducing boxes with labels from A to N and providing a detailed discussion of the key implementation steps for each of these boxes.
In this subsection, the essential information for the most critical implementation characteristics of the interactive flowsheets is provided to enable a general reader to grasp the overall technical innovation.
In the core matrix-based kMC simulations, reaction event per reaction event is sampled based on reaction probabilities as represented by means of a cumulative probability curve, starting from a sufficiently large number of initial molecules to ensure numerical convergence as shown in Supplementary Figure 16 for chemistry 1 selecting some key molecular characteristics. The update of the reaction probabilities after the occurrence of a reaction event requires the discrete updates of the reaction rates, which are directly related to (i) the chemistry platform selected defining the potential reaction types between the functional groups present and (ii) the reaction mixture composition defining the (reactant) concentrations related to these reaction types. In the present work, most simulation emphasis is on the organosilica network synthesis case ( Fig. 3 and available, which is e.g. the case if a structural defect is encountered that has been created due to a chemical reaction omitted in typical MD simulations, as illustrated for the organosilica case in Supplementary Table 5.
Due to the discrete nature of the developed modeling framework we can select the simulation times when we want to switch to the more detailed part regarding molecular specifications and configurations.
The more frequently this switch is performed the more detailed the modeling outcome becomes albeit at computational cost. The need for a switch is directly related to the overall kinetics as for instance deducible from the yield profile with a steep variation implying strongly alternating kinetics. In the present work, we have carefully analyzed the switching times to guarantee numerical convergence and sufficient molecular detail on the one hand and an acceptable computational cost on the other hand. For defect free structures, using supercomputer platforms as in the present work, we have a couple of hours as simulation time. With more structural defects formed the simulation time increases and we have with current supercomputer capacities a simulation time of a few days. Of course the MD input data need to be acquired as well and there we have simulations times on the scale of hours, days or weeks depending on the required potential field.
To highlight the general accuracy of the developed modeling framework we showcase in Supplementary In addition, in the Supplementary Discussion, all the required steps for the determination of the chemical and diffusional limitations related input parameters are covered for each chemistry. A general research strategy is followed in which smaller synthesis times or reaction systems with monofunctional analogues are considered to determine chemistry related kinetic parameters free from the impact of diffusional limitations. The larger synthesis times are subsequently employed to tune diffusional limitations related parameters with as input the previously determined chemistry related kinetic parameters. Hence, inherently parameter correlation is minimized also bearing in mind that specific experimental responses or variations within a selected response are considered. This is further illustrated in Supplementary