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

3.5

SCCMS Projects

326

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

Development of molecular simulation framework forstudying mesoscale liquid dynamics
Hayato SHIBA
Information Technology Center, University of Tokyo
Kashiwanoha, Kashiwa, Chiba 277-0882

I now work on a large-scale molecular simulation of the dynamics of liquids. This year
was in the phase of developing new code. The
achievements are described in two parts in thefollowing, in which simulations on the scale ofmillions of atoms (particles) are enabled.

1

Development for large-scale electrolyte simulators

Dilute electrolytes, water with a small amountof dissolved ions, are indispensable for functional properties in, for example, biological
function and batteries. For the purpose of
elucidating many-body properties produced bycorrelated ions, million-atom scale molecular
dynamics simulations may provide valuable insights by directly incorporating the cooperative eﬀects induced by electrostatics.
In this year, I have developed molecular simulation code for TIP4P/2005 water + ions,based on Framework for Developing Particle
Simulator (FDPS) [1], FDPS provides a librarythat automatically produces spatial decomposition by the recursive multi-section method
associated with tree construction for for fastcommunication between diﬀerent MPI processes, which means that the users can concentrate on implementing force ﬁelds and time development schemes. Importantly, tree method
ensures scaling with O(N log N ), bypassing thenecessity for using FFT that prevents parallel scalability. Currently, I have achieved the
speed of 12,000 steps per day for a systemwith 192,000 atoms. Further application tuning is ongoing, targeting 50,000 steps per day

on ISSP System B (Ohtaka).

2

Development simulators for longtime correlations in glasses

When a liquid is rapidly cooled down, it turns
into a glass with solidity, keeping the random
atomic conﬁgurations. About two years ago,
I have revealed power-law correlations in thelong-time behavior of the velocity autocorrelation function in such a supercooled liquidin two dimensions. In order to fully address
such a long-time correlations both in two andthree dimensions, it is necessary to calculate
time-dependent trajectories of each individual particle. This brings about diﬃculties in
MPI parallelization using the spatial decomposition, because particles are no longer boundto the memory of a speciﬁc CPU. In FY2020,
I have added some time-series analysis codefor shared memory environment to my existing distributed parallel code for a LennardJones binary mixture with Kob-Andersen parameters. The whole codes are also extended
to equip CUDA routines, enabling us to calculate the autocorrelations on shared memory
environment by accelerating the force-ﬁeld calculation on GPUs. I am currently using thisfor studying space and time correlations in 2D
and 3D glass-forming liquids with O(107 ) particles.

References[1] http://fdps.jmlab.jp327

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

Creating a Wannier function database toward material design
Kazuyoshi YOSHIMI
Institute for Solid State Physics, University of Tokyo, Kashiwa 277-8581
As a numerical scheme for non-empirical

ESPRESSO (QE) [4] and Wannier90 [5]. In

analysis of electronic

states in strongly

order to create a dataset, we have modified a toolcorrelated materials, the development of a

cif2qewan [6] made by T. Koretsune group at

calculation method based on first-principles

Tohoku Univ. to output the input files for QE

calculations has been actively studied in recent

calculation from cif files containing crystal

years. In this method, the effective Hamiltonian

structures. By downloading cif files from thedescribing the low-energy degrees of freedom isMaterials Project [7] under certain conditions

derived from the band structure obtained by theand using them with this tool, we performed

first-principles calculation [1]. The derived ab-comprehensive

initio effective Hamiltonian can be analyzed

functions and created a dataset.

with high accuracy using low-energy solvers

First, we performed exhaustive calculations forsuch as numerical exact diagonalization [2] andabout 1000 substances as a trial calculation. As

multivariate variational Monte Carlo methods

a result, the calculation was successful for about[3] to reveal the electronic states of real solids

400 materials, and about 10 Wannier functions

in a non-empirical manner. On the other hand,were generated for each material. In order tothese analyses require a high level of expertiseincrease the efficiency of data generation, it wasand experience in the generation of the Wannierfound that it could be reduced to about 10% byfunction, which is required in the derivation ofcontrolling the parameter values. We also foundthe effective Hamiltonian. If this task can bethat the bottleneck area, which takes a lot of

automated, it will be possible to automaticallycomputation time in Wannier90, could be

generate model analyses for various materials,

speeded up by MPI parallelization, and we

which will greatly expand the framework ofsucceeded in achieving a 10-times speedup aftermaterial design.

making changes.

In this project, we have developed a tool for

In future research, we plan to perform high-automatic generation of Wannier functions andthroughput calculations using the tools createdmade dataset for machine learning. The Wannier

and accelerated in this project to create datasets

functions can be calculated using the first-of Wannier functions, and also to try toprinciples

automatically generate Wannier functions using

calculation

softwareQuantum

328

calculationsofWannier

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

machine learning such as deep learning. This

[2] M. Kawamura, K. Yoshimi, T. Misawa et

research was done by collaboration with T.

al., Comp. Phys. Commun. 217, 180 (2017).

Koretsune, K. Kurita, and T. Oikawa of Tohoku[3] T. Misawa, S. Morita, K. Yoshimi et al.,Univ., T. Misawa of Beijing Academy ofComp. Phys. Commn. 235, 447(2019).

Quantum Information Sciences, K. Ido of[4] P. Giannozzi, O. Andreussi, T. Brumme et al.,Institute for Solid State Physics, and TOYOTA

J. Phys.: Condens. Matter 29, 465901(2017).

MORTOR CORPORATION.

[5] http://www.wannier.org[6]https://github.com/wannier-utils-

References

dev/cif2qewan[1] K. Nakamura, Y. Yoshimoto, Y. Nomura et[7] https://materialsproject.org

al., Comp. Phys. Commun. 261, 100781(2021).

329

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

Computational Investigation of Charge
Photogeneration in Organic Solar Cells by
Fragment-Based GW/BSE Method
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. We have recently
developed the large-scale GW method basedon the fragment molecular orbital method.
The fragmentation approximation for the total polarization function and the ∆COHSEX
approximation have been employed for eﬃcient evaluations of GW quasiparticle energies
from localized molecular orbitals. More recently, we have proposed the novel large-scale
excited-state methods based on the FMO, exicton model, and GW/BSE with Tamm-Dancoﬀ
approximation (TDA). In this method, theexcited-state Hamiltonian matrix elements at
the GW/TDA-BSE level is calculated in thebasis of fragment conﬁguration state functions
which describe intrafragment excitations or interfragment charge-transfer excitations. The
excited-state Hamiltonian is then diagonalized
to approximate the adiabatic excited states ofan entire system. We have conﬁrmed thatfor molecular clusters, the novel fragmentbased GW/TDA-BSE can reasonably reproduce the excited states obtained from conventional GW/TDA-BSE method. Our developments enables the applications of GW/TDABSE method to large molecular assemblies,
which will be useful for investigating charge
separation in condensed phase.

energy levels and charge-transfer states ofthe face-on and edge-on orientations of pentacene/C60 bilayer heterojunctions. FMOGW/BSE calculations were performed for the
local interface structures in the face-on andedge-on bilayer heterojunctions, which containapproximately 2000 atoms. Calculated energy levels and charge-transfer state absorption spectra are in reasonable agreements withthose obtained from experimental measurements. We found that the dependence ofthe energy levels on interfacial morphology is
predominantly determined by the electrostatic
contribution of polarization energy, while theeﬀects of induction contribution in the edgeon interface are similar to those in the faceon. Moreover, the delocalized charge-transfer
states contribute to the main absorption peakin the edge-on interface, while the face-on
interface features relatively localized chargetransfer states in the main absorption peak.
Our results indicate that optimizing the molecular orientation and interfacial morphologies isessential for improving the power conversion
eﬃciency of organic solar cells.

References[1] T. Fujita, Y. Noguchi, T. Hoshi, Materials, 13 (2020) 2728.

As recent applications [1], we explore the330Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo!"#$%&%'()*+,-#./)0+,/"%+1&)2#")3#.%4$5%#')6/7-+"(/+81/)
9+,,/"%/&:)!"#$%$&$')%';/&,%(+,%#'&!""#$"%#&'()*!"#+#!,-./.01,#2(23$4#"$+#567,.#(*'8!3#"$%!

!"#$%#&'()*+,,-(,.(/0120''"2013(45#6#(7028'"52&93()%2&#3(45#6#(:;:<=>?@3(A#B#0(#

C)CD3(45#6#(7028'"52&93(CE#"#623(45#6#(:;?<==F?3(A#B#0($/)CGH3(I9,&,(7028'"52&93(I9,&,(;@:<>JF:3(A#B#0(!"# $%&"# '()*+"*# (,-# (./"'# -0# /1-2+'+34#($"# 2%("1+%6# ')00"1'# 01-2# ($"# '"1+-)'# /$%'"#

5%($-*"# 2%("1+%6'# %//6+5%76"# 0-1# '-*+)28+-3#(1%3'0-12%(+-3# *)1+34# +('# 5$%14+34E*+'5$%14+349#

1"5$%14"%76"#7%(("1+"'9##H'+34# ($"# 0+1'(# /1+35+/6"'# *"3'+(.# 0)35(+-3%6#($"-1.G##,"#"</6-1"*#'.'("2%(+5%66.#+('#

4"-2"(1+5G#"6"5(1-3+5#/1-/"1(+"'G#/$%'"#'(%7+6+(.#%3*# *+00)'+-3# 2"5$%3+'2# 0-1# =%"# 2+41%(+-39# I(#+'# 0-)3*# ($%(# %(# $+4$# -1# 6-,# =%# 5-35"3(1%(+-3G#=%"#+-3'#/1"0"1#6-5%(+34#-36.#%(#($"#"*4"8'$%1"*#/1+'2G# ,$+6"# %(# !J@9BG# ($"# 2%("1+%6# 7"5-2"#
2-1"#'(%76"#,+($#($"#1%(+-#-0#($"#=%#+-3#%(#"*4"8#!"#$% &'% % ()*+*,-./0% 1"2234"5,% 6/57*44*4% .-% 23))%'$%1"*# =%## %3*# 0%5"8'$%1"*# =%$# /1+'2# -0#=%#K=%$ ! ;KL9# !$"3# %66# ($"# =%"# +-3'# %1"#

78./#",#%4-.-*49%",7)31",#%:.;%7/544",#%6/57*44%.,1%:<;% 6./.))*)% 6/57*44% .,1% .7-"=.-"5,% *,*/#0% 6/52")*%*"+3("15%6%("*# 01-2# =%>3?!G# ($"# 5-3'/+5)-)'#

52%>.%=.7.,70?654"-"=*%65)./5,%75+6)*@%5<-.",*1%

5$%34"# +3# ($"# 6%((+5"# 5-3'(%3(# "# 5%)'"'# %3#<0% 5<-.",*1% <0% :7;% AABC!% .,1% :1;% DE(FG%'+43+0+5%3(# 1"*)5"# -0# &-6)2"# )/# (-# ;@M9# !"#+*-8514$%H8*%1./I%<)3*%57-.8*1/5,4%/*6/*4*,-%-8*%/-+3("*# -)(# ($%(# ($"# $.71+*+N"*# -17+(%6'# -0# ?O'#

J>>%K,%4"-*4%L8*/*%65)./5,%25/+*19%.,1%-8*%1./I%;#"#'(%("'#%3*#>3O'#A$"##%3*#A$!!%"!#'(%("'G#/6%.#70.,% <.))4% 4-.,1% 25/% >.% =.7.,70% -/.7*% 1"2234"5,$%%3# "''"3(+%6# 1-6"# +3# ($"# "6"5(1-3+5# '(1)5()1"# -0#

H8*% 5-8*/% >.% .-5+4% ./*% 8"11*,% 25/% +5/*%($"# 2%("1+%69# !"# /1-&"*# ($"# /-6%1-3# 0-12%(+-3#7)./"2"7.-"5,$!"#%%3*#($"#/-6%1-3#+'#/1"0"1%76.#5-)/6+34#,+($#=%"#%&%5%35.E+-3# %(# ($"# ($+1*# 3"%1"'(# 3"+4$7-1# '+("'9#!"! #$!%&'"()$*+,+-(./0-+((P,-# "6"2"3(%1.# *+00)'+-3# /1-5"''"'G# +356)*+34#:;# 6%."1"*# -<+*"# =%!>3?!# "<$+7+(# %#/%1%66"6#%3*#51-''+34#/1-5"''"'G#%'#'$-,3#+3#Q+49#&-6(%4"# ,+3*-,# 01-2# ;9@# (-# A9BC# %3*# $+4$#

LG#5-)6*#-55)1#%3*#1"R)+1"#%3#%5(+&%(+-3#"3"14.#'/"5+0+5# 5%/%5+(.# )/# (-# ;;@# 2D$E49# F-,"&"1G#-0# S;A# 2"C# TBLU# 2"CV# %(# $+4$# =%#331

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

concentration and 273 meV (327 meV) at low

meV in the Na-rich and poor contents,

Na content by GGA+U (HSE06) method.

respectively. Therefore, the proposed cathodeCompared to olivine phosphate, the Na+ ion

can be expected to exhibit significantly fasterwould diffuse better in NaxMnO2 than in

Na-ion diffusion than that in popular cathodespolyanionic framework like LiFePO4 due to thesuch as NaVOPO4 (627 meV).[2]significantly lower activation energy.[1]2. NaxVPO4F
A new orthorhombic phase of tavorite-like
NaxVPO4F cathodes (0 ≤ x ≤ 1) applicable for
Na-ion battery was proposed using density
functional theory. We have explored the phasestability, electronic structures, electrochemical
properties, and diffusion mechanisms of Na
ions/vacancies of NaxVPO4F. No imaginary

Fig. 2: Possible elementary diffusion processes offrequency in phonon dispersion provides

an Na vacancy (vc) accompanying by a positive

evidence that the proposed material is thermallypolaron (pp) along the [010] direction instable. It is suggested that the material can be

NaVPO4F. O (red), F (blue), V (brown), P (grey),

synthesized from the reaction: NaF + VPO4 à

Na (cyan), and vc (black) are shown by balls, and

NaVPO4F because of the negative formation

polyhedral structural units of V1NNO4F2 (brown),

energy. The material NaVPO4F would undergo

V2NNO4F2 (green), AO4F2 (yellow) octahedra, andbi-phase reactions upon the Na deintercalation

NaO3F (purple) tetrahedra are shown. When vc

in the ranges of 0.125 ≤ x ≤ 0.5 and 0.5 ≤ x ≤ 1,

moves along a path from vcOA to vcOB via vcTand show high voltages of 3.80 V and the highsites, pp can jointly hop from A to C sites (parallel

theoretical capacities of 143 mAhg−1. The

process indicated by a green arrow), from A to B

polaron migration effect on the Na+ ion

sites (crossing process indicated by a blue arrow),

diffusion mechanism was precisely investigated.

or remain at a same site (single process).

Three elementary diffusion processes, including
single, crossing and parallel processes, are

Referencesindicated in Fig. 2. The parallel and crossing[1] H. D. Luong, V. A. Dinh, H. Momida, T.

processes require much higher activation

Oguchi, Phys. Chem. Chem. Phys. 22, 18219

energy for Na+ diffusion, implying a substantial(2020).

polaronactivation[2] H. D. Luong, V. A. Dinh, H. Momida, T.

barriers of Na ion diffusion gain 323 and 530

Oguchi, J. Alloys Compd. (2021) accepted.

migrationeffect.

The

332

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

Ab initio study toward abundant element
nanocatalysts with less precious metals
Tetsuya TAKETSUGU
Faculty of Science, Hokkaido University
N10W8, Sapporo, Hokkaido 060-0810

DFT computations with QuantumEspresso(QE) are used to gain chemical insights andto design novel catalysts with abundant elements such as hexagonal boron-nitride (h-BN)
for reducing costs and dependences on precious
metals used in various heterogeneous catalysts.
We have shown that gold clusters supportedon h-BN/Au(111) are highly active catalysts
for hydrogen evolution reaction (HER) [1]. It
is demonstrated that the most stable gold clusters on the h-BN/Au(111) support are not necessarily the most reactive; stabilization on thesurface gold clusters with the diﬀerent geometries and orientation with respect to the support can form a large variety of energeticallyfavorable H adsorption sites considerably promoting HER activity. The charge transfer between the supported gold clusters and the hBN/Au(111) surface promotes the catalytic activity of the Aun @h-BN/Au(111) system for
HER. Therefore, h-BN/Au(111) can be considered as promising support for gold clusters
as electrocatalysis [1].
Another application of BN is for oxidative
dehydrogenation (ODH) of light alkanes suchas olefin. Olefins are key feedstocks for rangeof commercially important chemicals such as
polymers, fibers, and their derivatives, which
are important for the energy and day to day
applications. Therefore, a deeper insight into
the mechanism and ways to explore catalysts
for eﬃcient oxidative dehydrogenation reaction is needed. We have proposed that oxygen functionalized h-BN with B-O-O-B sites

at the zigzag edges is active catalyst for ODH.
This suggestion is consistent with experimentally reported concentration of B and N atoms(high B:N ratio) to explore the ODH mechanism in a series of light alkanes [2]. DFT
calculations are performed to verify our proposal. We have explored the catalytic activity of the newly proposed oxygen functionalized edges of h-BN sheets for oxidative dehydrogenation of ethylbenzene. It is demonstrated that oxygen functionalized highly defective boron-nitride can demonstrate extraordinary catalytic activity for the above reaction[2].
We have elucidated the mechanism of thenovel ”heterocyclic segregation” reaction of
one pyrazine ring with two nitrogen atoms
into two quinoline rings with one nitrogen
each. Our DFT calculations showed that the
intramolecular ring-forming and -opening of
nitrogen-heterocycles are strongly aﬀected by
the initial hydrogen-substrate interaction [3].

References[1] M. Gao, M. Nakahara, A. Lyalin, and T. Taketsugu, J. Phys. Chem. C 125, 1334-1344 (2021).
[2] R. Han, J. Diao, S. Kumar, A. Lyalin, T. Taketsugu, G. Casillas, C. Richardson, F. Liu, C. W.
Yoon, H. Liu, X. Sun, Z. Huang, J. Energy Chem.
57, 477-484 (2021).
[3] K. Sun, K. Sugawara, A. Lyalin, Y. Ishigaki,
K. Uosaki, T. Taketsugu, T. Suzuki, and S. Kawai,
Angew. Chem. Int. Ed. 60, 9427-9432 (2021).

333

Activity Report 2020 / 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 applied the scheme to (Nd1−α−β−γPrα Laβ Ceγ )2 (Fe1−δ−ζ Coδ Niζ )14 B [1]. In thefirst step, we prepared 119 experimental samples, and measured the magnetization at 3-7
temperatures within a range from 300 K to473 K. The magnetization at zero temperature (µ0 M ) and Curie temperature (TC ) areevaluated from these data by Kuzmin’s formula. Independently, we performed systematic first-principles calculation for 2869 compositions using the KKR-CPA method. The
saturation magnetization at zero temperature

is calculated in the local spin density approximation. The intersite exchange couplings
are calculated using Liechtenstein’s method,from which the Curie temperature is evaluated in the mean-field approximation. We thenadopted the data assimilation, and obtainedµ0 M and TC at arbitrary composition. The
magnetization at arbitrary temperature is obtained using Kuzmin’s formula. Figure 2 shows
the magnetization at 0 K, 300 K and 400 K.
We see that the magnetization monotonicallydecreases with increasing Co concentration (δ)
at 0 K, whereas small amount of Co doping
enhances the magnetization at 400 K.

5
4
3

y

High-throughput first-principles calculation
is a powerful tool for developing functional materials. However, systematic errors containedin computational data often prevent us fromquantitative evaluation of physical properties.
In order to overcome this problem, we have developed a data assimilation method in which asmall number of experimental data are integrated with a large number of computational
data [1]. Figure 1 is a schematic illustrationof the method. Here, open circles denote accurate but a small number of data, whereas
open squares are a large number of data that
contain a systematic error. If we consider theformer data only, we cannot estimate the target variable accurately for large x where thedata do not exist, whereas accurate estimationis possible if we consider both the data. Basedon this framework, we have developed a practical scheme for evaluating finite-temperaturemagnetization.

true model (i)sample (i)OLS (i)assimilation (i)true model (ii)sample (ii)
OLS (ii)
assimilation (ii)

2
1
0
1
0.0

0.2

0.4

x

0.6

0.8

1.0

Figure 1: An example of the data assimilation.

334

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(Nd 1 − α − β − γ Pr α La β Ce γ ) 2 (Fe 1 − δ − ζ Co δ Ni ζ ) 14 B,
1.0

0.5

(α, ζ): (0.00, 0.00)

γ = 0.00

1.0

0.5

δ

δ = 0.00

δ

γ

0K

1.0

β = 0.00
1.8
1.7

0.5

1.6

β

0.0
0.0

1.0

0.5

β

0.0
0.0

1.0

1.0

1.0

0.5

0.5

0.5

δ

1.0

0.5

γ

1.0

1.5
1.4
1.3
1.2
1.1

0.5

β

0.0
0.0

1.0

0.5

β

0.0
0.0

1.0

1.0

1.0

0.5

0.5

0.5

0.0
0.0

0.5

β

1.0

0.0
0.0

δ

1.0

δ

γ

0.0
0.0

400 K

0.5

δ

γ

300 K

0.0
0.0

0.5

β

1.0

0.0
0.0

0.5

1.0

0.5

1.0

γ

γ

1.0

Figure 2: Magnetization of (Nd1−β−γ Laβ
Ceγ )2 (Fe1−δ Coδ )14 B.

References[1] Yosuke Harashima, Keiichi Tamai,
Shotaro Doi, Munehisa Matsumoto,Hisazumi Akai, Naoki Kawashima,
Masaaki Ito, Noritsugu Sakuma, Akira
Kato, Tetsuya Shoji, and Takashi Miyake,Phys. Rev. Materials 5, 013806 (2021).

335

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

Optimized design of magnetic materials based on the
integration between quantum beam experiment andfirst-principles calculation
Kanta Ono
Institute of Materials Structure Science, High Energy Accelerator Research Organization
Oho 1-1, Tsukuba 305-0801

We have been working on characterization and material design of rare-earth permanent magnets (REPMs) which is one of themost important materials for traction motorsof vehicles and robotics. Quantification of
REPMs has to be done on multiple lengthscales and energy scales, where our scope liesin the microscopic thermal equilibrium andoff-equilibrium magnetization reversal process
closer to the mesoscopic scales [1]. Since themain-phase ferromagnets in REPMs are typically made of multiple sublattices and eachsublattice contributes to the bulk magnetic
properties differently, it is an important problem to quantify the site-resolved contribution
from the microscopic measurements to understand which sublattice can be more improved [2]. Based on the measured data fromthe neutron diffraction experiments, the Rietveld analysis can give us such information
which so far has needed expertise of experienced practitioners. Still, unphysical solutions
can be given as a metastable solution in theoverall landscape of the cost function in thedata fitting. In order to quickly select a physically meaningful solution and eliminate human bias as much as possible, we have beendeveloping a data analysis technique combining Rietveld analysis and theoretical calculations [3, 4]. Ab initio electronic structure calculations for metallic ferromagnets can addresssite-specific magnetic moments based on the
crystal structure data yielded by the Rietveld

336

analysis. The output of the ab initio calculation can be fed back into the Rietveld analysis and such hybrid process between measured data analysis and ab initio calculations
can be iterated until the overall convergence isreached.
While the convergence in the hybrid analysis
is quite fast in the ThMn12 -type ferromagnets,
feasibility for the Nd2 Fe14 B-type ferromagnets
has been seen only recently with the new system B where the 128 cores can shorten the
ab initio part based on AkaiKKR in less thanan hour. This is in strong contrast to otherelectronic structure calculations for Nd2 Fe14 B
which can take weeks on older machines. With
the new hardware, the good ab initio calculation framework, and our new hybrid data
analysis methodology, further developments of
REPMs will be accelerated.

References[1] H. Tsukahara, K. Iwano, T. Ishikawa,
C. Mitsumata, K. Ono, NPG Asia Materials 12, 29, (2020).
[2] K. Saito, S. Doi, T. Abe, K. Ono, J. Alloys
and Compounds 721, 476 (2017).
[3] M. Matsumoto, T. Hawai, K. Ono,
Phys. Rev. Applied 13, 064028 (2020).
[4] T. Hawai, M. Matsumoto, K. Ono, in
preparation.

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

Evaluation of the polymer blend miscibility by using
chain-increment method with all-atom molecular
dynamics simulation
Kazuo YAMADA and Nobuyuki MATUBAYASIDivision of Chemical Engineering, Osaka University, Toyonaka, Osaka 560-8531
PE

PVDF, 423K

PVDF, 553K

PVA

PVP

70

80

70

80

0

(kcal/mol)

Polymer blend, a mixture of polymers, may

-2
-3

∆µi

incr

have distinct structural and dynamical features-1

from the component polymers and is employed

-4
-5
30

in a wide variety of applications including

40

50

60

monomer index i-5<u>i (kcal/mol)

home appliances, biomedical devices, andautomotive products. The mutual miscibility of

-10
-15
-20

polymer species is determined by the free

30

40

50

60

energy of mixing. The mixing free energy is

Fig 1: Incremental free energy ∆𝜇!"#$% and the average interaction energy

given by the changes in the (excess) chemical〈𝑢〉! against the index 𝑖 of the incremented monomer in the tagged

potentials between the pure and mixed states ofpolymer at 𝑖 = 30, 40, 50, 60, 70, and 80. The horizontal, dashed line isthe component species, and the chemical

the averaged value of ∆𝜇!"#$% or 〈𝑢〉! .

potentials reflect sensitively the effect of intra-

calculation has been conducted at 1 bar on theand intermolecular interactions in turn. To

supercomputer system at ISSP. We compute theapproach the chemical potential of a polymer

incremental

molecule, we formulated the method of chain

averaged interaction energy of the incremented

increment and provided a scheme to computemonomer with the surroundings 〈𝑢〉! . Figure 1

the free energy of incrementing a monomer

shows that 𝜇!"#$% stays constant within a few

with all-atom models.[1] In this work, we

tenths of kcal/mol for inner monomers,

extended the method and treated the single-

demonstrating the validity of the chain-component(PE),

increment method. 〈𝑢〉! is also constant against

polyvinylidene difluoride (PVDF), poly(vinyl

the index i of the incremented monomer,

alcohol) (PVA), polyvinylpyrrolidone (PVP) inshowing that the local environment around thetheir linear forms.

monomer is independent of the monomer

Results

location along the polymer chain.

melts

of

polyethylene

freeenergy ∆𝜇!"#$% andtheReferences

The number of monomer units is fixed at 100
for all the systems. All-atom molecular[1] K. Yamada and N. Matubayasi: Macromol.

dynamics

53 (2020) 775-788.

simulation

andfree-energy

337

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

Development of high-throughput calculation tools
and evaluation of magnetic properties in hardmagnetic materials
Tetsuya FUKUSHIMA
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

In this year, we developed (i) an automatic exhaustive calculation tool and (ii) an
atomistic spin Monte Carlo simulation tool
for permanent magnetic materials. Our automatic exhaustive calculation tool is basedon the all-electron Korringa-Kohn-Rostoker(KKR) Green’s function method (AkaiKKRcode). The KKR Green’s function method
can be combined with the coherent potential
approximation (CPA) which can treat configurational disordered systems, is quite useful for investigating the electronic structures,
magnetic properties, and transport propertiesof rare-earth transition-metal alloys. In the
CPA approach, the multiple scattering effectby disordered potentials is replaced with effective medium. Therefore, we can perform theelectronic structure calculations for disorderedmaterials without supercells and decrease thecomputation costs dramatically. Additionally,
since a Green’s function is directory obtainedin the KKR method, we can efficiently investigate finite temperature magnetism and transport properties by the linear response theory.
We applied this tool for and high entropy alloys
and rare-earth magnet Sm-Fe-N systems and,succeeded in constructing large-scale magnetic
material databases.
Atomistic spin Monte Carlo simulation wasused for investigating the temperature dependent magnetization reversal in permanent
magnets. The microscopic mechanism of thecoercivity at finite temperature is a crucial

Figure 1: Free energy landscape for a isolated
Nd2 Fe14 B particle by atomistic spin model
simulation. The size of the particle is 14 nm,
corresponding to 212,536 spins.
issue. Since the magnetization reversal is astochastic transition from a metastable state
to a stable state, we tried to capture the physical insight of the coercive force from the viewpoint of the free energy landscape. For the
calculation of the free energy landscape, wedeveloped a numerical simulation method fordescribing the magnetization reversal by nucleation, where highly efficient parallelization
scheme (around 1,000–3,000 MPI parallelization) for the replica exchange Wang-Landau
MC method and the reduction of systematic
errors by the 1/t method are implemented.
Applying the above method to the atomistic
spin model with a isolated Nd2 Fe14 B particle,
we performed the magnetization reversal anal-

338

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

long-ranged interaction in computational costO(N logN ) without any approximations. Thismethod can be combined with Monte Carlosimulation, leading to the quantitative analysis of the temperature dependent coercive
forces of permanent magnets from the free energy landscape calculations with the magnetic
dipole interactions (see Fig. 2). It is foundthat the influence of the magnetic dipole interaction (demagnetizing field) is strongly suppressed due to the surface thermal fluctuation.

References
Figure 2: Coercive energy of Nd2 Fe14 B atroom temperature as a function of the size ofthe domain. The blue (red) circles indicate thecalculation results without (with) the magnetic
dipole interactions, respectively. White circles
are the experimental results.

[1] Y. Toga, S. Miyashita, A. Sakuma, and T.
Miyake, Npj Comput Mater 6 (2020) 1.
[2] T. Hinokihara, M. Nishino, Y. Toga, and
S. Miyashita, Phys. Rev. B 97 (2018)
104427.

ysis for the effect of the thermal fluctuationand investigated the formation of the reverse
nucleus which is a origin of the coercive force(see Fig. 1). It was considered that the barrierof the magnetization reversal is determined by
the competition between the formation energyof the reverse nucleus and the Zeeman energy
by external magnetic field: however, we newly
proposed that the energy gain due to the expansion of the magnetization reversal region by
the domain wall motion is a key role. This fact
is quantitatively consistent with experimental
results, where the size of the reverse nucleus
dose not depend on external magnetic field.
In order to investigate the coercive forceof permanent magnets quantitatively, we needto consider the effect of the magnetic dipole
interactions (demagnetizing field). Althougha convolution integral method (fast Fouriertransform) is commonly used to incorporate
long-ranged interaction, it is not effective forpermanent magnet materials with low crystal
symmetry such as Nd2 Fe14 B. We have developed a modified stochastic cutoff method thatis applicable complex crystal structures with339

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

Novel chloride solid electrolytes for all solid-state sodium
metal battery
Masanobu NAKAYAMA
Department of Advanced Ceramics
Nagoya Institute of Technology, Gokiso-chi, Showa-ku, Nagoya, Aichi 466-8555
Using high-throughput computing and structure

which is a derivative of original BVFF, todatabase (Materials Project [1]), we evaluatedreproduce DFT-MD derived structural features,

exhaustively migration energies of Li, Na, andsuch as RDF and ADF.

Mg ion in inorganic solids. Figure 1 presents

approach is adopted to calculate the ion

the distribution of migration energies of Li, Na,conductive performance.

and Mg ions in oxides by using BVFF-

calculated Na ionic conductivities at room

percolation approach.[2,3] The results indicate

temperature for 86 compositions.

that various factors affect the ion conducting

verify reproducibility of present computation

performance. Thus, materials informatics (MI)and apply MI approach for thus obtained dataapproaches are adopted to clarify the crucial

near future.

factor for ionic conductivity in oxides. Among

7

various6

regression

approaches,

gradient

performance.
performance

However,
is

not

the

sufficient

Migration energy /eVboosting regression shows the best prediction
regression
from

the

diagnostic plots shown in Fig.2. Our tentative

Figure 3 shows
We plan to

Li
Na
Mg

5
4
3
2

analysis suggests that Li-O bond length and Li-

1

O-Li bond angles are important, i.e. structural

0
10

15

20

25

Cell volume per oxygen atom /Å

factors are more important than compositional
factors.

After that FF-MD

30
3

One of the reasons for rather poor

Fig. 1: Migration energies for Li, Na, and Mg

regression performance may stem from low

ions (indicated by orange, blue, and green plots,

reproducibility

respectively)ofuniversal

BVFF

parameters.[2] To improve, high-throughput

exhaustive BVFF approach.

FF approaches[4] are used for 86 Na containing
materials.

In

this

technique,calculated

DFT-MD

simulations are performed just for 2 ps, andthen metaheuristically optimize FF parameters,

340

using

automated

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

References[1] A Jain, S. P. Ong, G. Hautier, W. Chen, W.D.
Richards, S. Dacek, S. Cholia, D. Gunter, D.
Skinner,G.

Ceder,

K.A.

Persson,APL

Materials, 2013, 1, 011002.
[2] S. Adams, R. P. Rao, Phys. Status Solidi
Fig. 2 diagnostic plot of the gradient boosting

Appl. Mater. Sci., 2011, 208, 1746–1753

regression derived prediction function for Li[3] M. Nakayama, M. Kimura, R. Jalem, T.

migration energies in 711 Li-ion containing

kasuga, J. Appl. Phys. 2016, 55, 01AH05

oxides.

[4] R. Kobayashi, Y. Miyaji, K. Nakano and M.
Nakayama, APL Mater., 2020, 8, 081111

Figure 3 Histogram of room temperature Na
ion conductivities for 86 samples using highthroughput FF-MD approach.

341