M07 Computing Elastic Constants - Revision history

Yanmingw at 19:29, 23 June 2011

2011-06-23T19:29:46Z

← Older revision		Revision as of 19:29, 23 June 2011
Line 103:		Line 103:
	\|}		\|}

	If the only nonzero strain is <math>\epsilon_{11}</math>, then the ~~streain~~ energy per volume		If the only nonzero strain is <math>\epsilon_{11}</math>, then the strain energy per volume
	becomes <math>W = C_{11} \epsilon^2 /2</math>.		becomes <math>W = C_{11} \epsilon^2 /2</math>.
	We can determine <math>C_{11}</math> by ﬁtting <math>W(\epsilon_{11})</math> to a parabola as shown in Fig.2(a).		We can determine <math>C_{11}</math> by ﬁtting <math>W(\epsilon_{11})</math> to a parabola as shown in Fig.2(a).

Caiwei at 06:07, 25 October 2008

2008-10-25T06:07:32Z

← Older revision		Revision as of 06:07, 25 October 2008
Line 1:		Line 1:
	<H1 ALIGN="CENTER"><FONT SIZE="-1">Manual 07 for MD++ </FONT>		<H1 ALIGN="CENTER"><FONT SIZE="-1">Manual 07 for MD++ </FONT>
	<BR>		<BR>
	~~Calculating~~ Elastic Constants</H1>		Computing Elastic Constants</H1>
	<DIV>		<DIV>

Kwkang: /* Tcl version of MD++ script */

2007-12-03T22:35:34Z

Tcl version of MD++ script

← Older revision		Revision as of 22:35, 3 December 2007
Line 399:		Line 399:
	plot(strain11, sigma11, '.-',strain11, sigma22,'.-r',...		plot(strain11, sigma11, '.-',strain11, sigma22,'.-r',...
	strain12, sigma12, '.-g'), grid		strain12, sigma12, '.-g'), grid
	xlabel('Strain~~, \epsilon_{11}~~','fontsize',16),		xlabel('Strain','fontsize',16),
	ylabel('Stress (MPa)','fontsize',16)		ylabel('Stress (MPa)','fontsize',16)
	legend('\~~sigma_~~{11}~~','~~\sigma_{22}','\sigma_{12}',2)		legend('\epsilon_{11} vs. \sigma_{11}','\epsilon_{11} vs. \sigma_{22}',...
			'\epsilon_{12} vs. \sigma_{12}',2)

	figure(2)		figure(2)

Kwkang at 22:01, 29 November 2007

2007-11-29T22:01:19Z

Show changes

Kwkang: /* Using the potential energy */

2007-11-29T18:58:18Z

Using the potential energy

← Older revision		Revision as of 18:58, 29 November 2007
Line 113:		Line 113:

	{\|border="0" align="center"		{\|border="0" align="center"
	\|<math> B = \frac{1}{2}(C_{11} + 2C_{12}) </math>		\|<math> B = \frac{1}{3}(C_{11} + 2C_{12}) </math>
	\|}		\|}

Kwkang: /* Using the potential energy */

2007-11-29T01:15:57Z

Using the potential energy

← Older revision		Revision as of 01:15, 29 November 2007
Line 103:		Line 103:
	\|}		\|}

	If the only nonzero strain is <math>\epsilon_{11}</math>, then the streain energy ~~becomes~~ ~~<math>W = C_{11} \epsilon^2 /2</math>.~~		If the only nonzero strain is <math>\epsilon_{11}</math>, then the streain energy per volume
			becomes <math>W = C_{11} \epsilon^2 /2</math>.
	We can determine <math>C_{11}</math> by ﬁtting <math>W(\epsilon_{11})</math> to a parabola as shown in Fig.2(a).		We can determine <math>C_{11}</math> by ﬁtting <math>W(\epsilon_{11})</math> to a parabola as shown in Fig.2(a).
	The elastic constnat <math>C_{44}</math> can be determined in a similar way as shown in Fig.2(b).		The elastic constnat <math>C_{44}</math> can be determined in a similar way as shown in Fig.2(b).

Kwkang: /* Finding shear elastic constant for a non-tilting simulation cell */

2007-11-28T21:57:49Z

Finding shear elastic constant for a non-tilting simulation cell

← Older revision		Revision as of 21:57, 28 November 2007
Line 130:		Line 130:

	== Finding shear elastic constant for a non-tilting simulation cell ==		== Finding shear elastic constant for a non-tilting simulation cell ==
			In the above sections, we compute <math>C_{44}</math> by shearing the simulation cell. The
			simulation cell was rectangular before the shear, but becomes tilted after the
			shear. However, sometimes a simulation code allows to use only a rectangular
			simulation box and we can not tilt the box at all.
			In this case, how can we ﬁnd the shear elastic constant <math>C_{44}</math>?
			Notice that <math>C_{11}</math> and <math>C_{12}</math> can be always obtained with a rectangular cell,
			when the repeat vectors of the cell are aligned with the cubic axes of the crystal.
			With <math>C_{11}</math> and <math>C_{12}</math>C12 known, we can determine <math>C_{44}</math>
			using the orthogonal transformation rule for the 4th order tensor.
			When the repeat vectors are not aligned with the cubic axes, we should be able to compute
			linear combination of <math>C_{11}</math>, <math>C_{12}</math> and <math>C_{44}</math>
			while keeping the simulation cell rectangular.
			Remember that <math>C_{11}</math>, <math>C_{12}</math> and <math>C_{44}</math> are components
			of a 4th order tensor which transforms as

			{\| border = 0 align = center
			\| <math> C_{ijkl} = T_{mi} T_{nj} T_{sk} T_{tl} C_{mnst} </math>
			\|}

			where <math>T_{ij} = \mathbf{e}_i \cdot \mathbf{e}'_j</math> is relation between two coordinate systems
			{<math>\mathbf{e}_i</math>} and {<math>\mathbf{e}'_j</math>}.

			For example, if we have one coordinate system <math>\mathbf{e}_1</math> = [1 0 0],
			<math>\mathbf{e}_2</math> = [0 1 0] and <math>\mathbf{e}_3</math> = [0 0 1],
			and another coordinate system <math>\mathbf{e}'_1</math> = [1 1 0]/ \sqrt{2},
			<math>\mathbf{e}'_2</math> = [−1 1 0]/ \sqrt{2} and <math>\mathbf{e}'_3</math> = [0 0 1],
			then the matrix [T] is

			{\| border = 0 align = center
			\|<math> [T] = \begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} & 0 \\
			1/\sqrt{2} & 1/\sqrt{2} & 0 \\
			0 & 0 & 1 \end{pmatrix}</math>
			\|}

			In particular,

			{\| border = 0 align = center
			\|<math> \begin{align} C_{1111} &= T_{m1} T_{n1} T_{s1} T_{t1} C_{mnst} \\
			&= \sum_{i=1}^3 T_{i1} T_{i1} T_{i1} T_{i1} C_{11} \\
			& + 2(T_{11}T_{11}T_{21}T_{21} + T_{11}T_{11}T_{31}T_{31} + T_{21}T_{21}T_{31}T_{31})C_{12} \\
			& + 4(T_{21}T_{31}T_{21}T_{31} + T_{11}T_{31}T_{11}T_{31} + T_{11}T_{21}T_{11}T_{21})C_{44}.
			\end{align}</math>
			\|}

			Suppose we create a perfect crystal such that x, y, and z axes of the simulation
			cell are aligned with <math>\mathbf{e}'_1</math>, <math>\mathbf{e}'_2</math> and <math>\mathbf{e}'_3</math>.
			Then by straining the crystal along x axis,
			we can obtain <maht>C'_{1111}</math> = 180.4 (GPa) using the SW potential. Given <math>C_{11}</math> = 160.6,
			<math>C_{12}</math> = 80.5, the matlab script below gives <math>C_{44}</math> = 59.9(GPa) for silicon.

			<pre>
			% With given C11, C12, and C1111’, find C44 in cubic materials.
			% C11, C12, and C44 are the elastic constants in Cartesian
			% coordinate system.
			% C1111’ is the elastic constant in the new coordinate.

			syms C11 C12 C44
			C11n = 160.6; C12n = 80.5; C44n = 60.2;
			C1111prime = 180.40;
			xyz0 = [ 1 0 0; % original coordinate
			0 1 0;
			0 0 1 ];
			xyz1 = [ 1 1 0; % new coordinate
			-1 1 0;
			0 0 1 ];
			T = zeros(3,3);
			for i=1:3
			for j=1:3
			T(i,j) = dot(xyz0(i,:),xyz1(j,:)/norm(xyz1(j,:)));
			end
			end

			p=1; q=1; r=1; s=1;
			A = T(1,p)T(1,q)T(1,r)T(1,s) + T(2,p)T(2,q)T(2,r)T(2,s) + ...
			T(3,p)T(3,q)T(3,r)*T(3,s);
			B = T(1,p)T(1,q)T(2,r)T(2,s) + T(1,p)T(1,q)T(3,r)T(3,s) + ...
			T(2,p)T(2,q)T(3,r)*T(3,s);
			C = T(2,p)T(3,q)T(2,r)T(3,s) + T(1,p)T(3,q)T(1,r)T(3,s) + ...
			T(1,p)T(2,q)T(1,r)*T(2,s);
			CC = AC11 + 2BC12 + 4C*C44;
			cc = subs(CC, C11, C11n); cc = subs(cc, C12, C12n);

			disp(sprintf(’C44 = %f’,double(solve(cc-C1111prime,C44))))
			</pre>

	== Notes ==		== Notes ==

Kwkang: /* Using the potential energy */

2007-11-28T21:28:34Z

Using the potential energy

@@ Line 94: / Line 94: @@
 == Using the potential energy ==
 {|border="0" align="center"
-|<math> \mathbf{r} = a\mathbf{H}\cdot \mathbf{s} </math>.
+|<math> W = \frac{1}{2}\sigma_{ij}\epsilon_{ij} = \frac{1}{2}\epsilon_{ij} C_{ijkl} \epsilon_{kl}</math>.
 |}
 == Finding shear elastic constant for a non-tilting simulation cell ==

Kwkang: /* Using Virial Stress */

2007-11-28T21:15:10Z

Using Virial Stress

← Older revision		Revision as of 21:15, 28 November 2007
Line 88:		Line 88:
	here are consistent with the original paper.<ref>H. Balamane, T. Halicioglu, and W. A. Tiler, PRB '''46''' 2250-2279</ref>		here are consistent with the original paper.<ref>H. Balamane, T. Halicioglu, and W. A. Tiler, PRB '''46''' 2250-2279</ref>

	<gallery caption="Fig.1 ~~Visualization of FCC structure~~" widths="~~200px~~" heights="~~200px~~" perrow="3">		<gallery caption="Fig.1" widths="300px" heights="300px" perrow="2">
	Image:~~fcc~~.jpg\|(a) ~~Conventional image of FCC unit cell~~		Image:C11_C12.jpg\|(a)
			Image:C44.jpg\|(b)
	Image:al1X1X1.jpg\|(b) FCC unit cell in MD++. One corner atom and three face-centered atoms per cell
	Image:al3X3X3.jpg\|(c) FCC of 3X3X3 cells
	</gallery>		</gallery>

Kwkang: New page:

Manual 07 for MD++
Calculating Elastic Constants

Keonwook Kang and Wei Cai

2007-11-28T21:00:10Z

New page: <H1 ALIGN="CENTER"><FONT SIZE="-1">Manual 07 for MD++ </FONT> <BR> Calculating Elastic Constants</H1> <DIV> <P ALIGN="CENTER"><STRONG>Keonwook Kang and Wei Cai</STRONG></P> </DIV> <P ALI...

New page

<H1 ALIGN="CENTER"><FONT SIZE="-1">Manual 07 for MD++ </FONT>
<BR>
Calculating Elastic Constants</H1>
<DIV>

<P ALIGN="CENTER"><STRONG>Keonwook Kang and Wei Cai</STRONG></P>
</DIV>
<P ALIGN="CENTER">[[ Nov 28 ]], [[2007]]</P>
<P>
<BR><HR>

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(<math>C_{11}</math>, <math>C_{12}</math> and <math>C_{44}</math>) for the cubic crystal materials, ''e.g.'' silicon.
<BR><HR>

== 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: <math>C_{11}</math>, <math>C_{12}</math> and <math>C_{44}</math>. <math>C_{11}</math> is the proportional constant between the stress
<math>\sigma_{11}</math> and the strain <math>\epsilon_{11}</math>, and <math>C_{12}</math>
is the proportional constant between the stress <math>\sigma_{22}</math> and the strain <math>\epsilon_{11}</math>.
We can get <math>C_{11}</math> and <math>C_{12}</math> from one simulation by applying strain <math>\epsilon_{11}</math> and
monitoring two components of stress <math>\sigma_{11}</math> and <math>\sigma_{22}</math>.
MD++ script '''si_elastic.script''' is given below.

<pre>
# -*-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
</pre>

You can run the script by typing

$ bin/sw gpp scripts/si elastic.script

The above script ﬁrst 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 ﬁle, the stress components <math>\sigma_{11}</math> and <math>\sigma_{22}</math>
can be obtained and plotted as a function of the applied strain <math>\epsilon_{11}</math>. Fig. 1(a) plots the result from
running the script with Stillinger-Weber(SW) potential.
The slopes of these lines give <math>C_{11}=</math> 160.6 (GPa) and <math>C_{12}=</math> 80.5 (GPa).

To compute <math>C_{44}</math>, which is the proportional constant between <math>\sigma_{12}</math> and <math>\epsilon_{12}</math>,
we can simply insert the following lines before quit in the above script.

<pre>
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
</pre>

These commands shear the crystal in increments of <math>\epsilon_{12}=</math> 0.1%. Again, from the
log ﬁle, the stress component <math>\sigma_{12}</math> 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 <math>C_{44}=</math> 60.2(GPa). The values for C11 , C12 , and C44 obtained
here are consistent with the original paper.<ref>H. Balamane, T. Halicioglu, and W. A. Tiler, PRB '''46''' 2250-2279</ref>

<gallery caption="Fig.1 Visualization of FCC structure" widths="200px" heights="200px" perrow="3">
Image:fcc.jpg|(a) Conventional image of FCC unit cell
Image:al1X1X1.jpg|(b) FCC unit cell in MD++. One corner atom and three face-centered atoms per cell
Image:al3X3X3.jpg|(c) FCC of 3X3X3 cells
</gallery>

== Using the potential energy ==

{|border="0" align="center"
|<math> \mathbf{r} = a\mathbf{H}\cdot \mathbf{s} </math>.
|}

== Finding shear elastic constant for a non-tilting simulation cell ==

== Notes ==

<references/>

M07 Computing Elastic Constants - Revision history

Yanmingw at 19:29, 23 June 2011

Caiwei at 06:07, 25 October 2008

Kwkang: /* Tcl version of MD++ script */

Kwkang at 22:01, 29 November 2007

Kwkang: /* Using the potential energy */

Kwkang: /* Using the potential energy */

Kwkang: /* Finding shear elastic constant for a non-tilting simulation cell */

Kwkang: /* Using the potential energy */

Kwkang: /* Using Virial Stress */

Kwkang: New page:

Manual 07 for MD++ Calculating Elastic Constants

Manual 07 for MD++
Calculating Elastic Constants