20231109OpenMX Developersʼ Meeting @ ISSP

Atomic Stress Calculation in OpenMXTiCSiYoshinori Shiihara
Graduate School of Engineering, Toyota TechnologicalInstitute, Nagoya, Japan
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

1/27

Overview• Background• Atomic stress: formulation and Implementation
in OpenMX• Applications
De

•

•

Atomic stress distribution in multicomponent alloys and its origins•

Elastic heterogeneity in hetero-nanostructures••

De

TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

2/27

Background

TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN3/27

Materials Mechanics from ab initio
Ab initio calculation gives energy, force and stress based on quantum mechanics.
Ideal strength

Cell stressAtomic force
Total energy
Atomic structure of MgZn2(periodic boundary condition)

Cell stress

Wave function → electron density

Youngʼs modulus
Poisson rate

Cell strain

Stress-strain curve on cellStress distribution is not available in conventional methods.
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

4/27

Why atomic stress?
Elastic interaction between dislocation and solute atom
Strain energy
by dislocation

Pressure
Pressure bydislocation

Volume change by dislocation

Bulk modulus, Volumetric strain

Pressure bydislocation

Pressure by defectreflecting characterisity of defectVolumetricchange

Elastic interaction can be understoodbased on the electronic structure.

TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN5/27

Stress in continuum mechanics
Momentum conservation laws

Stress vector∫ ρadV = ∫ ρgdV + ∫ t dS

VDensity

Vacceleration Body force

S S
Surface force∫V ρr × adV = ∫V ρr × gdV + ∫ S r × tdSt n = σ(r)⋅ nSurface normalStress tensorρ a = divσ(r)+ ρ gf (r) = divσ(r)σ T (r) = σ(r)Divergence of stress is equal
to the forve acting on a material point.

Energy conservation lawForceDeformation∫ δ R(r)dV = ∫ f (r)⋅ δ u(r)dV = ∫ divσ(r)⋅ δ u(r)dV
WorkTOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN⎛ ∂F(r) ⎞σ(r) = ⎜⎟⎝ ∂ε(r) ⎠T
Stress field: strain derivativeof free energy density
6/27

Stress in atomic levelMomentum conservation lawikσ̂(k) = −∑ fi exp(ik ⋅ ri )Discretization Fourier transformationdivσ(r) = ∑ fiδ (r − ri )ii

Two-body force= −∑∑ fij exp(ik ⋅ ri )ii≠ jInverse Fourier
1
1σ = ∫ σ(r)dr =rij ⊗ fijtransformation∑∑VV
2V i j≠i
1σ(r)=rij ⊗ fijδ (r − ri )Cell stress∑∑1
2 i i≠ j
Basically, two-body force is not
Integrateandσi =r⊗f∑
ij
ij
averagedefined in ab initio calculation
Atomic stress 2Vi j≠iEnergy conservation law：
Cell stress

Total energy

1 ∂Etotσ=V ∂εCell strain

Atomic energy

1 ∂Eiσi =Atomic stress Vi ∂ε

TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

Inconsistent with the definition
by internal force?

Fig: Periodic boundary
condition cell

7/27

Definition of ab initio atomic level stress
Cell-strain derivative of atomic energy
Cell stress（Nielsen-Martin）

Depending on
atomic energy

Atomic stressAtom-wise decomposed
total energy

Total energy of cellE1 ∂Φσ αβ =
Ω ∂εαβCell volume

1 ∂ΦEiσ αβ , i =Ωi ∂εαβAtomic volume

Cell strain

Practical implementation LSMS- based methodof atomic stress

Cell strain

QMASOpenMX(ORNL)(AIST)(Univ. Tokyo)

Voronoi

Bader

Atomic orbital

Decomposedenergyʼs response
against cell strain

Stress decomposition should solve the non-uniqueness of decomposed stress.
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

8/27

Problem: Non-uniqueness of energy/stress density
Any function integrates tozero can be added to the total
stress/energy.

Total stress and total energy on DFT (under periodic boundary condition)
Total stresss ab =¶Etot= ò s ab (r )dr = ò s ab (r )+g (r)drWcWc¶e αβò g (r ) dr = 0 ò f (r ) dr = 0

Stress density*1

Total energy“Gauge-dependent problem” *1,2

Etot = ò e tot ( r ) dr = ò e tot ( r ) +f ( r ) drWc

Wc

Wc

*1Filippetti, et al., PRB, 61, 8433 (2000).
*2Chetty, et.al, PRB, 45, 6074 (1992).

Wc

Energy density*2

Ex. Kinetic energy density
Symmetric KED

2*ts (r ) =fÑy(r ) × Ñy i (r )åii2m i2
Generalized KED*1: t (r ) = ts (r ) + gÑ r (r )Stress componentProblem: find partial region where the gaugedependent term integrates to zero.
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN¶ ¶r (r )¶ra ¶rb¶ ¶r (r )òW¢ ¶ra ¶rb = 099/27

Solution: Stress decomposition based on atomic orbital
Locality can avoid the non-uniqueness problem.

Two KED distributions on potential-well model
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPANSim et al., JCP, (2003)

10 /27

Formulation and Implementation in OpenMX

TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

11 /27

Formulation of stress decomposition in OpenMX
Atomic stress is decompositionof total stress.

See https://www.openmx-square.org/tech_notes/tech23-1_0.pdf

Almost same formulation
Stress calculation is like force
calculation.
Atomic coordinate derivative w.r.t h axis

Strain derivative w.r.t gh component

Stress formulation is similar
to that of atomic force.
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

12 /27

Formulation of stress decomposition in OpenMX

TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

13 /27

Stress formulation: two-center integral (1)Two-center integral: Kinetic energy, etc.
AO center

Cell vector(R n )
ρ s,iα , j β φ iα (r − t i ) V̂ φ j β (r − t j − R n )Density matrix

Strain derivative

Potential operator

Basis function∂ (Rn )
ρ s,iα , jβ φiα (r − t i ) V̂ φ jβ (r − t j − R n )∂εγη=(R n )
∂ρ s,iα , jβ∂εγη(R n )+ ρ s,iα , jβφiα (r − t i ) V̂ φ jβ (r − t j − R n )Overlap term∂φiα (r − t i ) V̂ φ jβ (r − t j − R n )∂εγηr − t j − Rn → r

φiα (r − t i ) V̂ φ jβ (r − t j − R n )t ji,n = t j + R n − t i→ φiα (r) V̂ φ jβ (r + t i − t j − R n )t ji,nε = ( I + ε ) ⋅ t ji,n= φiα (r) V̂ φ jβ (r − t ji,n )TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

14 /27

Stress formulation: two-center integral (2)∂φiα (r) V̂ φ jβ (r − t ji,n )∂εγη) ∂,= ++ γ φiα (r) V̂ φ jβ (r − t ji,n ) .. t ηji,n* ∂t ji,n∂φiα (r) V̂ φ jβ (r − t ji,n )∂εγη) ∂,η= + γ φiα (r − t i ) V̂ φ jβ (r − t j − R n ) . tij,n* ∂ti-t ji,nε = ( I + ε ) ⋅ t ji,n∂tι= διη tη∂εγηt ij,n = t i − t j − R nForce component acting on atom i
Stress between atom i and j : sDecomposition rule
1/2 s for atom i
1/2 s for atom jTOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN15 /27

Stress formulation: integral over cellStress components from Hartree energy, exchange-correlation energy, etc.

∫ n(r)V (r)drVB

Strain derivative∂n(r)V (r)dr∫V∂εγη B= δγη ∫ n(r)Vdr + ∫Product rule

VB(in GGA)∂n(r)∂V (r)
V (r)dr + ∫ n(r)drVB ∂εVB∂εγη
γη

∫ f [n, ∇n ,n ]dr
VBXC

pcc

Strain derivative ∂ ∫ fXC [n,∇n,npcc ]dr∂εγη VB= δγη ∫ fXC dr + ∫Product & chain rules

VBTOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPANForce-like terms
appear as well

∂ fXC ∂n(r)∂ fXC ∂∇n(r)dr + ∫⋅
drVB ∂n(r) ∂ε
VB ∂∇n(r)∂εγηγη
16 /27

Applications•

Atomic stress distribution in multicomponent random alloys and its origins
Y. Shiihara, Y. Itai, I. Lobzenko, T. Tsuru, Frontiers in Materials, 9, 895626 (2022).

•

Elastic heterogeneity in hetero-nanostructures••

De

Partly, R. Hossain, H. Kimizuka, Y. Shiihara, and S. Ogata, Comput. Mater. Sci., 209, 111366 (2022).

TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN17 /27

Computational models

Y. Shiihara, Y. Itai, I. Lobzenko, T. Tsuru, Frontiers in Materials, 9, 895626 (2022).

n Atomic structures
fcc
CrMnFeCoNi types

bcc
VNbTaMoW types

TiZrHfNbTa + a types

Fig.: Atomic stress in CrFeMnCoNi
CrMnFeCoNi 60 atoms

VNbTaMoW 60 atoms

TiZrHfNbTa 60 atoms# of atoms: bi: 60, ter: 60,qua: 60 , qui: 60, sep: 60(10 variants ✕ 5 cases)# of atoms: bi: 60, ter: 60,qua 60: , qui: 60, sep: 60（10 variants ✕ 5 cases）# of atoms: bi: 60, ter: 60,qua 60: , qui: 60, sep: 60# of elements bi: binary, ter: ternary, qua: quarternary, qui: quinary, sep: septenaryn Ab initio calculation
OpenMX, Spin-polarized DFT , Norm-conserving pseudo potential, PBE-GGA,
k-grid: 0.022 Å, Electronic Temp.: 1500 K, FFT grid: 400 Ry
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN

Fig.: Atomic stress in VNbTaMoW

18 /27

Yield stress vs atomic stress dispersion
Oh, et al., Nature commun. (2019).

FCC
[MPa]

[MPa]BCC

Fig.: Friction stress vs lattice
distortion. Zhao et al. (2019)

Break line：Oh et al. (2019)

Fig.: Pressure distribution
in a random alloy

Fig: Standard deviation of atomic pressure vs experimentallyestimated friction stress (0K) in fcc random alloys.

TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPANFig: Standard deviation of atomic pressure vs yieldstress (experiment, room temp.) in bcc random alloys.

19 /27

How can tensile/compressive atomic stress appear?

Fig: Electronegativity and atomic volume of each element• Fcc: charge transfer，bcc: volume difference• Both should simultaneously affect!

Fig: Atomic pressure vs atomic quantities (charge transfer/volume difference)
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN20 /27

Regression by atomic-level quantities
Red bar: charge-related, blue bar: volume-related

fcc

Features：
E1: Charge transfer from thebulk
E2: Charge diff. from 1st NNave
V1: Vol. diff. from the bulk
V2: Vol. diff. from 1st NN ave.
V3: Vol. diff. from cell ave.
Fig: Regression accuracy of atomic pressurefor fcc and bcc alloys.

fcc
Fig: Feature importance of regression for fcc and bcc alloys

Clear difference: charge transfer (fcc), volume change(bcc)
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN21 /27

Atomic stress determined by 1st neighbors
Regressed by # of atoms for each element in 1st shell
fcc

bcc

We obtained rational explanation for the origin of atomic stress.
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN22 /27

Local elasticity calculation in hetero-nano structures
Elastic constants:

Hetero-nano struct.

StressFirst-principles method:

Elastic constant matrix

StrainTi
CTi
SiCTiTi
CTi

Cell elastic constantSiTiCTi

Change in cell stress

Given cell stress
Local elastic constant Change in local stress

MAX phase: Ti3SiC2M. Higashi, et al., Acta Materialia 161(2018) 161-170

Change in local strain••••••

Total energy - OK
Atomic force - OK
Cell stress - OK
Cell elastic constant - OK
Local stress - NG
Local elastic constant - NG

To calculate local elasticity, we need local stress and local strain.

TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN23 /27

Local elastic constant calculation
Local elastic constantsAtomicstress
Local volume
Voronoidecomposition

Atomic strain*

Atomic deformation gradient:minimize

Reference configuration

Atomic Green-Lagrange strain of atom i
Current configuration
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN*P. M. Gullet, Model. Simul. Mater. Sci. Eng., 16, 015001 (2008)24 /27

Computational settings
MAX phases
・Ab initio software: OpenMX(0001) stacking

Hardlayer

・Norm conserving pseudopotential

Ceramics Ti
Clayer・GGA

Ti

Metalliclayer

A = (Si, Ge, Sn, Al, Ga, In)
TiSoftlayer

Ceramics CTilayer

C

Metallic
layerA = (Si, Ge, Sn, Al, Ga, In)

Ti
Ceramics Clayer
Ti
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPANTi

differences in the local elastic constantsfor different A25 /27

MAX phases : Local elastic constants[GPa]

Partly, R. Hossain, H. Kimizuka, Y. Shiihara, andS. Ogata, Comput. Mater. Sci., 209, 111366 (2022).

A

layer

C11C12C13C33

C44

A

layer

C11C12C13C33

Al

soft

320

4246292 127

Si

soft

310

5786332 167

hard

440

7876355 148

hard

460100 148 474 195

soft

297

5154279 116soft

313

49

65

327 143

hard

420

76

27380 138

hard

395

60

79404 151

soft

300

55

45

227

soft

284

67

55

262 108

hard

398

82

62362 137

hard

360

7486412 149

GaIn

soft

97

GeSn

hard

soft

C44

Ti-A
Ti-C
Ti-C

hard

14 Gr.
13 Gr.

We obtained local elasitic constants correlating bond property.
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN26 /27

Summary and future plans
Summary• The formulation of atomic stress was shown and discussed.
• Orbital decomposition can avoid the non-uniqueness problem in atomic stress.
• In multicomponent random alloys, reasonable atomic stress was obtained.
• The local elasticity calculation is demonstrated in the case of MAX phases.
Future plans• Prepare papers describing the formulation(collaborative work with Prof. Shimada in Kyoto Univ.)• Publish my subroutines in Github as unofficial function of OpenMX.
• Revise the formulation based on pair interaction.

Thank you for your attention!
TOYOTA TECHNOLOGICAL INSTITUTE, NAGOYA, JAPAN27 /27