The 37th CMD workshop“Advanced study with ES-opt”

Meta-dynamics simulation with DFT
Koichi KusakabeGrad. Sch. Eng. Sci, Osaka University

1Authors of this instruction◼ Dr. Takahiro Ishikawa @ NIMS
◼ Prof. Hitose Nagara @ Tottori university

with contributions by K. Kusakabe.

2

A problem of conventional MD
CrystalCreation ofnucleuses ina final form

The nucleus formation cause
loss in stabilization energy
resulted in an energy barrier.

When the barrier energy is larger than
the thermal fluctuation, MD cannot
search for the free energy minimum!

The Meta Dynamics
R. Martoňák et al., Phys. Rev. Lett. 90, 75503 (2003)

V Gaussian-shaped potentials

V

gradually fill the potential well.

xxVg ( x) = W exp( − x(t ) − x(ti )2t i th : Stepping parameter
W : Weight of a Gaussian2h 2 )The motion allows
the system fall
into the minimum.

The Meta Dynamics

Application of the meta-dynamics for a simulation ofthe pressure induced structural transformation
Simulation cell


a3


a2
a1

Action of the simulation cellmakes the system escapefrom the potential well.
(A step of this action iscalled a meta-step.)
Between adjacent meta-steps,
the system moves following
the conventional MD.

As a result, search of the free-energy minimum
is done on a nine-dimensional phase space(reduced to six-dimensional space by utilizing
a symmetric matrix representation) specified
by ( , , ).

Flowchart of Meta-dynamics① Constant h MD

Starting from a minimal
as a potential well.
Filling the FES

MD-step

Already visited h･ New atomic position･ Pressure tensor④ Filling the FE wellt② Driving force Fto⑥ Total driving force③ Update of h⑤ Driving force F g
By adding Gausspotentials, thesystem goes out ofthe potential well.

Metadynamics Algorithm
R. Mortonak et al., Phys. Rev. Lett. 90, 75503 (2003)② FFto ijto ij= − G  htottijp : Internal pressure tensor

 1 = − G  h =  {[(p − pext )h ]ij + [(p − pext )h ] ji }1 −  ij  2 ⑥③tijtoF =F +Fthtot +1FFtg= h + h Ft⑤ G (h ) =  Wetgt −1tttt

F− h t −h t (i )tt −1t: Driving force from original Got: Driving force from Gaussian
gt : Number of metastepsth : Stepping parameter
22h 2

W : Height of a Gaussiant ( i )ttt (i )−h−hWtttt(i)t⑥ F g ij = − Gg  h ij =  2 (h − h )et ( i ) t h22h 2

Metadynamics simulation method in ESopt
Metadynamics parts are written inUNIX shell scripts and perl scripts.
Ab-initio MD part uses ESopt.

How to run the examples
on the workshop machines
will be given in separate documents.
(This time, it is not distributed, because it is underreconstruction.)

Very Simple Example
Impurity migration• H in Silicon
Si ---- 8 atoma in a cube andon diamond latice points
H ----- 1→ escape from a minimum energy point

1) Try a run with default parameters
2) Try runs with diffrent parametersMore Complicate Example
Pressure induced structural phase transition•Si ---- 8 atoms in a cubeon diamond lattice points•Apply the pressure→ Unit cell shape will changePredicting Crystal Structures• Very simple equation of motion

G (h) = F (h) + PV
where G (h) is the Gibbs free energy
and F (h) is the Helmholtz free enrgy
as functions of h = (a, b, c) with V = det(h).

Driving force for the change of hG1t−1t−1−= V {[(p − p )h ]ij + [(p − p)h ] ji}(1 − ij )hij2Gt
F =−,ht
Fh t +1 = h t + h tFt

Internal stress P are calculatedby an ab-initio MD pacckageIn metadynamics, history dependent Gaussian
potential is added to get G t (h),whereG (h) = G(h) + åWet- h-ht' 2/2(d h)2t '<t

between t +1and t, MD run at constant h, P, and T
long enough to allow reluxation.

• Using metadynamics to get information
how to synthesize new materials• High pressure is one of those methods• Nanoscience→Large sample size simulationis necessary(pararell computor)

CMD-WS 13 @ IIAS
2008.9.10

Determining the structure of compressed calciumin phase Vby the ab-initio molecular dynamics simulation

Graduate School of Engineering Science, Osaka University
T. Ishikawa, A. Ichikawa, H. Nagara, K. Kusakabe and N. Suzuki

Table of contents
1.Previous research of compressed calcium
2.Computational details
3.Candidate structure of Ca-V
4.Candidate structure of Ca-IV
5.SummarySuperconductuvity in Caat high Pressures• High Tc : 25K at 161GPa• Highest recordin elements（2nd is Li 20 K at 48GPa)
T.Yabuuchi,T.Matsuoka,
Y.Nakamoto, and K. Shimizu,
J.Phys.Soc.JPN, 75, 83703(2006)

Previous research of compressed calcium
Unknown structures [1]

fcc bcc

20 32sc (simple cubic)

Ca-IV

113

Ca-V

139

[GPa]

• Decrease of the coordination number with increase of pressure.
12(fcc)→8(bcc)→6(sc)• The simple cubic structure can exist in wide pressure range.
• Superconducting transition temperature Tc linearly increases
with increase of pressure. → 25K at 161GPa [2][1] T. Yabuuchi et al., J. Phys. Soc. Jpn. 74, 2391 (2005).
[2] T. Yabuuchi et al., J. Phys. Soc. Jpn. 75, 83703 (2006).

Computational details• Metadynamics code + PWscf code
• Ultrasoft pseudopotential, GGA(Perdew-Wang 91)• Valence electron → 3s, 3p, 4s• k-space integration → 8×8×8, Cut-off energy →24Ry• Update of the cell by the steepest descent method (Meta-step)h : Stepping parameter•Driving force

GO : Gibbs free energy•Guassian potential

W : Weight of a Gaussian•120GPa, 100 MD-steps/1 meta-stepResult of the conventional MD
Simple cubic[1] T. Yabuuchi et al.,J. Phys. Soc. Jpn. 74, 2391 (2005)
8 atomsg
139 GPagg

4.5Å

Intensity

Ca-VStructure obtainedby the conventional MDgg

[1]ggCa-IVgg[1]Theory120GPa

127 GPa
gStructure obtained
120 GPaby the conventional
MDba 85.9°

10c152025302θ[deg]aa=3.05Å, b/a=1.08, c/a=1.44

The diffraction pattern does not matchthe experimental patterns!

Result of the metadynamics simulation
Switch-on

Orthorhombicba 85.9°caa=3.05Å, b/a=1.08, c/a=1.44

Switch-off

New structure(Monoclinic)

New structure obtained by the metadynamics

1’2’b

122 2’ cV=2cum
1 1’caγ=92.66°a• Interpenetrating structure oftwo monoclinic lattices• Zigzag pattern along c axis

g : Peaks from a metal gasket↓: Diffraction peaks from the sample

The zigzag structure is thecandidate structure of Ca-V!

Exploration of the structure of Ca-IV
Final structureof Ca-V

Predicted structure
 of Ca-IVcc

q = (0,0,0.5)

q = (0,0,0.25)

cc

cO

Most distinct change at thetransition from Ca-IV to Ca-V!

Final structureof Ca-VcO

c/cO=2

c/cO=4The peak around 17°disappears
by changing q from (0,0,0.5) to (0,0,0.25)!

Structures of Ca-IV and Ca-V
Candidate of Ca-IV

4’
3’ 1’2’4
3 12ba4 4’
3 3’cIV=4cum

2 2’c

1 1’

Helical patternaCandidate of Ca-V

1’2’b122 2’ cV=2cum
1 1’caγ=92.66° Zigzag patternaComparison of the enthalpy curves

Calculated Tc for the zigzag structure
μ*=0.1, k-points=(32,32,32), q-points=(4,4,4)P[GPa]

λωlog[K]

Tc[K]

100

1.1113245.0719.87

1101.0524240.0118.04

1200.9909

250.8717.23

1300.9149

294.3517.74

1400.8932

314.7918.19

Summary• The predicted structure of Ca-IV is theinterpenetrating structure of two tetragonal lattices
with the helical distortion.
• The predicted structure of Ca-V is theinterpenetrating structure of two monoclinic lattices
with the zigzag distortion.
•We calculated Tc for the zigzag structure and obtained
high values which are almost equal to the experimentaldata.