M07 Computing Elastic Constants
Manual 07 for MD++
Calculating Elastic Constants
Keonwook Kang and Wei Cai
Elastic constants are physical properties of crystals to relate the mechanical response to the material deformation (i.e. stress and strain) within elastic regime. To calculate the elastic constants is one simple way of validating a potential model. There are two methods to compute elastic constants: one is to monitor (Virial) stress and the other is to monitor strain energy as a function of the applied strain. Conceptually, the two approaches are equivalent except that stress is a linear function of strain and the strain energy is a quadratic function of strain. In this manual, both methods will be explained to get three independent elastic constants(Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{11}} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{12}} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{44}} ) for the cubic crystal materials, e.g. silicon.
Using Virial Stress
If Virial stress is available to compute, the task of obtaining elastic constants is relatively easy. In crystals with cubic symmetry, there are three independent elastic constants: Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{11}} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{12}} and . is the proportional constant between the stress and the strain , and is the proportional constant between the stress and the strain . We can get Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{11}} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{12}} from one simulation by applying strain Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \epsilon_{11}} and monitoring two components of stress Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sigma_{11}} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sigma_{22}} . MD++ script si_elastic.script is given below.
# -*-shell-script-*-
# Calculate Elastic Constants of Si
# setnolog
setoverwrite
dirname = ~/Codes/MD++/runs/si_elastconst
#--------------------------------------------
#Create Perfect Lattice Configuration
crystalstructure = diamond-cubic
latticeconst = 5.430953e+00 #(A) for Si
latticesize = [ 1 0 0 3
0 1 0 3
0 0 1 3]
makecrystal
#-------------------------------------------------------------
#Conjugate-Gradient relaxation
conj_ftol = 1e-4 conj_itmax = 1000 conj_fevalmax = 1000
conj_fixbox = 1
eval relax eval finalcnfile = relaxed.cn writecn
#-------------------------------------------------------------
conj_fixbox = 1 #conj_monitor = 1 conj_summary = 1
input = [ 1 1 0.001 ] shiftbox relax eval
input = [ 1 1 9.99000999000999e-04 ] shiftbox relax eval
input = [ 1 1 9.98003992015968e-04 ] shiftbox relax eval
input = [ 1 1 9.97008973080758e-04 ] shiftbox relax eval
input = [ 1 1 9.96015936254980e-04 ] shiftbox relax eval
quit
You can run the script by typing
$ bin/sw gpp scripts/si elastic.script
The above script first generates a 3 × 3 × 3 perfect diamond cubic structure of Si. It then applies tensile strain in increments of 0.1% along [100] direction by the command shiftbox, and calculate the corresponding Virial stress at each strained state by the command eval. From the log file, the stress components Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sigma_{11}} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sigma_{22}} can be obtained and plotted as a function of the applied strain Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \epsilon_{11}} . Fig. 1(a) plots the result from running the script with Stillinger-Weber(SW) potential. The slopes of these lines give Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{11}=} 160.6 (GPa) and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{12}=} 80.5 (GPa).
To compute Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{44}} , which is the proportional constant between Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sigma_{12}} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \epsilon_{12}} , we can simply insert the following lines before quit in the above script.
incnfile = relaxed.cn readcn input = [ 1 2 0.001 ] shiftbox relax eval input = [ 1 2 0.001 ] shiftbox relax eval input = [ 1 2 0.001 ] shiftbox relax eval input = [ 1 2 0.001 ] shiftbox relax eval input = [ 1 2 0.001 ] shiftbox relax eval
These commands shear the crystal in increments of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \epsilon_{12}=} 0.1%. Again, from the log file, the stress component Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sigma_{12}} can be obtained and plotted as a function of the applied strain. The result from SW potential is plotted in Fig. 1(b), from which we obtain Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle C_{44}=} 60.2(GPa). The values for C11 , C12 , and C44 obtained here are consistent with the original paper.[1]
- Fig.1 Visualization of FCC structure
Using the potential energy
| Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathbf{r} = a\mathbf{H}\cdot \mathbf{s} } . |
Finding shear elastic constant for a non-tilting simulation cell
Notes
- ↑ H. Balamane, T. Halicioglu, and W. A. Tiler, PRB 46 2250-2279