Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

3.2

First-Principles Calculationof Material Properties

41

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Clarification of Microscopic Mechanisms of Semiconductor Epitaxial
Growth and Device-Interface Formation by Large-Scale Quantum-TheoryBased Computations
Atsushi Oshiyama
Institute of Materials and Systems for Sustainability, Nagoya University
Furo-cho, Chikusa-ku, Nagoya 464-8601
In the fiscal year of 2021, on the basis of theGaN, a premier material in optoelectronics is

total-energy electronic-structure calculations andnow emerging in power electronics due to its

molecular

superior

dynamics

simulations

within thephysicalproperties.

However,todensity-functional theory (DFT), we have studied

guarantee the expected high performance of thethe epitaxial growth of wide-gap power semi-

devices, it is necessary to forge high-quality

conductors, atomic and electronic structures of

epitaxial

semiconductor-insulator interfaces. The main

experimental efforts in this decade, the current

computational tools are our RSDFT (Real Space

quality of GaN crystalline films is unsatisfactoryDFT) code and RS-CPMD (Car-Parrinello Mole-

for replacing Si power devices in the market. To

cular Dynamics) code as well as VASP code.

unveilSpecifically, we have studied 1) the atom-scale

identification of the surfaces on which atomic

mechanism of GaN epitaxial growth [1], 2) thereactions take place is an essential prerequisite.

films

the

of

growthGaN.

Despitephenomena,

many

atom-scale

mechanism of the adatom diffusion on SiC
stepped surface [2], 3) the formation of danglingbond free interface of GaN and a gate insulator[3], and 4) the structural identification of thenano-meter-scale GaN growth surface using
Bayesian

optimization

[4].

We

have

also

developed 5) an order-N orbital-free DFT (OFDFT)

schemebyinventingFig.1: Top view of the stable 6x6 surface

kinetic-energy

structurefunctional using the neural-network [5]. The

discovered

by

the

present

Bayesian optimization (structure name#130 in Fig. 2). The black dashed line

below is the explanation of the issues 4) and 5)

indicates a 6x6 area. Brown, blue andabove.

white atoms correspond to Ga, N and H,
respectively.

Structural identification of GaN (0001) byBayesian optimization

42

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

We have identified the stable structure ofand more complex adsorbate-arrangements, andnanometer-scale GaN (0001) surfaces with Ga-then determination of surface phase diagrams for

and H-adsorbates, which are the fundamental

design of further preferred growth conditions.

basis for crystal growth modeling, by large-scale

Order-N OFDFT scheme with machine

density-functional calculations combined withlearning

machinelearning

Bayesian

optimization

A scheme in DFT in which the kinetic energy

technique. We have been able to reach a single

functional is expressed by the sum of single-

stable structure (Fig. 1) with satisfactorily low

particle orbitals squared is introduced by Kohn

mixing enthalpy by 130 trials based on Bayesianand Sham (KSDFT) and has been applied to a

optimization among

numerous

23C15=490,314

candidatenumberofapplications

withunprecedented success.

structures (Fig. 2). We have found that the

However, solving the Kohn-Sham equations for

obtained structure lacks any postulated highby humanall the occupied orbitals in the system is a

intuition, satisfies electron counting rule locally,

computational burden scaling with the system sizeand shows a complex adsorbate-rearrangement,

N as O(N3). The scheme with lower-order scalingleading to lower mixing enthalpy. The present

is highly demanded in materials science and also

scheme of Bayesian optimization combined within advancing DFT. One of the solutions in a

first-principle calculations paves a way toward

legitimate

identifying surface structures with a larger scale

functional theory (OFDFT) in which the kinetic

symmetry previously introduced

way is the orbital-free density-

energy functional is expressed as a functional of
the electron density ρ (r ) itself and an single
Euler equation rather than Kohn-Sham equations
is solved, being expected to be O(N) scheme.
We express the kinetic energy functional asT [ ρ ] = ∫ τ TF (r ) F ( s 2 , q; W )dr ,where τ TF is the Thomas-Fermi kinetic energy
density, the enhancement factor F is expressed as
Fig. 2: History of the exploration of thea gradient expansion with s = ∇ρ [2(3π 2 )1/3 ρ 4/3 ]surface structure by Bayesian optimization.
The objective is to minimize the mixing

and q = ∇ 2 ρ [4(3π 2 )2/3 ρ 5/3 ] . The enhancement

enthalpy.

factor is determined by the deep learning

The

black

dashedline

corresponds to the minimum value of theusing the neural network where the set ofbaseline models, and the red marker
corresponds to the most stable structure

weight parameters W is determined so as#130 in the present exploration.

to minimize the cost function L expressed43Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyoin terms of the kinetic energy functional

systems which include semiconductors,

derivative δ T δρ :metals, ionic insulators and some of their2=Lpolytypes: diamond, graphene, ds-Si, face-
1 Nt 1  δ T
δ T KS()(rn )  ,−r∑n
Nt n=1 2  δρδρcentered-cubic(fcc)-Si,

beta-tin

Si,

zinc-

functional

blende (3C)-SiC, body-centered-cubic(bcc)-

derivative is computed at each real-space

Li, fcc-Al, fcc-Cu, bcc-Na, and NaCl. The

mesh

the

structural properties, i.e., lattice constants

corresponding values by the Kohn-Shamand bulk moduli of all those materials

scheme. As a training data, we use theobtained by the present functional shows

functional

better agreement with the corresponding

where

the

pointskinetic-energy

andcomparedderivative

in

with

diamond-values by the Kohn Sham scheme.

structured (ds) carbon at some mesh points.

The order N behavior of the computation
time in the present OFDFT scheme isclearly shown in Fig. 3.

The group id codes for the achievements
above are k0042 and m0037.

Fig. 3: Computational time of KSDFT and

Related publications

OFDFT. Comparison of (a) t1 (subspace[1] K. M. Bui, K. Shiraishi and A. Oshiyama,diagonalization), t2 (conjugate gradientAppl. Surf. Sci. 557, 149542 (2021).

minimization and gram-schmidt diagonalization), and t3 (SCF mixing), and (b) total[2] K. Seino and A. Oshiyama, Appl. Surf. Sci.

time for the SCF iteration. The test

561, 149927 (2021).

system is 4H-SiC supercell with thenumber of atoms of 576 -4704. In both (a)[3] K. Chokawa, K. Shiraishi and A. Oshiyama,and (b), data points obtained by KSDFT

Appl. Phys. Lett. 119, 011602 (2021).

and OFDFT are plotted in red and blue

[4] A. Kusaba, Y. Kangawa, T. Kuboyama and A.

colors, respectively.

Oshiyama, Appl. Phys. Lett. 120, 021602 (2022).

The transferability of the kinetic-energy[5] F. Imoto, M. Imada adn A. Oshiyama, Phys.

functional obtained by the deep learning

Rev. Research 3, 033198 (2021).

with the training data of only diamond has
been examined by the application to 23 test

44

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles prediction of stability and functionality ofcomplex materials
Osamu SUGINO
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
In recent years, considerable progress has

systems will also be effective not only forbeen made in predicting material properties

electronic systems but also the hybrid solvation

using the density functional theory (DFT)models that we have long been studying.

simulations

and

the

DFT-basedmodel

In view of ability of meta-GGA in describing

simulations. In our group, we have developed a

the Mott insulators, this functional was used for

scheme to construct the exchange-correlation

our study on exotic materials such as the solid(xc) functional of the Kohn-Sham DFT byoxygens in the high-pressure phases and the Hg-machine-learning the many-body wave function

cuprate (HgBa2Can−1CunO2n+2+x) [4]. We foundof small molecules [1,2]. In our prior study [1],that the lattice parameters and the bandgaps

small molecular systems were simply used as the

become consistent with experiments with thislearning set, but the success was limited tofunctional, indicating that we are now at the startmolecular systems only; problems remain inline of quantitative research.

crystalline systems. To overcome this problem,

One of the open problems of DFT-based

physical conditions applicable to general

simulation is to incorporate the nuclear quantum

materials were introduced into the learning

effect. Although the path integral molecular

process, and the resulting xc functional within

dynamics (PIMD) simulation is an establishedthe meta-GGA format was found to significantly

first-principles scheme, it becomes difficult toimprove the applicability [2]. The improved xc

reducefunctional was implemented to an electronic

temperature is reduced and thus different

structure calculation package VASP and is readyscheme is needed. To overcome this problem forto distribute to supercomputer users. For further

the hydrogen diffusion on metal surfaces, we

progress, we focused on techniques for the DFT

have reformulated the quantum transition state

for classical systems (classical DFT) and, in thistheory that is based on a semi-classicalcase, low-density approximation of the n-

approximation using an apparently full-quantum

particle density was found to be effective [3]. We

formulation. The fluctuation effect was found to

consider that applying this method to electronic

be significantly reduced with the new method

45

the

statisticalfluctuation

as

the

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[5].

reduction reaction occurring thereon andWe have also investigated the effect ofconsidered the effect of liquid water as well.

nuclear quantum motion on the bandgap ofThese results were shared with experimentalists

insulators. This was done using the Allen-Heine-and helped them in planning the electrode

Cardona (AHC) theory and its extensions [6].

development.

Herein the Debye-Waller terms of AHC wasreformulated using the momentum operator, by

References

which some of the electron-phonon coupling[1] R. Nagai et al.: Phys. Rev. Res. 4 (2022)

terms can be added to infinite order.

013106.

The description of electrochemical interfaces[2] R. Nagai et al.: npj Compt. Mater. 6 (2020)is a challenging problem of first-principles

1.

simulations, on which our group has long been

[3] T. Yokota et al.: Phys. Rev. E 104 (2021)

studying. The difficulty of the simulation is due

014124.

to existence the electric dipole layer (EDL) that[4] A. Tatan et al.: arXiv cond-mat 2204.05719

can extend macroscopically. The thickness of(2022).

EDL, however, is reduced to several nm when[5] Y. Kataoka et al: unpublished.

the ions are densely distributed as in strong[6] K. Ishii et al.: Phys. Rev. B 104 (2021)

acidic/alkaline solutions, but methods for the245144.

simulation has been developed only for the[7]

solution side. Therefore, EDL existing in thecenter-dev/abICS

solid side of the interface remains to be achallenge; this is the case when semiconducting
electrode is used. To simulate the ZrOxN2-xelectrode, which is a candidate for Pt-freeelectrode of the next-generation fuel-cells, wehave applied the first-principles Monte Carlosimulation with abICS [7] and then sped it up
with the help of a machine-learned model
potential.

By

the

simulation,

we

could

successfully obtain the distribution of oxygen
vacancies and the nitrogen impurity atoms, andthe result was used to describe the EDL. We
further investigated efficiency of the oxygen

46

S.

Kasamatsu:

https://github.com/issp-

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Analyses related to atomic structures and atom
dynamics at complex structures such as surfaces,interfaces and defects
Satoshi WATANABE
Department of Materials Engineering, the University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-8656

1

Introduction

Keeping in mind the importance and diﬃculty of atomistic simulations on structuresand dynamics at surfaces, interfaces and defects toward development of novel informationand energy devices, we have been investigating
atomic structures and atom dynamics in various systems by first-principles calculations, insome cases combining them with multi-scale ormachine-learning approaches. In the following,
two of our recent results are described.

2

Band alignment of LiTi2 O4and SrTiO3Recently, it was found in thin-film solid-state
Li batteries using LiTi2 O4 thin films as positive electrodes, the spontaneous Li-ion migration between the solid electrolyte (Li3 PO4 ) andLiTi2 O4 occurs in the case of the conducting
substrate, Nb-doped SrTiO3 (001) [1]. Thatis, the LiTi2 O4 turned to the Li2 Ti2 O4 phasespontaneously.
To understand the above, we have investigated the band alignment of the
LiTi2 O4 /SrTiO3 superlattice model using
the density functional theory (DFT) calculations. Here, we adopted the SrTiO3 without
Nd doping to examine the band alignment
before the charge transfer. The calculated
projected density of states (pDOS) reveals
that the bottom of the conduction band islocated above the Fermi level EF in the central
region of the SrTiO3 layer. Since EF of the
Nb-doped SrTiO3 is considered to be locatednear the bottom of the conduction band,
this result suggests the electron transfer fromthe Nb-doped SrTiO3 to LiTi2 O4 . From thisresult and the absence of a potential gradientin the metallic LiTi2 O4 layer, we can saythat the transferred electrons charge on thesurface of LiTi2 O4 , which leads to a decreasein the electrochemical potential of Li ionsat the LiTi2 O4 surface. Consequently, thespontaneous migration of the Li ions occurs.
The present result suggests the possibility tocontrol the spontaneous Li-ion migration between the solid electrolyte and LiTi2 O4 by tuning the electrochemical potential diﬀerence via
the substrate selection.

3

Neural network potential
for multiple charge states ofdefects

In semiconductors and insulators, defects
and impurities sometimes take various charge
states depending on EF . To examine the point
defect behavior in such multiple charge states,
we have proposed a new scheme of neural network potential (NNP) [2]. In our scheme,
only minimal modifications to the conventional
scheme [3] is needed: One system charge (Qsys )node is added in the input layer. Here, Qsys

47

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyois defined as a constant value times a charge
state of the considered supercell divided by its
volume.
We have applied this scheme to wurzite-GaN
with a nitrogen vacancy [2]. For the charge
states of the vacancy, 0, 1+, 2+, and 3+
have been considered. By performing molecular dynamics simulations using the StillingerWeber potential, we have generated 4,550 pristine GaN structures and 4,550 structures witha N vacancy in total. Then the DFTdata havebeen obtained for these sturctures. Total number of the DFT data is 22,750, because calculations have been performed for 0, 1+, 2+ and3+ charge states in the cases of the structureswith a N vacancy.
The NNP constructed using our scheme andthe aforementioned data predict the total energies and atomic forces accurately for all the
charge states: The root-mean-square-errors of
the total energies and atomic forces were, respectively, 1.45 meV/atom and 63.8 meV/Å
for the training dataset and 1.44 meV/atomand 64.6 meV/Å for the test dataset. Moreover, it has reproduced the phonon band structures and themodynamic properties of the defective systems fairly well.

References[1] K. Nishio, S. Ichinokura, A. Nakanishi,
K. Shimizu, Y. Kobayashi, N. Nakamura,
D. Imazeki, R. Shimizu, T. Hirahara, S.
Watanabe, and T. Hitosugi: Nano Lett.
21 (2021) 10086.
[2] K. Shimizu, Y. Dou, E. F. Arguelles, T.
Moriya, E. Minamitani, and S. Watanabe:
arXiv:2203.16789.
[3] J. Behler and M. Parrinello: Phys. Rev.
Lett. 98 (2007) 146401.

48

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Quantum Simulations on Dynamical Heterogeneous
Catalysts
Y. J. WONG, J. I. ENRIQUEZ, H. H. HALIM, S. E. M. PUTRA, T. N. PHAM,
Y. HAMAMOTO, K. INAGAKI, I. HAMADA, and Y. MORIKAWA
Department of Precision Engineering, Graduate School of Engineering,Osaka University, 2-1 Yamadaoka Suita, Osaka 565-0871

In 2021-2022, we carried out theoretical investigation of chemical processes at surfaces
and interfaces, adsorption and reaction of polymeric formic acid on Cu(111), NO adsorption
and dissociation on Cu(100), oxygen vacancy
migration in SrFeO3−δ and Sr3 Fe2 O7−δ , dry reforming of methane on Co catalysts, and adsorption of CO2 on terrace, stepped and defectsites on Pt surfaces.
In this report, we report the dry reformingof methane on Co catalysts,[1]. The emissionsof greenhouse gases (CH4 and CO2 ) are rising
exponentially due to human activities and increased energy demand, contributing to global
warming and climate changes. To make useof CH4 and CO2 into renewable energy, catalytic dry reforming of methane (DRM) is aneﬃcient pathway to convert them into syngas(CO and H2 ). Currently, active transition metals known to be capable of catalyzing DRM aremainly group VIII elements such as Ru, Rh,
Ni, Co and Pt. Among them, supported Ni
could be the most studied catalyst due to its
relatively high activity and low cost. However,
the eﬃciency of Ni is widely known to be easily
suppressed by carbon deposition.
On the other hand, Co is particularly interesting due to its excellent catalytic activity
and considerably low carbon deposition. Experimentally, Co is more carbon resistant thanNi. Co could enhance the adsorption of CO2 ,which helps to inhibit the deposition of carbon
and enhance carbon removal.

To gain insight into the occurrence of carbon
deposition vs removal during DRM, we calculated the CO* formation reaction throughthe surface reaction of C* and O*, and theCC* formation reaction through the C-C coupling. The present direct quantitative comparison between the activation energies of the involved reactions would clarify the competition
between carbon deposition vs removal over different Co surfaces. In addition, we further
compared our results with those of Ni to reveal the cause of the lower carbon deposition
on Co than on Ni. The quantitative comparison of CO2 activation, C-C coupling, and CO*formation over Co and Ni could reveal the theoretical underlying mechanism of lower carbon
deposition on Co than Ni.
The overall energy diagram illustrating theminimum energy pathway for reactions between CO2 and C* on all the investigated Co
surfaces, which includes CO2 activation (green
region), O* and C* diﬀusion (blue region), andthe competition between CO* formation and
C-C coupling (red region), is summarized in
Figure 1. We consider the reaction of CO2
with two predeposited C*. As observed in Figure 1, the rate-determining step for the reaction between CO2 and C* on all the investigated Co surfaces is located in the red region, indicated by the energy needed for CO*
or CC* dimer formation. Two conclusions canbe drawn from evaluating the competition between CO* and CC* dimer formation on each49

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

surface. First, carbon graphitization is dominant on (111) and possibly on (110) because
the CC* dimer ﬁnal state is of lower or equal
energy r elative to its respective initial state,indicating that the CC* dimer is stable anddoes not dissociate into atomic C*. In addition, CC* dimer formation on these surfaces
has a lower activation energy than CO* formation. Even though the CO* ﬁnal state ismore stable than the CC* dimer, the stabilityof the CC* dimer compared to its initial statemight give rise to nonreversible carbon graphitization, which could eventually cover the active surfaces. Thus, on (111) and (110), CC*
dimer formation has a higher degree of preference than CO* formation. By combining thisresult with the fact that (111) is the most exposed surface of a Co nanoparticle, as indicated by the Wulﬀ constructions, we hereby
deduce that the origin of carbon deposition
that leads to catalyst poisoning during DRM
is caused by the high degree of preference for
C-C coupling than CO* formation, especiallyon (111). Second, on (100), (211), and (221),
CC* dimer formation is not favorable because
the CC* dimer ﬁnal state has a higher energy
level than its respective initial state, indicatingthat the CC* dimer is easily dissociated into
two atomic C*s. In addition, the CO* ﬁnal
state is more stable than the CC* dimer. The
combination of the above mentioned points implies a higher degree of preference for CO* formation on (100), (211) and (221), even though
the energy needed for CC* formation might
be lower than that of CO* formation. On theother hand, other surface reactions, such as
CO2 activation and surface diﬀusion of O* andC* have lower activation energies than those of
CO* and CC* dimer formation on each respective surface (Figure 1). The low energy needed
for CO2 activation implies that surface O*s are
suﬃciently supplied for CO* formation. In addition, the kinetic rate constants for CO* formation, CC* formation, and CC* dissociation
over all the investigated surfaces were calcu-

lated and discussed in Section 12, S.I.. The result corresponds well with the ﬁndings where
C-C coupling has the highest degree of preference on (111), followed by (110), however, is
not favorable on (100), (211), and (221).

Figure 1: Energy diagram showing the minimum energy pathways of reactions between
CO2 and two C* on FCC Co (111), (110),(100), (211), and (221).
The present study provides a comprehensive
clariﬁcation on the mechanism of DRM, specifically on the surface reaction of CO2 and C*
on the major ﬂat and step FCC Co surfaces,
namely, (111), (110), (100), (211), and (221).
and/or promoted Co nanoparticles for DRM
to restrict C-C coupling.

References[1] Y.-J. Wong, H. H. Halim; N. F.
Khairudin, T. N. Pham, S. E. M. Putra,
Y. Hamamoto, K. Inagaki, I. Hamada, A.
R. Mohamed, and Y. Morikawa: J. Phys.
Chem. C 125, 21902-21913 (2021).

50

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Understanding superionic conductivity in disordered systems
using machine learning potential molecular dynamics andMonte Carlo sampling
Shusuke KASAMATSU
Academic Assembly (Faculty of Science),
Yamagata University, Kojirakawa, Yamagata-shi, Yamagata 990-8560
The competition between order and disorderenvironments, Sc-VO-Sc, Sc-VO-Zr, or Zr-VO-

in materials leads to a wide variety of materials

Zr, are responsible for the hydration behavior.

functionalities such as electronic, ionic, andUsing abICS, we have clarified that the former

thermal conductivity/insulation and catalytic

two are active sites for hydration, while the Zr-

activity. However, because of the huge number

VO-Zrof degrees of freedom in disordered systems,

thermodynamic conditions and contributes little

successes in theory-driven design of suchto the hydration reaction [3].

site

seldomexistsunder

normalfunctionalities have been rather limited. In thisIn the case of amorphous systems, we

research project, we tackled this problem

considered the use of the neural network

through a combination of first-principles

potential (NNP) approach originally proposed

calculation and machine learning.

by Behler and Parinello [4]. Although the

For disordered crystalline systems, weneural network potential methodology seems toemployed our abICS (ab Initio Configuration

be a well-established approach at this point, we

Sampling) framework for deriving an on-lattice

have found that it is quite a challenge to apply

neural network model from first-principles

to many-component glass systems. This iscalculations [1] and used that to perform

because the NNP requires that the training set

replica-exchange Monte Carlo sampling. One

sufficiently

system that we have been examining is heavily

relevant configuration space, which will be

Sc-doped BaZrO3, which is known as ahuge for many-component glasses. We observe

promising perovskite proton conductor for solid

regularly that NNP molecular dynamics (MD)oxide fuel cells [2]. The oxygen vacancy is

calculations ‘fall apart’ and exhibit completely

known to be the active site for the hydration

unphysical structures after it wonders into

reaction, which is the key reaction to introduce

regions of previously unlearned configuration

protons into the system. However, it is yet

space. To tackle this issue, we have been trying

unclear which of the three types of vacancy

an active learning cycle approach where we51

coversthethermodynamicallyActivity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

take structures from the NNP-MD trajectory,
perform DFT calculations on those structures,
add them to the training set, and run the NNPMD again with the new NNP. As shown in
Fig. 1 for the AgI-As2Se3 superionic mixed
glass system, increasing the number of cycles
leads to improvement in the NNP predictions,although we have not been able to fully

Fig. 1: Improvement of NNP accuracy withconverge the low energy structures yet.

active learning iterations for AgI-As2Se3

We have also been examining the feasibilityof training on DFT data using high-costfunctionals such as the HSE06 hybrid. Since it
is not feasible to obtain a long enough MD
trajectory to use as the training set with HSE06,
we first performed first-principles MD usingtherelatively

low-cost

GGA

functional.

Afterwards, we took a subset of the structures

Fig. 2: Total structure factors from NNPs(~1500 structures) and recalculated them using

trained on various functionals and compared to

HSE06, which were then used as the training

experiment.

set in NNP training. This approach was applied
to amorphous GeO2, and we were able to obtain

References

structure factors in excellent agreement with

[1] S. Kasamatsu, ISSP Supercomputer Center

experiment (Fig. 2). This shows that NNP

Activity Report 2020, pp. 16-26.

technology is promising for realizing long-time[2] J. Hyodo et al.: Adv. Energy. Mater. J. Phys.

and large-scale MD simulations with theSoc. Jpn. (2020) 2000213.

accuracy of high precision (and high cost)[3] K. Hoshino, S. Kasamatsu, et al., inelectronic structure methods.

preparation.
[4] J. Behler: Int. J. Quant. Chem. 115 (2015)
1032-1050.

52

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Formation energy calcuation of defects in candidates of thermoelectricmaterials
KAZUHIKO KUROKI
Department of Physics, Osaka University
1-1 Machikaneyama, Toyonaka, Osaka, 560-0043, Japan

Burke-Ernzerhof (PBE) exchange-correlation functional and the projector augmented wave method.
We ﬁnd that the interstitial-site doping of AE =Ca, Mg or G3 = Sc, Y, and G3 = La, Y substitutional doping for the Ca site is favorable both for
CaZn2 X2 (X = As, P) in terms of energy stability. In particular, the formation energy of the La
substitutional doping is the lowest among the considered cases both for CaZn2 X2 (X = As, P) andis negative, so that La is expected to be substituted for the Ca site and provide electron carriers
spontaneously (Fig.1). We also ﬁnd that the formation energy for CaZn2 As2 is smaller than thatfor CaZn2 P2 , so that n-type doping is considered tobe relatively easier in CaZn2 As2 than in CaZn2 P2 .

CAAL2 SI2 -TYPE ZINTL PHASECOMPOUNDS CAZN2 X2 (X = AS,P)

Thermoelectric performance of many kinds of
CaAl2 Si2 -type Zintl phase compounds has beentheoretically estimated by ﬁrst-principles bandcalculations[1]. This work suggests that CaZn2 X2(X = Sb, As, P) exhibits large ZT values in theelectron-doped regime due to its high valley degeneracy of the electronic band structure. Indeed,this kind of 122 Zintl phase compounds have beensynthesized and their electronic structure have also
been investigated theoretically. However, these 122
Zintl phase compounds become p-type semiconductors because the cationsite vacancies are easily
formed under usual synthesis conditions.
Given the above background, in 2021 ﬁscal year,
we have explored the possibility of n-type (electron
carrier) doping of CaZn2 X2 (X = As, P) using ﬁrst
principles calcuation[2]. We have considered thepossibility of n-type doping into CaZn2 X2 for thefollowing two situations: interstitial-site doping of
the alkaline earth metals AE (= Mg, Ca, Sr, Ba)and group 3 elements G3 (= Sc, Y, La), and G3
substitutional doping for the Ca site. The former
is analogous to Mg3 (Sb,Bi,Te)2 , in which excess Mg
is introduced to the interstitial site of Mg3Sb2 so as
to realize n-type doping[3]. On the other hand, thelatter is another possible way to introduce electron
carriers utilizing the diﬀerence in the valence number of the ion (Ca 2+ and G3 3+) , and the relevant
compounds with full rare earth substitution including La, CaAl2 Si2 -type arsenides RE2/3 Zn2 As2 (RE= La, Ce, Pr, Nd, Sm), do exist in reality.

LAOXS2 (X = SB, BI)

In LnOX Ch2 , electron carriers are usually introduced by ﬂuorine substitutional doping for oxygen, or the valence ﬂuctuation of cerium ions for
Ln = Ce compounds. Since atoms in insulating
layers work as a charge reservoir in both ways, thecharge carrier can be introduced without a largechange in the electronic state of the conducting layers, which is also advantageous for studying theirtransport properties. In fact, our previous theoretical study showed that the thermoelectric propertiesof LnOX Ch2 can become very good, provided thatelectron doping by introducing defects makes thesystem metallic without largely changing the elec-

In order to describe dopants and defects in a
crystal under periodic boundary conditions, the supercell approach was adopted in our calculations.
The atomic positions in the supercell were relaxedin the structural optimization after the lattice constants and atomic positions of the primitive cellare optimized. To consider the possibility of chemical doping, its defect formation energy of charged
states was calculated using the supercells. To numerically evaluate the defect formation energy ofchemical doping in CaZn2 X2 (X = As, P), we have
performed ﬁrst-principles calculations within thegeneralized gradient approximation using VASP.
We perform DFT calculations with the Perdew-

FIG. 1. The formation energy of the G3 substitution for
the Ca site against the Fermi energy ∆EF = EF − EV
for CaZn2 As2 . The group 3 elements G3 (= Sc, Y, La)
are chosen as dopant D. The dotted line denotes the
position of the conduction band minimum[2].

53Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo2
placement of O by S) and OS are not stable while
VO (O vacancy) and VS can take place, where theformation energy of VO is higher in X = Sb thanin X = Bi (Fig.2). A characteristic feature is that
VSeq (Seq is the in-plane S) becomes much morestable in X = Sb than in X = Bi, due to the formation of an Sb2 dimer and the occupation of the
impurity energy levels. The formation energies ofcation defects, XS , and SX are found to be positive. We have also found that Fluorine likely replaces oxygen for both X = Sb and Bi. The defect formation energy of FO is negative for bothcompounds, while that for X = Sb is much higher
than X = Bi. Our study has clariﬁed the stability of several point defects and suggested that thein-plane structural instability is enhanced in X =Sb. We expect that, although further study is necessary, the knowledge obtained in the study shouldbe helpful for understanding/controlling the transport properties of LaOXS2 and related compounds
by impurity doping.

FIG. 2. Defect formation energies of anion point defectsfor (a) LaOSbS2 (X = Sb) and (b) LaOBiS2 (X = Bi).
The horizontal line (the Fermi level) is restricted to theenergy range between the band edges. The values of thecharge q of the defects are shown beside the line[6].

tronic structure[4]. However, it has recently beenfound that it is experimentally diﬃcult to enhancethe electrical conductivity in X = Sb compounds
against their robust insulating nature[5].
In 2021 ﬁscal year, in order to shed light on
the puzzling transport properties of the doped layered pnictogen dichalcogenides, we have systematically investigated the defect formation energy ofseveral point defects in LaOXS2 (X = Sb, Bi) using ﬁrst-principles calculation[6]. We used the PBE
parametrization of the generalized gradient approximation and the projector augmented wave method
as implemented in VASP.

[1] H. Usui and K. Kuroki, J. Phys. Soc. Jpn. 89 (2020)
124707.
[2] K. Nishiguchi, M. Ochi, C.H.Lee, and K. Kuroki, toappear in J. Phys. Soc. Jpn.
[3] H. Tamaki, H.K. Sato, and T. Kanno, Advanced
Materials 28 (2016) 10182.
[4] M. Ochi, H. Usui, and K. Kuroki, Phys. Rev. Appl.
8 (064020) (2017).
[5] Y. Goto et al., J. Phys. Soj. Jpn. 289 (2019) 38.
[6] M. Ochi and K. Kuroki, Phys. Rev. B 105 (2022)
094110.

We have found that anion replacements SO (re-54

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Studies on superconductivity in nickelates from a multiorbital viewpoint
KAZUHIKO KUROKI
Department of Physics, Osaka University
1-1 Machikaneyama, Toyonaka, Osaka, 560-0043, Japanplan to publish the results as soon as possible.

DESIGNING NICKELATESUPERCONDUCTORS WITH d8+δELECTRON CONFIGUATIONEFFECT OF RESIDUAL HYDROGEN IN
THE INFINITE LAYER NICKELATESUPERCONDUCTORSIn 2020 ﬁscal year (the last ﬁscal year), we theoretically designed unconventional nickelate superconductors with d8 electron conﬁguration[1]. We
considered d8 nickelates with K2 NiF4 structure,
where the apical oxygens are replaced by halogens or hydrogens. One of the best candidates waselectron-doped Ca2 NiO2 Cl2 . The key point wasthat the large crystal ﬁeld splitting of the orbitalenergy levels, induced by replacing the apical oxygens by halogens or hydrogens, is approximately
equivalent to large interlayer hopping in a bilayer
system. Applying the ﬂuctuation exchange approximation to the ﬁve orbital model of these materials, we showed that the maximum superconducting transition temperature of some of them may beeven higher than that of the high-Tc cuprates.
A major drawback of our proposal was that
Ca2 NiO2 Cl2 has never been synthesized previously,
and it may indeed be diﬃcult to synthesize. On
the other hand, a similar material Sr2 NiO2 Cl2 hasbeen synthesized experimentally[2], but its lattice
constant is too large compared to Ca2 NiO2 Cl2 , resulting in a smaller crystal ﬁeld splitting, and hence
the predicted superconducting Tc , if any, is low. In
fact, a high spin state is suggested in Sr2 NiO2 Cl2 ,which is unfavorable for our mechanism to work,
Given this background, in 2021 ﬁscal year, we studied the hydrostatic pressure eﬀect on Sr2 NiO2 Cl2 in
an attempt to mimic the situation of Ca2 NiO2 Cl2 .
The structural parameters under pressure has been
determined using VASP. We have obtained ﬁve orbital models through Wannierization, and studied
the model using the ﬂuctuation exchange approximation. This is still an ongoing research, and we

Our theory on desigining d8+δ nickelate superconductors may also be related to another newly
discovered superconductor (Nd,Sr)NiO2 [3]. In thissuperconductor, the nickel electron conﬁguration
is expected to be close to d9 , but this might beaﬀected, e.g., by the reduction process during thesynthesis of the material. Since the band structureof (Nd,Sr)NiO2 itself, due to the absence of apical oxygens, resembles that of the proposed mixed
anion-nickelates, (unintentional) reduction of theelectron number in (Nd,Sr)NiO2 may also result ina similar s±-wave superconductivity.
In 2021 ﬁscal year, we have studied the eﬀect ofresidual hydrogen using the adopting the supercell
approach, which contains multiple Ni atoms in aunit cell due to the presence of eccess hyrdogens
added at the apical site periodically. The structural parameters were determined using VASP. We
have so far constructed 10 and 20 orbital modelsof these systems through Wannierization. We plan
to analyze the possibility of superconductivity inthese multiorbital models.

[1] N. Kitamine, M. Ochi and K. Kuroki, Phys. Rev.
Res. 2, 042032(R).
[2] Y. Tsujimoto, C.I. Sathish, Y. Matsushita, K. Yamaura, and T. Uchikoshi, Chem. Commun. 50, 5915
(2014).
[3] D. Li et al., Nature (London) 572, 624 (2019).

55

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First principles studies of impurity-dislocation complexes in
GaNKenji SHIRAISHI
Institute of Materials and Systems for Sustainability,
Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8601
Line defects and points defects haveresults indicate that Mg-H complexes do notintensively studied so far [1, 2]. However,form near the screw dislocation, although screw

complexes composed of line defects and point

dislocations attract both Mg and H atoms.

defects have not been considered yet. In thisFigure 2 shows the atomistic structures of the

project, we have studied dislocation-impuritycomplexes in GaN by the first principles
calculations.

b=[0001]

Fig.2: Atomistic structures of a complex
composed of a GaN screw dislocationand Mg and H impurities.
complex composed of a GaN screw dislocationand Mg and H impurities. As shown in thisfigure, the local Mg-N-H structure which is

Fig.1: Atomistic Structures of acalculation model of the complex
composed of GaN screw dislocationand Mg and H impurities.

observed in bulk GaN is broken by the effect of
the screw dislocation. This indicate thatproperties of point defects are totally changed

Calculations were performed by VASP
package.

by the effect of line defects in the dislocation-

We prepared 800 atom unit cells

impurity complexes.

which contain a screw dislocation and Mg andReferences

H impurities as shown in Fig.1.

[1] S. Berber et al. : Phys. Rev. B 77, (2008)

It is well-known in bulk GaN that Mg and H

165405.

impurities tend to form complexes of the Mg-[2]

N-H local structure. However, calculated

645.

56

K. Matsunaga: Acta Mater. , 195, (2020)

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First Pricnciples Study of N Vacancies in Si3N4 for Flash
Memory Application
Kenji SHIRAISHI
Institute of Materials and Systems for Sustainability,
Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8601
Defects in Si3N4 are now used for theapplication of 3D NAND flash memories.
Among Defects in Si3N4, N vacancies are
considered as the physical origins of charge
traps used in memory functions [1, 2].
Recently, the charge transfer between charge
traps in Si3N4 has become serious problem inthe miniaturized 3D NAND flash memories.

Fig. 2: Formation energies of various
charged state of the N vacancy as afunction of system Fermi energy. -2, 0
and +2 charged states are stable.

However, the physical mechanisms of thischarge transfer is not clarified. In this study, weinvestigated the charge transfer mechanism byenergy.

the first principles calculations.

energies of a VN-H complex as shown in Fig. 2.

First, we calculated the formation energies

As shown in Fig. 2, stable charged states are -2,of a N vacancy (VN). The results are shown in
Fig. 1.

Next, we calculated the formation

0, and 2 in the Si3N4 Fermi energy. ThisAs shown in Fig. 1, stable charged

indicate that H diffusion from one VN tostates are -1, 0 and +1 in the Si3N4 Fermi

another VN causes charge transfer when VN
traps electrons.

References[1] K. Yamaguchi, A. Otake, K. Kobayashi, andK. Shiraishi, Tech. Digest of IEEE Int. Electron
DevicesMeeting (IEDM), Baltimore, MD, USA,2009, pp. 1–4.
Fig.1: Formation energies of various
charged state of the N vacancy as afunction of system Fermi energy. -1,
0 and +1 charged states are stable.

[2] K. Yamaguchi, A. Otake, K. Kamiya, Y.
Shigeta, and K. Shiraishi,Jpn. J. Appl. Phys.,50,(2011) 04DD05.

57

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Development of quantum algorithm for imaginary-time
evolution
Yu-ichiro MATSUSHITATokyo Tech Academy for Convergence of Materials and Informatics (TAC-MI),Tokyo Institute of Technology, Ookayama, Meguro-ku, Tokyo 152-8550
hybrid classical-quantum algorithm.

Recently, the development of materials

We have developed a so-called “nonlocal”

computation methods using quantum computers
has attracted a great deal of attention. In

approximationasamethodparticular, we have paid attention to theapproximate

unitary

operators,

imaginary-time evolution method, which hasdeveloped the method to realize quantum

been known for a long time in materials

imaginary-time evolution with a quantum

computation,

thecircuit shallower than the conventional methods

development of algorithms for its execution onof imaginary-time evolution. We constructed ana quantum computer. As a result, we haveactual quantum circuit and applied it to thesucceeded in developing an algorithm for

MAXCUT problem. For the calculations, wequantum imaginary-time evolution methods [1].

verified the effectiveness of the algorithm withOn a quantum computer, only unitary

a quantum circuit simulator. The results showedtransformations and observations must be usedthat the quantum imaginary-time algorithm

to achieve the desired computation. The

does work, and at least 94% of accuracy isimaginary-time evolution operator is a non-achieved for unweighted 3-regular MAXCUT

unitary transformation, and its realization on a

problems with ten vertices. At the same time,quantum computer is non-trivial. One approach

we have found that there is an accuracy limit

developed in this study is to realize thedue to the introduction of an approximate

approximatedunitary operator, and the algorithm does not

and

have

workedimaginary-timeonevolutionforandfinding
have

completely match the exact solution.

operator on a quantum computer by finding aunitary operator that best approximates thequantum state after the operation of theReferences

imaginary-time evolution operator. This method[1] H. Nishi: npj Quantum Information 7

uses a classical computer in finding the unitary(2021) 85.

operator, and its unique feature is that it is a58

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-Principles Molecular-Dynamics Study of Structural andElectronic Properties of Covalent Liquids and Glasses under
Pressure
Fuyuki SHIMOJO
Department of Physics, Kumamoto University, Kumamoto 860-8555
It is important for the investigation of theby using a so-called two-phase simulation

dynamic properties of covalent liquids andmethod based on FPMD [1]. Both results are inglasses under pressure to consider the stability

reasonable agreement with each other.

of materials states under various conditions

We

have

also

obtainedthemeltingfrom the viewpoint of free energy. For thistemperature of binary alkali NaxRb1-x mixture

purpose, the thermodynamic integration (TI)

[2] as a function of concentration x, and foundmethod based on the first-principles molecularthat the melting temperature minimum occurs

dynamics (FPMD) simulations is useful andat x = 0.2, being in agreement with theeffective. However, the system size andexperimental observation. These two findings

simulation time are limited due to the highfor diamond and alkali mixture demonstrate thecalculation cost of FPMD. To solve this cost

validity and usefulness of the current method.

problem, we have used machine-learninginteratomic potential based on artificial neuralnetworks (ANN) trained from the results of
FPMD simulations.
In this study, we have calculated the melting
temperature of diamond, one of the typicalcovalently bonded materials, using the TI
method with the ANN potential as a function ofpressure.

Figure

1

shows

the

Fig. 1: Pressure dependence of melting
temperature of diamond.

pressuredependence of calculated melting temperature.
The results of this study are indicated by the

References

red diamond symbols. The triangle up and[1] A. A. Correa et al., PNAS 103, 1204 (2006).

down symbols indicate the results of the[2] A. Irie et al., J. Phys. Soc. Jpn. 91, 045002

previous theoretical study, which carried out(2022).

the direct simulation of solid/liquid coexistence59Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Development of first-principles electronic-structure
and transport calculation code RSPACE andsimulations for device
Tomoya ONO
Graduate School of Engineering, Kobe University
Rokkodai-cho, Nada, Kobe 657-8501

at the interface.

Silicon carbide (SiC) attracts attention because it is one of the most promising widebandgap semiconductors for developing nextgeneration switching devices operating in highpower and high-frequency applications. However, the potential of 4H-SiC has not been
fully utilized owing to the low carrier mobilityof 4H-SiC-based metal-oxide-semiconductor
field-eﬀect transistors, which is caused by defects at the 4H-SiC/SiO2 interface. To reducethe density of interface defects and increase
the channel mobilities, post-oxidation annealing using N2 O, NO, or N2 gas has been proposed. In this study, we propose the interface
atomic structures in which N atoms are incorporated by replacing C atoms, and we investigate the total energies and electronic structures of the nitride layers by density functional
theory (DFT) calculation. Universal atomicscale models describing 4H-SiC with a high-Natom-density layer in arbitrary crystal planes
are proposed. To explore the most preferable
crystal plane of 4H-SiC for N-atom incorporation, we study the stability and electronic
states of the nitride layers on a 4H-SiC surface.
We find that the nitride layers grow along thea face, which is consistent with the conclusion
obtained using the bulk model[1]. Any defectstates in the bandgap of 4H-SiC are not generated in our structure. The formation energyof the nitride layers at the topmost layer ofthe interface is smaller than that at the second
layer, indicating that the N atoms accumulate

In our model, the surface of the 4H-SiC substrate is terminated with OH groups to imitate
the SiO2 layer and the other side of the substrate is terminated with H atoms. The rectangular supercells of 10.1 Å × 10.7 Å × 26.3
Å for the a-face model, 10.1 Å × 12.3 Å × 27.3
Å for the m-face model, and 12.3 Å × 10.7
Å × 26.8 Å for the Si(C)-face model were employed. The z-axis is taken to be the direction
perpendicular to the surface. The numbers of
atoms for the a-face, m-face, and Si(C)-facemodels are 176, 208, and 208. We considered
the modification incorporating four Si vacancies (VSi s) and 16 N atoms at C sites (NC s).
The cases in which VSi s are arranged parallel to the surface were investigated to evaluate
the anisotropy of the formation energy of thenitride layer. In our model, the areal N-atomdensity is on the order of 1014 cm−1 , which isconsistent with the experimental results. Although there are two inequivalent lattice sitesof 4H-SiC, i.e., h (hexagonal) and k (quasicubic) sites, our previous studies have shownthat nitride layers tend to grow at the k site.
Thus, we consider the cases in which VSi s arearranged at the k site. The interfaces with
the nitride layer at the second layer are also
investigated for all the faces. For the firstprinciples calculation, we used the RSPACE
code[2, 3], which is based on the real-space
finite-diﬀerence approach within the frameworks of DFT and is developed in this project.

60

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

We employed the local density approximationof the DFT to describe the exchange and correlation eﬀects. Electron–ion interactions are approximated by the projector-augmented wave
method. We adopt 6 × 6 × 1 Monkhorst-packk-point meshes including a Γ-point in the Brillouin zone. The real-space grid spacing is chosen to be ∼ 0.21 Å. The structural optimization is performed until the residual forces are
smaller than 0.001 Hartree/Bohr radius.
We assumed the interaction of the SiC substrate and the arriving NO molecule at the interface as

Table 1: Formation energy of nitride layers defined by Eq. (1).
face
a
a
m
m
Si
Si
C
C

SiC(w/o mod) + 24NO

layer
top
2nd
top
2nd
top
2nd
top
2nd

Eform (eV)
−7.43
−7.10
−7.13
−7.07
−7.07
−7.01
−7.14
−7.01

annealing generates the nitride layers immedi→ SiC+ 4SiO2 + 16CO + 4N2 ,ately below the SiO2 layers without any transition layers. The diﬀerence in the formation
with SiC(w/o mod) (SiC(w/ mod) ) being the interenergy between the nitride layers at the topface without (with) the modification incorpomost and second layers is the largest for the
rating four VSi s and 16 NC s. The formationa face. This result agrees well with the obserenergy of the above interaction for generating
vation that the formation energy of the nitride
one VSi is obtained aslayer along the a face is the smallest. The local(w/ mod)densities of states of the interfaces reveal that
Eform = Etotal/4 + E(SiO2 )(w/o mod)+4µCO + 2µN − Etotal/4 − 6µNO , (1) there are no gap state at the interface. These
results imply that NO annealing after the oxwhere E(SiO2 ) is the total energy of a SiO2 idation can contribute to the reduction in the(w/o mod)unit in a bulk of quartz SiO2 . Etotal
and density of interface defects by forming the ni(w/ mod)Etotalrepresent the total energies of the tride layer immediately below the SiO2 layerinterface without and with N-atom incorpo- of the interface.
ration, respectively. In addition, µNO , µCO ,and µN are the chemical potentials of a NO
References
molecule, a CO molecule, and N atoms in a N2
molecule, respectively. We set the temperature[1] M. Uemoto, N. Komatsu, Y. Egami, andat 1000 K and the partial pressure of NO gas
T. Ono, J. Phys. Soc. Jpn. 90, 124713(pNO ) at 1 atm. The partial pressures of CO(2021).
−5(pCO ) and N2 (pN2 ) gases are set at 10 atm[2] T. Ono and K. Hirose, Phys. Rev. Lett.
and 0.25 × 10−5 atm.
82, 5016 (1999).
The formation energies of the nitride layers
are shown in Table 1. We find that all the re[3] K. Hirose, T. Ono, Y. Fujimoto, andactions are exothermic and the nitride layers
S. Tsukamoto, First Principles Calculagrowing along the a face are the most stabletions in Real-Space Formalism, Electronic
among those along the a, m, C, and Si faces.
Configurations and Transport Properties
The formation energy of the nitride layer at theof Nanostructures (Imperial College, Lontopmost layer is smaller than that at the secdon, 2005).
ond layer, resulting in the localization of the
N atoms at the interface. Note that the NO(w/ mod)61

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Analysis of Surface Reactions in Atomic Layer Etching
Processes
Satoshi HAMAGUCHI
Center for Atomic and Molecular Technologies,
Osaka University, 2-1 Yamadaoka, Suita, Osaka 567-0871
Thermal atomic layer etching (ALE)First, we examined the interaction of

can be used for precise and damageless−diketones with a metallic Ni surface,etching. For example, it can be used tousing STATE [6,7]. STATE is a simulation

fabricate magnetic tunnel junction (MTJ)

code based on density functional theory

stacks of magnetoresistive random-access(DFT) with pseudo-potentials and a plane

memory (MRAM) devices [1]. The aim ofwave basis set. It has been found that the

this

the

cleavage of a C-F or C-H bond of amechanisms of thermal ALE of metal in−diketone is more likely to take place

general, whose cycle consists of anthan the deprotonation on a metal surface.

oxidation step and an exposure step to

The surface roughness can also help to−diketones,

form

study

is

to

understand

such

as;

hexafluoroacetylacetone

C-Nibonds,

which

were

also

(hfacH),observed in experiments. The reaction

trifluoroacetylacetone(tfacH),andmechanisms that we observed in theseacetylacetone

molecules.

Forsimulations(acacH)have

revealed

the

self-

instance, thermal ALE of Ni with hfacH

limiting nature of thermal ALE for metals

starts with the deprotonation of hfacH

with −diketone gases in general, which ismolecules after the oxidation of the Ni

consistentsurface.

observations [5,8].

By

increasing

thesurfacewithearlierexperimental

temperature in the range of 300 ~ 400 °C,

Second, we examined the interaction ofvolatile nickel complex Ni(hfac)2 and−diketones with oxidized metal surfaces,

water H2O molecules are formed on theusing TURBOMOLE program [9]. It has

surface and remove the NiO layer. In thisbeen found that, as −diketones approachprocess, once a new metallic Ni surface

anappears, no further etching occurs, which

deprotonate. In this way, volatile metal

ensures the self-liming nature of thecomplexes and H2O molecules can bethermal ALE process [2-5].

formed62oxidized

when

metal

surface,−diketone

theymolecules

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

interact with an oxidized metal surface.

measurement. The results indicate that

The reaction energies of such interactions

carbon contained in the incident CF3+ ionshave been evaluated from the simulations

remove oxygen of Y2O3 more efficiently.

and we found that etch oxidized metal has

Furthermore, high-energy incident CF3+

its own desorption energy or temperatureions promoted the formation of Y-F bonds

to volatilize its metal complex. Therefore,

near the surface region of Y2O3 films more

the applied temperature can be used as aefficiently than low-energy incident CF3+

self-limitingtargetions. Our simulation results are found tocontains different materials. In addition,

be consistent with these experimental

hfacH is found to get more stably bonded

observations.

conditionifthethan acacH or tfacH with an oxidized
metal surface because the negatively

REFERENCES

charged F atoms in CF3 are repelled by O

[1] G. S. Oehrlein and S. Hamaguchi,

atoms of the surface, which forces a

Plasma Sources Sci. Technol. 27, 023001

deprotonated hfac to adsorb on the(2018).

surface nearly vertically. This does not[2] A. H. Basher et al., JJAP 59, 090905

happen to acacH, whose CH3 groups tend

(2020).

to interact directly with the metal oxide.

[3] A. H. Basher et al., JVST A 38, 052602

As a result, hfacH is considered more(2020).

suited for ALD processes than other β-[4] A. H. Basher et al., JVST A 38, 022610

diketones.

(2020).

In addition, we also examined the[5] T. Ito, et al., AVS 65th (2018).

interaction of Y2O3 with incident F+ and[6] I. Hamada, Physical Review B 89 (12),CF3+ ions, using ion beam experiments

121103, (2014).

with XPS and DFT simulations. It has[7] Y. Morikawa et al., Phys. Rev. B, 69,

been found that the penetration depth of

041403 (2004).

F atoms in Y2O3 exposed to F+ or CF3+ ion[8] H. L. Nigg and R. I. Masel, JVST A 17,

injection depends essentially on the3477, (1999).

incident energy and hardly depends on[9]

the ion dose once a sufficient amount ofhttps://www.turbomole.org.

fluorine is provided to the surface andsaturate. We also obtained the atomic
concentration depth profiles by XPS

63

TURBOMOLE

V7.5

2020,

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles study on atomic and electronic
structures of graphene/SiC interfaces
Yoshiyuki EGAMI
Division of Applied Physics, Faculty of Engineering, Hokkaido University
Kita 13, Nishi 8, Kita-ku, Sapporo, Hokkaido 060-8628

Graphene is in the limelight as a post-siliconmaterial in next-generation devices due to its
unique electronic and mechanical properties
and conductivity. In recent years, the epitaxial growth in the silicon carbide (SiC) pyrolysis method has attracted attention as one ofuseful methods for producing large-scale andhigh-quality graphene. However, in the pyrolysis method, the bottom graphene layer iscovalently bonded to the Si atom of the SiC
substrate. To separate the bottom layer fromthe substrate, atomic intercalation technique
is employed. As the result, the bottom layer
becomes an ideal single-layer graphene, andthe Dirac cone appears in the electronic band
structure. On the other hand, atomic vacancies are formed in the graphene during the epitaxial growth process. Since the intercalants
can be bonded to the dangling bonds, which
are aﬀected on the energy band structure andthe position of the Dirac point of graphene.
In this study, ﬂuorine atoms are adopted asintercalants, and the eﬀect of the bond between a ﬂuorine atom and a dangling bondon the electronic state of graphene is investigated using the ﬁrst-principles calculation
code: RSPACE[1, 2]. We investigate the difference in the eﬀects of the number of ﬂuorine
atoms adsorbed on the single vacancy in thegraphene on the electronic state. It is foundthat, in the model where two ﬂuorine atoms
are bonded to the defect, the dangling-bondstate is hybridized with the conduction band
bottom of SiC substrate. Moreover, the conju-

gated π state of graphene is extended over theenergy range in which the bandgap of SiC substrate is located (Fig. 1). On the other hand,in the system produced by hydrogen intercalation, the dangling bond state is in the middleof the band gap of the SiC substrate and extends to the inside of the substrate, which cancause a large leakage current.
This work has been performed on System Bof the Supercomputer Center, the Institute for
Solid State Physics, the University of Tokyo.

Figure 1: Local density of states of graphenewith single vacancy on ﬂuorine-intercalated
SiC substrate. Two of dangling bonds are terminated by ﬂuorine atoms.

References[1] K. Hirose et al.: First-principles calculations in real-space formalism (Imperial
College Press, 2005).
[2] Y. Egami, S. Tsukamoto and T. Ono:
Phys. Rev. Res. 3, 013038 (2021).

64Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

High-Throughput Polymer Electret Material Screening Based on
Solid-State Quantum Chemical Calculation and Deep Learning
Yucheng ZHANG, Zetian MAO, and Yuji SUZUKI
Department of Mechanical Engineering
The University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656
Electret is a dielectric material withionization potentials (IPs) by the density

quasi-permanent charges and can trap

functional theory (DFT) with polarizable

charges stably for decades. Vibrational

continuum model (PCM) correction agree

electret energy harvester (VEEH) canwell with the trends of the surface charge

generate electricity efficiently from kinetic

density and charge stability measured byenergy via the electrostatic induction. The

the wet experiments.

VEEH with bipolar (both positively andMost recently, by employing the de novo

negatively charged) electret leads to fourmolecule generation and the functionaltimes output power theoretically, thusgroup enrichment analysis, we successfully

specifying the importance of designing

proposed the electret CTX-A/APDEA withhigh-performance bipolar polymer electrets.

high bipolar charging performance via

CYTOP (Cyclic Transparent Optical

combining CTX-A with N-(3-aminopropyl)

Polymer, AGC Chemicals) is one of thediethanolamine [2]. This result encouraged

best polymer electret materials, with three

us to move further to search new amines incommercialized types displayed in Fig. 1.

the open-source database PubChem, aimed

Recently, we utilized quantum chemical

at designing better bipolar polymer electret

analysis for evaluating the charge trap of

based on CYTOP CTX-A after the

CYTOP electrets [1]. It is found thatdehydration reaction.

although the repeat-unit-number of CTX-A
is above 1000, the trapped charge islocalized at the amide bond, which isformed by dehydration reaction between
the carboxyl end group and the amine
introduced. It is also shown that thecomputed electron affinities (EAs) and

Fig. 1. Chemical structure of CYTOP series.

65Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Fig. 2 illustrated the prediction results on
the established dataset with one kind ofdeep learning models named MEGNet [3].
The mean absolute error (MAE) of IP is as
small as 0.326 eV, with computational costdecreased significantly. Therefore, thepresent approach will be employed forscreening bigger database for new highperformance bipolar polymer electrets.
The geometry of the molecule is relaxed
Fig. 2. Prediction of PCM-DFT IPs forunder the open-shell condition and then the

training and test sets using MEGNet.

energy of both open-shell and closed-shell
system are calculated and compared. All

Although solid-state quantum chemical

the quantum chemical calculations arecalculation based on PCM-DFT can give

made with the software GAUSSIAN [4] at

quantitative EA and IP values, the average

the DFT level with the CAM-B3LYP

computation time is 41.4 CPU hours perfunctional and 6-31+G(d,p) basis set, whilemonomer molecule and even as high as

PCM is combined for solid-state analysis

2500 CPU hours per hexamer molecule. To

due to the solvation effect correction.

reduce such heavy computational cost, we
employdeeplearning

modelsThis work was partially supported byto

JST CREST Grant Number JPMJCR15Q3

approximate the EAs and IPs efficiently.

and JPMJCR19Q1. The calculations are

To prepare the labeled data, 5522 amines

performed with SGI ICE XA ISSP systemwere initially sampled from the PubChem.

B and every calculation is done with 1 node

Thereafter, the virtual dehydration reaction(24 CPUs).

was conducted and PCM-DFT computation

References

was automatically made. Subsequently, a[1] Kim S., et al., J. Phys. Chem. B, Vol.

dataset mapping the PCM-DFT IPs and124, No. 46, pp. 10507-10513 (2020).

molecular structures was establised and[2] Zhang, Y., et al., Appl. Phy. Lett., 118,

splited into training-validation-test sets

223904 (2021).

with the ration of 80%, 10%, and 10%,[3] Mao, Z., et al., IEEE Conf. Electr. Insul.

where the training set was used to teach theDielectr. Phenom., pp. 667-670 (2021).

deep learning models and the performance[4] M. Frisch et al., “Gaussian 09, Rev. A.

of the models were tested on the test set.

1,” Gaussian, Inc. Wellingford, CT, (2009).

66

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles calculations on muon in solids
Mineo SAITO
Institute of Science and Technology, Kanazawa University
Kakuma, Kanazawa, 920-1192 Japan

Hydrogen plays an importantdensity functional theory is 1.36 eV, whichrole in semiconductor materials;

is close to the experimental value. On the

Well known examples are hydrogen

other hand, the value calculated based on

passivation and carrier injection. The

the GGA is small (0.71 eV). The Simuon spin resonance is an effective tool tomonovacancy is found to form a triangle

determine the dynamics of hydrogen inof the C2v symmetry and to induce a mid-

semiconductors since muonium consisting

gap level having large amplitudes at theof a positive muon and an electron bears

equivalent two atomic sites.

similarities

tohydrogen.

hyperfinestructures

Observed

provideWe findusefulthat the monohydrogen-vacancy pair induces the gap level.

information on the atomic site ofOnmuonium, and thus that of hydrogen. In a

dihydrogen passivates the defect level. It

previous study, we clarified that negative

is suggested from out calculation thatFermi contact terms are caused by thetrihydrogen-vacancy pair leads to n-type

electron correlation and observing thedoping. We conclude that hydrogenation

negative value suggests the formation ofplays various roles in barium silicide,resonance bond of hydrogen [1].

which depends on the concentration of

the

other

hand,

we

find

that

One of the scientifically importanthydrogen. The muon spin resonancenew targets is hydrogen impurity inmeasurement is expected to confirm thebarium silicide (BaSi2), which is a

atomic geometry of the monohydrogen-candidate for solar cell materials. It isvacancy pair.

known that hydrogen passivates the defectlevel of Si vacancy. This crystal includes[1] M. N. Manal, S. Minami, F. Ishii and

Si tetrahedrons and electrons transfer from

M. Saito, Jpn. J. Appl. Phys., 58 (2019)barium to silicon. We first confirm that

081008.

this crystal has an indirect gap and the gap
value calculated based on the hybrid

67

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

GW + BSE calculation of RIXS spectraby the all-electron mixed basis program, TOMBO
Kaoru OHNO
Graduate School of Engineering, Yokohama National University79-5 Tokiwadai, Hodogaya-ku, Yokohama 240-8501
The initial state of X-ray emission spectroscopy (XES) and resonant inelastic X-ray emission spectroscopy (RIXS) is a highly excited
eigenstate with a deep core hole after the X-ray
photoelectron spectroscopy (XPS) and X-ray
photoabsorption spectroscopy (XAS) processes,
so that the XES and RIXS calculation offers agood example of extended quasiparticle theory(EQPT) [1] associated with many-body pertur-

Fig. 1: Calculated (above) and experimental [3]bation theory. Moreover, to overcome the basis(below) RIXS spectra of a CH3OH molecule.

set incompleteness problem, we adopt the all-theory [4], we gave a justification and comment

electron mixed basis program, TOMBO, whichof applying the method relying on time-depen-

uses both plane waves and numerical atomic

dent density functional theory as well as theorbitals as basis functions. We have success-one-shot GW + BSE approach to this problem.

fully applied this approach in our previous XES[2] study. Here we applied the standard one-

References

shot GW + Bethe--Salpeter equation (BSE)[1] K. Ohno, S. Ono, and H. Isobe: J. Chem.

approach in MBPT to this problem on the basis

Phys. 146 (2017) 084108.

of EQPT and analyze XES and RIXS energies[2] T. Aoki and K. Ohno: Phys. Rev. B 100

for CH4, NH3, H2O, and CH3OH molecules.

(2019) 075149.

Figure 1 shows the RIXS spectra of a[3] A. Benkert, F. Meyer, D. Hauschild, M.

CH3OH molecule calculated by our method

Blum, W. Yang, R. G. Wilks, M. Bär, F. Reinert,

together with the experimental data [3]. Both

C. Heske and L. Weinhardt,: J. Phys. Chem. A

results coincides very well without introducing

120 (2016) 2260.

any empirical fitting parameter.

[4] T. Nakashima, H. Raebiger, and K. Ohno:
Phys. Rev. B 104 (2021) L201116.

Moreover, according to extended Kohn- Sham

68

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Theoretical study of electric double-layer structure of MXene
electrode dependent on surface functional groups
Atsuo YAMADA
Department of Chemical System Engineering,
The University of Tokyo, Bunkyo-ku, Hongo, Tokyo 113-8656
Aqueous electric double-layer capacitor

while electrolytes were described using the(EDLC) is promising candidate for large-scale

RISM method. All the ESM-RISM calculations

stationary storage due to non-flammability, lowwere conducted with Quantum Espresso code.

materials cost, and non-toxicity of aqueous

The calculated EDL capacitance increased

electrolyte. Transition metal carbides/nitrides

with the increase in the atomic number of(MXenes) are appealing electrode materials for

halogen atoms (F < Cl < Br < I) (Fig.1a). The

aqueous EDLCs owing to their high specificincreased EDL capacitance is due to the biased

capacitance (> 300 F/g) in conventional

distribution of excess electrons in the electrodeaqueous electrolyte. Herein, we investigated thetoward the electrolyte, thus reducing the EDL

dependence of EDL capacitances on the surface

thickness (Fig. 1b). The induction of electrodefunctional groups of the MXene electrodes [1].

electrons into the EDL region is facilitated by

The ESM-RISM method was adopted tothe smaller negative charges of surface halogen

evaluate the capacitance of the MXene

atoms with lower electronegativity. Therefore,

electrodes. The simulation cells are composed

lowering the valence numbers of the MXeneof MXene Ti3C2T2 with different surface

surface atoms should enhance the capacitance.

functional groups (T = F, Cl, Br, and I) and theReferences

1 M KOTf/H2O electrolyte. Here, MXene[1] T. Shimada N. Takenaka, Y. Ando, M. Otani, M.

electrodes were treated quantum-mechanicallyOkubo, A. Yamada, Chem. Mater., 34, 2069 (2022).

using the density functional theory method,

Figure 1 (a) Calculated capacitances per unit surface area of MXene. (b) Spacial distributions of excess electronsof the electrode (red) and counter positive charges of the electrolyte (blue) at the electrode/electrolyte interface.

69

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Elucidation of mechanism and theoretical prediction ofsuper-functional materials using high-density hydrogen
Shinji TSUNEYUKI
Department of Physics, School of Science,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-Ku, Tokyo 113-0033
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Hydrogen, a ubiquitous element, is present as an

exhibit without experiments. In this project, we

impurity in various materials and is known toapplied the methods to elucidate the mechanism

significantly change materials' electrical andof the photo-induced metallization of Yttrium

mechanical properties by assuming different

oxyhydride recently discovered by a Tokyo Techchargegroup experimentally [1].

statesdependingontheatoms'

environment. In recent years, record-breaking

Insuperconductors and superionic conductors havesignificantly reduced the electrical resistance ofbeen discovered in which protons and hydrides

an insulator YOxHy epitaxial thin film, and theare major components. One of the most common

metallic state was maintained on a scale ofproblems in studying these materials is theseveraldifficulty in experimentally quantifying and

YOxHy on glass substrates exhibits a one-digit

observing the position of hydrogen, which has

decrease in electrical resistance upon solar

no inner-shell electrons.

irradiation, while the epitaxial thin film has

Electronic structure calculations based onachieved a resistance reduction of more thandensity-functional theory and first-principles

three orders of magnitude. Furthermore, UV

molecular

allow

laser irradiation of the epitaxial thin film

quantitatively predicting materials' structure andreduced the electrical resistance by more thanphysicalseven orders of magnitude.

dynamics

properties

methods

without

empirical

theirexperiment,

days.

light

Conventional

irradiation

polycrystalline

parameters. By these methods, the energy of the

To explain these results, we constructed asystem and the forces acting on the atoms can be

structural model based on high-resolution

calculated using atomic configurations as input

measurements of the local structure anddata. In principle, it is possible to predict how

chemical composition, and calculated themuch and where hydrogen atoms are stable in aelectronic structure based on this model. As acrystal and what physical properties they will

result, we found that the hydrogen atoms change

70

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

their position and the charge state, leading to thefound that the ion conductivity of the high-T

long-life metallization [1,2].

phase is enhanced one order of magnitude by the

Another topic in this project is the study ofconcerted motion of Li ions and another one

ionic superconductor LiCB9H10. We first appliedorder of magnitude by the rotation of thethe data assimilation method for structure

CB9H10- complex [4].

prediction [3] to find no significant difference
between the high-T and low-T structures. From

References

a first-principles molecular dynamics simulation,[1] Y. Komatsu et al.,, Chem. Mater. 34, 3616

we found that the motion of Li ions strongly(2022).

correlates with the rotation of the CB9H10-[2] R. Sato et al., to be submitted.

complex. Finally, by constructing neural-[3] N. Tsujimoto, D. Adachi, R. Akashi, S.

network interatomic potentials, we calculated

Todo and S. Tsuneyuki, Phys. Rev. Materials 2,

the ion conductivity by a large-scale molecular

053801 (2018).

dynamics simulation in good agreement with

[4] R. Sato et al., submitted.

experiments. Through these simulations, we

71

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Reduction of Rare Metals in Fuel Cell Catalysts
and Hydrogen Permeable Membrane
Norihito Sakaguchi
Center for Advanced Research of Energy and Materials, Faculty of Engineering,
Hokkaido University, Sapporo, Hokkaido 060-8628

We investigated the catalyst adsorption

diffusion and coarsening. In addition, we also

and diffusion properties on light-element

revealed that the Cl- ion adsorption enhances

doped graphene supported on the magnetic

the anisotropic catalyst nanorod growth.[1]substrates and hydrogen dissolution and

We

also

investigatedthehydrogendiffusion properties in the amorphous ceramic

dissolution and diffusion properties in thematerials, with the aid of the first-principles

amorphous ceramic materials. We considered

calculation based on the density functional

the Al2O3 as the hydrogen permeation

theory (DFT) and machine learning potentials.

membrane.

Since

DFT

calculationsofAt first, we investigated the stabilities ofamorphous materials require quite large

Pt atoms on light-element doped graphenecomputing costs, we constructed the machine

supported on the magnetic substrates. We

learning

considered Fe and Ni as the magnetic

interatomic potentials obtained with DFT

substrates. We performed the total energy andcalculations. We performed melt-quenching

electronic structure calculations using The

technique to obtain the amorphous structure.

Vienna Ab initio simulation package (VASP).

We performed 2000 DFT calculations with

We revealed that the adsorption of Pt atoms

various hydrogen positions in amorphouson light-element doped graphene with flat

Al2O3 to construct the database for machine

structures

becomes

supportingon

potential

toreproducethestrongerthroughlearning. The obtained machine learningmagnetic

substrates.

potentials reproduce the stable adsorption

However, in most cases of steric light-

structures, however, cannot reproduce theelement

diffusion

doped

the

graphene,themagnetic

barriers

obtained

with

DFT

substrates weaken the Pt adsorption on thecalculations

light-element doped graphene. We also found

Referencesthat the magnetic substrates reduce the[1] M. J. Saw, M. T. Nguyen, Y. Kunisada, T.

diffusion barriers of Pt atoms. Therefore, theTokunaga, T. Yonezawa, ACS Omega 7

magnetic substrates enhance the catalyst(2022) 7414.

72

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Reduction of Rare Metals in Formic Acid Decomposition
Catalysts and Oxygen Storage Materials
Yuji KUNISADA
Center for Advanced Research of Energy and Materials, Faculty of Engineering,
Hokkaido University, Sapporo, Hokkaido 060-8628

We investigated the catalyst adsorption

graphene. This is because the active dangling

and diffusion properties on light-elementbond forms strong chemical bonds with

doped graphene supported on the magnetic

magnetic surfaces, hence, Cu atoms cannot

substrates with the aid of the first-principles

form stable chemical bonds with the light-

calculation based on the density functional

element doped graphene. We also investigatedtheory (DFT). Besides, we performed electron

the diffusion properties of Cu atoms on light-

energy loss spectroscopy (EELS) simulations

element doped graphene supported on theof heteroatom-substituted brownmillerite-typemagnetic substrates. We found that theoxides to characterize the EELS spectra

magnetic substrates reduce the diffusion

obtainedbarriers of Cu atoms. In addition, we also

withtransmission

electronrevealed the effect of Cl- ions on the catalyst

microscopy.
At first, we investigated the stabilities ofnanorod growth.[1]Cu atoms on light-element doped graphene

We performed EELS spectra of the

supported on the magnetic substrates. We

partially

considered Fe and Ni as the magnetic

considered Sr dopants, which show promising

substrates. We performed the total energy andoxygenelectronic structure calculations using The

obtained EELS spectra show good agreement

Vienna Ab initio simulation package (VASP).

with the experimental results.

We installed a parallelized VASP with Intel®

References

MPI Library and Intel® Math Kernel Library.

[1] M. J. Saw, M. T. Nguyen, Y. Kunisada, T.

We revealed that the adsorption of Cu atoms

Tokunaga, T. Yonezawa, ACSOmega 7

on light-element doped graphene becomes(2022) 7414.

stronger

with

the

magnetic

substituted

absorption

Ca2AlMnO5.

properties.[2]

We

The

substrates.

[2] K. Tanahashi, Y. Omura, H. Naya, K.

However, in cases of some light-element

Miyazaki, G. Saito, Y. Kunisada, N.

dopants, the magnetic substrates weaken the

Sakaguchi, T. Nomura, ACS Sustain.

Cu adsorption on the light-element doped

Chem. Eng. 9 (2021) 9317.

73

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles calculation of energy materials
Fumiyuki ISHII, Naoya YAMAGUCHI, Rifky SYARIATI,
Hikaru SAWAHATA, Hana Pratiwi KADARISMAN, Kaiki SHIBATA,
Sefty YUNITASARI, Yaotang ZHANG, Vierta SARASWATI,
Takahiro TOMINAGA, Yedija Yusua Sibuea TEWENG, Juhri HENDRAWAN
Nanomaterials Research Institute, Kanazawa University
Kanazawa, 920-1192, Japan

1. Topological thermoelectric materials
Thermoelectric conversion based on theanomalous Nernst eﬀect has received a lot of
attention as ANE achieves more eﬃcient heatrecovery and heat recovery from more heatsources than those based on the Seebeck effect. For the anomalous Nernst eﬀect, thetopological nature of the wave functions, andthe Berry curvature, a ﬁctitious magnetic ﬁeld,
is essential. We have studied the anomalous
Nernst eﬀect of iron-based ferromagnetic materials such as Fe3 X(X=Al, Ga)[1] and FeCl2 [2]by ﬁrst-principles density functional calculations. To evaluate the anomalous Nernst coefﬁcient, a detailed Fermi energy dependence of
the anomalous Hall conductivity is required.
We have implemented a ﬁnite diﬀerence algorithm in OpenMX[3] to calculate the Berrycurvature of metallic system given by the averaged local Berry phase of the momentum
space[4]. This implementations open up the
possibility of eﬃcient thermoelectric material
design based on high-throughput screening.
2. Photovoltaics
Material property data was extracted fromabout 1500 papers on organic thin ﬁlm
non-fullerene solar cell materials, and highperformance material conditions were derived
based on data statistical methods[5]. One of
the conditions obtained was the empiricallyknown ”open-circuit voltage (Voc) of 0.8-0.9 V
is appropriate”, and the database we createdrediscovered this rule of thumb. In addition, it
was suggested that the sulfur atom content is
an important factor for high performance photovoltaics materials.
We have performed ﬁrst-principles calculations of the topological phase transition instrained inorganic halide perovskite CsPbI3 .
For the experimental lattice constant, thebandgap decreases from 0.15 eV to 0 eV
when a strain of c/a > 1 is applied. Increasing the strain will gradually open the
bandgap. We calculated the topological Z2 invariants that characterize the topological insulating phase. From the results, the bulk crystal CsPbI3 showed a strain-induced topologicalphase transition, changing from a trivial insulator to a topological insulator phase. This
topological phase transition of the strained inorganic halide perovskite has the potential tobe applied to new energy devices.
3. Graphene and diamond
Applying the supercell scheme to the study ofalloys and various perturbed system including
impurities and vacancies raises band folding
problems. The unfolding method is expected
to obtain the meaningful energy-momentum
relationships by avoiding it. The method proposed by Chi-Cheng Lee et al proved[6] very
useful for calculating the electronic structureof a vacancies system with OpenMX code. We
have calculated and observe the graphene ﬂake
size dependence of Dirac cone in the unfolding74

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

band structures.
Negatively charged diamond nitrogenvacancy centers are promising for many
applications, such as high-resolution sensors
for electric ﬁelds and quantum computer
Qubits . Density functional calculations areoften performed in supercell systems. However, as the number of supercells increases,
the ﬁrst Brillouin zone (FBZ) shrinks and theband becomes more complex. An unfolding
band method calculation was performed to
study the electronic properties of the nitrogenvacancy center diamond using the OpenMX
code. Our results show that the unfoldedband of the bulk diamond supercell and thenitrogen vacancy center of the diamond havefewer lines than those of the traditional folded
band. The shape of the unfolded energy band
matches the primitive cell band. The dedailed
electronic band structure of diamond withnitrogen vacancy center are analyzed.

that has been observed in experiments and 3
brand new structures, and then performed detailed ﬁrst-principles analysis of the elctronic
structures to predict their physical properties.
We calculated the Z2 topological invariant andspin-textures of the new structures. Amongthe 3 new structures, we found that two structures with centro symmetry were strong topological insulators, and one structure without
centrosymmetry had wide band gap and strong
Rashba splitting.

References[1] A. Sakai, S. Minami, T. Koretsune, T.
Chen, T. Higo, Y. Wang, T. Nomoto, M.
Hirayama, S. Miwa, D. Nishio-Hamane,
F. Ishii, R. Arita and S. Nakatsuji, Nature 581, 53(2020).
[2] S. Rifky, S. Minami, H. Sawahata, and F.
Ishii, APL Materials 8, 041105 (2020).

The Al2 O3 /diamond(111) interface containing three types of impurities was calculated.
We predicted three diﬀerent
Al2 O3 /diamond(111) interface structures andinvestigated their band edge changes. Bandedges are evaluated by calculating the local
electron density. Compared to diamond bulk,
the band of Al2 O3 bulk is generally shifted
downward. The band oﬀsets for dangling
bonds, C-H bonds, the ideal O-terminal interface are estimated. Our research helps toassess the impact of interface structure on device performance and provide a reference for
improving experimental methods.

[3] T. Ozaki et al., http://www.openmxsquare.org/[4] H. Sawahata, N. Yamaguchi, S. Minami,
and F. Ishii, arXiv:2204.05949.
[5] N. Yamaguchi, H. Sano, H. Sawahata, M.
Nakano, T. Taima, F. Ishii, M. Karakawa,
Jpn. J. Appl. Phys. 61 030905 (2022), selected for Spotlights 2022.
[6] C-C. Lee, Y. Yamada-Takamura, and T.
Ozaki, J. Phys.: Condens. Matter 25345501 (2013).

4. Crystal structure prediction for spintronics materials
Bismuth tellurium iodide (BiTeI) has attracted
great interest as a potential material for spintronics device. We have performed crystal
structure prediction for BiTeI crystals by genetic algorithms and density functional calculations. We found 4 BiTeI crystal structureswith high stability , including one structure75

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Development of crystal structure prediction methods using
machine learning
Tomoki YAMASHITATop Runner Incubation Center for Academia-Industry Fusion,
Nagaoka University of Technology
1603-1 Kamitomioka-machi, Nagaoka, Niigata, 940-2188
We

We search for solid electrolyte materials that

have

also

studiedonalgorithm

can be synthesized in the systems of Li-Mg-

development for CSP. Look Ahead based on

SiO4 using crystal structure prediction (CSP) by

Quadratic Approximation (LAQA), which is

CrySPY [1]. We investigate whether it is

one of the selection-type algorithms wepossible to synthesize structures with high ionicpreviously developed [3], can control the

conductivity by changing the ratio of Li to Mg

optimization priority of the candidates. We

from

CSP

have proposed an improved score of LAQA,

simulations with random search. In the CSP

where the stress term is added to overcome thesimulations, first-principles calculations were

drawbacks

carried out using the VASP code [2], and asimulations by this improved algorithm are

generalized

performed

Li2MgSiO4.

We

gradient

performed

approximation

was

of

to

the

previous

investigate

score.

theCSP

searchingemployed for exchange-correlation functional.

efficiency for typical materials such as Si (16The compositions calculated in this study were

atoms), Al2O3 (10 atoms), NaCl (16 atoms), andLi2MgSiO4,

andSrCO3 (20 atoms). We used the CrySPY code

Li2.5Mg0.75SiO4. The number of searching trials

for the CSP simulations. For the system offor each composition was 200, which were

SrCO3, we employed the VASP code to

conducted using the ISSP supercomputer. Ourevaluate the energy and to carry out the local

calculated formation energies showed thatstructure

Li2MgSiO4 can be synthesized, however,supercomputer. For the other systems, theLi1.5Mg1.25SiO4are

QUANTUM ESPRESSO was also used. The

unlikely to be synthesized, leading to phasetotal-energy calculations were performed using

separation. These results are consistent with ourthe density functional theory with the projector-

experiments. Therefore, it is difficult toaugmented wave method. The generalized

synthesize solid electrolytes in the system ofgradient approximation was used for exchange-

Li-Mg-SiO4 except for Li2MgSiO4.

correlation functional. Figure 1 illustrates theLi1.5Mg1.25SiO4,

and

Li2.5Mg0.75SiO476

optimization

using

the

ISSP

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

number of total optimization steps required toReferences

find the most stable structure, obtained by the[1] T. Yamashita, S. Kanehira, N. Sato, H. Kino,

CSP simulations using random search (RS),H. Sawahata, T. Sato, F. Utsuno, K. Tsuda, T.

Greedy, and LAQA. The searching efficiencies

Miyake, and T. Oguchi: Sci. Technol. Adv.

of LAQA with the weight parameter of stress

Mater. Meth. 1 (2021) 87.

wS of 10 (red bars in Fig. 1) are not always the[2] G. Kresse and J. Furthmüller: Phys. Rev. B

best, however, it demonstrates high searching

54 (1996) 11169.

efficiency on various data sets. In other words,[3] K. Terayama, T. Yamashita, T. Oguchi, andthe introduction of the stress term makes

K. Tsuda, npj Comut. Mater.: 4 (2018) 32.

LAQA more robust and versatile [4].

[4] T. Yamashita and H. Sekine: Sci. Technol.
Adv. Mater. Meth. 2 (2022) 84.

Fig. 1: Number of total optimization steps required to find the most stable structure.

77

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Development of Exciton Maps for Identifying Active/Inactive
Thermally Activated Delayed Fluorescence Molecules
Yoshifumi NOGUCHI
Department of Applied Chemistry and Biochemical Engineering, Graduate School ofEngineering, Shizuoka University, Johoku 3-5-1, Hamamatsu, Shizuoka 432-8561, Japan
electroluminescencethat our method successfully reproduced thequantum efficiency in the third-generation

spectra. By sandwiching an arbitrary operator

organic light-emitting diodes is based on thewith the exciton wave functions, we estimatedmechanism of thermally activated delayed

the expected values such as the exciton size,

fluorescence (TADF), in which triplet spinelectron-hole separation distance, electron (orexcitons are converted to singlet spin excitonshole) delocalization, and exciton binding

via reverse intersystem crossing (RISC). Theenergy based on two-particle picture andexisting energy barrier between a singlet spin

proposed an exciton map [2]. We revealed that

exciton and a triplet spin exciton, commonly

the TADF active/inactive molecules could be

named as singlet-triplet splitting (Δ!" ), shouldclassified well in the exciton map by choosing

be almost zero for TADF to occur. In recent

the exciton binding energy as the vertical axis

experimental

certainand the ratio of electron and hole delocalization

carbazole benzonitrile (CzBN) derivatives, inas the horizontal axis. Indeed, we concluded

which the carbazole moiety acts as an electronthat there are two kinds of TADF active

donor and the benzonitrile moiety acts as anmolecules existed based on the following: first,

electron acceptor, can be TADF active andthe ratio of electron and hole delocalization wasinactive: nevertheless, their Δ!" are typically

proportional to the exciton binding energy, andabout 0.2 eV. However, we must reconsider thesecond, the hole state spatially localized moremechanism of TADF to completely understandthan the electron.

A

high

internalstudies,

however,it. We have studied the mechanism of TADF
using the first-principles GW+Bethe-Salpeter

References

method and exciton analysis method with

[1] T. Hosokai, H. Noda, H. Nakanotani, T.

exciton wave functions [2]. Before discussing

Nawata, Y. Nakayama, H. Matsuzaki, and C.

the optical properties of CzBN, we compared

Adachi, J. Photonics Energy, 8, (2018) 032102.

the simulated photoabsorption spectra with the[2] Yoshifumi Noguchi, J. Chem. Phys., 155,

available experimental spectra and confirmed(2021) 204302.

78

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Maximization of interface magnetoelectric couplingwith Bayesian optimization
Yoshihiro GOHDA
Department of Materials Science and Engineering, Tokyo Institute of Technology
J1-3, Nagatsuta-cho 4259, Midori-ku, Yokohama 226-8502, JapanTo reduce energy consumption in magnetic
devices, the control of magnetism by the electric current through the magnetic ﬁeld should
be replaced with that by the electric ﬁeld. Interface multiferroicity is promising, because
the magnetoelectric (ME) coupling is muchlarger than that of single phase multiferroics [1]. In addition to choosing optimal materials for ferromagnetic and ferroelectric phases,
the interface modiﬁcation is also a choice to enhance the ME coupling.

density of states close to the Fermi level [2]. In
contrast, the following three cases exhibit tiny
ME coupling. Insuﬃcient occupancy causes
stronger ionicity reducing interaction with Ti.
Excess occupancy makes the hybridized states
too deep. In addition, 4d electrons are too itinerant to have suﬃcient hybridization with Ti.
In addition, the enhancement of the ME
coupling due to the change in the magnetic
anisotropy has also been studied. Changes in
the magnetic anisotropy of a whole ferromagnet by the exchange coupling [5] are further enhanced by the strain eﬀect at the interface [3].

In this project, we performed large-scaleﬁrst-principles calculations of multiferroic interfaces [2, 3, 4].
We explored better
insertion of transition-metal atomic layers
into Fe3 Si/BaTiO3 (001) interfaces using theBayesian optimization.
As is shown in
Fig. 1(a), the Bayesian optimization improved
the interface ME coupling by 55% that is obviously much better than the random search.

(a)

(b)

Figure 1:(a) The Bayesian optimization of
Fe3 Si/BaTiO3 interfaces with inserted atomic layers. (b) The change in the interface magnetic moments.

As a by-product of the Bayesian optimization, we successfully provided an explanation of the mechanism for the multiferroicity enhancement [2]. Figure 1(b) shows thedependence of the change in the interface
magnetic moment on the exchange-correlation
functional and the electron localization controlled by the Hubbard U parameter. The
large ME coupling due to the change in themagnetic moment in V comes from the nearly
half occupied 3d bands causing peaks in the[1] T. Taniyama, J. Phys.: Condens. Matter 27,
504001 (2015).
[2] R. Costa–Amaral and Y. Gohda, Phys. Rev.
Appl. 15, 064014 (2021).
[3] S. Fujii, Y. Gohda, et al., NPG Asia Mater., inpress.
[4] S. Tsuna, R. Costa–Amaral, and Y. Gohda,
submitted.
[5] K. Fujita and Y. Gohda, Phys. Rev. Appl. 11,024006 (2019).

79

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-Principles Design of
Spin-Polarized Bandstructure at MoTe2/Bi2Se3 Interface
Kunihiko Yamauchi
Institute of Scientific and Industrial Research (ISIR-SANKEN), Osaka University, 8-1
Mihogaoka, Ibaraki, Osaka 567-0047

andmay merge. We show that MoTe2 and Bi2Se3

topological properties in various materials bylattices match well and in-plane strain provides

means of the first-principles density functionala handle to manipulate the spin-polarized band

theory (DFT) calculations aiming at thealignment near the Fermi energy. This interface

spintronics application.

shows the various spin degree of freedom, and it

We

are

studying

magnetism

Since last year, we have investigated the

may be useful as an advanced spin generator [2].

magnetic stability in 2D magnetic material VI3

In this year, we also studied the cathodemaking a comparison with CrI3 using the VASP

properties of sodium transition-metal oxides

code at the ISSP ohtaka supercomputer. In VI3,

(such as Na2NiTiO4) by means of DFT

it turned out that under the trigonal crystal fieldcalculations and the machine-learning approachwith D3d point group, the t2g orbital state is split

to reduce reliance on Li-ion batteries and aim at

into a1g and eg’ states. The a1g orbital state plays

the sustainable development goals. We found aan important role in inter-layer magnetic

strong relation between the cation ionic radii andexchange interaction in a bilayer structure. By

the structural stability. This finding may lead tousing Wannier90 code that is interfaced with thethe future industrial application to the sodium-

VASP code, we estimated the hopping integral

ion rechargeable battery [3].

and discussed the inter-layer magnetic exchangeinteraction

based

on

the

References

virtual-hopping

mechanism. This finding may pave a way to the[1] T. P. T. Nguyen, K. Yamauchi*, T. Oguchi,

spintronics application with 2D materials [1].

D. Amoroso, and S. Picozzi, Phys. Rev. B 104,

Next, we explored a heterostructure that

014414 (2021).

consists of transition-metal dichalcogenide and[2] K. Yamauchi, R. Shimazu, and T. Oguchi, J.

3D topological-insulator to investigate thePhys. Soc. Jpn. 91, 044705 (2022).

hybrid spin state in which the out-of-layer spin[3] K. Yamauchi, H. Momida, and T. Oguchi, J.

polarization and the helical spin polarization

Phys. Soc. Jpn. 91, 034704 (2022) .

80

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles lattice-dynamics calculations
for negative-thermal-expansion materials
Yasuhide Mochizuki,1 Soungmin Bae,2 and Hannes Raebiger2
1

Department of Materials Science and Engineering,Tokyo Institute of Technology, 2-12-1 O-okayama, Meguro-ku, Tokyo
2

Department of Physics,

Yokohama National University, Yokohama 240-8501
Negative-thermal-expansion(NTE)materials have attracted considerable attention
because they contribute to sustaining largely
expanding materials and also because of theirfascinating

andcounter-intuitive

thermalproperties. The well-known and representative
NTE materials are Si [1], ReO3 [2], ScF3 [3]ZrW2O8 [4] and KZr2P3O12 [5]. Since the reportsof negative thermal expansion in various
materials in the past few decades, it has been
widely known that the mechanism for NTE is thetransverse acoustic phonon modes, which havenegative Grüneisen parameters. However, thedetailed mechanisms for the NTE properties are
uncovered ever since now. Here, we report ourprogress for elucidating the microscopic origin

Fig 1. Calculated phonon bands for (a) Si, (b)of NTE in Si, ReO3, and KZr2P3O12.
Our

first-principles

calculations

ReO3, and (c) KZr2P3O12 and Grüneisen

wereparameters for (d) Si, (e) ReO3, and (f)performed by using the projector augmented-

KZr2P3O12. The red and blue respectively

wave method [6] as implemented in VASP

indicate the degrees of positive and negative[7]. The phonon band structures were derived

values of Grüneisen parameters.

from the calculated force constants using
calculation results of phonon band with

PHONOPY [8].

Grüneisen parameters for Si, ReO3, andAs shown in Fig. 1, we present the81Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of TokyoKZr2P3O12. From Figs. 1(a, b, d, e), one can seethat the zone-boundary phonon modes in Si and
ReO3 generate negative Grüneisen parameters,
triggering the NTE. In Si, the most importantphonons for realizing NTE are the atomic
vibration eigenvectors which transform as
irreducible representations X4 and L3–. On theother hand, in ReO3, the most important phonons
for the NTE are the eigenvectors transforming as+

Fig 2. Theoretical and experimental thermal–

irreducible representations M2 and R5 . The+

expansion coefficients for KZr2P3O12. The blue–

phonon modes of M2 and R5 correspond to theand red colors respectively indicate the results

ReO6 octahedral rotations. Note that these zone-

from calculations and experiments. The redboundary phonon modes of Si and ReO3 are low-

broken line is derived by least-square method.

energy transverse acoustic phonon modes. It is

References

also noteworthy that transverse acoustic phonon[1] S. Biernacki and M. Scheffler, Phys. Rev. Lett.

modes do have crystal momentums, which

63, 290 (1989).

clearly indicate that they must be Brillouin-[2] T. Chatterji et al., Phys. Rev. B 78, 134105

zone-boundary phonon modes. Very intriguingly,(2008).

however, as shown in Figs. 1(c, f), a Brillouinzone-boundarymode

transforming[3] B. K. Greve et al., J. Am. Chem. Soc. 132,as

15496 (2010).

irreducible representation Γ2+ in KZr2P3O12[4] C. Martinek and F. Hummel, J. Am. Ceram.

contributes to NTE the most. These results

Soc. 51, 227 (1968); T. A. Mary et al., Sciencesuggest that not all negative Grüneisen

272, 90 (1996).

parameters are from transverse acoustic phonons.

[5] W. Miller et al., J. Mater. Sci. 44, 5441

Further, to confirm the validity of our(2009).

calculation results for KZr2P3O12, we also

[6] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).

prepared the polycrystal of KZr2P3O12 and[7] G. Kresse and J. Furthmüller, Phys. Rev. B

examined the thermal expansion coefficient byconducting

X-ray

diffractionin54, 11169 (1996); G. Kresse and D. Joubert, ibid

relevant

59, 1758 (1999).

temperatures. As shown in Fig. 2, our[8] A. Togo and I. Tanaka, Scr. Mater. 108, 1

calculations are qualitatively consistent with the(2015).

experimental results.

82

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Analyses on electronic structure and magnetoelectric effect inhigh-performance spintronics and magnetic materials
Tatsuki ODA1,2, Ko HYODO2, Daiki YOSHIKAWA2, Daichi TAKAGI2, Masao OBATA1,2
1

Institute of Science and Engineering, Kanazawa University, Kanazawa, Ishikawa 920-1192
2

Graduate School of Natural Science and Technology, Kanazawa University, Kanazawa,
Ishikawa, 920-1192
We have studied several systems in high-was evaluated using the total energy (TE) andperformance spintronics and magnetic materials

grand-canonical force theorem (GCFT) [1,2].

on electronic structure and magnetoelectric

Our self-consistent calculations revealed that theeffect. In these investigations, we took someintra-c-plane AFM configuration was stabilized

attempts for improving physical accuracy oras the ground state for both of the materials. The

finding new candidates in the materials of

MCAE is obtained to be negative, indicating a

spintronics application. Here, we report thefavor of in-plane anisotropy. The absolute value

developments of (A) investigation on magneticof MCAE tends to decrease as the ratio of latticeanisotropy energy (MAE) in antiferromagnet

constants c/a increases. The SMAE is always(AFM) L10-MnX (X=Ir, Pt), (B) optimization inpositive, reflecting the AFM configuration. In

a computational code for quasi-particle self-

MnPt, we obtained an anisotropy transition inconsistent GW (QSGW) method, and (C)the investigated range of strain, namely, theinvestigation on electronic structure in magnetic

whole MAE became positive (perpendicular toshape memory alloy (MSMA) Ni2MnGa.

c-plane) at larger c/a, from the in-plane(A) MAE in AFM MnX (X=Ir, Pt)

anisotropy. Through this investigation, theWe investigated the control of magnetic

approach GCFT was confirmed to be an efficient

properties by structural strain using the first-

approach to analysis of MAE.

principles calculations based on spin-density(B) QSGW method

functional theory (SDFT) in L10-type Mn alloy;

It is well known that the DFT suffers from

MnIr and MnPt. Contributions to the MAE weremany deficiencies such as small band gap,

calculated from both the magnetocrystalline

electron over binding property, etc.. To improveanisotropy energy (MCAE) arising from thethese deficiencies, extension methods of DFT

spin-orbit coupling and the shape magnetic

such as DFT+U have been developed. However,anisotropy energy (SMAE) arising from thethey often have a semiempirical character, andmagnetic dipole interaction energy. The MCAE

are not universal and parameter-free. In the other83

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

decrease in parallelization efficiency was foundto be due to the communication time among PEs.
(C) Electronic structure in Ni2MnGa
Recent results on Ni2MnGa show that themodulatedmartensite

phase

is

more

energetically stable than the non-modulated(NM) phase, supporting an importance ofelectron-correlation effect [4]. The GGA+U

Fig. 1: Speedup factor in the new code.

reveals that the stability of the austenite andhand, QSGW is a universal approach to themartensite phases is dependent on U. We

reliable prediction of electronic structure [3].

investigated the electronic structure of MSMA

The idea of QSGW is an optimization of the one-Ni2MnGa using QSGW [4,5]. The Ni eg orbitals

body Hamiltonian by performing a self-

in the cubic (austenite) phase, which lead toconsistent

many-bodymartensite phase transition by band Jahn-Teller

approach with GW approximation, where theeffect, were found to locate on the Fermi level.

self-consistency is constructed to minimize the

The NM and 4O structures of Ni2MnGa were

perturbation. QSGW has already had success incalculated with QSGW. The peak of density ofthe estimation of band gaps and other properties

states near the Fermi level for the minority-spin

for a wide range of materials. The requirement

channel was located above the Fermi level into be paid is computational cost. QSGW needs

NM and below the Fermi level in 4O. Thisroughly 100-1000 times more expensive in

difference might be related to the presence ofcomputational time than DFT. We parallelized

nanotwin double layer in modulated martensite

the part of calculation for the polarization

phases.

function (PF) in GW approach, which is one ofReferences

the most expensive (time-consuming) parts in[1] I. Pardede et al., J. Magn. Magn. Mater., 500,

QSGW. As a result, the computational time of

166357 (2020).

PF was reduced from 11 hours to 1 hour in the[2] I. Pardede et al., Crystals 10, 1118 (2020).

application to 4O structure (16 atoms in the unit

[3] T. Kotani, J. Phys. Soc. Jpn. 83, 094711

cell) of Ni2MnGa. The parallelization efficiency,(2014); ecalj: https://github.com/tkotani/ecalj/.

as shown in Fig. 1, was almost ideal when the[4] M. Zelený et al., Materials&Design, 209,

number of processor elements (PEs) was small,

109917 (2021).

but decreased as the number of PEs increases;[5] M. Obata et al., submitted, (2022).

calculation

in

the

with an efficiency of 60% at 30 PEs. This

84

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Prediction of new crystal structures based on the densest ternary sphere packings

Ryotaro Koshoji and Taisuke Ozaki
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Densest unary sphere packing was proven to

P6/mmm symmetry, and if a line of small spherebe the FCC structure in the 2000s, but thein the (4-3-1) structure is replaced of one small

densest binary and ternary sphere packings havesphere, we can get the crystal of 𝑆𝑚𝐶𝑜5. Note

not yet been identified analytically. To discoverthat we name a PDTSP consisting of 𝑙 small, 𝑚

the unknown putative densest packings, we havemedium, and 𝑛 large spheres per unit cell as (𝑙-created an algorithm to explore them, and we𝑚-𝑛) structure, and the small, medium, and large

successfully found the 12 putative densest

spheres in figures are represented by pink,binarywithyellow, and purple balls, respectively. The (9-7-confirming the known PDBSPs with a small3), (10-6-3), and (13-3-1) structures can becorrection of the phase diagram [1].

derived from the Perovskite structure, and they

spherepackings(PDBSPs)

On top of the success of the discovery, we

have the 𝑃𝑚3̅𝑚 symmetry. Remarkably, an

planned to explore the putative densest ternaryedge of unit cell of the (9-7-3) structure issphere packings (PDTSPs). If we only search theconstituted by a line of four medium spheres.

periodic packings of 𝑛 different sizes of spheres,

The (13-2-1) has the 𝐹𝑚3̅𝑚 symmetry. Largethe packing fraction at a given composition ratio

spheres constitute the FCC structure without

is maximized by filling the space with not morecontact, and a tetrahedral site is occupied by one

than 𝑛 kinds of packing structures. Since no

medium sphere. The medium and large spheres

seminal study seems to construct the phaseare surrounded small spheres constituting semi-

diagram, we made the ternary phase diagrams toidentify the PDTSPs that are denser than any
other combination of packing structures. The
example of the ternary phase diagram is shownin Fig. 1.
As a result, we have discovered the 59 kindsof PDTSPs [2,3]. The newly found PDTSPs
have unique structural properties, and some of

Fig.1: An example of ternary phase diagram

them have high symmetry. The (4-3-1) has the85Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

regular polyhedrons. This structure is predicted

hydrides that have the isotypic crystal structures

to become the PDTSP at the several radius ratios

to the (13-2-1) and (13-3-1) structures, but thesuch as 0:44 : 0:64 : 1:00, at which small andcalculations have not yet been finished.

medium spheres would be too large to be placed

The computational cost of the structuralin the tetrahedral and octahedral sites if large

optimizations is low, so all the process can be

spheres comprising the FCC structure withfinished in one processor, so the exploration ofcontacts. Finally, large spheres in the (16-2-2)

the structures was effectively parallelized.

structure are surrounded by the network

References

structures of small and medium spheres.
The six exemplified PDTSPs might be used

[1] R. Koshoji, M. Kawamura, M. Fukuda and T.

as the structural prototypes for crystals under

Ozaki: Phys. Rev. E 103 (2021) 023307.

pressure. For example, the (13-2-1) and (13-3-1)[2] R. Koshoji and T. Ozaki: Phys. Rev. E 104

structures might be realized by metal hydrides(2021) 024101.

with replacing small spheres of hydrogen atoms.

[3] R. Koshoji and T. Ozaki: arXiv:2110.15505

We have been exploring the quartic metal

Fig. 2: (4-3-1)

Fig.3: (9-7-3)

Fig.4: (10-6-3)

structure

structure

structure

Fig.5: (13-2-1)

Fig.6: (13-3-1)

Fig.7: (16-2-2)

structure

structure

structure86

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Low Thermal Conductance Generated by van Der Waals
Interaction
Masato OHNISHI
Department of Mechanical Engineering,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan

Clathrate compounds are promisingphonon (SCP) theory [1] and unified

candidates

theory [2]. The SCP theory is used toforthermoelectricmaterials in terms of the “electron-

considercrystalanharmonic potential and to obtainphonon-glass”

concept.

effectsofforth-order

Clathrate compounds are composed

temperature-dependentof guest atoms encapsulated ininteratomic force constants (IFCs).

cagelike structures. It is expected

Using the unified theory, one canthat electrons can transport in theconsider the coherent contribution,framework

lessthat of interbranch tunneling, toforheat transport as well as thephonons because of guest atoms

population contribution, which isleading

phononconsidered also in the context of

scattering. Particularly, in type-I

Peierls Boltzmann transport theory.

clathrate compounds, guest atoms

We

vibrate in a strong anharmonic

package [1] to perform the above

potential because they are composed

approaches and VASP to calculate

of a Weaire-Phelan structure, which

forces on atoms.

scattering

with
rate

to

much
than

that

significant

divides a space with the maximum
volume

with

the

same

have

employedharmonic

ALAMODE

While phonon dispersion haveimaginarycross-(negative)

frequencies

sectional area (Fig. 1(a)). In thiswhen guest atoms have negative

study,

harmonic

wehavestudies

type-I

potential

atlow

Ba8Ga16X30 (BGX, X = Si, Ge, Sn) totemperatures, negative frequencies

gain insights of the effect of phonon

becameanharmonicity on heat transport.

temperatures, as shown in Fig. 1(b).

Phononproperties

positiveatfinitewere

For type-I BGSn, harmonic potential

calculated with the self-consistentwas positive above 100 K, and87

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Fig. 1 Phonon properties of type-I clathrate compounds. (a) Crystal structures of typeI clathrate. (b) Phonon dispersion of Ba8Ga16Sn30 (BGSn) without considering
renormalization (grey), at 300 K (blue), and at 600 K (red). (c) Temperature-dependent
thermal conductivities of BGSn. Population (blue), coherent (orange), and electron(green and red) contributions were considered.

thermalconductivities

could

be

temperature-dependence of type-I

calculated at this temperature range.

clathrate at low temperatures (≈ 10

The

thermalK) have been widely studies, we haveall

revealed

calculated

conductivities

includingtheintriguing

phononcontributions (population, coherent,

properties at high temperatures inand electron contributions) were in

type-I clathrate [3].

good agreement with experimental
data, as shown in Fig. 1(c). To gain

References

insights into the effect of phonon[1] T. Tadano and S. Tsuneyuki, J. Phys. Soc. Jpn.

anharmonicity,

87, 041015 (2018).

we

have

also

analyzed mode dependent phonon[2] M. Simoncelli et al., Nat. Phys. 395, 1 (2019).

properties at high temperatures in[3] M. Ohnishi et al., arXiv:2103.00413.

details (> 100 K). While a unique

88Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

A Unified Understanding of Thermodynamic
Properties in Icy Planets and Icy Moons
Yu KOMATSU†‡ , Terumasa TADANO ¶ , Yasunori HORI†‡†
Astrobiology Center, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan.
‡
National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan.
¶
Research Center for Magnetic and Spintronic Materials,
National Institute for Materials Science, Tsukuba 305-0047, Japan.
The high-pressure properties of hydrogensulfide have drawn considerable attention interms of its superconductivity and metalization [1]. Astronomical observations have detected H2 S from the atmosphere of icy planets(ice giants) like Neptune and Uranus [2, 3], andalso in plumes coming from the interior oceanof Enceladus, a Jupiter icy moon [4]. Thesefindings provide insight for estimating the interior of icy bodies and their physical and chemical origins. Although H2 S would be one of thetypical minor components of ice mantle in icy
planets and satellites inside and outside oursolar system [5], the phase diagrams and properties of H2 S and its mixture have not yet been
revealed, especially at high P -T conditions.
Our project aims to elucidate possible mixing states and thermodynamic conditions in
the icy bodies by ab initio calculations. We
focus on the properties of H2 S as the first onset from the interest above and perform density functional molecular dynamics (DFT-MD)
simulations for H2 S in the temperature rangeof 500 ∼ 10,000 K and the density range of 1
∼ 5 g/cm3 using the Quantum ESPRESSO
package. More than 10 ps for MD simulations
are allowed after thermal relaxation with a 2 ×
2 × 2 supercell containing 32 H2 S molecules,
and the ion temperature is controlled in the
N V T ensemble with an energy cutoff of 80 Ry.
Pressures in our calculations were up to ∼250 GPa, corresponding to the values at the

Figure 1: Radial distribution function gHS for
Pmc21 structure.

deep mantle of ice giants [6]. Figure 1 presents
the radial distribution function between H andS atoms for Pmc21 phase of H2 S. The fluid
phase emerged at higher temperatures (e.g.,∼3,000 K with 3 ∼ 5 g/cm3 ), and at lower
densities like 1 g/cm3 , the phase appeared atmore higher temperatures. Moreover, diffusion coefficients in Figure 2 implicate the transition from the solid phase to the superionicand then to the fluid, like H2 O and NH3 [6].
The difference in diffusion between hydrogenand sulfide becomes more prominent at higher

89

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[6] Cavazzoni, C.; Chiarotti, G. L.; Scandolo,S.; Tosatti, E.; Bernasconi, M.; Parrinello,
M. Science 1999, 283, 44–46.

Figure 2: Diffusion coefficient of H and S for
Pmc21 structure.
temperatures (all the conditions are not plotted in Figure 2). Thermodynamic quantities,e.g., free energy, are being evaluated, and an
equation of state (EOS) including H2 S is being constructed, which are significant for understanding conditions realized in icy bodies.
The computation was partly carried out using the computer resource offered under the
category of General Projects by the ResearchInstitute for Information Technology, Kyushu
University.

References[1] Li, Y.; Hao, J.; Liu, H.; Li, Y.; Ma, Y. The
Journal of Chemical Physics 2014, 140,
174712.
[2] Irwin, P. G. J.; Toledo, D.; Garland, R.;
Teanby, N. A.; Fletcher, L. N.; Orton,
G. A.; Bézard, B. Nature Astronomy 2018,
2, 420–427.
[3] Irwin, P. G. J.; Toledo, D.; Garland,R.; Teanby, N. A.; Fletcher, L. N.; Orton, G. S.; Bézard, B. Icarus 2019, 321,
550–563.
[4] Waite, J. H., Jr; et al. Nature 2009, 460,
487–490.
[5] Johnson, T. V.; Mousis, O.; Lunine, J. I.;
Madhusudhan, N. The Astrophysical Journal 2012, 757, 192.

90

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles transport calculations for single-molecularjunctionsTatsuhiko OHTO
Graduate School of Engineering Science,
Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka 560-8531
Single-molecule devices attract attentions

the MD trajectory. We observed that the peak top

because of possible bottom-up construction,

value largely depends on the molecular junction

potential

structures.

fordeviceminiaturization,andutilization of unique electrical responses oforganic materials. Since it is not possible todirectly observe the bridging structure of themolecule between electrodes, the presence of themolecule should be confirmed through theconductancehistogram.

The

conductancehistogram is normally broad, ranging orders of
conductance values, while conventional firstprinciples transport calculations are performed

Fig. 1: Snapshots of the molecular dynamics

only for the most stable structure, which is

trajectory.

assumed to be the peak top value of the
conductance histogram. A combined molecular
dynamics and transport calculation is required toreproduce the broad conductance histogram anddiscuss the structure-property relationship.
In this work, a reactive force field to simulate
the break junction procedure has been developed.

Fig. 2: Calculated conductance histograms for

Three junction structures (on top, bridge, and

Au-BDT-Au molecular junctions.

hollow) were prepared and 200 break junctions
for each structure were simulated using theReferences

molecular dynamics (MD) method. We further[1] A. R Rocha et al. Phys. Rev. B 73, 085414

performed first-principles transport calculations(2008).

using the SMEAGOL codes for the snapshots for91

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles study ofquantum transport in nanostructuresNOBUHIKO KOBAYASHI
Department of Applied Physics, University of Tsukuba
1-1-1 Tennodai Tsukuba Ibaraki 305-8573

An atomic-scale understanding of quantum
transport properties in materials is an important issue in materials science and nanotechnology. The aim of this project is to reveal
and predict quantum transport properties inmaterials from first-principles. We have developed the nonequilibrium Green’s function(NEGF) method, and the O(N) time dependent wave-packet diﬀusion (TD-WPD) methodon the basis of the density functional theory(DFT). By using these methods, we have investigated charge, head and spin transport properties of materials. [1]We have developed a Simulation code for
Atomistic Kohn-sham Equation (SAKE) for
the ab-initio electron transport calculation
based on the DFT and NEGF formalism. [2]We have applied the method to analyses ofthermoelectricity of Heusler compounds. The
electronic structures and thermoelectric properties are analyzed, and high performance systems are designed.

Figure 1: Schematic pictures of the C8 DNBDT monolayer organic field-eﬀect transistor.
attention, since the ionic liquids induce an intense electric field at a semiconductor interface, resulting in a high concentration of charge
carriers. We elucidated that strong electricfields and highly doped charges induce significant modulation eﬀects on the electronic band
structure and the mobility.

We developed the O(N) TD-WPD method
for the quantum transport calculation of hugesystems of up to 100 million atoms a decadeago. It is possible to calculate the conductivity and the mobility of the system withmicron-order lengths at room temperature at
the atomistic levels. We have recently expand
the method to analyze thermoelectric properties, and demonstrated Seebeck coeﬃcient andpower factor calculations for a huge system of
100 million atoms.

References[1] N. Kobayashi, H. Ishii, K. Hirose, in 3D Local Structure and Functionality Design of Materials eds H. Daimon, Y. C. Sasaki (World
Scientific, 2018) p133.

Furthermore, we have analyzed the gate induced modulation eﬀects of electronic statesin monolayer organic field-eﬀect transistor as
shown in Fig. 1. [3] Organic field-eﬀect transistors with ionic liquids have attracted much[2] H.Takaki, N.Kobayashi, K.Hirose,
J. Phys.: Condens. Matter. 32 325901(2020)[3] H. Ishii, N. Kasuya, N. Kobayashi, K. Hirose,
S. Kumagai, S. Watanabe, and J. Takeya,
Appl. Phys. Lett. 119 223301 (2021).
92

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles study on charging-induced structure phasetransitions of thin film oxides
Takashi NAKAYAMA
Department of Physics, Chiba University
1-33 Yayoi, Inage, Chiba 263-8522
HfO2 has various crystal structures dependingon the growth conditions. Although the orthorhombic (Ortho) HfO2 is thermodynamicallyunstable, it can be grown on metal electrodesubstrates such as TiN and Pt and has giant
ferroelectricity, thus being expected as a promisingmaterial for future memory devices. However, it
has not been clarified why such unstable Orthophase is realized in experiments. In this work, we
study the stabilization mechanism of Ortho HfO2
phase using the first- principles calculations.
Since the ferroelectric phase of representative
BaTiO3 is normally realized by the substrateinduced strain, we first consider the effect of strain.
Fig.1 shows formation energies of stablemonoclinic (Mono) and unstable Ortho phases inbulk and thin-film forms as a function of thecross-section area, i.e., biaxial strain. The energy
increase in thin films occurs due to the surface
energy loss. As seen in both figures, by applying
tensile or compressive strain, the Mono phaseincreases the energy compared to the Ortho phase.
This is because the Mono phase is stabilized bylowing its crystal symmetry and the deformationfrom such stable structure promotes large energy
loss. However, the Ortho phase always has larger
energy than the Mono one. Thus, the strain is not asingle origin of the Ortho-phase formation.
Next, we consider the effect of oxygen vacancy(Vo), which is a representative defect and is often
generated in oxides. Vo produces an electronoccupied state below the conduction bands, thus
stabilized by the positive ionization. Fig.2 showsthe energy difference between Ortho and Monophases as a function of Vo ionization charge. From
this figure, we know that the Ortho phase becomesFig.1. Calculated formation energies of Mono andOrtho HfO2 in bulk (left) and thin-film (right)
forms under (111) biaxial strain, relative to themost stable Monoclinic bulk phase, as a functionof an area of formula unit.

Fig.2. Calculated difference in energy between
Ortho and Mono phases in bulk forms, as afunction of the ionization charge of oxygen
vacancy (Vo). The zero energy corresponds to the
Mono phase.

the most stable when the 4-5% Vo’s exist and they
are +2 charged. However, considering that theformation energy of Vo, is relatively large around
4-5eV, the assumed composition 4-5% is still large.
Thus, the ionization of Vo is also not a single originof the Ortho-phase formation.

93

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Fig.4. Calculated local-potential profile around
Ortho HfO2 (001) thin film. To correct theperiodic boundary condition of the repeatedslab geometry, the artificial electrodes areinserted in vacuum region so as to keep thehorizontal potential.

calculation method using the thin film geometry
should be developed. Thus, we try to calculate thepolarization charge using the conventional
repeated-slab unit cell.
Due to the cancellation of electric fields causedby the front and back surfaces of the film, thepotential should be flat outside the ferroelectric film.
To realize such situation, we arrange the artificial
electrodes (dipole) in the vacuum region in the unit
cell, as shown in Fig.4. Then, we calculate thepolarization charge using the Ortho HfO2 (001) thinfilm. The calculated polarization charge is about
53μC/cm2, which agrees with the values by the
Berry phase method around 40-50μC/cm2 and theobserved values around 40μC/cm2. In the future,
we will use this method and study how thepolarization changes with varying the substrate.
All these calculations were performed using thexTAPP, VASP, and our original pspwf codes. In
order to realize the calculations for the present
thin-film systems, because the system is made of alarge number of atoms (1000-2000 atoms), theadvanced computing facility having multi-task andhigher-speed CPU (more than 1280 cores × 2.0GHz), larger-size memory (around 256GB), andlarger-size strange (more than 2 TB) is indispensable. These conditions are realized only by thepresent ISSP supercomputing systems.

Fig.3. Calculated difference in energy between
Ortho and Mono phases in thin-film forms,
with/without 1-4% strain and 3% charged
oxygen vacancy (Vo). The cases of (a) (111)tensile strain and (b) (001) compressive strain
are shown.

Then, we consider the case when both the strain
and the Vo ionization coexist in HfO2. Fig.3 showsthe energy difference between Ortho and Monophases for the cases with/without strain and Vo
charging. Figures 3(a) and 3(b) respectively
correspond to the cases of (111) tensile strain and(001) compressive strain. It can be seen that boththe strain and the Vo charging work cooperatively
to stabilize the Ortho phase. In these ways, we canconclude that not only the substrate-induced strain
but also the ionization of Vo are important factors
to grow and stabilize the ferroelectric Ortho phase.
Finally, we consider the ferroelectricity of Ortho
HfO2. It is well known that the polarization chargeof ferroelectric material can be calculated by the
Berry phase method using the bulk unit cell. The
result by this method corresponds to the freestanding ferroelectric film located isolate in thevacuum. However, in experiments, such material isgrown on some substrates. In order to study theeffect of the substrate on the ferroelectricity, the94Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Elucidation of surface structure and reactive property of (100)exposed CeO2 nanocubes through first-principles calculation
Akira Yoko
Advanced Institute for Material Research,
Tohoku University, Katahira 2-1-1, Aoba-ku, Sendai, Miyagi 980-8577
CeO2 nanocubes with (100) facet exposed

Cr2O3 and CeO2 [3]. The Monkhorst-Pack k-have been synthesized by supercritical method,point mesh of 4×4×4 was used for the unit celland it was found to have high oxygen storageofcapacity. In this study, the stabilization of the

conducted with supercell of 2 × 2 × 2 extension(100) facet by organic modification was studied(96 atoms).The slab model exposing (100) facet

by using first-principles calculations.

includes 60 atoms, and that for (111) facet

Another approach to improve oxygen

CeO2.

The

defectsimulations

were

includes 108 atoms. The Monkhorst -Pack k-

storage capacity of CeO2 is metal doping into

point mesh for the slab model was 4×4×.

the structure. In this study, Cr-doping into CeO2

Figure 1a and 1b show CeO2 (100) and (111)

was also studied. The structural effects bysurface structure. The (111) facet of CeO2 isdoping were investigated.

nonpolar, contrary to the polar surface of (100).

In this study, first principles simulations

For the (100) facet model, half of the oxygen

based on density functional theory (DFT) was

was moved to the other side to exclude polarity

conducted. Plane wave basis DFT simulations

keeping the CeO2 composition. The surface

were conducted using VASP code. Structures ofenergy obtained from these calculations were

CeO2 at surface and bulk were calculated with

1.46 J/m2 for (100) and 0.78 eV/m2 for (111).

generalized gradient approximation proposed

The CeO2 (111) is more stable than CeO2 (100).

by Perdew, Burke, and Ernzerhof (PBE) [1].

The instability of (100) is the origin of high

The

the

activity as catalysts. When CeO2 is synthesized

pseudopotentials were 5s25p64f15d16s2 for Ce,

with organic surfactant, it is known that (100)3s23p63d54s1 for Cr, and 2s22p6 for O. The

facet exposes. In Figure 1c, structure of

DFT+U method introduced by Dudarev et al.

carboxylic acid attachment on CeO2 (100) is[2]was used to treat electron localization. The

shown. When calculating carboxylic attachment,

parameter U-J was set to 3.5 eV for Cr3d, 5.0

dipolar correction implemented in VASP waseV for Ce4f, and 5.5 eV for O2p states,used. After the deprotonation of the carboxylic

respectively, based on previous studies ofacid one oxygen of COO locates between twovalenceconfigurationsof95

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Ce atoms and the other oxygen locates on one

exists as 4+. Figure 2c shows the structure that

Ce atom. The adsorption energy was found to

one of the oxygen coordinating Cr (O-1 inbe -26.6 kJ/mol. The apparent surface energy

Figure 2b) was removed. As a result of the

includingvacancy introduction, local structure around Cr

theadsorptionenergyofthechanged to six-coordination, and Cr became 3+

surfactant is lower in (100) than (111).

state, and the oxygen vacancy formation energy
was 0.99 eV. The significant decrease inoxygenvacancyformationenergywasobserved. Also in this case, Ce3+ exits in thenext nearest neighbor of oxygen vacancy.
Further oxygen removal (O-2 in Figure 2c)occurred with lower oxygen vacancy formationenergy as 0.65 eV. The local structure and thevalence of Cr did not change, and two more
Ce3+ are formed at the next nearest neighbor of

Fig. 1 Surface structures of CeO2. a) CeO2
(100),b) CeO2 (111),

the oxygen vacancy. At the vicinity of Crc) Carboxylic aciddopant oxygen vacancy can be formed easily.

adsorption on CeO2(100).

Moreover, an oxygen atom, which is not at theThe bulk calculation was conducted tonext of Cr surrounded by four Ce4+ (O-3 inobtain the oxygen vacancy formation energy in

Figure 2d) was removed. As a result, thethe Cr-doped system. The supercell consisting

structure shown in Figure 2e was obtained. For

96 atoms were used. Figure 2 shows the local

the three oxygen vacancies, one Cr3+ and 5

structureThe

Ce3; was formed. Five Ce3+ locates not at the

information about the structures is summarizednearest neighbor (NN) of oxygen vacancy but

in Table 1. Figure 2a shows the pure CeO2,

at the next nearest neighbor of the oxygen

which includes an oxygen vacancy. Two Ce3+

vacancies. The oxygen vacancy in this case wasions exist in the next nearest neighbor (NNN)1.20 eV, and significantly lower than pureof the oxygen vacancy. The oxygen vacancy

CeO2 although the local structure was similar.

formation energy was 2.00 eV. Figure 2b shows

This is not the oxygen vacancies around dopant.

the structure of Cr doped CeO2, which one of

Doping affects whole the structure occurring

the Ce are substituted with Cr. The doping

distortion, which results in the lower oxygen

concentration is corresponding to 3 at%. In thisvacancy formation energy at far position fromcase four oxygen atoms coordinates Cr, and Cr

the dopant.

aroundoxygenvacancies.

96Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Figure 2
substituted

Structure of Cr 3% doped CeO2 a)Non-doped CeO2 with an oxygen vacancy, b) Cr4+
CeO2,c) ～ e)

Cr

doped

97

CeO2

with

1-3

oxygenvacancies

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Table 1 Structural information about Cr 3% doped CeO2Bader charge, |e|
Structure

Ce3+

Ce4+

Cr3+

Non-doped CeO2
with 1 VO+2.11 +2.44 n/a(Figure 6A)
Cr4+ doped CeO2
n/a
(Figure 6B)+2.45 n/a

Cr4+

O

VO
positionsE(distancefrom Vo
[eV]
Cr)

-1.28
to
n/a
-1.20
-1.24
+2.00 to
No VO
-1.16

Cr3+ doped CeO2
with 1 VO+2.15 +2.44 +1.86 n/a(Figure 6C)
Cr3+ doped CeO2
with 2 VO+2.15 +2.43 +1.85 n/a(Figure 6D)
Cr3+ doped CeO2
with 3 VO+2.13 +2.42 +1.85 n/a(Figure 6E)

Ce3+ positions
relative to VO

2.00

NNN (2Ce3+)

n/a

n/a

-1.26
O-1 in Figure 6B
0.99
to
(2 Å)
-1.17
O-1 in Figure 6B
-1.27
(2 Å)
to
O-2 in Figure 6C
-1.20
0.65
(2 Å)
O-1 in Figure 6B
(2 Å)
-1.37
O-2 in Figure 6C
to
(2 Å)
-1.20
O-3 in Figure 6D
1.20
(6 Å)

NNN (1Ce3+)

NNN (3Ce3+)

NNN (5Ce3+)

発表（投稿）論文等
1．Hao, X., Yoko, A., Inoue, K., Xu, Y., Saito, M., Chen, C., Seong, G., Tomai, T., Takami, S.,
Shluger, A.L., Xu, B., Adschiri, T., Ikuhara, Y.
“Atomistic Origin of High-Concentration Ce3+ in {100}-Faceted Cr-Doped CeO2 Nanocrystals”
Acta Materialia, 203, 116473 (2021)
DOI: 10.1016/j.actamat.2020.11.015
References[1] J.P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys.
Rev. Lett., 77 (1996), p. 38650.
[2] S. Dudarev, G. Botton, S. Savrasov, C. Humphreys, A. Sutton, Electron-energy-loss spectra andthe structural stability of nickel oxide: An LSDA+ U study, Phys. Rev. B, 57 (1998), p. 1505
[3] P.R. Keating, D.O. Scanlon, B.J. Morgan, N.M. Galea, G.W. Watson, Analysis of intrinsic
defects in CeO2 using a Koopmans-like GGA+ U approach, J. Phys. Chem. C, 116 (2012), pp.
2443-2452

98

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Magnetocaloric Eﬀect of Transition-Metal Alloys
Hung Ba TRAN1 and Tamio OGUCHI1,2
1
SANKEN, Osaka University, Mihogaoka, Ibaraki, Osaka 567-0047, Japan
2
CSRN, Osaka University, Machikaneyama, Toyonaka, Osaka 560-8531, Japan

3.5

"'

m (!B/Cr)"%!!"# ""#$

"$
"!
"#

0.5
00

H3456!"7
H3456$"7
H3456%"7
H3456&"7
H3456'"7"&
ï6S*"+,-./012!""#$%!#&$%!"#&$%'"#&$%

3
2.5
2()*+,#!!"# --#"#
1.5 *./0(,#! --#$
!"#
1

50

!"#(

100
T (Q)

150

200

ï!

"#

"'

)567.ï89:;5!-.7<5=
1,T3>K"ï,K"0T83#T

ï6S*"+,-./012

"%

!"#&

"$# "&# "(# ")# "!## "!$# "!&# "!(# "!)# "$##
T"+12

H3456!"7
H3456$"7
H3456%"7
H3456&"7
H3456'"7"&!"#')*+!,-./01#23

We study the magnetocaloric eﬀect oftransition-metal alloys by combining ﬁrstprinciples calculations and Monte-Carlo simulations. Several types of the magnetocaloric
eﬀect such as conventional (direct) magnetocaloric eﬀect, giant magnetocaloric eﬀect,
anisotropic magnetocaloric eﬀect, and inverse
magnetocaloric eﬀect can be considered by using our in-house program [1–6]. The peakin magnetization when applying an intermediate magnetic ﬁeld along the hard axis of
CrI3 , a well-known 2D material for the spintronic application, is successfully reproduced
by using our program as shown in Fig. 1.
The origin comes from the competition between magnetic anisotropy and external magnetic ﬁelds, where the temperature dependence of magnetic anisotropy energy plays
the role. The ﬁnite magnetic anisotropy atslightly higher than Curie temperature is reproduced, where the magnetization anisotropy
and anisotropic magnetic susceptibility is the
crucial factor. The high accuracy in the estimation of anisotropic magnetic entropy change
in CrI3 paves the way to predict the magnetocaloric properties of new materials. In thenext step, we have some ideas to enhance theaccuracy and predict some novel propertiesof magnetic materials within the theoreticalframework.

!!"# ""%

"$

!"#%

"!

!"#$
!"

"#

!("

!$""

T!,43

!$("

!%""

ï!

"#

"$# "&# "(# ")# "!## "!$# "!&# "!(# "!)# "$##
T"+12

Figure 1: Magnetization, magnetic anisotropy energy, and anisotropic magnetocaloric eﬀect of CrI3
[6].
[2] H. B. Tran, T. Fukushima, Y. Makino,
and T. Oguchi, Solid State Commun. 323,
114077 (2021).
[3] H. B. Tran, T. Fukushima, H. Momida, K.
Sato, Y. Makino, and T. Oguchi, Comput.
Mater. Sci. 188, 110227 (2021).
[4] H. B. Tran, T. Fukushima, H. Momida,
K. Sato, Y. Makino, and T. Oguchi, Direct and Inverse Magnetocaloric Eﬀects
in FeRh Alloy: A Theoretical Study,
http://dx.doi.org/10.2139/ssrn.3799705.
[5] H. B. Tran, H. Momida, Y. Matsushita,
K. Sato, Y. Makino, K. Shirai, and T.
Oguchi, Phys. Rev. B 105, 134402 (2022).
[6] H. B. Tran, H. Momida, Y. Matsushita,
K. Shirai, and T. Oguchi, Acta Mater.
231, 117851 (2022).

References[1] H. B. Tran, T. Fukushima, K. Sato, Y.
Makino, and T. Oguchi, J. Alloys Compd.
854, 157063 (2021).

99

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo!"#$%&'$(#)*$+",&#+$%*'-,+-+"./*/,0/*-1,2#-/*%.,30-$%*'-+",$+"$0"+%*'-,4,*)5"*$*%,/'"6+%*'-,)'2#" 7!"#$%&'()&*&!"#$%$&$'()*+(,*-%.(,$/$'(012#%3#4
51'(6"%7'+#%$2(*)(5*82*4(9/#1%:/;"*;1/4(9/#1%:/4(<1%=/(>??;@A@B!"# $%&"# '()*+"*# ","-(./-$"0+-%, ."%-(+/1'2'("0' )'+13# *"1'+(2# 4)1-(+/1%,# -/05+1"*#
6+($# +07,+-+(# '/,&%(+/1# 0/*",8# -%,,"*# 9:;<=>:; '+0),%(+/1' ?@ABCD#E$+'#2"%.8#6"#4/-)'"*#/1# ($"# %*'/.7(+/1# '(.)-()."'# /4# F(G@@@H<I!JG'""# K+3D# @%HD E$"# )1*".'(%1*+13# /4# ($+'#
0"(%,L6%(". +1(".4%-"# '2'("0 +'# +07/.(%1(# 4/.#5%'+-# ","-(./-$"0+-%,# ."%-(+/1'# ')-$# %'#/M23"1L$2*./3"1#."*)-(+/1L"&/,)(+/1#."%-(+/1'D
K+3D# @%# '$/6'# ($"# 0/'(# '(%5," '(.)-()."8#
F(G@@@H# NOPQNOP<I!J# +1-,)*+13# R< %1*# S<
0"05"."*# .+13'8# 6$+-$# +'# -/1'+'("1(# 6+($#"M7".+0"1(D K+3D#@5#'$/6'#($"#7%+.#*+'(.+5)(+/1#4)1-(+/1'# /4# ($"# /5(%+1"*# '(.)-()."8# 6$+-$#'$/6' ($"# F(<J# %1*# JAJ 7"%T' /4# UDO# V %1*

K+3D#@Y#G%H#X*'/.7(+/1#'(.)-()."#/4#F(G@@@H<I!JD#

UDR# VD# E$+' '0%,, JAJ# *+'(%1-"# -/),*# +1*)-"#

G5H#F%+.#*+'(.+5)(+/1#4)1-(+/1#/4#F(G@@@H<I!JD#($"# ."*'$+4(# /4# JI# &+5.%(+/1%,# 4."W)"1-+"'D# !"#%,'/# -/1*)-("*# ($"# %1%,2'+'# /4# &+5.%(+/1%,#

References4."W)"1-+"'# /4# ($"# /5(%+1"*# '(.)-()."D# X'# %#?@C# ;D# J(%1+# %1*# JD# :)3+1/Y# F$2'D# ="&D# Z# !"."'),(8 -%,-),%("* &+5.%(+/1%,# 4."W)"1-+"'# 6"."#

GU[[\H#@@RB[SD#+1#3//*#%3.""0"1(#6+($#"M7".+0"1('D I/6"&".8#?UC#:D#]+'$+$%.%#%1*#;D#J(%1+Y#F$2'D#="&D#Z##$($"# %*'/.7(+/1# "1".32# /4# ($+'# '(.)-()."# +'# ,"''#

GU[@SH#@@RBUPD($%1# ($"# %-(+&%(+/1# "1".32# /4# *"'/.7(+/1# 4./0#?OC# ^D# I%.)2%0%8# ED# >T"'$/_+8# %1*# ;D# J(%1+8#"M7".+0"1(8# 6$+-$ +1*+-%("' ($%(# ($"# F(<I!J#F$2'D#="&D#;%(".D#% GU[@`H#[PR`[@D+1(".%-(+/1# '$/),*# 5"# +07./&"*# 52# 0/."#?BC#^D#I%.)2%0%8#ED#>T"'$/_+8#%1*#;D#J(%1+Y#^D#'/7$+'(+-%("*#*"1'+(2#4)1-(+/1%,D

F$2'D#a$"0D#a#&%% GU[@`H#P`[BD

100

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Study of Efficient Training Data Generation Method
for Constructing Artificial Neural Network Force Field III
Kohei Shimamura
Department of Physics,
Kumamoto University, 2-39-1 Kurokami, Chuo-ku, Kumamoto 860-8555
Machine Learning Interatomic Potential

potentials is closely related to pressure. Thus, it(MLIP) not only achieves high accuracy byhas been found that without pressure training,learning

obviousfromfirst-principles

molecular

dynamics (MD) data, but also has low

errors

can

occur

in

thermalconductivity [2].

computational cost. Active Learning is the

Since the current MLIPs are trained not only

driving force behind the recent development ofwith potential energy but also with atomic force

MLIPs. This approach primarily uses theand pressures, the cost function is composed of

uncertainty in the output of MLIPs to discovertheand generate unlearned data for efficiently

training, adjustable coefficients are therefore

constructing robust potentials. In this study, we

introduced into these terms having different

continued to work on improving Active

units and data size. Although previous studies

Learning. The MLIP using artificial neuraldidnetwork (called ANN potential) was used. In

demonstrated that the coefficient adjustment

the previous studies, the use of standard

was a significant effect on the accuracy of ANN

deviation of atomic force is consideredpotentials [3].

effective for finding unlearned data. We havecorrespondingnotemphasize

terms.

theFor

efficiently

importance,weNo matter how efficiently untrained data canalso found that the method worked well.

be obtained, it becomes meaningless if theOn the other hand, we found that the MLIPs

training is inadequate. As the application oftraining method needs to be improved. The

MLIPs expands, the training of pressure and thecurrent standard training method does notcoefficient adjustment would become essential.

guarantee the accuracy of pressure required for

References

MD. We have demonstrated that the training of[1] A. Irie, et al.: J. Phys. Soc. Japan 90,

pressure is essential for ANN potentials to

094603 (2021).

reproduce the solid-liquid phase transition [1].

[2] K. Shimamura, et al.: Chem. Phys. Lett. 778,

Pressure training is also essential for thermal

138748 (2021).

conductivity calculations using MLIPs. The[3] A. Irie, et al.: J. Phys. Soc. Japan 91,

heat flux in the Green-Kubo formula for ANN

045003 (2022).

101

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Structure stability and electronic structures ofpoint defects in wide-gap oxide materials
Hiroyoshi MOMIDA
Institute of Scientific and Industrial Research (SANKEN), Osaka University
8-1 Mihogaoka, Ibaraki, Osaka 567-0047

The wide-gap oxide materials such as Al2 O3
and Ga2 O3 have been attracted great interests
in applications for power electronics devices.
In this work, we have studied the structuraland electronic properties of oxygen defects (vacancies or interstitials) in several amorphousoxides by the first-principles calculations.
The amorphous model structures of a-Al2 O3
and a-Ga2 O3 are generated by the firstprinciples molecular dynamics simulations andthe melt-and-quench techniques. Figure 1
shows an example of the generated model
structures using the 120-atom supercells. The
mass densities are set to be 3.3 g/cm3 for aAl2 O3 and 4.5 g/cm3 for a-Ga2 O3 , which arein their respective experimental ranges and arelower than those of the α crystal phases (4.0
and 6.4 g/cm3 for Al2 O3 and Ga2 O3 , respectively). The O coordination numbers around
Al are four (60.4%), five (29.2%), and six(10.4%) in a-Al2 O3 , and those around Ga are
four (66.7%) and five (33.3%) in a-Ga2 O3 .

(a) a-Al2O3O(b) a-Ga2O3

Al

Ga

0
O Al Oi

Figure 2: Density of states (DOS) of a-Al2 O3
with the excess O (O0i ). Total (black line withleft scale) and O0i (red line with right scale)DOSs are shown. The inset shows the local
structure near O0i , forming the peroxy linkage.
The O defects are introduced by adding andremoving O (interstitials Oqi and vacancies VOq )in the pristine oxide models, considering thecharged states q of defects. In the charge neutral state (q = 0), it is found that O0i tends
to form the O–O peroxy linkage, giving theπ ∗ - and σ ∗ -type O2 -like states near the valenceband top and the conduction band bottom regions, respectively (Fig. 2). It is also foundthat such the O2 -like states show interesting
behavior by changing q, and similar properties
of Oqi are confirmed in a-Ga2 O3 and in several
oxide crystals. We also study electric properties due to VOq in oxides to discuss ReRAM [1].

O

References
Figure 1: Amorphous structures of (a) a-Al2 O3and (b) a-Ga2 O3 120-atom supercell models.

[1] H. Momida and T. Ohno: Appl. Phys. Lett.
117 (2020) 103504.

102

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Theoretical study for magnetic-nano-particles formationfrom first-principles calculations
Yasutomi TATETSU
Meio University, Biimata, Nago, Okinawa 905-8585
There is a growing interest in searching novel materials that are not only thermodynamicallystable, but also thermodynamically unstable for applications. Chemical synthesize techniques, which
control nano particles precisely, are being developed in order to find new functional materials fortechnologies. Recently, Teranishi’s group at Kyoto University found a new behavior of crystal
structure transformation for an ionic nanocrystal [1], and a new Z3-type Fe(Pd,In)3 crystal structure,
which was synthesized by introducing a third element of In into L12-FePd3 [2].
In this report, we performed first-principles non-collinear calculations in order to explain thehigher coercivity value in the Z3-Fe(Pd,In)3 structure. We used the ab-initio code, OpenMX [3], forcalculating the magnetic anisotropy energies of L12-(Fe,In)Pd3 and Z3-Fe(Pd,In)3 structures andcompared these energies to find easy axes. In each calculation, the spin directions of each atom inthese systems were set from [001] to [100] when 𝜑 = 0 and from [100] to [110] when 𝜃 = 0. We
found that the magnetic anisotropy energies show −1.38 μeV/atom for the L12-(Fe,In)Pd3 structureand −0.213 meV/atom for the Z3-Fe(Pd,In)3 structure. The easy axis for Z3-Fe8(Pd20,In4) is the c axis,
as shown in the figure. The tendency which the Z3-Fe(Pd,In)3 structure has higher coercivity values
than the L12-(Fe,In)Pd3 structure is consistent with the experimental measurement. We assume that
added In atoms might be one of the reasons toincrease the magnetic anisotropy of the Z3Fe(Pd,In)3 structure. Further theoretical studies areneeded to explain the mechanism of the coercivity
increase in the Z3-Fe(Pd,In)3 structure.

References[1] Z. Li, M. Saruyama, T. Asaka, Y. Tatetsu, and T.
Teranishi, Science, 373, 332 (2021).
[2] K. Matsumoto et al., Nat. Commun., 13, 1047

Figure: Total energy analysis for the Z3Fe(Pd,In)3 structure depending on the spin

(2022).

directions.

[3] http://www.openmx-square.org

103

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Evolutionary algorithm for simulation of fastchemical reaction process
Koichi KUSAKABE
Graduate School of Science, University of Hyogo
Kouto, Kamigori, Hyogo 5-8581

1

Introduction

One methodology in the quantum materials informatics, based on density functional theory
to evaluate many-body correlation eﬀects, hasrecently been used to determine the electronic
phase diagram of oxides [1], to design hightemperature superconductors [2], and spinelectronic devices based on atomic layer material [3]. Materials design of measurementbased quantum computing devices [4, 5], andthe elucidation of unknown photo-induced
phase transitions based on measured data [6]were performed. We have also designed an
activation method of catalytic properties andevolutionary operations for detonation processes.

2

nanographene catalyst, in which multiple active sites are involved in the reaction, conﬁrm that propylene generation from propane
on our catalytic structure has the eﬀect of signiﬁcantly lowering the endothermic reaction
heat.
At the same time, in-situ AI analysis of timedependent small-scale data during the synthesis of catalytic materials was requested toelucidate the experimental procedures. Theoretical analysis of emission spectra during
the material conditioning (activation processby high-energy particle irradiation) was performed. As a result, it became possible to develop a measurement method to support insitu observation of molecular desorption by
adopting machine-learned classiﬁcation using
ResNet.

Alkane dehydrogenation on
3
nanographene catalyst

To support designed activated nanographene
catalysts for dehydrogenation of alkane, several dehydrogenation reaction pathways were
determined based on density functional theory
calculations. The catalytic eﬀect of vacancy
created graphene surfaces is simulated using
the alkane decomposition reaction as an example.
The eﬀect of the atomic defects was investigated, and the eﬀect of lowering the reaction energy barrier due to substitution by nitrogen and other group V elements was conﬁrmed. Furthermore, the characteristics of the

Algorithms for elucidationof detonation reaction paths

We are also promoting the elucidation of thedetonation reaction mechanism by combining
evolutionary algorithms, which reveals oxidation and reduction reactions in the TNT-RDX
composite that produces the detonation diamond. An evolutionary algorithm systematically generates intermediate states of chemical
reaction processes.
Stable crystals initiate rapid chemical reactions triggered by mass transport. We introduced mass transport into the evolutionary algorithm. The basis for the continuation of the104

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

reaction was found by structural optimization
calculations, ﬁrst-principles molecular dynamics calculations, and the migration paths ofradicals within the crystal.
Consider that radical species are supplied toadjacent material phases in the processes ofredox reactions. To determine the initiation
process of this chain reaction, we applied theaccelerated chemical reaction simulation using
the generational operation of the evolutionary
process. We applied this method to oxidationof RDX and reduction of TNT, and conﬁrmedthat our acceleration method works
The initiation point of the detonation process is considered to be the formation of OH
radicals by electron transfer excitation. The
OH radicals pull out hydrogen atoms andrapidly destabilize RDX. RDX then decomposes to form NO2 . (Figure 1 (a), (b))
Hydrogen radicals are supplied to TNT fromthe destabilized RDX. When a certain amountof hydrogen radicals are added, the nitro
groups of TNT are decomposed to generatewater, resulting in nitroso groups. This nitrosogroup functions to join adjacent carbon skeletons. (Figure 1 (c)) From these carbon skeleton linkages, the formation of clusters withlarger carbon skeletons may follow.
By conﬁrming that the supply of hydrogenradicals to TNT and OH radicals to RDX initiate rapid reactions, and by successfully drawing the initial reaction cycle, we conﬁrmed that
our method has potential for providing insight
into real-time evolution problems.

Acknowledgement
The author is grateful to Dr. T. Ishikawa, Mr.
Y. Makino, and Prof. S. Sakamoto for theircollaboration. This work is partly supported
by collaborative joint research in 2021 on “Research of detonation soot as a non-oxidative
dehydrogenation catalyst” and “Study of detonation reaction mechanisms combining molecular dynamics and evolutionary algorithms”

with Daicel engineering science joint researchin Osaka university.

Figure 1: Molecular structures arising from operations of evolutionary algorithms. (a) RDX,
which is destabilized by the addition of hydroxyl radical, generates (b) a water moleculeand simultaneously releases of NO2 . (c) TNT
destabilized by hydrogenation can generate
carbon clusters bridged by NO groups.

References[1] S. Teranishi, K. Nishiguchi, S. Yunoki, K.
Kusakabe, J. Phys. Soc. Jpn., 90, 094707(2021).
[2] S. Teranishi, K. Nishiguchi, K. Kusakabe,
J. Phys. Soc. Jpn., 90, 054705 (2021).
[3] H. Harfah, Y. Wicaksono, G. K. Sunnardianto, M. A. Majidi, K. Kusakabe,
Nanoscale Adv., 4, 117 (2022).
[4] N. Morishita, Y. Oishi, T. Yamaguchi,
K. Kusakabe, Appl. Phys. Express 14,
121005 (2021).
[5] N. Morishita, K. Kusakabe, Phys. Lett. A,
408, 127462 (2021).
[6] M. Nagai, Y. Higashitani, M. Ashida,
K. Kusakabe, H. Niioka, A. Hattori, H.
Tanaka, G. Isoyama, N. Ozaki, Researchsquare, DOI:10.21203/rs.3.rs-130295/v1(Preprint).

105

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Study on Hydration Process of Cellobiose by First Principles Calculations
Tomomi Kondo1,2, Daisuke Akazawa1,2, Takehiko Sasaki1，Motoyuki Shiga2
1 Department of Complexity Science and Engineering, Graduate School of Frontier

Sciences, The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8561
2 Center for Computational Science and E-Systems, Japan Atomic Energy Agency

148-4 Kashiwanoha Campus, 178-4 Wakashiba, Kashiwa, Chiba, 277-0871

Cyclizationof

30 water molecules, and acidic molecules suchpolyalcohol in hot water has been proven to be

as HCl added under the same thermodynamic

efficient method for utilization of biomass

conditions as in Experiment [1] (20 MPa, 523 K).

materials [1]. For this purpose we proposed a

CV is defined for each product as shown in Table

new

1. Fig. 2 shows the illustration for 1,4-AHSO3

methoddehydrationofreactionrefinement(improved

metadynamics method: Refined MTD),

andas one example. The dihedral angle φ between

applied it to cyclization dehydration of 2,5-

the hydroxyl group and the main chain, theHexanediol (HDO) [2]. Based on these results

difference d in the O-C distances of leaving and

we studied cyclization dehydration process of d-bonding, and the proton coordination number nsorbitol (SBT) yielding 1,4-Anhydro-d-sorbitolof the hydroxyl group were set according to the(1,4-AHSO),

2,5-Anhydro-d-sorbitol(2,5-

reaction path. Reaction processes in HCl(aq),AHSO)

1,5-Anhydro-d-sorbitol(1,5-

H2CO3 (aq) and pure water were examined.

andAHSO) as shown in Fig. 1.

Fig. 2 Illustration for CV of 1,4-AHSO3Table 1 Definitions of CVs
Fig. 1 Cyclization dehydration process of d-sorbitol.

Calculations were conducted using software
PIMD [3] working with DFTB or VASP. The

According to DFTB based calculations overcalculation system consists of 1 SBT molecule,3 ns dehydration reaction was observed and the106Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

transition state was found on the free energy

barriers for two of the pathways leading to thecontour map at CV with φ =155° , d=0.0 andfive-membered ring products, 1,4-AHSO and

n=2.1. Molecular dynamics calculations were

2,5-AHSO were found to be lower than those

conducted around this transition state so that theleading ot the six membered ring product, 1,5-

real time mechanism is visualized. After proton

AHSO. This is in agreement with experiments

is bonded on the oxygen atom (Fig. 3 ①) thewhere only the five-membered ring products

dehydration (Fig. 3 ②) and the formation of O-

were detected. The reaction pathway to the five-

C bond (Fig. 3 ③) occur almost simultaneously,

membered ring products goes through a

indicating that the five-membered ring is formed

metasbable state of SBT in the confined

via SN2 process.

conformation. These results were published in➀ 0 ps➁ 0.145 ps[4].
n = 2.1

We also calculated ab initio molecular
dynamics (MD) simulations for cellobiose
aqueous solution using PIMD and CP2K

➂ 0.240 ps➃ 0.498 ps

program and the subsequent C-K edge XAS
spectra simulation using CP2K. Experimental
XAS spectra can be interpreted by watermolecules interacting with the cellobiose
molecule [5].

Fig. 3 Reaction process from SBT to 1,4-AHSO

References[1] A. Yamaguchi et al., Green Chem., 13,873 (2011).
[2] T. Kondo et al., J. Comput. Chem., 42, 156(2021).
Fig. 4 Dehydration mechanism for SBT.

[3]http://ccse.jaea.go.jp/ja/download/pimd/in
dex.jp.html.

The mechanism is summarized in Fig. 4. The[4] T. Kondo et al., J. Comput. Chem., 42,

simulation results predict that the reaction is1783(2021).

faster in HCl solution than in carbonic acid[5] D. Akazawa et al., J. Chem. Phys., 156

solution, due to the acid-catalyzed mechanism.

044202 (2022).

In most of the calculated conditions, the reaction

107Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Ab initio calculation for structural stability oftransition-metal-oxide superlattice
Kazuma NAKAMURA
Quantum Physics Section, Kyushu Institute of Technology
1-1 Sensui-cho, Tobata, Kitakyushu, Fukuoka, 804-8550

0.20

D H (eV/f.u.)

In this study, we investigate the structural stability of transition-metal-oxide superlattices using first-principles calculations.
We consider superlattices built by polar andnon-polar transition-metal oxides. The polar type is a transition-metal oxide in which
the charge neutrality of the constituent parts
is not maintained, but the entire charge neutrality is maintained. For example, in a perovskite LaAlO3 , it is made as the stacking
structure consisting of a LaO plane charged
with +1 valence and an AlO2 plane charged
with −1 valence. On the other hand, in the
non-polar type, the charge neutrality is maintained including the partial system; for example, in SrTiO3 , it consists of charge-neutral
SrO and TiO2 planes. We consider 10 typesof superlattices consisting of SrTiO3 , SrVO3 ,BaTiO3 , LaAlO3 , CaCuO2 . In the present report, we focus on SrTiO3 /SrVO3 (STO/SVO)
and CaCuO2 /SrTiO3 (CCO/STO) superlattices. Also, we write the superlattice as An /Bn
where A and B specify materials, and n is their
number. We used xTAPP [1](a) STO/SVO

0.5

DH
DHI
DHS

0.15
0.10

(b) CCO/STO
DH
DHI
DHS

0.4
0.3
0.2

0.05

0.1
0.0

0.00
1

2

3

4

5
Number of layers n in [An/Bn]

1

2

3

4

5

Number of layers n in [An/Bn]

Figure 1: Calculated standard formation
enthalpy ∆H(n) (black lines) for the (a)(STO)n /(SVO)n and (b) (CCO)n /(STO)n superlattices. Red and blue lines are an interfacial contribution ∆HI (n) due to the interface
formation, and a strain-storage one ∆HS (n)due to the lattice mismatch.
and

 

∆HS (n)= E(A, akn ) − E(A) + E(B, akn ) − E(B) , (3)respectively. ∆HI (n) is an interface energy
due to the interface formation, and ∆HS (n) is
a strain-storage energy due to the lattice mismatch. Also, ak is an in-plane lattice parameter of the superlattice.
Figure 1 shows our calculated standard forTo discuss the structural stability of the su- mation enthalpy ∆H(n) (black lines), as wellperlattice, we calculate the n dependence of as its decomposition ∆HI (n) (red lines) andthe following standard enthalpy of formation ∆HS (n) (blue lines). Panels (a) and (b) are
the results for STO/SVO and CCO/STO su∆H(n) as
perlattices, respectively. The case ∆HI (n) >∆HS (n) indicates that superlattice does not∆H(n) = ∆HI (n) + ∆HS (n),(1)
develop epitaxially; the STO/SVO superlattice
is stable with the epitaxial growth, while thewhere ∆HI (n) and ∆HS (n) are defined as
CCO/STO superlattice is unstable.
About the optimized structure of the superE[An /Bn ] ∆HI (n)=− E(A, akn )+E(B, akn ) , (2) lattice, the following trends were observed:n108

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(a) Before relaxation(b) After relaxation

(a)

(b)
−

E
+

E=0
SrO
: Sr, Ba, La: Ti, V, Al

TiO2
Ca
CuO 2

E

:O

E=0

: Ca

: Cu

−
STO
Insulator
+
CCO
Metal
−

E
+

Figure 2: Ferroelectric displacement occurredin the CCO-based superlattice.
1. In the superlattices composed only of theperovskite-type oxides, an in-plane lattice
constant ak shows a tendency to converge
to an averaged value of the two materials
that make up the superlattice.

Figure 3: Left: Calculated effective potential for electrons in the superlattice, along thestacking direction z. Right: A schematic figure
for the resulting electric field.
Table 1: Interfacial charge in the CCO-based
superlattice. Unit of e.

2. The average c length of of each block constituting the superlattice basically reflects
the bulk parameters.

Interfacial charge
n

3

4

5

CCO/STO 0.33 0.30 0.28 0.27
CCO/BTO 0.33 0.33 0.31 0.29
CCO/LAO 0.38 0.31 0.28 0.26
LAO/SVO 0.25 0.20 0.19 0.16

3. In the superlattice of polar and insulating
non-polar transition-metal oxides, a significant ferroelectric displacement occurs
in the non-polar side (Fig. 2).
To understand the tendency of the standard
formation enthalpy (Fig. 1) and the ferroelectric displacement in the CCO-based superlattice (Fig. 2), we calculated an effective potential for electrons in the superlattice along a
stacking direction z. Figure 3 (a) displays results for the CCO/STO superlattice. From thefigure, we see that a potential gradient is occurring in the STO side (blue arrow), leading to
the electric field from TiO2 /Ca to CuO2 /SrO
interfaces. From the behavior of this effective
potential, the CCO/STO superlattice can be
regarded as dielectric inserted into a capacitor [Fig. 3 (b)]. With this view, the interfacial
charge can be estimated from the classical elec-

2

tromagnetism as
s
Q=

ǫ
ǫ0



ǫ0

dv(z)
∆S.
dz

(4)

Table 1 shows our estimated interfacial
charges for the CCO-based superlattice, weshow an n dependence. The resulting charge
is near 0.3 e, which gives the electric field of
0.87 VÅ−1 .

References[1] http://xtapp.cp.is.s.u-tokyo.ac.jp

109Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Study on Removal Mechanism of Single-Crystalline Si
Planarized by Catalyst-Referred Etching in Pure Water
Pho Van BUI, Taku HAGIWARA, Kiyoto KAYAO, Daisetsu TOH, Kouji INAGAKI,
Yoshitada MORIKAWA
Graduate School of Engineering,Osaka University, 2-1 Yamada-oka, Suita, Osaka 565-0871
There is a growing demand for ultra-

Our recent studies indicated that the

precision optical components for scientific andremoval mechanism of SiC and SiO2 etchingindustrial applications, especially in extreme

via CARE is a hydrolysis reaction, in which theultraviolet (EUV) and X-ray regimes. Using

catalyst

short-wavelength light, scientific imaging ofstabilization of the hypervalent state, increasing

cutting-edge materials/biological samples andthe rate of the reaction [3]. Thus, the purpose ofnanoscale lithography has become possible

the proposed research is to clarify the removal

thanksopticalmechanism of Si etching via CARE using first-

components. Si is the material of choice forprinciples calculations. Additionally, in Si

optical components, thanks to its attractive

polishing, the effect of self-oxidation byproperties,

to

the

such

ultra-precisionwater

dissociationandhighercleanliness,

oxygen and water is still not well understoodworkability.

A highly

due to the instability of Si under ambient

ordered surface with a root mean squareconditions. Thus, the oxidation effects ofroughness at the level of several tens ofoxygen and water are also taken into account.

machinability,

and

as

assists

picometer is greatly desired for the highest

The mechanistic removal pathway of CARE

reflectivity and the lowest unwanted scattering

for Si(111) has been clarified and shown in Fig.

[1].

1. In the pathway with the Pt catalyst, at first, a

In this study, we apply catalyst-referredwater molecule is dissociatively adsorbed on Pt.

etching (CARE), an abrasive-free polishing

In the next step, the adsorbed OH is transferred

method, to the planarization of a crystalline Si

to the targeted Si, forming the five-foldsurface, using Pt as the catalyst and pure water

coordinated state in the metastable state (MS).

as the etching solution [2]. An atomicallyAt this state, the Si-Si back-bond is elongatedsmooth surface with sub-Angstrom root-mean-

and weekend. Finally, the Si─Si back-bond issquare roughness is achieved on the Si

broken by a proton transfer at the final state

substrate.

(FS). The activation barrier is reduced from 1.0

110

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of TokyoeV (without the Pt) to 0.6 eV (with the Pt).

The obtained results indicated that the

Similar results are observed on a Si(100).

oxidation occurs readily in the air at a step-edgesite but not at a terrace site. The obtained
results are consistent with the previous study[4]. The obtained results also explained why thewettability of the Si surface could last forseveral hours in the ambient conditions in theprevious experimental study [5].

Fig. 1: Mechanistic reaction pathways of the Ptcatalyzed water etching of Si(111).
To investigate the effect of oxidation for an
HF-treated Si(100) surface, Si(100) with an Hterminated surface was used. Si-H and Si-OH
exchange reaction is investigated, as shown in
Fig. 4: Oxidation by an oxygen molecule at aFig. 2. The activation barrier is 0.75 eV.

step-edge and a terrace site.
The study has clarified the mechanistic
removal pathway of CARE for Si in pure water
Fig. 2: Oxidation effect of Si(100) by waterand the oxidation effect of a Si surface by

To investigate the effect of oxidation by

oxygen or water. Oxygen plays a main role inoxygen at the terrace and step-edge sites,

surface oxidation. Moreover, the combinedSi(111) is employed, as illustrated in Fig. 3.

oxidation effects of oxygen and water might be

The corresponding activation barriers are 0.46

an important factor to be considered in futureand 0.15 eV, respectively (Fig. 4).

studies.

References[1] R. J. Noll: Opt. Eng. 19 (1980) 249.
[2] H. Hara et al.: J. Electron Mater. 35 (2006) L11.
[3] P.V. Bui at al.: Jpn. J. Appl. Phys. 57 (2018)055703.
[4] F. A. Soria, E. M. Patrito, and P. Paredes-Olivera:

Fig. 2: Si(111) oxidation effects by oxygen at

J. Phys. Chem. C 116 (2012) 24607.

the terrace and step-edge sites (top-view).

[5] M. Niwano at al.: J. Appl. Phys. 76 (1994) 2157.

111

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Ab initio molecular dynamics study on ion-conduction
mechanisms of perovskite-type oxide
Satoshi OHMURA
Department of Environmental and Civil Engineering,
Hiroshima Institute of Technology, Miyake 2-1-1, Saeki-ku, Hiroshima 731-5193
We have investigated the mechanism ofcrosswise position from the potential valleys ofphotoinduced oxygen dynamics in cobalt

the ground state via intersection of the potential

double-perovskite

densityenergy curves in the excited and ground state atfunctional theory (DFT) based quantum-~1 Å. Figure 2 shows the ionic configurations

molecular-dynamics (QMD) simulations.

of the oxide ion (O1) as it moves from thecrystalusing

We adopted the LR-TDDFT for calculating

diagonal to the crosswise position. The oxide ion

light-absorption spectra, which are calculated as

is subjected to the repulsive and attractive forces

folded oscillator strengths

from the other oxide and cobalt ions,
respectively.

Fig. 2. Ionic configurations during
the travel of the oxide ion
Fig. 1. Calculated oscillator strength

Referencesof EuBaCo2O5.39 crystal[1] M. Hada, S. Ohmura et al.: Carbon. 18(2021) 612.

We found that the photoexcitation with an[2] M. Hada, S. Ohmura et al.: Appl. Mater.

energy exceeding 3 eV provides the oxide ion

Today 24 (2022) 101167.

with excess energy, arrowing it to pass throughthe barrier. The oxide ion may travel to the112

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Atomic structure of metal-oxide interfaces and oxide
grain boundaries
Kazutoshi INOUE1 , Mitsuhiro SAITO2 , and Yuichi IKUHARA1,2,∗
1Advanced Institute for Materials Research, Tohoku University, Sendai, Miyagi, 980-8577
2
Institute of Engineering Innovation, The University of Tokyo, Yayoi, Tokyo, 113-8656

The atomic structures of SiO2 /4H-SiC(0001) [Si face] and (0001̄) [C face] interfaces were analyzed by ﬁrst-principles molecular dynamics (MD) calculations combinedwith aberration-corrected scanning transmission electron microscopy (STEM). MD simulations were performed using the Vienna
ab initio Simulation Package (VASP). The
generalized gradient approximation (GGA)in the form of the Perdew–Burke–Ernzerhof(PBE) exchange-correlation functional was implemented. A kinetic energy cutoﬀ of 450 eV
was used for the expansion of the charge density, and all calculations were performed usinga Brillouin zone sampling with 2 × 2 × 1 kmesh. The cell sizes are 1.070 nm (∼ 4 unit
cells of SiC) along the [11̄00] axis and 0.926
nm (∼ 3 unit cells of 4H-SiC) along the [112̄0]axis where 6 Si-C bilayers along the [0001] axis
were involved as the 4H-SiC substrates. SiO2
structures were constructed by randomly distributed 30 Si and 60 O atoms on each faces
as the initial models. The opposite side of the
4H-SiC {0001} slabs was passivated by hydrogen atoms to eliminate the eﬀects of dangling
bonds, and a vacuum slab with a thickness of
1.5 nm along the [0001] axis was embedded
into the outer side of SiO2 amorphous layer toavoid unwanted interactions. The supercell includes 246 atoms in total, including 102 silicon,
72 carbon, 60 oxygen and 12 hydrogen atoms.
The initial model was heated at 3500 K for 15.0
ps, cooled at −355 K/ps, and then relaxed at
300 K for 9.0 ps. Finally, the atomic positions

were fully relaxed until the maximal force oneach atom decreased to below 0.05 eV/Åandthe total energy converged to less than 1.0 ×
10−4 eV/atom.
We observed the interfacial O atoms on theamorphous SiO2 layer on the Si face show clear
atomic ordering with a rigid O-Si bridge structure across the SiO2 /4H-SiC interface, whichindicates a slow oxidation rate. The C faceinterface, on the other hand, was atomicallyrough and formed many dangling bonds withhigh oxidation rate. We concluded that theformation of a stable and ﬂat oxidation front
by O atom ordering and the suppression of interface defects or residual C will be useful for
designing high-performance 4H-SiC MOSFET
devices [1].
Polycrystalline ZnO is widely used andits properties are largely inﬂuenced by grainboundaries (GBs). We systematically investigated [0001]-symmetrical tilt GBs in polycrystalline wurtzite ZnO by ﬁrst-principles
calculations and atomic-resolution aberrationcorrected STEM. The Zn sites form an octahedron capped by Zn tetrahedra, with O atoms
located at the center, forming a Zn octahedron
which can be regarded as a hexagon viewed
along the c = [0001]-axis. The misorientationof ZnO [0001]-symmetrical tilt GBs is deﬁned
by the dihedral angle 2θ between the (112̄0)planes in adjacent grains, where θ is the dihedral angle between the (112̄0) and (p q p + q 0)GB planes.
The stable structures of the Σ7 (1230) at

113

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

21.79◦ , Σ13 (257̄0) at 27.80◦ and Σ13 (1340)
at 32.20◦ CSL GBs were independently simulated. The upper crystal was shifted toward
the lower crystal by 0.01 nm along the x-, y, and z-directions. Static lattice calculations
were carried out with General Utility Lattice
Program (GULP) using the Buckingham-form
two-body interatomic potential for ZnO. The
local minima were then optimized in the ﬁrstprinciples calculations with the plane-wave basis projector-augmented wave method included
in VASP. The GGA was used for the exchangecorrelation potentials in the PBE form and thecutoﬀ energy was 400 eV for all plane-wave basis sets. All ionic positions as well as the supercell volume were relaxed such that all forceswere < 1.0×10−5 eV/Å, and the residual force
on each relaxed atom was < 0.05 eV/Å under the constant pressure of 0 Pa and constanttemperature of 0 K. The Brillouin-zone integrations were performed over a 4 × 2 × 4 kpoint mesh generated by the Monkhorst-Packscheme.
The GB and bulk structural units foundin this study are similar to capped trigonal
prisms (CTPs) and octahedra, corresponding to octagons and hexagons viewed along
the c-axis, respectively.
After relaxation,
the stable atomic structures were obtained
for the Σ7 (123̄0) GB at 21.79◦ concave octagons are aligned periodically and are ﬁlled
in by two hexagons. On the other hand,the stable atomic structure for the Σ13 (134̄0)
GB at 32.20◦ was alternative arrangement ofa hexagon and a concave octagon in a period. The stable atomic structure calculated
for the Σ13 (257̄0) GB is composed of four
hexagons alternately arranged with three octagons, which is reportedly a one-to-one combination of the Σ7 (123̄0) GB at 21.79◦ and theΣ13 (134̄0) GB at 32.20◦ . These hexagons andconcave octagons are fundamental structuralunits to describe [0001]-symmetrical tilt GBs
in wurtzite ZnO, and we found that the arrangement of structural units changes aroundthe critical angle of 32.20◦ , where the (134̄0)
GB at 32.20◦ is considered a singular GB. The
zigzag arrangements are achieved because of
the octagons inclined at ±30◦ to the GB plane
above 32.20◦ , exhibiting the same core structure as the Σ7 (1230) GB at 21.79◦ [2].

References[1] M. Saito, H. Li, K. Inoue, H. Matsuhata,
Y. Ikuhara: Acta Mater. 221, 117360(2021).
[2] K. Inoue, J.Y. Roh, K. Kawahara, M.
Saito, M. Kotani, Y. Ikuhara: Acta
Mater. 212, 116864 (2021).

114

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Structure and Property of Basalt Melt and Glass
Toshiaki IITAKA
RIKEN Center for Computational Science,
2-1 Hirosawa, Wako, Saitama, 351-0198 JAPANKnowledge of the macroscopic physicaltemperatures of 2200 K (0 to 82 GPa) and 3000

properties of silicate melt (magma) under highK (40–70 GPa). A key finding is that, althoughtemperature and high pressure is extremely

the density and coordination numbers around Si

important for understanding the phenomena ofand Al increase with pressure, the Si–O and Al–

terrestrial planets such as the formation of the

O bonds become more ionic and weaker. The

early Earth, the transport of volatile elements,

temporal atomic interactions at high pressureand volcanic eruptions.

are fluxional and fragile, making the atoms

We study the influence of the microscopic

more mobile and reversing the trend in

structures (structural factors, radial distribution

transport properties at pressures near 50 GPa.

function,

network

The reversed melt viscosity under lower mantle

topology, etc.) of basalt melts and basalt glasses,

conditions allows new constraints on thewhich are related to the history of earth

timescales of the early Earth’s magma oceansformation, on their macroscopic physicalandproperties

viscosity,

explanation for the horizontal deflections ofthermal conductivity, etc.), paying particularsuperplumes at ~1000 km below the Earth’s

attention to the medium-range structure, trace

surface [1, 2].

coordination(density,

number,diffusivity,

also

provides

the

first

tantalizing

References

elements, and entropy effect (temperature[1]

effect).

A. Majumdar et al., Nature

Transport properties like diffusivity andCommunications 11, 4815 (2020).

viscosity of melts dictated the evolution of the

http://dx.doi.org/10.1038/s41467-

Earth’s early magma oceans. We reported the020-18660-w

structure,[2]

density,diffusivity,

electrical

S. Q. Feng et al., Phys.

conductivity and viscosity of a model basaltic

Chem. Miner. 48, 41 (2021).

(Ca 11 Mg 7 Al 8 Si 22 O 74 )http://dx.doi.org/10.1007/s00269-

melt

from

first-

021-01165-3

principles molecular dynamics calculations at

115

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First principles calculations of two-dimensional ordering
process at organic/inorganic materials interface
Tomoya Shiota,1 Yuuichi Orimoto,2 Yuriko Aoki2
1

Department of Molecular and Material Sciences, Interdisciplinary Graduate School ofEngineering Sciences, Kyushu University, Fukuoka 816-8580, Japan.
2

Department of Material Sciences, Faculty of Engineering Sciences,
Kyushu University, Fukuoka 816-8580, Japan

We have studied dynamic processes at the

500K. The kinetic Monte Carlo simulation atinterface between the organic and inorganic

500K using the energy barriers obtained by thesemiconductor materials by means of the first-

DFT calculations generally reproduced theprinciples calculations and the kinetic Monteexperimental observation.

Carlo simulations. All the first-principles

Toward a more accurate simulation, we

calculations were conducted within the density

investigated the flip-flop motions of surfacefunctional theory (DFT) via the Vienna ab initio

buckled dimers both for the clean Si andsimulation package (VASP), version 5.4.4 in theGe(001) as elementary processes because their

ISSP's massively parallel supercomputer system.

flip-flop motions are well known for their

In this project, a system of 3,4,9,10 - perylene

dynamic properties on both surfaces. It was(PTCDA)

found that the flip-flop motion is not a simple

molecules adsorbed on a Ge(001) surface wasseesaw motion but a two-step motion. In

employed as a prototype of the organic-

addition, by considering the flip-flop motion, theinorganic semiconductor interface. Previously,

diffusion barriers of a PTCDA molecule on

we found that a solely adsorbed PTCDA

Ge(001) were reduced by about 20%.

tetracarboxylic

dianhydride

molecule is not mobile at 500K but is mobile at
700K. Furthermore, it was found that the
PTCDA molecules diffuse concertedly with the
dynamics of Ge(001) surface, involving the flipflop motion of the Ge buckled dimers. [1]First, we explored the diffusion process ofinteractions. As a result, it was found that theFig. 1. (a) Adsorption of PTCDA on Ge(001)surface. (b) The flip-flop motion of the Si andGe(001) surfaces.

diffusion mechanism of PTCDA molecules

References

under the intermolecular interactions is similar[1] T. Shiota, W. Mizukami, H. Tochihara, K.

to that of the lone adsorbed molecule, and the

Yagyu, T. Suzuki, and Y. Aoki, J. Phys. Chem. C

PTCDA molecules were found to be mobile

124, 24704 (2020)

PTCDA molecules under the intermoleculareven under the intermolecular interactions at

116

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of TokyoRSDFT calculation of atomic displacement captured byenergy dissipation channel of noncontact atomic force
microscope
Toyoko ARAI1, Takahiro Kato1, Jun-Ichi Iwata2, Atsushi Oshiyama3
1

Graduate school of Natural Science and Technology, Kanazawa University, Kanazawa,
Ishikawa, 920-1192
2

Institute of Innovative Research, Tokyo Institute of Technology, Tokyo 152-8550
3

IMaSS, Nagoya University, Nagoya, Aichi 464-8601

In the last two decades the conservative andbackbonds of the Si adatom and subsequent

non-conservative (dissipative) force interactions

bond formation in an alternative atomic

between a tip and a sample in proximity haveconfiguration of the Si surface atoms, including

been intensively examined by non-contact

the Si atom at the apex of AFM tip as shown inatomicIn

Fig. 1. In this study, to explicate the mechanism

particular, the channel of nc-AFM to measureof the phenomena, the forces between Si tip andthe energy dissipation through the non-Si(111)-(7×7) surface were calculated using

conservative

been

real-space density functional theory (RSDFT)

regarded to have great potential to explore nano-[2]. RSDFT is a first-principles program

mechanical phenomena; the dissipation channel

developed by Oshiyama and Iwata, which uses agives the change in the amount of mechanicalreal-spaceenergy to maintain the constant oscillation

pseudopotential method.

force

microscopy

force(nc-AFM).

interactions

hasdifference

methodandaamplitude of an AFM cantilever. We had

First, we optimized the structure of a Si(111)-

experimentally found that, when the nc-AFM(7×7) reconstructed surface consisting of 298 Si

image showed the high resolution for theatoms, the back surface of which was terminatedindividual adatoms on a Si(111)-(7×7) surface,with 49 H atoms, and the structure of a [001]-

the simultaneously obtained energy dissipation

oriented Si tip by RSDFT; the tip consisted of 30

signal increased in close proximity over the

Si atoms having a Si dimer at the tip apex, thehollow sites surrounded by a Si adatom and a Si

back side of which was terminated with 18 H

rest atom with a dangling bond per each [1].

atoms, named “dimer tip” [3]. Next, we

Based on the experimental results, the dynamic

calculated the structure changes and the forces

atomic processes responsible for the dissipation

acting between the tip and the Si surface during

had been discussed in regard to breaking of the

tip approaching and retracting over the hollow

117

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Fig. 1 Model of backbond breaking andformation during a tip approach and retraction
cycle with a hysteresis loop in the forcedistance curve, which is estimated fromexperimental results.

distance curves with hysteresis loop weresite at the separations between the tip apex atom

obtained by RSDFT calculations. The greyed

and the Si adatom on the surface from 5.0 to 2.0

area

Å, as shown in Fig. 2. On the approaching at 3.5

corresponding

Å, one of the three backbonds of the adatom wasdissipation energy, was about 1.6 eV. This value

broken, and the adatom was pulled toward thewas larger than the value of about 0.4 eV thattip. At 2.2 Å, the force between the tip and thewas experimentally measured in ref. １. Thissurface turned to repulsive. At 2.05 Å, the tipmeans that the initial structure of the tip should

apex atom bonded with the rest atom under abe reconsidered to clarify the mechanism for thestrong attractive force. Starting from thedissipation energy of the tip–sample system.

Fig. 2 Force–distance (separation) curve withthe optimized structures of the Si tip andSi(111)-(7×7) surface, calculated by RSDFT.

surroundedthetocurves
the

in

Fig.

2,

nonconservative

structure of the tip and the Si surface, we
calculate the optimized structure by retracting

References

the tip by 0.1 Å step, which was the separation[1] T. Arai, R. Inamura, D. Kura, and M.

between the H-terminated Si planes of the tip

Tomitori, Phys. Rev. B 97, 115428 (2018).

and the Si surface. Up to the separation of 2.8 Å,[2] J.-I. Iwata, D. Takahashi, A. Oshiyama, T.

the attractive force became stronger, and at 2.9

Boku, K. Shiraishi, S. Okada, and K. Yabana, J.

Å the structure and the force almost reverted toComput. Phys. 229, 2339 (2010).

those at 2.9 Å on the approaching. When being[3] P. Pou, S. A. Ghasemi, P. Jelinek, T.

further retracted, the structure and the force were

Lenosky, S. Goedecker and R. Perezet,

the same as those on approaching. The force-

Nanotechnology 20, 264015 (2009).

118

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Elucidation of Photocatalytic Reaction Mechanism by Ab
Initio Calculations
Tatsuya JOUTSUKA
Institute of Quantum Beam Science, Graduate School of Science and Engineering, Ibaraki
University, Hitachi, 316-8511 Ibaraki, Japan

Frontier Research Center for Applied Atomic Sciences, Ibaraki University, 162-1
Shirakata, Tokai, Ibaraki 319-1106, Japan.
We have studied the reaction mechanism ofanalysis, [4] and compute the band structure in

photocatalytic oxidation in TiO2 [1] using cp2k
program packagetitania and the other photocatalysts.

(MPI parallelization withmultiple threads). We focused on the energetics
and dynamics during the polaron migration.
This is a collaborative work with the Prof.
Satoshi Yamauchi at Ibaraki University, who
conducted the LPCVD experiment.
This experimental and theoretical study
clarifies

the

facet

dependence

of

the

photocatalytic reaction at the anatase TiO2
surface. The <112>-oriented anatase TiO2 layer
is deposited on the Ru(0001) substrate at 360°C

Fig. 1: Molecular mechanism of hole migration.

by low-pressure chemical vapor deposition.
The deposition rate is three times higher thanReferencesthat of the multi-orientation layer on the Pylex[1] T. Joutsuka, H. Yoshinari, and S. Yamauchi:

glass. The density functional theory (DFT)

Bull. Chem. Soc. Jpn., 94 (2021) 106–111.

calculations using the constrained DFT method(Selected Paper, Inside Cover, Open Access).

[2-3] and the hybrid functionals show that the[2] T. Joutsuka and K. Ando: J. Phys. Chem. B,(112) surface stabilizes the adsorbed water

124 (2020) 8323–8330.

molecule most strongly.

[3] T. Joutsuka and K. Ando: Chem. Lett., 50

We are extending this study to analyze thedynamical

properties

by

the

(2021) 1325–1328.

constrained[4] T. Joutsuka and K. Ando: J. Phys. Chem. A,

DFT[1-3] and the dynamical reaction path

115 (2011) 678–684.

119

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Theoretical Analysis of Metal Film Growth Mechanism in
Chemical Vapor Deposition
Tatsuya JOUTSUKA
Institute of Quantum Beam Science, Graduate School of Science and Engineering, Ibaraki
University, Hitachi, 316-8511 Ibaraki, Japan

Frontier Research Center for Applied Atomic Sciences, Ibaraki University, 162-1
Shirakata, Tokai, Ibaraki 319-1106, Japan.
We have studied low-pressure chemical

catalytic performance by chemically modifying

vapor deposition of Cu on Ru substrate using

the catalytic reaction field.

CuI using cp2k and VASP program packages(MPI parallelization with multiple threads). In
this year, we also studied the mechanism of Cu
nanoparticles by reducing spinel type oxides.
This

is

a collaborative work withtheexperimentalists [1-3] (Prof. Satoshi Yamauchi
and Assistant Prof. Shohei Tada at Ibaraki

Fig.

University), who conducted the experiment.

nanoparticles by H2 reduction of inverse spinel

In this study, we elucidated the molecular1:

FormationmechanismofCutype Cu-doped MgAl2O4.

mechanism of Cu deposition using copper(I)iodide (CuI) on Ru substrate by low-pressure

Referenceschemical vapor deposition through density[1] T. Joutsuka, S. Yamauchi: Chem. Phys. Lett.

functional theory (DFT) molecular dynamics

741 (2020) 137108.

(MD) simulations by cp2k. In addition, the[2] T. Nishikawa, K. Horiuchi, T. Joutsuka andcharge and chemical environment around Cu

S. Yamauchi: J. Cryst. Growth, 549 (2020)

atoms in spinel type oxide Cu-doped MgAl2O4.

125849.

Figure 1 shows the proposed formation[3] T. Joutsuka, S. Tada, et al.: submitted.

mechanism of Cu nanoparticles after reducing[4] T. Joutsuka and K. Ando: Chem. Lett., 50

inverse spinel-type Cu-doped MgAl2O4.

(2021) 1325–1328.

Based on this study, we are developing the[5]computational method to compute the surface

T.

Joutsuka:submitted

10.26434/chemrxiv-2022-9t6vf).

acidity on solid oxides, [4-5] and improving the

120

(DOI:

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Prediction of properties of organic ferroelectrics andpiezoelectrics by ﬁrst-principles calculation
Shoji ISHIBASHI
Research Center for Computational Design of Advanced Functional Materials (CD-FMat),
National Institute of Advanced Industrial Science and Technology (AIST)
Tsukuba, Ibaraki 305-8568

Among hydrogen-bonded systems consistingof small organic molecules, some shows ferroelectricity or antiferroelectricity [1]. They
contain neither toxic nor rare elements andare expected as environmentally friendly materials.
Spontaneous polarization, which
is one of the most important properties
for the ferroelectrics, can be predictedby electronic-structure calculations together
with the Berry phase theory.
We calculated polarization values for seven representative hydrogen-bonded ferroelectrics (croconic acid (CRCA), 2-phenylmalondialdehyde(PhMDA), 3-hydroxyphenalenone (HPLN),
cyclobutene-1,2-dicarboxylic acid (CBDC),
2-methylbenzimidazole (MBI), 5,6-dichloro2-methylbenzimidazole (DC-MBI), and 3anilinoacrolein anil (ALAA)) and the obtained
results are in excellent agreement with the experimental results [2]. These calculations were
performed using experimentally obtained crystal structures with only hydrogen positions
computationally optimized. It is well known
that hydrogen positions are diﬃcult to accurately determine by x-ray diﬀraction analysis. Accurate experimental structure data are
not always available. In such cases, computational optimizations for all of the crystal structure parameters are necessary. Nonetheless,the van der Waals interaction, which is often
crucial in organic solids, cannot be accurately
described by the local density approximation(LDA) or the generalized gradient approxima-tion (GGA). In the present study, we apply
van der Waals density-functional theory (vdWDFT) calculations to reproduce crystal structure parameters of the seven compounds mentioned above. For the vdW-DFT method, weused two forms: the van der Waals densityfunctional consistent-exchange (cx) method [3]and the revised Vydrov–van Voorhis (rVV10)
method [4]. Calculations were performed using
the QMAS code.
Figure 1 shows the deviations of the calculated structure parameters (lattice parameters a, b, and c, and hydrogen-bond lengthd, namely, the O· · ·O or N· · ·N distance) (a)from room-temperature experimental resultsand (b) from 0-K extrapolated values, as wellas (c) their averages and bounds [5]. The calculated results obtained by the LDA and thePBE version of the GGA are also shown. The
experimental structures are at room temperature, whereas the calculation results correspond to 0 K. Unfortunately, only a limited
number of low-temperature structures havebeen reported and the 0-K values were extrapolated using them to prepare Fig. 1 (b). The
vdW-DFT results show good agreement withexperimental results, implying that an important step for the computational materials design of organic ferroelectrics has been achieved.
Spontaneous polarization values obtained by
the vdW-DFT calculations are again in good
agreement with the experimental results. The
success in reproducing the structure parame-121

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo
25(c) Overall(b) 0 K Extrapolated(a) Room Temperature

cx

20

rVV10
LDA
PBE

Deviation (%)

15

10

5

0

-5

-10

abcd

abcd

CRCA

abcd

abcd

HPLN
PhMDA

abcd

abcd

MBI
CBDC

abcd

abcd

ALAA

CRCA

DC-MBI

abcd

abcd

abcd

HPLN
PhMDA

abcd

abcd

MBI
CBDC

abcd

RT

0 K

ALAA
DC-MBI

Figure 1: Deviations of calculated structure parameters (a) from room-temperature experimental
results and (b) from 0-K extrapolated values, as well as (c) their averages and bounds (S.
Ishibashi, S. Horiuchi, and R. Kumai, “Hydrogen-bonded single-component organic ferroelectrics
revisited by van der Waals density-functional theory calculations”, Phys. Rev. Mater. 5, 094409(2021), DOI: 10.1103/PhysRevMaterials.5.094409).
ters also enables the simulation of electromechanical responses. Direct and converse piezoelectric coeﬃcients are evaluated for CRCA,
PhMDA and DC-MBI and the resultant values show acceptable agreement with the experimental values if possible objective factors
are considered [5].

Organic antiferroelectrics also are promising. For example, their applications as highpower electrical energy storage materials are
expected utilizing electric-ﬁeld-induced phasechanges. Usually, they consist of two or morepolar subunits. For squaric acid (SQA), we
have demonstrated that the polarization of theelectric-ﬁeld-induced phase can be estimatedby proper summation of the subunit polarization vectors [6]. By applying this procedure,
we have successfully predicted polarization values for the electric-ﬁeld-induce phases of several organic antiferroelectrics and the obtained
results are in excellent agreement with experimental results [7].

References[1] S. Horiuchi and S. Ishibashi, J. Phys. Soc.
Jpn. 89, 051009 (2020).
[2] S. Horiuchi, K. Kobayashi, R. Kumai,
and S. Ishibashi, Nat. Commun. 8, 14426(2017).
[3] K. Berland and P. Hyldgaard, Phys. Rev.
B 89, 035412 (2014).
[4] O. A. Vydrov and T. Van Voorhis, J.
Chem. Phys. 133, 244103 (2010); R.
Sabatini, T. Gorni, and S. de Gironcoli,
Phys. Rev. B 87, 041108(R) (2013).
[5] S. Ishibashi, S. Horiuchi, and R. Kumai,
Phys. Rev. Mater. 5, 094409 (2021).
[6] S. Horiuchi, R. Kumai, and S. Ishibashi,
Chem. Sci. 9, 425 (2018).
[7] S. Horiuchi and S. Ishibashi, Chem. Sci.
12, 14198 (2021).

122

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Global search for the stable structures of two-dimensional
materials based on Gaussian process regression
Yuji HAMAMOTO
Department of Precision Engineering, Graduate School of Engineering,Osaka University, Yamada-oka, Suita, Osaka 565-0871
Silicene, a silicon analog of graphene has

algorithms. As an example, we here apply

drawn a great deal of attention not only as a

GOFEE to the global search for the stablepromising candidate for topological insulatorsstructures of silicene supported on the Ag(111)from a scientific point of view, but also as ansurface. Our results reveal that GOFEE not only

atomically thin semiconductor material from an

reproduces successfully the stable structures

application point of view. On the Ag(111)

identified so far, but also predicts the existence

surface, however, the sp3-like buckling nature ofof a variety of metastable structures which are

its honeycomb lattice leads to the emergence of

too subtle to determine experimentally.

a varieties of phases characterized by different
patterns of buckling structures as compared withthe freestanding structure predicted theoretically[1], hindering the elucidation of the whole
picture of the electronic and structural propertiesof silicene supported on solid surfaces.
First-principles calculation based on thedensity functional theory (DFT) is a powerful
technique to predict stable structures of newmaterials, but its performance depends on the

Figure: Obtained structures corresponding toinitial guesses, and the global search for stablethe 4×4-a and b (top left and right) and thestructures is computationally demanding in most

Ö13´Ö13-a and b (bottom left and right) phases

cases. To tackle the difficulty in the globalof silicene on the Ag(111) surface.

structure search by the DFT calculation, we

References

adopt an efficient machine-learning technique[1] K. Takeda and K. Shiraishi, Phys. Rev. B

named global optimization with first-principles

50, 14916 (1994).

energy expression (GOFEE) [2], which is based[2] M. K. Bisbo and B. Hammer, Phys. Rev.

on Gaussian process regression and evolutionary

Lett. 124, 086102 (2020).

123

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles study on the origin of phase stabilityof Mg–Zn–Y alloys with long-period stacking order
Takao TSUMURAYA
POIE, Kumamoto University, Kumamoto 860-8555, Japan
Magnesium Research Center, Kumamoto University, Kumamoto 860-8555, JapanThe design of alloys based on magnesium(Mg) has attracted much attention since Mg
is the lightest substance element for structural metals. A dilute Mg alloy containing
one at.% of Zn and two at.% of Y showed
high-strength properties with a tensile yieldstrength of ∼600 MPa [1]. This strength iscoupled with the appearance of a unique crystal structure called long-period stacking order(LPSO). In the LPSO-based Mg alloy, soluteatoms of Zn and Y are enriched in a few layers of the (0001) plane of the hexagonal closepacked (HCP) Mg matrix [2]. Shockley partial
dislocation (stacking faults) occurs in the solute atoms enriched layer, which orders periodically along the [001] direction of the hexagonallattice. Furthermore, structural analysis using scanning transmission electron microscopy(STEM) revealed that L12 -type clusters of solute elements (Zn and Y) are embedded in the
Mg matrix near the stacking fault [3].
In this study, the electronic origin of thephase stability of Mg–Zn–Y based alloys with
LPSO is studied based on ﬁrst-principles
density-functional theory (DFT) calculations.
In particular, we aim to understand the role of
the Zn atom in 18R–Mg58 Zn7 Y8 alloy, and thepossible realization of Zn vacancy formationis studied using a convex hull energy diagram.
The heat of formation as a function of the number of Zn vacancies is calculated, and we show
how the electronic structure near the Fermi
level aﬀects structural stability. In the present
DFT calculations used a pseudopotential tech-nique with plane-wave basis sets adopting
the projected augmented-wave method implemented in Quantum ESPRESSO. The
exchange-correlation functional is the generalized gradient approximation proposed by
Perdew, Burke, and Ernzerhof.
For calculating the heat of formation, thenumber of Zn atoms in Mg58 Zn7 Y8 is deduced
using the ﬁxed number of Mg and Y atoms
to understand whether the Mg-Zn-Y alloy includes Zn vacancies that can be formed fromthe elemental metals of Mg, Zn, and Y. The
obtained scatter plot of the heat of formation
as a function of the Zn composition ratio indicates that all phases, including Zn vacancies,
are thermodynamically unstable and decomposed into Mg58 Zn7 Y8 (x = 0) and elemental metals of 58Mg+8Y phases. Nevertheless,the hull distance decreases as vacancies (x) decrease. Therefore, we consider the state inwhich Zn7 Y8 clusters with an L12 structure are
already formed and plot the convex hull diagram. Zn vacancies with x = 1 and 3 are on theconvex hull. These states should be observed
below the decomposition temperature. Therefore, we conclude that the Zn atoms in the
LPSO alloy are stable even if they are about
half of the number of Y atoms.
In general, the density of states (DOS) for
HCP metal is characterized by a deep valley
or dip near the Fermi level (EF ), where EF isat the lowest position of the valley. The depthof the valleys in the DOS represents the magnitude of splitting the prominent peaks of the124Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

DOS between those in occupied and unoccupied states. The width of the valley or dip
appears as a diﬀerence in the enthalpy of formation (or cohesive energy) between diﬀerent
metals [4, 5]. By focusing on the partial DOS
of Zn and its neighboring Mg, we ﬁnd that theformation of the bonding state between Zn andthe surrounding Mg atoms has a crucial role instabilizing the LPSO structure.

References[1] Y. Kawamura, K. Hayashi, A. Inoue, andT. Masumoto, Mater. Trans. 42 1172
(2001).
[2] E. Abe, Y. Kawamura, K. Hayashi, and
A. Inoue, Acta Metall. 50 3845 (2002).
[3] D. Egusa and E. Abe: Acta Mater. 60 166(2012).
[4] S. Inoue and J. Yamashita: J. Phys. Soc.
Jpn 35 677 (1973).
[5] T. Tsumuraya, I. Watanabe, and T.
Sawaguchi: Phys. Rev. Research 3 033215(2021).

125

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Metastability of 2D and 3D binary ordered alloys
Shota ONO
Department of Electrical, Electronic and Computer Engineering, Gifu University
Gifu 501-1193

It is well-known that the fcc (hcp) metals
in the hcp (fcc) structure are dynamically stable, while the fcc (bcc) metals in the bcc (fcc)
structure are unstable [1]. To understand themetastability relationship between diﬀerent
crystal structures, we calculated phonon dispersions for 2D and 3D systems using densityfunctional perturbation theory implemented in
Quantum ESPRESSO [2]. The calculations
were done using system B and C.
In our previous work, we showed that if theelementals metals in the 2D hexagonal structure is dynamically stable, then the fcc and hcp
strcutures are also dynamically stable, showing that 2D structures serve as building blocks
for constructing the 3D structures [3]. In thepresent work, we extend this investigation tobinary ordered alloys of CuX with X being 46
metallic elements [4, 5]. We showed that if the
CuX in the AB layer structure is dynamicallystable, then that in the AB and ABC stacking structures is also dynamically stable (seeFig. 1 for the crystal structures). In addition,
we calculated phonon dispersions for the CuX,
AgX, and AuX in the AB layer structure, andfound that 48 out of 135 ordered alloys are dynamically stable, predicting that those alloys
in the AB and ABC stacking structures arealso dynamically stable [5]. The present work
has shed light on the 2D structures as a building block for computational material design.
We also predicted potential substrates thatallow immiscible metals (Pb and Sn) to form
the surface alloys by using density-functional
theory [6]. The prediction of their surface geometry is left for future work.

!"#.&/01#,23'(.04#56*0/'672-!"#$%&'()*+#,"!-!"8#$%&'()*+#,9:"-Figure 1: Schematic illustration of 2D and 3D
strcutures.

References[1] G. Grimvall et al., Rev. Mod. Phys. 84(2012) 945.
[2] P. Giannozzi et al., J. Phys.: Condens.
Matter 29 (2017) 465901.
[3] S. Ono: Phys. Rev. B 102 (2020) 165424.
[4] S. Ono: Sci. Rep. 11 (2021) 14588.
[5] S. Ono: Phys. Rev. Materials 5 (2021)
104004.
[6] S. Ono, J. Yuhara, and J. Onoe: Chem.
Phys. Lett. 776 (2021) 138696.

126

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Ab Initio Molecular Dynamics Study on Static Structure of
Glass Materials
Akihide KOURA
Technical division,
Kumamoto University, Kurokami, Chuo-ku, Kumamot,o 860-8555

We have studied a condensed matter system

based on the AIMD simulations of liquid andwhich are crystal, liquids, and amorphous. In

amorphous V2O5. Figure 1 shows the root mean

this project, we focused on the static structure

square errors (RMSE) and partial radialof amorphous V2O5 state. The ab initio

distribution functions 𝑔𝑔𝛼𝛼𝛼𝛼 (𝑟𝑟) (𝛼𝛼, 𝛽𝛽 ∈ V, O) of

molecular dynamics (AIMD) simulations wereliquid V2O5. According to our previous study

performed with our original QXMD code [1].

[3], considering not only energy but also atomic

In order to investigate the static structure offorce and virial tensor improves ANN potential

amorphous V2O5, we employed the machine

better. Based on the learning condition, our MD

learning method using artificial neural network

simulation using the ANN potential basically(ANN) to construct the interatomic potential [2].

represent the FPMD simulations except for the

The Aenet package [2] was used for training

unreasonable peak of 𝑔𝑔OO (𝑟𝑟) at approximately

𝑟𝑟 = 0.4 Å. In the future work, we develop thetraining method by other technics, e.g., active

learning from the result of MD simulations
with ANN potential for liquid systems.

References[1] F. Shimojo, S. Fukushima, H. Kumazoe, M.
Misawa, S. Ohmura, P. Rajak, K. Shimamura, L.
Bassman, S. Tiwari, R. K. Kalia, A. Nakano,and P. Vashishta: SoftwareX 10 (2019).
[2] N. Artrith and A. Urban: Comput. Mater.
Sci. 114 (2016) 135.
[3] A. irie, S. Fukushima, A. Koura, K.
Fig. 1: (Left) Root mean square errors and(right) partial radial distribution functions.

Shimamura, and F. Shimojo, J. Phys. Soc. Jpn.
90, 094603 (2021).

127

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Atomistic simulations for excited materials by strong laser
Yasushi Shinohara1,2 , Mizuki Tani3
1

Photon Science Center,

The University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656.
2

NTT Basic Research Laboratories,

NTT Corporation, 3-1 Morinosato Wakamiya, Atsugi, Kanagawa 243-0198.
3

Department of Nuclear Engineering and Management

The University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656.
We have performed investigations thatelectron occupation shows symptom of multiple

describe nonequilibrium electron-ion dynamics

Rabi-flopping to higher energy bands for bulk Al.

realized by strong laser fields irradiation.

We

have

developedasemiclassicalComplex pathways accompany these dynamics

theoretical framework for electron-ion coupledalong different phases depending on time anddynamics under a strong laser field [3]. Free-

positions.

andelectronic behavior of metal is accurately

nonequilibrium characteristics prohibit direct

reproduced with more than 50-times lightermeasurements for the intermediate processes ofcomputational

the dynamics. Our investigation is aimed tolightweight cost allows us to simulate spatially

understand what processes could happen in thelarger systems, including non-uniform structures,

laser-excited materials, including insulators andsuch as surfaces or impurity, and electron-ion

metals, via theoretical simulations.

coupled dynamics for the laser-excited material

These

multi-physicalWe have accumulated simulated electronic

cost

than

TDDFT.

The

response.

structure data just after ultrashort laser pulse byusing SALMON code [1] based on time-

References

dependent density-functional theory [2]. This

[1] Erich Runge and E. K. U. Gross, Phys. Rev.

data is the key quantity of the laser-excited

Lett. 52 (1984) 997.

dynamics' initial process, which triggers whole[2] M. Noda, et al., Comp. Phys. Comm. 235

subsequent processes. We reported Rabi-type(2018) 356.

absorption oscillation in bulk Al with resonant[3] Mizuki Tani, Tomohito Otobe, Yasushi

excitation by increasing the intensity of the

Shinohara, and Kenichi L. Ishikawa, Phys. Rev.

electric field [3]. By analyzing the occupation

B 104 (2021) 075157.

evolution as a function of the intensity, the128

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Study on structural elementary excitations at
semiconductor surfaces and interfaces
Hiroyuki KAGESHIMA
Graduate School of Natural Science and Technology, Shimane University
1060 Nishi-Kawatsucho, Matsue, Shimane 690-8504

In this project, we have been focused on
physical properties of structural elementary
excitations of semiconductor surfaces and interfaces. In this year, we have focused on
the physical properties of defects in SiO2 .
The calculations were performed based on theﬁrst-principles calculation. Program package
PHASE/0 was employed. [1]Defects in SiO2 have been widely studied
because they are the source of device performance degradation for MOS ﬁeld-eﬀect transistors (MOSFETs), which is widely used forthe integrated circuits. The gate insulator of
MOSFET is made of silicon oxide SiO2 . Thosestudies revealed that O vacancies, one of typical defects, can trap charges, assist gate leakcurrent, and also lead to the formation of eternal gate leakage paths. They also revealed
the role of post oxidation annealing. However,recent progress of integrated circuit technology have forced the MOSFET device structurefrom planar to three dimensional. If we consider the vertical body-channel-MOSFET withsilicon oxide gate insulator, the oxide can bestrained because of the device formation process. The interface oxide might be compressed
and the surface oxide be expanded due to thethermal oxidation of silicon for the gate formation process. The compressive and expanded
strain are expected as much as 3% in the worst
case. Such strain should aﬀect on the physicalproperties of defects in silicon oxide gate insulator

and studied strain dependence of its formationenergy and diﬀusion barrier. Based on theseresults, in this year, we have studied stability of various defects including the O vacancy,
the charging eﬀects on them, and their strain
dependences. Some of the vacancies we havestudies are fully terminated by H atom or OH
group, and other are partially terminated. We
modeled the host SiO2 by alpha-quartz crystal with supercell containing 24 Si and 48 O
atoms.
First, we have found a novel O vacancy defect which has only two -OSi bonds. This defect is unique because the widely known O vacancy defect has three -OSi bonds. This new
O vacancy defect is stable because it accompanies a three-Si-member ring within the model,
while ordinary alpha-quartz just contains only
six-Si-member rings. Next, we have found thata dangling bond of Si can surely assist thegate leakage current because they show charging levels in the middle of band gap. Then, we
have found that the compressive strain makes
it easier to release the terminators H or OH,
while it does not signiﬁcally change the charging level positions. From these results, we have
conﬁrmed that strain in the gate insulator silicon oxide surely enhances the device performance degradation. Careful device formation
process must be necessary.

References[1] https://azuma.nims.go.jp

Therefore, we ﬁrst focused on the O vacancy

129

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Ab-initio molecular-dynamics calculations for moleculemodified nanointerfaces
Kenji TSURUTA, Masaaki MISAWA, Narumi MACHIDA,
Kyosuke SATO, and Satoshi OHMURAa
Graduate School of Natural Science and Technology, Okayama University, Okayama 700-8530aDepartment of Civil and Environmental Engineering, Hiroshima Institute of Technology,
Hiroshima 731-5193
Surfaces and/or interfaces of oxide materials
are often modified by other molecules in order
not only to stabilize their structure against
deterioration but also to control surface
morphology. Calcium carbonate (CaCO3),
which is widely used as filler in resins, plastics,
and rubber materials to improve their physical,
mechanical, and thermal properties, is a typicalexample of such phenomena. To control particle
morphology and to improve affinity with basematerials, calcium carbonate is often surfacetreated with organic compounds. This activatedcalcium carbonates (ACCs) not only exhibit
excellent application properties, but also can beapplied as eco-friendly inorganic-organic hybridmaterials. It has been reported [1] that surface
deposits with ~5 nm thickness exist on the (104)surface of calcite covered with stearic acid (SA).
Since this thickness corresponds to about twice
the length of a SA molecule, it is suggested thatSA bilayers may be formed on the calcite (104)surface, but the detailed adsorption state has notbeen clarified to date.
In this study, we performed first-principles
calculations based on density functional theory
using the Vienna Ab initio Simulation Package(VASP). We created a simulation model for three
types of SA-SA and two types of calcite-SAinterfaces as shown schematically in Fig. 1 A-E.
First-principles molecular dynamics (FPMD)
simulations with canonical ensembles at room
temperature were performed for these models,
and the presence or absence of adsorption wasverified by observing changes in the interface
structure in each model. We also verified thestability of the interface structures by analyzing
the energetic properties in the relaxed system.
According to the results of the FPMD
simulations, we considered that the interface
structure in the A, B, and E models are
intrinsically unstable because the interatomic
distance of the bonding part in each modelFig. 1 Schematic illustration of the simulation models
for SA-SA (A-C) and calcite-SA (D, E) interfaces.

gradually increased at room temperature . On theother hand, such behavior has not been observedin the C and D models. Comparing of thepotential energy of the C model with that of two
isolated SA molecule, it was found that theenergy of the C model was 0.4 eV lower. This
result indicates that the C model is energeticallyfavorable and stable interface structure may beformed with an adsorption with C model. Fromthese results, it is considered that despite the firstSA layer is favorable with D model, somemolecules on the first layer may be inverted (or"alternately aligned") in order to adsorb the
second layer with C model in the SA bilayer on
calcite (104) surface[2]. Further computation
with larger system size is necessary to validate
this hypothesis. Development of an effective
interatomic potential based on a Deep-Learning(DL)-based molecular design methods [3] is inprogress for simulations of the bilayer system.
References[1] M. Ricci et al., Langmuir 31, 7563 (2015).
[2] N. Machida, M. Misawa, Y. Kezuka, and K.
Tsuruta (to be published).
[3] K. Sato and K. Tsuruta, Mater. Sci. Forum
1016, 1492 (2021).

130

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles Calculations of Magnetic Insulators
Under Uniform Electric Fields
Naoya YAMAGUCHI
Nanomaterials Research Institute (NanoMaRi), Kanazawa University
Kakuma-machi, Kanazawa, Ishikawa 920-1192
This year, we tried to focus on electric ﬁeld
eﬀects of magnetic insulators. We have developed a linear-combination-of-pseudo-atomicorbital (LCPAO) scheme of a ﬁnite electricﬁeld method based on the modern theory ofpolarization [1] and implemented it to the
OpenMX code [2], a ﬁrst-principles calculation code based on the density functional theory (DFT). Through the scheme, one can consider electric ﬁeld eﬀects even for periodic systems such as bulk insulators. We improved
our implemented code this year. The socalled egg box eﬀect causes numerical errors
due to the real space grid, especially in the
LCPAO method, and it disturbs evaluation of
the forces, dipole moment and stress tensor.
It can be resolved by the grid cell sampling
to repeat the evaluations for the shifted grid
several times based on the frozen LCPAO density matrix, and we implemented the functionin the previous year. We found that a solution was the “one-shot” grid cell sampling, thatis, the forces under the ﬁne grid were evaluated from the density matrix obtained from therough grid without the grid shifts. We calculated electronic and static dielectric constants
and Born eﬀective charges of III-V and II-VI
semiconductors and group IV insulating materials. The calculated values coincided withthe reported computational values even in the
GGA case that the density gradient is sensitive to the egg box eﬀect, and they were anagreement with the experimental values. In
addition, we also debugged the implemented
code for the non-collinear spin density functionals applicable to materials that spin-orbit
interaction (SOI) is important. We also foundthat the SOI has a non-negligible inﬂuence on
the dielectric constants of cases including theheavy elements. For the magnetic insulator
cases, we conﬁrmed that the debugged code
worked well for the antiferromagnetic NiO fordiﬀerent U parameters based on the DFT+U
scheme through the comparison to the previous
study. For some insulators with noncollinear
spins, the self-consistent-ﬁeld (SCF) calculations were diﬃcult to get converges because
the careful mixing during the SCF calculations
were essential. Our implementation will become a powerful tool to expand subjects of researches on materials in investigating the electric ﬁeld eﬀects.
We also developed the new evaluation
scheme for the anomalous Hall and Nernst conductivity without constructing the Wannierfunctions or explicit evaluation of the Kubo
formula [3]. We reproduced the intrinsic contribution of the typical metals such as bcc iron.
Our scheme has advantages of eﬃciency andstability: it is the non-iterative scheme unlike
Wannierization and can avoid numerical instability such as the k-space Dirac monopole. The
scheme is more appropriate for not only large
scale cases but also high-throughput calculations than a Wannier function method.

References[1] N. Yamaguchi and F. Ishii,
https://doi.org/10.48550/arXiv.2203.10441.
[2] http://openmx-square.org.
[3] H. Sawahata, N. Yamaguchi, S. Minami and
F. Ishii,
https://doi.org/10.48550/arXiv.2204.05949.

131

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Optical conductivity of Ca5Ir3O12 by first-principles calculations
Shoya Kawano
Department of Electrical and Electronic Engineering, Kyushu Institute of Technology,
Kitakyushu, Fukuoka, 804-8550

We have studied the title compound

Perdew-Burke-Ernzerhof

type.

Γ-point-

Ca5Ir3O12 which have attracted attention due tocentered 5×5×13 and 10×10×25 k-points

the its strong spin-orbital interaction (SOI) [1].

samplings is set to be compared.. We use theSuch strong SOI is affected to the physicalenergy cutoff of 96 Ry for the wave function, aproperties.

anconvergence criterion for the electronic self-

antiferromagnetic ordering below 7.8 K, aconsistency loop of 10−10 Ry, and the Gaussian

second order phase transition at 105 K and non-smearing method with a smearing width of

linear conductivity along c-axis and so on. We

0.002 Ry. We included the SOI effect. In

performed synchrotron radiation-based infrared

RESPACK,(SR-IR)

its

calculated within random phase approximation.

electronic property for single crystal of

For this calculation, we set cutoff energy for

Ca5Ir3O12.

wepolarization function of 6 Ry. Smearing value

observed reflectivity and optical conductivity.

used in tetrahedron calculation is set to 0.01 eV.

Then, we investigated reflectivity and optical

Lattice parameters of Ca5Ir3O12 is referred to

conductivity of Ca5Ir3O12 by first-principles

previous report [5].

Thiscompoundspectroscopy

From

this

to

shows

revealthemeasurement,

calculations to discuss about our experimental

polarization

functionisFig 1 shows band structure of Ca5Ir3O12
from 5 × 5 × 13 and 10 × 10 × 25 k-point

results.
We

theperformed

densityfunctional

sampling. E = 0 is fermi level. Both results arecalculations by using Quantum ESPRESSO

consistent with previous report []. As showncode [2]. These results were used to calculate

previous report, this material is metal on thisreflectivity and optical conductivity by usinglevel of density functional calculation. ThisRESPACK code [3]. In Quantum ESPRESSO,

metallic trend is different from experimental

we use the norm-conserving pseudo potentials

observation which shows as semiconducting

generated by the code ONCVPSP (Optimized

electrical properties. This is limitation of thisNorm Conserving Vanderbilt PSeudopotential),level calculation. From this result, we use 5×5

and are obtained them from the PseudoDojo [3].

×13 k-point sampling for polarization function

The exchange-correlation function is used

calculation because this calculation is time132

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

consuming.
Fig. 2 shows the calculated reflectivity of
Ca5Ir3O12 parallel to c-plane and c-axis. Thisshows reflectivity R = 1.0 at  = 0 because of
metallic behavior in our results. The broad peakis shown around 5000 cm-1.
Fig.

3

shows

the

calculated

opticalFig. 1 Band structure of Ca5Ir3O12 fromconductivity of Ca5Ir3O12 parallel to c-planedifferent k-point sampling. E = 0 is fermi
and c-axis. These have two peaks at about 1000
level.

cm-1 and about 5000 cm-1 for parallel to c-planeand c-axis. We considered that the peak of
about 1000 cm-1 is shown due to the metallic
behavior in the calculation. Experimental result
has the peak at about 4000 cm-1. This seems to
be good agreement with our 5000 cm-1 peak.
However, we note that the calculation andexperiment have difference of metal or
semiconductor.

Fig. 2 Calculated reflectivity of Ca5Ir3O12
We need to more discussion by using theseparallel to c-plane and c-axis.
results, and also improvement such as metallic
behaviour in the calculation. Using ISSP
Supercomputer, we achieved the first principles
calculation of optical properties such as thecompound Ca5Ir3O12.

References[1] K. Matsuhira et al. J. Phys. Soc. Jpn. 87,(2018).
[2] https://www.quantum-espresso.org
Fig. 3 Calculated optical conductivity of[3] https://sites.google.com/view/kazuma7k6
Ca5Ir3O12 parallel to c-plane and c-axis.
[4] http://www.pseudo-dojo.org[5] M. Charlebois et al. Phys. Rev. B 104,
075153, (2021)

133

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Critical exponent of metal-insulator transition in
doped semiconductors and spin ordering
Yosuke HARASHIMA
Center for Computational Sciences, University of Tsukuba
Tennodai, Tsukuba, Ibaraki 305-8577
Devision of Material Sciences, Nara Institute of Science and Technology
Takayama-cho, Ikoma, Nara 630-0192

Impurity doping have an important eﬀect
on the conduction properties of semiconductors. Evaluating electrical conductivity at azero temperature, the material becomes an insulator at low impurity concentrations, a metal
at high impurity concentrations, and a metalinsulator transition occurs at some concentration in between. This metal-insulator transition is known to occur in many materials regardless of semiconductor or impurity type andis a general and fundamental property of quantum transport.
Electrons in semiconductors are scatteredby randomly aligned impurities and the phaseshifts of wavefunctions are distributed randomly. When those waves overlap, they cancel
each other out and a localized state appears.
This is Anderson localization, and the metalinsulator transition mentioned above is considered to be the Anderson transition between theAnderson localized phase and the delocalized
phase. Continuous transitions such as the Anderson transition can be classiﬁed according to
the critical exponent. The critical exponent νcharacterizes the divergence of the correlation
length ξ ∝ |n − nc |−ν at the transition point(n and nc are impurity concentrations and thecritical concentration). It is universal and independent of the details of the system. The
critical exponent represents the universality ofthe phase transition, therefore, deriving its exact value is an essential problem.

Two diﬀerent values for this critical exponent have been reported from experiments,ν ≈ 0.5 and ≈ 1.0, but the origin of this discrepancy has not been clariﬁed. This problem
is called the exponent puzzle. To solve thisexponent puzzle, the actual situation in the
doped semiconductor needs to be simulated
from ﬁrst-principles.
In this study, we consider a system of randomly distributed P impurity donor ions in
Si. For the analysis of its electronic structure, we used a program we implemented that
self-consistently solves the Kohn-Sham equation based on density functional theory. Experiments suggest that the spin order changes
near this metal-insulator transition point, andthus, the present calculations take into accountthe spin degrees of freedom.
We generated random impurity conﬁgurations. The impurity concentration is ranged
from ≈ 0.8 × 1018 cm−3 to ≈ 2.1 × 1018 cm−3 ,and the system size is taken by 228.6, 285.8,
342.9, and 400.1 Å. The critical concentration has been predicted by convolutional neural network [1]. We sampled ≈ 180, 000 data inthe vicinity of the predicted concentration, andsolved the electronic state self-consistently.

References[1] Y. Harashima, et al.: J. Phys. Soc. Jpn.
90 (2021) 094001.

134

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Exploration of exotic surface sites for catalyst informatics
Yoyo HINUMA
National Institute of Advanced Industrial Science and Technology (AIST)
1-8-31, Midorigaoka, Ikeda, Osaka 563-8577

Surface point defects of metal oxides, forthat are metastable against macroscopic facet

instance O vacancies, have a dominant effectformation and by comparing cleaved surfaces

on heterogeneous catalysis. The Mars-Van

and surfaces suggested by a genetic algorithm

Krevelen mechanism is one of the most(GA) for promising orientations is proposedfrequently encountered catalytic process. In one

and demonstrated using 34 orientations of β-

example, O vacancies on a metal oxide catalyst

Ga2O3

surface act as reaction sites. The energy

terminations

minimum required to remove O from a surface,

accessible are found for β-Ga2O3 and θ- Al2O3,

which is denoted as the surface O vacancy

respectively, where the highest surface energy

formation energy (EOvac), can be used to

was roughly twice of the lowest. The lowest

rationalize and predict catalytic performance insurface O vacancy formation energy (EOvac) insuch a catalytic process. Calculation of EOvac

an accessible surface is 3.04 and 5.46 eV in therequireswith(101) and (20-1) terminations for β-Ga2O3 andsufficient spacing between O vacancies, henceθ-Al2O3, respectively, where the decrease insome estimation of EOvac from less costly

EOvac from the most stable surface is 1.32 andcalculations, such as slab-and-vacuum model

1.11 eV, respectively. The EOvac in accessible

calculations with minimum cell size and evensurfaces showed a good correlation withbulk calculations, will be effective in screening

descriptors

materials for a given purpose.

environment, suggesting that exploiting surfaceaslab-and-vacuum

model

Computational exploration of previously

andθ-Al2O3.

Seven

considered

of

the

local

and

six

experimentally

coordination

O in an unfavorable environment in an

unknown reactive sites is a powerful strategy

accessible

for emergence of new catalytic reactions.

vacancy-related catalyst performance even in

Exoticmaterials that do not show reactivity on thesurfacescanbetheoreticallyinvestigated, but there are very few, if any,
computationalmodelsofhightermination

most stable surface.

indexorientations that considers reconstruction of thesurface. A workflow to efficiently obtain a setof accessible terminations by removing those

135

would

enhanceOActivity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Analysis of the electronic structure of ferromagnetic shapememory alloy and ferrimagnetic spinel oxide
Masao OBATA1,2
1

Institute of Science and Engineering, Kanazawa University, Kanazawa, Ishikawa 920-1192
2

Graduate School of Natural Science and Technology, Kanazawa University, Kanazawa,
Ishikawa, 920-1192
We have investigated the electronic andvector length in Fermi surface nesting. The

magnetic structure of the magnetic shapeobtained nesting vector from the generalized

memory alloy of Ni2MnGa and ferrimagnetic

susceptibility

spinel oxide of NiCo2O4. In these materials, theexperimentally observed crystal structure at low

LDA/GGA

showedtemperatures, such as 10M or 14M long-periodinconsistency with experimental facts, such as

modulated structures [2]. We also estimated thethe period of martensite modulation due toeffect of Mn's on-site Coulomb interaction on

Fermi surface nesting and the properties of

the tetragonal phase [3].

electrical resistivity, respectively. Since these(B) NiCo2O4

discrepancies may originate from electron

The NiCo2O4 spinel has drawn much attention

correlation effects, we employed the quasi-

for its material application due to its low cost,

particle self-consistent GW (QSGW) approach,high electronic conductivity, redox reaction, andwhich is one of the beyond LDA/GGA methods,

high electrochemical activity. However, its

to investigate the electronic structure of thoseelectronic structure is not well understood. We

materials.

investigated the inverse spinel β structure by

level

The

investigation

first-principles

calculationisconsistent

with

the

package ecalj [1] is used for this study.

QSGW. We found that the eg orbital of Co(A) Ni2MnGa

tetrahedral site appears on the Fermi-level and

The Ni2MnGa ferromagnetic shape memory

splits due to the uniaxial crystal symmetry,alloys, which have magnetic and shape memory

stabilizing 3z2 and destabilizing x2-y2 orbitals.

properties, have been attracted as magnetic

References

actuator materials. We applied the QSGW to the[1] https://github.com/tkotani/ecalj/

high-symmetry (cubic) and[2] M. Obata, T. Kotani, T. Oda, submitted.

low-symmetry(tetragonal) phases of Ni2MnGa. We found that[3] M. Zelený, P. Sedlák, O. Heczko, H. Seiner,

the Fermi surface of the cubic phase expands

P. Veřtát, M. Obata, T. Kotani, T. Oda, L. Straka,

compared to GGA, resulting in the reduction ofMater. Des., 209, 109917 (2021).

136

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles calculations of complex metallic alloy
surfaces
Kazuki NOZAWA
Department of Physics and Astronomy, Kagoshima University
1-21-35, Korimoto, Kagoshima 890-0065

We studied the adsorption of pentacene on
the 2-fold surface of an icosahedral Ag-In-Yb
quasicrystal [1](QC). An STM experiment reported that pentacene molecules adsorbed in aquasiperiodic manner on this surface[2].
The QC surface has approximated by a diskshaped atomic cluster with a radius of 2.2 nm(Fig. 1). Atomic positions of the QC surface
has extracted from the structural data of an
X-ray diﬀraction experiment for the isostructural Cd5.7 Yb quasicrystal[3], where Ag and
In atoms occupy the Cd sites in the Ag-In-Yb
QC. In the calculations, the molecular axis ofpentacene is constrained to be parallel to thesurface, and all atomic positions, including thesubstrate atoms, are ﬁxed. Calculations arecarried out with the VASP code mainly on thesystem-B as bulk jobs using 72 nodes.

Fig. 2(a)). Comparably stable orientational
conﬁgurations are found when the molecularaxis is aligned to two 5-fold axes on the surface.
Obtained orientational dependence of the adsorption energy well explains the result of statistical analysis from an STM experiment[2].

(a) 2f-90°

First, we carefully checked the cluster thickness dependence of the adsorption energy. And
we conﬁrmed that the thin ﬁlm composed onlyof the top layer atoms well describes the interactions with pentacene. Thus, we describe theresults below using the thin ﬁlm as the substrate.
The adsorption energy of pentacene shows
strong orientational dependence on the molecular axis of pentacene. Figures 2(a) and (b)show the adsorption energies of pentacene withits molecular axis aligned to two diﬀerent 2fold axes, 2f-90◦ and 2f-180◦ , which are parallel to the surface and perpendicular to eachother. It is clearly shown that pentacene favors one of them (2f-90◦ direction shown in(b) 2f-180°

Figure 1: Adsorption energy of pentacene withthe molecular axis parallel to the 2-fold (a) 90◦
and (b) 180◦ directions.
We also studied the adsorption of oxygen
atoms on the (111) surface of the PdZn inter-137

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

metallic compound. We used the VASP codewith a slab model consisting of six atomic layers. Using a 2x2 unit cell containing eightatoms in the top layer, we investigated theadsorption structures for various oxygen coverages between 0.125 - 1.0 monolayer (ML).
Calculations were mainly carried out on thesystem-C as ﬂat MPI jobs with 4-9 nodes. At
high oxygen coverages, we found that Zn atoms
form a Zn-O layer with a hexagonal ring reminiscent of the Wurtzite ZnO. PdZn is known
as a candidate catalyst for the methanol steam
reforming (MSR) because of its higher thermal
stability than the traditional Cu-based catalyst. It has been reported that a ZnO patch isformed on the PdZn catalyst during the MSR
reaction, and considered that the high CO2 selectivity of PdZn in MSR is due to the formation of the ZnO patch[4]. However, there
is still no report of the ﬁrst-principles calculations about the formation of the Zn-O layer.

References[1] M. Sato, Master Thesis (Kagoshima University, 2022).
[2] A. Aloﬁ et al., Mat. Trans. 62, 312(2021).
[3] H. Takakura et al., Nature Matter. 6,58(2007).
[4] M. Friedrich et al., Angew. Chem. 52,
4389 (2013).

138

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Atom exchange of fcc metal surfaces
Tsuneo FUKUDA†
Department of Physical Electronics and Informatics, Graduate School of Engineering,Osaka City University, 3-3-138 Sugimoto, Sumiyoshi-Ku, Osaka 558-8585

The

migration ofan adsorbed atom

vacuum layer for both sides for both the(adatom) impinged from gas phase to a surface

Ni(110) and Ni(100) surfaces. After fully

is one of the fundamental processes for crystal

relaxed atomic structures of slabs for pristine

growth. In some adatom-surface combinations,surfaces, Ni or Cu adatom was located andinstead of migration of an adatom on theadsorbed surfaces were again relaxed bysurface, an adatom may be embedded in the

keeping atomic positions of bottom 4 layers

surface layer by displacing the surface atom,

fixed. Electronic and ionic iteration parameters

and the released atom migrates on the surface,

are listed in Table 1. The NEB method waswhich is often called "displacive adsorption".

performed for 5 images with both the initial and

Such

and

final structures being fixed (IMAGE=4). Near

surface atoms are frequently reported by using

the potential maximum, the migration potential

field ion microscopy (FIM) [1].

was further evaluated by dividing these images.

substitutions

betweenadatoms

We have investigated adatom migration andadatom-surface

atom

substitution

The obtained migration potential barriers

fromare listed in Tables 2 and 3. For the Ni(110)theoretical point of view. We chose the Ni(100)surface, a simple surface diffusion of adatom

and Ni (110) surfaces as model systems andalong [110] direction, path A in Fig. 1(a), had ainvestigated the migration and the substitutionof Ni and Cu atoms based on first-principles
calculations. We employed a popular firstprinciples calculation code VASP, and thenudged elastic band (NEB) method was applied
Fig. 1: Adatom migration pathways for (a) Ni(110) andNi(100) surfaces. The unit cells used in the structural
calculation are shown by the dotted lines. Initialadatom position is indicated by blue circle, and path A,
B, and C on the Ni(110) surface and path A and B on
the Ni(100) surface are indicated by purple, green, andyellow, respectively.

to the evaluation of adatom migration andexchange potentials [2].
Surface models are shown in Figs. 1. We
employed 7-layer slab model with ~0.8 nm†

Present address: Division of Electronics and Physics, Graduate School of Engineering, Osaka
Metropolitan University, 3-3-138 Sugimoto, Sumiyoshi-Ku, Osaka 558-8585

139

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Table 2 Migration potentials of Ni adatom
corresponding paths shown in Figs. 1

Table 1 Parameters for present calculation
VASP Version
slab geometry
k-points
potential
cut off energy
Fermi level smearing
convergence condition

SCF mixing parameters

5.3.5
Ni(110): 4×2 unit 7 layerNi(100): 3×3 unit 7 layerNi(110): 12×16
Ni(100): 12×12
PBE [3]384.08 eV
Methfessel & Paxton
 = 0.12 eVElectronic:10-6 eVIonic: 0.01 eV/Ang.
AMIX = 0.02
BMIX = 0.0001
AMIX_MAG = 0.08
BMIX_MAG = 0.0001(linear mixing)

Time period

POTIM = 1.0

Number of core576 (192/image)face
Ni(110)Ni(100)path
Path APath B
Path C
Path APath B

Potential energy [eV]
0.321
0.444
0.443
0.728
1.105

Table 3 Migration potentials of Cu adatom
corresponding paths shown in Figs. 1
face
Ni(110)Ni(100)path
Path APath B
Path C
Path APath B

Potential energy [eV]
0.309
0.426
0.428
0.523
1.080

report.[4]In conclusion, the present calculation

lowest energy. Atom exchange processes suchshowed that in addition to the surface migrationas path B and C have higher energy than path A.

of adatom such as Ni and Cu, atom exchange

Experimentally, only path C was seen by FIM

may be possible on the Ni(110) surface, but

study [1]. Although the present calculation ofslightly higher energy, i.e., higher temperature

the energy barrier of path C gave 0.36 mV

is needed for replace surface atom with adatom

lower than that of path B, but we failed toon the Ni(100) surface.

reproduce the experimentally observed large

References

difference in the exchange probability between[1] G. Antezak and G. Ehrlich, Surf. Sci.

path B and C [1]. This seems to indicate theReports 62 (2007) 39.

thermal effect, i.e., collective motion of[2] H. Jonsson, G. Mills and K. W. Jacobsen,

adatoms and surface atoms is crucial to theCondensed Phase Simulations, ed. B. J. Berne,

atom exchange process and molecular dynamic

G. Ciccotti and D. F. Coker (World Scientific,

study should be required for understanding this1998).

discrepancy.

[3] J. P. Perdew, K. Burke, and M. Ernzerhof,
Phys. Rev. Lett. 77 (1996) 3865.

For the Ni(100) surface, adatom diffusion
along the [110] direction, path A in Fig. 1(b),[4] B. Müller, B. Fischer, L. Nedelmann, A.

again had a lower energy than atom exchange

Fricke, and K. Kern, Phys. Rev. Lett. 76 (1996)
2358.

process path B, in agreement with the previous

140

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Development of the ﬁrst-principles method for
calculating Tc based on density functional theory forsuperconductors
Ryosuke AKASHI
Department of Physics, The University of Tokyo

The ﬁrst-principles prediction of the superconducting transition temperature Tc of metals
is a challenging goal in the condensed matterphysics. For phonon-mediated superconductors, accurate calculation schemes is now available on the basis of the Eliashberg theory anddensity functional theory for superconductors(SCDFT) [1, 2]. In this year, we conducted adevelopment of a ﬁrst-principles method withthe latter theory.
In the conventional superconductors, thephonon-mediated attraction is the dominantsource of the electronic pairing interaction,
where the electron-electron Coulomb repulsion
has only a secondary role that it partially cancels the attraction. On the other hand, as the
Coulomb interaction becomes stronger, there
emerges a chance for it to aﬀect the superconductivity in a diﬀerent way. A prominent example is the plasmon-mediated pairing
mechanism [3, 4]: the frequency dependenceof the screened Coulomb interaction enhances
the pairing. Another is the spin-ﬂuctuation
mechanism, where ferromagnetic ﬂuctuationstrongly suppress the singlet pairing [5, 6, 7, 8].
We have developed an SCDFT based method
to calculate Tc with the plasmon and spin ﬂuctuation eﬀects included. In this project, wemade further progress along this line. In the
SCDFT, the superconducting Tc can be evaluated by solving the following gap equation∆nk = −Znk ∆nk−

1∑tanh[(β/2)En′ k′ ]Knkn′ k′∆n′ k′ .(1)2 n′ k′En′ k′The kernels K and Z represent the eﬀectof the electron-phonon and electron-electron
Coulomb interactions. By developing new
forms of K we can in principle incorporate
novel interaction eﬀects. β is the inverse temperature and Enk is the Kohn-Sham energy
eigenvalue.
First, we constructed a formalism of theexchange-correlation functional in the SCDFT
based on the adiabatic connection formula [9].
We formally conﬁrmed that this formalism reproduces the previous one.
Next, we extended this formalism and derived a formula for the exchange-correlation
kernel beyond the local density approximation. As this is numerically too demanding,
we are continuing exploration of eﬃcient numerical implementation of the newly developed
formula.
We also tested an implementation of the
SCDFT-based Tc calculation in Superconducting Toolkit [10] to see if the previously
published values [8], calculated with anotherin-house code, are well reproduced. The calculated values (Table 1), though fairly agreeing
with the previously published ones, were generally larger. The overestimation is likely attributed to a subtle implementation details likehow to numerically integrate the internal frequencies. These calculation were performed in
System B (Ohtaka) with parallelization within
a single node. The elapsed time for each Tc

141

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

was approximately less than 30 minutes.
Table 1: Superconducting transition temperature Tc (K) calculated with the SCDFT gap
equation Eq. (1) using Superconducting
toolkit, with and without the spin ﬂuctuation eﬀects. In RPA and LDA, the exchange
correlation kernels entering the K formula [7]are ignored and incorporated, respectively.
RPA
LDA
w.o. SF w. SF w.o. SF w. SF
20.5
2.8
18.1
2.5
V
Al
1.2
0.6
0.9
0.4

References[1] M. Lüders, et al., Phys. Rev. B 72 (2005)
024545.
[2] M. A. L. Marques, et al., Phys. Rev. B 72
(2005) 024546.
[3] Y. Takada, J. Phys. Soc. Jpn. 45 (1978)
786.
[4] R. Akashi and R. Arita, Phys. Rev. Lett.
111 (2013) 057006.
[5] N. F. Berk and J. R. Schrieﬀer, Phys. Rev.
Lett. 17 (1966) 433.
[6] F. Essenberger, A. Sanna, A. Linscheid,
F. Tandetzky, G. Profeta, P. Cudazzo,
and E. K. U. Gross, Phys. Rev. B 90
(2014) 214504.
[7] M. Kawamura, Y. Hizume, and T. Ozaki,
Phys. Rev. B 101 (2020), 134511.
[8] K. Tsutsumi, Y. Hizume, M. Kawamura,
R. Akashi and S. Tsuneyuki: Phys. Rev.
B 102 (2020) 214515.
[9] D. C. Langreth and J. P. Perdew, Solid
State Commun. 17 (1975) 1425.
[10] M. Kawamura, R. Akashi and S.
Tsuneyuki, Phys. Rev. B 95 (2017)054506 .

142

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Calculation of photoelectron angle distribution in thephotoexcitation process on the organic molecules adsorbed
surface
Kaori Niki
Graduate school of science
Chiba University, Yayoi-cho, Inage-ku, Chiba 263-8522
In recent years, it has been reported electronic

Burke-Ernzerhof (PBE) formalism. We used aand physical properties of organic molecules

unit cell with lattice constants a = b = 3.53 Å,

which interact with the surface using Wavec = 20.0 Å, α = β = 90° and γ = 120°, and plane-

numberPhotoelectron

wave basis set with the projector augmented

Spectroscopy. We have investigated thewave (PAW) approach with an energy cutoff ofcopper-380 eV. The Brillouin zone integration was

monolayer substrate system. Monolayer TiSe2

performed on a Monkhorst-Pack of 18×18×1

have two typically phases depending on thegrid of k-points. The energy convergencetemperature. A (1 × 1) structure (normal phase)

criterion was chosen for the self-consistency

is shown at room temperature and a (2 × 2)

cycle was 1×10-7 eV. To account for Fermi

superstructure (CDW phase) is realized undersurface broadening, temperature parameter σ

the 200 K. Each phase has band gap of 83 meV

for Methfessel-Paxton smearing was appliedand 153 meV, respectively [1, 2]. We have[5].

adopted this interesting thin film as the

We investigated a local Coulomb repulsion ofsubstrate and reported the structural changes

Ti 3d electrons and the effect of the electronic

depending on the temperature in the previous

smearing parameter σ on the band gap. We

year.

In this year, we investigated thecalculated band gaps with different sets ofelectronic state of a (1 × 1) monolayer TiSe2.

Coulomb repulsion U and σ. Fig. 1(a) show a

Based on these results, we optimized thenegative bandgap without considering U andstructure of CuPc adsorbed on TiSe2.

σ. On the other hand, fig. 1(b) and (c) have

1. Band calculation of TiSe2 monolayer

positive bandgaps close to 0. It is found that

We performed density functional theory

both σ and U increase the gap between the top(DFT) calculations using the Vienna Ab Initio

of Se 4p band (Γ, red point) and the bottom of

Software Package (VASP) version 5.4.4 [3,4].

Ti 3d band (M, blue point). However, the rise

The

were

in σ push up the upper end of the Se4p band

Gradient

above the fermi surface, so Fig. (b) is also

Approximation (GGA) within the Perdew-

excluded. As a result, we find the parameter set

space

(Ⅱ)

exchangedescribed

byresolvedphthalocyanine

correlation
the(CuPc)

effects

Generalized143

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(U = 3.9 eV and σ = 0.01 eV) which reproduce

structure of the monolayer TiSe2 with lateralexperimental result [6]. We used these

dimensions of 693 Å2. For Cu atom of CuPc,

parameters for optimization of CuPc / TiSe2

U value was adjusted to 4.0 eV which could

structure.

reproduce the qualitative orbital ordering ofisolated CuPc. In addition, CuPc had an
unpaired electron, thus we considered spinpolarization. Ionic positions were optimized
until the residual force on each ion was lessthan 1×10-2 eV/Å. When starting with arandom wave function, it took about 12 hours
to converge one electronic state using 5 nodes
with total 144 cores. From this structure
optimization calculation, it is found that theadsorption distance between molecule andsubstrate is 3.16 Å, resulting in a flat molecule
structure.
We have reached the stage where it is possible
to calculate the Wave number space resolved
Photoelectron Spectroscopy. We are planningto perform this calculation next year.

References[1] J.-P. Peng et. al., Phys. Rev. B 91, 121113(2015).
[2] P. Chen et. al., Nat. Commun 6, 8943
(2015).
[3] G. Kresse and J. Furthmuller, Phys. Rev. B,
54, 11169 (1996).
Figure 1 Band dispersion of the (1×1)[4] G. Kresse and J. Hafner, Phys. Rev. B, 49

monolayer TiSe2 with different electronic

14251 (1994).

smearing parameter σ and effective Hubbard[5] M. Methfessel and A. T. Paxton, Phys. Rev.

parameter U for Ti.

B 40, 3616 (1989).

2. Structural optimization of CuPc/TiSe2[6] R. Bianco et.al., Phys. Rev. B 92, 094107

We expanded above model for an 8 × 8(2015).

144

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Physical properties of layered graphene
Yoshitaka FUJIMOTO
Depertment of Physics, Tokyo Institute of Technology
Oh-okayama, Meguro, Tokyo 152-8551

Graphene is one of the most promissing materials for electronics devices including ﬁeld effect transistors and molecular sensors. Bilayerof graphene would have rich properties because electronic structures of graphene change
depending on how two layers interact witheach other. Here, we study the stackingpattern eﬀects on the electronic transport andmolecular adsoprtion of bilayer graphene using a ﬁrst-principles density-functional calculation [1, 2, 3].
It has been reported that the bilayer
graphenes with three diﬀerent stacking patterns (AA, AA′ , and AB) have been synthesized and observed experimentally. Table 1
lists the adsorption energies and the binding
distances between the B atom and the adsorbed molecule for various molecules on the
B-doped bilayer graphenes with three diﬀerent
stackings. For all stacking patterns, CO andCO2 molecules are adsorbed not strongly but
weakly, whereas NO and NO2 molecules canbind strongly with the large adsorption energies and the short binding distances. O2 andN2 molecules weakly bind with small adsorption energies. Therefore, only NO and NO2
molecules in air are adsorbed with chemicalbonds on the B-doped bilayer graphene without depending on the stacking patterns.
We study electronic transport of the pristine
bilayer graphene for the three diﬀerent stackings. The conductance for the AB stacking has
the linear dispersion near E = 0 eV as in thecase of the monolayer. On the other hand, unlike monolayer, the slope of the conductancecurve between 0 eV and ∼ ±0.25 eV is diﬀerent from that between ∼ ±0.5 eV and ±1.0
eV, respectively. Due to the diﬀerence of theslope of the conductance curve in diﬀerent energy ranges, the transport properties induced

145

by the adsorption of molecules would be diﬀerent from that of the monolayer. The ﬁnite conductivity in the AB-stacked bilayer graphene
was reported theoretically and observed experimentally. The conductance of the AA-stacked
bilayer graphene seems to have constant value
from the energies E = −0.25 eV to +0.25 eV,
and it increases linearly below −0.25 eV andabove +0.25 eV. Unlike the AA-stacked bilayer
graphene, the conductance of the AA′ stacking
shows asymmetric behaviors near the Fermi
energy. Thus, it is interesting that the conductances of the bilayer graphenes with AA,
AA′ and AB stackings are diﬀerent from one
another.
We study the molecular adsorption eﬀects
on the transport of the B-doped bilayer
graphenes for three diﬀerent stackings, anddeﬁne ∆G(E) = (Ga − Gb )/Gb × 100(%),
where Ga and Gb are the conductances of the
B-doped bilayer graphene with and without
molecules, respectively. When NO and NO2
molecules are adsorbed on the AB-stacked bilayer graphene, the conductances at the energy E = 0 eV change with ∆G(0) = +8.2%
and +3.0%, respectively. For the AA stacking,the adsorptions of NO and NO2 molecules can
induce ∆G(0) = +2.1% and −17.4% in the
conductances, respectively. In the case of the
AA′ stacking, there are sizable conductancechanges of −17.2% and −47.7% for the adsorptions of NO and NO2 molecules, respectively.
In summary, the stacking-dependent behaviors of the transport and the adsorption
properties of molecules of bilayer graphene
have been investigated for various molecules
including environmentally polluting or toxic
molecules and common molecules (NO, NO2 ,CO, CO2 , O2 , and N2 ) using a ﬁrst-principles
density-functional study. Irrespective of stack-Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Table 1: Adsorption energy Ea (eV) and distance d (Å) between molecule and B atom for eachmolecule adsorbed on B-doped bilayer graphenes with AA, AA′ , and AB stackings.
Stacking
NO
NO2
CO
CO2
O2
N2
AA
Ea −1.20 −1.12 −0.12 −0.04 −0.16 −0.28
d
2.15
1.60
2.90
2.87
1.76
2.96
AA′
Ea −1.22 −1.15 −0.12 −0.04 −0.19 −0.28
d
2.15
1.60
2.90
2.87
1.75
2.96
AB
Ea −1.19 −1.13 −0.12 −0.03 −0.20 −0.27
d
2.14
1.60
2.89
2.84
1.75
2.93

Table 2: Percentage of changes in the conductance at the Fermi energy E = 0 [∆G(0)] foreach molecules adsorbed on B-doped bilayer
graphene with AA, AA′ , and AB stackings.
The changes in the conductance are calculatedwith respect to the B-doped bilayer graphene.
AA
AA′
AB
NO
2.1
-17.2
8.2
NO2
-17.4
-47.7
3.0

ing patterns, only NO and NO2 moleculesin air can bind strongly on B-doped bilayer
graphene. It is found that the conductancecurves of bilayer graphenes change depending on the stacking patterns and the presenceof the NO/NO2 molecule, and the variationof the conductances induced by molecular adsorptions could be used by gas sensing devices
such as ﬁeld-eﬀect transistors under low bias
voltages.

References[1] Y. Fujimoto and S. Saito, Appl. Surf. Sci.
Adv. 1, 100028 (2020).
[2] Y. Fujimoto, Mod. Phys. Lett. B 35,
2130001 (2021).
[3] Y. Fujimoto and S. Saito, J. Electrochem.
Soc. 169, 037512 (2022).

146

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Computational materials design of Ag and Cu chalcogenide
based thermoelectric materials
Ho Ngoc Nama) and Kazunori Sato a, b)
a)b)

Graduate School of Engineering, Osaka Univ., Yamadaoka, Suita, Osaka 565-0871

CSRN, Graduate School of Engineering Science, Osaka Univ., Toyonaka, Osaka 560-8531
The development of compact, wearable, andrespectively [2]. The effects of the doping ofespecially self-charging personal electronics

transition metal impurities were investigatedhas become more feasible in recent years. As aand the results were consistent with thecontinuous power supply for an uninterrupted

experimental observations [2]. We have also

operation, it is desirable that these electronic

found that experimental transport properties

devices can convert directly the waste heat ofwere reasonably reproduced by using thethe human body or environment into electricity

electron-phonon coupling method, showing thevia flexible thermoelectric (TE) generators. A

suitability of this scattering model in estimating

recent study discovered the extraordinary

the relaxation time [3].

ductility of the inorganic semiconductor silver

In addition to Ag chalcogenides, we havesulfide (Ag2S) [1] paving the way for flexible

also applied our calculation procedures to Cu

TE developments. In this study, we use first-

sulfide systems and found that the defectprinciples calculations combined with theformation energy and transport properties canBoltzmann transport theory to investigate thebe simulated by using the acanthite-likeelectronic and transport properties of Ag andstructure for this material. Based on these

Cu chalcogenide materials.

knowledges, we will extend our computational

First of all, we have chosen the exchange -materials design for (Ag1-x, Cux) (S1-y, Sey)2

correlation functional SCAN+rVV10 in orderchalcogenides to realize flexible TE materials.

to reproduce band gap energy and lattice

References

constants of a-Ag2S reasonably. According to[1] X. Shi et al., Nat. Mater. 17 (2018) 421.

the calculations of the formation energy of[2] H. N. Nam et al., Phys. Chem. Chem. Phys.

typical intrinsic defects, it is found that the23 (2021) 9773.

intrinsic defects greatly affect the system’s[3] H. N. Nam et al., Appl. Phys. Lett. Accepted.

conductivity where Ag vacancy and interstitial[4] H. N. Nam et al., Phys. Rev. B105 (2022)

Ag

075205.

act

asp-typeandn-typedefects,

147

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Thickness-Dependent Electronic Properties ofAtomically Thin Boron Monosulfide
Masayuki TOYODA
Department of Physics, Tokyo Institute of Technology
Ookayama, Meguro, Tokyo 152-8551

Rhombohedral boron monosulfide (r-BS) isa meta-stable phase of binary boron-sulfersystem obtained under high-temperature andhigh-pressure condition [1]. Recently, Kondo
et al. have succeeded to exfoliate atomicallythin layers of r-BS from its bulk.
We have studied the electronic structures ofthin layers of r-BS using the first-principles
density functional theory calculations. The
calculations have been performed using themassively parallel computers in ISSP, University of Tokyo.
From the electron localization function(ELF) analysis, we have found that B and S
atoms are covalently bonded in each S-B-B-S
unit layer (Fig. 1). The valence-band maximum (VBM) states are consisted of hybridized
B-pz and S-pz orbital states. By the stackingof layers, the VBM states exhibit splitting, resulting in the systematic narrowing of the bandgap (Fig. 2). This band-gap narrowing has also
been observed in the optical experiement [2]and is suggesting the potential application of
r-BS as 2D semiconductos with tunable band
gap.

Figure 1: Electron localization function along
the atomic bonds in (a) r-BS, (b) h-BN, and(c) MoS2 .

References[1] T. Sasaki, H. Takizawa, K. Ueda, and
T. Endo, phys. stat. sol. (b) 223, 29(2001).
[2] H. Kusaka, R. Ishibiki, M. Toyoda, et al.,J. Mater. Chem. A 9, 24631 (2021).

Figure 2: Calculated band gap values of thinlayer (1 < N < ∞) and bulk (N = ∞) r-BS.

148

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Large-Scale Molecular Dynamics Simulations on NonEquilibrium Phenomena Using First-Principles Calculationand Machine Leaning
Masaaki MISAWA
Faculty of Natural Science and Technology, Okayama University
Artificial neural network (ANN) potential,

movement of atoms occurs during the elastic-

which is an interatomic potential constructed bymachine-leaning,a

In order to improve the predicting ability ofpromising method to achieve extra-large-scale

ANN-potential, we introduce a more accurate

molecular dynamics (MD) simulation withtraining method that uses not only potential

first-principles accuracy [1]. Application of thisenergy but also atomic force and pressure as

ANN-MD to far-from-equilibrium phenomena,reference data [4]. Using this improved ANN-such

pressure-induced

potential, it was succeeded that reproducing thetransformation, is important to understandelastic-to-plastic transition behavior and plastic

structural properties of materials under extreme

deformation of α-quartz with an accuracy closeconditions. For this purpose, we have tried to

to FPMD simulation. Thus, the ANN-potential

perform ANN-MD simulation of shock-induced

trained with energy, force, and pressure is a

structural transformation of silica.

potentially powerful tool to investigate wideasattracts

fracture

andattentionas

to-plastic transition process.

To compute shock-compression behaviors

range of far-from-equilibrium phenomena andwithin the framework of MD method, thewill also provide useful information for solid

multiscale-shock technique [2] was employedearth physics field in future.

in our simulations. Potential energy of shockcompressed

α-quartz

obtainedbyfirst-

References

principles MD (FPMD) method is used as a[1] J. Behler and M. Parrinello, Phys. Rev. Lett.

reference data for training of ANN-potential.

98, 146401 (2007); N. Artrith and A. Urban,

As a result of ANN-MD simulation for elastic

Comput. Mater. Sci. 114, 135 (2016).

shock-wave
behaviorsofregion,elastic-deformation[2] E. Reed et al., Phys. Rev. Lett. 90, 235503

α-quartzwas(2003).

successfully

reproduced with high-accuracy [3]. On the[3] M. Misawa et al., J. Phys. Chem. Lett. 11,

other hand, for plastic shock-wave region, the4536 (2020).

ANN-potential[4] K. Shimamura et al., J. Chem. Phys. 153,has

completely

failedon

234301 (2020).

prediction of both structure and energy, becauseof that irregular structure appears and rapid149

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Defect formation energy and band structure of Pb-based
topological insulators
Yuki TOKUMOTO
Institute of Industrial Science,
The University of Tokyo, Komaba, Meguro-ku, Tokyo 153-8505

PbBi2Te4 topological insulator (TI) is

(DFT)

using

Viennaab-initio

simulationconsidered as one of the promising candidates

package (VASP) version 6.2.1. We adopted

for future spintronics material with a large spin

generalized gradient approximation (GGA)

current density. However, experimental studies

proposed by Perdew, Burke, and Ernzerhofof its surface transport properties have not yet(PBE) as the exchange correlation energy

been reported, owing to a lack of bulk

functional. For all calculations, we included

insulating crystals. So far, we achieved bulk

spin-orbit coupling (SOC). We obtained band

insulating states in Pb(Bi,Sb)2Te4 topological

gap of 115 meV for PbBi2Te4, which is muchinsulators by tuning compositional ratio of Bi

lower than experimentally reported value ~230

and Sb [1]. Although such control ofmeV. In order to estimate band gap morecompositional ratio leads to a change in theaccurately,

type and concentrations of intrinsic point(HSE06) and HSE03 hybrid functional was

defects, the Fermi level is not necessarily well

employed. As a result, band gap was further

tuned. For further success in achieving theunderestimated. We found that the HSE is notbulk-insulating state, not only the reduction ofeffective for accurate estimation of band gap of

absolute amount of charged defects but also thenarrow-gap semiconductors.

theHeyd-Scuseria-Ernzerhof

enhancement of the activation energy is

In 2022, we plan to calculate defectrequired. Recently, it was reported that Snformation energies using GGA-PBE including

doping for Bi2Te2Se TI is effective for lowering

SOC.

carrier concentrations and activated resistivities,
with a localized impurity band that acts as aReferencescharge buffer occurring inside the bulk band[1] Y. Hattori, Y. Tokumoto, and K. Edagawa:gap.

Phys. Rev. Mater. 1 (2017) 074201.

We carried out band structure calculation of
PbBi2Te4 based on density functional theory

150

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Electronic structure of the interfaces
between diamond and alkali metals
Toshihiro SHIMADADivision of Applied Physics, Faculty of Engineering
Hokkaido University, Kita-13 Nishi-8, Kita-ku, Sappooro, Hokkaido 060-8628
The interaction between diamond surfaces

functional. The

initialstructure

for

the

and liquid alkali metals has been little

optimization was created considering the bulk

investigated. We observed the etching process

lattice constants of the Li, Na, and K. 3 × 3 × 4

upon heating and the change in fluorescence ofunit cells of the diamond with 4 layers of alkali

DFT

metal atoms on both sides were employed. We

calculations using the ISSP supercomputer

tried several starting structures do the alkali

were performed to investigate the electronic

metals (lattice constant, orientation and shift),structure,theand the final structures were almost the same.

experimental results can be explained by the

The results are shown in Fig. 1. An upwarddifference in band bending.

shift in the density of states of interfacial

nitrogen

impurities

and

it

(NV-

was

centers).

shown

that

Experimental observation revealed that the

carbon is observed, and it is stronger at theetching of diamond occurred in heated lithium.
K and Na did not show etching.

The

fluorescence

was

of

the

NV-

centerinterface with (c) K.

significantly reduced only in contact withmolten K. This is understandable if we consider
that the negatively charged NV- that emits
fluorescence is lost due to upward band
bending. Although the extent of band bending
dependsonthedopingstateofFig.1 Density of states of carbon atthesemiconductor, the amount and direction of

diamondband bending can be estimated by determining

materials [1].

the positional relationship between the Fermi

Referenceslevel of the undoped semiconductor (in this[1] H. Takehana et al. Carbon 182 (2021) 585.

case diamond) and the metal using DFT[2] https://www.vasp.at

calculations. We used VASP[2] with PBE

151

interfaces

withvarious

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Automatic generation of first-principles effective models based on MLO
Hirofumi SAKAKIBARA
Advanced Mechanical and Electronic System Research Center (AMES),
Tottori University, Tottori, Tottori 680-8525

To thematerial informatics,

structure. Historically, we already have theautomatic generation of data is the central part

technique of Maximally Localized Wannierof calculations. According to the density

Functionfunctional theory (DFT) and local density

contains several difficulties or automation

approximation

calculateespecially for the system spanned by s- and p-

electronic structure, which generally controls

orbitals like semi-conductors. Therefore, fortheproperties,

calculating mixed system (superlattice and sosystematically without any empirical inputs.

on) containing not only the localized electrons

The DFT calculation itself is not so much

such as d-electron but also (s- or p-) delocalized

demanding in computational resources as longelectrons, MLWF is not a perfect method.

most

goal

of

(LDA),

part

of

we

can

material(MLWF)[1],

however,

MLWF

as we consider tens of atoms per unit cell.

For the goal of automation, we have already

Therefore, DFT calculations have been taken

constructed a new method named Muffin-tin

advantage of material informatics.

Like Orbital (MLO). MLO is given as a linear

The band structure within few eV measured

combination of atomic basis functions, which isfrom Fermi energy is effective for material

used for the first-principles calculations based

properties at low temperature. In other words,

on LMTO algorithm. Thanks to the PMT

wholeexcessmethod [2], we can take into account the effects

information. Even the drastic development ofof extended electron orbital such unoccupied

computational power, the computational coststates can be considered in the localized basis.

beyond the DFT calculations, for example

Here, PMT is a mixed basis function for first-

many body calculations, are still sometimes

principles calculation between Plane-wave anddemanding.

new

Linear Muffin-Tin orbital. By a kind offunctional material such as magnets and/or

weighted sum method, we can take into accountsuperconductors driven by electronic many

the effect of unoccupied orbital to the low-body effects, it is convenient that we canenergyautomatically construct a model which only

numerically.

band

structure

Therefore,

containsforminingcontains the low-energy part of electronicmodels’senergyeigen

valuesFor the benchmark materials, we have152

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

free electron’s Bloch orbital.

chosen four simple samples: Si(already shown
Fig.1 of the 2020’s report), SiO2(cubic),

For the future work, we have to simplify and

SrTiO3, and graphene. For all samples, thesophisticate MLO algorithm to save theautomation of model construction is succeeded.

computational times and to adopt for more

For SrTiO3, components described by plane-larger system such as superlattice, heterowave (namely, interference of unoccupied

junctions, and surfaces.

orbitals to the low-energy model) bases is few.
This can be interpreted that LMTO is asufficiently accurate basis function for SrTiO3.
For graphene, we can obtain a two-orbital
model which contains only C-pz’ π -bondReferences

orbitals forming the Fermi surface (or, Dirac[1] N. Marzari and D. Vanderbilt, Phys. Rev.

cone). Surprisingly, it is revealed that such π-

B 56, 12847 (1997).

bonds contains much plane-wave components.

[2] T. Kotani, H. Kino, and H. Akai, J. Phys.

For SiO(cubic), we find that an empty sphere inSoc. Jpn. 84, 034720 (2015).

a relatively large interstitial region is necessary
to reproduce the extended orbital shaping like a

153

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Theoretical Study on Electronic Properties
in New Nanoscale Surfaces and Interfaces
Katsuyoshi KOBAYASHI
Department of Physics, Ochanomizu University
2-1-1 Otsuka, Bunkyo-ku, Tokyo 112-8610

References[1] M. Bonilla et al.:(2018) 289.

Nat. Nonotech. 13[2] K. Sumida et al.: Phys. Rev. Materials 6(2022) 014006.

154

Energy (eV)

(a)

Γ

M

Γ

M

(b)
Energy (eV)We theoretically studied the electronic
structures of VSe2 monolayers.
Recentlymonolayers of transition-metal dichalcogenides(TMD) are extensively studied. Much attention is paid to magnetic properties of TMD
monolayers in recent years. Room temperature
ferromagnetism was reported for VSe2 monolayers in 2018 [1]. However, there is a controversy over its emergence. A recent experiment suggested that VSe2 monolayers are on
the verge of ferromagnetism [2]. We performed
a density-functional calculation on this system.
We used the program package VASP for thedensity-functional calculation. Figure 1 shows
the calculated result. Figure 1(a) and (b) show
band structures of the VSe2 monolayer in nonmagnetic (NM) and antiferromagnetic (AFM)
states, respectively. The magnetic structureof the AFM
state is shown in Fig. 1(c). We√√
used a 3 × 3 supercell of the primitive cellin the calculation. Γ is the center in the twodimensional Brillouin zone. M is the point at
the half of a primitive reciprocal lattice vector.
We found that the band structures of NM andAFM states are very similar.
In addition to the VSe2 monolayers we performed calculations of electronic structures of
V5 Se8 . We unfolded the calculated band structures in order to compare with experimental
band structures obtained from angle-resolved
photoelectron spectroscopy.

(c)Figure 1: Band structures of a VSe2 monolayer in nonmagnetic (a) and antiﬀeromagnetic(AFM) (b) states. The magnetic structure of
the AFM state is shown in (c). Red and green
circles show V and Se atoms, respectively. Arrows express the spin structure.

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles calculation of graphitization of diamondsurface and its exfoliation process
Kouji Inagaki
Graduate School of Engineering, Division of Precision Engineering and Applied Physics
Osaka University, Yamada-oka 2-1, Suita, Osaka 565-0871
Diamond has a wide bandgap electronic
property

and

high

mechanicalatoms (shown in (a) and (b)). This is carried out

strength;

by artificially adding a constant value of 0.002

application of diamond to semiconductor

Bohr to z axis of the terminal atoms in every

electronic devices or to mechanical cutting tools

quenched MD step with optimizing along x andare expected [1]. In electronic device appli-

y directions for the atom and other surface atoms.

cations, polishing is required to produce a flat

This procedure is repeated until saturated. The

surface, but diamond is the hardest material andbond breaking energy 2.1eV is obtained by theit is not easy. On the other hand, in toolenergy difference between (b) and (c), while cor-applications, a serious wear property in cutting

responding value of pristine surface is muchiron-based materials is reported. By utilizing thislarger energy of 3.2 eV. This indicates C-C bondwear phenomenon, it may be possible to realize

is effectively weakened by a Fe atom.

a novel diamond surface flattening method.

It is concluded that a Fe atom can play a

Currently, two possible causes of wear arecatalytic role in peeling graphite layer formed on

proposed: graphitization of diamond surface,diamond surfaces. Further investigations will be

and diffusion of C atoms into metal. This issue

required to evaluate metal element, atomic

remains unclear, particularly in the mechanism

environment effect, etc.

how strong C-C bonds break.
In this study, we investigated weakening of
C-C bond by the effect of Fe atom to confirm(a) Pristinegraphite layer peeling. Figure 1 (a) and (b)(b) Fe inserted (c) Peeled

Fig. 1: Peeling of graphite layer formed on dia-shows side views of the surface structures ofmond (221). Bright, dark and small atoms indi-

diamond (221) with and without Fe atom, andcate C, Fe and H atom. Atoms directed by arrow(c) peeled structure with Fe atom. (Peeled inare pulled-up to obtain peeled structures.

pristine structure is also perfomed but not

Reference

shown.) A DFT-based STATE code is employed[1] J.Wang, G.Zhang, N.Chen, M.Zhou, Y.Chen

to obtain stable structures and energies. Peeled

Int. J. Adv. Manufac. Technol. 113 (2021) 3027.

structures are obtained by pulling-up terminal

155

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

AB type 2D materials searchby high-throughput DFT calculations
Masahiro FUKUDA
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Recent rapid research progress of twodimensional (2D) materials have revealed that
2D materials are expected to be goodcandidates for future applications and devices.
In our previous research, we construct a
structure map [1] and a web-based database [2]of AB2 type monolayers by high-throughput
DFT calculations using OpenMX [3].
In this project, a structure map for AB typemonolayers of 2016 compounds which are all

Fig.1:

the combinations of 63 elements selected from

Examples

of

basic

planarstructures of AB type monolayers.

the periodic table (except for noble gases,

distortion, buckling and so on.

lanthanoids and actinoids) is constructed byhigh-throughput calculations based on the DFT

The genealogical classification helps us to(GGA-PBE) with the symmetry-unconstrained

understand phase transitions caused by latticegeometry optimizations and variational cell

matching with a substrate, stacking layers,

optimizations starting from several structures

external bending, tensile and compressive

such as ferromagnetic planar square lattice andstresses, and so on, when we invent newplanar honeycomb lattice structures as initialdevices.

states. We classified the converged structuresby planarity, the number of atoms in a unit cell,

References

2D Bravais lattice and space-group symmetry.

[1] M. Fukuda, J. Zhang, Y.-T. Lee, T. Ozaki,

The classification makes it clear that most ofMater. Adv., 2, 4392 (2021).

the AB type monoatomic layers are understood[2] Structure map and database of AB2 type 2D

as basic planar structures and their derivationmaterials (www.openmx-square.org/2d-ab2/).

patterns such as the shape of the primitive cell,[3] OpenMX, http://www.openmx-square.org

156

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Charge storage mechanism of Ti3C2Ch2 (Ch = O, S, Se, Te)
Kosuke Kawai, Masashi Okubo
Department of Electrical Engineering and Bioscience, School of Advanced Science andEngineering, Waseda University, 3-4-1 Okubo, Shinjuku-ku, Tokyo 169-8555

We have investigated the charge storage
mechanismsoftransition

metalof Ti27C18Ch18 (Ch = O, S, Se, Te), all Al atoms

carbidesand a Ti3C2 monolayer were removed from

(MXene) based on density functional theory

Ti6Al2C4, which was obtained from the(DFT) calculations. This year, we focused onMaterials Project (ID: mp-3747), and expanded

the dependence on termination function group

to a superstructure of 3×3×1 (Fig. 1a).

of a Ti3C2Tx monolayer.

Structure optimization calculations were carried

DFT calculations were performed with theout before and after adsorption of Zn atom on

Vienna ab initio simulation package (VASP).

A–C sites (Fig. 1b). Bader charge analysis was

The projected augmented wave method and aperformed using the script provided byplane basis set implemented in the VASP code

Henkelman’s group.3

were utilized with a cutoff energy of 520 eV

Adsorption at A and B sites are preferred toand k-point sampling on a 4×4×2 grid.1 The

that at C site (Fig. 2a) due to more coordination

generalized gradient approximation of thenumber with Ch atoms (3 for A and B site vs. 1

Perdew–Burke–Ernzerhof exchange-correlation

for C site). Averaged adsorption energy, whichfunctional was used.2 To prepare the slab model

would

be

inversely

proportional

Fig. 1 (a) The model structure of Ti27C18Ch18 (Ch = O, S, Se, Te) and (b) absorption
sites on the Ti27C18Ch18 monolayer. These structures were described using VESTA.4

157

to

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

charge/discharge rates, was in the order of Te <References

Se < O < S. Net charge transfer through Zn[1] G. Kresse, J. Furthmüller: Phys. Rev. B 54

adsorption increased in the order of Te < Se < S(1996) 11169−11186.

< O (Fig. 2b), which likely arises from more

[2] J. P. Perdew, K. Burke, M. Ernzerhof: Phys.

electronic interaction with more electronegative

Rev. Lett. 77 (1996) 3865−3868.

atoms. These results imply that Ch = O is the[3] W. Tang, E. Sanville, G. Henkelman:most preferable to achieve both a large

J. Phys.: Compute Mater. 21 (2009)

capacitance and fast charge/discharge rates

084204.

simultaneously in electrochemical capacitors[4] K. Momma, F. Izumi: J. Appl.

using Zn2+ ion as an adsorbed species.

Crystallogr. 41 (2008) 653−658.

Fig. 2 (a) Adsorption energy and (b) net
charge transfer for Zn atom on the
Ti27C18Ch18 (Ch = O, S, Se, Te) monolayer.

158

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Diffusion behavior of hydrogen in tungsten-rhenium alloy
Kazuhito Ohsawa
Institute for Applied Mechanics, Kyushu University,
Kasuga-koen, Kasuga, Fukuoka 816-8580

Introduction

Simulation Method

Tungsten (W) and its alloys are plausible

The stable positions of H in the W-Re alloys

candidates for plasma facing materials (PFMs)

and their energy levels are calculated in terms ofused in fusion reactors because of their excellent

first-principle calculations based on density

properties, e.g. high melting point, low hydrogen

functional theory. We use Vienna ab-initio(H) solubility, and toughness for irradiation.

simulation package (VASP). Large simulation

However, a large amount of H and H isotopes

cell composed of 5x5x5 bcc lattice (250 atoms)

are retained in vacancy (V) and vacancy typeare used to avoid the effects of periodic

lattice defects in the W specimen nucleatedboundary condition imposed on the simulation

under the irradiation circumstance. In particular,

cell. The cut-off energy of plane wave is 350eV.

tritium (T) retention in the PFMs is a serious

Re density is assumed to be low (less than 0.8

problem for safety operation of fusion reactors.

atomic percent) in the W-Re alloys. Besides,

By the way, rhenium (Re) is created in the Wenergy profiles of H along the migration paths in

specimen by nuclear transmutation in thethe W-Re alloys are estimated by nudged elastic

neutron irradiation circumstance in fusion

band (NEB) method.

reactors. The effects of Re addition on the W

We investigate energy profiles of H in two

specimen were investigated. Fortunately, thetype of W-Re model alloys, including one andmaterial properties are improved by the Re

two Re atoms, as shown in FIG. 1 and 2,

addition. For example, the H retention isrespectively.

decreased because

vacancy

nucleationissuppressed by the presence of Re. Besides, theResults

plasticity is also improved. Therefore, theH atom in the pure W and W-Re alloy is suited

application of W-Re alloys are positively

at tetrahedral interstitial site (T-site). Migrationconsidered for the PFMs instead of pure W.

energies are estimated as barrier heights for H to

In the present work, properties of H in W-Re

move to the neighboring T-sites. In pure W

alloys is investigated. In particular, diffusion

lattice, migration energy for H is estimated to be

behavior of H in the W-Re alloys.

about 0.21eV.

159

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

In FIG. 1, T-sites (a-i) in the vicinity of a Re

from the two Re atoms. So, the migration

atom and migration path of H from T-site (a) toenergies of H to move the neighboring T-site are(i) are shown. T-site (a) and (b) situated nearest

estimated to be 0.21-0.22 eV, which is almost

to the Re atom are geometrically equivalent.

same as that in pure W.

Energy profiles of H at the T-sites and along the
migration path A are shown in FIG. 3. H located

Summary and Discussion

at T-site (a) and (b) are about 0.1 eV unstable,

According to H energy profile, H will keepcompared with that located far enough away

away from Re in the W-Re alloys. Therefore, thefrom the Re. Therefore, T-sites close to Re arehigher Re density becomes, the lower H

supposed to be unfavorable for H.

solubility becomes. If Re density is low, theFIG. 2 shows likely migration path of H in the

migration energy of H in the W-Re alloy is

W-Re alloy including two Re atoms. The path isestimated to be about 0.21eV, which is similar toassumed to avoid T-sites close to the Re atoms.

that in pure W. However, the room available for

FIG. 4 shows energy profile of H along the

H diffusion becomes smaller by the Re addition.

migration path B from T-site (1) to (12). The

So, H diffusivity in W-Re alloy would be lower,

migration path B keeps an appropriate distancecompared with that in pure W.

FIG. 1: Schematic view of W-Re alloy
including one Re and H migration path A.

FIG. 2: Schematic view of W-Re alloy
including two Re and H migration path B.

FIG. 3: Energy profile of H along the path A in FIG. 1.

FIG. 4: Energy profile of H along the path B in FIG. 2.

160

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Clarification of atomistic mechanism application of process
design for adhesion interface between metal and plasmatreated fluoropolymers using first principles calculation
Yuji Ohkubo, Misa Nishino, and Kouji Inagaki
Graduate School of Engineering,Osaka University, 2-1 Yamadaoka Suita, Osaka 565-0871
Recently, increasing frequency of the digital

Our calculations were performed using thesignals is necessary because the amount of

STATE (Simulation tool for atom technology)

information

increasing.

code with plane-wave basis and ultra-soft

However, the higher the frequency, the greater

pseudopotential[2, 3]. The wave function and thetheTherefore,charge density of cutoff energies were 36 andpolytetrafluoroethylene (PTFE) which has low

400 Ry, respectively. Generalized gradientrelative dielectric constant and dielectric loss

exchange-correlation

tangent is considered as suitable as dielectrics ofPredew-Burke-Ernzerhof (PBE) was used forhigh frequency printed wiring boards. However,

the exchange-correlation functional and DFT-PTFE has low adhesion property, so adhering

D2 was used for the dispersion correction

PTFE to Cu is difficult. We have achieved

method[2, 3]. Cu surface was modeled Cu(111)

generation of oxygen-containing functional

4×4 slab composed of 3 atomic layers. As thegroups and adhesion strength of 0.94 N/mm

model of as-received PTFE, CF3–CF3 was

between PTFE and Cu when PTFE was He-

calculated. In addition, CF2=CF2, CF3–CF2–OO・,

plasma-treated at over 200℃ under atmospheric

CF3–CF2–OH, and CF3–CF=O were calculatedpressure[1]. However, the factor improving

as the model of plasma-treated PTFE. These

adhesion property, such as the chemical reaction

were the models of molecule containing oxygen-

at PTFE/Cu interface and the value of adsorption

containing functional groups observed by C1s-

energy has not been investigated. Therefore, in

XPS spectrum of the plasma-treated PTFE

this study, in order to explore functional groups

surface. Cu surface, fluorinated molecules, andcontributing to improvement in adhesion

adsorption systems were optimized then total

property, the reaction between Cu(111) surface

energies were calculated. The adsorption energy

and fluorinated molecules was analyzed and thewas defined as the difference between the sum

adsorption energy was calculated by means ofof total energy of Cu surface and fluorinatedfirst-principles calculations.

molecules and total energy of adsorption system.

communicationtransmission

loss

is

is.

161

functional

based

on

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Fig. 1shows theoptimized

Table 1 Value of adsorption energies on
Cu(111).

atomicgeometries for CF3–CF3. CF3–CF3 was notbonded to Cu surface and the value of adsorptionenergy was 0.25 eV in case of both on top siteand on bridge site. This is because C–F bonding
energy of 5.02 eV[4] is large and C–F bond isstable. After optimization for CF3–CF2–OO・, O–

COO・, the adsorption energy was the highest

O bond in CF3–CF2–OO・ was broken and CF3–

and the molecule was stable.

CF2–O・ and O were bonded to Cu surface,In this study, the functional groups

respectively. In this case, the sum of adsorption

contributing to improvement in adhesionenergies of CF3–CF2–O・ and O was calculated

property were explored by calculating theas the adsorption energy of CF3–CF2–OO・.

adsorption energies by means of first-principles

Therefore, the calculation of adsorption energy

calculations. From the result of calculations, theof CF3–CF2–O・ which is broken O–O bond atadsorption energy of molecules containing

initial state was added. Fig. 2 shows the

COO・ was the highest. Therefore, generating the

optimized atomic geometries for CF3–CF2–OO・

functional groups containing COO・ is possible

and CF3–CF2–O・. The adsorption energy wasto generate PTFE surface having high adhesion

3.10 eV, and the highest adsorption energy wasproperty. However, the calculations about Cu

obtained. Table 1 shows the value of adsorption

surface having different miller indices areenergies by calculation. The value of adsorption

needed because Cu surface has not only Cu(111),

energies of CF3–CF3, CF2=CF2, CF3–CF2–OO・,

but also Cu(100) and Cu(110). In addition, theCF3–CF2–O・, and CF3–CF=O was 0.25 eV, 0.68

calculations copper oxide such as CuO andeV, 4.73 eV, 3.10 eV, 0.60 eV, and 0.74 eV,

Cu2O are also needed because Cu surface isrespectively. In case of molecules containing

likely to be oxidized.

References[1] M. nishino et al., INTERFINISH2020
Abstract, C-SEP06-004 (2020) 210.
[2] K. Toyoda et al., J. Chem. Phys. 132 (2010)

Fig. 1 Atomic geometries for CF3–CF3
model (a) on top site and (b) on bridge site.

134703.
[3] K. Toyoda et al., J. Nanosci. Nanotechnol. 11(2011) 2836–2843.
[4] H. O. Pritchard et al., Chem. Rev. 55 (1955)

Fig. 2 Atomic geometries for (a) CF3–CF2–
OO・ model and (b) CF3–CF2–O・ model.

745–786.

162

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles electronic structure calculation of
double-perovskite photocatalyst
Kazume NISHIDATE
Department of Systems Innovation Engineering, Graduate School of Science and Engineering,
Iwate University, 4-3-5 Ueda, Morioka, Iwate, 020-8551

Double-perovskite oxides, A2 B ′ III B ′′ V O6 ,have been extensively studied due to their intriguing physical and chemical properties originating in their mixed valence nature[1]. Above
all, Ba2 B ′ III B ′′ V O6 was found to show an efﬁcient photo catalytic activity to dissolve water molecules into oxygen and hydrogen gases.
While TiO2 has been known to be an practical photo catalyst, it uses only the ultraviolet
rays which is less than 3 % of the sunlight pouring on the earth surface. On the other hand,the catalytic function of Ba2 B ′ III B ′′ V O6 double perovskite is activated under visible lightand hence has a great potential as a photocatalyst in daily use. In this work, we study
the surface electronic properties of Ba2 PrBiO6
double-perovskites. We used the DFT code of
Quantum ESPRESSO throughout our study.
PAW potentials were employed to represent
the core region and the PBE functional wasapplied for exchange correlation. We also usedthe HSE06 hybrid functional to evaluate theband gap values. We constructed slab system with three units along the z-axis exposing(001) surface. The vacuum region was taken to
be 20 Å height. Where two types of surfaces
can be exposed as the (001), namely, the BaO polar-surface and the Pr-Bi polar-surface.
We calculated the surface band structures andidentiﬁed the band gaps as well as their work
functions for these polar-surfaces.

1). Top components of the valence band arecomprised of oxygen orbitals, while the isolated bottom part of the conduction band isderived from the Bi orbitals. The surface bandgap value was estimated to be 1.93 eV within
the HSE hybrid functional level. The band gap
value of the bulk system with the HSE is foundto be 3.70 eV, therefore we think that the bulk
band gap is signiﬁcantly modiﬁed at the surface region. However, experimentally obtained
optical band gap value is 1.07 eV and which isstill half of the surface band gap obtained inthis study. We think the origin of the discrepancy is the B-site disordering (Pr↔Bi). Now,
we have been studying the eﬀect of the B-site
disordering by evaluating the band gap of thesystem.

Figure 1: Ba-O surface band structure.

References
As an example, we show the Ba-O surface band structure calculated within the
PBE exchange-correlation functional (Figure[1] S. Vasala and M. Karppinen: Prog. in
Solid State Chem. 43 (2015) 1.

163

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Electronic structure calculation for a transition-metalmonoxide using a four-atom-configuration model
Takayuki MAKINO
Research Center for Development of Far-Infrared Region,
University of Fukui, Bunkyo, Fukui, 910-8507

We have studied an electronic property of atransition-metalmonoxide.

Recently,

Ry which required a 6×6×6 k-point mesh.

thisAdopted pseudopotentials are based on themonoxide was reported to have an indirect gap

projector augmented-wave method.

whose energy is ca. 0.2 eV [1]. It is well known
that these compounds have been poorly

Introduction of Hubbard U – It should be

described by the density functional theory

noted that in addition to the Hubbard U on the d(DFT) with using commonly used approximate

orbital (Ud) in the transition-metal element, wefunctionals. It was very difficult for us to

introduced that (Up) on the p orbital of the

reproduce this insulating ground state of thisoxygen atom. An additional introduction of Upmonoxide, only with the simple DFT+U

is believed to be useful to open the gap (i.e.,

calculation. The word “simple” means anleading to a correct insulating ground state) due

application of only the on-site Hubbard energy

to the existence of the exchange-inducedon the d electrons. In this work, we explain

splitting. The value of Up (ca. 21 eV) has been

about several extensions [2] that we learned

determinedfrom literatures [2] and used to reproduce theresponse approach. On the other hand, if we

insulating ground state. We used the pwscf

tried to decide the value of Ud self-consistently,

codes included in the Quantum ESPRESSO

we could not reproduce the insulating ground(QE) package [3]. We performed the following

state theoretically. Therefore, as a compromise,jobs to obtain the theoretical electronic andwe set Ud of ca. 4.4 eV. We did not adopt thevibronic properties.

magnetic exchange J parameter for this

by

the

self-consistentlinearcalculation.
Self-Consistent Field Calculation – No matter
what physical quantity is calculated; one must

Spin-polarized DFT+U – The primitive cell offirst determine the electron density. We used a

the rock-salt structure is composed of two

50 Ry energy cutoff for the electronic wave

atoms (one is anion, while the other is cation).

functions and a small smearing width of 0.02

To implement the antiferromagnetic ordering,164

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

we adopted four-atom-configuration including

points. This is qualitatively consistent with thetwo cations with up and down spins. By

experimental observation [1]. The ground state

adopting this configuration, we can expect anobtained in the four-atom cell has slightly lower

effective lowering of the symmetry which may

energy per the cation-anion pair compared to

be useful to open the gap.

thenonmagnetic

andmetallic

two-atomprimitive cell, and thus the antiferromagnetic
Structural optimization – We also optimized

insulating ground state is energetically favored.

the structure by using the Broyden–Fletcher–

The result of the density-of-states calculation

Goldfarb–Shanno method. This turned out to besuggests that both the valence and conduction

essential to make the Fermi energy lie in the

bands are mainly composed of the d electron,

splitting region of the d electrons. Otherwise,

consistent with the narrow Mott-Hubbard gap

one cannot reproduce the insulating ground[1].

structure with the experimental lattice spacings.
atom
1
Tr[ns(na)] (up, down, total) =2.48039 1.01275 3.49314
spin 1eigenvalues:
0.292 0.292 0.448 0.449 1.000⁝
spin 2
eigenvalues:
0.151 0.151 0.192 0.259 0.259⁝

This resulted in the structure with the symmetry
C2h (2/m). It is known that the theoreticalground state depends on the starting eigenstates(occupation) during the spin-polarized DFT+U
calculations.

In

some

of

transition-metalmonoxides such as cubic FeO and CuO, thisFig. 1: An example of the orbital level output

adjustment is essential to reproduce theafter the initial eigenvalue settings are reflected.

insulating ground state [2]. This is also true inour case. The example of the output is shown inReferences

Fig. 1.

[1] K. Kaminaga et al., Appl. Phys. Lett. 108, 1
Band

dispersion

andelectronic(2016).

DOS[2] B. Himmetoglu et al., Phys. Rev. B. 84,

calculations –It should be noted that all of the115108 (2011).

above-mentioned extensions such as non-zero[3] P. Giannozzi et al., J. Phys.: Condens. Matter

Up, the spin polarization, structure optimization,

21, 395502 (2009).

and occupational adjustments turned out to be
essential to open an indirect gap (0.55 eV) in
our case with a transition between L and Γ

165

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Interactions in transition-metal-phthalocyaninemonolayer weakly bound to metal surfaces
Tetsuya ARUGA
Department of Chemistry, Graduate School of Science, Kyoto University
Kitashirakawa-Oiwakecho, Sakyo, Kyoto 606-8502

Phthalocyanine is a π-conjugated planarmolecule with a size of ∼1.2×1.2 nm2 . Two H
atoms at the center of the molecule can be substituted by a variety of transition metal ions
(M 2+ ), yielding a transition metal phthalocyanine (M Pc). The M 2+ ions are coordinatedby lone pairs of four N atoms arranged in aplanar tetragonal geometry. The adsorptionof M Pc on transition metal surfaces [1, 2] hasbeen widely studied experimentally.
We recently have studied experimentally theadsorption and ordering of metal phthalocyanine molecules on a bilayer indium ﬁlm fabricated on Si(111) surface [3, 4, 5] , on which theinteraction of M Pc with the substrate is expected to be weaker than the transition-metals
substrates. Our aim is to establish nearlyfree-standing transition metal phthalocyaninemonolayers eﬀectively decoupled from the substrate.
The experimental studies of M Pc adsorption on metal surfaces commonly show thatthe phthalocyanine molecules form monolayerswith long range order, in which the symmetry
axes of molecules are aligned at a particular
azimuthal angle which is not parallel to theaxes of the translational symmetry of the twodimensional lattice the molecules form. Thishas sometimes been ascribed to the direct interaction of nitrogen atoms (imide N) of themolecule with the substrate metal atoms.
Our computational result suggests that theelectronic interaction between phthalocyanineand In substrate is indeed week due to a large

van der Waals gap between the molecular plane
and the topmost In atomic layer of the substrate and further to the weak electronic interaction of the π orbitals of the molecules withthe sp electrons of the In bilayer.
In order to highlight the intermolecular interaction, we studied the “free-standing” M Pc
monolayer by ﬁrst-principles calculation by
VASP. The molecules were assumed to be ordered in a nearly square translational symmetry which exactly is the one observed on
In/Si(111). Note that the characteristics of themolecular arrangement is more or less similar
on other metal surfaces studied. The total energy of the molecular ﬁlm was calculated as afunction of the azimuthal rotation angle of themolecule. It was found that the total energy
is minimized at an azimuthal angle close tothat observed experimentally for the molecule
adsorbed on many metallic surfaces. This implies that the molecular conﬁguration in the
M Pc monolayer is governed by the intermolecular interaction. We have examined several
possibilities of the physical origin of the inplane intermolecular interaction between M Pc
molecules and succeeded in identifying the interaction that stabilizes the free-standing M Pc
monolayer.
We also studied the d electronic structureof M 2+ ions in the metal phthalocyanine adsorbed on In/Si(111) surface to probe the interaction with the substrate sp electrons. A
good correspondence between the density ofstates measured by photoemission experiments

166

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

with those obtained by the calculation yielded
an insight into the interaction between central
metal ions and metallic substrate in this system.

References[1] H. Koshida, Y. Takahashi, H. Okuyama,
S. Hatta, and T. Aruga, e-Journal of
Surface Science and Nanotechnology, 20(2022) 25.
[2] H. Koshida, H. Okuyama, S. Hatta, T.
Aruga, Y. Hamamoto, I. Hamada, Y.
Morikawa, J. Phys. Chem. C 124 (2020)
17696.
[3] S. Terakawa, S. Hatta, H. Okuyama, T.
Aruga, Phys. Rev. B 105 (2022) 125402.
[4] S. Terakawa, S. Hatta, H. Okuyama, T.
Aruga, Phys. Rev. B 100 (2019) 115428.
[5] J. W. Park and M. H. Kang, Phys. Rev.
Lett. 109 (2012) 166102.

167

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Electronic Structure and Fermiology of d-electroncompounds
Isaki Kinjo, Juno Iraha, Yasutomi Tatetsu1 , and Takahiro Maehira2
Graduate school of Engineering and Science, Nishihara, Okinawa 903-0213
1
Center for Liberal Arts Education, Meio University, Nago, Okinawa 905-8585
2
Faculty of Science, University of the Ryukyus, Nishihara, Okinawa 903-0213

represents that WC is a compensated metal.
We move to discuss the Fermi surface of WC.
Three bands cross the EF , therefore WC has
three Fermi surfaces in Figure 2. The color on
the Fermi surface illustrates the distribution of
the Fermi velocity components. Fermi surfaces
for WC centered at the A, K, and M points,
respectively. The 5th band has one dog’s bonelike hole sheet centered at the each Γ point.
The Fermi surfaces were calculated with FermiSurfer [3].

eVThe monocarbides of transition elementsfrom the group-IVB, -VB, and -VIB columnsof the Periodic Table are well known for theirgreat hardness and strength, high melting temperatures, metal-like thermal and electricalresistivities, nonstoichiometry, and superconducting properties [1, 2]. Most experimental
and theoretical eﬀorts to understand the fundamental properties of these technologicallyimportant materials have focused thus far on
the group-IVB and -VB carbides, all of which
crystallize with the cubic rocksalt structure.
The compound WC has two phases: a low temperature phase is hexagonal WC and a hightemperature phase above 2525◦ C is cubic WC.
We studied simple hexagonal WC with a spacegroup of P6m2 which contains one molecule in
the unit cell. The lattice parameters a and care 2.906 and 2.837 Å, respectively [1]. The
calculation for the energy band structures was
carried out by using the OpenMX code(http://www.openmx-square.org/). First, we discuss the calculated results for WC as shown in
Figure 1, in which we depict the energy band
structure along the symmetry axes in the Brillouin zone, which is shown in Figure 1, in theenergy region from -20.0 to 20.0 eV. The Fermi
level EF is at 0.0 eV. for WC and indicated bydashed lines in Figure 1. With respect to theenergy band structure near EF , we emphasizethat there is always hybridization between the
W 5d and C 2p states in WC. The total number of holes is equal to that of electrons, which

20
19
18
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
0
-1
-2
-3
-4
-5
-6
-7
-8
-9
-10
-11
-12
-13
-14
-15
-16
-17
-18
-19
-20

Γ

M

K

Γ

A

L

H

A,L

M,K

H

Figure 1: The energy band structure for WC,
the symmetry points and axis for the simple
hexagonal-type structure. EF indicates the position of the Fermi level.
ReO3 -type structures can be described as
ABX3 perovskites in which the A-cation site
is unoccupied. They therefore have the general composition BX3 , where B is normally acation and X is a bridging anion. The chem-

168

eV

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Figure 2: The Fermi surface for WC centeredat Γ point. Colors indicate the Fermi velocity
components on the Fermi surface.

25
24
23
22
21
20
19
18
17
16
15
14
13
12
11
10
9
8
7
6
5
4
3
2
1
0
-1
-2
-3
-4
-5
-6
-7
-8
-9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25

Γ

X

M

Γ

R

M

R

Figure 3: The energy band structure for ReO3 .
EF indicates the position of the Fermi level.
8

T otalRe- s
O 1- s
O2- s
O3- s
Re- p
O 1- p
O2- p
O3- p
Re- d
O 1- d
O2- d
O3- d
Re- f

7

6

5
D(E)(S tates /cell eV .)ical diversity of such structures is very broad,
ranging from simple oxides and fluorides, suchas WO3 and AlF3 , to more complex systems
in which the bridging anion is polyatomic, as
in the Prussian blue-related cyanides such as
Fe(CN)3 and CoPt(CN)6 . We found that theenergy bands near the Fermi level are mainly
because of the hybridization between the Re 5d
and O 2p electrons in Figure 3. Since the lowest twelve bands are fully occupied, as shownin Figure 3, the next three bands are partially occupied, while higher bands are empty.
This compound is compensated metal becauseof having two chemical units in the unit cell.
The calculated total density of states (DOS)for ReO3 is shown in Figure 4, in the energy
region from -25.0 to 25.0 eV. The total densityof states for ReO3 is shown by the solid blackline and the s, p, d and f states are shown by
the solid color lines, respectively. This figuretells us that the DOS of the d states is higher
than that of the f states at the Fermi level.
The Fermi surface of ReO3 is found to consistsof two hole sheets and two electron sheets.

X

4

3

2

1

0

- 20

- 10

0

10

20

Energy ( eV .)Figure 4: The density of states for ReO3 . EF
indicates the position of the Fermi level shownby the dashed line.

References[1] L. E. Toth:Transition Metal Carbides andNitrides., Academic, New York(1971).
[2] E.K. Storms:The Refractory Carbides.,Academic, New York(1967).
[3] M. Kawamura: Phys. Rev. B 89, 094515(2014).

169

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

First-principles Study of Defects of Magnesium Alloys
Daisuke Matsunaka
Department of Mechanical Systems Engineering,
Shinshu University, 4-17-1 Wakasato, Nagano 380-8553
atom, consistent with the previous DFT study.

Magnesium (Mg) has been of increasingly
viewpoint,

Several magnesium alloys include long

because of its low density and relatively highperiod stacking ordered (LPSO) structures

specific strength.

which consist of periodic arrangement of basalinterest

from

the

engineering

For molecular dynamics simulation widely

stacking faults and enrichment of solute atoms

used to obtain atomic-level information onin the vicinity of the stacking faults. The solute-strength and deformation of materials, anenriched layers give rise to a significant massinteratomic potential is very important. In thischange as well as local stiffness due to the L12-

study, a highly accurate interatomic potential oftype solute clusters, which affect the lattice

pure magnesium was developed based onvibration properties of the LPSO structures.

artificial

The

Our previous study showed that for a heavyreference data for training the ANN potential

enough mass change, a phonon-mediatedwere systematically constructed using density

interaction of the solute-enriched layers isfunctional theory (DFT) calculations. The

attractive and stabilizes the periodic stacking ofdeveloped ANN potential can accurately

the solute-enriched layers. In this study, we

estimate lattice constants and elastic moduli.

analyze phonon states of LPSO magnesium

Additionally,

generalized

stacking

fault

alloys, using first-principles calculations. We

energiesslipinmagnesium

calculate phonon density of states and phonon

calculated with the ANN potential are also in

dispersion based on the density functional

agreement with the DFT results. We also

perturbation theory.

developed the ANN potentials for several kinds

Referencesof magnesium alloys; Mg-Y, Mg-Ca, Mg-Al,[1] T. Nagata: Master Thesis, Shinshu Univ.

Mg-Zn. The ANN potentials of magnesium(2022).

alloys[2] D. Matsunaka, K. Oboso, Y. Shiihara and H.

neuralof

network

successfullysystems(ANN)[1].

described

interactions

between a basal dislocation and each solute

Mori: Proceedings of CMD2021, 204.

170