OriginalPaper
Some Applicaitons of Molecular Statics Method to Lattice Defects Problems
Yutaka TAKAHASIlI
(Department of MechanicalEngineeI･ing)
(Received Novemberl,1993)
New Tnethod to calculate equilibrium state of a.tom system
bonded with pair‑POtential was developed. The met,hod was
applied to mechanical property problems such as uniaxial
t,enSion and,laしt,ice defect,prOblems such as dislocationin t,WO‑
dimensionall,ennard‑Jones crystals.These pI･Oblems were able to
be soIvable by the continuum elast,ic t,heory and t,he accordanc(ユ
and dit‑ference were discussed.
Key Words: equjlibrium sLructure,pajr‑pOLential,Principle of
minimum potentialenergy,Ⅰ,ennard‑Jones crystal
1.Introduction
Material̲S are COnVent′ionally tr･eated by conLinuum theory, howevel･,they
act,ually consist of a number of atoms. Because bonding among atoms doII】inates
macroscopic mechanical pr･Operties such as (チ1ast,ic behavioI･, plastic
deformation and fracture toughness,many Studies in the atomic scale has been
performed by means of computer simulation to understand and combine
macroscopic and microscopic worldin ma†.erials. As calculuses for such
simulation, many kinds of methods wer･e already developed such as 'molecular
dynamics methodl)'and'Monte Carlo method 2))' 'car‑Parrinello's method
3)'･is thelatest,trendin this field. Principles of these met,hod are based
on Newton's equation, Stat,istical mechanics and quantum mechanics,
respect.ively.
Recently, the present aut,hor has developed new method to calculat,e
equilibrium state at zero Kelvin based on the principle of minimum of
potentialenergy4)I Of which conceptis wellknownin structuralmechanics.
The method was termed
asIJnOlecular statics methodI.Although Isimulated
annealing method5)IIImodified MI)methodIandIsteepesrt‑descent
method,have
been frequently applied to such equilibrium state calculations,the new
methodis superior to them with respect to convergency and calculation time.
In this reportISOme applications were presented to show the validity of
the new method and the applicability to defect problems.
2.Theory
Formulation of the methodis briefly denoted,Since detailwas discussed
in my previous paper･In a word,itis quite similar to finite element
method.
⊥土⊥.̲̲Pri坤IpOtent主旦L̲旦旦邑E已又
Force balanceis equalto a state where potentialenergy His minimized;
H = ◎ ‑ Wext(◎:internalpotentialenergy,Wext:eXternalwork)(1).
The state corresponds to equilibrium state at OK thermodynamically because H
is enthalpy and Hisidenticalwith free energy at OK.Ifinteraction between
atomsis expressed by pair‑POtential￠･(r),Eq.(1)is written as
H = ∑￠(=r(k)+u(k)‑r(1)‑u(1)‖)‑ ∑(u(k)･F(k))
k〉1 k
(2)
Where r(k)I u(k) and F(k)areinitialpositionI displacement vector and
applied force of k‑th
atom,and(Ⅹ･y)isinner product of vectors x and y.
▼2
Talor ex ansion a roximat迦
Note that ￠is function depending on distance r only,gradient and Hesse
matrix of ￠ are written as
￠ = ￠'‑
1txIl
llxtl￠''‑￠' ￠, (∂2￠/∂Ⅹi∂xj)=
【Ⅹ=
〔……j
llxll3 1lxll
,E:unit matrix].
Thus,◎in Eq.(1)is approximately expressed by
r(k)‑r(1)
◎=∑【￠(11r(k)‑r(1)‖)+￠(‖r(k)‑r(1)‖)(
k〉1 llr(k)‑r(1)ll
1
+ ‑(u(k)‑u(1)･Skl(u(k)‑u(1)))】,
2
u(k)‑u(1))
llr(k)‑r(1)‡l￠‑'‑￠
(Skl= ((ri(k)‑ri(1))(rj(k)‑rj(1)))
ltr(k)‑r(1)tl3
using Taylor series of second order.Eq･(2)is simply rewritten as
I H = α一(b･u)+‑(u･Su)
2
in the quadratic form of u =
髪 望遠R相通画届空軍 冨 圃醐通日冨岡
(5)
(6)
In Eq.(6),finding u such that His minimizedis equivalent to aH/∂u=0･
Therefore,following equationis obtained
Su = b. (7)
The linear
equation is soIvable for given boundary condition using
conventional calculuses such as Gauss's sweep‑Out method or conJugate
gradient method.
Because above explanation was derived from energy principle,
relationship between the present method and FEMis not clear･ According to
formulation base on force balance, b and u correspond to array vectors of
force and displacement and Sis overa11stiffness matrix combining b and u
linearly.
3.Examples
In following calculations, tWO‑dimensionalLennard‑Jones crystal was
treated.The pair‑pOtentialis given by
￠(r)= Ⅴ【(a/r)12‑2(a/r)6】 (8)
where two parameters(V, a)are potentialdepth and bonding distance between
two at,OmS,reSPeCtively.Elastic constants such as YoungIs modulus E,Poisson
ration v and shear modulus G are
E = 95.983 Va‑2,V = 0.3333,G = 35.842 Va 2
in the two‑dimensionalsystem fromlattice statics theory4)･
(1)■Uniaxial‑Lensile prou旦担
Results of elastic deformation when one‑directional tensile load is
applied is described first. Because finite rectangular‑Shape specimen was
treated in the previous paper4)andlocal inhomogeneous deformation and
stress concentration near surface and corners wasin problem,finite specimen
was calculated to prevent surface effect and edge effectin this calculation･
(1) Atoms were arrangedin hexagonalshape. Periodic boundary condition was
x‑and y‑directions of rectangular unit cell(size of ce｣1:1JXXLy,
number of atoms‥ n=576,inset of Fig.1).
(2)Equilibrium state was calculated for severalsets of(Lx,I,y).Intemal
potential◎is also calculated foI･the equilibrium staLe.Noload condiLion
COrreSpOnds to a state where￠is minimized for Lx and Ly.
(3)◎is approximately expressed by quaTtLic function oflongitudinalsLrain
Eyy and transversaユ strain Exx usingユeast‑Squar･e method for severalsets of
(fxx,ぎyy)=ぎ★ぃく3%,!ぎ,y蓼く3%ト
①=￠(どxx,どyy)
(9)
(4)Relations among tensile stress oyy aIld strains Exx,fyy WaS Obtained from
pt･inciple of vil･tualwork･If displacement6Exx and6Eyy are given,increment
Of q)is equalt,0 Virtualwork.
(Lx〔アy)(Lyy6ttyy)=(D(Exx+6Exx,Eyy+̀了Eyy)‑①(Exx,Eyy)
∂◎/∂.Fxx = 0
0yy =(∂◎/∂Eyy)/LxLy.
tn Fig.l,Open a11d closed
Ci上･Cles show calculated results
()f‑ uyy‑Eyy and Exx‑Eyy
T'elaLions frol11Eqs.(11) and
(12)･They are perfectlyin good
agreement with theoretical
Cu上●V(‑S Oflattice statics t.heory
draw‖ by solid lines in the
figure‑ Although Lhe exampleis
Very primary and notinteresting
for a problem of structural
englneering, it is confirmed
from the results that the
present method is valid for
Statistical equilibrium
Calculation.
z･ロ○>､丈b､ss2)SP巴コPむ∝
…U
ご馬上SO 5 0
0
∩)
∩〜
0
0
Fig･1 Stress‑Strain(oyy‑ぎyy)and
transversalstrain‑longitudinal
Strain(Exx‑Uyy)curvesin
iniaxialloading two‑dimensional
Lennard‑JorleS CryStal.The Number
Of atomsin a rectangular cell
is 576.
坤生臼旦旦̲豊里旦￡gヱ
Dislocations are the origin of plastic flow･Macroscopi,C deformationis
naturally achieved in consequence of collective motion andinteractions of
to know deformation.
A dislocation accommodates strain field characteri2;ed by direction of
dislocat.ion line and Burgers vector.Isotropic elastic theory has been well
used because the strain fieldis analytica11y soluble. The solution is,
however,incorrect physically because strainisinversely proportional to
distance r from dislocationline andit diverges at r=0. Thus,nOt COntinuum
approximation but discrete treatment ofindividualatoms are required for
discussion of dislocation core reglOn.
(a) Unit
｢Cell｢
ny=7[
ny
=x十1 甜壬弓Fauけed
盛葦灘迂葦丑葦PLane
ミミ票請甑だごご;::::ミ:ミm翳諮;岩
000C=⊃000=====000̀⊃(⊃00 000(⊃000=====00(⊃00000 00000000=====0000000
(⊃000000=====00000000 00000000=====0000000
0000〔)00=====000CI0000
nx=15
1∴.̲｣
Lx=15.Oq
(b)
毒頚蓋蒙
｣
Lx=15.6G
Fig.2 Procedure of dislocation core structure calculation;(a)initial
atom arrangement. Two crystals with differentlattice planes are
adhered.(b)Shrinkage of core occurred after severalit,erations.
(1) First, COre StruCture WaS Calculated by the present method.It is
explainedin Fig. 2 using smallspecimen.In Fig.2(a),tWO perfect,CryStals
were adhered where the number oflattice planes(nx)was not same. Because
extra planes wereinserted from top surfacein this case and Burgers vector b
is parallel to x‑direction, a dislocationlays on horizontal plane. Core
region was, however, Widely extended on the faulted plane and not clear in
the figure.
Periodic boundary condition wasimposed to x‑direction to avoid free
Surface effect.If not, because mirror force was applied to core from
Surface,dislocation slipped andit is sunk at the surface.
Asiterations were repeated untilforce balance was achieved, Shrinkage
of faulted plane occurred.In converged state of Fig. 2(b),PentagOnalatom
arrangement appeared which was dislocation core. One can always see such a
picture in a text book of dislocation theory as Bragg's bubble mode16).
Therefore,bubble bath was reproduciblein computer.Additionally,itis very
advantageous that any pair‑pOtential and any boundary condition were
availablein the present method.
(2) Secondly, dislocation energy was
discussed.Similarly,equilibrium of a
dislocation was calculated for n=840
atoms (Fig. 3). On condition that
Periodic boundary condition wasimposed
to x‑direction, Calculations were
Performed making ce111ength (Lx) a
●●●●●●●●l
●●●●●●●●●
l■■■■l■l■●l==●1
==●●●●●●●●
====●1
====●l
●●●●●●●●●●●●●
====●l
●●●●●●●●●●●●●
====●l
●●t■■●●●●●●●●
●●●●●●●●●●●●■
●●●●●●●●●●●●
====l
■■■■■■■■●●●●●●
■●●●●●●●●●●●l
==■●●●●●●
====l
==1■1■■l■■l■■ll====l
●●●●●●●●●●●●
l●●●●●●●●●●●l
free parameter and minimizing ￠.Because 〟㍊㌫〃㌫
X‑directional tensile load was applied
in t,his case, unloading condition was a
●●●
l●●1
●●●
====l
●●==■●●●●●
====1
==●●●●●●●
====1
==■●●●●●●●
●●●●●●●●●●●●l
Unit Cell
l=====●l
●●●●●●●●●●●●●●●●●
l====●●●●1 l======●1
●●■t●●●●●●●●●●●●●●●●
l======●l
====●●●●●●●●
l======●l
●●●●●●●●●●●●●●●●●●●●
l======●l ●1■l■■===l l==●●●●●●
====l l●●●●●●●●●●●●
===●●●●l l==●●■●●●
==l■●
l==●1
■●●●●●●●
■●●●●●●●●●●●●●●●●●●●
l======●l
======■●
l●●●l==l■●l■●l■t■■■■■■■■■■■l===●●●●●●●●●●
l==■●■●■●●■=●1
=●●●●■●■●■●●■●■●●●
l●■=====●●●1
●●●●●●●●●●●●●●●●●●●●
l●●●●●●●●● ●●●●●●
==●■==●l l●●●●●●●●●●●●●●●
●●●l■l■1■1■■●●●l■l■■l■■■
l●●●●●●●
==●l l●●●●●●●
●●●●●●●
l●●●●●●l
●●●●●●●●
l●●●●●●●l
●■●●■●●●
l●●●●●●●●●●●●●
====●l l====●
===●｣■l■l■l====●l l●●●●●●●●●●●●●
====■l●●●●●●●●●●●●●●●
l●●●●●●●●●●●●●●■
●●●●●●●●●●●●●●●
l●●●●●●●●●●●●●●l
●●●･===●●●
●●●●●●●●
l■==1
==･●●
==●1■●●●●●
l===●●●●
●■●●｣■●●●●●●●●
ll==●●●●●●●
=●●●●
l●●●●●●
●●●●●●●
l●●●●●●
l●●●●●●●
==●l
lt●●●●●●●●●l
●●●●●●●●●●●
l===●1 l●●●●●●■■■■■■●●
●●●●●●●●●●●●
l●●●●●●●●●●●●
●●●●●●●● ●●●
l=●■●■●■ ●●●
●●●●●●●●●●●●
l●●●●●●●●●●●●
●●●●●●●●●●●●
l●■==●■●
====1 1●●●●●●●●●●●
State Where ◎is minimi2;ed from Eq.(11).
Fig.3 Equilibrium atom
By comparing ￠ of perfect crystal with
arrangement of
the minimum value, dislocation energy r
dislocation in
was 25.1V or O.70Gb2 because the
n=840 Lennard‑
magnitude of Burgers vector b was a and
Jones crystal.
shear modulus G was 35.842Va‑2.
The result comparatively agreed with a value calculated by elastic
theory.According toisotroplC elastic theory,Strain energy of a dislocation
is estimated to be r=Gb2 7).
4.Conclusion
New method to calculate a static state in discrete atom systems was
applied to some probleTnS Of elastic deformation andlattice defect.The state
COrreSpOnds to force‑balancingin structural engineering and equilibrium
State in thermodynamics. Above examples may prove the usefulness of the
method for the study of atomic‑SCale simulations of static behavior.
References
l)L.Verlet:Phys.Rev.,.巳坦,98‑104(1967).
2)N.Metropolis,A.Rosenbluth,M.Rosenbluth,A.Teller and E･Teller:J･CheT･Phys･･Bi･831‑836(1953)･
3)R.Car and M.Parrlnello:Phys.Rev.Lett.,屋旦,2471‑2479 (1985)̀.
4)Y.Takahashi:J.Jpn.Inst.Metall.,51,597‑603(1993)【inJapanese】.
5)S.Kirpatrick,C.D.GelattJr.and M.P.Vecchi:IBM Res.
Rep.,RC‑9335(1982).
6)W.L.Bragg: Proc.Roy.Soc.London,全土旦旦,474‑478(1947), 皇皇旦旦,171‑174(1949).
7)H.Si2:uki:Introduction of Dislocation THeory,Agne,Tokyo(1980), pp.69‑77【inJapanese】.
