Elastic constants and homogenized properties of two monoclinic structures (gypsum and tobermorite) were investigated by first-principles method. The gypsum (chemical formula of CaSO4•2H2O) is an evaporite mineral and a kind of hydration product of anhydrite. Besides, the 11 Å tobermorite model (chemical formula: Ca4Si6O14(OH)4·2H2O) as an initial configuration of C-S-H structure is commonly used. Elastic constants are calculated based on density functional theory (DFT), which can also contribute to provide information for investigating the stability, stiffness, brittleness, ductility, and anisotropy of gypsum and tobermorite polycrystals. In addition, based on elastic constants (13 independent constants) of the monoclinic gypsum crystal, the elastic properties of polycrystals are obtained. The bulk modulus B, shear modulus G, Young’s modulus E, and Poisson’s ration ν are derived. Therefore, it is fairly meaningful to study the elastic constants to understand the physical, chemical, and mechanical properties of two monoclinic structures. Elastic constants can be used as the measure criterion of the resistance of a crystal to an externally applied stress. The calculated parameters are all in excellent agreement with reference.
- DFT calculation
- single crystal
- nano scale
- elastic constants
- homogenized moduli
The density functional theory (DFT) is commonly used to study the crystal structure, lattice energy, the equation of state, the electronic bandgap, and vibration spectra properties . Based on the kinetic energy density functional of Thomas  and the exchange-correlation effects of Dirac , DFT has been greatly developed by Kohn and Sham (KS) , who have established the fundamental approximation theorem on the functional status to describe real systems by electronic structure calculations. The eigenvalues of KS equations have no physical meaning, and the ionization energy is in the opposite state direction . Moreover, one proposed approach is to introduce the eigenstates to calculate multi-body (many-body calculation) on the basis of Monte Carlo calculations  and perturbation theory . The calculation of elastic constants is preceded by full geometry optimization and the stress tensor calculation of a number of distorted structures at the atomic scale. Polycrystalline structure constituted by a single crystal structure contains a variety of information (e.g., orientation) and the properties of a single crystal, such as anisotropy. Within the mechanics of typical crystals structures, the transition from the micro- to the meso-scale (homogenization) and vice versa (localization) can be estimated. Homogenization is an idealized description of a statistical distribution inside the actual heterogeneous material. Once the continuity model is admitted, the concept of homogeneity is deduced from it . For quasi brittle materials, Zhu et al.  have formulated the anisotropic model in the framework of Eshelby-based homogenization methods. X-rays diffraction measurement is one of the stress assay test methods in physics field, of which the stress is actually determined by the strain . Diffraction-based stress analysis depends critically on the use of the correct diffraction elastic constants . X-ray method to test the material stress and to obtain elastic constants  is commonly based on the Reuss model . Elasticity of single crystal and mechanical properties of polycrystalline material have been closely integrated. Various calculations methods are compared to determine homogenized moduli of the polycrystalline material composed of a single crystal, for example, the certain stress of Reuss model  and the certain strain of Voigt model .
DFT as a first-principles theory and a solid band theory in quantum mechanics has own a great success in linking physical properties and molecular structure, the calculation with exact accuracy but for low computational efficiency for macromolecular structure, which can be used to calculate elastic constants of anisotropic crystals, the monoclinic gypsum, and tobermorite crystals, for example. The chemical formula of gypsum is CaSO4•2H2O, which is an evaporite mineral and a kind of hydration product of anhydrite (chemical formula: CaSO4). Moreover, the 11 Å tobermorite model (chemical formula: Ca4Si6O14(OH)4·2H2O) as an initial configuration of C-S-H structure is commonly used. Since Young’s modulus parameters of gypsum and C-S-H are important to the multi-scale model , elastic constants of the gypsum crystal are investigated. The crystal is monoclinic, with 13 independent constants. For the homogenization of elastic deformation, especially for polycrystalline structures, the traditional Reuss-Voigt-Hill method is used to calculate the elastic moduli of monoclinic structures. Based on the ab initio plane-wave pseudopotential density functional theory method mentioned earlier, we focus on the monoclinic crystals to estimate their homogenized elastic moduli.
2. Theoretical calculation by density functional theory (DFT)
Despite the above advantages of DFT, however, the resolution of a system by Kohn-Sham equations involves difficulties due to an infinite number of electrons. These electrons maybe changed under an effective potential generated by an infinite number of cores or ions.
2.1. Equation of the theoretical approximate solution
From a microscopic point of view, Schrödinger equation describing a periodic crystal system composed of atomic nuclei n in mutual interaction and electron spin σi is positioned
Hamiltonian, in simple cases, consists of five terms: the kinetic energy of the electrons and nuclei, and the various interactions between them.
The possible analytical representation and resolution of such a problem become a difficult task due to the limited memory of the computer tools. However, it is possible to reformulate the problem using appropriate theorems and approximations.
The fundamental principle approaches of mean field theory, in particular the DFT, are that any properties of an interacting particle system can be considered as a functional density in the ground state of the system n0(r). Besides, the scalar function of the position n0(r) essentially determines the wave functions of the system at the ground state and the excited states. Electronic and mechanical properties of a periodic crystal refer to solid state physics, quantum mechanics, and crystallography.
The crystalline ion movement of the electron is as
According to the Born-Oppenheimer or adiabatic approximation , the dynamics of the system (electrons and nuclei) is described. The electrons are assumed to react instantly to ionic motion. In electronic coordinates, the nucleus positions are considered as immobile external parameters.
where the last term of the Hamiltonian is constant and has been introduced in order to preserve the neutrality of the system and avoid the divergence of the eigenvalues. Clean the ground state of the system for fixed nuclear positions, total energy is given by the formula:
This energy has a surface in the space coordinates that is said to be ionic Born-Oppenheimer surface. The ions move according to the effective potential energy, including Coulomb repulsion and the anchoring effect of the electron, which are as follows:
The dissociation degrees of freedom of electrons from those of nucleons, obtained through the adiabatic approximation, are very important, because if the electrons must be treated by quantum mechanics, degrees of freedom of ions in most cases are processed in a conventional manner.
This theorem/approach of Hohenberg and Kohn tries to make an exact DFT theory for many-body systems. This formulation applies to any system of mutually interacting particles in an external potential
The total energy of the ground state of a system for interacting electrons is functional (unknown) of the single electron density
As a result, the density n0(r) minimizing the energy associated with the Hamiltonian (9) is obtained and used to evaluate the energy of the ground state of the system.
The principle established in the second theorem of Hohenberg and Kohn specifies that the density that minimizes the energy is the energy of the ground state
Because the ground state is concerned, it is possible to replace the wave system function by the electron charge density, which therefore becomes the fundamental quantity of the problem. In principle, the problem boils down to minimize the total energy of the system in accordance with the variations in the density governed by the constraint on the number of particles
2.2. The approximation approach of Kohn-Sham
The approach of Kohn-Sham system substitutes the interacting particles, which obeys the Hamiltonian in Eq. (3), by a less complex system easily solved. This approach assumes that the density in the ground state of the system is equal to that in some systems composed of non-interacting particles. This involves independent particle equations for the non-interacting system, gathering all the terms complicated and difficult to assess, in a functional exchange-correlation .
T is the kinetic energy of a system of particles (electrons) independently (non-interacting) embedded in an effective potential which is no other than the real system,
The Hartree energy or energy of interaction is associated with the Coulomb interaction of the self-defined electron density.
Solving the auxiliary Kohn and Sham system for the ground state can be seen as a minimization problem while respecting the density n(r). Apart from orbital function TS, all other terms depend on the density. Therefore, it is possible to vary the functions of the wave and to derive the variational equation:
With the constraint of orthonormalization , this implies the form of Kohn-Sham for Schrödinger equations:
represents the eigenvalues, and is the effective Hamiltonian H,
Eqs. (16)–(18) are known equations of Kohn-Sham, the density n(r) and the resulting total energy EKS. These equations are independent of any approximation on the functional EXC(n), resolution provides the exact values of the density and the energy of the ground state of the interacting system, provided that EXC(n) is exactly known. The latter can be described in terms of Hohenberg Kohn function in Eq. (8)
or more precisely,
This energy is related to potential exchange-correlation .
For the exchange-correlation functional, the only ambiguity in the approach of Kohn and Sham (KS) is the exchange-correlation term. It is subject to functional approximations of local or near local order of density that said energy EKS can be written as
where εxc([n], r) is the exchange-correlation energy per electron at point r, it depends on n(r) in the vicinity of r. These approximations have made enormous progress in the field.
The approximation of the local density (LDA)
The use of the local density approximation (LDA) in which the exchange-correlation energy is another integral over all space, assuming that εxc([n], r) is the exchange-correlation energy per particle of a homogeneous electron gas of density nE22
This approximation has been particularly checked to deal with non-homogeneous systems.
The generalized gradient approximation (GGA)
The generalized gradient approximation (GGA) involves the local density approximation providing a substantial improvement and better adaptation to the systems. This approximation is equal to the exchange-correlation term only as a function of the density. A first approach (GEA) was introduced by Kohn and Sham then used by the authors of Herman et al. .
This notion of GGA is the choice of functions, which allows us a better adaptation to wide variations so as to maintain the desired properties. The energy is written in its general form :E23
2.3. Parameters of Bloch theorem and Brillouin zone
Different states of the Schrödinger equation for an independent particle in a system. By Kohn and Sham equations, responding to eigenvalue equation is as:
where the electrons are immersed in an effective potential
The effective potential has the periodicity of the crystal and may be expressed using Fourier series, in a periodic system:
Gm is the reciprocal lattice vector:
where is the volume of the original mesh.
As the translational symmetry, it is that states are orthogonal and conditioned by the limits of the crystal (infinite volume). In this case, the Eigen functions of KS are governed by the Bloch theorem: they have two quantum numbers: the wave vector k in the Brillouin zone (BZ) and the band index i, and this can be expressed by a product of a plane wave exp.(ik, r) and a periodic function:
where R is the vector of direct space defined by ai with and Ni is the number of primitive cells in each direction (in the case of perfect crystal).
Solving Eq. (24) is equivalent to increase the periodic function , in a database-dependent functions points :
where is the wave function developed in a space of infinite dimensions; this means that j should be in principle infinite. But, in practice, we work with a limited set of basic functions, which imply that the description of will approximate. That the selected database simply solves the system:
where each point is a set of k eigenstates, the label having i = 1, 2, … obtained by diagonalization of the Hamiltonian in Eq. (29).
It is necessary to integrate the points k in the Brillouin zone. For a function fi(k) where i defines the band index, the average value is
is the cell volume of the original mesh in the real space and of the cell volume of the Brillouin zone are determined using a sampling points k. Several election procedures exist for these points. Particularly those of Baldereschi , Chadi and Kohen , and Monkhorst and Pack  are the most frequently used.
3. Elastic constants and homogenized moduli of monoclinic structure
According to the crystal theory , any crystal lattice system contains six independent variables, namely the cell side length a, b, and c; unit cell angle α, β, and γ. Generally, the crystal under a certain deformation, temperature, and pressure can be described by the corresponding six-dimensional deformation tensor. The temperature and pressure will cause cell-deformed configuration tensor as
where and separately represent the cell configuration tensor and the deformation tensor at temperature T (K) under the pressure P (GPa).
3.1. Calculation of elastic constants for single crystal structure
A multi-particle electronic structure satisfies the Schrödinger equation. As in , Kohn-Sham equation as an approximation to simplify Schrödinger equation is described. For crystal composed by vibrator with the vibration frequency wi, the total Helmholtz free energy is
Helmholtz free energy can be calculated for all the thermodynamic quantities. DFT-QHA (quasi-harmonic approximation) is a precise calculation method to calculate thermodynamic properties of solid materials elastic constants and Debye temperature with the accurate predictions.
According to the theory of elasticity, under the isothermal strain, the elastic modulus of Helmholtz free energy can be described by the form of the Taylor expansion, of which the coefficients of the polynomial are the elastic coefficient:
where ,, and are the coefficients of Lagrange deformation tensor, is the isothermal first-order elastic coefficients, and is the Helmholtz free energy.
The components of the stress tensor can be extracted by , after the applied strain, the total energy variation of the system can be expressed as
The second-order elastic coefficients can be obtained by the coefficient of the second-order Taylor expansion of Helmholtz free energy with the strain,
Here, strain and thermodynamics deformation are symmetric. There is only six independent deformation tensor in the nine-dimensional deformation tensor. LCEC is a second-order linear combination of independent elastic coefficients corresponding to Helmholtz free energy coefficient under some deformation mode [30, 31]. For all directions under monoclinic crystals, if a strain is added, the corresponding simultaneous equations can be solved to determine all elastic coefficients.
3.2. The energy-volume relationship of the monoclinic crystal
Deformation tensors to calculate independent Cij constants of monoclinic crystal are listed in Table 1.
|Deformation tensor||ΔE-V relation of LCEC||LCEC|
|c11+ c22 + 2c12|
|c22 + c33 + 2c23|
|c11 + c33 + 2c13|
|c44 + c55 + 2c45|
|c11 + c66 + 2c16|
|c22 + c66 + 2c26|
|c33 + c66 + 2c36|
|c11 + c22 + c33 + 2c12 + 2c13 + 2c23|
|c11 + c22 − 2c12|
|c11 + c33 − 2c13|
|c22 + c33 − 2c23|
|c11 + 4c55 + 4c15|
|c11 + 4c55 − 4c15|
|c22 + 4c55 + 4c25|
|c33 + 4c55 + 4c35|
For monoclinic crystal, elastic constants include C11, C22, C33, Cl2, C13, C23, C44, C55, C66 Cl5, C25, C35, and C46; the strain energy-volume relation and elastic moduli of monoclinic symmetry based on E-V method can be obtained. The calculated E-δ points are fitted to second-order polynomials E(V, δ). For all strains, different strain forms δ are taken to calculate the total energies E for the strained crystal structure. By applying a series of δ strain amplitude, the independent elastic constants of monoclinic crystal by these simultaneous ΔE-δ equations can be obtained.
3.3. Homogenization of monoclinic polycrystals by RVH estimation
Stress-strain relation in an orthotropic monoclinic crystal can be defined by the independent elastic stiffness parameters :
where σ represents the normal stress and shear stress in each direction (unit: nN/nm2); ε and γ are the normal strain and shear strain in each direction, respectively.
The homogenized elastic properties of polycrystals can be calculated, of which elastic moduli and Poisson’s ratio can be obtained by calculating Voigt and Reuss bounds and averaging term as 
For monoclinic crystal structure, elastic constants include C11, C22, C33, Cl2, C13, C23, C44, C55, C66 Cl5, C25, C35, and C46. The criteria for mechanical stability are given by Wu :
Young’s modulus and Poisson’s ratio can be rewritten based on the Voigt-Reuss-Hill approximation . In terms of the Voigt-Reuss-Hill approximations , MH = (1/2)(MR + MV), M refers to B or G. Thus, Young’s modulus E and Possion’s ratio μ are obtained as
Then, Voigt-Reuss-Hill average  will be determined, and Young’s modulus can be calculated.
4. Modeling and homogenized elastic moduli of gypsum structure
4.1. Nanoscale modeling of monoclinic crystals
4.1.1. Nanoscale modeling of monoclinic gypsum crystal
The gypsum morphology is monoclinic, and the initial lattice is as a = 5.677Å, b = 15.207Å, c = 6.528Å, α = β = 90°, and γ = 118.49°, its structure is monoclinic with space group I 2/a .
In Figure 1, the gypsum crystal can be summarized as follows: (1) the two hydrogen atoms of water molecules formed weak hydrogen bonds with the O atoms of Ca and S polyhedra; (2) a stacking sequence of CaO8 and SO4 chains in the (010) plane alternates with water layers along the b-axis; and (3) in (010) plane, the sulfate tetrahedra and CaO8 polyhedra alternate to form edge-sharing chains along  and zigzag chains along  direction  (Table 2).
|Atom||x||y||z||Occupancy rate||Uiso or Ueq|
4.1.2. Nanoscale modeling of monoclinic 11 Å tobermorite crystal
Hamid model  as the 11 Å tobermorite (formula: Ca4Si6O14(OH)4·2H2O) as an initial configuration is commonly used. The morphology is monoclinic, and the initial lattice is : a = 6.69 Å, b = 7.39 Å, c = 22.779 Å, α = β = 90°, and γ = 123.49°, space group P21. Modeling of 11 Å tobermorite is shown in Figure 2.
In Figure 2(a), the 11 Å tobermorite crystal can be summarized as follows: (1) the structure is basically a layered structure. (2) The central part is a Ca-O sheet (with an empirical formula: CaO2, of which the oxygen in CaO2 also includes that of the silicate tetrahedron part). (3) Silicate chains envelope the Ca-O sheet on both sides. (4) Ca2+ and H2O are filled between individual layers to balance the charges. The infinite layers of calcium polyhedra are parallel to (001), with tetrahedral chains of wollastonite-type along b and the composite layers stacked along c and connected through the formation of double tetrahedral chains . Atomic coordinates and displacement parameters are seen in Table 3.
|Atomic species||X||Y||Z||Occupancy rate||Uiso or Ueq||Atomic species||X||Y||Z||Occupancy rate||Uiso or Ueq|
4.2. Initial conditions and elastic constants of monoclinic crystals
4.2.1. Initial conditions and elastic constants of gypsum
The initial conditions are as follows: the pressure region of 0–1 GPa is used. Besides, a plane-wave basis set and ultrasoft pseudopotentials using GGA are used with a plane-wave cutoff energy of 400 eV. Brillouin zone is 6 × 6 × 4. Self-consistent convergence of the total energy per atom is chosen as 10−4 eV. Elastic constants of monoclinic gypsum crystal under 0–1.0 GPa are shown in Figure 3.
From Figure 3, elastic constants at 0 GPa are given as c11 = 82.464 GPa, c12 = 34.751 GPa, c13 = 33.643 GPa, c15 = −1.987 GPa, c22 = 63.046 GPa, c23 = 34.920 GPa, c25 = −8.071 GPa, c33 = 57.549 GPa, c35 = −3.054 GPa, c44 = 20.863 GPa, c46 = −4.688 GPa, c55 = 28.062 GPa, and c66 = 28.556 GPa. It is found that the oxygen atom of the water molecule did not change its position or occupancy under pressure conditions. A simple pressure increase at an ambient temperature cannot induce dehydration because of the unchange of water molecular in the gypsum structure within pressure range .
Elastic constants of gypsum crystal model based on DFT are calculated, and parameters are detailed in Table 4.
4.2.2. Initial conditions and elastic constants of tobermorite
Initial conditions of tobermorite are quite the same with that of gypsum crystal. Elastic constants of 11 Å tobermorite crystal under 0–1.0 GPa are shown in Figure 4. Elastic constants are shown in Table 5.
A comparisonal results of Shahsavari  are provided. Elastic constants at 0 GPa are as follows: c11 = 106.63 GPa, c12 = 50.37 GPa, c13 = 41.09 GPa, c15 = −3.50 GPa, c22 = 131.67 GPa, c23 = 22.78 GPa, c25 = −0.78 GPa, c33 = 71.45 GPa, c35 = −0.83 GPa, c44 = 26.03 GPa, c46 = −0.02 GPa, c55 = 27.61 GPa, and c66 = 45.26 GPa. Thus, elastic modulus can be homogenized to compare with the results of LD C-S-H phase in nano-indentation test by Vandamme and Ulm .
4.3. Homogenized elastic moduli of typical monoclinic structures
4.3.1. Elastic modulus of monoclinic gypsum structure
Based on elastic constants, the elastic moduli of gypsum at 0 GPa are verified and averaged in Figure 5.
As gypsum shows anisotropic compressibility along three crystallographic axes with b > c > a below 5 GPa , the pressure region of 0–1.0 GPa is used to verify whether the performance of model under low pressure is stable. Mechanical moduli of gypsum polycrystalline are listed in Table 6.
|Pressure (GPa)||Gv (GPa)||Bv (GPa)||Gr (GPa)||Br (GPa)||B (GPa)||G (GPa)||E (GPa)||μ|
As an acoustic method  and mechanical properties  have been investigated, according to elastic constants of gypsum crystal , elastic moduli by experiment can be calculated, as shown in Table 6. Elastic moduli are as follows: Gv = 22.146 GPa, Gr = 19.705 GPa, Bv = 45.521 GPa, Br = 43.822 GPa, B = 44.672 GPa, G = 20.926 GPa, E = 54.299 GPa, and μ = 0.2974. These results are close to the plane-strain value of Young’s modulus by reference  E = 50 GPa, μ = 0.45. By comparison of gypsum crystal and CH crystal, axial moduli of gypsum in x, y, and z directions are 57.75, 37.22, and 34.91 GPa, while axial moduli of Ca(OH)2 in x, y, and z directions are 93.75, 93.75, and 42.39 GPa, showing that gypsum crystal is much less anisotropic than hydrogen-bonded layered Ca(OH)2 structure .
4.3.2. Elastic modulus of monoclinic tobermorite structure
Elastic moduli at 0 GPa are verified and averaged as Gv = 32.815 GPa, Bv = 59.803 GPa, Gr = 29.908 GPa, Br = 54.276 GPa, E = 79.512 GPa, and μ = 0.268. Young’s modulus is about 79.512 GPa by Reuss-Voigt-Hill estimation, which is close to the simulation result of 89 GPa  by Pellenq and result of 78.939 GPa  by Shahsavari. Mechanical moduli by different methods are listed in Table 7.
|Pressure (GPa)||Bv (GPa)||Br (GPa)||Gv (GPa)||Gr (GPa)||B (GPa)||G (GPa)||E (GPa)||μ|
However, these values considering the ordered Si-chain at a long range are far away from the nano-indentation experiment performed on the C-S-H phase [40, 46]. It confirms the absence of order at a long range in this phase and that the up-scaling to polycrystals cannot be done with the tobermorite model. Modeling of C-S-H structure with disordered Si chain should be fairly considered.
Elastic constants of gypsum and tobermorite structures under a certain pressure region are calculated by DFT method, which has a certain value for both application and reference. Results are as follows:
For monoclinic gypsum and tobermorite crystals, elastic coefficients are obtained in 0–1-GPa pressure range to verify the reliability of the model by comparing other literatures.
Elastic constants of gypsum single crystal at 0 GPa are given as follows: c11 = 82.464 GPa, c12 = 34.751 GPa, c13 = 33.643 GPa, c15 = −1.987 GPa, c22 = 63.046 GPa, c23 = 34.920 GPa, c25 = −8.071 GPa, c33 = 57.549 GPa, c35 = −3.054 GPa, c44 = 20.863 GPa, c46 = −4.688 GPa, c55 = 28.062 GPa, and c66 = 28.556 GPa.
Elastic constants of 11Å tobermorite single crystal at 0 GPa are as follows: c11 = 106.63 GPa, c12 = 50.37 GPa, c13 = 41.09 GPa, c15 = −3.50 GPa, c22 = 131.67 GPa, c23 = 22.78 GPa, c25 = −0.78 GPa, c33 = 71.45 GPa, c35 = −0.83 GPa, c44 = 26.03 GPa, c46 = −0.02 GPa, c55 = 27.61 GPa, and c66 = 45.26 GPa.
Young’s modulus of gypsum is about 54.299 GPa. Elastic moduli at 0 GPa are as follows: Gv = 22.146 GPa, Gr = 19.705 GPa, Bv = 45.521 GPa, Br = 43.822 GPa, E = 54.299 GPa, and μ = 0.297.
Young’s modulus of 11Å tobermorite is about 79.512 GPa. Elastic moduli at 0 GPa are as follows: Gv = 32.815 GPa, Bv = 59.803 GPa, Gr = 29.908 GPa, Br = 54.276 GPa, E = 79.512 GPa, and μ = 0.268.
Structural, elastic properties of monoclinic crystals are investigated, and Cij determination is given by DFT method. Reuss-Voigt-Hill estimation has been used for polycrystal structures and can be seen as an intermediate step in the homogenization of elastic properties.
The authors greatly acknowledge the financial support for this work provided by the China Scholarship Council (CSC) and the support of start-up foundation of Xi’an Shiyou University. Thanks to Qiufeng Wang for her proofreading.