CMD-WS

STATE-Senri

STATE-Senri
Ikutaro Hamada, Kouji Inagaki,
Yuji Hamamoto,
Hidetoshi Kizaki, and Yoshitada Morikawa
Department of Precision Science and Technology,
Osaka University
morikawa@prec.eng.osaka-u.ac.jp
Osaka University

CMD-WS

STATE-Senri

STATE-Senri(Simulation Tool for Atom TEchnology)• Density Functional Theory
LDA8GGA8GGA+U8vdW-DF• Ultrasoft pseudopotential• Plane wave basis set• Iterative diagonalization
Davidson method, RMM-DIIS method• Broyden charge density mixing• Fortran90+MPI+OpenMP: Very good parallelization
efficiency.
--> Applicable to wide range of materials.
Osaka University

STATE-Senri

STATE-Senri(Simulation Tool for Atom TEchnology)• Structure Optimization• Finite Temperature Molecular Dynamics• Reaction Path Search• ESM to calculate electro-cheimcal
reactions and electric field effects.
• Vibrational Spectra (IR, HREELS)8Electronic States (XPS, UPS)8STM
Image calculations
Osaka University

CMD-WS

CMD-WS

STATE-Senri

Pseudopotential Method• Core ElectronsQ Tightly bound to nuclei--> Localized close to nuclei.
• Valence ElectronsQ Weakly bound to nuclei--> Delocalized around neighboring atoms.
• Physical and chemical properties of materials
should be mainly governed by the valence
electron states.
--> Necessary to reproduce valence electronstates accurately.
Osaka University

CMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Concept of Pseudopotential

Advantages of the pseudopotential method• AtomQScattering potential of electrons.
• Core electrons are localized within spheres of radius rc• If the scattering properties of pseudopotentials outside ofradius rcare the same as those of true atomic potentials,
then the valence electron states should be reproduced bypseudopotentials exactly.

• It is as accurate as an all-electron full-potential
linearized augumented plane wave (FLAPW)
method for valence states.
• It is computationally fast because it does not
calculate core states explicitly.
• Required number of effective digits can be smaller
compared with all-electron methods because it
does not treat core states of which energies havelarge absolute values.
• Pseudopotentials and pseudowave functions aresmooth and can be expanded by plane wave basis
sets efficiently.

Electrons

Electrons
Pseudopotentials which haveexactly the same scattering
properties for valence electrons.

True atomic potentials
Osaka University

CMD-WS

STATE-Senri

Osaka University

CMD-WS

STATE-Senri

Norm-conserving pseudopotentials

Costruction of pseudo wave function

To construct norm-conserving pseudopotentials, we first
perform all-electron calculation for an isolated atom andobtain energy eigenvalues εl, wave functions ψ l (r), and,self-consistent effective potential V AE ( r). l is an angularmomentum quantum number of a valence state.
€satisfies the following Kohn-Sham equation.

Next, we construct pseudowave function ϕ l ( r) by smoothingtrue wave function ψ l (r) inside the cutoff radius of rcl.
Usually, we expand ϕ l ( r) by polynomial functions inside the rcl .

(T€+ V AE (r€))ψ l (r)€= εlψ l (r),T =−1 d 2 d l( l + 1)
r+,
2r 2 dr dr
2r 2
∞

2

ψ l ψ l = ∫ 0 ψ l ( r) r 2 dr = 1.
Osaka University€€€ for r ≥ r

ψ l ( r)ϕ l ( r) = M

∑c r

cl

€
for r ≤ rcl

2t +l

l,2t

t= 0

Two more constraints are imposed on ϕ l ( r.)• They are nodeless inside rcl.
• The€norm-conserving condition:
2
r
2 2
€r
ψ
r
r
dr
=
ϕ
r
∫0 l ( )∫ 0 l ( ) r 2dr.
cl

€

cl

The norm-conserving€ condition is very important toensure the accuracy of pseudopotentials.
Osaka University€

CMD-WS

STATE-Senri

CMD-WS

Construction of pseudopotentials

Norm-conserving Pseudopotential•

The norm-conserving pseudopotentials are constructedby imposing the following three conditions.
1. They reproduce correct eigenenergies ε l for valence
electrons.
2. Pseudo-wave functions coincide with all-electron wave
functions outside the cutoff radii rcl.
3. Pseudo-wave functions have correct norm inside thecutoff radii rcl.

By solving the Kohn-Sham equation inversely, pseudopotential Vl(r) which acts pseudowave function ϕ l ( r).
(ε − T )ϕ l ( r) .
Vl ( r) = l
ϕ l (r)€ dependent andPseudopotentials become angular momentum
therefore, they are non-local. Usually, angular momentum
dependence are taken into account only up to l=2 and pseudopotentials for€l>2 are replaced by a common potential Vloc(r).

rc

4π ∫ ψ

l= l max ,m= lpsV =∑ Y (V (r) − V (r)) Y
lm

l

loc

lm

0

+ Vloc ( r)l,m= 0Osaka€ University
STATE-Senri

The norm-conserving condition
The norm-conserving condition can be related to logarithmic
derivatives of wave functions by the following equation.
2
2
rc
rc
2 d d
1
2
ψ
r
r
dr
=
φ
r
(
)
(
)
∫0 l
∫ 0 l r 2dr = − 2 (rφl (r)) dε dr ln φl (r)rc
The logarithmic derivatives of wave functions can be related
to phase shifts of scattered wave functions by sphericalpotentials.
k j l ' ( krc ) cosδ l − n l ' ( krc ) sin δ ldln φ l ( r) =,drj l ( krc ) cosδ l − n l ( krc ) sin δ lrc[

δ l (ε )Incident wave
€Osaka University

2

rc

ps
l

2

(r) r dr = 4 π ∫ 0 ψ ( r) r 2 dr
2

2

d d
lnψ lAEdε drrc

Osaka University

CMD-WSk = 2ε ,AEl= −2π ( rψ lAE )where, |Ylm> are projectors onto states with angularquantum numbers of l and m momentum.

€

STATE-Senri]d

CMD-WS€

STATE-Senri

The norm-conserving condition
A pseudo-wave function ϕ l ( r) is constructed in such a
way that it has the correct scattering property at the
eigen energy el outside the cutoff radius rcl. By€
imposing the norm-conserving
condition, thescattering property becomes correct not only at the
eigen energy el but also at other energies e up to thefirst order of e-el. Therefore, the norm-conserving
condition is very important to ensure the accuracy ofpseudopotentials.

Scattered wave
Osaka University

CMD-WS

STATE-Senri

CMD-WS

Norm-conserving pseudopotential

Norm-conserving pseudopotential	 • The logarithmic derivative of pseudo-wave function agrees
with that of true wave function in wide range of energy.

Si 3s
ψ (r)

STATE-Senri

	
 

.
	
 

y/%1

-

) .	
 

7

. .

	
 

rc

. (- .	
 1 hartree = 27.211eV
Osaka University

Osaka University

CMD-WS

STATE-Senri)

0 8 9E

CMD-WS5 B 88εl
STATE-Senri

Plane Wave Basis Set

Cutoff Wavenumber (Cutoff Energy)• Pseudopotentials and pseudo-wave functions can be
expanded by plane waves efficiently because they
are quite smooth.

The number of reciprocal points G used in plane wave expansion is,in principle, infinite. However, in practice, it should be finite and thesummation with respect to G should be cutoff. Usually, we cutoff G
in such a way that the norm of k+G is smaller than a cutoff wave
number Gmax.

k + G ≤G max

ψ i,k =

k +G ≤Gmax∑ c i,k + G k + G ,ψ ik (r ) =G

r k+G =

k + G k'+G' = δkk'δ GG'
Osaka UniversityG1
exp(i(k + G) • r ),ΩaΩa : Unit cell volume∑€

c ik +G

1
exp(i(k + G) • r ),Ω

cik is a Fourier transofrm of yik(r).
1c ik +G =∫ exp(−i(k + G) • r )ψ ik (r )dr .
Ω Ω
2p/Gmax corresponds to the spacial resolution of the basis set. Therefore,
the required Gmax depends on the steepness of pseudo-wave functions.
Gmax2 is called “Cutoff Energy”.
€
Osaka University

CMD-WS

STATE-Senri

CMD-WS

AB5 B 5.

)

Advantages of Plane Wave Basis Set

)
(

(

(
1 5 5.

((

35BB 68

0 8 9E 25 B 88

Cutoff Wavenumber: Example Si((
1 5 5.

STATE-Senri((

By using Gmax=3.5 aB-1, the total energy is converged
within 0.1 eV and the lattice constant is within 0.001
nm. In general, energy differences among different
structures converge better than the absolute values of
total energies.

• The size of a plane wave basis set can be enlarged
systematically by increasing the cutoff Gmax.
--> No Over Completeness Problem!
• The whole space can be covered with the same accuracy.
--> No Basis Set Superposition Error!
• Atomic forces can be easily calculated because of theHellmann-Feynman theorem for the derivative of thetotal energy with respect to atomic positions.
FI= −

d
ET
dRI

= −∑ f ik ψ ik
ik

1 Hartree = 27.2116 eV∂d(Vlocion + VNL ) ψ ik − dR E EW
∂RII• Very efficient FFT can be used.

Osaka University

Osaka University€
CMD-WS

STATE-Senri

CMD-WSUltrasoft Pseudopotential

Problem in Norm Conserving Pseudopotentials

4A8C7
y

• Nodeless wave functions such as O2p and Cu 3d.
àWave functions have large amplitude close to nuclei.
àTo satisfy the norm-conserving condition, pseudo-wave
functions should also have large amplitude near nuclei.
àThe required number of plane waves to expand steep
pseudo-wave functions becomes large and calculations
become very heavy.

STATE-Senri.
• By removing the norm-conserving condition,
pseudo-wave functions can be very smooth.
à Ultrasoft Pseudopotential
Osaka University

Osaka University

CMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Generalized Norm-conserving Condtion
4π ∫ ψ
0

AE
l

2

Cu 3d

(r) r 2 dr = ψ lAE ψ lAE r = ψ lPS S ψ lPS r= −2π ( rψ lPS )

c

2

c

d d
lnψ lPS ,dε drrcGeneralized Norm-conserving Condition

S = 1+ ∑ qnm β n β m ,Y'(r)/ Y(r)rc

Ultrasoft Pseudopotential

ψ lPS β m = δlm .

)nm

H ψ lPS = εl S ψ lPS .

0 8 9E 25 B 88
Generalized Eigen Value Problem

Osaka University

The accuracy of pseudopotentials can be improved
further by employing multiple reference states.
Osaka University

CMD-WS

STATE-Senri

3E-trilemma

CMD-WS

Important roles of Interfaces!Economy!Energy!Environment # "  "

Solve the trilemma by Quantum Simulations.
Osaka University

STATE-Senri

2018/09/18-19

CMD-WS

Heterogeneous Catalyst

STATE-Senri

CMD-WS

STATE-Senri

Heterogeneous Catalyst

H2O dimer
H-bondNO dissociation by H2O
Thanh Ngoc Pham et al.,JPhysChemC.(2018).

HCO
CO2 hydrogenation on Cu surfaces
Fahdzi Muttaqien, et al., J. Chem. Phys. (2014),
J. Chem. Phys. (2017) , Chem. Comm. (2017),
Nature Chem. (2019).

Etching of SiC by Pt catalyst
B.V. Pho et al.,Appl. Phys. Lett., (2015).
Suppression of Pd Sintering
I. Hamada et al., J. Am. Chem. Soc.,(2011).

4

4

7

6

8

8

Suppression of Pd Sintering
I. Hamada et al., J. Am. Chem. Soc.,(2011).
99HCO
CO2 hydrogenation/decomposition on Cu
G.C. Wang et al.,J. Phys. Chem. B, (2006).
Fahdzi Muttaqien et al., Chem. Comm. (2017).

Etching of SiC by Pt catalyst
Bui Van Pho et al.,Appl. Phys. Lett., (2015)2019/08/26 J.J.Appl.Phys (2018).

H2O dissociation on MgO/Ag,
H.-J. Shin et al., Nature Materials (2010).

Osaka University

Chemicalat Water/Solid
Chemical Reactions
Reactions at
Water/Solid Interfaces
Interface

CMD-WS

STATE-Senri

Methanol Oxidation on Pt ElectrodeJ. Herron et al., PNAS USA (2016).
Hydrogen Evolution Reaction
M. Otani et al., J. Phys. Soc. Jpn (2008).
PCCP (2008), JPCC (2008).

H2O dissociation on MgO/Ag,
H.-J. Shin et al., Nature Materials (2010).

Osaka University

CO2 adsorption/dissociation on Cu surfaces
Fahdzi Muttaqien, et al.,J. Chem. Phys. (2014), J. Chem. Phys. (2017).

Water Dimer / Cu(110)
T. Kumagai et al.,Phys. Rev. Lett., (2008),
Phys. Rev. B (2009).

CMD-WS

STATE-Senri

Input file: nfinp_1• Ictrl : Should be always “ 0 0 0 0 0 0 “• GMAXThe cutoff wave number of plane wave basis set for wave functions.
• GMAXPThe cutoff wave number of plane wave basis set for charge density.
These cutoff numbers depend on pseudopotentials.
ExampleQ
pot.Al_pbe1, pot.Si_pbe1, pot.S_pbe1, pot.Cl_pbe1, pot.Br_pbe1,
pot.Ge_pbe1,
GMAX = 3R4, GMAXP = 10
pot.H_lda3, pot.C_pbe3, pot.O_pbe3,
GMAX = 5.0 , GMAXP = 15.0
pot.H_lda1, pot.C_pbe1, pot.O_pbe1,
GMAX = 5.5, GMAXP = 20.0
It isUniversity
necessary to check the dependence of total energies on cutoff.
OsakaCMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Input file: nfinp_1

Input file: nfinp_1• Ictrl : Should be always “ 0 0 0 0 0 0 “• GMAXThe cutoff wave number of plane wave basis set for wave functions.
• GMAXPThe cutoff wave number of plane wave basis set for charge density.
These cutoff numbers depend on pseudopotentials.
ExampleQ
pot.Al_pbe1, pot.Si_pbe1, pot.S_pbe1, pot.Cl_pbe1, pot.Br_pbe1,
pot.Ge_pbe1,
GMAX = 3R4, GMAXP = 10
pot.H_lda3, pot.C_pbe3, pot.O_pbe3,
GMAX = 5.0 , GMAXP = 15.0
pot.H_lda1, pot.C_pbe1, pot.O_pbe1,
GMAX = 5.5, GMAXP = 20.0
It isUniversity
necessary to check the dependence of total energies on cutoff.
Osaka

Osaka University

CMD-WS

CMD-WS

STATE-Senri• NTYPE : Number of atomic types in a unit cell
• NATM = NATM2
The number of atoms in a unit cell
• NUM_SPACE_GROUP
Specify the symmetry of the system by the group space number.
• TYPE_OF_BRAVAIS_LATTICE
Specify the type of Bravais lattice.
0 SIMPLE
1 BODY_CENTER
2 FACE_CENTER
3 A_CENTER
4 B_CENTER
5 C_CENTER
6 RHOMBOHEDRAL_CENTER

STATE-Senri

Unit Cell

K-points

ONSpecify the lengths of the lattice vectors and the angles among
them.
a

b

c

a

b

g

kx

ky kz 1ca ba
PNSpecify the Cartesian coordinates of the lattice vectors. g b
Cartesian
ax ay
bx by
cx cy

Osaka University

az
bz
cz

Osaka University

1

1 : KNX, KNY, KNZ, K-POINT_SHIFTCMD-WS

STATE-Senri

CMD-WS

Atomic coordinates

Input file: nfinp_1

NCORD: Unit of coordinates
0
Unit of Primitive Cell
1
Cartesian coordinate in Bohr radius
2
Unit of Conventional Cell
NINV: Obsolete
x
x
x

y
y
y

z
z
z

STATE-Senri

•
•
•
•
•
•

IATOMN : Atomic Number of atom type 1
ALFA: Obsolete
AMION: Atomic Mass
ILOC: Obsolete
IVAN: Obsolete
ZETA1: Initial spin polarization of atom

IWEI IMDTYP ITYP
IWEI IMDTYP ITYP
IWEI IMDTYP ITYP…
IWEI : Obsolete and should be always 1
IMDTYP = 0: Fix atomic position1: Move atomic position
ITYPE: Atom type number
Osaka University

Osaka University

CMD-WS

STATE-Senri

CMD-WS

Input file: nfinp_1• ICOND
How to prepare the initial wave functions.
0 Q Generate from scratch
1 Q Read from “restart.data” which has the latest date.
• INIPOS
How to prepare the initial atomic positions.
0 Q Read from nfinp_1.
1 Q Read from “restart.data” which has the latet data.
2 : Read from “GEOMETRY” which can be modified by text
editor.
• INIVEL
Specify the atom velocities actually used;
0 à Initialized, 1à Read from restart file, 2à Read fromGEOMETRY file
Osaka University

STATE-Senri

Input file: nfinp_1• ININOS
Specify the thermostat variables actually used;
0 Q Initialized
1 Q Read from “restart.data” which has the latest date.
• INIACC
Reset accumulators and the counter of MD steps.
0 Q Initialize.
1 Q Read from “restart.data” and continue accumulating.
• NMD1
Maximum number of iterations for SCF calculations.
• NMD2
Maximum number of molecular dynamics steps.
• CPUMAX
Maximum CPU time in seconds.
Osaka University

CMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Input file: nfinp_1

Structure Optimization and Molecular Dynamics• Way-of-mix
1 Simple
2 BROYD1
3 BROYD2
4 DFP
5 PULAY
6 BLUGEL• Mix-what: 1 Charge Mixing, 2 Potential Mixing• KBXMIX
The maximum total number of iterations to be used for mixing.
• MIX_ALPHA
Fraction of new charge density in charge density mixing. 0.3R0.8.

DTIO:• Increment of time for each step of molecular dynamics simulation.
(O a.u. =2.4x10-17seconds)• The time step should be carefully chosen to assure the accuracy ofmolecular dynamics. It depends on atomic types, maximum forces
acting on atoms, and integration scheme of the Newton equation of
motion.
• For QUENCHING and NEB schemes, the time step should be 50R
300a.u.
• For GDIIS scheme, the time step should be 400R800a.u.
• For NORMAL_MODE, it is not used.
• For FINITE_TEMPERATURE_MD and if the system includes
hydrogen, the time step should be about 20a.u. If hydrogen is notincluded, the time step should be 50R100a.u.

Osaka University

Osaka University

CMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Structure Optimization and Molecular Dynamics

Input file: nfinp_1

IMDALG: Algorithms of MD
1 Newtonian equations of motion
2 Quenching Molecular Dynamics
This scheme works well when forces acting on atoms arelarger than 0.02 hartree/aB.
3 Vibration normal mode analysis.
4 GDIIS: Generalized Direct Inversion of Iterative
Subspace
This scheme becomes efficient when the maximum force
is smaller than about 0.02 hartree/aB.
6 Nudged Elastic Band Method
Method for reaction path search proposed by Mills,
Jonsson, Schenter, Surf. Sci. 324 305 (1994).

• EDELTA: Energy criteria for convergence [ in Hartree/atom]
Recommended values: 0.1d-08  0.1d-09• WIDTH: Scheme of Brillouin Zone Sampling
WIDTH < -10.0 Tetrahedral scheme-10.0 < WIDTH < 10.0 Special k-points scheme is used andthe absolute value of WIDTH, namely, |WIDTH| means the energy
broadening around the Fermi level• FORCCR: Criterion for force in MD [Hartree/Bohr radius]
If the maximum force is smaller than FORCCR, the molecular
dynamics is terminated.
Recommended value: 1.d-03 [Hartree/Borh radius]• XCTYPE
Selection of exchange-correlation function.
GGAPBE QGeneralized Gradient Approximation by Perdew, Burkeand Ernzerhof, Phys. Rev. Lett. 77 3865 (1996). (Recommended)
LDAPW91 : Local Density Approximation.

Osaka University

Osaka University

CMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Input file: nfinp_1

Input file: nfinp_1• NSPIN
0 QNo spin polarization.
1 : Spin polarized calculation.
• DESTM
Energy width for STM-imaging.
• KEG
Number of eigenvalues.

• IMSD
Choose optimization method.
2 QDavidson method.
1 : Residual-vector Minimization Method.
Davidson scheme is quite stable but computationally less efficientcompared with RMM scheme. RMM is very fast but it is less stablecompared with Davidson scheme. When you start calculations fromscratch wave functions, you should use Davidson scheme for first 510 iterations and switch sot RMM scheme.
• EDELTA
Convergence criteria for wave function optimization.
Recommended values: 0.10D-7R0.10D-8

Osaka University

Osaka University

CMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Reaction Path Search&JOB_TYPE
NEB&END_JOB_TYPE
NEB: Nudged Elastic Band Method
Mills, Jonsson, Schenter, Surf. Sci. 324 305 (1994).

Reaction Path Search&JOB_TYPE&END_JOB_TYPE

NEB

NEB: Nudged Elastic Band Method
Mills, Jonsson, Schenter, Surf. Sci. 324 305 (1994).

Saddle Point
MTransition StateN?DIAL<FCKBGI;Osaka University

Osaka University

;

:

CMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Reaction Path Search&JOB_TYPE&END_JOB_TYPE

Reaction Path Search

NEB&JOB_TYPE&END_JOB_TYPE

NEB: Nudged Elastic Band Method
Mills, Jonsson, Schenter, Surf. Sci. 324 305 (1994).

Osaka University

NEB: Nudged Elastic Band Method
Mills, Jonsson, Schenter, Surf. Sci. 324 305 (1994).

Osaka University

CMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Reaction Path Search&JOB_TYPE&END_JOB_TYPE

NEB

Reaction Path Search

NEB

NEB: Nudged Elastic Band Method
Mills, Jonsson, Schenter, Surf. Sci. 324 305 (1994).

Rn+1

7+!.!.!.⊥Rn
F∥=M-∇VN∥

Rn-1
JEH@>;JEH@>=Osaka University

;



9



JEH@>;JEH@>=Osaka University

;



9



CMD-WS

STATE-Senri

CMD-WS

Untar state files

STATE-Senri

Submitting a job• Copy a state file to your directory.

• Submit a job

cp ~localuser/cmd/morikawa/state-senri.tar.gz .

> mpirun –np 1 ./STATE < nfinp_1 > nfout_1 &• Plot density of states (DOS)• Unzip and untar state files.

> gnuplot>> plot “dos.data” u 1:2 w l↑“el”> gzip –d state-senri.tar.gz> tar xvf state-senri.tar .

Osaka University

Osaka University

CMD-WS

STATE-Senri

CMD-WS

STATE-Senri

Cl on Al(100)Submitting a joby• Plot density of states (DOS)T> gnuplot>> plot “dos.data” u 1:2 w l , “” u 1:3 w l
↑
“el”

B
H

Adsorption Sites
T: Top
B: Bridge
H: Hollow
xFirst-layer Al
Second-layer Al

Osaka University

Osaka University

CMD-WS

STATE-Senri

CMD-WS

Cl on Al(100)yB

B

STATE-Senri

Loging to cmd2

Adsorption Sites
T: Top
B: Bridge
B H: Hollowx>>ssh –Y astud??@cmd2.phys.sci.osaka-u.ac.jp

Account ID
astud??
Password
cmdaastud??

?? = 01

24

First-layer Al
B

B

B
Second-layer Al

Osaka University

Osaka University

CMD-WS

STATE-Senri

copy state-senri tar file
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>
>>

cp /home/CMD/teac04/state.tar.gz .
tar xvfz state.tar.gz
cd STATEls –lexample: example files
external: external libraries
gncpp: pseudo-potential database
README: README filesrc
: source codes
test: test files
tools: tool programs

Osaka University

CMD-WS

STATE-Senri

Installation of mpich and libfftw3-dev24%/ "03(&3*.23",,-0*$)24%/ "03(&3*.23",,,*#''36%&5Osaka University