Al-driven cement functionality by manifold structuring & disorder

Abstract

Cementitious calcium–silicate–hydrate (C–S–H) gel structures with differing aluminium (Al) content and structuring were characterised by a series of ab initio molecular dynamics (AIMD) experiments. Differing Al-binding modes and coordinations (4, 5, 6) were incorporated into dynamic simulations of various C–S–H bulk structures (layered, non-layered, anhydrous porous, hydrated porous) to resolve Al's configurational preference as well as any effects on Si/Al chain structure, including Lowenstein type ‘Al-avoidance’ (preference for non-neighbouring Al atoms) and Pore sizes. Emphasis was placed on modelling fault lines in C–(A–)S–H structures, where local order is thought to be lost. These regions tend to have high concentrations of micro stresses that potentially induce crack propagation and would thus possibly benefit from Al-toughening. Principal aspects of Loewenstein's Al avoidance are not operative in the majority of cases, with alternating Al atoms increasing chain length and branching giving rise to larger pores. O–Al–O and Al–O–Al bond angles were shown to be more flexible than their Si counterparts (O–Si–O, Si–O–Si), the former Al-units effectively acting as atomic hinge points. Al-coordination was shown to shift dynamically throughout each simulation, generally occupying tetrahedral Al-IV geometry with some preference for penta-coordinated Al-V at Q₂ and Q₃ sites. Loewenstein's principle of opposite coordination modes for adjacent Al are observed; specifically, Al neighbours will seldom hold the same coordination geometry. Through this approach we have demonstrated structuring of C–(A–)S–H containing stable hydrated pores, in agreement with experimental trends arising from empirical NMR and small angle neutron scattering (SANS) measurements in established works.

---

1. Introduction

Production of modern cement is currently responsible for ∼5–8% of global CO₂ emissions and rising, in addition to it suffering from low durability and longevity as well as issues with creep, shrinkage and chemical instability.^1,2 These shortcomings must be addressed to evolve production of more sustainable cement blends. Ancient Roman cements are much tougher, having stood the test of time for more than 2000 years with no steel reinforcement (rebar). While this is in part attributed to the exceptional architecture and compaction methods used at the time, toughness also arises from content of pozzolana in the cement blends and phase compositions resulting from heating–cooling and composition, as in other materials.^3–8 This pozzolana is a natural alumino-siliceous material with raised aluminium (Al) content that is “a siliceous or siliceous and aluminous material which, in itself, possesses little or no cementitious value but will, in finely-divided form and in the presence of moisture, react chemically with calcium hydroxide at ordinary temperatures to form compounds possessing cementitious properties”.^7,9

Functional shortcomings of Portland cement originate at the atomistic level, within the calcium–silicon–hydrate (C–S–H), the main binding phase that essentially cements all constituent component particles together and thus yields the bulk cement's mechanical properties. Al-containing cements, including Roman cements, have Al-enriched C–S–H phases, or C–A–S–H, the latter being inherently more amorphous and known to help reinforce underwater concrete.¹⁰ C–A–S–H forms ‘platy-intergrowths’ of mesoscale sizing, that act to obstruct crack propagation along otherwise stark interfacial zones.^4,5 Even C–S–H defined phases have ∼2.5 wt% Al, with fractions up to Al/Si ≈ 0.1 (10% Al), substituting for Si mainly at the bridging (Q₂) position.¹¹ Forcefield models of C–A–S–H have shown that minimal Al uptake increases and maintains mean chain length even as tensile loading is applied.¹² The correlation between Al-content and tensile strength justifies the lack of rebar in ancient roman cements, which the latter would otherwise serve to relieve tensile loading. Recent work has also revealed the mechanism by which Al eases silicate chain oligomerisation, employing a molecular cluster approach.^13,14 Maximising toughening via Al-substitution (e.g. of Si) requires knowledge of the best configuration (position and Al–Al ordering) and conformation (local Al structuring, coordination and ‘behaviour’ with respect to pores and neighbouring atoms) which has prompted further studies on C–S–H gels with raised aluminium contents.^15–17 Loewenstein's law of Al avoidance must also be considered, which states that there is a thermodynamic preference for Al–O–Si–O–Al (Al avoidance) versus Al–O–Al–O–Si (Al-proximity).¹⁸ To date, this rule is widely accepted in most aluminosilicates with particular exceptions in calcium-alumino-silicates.¹⁹ A study by Pegado et al. employed CP2K's QUICKSTEP simulation module to very accurately depict a C–S–H model based on the mineral tobermorite and perform static energy calculations with mono and disubstituted Al atoms²⁰ Although Tobermorite at-best represents an idealised and moderately representative scaffolding of C–(A–)S–H's true structure (which includes quite a bit of amorphous and disordered structuring), they revealed Al-avoidance was in fact favourable, yet the change of relative energy (ΔE) plateaued at large Al–Al distances. In conclusion, the authors stressed the need for added complexity in the models, particularly Al-coordinations > 4 (Al-5, Al-6) as well as the modelling of interlayer and edge regions.

Given Al's integral role in chain formation, this current work employed dynamic simulations to determine how Al-avoidance contributes to the final structuring of the Si/Al chain network as well as seeing if the energy penalty of Al-avoidance emerges as a distinct limiting factor on C–(A–)S–H structures.

1.1. Order–disorder

While it is known that C–S–H is a semi-crystalline phase, its exact atomistic structure remains unresolved. As C–S–H is the combination of several hydrated phases and not their respective powder precursors, there is little homogeneity, rendering typical diffraction techniques less effective in resolving long range structure. Nevertheless, X-ray diffraction and microscopy have been used to draw a comparison between C–S–H and the crystalline mineral Tobermorite. This mineral forms the basis of the Feldman and Sereda model for C–S–H, which in turn has been adapted and evolved using modern spectroscopic and computational techniques.²¹ This model originally evolved from the Powers model (based on cellulose) that describes C–S–H as layers of SiO₄ tetrahedra trapping tubes of interlayer and intralayer water and Ca.²² C–A–S–H have shown maximum Al-uptake at low Ca/Si ratios (Ca/Si ∼ 1.0), while SANS experiments by Allen et al. have been very clear in resolving that the same ratio for C–S–H is 1.7 and the density of C–S–H is approximated by the same study to ∼2.6 g cm⁻³.²³ SANS and SAXS have little to no resolution down the c/z axis of hydrated cement due to both the sheer thickness and the likelihood of any layering being badly aligned down this axis.

The standard application of these techniques are also not able to resolve sub-nanometre structure and so the data from this study is combined with the evolved Feldman–Sereda model to give the “globule” model for C–S–H. Here, C–S–H can almost be thought of as crystalline cement powder particles dissolved in solution. Each (less than spherical) globule houses the aforementioned Ca + Si layers, separated by interlayer water which is said to enter irreversibly during adsorption. The water fills and expands the space between layers uniformly, while maintaining a stark boundary between Ca/Si layers and water itself, resulting in the “tubes” of water described above. The globules then flocculate, trapping exo-globular water in nanopores (Fig. 1).

Fig. 1 Possible morphologies of C–S–H globules with layered, expanded layered and non-layered structuring, reflected as a function of density. The layers comprise combined Ca and Si phases (green) while a disordered approach involves a spread of Si chains (blue) with Ca ions (yellow) interspersed between the chains. Bound water (blue pattern) is trapped in the interlayer space (order/left), else in rounded or oblong gel pore (disorder/right).

The density can be used as an indicator of order–disorder transitions in cement, with less ordered layers packing inefficiently leading to a less dense C–S–H globule and vice versa. Characterisation with ¹H-NMR on white (non-ferrous) cement pastes revealed a density of ∼1.9 g cm⁻³ for C–A–S–H gels.²⁴ This aligns with a multitechnique analysis of Al-uptake in C–A–S–H gels by Kapeluszna et al. which found that increasing Al-content widens the basal spacing and decreased crystallinity.¹¹

Pellenq et al. have evolved their models by factoring in defects to account for C–S–H's glassy character, stating that when averaged across all atoms, “… C–S–H is closer to a glass than a perfect crystal.”²⁵ However, their study did highlight some configurational ordering of silicate chains along with a “non-random arrangement of Ca atoms”. Neutron scattering structure factors, and by extension partial structure factors, were used to compare glass, C–S–H and crystal phases and revealed that the distribution of water and hydroxyl groups was amorphous. The implication here is that hydrated species lead to a much more disordered system. In any case this study utilised a broad approach with various topologies and morphologies (order vs. disorder) as starting points for dynamic simulations of the C–(A–)S–H phase.

1.2. Chains and pores

A common property of all C–S–H models are silicate chains of SiO₄⁻ tetrahedra that increase in length with age and temperature. These chains act as the foundation, with other constituent phases consolidating between and around them. The length of the Si chains follows a 3n − 1 rule (n = 1, 2, 3…) with a dreierkketten arrangement where the chain more-or-less repeats every three units. NMR data produced by Bauchy et al. confirmed that Q_n ratio of Si units in C–S–H averages to: Q₀ ≈ 10%, Q₁ ≈ 23.3%, Q₂ ≈ 66.7%, Q_3–4 < 0.1%.²⁶ The exact spatial distribution of these chains with respect to each other is unknown and requires extensive computational modelling coupled with experimental data to resolve.

Pores are another structural property to be incorporated into a model evolving towards an ever more realistic C–(A–)S–H structure, and present a challenge in terms of computational modelling and accounting for the interplay between cementitious (alumino-)silicate (cement) chains and H₂O, within the C–(A–)S–H matrix. Non-saturated hydrated cement pastes have been shown to reach total porosities of 36.6% and over time these pores generally shrink as hydration proceeds.²⁷ Kim et al. explored the influence of hydration degree by synthesising OPC pastes with water to cement ratios (w/c) ranging from 0.45 to 0.6 using cement that had already been cured for 91 days (achieved full hydration).²⁸ Using mercury intrusion porosimetry (MIP) they found that adding further water to the cement in fact increased both porosity and compressive strength.²⁹ While pores are normally associated with weakness in cement due to their active ‘allowance’ and promotion of chemical attack, corrosion, leaching and expansion, the role of pore networks are nuanced and yet unresolved at the local order. Pores are classified into two types: capillary pores >10 nm in diameter and Gel pores ∼0.5–10 nm.^30,31 Water is also classified separately with 1.3–1.5 H₂O/Si making up bound water which refers to water that is chemically bound through electrostatic interactions.³² Non-bound is merely physically trapped or adsorbed, often inside gel pores, with the exact amount left unclear. Gartner et al. attribute the presence of non-bound water to it being a by-product of Si oligomerisation.³³

Given that Bauchy et al. have already established a good degree of disorder inside cement, we actioned the modelling of C–S–H that fits the experimental data outlined above and allows us to investigate chain and pore length/volume and distribution with and without Al in ordered and disordered starting configurations.

2. Methods

2.1. Model construction

A multivariate analysis on various C–S–H models was conducted through a reductionist, ‘ground-up’ approach using the packmol programme package³⁴ to fill a unit cell/box for periodic boundary condition (PBC) calculations. Geometry optimisations using CP2K's QUICKSTEP with the BlYP method and triple-zeta TZVP-MOLOPT-SR-GT basis set, on all the starting configurations were performed to avoid any spurious effects; including chemically illogical structuring such as incomplete bonds, for example with ‘dangling’ or ‘hanging’ O-atoms which can occur in hydroxyl group containing systems (erroneous terminating with an –O group), else a final structure containing improbable multiplicity (e.g. 2 or unexpectedly higher) due to unpaired electrons.

Initial box dimensions were set to a = 26, b = 26, c = 26 Å for the disordered simulations and a = 35, b = 35, c = 35 Å for the ordered simulations. These orthogonal supercells allowed the volume of each cluster to vary freely as we tested the density limits of each cluster by subjecting them to 1000 bar over 125 ps in 0.25 fs increments via the NPT ensemble. At the end of the NPT simulations, the box dimensions had nearly halved in all directions culminating in cluster densities of ∼2.6 g cm⁻³. A series of 10 ps NVT runs (target T = 300 K, P = 1 atm) were completed for each model, to relax structures and remove spurious effects post-NPT, followed by a 25 ps NVE production run, in 1fs increments.

Fig. 2 outlines the variables that each model accounts for. This set of 246 atom models covered composition, configuration and conformation. At the composition level, the following three stoichiometries were explored: 0% Al (Si only), 8.3% Al – which represents the accurate Al/Si ratio of ∼1/12 and 25% Al – an extreme used to test an upper limit of Al's action in C–(A–)S–H. The percentages used here correspond to the proportion of tetrahedral units (24 possible) that are occupied by Al (see formulae, Fig. 2). From these three ratios, there are at least four possible conformations. The first is the ordered mode (Fig. 2, left) – whereby the layered-type structure seen in the literature are replicated and pores are modelled using interlayer water.

Fig. 2 Model construction outlining the three general variables. (1) Composition – stoichiometry varying with increasing Al amount. (2) Order/disorder – the topology of the model varying between the classical layered approach for C–S–H seen in literature against the “Disordered” model which represents a more homogenous blend of phases. Pore vs. hydrated pore (3) Al avoidance, resolving the influence of (non-)neighbouring Al atoms.

This classical interpretation is contrasted against three disordered modes (disordered, disordered pore, disordered hydrated pore). The silicate chains in each model are built upon the same Q_n-ratio emerging from the NMR data (Q₀ ≈ 10%, Q₁ ≈ 23.3%, Q₂ ≈ 66.7%).²³ Each cluster has 2 silicate monomers, 6 dimers and 2 pentamers, which Al are substituted into. The two modes of Al-substitution can be categorised in terms of Al avoidance: spatial and chain avoidance. These modes are highlighted in the third section of Fig. 2.

Models with 8.3% Al contain only two Al atoms, here Al-avoidance is defined in terms of real space in C–S–H, instead of being constrained to aluminosilicate chains. Clusters with spatial Al-avoidance have Al substituted into the central Q₂ position in both pentamer chains, while Al proximity sees these atoms adjacent to each other in the same chain. By contrast, 25% Al models investigate Al-avoidance in the manner originally implied by Loewenstein, in that Al–O–Si–O–Al is normally favoured over Al–O–Al–O–Si.^35,36 Models with Al avoidance are henceforth indexed with L_w (Loewenstenian) while models with Al in proximity to itself are left unlabelled. Another aspect to Al avoidance is that when two substituted Al atoms share an oxygen bridge, one must adopt a higher coordination number in accordance with the Pauling radius ratio rule.¹⁸

The mass (m) of each cluster is ∼6.58 × 10⁻²¹ g and with C–S–H having density (ρ) ≈ 2.6 g cm⁻³ when Ca/Si = 1.7, produces a target unit cell volume (V) of ∼2.531 × 10⁻²¹ cm⁻³ or 13.6 ų, where V = m/ρ. Volume may be quantified via simple multiplication of 3 chosen orthogonal axes else more rigorous means of evaluation (Section S1). This resulted in a pore that demonstrates changing structure, yet does not collapse in the NPT or NVE run trajectories. Models of the same nature were also constructed at quadruple the scale (984 atoms) and subjected to a similar simulation protocol of relatively short NPT and NVE runs. These larger models could be used to demonstrate whether the coalescing of layers in ordered C(A–)–S–H models was a conserved property, regardless of scale. The overall objectives of the AIMD simulations are summarised below:

• To confirm the adherence to the Al avoidance rule in C–(A)–S–H (configuration).

• To resolve what influence Al-ordering has on C–(A–)S–H chains, if any (configuration).

• To inform on possible mesoscopic ordering in C–(A–)S–H via highlighting if ordered conformations persist and comparing their final topologies with their disordered counterparts. (conformation).

• To investigate factors controlling mean chain length.

• To investigate the stability and shape of empty and hydrated pores over time (conformation).

• To produce a model for C–(A–)S–H that aligns with the SANS-resolved density and NMR-determined ratio of Al/Si chains.

2.2. Computational details

Towards investigating the time-dependent stabilities and propensity for dynamical interconversions of Al coordination (4 ↔ 5, 5 ↔ 6, 4 ↔ 6) at standard temperature and pressure (STP), selected geometry-optimised nanoclusters were subsequently investigated with ab initio molecular dynamics (AIMD) simulations, employing the CP2K code. O–Si–O, O–Al–O, Si–O–Si and Al–O–Al bond-angles were also characterised towards resolving angular populations and flexibility. CP2K uses a hybrid Gaussian plane-wave approach to implement the DFT method.^37,38 Relevant exchange–correlation was accounted for using the Becke–LeeYang–Par (BLYP), functional for its success in accurately predicting dynamic and electronic properties in aluminosilicate glasses;³⁹ whilst also presenting a good methodological match with the preliminary electronic structure determinations. Geodecker–Teter–Hutter pseudopotentials were used to mitigate expensive determinations of core configurations and all atomic species were represented using the double-zeta valence polarised MOLOPT basis sets.⁴⁰ Plane-wave kinetic energy cutoff was set to 250 Ry. Microcanonical (NVE) ensemble (involving constant number of particles N, volume V and energy E) were used to run the AIMD simulations over 100 ps trajectories, within which most of the local order phenomena occur as the gel anneals. These first principles (ab initio) dynamical simulations with NVT and NVE ensured to account for atomistic and quantum chemical aspects, within a dynamical context, and under standard conditions (T = 300 K, P = 1 atm). These AIMD models represent a compromise between accuracy and the limits of reasonable computing resources (i.e. high-end workstations, else non-specialist HPC allocations); accessible by the majority in the wider physical, materials, engineering and life sciences.

Towards further affirming the structures as being representatively stable, subsequent geometry-optimisation and frequency analyses confirmed the pores as being stables, with the structures residing at minima on their respective PEHSs using ab initio electronic structure computations, detailed in Section S2.

3. Results

3.1. Order disorder transitions

The transition shown in Fig. 4 represents a key finding, that under dynamic simulations, any layers present begin to coalesce and eventually all significant layering disappears into the fog of disorder. Originally this was thought to be a by-product of compression under NPT, hence to test the consistency of these layers at constant volume, the simulation was restarted with NVE only.

Fig. 4 shows how even under NVE, models of a water layer/hydrated pore is integrated into the bulk almost immediately and is no longer distinguishable from the bulk structure. This could also have been a limitation of the size and packing of the relatively small (246 atom) models, so the same test was performed using a larger model set. The bottom half of Fig. 4 shows one example of various attempts at packing these layers in accordance with their literature spacing, All attempts at layering failed in a similar fashion. Even under limited time scales, layers would begin to mix and result in a cluster that was indistinguishable from their disordered counterparts. Regardless, the ordered models that began as such in Fig. 2 were still subject to the full simulation protocol. Curiously, the final product of these layered simulations resembled a suggestion for the microstructure of C–S–H globules by Chiang et al.³⁹ While the internal ordering is not present the overall topology resembles the cylindrical model that was extrapolated from their measurements, with radius 1 nm. The timescale for these NPT runs equilibrate within the first ∼10th of the simulation (see Fig. S3), while results from production runs (NVE) are observed within the first half, and are maintained for the rest of the 25 ps. The remaining simulation time was prescribed as a ‘buffer’ for validation and any further changes in pore size, chain distribution etc. The structures were further confirmed as being stable, including pore robustness from ab initio electronic structure computations (Section S2).

3.2. Dynamic structure of pores

Generally empty pores made through the hybrid approach in Fig. 3 would lose half their volume and collapse well into the first third of the NVE run. Empty pores detected in experiment are a drying product leftover from non-bound water evaporating out of the C–S–H gel, as opposed to being an intrinsic structural feature of the gel. Otherwise, any gel pores in C–S–H, if they are to persist for any significant amount of time, must contain water. As shown in Fig. 4, pores manifesting as water layers quickly collapse from the beginning of the simulation. Water in these layers is not chemically bound to the rest of the cluster, hence dissolution of these pores upon drying is expected.

Fig. 3 Attempts to model disordered pores in AIMD simulations by 3 differing approaches ((1) retroactive, (2) proactive, (3) hybrid). The retroactive approach was successfully (Section 4), as was the hybrid one. The proactive approach 2 failed.

Fig. 4 Loss of short-range layering under NVE simulation (constant volume) in 246 atom clusters (top) and 984 atom clusters (bottom). Final structuring of these two systems is shown after 10 and 8 ps, respectively.

Instead, the simulations show how pores that persist over longer time periods are more likely spherical collections of non-bound water. While this water is not chemically bound and incorporated into the Si/Ca phases, the position of the pore is what renders escape difficult, given that the pore is embedded in the centre of the gel. Fig. 5 shows these hydrated pores maintaining their shape throughout the simulation. Pore size in dynamic runs also seemed to be dependent on Al/Si chain structure. Outside of flexing to accommodate larger pores, these chains could also form hydrogen bonds to leech water molecules out of the pore. By measuring the x, y and z dimensions of the voids left behind after each simulation, Table 1 lists the relative volume changes for each pore.

Fig. 5 Simulation of water filled pores. Atoms are greyed out (right side) to highlight the intrapore water molecules which help maintain pore size throughout the simulation.

Table 1 Pore volumes and relative changes in volume (nm³) and the percentage change, as determined from inter-atomic separation of the inside extremities of the pores. The Al-content of the models are distinguished as follows: 0% Al, 8.3% Al, 25% Al. Underlined values denote minimal volume changes. Model nomenclature and indexing follows: D = disorder, P = pore, H = hydrated, L_w = Lowenstein

Model	Vpore/nm^3		Δ
	Initial	Final	nm^3	%
0% Al (Al : Si = 0 : 12)
DP	0.292	0.119	−0.17	−59.3
DH	0.416	0.260	−0.16	−37.5
8.3% Al (Al : Si = 1 : 11)
DP Lw	0.331	0.152	−0.18	−54.2
DP	0.285	0.161	−0.12	−43.3
DH Lw	0.361	0.344	−0.02	−4̲.9̲
DH	0.285	0.244	−0.04	−14.2
25% Al (Al : Si = 3 : 9)
DP Lw	0.244	0.125	−0.12	−48.7
DP	0.252	0.127	−0.13	−49.8
DH Lw	0.363	0.397	+0.03	+9̲.4̲
DH	0.410	0.243	−0.17	−40.7

---

The primary factor affecting pore size is reaffirmed as water content. Secondary to this is the role aluminium plays. Al substitution at Q₂ creates a fulcrum for the rest of the Si/Al chain to flex around. At 8.3% Al, the Al avoidance model (L_w) has both Al atoms at the Q₂ position of each pentamer, meaning both chains have a hinge at the midpoint. As a result, the pore in this model only loses 4.85% of its volume, while its Al proximity counter part loses 14.22%. These figures are attributed to the shape of the chains and their proximity to the pore.

One of the two pentamer chains is slightly closer to the pore, serving as the boundary between the pore and the rest of the cluster. In both Al 8.3% cases, the pentamer that runs closest to the pores contains Q₂ Al, and thus flexes around the void while remaining fairly parallel to the pore surface (i.e. not intruding on the pore space). By comparison, the same hydrated pore model with 0% Al sees a loss of 37.52% of its volume, as this pentamer breaks. This cleavage of the pentamer leaves a smaller Si dimer which is free to move into the pore and form hydrgoen bonds with the water there, reducing pore volume. Al avoidance is further shown to facilitate larger pores at high Al (25% Al models), with the contrast between Al proximity yielding a 40.67% reduction versus and increase of 9.36% in volume for Al avoidance. The slight increase in pore size here can be attributed to the fact that the original pore volume is comparatively low (0.363 vs. 0.410 nm³) due to compression from the NPT phase. Once the structure is relaxed the pore size balloons with the help of flexible Al/Si chains slightly increasing its size to 0.4 nm³. The second column in Table 1 corresponds to the pore sizes at the beginning of the NVE run, whereupon the clusters had already been compressed and achieved the correct density. Hydrated Pores contained 4 water molecules which resulted in larger pores that generally kept more of their original volume. Empty pores, regardless of Al substitution or ordering, will lose at least half their volume. Having established the density and presence of hydrated pores, what remains is to confirm which of the models, if any, still conform to the Q_n ratio for C–S–H chains given by NMR (Table 2).

Table 2 Summary of Al coordination and Q-number for differing Al-content (0, 8.3, 25%), respectively. All models retain NMR-determined values of Q-number ratioing of Al, as follows: Q(0 : 1 : 2) = 8.3 : 66.67 : 25.0. Again, with D = disorder, P = pore, H = hydrated, L_w = Lowenstein

Model	Al-coord’n	Al Q# ‘count’
0% Al (Al : Si = 0 : 12)
D	N/A	N/A
DP	N/A	N/A
8.3% Al (Al : Si = 1 : 11)
DLw	1(Al4) + 1(Al5)	2(Q2)
DPLw	1(Al4) + 1(Al5)	2(Q2)
DH	2(Al4)	2(Q2)
25% Al (Al : Si = 3 : 9)
DPLw	2(Al4) + 4(Al5)	3(Q1) + 3(Q2)
DHLw	3(Al4) + 3(Al5)	3(Q1) + 3(Q2)
DH	5(Al4) + 1(Al5)	6(Q2)

---

3.3. Al-avoidance influence on chain configuration

There is no direct trend between Al avoidance and energy, rather the issue affecting C–S–H energy is the sum of many individual parts. The order disorder transitions outlined above are one such example where ordered models lose their layering and become difficult to distinguish from the disordered models. Anhydrous pores also collapse leaving them indistinguishable from the non-porous disordered models. Given that the morphology and general structure of these models (ordered, disordered and disordered + pore) should be very similar, we expect their energies to be similar. Instead, the energies vary by hundreds of kJ mol⁻¹. Even given the unreliability of energy values taken from AIMD over free energy, these differences are too high to be in the order of coordination shifts or even the breaking/making of a few bonds. We can instead explain the trends in energy by looking specifically at the chains that arise because of Al avoidance, coordination etc.

Bond lengths in all models barely deviate over time. Fig. 6 shows a radial pair distribution function comparing all bond lengths for hydrated pore models with Al avoidance and with no Al. The three distributions layered over each other reflect the nature of the remaining models, in that average bond distances for Al–O, Si–O and Ca–O pairs were consistent and independent of Al avoidance. Ca–O pairs have been shown to float around the 2–3 Å range, given that their interactions are largely ionic in nature. The Ca–O lengths are dependent on oxygen availability, since O²⁻ ions are the only source relative negative charge and are craved by every other positively charged ion (Si⁴⁺, Al³⁺, H⁺) among which Ca²⁺ has the lowest charge density and thus has the ‘loosest grip’. Radial distribution functions including Al/Si chains normally show a shoulder around 1.8 Å, where Al–O bonds stretch to accommodate the shorter, more rigid, Si–O bond lengths that neighbour them. The real show of flexibility is not found in bond length but rather angular flexibility and its influence on chain length.

Fig. 6 Radial pair distribution functions for 33%, 9% and 0% Al hydrated pore models with Al avoidance (where applicable).

Fig. 7 contains the angular distribution functions for Al–O–Al, O–Al–O, O–Si–O, Si–O–Si as well as Al–O–Si hinges for every model. An easy means of interpretation is using O–Al–O and O–Si–O as indicators of local coordination, while Al–O–Al and Si–O–Si are the hinges that the control the shape and flexibility of the chain at large. We can see how the chain hinges are dependent on local coordination by relating the rigidity of Si, indicated by a sharp 110° peak (Fig. 7, top right), to the same behaviour in Si–O–Si hinges (Fig. 7, top left). As we move down to Al–O–Si hinges, the effect of Al atoms both in small and large quantities is apparent by contrast to 0% Al. These new Al–O–Si hinges show a more amorphous spread of angular allowance in cases with 25% Al, where multiple Al atoms are substituted into the same chain. This flexibility is derived from O–Al–O angles that dictate Al coordination. Preferences for Al-4 and Al-5 are present, with the main peaks in the range of 110° (Al-4) and 90° (Al-5). Where Al units neighbour each other, and thus are viable as models for C–S–H short range order. The model ‘disorder Al hydrated pore L_w’ represents a success at modelling Al substituted C–S–H with a hydrated pore, while meeting the NMR and SANS data requirements. Fig. 8 shows every 246 atom C–S–H cluster that was simulated. Each cluster began the simulation with two pentamer chains which evolve depending on Al coordination. These chains are highlighted, and their lengths given in blue. Interestingly, the remaining Si dimers and monomers would not oligomerise independent of Al, further verifying trends seen Al assisted oligomerisation.¹¹ Chains with no Al can be seen cleaving in the top row of Fig. 8. Al 0% DH shows a hydrated pore pushing against an Si pentamer, leading to cleavage. As far as Al-position in the structuring, Al–O–Al moieties have been zeolites, yet the more relevant trait of Al-avoidance seen here is its effect on chain growth.⁴¹ Al makes for longer chains with higher degrees of branching, whilst suggested as a cause for kinetic energy destabilisation Al proximity causes the chains to coil and fold in on in themselves. Comparing 25% Al order L_w to non-L_w gives mean chain lengths of 5.3 and 5, respectively. Al avoidance leads to longer chain lengths due to the Al units being spread out to the ends of the chain. These units serve as better chain mergers.¹¹ In a chain of length n = 5, having chain mergers at n = 1, n = 3 and n = 5 means Si monomers and dimers can merge at both termini as well as above the centre (Q₃ branching Al). An example of this is the comparison between 25% Al L_w and non-L_w, where one pentamer merges with part of the other pentamer, creating a branched heptamer (7 + 4). The non-L_w (Al proximity) counterpart exhibits similar behaviour but instead coils in on itself due to the flexibility of the Al hinges that are clustered at the centre of the chain, taking on a ring like conformation. This is a direct example of how Al ordering can potentially shift otherwise symmetrical silicate chains, to being branching chains with raised disorder, again echoing findings in zeolite systems.⁴² As for the implications for cement, our model aligns with the local properties of C–(A)–S–H which is shown to have higher chain branching (Q₃ Al) than tobermorite, the former being mechanically tougher, the latter being stronger and more brittle. These phases differ very little in composition, yet their configuration as a function of chain connectivity is what leads to the difference in mechanical properties.¹¹ The argument that the disorder, branching and thus variance in chain order is reserved only for Al substituted systems, is muddied when in practical C–S–H examples. The supposed non-Al substituted C–S–H has been shown to uptake a minimal amount of Al and can achieve an Al/(Al + Si) ratio of 0.03. Instead, there seems to be some level of short-range disorder in 0% Al C–S–H, which is further enhanced on addition of Al. As with previous studies we present a theoretical scaffold that potentially connects the nanostructure to the mechanical behaviour.⁴³

Fig. 7 AIMD simulation results on the time-dependent structuring and dynamical angular flexibility resolved from NVE-determined angular distribution functions (ADFs) for O–Si–O and Al–O–Al hinges. Additional ADFs for Si–O–Si, Si–O–Al and Al–O–Al are disseminated in Fig. S2. D = disorder, O = order, S = no Al, P = pore, hydrated pore L = Al avoidance N = no Al avoidance. Top right shows example illustrations of DHL_w, DH, O and OL_w.

Fig. 8 Al/Si chains highlighted (transparent ‘goldish’ oval overlays) in each of the C–S–H clusters, sorted by conformation (columns) and Al amount/avoidance (rows). Blue values used to highlight the new lengths of the two chains that began as pentamers. E.g., 5 indicates a chain length of n = 5, branching is shown using +, where 5 + 1 is n = 5 chain with one branching unit.

4. Conclusions

The dynamic C–S–H simulations were successful in suggesting a new model for C–S–H, with and without Al, that conforms to the experimental data informing chain length, connectivity and overall density. The success in modelling a hydrated gel pore of 0.4 nm³, further highlights the importance of Al in sustaining pore volume. The Al avoidance rule is adhered to in terms of Al coordination and leads to longer chains.

There is no direct trend between Al avoidance and energy, as this property is buried amongst all the other variables that contribute to the shape and length of Al/Si chains, which in turn influence the morphology of the overall cluster. Ordered conformations for C–S–H default to their disordered counterparts, as layering does not persist in an unconstrained simulation (supercell). While these order disorder transitions establish a general idea of the mesoscopic conformation, the conformation of individual Al/Si chains must be indexed to find favourable modes. In theory, cataloguing these chains would lead to methods of controlling their shape and orientation, which would affect the integrity of either the C–S–H or C–A–S–H structures.

The AIMD method is suitable for these systems over the trajectories of the NVE and NPT runs, due to its accuracy of evaluating the nature of the chemical bonding (including making/making), in particular the dynamic bond-angle flexibility. This allowed for the models to be characterised without relying on forcefields (reactive or otherwise) built on crystal structures and simplified C–S–H models, which may not have revealed the order–disorder transitions in these nano-particle cement models. The trends uncovered in this work provide the groundwork for building larger C–(A–)S–H clusters, into mesoscale sizes (100s of nm and larger. As envisioned in Fig. 4, future studies could evolve such mesoscale models with their short range order (SRO) built-up from (semi-)amorphous Al-substituted cements, where we anticipate similar structural and functional trends remain consistent.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Additional methodological and data is available, including for the computational models presented in the work. See DOI: https://doi.org/10.1039/d5ma00244c.

Acknowledgements

All authors thank MEH Srl (Rome, Italy), MEH-Canada (Toronto & Calgary, Canada), and GIOCOMMS (Toronto, Canada) for providing support to this work and the wider scope of projects ongoing, as well as Modern Age Plastics Inc. (MAPI), Toronto, Canada for their continued support (2020-present) for fabrication of equipment to advance empirical industrial tests on related material systems. The European Commission is thanked for co-funding the BUCK$$$ project under the joint call CETPartnership 2022 (Grant Agreement No. 101069750). Further, STFC and RAL-ISIS are acknowledged for funding related activities contributing to this work (RB2310655, RB1710447 and RB1710444). MS, DDT and GAC acknowledge Leverhulme Trust (RPG-2023-239) for resources supporting projects advancing materials characterisations by both computational and empirical means. MS acknowledges PRINN (Italy) for computational support. DHS acknowledges the support of dgNetrix, Canada (award 202408). RYY acknowledges the support of the Provost's office at the University of British Columbia (UBC), Vancouver, Canada as well as NSERC, Canada (RPGIN04598). KVT thanks McMaster University for continued personal support to advance project work between Canada and Italy. Finally, all authors are grateful for the efforts by the above acknowledged to continue forming the growing synergy between academia, science and industry to generate impact in Canada, the UK and Italy.

References

1. D. Hodgson, IEA Cement – Analysis, 2020.
2. T. W. Bremmer, et al., Environmental aspects of concrete: problems and solutions. In All-Russian conference on concrete and reinforced concrete, 2001.
3. P. Brune, R. Perucchio, A. R. Ingraffea and M. D. Jackson, The toughness of imperial roman concrete. In Proceedings of the 7th International Conference on Fracture Mechanics of Concrete and Concrete Structures, Jeju Island, Korea, 2010, vol. 52, pp. 23–28.
4. K. Shen, et al., Revisiting Ancient Roman Cement: The Environmental-Friendly Cementitious Material Using Calcium Hydroxide-Sodium Sulfate-Calcined Clay, ACS Sustainable Chem. Eng., 2023, 11(13), 5164–5174.
5. N. J. Delatte, Lessons from Roman cement and concrete, J. Prof. Issues Eng. Edu. Pract., 2001, 127(3), 109–115.
6. F. Massazza, Pozzolanic cements, Cem. Concr. Compos., 1993, 15(4), 185–214.
7. H. Liu, et al., From molten calcium aluminates through phase transitions to cement phases, Adv. Sci., 2020, 7(2), 192209.
8. R. A. Davies, et al., Geometric, electronic and elastic properties of dental silver amalgam γ-(Ag3Sn), γ1-(Ag2Hg3), γ2-(Sn8Hg) phases, comparison of experiment and theory, Intermetallics, 2010, 18(5), 756–760.
9. V. H. Dodson, Pozzolans and the pozzolanic reaction, Concrete admixtures, Springer US, Boston, MA, 1990, pp. 159–201.
10. M. D. Jackson, et al., Material and elastic properties of Al-tobermorite in ancient Roman seawater concrete, J. Am. Ceram. Soc., 2013, 96.8, 2598–2606.
11. M. D. Jackson, et al., Unlocking the secrets of Al-tobermorite in Roman seawater concrete, Am. Mineral., 2012, 98(10), 1669–1687.
12. E. Kapeluszna, et al., Incorporation of Al in CASH gels with various Ca/Si and Al/Si ratio: Microstructural and structural characteristics with DTA/TG, XRD, FTIR and TEM analysis, Constr. Build. Mater., 2017, 155, 643–653.
13. M. S. Salha, R. Y. Yada, D. H. Farrar, G. A. Chass, K. V. Tian and E. Bodo, Aluminium catalysed oligomerisation in cement-forming silicate systems, Phys. Chem. Chem. Phys., 2023, 25(1), 455–461.
14. K. V. Tian, M. Z. Mahmoud, P. Cozza, S. Licoccia, D. Fang, D. Di Tommaso, G. A. Chass and G. N. Greaves, Periodic vs. molecular cluster approaches to resolving glass structure and properties: Anorthite a case study, J. Non-Cryst. Solids, 2016, 138–145.
15. E. L’Hôpital, et al., Incorporation of aluminium in calcium-silicate-hydrates, Cem. Concr. Res., 2015, 75, 91–103.
16. I. Lognot, et al., NMR and infrared spectroscopies of CSH and Al-substituted CSH synthesised in alkaline solutions, Nuclear magnetic resonance spectroscopy of cement-based materials, Springer Berlin Heidelberg, 1998, pp. 189–196.
17. X. Pardal, Experimental study of Si–Al substitution in calcium-silicate-hydrate (CSH) prepared under equilibrium conditions, Cem. Concr. Res., 2009, 39(8), 637–643.
18. W. Loewenstein, The distribution of aluminum in the tetrahedra of silicates and aluminates, Am. Mineral., 1954, 39(1–2), 92–96.
19. R. E. Fletcher, S. Ling and B. Slater, Violations of Löwenstein's rule in zeolites, Chem. Sci., 2017, 8(11), 7483–7491.
20. L. Pegado, C. Labbez and S. V. Churakov, Mechanism of aluminium incorporation into C-S-H from ab initio calculations, J. Mater. Chem. A, 2014, 2(10), 3477–3483.
21. R. F. Feldman and P. J. Sereda, A new model for hydrated Portland cement and its practical implications, Eng. J., 1970, 53(8–9), 53–59.
22. T. C. Powers, Structure and physical properties of hardened Portland cement paste, J. Am. Ceram. Soc., 1958, 41(1), 1–6.
23. A. J. Allen, J. J. Thomas and H. M. Jennings, Composition and density of nanoscale calcium–silicate–hydrate in cement, Nat. Mater., 2007, 6(4), 311–316.
24. A. C. A. Muller, K. L. Scrivener, J. Skibsted, A. M. Gajewicz and P. J. McDonald, Influence of silica fume on the microstructure of cement pastes: New insights from 1H NMR relaxometry, Cem. Concr. Res., 2015, 74, 116–125.
25. M. Bauchy, M. A. Qomi, F. J. Ulm and R. M. Pellenq, Order and disorder in calcium–silicate–hydrate, J. Chem. Phys., 2014, 140(21), 214503.
26. R. J. M. Pellenq, A. Kushima, R. Shahsavari, K. J. Van Vliet, M. J. Buehler, S. Yip and F. J. Ulm, A realistic molecular model of cement hydrates, Proc. Natl. Acad. Sci. U. S. A., 2009, 106(38), 16102–16107.
27. E. Gartner, I. Maruyama and J. Chen, A new model for the C-S-H phase formed during the hydration of Portland cements, Cem. Concr. Res., 2017, 97, 95–106.
28. F. Avet, E. Boehm-courjault and K. Scrivener, Cement and Concrete Research Investigation of C-A-S-H composition, morphology and density in Limestone Calcined Clay Cement (LC3), Cem. Concr. Res., 2019, 115, 70–79.
29. Y. Kim, K. Lee, J. Bang and S. Kwon, Effect of W/C Ratio on Durability and Porosity in Cement Mortar with Constant Cement Amount, Adv. Mater. Sci. Eng., 2014 DOI:10.1155/2014/273460.
30. L. Li, H. Zhang, X. Guo, X. Zhou, L. Lu and M. Chen, Pore structure evolution and strength development of hardened cement paste with super low water-to-cement ratios, Constr. Build. Mater., 2019, 227, 117108.
31. H. M. Jennings, J. W. Bullard, J. J. Thomas, J. E. Andrade, J. J. Chen and G. W. Scherer, Characterization and Modeling of Pores and Surfaces in Cement Paste, J. Adv. Concr. Technol., 2008, 6(1), 5–29.
32. R. R. Lloyd, J. L. Provis, K. J. Smeaton and J. S. J. Van Deventer, Microporous and Mesoporous Materials Spatial distribution of pores in fly ash-based inorganic polymer gels visualised by Wood’ s metal intrusion, Microporous Mesoporous Mater., 2009, 126(1–2), 32–39.
33. H. M. Jennings, A model for the microstructure of calcium silicate hydrate in cement paste, Cem. Concr. Res., 2000, 30, 101–116.
34. A. Tossell, A theoretical study of the molecular basis of the al avoidance rule and of the spectral characteristics of al-o-al linkages, Am. Mineral., 1993, 78(9–10), 911–920.
35. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, Packmol: A package for building initial configurations for molecular dynamics simulations, J. Comput. Chem., 2009, 30(13), 2157–2164.
36. S. K. Lee and J. F. Stebbins, The degree of aluminium avoidance in aluminosilicate glasses, Am. Mineral., 1999, 84(5–6), 937–945.
37. M. Filatov, W. Thiel, B. M. Filatov, W. Thiel and C. Zu, Density functional, 2010, vol. 8976.
38. J. VandeVondele and J. Hutter, Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases, J. Chem. Phys., 2007, 127(11), 114105.
39. W. L. Li, et al., Optimized pseudopotentials and basis sets for semiempirical Density Functional Theory for electrocatalysis applications, J. Phys. Chem. Lett., 2021, 12.42, 10304–10309.
40. W. S. Chiang, E. Fratini, P. Baglioni, D. Liu and S. H. Chen, Microstructure determination of calcium-silicate-hydrate globules by small-angle neutron scattering, J. Phys. Chem. C, 2012, 116(8), 5055–5061.
41. A. V. Larin, The Loewenstein rule: the increase in electron kinetic energy as the reason for instability of Al–O–Al linkage in aluminosilicate zeolites, Phys. Chem. Miner., 2013, 40, 771–780.
42. T. Takaishi, Ordered distribution of Al atoms in the framework of analcimes. Journal of the Chemical Society, Faraday Transactions, 1998, 94.10, 1507–1518.
43. N. Li, D. Zhou, J. Wan, Z. Zhang, F. Jia, J. Ye, X. Zhi and W. Chen, Atomic-scale insights into mechanical properties of calcium silicate hydrates: role of hydrogen bond networks and bond order distributions, RSC Adv., 2025, 15(13), 10365–10377.

---