Realistic atomic structure of fly ash-based geopolymer gels: Realistic atomic structure of fly ash-based geopolymer gels: Insights from molecular dynamics simulations Insights from molecular dynamics simulations

Geopolymers, synthesized through alkaline activation of aluminosilicates, have emerged as a sustainable alternative for traditional ordinary portland cement. In spite of the satisfactory mechanical performance and sustainability-related benefits, the large scale acceptance of geopolymers in the construction industry is still limited due to poor understanding of the composition-property relationships. Molecular simulation is a powerful tool to develop such relationships, provided the adopted molecular structure represents the experimental data effectively. Towards this end, this paper presents a new molecular structure of sodium aluminosilicate hydrate geopolymer gels, inspired from the traditional calcium silicate hydrates gel. In contrast to the existing model — where water is uniformly distributed in the structure — we present a layered-but-disordered structure. This new structure incorporates water in the interlayer space of aluminosilicate network. The structural features of the new proposed molecular structure are evaluated in terms of both short-and medium-range order features such as pair distribution functions, bond angle distributions, and structure factor. The structural features of the newly proposed molecular structure with interlayer water shows better correlation with the experimental observations as compared to the existing traditional structure signifying an increased plausibility of the proposed structure. The proposed structure can be adopted as a starting point towards realistic multiscale simulation-based design and development of geopolymers.

Presence of inorganic three dimensional alkali aluminosilicate network structure is the primary contributor to such excellent characteristics of geopolymers 7,[26][27][28][29][30][31] .As the significance of geopolymers gained more attention due to its environmental benefits, significant works have been focused on studying the local structure of alkali aluminosilicate gels (such as NASH gel) using various advanced materials characterization techniques such as x-ray diffraction (XRD) and nuclear magnetic resonance (NMR) spectroscopy [32][33][34][35] .In order to obtain detailed comprehensive insight on the local structural characteristics as well as the mechanical properties of such aluminosilicate gels, molecular dynamics (MD) simulations have been carried out 7,[36][37][38][39] .In these MD studies 7,[36][37][38][39] , the atomic structure of NASH gel has been constructed by traditional melt-quenching method 40,41 to obtain the amorphous sodium aluminosilicate glass (NAS) structure.Finally, water is adsorbed randomly into the 3D glass (NAS) skeleton to obtain the wet NASH structure.Water in the NASH structure promotes diffusion and dissociation of alkali species 42- 44 and may also control the thermal expansion behavior as in the case of calcium silicate hydrate gels in traditional hardened cement paste 45 .Hence, distribution pattern of water in the NASH structure-in particular, any possibility of presence of interlayer water in the NASH structure-is likely to influence the mechanical properties and durability characteristics of geopolymers.The importance of water at the interlayer region in cementitious gels can be traced back to the work by Feldman and Sereda 46 , where the model for hardened cement paste was developed with interlayer water.The assumption was made based on immense value of surface area in the cementitious gels which can be explained from the presence of gel water held in the interlayer similar to clay 47 .Their assumption was later confirmed by the experimental results of measurement of helium flow in hydrated Portland cement paste and hydrated tricalcium silicate 48,49 .While all the previous studies adopted random adsorption of water into the 3D NAS structure, the current study evaluates, for the first time, the possibility of presence of interlayer water in the NAS skeleton using molecular dynamics simulations.Towards that end, a vacuum layer is introduced into the 3D dry NAS structure followed by adsorption of water into the created interlayer space using a Grand Canonical Monte Carlo (GCMC) simulation 50 to obtain a realistic 3D structure of NASH with interlayer water.The reliability of molecular simulation depends significantly on the ability of the inter-atomic potential to model the simulated structure 40,[51][52][53] .In this study, ReaxFF potential 54 has been adopted, which accounts for the potential energy based on the local structure around each atom 55 and has been shown to yield promising results towards realistic interactions of ordered and disordered structures 36,[56][57][58][59][60] .The structure of NASH gel with interlayer water is likely to address the leaching issue in geopolymers in a more efficient manner (further investigation is needed to confirm this) as compared to the structure with random distribution of water where significant caging effect prohibits realistic representation of durability characteristics of geopolymers.The characteristics of the proposed NASH structure are compared with experimental observations as well as with the features of traditional simulated structure of NASH without any interlayer water with a view to evaluate the possibility of the existence of interlayer water in the NASH structure.

SIMULATION METHODOLOGY
This section presents the model construction procedure as well as various measures of characterization of atomic structure adopted in this study to evaluate the proposed atomic structure of NASH.

Model Construction
.Here Si/Al ratio and Na/Al ratio are considered 3 and 1 respectively.The unit cell is replicated three times along x and y axes and duplicated along z axis to generate the simulation box.The initial structure, thus obtained, is melted at 4000 K and 0 atm in the isothermal isobaric (NPT) ensemble for 400 ps to make sure that there is a memory loss of the atomic structure from the initial configuration.The structure is then quenched linearly from 4000 K to 300 K with a cooling rate of 1K/ps in NPT ensemble so as to obtain the glass structure.Note that a similar cooling rate is adopted by Oey et al. 63 The glassy structure, thus formed, is then relaxed at zero pressure and 300 K for 400 ps in NPT ensemble, followed by NVT ensemble for 400 ps.For statistical averaging, the structure is run for another 200 ps.In order to create an interlayer space for water adsorption, a vacuum layer of thickness 9 Å is introduced in the z-direction at the midthickness region.Layer of similar thickness has been adopted for modelling the atomic structure of calcium silicate hydrate (CSH) 58,64,65 .To hydrate the structures, a Grand Canonical Monte Carlo (GCMC) simulation 50 is performed in the grand canonical ensemble (VT).In this study, the chemical potential  for the fictitious water reservoir is set to 0 eV at a temperature of 300 K to provide unlimited supply of water.The simulation is continued till the structure is saturated.Figure 2 shows the equilibrium run in the  The saturated structure, thus obtained, is then equilibrated for 400 ps in NPT ensemble and the statistical average is calculated for another 200 ps in NVT ensemble to obtain the final NASH structure with interlayer water as shown in Figure 3(a).For comparison purpose, a traditional NASH structure with random distribution of water (without any interlayer water), similar to the structures reported in literature 7,[36][37][38][39] is also generated.In this case, traditional method of assigning random position to all the atoms (Si, Al, Na and O) in the bounding box is followed keeping the values of the Si/Al ratio and Na/Al ratio 3 and 1 respectively.Hereafter, similar melting-quenching procedure is followed as explained earlier (refer to Figure 1) to obtain amorphous NAS structure.Afterwards, the water is adsorbed into the NAS structure using GCMC 50 as explained earlier.In this case, water gets adsorbed into the open spaces inside the NAS structure.The NASH structure, thus generated for comparison, is shown in Figure 3(b).The MD simulation is performed using an open-source software, LAMMPS 66 .Verlet method 50 is applied to integrate the trajectories with time step of 0.5 fs.ReaxFF potential 54 is used for the inter-atomic interactions along with a charge equilibrium (qeq) approach 67 .The simulated bulk density for the proposed and traditional NASH structures are 1.825 g/cm 3 and 2.12 g/cm 3 respectively which are in line with values reported in the literature 38,68,69 .

Interatomic Potential
The ReaxFF potential is used in this study for the atomic interaction in the molecular structure based on the reactive force field approach 59 .In ReaxFF, the bond order formation is implemented using the parameters that are tested in the quantum chemical data such that the bond forming and breaking processes in a chemical reaction can be achieved.The ReaxFF includes Coulombic interactions, nonbonded van der Waals, and bonded interactions.ReaxFF has been successfully implemented for MD simulations and the features of the simulated atomic structures (using ReaxFF) have been shown to correlate well with the experimental observations 58,[70][71][72][73][74][75][76] .The ReaxFF file has been added as a supplementary material.

Structural Characterization
This section describes the structural characterization measures adopted to evaluate short-range and medium-range order of the NASH structure.

Pair distribution function
To extract the in-depth information on the local structure of the material, the pair distribution function (PDF) 50 is adopted.PDF represents the probability of finding an atom at a given distance 'r' from the reference atom.Mathematically, it is the ratio of the local density of atoms with respect to the global density of the atoms expressed as a function of distance from an atom.In particular, neutron PDFis given as 53,58,77 : where   is the fraction of  atoms ( = Al, Si, Na, or O),   is the neutron scattering length of the species, and   are the partial PDFs.Using the partial PDF, coordination number of each of the species can be obtained by computing the number of neighbors within the first coordination shell of the respective atoms.The cutoff distance is taken from the first minimum of the respective partial PDFs.The PDF is mainly applicable for investigating the properties of the structure in short range order (< 3 Å).

Structure Factor
The structure factor is adopted here to evaluate the medium-range properties of the atomic structure.
The partial structure factors are calculated from Fourier transformation of the partial PDF   () as shown in Equation 2.
where is the scattering vector,  0 is the average atom number density, and  is the maximum value of the integration in real space, which is set to half of the size of one side of the simulation cell.The total structure factor is calculated as: where   and   are the fractions of atoms and   and   are neutron scattering lengths, for elements i and j, respectively.The partial PDF and structure factor are plotted by taking statistical average over 100 frames at 300 K.

Local structure and short-range order
In this section, the short-range order (< 3 Å) of the NASH structure is evaluated through the total pair distribution function and partial pair distribution function.

Total pair distribution function
Figure 4 shows the computed total PDF of NASH structure with interlayer water.In addition, the experimentally obtained PDF using neutron diffraction (plotted against secondary/right y axis) as well as simulated PDF for the NASH structure without interlayer water (random distribution of water) is also plotted in Figure 4 for comparison.Here, first sharp peak at 1.0 Å corresponds to D-O interactions (where D refers to deuterium).In the experimental study 32,33,35 , heavy water was added to prepare the geopolymer binder.It needs to be noted that the first peak at 1.0 Å obtained here correlates very well with the experimental data 32,33 as well as simulations without interlayer water.In order to have a meaningful comparison of the computed pair distribution with the neutron experimental data, the PDF is broadened using the methodology described by Wright et al. 78 .Figure 5 ](4) where (  ) is the experimental total PDF.These factors are calculated over the range from 1 Å to 10 Å and the computed ℛ  for NASH with interlayer water and without interlayer water are found to be 5.79% and 5.87% respectively.Hence, in terms of PDF both the structures show good agreement with the experimental data (since the value of ℛ  is below 9 % ).To shed more light, the degree of agreement with the experimental data is also presented in terms of structure factor in a forthcoming section for both the NASH structures presented here.

Partial pair distribution function
In this section, the partial PDF is analyzed.we observe a minor peak at 2.0 Å corresponding to the partial Na-O PDF of both NASH with and without interlayer water.In particular, the peak is more pronounced for NASH structure with interlayer water.
This could be attributed to the clustering of Na atoms around the water molecules in addition to the formation of non-bridging oxygen atoms.To shed more light on this, Table 1 reports the percentage of oxygen molecules connected to various bond angle configurations for a cutoff distance of 2.12 Å (obtained from the minima after the minor peak of Na-O PDF).It is evident from Table 1 that the fraction of oxygen   The partial PDFs of Na-Na, Si-Na, Al-Na, and Na-H are reported in Figure 7(a), (b), (c) and (d), respectively, so as to evaluate the influence of interlayer water on cation interactions.It is observed that the width and position of the peaks corresponding to Na-Na, Si-Na, Al-Na, and Na-H interactions are not significantly affected by the incorporation of interlayer water, thus justifying the possibility of interlayer water in the NASH structure.To shed more light on this, the forthcoming sub-section evaluates the influence of the incorporation of the interlayer water in terms of the bond angle distributions.

Bond angle distributions (BAD)
A partial pair distribution in terms of bond angles is presented herein.Figure 8 shows the intra-tetrahedral bond angle distributions.Figures 8(a The BAD for intra-tetrahedral O-Si-O (Figure 8(a)) bond angle shows a peak at 107 o for all the structures, which is comparable to the value reported for silicate glass 51 .It should also be noted here that the presence of water does not influence the intra-tetrahedral O-Si-O BAD.Similar observations are also found in the literature 36 This suggests that the Si tetrahedra is not affected or distorted by the presence of water molecules.In the case of O-Al-O bond angle (Figure 8 Figure 9(a) shows the Si-O-Si bond angle distribution.While the dominant peak is observed at the same position (135 o ) for dry NAS glass as well as the NASH structure with interlayer water, the width of the peak increased slightly for NASH structure with interlayer water reflecting slight increase in the degree of disorder in the system with water-incorporation, which is expected.Moreover, the width of the peak increases significantly as compared to dry NAS when water is adsorbed into the entire NASH system without any interlayer water.The difference in the width can be explained from the fact that the access to the Si-O-Si bond for water is limited in the NASH structure with interlayer water resulting in limited interactions between water and Si-O-Si bonds.This results in a slight increase in width of the peak as observed in that case.On the other hand, the interactions between Si-O-Si bonds and water increases significantly for the NASH structure without interlayer water since water gets better access to the whole of Si-O-Si bonds in the system.This explains the higher degree of disorder and presence of significantly wider Si-O-Si peaks in NASH system without interlayer water The shoulder peaks observed at 110 o (approx.) for NAS glass can be attributed to the Si-Si steric repulsion 53,71,79 .Presence of both these peaks are in line with the observations reported in the literature 80 .While no significant change in the width of the major peak is observed when water is incorporated in case of both the NASH structures, the width of the secondary peak (at approximately 100 o ) increases.The increase in width of the secondary peak is higher for the NASH structure without interlayer water as compared to the NASH structure with interlayer water as expected, similar to the trends observed in the Si-O-Si bond angle distributions for the difference in the distribution of water in these two NASH structures as Overall, the good correlation of the short-range order of the NASH structure with the experimental data validates the plausibility of the NASH structure with interlayer water proposed in this study.To elaborate more on the plausibility of the proposed NASH structure with interlayer water, the forthcoming section assesses the medium range order of NASH in terms of structure factor.

Medium range order
The structure factor represents correlations at larger distances (between 3 Å to 10 Å) as compared to PDF which analyzes structure correlations at smaller distances (≤ 3 Å).In this section the neutron structure factor is computed using Equation 3 and reported in Figure 10.From Figure 10, it is clear that the structure factor for NASH with interlayer water contains four distinct peaks, the position of which correlates very well with the experimental observations.In particular, the first sharp diffraction peak (FSDP) signifies the extent of correlation in the medium-range structure of glassy/disordered structures 52,53,77 .The position of FSDP for the proposed structure exhibits excellent match with the experimental observation signifying correlation the medium-range structure.This is in contrast to the NASH structure without interlayer water, where the first peak is notably shifted to a higher wave vector.Moreover, the position of the first two peaks in case of NASH with interlayer water exhibits better correlation with experimental observation as compared to the traditional NASH structure without interlayer water.This difference in FSDP between the two structures presented here originates from their nature of water distributions (refer to supplementary materials document for more details).To shed more light, the degree of agreement between the simulated and experimental structure factor is quantified using the Wright factor (ℛ  ) 78 (equation 4).While the NASH structure with interlayer water exhibits a value of 9.3% (ℛ  ), the structure with random dispersion of water without any interlayer water shows a significantly higher value of 13.8% (ℛ  ).Hence, it can be concluded that the proposed structure with interlayer water shows better correlation with the experimental data as compared to the one without interlayer water.Overall, a good match between the experimental and simulated structure factor peaks confirms the plausibility of the NASH structure with interlayer water, presented in this paper.

Figure 1 :
Figure 1: Schematic representation of the model construction procedure (Colors scheme -Al: Green, Si: Blue, Na: Yellow, O: White, and H: Red) The proposed atomic structure of the NASH is derived in this study from a unit crystal cell of albite mineral 61 (chemical composition: NaAlSi3O8) with unit cell parameters (a, b, c) = (4.803Å, 8.327 Å, 8.78 Å)62 .Here Si/Al ratio and Na/Al ratio are considered 3 and 1 respectively.The unit cell is replicated three simulation.Here, the equilibrium is reached after approximately 2 x 10 5 steps beyond which the number of adsorbed water molecules fluctuates around a constant value.

Figure 2 :
Figure 2: Number of water molecules adsorbed with increasing GCMC steps

Figure 3 :
Figure 3: Atomic structure of: (a) NASH gel with interlayer water, (b) NASH gel with random dispersion of water into the structure.
The following three peaks correspond to Si-O, Al-O and O-O interactions, respectively, which also show good correlation with experimental data32,33 as well as simulations without interlayer water.Such close match in the position of the peaks suggests that first neighbor arrangements of the proposed atomic structure of NASH with interlayer water is in good agreement with the experimental data signifying the possibility of presence of interlayer water in the NASH structure.In order to confirm this further, the forthcoming sub-section evaluates X-O partial PDFs (X = Al, Na, Si, H) of the proposed structure with respect to the experimental data as well as traditional NASH structure with random distribution of water for detailed insight.

Figure 4 :
Figure 4: Neutron pair distribution function (PDF) predicted for NAS glass, NASH gel with interlayer water, and NASH gel with random dispersion of water into the structure

Figure 5 :
Figure 5: Neutron pair distribution function (PDF) predicted for (a) NASH gel with interlayer water, and(b) NASH gel with random dispersion of water into the structure (PDF is broadened using the methodology described by Wright et al.78 ) Furthermore, the Wright factor is computed to compare the degree of agreement between the computed pair distribution function (PDF) and the experimental data for both the structures.The Wright factor, ℛ  is expressed as78

Figure 6 (
a), (b), (c) and (d) plots the Al-O, Na-O, Si-O, and H-O interactions respectively.No significant difference in the position and width of the peaks for the X-O partial PDFs (X = Al, Na, Si and H) are observed with incorporation of interlayer water.Major peaks for Al-O, Na-O, Si-O, and H-O interactions are observed at 1.8 Å, 2.3 Å, 1.62 Å and 1 Å, respectively, and the positions of the peaks are in good agreement with the values reported in the literature[32][33][34][35] .Interestingly, atoms contributing towards Si-O-Si, and Si-O-Al bonds decreases and significant fraction of oxygen atoms connecting to hydrogen atoms (X-O-H where X: Si, Al, H) is being observed when water is adsorbed into the NAS structure.Hence, the presence of the minor peak at 2.0 Å can be primarily attributed to the incorporation of water in the structure.A detailed discussion on the bond angles is presented in a later section.

Figure 6 :
Figure 6: Partial pair distribution function (PDF) of (a) Al-O, (b) Na-O, (c) Si-O, (d), H-O for NAS glass, NASH gel with interlayer water, and NASH gel with random dispersion of water into the structure

Figure 7 :
Figure 7: Partial pair distribution function (PDF) of (a) Na-Na, (b) Si-Na, (c) Al-Na, and (d) Na-H for NAS glass, NASH gel with interlayer water, and NASH gel with random dispersion of water into the structure

Figure 8 :
Figures 8(c) and (d) shows the BAD for Al-O-Al corresponding to coordination numbers 4 and 5 respectively for Al atoms.
(b)) two peaks were observed.This could be attributed to the presence of five-coordinated Al atoms in the structure as shown in Figure 8(d).While the peak at 107 O (O-Al-O BAD) is contributed by both four and five-coordinated Al atoms (Figures 8(c) and 8(d)), the minor peak at 162 O stems from the presence of five coordinated Al atoms only (Figure 8(d)).However, for the O-Al-O BADs, no significant changes are observed between NAS, NASH without interlayer water, and NASH with interlayer water signifying that the intra-tetrahedral O-Al-O bonds are not affected by the incorporation of water in the structure, similar to the trends observed for O-Si-O BAD.

Figure 9 :
Figure 9: Bond angles distribution of (a) Si-O-Si, (b) Si-O-Al, and (c) Al-O-Al for dry NAS glass, NASH gel with interlayer water, and NASH gel with random dispersion of water into the structure.

Figure 9 (
Figure 9(b) shows the bond angle distribution for Si-O-Al.The dominant peak is observed at approximately 130 o for all the three structures.A second peak is also observed at approximately 100 o .
explained earlier.In the case of Al-O-Al bond angle (Figure 9(c)), a bimodal distribution of the angles is observed which is in line with the observations reported elsewhere 80 .Such bimodal distribution in Al-O-Al could be attributed to the presence of five coordinated Al atoms and edge sharing tetrahedra of Al.Five coordinated Al atoms have been shown in the literature 81 to yield stable structure for Al.While the presence of five coordinated Al atoms in all the three structures is confirmed in this study as explained previously in this paper (refer to Figure 8(c) and 8(d)), R.I.N.G.S code 82 is employed to evaluate the presence of edge-shared tetrahedra.It is observed that for both NAS glass and NASH without interlayer water, only corner-shared tetrahedra is present whereas, in the case of NASH structure with interlayer water, 2 (10%) edge-shared and 18 (90%) corner-shared tetrahedra are observed.This explains the significant difference in the Al-O-Al bond angle distributions observed between the two NASH structures presented in the paper.It should be noted that the Al-O-Al angles in these structures exist only in a few quantities close to 1%.This is in agreement with the Loewenstein rule 83 which suggests that heteropolar Al-O-Si angles are preferred at the expense of homopolar Al-O-Al angles in aluminosilicate minerals.

Figure 10 :
Figure 10: Neutron structure factor predicted for NAS glass, NASH gel with interlayer water, and NASH gel with random dispersion of water into the structure

Table 1 :
percentage of oxygen molecules connected to various bond angle configurations for a cutoff distance of 2.12 Å in the NAS, NASH without interlayer water, and NASH with interlayer water