Calculation of the shear moduli

If an external force acting on a body or if one part of the body applies force on another part, it is said that the body is in the state of stress. Stress is defined in units of force per unit area $\sigma=F/A$ and can be characterized in general case by the stress tensor. Strain describes the state of deformation of a solid; there is dilatational strain which changes the volume, but not the shape and deviatoric strain which in contrast changes the shape, but not the volume.

For a field of deformations $\vec u(x,y,z)$, we can define a symmetric strain tensor in the following way:

$$\epsilon_{\alpha \beta}=\frac{1}{2}\left(\frac{\partial u_\a... ...ial u_\beta}{\partial x_\alpha}\right ),~~\alpha,\beta~=~1,2,3$$ (10.1)

Hook discovered that the strain $\epsilon$ is proportional to the stress $\sigma$:

$$\epsilon=S\sigma~~or~~\sigma=C\epsilon$$ (10.2)

These relation between the stress and the strain could be generalized for an anisotropic solid, in the terms of tensors:

$$\epsilon_{\alpha \beta}=\sum_{ \gamma \delta} S_{\alpha \beta \gamma \delta}\sigma_{ \gamma \delta}$$ (10.3)

or

$$\sigma_{\alpha \beta}=\sum_{ \gamma \delta} C_{\alpha \beta \gamma \delta}\epsilon_{ \gamma \delta}$$ (10.4)

The fourth-rank tensor $C_{\alpha \beta \gamma \delta}$ has $3x3x3x3=81$ components, but thanks to cubic symmetry (of simple cubic, fcc and bcc lattices) there are only three independent components of the elasticity tensor $C_{\alpha \beta \gamma \delta}$.

It is customary to reduce the number of subscripts by means of an abbreviated notation for pairs of coordinate directions, as follows [67]:
in the full tensor notation the indices of $C_{\alpha \beta \gamma \delta}$: 11, 22 , 33, 23 32, 13 31, 12 21,
while in the abbreviated notation they are $C_{\alpha \beta}$: 1, 2, 3, 4, 5, 6
and there is a correspondence:

$$11 \rightarrow 1=> C_{1111} \equiv C_{11}$$ (10.5)

$$22 \rightarrow 2=> C_{1122} \equiv C_{12}$$ (10.6)

$$33 \rightarrow 3=> C_{1133} \equiv C_{13}$$ (10.7)

$$23,32 \rightarrow 4=> C_{3223}=C_{2323} \equiv C_{44}$$ (10.8)

$$13,31 \rightarrow 5=> C_{3113}=C_{1313} \equiv C_{55}$$ (10.9)

$$21,12 \rightarrow 6=> C_{1221}=C_{2121} \equiv C_{66}$$ (10.10)

In these notation the three independent components of the elasticity tensor [67]:

$$C_{11},~C_{12},~C_{44}$$

Of course, one can use linear combinations of these elastic constants, to express other physical quantities such as: bulk modulus: $B=C_{11}+\frac{1}{2}C_{12}$ and shear modulus: $C'=\frac{1}{2}(C_{11}-C_{12})$ The free energy of a distortion can be expressed in powers of the strain tensor components. Assuming that the distortion of a crystal lattice too small to break the four-fold cubic symmetry and neglecting high order terms we can write the free energy as:

$$U_{el}=\sum_{\alpha \beta}C_{\alpha \beta}\epsilon_{\alpha \... ...on_{11})^2 +\frac{1}{2}C_{44}(\epsilon_{11}^2+\epsilon_{13}^2)$$ (10.11)

On the base of elasticity theory, Born [68] derived the general conditions for stability of a crystal lattice:

$$B>0,~ C'>0~ and~ C_{44}>0$$ (10.12)

and showed that the melting temperature could be found from the condition $C_{44}(T_m)=0$. Therefore, calculation of these elastic coefficients is very important in studying the melting transition. The elastic constants $C_{11},C_{12},C_{44}$ are calculated by means of fluctuation formulas obtained by Ray and Rahman [57,58]. It was shown that the elastic coefficients $C_{\alpha \beta \gamma \delta}$ are related to fluctuations of the stress tensor

$$<\sigma_{\alpha \gamma}\sigma_{\beta \rho}> -<\sigma_{\alpha \rho}><\sigma_{\beta \gamma }>$$ (10.13)

and the ensemble average of the Born term $B_{\alpha \beta \gamma \rho }$. The elastic coefficients are calculated in the following way:

$$\begin{array}{l} C_{\alpha \beta \gamma \rho }=\left(-{V}/{k... ...\gamma }> \right) +<B_{\alpha \beta \gamma \rho }> \end{array}$$ (10.14)

where $\sigma_{\alpha \beta}$ is the stress tensor, $\delta_{\alpha \beta}$ is Kroneker delta, and $<>$ is ensemble averaging, $N,V,T$ are the volume, the number of atoms and the temperature, respectively.

The Born term is defined as:

$$B_{\alpha \beta \gamma \delta }=\frac{1}{2V}\sum_{i,j:~i \ne... ...ac{r_{ij\alpha}r_{ij\beta}r_{ij\gamma}r_{ij\delta} }{r^2_{ij}}$$ (10.15)

here $\alpha \beta \gamma \delta = 1,2,3$ are Cartesian indices and $i,j=1,2,..,N_{neigh}$ numerate the neighbor atoms, and $r_{ij}=\vert\vec r_i-\vec r_j\vert$ is distance between the atoms $i$ and $j$. Potential energy is represented by the Finnis - Sinclair potential:

$$U_{pot} =\frac {1}{2} \sum_{i,j:~i \neq j}^N \Phi(r_{ij}) +\sum_{i=1}^N F(\rho_i)$$ (10.16)

and therefore the Born term is given by the formula:

$$B_{\alpha \beta \gamma \delta }=(B1)_{\alpha \beta \gamma \d... ...alpha \beta \gamma \delta }+(B3)_{\alpha \beta \gamma \delta }$$ (10.17)

where:

$$(B1)_{\alpha \beta \gamma \delta }= \frac{1}{2V}\sum_{i,j:~i... ...ac{r_{ij\alpha}r_{ij\beta}r_{ij\gamma}r_{ij\delta} }{r^2_{ij}}$$ (10.18)

$$(B2)_{\alpha \beta \gamma \delta }=\frac{1}{2V}\sum_{i,j:~i ... ...ac{r_{ij\alpha}r_{ij\beta}r_{ij\gamma}r_{ij\delta} }{r^2_{ij}}$$ (10.19)

here

$$\rho_i=\sum_j^{N_{neigh}}\rho(r_{ij})$$ (10.20)

and finally:

$$(B3)_{\alpha \beta \gamma \delta } =\sum_i^{N} F''[\rho_i]g_{i\alpha \beta}g_{i\gamma \delta}$$ (10.21)

while the function $g_{i\alpha \beta }$ defined as:

$$g_{i\alpha \beta }=\frac{1}{2}\sum_{j:~j \neq i}^{N_{neigh}} \frac{\rho'(r_{ij})r_{ij\alpha}r_{ij\beta} }{r_{ij}}$$ (10.22)

The main advantage of this method is the fast convergence of the stress fluctuation term to its equilibrium value.

In order to calculate the shear elastic modulus using the above mentioned formula two steps are required. The first step is to find the zero-stress reference matrix for the computational cell $H^0_{\alpha \beta}$ and the second step is to run the NVT MD simulation to calculate the elastic coefficients, using stress tensor fluctuations and average value of the $<B_{\alpha \beta \gamma \rho }>$ term.