Interparticle Interactions in Colloidal Systems: Toward a

Jul 14, 2017 - We introduce a simple yet effective mechanism to capture the sliding behavior of colloids with surface roughness. ... Molecular Simulat...
3 downloads 21 Views 8MB Size
Subscriber access provided by UNIV OF DURHAM

Article

Interparticle interactions in colloidal systems: towards a comprehensive mesoscale model Saeed Masoumi, Hamid Valipour, and Mohammad Javad Abdolhosseini Qomi ACS Appl. Mater. Interfaces, Just Accepted Manuscript • DOI: 10.1021/acsami.7b08465 • Publication Date (Web): 14 Jul 2017 Downloaded from http://pubs.acs.org on July 18, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

ACS Applied Materials & Interfaces is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

Interparticle Interactions in Colloidal Systems: Towards a Comprehensive Mesoscale Model

Saeed Masoumi,



Hamid Valipour,



and Mohammad Javad Abdolhosseini Qomi

∗,‡

†Centre for Infrastructure Engineering and Safety (CIES), School of Civil and Environmental Engineering, UNSW Australia, UNSW Sydney, NSW 2052, Australia

‡Advanced Infrastructure Materials for Sustainability Laboratory (AIMS Lab), Department of Civil and Environmental Engineering, Henry Samueli School of Engineering, E4130 Engineering Gateway, University of California, Irvine, Irvine, CA 92697-2175 USA. E-mail: [email protected]

Abstract Intermolecular interactions control the collective behavior of colloidal clusters. The overwhelming majority of literature focuses on cohesive attributes of intermolecular forces as it governs the jamming process. Yet, the overlooked sliding component plays a critical role in the slow relaxation dynamics of colloidal aggregates and the emergence of discontinuous shear thickening in dense suspensions. Here, we use crystalline calciumsilicate-hydrate (C-S-H) as a model system to explore synergistic cohesive and sliding interactions. We use the free energy perturbation approach to reconstruct potential of mean force prole between two nite-sized nanolayers in an aqueous environment. We show that sliding free energy barriers decay exponentially with increasing separation distance. The characteristic length scale of the decay is related to the interface corrugation. We introduce a simple yet eective mechanism to capture the sliding behavior of colloids with surface roughness. Moreover, we develop a global free energy landscape model by juxtaposing cohesive and sliding interactions. This model enables 1

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

us to measure the height of energy barriers which is essential to elucidate deformation mechanism and dynamics of colloidal aggregates. For cohesive colloids, our approach predicts a sublinear relation between applied normal and shear stress at the onset of sliding that is in contrast to Amontons' laws of friction. We demonstrate that the sublinear trend is due to the cohesion and nature of soft contact at the nanoscale. The proposed framework provides a new route to draw a more realistic picture of intermolecular interactions in nanoparticulate systems such as geomaterials, cementitious systems, complex colloidal assemblies and dense suspensions. keywords: free energy landscape, cohesive interaction, sliding free energy barrier, mesoscale potential of mean force, sublinear friction, corrugated surface

Introduction The intermolecular interactions govern the jamming, 1,2 rheology, 3,4 and dynamical behavior 5 of colloidal systems. These intermolecular interactions are conventionally studied in terms of normal and tangential surface forces. The normal interactions control the dispersion and agglomeration of colloidal particles via attractive and repulsive intermolecular forces and hence are responsible for jamming phenomenon. 69 The tangential component of intermolecular interactions are studied in the context of friction. Recently, researchers have directly observed that frictional forces are at the center of discontinuous shear thickening phenomenon in dense suspensions. 10 This motivates a deeper understanding of the relation between friction and normal forces in colloidal systems at the mesoscale. 11 In non-cohesive colloidal contacts, these two forces are related linearly following Amontons' rst friction law. 12 However, observations of nonlinear behavior in cohesive systems disproved the continuum description of friction at the molecular level and undermines its application to study colloidal systems at the mesoscale. 1315 Despite the progress in the understanding of intermolecular and surface forces, interparticle sliding forces have been ignored until only recently in mesoscale simulations 16,17 and are currently limited to the continuum-level description of friction forces via 2

ACS Paragon Plus Environment

Page 2 of 42

Page 3 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

Amontons-Coulomb laws. In this paper, we aim at studying the interaction between colloidal nanoparticles by considering synergistic cohesive and sliding contributions using free energy perturbation (FEP) approach. We use calcium-silicate-hydrates (C-S-H), the main binding phase in cementitious systems 1827 merely as a model system in this study. Even though the presence of C-S-H in the built environment dates back to two millennia ago in Roman concrete, 28 it was not until the last few years that researchers found that the nanogranular structure of C-S-H possesses cohesive-frictional attributes. 26,2931 This makes C-S-H a truly ideal system to study combined cohesive and sliding eects in colloidal systems. While the cohesive properties of C-S-H have been eectively studied using primitive models, 32,33 such models fail in the limit of small distances and high surface charge densities. In these cases, the local eects are signicant, and properties of conned solution becomes strikingly dierent from that of bulk. 3436 This encouraged recent full atomistic simulations studies of cohesive properties of C-S-H via grand canonical Monte Carlo 37 and molecular dynamics (MD). 38 Despite these eorts, a wholistic picture of free energy landscape (FEL) in cohesive-frictional colloids remains elusive in all colloidal systems.To this end, this paper aims to pave the route to establish a comprehensive approach toward modeling cohesive-frictional colloids at the mesoscale. The present paper is structured as follows: First, we study the origin of cohesive interactions between nite-sized crystalline C-S-H nanolayers along predened interaction paths using the FEP approach. We further investigate the possibility of describing normal interactions by modifying the long-established Gay-Berne potential. 3941 Next, we explore fundamental mechanisms of resistance against sliding and describe the sliding energy barriers in the context of nanoscale surface roughness. Furthermore, we elaborate on our proposed mathematical framework to describe sliding interactions. Finally, we reconstruct the global FEL by superimposing the normal and sliding interactions and subsequently employ it to better understand cohesive-frictional nature of colloids upon application of an external stress.

3

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Molecular Models and Simulation Details C-S-H at the low calcium-to-silicon ratio (C/S) has a crystalline structure that is akin to that of Tobermorite. We use the crystal structure of Hamids's 11 Å Tobermorite unit cell 42 as a starting point. We curtail the large tilt factor in the unit cell by applying a hexagonal-toorthogonal transformation. Then we produce a supercell of 2×2×1 Tobermorite. This leads to a layer with nal lattice parameters of a=26.54 Å, b=48.87 Å, c=12.20 Å, α = 82.5◦ , β =

90◦ , γ = 90◦ . To take into account the nite-size eects, we remove the periodic boundary condition and subsequently correct ending sites such that silicon atoms are tetrahedrally coordinated (see Supporting Information for further details). The extra charge , as a result of correcting coordination number of the ending sites, is balanced by introducing calcium counterions. This manipulation increases the calcium-to-silicon ratio (C/S) of the layer to 1.05. We subsequently place two crystalline C-S-H nanolayers, consist of 1968 atoms, inside a large simulation box and saturate the medium with simple point charge (SPC) water molecules(see inset of Figure. 1 a, c). 34,43 We keep sides of the interlayer pore open for realistic exchange of counterions in the aqueous medium. To reduce the number of degrees of freedom in the model and hence the computational cost associated with FEP calculations, we constraint atomic vibrations inside nanolayers, as entropic contributions of such vibrations are not a signicant part of intra-molecular interactions. We also constrain the internal vibrational degrees of freedom in water molecules using SHAKE algorithm. 44 However, we allow charge-balancing calcium divalent cations and water molecules to move during MD simulations freely. Transferable CSH-FF potential that operates based on partial atomic charges describes the interatomic interactions in the system. 45 Our approach to calculating free energy dierence comprises a series of parallel MD simulations at dierent center-to-center distances along a predetermined interaction path. After the MD simulation reaches equilibrium in a given state, we begin an uncorrelated simple overlap sampling (SOS) with forward-backward perturbations. In SOS approach, one ideates an imaginary intermediate state and performs perturbations with regards to it. 4

ACS Paragon Plus Environment

Page 4 of 42

Page 5 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

We assume a linear variation of the system Hamiltonian between two successive states. To ensure proper convergence of FEP calculations, we choose a perturbation step size of 0.125 Å. In canonical ensemble (NVT), we can compute the free energy dierence between two successive states, i and i+1, as, 46

∆Ai,i+1

  hexp(−β∆Ui,i+1 /2)ii 1 = − ln β hexp(−β∆Ui+1,i /2)ii+1

(1)

where A corresponds to the so-called Helmholtz free energy, β = 1/kB T , kB , T are respectively Boltzmann constant and the absolute temperature. Also, ∆Ui,i+1 is the total potential energy dierence between states i + 1 and i. 38 To properly sample important regions, we build a catalog of MD trajectories and potential energies in an NVT ensemble at T = 300 K using Nosé-Hoover thermostat. 47 At each state, we run a 2 ns-long MD simulation, in which we take the rst 0.5 ns as relaxation period and use the rest for sampling purposes. This provides 1500 frames and 3000 associated perturbations to calculate the free energy dierence in eq. 1. We readily calculate the PMF by sequentially summing the corresponding free energy dierences between successive states with respect to an arbitrary reference point. As discussed above, FEP technique in conjunction with full atomistic simulations allows us to acquire cohesive-sliding interactions between two colloids. However, considering Np as the average number of perturbation states in a given degree-of-freedom (DoF) while keeping the other DoFs xed, this corresponds to Np6 independent MD simulations to reconstruct FEL. For a reasonable number of intermediate perturbation states, Np ≈ 100, such FEL reconstruction requires 1012 of long statistically-independent MD trajectories, which proves to be computationally infeasible. Consequently, we reconstruct the FEL by interpolating a parametric analytical function to a set of free energy proles extracted from FEP simulations along independent predened coordinates. Furthermore, we verify the robustness of our proposed reduced-order energy landscape model by comparing its predictions against other interaction paths.

5

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 42

Normal Component of Free Energy We can describe the relative spatial arrangement of any two rigid layers in terms of their center-to-center distance vector, r 12 , and three relative Euler angles, ω12 = (α12 , β12 , γ12 ). This yields a six DoF model for normal interactions in the absence of surface roughness eects as follows,

Anrel = Anrel (r 12 , ω12 )

(2)

where Anrel is the normal free energy component in the global FEL. To construct Anrel , we perform FEP in face-to-face (FTF) and edge-to-edge (ETE) directions. Since interactions between layers become increasingly weaker at large center-to-center distances, we assign a value of zero to free energy at largest distance in our analysis and determine the free energy with respect to this reference state. Fig. 1.(a) shows the calculated free energy prole required to bring two nite-sized Tobermorite layers toward each other in FTF scenario as depicted in the inset. Recently, Masoumi et al. 38 demonstrated that intermolecular forces between innite crystalline C-S-H layers in FTF interaction are governed by three distinct terms, namely entropic solvation, combined attraction, and steric repulsion. The oscillatory entropic solvation contributions are attributed to characteristics of nano-conned water and in particular, its ordering adjacent to the hydrophilic surface of Tobermorite. 38,48 The presence of ordered water molecules at the surface of hydrophilic layers arises solvation forces. Thermodynamic properties of these water molecules are dierent from that of bulk and consequently, continuum theories fail to capture their realistic contribution. We discuss the origins of attractive and repulsive interactions in the light of non-bonded interactions, i.e. the sum of Lennard-Jones and electrostatic energies, between calcium-silicate layers and charge-balancing calcium counterions. As illustrated in Fig. 1.(b), the attractive interactions emanate from the interplay between the negatively charged calcium-silicate layer and the charge-balancing calcium counterion surrounding the other layer. Besides, the repulsion occurs between similarly charged species, either calcium counterions or calcium-silicate lay-

6

ACS Paragon Plus Environment

Page 7 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

ers. At short distances, region (I) of Fig. 1, colloidal surfaces overlap and a strong repulsion is induced as surfaces come to contact. 9 This can be observed from the fast growing repulsion in the Layer-Layer interaction. On the other hand, owing to the corrugated pattern of the surface, interlayer calcium atoms condensate on the surface without actual contact to each other, and this leads to a slow growing repulsion in Ca-Ca interaction. Also, CaLayer attraction grows mildly due to van der Waals interaction and counterion condensation eect. 4951 . In region (II), the Ca-Layer attraction outweighs the repulsive contributions, which explains the dominant cohesive behavior of between layers. This trend extends in the region (III) with an exclusion of the counterion condensation eect. We also conduct same FEP strategy to work out the OMF in ETE direction. Fig. 1.(c) depicts ETE interaction per area between Tobermorite layers. For the sake of comparison, it is worth mentioning that for up-scaling, the FTF interaction usually scales with the square of length [L]2 while the ETE interaction scales with length [L]. In total, FTF interaction in Tobermorite with realistic size is signicantly stronger than ETE interaction. 52 Having access to FTF and ETE free energy proles, we aim to reconstruct the cohesive FEL for nite-sized layers via an analytical function. We note that our nanolayers are parallelepiped in which the layer thickness is considerably smaller than the other two dimensions. Therefore, we can assume oblate ellipsoids to be an acceptable approximation of our nitesized Tobermorite nanolayers. Gay-Berne (GB) potential presents a framework to describe a 6 DOF interparticle interactions between ellipsoidal objects and has been extensively used for coarse-grained computations at the mesoscale. 3941 Inspired by GB approach, we propose the following functional form,

Anrel (r 12 , ω12 ) = η12 χ12 Anrel

7

ACS Paragon Plus Environment

(3)

ACS Applied Materials & Interfaces

2

100

Free Energy Perturbation Analytical Model 50 2

E (KJ/mol.Å )

2

(KJ/mol.Å )

0

A

dcc

n

rel

-2

0

-50

(I)

(a)

12

14

(KJ.mol.Å )

2

rel

0

30

32 dcc(Å)

34

18

(d)

ζ=6Å ζ=4Å

-2 28

16

θ

-0.4

(c)

(III) 14 dcc(Å)

2

n

dcc

12

4

0

-0.2

(II)

10

A

rel

(b)

Free Energy Perturbation Analytical Model

2

(KJ/mol.Å )

18

16 dcc(Å)

0.2

n

Ca-Ca Ca-Layer Layer-Layer

Condensation Region

-4

A

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 42

0

Free Energy Perturbation Analytical Model a = Req Analytical Model a = 1.7 Req

ζ=2Å

5

10

15 θ (°)

20

25

30

Figure 1: Normal interactions between two nite-sized Tobermorite layers based on free energy perturbation calculations. (a) The potential of mean force for face-to-face interaction. The solid line is the tting exercise based on eq. 3. The inset shows a schematic representation of Tobermorite layers. Red, purple and green circles respectively indicate oxygen, intralayer and interlayer calcium atoms. The silicon-oxygen bonds are shown as yellow rods. (b) Decomposition of the internal energy content between Tobermorite layers in terms of non-bonded interactions between interlayer calcium atoms and calcium-silicate layers at varying center-to-center distance. (c) Testing the tted normal potential of mean force when two layers approach each other in the edge-to-edge direction. The inset shows the schematics of the approaching scenario. (d) Testing the delity of the tted normal potential of mean force for the out-of-plane angular interaction of layers. The inset schematically shows the out-of-plane rotation. The analytical model in eq. 3 with two equivalent ellipsoidal radii is compared with free energy perturbation calculations at dierent rotation angles and center-to-center distances.

8

ACS Paragon Plus Environment

Page 9 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

where η12 denes the strength of the interaction based on the relative orientation of platelets,

 ν2 2s1 s2 η12 = det (G12 ) p si = (ai bi + c2i ) (ai bi ) 

(4)

G12 = RT1 S 21 R1 + RT2 S 22 R2 where Ri is the transformation matrix between the local and global coordinates, S i =

diag(ai , bi , ci ) is the structure matrix of particle i with radii of (ai , bi , ci ), and ν is an empirical tting constant. For an oblate particle we have ci < bi ≤ ai . The other factor in eq. 3,

χ12 , takes into account the anisotropic relative well-depth of the FTF ( ai ) and ETE (bi , ci ) interactions in between particles i = 1, 2 considering their relative orientations,

 T −1 µ χ12 = 2ˆ r 12 B 12 rˆ12 B 12 =RT1 E 1 R1 −1

−1

+

(5)

RT2 E 2 R2

−1

ˆ 12 = r 12 /|r 12 |, E i = diag(ai µ , bi µ , ci µ ), and µ is an empirical constant. The where r distant-dependent nature of interparticle interactions in eq. 3 is reected in Anrel (r 12 ). As discussed earlier, we can decompose the intermolecular interactions between two nano-sized Tobermorite layers into attractive, repulsive and entropic solvation terms. We nd that (12,6) Lennard-Jones potential along with a decaying oscillatory term can adequately t our FEP data. Therefore, we use this mathematical form as it has well-described features for a 6 DOF interaction between two layers. 52 The decaying oscillatory term describes the presence of solvtion forces. 8,9,53,54 Therefore, we put forward the following functional form for the distance-dependent interaction,

r) = 4a1

Anrel (

(

a2 h12 + a2

12

 −

a2 h12 + a2

6 )



+ a3 e

h12 a4

 sin

2π (h12 + a5 ) a4



(6)

where a1 is related to the well-depth of free energy, a2 is a geometrical parameter depending 9

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 42

on the shape and coordinates of the layers, a3 is the magnitude of entropic solvation contribution to the total free energy, a4 is the periodicity of oscillations due to entropic eects of liquid ordering adjacent to layers, and a5 is the delay parameter in the oscillatory behavior of solvation forces. Also, h12 = |r 12 | − σ12 is the anisotropic distance based on the closest approach of particles. 40 For an ellipsoidal particle, σ12 is dened as,



σ12

1 ˆ12 = rˆT12 G−1 12 r 2

− 12

(7)

where the Lennard-Jones term in eq. 6 is zero at |r 12 | = σ12 . Here, we t the analytical model in eq. 6 to our FTF and ETE free energy proles. Based on FEP results in Figs. 1, we take values of 11.2 Å and 28.6 Å for σ12 and ci = 1 and ai = bi = 0.08. for relative energy well-depth in FTF and ETE directions, respectively. In our tting procedure, we place the majority of tting weight on capturing FTF free energy prole as it governs the cohesion in cementitious systems. The empirical coecients µ = 1 and ν = 8 are kept constant during the tting process. Our tted coecients for a particle with a surface area of S in eq. 6 are

η12 × a1 /S = 4.88 KJ/mol.Å2 , a2 = 5.91Å , η12 × a3 /S = −0.2 KJ/mol.Å2 , a4 = 2.6 Å, and a5 = −0.73 Å. To nd a1 and a3 , one needs to prespecify dimensions of ellipsoidal particles. Figs. 1 demonstrate the prociency of the unied analytical model in eq. 3 to describe free energy prole in FTF and ETE directions, respectively. The analytical formula in eq. 3 provides access to the normal interaction between two nearly oblate Tobermorite layers. To examine the appropriateness of our proposed model to predict PMF along an arbitrary path for which it is not trained for. Toward this end, we perform angular perturbation of nanolayers in FTF case to work out its PMF. This is especially crucial since the rotation of layers is a shape dependent problem. We dene the maximum angle of rotation as the characteristic angle in which the layers come to contact and is always θmax ≥ arcsin(ζ/l), where ζ is the interlayer separation distance in FTF case, and l is the farthest distance of the body from the center of rotation and perpendicular to

10

ACS Paragon Plus Environment

Page 11 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

the axis of the rotation. We extend our supercell to a 4 × 2 × 1 model to have a decent oblate approximation. Then, we conduct an angular perturbation of layers positioned on top of each other at separation distances ζ = 2, 4 and 6Å (ζ = 0 corresponds to the equilibrium state). Due to small entropic eects that arise from angular perturbations and to reduce the computational cost, we model water molecules implicitly using dielectric continuum theory. Fig. 1.(d) shows the normalized PMF of angular perturbation along with the predictions of eq. 3 for an oblate particle with radii of ( a,a,c), where 2c = 12.1Å is the FTF equilibrium p distance in Fig. 1.(a). We consider two cases of a = Req and 1.7Req with Req = S/π . As it is evident, GB potential captures the contact between particles. The PMF for angular perturbation does not show a signicant minimum and it rises up fast at θ ≥ θmax . From this, we conclude that the center-to-center distance plays the most crucial role in the strength of the angular interactions, especially for small and medium range rotations. What GB is not accurately predicting is the angle in which contact occurs. For the case where a =

Req , GB predicts a larger angle of contact than simulations. This is primarily due to the distinction between the shape of the platelets in the simulations and oblate approximation in GB potential. On the other hand, for a = 1.7Req , GB predicts the correct value of contact angle, however, we emphasize that the equivalent surface area is compromised and care should be taken to use such approximation for upscaling purposes.

Sliding Component of Free Energy In essence, the total intermolecular interaction is a synergistic eect of cohesive and sliding interactions. While the former is mainly associated with non-bond interactions and entropic eects, the later is more inuenced by the surface chemistry and morphological aspects. Thus, it becomes crucial to consider sliding interactions when studying the FEL of colloids with atomically uneven surfaces. In this section, we study the free energy barrier associated with the onset of relative transverse displacement (sliding) between stacked layers of Tober-

11

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

morite. The free energy calculated in this scheme is merely a function of relative coordinates of layers and gives information about materials resistance against sliding. Contrary to the normal case-study, restricting the motion of atoms in layers would adversely aect the calculation of energy barriers. 38 Henceforth, we allow atoms at surfaces to move freely during MD simulations. The FEP procedure used to calculate the sliding energy barriers is analogous to that of normal interactions except from smaller perturbation step size of 0.0625 Å and periodic boundary condition for layers, which ensures a smoother free energy variation. The morphology of the surface of colloids can be characterized by corrugation spacing. For crystalline C-S-H, there are two principal directions with periodic nature in parallel ( x1 ) and perpendicular ( x2 ) to silicate chains as shown in the inset of Fig. 2a. 38 The periodicity length along x1 and x2 directions is 7.3 Å and 11 Å respectively. Also, we note that the atomic structure of each periodic block is centrosymmetric. To reconstruct the sliding FEL within each periodic block, we rst calculate free energy in the two principle directions. As outlined before, the sliding takes more eect from structural motifs, and therefore it is reasonable to consider such interactions at dierent interlayer separation distances, i.e.,

ζ . Furthermore, we emphasize that sliding simulations at small separation distances require enforcing bridging sites of opposite layers to intersect each other, which can induce unrealistic bond breakage in our nonreactive force eld. Henceforth, we study the sliding at greater separation distances, i.e. ζ ≥ 0.75Å. The insets in Fig. 2a show PMFs along x1 and x2 at ζ = 0.75, 1.3, and 3.9 Å . We choose the transverse equilibrium point ( x1 = 0 and

x2 = 0) as the reference free energy, and we sequentially sum the free energy dierences with respect to it.As it is evident from Fig. 2a, the increase in the separation distance has a considerable impact on the PMF. Also, regardless of the direction and separation distance, the PMF follows the same trend, that is trigonometric-like behavior. The insets of Fig. 2.(a) also indicate that for greater separation distances, the sliding free energy barriers (SFEB) decreases.

12

ACS Paragon Plus Environment

Page 12 of 42

Page 13 of 42

5

x1 (Å)

0

x1

2 t

3 2 1

(a)

0 0

x2

xx11 xx22 Exponential Fit 2

4 ζ (Å)

2

6

4

ζ= 0.75 Å

4

Asrel (KJ/mol.Å2)

4 Asrel (KJ/mol.Å2) rel(KJ/mol.Å

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

3

ζ= 1.3 Å

2 1

ζ= 3.9 Å

0 4 3 2 1 00

ζ= 0.75 Å ζ= 1.3 Å

ζ= 3.9 Å

2

4

6

x2 (Å)

6

8

10

8

(b)

(c)

Figure 2: Sliding interactions between corrugated crystalline C-S-H layers. a) Training sliding free energy barriers at various center-to-center distances, parallel (x 1 ) and perpendicular (x2 ) to silicate chains. The graphical inset shows the direction of silicate chains with respect to the edges of the nite-sized layer. Inset plots compare calculated free energy proles along x1 and x2 directions at three center-to-center distances with the analytical transverse energy landscape model in eq. 13. Note that free energy proles are shifted to a baseline at zero. b) Schematics of the nite silica chains in which red and yellow particles respectively indicate oxygen and silicon atoms. A set of protruding silicate tetrahedron are designated by blue circles. c) The protruding sites are structural motifs that present energy barriers against transverse sliding. The intensity of colors represents the strength of free energy barrier. The maximum barrier takes place at positions where protruding sites of a layer, shown with blue circles in (b), coincide with the protruding sites of the other layer. A general sliding free energy barrier can be reconstructed along an arbitrary path using the free energy barrier along x1 and x2 directions.

13

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 42

As we seek to dene an analytical model for SFEB, we rst measure the dependency of free energy barrier to the separation distance ζ . Fig. 2a depicts the change in the maximum SFEB in x1 and x2 directions as a function of the separation distance. We note that the dierence between SFEB along x1 and x2 directions is negligible. From here we use a unied mathematical expression for both directions. The maximum SFEB shows an exponential dependence on ζ . A similar trend has been corroborated by the experimental study of colloidal systems by Feiler et al. 55 Furthermore, the reduction in the amplitude of stick-slip friction as a result of small normal displacements has been observed for various interfaces in experiments. 15 The current results show that sliding at the nanoscale has a soft contact character. Henceforth, we t SFEB results in Fig. 2.(a) with a decaying exponential function as follows,

ζ a7 Am rel (ζ) = a6 e −

(8)

where a6 is the height of SFEB at the zero separation and is related to the number of contact points between two layers and will be adjusted for layers with dierent surface roughness (e.g. crystalline C-S-H with higher C/S). Also, a7 is the characteristic contact length. The tting exercise yields a6 = 8.5 KJ/mol.Å2 and a7 = 1 Å. Eq.8 indicates that as colloids approach each other during the course of sliding, the eect of substrate corrugation on SFEB becomes more pronounced. To further generalize eq. 8 to encompass sliding interaction between two nite-sized particles with arbitrary shape, relative distance and spatial orientations, we propose the following form:

Agrel

1 (r 12 , ω12 ) = Sint

Z

g(S)Am rel (ζ)dS

(9)

S21

where layer 2 with surface area S2 is sliding on top of layer 1 with surface area S1 . S21 is the normal projection of the second layer's area on the rst layer's mid-plane. Sint is the

14

ACS Paragon Plus Environment

Page 15 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

overlapping area of S21 and S1 . Moreover, we dene g as the the overlap function,

g(S) =

   1, ∀ dS ∈ Sint

(10)

  0, ∀ dS 6∈ Sint For a special case of an oblate particle with radius a and S21 ∈ S1 , eq. 9 leads to the following rotation-dependent analytical form,

Agrel

2πaa6 a7 cot (θ) (h12 ez , θ) = I1 Sint



a sin (θ) a7

 e



h12 a7

(11)

where I1 is the modied Bessel function of rst kind and θ is the out-of-plane angle of rotation illustrated in the inset of Fig. 1d. Eq. 11 has an indeterminate form at θ = 0, which after limit evaluation yields eq. 8. The next step toward deriving a comprehensive analytical model is to describe the SFEB as a function of any point in the x1 -x2 plane. Here, we introduce a simple yet realistic notion to explain the SFEB mechanism. We consider the source of energy barriers to be located at positions where protruding sites on the colloids surface coincide in the course of sliding. We designate these positions, as a source of repulsion in sliding, with blue circles in Fig. 2.(b). Then, we take the SFEB of any arbitrary point P to be proportional to its distance from the closest SFEB's prole d1 and d2 respectively along x1 and x2 directions, see the schematic in Fig. 2.(c). The color intensity in Fig. 2.(c) represents the magnitude of SFEB. Here, we propose that the move in any arbitrary path in x1 − x2 plane can be dened as a function of SFEB along principal axes and the distance from the closest periodic length. We display the border of the periodic block on Tobermorite surface with a solid rectangle. Furthermore, we show the primitive cell inside the periodic block via a dashed rectangle in Fig. 2.(c). Here, we performed FEP only inside the primitive cell and construct the entire SFEB in the periodic block using symmetry operations. Inspired by the model presented in Fig. 2, we propose an analytical expression for the trend of SFEB

15

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 42

in x1 − x2 plane as follow:

 Atrel (r 12 ) ∝ 

1 − cos



2

2πx1 Lp1

  

1 + cos



2πd2 Lp2



2



1 − cos

+



2πx2 Lp2

2

 

1 + cos



2



2πd1 Lp1

  (12)

where d1 and d2 are distances from closest rectangle sides, x1 = r 12 .ex1 and x2 = r 12 .ex2 .

Lp1 and Lp2 are the periodic lengths along x1 and x2 . See supporting information for further details about eq. 12. In the light of the above discussions, we construct SFEB as a multiplication of the maximum intensity of SFEB and its trend, and consequently propose the following analytical expression for sliding free energy prole,

Asrel (r 12 , ω12 ) = Agrel (r 12 , ω12 ) Atrel (r 12 )

(13)

Interestingly, simultaneous Atomic Force Microscopy (AFM) and Scanning Tunneling Microscopy (STM) measurements reported by Ternes et al. 56 determines the energy barrier during transverse movement of adsorbed atoms on various substrates. They found the decaying exponential dependence of energy barrier on separation distance. Moreover, they reported a trigonometric-like trend for their sliding energy barrier, which is consistent with the current FEP calculations. Obviously, for dierent crystalline structures, eq. 12 can be tailored and replaced in eq. 13 without loss of generality. we assess the transferability of eq. 13 by comparing FEP calculations along the diagonal path of the periodic rectangle for zero tilt angle (θ = 0) at ζ = 1.3Å (See Supporting Information for comparison of the predicted PMF by eq. 13 and FEP along the complex diagonal path).

Global free energy landscape and Its Application So far, we have developed an analytical expression for variation of free energy in normal and sliding contributions. The amount of work needed to displace a layer from its equilibrium

16

ACS Paragon Plus Environment

Page 17 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

state to a specic separation distance should be equal to the height of corresponding normal free energy barrier. On the other hand, when layers are farther from each other, less work is required for sliding. Therefore, the global FEL is aected by the interplay between two mechanisms, namely normal and sliding free energy barriers. Since the normal PMF would have a dierent trend for various value of x1 and x2 , we consider that eq. 6 is part of a bigger picture and needs to be extended so that it reects the inuence of initial position in

x1 -x2 plane. In addition, in the previous section, we took the origin of SFEB equal to zero for various separation distances, we stress here that in the global sense, the initial value of the free energy varies with the separation distance. Therefore, we can reconstruct the global FEL, Aglb rel , by superposing the normal and sliding free energy proles, n s Aglb rel (r 12 , ω12 ) = Arel (r 12 , ω12 ) + Arel (r 12 , ω12 ) + c

(14)

where we take c = 4.95 KJ/mol.Å2 such that at the equilibrium state global free energy barrier pertains to zero and accordingly any point in the phase space takes a relative value with respect to this reference point. Eq. 14 packs full atomistic FEP calculations into a simple analytical function. We can further explore applications of such reduced-order free energy barrier function to understand mechanisms of deformation at the coarse-grained level. In what follows, we restrict our attention to interactions between two parallel layers,

ω12 = (0, 0, 0), as it is energetically the most favorable agglomeration scenario due to the great FTF energy barrier. The question here is what is the mechanism for sliding and whether it involves normal jumps in the absence of external load? It would be of particular interest to see whether the hypothesis of occurrence of normal jump during possible sliding is valid or not. 38 To this end, we depict the global free energy landscape as a function of interlayer separation at the saddle point located at x1 = 1.85Å and x2 = 2.75Å . Fig. 3 shows the competition between increasing free energy due to partial rupture in the normal direction and decreasing SFEB as a result of the further distance between protruding sites on the surfaces.

17

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 42

The minimum at ζ = 0.6Å is the most probable interlayer separation that sliding occurs. This behavior meets our expectation of normal jumps during sliding as it is energetically more favorable. These results indicate that during transverse displacement of a layer, there are extra normal forces induced into the system. In the absence of external resistance, a normal displacement is expected to satisfy the equilibrium condition, i.e. the sum of forces become zero. Alternatively, if normal displacement of layers is constrained (displacement control sliding), which is the more probable case for multilayer material systems, normal stress arises. The inset of Fig. 3 illustrates 3-dimensional FEL in x1 − x2 plane at optimum separation distance for sliding. The presence of saddle point in between local minimums represents the energy barrier during the transition between neighboring basins. The transition between two neighboring free energy basins comes about through saddle points. Since the energy barriers are much greater than thermal uctuations, this transition is rather a rare event. The required time for such transition is proportional to the height of eective energy barrier in the presence of external force, F, and can be described via Bell model, 57

 t = t0 exp

E0 − γF kB T



(15)

wherein its original form, t0 is the reciprocal of a natural frequency of atomic vibration in the solid, E0 is the energy barrier with respect to the minimum well depth, and γ is a constant that takes into account the structure of the solid and can be determined empirically. In eq. 15, the height of energy barrier and its position is assumed to be xed. γ is taken as a distance between the position of the minimum ( x∪ ) and energy barrier ( x∩ ). This distance in eq. 15 is xed and equal to ∆x0 such that at F = E0 /∆x0 , the eective energy barrier disappears. Before advancing further in the quantitative study of the transition time, we examine the eective energy barrier landscape. It is well-known that the energy landscape twists upon applying external stress. 58 The governing equation of eective FEL

18

ACS Paragon Plus Environment

Page 19 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

Figure 3: The variation in the height of the saddle point versus interlayer separation distance. The competition between the rise in normal free energy along ζ and fall of sliding free energy at saddle point in x1 −x2 plane results in an oscillatory behavior. The minimum of the graph shows the optimum interlayer separation distance for sliding concerning global free energy.

19

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 42

in the presence of external stress eld is,

Grel (r 12 , ω12 ) =

Aglb rel

1 (r 12 , ω12 ) − S

Z Z σ : d Ω

(16)



where Grel is the so-called Gibbs free energy, σ and  are respectively the external stress and strain tensors in the domain Ω. Here, our stress tensor includes τx1 ζ and τx2 ζ . Also, we take into account the work of normal force to be exactly equal to the normal energy barrier at any given separation distance. Fig. 4 shows the 2D color map of the energy barrier landscape in x1 − x2 plane within a periodic block explained in Fig. 2.b,c. As it is evident, increasing forces tilts the free energy landscape. In the event of applying an external stress, the energy barrier bends in the direction of stress. Looking at the energy barrier colormap, it is interesting to note that an increase in separation distance does not guarantee the reduction in sliding free energy barrier, for instance, compare the second and third columns in Fig. 4, which further conrms the results shown in Fig. 3. It is known that the magnitude of the external stress and the shape of FEL alters not only the height but also the position of the saddle point. 59 Therefore, such alterations should be considered when studying transition between two adjacent free energy basins. Back to our discussion of transition time between states, we take into consideration the fact that external force eld aects both the magnitude of the eective energy barrier and the relative position of the energy barrier. Therefore, we rewrite the transition time as a function of force via Bell's equation at an updated path toward the energy barrier which buries in itself the new position of x∩ and reduced height of eective energy barrier as follow,

t∗ (τ ) = exp {A∗b (1 − W ∗ /A∗b )} where, t∗ = t/t0 , A∗b =

Ab (x∩ (τ )) , kB T

and W ∗ =

τ.(x∩ (τ )−x∪ (τ )) . 2kB T

(17)

For a colloidal system with nano-

metric sizes, t0 would have a dierent interpretation than atoms. Moreover, to investigate the behavior of the transition time, we take kB T = 1 KJ/mol.A2 to facilitate the calcula20

ACS Paragon Plus Environment

Page 21 of 42

ζ =0Å

ζ = 0.75 Å

ζ = 1.3 Å

(GPa) τx1ζ = 0 τx2ζ = 0 8 5

τx1ζ = 0.8 τx2ζ = 0

2 -1 -4

τx1ζ = 0 τx2ζ = 0.8

Grel (KJ/mol.Å2)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

-7

x1 τx1ζ = 0.8 τx2ζ = 0.8 x2 Figure 4: Tow-dimensional colormap of the Gibbs free energy landscape in a periodic block as a function of the interlayer spacing and external shear stresses applied along x 1 and x2 directions. The application of external force bends the energy landscape along the direction of the applied shear stress.

21

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tions. Fig. 5.a illustrates the distribution of dimensionless transition time in the presence of external stress. At weaker force elds the transition requires a longer time to take place, whereas at stronger force elds t∗ falls below the value of 1 which implies that the eective energy barrier has completely vanished and transition occurs almost immediately. The access to the global FEL also allows us to investigate the interplay between the normal and shear stress components at the onset of sliding. This is conventionally studied in the context of macroscopic friction. Here, we are concerned with friction in colloidal systems. In what follows, we assume that thermal uctuation contribution in sliding is negligible. Also, we study onset of frictional forces by moving across FEL and therefore no bond breakage or energy dissipation is being considered as our force-eld is nonreactive. To study the onset of sliding at a prespecied normal stress, we nd the shear stress where the Hessian matrix of Gibbs free energy possesses at least one negative eigenvalue. This point corresponds to states in which the colloidal assembly becomes unstable, i.e. yielding stress. We subsequently study the relation between the applied normal stress ( σn ) and corresponding shear stress prior to the failure ( τcr ). The applied external normal stress during sliding is counterbalanced by two internal stress components. The rst component comes from the normal interaction introduced via Anrel and the second part is due to the induced repulsion as the protruding corrugations are forced against each other during the transverse displacement. In order to dierentiate between these two reactions, we introduce a coefcient to the normal interaction αAnrel and adjust the contribution weight by changing α. Therefore, the externally applied stress is going to be carried as ασnn + σns = σn , where σnn and σns are respectively, induced normal stress due to rst and second above-mentioned interaction components. The inset of Fig. 5b schematically illustrates these contributions. As it is evident from Fig. 5b, in the absence of normal interaction ( α = 0), τcr − σn relationship is reminiscent to that of non-cohesive granular materials. In this case, τcr − σn relation is linear and the cohesion, i.e. the shear stress at zero conning pressure, is zero. As we increase α, the cohesion increases and the friction angle at σn = 0 decreases simultaneously.

22

ACS Paragon Plus Environment

Page 22 of 42

Page 23 of 42

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Applied Materials & Interfaces

More importantly, the sublinear trend in τcr − σn relation emerges as we intensify the normal contribution. The intersection of τcr − σn lines is the same for all values of α and corresponds to σnn = 0. At this point, contributions of attractive and repulsive terms in σnn cancel each other out. Prior to this point ( |σn | ≤ |σn0 |), the attractive forces are dominant, while beyond it (|σn | ≥ |σn0 |) repulsive forces due to the contact play the major role in tuning the τcr − σn relation. To evaluate the contribution of adhesion on τcr , we adjust the strength of attractive forces (mainly van der Waals interactions) in the normal contribution. To achieve this, we specify a parameter β as in for β = 0 there is no cohesion in the material and for β = 1 the attraction ∗ ) as a function is fully restored. Fig. 5c depicts the normalized critical shear stress ( τcr

of β in the absence of external normal stress. For small values of β , the repulsion due to sliding results in large separation distance and therefore the sliding takes place without any ∗ ∗ signicant resistance ( τcr ≈ 0). As the van der Waals attraction increases, τcr becomes more

signicant. This result further supports the contribution of attractive forces in the trend observed in Fig. 5b for ( |σn | ≤ |σn0 |) region. The results in Fig. 5b,c mark the contribution of attractive forces at the onset of sliding for the same corrugation morphology. In order to investigate the eect of corrugation spacing, we modulate the periodicity length along x1 direction and study its eect on the critical shear stress. We dene variable corrugated spacing as Lp = n × Lp1 and adjust it by varying n parameter. As it is evident from Fig. 5d, for smaller values of n, the surface is rough and the critical shear stress tends to larger values. This implies the presence of a interlocking phenomena during the sliding process. We realize that the increased critical shear stress for smaller spacings is based on the assumption that the protruding sites on surfaces can penetrate each other, which might become invalid at the limit of very small spacings. In reality, the applied stresses should distort the surface through an irreversible thermodynamics process to overcome such interlocking situation. Such distortion results in a dissipation of frictional heat. 60 For greater corrugation spacings, the surface is smoother and less force is

23

ACS Paragon Plus Environment

ACS Applied Materials & Interfaces

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 42

required for the onset of sliding. The increase in the average and uctuations of the friction stress as a result of larger spacing has also been reported for atomically stepped surfaces. 61 Overall, these results indicate that the sublinear trend stems from the presence of cohesive interaction 14 and the onset of sliding is highly correlated to the morphological aspects of the interface.

Comparison of Mechanical Properties with Experiments To further validate our results, we calculate mechanical properties and compare them with available experimental data from the literature in Table 1. We calculate elastic modulus as:

 Ezz = d0

∂ 2 hAnrel i ∂ζ 2



(18) d0

where d0 is the FTF equilibrium distance and is 12.1 Å in our simulation. We nd Ezz =

77.6GP a, which is in a good agreement with experimental results considering the perfect atomic structure studied here. Also, we compute surface free energy as:

 γs =

hAnrel i − 2



(19) d0

The value we work out for surface free energy of Tobermorite is 0.42 (J/m 2 ) that is close to the experimental report of 0.32-0.4 62,63 (J/m2 ) and lower than previous MD simulations of 1.15 (J/m2 ) using topological constraint theory. 64,65 Also, we compute the the cohesive pressure (Pco ) as the maximum force per area in FTF interaction (see mean force graph in Supporting Information). We calculate mean force at any separation distance as follow,

hF ii =

h∆Aii,i+1 ∆λ

(20)

where ∆λ is the perturbation step size. From FEP results we nd P co =5 GPa which is

24

ACS Paragon Plus Environment

Page 25 of 42

t* (arb. Units)

2

σn = ασnn +σns

σn = σn n

τcr

20 16

1

12

τcr (GPa)

τx2ζ (GPa)

24

8 4