Structure, hydration, and chloride ingress in C-S-H: Insight from DFT calculations

The structure of Calcium-Silicate-Hydrate (C-S-H) and the effect of variations in its water content have been investigated using density functional theory (DFT) calculations. Trends for calculated densities as a function of hydration are in good agreement with experimental values, and in line with what is found using molecular mechanics in the literature. While we observe very little variation in SieO and CaeO bond lengths between different structures, structural diversity is otherwise great, in accordance with experimental observations, as we see no obvious correlation between structural features and material system stability. A mapping of energetics of hydroxyl substitution with chloride reveals, unsurprisingly, that chloride preferentially coordinates to calcium. More specifically, it was found that the most stable sites for chlorine substitution involves at least two adjacent calcium atoms. Computed chloride substitution energies indicate that the C-S-H phase may bind chloride from aqueous solution, potentially influencing chloride diffusion in concrete.


Introduction
Portland cement has been in use for almost 200 years and constitutes the great majority of the global annual cement production, which stands at some 4 billion tons [1]. The main component of hydrated cement paste is calcium silicate hydrate (C-S-H). Thus, its properties are important in determining concrete behaviour. More than 40 C-S-H phases have presently been characterized crystallographically [2]. The C-S-H phase encountered in concrete is a layered, mostly amorphous but also nanocrystalline material including a complex network of interconnected, water filled pores. The pores contribute to the elemental transport properties of the concrete mix and are categorized according to size: Macro pores (> 10 μm), capillary pores (≈10 nm-10 μm), and gel pores (< 10 nm) [3]. C-S-H is hence a material with structure on multiple scales. In the present contribution, we consider the atomic scale.
The vast present répertoire of solid-state characterisation techniques can give very detailed information on simpler materials, but in the case of cement, its semi-crystallinity, variable stoichiometry, and structural diversity lead to a need for further information to determine its structure on the atomic scale [4,5]. C-S-H is found to have a Ca/Si molar ratio that varies between 1.2 and 2.1, even in the same sample at a spatial resolution of 100 nm [6]. Richardson [7] presented a perfectly ordered model of C-S-H bulk structure fully based on X-ray diffraction experiments. Chen et al. [8], studied C-S-H structure using MASNMR experiments, affording results that are complementary to those found with diffraction techniques. Among other things they report mean silicate chain lengths as a function of Ca/Si ratio and note that this dependence varies with synthesis method of the tested material. Mohamed et al. [5] summarize that the structure of C-S-H may be regarded as defective, nanocrystalline tobermorite with missing bridging silicate chains. The latter lead to the modest average length of silicate chains observed experimentally, and to deprotonated silanol groups. Subsequently, this charge is compensated by additional calcium ions and Ca-OH groups.
Most concrete structures are reinforced with steel rebar for enhanced structural integrity. The carbon steel used is stable in the alkaline environment of a pure hydrated C-S-H but will corrode if the chloride concentration in the solution in equilibrium with the C-S-H exceeds a certain threshold and induces passive film breakdown. When chloride reaches the steel rebar, irreversible corrosion damage begins [3,9] followed by cracking, spalling and failure. In a saline environment, corrosion of the steel rebar due to chloride ingress is one of the major causes of degradation of reinforced concrete. The thermodynamic driving force for the chloride diffusion is concentration gradients [3]. In saline environments, the chloride concentration is higher at the outer surface of the concrete structure, and chloride ions will move from the surface into the water filled pore structure of the C-S-H. Thus, the kinetics of macroscopic chloride diffusion in C-S-H is assumed to depend largely on the pore structure [10,11]. The degree of hydration [12,13] may also affect chloride transport, which is assumed to occur mainly in the larger pores of C-S-H, i.e. macro pores and capillary pores [3]. Chloride in C-S-H are generally divided into two different types; "free" and "bound" chloride. The "free" chloride ions are found in the aqueous phase filling the pores, and their relative mobility makes them relevant for corrosion. The latter type is chemically bound chloride species. As the latter can be released over time, both types may contribute to corrosion damage [14]. The focus of this work is on bound chloride, and on processes on the electronic and atomistic scales.
Atom-scale modelling has been applied to cement in recent years [5,[15][16][17][18][19][20][21][22][23][24][25], mostly based on empirical force fields. The database cemff [19] contains a thorough comparison and validation of interatomic potentials for cementitious systems, allowing researchers to choose the most appropriate potential for their application. Interatomic potentials are generally valid for the problem for which they are developed; their transferability to related problems may be questionable. ClayFF [26] and its re-parameterization CSH-FF [27] are the most commonly used potentials for tobermorite and C-S-H gels. They perform well for mechanical property calculations, but are not validated for surfaces or interfacial properties. They do not describe chemical bond breaking and are therefore unable to capture, for instance, the splitting of water and formation of Ca-OH in C-S-H. However, ReaxFF potentials [28], while being more computationally expensive, contain terms designed for the description of formation and dissociation of bonds.
Churakov has investigated aspects of the structure [22] and water sorption [23] of 11 Å tobermorite using DFT-based molecular dynamics simulations. Furthermore, Churakov et al. [24] studied the intrinsic acidity of the C-S-H surface towards the water phase in a pore using DFT-based molecular dynamics and thermodynamic integration. A follow up study [25] describes aluminium/silicon exchange in the pore surface through a similar approach. In their molecular dynamics study, Pellenq and co-workers [15,16] used C-S-H models based on the structure of tobermorite in conjunction with ReaxFF potentials, addressing the effect of Ca/Si ratio on the structure and mechanical properties of C-S-H. Kovačević et al. [17,18] applied reactive force fields and semiempirical quantum chemical methods to study C-S-H. Three different C-S-H model structures were optimized, and molecular dynamics simulations were carried out to investigate volume changes as function of water content. In the works by Pellenq et al. and Kovačević et al., unit cells of sizes varying from roughly 20 × 20 × 20 Å 3 to 30 × 30 × 36 Å 3 were used. Mohamed et al. [5] performed both density functional theory (DFT) and forcefield-based simulations to establish an atomic structure of C-S-H, using smaller unit cells of about 7 × 7 × 28 Å 3 . They also proposed a building block scheme for creating C-S-H models for a range of Ca/Si ratios in order to accelerate identification of atomistic structures for C-S-H of different stoichiometries and to allow fast calculations of key characteristics.
Atomistic modelling has been used to study chloride diffusion in C-S-H [21,[29][30][31]. Kalinichev et al. [31] used MM/MD to study the interaction of chloride with surfaces relevant to hydrated cement in one of the earliest applications of their developed ClayFF force field. They did not observe any surface-bound chloride ions, in contrast to the other phases under investigation (portlandite, Friedel's salt, and ettringite). Hou et al. [21] studied the transport of chloride ions in capillary pores of C-S-H. They used the same methods and found that the tobermorite surface along pores repulses chloride ions and absorbs calcium ions. Pan and co-workers found a quite strong adsorption of chloride to the surface of portlandite, whereas no detectable binding were found for the surface of tobermorite or jennite [29]. Later, Zhou et al. [30,32] used a combination of molecular mechanics/molecular dynamics (MM/ MD) with the ClayFF potential and nuclear magnetic resonance (NMR) experiments, and found that binding of chloride ions to the tobermorite surface was greatly enhanced with increasing number of calcium ions in the vicinity. This effect was much stronger for calcium present in the C-S-H structure as opposed to free calcium present as ions in the water/ torbermorite interphase. However, they did not observe chloride ion diffusion into the bulk C-S-H-structure.
No clear consensus can be found in the literature, neither on the atomic structure of C-S-H nor on the mechanism of chloride ingress. Simulations based on first principles calculations may provide valuable insight exceeding what is available solely from the variety of experimental results and classical atomistic simulations. Earlier computational studies of chloride in C-S-H have been based on non-reactive classical forcefields, limiting the possibility to include the interaction with surfaces. Furthermore, mechanistic insight into the possible interaction of chloride within the gel pores is, to our knowledge, missing.
In the present work, the structure of C-S-H, its volume as a function of water content, and the interaction of chloride with C-S-H are investigated through density functional theory (DFT) calculations. A quantum mechanical approach has been chosen, despite being more computationally demanding, to avoid the approximation of interatomic potentials used in most literature to date. In the present contribution, results on the structure and volume correlations depending on water content in C-S-H models using DFT calculations are presented. Sampling of structures prior to DFT calculations were performed with classical MM/MD using the reactive ReaxFF force field. Further, chloride absorption properties within these C-S-H models are discussed, based solely on relative thermodynamic stability of different chloride sites within the material.

Models for the C-S-H phase
To construct the models used in the simulations the procedure of Kovačević et al. [18] was followed. The experimental structure determined for the 11 Å tobermorite by Merlino et al. [33] was used as starting point with lattice parameters of a = 6.732, b = 7.369, c = 22.680, α = β = 90.0° and γ = 123.18°, corresponding to a 2 × 1 × 1 supercell of the tobermorite structure, was employed. The models consist of two intralayers represented by calcium and silicate chains, and two interlayers with bridging silicon atoms and calcium and water. As illustrated in Fig. 1, there are two types of silicate tetrahedra in the intralayers; bridging tetrahedra (BeSi) and paired tetrahedra (PeSi).
A C/S (calcium/silicon) ratio of 1.647 is achieved by first increasing the number of interlayer calcium atoms followed by randomly removing silicon and calcium atoms. Two different C-S-H models were created, model A and model B, observing the following procedure (see • Model B: randomly remove most bridging silica tetrahedra (denoted BeSi in Fig. 1) and one paired silica tetrahedron (denoted PeSi in Fig. 1) 4. Add 2 protons to replace each removed silicon atom, forming water or hydroxyl groups. 5. Leave some random bridging silica tetrahedrons (the remaining bridging silica tetrahedra are not directly connected to another bridging silica tetrahedron) 6. Randomly remove calcium atoms in interlayer to reach the desired C/S ratio. For Model B one intralayer calcium atom close to the vacancy from the removed intralayer silicon is removed. 7. Add water randomly in interlayers to desired hydrate to silica (water/Si) ratio. 8. Perform MM/MD simulations as described below. 9. Optimize resulting structures using periodic DFT as described below.
Starting structures for MM/MD runs were prepared in the following way: A set of four different C-S-H structures were created for each of the models A and B, and for each of these 8 structures representing three levels of water content were investigated (water/Si = [1.5, 1.8, 2.1]), giving a total of 24 different initial structures. For each of these structures, three full computational procedures (as defined above) were performed differing only in random generator seeds for initial velocities in the MD calculations. The total number of parallel runs were therefore 72, as shown schematically in Fig. 2. The chemical formula at the lowest water level (water/Si ratio of 1.5) for both models A and B is Ca 28 Si 17 O 88 H 52 . The total number of water molecules (before optimization) within the simulation cell is 26 at this level, and 5 water molecules are added to the unit cell for each increase in water/Si giving a final chemical composition of Ca 28

Computational details
To optimize the generated structures and achieve a good starting point for the DFT calculations, molecular dynamics simulations where performed prior to the DFT relaxations for identifying the most stable structures. The molecular dynamics simulations where performed using the ReaxFF potential [28] parameterized for Ca-Si-O-H [19] as implemented in the LAMMPS simulation software package [34]. The DFT calculations at the PBE-GGA level [35] were performed using the Vienna ab-initio simulation package (VASP) [36,37]. Based on initial convergence tests a plane-wave energy cut-off of 700 eV was chosen, and a Γ-centered k-point grid of 2 × 4 × 1 was used to sample the Brillouin zone. Van der Waals corrections, using the zero damping DFT-  D3 method of Grimme, [38] were included. The computed results were analysed using ASE [39] software, while Ovito [40] was employed for visualization.
Models for the C-S-H phase were optimized as follows: • A series of simulated annealing MM/MD (ReaxFF) calculations, starting from different random seeds, were performed for each of the 72 starting structures according to the following protocol: 1. Energy minimization of the original structure and volume to a force tolerance of 2.5 · 10 −4 eV/Å. 2. NPT molecular dynamics at 400 K for 1000 ps using a timestep of 0.5 fs. The temperature was chosen after empirically determining that any higher temperature led to disruption of the C-S-H structure. 3. Molecular dynamics cooling down the system from 400 K to 0.1 K over 250 ps at a timestep of 0.5 fs. 4. Energy minimization of the final structure and volume to a force tolerance of 2.5 · 10 −4 eV/Å.
• The structures and volumes were then optimized to a force tolerance of 0.02 eV/Å with DFT as follows: 5. Ionic relaxation. 6. Volume relaxation to an electronic energy tolerance of 10 −5 eV/ Å. 7. Ionic relaxation to a force tolerance of 0.1 eV/Å. 8. Volume optimization [41]. 9. Ionic relaxation to a force tolerance of 0.02 eV/Å. 10. If convergence, i.e. forces of 0.02 eV/Å, is not achieved after step 9, the steps 8 and 9 of the procedure are repeated until convergence.
A referee commented on the issue of pH in our model. As the latter represents bulk solid, the concept of pH is not well defined here. On the other hand, the model allows flexibility in the description of oxygen as water, hydroxyl, or oxide as it is based on first principles. To evaluate the relative energetic stability of the various C-S-H structures at the various water levels, the following metric was chosen: where E CSH+(n+m)H 2 O is the total energy of the C-S-H system with (n + m) water molecules, E CSH+mH 2 O is the total energy of C-S-H system with the least amount of water where we used the most stable of the CSH + mH 2 [42]) is the experimentally determined heat of evaporation of water. We chose this approach since liquid water is the practically relevant reference for applications of concrete. The choice of reference is, of course, irrelevant for comparison of systems of identical stoichiometry; however, comparing systems with differing water content, the liquid water reference leads to much smaller energy differences than what would result from referencing gaseous water. A more self-consistent procedure would include calculating the energy of vaporization from first principles, but the present state of the art does not, to our knowledge, allow this with any degree of accuracy that would be useful in the present context. Chloride sorption was investigated for all of the optimized C-S-H models with water/Si ratio of 1.5 and 1.8. To maintain charge neutrality, chloride was introduced by substitution of a hydroxyl group by a chloride ion. For each of the 48 model systems, each hydroxyl group in turn was changed to chloride by positioning the chloride atom where the oxygen used to be, and the structure was then optimized. Another set of structures was generated by placing chloride at a sterically accessible site within a distance between 3 and 7 Å from the originally removed OH group, followed by optimization. Some resulting duplicates were discarded. The remaining C-S-H structures containing chloride were optimized with no restrictions placed on the ion positions. These structures, each with one chloride ion within the computational cell, were analysed further.

Structure
The relative stability ΔE of the different C-S-H structures at the lowest water content investigated (water/Si ratio of 1.5) are listed in Table 1 with corresponding structural properties. The choice of water/ Si ratio can be debated, as higher ratios are generally observed experimentally. The choice of hydration level for the structural investigation was made considering the results reported by Allen et al. [43] who found a water/Si ratio of 1.8 for C-S-H paste, decreasing to 1.4 upon drying in vacuum. As our model does not describe macro or capillary pores, from which water presumably is removed by the drying process, a ratio of 1.5 was considered appropriate. Further remarks on relative stabilities can be found in Section 3.2 on swelling. Our results are based on calculations of the different C-S-H structures generated for model A and model B at a water/Si ratio of 1.5, and the results for the most stable parallel simulated annealing calculation is reported in the table. During the initial optimization using molecular dynamics, the interlayer atoms undergo rather large structural changes. For model A, minor changes are found for the intralayer atoms where the SieO and CaeO chains are preserved. The water molecules remain in the interlayer for this model. For model B, on the other hand, one silicon and one calcium atom move from the intralayer, making room for hydroxyl groups. During preparation of the starting structures for DFT calculations, after the water addition step, more OH groups are formed through dissociation of water upon the initial optimization using molecular dynamics with ReaxFF. This is expected, due to the presence of unsaturated SieO and Ca species. One hydrogen atom of the H 2 O molecule attaches to the SieO species of the intralayer silicate chains, Table 1 Overview of structural properties of the energetically favoured C-S-H models from the three parallel simulated annealing runs at water/Si = 1.5.

Model
Ca whereas the OH group bonds to calcium in the interlayer. The number of OH groups is in most cases reduced somewhat after the DFT calculations, showing that ReaxFF favours splitting of water to a greater extent than DFT. Taking the limited transferability of forcefield-based approaches into account, one should place larger confidence in the DFT results.
As the models A and B have identical stoichiometry, the energy per simulation cell can be directly compared. Considering all the C-S-H structures investigated, one of the model B structures is the energetically most favourable at a water/Si ratio of 1.5 (this appears to change upon further hydration, vide infra). The three most stable (within 0.3 eV) model A C-S-H structures vary only in the positioning of the water molecules. It should be noted that while the overall structure is charge neutral, the charge distribution within the unit cell varies. Mohammed et al. [5] discussed the charge distribution in different C-S-H structures, and observed large charge fluctuations in the structures reported by Kovačević et al. [18]. They noted that the different building blocks of their structures had similar energies despite having different charge distributions, confirming the existence of charge fluctuations of the C-S-H structures. However, the magnitude of these charge fluctuations and the length scales over which they are compensated remain to be studied. The charge distribution in the models studied here can be estimated from the Ca/Si distribution in the two interlayers in the simulated unit cell as shown in Table 1. In model A, only bridging silicon tetrahedra are removed, while in model B one bridging silicon tetrahedron less is removed but one intralayer silicon is also removed. Hence model B has one more silicon atom in the interlayer compared to model A. Likewise, there is one more calcium atom in the interlayers in model B compared to model A. As a result, the charge distribution between intralayers and interlayers differ between the models. Two of the model A structures have the same number of calcium atoms (6) in both interlayers, but one silicon only in one of the interlayers. The other two models A have an uneven distribution with 5 and 7 calcium atoms in the two interlayers, respectively. Also, in one of the models the interlayer silicon is left in the interlayer with most calcium, while in the other model it is opposite. For model B there is variation in the distribution of calcium and silicon atoms in the interlayers, with ratio 5/8 or 6/7. The two models B with calcium distribution 5/8 and one of the two models with distribution 6/7 have the two remaining silicon atoms in the interlayer with 5 calcium. The last model has one silicon in each interlayer. Variations in the number of cations within the layers will influence the number of anionic species present and this may have impact on stability. No systematic effects from this variation can be identified here with regard to energetic stability. For instance, for model A, the same number of calcium in each intralayer, gives structures both at lowest and highest relative energy. As a corollary to this observation, we note the great variety of apparently stable structures found for calcium silicates, both as natural minerals and synthesized materials [2].
The characteristics of the C-S-H structures have been analysed to determine favourable configurations of calcium and hydroxyl groups (Fig. 3). Table 1 lists the structural and energies of the calculated C-S-H and Fig. 4 shows the corresponding structural properties of calcium and silicon for all C-S-H models investigated. The cut-off distance for CaeO coordination is set to 2.9 Å in accordance with Richardson [7]; for SieO, 2.0 Å was chosen. The interlayer calcium atoms mainly coordinate to 6 other atoms, albeit with a significant admixture of 5-and 7-coordination. The higher coordination of the intralayer calcium atoms is expected due to the higher degree of structural order within the intralayer, the structure of which is similar to that of the original tobermorite. Without exception, the silicon atoms coordinate to four oxygen atoms. This is in agreement with experimental findings with the average CaeO and SieO coordination numbers of 6.1 and 3.9, respectively, experimentally determined for samples with C/S ratio of 1.6 [44]. Fig. 4 shows the distribution of bonding distances of calcium, bonded to silicon and different oxygen species. The average SieO bond length is 1.65 Å, in agreement with experiments [44]. The majority of the CaeO bonds are within 2.3-2.5 Å, with an average bond length of 2.44 Å. This is in excellent agreement with experiments on crystalline calcium silicate hydrates and related phases [7,44,45] as well as theoretical findings [5,17]. It is also similar to the CaeO bond length found for portlandite (Ca(OH) 2 ), which is 2.37 Å [46]. There is a larger number of interlayer calcium atoms with shorter CaeO bond lengths compared to intralayer calcium atoms. This can be associated with the higher coordination number of the intralayer calcium atoms, which is mainly 6 or 7. The structural properties of model A at higher hydration level (water/Si ratio of 1.8 and 2.1) does not change significantly as shown in Figs. S1-S4 and Tables S1-S2 in the supplementary information.
The results reveal that calcium atoms in the interlayer are often linked to another calcium connected through OH groups forming Ca (OH) x Ca structures. Fig. 5 exemplifies the local environment of Ca-OH structures formed in the interlayers. Also Mohamed et al. [5] found Ca (OH) 2 to be the most stable structure in the interlayer. However, the combination into more complex Ca(OH) x Ca structures in the interlayer has, to our knowledge, not been described earlier. These large Ca (OH) x Ca moieties seem to occur more frequently in the energetically more stable configurations. Ca in the interlayer coordinates also to the oxygen atom of H 2 O, but to a lesser degree and with longer bonds. Ca in the intralayers coordinates exclusively to O bonded to silicon in the intralayer, staying in the initial tobermorite silicate chains. We also observe some coordination from silicon to OH groups just on the border between the intra-and interlayers.
Hydroxyl groups account for 21% of the oxygen positions coordinated to calcium for the structure presented here (Ca/Si ratio of 1.65). This corresponds well with experimental results from inelastic neutron scattering and NMR correlation [4,47]. Thomas et al. [4] found that around 23% of the charge of Ca 2+ ions are compensated by OH − groups for a Ca/Si ratio of 1.7. For silicon, the corresponding amount of Si-OH groups are 7% on average in the present work, whereas experimental results indicate that this value should be almost negligible for C-S-H structures with high Ca/Si ratios [48,49].
In conclusion regarding the structure of the C-S-H phase in cement, we recall from the Introduction that a large number of phases consisting of calcium, silicon, and oxygen have been characterized experimentally. They vary widely both in structure and in stoichiometry, and we take this as evidence for great structural diversity available within a narrow range of energies of formation. The present computational results corroborate this, as we find few trends in structure vs. stability for different model systems.

Swelling
The density of C-S-H with a Ca/Si ratio of 1.65 was investigated at three different hydration levels, i.e. water/Si ratios of 1.5, 1.8 and 2.1. Change of density and related swelling are of great interest, especially very close to the steel rebar, due to possible accumulation of local stress and enhanced formation of previously mentioned nanoscale pores and cracks. Thus, swelling might be considered as a damage accelerator and quantification becomes of great interest. Fig. 6 shows the energy per atom and the density for the different C-S-H structures. At water/Si ratio of 1.5, one of the C-S-H structures of model B is energetically favoured. However, at higher water content model A is favoured. On average, model A is favoured for all water contents. The optimized structures of the most stable structure for models A at all water contents are shown in Fig. 7. This result may be seen in light of the MASNMR results of Chen et al. [8] which indicate that our structure A is the closest model to the experimental situation, as it does not involve removal of pair sites. There is an increase in stability gained by adding more water as judged by the decreasing ΔE. The relative energies of systems with different water content reported here are calculated using a semiempirical method of potentially modest accuracy in comparing different hydration levels. Higher water/Si rations (1.8, 2.1) are more stable experimentally, and appear to be more stable here as well. On the other hand, our model does not describe capillary and macro pores, and should hence show a lower water content. We have presently no satisfactory solution to this dilemma. The densities as a function of water content should, however, not be affected by this, as they are computed entirely from first principles.
The densities of the most energetically stable structures are given in Table 2. They decrease with increasing water content, capturing the experimental trend [43,50,51]. The predicted average densities are lower than the experimentally obtained values. A density of 2.604 g/ cm 3 was found experimentally for C-S-H with composition (CaO) 1.7 (SiO 2 )(H 2 O) 1.80 [43], whereas the comparable average densities obtained here at water/Si ratio of 1.8 are 2.46 g/cm 3 and 2.44 g/ cm 3 for models A and B, respectively. These values are close to previous MD calculations using ReaxFF, where an average density of 2.49 g/cm 3 was obtained for C-S-H a similar water/Si ratio [18]. In Fig. 8, we compare the present results with earlier computational and experimental data. While force-field based results from the literature are in good agreement with the present DFT data, experiments appear to give densities that are from 5% to 8% higher.
Our results show large variations, especially in density. This may be an effect of the limited system size. Model B has a slightly lower density than model A, on average, for the two situations with lowest water content. At the highest water/Si ratio investigated, this relation is reversed.

Chloride sorption properties
Chloride sorption properties were investigated for the C-S-H structures with water/Si ratio of 1.5 and 1.8. The most stable configuration for each structure was investigated where one OH group were replaced by chloride as outlined in the Computational Details Section above.
The energetics of chloride sorption were calculated through a semiempirical approach akin to a Haber cycle to generate an estimate for a realistic reaction:    It may be expressed as the sum of the following three reactions: (1) CSH-OH − + NaCl (s) ⇔ CSH-Cl − + NaOH (s) (2) NaOH (s) ⇔ Na + (aq) + OH − (aq) (3) Na + (aq) + Cl − (aq) ⇔ NaCl (s) While the three reactions above may not take place in the system under study, adding them to arrive at a total reaction is still appropriate under the assumption that the energy is a thermodynamic state function. The chloride sorption energies given in Fig. 8 have been calculated according to where E CSH−OH+Cl is the total energy of the C-S-H model with a hydroxyl group substituted by chloride. E CSH is the total energy of the C-S-H model before substitution. E NaCl and E NaOH are the energies of solid NaCl and NaOH, respectively, as calculated with DFT. ΔE NaCl (0.04 eV) and ΔE NaOH (−0.46 eV) are the experimentally determined enthalpies of solvation [42]. Negative values of ΔE Cl indicate that the chloride substitution is energetically favourable. Calculating a Haber cycle in this semiempirical fashion cannot be relied upon to produce very accurate energies. Moreover, we disregard the effects of entropy change and deviations from unit concentrations/pressures. However, the most favourable substitutions show reaction energies as exothermic as ≈−0.85 eV, suggesting that the substitution of hydroxide with chloride at the most favourable sites in dense cement paste is exothermal.
The energetics of chloride sorption versus coordination number for models A and B are shown in Fig. 9 and structural properties are illustrated in Fig. 10. The most stable configurations of chloride are when chloride bonds to calcium atoms within the interlayer. The favoured sorption energies of chloride for model A are −0.85 eV and −0.64 eV for water/Si ratios of 1.5 and 1.8, respectively. For model B, the corresponding energies are −0.33 eV and −0.46 eV, respectively. Overall, the computations reveal that 17% of the structures investigated here form chloride species that are stable with respect to aqueous NaCl. This indicates that chloride sorption into the cement gel may play a role in the greater picture of chloride sorption in concrete. Fig. 9 displays the structural properties calculated for chloride in C-S-H, where chloride neighbouring calcium was considered within a radius of 3.5 Å, whereas a radius of 2.9 Å for silicon was employed. Chloride coordinates to up to five neighbouring atoms (Ca, Si), with chloride linked to 2 or 3 neighbouring atoms being the most abundant configuration within the structures investigated here, as illustrated in Fig. 9. It is evident that there is a large variation in the Cl-X distances (X being Si or Ca). CleCa bond lengths vary from around 2.5 Å to 3.5 Å. The majority of the energetically stable chloride configurations exhibit CleCa bond lengths between 2.7 Å to 2.9 Å. For comparison, the CleCa distance in solid CaCl 2 is 2.75 Å [52].
Chlorine atoms linked with silicon is highly endothermic in all cases (Fig. 9), suggesting that chlorine will not interact with silicon in the intralayer of C-S-H. This is in accordance with chemical intuition as well as experimental findings; no interaction between silicon and chloride was observed using NMR [53] or in previous MM/MD simulations [21,29]. Moreover, no stable bonding of chloride has been found close to the (silicon-terminated) tobermorite surface itself.
The local environment of the favoured chloride sorption configurations of model A and B are illustrated in Fig. 10, 11 and 12, where chloride is well within van der Waals distance of several calcium atoms. The most energetically favoured chloride configurations arise when chloride has two to four calcium atoms within its coordination sphere. This suggests that chloride prefers to cluster with calcium within the C-S-H interlayers. There is no obvious correlation between chloride sorption energy and number of neighbouring calcium atoms beyond two. It has previously been shown theoretically that calcium and chloride may form stable clusters in C-S-H phases [32,54], indicating a possible retention mechanism of chloride in the cement paste, which could potentially inhibit chloride diffusion. Quantification of this effect could be an interesting topic for further studies.

Conclusions
The structure of C-S-H and the effect of varying water content on its density have been investigated using density functional theory calculations. Results for densities as a function of hydration are in good agreement with experimental values, and in line with what has been found with molecular mechanics in the literature. While we observe very little variation in SieO and CaeO bond lengths between different structures, structural diversity is otherwise great. This is in accordance with experimental observations, as we see modest correlation between structural features and stability. A mapping of energetics of hydroxyl substitution with chloride reveals, unsurprisingly, that chlorine preferentially coordinates to calcium. More specifically, the most stable sites for chlorine substitution involve at least two calcium atoms. Computed chloride substitution energies indicate that the C-S-H phase will bind chloride from aqueous solution, potentially influencing chloride diffusion in concrete. Quantification of this effect and its impact on chloride transport is still lacking; a more comprehensive The corresponding lighter colours illustrate the distribution of the energetically stable chloride configurations. The bond length cut-off distance for CleCa coordination is set to 3.5 Å and the cut-off distance of CleSi is set to 2.9 Å.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.