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

3.5

SCCMS Projects352

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

Developing graph neural networks that helps
predicting the molecular dynamics of glassy liquids
Hayato SHIBA
Information Technology Center, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-0882

In the past two decades, structure–dynamics
correspondence in glass-forming liquids havebeen hot topic. Along with the drastic slowdown of atomic motion accompanying the glasstransition, localized domains of les that rearrange more preferentially than others growup. Structural analyses of local geometric orders and of real-space normal modes of vibrations have clariﬁed the existence of correlations between the structure and the dynamics. It has remained, however, unclear to what
extent the dynamics can be dictated by the
structure. Recently, machine learning methodologies, in particular, graph neural networks(GNNs), enabled us to predict the slow dynamics far better than these methodologies. The
GNNs have turned out to be a good ﬁt formaking predictions more accurately than ever,and thus, open a new paths to further reveal
the structure-dynamics relationship.
In FY2022, we have engaged in improving
our model, BOTAN (BOnd TArgeting Network), which we have proposed in FY2021.
In BOTAN, a decoder for edge features is additionally implemented to the previous GNN
model to enable itself to learn relative motion
between neighbor pairs (edge target) in addition to self-motion of particles (node target).
This model has tens of thousands of weightparameters in it, which is the reason why it
can capture the complex relation between the
structure and the dynamics. However, in themeawhile, it is diﬃcult to train such a large
model, especially when edge features come in

0.03ταBapst GNN
FY2021
FY2022

T = 0.44
10-1

100

τα

101

102

t

103

104

105

Figure 1: Comparison of pearson correlationcoeﬃcients between diﬀerent GNN models:the original GNN model and of BOTAN (inFY2021 and FY2022).
as the training target. In order to ﬁx such
shortcomings, we have introduced a set of pretrained weight parameters that reﬂect the factthat closer pairs are more likely to separate
in a short-time range. This improvement has
updated the best predictive performance evenon the time-scale of glassy structural relaxation, see Fig. 1 in which the Pearson correlation coeﬃcients representing the predictive
accuracy (τα denotes the structural relaxationtime). BOTAN is now one of the state-ofthe-art models for glassy dynamics relaxation,
as well as other models that have also appeared in FY2022. The full article introducing his state-of-the-art model of BOTAN has
been published [1], together with training code
implemented on PyTorch Geometric [2].
In FY2022, we have started our next trial to353

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

develop a new application of BOTAN: learning
the dynamics of a glassy liquid under quasistatic shear at zero temperature. We havesimulated using so-called the Kob-Andersen binary glass-forming mixture and computed thedataset for training by using System B at ISSP,
and now we are trying to apply BOTAN for
the prediction of rearranging hotspot from theinitial structure.

References[1] H. Shiba, M. Hanai, T. Suzumura, andT. Shimokawabe: J. Chem. Phys. 158,
084503 (2023).
[2] https://github.com/h3-Open-BDEC/
pyg botan

354

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

A systematic ab initio study of quasi-one-dimensional
molecular conductors TM salts
Kazuyoshi Yoshimia , Takahiro Misawaa,b , Takao Tsumurayac , and Hitoshi Seod,eaInstitute for Solid State Physics, University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
b
Beijing Academy of Quantum Information Sciences, Haidian District, Beijing 100193, Chinac
POIE, Kumamoto University, 2-39-1 Kurokami, Kumamoto 860-8555, Japand
Condensed Matter Theory Laboratory, RIKEN, Wako, Saitama 351-0198, Japane
Center for Emergent Matter Science (CEMS), RIKEN, Wako, Saitama 351-0198, JapanThe quasi-one-dimensional organic conductor TM2 X (TM=TMTSF, TMTTF, X: monovalent anion) has been extensively studied experimentally and theoretically due to its diverse ground states induced by changes in pressure and molecular substitution[1, 2]. Whilethe importance of electronic correlation effects
has been pointed out in these studies, a systematic evaluation of electronic correlation effects has not been carried out due to the difficulty of quantitatively estimating electronicinteractions in molecular systems.
In this study, we derive effective models
using first-principles calculations as a starting point and the constrained random phaseapproximation method. The analysis of organic conductors involves molecules with many
atoms as a unit, which requires significant
computational effort and memory to derive
effective models. We addressed this problem by using the recently developed software
RESPACK[3] and exploiting its parallelization
feature to take advantage of large computing resources. Using the first-principles software Quantum ESPRESSO[4] and RESPACK,
we systematically derived effective models for
TM2 X and analyzed the electronic states using a many-variable variational Monte Carlosolver, mVMC [5].
Based on the obtained results, we discuss
the material dependence of the electronic cor-

relation effects in the TM molecular systemand clarify the relationship between the different ground states observed in the experimentand the electronic correlation effects. We showthat the correlation effect in this system can becharacterized by the screened on-site Coulombinteraction U relative to the intrachain transfer integrals. This parameter is material dependent and we show that it plays a key rolein understanding the spin and charge ordering
phenomena in the system[6].
We thank R. Kato, S. Kitou, and T. Nakamura for fruitful discussions.
This work
was supported by MEXT/JSPJ KAKENHI
under grant numbers 19K03723, 19K21860,
20H04463, 21H01041 and 22K03526.
KY
and TM were supported by Building of Consortia for the Development of Human Resources in Science and Technology, MEXT,
Japan.
The input and output files for
RESPACK and mVMC are uploaded to ISSP
data repository[7].

References[1] D. Jerome, Chem. Rev. 104, 5565 (2004).
[2] K. Yoshimi, H. Seo, S. Ishibashi, and S.
E. Brown, Phys. Rev. Lett. 108, 096402(2012), and references therein.

355

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[3] K. Nakamura, Y. Yoshimoto, Y. Nomura
et al., Comp. Phys. Commun. 261, 100781(2021).
[4] P. Giannozzi, O. Andreussi, T. Brumme
et al., J. Phys.: Condens. Matter 29,
465901 (2017).
[5] T. Misawa, S. Morita, K. Yoshimi et al.,Comp. Phys. Commun. 235, 447 (2019).
[6] K. Yoshimi, T. Misawa, T. Tsumuraya,
and H. Seo, arXiv:2210.13726.
[7] https://isspns-gitlab.issp.utokyo.ac.jp/k-yoshimi/tm-salts

356

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

Development of an effective model estimation tool by
Bayesian optimization
Kazuyoshi Yoshimi
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

Target vs Estimated magnetization curve
0.5

TargetEstimated (BO_NO13)

0.4
0.3
Sz

Analysis based on effective models is useful to understand properties of materials. The
reliability of an effective model is often evaluated by comparing it with experimental data,
but this requires trial and error due to the existence of the large number of model parameters.
We have developed a tool called Bayesian optimization tool of Effective Models (BEEMs)[1] to estimate effective models as an inverse
problem using experimental data and Bayesian
optimization with Gaussian processes [2].
In BEEMs, the quantum lattice model solverHΦ[3] is used as the forward solver and themagnetization curves are calculated based on
the given Hamiltonian. The cost function isthe difference between the obtained magnetization curve and the target magnetization
curve. We use the Bayesian optimization library PHYSBO[4] to find the next candidate
point of the Hamiltonian that may have asmaller cost function.
BEEMs is available from the GitHub pageand can be installed by the following procedures:

0.2
0.1
0.0
0.0

0.5

1.0

1.5
H

2.0

2.5

3.0

Figure 1:
The magnetization curve at(J2 , J3 ) = (0.1, 0.02) (represented by red symbols).The proposed magnetization curve using
BEEMs (represented by the blue line).

is a tool for plotting magnetization at specified
candidate points. output_score is a tool for
displaying the score and the candidate point
parameters up to the specified rank.

As a demonstration, we performed model
parameter fitting for S = 1/2 Heisenberg spin
chain with 12 sites under periodic boundary
condition. The Hamiltonian is given by H =P12i=1 [J1 Si · Si+1 + J2 Si · Si+2 + J3 Si · Si+3 ].
$ git clone
For simplicity, we set J1 = 1 and calculated
https://github.com/k-yoshimi/BEEMs.git
the magnetization curve at J2 = 0.1 and$ cd BEEMs
J3 = 0.02 as the reference data as shown in$ pip3 install .
Fig. 1. Figure 2 shows the mean squared error
Then, you can use beems as a python pack- (MSE) between the reference data and theage. Post tools gauss_fit, mag_plt, and calculated data. Up to the fifth step, ranoutput_score are also available. gauss_fit dom sampling is performed even if Bayesianis a tool for constructing a score function us- optimization is selected. It is seen that aftering a Gaussian process based on the results performing Bayesian optimization, the cost
obtained by Bayesian optimization. mag_plot function reaches almost zero in the 13th step,357

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

0.008

BO (from 5 steps)

0.5

Random sampling

0.4

0.006
0.004
0.002
0

MSE

0.3

J3

Mean Squared Error (MSE)

0.01

0.2
0.1

0

5

0

10 15 20 25 30 35 40 45 50
Steps

Figure 2: Mean squared error(MSE) between
the reference data and the calculated dataobtained by the Bayesian optimization (BO)and random sampling, respectively. The green
shaded area indicates the region where the random sampling is performed.
while random sampling is underperformed
comparing to Bayesian optimization evenafter 50 steps. The optimal parameters foundby Bayesian optimization for the best score
are J1 = 1.0, J2 = 0.1, and J3 = 0.025. We
created a trained model using gauss_fit.py
to estimate MSE obtained with the Gaussian
process after 50 trials. As a result, we foundthat the region where (J2 , J3 ) is relatively
small has a lower score as shown in Fig. 3.
Using such a map, it is possible to check theaccuracy of the estimated parameters. For
example, we could explore approaches such as
refining the grid size in the low score region.
In summary, we have developed BEEMs: aneffective model estimation tool by combiningHΦ and PHYSBO. In version 1.0, it is possible
to estimate the effective model based on themagnetization curve. In the future, we aim to
improve the calculation speed, as currently all
magnetic fields have to be calculated. We plan
to add a mode that can calculate the magnetization curve by splitting the calculation by Sz
when Sz is conserved. We also plan to extend
the tool to estimate other physical quantities
if there is demand. We also plan to provide anexpert level manual and allow users to define
their own classes for use.

0

0.2 0.4 0.6 0.8

J2

1

0.018
0.016
0.014
0.012
0.01
0.008
0.006
0.004
0.002
0

Figure 3: The color map of MSE estimated by
the trained model using gaussian process.
We would like to acknowledge T. Misawa,
R. Tamura, and Y. Yamaji for their contribution to code development and useful discussion. KY was supported by JSPS KAKENHI
Grant Number 21H01041.

References[1] https://github.com/k-yoshimi/BEEMs[2] C.E. Rasmussen and C.K.I. Williams,
Gaussian Processes for Machine Learning,
The MIT Press (2005).
[3] M. Kawamura, K. Yoshimi, T. Misawa
et al., Comp. Phys. Commun. 217, 180(2017).
[4] Y. Motoyama, R. Tamura, K. Yoshimi
et al., Comp. Phys. Commun. 235, 447(2019).

358

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

Exciton Properties in Organic Optoelectronic Devices
from Large-Scale Electronic Structure Calculations
Takatoshi Fujita
National Institutes for Quantum Science and Technology
Anagawa, Inage-ku, Chiba 263-8555

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.
We have developed the novel large-scale excited-state methods [1] based on the fragment molecular orbital(FMO) method, exicton model, and GW and
Bethe-Salpeter equation (BSE) method with
Tamm-Dancoﬀ approximation (TDA). In thismethod, the excited-state Hamiltonian matrix
elements at the GW/TDA-BSE level is calculated in the basis of fragment conﬁguration
state functions which describe intrafragment
excitations or interfragment charge-transfer
excitations. The excited-state Hamiltonian isthen diagonalized to approximate the adiabatic excited states of an entire system. The
novel fragment-based GW/BSE method allows for the eﬃcient calculation of delocalized excited-states in realistic molecular aggregates. In addition, the excited-state Hamiltonian derived this method can be used forreal-time simulations in combinations with the
wavepacket prorogation method.

model Hamiltonian for single-excitation anddouble-excitation manifolds was derived on thebasis of fragment-based excited-state calculations within the GW approximation and theBSE. After elucidating the energetics of theelectron–hole separation and examining linear absorption spectrum, we investigated thequantum dynamics of exciton and charge carriers in comparison with the pump–probe transient absorption spectra. In particular, we
introduced the pump–probe excited-state absorption (ESA) anisotropy as a spectroscopic
signature of charge carrier dynamics after exciton dissociation. We found that the charge
separation dynamics can be probed by thepump–probe ESA anisotropy dynamics after
charge-transfer excitations. The present study
provides the fundamental information for understanding the experimental spectroscopy signals, by elucidating the relationship between
the excited states, the exciton and charge
carrier dynamics, and the time-resolved spectroscopy.

We have developed a computational method
for excited-state dynamics and time-resolved
spectroscopy signals in molecular aggregates [2], on the basis of ab initio excitedstate calculations. As an application, we
consider the charge separation dynamics andpump–probe spectroscopy in the amorphous
P3HT/PCBM blend. To simulate quantum
dynamics and time-resolved spectroscopy, the[1] T. Fujita, Y. Noguchi, J. Phys. Chem. A
125 (2021) 10580.

References[2] T. Fujita, T. Hoshi,
https://doi.org/10.26434/chemrxiv-2023773xs.

359

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

Development of an efficient method for calculating
miscibility of long-chain polymer blends using molecular
dynamics simulations
Kazuo Yamada and Nobuyuki Matubayasi
Division of Chemical Engineering, Graduate School of Engineering Science,
Osaka University, Toyonaka, Osaka 560-8531
We have studied the mutual miscibility ofof a theory of solutions (energy-representation

polymer species in polymer mixtures using all-method). In our method, we compute theatom molecular dynamics simulations. The

monomer-wise chemical potential ∆𝜇!"#$% formutual miscibility of polymers is determined

the incremented monomers in the taggedby the free energy of mixing given by thepolymer and the total free-energy of solvationdifference between the chemical potentials of∆𝜇 of the polymer is the sum of ∆𝜇!"#$% . This

the pure and mixed states of the constituent

year, we focused on the computation of the freespecies. Since the chemical potential isenergy of solvation in polymer mixture systems.

sensitive

tointra-

andintermolecularinteractions, an atomistic treatment of polymer

Results

mixture systems is desired to evaluate free

The free energy was computed for the mixing

energyatomisticof polyethylene (PE) and polyvinylidene

calculation of the chemical potential of an

difluoride (PVDF) mixture systems with theentire polymer chain is a difficult task withmolar ratio 1:0, 3:1, 1:1, 1:3, and 0:1. The

current computer capabilities because of thenumber of the monomer unit n and polymer

enormous number of intramolecular degrees of𝑛&'( are fixed 100 and 200, respectively. All-

freedom in polymer chain. To solve thisatom molecular dynamics simulation and free-problem, we developed a chain increment

energy calculation has been computed bymethod that approaches the chemical potential

supercomputer system at ISSP. Figure 1 showsof polymers using an all-atom model [1]. The

the monomer-wise chemical potential ∆𝜇!"#$%

intermolecular interaction of a polymer ofand averaged interaction energy 〈𝑢〉! for eachinterest with the surrounding molecules ismonomer i = 20, …, and 80 by using our

introduced sequentially for the monomers

scheme [1]. It was seen that both ∆𝜇!"#$% and(segments), and the free energy for turning on〈𝑢〉!

the interaction is treated within the frameworkmonomer in the inner part of the polymer chainof

mixing.

However,360

are insensitive to the location of theActivity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyoand depends on the mixing ratio of the polymer

polymeric systems, atomistic calculations ofsystems. The standard error of ∆𝜇!"#$% at i eachsolvationof i = 20, ... , and 80 is less than ∼0.1 kcal/mol

computationallyfor both PE and PVDF at all the molar ratios inincrement

Figure 1.

approaches free energy of mixing by focusingPE∆µincr(kcal/mol)
i

1:0

1:1

1:3

0:1

are

expensive.

methodweextremely
The

have

chain

developed

monomers turn on their interactions with

−0.5
−1.0

surrounding molecules. Combined with energy-

−1.5
−2.0

representation method, the pair-wise solvation

−2.5
−3.0

<u>i (kcal/mol)energieson the solvation free energy as repeated

PVDF

3:1

free

free energy can be calculated in computational

−5.0
−6.0
−7.0

time

−8.0
−9.0

that

simulations.

isfeasiblefor

all-atomMD

Extendingthisapproachto20 40 60 80 20 40 60 80 20 40 60 80 20 40 60 80 20 40 60 80

polymer mixtures and applying it to themonomer index i

Figure 1: Incremental free energy ∆𝝁𝐢𝐧𝐜𝐫
and the average𝒊

PE/PVDF system demonstrates that polymer

interaction energy 〈𝒖〉𝒊 against the indices i of the incremented

mutual miscibility evaluations can be computed

PE and PVDF monomers in the PE/PVDF mixture systems.

in an all-atom model.

Figure 2 shows the mixing free energy per

0.3∆mixG/N (kcal/mol)monomer ∆&") 𝐺/𝑁 of the PE/PVDF mixture
versus the molar fraction of PVDF. According
to Figure 1, ∆𝜇!"#$% of PE is within 0.1 kcal/mol
over the entire range of mixing ratios, and themole fraction dependence of

0.2

0.1

∆&") 𝐺/𝑁 is0dominated by that of PVDF. Also, the sign of0

1/4∆&") 𝐺/𝑁 is positive for each mixing ratio
examined since ∆𝜇!"#$%

1/2

3/4

1

mole fraction of PVDF

of PVDF becomes

Figure 2: Free energy of mixing per monomer ∆𝐦𝐢𝐱 𝑮/𝑵 for

negative for larger mole fractions of PVDF.

the PE/PVDF system at mole fractions of 0, 1/4, 1/2, 3/4, andFigure 2 thus shows that PE and PVDF are

1 for PVDF. Since the total number of polymer chains in theimmiscible with each other, which is consistentsystem is 200 and the degree of polymerization n is 100 for PE

with experimental reports.

and PVDF, the total number of monomers N is equal to 20000.

References
Conclusion[1] K. Yamada and N. Matubayasi: Macromol.

When solute degrees of freedom are large, as in

53 (2020) 775-788.

361

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

Fugaku Battery and Fuel Cell Project
Osamu Sugino
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
chose

Fiscal 2022 was the final year of thethosecontaining

relatively

small

Fugaku-denchi Project [1], in which weendothermic elementary reactions only. The

participated in the study of the oxygen

hydration of OH* was suggested to be the rate

reduction reaction (ORR) occurring on an oxide

determining step under the assumption of minor

electrocatalyst, i.e., zirconium oxynitride. This

contribution of the details of activation barrier.

defective zirconia is known to produce high

The corresponding reaction heat is about 0.4 eV,

voltage comparable to platinum under open

which is comparable to that on the platinum

circuit conditions, but the reaction mechanism

surface. Importantly, a similar value wasis still unknown. Due to the lack of detailed

obtained for the surface without defects,

experimental information, we have performed asuggesting that the defects play a minor role inthought experiment based on reliable first-

determining

principles simulations. We first performed a

Experimentally, however, the defects werefirst-principles Monte Carlo simulation with

reported to enhance the reaction current density.

abICS [2] to find a distribution of nitrogen

A possible explanation for the apparent

impurities (NO) and oxygen vacancies (VO).

discrepancy is that defects contribute to

The results were used to create a model where

increasing the carrier density without changing

NOs and VOs are located at a subsurface andthe activity as much. We are now examining

surface site, respectively. We then investigatedother possibilities as well by using advanced

the free-energy profile for ORR occurring on

exchange-correlation functionals, increasing thethe model [3].

number of possible structures, and consideringtheactivityof

ORR.

Since there are many possible reaction sites

the activation energies. Making a reliable

for oxide surfaces, we examined all possible

machine learning potential will be an important

adsorption structures of the ORR intermediates

step for that purpose.

and then used a quick and easy method to

We are grateful for the opportunity

determine the reaction pathways based on theprovided by the ISSP supercomputer systems to

structural similarity of the intermediates.

study the ORR intermediates in detail, which

Among the resulting free energy profiles, wewould otherwise be too challenging.

362

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

References

Vacancy on the Oxygen Reduction Reaction on[1] https://www.nims.go.jp/fugaku-

the Tetragonal Zirconia (101) Surface”, S.

denchi/en/index.html

Muhammady, J. Haruyama, S. Kasamatsu, and[2] https://www.pasums.issp.u-

O. Sugino, J. Phys. Chem. C126, 15662–15670

tokyo.ac.jp/abics/(2022).

[3] “Effect of Nitrogen Doping and Oxygen363

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

Materials design of high performance magnetic materials
Tetsuya FUKUSHIMA
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

In this year, we have developed an all electron first-principle calculation code, where afull-potential Korringa-Kohn-Rostoker (KKR)
Green’s function method is implemented, forthe design of permanent magnets and spintronics materials.
It is well known that the local density
approximation (LDA) and generalized gradient method (GGA) are poor approaches forstrongly correlated materials. In particular,
LDA and GGA give unrealistic results for localized f -electron systems. When the electron
configurations are not f 7 or f 14 , the occupied
and unoccupied orbitals split off due to thecorrelation energy. However, in the LDA andGGA calculations, the f orbitals do not split
and are pinned on the Fermi energy, leading
to significantly overestimated cohesive energy.
One of the origin of this problem is the error for self-interaction energy. The Coulomb
and exchange terms in a same orbital should
be canceled out, because the energy of electroninteracts with itself (self-interaction energy) iszero. This is satisfied in Hartree-Fokc approximation; however, LDA and GGA lead to finite
self-interaction energy.
We implemented the self-interaction correction (SIC) method to our full-potenital KKR
Green’ function code and applied it rare earth
compounds. In our code, the spherical electron densities of the 4f orbitals are extracted
from the Green’s function, and then the spherical components of the SIC potential are calculated. This means that the our SIC potential
is not full-potential. The full-potential SIC isneeded to accurately describe the orbital split-

(a)

(b)

(c)

(d)

-1.5

Figure 1: Calculated density of states of
NdCo5 and SmCo5 . (a,c) and (b,d) correspondto LDA and SIC results, respectively.

tings reflecting from the geometric shape of the
f orbitals. The SIC potential is constructedby using the exchange correlation potential parameterized by Moruzzi, Janak, and Williams.
In order to investigate the effect of SIC on
the electronic structure, we compared the density of states of NaCo5 and SmCo5 calculated
by LDA and SIC (see figure 1). We used theexperimental lattice constants in these calculations. The space group is P 6/mmm. As
shown in Fig. 1(a,c), in the LDA calculation,
the partially occupied f states are located on
the Fermi level. On the other hand, SIC lead tothe splitting unoccupied and occupied f states.
Note that the LDA+U method is often used tocorrect for such LDA errors, but the value ofon-site Coulomb interaction (U ) is an empirical parameter and is not suitable for materials
screening and design.
Figure 2 shows the density of states of Ce2 O3364

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

Density of states (1/eV)

30
25

KKR-LDA
Ce-4f

20
15
10

O-2p

5
0
30

Density of states (1/eV)

235101 (2007)

Total

25

-6
-4
-2
0
2
4
Energy relative to the FermiTotal
level (eV)KKR-SIC

20

Ce-4f

15
10

O-2p

Ce-4f

5
0

-6
-4
-2
0
2
4
Energy relative to the Fermi level (eV)Figure 2: Density of states of Ce2 O3 calculated
by LDA and SIC.
calculated by LDA and SIC. The space groupis P¯3m1 and the antiferromagnetic state is
assumed in the calculations. It is thoughtthat the valence of the Ce ion in Ce2 O3 is
3+ and the electron configuration f 1 is realized. In the SIC calculation, the occupied fstates containing one electron are pushed down
to energetically lower region, and are locatedaround -3.0 eV. The unoccupied f states lie
around 1.0 eV to 3.0 eV. Pourovskii et al. performed the electronic structure calculations for
Ce2 O3 on the basis of Linear Muffin-Tin Orbital Method with dynamical mean field approximation (DMFT) [1]. Our results obtained
from the SIC calculations are quite consistent
with their DMFT calculations.

References[1] L. V. Pourovskii, B. Amadon, S. Biermann, and A. Georges, Phys. Rev. B 76,

365

Activity Report 2022 / 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, National Institute of Advanced Industrial Science and Technology
Umezono, Tsukuba, Ibaraki 305-8568

The calculation is based on density functional theory in the local density approximation. We use AkaiKKR which is based on
the Korringa-Kohn-Rostoker Green function
method. The open-core approximation is used
for the f electrons in Nd and Sm. Randomness caused by doping is treated by coherent
potential approximation.
In order to investigate the stability, we
consider R1−ζ1 −ζ2 −ζ3 −ζ4 Zrζ1 Feζ2 Coζ3 Tiζ4 . In
the four-dimensional space (ζ1 , ζ2 , ζ3 , ζ4 ), we
compare the total energy for the 1-12, 217, and unary phases, and identified stablephases. Figure 1 shows the phase diagramof the Sm system on the (Sm,Zr)(Fe,Co)12
plane. We then compare the total energy difference between (R1−α Zrα )(Fe1−β Coβ )12−γ
Tiγ and stable phases, and calculate thehull distance. Figure 2 shows the result for(Sm1−α Zrα )(Fe1−β Coβ )12 . We see that thehull distance is not monotonic as a functionof α and β. We also see that the hull distanceis small for high Zr concentration.
1.0 (1)0.8
Zr concentration

RFe12 -type compounds with the ThMn12
structure are potential main phase for rareearth magnets. Because of high Fe content,high magnetization and high Curie temperature can be expected. However, RFe12 is
known to be unstable, and there have been
many trials for stabilizing the ThMn12 structure by chemical doping. In this study, weinvestigate the stability of (R,Zr)(Fe,Co,Ti)12(R=Y, Nd, Sm) against the 2-17 and unaryphases by first-principles calculation [1].

(4)

(5)

(6)

0.6
0.4
0.2
0.0
0.0

(2)

(7)

(8)

(3)

0.2

0.4
0.6
Co concentration

0.8

1.0

Figure 1:
Phase diagram on the(Sm,Zr)(Fe,Co)12 plane.
(1) (S,Zr)Fe12 ,Zr(Fe,Co)12 , Zr2 Fe17 and unary Fe; (2)(Sm,Zr)Fe12 , Sm2 (Fe,Co)17 , and unary Sm
and Fe; (3) (Sm,Zr)Fe12 , Sm2 (Fe,Co)17 , andunary Fe; (4) (Sm,Zr)Fe12 , Zr(Fe,Co)12 ,Sm2 (Fe,Co)17 , and unary Fe; (5) Zr(Fe,Co)12 ,Sm2 (Fe,Co)17 , and unary Fe; (6) Zr(Fe,Co)12 ,Sm2 (Fe,Co)17 , Zr2 Co17 , and unary Fe; (7)Sm2 (Fe,Co)17 , Zr2 Co17 , and unary Fe; (8)Sm2 Co17 , Zr2 Co17 , and unary Fe and Co.

References[1] Taro Fukazawa, Yosuke Harashima,
Hisazumi Akai and Takashi Miyake:Phys. Rev. Mater. 6, 054404 (2022).

366

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

0.030

1.0
Zr concentration

0.025

0.010

0.8

0.020

0.6

0.015
0.4

0.0
3

0

0.005

40

0.0
0.0

0.010

20

0.0

0.2

0.0

0.2

0.4
0.6
0.8
Co concentration

1.0

0.000

Figure 2:
Hull distance in eV for(Sm1−α Zrα )(Fe1−β Coβ )12 with the ThMn12
structure.

367

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

Data driven analysis for impurity eﬀects on
photocatalysts
Mikiya FUJII and Yosuke HARASHIMADivision of Materials Science, Nara Institute of Science and Technology
Takayama-cho, Ikoma, Nara 630-0192

Photocatalytic generation of hydrogen fromsolar energy (green hydrogen) is one of the key
elements for solving resource, energy, and environmental problems. The purpose of thisstudy is to obtain guidelines for controlling theproperties of new photocatalytic compounds
and impurities to improve the yield of hydrogen production in response to visible light,
which is still not suﬃcient. For this purpose, a combined analysis of ﬁrst-principles
calculations and informatics will be performed.
While materials development has been conducted within a limited search space basedon experience, this study will develop a newresearch in photocatalytic materials development by utilizing information science to predict and develop new materials with a greatly
expanded search space. [1, 2] In addition tosolving energy problems through hydrogenproducing photocatalysts, this research will
have a spillover eﬀect on the synthesis ofcarbon-neutral petrochemical feedstocks from
CO2 and green hydrogen.
Process informatics for calcination processwith dynamical Monte Carlo simulations:
The hydrogen evolution rate (HER) of photocatalytic materials depends strongly on theircalcination processes. The calcination process will determine an amount of product compound, defect, and surface area. The HER
is considered to be highly dependent on thosefeatures, however, experimental analysis to detect the features is time-consuming. We developed a dynamical Monte Carlo (DMC) simu-

lation code to predict the features, and constructed an HER prediction model using the
calculated features. The DMC code is written in C++ and the HER prediction code isin Python. The DMC code is based on formation energies of product compounds, defectenergies, and surface energies. We will use
the ﬁrst-principles calculation code, Vienna Ab
initio Simulation Package (VASP), to calculatethese quantities. In Ref. [3], quantum chemical
calculations were used as explanatory variablesin the model to predict experimental values in
ﬂow copolymerization [4]. It was found that
information on the so-called transition statesof intermediates is important in predicting the
products of chemical reactions. This is a kindof data assimilation [5], which we will incorporate into our future method of calculatingintermediates in crystal formation by processsimulation.
Highthroughput search of high entropy catalyst:
Platinum is known for its high catalytic performance, but it is very rare, so there is aneed to develop alternatives. We performed asystematic search for high-entropy alloys thatwould mimic platinum. The electronic structures were calculated by using ﬁrst-principles
code, AkaiKKR, which is based on the KKRGreen’s function method. It enables us toconsider nonstoichiometric chemical compositions by the coherent potential approximation. We considered up to a quaternary system
which consist of transition metal elements andperformed ﬁrst-principles calculations on more

368

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

than 100,000 samples.

References[1] M. Fujii: JSAP Review, 2022, 220416(2022).
[2] 藤井幹也: 応用物理, 91, 688 (2022).
[3] S. Takasuka, et al.: Digital Discovery, accepted.
[4] A. Wakiuchi, et al.: Macromol. Mater.
Eng. 2200626 (2022).
[5] Y. Harashima: (invited talk) The 24th
Annual Meeting of Japan Society of Theoretical Chemistry, May 18th 2023.

369

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

Computational and Data Materials Science Study for ET
Revolution by Developing Next-Generation Battery and
Fuel Cell
Yoshitaka TATEYAMANational Institute for Materials Science,
Namiki 1-1, Tsukuba, Ibaraki 305-0044
We clarified a long-standing problem of

We further addressed the speed-up of theatomistic insight of the dopant effect on the ionic

conductivity calculations by developing a new

conduction around grain boundary, by extensive

method of chemical charge diffusion non-first-principles MD calculations of diffusion

equilibriumcoefficients[1].

demonstrated that this new method surely

MD(CCD-NEMD),

andaccelerates the sampling of mean-squaredisplacement[3].

Fig. 1: Calculated self-diffusion coefficients of
Li7La3Zr2O12 solid electrolyte.
Fig. 3: (left) Definition of flux and color
charges to the Li ions and anion units in the
Li10GeP2S12 model system. (right) Schematic
picture about the speed-up of the evaluation ofcorrelated ionic conductivity.

We also developed a high-throughput firstprinciples MD sampling flow for the ionic
conductivity

calculations,andproposedpotential solid electrolytes with higher ionic

References

conducitivy by applying the flow to Na-ion[1] B. Gao, R. Jalem, and Y. Tateyama, J. Mater.

sulfide solid electrolytes[2].

Chem. A, 10, 10083-10091 (2022).
[2] S.-H. Jang, Y. Tateyama, and R. Jalem, Adv.
Funct. Mater. 32, 2206036 (2022).
[3] R. Sasaki, B. Gao, T. Hitosugi, and Y.

Fig. 2: Schematic picture of the highthroughput screening flow for the search forsolid electrolyte materials with the higher ionic
conductivities.

370

Tateyama, npj Comput. Mater. 9, 48 (2023).

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

Theoretical study of hydrogen-evolution
semiconductor photocatalysts using first-principles
calculations
Seiichiro L. TEN-NO
Graduate School of System Informatics, Kobe University
Rokkodai-Cho, Nada-Ku, Kobe 657-8501

Planar Defect in Layered Perovskite
Photocatalyst Y2 Ti2 O5 S2
Layered perovskite Y2 Ti2 O5 S2 (YTOS) isa promising semiconductor photocatalyst with
an electronic structure suitable for overall water splitting under visible light. Similar toother photocatalysts, structural defects during
synthesis should be controlled to promote photocatalytic performance [1, 2]. Very recently,
transmission electron microscopy (TEM) haverevealed extremely large planar defects composed of S–Mg–S layers in YTOS synthesizedusing the flux method. It is not clear how sucha characteristic planar defect makes an eﬀect
on photocatalytic activities in YTOS.
Motivated by this observation, we investigated the planar defect structure and electronic structure using first-principles calculations based on the density functional theory(DFT) [3]. The DFT calculations were performed within the PBE+U (Perdew–Burke–
Ernzerhof functional with Hubbard U correction) exchange-correlation functional using the
projector augmented wave (PAW) method, as
implemented in the Vienna ab initio simulation package (VASP). Here we set U = 7.5 eV
for the 3d orbitals at the Ti sites to reproduce
the experimental band gap, 1.9–2.0 eV [1].
We have microscopically determined the planar defect structure using the structural optimization, where a Mg atom is centered in
the ab-square intersecting the defect S layers(Fig. 1). The obtained defect structure is consistent with the ADF–STEM image: for example, the distance between the Y layers of thedefect model, 9.93 Å, is in quantitative agreement with the experimental value, ca. 10.00
Å.
(b)

Y layerr

Y
Ti
O
S

(a)
Defect-S

9.93 Å

Defect-Mg
Defect-S

Y layerMg(001)

Figure 1: The optimized structure of (a) the
conventional unit cell and (b) supercell with
the MgS2 planar defect.
To understand the electronic structures, wehave investigated the projected density ofstates (PDOS) and band structure shown in
Fig. 2. The electronic structure shows that the371

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

defect bands appear within the original bandgap. The large dispersion of the bands reflects
the planar (two-dimensional) nature of the defect in S and Mg with a high state density,
in contrast to point defects. It is expectedthat such defect bands considerably increase
the electron–hole recombination rate, enhancing the interband transitions of the defect valence and/or defect conduction bands via radiative or nonradiative transitions, which is anundesirable factor for the eﬃciency of the overall water splitting of YTOS.

Mg-poor

60

Total
Y−4d
Ti−3d

(a)

PDOS (states/eV)evant compounds in the synthesis of YTOS.
Fig. 3 shows Ef as a function of the chemical potential of S relative to bulk precipitate(∆µS ), where A, B, and C denote the representive points of the phase diagram. The evaluation of the formation energy suggests that Ef
increases with decreasing ∆µS and/or ∆µMg .
Therefore, the optimum conditions for the S
and Mg environments must be optimized forthe synthesis of YTOS,which suppresses theformation of the planar defect.

50

O−2p
S−3p
Mg−3s

40

Mg-rich30
Defectlevels20

S-poor

S-rich

10
0
−6

−4

−2

0

2

4

6

E−EV (eV)Figure 3: Formation energy of the MgS2 planar defect Ef against ∆µS . The red, blue, andyellow lines represent Ef at points A, B, andC of the phase diagram, respectively.

6(b)E−EV (eV)

4
2
0
−2

References

−4
−6
Γ

X

M

ΓZ

R

A

Z

k−pathFigure 2: (a) PDOS for YTOS with the planar defect. (b) Band structure with the planar
defect. The defect S-3px,y , S-3pz , and Mg-3s
weights are colored in cyan, purple, and orange, respectively. The energy (E) for eachcase is measured from the valence band maximum (EV ).

[1] Q. Wang et. al.: Nat. Mater. 18, 827(2019).
[2] Z. Pan et. al.: Res. Chem. Intermed. 47,
225 (2021).
[3] M. Nakabayashi, K. Nishiguchi, X. Liang,
T. Hisatomi, T. Takata, T. Tsuchimochi, N. Shibata, K. Domen, and S.
L. Ten-no: J. Phys. Chem. C (in press:
https://doi.org/10.1021/acs.jpcc.3c00820).

To consider the stability of the planar defect,
we have finally evaluated the defect formationenergy (Ef ). The chemical potentials µi were
determined by the phase diagram of the rel-372