Open access peer-reviewed chapter

Elastic Constants and Homogenized Moduli of Monoclinic Structures Based on Density Functional Theory

Written By

Jia Fu

Submitted: 14 June 2017 Reviewed: 08 November 2017 Published: 20 December 2017

DOI: 10.5772/intechopen.72301

From the Edited Volume

Density Functional Calculations - Recent Progresses of Theory and Application

Edited by Gang Yang

Chapter metrics overview

1,576 Chapter Downloads

View Full Metrics


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

1. Introduction

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 [1]. Based on the kinetic energy density functional of Thomas [2] and the exchange-correlation effects of Dirac [3], DFT has been greatly developed by Kohn and Sham (KS) [4], 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 [5]. Moreover, one proposed approach is to introduce the eigenstates to calculate multi-body (many-body calculation) on the basis of Monte Carlo calculations [6] and perturbation theory [7]. 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 [8]. For quasi brittle materials, Zhu et al. [9] 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 [10]. Diffraction-based stress analysis depends critically on the use of the correct diffraction elastic constants [11]. X-ray method to test the material stress and to obtain elastic constants [12] is commonly based on the Reuss model [13]. 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 [13] and the certain strain of Voigt model [14].

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 [15], 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

and assumes that the electron mobility (φ) does not depend on the speed nuclei but on their positions.

According to the Born-Oppenheimer or adiabatic approximation [16], 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

, where the Hamiltonian is written as.


DFT and its founding principle are summarized in two theorems, first introduced by Hohenberg and Kohn [17], which refer to the set of potential

and the density minimizing of Eq. (5).

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

. In this stage, the DFT can reformulate the problem rather than solve an uncertain functional FHK(n).

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 EXCn.


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 φiφj=δi,j, this implies the form of Kohn-Sham for Schrödinger equations:


εi represents the eigenvalues, and HKS 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 Vxc=Excnr.

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.

  1. The approximation of the local density (LDA)

    The use of the local density approximation (LDA) in which the exchange-correlation energy ExcLDAn 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 n


    The exchange term Exbomnr can be expressed analytically, while the correlation term is computed accurately using the Monte Carlo by Ceperley Alder [18] and then set in different shapes [19].

    This approximation has been particularly checked to deal with non-homogeneous systems.

  2. 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. [20].

    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 [21]:


    where εxbom is the exchange energy of an unpolarized density n(r) system. There are many forms of FXC, the most used are those introduced by Becke [22] and Perdew [23, 24].

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 Ωsell 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 i123 and Ni is the number of primitive cells in each direction (Ni in the case of perfect crystal).

Solving Eq. (24) is equivalent to increase the periodic function ui,kr, in a database-dependent functions points k:ϕjkrj=1Nbask:


where ϕjk 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 ϕjk 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


Ωcell is the cell volume of the original mesh in the real space and 2πd/Ωcell 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 [25], Chadi and Kohen [26], and Monkhorst and Pack [27] are the most frequently used.


3. Elastic constants and homogenized moduli of monoclinic structure

According to the crystal theory [28], 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 axP,TayP,TazP,TbxP,TbyP,TbzP,TcxP,TcyP,TczP,T and α1PTα6PTα5PTα6PTα2PTα4PTα5PTα4PTα3PT 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 [29], 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 ηij,ηkl, and ηmn are the coefficients of Lagrange deformation tensor, cijklTis the isothermal first-order elastic coefficients, and FηijT is the Helmholtz free energy.

The components of the stress tensor can be extracted by σi=j=16cijεj, 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 LCECLCEC
e=δδ0000ΔEV0=c112+c12+c222δ2c11+ c22 + 2c12
e=0δδ000ΔEV0=c222+c23+c332δ2c22 + c33 + 2c23
e=δ0δ000ΔEV0=c112+c13+c332δ2c11 + c33 + 2c13
e=000δδ0ΔEV0=c442+c45+c552δ2c44 + c55 + 2c45
e=δ0000δΔEV0=c112+c16+c662δ2c11 + c66 + 2c16
e=0δ000δΔEV0=c222+c26+c662δ2c22 + c66 + 2c26
e=00δ00δΔEV0=c332+c36+c662δ2c33 + c66 + 2c36
e=δδδ000ΔEV0=c112+c222+c332+c12+c13+c23δ2c11 + c22 + c33 + 2c12 + 2c13 + 2c23
e=δδδ2/1δ2000ΔEV0=c112c12+c222δ2c11 + c22 − 2c12
e=δδ2/1δ2δ000ΔEV0=c112c13+c332δ2c11 + c33 − 2c13
e=δ2/1δ2δδ000ΔEV0=c222c23+c332δ2c22 + c33 − 2c23
e=δ0002δ0ΔEV0=c112+2c15+2c55δ2c11 + 4c55 + 4c15
e=δ0002δ0ΔEV0=c1122c15+2c55δ2c11 + 4c55 − 4c15
e=0δ002δ0ΔEV0=c222+2c25+2c55δ2c22 + 4c55 + 4c25
e=00δ02δ0ΔEV0=c332+2c35+2c55δ2c33 + 4c55 + 4c35

Table 1.

Deformation tensors to calculate independent elastic constants of monoclinic crystal [30, 31].

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 [32]:


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 [32]

GR=15{ 4 [ (c33c55c352)(c11+c22+c12)+(c23c55c25c35)(c11c12c23) +(c13c35c15c33)(c15+c25)+( c13c55 c15c35 )(c22c12c23c13)+(c13c25c15c23)(c15c25)+f ]/Ω+3[ g/Ω+(c44+c66)/(c44c66c462) ] }1E38
BR=Ω [ (c33c55c352)(c11+c222c12)+(c23c55c25c35)(2c122c11c23) +(c13c35c15c33)(c152c25)+(c13c55c15c35)(2c12+2c23c132c22)+2(c13c25c15c23)(c25c15)+f ]1E40

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 [32]:


Young’s modulus and Poisson’s ratio can be rewritten based on the Voigt-Reuss-Hill approximation [33]. In terms of the Voigt-Reuss-Hill approximations [34], 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 [32] 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 [35].

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 [100] and zigzag chains along [001] direction [36] (Table 2).

Figure 1.

Modeling of gypsum crystal. (a) Gypsum structure [36] along [001]; (b) the real cell; (c) in x-direction; (d) in y-direction; and (e) in z-direction.

AtomxyzOccupancy rateUiso or Ueq

Table 2.

Atomic coordinates and displacement parameters of gypsum [36].

4.1.2. Nanoscale modeling of monoclinic 11 Å tobermorite crystal

Hamid model [37] 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 [37]: 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.

Figure 2.

Modeling of 11 Å tobermorite crystal. Silicate chains, calcium octahedral, and oxygen atoms are shown as yellow tetrahedra, green spheres, and red spheres. (a) 11 Å Tob monoclinic crystal; (b) in x-direction; (c) in y-direction; and (d) in z-direction.

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 [38]. Atomic coordinates and displacement parameters are seen in Table 3.

Atomic speciesXYZOccupancy rateUiso or UeqAtomic speciesXYZOccupancy rateUiso or Ueq

Table 3.

Atomic coordinates and displacement parameters of 11 Å tobermorite [38].

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.

Figure 3.

Gypsum monoclinic crystal under pressure 0–1.0 GPa by DFT. (a) Relative change of a, b, c, and V and (b) elastic constants.

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 [36].

Elastic constants of gypsum crystal model based on DFT are calculated, and parameters are detailed in Table 4.


Table 4.

Elastic coefficient Cij (GPa) of gypsum by DFT.

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.

Figure 4.

11 Å tobermorite monoclinic crystal under pressure 0–1.0 GPa by DFT. (a) Relative change of a, b, c, and V and (b) elastic constants.


Table 5.

Elastic coefficient Cij (GPa) of 11 Å tobermorite by DFT.

A comparisonal results of Shahsavari [39] 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 [40].

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.

Figure 5.

Elastic moduli of gypsum crystal under pressure 0–1.0 GPa.

As gypsum shows anisotropic compressibility along three crystallographic axes with b > c > a below 5 GPa [44], 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)μ
Reference [43]26.533339.255624.807739.238125.670539.246963.22650.2315

Table 6.

Mechanical moduli of gypsum polycrystalline by different methods.

As an acoustic method [41] and mechanical properties [42] have been investigated, according to elastic constants of gypsum crystal [43], 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 [44] 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 [42].

4.3.2. Elastic modulus of monoclinic tobermorite structure

Based on elastic constants of 11 Å tobermorite crystal using GGA calculation method by DFT, bulk modulus B and shear modulus G are separately calculated by Eqs. (37)(50) (Figure 6).

Figure 6.

Elastic moduli of 11 Å tobermorite crystal under pressure 0–1.0 GPa.

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 [45] by Pellenq and result of 78.939 GPa [39] 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)μ
Reference [191]54.213351.697634.156028.916852.955531.536478.93910.2516

Table 7.

Mechanical moduli of 11Å tobermorite polycrystalline by different methods.

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.


5. Conclusions

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:

  1. 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.

  2. 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.

  3. 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.

  4. 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.

  5. 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.


  1. 1. Messaoudi IS, Zaoui A, Ferhat M. Band-gap and phonon distribution in alkali halides [J]. Physica Status Solidi B. 2015;252(3):490-495. DOI: 10.1002/pssb.201451268
  2. 2. Thomas LH. The calculation of atomic fields, [C]. Mathematical Proceedings of the Cambridge Philosophical Society. Cambridge University Press. 1927;23(05):542-548. DOI: 10.1017/S0305004100011683
  3. 3. Dirac PAM. Note on exchange phenomena in the Thomas atom [J]. Mathematical Proceedings of the Cambridge Philosophical Society. 1930;26:376-385. DOI: 10.1017/S0305004100016108
  4. 4. Kohn W, Sham LJ. Self-consistent equations including exchange and correlation effects [J]. Physical Review. 1965;140(4A):A1133. DOI: 10.1103/PhysRev.140.A1133
  5. 5. Levy M, Perdew JP, Sahni V. Exact differential equation for the density and ionization energy of a many-particle system [J]. Physical Review A. 1984;30(5):2745. DOI: 10.1103/PhysRevA.30.2745
  6. 6. Foulkes WMC, Mitas L, Needs RJ, et al. Quantum Monte Carlo simulations of solids [J]. Reviews of Modern Physics. 2001;73(1):33. DOI: 10.1103/RevModPhys.73.33
  7. 7. Aulbur WG, Jönsson L, Wilkins JW. Quasiparticle calculations in solids [J]. Solid State Physics. 1999;54:1-218. DOI: 10.1016/S0081-1947(08)60248-9
  8. 8. Luther DIT. Homogenization of Damaged Concrete Meso-Structures Using Representative Volume Elements-Implementation and Application to Slang, Doctoral dissertation. Weimar Germany: Bauhaus-University; 2005
  9. 9. Zhu QZ, Kondo D, Shao JF. Micromechanical analysis of coupling between anisotropic damage and friction in quasi brittle materials: Role of the homogenization scheme. International Journal of Solids and Structures. 2008;45(5):1385-1405. DOI: 10.1016/j.ijsolstr.2007.09.026
  10. 10. Behnken H, Hauk V. Berechnung der röntgenographischen Elastizitäts-konstanten (REK) des Vielkristalls aus Einkristalldaten für beliebige Kristallsymmetrie[J]. Zeitschrift für Metallkunde. 1986;77:620-626
  11. 11. Gnäupel-Herold T. A software for diffraction stress factor calculations for textured materials[J]. Powder Diffraction. 2012;27(02):114-116. DOI: 10.1017/S0885715612000267
  12. 12. Kneer G. Die elastischen Konstanten quasiisotroper Vielkristallaggregate[J]. Physica Status Solidi B. 1963;3(9):K331-K335. DOI: 10.1002/pssb.19630030924
  13. 13. Reuss A. ZAMM – Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik. 1929;9(1):49-58. DOI: 10.1002/zamm.19290090104
  14. 14. Voigt W. Lehrbuch Der Kristallphysik Teubner, Leipzig 1910; Reprinted (1928) with an Additional Appendix. Leipzig, Teubner, New York: Johnson Reprint
  15. 15. Kamali-Bernard S, Bernard F. Effect of tensile cracking on diffusivity of mortar: 3D numerical modelling [J]. Computational Materials Science. 2009;47:178-185. DOI: 10.1016/j.commatsci.2009.07.005
  16. 16. Born M, Oppenheimer R. Zur quantentheorie der molekeln [J]. Annalen der Physik. 1927;389(20):457-484. DOI: 10.1002/andp.19273892002
  17. 17. Hohenberg P, Kohn W. Inhomogeneous electron gas [J]. Physical Review. 1964;136(3B):B864. DOI: 10.1103/PhysRev.136.B864
  18. 18. Ceperley DM, Alder BJ. Ground state of the electron gas by a stochastic method[J]. Physical Review Letters. 1980;45(7):566. DOI: 10.1103/PhysRevLett.45.566
  19. 19. Perdew JP, Zunger A. Self-interaction correction to density-functional approximations for many-electron systems [J]. Physical Review B. 1981;23(10):5048. DOI: 10.1103/PhysRevB.23.5048
  20. 20. Herman F, Van Dyke JP, Ortenburger IB. Improved statistical exchange approximation for inhomogeneous many-electron systems [J]. Physical Review Letters. 1969;22(16):807. DOI: 10.1103/PhysRevLett.22.807
  21. 21. Perdew JP, Burke K. Comparison shopping for a gradient-corrected density functional [J]. International Journal of Quantum Chemistry. 1996;57(3):309-319. DOI: 10.1002/(SICI)1097-461X(1996)57:33.3.CO;2-A
  22. 22. Becke AD. Density-functional exchange-energy approximation with correct asymptotic behavior [J]. Physical Review A. 1988;38(6):3098. DOI: 10.1103/PhysRevA.38.3098
  23. 23. Perdew JP, Wang Y. Accurate and simple analytic representation of the electron-gas correlation energy [J]. Physical Review B. 1992;45(23):13244. DOI: 10.1103/PhysRevB.45.13244
  24. 24. Perdew JP, Burke K, Ernzerhof M. Generalized gradient approximation made simple[J]. Physical Review Letters. 1996;77(18):3865. DOI: 10.1103/PhysRevLett.77.3865
  25. 25. Baldereschi A. Mean-value point in the Brillouin zone [J]. Physical Review B. 1973;7(12):5212. DOI: 10.1103/PhysRevB.7.5212. DOI: 10.1103/PhysRevB.13.5188
  26. 26. Chadi DJ, Cohen ML. Special points in the Brillouin zone [J]. Physical Review B. 1973;8(12):5747. DOI: 10.1103/PhysRevB.8.5747
  27. 27. Monkhorst HJ, Pack JD. Special points for Brillouin-zone integrations [J]. Physical Review B. 1976;13(12):5188
  28. 28. Shao TJ, Wen B, Melnik R, et al. Temperature dependent elastic constants for crystals with arbitrary symmetry: Combined first principles and continuum elasticity theory [J]. Journal of Applied Physics. 2012;111:083525. DOI: 10.1063/1.4704698
  29. 29. Bauernschmitt R, Ahlrichs R. Stability analysis for solutions of the closed shell Kohn–Sham equation [J]. The Journal of Chemical Physics. 1996;104(22):9047-9052. DOI: 10.1063/1.471637
  30. 30. Catti M. Calculation of elastic constants by the method of crystal static deformation [J]. Acta Crystallographica Section A. 1985;41:494-500. DOI: 10.1107/S0108767385001052
  31. 31. Ting TCT. Anisotropic Elasticity-Theory and Applications [M]. Oxford: Oxford University Press; 1996
  32. 32. Wu ZJ, Zhao EJ, Xiang HP, et al. Crystal structures and elastic properties of superhard IrN2 and IrN3 from first principles [J]. Physical Review B. 2007;76(5):054115. DOI: 10.1103/PhysRevB.76.054115
  33. 33. Hill R. The elastic behaviour of a crystalline aggregate[J]. Proceedings of the Physical Society. Section A. 1952;65(5):349. DOI: 10.1088/0370-1298/65/5/307
  34. 34. Raabe D. Computational Materials Science: The Simulation of Materials Microstructures and Properties. Weinheim: Wiley-VCH; 1998
  35. 35. Knight KS, Stretton IC, Schofield PF. Temperature evolution between 50 K and 320 K of the thermal expansion tensor of gypsum derived from neutron powder diffraction data[J]. Physics and Chemistry of Minerals. 1999;26(6):477-483. DOI: 10.1007/s002690050
  36. 36. Comodi P, Nazzareni, et al. High-pressure behavior of gypsum: A single-crystal X-ray study[J]. American Mineralogist. 2008;93(10):1530-1537. DOI: 10.2138/am.2008.2917
  37. 37. Hamid SΑ. The crystal structure of the 11Å natural tobermorite Ca2.25[Si3O7.5(OH)1.5]·1H2O[J]. Zeitschrift für Kristallographie-Crystalline Materials. 1981;154(1–4):189-198. DOI: 10.1524/zkri.1981.154.3-4.189
  38. 38. Merlino S, Bonaccorsi E, et al. The real structure of tobermorite 11Å normal and anomalous forms, OD character and polytypic modifications[J]. European Journal of Mineralogy. 2001;13(3):577-590. DOI: 10.1127/0935-1221/2001/0013-0577
  39. 39. Shahsavari R, Buehler MJ, Pellenq RJM, et al. First-principles study of elastic constants and interlayer interactions of complex hydrated oxides: Case study of tobermorite and jennite[J]. Journal of the American Ceramic Society. 2009;92(10):2323-2330. DOI: 10.1111/j.1551-2916.2009.03199.x
  40. 40. Vandamme M, Ulm FJ. Nanondentation investigation of creep properties of calcium silicate hydrates[J]. Cement and Concrete Research. 2013;52:38-52. DOI: 10.1016/j.cemconres.2013.05.006
  41. 41. Haussuhl S. Elastische und Thermoelastische Eingenschaften von CaSO4.2H2O (Gips)[J]. Zeitschrift für Kristallographie. 1960;122:311-314. DOI: 10.1515/zkri-1965-1-628
  42. 42. Watt JP. Hashin-Shtrikman bounds on the effective elastic moduli of polycrystals with monoclinic symmetry[J]. Journal of Applied Physics. 1980;51:1520-1524. DOI: 10.1063/1.327803
  43. 43. Meille S, Garboczi EJ. Linear elastic properties of 2D and 3D models of porous materials made from elongated objects[J]. Modelling and Simulation in Materials Science and Engineering. 2001;9(5):371. DOI: 10.1088/0965-0393/9/5/303
  44. 44. Huang E, Ku J, Lin J, Hu J. Pressure-induced phase transition in gypsum[J]. High Pressure Research. 2000;17:57-75. DOI: 10.1080/08957950008200306
  45. 45. Pellenq RJM, Lequeux N, Van Damme H. Engineering the bonding scheme in C-S-H: The iono-covalent framework [J]. Cement and Concrete Research. 2008;38(2):159-174. DOI: 10.1016/j.cemconres.2007.09.026
  46. 46. Miller M, Bobko C, Vandamme M, Ulm F-J. Surface roughness criteria for cement paste nanoindentation[J]. Cement and Concrete Research. 2008;38:467-476. DOI: 10.1016/j.cemconres.2007.11.014

Written By

Jia Fu

Submitted: 14 June 2017 Reviewed: 08 November 2017 Published: 20 December 2017