Transcript Invariant MD w/ Variable Cell Shape R. Wentzcovitch U

Invariant MD w/ Variable Cell Shape
R. Wentzcovitch
U. Minnesota
Vlab Tutorial
-Simulate solids at high PTs
-Useful for structural optimizations
-Useful for structural search (shake and bake)
-Various fictitious Lagrangian formulations
Fictitious molecular dynamics
(N,E,V)
H. C. Andersen (1978)
dU
mI RI  
dRI
(N,H,P)
dU
mI RI  
 fI
dRI
d (U  PextV )
WV  
F
dV
Variable Cell Shape MD
h2
h1
hij (t )
i=vector index
j=cart. index
dU
mI RI  
 fI
dRI
d (U  PextV )
Whij  
 Fij
dV
Anderson’s Fictious MD (HPN ensemble)
Anderson’s variable volume
fixed shape
constant pressure MD

L  K U
K at
(Anderson, J. Chem. Phys 72,2384(1980))
E  K  U  cte
U at
K cell
U cell
N
N
m
W 2
2/ 3
i T
 V  si .si   ij  rij   V  PV
2
i 1 2
i 1 i  j i1
N
LAnd
ri  V 1/ 3 si
V (t )  Cell volume
ri  V 1/ 3 si
E   Kat  Kcell  U at   Ucell  cte  " H "
The ensemble (trajectory) averages produce the HPN ensemble averages
Fictious MD (continue…)
Parrinello/Rahman variable cell shape MD
(Parrinello and Rahman, J. Appl. Phys 52, 7182 (1981))
 
N N
mi T
W
  si .gsi  ij  rij   Tr hT h  PV
2
i 1 2
i 1 i 1
N
LPR
Applying Lagrange’s equation
1
si  
mi

i 1
ij  rij 
rij
1
h    P  
W
   rij  T
1 N
1 N
T
   mi vi vi  
rij rij
V i 1
V i 1 j i rij

g  hT h
L
 L

t q
q
1
s

s

g
 i j  gsi

h  a, b , c
ri  hsi
" vi  ri  hsi "

  b  c , c  a, a  b
V
 h
I
2
T
2 T
 h1
V

- KLatt in PR-VCSMD is not uniquely defined
2a
a
a
a
a 0
h

0
a


a
h
0
v
h
0
 v 0
h

0
v


K Latt  Wv
2
equivalent
equivalent

K Latt
a

a
v

v
3 2
 Wv
2
The trajectory is not uniquely defined.
It does not depend only on the initial conditions.
Solution: use strain ε instead of h as dynamical variable
N N
mi T
W
  qi .dqi ij  rij   Tr  T   PV
2
i 1 2
i 1 j i
N
LInv
h  1    ho

ho  ao , bo , co
ε is strain

ri  (1   )qi
d  1    1   
T
V
T 1
    P 1   
W
1
qi  
mi
N
ij  rij 
i 1
rij

1
q

q

d
 i j  dqi
Invariant dynamics I
Wentzcovitch, PRB 44,2358 (1991)
Alternative form of LInv-I in terms of h and s:

LInv
fo  ( To  o )
h

N N
mi T
W
  si .gsi ij  rij   Tr hT f o h  PV
2
i 1 2
i 1 j i
N
V
  P   f o1
W
ij  rij 
1
si  

m
i i 1
rij
 s  s  g
i

 o  bo  co , co  ao , ao  bo
with
j
1

gsi
Final observation:
Inv  I
And
K Latt
 K cell
Solution: K
Inv  II
Latt
In the limit of variable V-only


W
W 2
T
 Tr hfh ~ V with f  ( T  )
2
2
Eq. of motion given by Eq. 9 in PRB 44, 2358(1991)
Fluctuations in the cell edges lengths of fcc
X-tal of Ar initially placed away from Veq.
Beeman integration algorithm
dt= 10 fmt (1 a.u. = 2.5 x 10-17s (in Ry))
Mi = 39 mp
W= 35 mp in (a); W= 0.0007 mp/ao3 in (b)
Rc= 10 ao
2a
a
a
a 0 0 


h  0 a 0 
 0 0 2a 


~
B
WV
Wentzcovitch, PRB 44,2358 (1991)
fcc
Potential
energy
isosurfaces
d
θ
sc
d
bcc
fcc
d
bcc
fcc
a b
a 0
1
h  a a
2
0 a
a
d
2
c
a

0
a 
Wentzcovitch, PRB 44,2358 (1991)
sc
Basins of
attraction if we
use a b and c
in the MD
bcc
Basins of
attraction if we
use a b and
c  a  b  c
in the PR-MD
Typical Computational Experiment
(Wentzcovitch, Martins, and Price, PRL 1993)
Damped dynamics (Wentzcovitch, 1991)
 ~ (   PI )
r ~ f int  f ( r, )
P = 150 GPa
hcp to bcc transition in Mg
(Wentzcovitch, Phys Rev. B 50, 10358 (1994))
u=1/6 or 1/3
(0001)
u=1/4
(110)
1
Atoms at  (ua  ub  c )
4
Distortion of the (0001) plane of the hcp structure into the (110) plane
of the bcc structure. Arrows indicate atomic displacements.
~150 K
Enthalpy barrier separating the
hcp from the bcc phases
at P=35 GPa at T=0K.
u=1/6 ↔ hcp
u=1/4 ↔ bcc
Ideal phase boundary (solid)
and blurry cause by hysteresis
(dashed). Phase transitions will
be simulated at the points
marked by dots and error bars
(undertainties in P and T).
Exp. PT = 45-55 GPa
at 300 K
u=1/6
u=1/4
u=1/6
u=1/4
hcp to bcc transition
u=1/4
Time evolution of the internal parameters
u’s, and angles and lengths of simulation
cell vectors.
Simulation w/ 16 atoms only
T = 700 K
P = 72 GPa
dt = 6 fts
W=0.02
mat=24.3 mp
Θab = 70.53o
Θab = 60o
u=1/6
u=1/4
bcc to hcp transition
u=1/6
u=1/4
u=1/4
Time evolution of the internal parameters
u’s, and angles and lengths of simulation
cell vectors.
Simulation w/ 16 atoms only
T = 500 K
P = 12 GPa
dt = 6 fts
W=0.02
mat=24.3 mp
Θab = 70.53o
Θab = 60o
MgSiO3 Perovskite
----- Most abundant constituent in the Earth’s lower mantle
----- Orthorhombic distorted perovskite structure (Pbnm, Z=4)
----- Its stability is important for understanding deep mantle (D” layer)
Crystal structure of post-perovskite
Tsuchiya, Tsuchiya, Umemoto, Wentzcovitch, EPSL, 2004
120 GPa
Exp
131
132
113
004
8
10
12
2 theta (deg)
Lattice system:
Bace-centered orthorhombic
Space group:
Cmcm
Formula unit [Z]: 4 (4)
Lattice parameters [Å]
a: 2.462 (4.286)
[120 GPa]
b: 8.053 (4.575)
c: 6.108 (6.286)
3
Volume [120 GPa] [Å ]:
121.1
(123.3)
( )…perovskite
042
041
040
111
002
021
020
6
a
Calc
110
b
c
023/130
022
Intensity (arbitrary unit)
Pt
 = 0.4134 Å
14
16
Ab initio exploration of post-perovskite phase in MgSiO3
- Reasonable polyhedra type and connectivity under ultra high pressure -
SiO4 chain
Perovskite
SiO3 layer
SiO3
Mg
SiO3
Mg
SiO3
MgSiO3
Structural relation between Pv and Post-pv
Tsuchiya, Tsuchiya, Umemoto, Wentzcovitch, EPSL, 2004
Perovskite
θ
b
c
a
b’
c’
a’
Post-perovskite
Deformation of perovskite under shear strain ε6
Conclusions
-VCSMD is very useful for structural optimizations when
the dynamics has the correct symmetry properties
(invariant dynamics)
- It is capable of simulating a phase transition when
one knows how the transformation occurs
- There is unavoidable hysteresis associated with
the simulation, which makes the simulation often
unfeasible
-Alternative approaches for obtaining phase boundaries
by computations will be discussed throughout the course
Practice
(Go to http://www.msi.umn.edu and navigate to the tutorial web site…
…to … software. You will use VCSMD today. Click and download program,
Input, and instruction.)
Some Instructions for Lind24-Lab
1)
OpenDX is a visualization software you may use. To enable access to OpenDX:
module load soft/opendx
module initadd soft/opendx
The first line enables the software for the current session, the second for every future session. Every
user will need to type those two lines, but once they do, the software will be permanently enabled for
your individual accounts.To launch the software, type 'dx'.
2) xmgr is a basic plotting software available in Linux. To launch it type ‘xmgr'.
3) The command for compiling fortran a code is 'f77'. It's part of the GCC 3.3.5 package built into Linux.
4) You can SSH to MSI machines. They are on a different network and use a different account, so you
will need to incorporate that into the command. For example, if your username is 'user' and the computer
is 'altix.msi.umn.edu', you would need to type ‘ssh [email protected]'.
5) They machines called lind24-01.itlabs.umn.edu, lind24-02.itlabs.umn.edu, etc, all the way up to
lind24-40.itlabs.umn.edu. Both OpenDX and Xmgr are graphical, so you'll need to enable X Forwarding
for the SSH connection if you're logging in remotely. Usually this can be done by adding the '-XY' flag to
your SSH command in Unix.
Run1
Test: md of Ar atom in fcc cell
(title)
nd
(calc)
s n
(ic,iio)
11.000000
(alatt)
1 1 1
(nsc)
1.000000
0.000000
0.000000
(avec)
0.000000
1.000000
0.000000
0.000000
0.000000
1.000000
0.00100 0.00000
(cmass, press)
1
(ntype)
4 Ar 40.00000
(natom,nameat,atmass)
0.000000
0.000000
0.000000
(rat)
0.500000
0.500000
0.000000
0.000000
0.500000
0.500000
0.500000
0.000000
0.500000
40.000000
(rcut)
5 5 5
(ncell)
1000 1110 10
(nstep,ntcheck,ntimes)
000.00000 0.00100 200.00000
(temp,ttol,dt)
~
Run2
Decrease step size by ½ and increase
# of steps by 2
Run3
Test: md of Ar atom in fcc cell
(title)
nd
(calc)
s n
(ic,iio)
11.000000
(alatt)
1 1 1
(nsc)
0.500000
0.500000
0.000000
(avec)
0.000000
0.500000
0.500000
0.500000
0.000000
0.500000
0.00100 0.00000
(cmass, press)
1
(ntype)
1 Ar 40.00000
(natom,nameat,atmass)
0.000000
0.000000
0.000000
(rat)
40.000000
(rcut)
9 9 9
(ncell)
2000 2110 10
(nstep,ntcheck,ntimes)
000.00000 0.00100 100.00000
(temp,ttol,dt)
~
Run4
Adjust cell mass to get same
period of oscillation
Run5
Test: Optimization under pressure (fcc)
(title)
nm
(calc)
s n
(ic,iio)
11.000000
(alatt)
1 1 1
(nsc)
1.000000
0.000000
0.000000
(avec)
0.000000
1.000000
0.000000
0.000000
0.000000
1.000000
0.00100 0.00000
(cmass, press)
1
(ntype)
4 Ar 40.00000
(natom,nameat,atmass)
0.000000
0.000000
0.000000
(rat)
0.500000
0.500000
0.000000
0.000000
0.500000
0.500000
0.500000
0.000000
0.500000
40.000000
(rcut)
6 6 6
(ncell)
100 1110 10
(nstep,ntcheck,ntimes)
000.00000 0.00100 500.00000
(temp,ttol,dt)
~
Run6
Test: Optimization under pressure (hcp)
(title)
nm
(calc)
s n
(ic,iio)
9.000000
(alatt)
1 1 1
(nsc)
1.000000
0.000000
0.000000
(avec)
0.500000 s 0.750000
0.000000
0.000000
0.000000
1.633000
0.00100 0.00000
(cmass, press)
1
(ntype)
2 Ar 40.00000
(natom,nameat,atmass)
0.000000
0.000000
0.000000
(rat)
t
1.000000 t 1.000000
0.500000
40.000000
(rcut)
9 9 9
(ncell)
100 1110 10
(nstep,ntcheck,ntimes)
000.00000 0.00100 500.00000
(temp,ttol,dt)
~
Run7
Test: MD of 32 atoms at 200K
(title)
md
(calc)
s n
(ic,iio)
10.000000
(alatt)
2 2 2
(nsc)
1.000000
0.000000
0.000000
(avec)
0.000000
1.000000
0.000000
0.000000
0.000000
1.000000
0.00100 0.00000
(cmass, press)
1
(ntype)
4 Ar 40.00000
(natom,nameat,atmass)
0.000000
0.000000
0.000000
(rat)
0.500000
0.500000
0.000000
0.000000
0.500000
0.500000
0.500000
0.000000
0.500000
40.000000
(rcut)
3 3 3
(ncell)
1000 100 10
(nstep,ntcheck,ntimes)
200.00000 0.00100 200.00000
(temp,ttol,dt)
Run8
Test: MD of 32 atoms at 2000K
(title)
md
(calc)
s n
(ic,iio)
10.000000
(alatt)
2 2 2
(nsc)
1.000000
0.000000
0.000000
(avec)
0.000000
1.000000
0.000000
0.000000
0.000000
1.000000
0.00100 0.00000
(cmass, press)
1
(ntype)
4 Ar 40.00000
(natom,nameat,atmass)
0.000000
0.000000
0.000000
(rat)
0.500000
0.500000
0.000000
0.000000
0.500000
0.500000
0.500000
0.000000
0.500000
40.000000
(rcut)
3 3 3
(ncell)
1000 100 10
(nstep,ntcheck,ntimes)
2000.00000 0.00100 100.00000
(temp,ttol,dt)