3.5

SCCMS Projects

359

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

First-Principles Phase Field Mapping
Kaoru Ohno
Department of Physics, Faculty of Engineering Science,
Yokohama National University, Tokiwadai, Hodogaya-ku, Yokohama 240-8501
In order to predict microstructures of alloys

effects. We succeeded in applying the first-

without relying on any empirical or experi-principles phase field method to this systemmental parameter, we have developed the first-

without using any thermodynamic parameter

principles phase field method by combining

and identifying the element partitioning

cluster expansion theory, potential renormali-

effect, i.e., enrichment of V in the β grainzation theory and density functional theory andboundary and Al in the α phase at 977 ºC.

applied it to NiAl alloys that are used for

Figs. 1(a) and (b) show the resulting V and Al

example in jet engine turbine blades at highconcentrations, while Figs. 1(c) and (d) showtemperature [1]. This year, we focused on Ti-

the corresponding experimental images [3],6wt%Al-4wt%V (Ti64) alloy, which is one ofwhich look very similar to each other. As a

the most used Ti alloys in industry in a broadconclusion, the present method has an ability

range from aerospace to medical applications.

to predict complex microstructure of alloys.

The mechanical properties and the phasetransformation of Ti64 are strongly influen-

References

ced by the microstructure, which is reliant on[1] S. Bhattacharyya, R. Sahara, and K. Ohno,

the history of processing, thermal treatment

Nature Communications 10, 3451;1-10 (2019).

procedures and alloying element partitioning[2] T. N. Pham, K. Ohno, R. Kuwahara, and S.
Bhattacharyya, J. Phys.: Cond. Mat. 32,

(a)

(c)

264001;1-9 (2020).
[3] S. Huang, et al. J. Alloys Compounds. 791,
575-585 (2019).

V

V(ex)

(b)

(d)

Fig. 1: (a) V and (b) Al concentrations obtainedby the present method [2], and (c), (d) the
corresponding experimental images [3].

Al

Al(ex)

360

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

Energy conversion --chemical energy-Osamu SUGINO
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581

Challenge

electrochemical

provide a hint to the hitherto

interface study is to establish a

unexplained structure of H/Pt(111).

quantum theory of electrode and to

It is also important to utilize thepredictofanultimately

efficientsimulation for electrocatalyst design.

energy conversion. To this end we

Platinum alloys are used in thestudied two topics, manifestation ofexisting fuel cells, but they are notnuclear quantum effect on H/Pt(111)

activeand environmental effect on theelectrocatalysts,

oxygen reduction reaction (ORR) on

oxides are attracting attention, buttransition metal oxides. The nuclear

thequantum effect was investigated byunderstood. In collaboration with

applying the first-principles pathexperimental groups, we calculated

integral

dynamics

the correlation of the adsorption(PIMD) under different coverageenergies of the ORR intermediates,

conditions [1]. The hydrogen atoms

O2H, OH, and O; the correlation isaredelocalized

known to be a descriptor of thequantum mechanically under low

activity. For most materials, thecoverageincreasingcorrelation curve does not cross thecoverage, the interatomic repulsion

point for ideal activity, but TiO2 andis

ZrO2 are found to cross it whenmolecular

known

found

to

but,
to

be
with

inducelocalization

enough.

mechanism

bettertransition

metal

has

some

not

noble

been

within an adsorption site; atop, fcc,

doped

or hcp.

The interaction induced

dopants [2]. The activity versus

localization, which had been showndefects relation is consistent withto exist in a low index noble metal

experimental findings. From thesesurfaces, was captured on platinum

results, the joint research group is

and this phenomenon is found tonow planning to realize an ideal

361

with

For

metal

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

catalyst.

C123 19486 (2019).

References[1] L. Yan et al. Phys.

Rev.

B101,

165414 (2020).
[2] Y. Yamamoto et al. J. Phys. Chem.

362

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

Molecular Dynamics Simulation Studies ofElectrochemical Properties of Ionic Liquid Electrolytes
Hayato SHIBA1,2 and Patrick A. BONNAUD1
1

Institute for Materials Research, Tohoku University, Katahira, Aoba-ku, Sendai 980-8577
2

Information Technology Center, The University of Tokyo, Kashiwa-no-ha,
Kashiwa, Chiba 277-8589

Ionic liquids (ILs) are salts that are liquid at

We have observed that ions migrate toward

the good reason that they are composed only ofelectrodes of opposite charge upon an appliedions, they are expected to be good candidates

voltage. The fluid structure in vicinal layers

for electrolytes in supercapacitors. They may

that span over roughly 1 nm with respect to theexhibit fast charging times, excellent powerelectrode surface is more pronounced than in

performances, and good aging properties. In

the neutral state. In the core of the gap between

this project, we focused on capacitors based onelectrodes, the higher the applied voltage, theproperties

more

of

Electrical

DoubleLayerstructured

thefluid.

Ions

formedCapacitors (EDLC), in which ionic liquid arealternative layers of cations and anions. Fromconfined between a pair of planar electrodes.

the average total charge sampled on one of theWe employed molecular simulation to study

electrodes, we computed the total gravimetric

underlying molecular mechanisms of electrical

capacitance and observed a fast increase at low

properties in EDLC. we selected an ionic pairapplied voltages before reaching a nearly

widelyconstant behavior beyond 1V. Those results set

studiedexperimentally:

methylimidazolium/

1-butyl-3-bis[(trifluoromethyl)acornerstone

for

the

understandingofsulfonyl] imide. Then, we inserted those ions inmolecular mechanisms underlying electricala simulation box made of two electrodes

properties of EDLC.

separated by a slit-shape gap of 4 and 6 nm up
to a mass density corresponding to given
thermodynamics conditions (300 < T < 400 K
and p~1 atm). We employed LAMMPS in
combination with the constant potential method
to simulate EDLC at the molecular scale, byusing a subroutine by Zhenxing Wang et al.1
Fig. 1: Snapshot of a molecular model for
simulations at constant potential.
1 https://github.com/zhenxingwang/lammps-conp363

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

Development of accuracy verification method on
analytic continuation results using cross validation
Kazuyoshi YOSHIMI1, Yuichi MOTOYAMA1, Yoshinori NAKANISHI-OHNO2
1Institute for Solid State Physics, University of Tokyo, Kashiwa 277-8581
2Graduate School of Arts and Sciences, University of Tokyo, Tokyo, 153-8902

Recentdevelopments

incomputationalverification on the results obtained by SpM

sciences have made it possible to construct andusing cross-validation.

analyze effective models for realistic materials.

In SpM, artificial oscillations appear in the low

Quantum Monte Carlo (QMC) calculations are

frequency range due to the cutoff of the bases.

one of powerful methods for finite temperature

To solve this problem, we focused on the Pade

analysis of effective models. To compare withapproximation,

experimental results such as excitation spectra,

analytic continuation in the low-frequency

it is necessary to connect the results of imaginary

region. In our method (SpM-Pade), the costtime simulations in QMC method to real time byfunction which becomes large when themeans of a numerical analytic continuation.

difference between the spectrum and that by

However, there is a problem that small statistical

Pade becomes large in the region where thenoise of the input data has a great influence onspectrum by Pade is stable to noise. We

which

gives

high-precise

the results of the analytic continuation, and there
has been no decisive method to enable highprecision spectral analysis until now.
Recently, we have proposed a method toautomatically select bases which are insensitive
to noise using sparse modeling and to achieve a
numerically stable analytic continuation against

Figure 1 Reconstructed spectrum by the Padenoise (SpM) [1-3]. In addition, SpM can easily(left), the SpM (center), and the SpM-Pade

add physically required conditions as constraints,(right) method. The line denotes the mean value

such as spectral non-negativity and sum rule. In

and the shaded region denotes the standard

this work, we have developed a new method todeviation from 30 independent samples. The

improve the validity of the results obtained byexact spectrum is shown as the black dashed

SpM. We also developed software to enable the

line.

364

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

examined the validity of SpM-Pade using thedeviation. Actual data, such as QMC data, havetest data obtained by adding noise to thea correlation between imaginary time data. In theartificially given spectrum and performing ancurrent simulation, we have ignored these

inverseeffects. We plan to apply SpM-Pade to the QMCtransformation

of

the

analytic

continuation. In order to see the sensitivity todata and investigate its applicability.

noise, the simulations were done for 30

For the cross-validation program, we develop

independent samples. In Fig. 1, we show thea program for large-scale parallel computers soreconstructed spectrum by the Pade (left), thethat it is designed not only to be used for analytic

SpM (center), and the SpM-Pade (right)

continuation but also for various methods. Now,

methods, respectively. It is seen that using Pade

we confirmed that it worked on the ISSP

approximation, the spectrum (solid line) is

supercomputer for the analysis in which thealmost same as the exact one (dashed line)

reduction of the noise of the image was donearound ω=0 and the standard deviation (shadedusing sparse modeling. We plan to release thisarea) is small. With increasing ω , both theprogram as open source software when it can bedifference between the spectrum and the exact

used for various solvers, including the program

one and the deviation become large toward theof SpM [2].

second peak aroundω=2 as expected. In the caseof SpM, the oscillations appear around ω =0

References

and the second peak becomes broad compared to[1] J. Otsuki, M. Ohzeki, H. Shinaoka, and K.

the exact one. However, the small variance

Yoshimi: Phys. Rev. E 95 (2017) 061302.

indicates that this method is insensitive to noise.

[2] K. Yoshimi J. Otsuki, Y. Motoyama, M.

In the case of SpM-Pade, due to introducing the

Ohzeki, and H. Shinaoka: Comput. Phys.

cost function, the oscillations around ω=0 is

Commun. 244 (2019) 319.

suppressed. On the other hand, the noise[3] J. Otsuki, M. Ohzeki, H. Shinaoka, and K.

tolerance is still strong as seen from the small

Yoshimi: J. Phys. Soc. Jpn. 89 (2020) 2061.

365

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

Ab initio study of quantum spin liquid in molecular
conductors
Takahiro Misawa
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

It has been experimentally reported that theX[Pd(dmit)2 ]2 molecular conductors exhibit
various phases such as antiferromagnetic, spinliquid, and charge-ordered states depending on
the choice of cation (X) [1]. In order to elucidate the origin of this variety of phenomena,detailed analyzes of band structures and transfer integrals based on the extended Hückel
method and first-principles calculations havebeen performed to quantitatively evaluate thematerial dependence of electronic states [2].
On the other hand, systematic analyses incorporating the electron correlation effect, whichplays an important role in the emergence ofvarious phases, have not been carried out due
to the lack of efficient tools and the needfor large computational resources for analysis
in molecular conductors. We note that thederivation of the low-energy effective Hamiltonian including the interaction parameters is
carried out only for EtMe3 Sb[Pd(dmit)2 ]2 [3],which has the quantum spin liquid groundstate.
In this project, using the open-source software package RESPACK [4, 5], we systematically derive effective Hamiltonians for all
available nine β 0 -X[Pd(dmit)2 ]2 (X = Me4 Y ,EtMe3 Y , Et2 Me2 Y and , Y = As, Sb, andP) [6]. In the part of obtaining the global band
structures, we use Quantum ESPRESSO [7].
Furthermore, using the quantum lattice model
solver HΦ [8, 9], we perform the numericallyexact analyzes for the obtained low-energy effective Hamiltonians and show that the mag-netic properties are well reproduced including
the quantum spin liquid behavior observed in
X=EtMe3 Sb. We also show that both the geometrical frustration and the off-site interactions play key role in stabilizing the quantum
spin liquid.
In addition to the ab initio study for themolecular conductors, we perform the ab initio
calculations for the high-Tc cuprates [10, 11].
We also analyze the quantum transport phenomena in the Weyl semimetals [12] and the
topological semimetals [13] using the realtime evolution of the quantum systems. Furthermore, using the many-variable variational
Monte Carlo method [14, 15], we analyze thequantum spin-nematic phase [16], the Kitaev
model under magnetic fields [17], and the dynamical properties of the correlated electron
systems [18].

References[1] R. Kato, Chem. Rev. 104, 5319 (2004).
[2] K. Kanoda, R. Kato, Annu. Rev. Condens. Matter Phys. 2, 167 (2011).
[3] K. Nakamura, Y. Yoshimoto, M. Imada,
Phys. Rev. B 86, 205117 (2012).
[4] https://sites.google.com/view/kazuma7k6r[5] K. Nakamura, Y. Yoshimoto, Y. Nomura,T. Tadano, M. Kawamura, T. Kosugi,

366

Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of TokyoK. Yoshimi, T. Misawa, Y. Motoyama,
arXiv preprint arXiv:2001.02351.
[6] T. Misawa, K. Yoshimi, T. Tsumuraya,
arXiv preprint arXiv:2004.00970.
[7] P. Giannozzi, O. Andreussi, T. Brumme,
O. Bunau, M. B. Nardelli, M. Calandra,
R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, et al., Journal of Physics: Condensed Matter 29, 465901 (2017).
[8] M. Kawamura, K. Yoshimi, T. Misawa,
Y. Yamaji, S. Todo, N. Kawashima, Comput. Phys. Commun. 217, 180 (2017).
[9] https://www.pasums.issp.u-tokyo.
ac.jp/HPhi/en/[10] M. Hirayama, T. Misawa, T. Ohgoe,
Y. Yamaji, M. Imada, Phys. Rev. B 99,
245155 (2019).
[11] T. Ohgoe, M. Hirayama, T. Misawa,K. Ido, Y. Yamaji, M. Imada, Phys. Rev.
B 101, 045124 (2020).
[12] T. Misawa, R. Nakai, K. Nomura, Phys.
Rev. B 100, 155123 (2019).
[13] T. Misawa, K. Nomura, Scientific Reports
9, 19659 (2019).
[14] T. Misawa, S. Morita, K. Yoshimi,
M. Kawamura, Y. Motoyama, K. Ido,
T. Ohgoe, M. Imada, T. Kato, Comput.
Phys. Commun. 235, 447 (2019).
[15] https://www.pasums.issp.u-tokyo.
ac.jp/mvmc/en/[16] T. Hikihara, T. Misawa, T. Momoi, Phys.
Rev. B 100, 214414 (2019).
[17] K. Ido, T. Misawa, Phys. Rev. B 101,
045121 (2020).
[18] K. Ido, M. Imada, T. Misawa, Phys. Rev.
B 101, 075124 (2020).

367

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

Unified Photonic-Electronic Devices
Kazuhiro YABANA
Center for Computational Sciences,
University of Tsukuba, Tsukuba 305-8577
now possible to calculate the ground state andWe develop a first-principles computational
dynamics

electron dynamics of systems of more than a

induced by ultrashort laser pulses based on

few thousand atoms. We intend to make thetime-dependent

new version public in early fiscal year of 2020.

methodtoinvestigate

electrondensityfunctional

theory(TDDFT) in real time. We develop the code
SALMON (Scalable Ab-initio Light-Matter

References

simulator for Optics and Nanoscience) [1] and[1] M. Noda, et.al, “SALMON: Scalable Ab-

make it open to the public at our website,

initio Light–Matter simulator for Optics andhttp://salmon-tddft.jp.

Nanoscience”, Comp. Phys. Comm. 235, 356

SALMON is developed as a unified(2019).

software that include traditional computational[2] K. Yabana, T. Sugiyama, Y. Shinohara, T.

methods

theOtobe, G.F. Bertsch, “Time-dependent density

finite-functional theory for strong electromagnetic

difference time-domain (FDTD) method, andfields in crystalline solids”, Phys. Rev. B85,

linear response calculation for susceptibities

045134 (2012).

based on TDDFT. There are two options for the[3]

coupling: macroscopic[2] and microscopic [3].

Yabana, “Time-dependent density functional

In addition, Ehrenfest molecular dynamics cantheory for interaction of ultrashort light pulse

be combined in the multiscale simulation [4].

with thin materials”, Phys. Rev. B98, 245147

oflight-matterelectromagnetism

analysisinteraction,
usingThis year, we have concentrated on the

S. Yamada, M. Noda, K. Nobusada, K.

(2018).

development of ver. 2 of SALMON. In the new[4] A. Yamada, K. Yabana, “Multiscale time-

version, we intended to improve the readability

dependent density functional theory for aand extensibility of the code by using structuresunified description of ultrafast dynamics:of Fortran extensively. We also paid much

Pulsed light, electron, and lattice motion ineffort to improve the efficiency and scalability

crystalline solids”, Phys. Rev. B99, 245103

for large-scale calculations. Using the code, it is(2019).

368

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

Dynamical DMRG study of magnetic excitations inmagnetization plateaus of a frustrated spin ladder
Takami TOHYAMA
Department of Applied Physics, Tokyo University of Science, Tokyo 125-8585

Magnetization plateau (MP) emerges inquantum spin systems due to spontaneously
breaking of translational symmetry. The broken symmetry can induce reconstruction ofelementary excitations, but its microscopic
mechanism and reconstructed quasiparticle
in MP phases have remained unclear. We
theoretically study magnetic excitations in
the MP phases of a frustrated spin ladder(FSL) by using the dynamical density-matrix
renormalization-group (DDMRG) method [1].
The Hamiltonian of FSL is deﬁned as H =∑H⊥ + H∥ + HZ with H⊥ = J⊥ N
i=1 Si,1 · Si,2 ,∑ ∑∑H∥ = 2η=1 Jη i 2j=1 Si,j · Si+η,j , and HZ =∑z , where S−h i,j Si,ji,1 (Si,2 ) is the S = 1/2
spin operator on ith rung in the upper (lower)
chain. Exchange energies of the ﬁrst-neighborbond in a leg, the second-neighbor bond in aleg, and the ﬁrst-neighbor bond in a rung, aredenoted by J1 , J2 , and J⊥ , respectively. We
set J1 /J⊥ = 0.2 and J2 /J⊥ = 0.65. Changing magnetic ﬁeld h, we obtain MPs at m =M/Ms = 1/3, 1/2, and 2/3, where M (Ms ) ismagnetization (saturation magnetization).
We calculate dynamical spin structure factor(DSSF) using DDMRG, where the correction
vector is expanded by the Legendre polynomialwith a Gaussian averaging [2] and a broadening factor is replaced by the width of the Gaussian, for which we set 0.02J⊥ . The system size,
N = 48 rungs, is suﬃcient to discuss dynamical behaviors. The number of states kept in
the DDMRG is m = 600, leading to truncation
error less than 1 × 10−4 .
Figure 1 shows DSSF for momentum along

Figure 1: DSSF for qy = π in the m = 1/3 MP.
rung qy = π as a function of momentum alongleg qx in the case of m = 1/3 [1]. A dispersive feature with zero-energy excitation at
qx = 2π/3, indicating a period with three times
the original unit-cell length in real space.
We also perform analytical approaches withthe perturbation theory to obtain an intuitive
view of magnetic excitations. A comparison
between numerical and analytical results indicates the presence of a reconstructed quasiparticle originating from spontaneously broken
translational symmetry, which is realized asa collective mode of the spin trimer called atrimeron [1].

References[1] K. Sasaki , T. Sugimoto, T. Tohyama, and
S. Sota: Phys. Rev. B 101 (2020) 144407.
[2] S. Sota and T. Tohyama: Phys. Rev. B
82 (2010) 195130.

369

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

Conversion and storage of energy−fuel cells andsecondary batteries: Research and development offundamental technologies of battery simulators.
Susumu OKAZAKI
Department of Materials Chmistry, Nagoya University
Furo-cho, Chikusa-ku, Nagoya 464-8603

The goal of our project is to develop thebasic technology of the whole battery simulator. One of the key techniques is molecularlevel design of polymer membranes controllingtransportation of protons and ions across themembrane with proper stiffness resistant tomechanical deformation under external stress.
Such membranes are widely applicable to thefuel cells used in the industrial products.
We previously performed fully atomistic
molecular dynamics (MD) calculations of hydrated perfluorosulfonic acid (PFSA) ionomers
composed of a hydrophobic polytetrafluoroethylene backbone with hydrophilic side
chains terminated by sulfonic acid, as a modelof proton exchange polymer electrolyte membrane of fuel cells[1, 2, 3]. On the system
B and C, we further performed a series of
MD calculations of different equivalent weightPFSA membranes containing H2 molecules atdifferent water uptakes[4] (examples are shownin Figs. 1 (A)–(C)). They explored atomistic
detail of H2 permeation through the membranes. The local semicrystalline structure of
the PFSA polymer and the morphology of thewater clusters in the membrane were found toaffect H2 permeation.
Moreover, new algorithms relating to thefast multipole method which enable to treat
two dimensional periodic boundary conditions(slab geometry)[5] and rectangular unit cellwith an anisotropic cell-partitioning[6] weredeveloped. They make it possible to perform

370

Figure 1: Schematic images of hydrated equivalent weight PFSA membranes with (A) 6, (B)16, and (C) 30 wt % of water uptakes[4].
MD calculations for more general calculation
systems by the software MODYLAS.
These results will contribute to the development of higher performance fuel−cells andsecondary batteries and to the realization of
the entire battery simulator in the near future.
References[1] A. T. Kuo, W. Shinoda, S. Okazaki : J.
Phys. Chem. C 120 (2016) 25832.
[2] A. T. Kuo, S. Okazaki, W. Shinoda : J.
Chem. Phys. 147 (2017) 094904.
[3] A. T. Kuo, K. Takeuchi, A. Tanaka, S.
Urata, S. Okazaki, and W. Shinoda :Polymer 146 (2018) 53.
[4] K. Takeuchi, A. Kuo, T. Hirai, T. Miyajima, S. Urata, S. Terazono, S. Okazaki,
and W. Shinoda: J. Phys. Chem. C, 123(2019) 20628.
[5] N. Yoshii, Y. Andoh, S. Okazaki: J. Comput. Chem., 41 (2020) 940.
[6] Y. Andoh, N. Yoshii, S. Okazaki: J. Comput. Chem., 41 (2020) 1353.

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

Development of Fragment-Based GW/BSE Method
and Application to Organic Optoelectronic Materials
Takatoshi Fujita
Institute for Molecular Science
Myodaiji, Okazaki, Aichi 444-8585

Accurate calculations of electronic states
are essential for computational studies of organic materials which are directed toward understanding of fundamental processes in organic electronic devices. Many fundamental
processes governing the device operation takeplace at the variety of interfaces, such as metalorganic and organic donor-acceptor interfaces.
Although developments of quantum chemistry
or ab initio program packages have enabled
one to routinely run a calculation for an isolated molecule or a molecular crystal, performing such a calculation for interface systems isstill a challenging issue.
We have recently developed the large-scale
GW method [1, 2] which can be applied
to large and disordered molecular systems.
Our implementation is based on the fragment molecular orbital (FMO) method, inwhich an entire system is ﬁrst divided into
many subsystems refereed as fragments, andtotal energy or physical properties are approximated from molecular orbital calculations
for fragment monomers, dimers, and optionally, trimmers. In the FMO-GW method,
the polarization function of an entire system
is approximated from molecular orbital localized within fragment monomer or dimers.
In addition, the ∆COHSEX approximation
is employed, in which dynamically-screened
Coulomb potential is explicitly evaluated for atarget fragment, while the statically-screened
Coulomb potential of entire system is evaluated at the static COHSEX level. Basedon these development, we have performed thelarge-scale GW calculation for systems whichcontain more than 1,000 atoms, such as theorganic semiconductor thin ﬁlm [1] and thedonor/acceptor interface [2]. In addition, we
have also developed the GW/Bethe-Salpeter
equation method which reasonably describes
excited states in polarizable molecular environments.
As applications of the FMO-GW, we have
explored the energy levels and charge-transfer
excited states in pentacene/C60 bilayer heterojunctions. In particular, we have investigated the dependence of interfacial morphologies on the electronic states and the eﬀects ofpolarization and delocalization eﬀects. Calculated energy levels and the excitation energies of charge-transfer states are in reasonableagreement with those estimated from ultraviolet photoelectron spectroscopy and from external quantum eﬃciency measurements. We
have found that the electrostatic contributionof polarization energies, which arises from pentacene quasrupole moments, governs dependence of interfacial morphologies on the electronic states.

References[1] T. Fujita, Y. Noguchi, Phys. Rev. B 98
(2018) 205140.
[2] T. Fujita, Y. Noguchi, T. Hoshi, J. Chem.
Phys. 151 (2019) 114109.

371

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

Data-driven materials design for new functional highentropy alloys
Tetsuya Fukushima
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

Recently, the constructions of materials
databases using first-principles calculations
have been actively carried out.
Combining such materials databases and data-mining
technique, one can not only analyze the mechanisms of physical phenomena, but also design new functional materials. Actually, there
are several large materials databases, such as
Materials Project, Open Quantum Materials
Database (OQMD), and Novel Materials Discovery (NOMAD), which researchers can freely
access. However, these inorganic materials
databases are only for compounds with stoichiometric compositions and do not contain
information about configurational disorderedsystems
In this year, in order to construct the materials databases for configurational disorderedsystems, we have developed an automatic highthroughput calculation method on the basis of
the Korringa-Kohn-Rostoker (KKR) Green’s
function method. There are several merits
in the KKR Green’s function method, compared to other first-principles approaches. The
most important advantage is the ability to perform the calculations of alloys and impurity
doped systems, combining with the coherent
potential approximation (CPA). Additionally,
since the Green’s function is directory calculated in this method, physical quantities, e.g.,magnetic exchange interactions and transport
properties, can be eﬃciently calculated by thelinear response theory.
Automatic high-throughput calculations by

372

first-principles approaches are not easy tasks,
because we need to appropriately control many
numerical parameters and self-consistent procedures. For example, in the KKR Green’s
function method, the energy integration of theGreen’s function is performed in the complex
energy plane to obtain the electron density.
This complex energy contour must cover thevalence bands and depends on systems. Our
python interface for the high-throughput calculations enable to determine the numericalparameters and to manage the self-consistentprocedures completely automatically. For theelectronic structure calculations, the AkaiKKR
program package is employed.
The automatic high-throughput calculations
have been demonstrated for quaternary highentropy alloys with BCC and FCC solid solution phases, where 4 principal elements havethe same atomic concentration and are randomly distributed in the crystals. The numberof the target elements are 38 as follows: Al, Si,
Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, Zn, Ga, Ge,
Y, Zr, Nb, Mo, Tc, Ru, Rh, Pd, Ag, Cd, In,
Sn, Hf, Ta, W, Re, Os, Ir, Pt, Au, Hg, Tl, Pb
and Bi. We randomly choose the 4 elementsfrom the above target elements and constructa quaternary high entropy alloy. Since boththe BCC and FCC cases are considered, thetotal number of the systems is 147,630. All
calculations are done by the ISSP supercomputer, system B.
We succeeded in automatically converging
99.99% of the systems, i.e., 147614/147630.

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

Figure 1 shows a part of our high-throughput
calculations, in which the magnetic moments
and Curie temperatures for the quaternaryhigh entropy alloys with the BCC phase areplotted. The Curie temperatures are calculated by the magnetic exchange exchangeinteractions and mean-field approximation.
Such database is quite useful for the design of
new functional materials. For example, in thecases of high-performance soft magnets, thehigh magnet moments and Curie temperatures
are needed. Using the database, we can screen
the candidate systems immediately. Here, it isconcluded that FeCoXY, FeCoNiY, and MnFeCoY systems are good candidates for the highperformance soft magnets

Figure 1: Automatic high-throughput calculations for quaternary high entropy alloys with
the BCC phase.

References[1] T. Fukushima, H. Shinya, A. Masago, K.
Sato and H. Katayama-Yoshida: Appl.
Phys. Exp. 12 (2019) 063006.
[2] S. Yamada, S. Kobayashi, A. Masago, L.
S. R. Kumara, H. Tajiri, T. Fukushima,
S. Abo, Y. Sakuraba, K. Hono, T. Oguchi,

373

and K. Hamaya: Phys. Rev. B 100 (2019)195137.
[3] A. Masago, H. Shinya, T. Fukushima, K.
Sato, and H. Katayama-Yoshida: AIP
Advances 10 (2020) 025216.
[4] T. Nunokawa, Y. Fujiwara, Y. Miyata,
N. Fujimura, T. Sakurai, H. Ohta, A.
Masago, H. Shinya, T. Fukushima, K.
Sato, and H. Katayama-Yoshida: J. Appl.
Phys. 127 (2020) 083901.

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

Electron Theory on Secondary-Battery Materials
Hiroki KOTAKA1 , Motoyuki HAMAGUCHI2 , Huu Duc LUONG2 ,Hiroyoshi MOMIDA1,2 , Tamio OGUCHI1,2
1
ESICB, Kyoto University, Katsuragoryo, Kyoto, Kyoto 615-8245
2
ISIR, Osaka University, Mihogaoka, Ibaraki, Osaka 567-0047

Microscopic mechanism of charging anddischarging reactions in several battery systems is studied by first-principles calculations to explore novel secondary battery materials. In this year, we further extended
our previous studies on Na/NaFeSO4 F and
Li/Li2+2x Mn1−x Ti1−x O4 systems.
Recently NaFeSO4 F has been expected to
be a candidate cathode material for the
next-generation Na secondary batteries. In
this first-principles study [1], the structural stability and voltage-capacity profiles of
Na/NaFeSO4 F system are investigated. Calculated total energy of NaFeSO4 F is comparedwith those of NaF and FeSO4 and found to beslightly stable to the phase separation. X-ray
absorption spectroscopy (XAS) spectra at theK-edges of Fe, Na, and F for NaFeSO4 F andFeSO4 F in the fully discharged and charged
states and of NaF and FeSO4 are computedwith the Fermi golden rule within the electric
dipole and quadrupole approximation. The
XAS spectra at Fe K-edge reveal that the discharged state in NaFeSO4 F and the charged
state in FeSO4 F are dominated by Fe2+ (d6 )and Fe3+ (d5 ), respectively. It is also foundthat the XAS spectra near K-edge of F in addition to Fe provide crucial information concerning the local structure in the relevant phases
existing during the reaction processes.
Li-excess cation-disordered rock-salt oxides
Li2+2x Mn1−x Ti1−x O4 have recently attracted
much interests as high-capacity cathodes associated with two-electron reactions. In these oxides, cations are disordered on the octahedral

sites coordinated by six oxygen ions and theoctahedra are expected to be stable during discharging and charging reactions. To clarify themechanism of the reactions, cathode propertiesof the oxides with the Li-excess of 0≤x≤0.3
are investigated by using first-principles calculations [2]. To evaluate the structure stability and reaction equations of the cathodematerials for each x, formation enthalpies of
Li2+2x−y Mn1−x Ti1−x O4 are calculated as afunction of Li-removal amounts y considering
several Li configurations. To model the oxides
with disordered cations and several Li concentrations x and y, supercells including 80 atomsof rock-salt lattice are assumed with the special quasi-random structure method. Basedon the estimated reaction equations, voltagecapacity profiles are obtained for each x andcompared with experiments. By analyzing theelectronic structures, roles of cation and anion
redox reactions depending on y as well as x
are discussed in detail. We discuss energy stability of O2 -release from the cathode during
charging/discharging processes by calculating
vacancy formation energies of O and Li. The
best x values in terms of voltages, capacities,
and stability against the O2 -releases are proposed.

References[1] H, Momida et al., J. Phys. Soc. Jpn. 88,
124709 (2019).
[2] M. Hamaguchi, et al., Electrochimica Acta
330, 135286 (2020).

374

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

Ab initio studies toward functional nanomaterials
based on abundant elements
Tetsuya TAKETSUGU
Faculty of Science, Hokkaido University
N10W8, Sapporo, Hokkaido 060-0810

Ab initio studies on various nano materials
are important to understand and design new
functional materials with abundant elements
for reducing costs and dependences on precious
metals used in many important aspects in ourdaily lives and therefore toward sustainable society. Heterogeneous catalysts rely mostly on
the expensive and rare metals such as Pt, Rh,
Pd, or Ru in automotive gas exhausts control,
water gas shifts, fuel cells, as well as combustive decomposition of ammonia. Platinum andpalladium are also used for hydrogen based energy strategies. Due to their high cost andlimited amounts, to reduce or even replacethese metals is emergent issues in industry.
To that end, we perform density functional
theory (DFT) computations under the periodic boundary conditions together with projector augmented wave method with VASP using ISSP supercomputers to gain chemical insights to understand the chemical mechanisms
in conventional catalysts as well as to designnovel catalysts with abundant elements.
We have developed a surface adsorption
model calculation database toward applicationto activity prediction of heterogeneous catalysts [1], where band calculations with SIESTA
are performed for various pure and alloy metal
surfaces adsorbed with molecules and atoms
relevant to important catalysis such as water gas shift. As a new functional nanomaterial, we have studied structural and electronic properties of borophene, an artificial
two-dimensional materials, whose structuraland electronic properties are not typical for
their bulk counterparts [2]. Scanning tunneling microscopy and density functional theory
calculations show that this structure forms asa single phase on iridium substrate in a wide
range of experimental conditions and maybe
then decoupled from the substrate via intercalation. We foresee that the adjacent borophene
sites with opposite electron doping might exhibit catalytic activity or facilitate highly regular adsorption of large molecules, metal clusters, and other objects interesting for catalysis,
sensors, and nanotechnology.
The excited states dynamics and spectroscopy are also important for the study on
functional nanomaterials. For spectroscopy,
the roles of silver nanocluster and plasmonic
nanocavity for surface- and tip-enhanced Raman spectroscopies, respectively, are investigated [3,4]. We are further studying electronic excitations of small molecule excitedby spatially nonuniform electric field, namely
the near-field, by using the real-time timedependent density functional theory codessuch as SALMON.

References [1] M. Kobayashi, H. Onoda,
Y. Kuroda, and T. Taketsugu, J. Comput.
Chem. Jap., 18, 251-253 (2019). [2] N.
A. Vinogradov, et al, ACS Nano, 13, 1451114518 (2019). [3] T. Tsuneda, T. Iwasa, andT. Taketsugu, J. Chem. Phys., 151, 094102(2019). [4] R. B. Jaculbia et al., Nature Nanotechnology, 15, 105-110 (2020).

375

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

Development of high-performance permanent
magnets by large-scale simulation and data-driven
approach
Takashi MIYAKE
CD-FMat, AIST
Umezono, Tsukuba, Ibaraki 305-8568

We have studied magnetic properties andstability of rare-earth magnet compounds using fist-principles calculation combined withmachine-learning techniques. The RFe12 -type
compounds having the ThMn12 structure haveattracted attention in recent years, because
high saturation magnetization can be expectedfrom their high iron concentration. NdFe12 Nx
and Sm(Fe,Co)12 films have successfully synthesized a few years ago, and it was reported experimentally that these compounds
have higher saturation magnetization, higher
anisotropy field and higher Curie temperaturethan Nd2 Fe14 B which is the main phase of theneodymium magnet. However, they are thermodynamically unstable, and partial substitution of Fe sites are necessary to stabilize a bulk
phase. The choice of rare-earth element also
aﬀects stability [1].
We calculated the magnetization, Curie
temperature and formation energy of(R1−α Zα )(Fe(1−β)(1−γ) Coβ(1−γ) Tiγ )12 (R=Nd,
Sm, Y, Z=Zr, Dy) based on density functional theory using Akai-KKR code. Nonstoichiometric composition is treated by thecoherent potential approximation. Intersite
magnetic exchange couplings are calculatedby Liechtenstein’s method, from which theCurie temperature is evaluated by solving thederived Heisenberg model in the mean-fieldapproximation.
To reduce the systematic
error in the formation energy obtained by the
KKR-CPA, data integration method (Fig.1)

is adopted, where the formation energies of
stoichiometric systems calculated by QMAS
are used to reduce the systematic error.
In order to determine the optimal composition eﬃciently, we adopt Bayesian optimization. We performed 1000 independent runs toevaluate the success rate statistically. Figure
2 shows the success rate of finding top 10 systems out of 3630 chemical compositions within
50 trials. The target variables are saturation
magnetization, Curie temperature and formation energy. The success rate is over 95 % if
we appropriately chose the descriptor, which
is much higher than that by the random sampling (12.9%).

References[1] Yosuke Harashima, Taro Fukazawa, andTakashi Miyake: Scripta Materialia 179(2020) 12.
[2] Taro Fukazawa, Yosuke Harashima,
Zhufeng Hou and Takashi Miyake: Phys.
Rev. Mater. 3 (2019) 053807.

376

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

Figure 1: Data integration method, where systematic errors contained in a large number ofdata for non-stoichiometric compositions arecorrected by a small number of more accurate
date for stoichiometric compositions.
[2].

100

Success rate [%]

80

µ 0M
TC
ΔE

60

40

20

0
RS BO #1 #2#3#4(without β nor γ)

#5

#6

#7

#8

#9

#10 #11(with β and γ)

β

γ

β, γ

(β, γ only)

Figure 2: Success rate of finding the systemwith the top 10 values of magnetization, Curie
temperature and formation energy out of 3630
candidates within 50 trials [2].

377

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

Multiscale simulations on complex multiphase flows
Youhei Morii, Takahiro Murashima, and Toshihiro Kawakatsu
Department of Physics, Graduate School of Science, Tohoku University, Sendai 980-8578
Wedeveloped

techniques(MSS)

multiscale

simulation

it (See Fig.1).

[1]

multiscale

drives the macroscopic flow is calculated using

andasimulation platform for complex fluids (MSSP)The instantaneous stress that

this microscopic simulator.

[2,3], under the support from Post-K project

We applied this technique to a viscoelastic

and HPCI. These MSS and MSSP are aimed topolymeric fluid flowing around a cylindricalapply for complex flow phenomena such as

obstacle that generates Karmann’s vortex street.

multiphase flows and viscoelastic flows whose

Figure 2 shows a comparison between (a) aconstitutive equations are not a priori known.

Newtonian fluid and (b) a Maxwellian fluid

In such systems, a simple model constitutivethat contains set of dumbbells as a model of

equation is not in general available, andpolymer solution. One can confirm the effecttherefore one must obtain the stress tensor fromof viscoelasticity coming from the polymer

the flow history of the fluid. As such a flow iscomponent on the Karman’s vortex street.

often non-Markovian, Lagrangian description isquite suitable to trace the flow history.
Smoothed particle hydrodynamics (SPH) is atypical simulation technique of fluid flows
based on Largangian description. Our MSSP is

Fig.2 Karman vortex street for (a) abased on this SPH method, where each SPH

Newtonian fluid and (b) a Maxwellian

particle contains a microscopic simulator inside

fluid.

References[1] T. Murashima, S. Urata, S. Li, Eur. Phys. J.
B, 92, 211 (2019).
[2] 森井洋平, 川勝年洋, 第 68 回高分子討論会, 福井大学 文京キャンパス, 2019/9/25.
[3] Youhei Morii, Toshihiro Kawakatsu, 10th

Fig.1 A schematic illustration of MSSP.

International

Conference

ConsortiumonofComputational

the

Asian

Materials

Science (ACCMS-10), Hong Kong, 2019/7/24.

378

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

Order-disorder transition in the superhigh-pressure phase of
Mg2SiO4
Koichiro UMEMOTO
Earth-Life Science Institute,Tokyo Institute of Technology, 2-12-1 Meguro-ku, Ookayama, Tokyo 152-8550
transitions

MgSiO3 bridgmanite with the orthorhombic

andthermodynamic

quantities

perovskite structure is a major constituent ofcalculated for these phases have been alreadythe lower mantle of the Earth. It takes the post-usedperovskite (PPV) structure near the core-mantle

dynamics of super-Earths [5].

innumerical

modelingofmantleboundary pressure (~125 GPa). The PPV phase

In this study, we predicted order-disorderis the final form of MgSiO3 in the Earth.

transition in the cation sublattice of tetragonal

However, at much higher pressures in “super-

Mg2GeO4 by first principles [4]. Mg2GeO4 is a

Earths” which are exoplanets expected to be

candidate of low-pressure analog of Mg2SiO4.

terrestrial with masses much larger than the

The dissociation of MgGeO3 PPV into pyrite-

Earth, MgSiO3 PPV should undergo “post-PPV”

type GeO2 and I-42d-type Mg2GeO4 wastransitions. So far, first principles studies havepredicted to occur at ~170 GPa [5]. Therefore,predictedof

Mg2GeO4 is a useful system to study the post-

MgSiO3 PPV: MgSiO3 PPV → I-42d-type

PPV transitions experimentally. In Mg2GeO4

Mg2SiO4 + P21/c-type MgSi2O5 → Mg2SiO4 +and Mg2SiO4, local oxygen arrangements

Fe2P-type SiO2 → CsCl-type MgO + SiO2 [1,2].

around Mg and Ge atoms are very close to each

When MgSiO3 coexists with NaCl-type MgOother. Therefore, configuration entropy isor pyrite-type SiO2, they were predicted to

expected to induce the order-disorder transition

recombine into Mg2SiO4 or MgSi2O5 [3].

at high temperature.

thethree-stage

dissociation

post-PPVThe order-disorder transition temperaturetransitions are extremely high, >~500 GPa, and(Tc) is given by the peak temperature of thehave been still difficult to be achieved

heat capacity which is calculated from theexperimentally. These predictions were based

partition function with respect to cationon the lowest-enthalpy phases under ultrahighconfigurations. For Mg2GeO4, we constructed apressures. Temperature effects by phonon weresupercell consisting of 56 atoms and generated

taken into account by using the quasi-harmonic

cation configurations until the convergence ofapproximation

Tc was achieved. Then we predicted that theTransitionpressuresof(QHA).

these

These

post-PPV379

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

order-disorder transition should occur at ~3000

Institute for Solid State Physics. This

K and 200 GPa; this pressure-temperature

research was supported by MEXT as

condition is now achievable by the diamond-“Exploratory

anvil-cell experiments. We found that Tc

computer” (Challenge of Basic Science –increases with pressure. Across the transition,

Exploring

the symmetry of Mg2GeO4 increases from

Physics and Multi-Scale Simulations)

Challenge

Extremesonthrough

Post-K

Multi-

tetragonal (I-42d) in the ordered phase to cubic(I-43d) in disordered one. In fact, the crystal
structure of the disordered cubic phase isReferences

identical to that of Th3P4. We also clarified that

[1] S. Q. Wu, et al.: J. Phys.: Condens. Matterthe effect of phonon on this transition is very

26 (2014) 035402.

small within the QHA.

[2] H. Niu et al.: Sci. Rep. 5 (2015) 18347.

Since Mg2GeO4 is the low-pressure analog[3] K. Umemoto et al.: Earth Planet. Sci. Lett.

of Mg2SiO4, the order-disorder transition

478 (2017) 40.

predicted in this study is expected to occur also[4] K. Umemoto and R. M. Wentzcovitch: inin Mg2SiO4 and will play an important role inpreparation.

modeling interiors of super-Earths.

[5] A. P. van den Berg, et al.: Icarus 317 (2019)
412.

The first principles calculations in this study
were performed using the Quantum-Espresso[6] K. Umemoto and R. M. Wentzcovitch: Phys.

packageRev. Mat. 3 (2019) 123601.

(https://www.quantum-espresso.org)

on the supercomputer (system B) at the380

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

Chain-Increment Method for Approaching the Chemical
Potential of a Polymer with All-Atom Model
Kazuo YAMADA and Nobuyuki MATUBAYASIDivision of Chemical Engineering, Graduate School of Engineering Science,
Osaka University, Toyonaka, Osaka 560-8531
We have computed the free energy of thepolymer as a sum of ∆𝜇#$%&' over the monomers

dissolution of the polymer in polymer binary

in the tagged polymer. In the present report, we

systems [1] using all-atom molecular dynamics

focus on the computation of the incremental free

simulations. An atomistic computation of freeenergy in polymer-melt systems [1].

energy is a challenging task, however, when a

The polymer species treated are polyethylene

polymer molecule is treated simultaneously as a(PE),whole. This is because a polymer is highly

methacrylate) (PMMA), and polyvinylidene

flexible and structurally diverse due to large

difluoride (PVDF) in their linear forms. All-

number of intramolecular degrees of freedom.

atom molecular dynamics simulations and free-

Given the fundamental and practical importanceenergy calculations have been implemented onof the chemical-potential calculation of polymer,

the supercomputer system at ISSP. The free

therefore, an elaborate scheme needs to beenergy of chain increment is shown in Fig. 1 at

established which exploits a specific nature of

the number of monomer units of 100. The

the polymeric structure.

incremental free energy ∆𝜇#$%&' is obtained forpolypropylene(PP),

poly(methyl

We develop a chain-increment method toeach monomer by using the scheme in [1], and it

approach the chemical potential of a polymer

is found that the incremental free energy stays

withall-atom model. The intermolecularconstant within 0.1 kcal/mol except at terminals.

interaction of a polymer of interest with the

The standard deviation is within 0.2 kcal/mol insurrounding

introduced

the inner part of the polymer, and the averaged

sequentially for the monomers, and the free

free energy of chain increment is insensitive toenergy for turning on the interaction is treatedthe chain length. All-atom computation is thus

withinenergy-feasible for the free-energy analysis of a polymer

representation theory of solutions. In our method,

system at the degree of polymerization of

we compute the free energy of chain increment

several tens or more. It is also seen that the∆𝜇#$%&' along the monomers in the tagged

preference order of ∆𝜇#$%&' agrees with that ofpolymer as well as the total free energy Dµ of thethe average interaction energy <u>i of thethemolecules

framework

is

of

the

381

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

PE

PP

PMMA

PVDF

PE

PMMA

0

∆µ (kcal/mol)∆µincr(kcal/mol)
i

-1
-2
-3
-4

-100
-50

-200
-300

-100

-5

-400
25

-6
0

<u>i (kcal/mol)

0

20

40

60

80

50

100

25

50

100

degree of polymerization N

100

-5

Fig 2: Total solvation free energy Dµ against the-10

degree of polymerization N at 25, 50, and 100 for
PE and PMMA. Dµ is estimated from ∆𝜇#$%&' at i =-15

1, 5, 10, … (N–10), (N–5), and N. The dashed lineis the least-square fit to the form of Dµ ∝ N, and

-20
0

20

40

60

80

100

the slope is −1.0 and −3.6 kcal/mol for PE andmonomer index i

Fig 1: Incremental free energy and the average

PMMA, respectively.

interaction energy against the index i of theincremented monomer in the tagged polymer. The

degree of polymerization of several tens or more.

degree of polymerization is 100 in this figure, andAtomistic computation of the chemical∆𝜇#$%&' and <u>i are plotted against the monomer

potential is of much computational demand for aindex at i = 1, 5, 10, 15, …, 90, 95, and 100. The

polymer when the interaction with the solvent ishorizontal, dashed line is the averaged value ofturned on simultaneously throughout the whole∆𝜇#$%&' or <u>i over i = 5, 10, 15, …, 90, and 95.

polymer. The method of chain increment is analternative

schemethat

introducestheinteractions sequentially along the polymer

incremented monomer with the surroundings.
We further compute the total solvation freechain. This work formulates an all-atom method

energy ∆𝜇 as a function of the degree of thefor chain increment by combining molecularpolymerization N for PE and PMMA, by using

simulation with the energy-representation theory

the incremental free energy of solvation ∆𝜇#$%&' atof solutions, and it is demonstrated that the freei = 1, 5, 10, … (N–10), (N–5), and N. In Fig. 2,

energy of chain increment can be obtained at aa good linearity of Dµ to N is observed. The free

precision of a few tenths of kcal/mol.

energetics of a long polymer is thus determined

Referencesby the chain-increment contribution ∆𝜇#$%&' , and[1] K. Yamada and N. Matubayasi: Macromol.

Figs. 1 and 2 demonstrate that ∆𝜇#$%&' can be

53 (2020) 775-788.

assessed from an all-atom simulation with the382Activity Report 2019 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Study on thermodynamic properties and NOx redox
reactions of ternary alloy nanoclusters by O(N)density functional calculations
Shih-Hsuan Hung and Taisuke Ozaki
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

We perform first-principles calculations andMonte Carlo sampling to investigate the structures of ternary PdRuM (M=Pt, Rh, or Ir)nanoparticles (NPs) with respect to three different spherical shapes. The calculation is carried out using OpenMX (Open source package
for Material eXplorer) software package.

1

Introduction

Metallic nanoparticles (NPs) are a promising
functional material in many aspects, such asenergy conservation, catalyst, and elements
storage. [1] To achieve better performance ofparticular functions, many experiments synthesize modified NPs to adjust their properties
physically and chemically. Bimetallic PdRu
NPs have been reported for providing better
catalytic reactions, but the drawback is suffered from the loss of metal element duringcatalytic cycles. [1, 2] Therefore, PdRuPt,
PdRuRh, and PdRuIr NPs or alloy have beenstudied and demonstrate their promising ability of NO reaction. [4, 3] Therefore, we would
like to investigation these NPs and provide
deeper insight in atomistic aspect.

2

ResultTo compare different morphology influencing different composition ternary NP, we select three different spherical shapes for theNPs. There are hexagonal close-packed (hcp),

truncated-octahedral (fcc), and icosahedral(fcc) shapes with 57, 55, and 55 atoms, respectively. To achieve higher coverage of energy survey, each nanoparticle is sampled 30
times using Monte Carlo sampling method,
thus there is 270 sampling in total. The statistical characteristic implies the position trendof individual atom for respective nanoparticles (detail is in published paper). By following the trend, we manually construct the bareternary NPs and find the most stable one. Figure 1 shows the most stable ternary NPs in
our calculation. The calculation shows thatthe atomic arrangements depend on the surface formation energy of individual elements.
In other words, the element possesses smaller
surface formation energy dominating the surface sites. The species have greater surface formation energy immerse in the NPs as bulk-likestate. Figure 1 demonstrates Ru atoms are always as a core in NPs due to its larger surface formation energy (γ = 1.31 eV/atom) offcc (111) surface. In addition, there are twoparticular shells are found in the investigation.
One configuration is binary solid solution alloy
and the other one represents that two elements
segregate from each other on the shell. The
reason of these two particular arrangements is
because the NPs tend to minimize its surface
strain, compared to lattice constant of Pd fcc
bulk (3.79 Å). For the NPs with solid solution shell, the Pd-M bond length decreases by
3.8 % in average. Meanwhile, the bond reduc-383

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

Figure 1: The most energetically stable
ternary PdRuM NPs in the calculations. The
unit of Ef is eV/atom.
tion of 3.0 % in average has been found for theconfigurations with segregation shell. In addition, the Ir-Ir bond length shows a significant
reduction of 6.2 %. By modeling the single
layer hexagonal as shell structure, the calculation shows a good agreement with our NPs
calculation. The reduced bond length of alloy
decreases the total calculated energy.
Oxidized NPs is an other problematic issueof catalyst. Due to the oxidation formationenergy variance, the elements of the NPs may
migrate from site to site. The oxygen adsorption on surface calculation demonstrates thatthe preference of oxygen adsorption (eV/atom)
on fcc hollow site is found to be Ru (-3.65)→ Rh (-2.78) → Ir (-2.47) → Pd (-2.06) →Pt (-1.84). Besides the single oxygen adsorption on cleaved surface, we examine the oxygen
molecule adsorption on pure Pd NPs. The calculation demonstrate taht the oxidation energy
is as a function of O-O distances. The increasing O-O distance up tp 5.0 Å decreases theoxidation formation energy. Figure 2 demonstrates the most stable oxidized ternary NPs
and their oxidation formation energy what we
can find. The result shows that the Ru atoms
migrate outward to the surfaces of NPs to beoxide. Pt, Pd and Pd atoms penetrate inwards
as a core of the ternary PdRuPt, PdRuRh, and
PdRuIr NPs, respectively. By calculating theNPs formation energy with respect to chemical
potential of oxygen molecule, the free energy
indicates hcp ternary NP is thermodynami-

Figure 2: The most energetically stable structures of high oxygen coverage on hcp, tOh, andIh ternary NPs.
cally unstable under oxygen-rich condition.

3

Conclusion and Discussion

In summary, we have provide deeper insight
into the most stable structures with three different morphologies and two oxidation statuses. The investigation shows the atoms areorganized in specific arrangement instead ofsolid solution alloy. In addition, the oxidation
condition can facilitate the atomic migrationin the ternary NPs. Although the scale of thestudy is not satisfying to interpret real situation, it can be a pilot investigation for real
scale NPs calculation and provide some possible configurations in advance.
References[1] Kusada Kohei et al.: J. Am. Chem. Soc.
136 (2014) 1864–1871[2] Wu Dongshuang et al.: Phys. Chem.
Chem. Phys. 14 (2012) 8051–8057[3] Sarker Md Samiul et al.: J. Mater. Res.
29 (2014) 856–864[4] Shang Changshuai et al.: Nano Res. 11
2018 4348–4355

384