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

3.4

Cooperative Phenomena
in Complex Macroscopic Systems

213

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

Development of new tensor-network-based methods andstudy of critical phenomena in quantum/classical magnets
Naoki KAWASHIMA
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Tensor-network is now widely recognized as

approximation produces a qualitatively correcta new framework/language for developing

phase diagram, the 2D iPEPS calculation

numerical methods for otherwise intractableyielded a quantitatively accurate one. We also

many-bodytheobserved that iPEPS calculation became slower

structure of the problem from the information-in pinning down the direction of the magnetic

scientific viewpoint. In SY2022, we developed

easy axes as one approaches the phasea few new numerical methods

for studying

boundary --- an indication of the ``easy-axis

quantum spin system in two spatial dimensions

softening''. In fact, the low-energy effective

and classical statistical-mechanical models

theory in the reduced Hilbert space restricted inbased on tensor-network representation. We

the lower atomic levels indicate the asymptotic

then applied them to some typical examples.

U(1) symmetry on the phase boundary, which isproblems

anddescribing

In [1], we studied the square-lattice S=2

consistent with the above-mentioned softening

Heisenberg model with cubic anisotropy. The

phenomena. From these results together withmagnetic anisotropy is quite ubiquitous --- itrenormalization

can be found in virtually every magnet thoughconcluded that the emergence of the continuousin many cases its effect is very weak. However,symmetry

the true nature of critical phenomena can beboundary. Therefore, this criticality belongs todramatically altered by such a weak symmetry

the same universality class as the 3D classicalbreaking field. In spite of the universal nature

XY model. The iPEPS calculation of theof the anisotropy, it is hard to incorporate in thecorrelation functions and other quantities

numerical calculation in an accurate way. We

support this interpretation. Thus, our iPEPS

studied the simplest lattice spin model for the

calculation provided the first systematic andtwo-dimensional quantum ferromagnet withcorrect

cubic anisotropy by means of mean-fieldapproximation, of the phase boundary and theanalysis and tensor network calculation. It

critical phenomena of this simple quantum spin

turnedsystem.

out

that

while

the

mean-field

214

group

becomesanalysis,considerations,exact

beyond

on

the

the

we

phase

mean-fieldActivity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of TokyoIn [2], we studied the phase transitions ofarbitrary two-dimensional lattice can be dealt

the two-dimensional q-state Potts model. While

with by deforming it to a non-uniform square

the universality class of the model is welllattice, though the numerical efficiency of such

established, the scaling function of the Binder

a treatment strongly depends on the original

parameter had not been accurately obtained

problem.

before. In particular, we focused on the nonmonotonic behavior of the Binder parameter.

References

We discovered that for q=3 its non-monotonic[1] Wei-Lin Tu, Xinliang Lyu, S. R. Ghazanfari,behavior as a function of the temperature is not

Huan-Kuang Wu, Hyun-Yong Lee, Naoki

the non-universal properties. Instead it reflects

Kawashima: arXiv:2204.01197.

a bump in the scaling function. Therefore, the[2] Hiroshi Watanabe, Yuichi Motoyama,

non-monotonicity is the universal property of

SatoshiMorita,the Potts model.

Progress

of

and

Naoki

Kawashima,

theoreticalandexperimental

physics 2023, 033A02 (1-19) (2023).

In [3], we developed an open source
software called TeNeS. It is a program package[3]for obtaining ground states of two-dimensional

Kazuyoshi Yoshimi, Satoshi Morita, Takeo

quantum lattice problems. It is based on the

Kato,corner transfer matrix renormalization group

Physics Communications 279 108437 (1-15)(CTMRG) method. In principle, it is designed(2022).

for an arbitrary two-body Hamiltonian and an

215

Yuichi

Naoki

Motoyama,

TsuyoshiKawashima:

Okubo,

Computational

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

Large-scale Molecular Simulation of Soft Materials using AllAtom and Coarse-Grained Models
Wataru SHINODA
Research Institute for Interdisciplinary Science,
Okayama University, 3-1-1 Tsushima-naka, Okayama 700-8530
We illustrate here a simulation study onof peptides in different initial configurations.

action of antimicrobial peptide (AMP), melittin,

One system showed lipid extraction and poreon lipid membranes as an example of ourformation by melittin peptides without inducing

research works using massive parallel computer

membrane curvature, while another showed

simulations. In this study, we conducted a series

large pores in a significantly curved membrane.

of coarse-grained molecular dynamics (CG-(Figure 1) This result indicated that, dependingMD) simulations to investigate the complicatedon the melittin arrangements on the membraneactions of melittin in lipid membrane systems. In

at high melittin concentration, pores werea membrane system with a lower peptide/lipid

formed either by relaxation of bending energy(P/L) ratio than the threshold value for melittin

due to large curvature or by the peptide insertion

pore formation, no pores were found during theinto the thinned membrane due to lipid

simulation, although the peptides formed dimers

extraction.

on the membrane. At a higher P/L ratio than theThrough this study, we have gained somethreshold value, we successfully observed pore

insight into the behavior of melittin, an AMP, on

formation caused by the insertion of melittin

lipid membranes. The membrane pore formation

peptides through their N-terminals. During theprocess of antimicrobial activity cannot bepore formation, it was found that membrane

induced by a single melittin peptide, even at P/Lcurvature was locally induced in the vicinity ofratios above the concentration required for pore

the peptide aggregate. These findings wereformation. This was energetically evidenced byconsistent with experimental observations andthe free energy profile of a single melittinsuggested melittin action. Additionally, we

insertion into the membrane. These results were

performed CG-MD simulations for large-scalealso supported by 5 µs CG-MD simulations.

membrane

peptidePore formation is always induced by cooperative

concentrations twice. Interestingly, two different

actions such as induction of melittin membrane

membrane pore formation processes were

curvature and extraction of lipids from theobserved in the systems with the same number

membrane. These two events are the key factors

systems

at

high

216

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(+*(# ()$>>")# 4-)"# 2-),*($-.@ C--4")*($?"# *1($-.#

2$./$.>'# -2# (+$'# '(&/0# <$%%# 1-.()$3&("# (-# (+"#

-11&)'#<+".#*>>)">*("'#-2#,"%$(($.#1-,4-'"/#-2#/"?"%-4,".(#-2#*.($,$1)-3$*%#*>".('#(+*(#,$,$1#*(#%"*'(#2-&)#-)#,-)"#4"4($/"'#*)"#2-),"/#-.#(+"#

678# *>*$.'(# ?$)*%# *./# 3*1(")$*%# $.2"1($-.@ !"#,",3)*."@ 6%(+-&>+# *>>)">*($-.# -2# 1*($-.$1#+-4"# (+*(# (+"# CD# 2-)1"# 2$"%/# 54Q8AC69#

4"4($/"'# *44"*)'# (-# 3"# /$'*/?*.(*>"-&': (+"#/"?"%-4"/ $.#(+$'#)"'"*)1+#<$%%#+"%4#&'#(-#3"((")#'1)"".$.># "22"1(' -2# "%"1()-'(*($1# $.(")*1($-.'&./")'(*./ (+"#*1($?$($"'#-2#678@$./$1*("# (+*(# *>>)">*($-.# ,*0 -11&)# -.# (+"#,",3)*.":# *' $.# (+"# CDB7E# '$,&%*($-.'@#!"#"$"%&"'

P-<"?"):# '4-.(*."-&'# *''",3%0# -2# ,"%$(($.#

RLS# T@# 7$0*U*=$# *./# !@# Q+$.-/*: V$-1+$,@#

4"4($/"'# /-"'# .-(# ."1"''*)$%0# )"'&%(# $.# (+"#

V$-4+0'@#61(* B V$-,",3)*."':# !"#$: LWXYNN2-),*($-.# -2# 4-)"'# 3"1*&'"# (+"# 2)""# ".")>0#5Z[ZZ9@3*))$")#)"M&$)"/#2-)#4-)"#2-),*($-.#$'#',*%%@ I+"#

K$>&)"#L\#J3'")?"/#,",3)*."#4-)"#2-),*($-.#
4)-1"''"'# $.# *# %*)>"# ,",3)*."# 4*(1+# *(#8FG]NFLZW@217

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

Crossover of ﬂuctuating surfactant adsorbed interface
Shunta Kikuchi and Hiroshi WATANABE
Department of Applied Physics and Physico-Informatics, Keio University
Yokohama, Kanagawa 223-8522, Japan

We investigated the mechanical properties ofa surfactant adsorbed at the interface between
two liquids. The restoring force of the membrane is largely due to interfacial tension andelasticity. These are explained by the elastic
energy model by Helfrich with the following
equation [1].
|h(q)|2 =kB T,
2γq + κq 4

where q is the wavenumber, h(q) is the Fouriercomponent of the ﬂuctuation of the membrane,
γ is the surface tension, and κ is the bending
rigidity. The biomembrane is composed of aphospholipid bilayer, and its restoring force isthought to be either interfacial tension or elasticity, but the observed ﬂuctuations are dominated by elasticity, with no eﬀect of interfacial
tension. We therefore investigated ﬂuctuations
in monolayer membranes to determine how theeﬀects of interfacial tension and elasticity contribute to the resilience of the membranes.
We represented two kinds of liquids by
preparing two kinds of atoms. The same typeof atoms interacted with the Lennard-Jones
potential, and the diﬀerent type of atoms interacted with the WCA potential to reproduce a system in which the two liquids arephase-separated. We modeled the surfactant
as a diatomic molecule and prepared an initial state in which the surfactant is adsorbed
at the interface. Then, molecular dynamics
was used for time evolution to observe the ﬂuctuations at the surfactant-adsorbed interface.
LAMMPS was used for time evolution.
Our investigation revealed a signiﬁcant de-pendence of the interfacial tension within thesystem on the surfactant’s natural length.
Modifying the natural length could alter theinterfacial tension from ﬁnite to nearly zero
while maintaining the system’s temperature
constant. The ﬁndings demonstrate that theinﬂuence of q 2 is discernible in the ﬁlm ﬂuctuations when large values characterize interfacial tension. Conversely, a transition from q 2
to q 4 becomes evident in the Fourier spectrum
when interfacial tension approaches near-zerovalues [2]. This crossover implies that whensurfactant molecules are suﬃciently adsorbed
at the interface, and the interfacial tension isnearly zero, the bending stiﬀness dominates in
the restoring force. In contrast, it plays thedominant role when the interfacial tension possesses ﬁnite values.

References[1] W. Helfrich, Z. Natuforsch C28, 693(1973).
[2] K. Shunta and H. Watanabe, J. Chem.
Phys. 158, 124901 (2023).

218

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

Correlation between the structure of charged
phospholipid bilayer membrane and ion distribution
Yuji HIGUCHI
Research Institute for Information Technology, Kyushu University
Motooka Nishi-ku, Fukuoka 819-0395

We have studied soft matter by molecularsimulations and revealed the self-assembly processes of amphiphilic molecules [1], water dynamics on lipid bulayer [2], and fracture processes of semicrystalline polymers [3]. This
year, we studied the ion distribution on acharged surface as a model of biomembranes
by coarse-grained molecular dynamics simulation using GROMACS [4] on systems B and C.
The electrostatics in soft matter systems isrelevant to its self-assembly structures anddynamics, and therefore, essential to be understood. In particular, ion distributions ona charged surface as a model of biomembranes containing charged phospholipids havebeen extensively studied. Although theoreticalstudies have revealed the ion distributions perpendicular to the homogeneous charged surface, the ion distributions parallel to the heterogeneous charged surface and their dynamics
are still unclear. Our coarse-grained molecular dynamics simulations show the dynamicalchanges of ion distribution proﬁles on the lipid
membrane during phase separation. Figure 1
shows a typical snapshot of phase separationof neutral and negatively charged lipids with a
NaCl concentration of 100 mM, and the cation
distribution near the lipid membrane at 2.0 µs.
The cations following negatively charged lipids
in microsecond order indicate rapid redistribution of cations parallel to the membrane. We
compare the simulation results with the modiﬁed Poisson–Boltzmann theory and reveal theimportance of the ion-ion interactions in an
aqueous solution. These results will contributeto the future development of the theory.

Figure 1: Phase separation and cation distribution of the charged and neutral lipids mixture at 2.0 µs with a NaCl concentration of 100
mM. The snapshot of the top view of the membrane and cation density (number/nm3 ) near
the surface.

References[1] R. Wakabayashi, R. Imatani, M. Katsuya,
Y. Higuchi, H. Noguchi, N. Kamiya, and
M. Goto: Chem. Commun. 58 (2022) 585.
[2] Y. Higuchi, Y. Asano, T. Kuwahara, and
M. Hishida: Langmuir 37 (2021) 5329.
[3] Y. Higuchi: Phys. Rev. E 103 (2021)
042502.
[4] M. J. Abraham, T. Murtola, R. Schulz, S.
Páll, J. C. Smith, B. Hess, and E. Lindahl,
SoftwareX 1-2 (2015) 19.

219

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

Constructing a Dataset of Thermal and Dielectric Properties
of Polymers Using Molecular Dynamics Simulations
Yuta YOSHIMOTO
Department of Mechanical Engineering,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656
We have constructed a dataset of thermal

GHz, validating our MD simulation protocol.

and frequency-dependent dielectric properties

Currently, virtual polymers in the PI1M

of polymers using high-throughput molecular

database have been screened using machine

dynamics (MD) simulations. More than 500

learning-aided prediction models of CLTEs andpolymers are constructed from the repeating

dielectric constants to further extend the dataset.

units of previously synthesized polymers [1,2]and virtual polymers in the PI1M database [3]Referencesusing RadonPy [4], an open-source Python[1] L. Chen et al.: npj Comput. Mater. 6(2020) 61.

library. As shown in Fig. 1(a), the MDcalculated values of coefficients of linear[2] S. Otsuka et al.: Int. Conf. on Emerging

thermal expansion (CLTEs) show reasonable

Intelligent Data and Web Technologies

agreement with the available experimental(2011) 22.
[3] R. Ma and T. Luo: J. Chem. Inf. Model. 60

values. Additionally, as shown in Fig. 1(b) and(2020) 4684.

(c), the MD-calculated dielectric constants also
exhibitreasonableagreementwith[4] Y. Hayashi et al.: npj Comput. Mater. 8

the(2022) 222.

experimental counterparts at 10 MHz and 1

Fig. 1: Comparison of the MD-calculated values of coefficients of linear thermal expansion (a),dielectric constants at 10 MHz (b) and 1 GHz (c) with the experimental counterparts [1,2].

220

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

Critical Phenomena of Active Brownian Particles
with Uniaxial Self-propulsion
Hiroyoshi Nakano
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

Active matter is a distinct class of
nonequilibrium many-body systems that show
unique collective phenomena such as ﬂocking,motility-induced phase separation, giant number ﬂuctuations and active turbulence. Manyexperimental realizations have been reported,
ranging from swimming bacteria to synthetic
colloidal particles, which motivate extensive
investigations of properties of active matter.
Active Brownian particles (ABPs) are one of
the minimal model modeled on biologial motions [1]. The motion of the particles basicallyconsists of three steps: (i) translational random walk, (ii) random rotation of the particle’s orientation, and (iii) propelling movement
along the particle’s orientation. The ABPs undergo a phase separation even when the particles interact with each other via only excluded
volume interactions. Its origin is sought forthe self-propulsion. This type of phase separation is called motility-induced phase separation (MIPS) [2]. The speciality of the origin of
the MIPS motivates the studies of whether themacroscopic behaviors of the MIPS are diﬀerent from those of the phase separation inducedby the intrinsic attractive force, such as theequilibrium liquid-gas transition.
Related to this problem, we focused on the
ABPs with uniaxial self-propulsion [3]. Thismodel is the ABPs that are reﬁned so that theparticles prefer to be oriented along a certain
axis. We particularly examined the propertiesof critical phenomena by using the ﬁnite-size
scaling method.

The result is summarized in Fig. 1. Figure 1(a) presents typical snapshot of the phaseseparated state, and Fig. 1 (b) presents thephase diagram. The phase-separated state appears at the large self-propulsion force F0 . The
critical point is located at the triangle point in
Fig. 1 (b). The binder ratios U for four system sizes are plotted in Fig. 1 (c), which showsthat the binder ratios U cross at the uniquepoint. It implies that this model undergoes
the second-order phase transition. The criticalexponent is estimated as
β = 0.35(4) , νx = 0.65(6).

(1)

Surprisingly, this result is quite close to thecritical exponent of the three-dimensional Isingmodel even though this model is in twodimensions.
Our sophisticated simulations demonstrated
that this result can be understood from thesimilarities between the ABPs with uniaxial
self-propulsion and randomly driven lattice gas(RDLG). For the RDLG, the particles are subjected to the random driving forces along one
axis. It is well known that the universality
class of the RDLG is the same as that of thetwo-dimensional Ising model with the dipolar
long-range interaction, which is close to that of
the three-dimensional Ising model. Our simulations showed that the ABPs with uniaxial
self-propulsion exhibits the same macroscopic
behaviors as the RDLG in the homogeneousstate and the phase-separated state. Based on
this result, we consider that the critical point

221

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

Figure 1: (a) Typical snapshot of phase-separated state. (b) Phase diagram of the ABPs withuniaxial self-propulsion. (c) Binder ratios for four system sizes.
of the ABPs with uniaxial self-propulsion belongs to the same universality class as that of
the RDLG [4].
In summary, we performed the simulationsof the two-dimensional ABPs with uniaxial
self-propulsion. We showed that this model
undergoes the MIPS and that the critical behaviors is quite close to the three-dimensional
Ising model. Our result suggested that thecritical points presumably belongs to the same
universality class as that of the RDLG.

References[1] Y. Fily and M. C. Marchetti: Phys. Rev.
Let. 108, 235702 (2012).
[2] M. E. Cates and J. Tailleur: Annu. Rev.
Condens. Matter Phys. 6, 219-244 (2015).
[3] H. Nakano, and K. Adach: in preparation.
[4] K. Adachi, K. Takasan, and K.
Kawaguchi:
Phys. Rev. Research 4013194 (2022).

222

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

Numerical Study on Novel Spin Nematic PhaseTôru SAKAI1,2 , Hiroki NAKANO1 , Rito FURUCHI1 ,and Tokuro SHIMOKAWA3
1
Graduate School of Science, University of Hyogo,
Kouto, Kamigori, Hyogo 678-1297, Japan
2
National Institute for Quantum Science and Technology,
SPring-8, Kouto, Sayo, Hyogo 679-5148, Japan
3
Theory of Quantum Matter Unit,
Okinawa Institute of Science and Technology Graduate University,
Onna-son, Okinawa 904-0495, Japan1

Field-induced spin nematic the frustration. This mechanism also had been
Tomonaga-Luttinger liquid prposed for S = 1 antiferromagnetic chain
with the single-ion anisotropy[2].
of the S=1/2 spin laddersystem with anisotropic fer2
Spin Nematic Liquid ofromagnetic rung interaction
the S=1/2 Distorted Diamond Spin Chain in MagThe S = 1/2 quantum spin ladder system withthe anisotropic ferromagnetic exchange internetic Field
action on the rung under magnetic ﬁeld is investigated using the numerical diagonalizationand the density matrix renormalization group(DMRG) analyses[1]. It is found that thenematic-spin-correlation-dominant TomonagaLuttinger liquid appears in some high magnetic ﬁeld. It is included in the TomonagaLuttinger liquid phase where the two-magnon
bound state is realized.
For some suitable parameters, after the ﬁeld-induced phasetransition from this two-magnon TomonagaLuttinger liquid phase to the single-magnon
Tomonaga-Luttinger liquid one, the re-entrant
transition to the two-magnon phase would possibly occur, as the magnetization curves by
DMRG indicate. Several phase diagrams withrespect to the coupling anisotropy, the magnetization and the magnetic ﬁeld are presented.
The present result is a proposal of the candidate system which exhibits the spin nematic
phase without the biquadratic interaction or223

The spin nematic order has attracted a lot ofinterest in the ﬁeld of magnetism. It is a kindof multipole order of spins. The previous theoretical and numerical studies[3] predicted thatthe spin nematic order would be induced by
the frustration of the ferromagnetic and antiferromagnetic exchange interactions, or thebiquadratic interaction. The spin nematic order is characterized by the long-range four spincorrelation and the two-magnon bound state.
The previous numerical diagonalization andthe ﬁnite-size scaling study had indicated thata two-magnon bound state can occur in the
S=1 antiferromagnetic chain with the singleion anisotropy under magnetic ﬁeld[4].The recent analysis of the critical exponents of thespin correlation functions[5] suggested that
this two-magnon bound state includes the spin
nematic liquid phase, as well as the spin density wave (SDW) liquid one. The phase di-Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

agrams with respect to the anisotropy andthe magnetization were obtained by the numerical diagonalization of ﬁnite size clusters.
In the present study[6], the S=1/2 distorted
diamond spin chain is investigated using thenumerical diagonalization of ﬁnite-size clusters. This system is a typical frustrated system. The recently discovered candidate material of this system, K3 Cu3 AlO2 (SO4 )4 called
alumoklyuchevskite[7], includes the ferromagnetic interactions. Thus we investigated the
distorted diamond chain with the ferromagnetic interactions, as well as the coupling
anisotropy. As a result we found that a twomagnon bound state appears in the magnetization process. We also discus about the spin
nematic liquid behaviors.

References[1] T. Sakai, R. Nakanishi, T. Yamada, R.
Furuchi, H. Nakano, H. Kaneyasu, K.
Okamoto and T. Tonegawa, Phys. Rev. B
106 (2022) 064433.
[2] T. Sakai, H. Nakano, R. Furuchi and K.
Okamoto, J. Phys.: Conf. Ser. 2164 (2022)
012030.
[3] T. Hikihara et al., Phys. Rev. B 78 (2008)
144404.
[4] T. Sakai, Phys. Rev. B 58 (1998) 6268.
[5] T. Sakai, H. Nakano, K. Okamoto and R.
Furuchi, J. Phys.: Conf. Ser. 2164 (2022)
012030.
[6] T. Sakai, H. Nakano, R. Fruchi andK. Okamoto, AIP Advances 13 (2023)
015313.
[7] M. Fujihara et al., Sci. Rep. 7 (2017)
16785.

224

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

Magnetization Plateau with Spontaneous Symmetry
BreakingTôru SAKAI1,2 , Hiroki NAKANO1 , Rito FURUCHI1 ,and Tokuro SHIMOKAWA3
1
Graduate School of Science, University of Hyogo,
Kouto, Kamigori, Hyogo 678-1297, Japan
2
National Institute for Quantum Science and Technology,
SPring-8, Kouto, Sayo, Hyogo 679-5148, Japan
3
Theory of Quantum Matter Unit,
Okinawa Institute of Science and Technology Graduate University,
Onna-son, Okinawa 904-0495, Japan1

S=2 Quantum Spin Chain
Translational
Symme- 2
try Broken Magnetization
with the Biquadratic Exchange Interaction
Plateau of the S=2 Antiferromagnetic Chain with The S = 2 quantum spin chain with the singleAnisotropies
ion anisotropy D and the biquadratic exchange

The magnetization plateau of the S=2 antiferromagnetic chain with interaction and singleion anisotropies is investigated using the numerical diagonalization of ﬁnite-size clusters
and some size scaling analyses[1]. The previous
level spectroscopy analysis indicated that two
diﬀerent magnetization plateau phases appear
at half of the saturation magnetization. One isdue to the large-D mechanism and the other isdue to the Haldane one. In the present study
the phase diagram is extended to wider regionof the anisotropies. As a result we ﬁnd anotherhalf magnetization plateau phase, where thetranslational symmetry is spontaneously broken .

225

interaction is investigated using the numerical
diagonalization of ﬁnite-size clusters and thelevel spectroscopy analysis. It is found thatthe intermediate-D phase corresponding to thesymmetry protected topological (SPT) phaseappears in a wide region of the ground statephase diagram. We also obtain the phase diagram at the half of the saturation magnetization which includes the SPT plateau phase[2].

References[1] T. Yamada, R. Nakanishi, R. Furuchi, H. Nakano, H. Kaneyasu, K.
Okamoto, T. Tonegawa and T. Sakai,arXiv:2208.10983; to appear in JPS Conf.
Proc.
[2] T. Sakai, T. Yamada, R. Nakanishi, R.
Furuchi, H. Nakano, H. Kaneyasu, K.
Okamoto and T. Tonegawa, J. Phys. Soc.
Jpn. 91 (2022) 074702.

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

Simulation of Quantum Many-Body Systems by
Tensor Network and Sampling
Synge TODO
Department of Physics, University of Tokyo, Tokyo, 113-0033
Institute for Solid State Physics, University of Tokyo, Kashiwa, 277-8581

We developed novel numerical methods by
combining quantum Monte Carlo, tensor network algorithms, and other optimization techniques. We investigated topological quantum
phase transitions in strongly correlated manybody systems and optimization of quantum operations in quantum computing.
Tensor Renormalization Group and Markov
Chain Monte Carlo Methods
The real-space renormalization groupmethod using tensor networks has recently
been widely used as a numerical method
for many-body spin systems. The tensor
renormalization group method can eﬃciently
compute physical quantities for large classical/quantum systems. By extending thetensor network to a form in which tensors
are placed not only on vertices but also onlines connecting vertices, we have increased
the accuracy of the real-space renormalization
method by about 100 times compared to existing numerical methods that require a similar
amount of computation time. This method
applies to general tensor networks. We arealso working on a new sampling method thatcombines tensor network representation andMarkov chain Monte Carlo methods.
Tensor Renormalization Group for Lattice
Fermion Systems
Quantum many-body systems involving lattice fermions appear in various ﬁelds ofphysics. For example, lattice quantum chromodynamics is typical in particle physics, amany-body system describing the interaction

between quarks via gluons. Strongly correlated electron systems in condensed matterphysics are another example of a many-bodysystem consisting of lattice fermions. We develop a method to calculate path integralsof lattice fermion systems. We represent thepath integral of a lattice fermion system bya tensor network and perform approximate
network contraction using the tensor renormalization group to compute the path integral in the inﬁnite volume limit. We extended the bond-weighted tensor renormalization group (BTRG) to handle fermionic systems and performed empirical calculations using 2-dimensional Wilson fermions [1].
Worm Algorithm with Geometric Assignment
MethodIn Monte Carlo methods, there are cases
where there are constraints on the “states”
to be summed, and it is diﬃcult to sample
while satisfying the constraints (e.g., satisﬁability problems). The worm algorithm is amethod for eﬃcient sampling by expanding
the state space in such situations. A worm isa kink (point) that breaks a constraint, andthe idea of the worm algorithm is that theconstraint is violated once and then reconciled later. This Monte Carlo method introduces a constraint-breaking worm (kink), andsampling is performed by stochastically moving the worm. The conventional worm algorithm moves the worm almost wholly at random. We proposed a guideline for moving theworm for eﬃcient computation—moving the226

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

Figure 1: Finite-temperature phase diagram ofsquare-lattice quantum dimer model.
worm straight ahead as much as possible—and
conﬁrmed a 25-fold improvement in computational eﬃciency compared to the conventional
worm algorithm [2].
Finite-Temperature Phase Transition of the
Quantum Dimer Model
The quantum dimer model was proposed by
Rokhsar and Kivelson in 1988 as a low-energy
eﬀective model for frustrated magnetic materials. Although the Hamiltonian of the quantum
dimer model has no negative sign problem,
Monte Carlo simulations have been complicated because of the strong geometric restriction on the dimer conﬁguration. We developed
an eﬃcient quantum Monte Carlo method thatcombines a cluster update method based on
the stochastic series expansion (SSE) with afast operator search method using priority
queues. By choosing the appropriate order parameter and ﬁnite scaling, we obtained a ﬁnitetemperature phase diagram of the quantitative
phase diagram of the quantum dimer model.
The critical exponent at the critical point is
consistent with the Kosterlitz-Thouless transition of the classical dimer model. The slope
in the classical limit is also compatible withthe critical temperature of the classical model.
The slope of the phase boundary decreases
as the quantum ﬂuctuations increase, suggest-

ing that the critical temperature becomes zeroat the Rohksar-Kivelson (RK) point, i.e., thecolumnar dimer phase at absolute zero extends
to the RK point without intermediate phases.
Continuous-Space Path-Integral Monte Carlo
The continuous-space path-integral Monte
Carlo method is the only method that cannumerically calculate physical quantities suchas the superﬂuid density and condensation
rate of interacting Bose particle systems. The
worm algorithm allows for simulations in the
grand canonical ensemble of superﬂuid helium
in systems of several hundred particles. However, implementing the algorithm is complex,
involving seven diﬀerent update procedures.
We developed more eﬃcient and simpler update methods in the canonical ensemble. Advanced Markov-chain Monte Carlo methods,
such as Hamiltonian Monte Carlo and event
chain Monte Carlo, have been used to reduce
the autocorrelation time of the simulation further. In addition, the Wang-Landau algorithm
allows the hyperparameters to be automatically adjusted during the simulation, whereas
the conventional worm algorithm has many hyperparameters that signiﬁcantly aﬀect the performance of the algorithm.

References[1] Shinichiro Akiyama, Bond-weighting
method for the Grassmann tensor renormalization group, J. High Energy Phys.
11, 030 (2022).
[2] Hidemaro Suwa, Lifted directed-worm
algorithm, Phys. Rev. E 106, 055306(2022).

227

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

Novel phenomena in frustrated magnets
Tsuyoshi OKUBO
Institute for Physics of Intelligence, The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 133-0033

Frustrated interactions often induce novel
phenomena in spin systems. One striking effect when we consider the quantum ﬂuctuation may be the stabilization of quantum spin
liquid without having long-range magnetic order [1, 2]. The ground state of the Kitaev
model on the honeycomb lattice is one of thefundamental quantum spin liquids in the twodimension[3, 4]. In the Kitaev model, the interactions are Ising type, Siγ Sjγ , and the spincomponent γ = x, y, z is determined from thedirection to the neighboring site. Another example is the Kagome lattice Heisenberg model[1, 2]. It has been widely discussed on thepossibility of the gapped Z2 spin liquid or thegapless U (1) spin liquid in the Kagome lattice
Heisenberg antiferromagnets, e.g., see Ref. [5]and references therein.
In this year’s project, we investigated ﬁnite temperature properties of frustrated spinmodels in the two-dimension, particularly concerning the thermal Hall conductivity. We
mainly considered the Kagome lattice Heisenberg model with Dzyaloshinskii-Moriya (DM)interaction and compared its thermal Hall conductivity with the Kitaev model. To numerically calculate the thermal Hall conductivities
at a ﬁnite temperature, we consider severalﬁnite-size clusters with open edges and represent their density matrix as a matrix productoperator (MPO).
We calculate the MPO representing a density matrix at a given temperature by the exponential tensor renormalization grope (XTRG)
algorithm [6]. In XTRG, we calculate the density matrix at an inverse temperature β, ρ(β),

through the relationship ρ(β) = ρ(β/2)ρ(β/2)
starting from the initial condition ρ(β0 ) =exp(−β0 H) ≃ 1−β0 H for the suﬃciently small
β0 , where H is the Hamiltonian of the system.
Once we have an MPO representation of a density matrix at β/2, we can approximate ρ(β)
through a standard tensor network compression algorithm to multiply two MPOs. Whenthe bond dimension of the MPO is D, boththe computation cost and the memory consumption of the XTRG algorithm scale O(D4 ).
Although they are higher than the costs forthe ground state, which scale O(D3 ), we canstill calculate ﬁnite temperature properties by
XTRG.
By the XTRG algorithm, we successfully
calculated the temperature dependence of theenergy current localized at the edge of thecylindrical system for the Kagome lattice
Heisenberg model with the DM interaction under a magnetic ﬁeld. We can calculate thethermal Hall conductivity by taking the temperature derivative of the energy current. In
this system, the direction of the thermal Hall
current depends on the sign of the DM interaction. Interestingly, in the XTRG calculation,
its sign is opposite to the one calculated by
Schwinger Boson mean-ﬁeld theory [7]. We
also found that similar sign inversion in thethermal Hall current is also observed in thecase of the Kitaev model at high temperature:
the sign is opposite to the expectation fromthe linear spin wave theory, although the signat the low-temperature value is consistent withthe expectation from the analysis of the Kitaev
spin liquid. It will be a future problem to un-

228

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

derstand the origin of the diﬀerence in signs inhigh-temperature regions.

References[1] L. Savary and L. Balents, Rep. Prog.Phys.
80 (2017) 016502.
[2] Y. Zhou, K. Kanoda, and T.K. Ng,
Rev. Mod. Phys. 89 (2017) 025003.
[3] A. Kitaev, Ann. Phys. 321 (2006) 2.
[4] Y. Motome and J. Nasu, J. Phys. Soc.Jpn.
89 (2020) 012002.
[5] Y.C. He, M. P. Zaletel, M. Oshikawa,and F. Pollman, Phys. Rev. X, 7 (2017)031020.
[6] H. Li, D.-W. Qu, H.-K. Zhang, et al.
Phys. Rev. Research 2 (2020) 043015.
[7] H. Doki, M. Akazawa, H.-Y. Lee,J. H. Han, K. Sugii, M. Shimozawa,N. Kawashima, M. Oda, H. Yoshida,
and M. Yamashita, Phys. Rev. Lett. 121(2018) 097203.

229

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

Theoretical design of novel proteins for medical and industrial
applications
Nao SATO1, Hirotaro SHIMAMURA2, Rina AOYAMA1, Mizuki TERANISHI1,
Koji OOKA3, Runjing LIU1, Kaho KOTANI1, Masataka YOSHIMURA1,
Yuuki HAYASHI1,4, Shunji SUETAKA1, and Munehito ARAI1,2
1
2
4

Department of Life Sciences, Graduate School of Arts and Sciences,

Department of Physics, Graduate School of Science, 3College of Arts and Sciences, andEnvironmental Science Center, The University of Tokyo, Komaba, Meguro, Tokyo 153-8902
Proteins are widely used in medicine and

12 designed peptides, we experimentally foundindustry as therapeutic agents for variousthat a mutant of the MLL TAD peptide had thediseasesindustrial

highest binding affinity for KIX among all

production. However, these proteins have beenpreviously reported inhibitors targeting thedeveloped

experiments,

MLL site of KIX [2]. Thus, our approach may

which require much time and cost. Therefore,

be useful for the rational design of helicalthere is an urgent need to develop efficient andpeptides that inhibit protein-protein interactions

theoreticalinvolved in various diseases.

and

as

throughmethods

enzymes
extensive

for

for

designing

novel

proteins. Here, we theoretically designed a

The programmed death-1 (PD-1) protein on

variety of novel proteins using the Rosetta

the surface of T cells inhibits excessive

software suite [1] and experimentally verified

inflammatory activity of T cells by binding to

the validity of our theoretical designs. In

programmed death ligand-1 (PD-L1) expressed

addition, we constructed a theory that canon antigen-presenting cells. However, reduced

predict the folding mechanism of proteins, toPD-1 function can cause abnormally activatedelucidate the principles of protein design.

T cells to damage normal cells, leading to

The KIX domain of CREB-binding protein

autoimmune diseases. To prevent the T-cellis associated with leukemia, cancer and various

activation by turning on the immune checkpoint,

viral diseases and has attracted considerable

we theoretically designed PD-1 agonists that

attention in drug discovery. We rationallytightly bind the PD-L1 binding site of PD-1.

designed a KIX inhibitor by performing

Using

theoretical saturation mutagenesis to search forsaturation mutagenesis on the fragment of

the peptides expected to bind KIX more tightly

human PD-L1 and selected the mutants that arethan the original partner proteins. Among thepredicted to have high affinity for PD-1.

230

Rosetta,we

performed

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

Experimental verification showed that we have

we are developing a universal method that can

successfully designed the PD-1 binders withgenerally improve the activity of various

higher affinity than the wild-type PD-L1, which

enzymes. We have applied this method to

may be useful as inhibitors of autoimmune

enzymes such as dihydrofolate reductase anddiseases.

aldehyde deformylating oxygenase. Verification

The

novel

coronavirus(SARS-CoV-2)of this method is underway.

infects human cells when its receptor-binding
domain

(RBD)

ACE2.

proteins fold into specific three-dimensional

Antibodies that inhibit this binding have been

structures is necessary to understand the designused as antiviral drugs. However, existing

principles of proteins. To solve this problem,

drugs may become less effective as new

we are developing a statistical mechanics

variants of SARS-CoV-2 with mutations in thetheory that can draw the folding free-energy

RBD emerge. Here, we theoretically modified

landscapes of any proteins [3]. The program isexisting antibody drugs to bind tightly to thewritten

RBD

successfully

mutant

binds

of

to

human

Elucidating the mechanisms by which

SARS-CoV-2

and

in

C++

and

calculated

Python.
the

We havefreeenergyexperimentally showed that the modified

landscapes of folding reactions for all-α and all-antibodies were able to bind tightly to the RBD

β proteins in agreement with experiments. Our

mutant. This approach is expected to be useful

method can be used to draw the free energy

in the future as new viral variants emerge.

landscapes of proteins associated with folding,Formationofaternary

complexbybinding, and catalysis, and to develop new

interleukin 33, its receptor ST2, and IL-1

protein design methods based on the prediction

receptor accessory protein (IL-1RAcP) induces

of protein dynamics.

allergic diseases. To inhibit the complex
formation, we have theoretically designed small

References

proteins that tightly bind ST2 but inhibit the[1] J. K. Leman, et al.: Nature Methods, 17,binding of IL-1RAcP. Experiments to verify the

665 (2020).

designs are ongoing.

[2] N. Sato, S. Suetaka, Y. Hayashi, and M.
Arai: Scientific Reports, 13, 6330 (2023).

Enzymes are widely used in the industrial
production of useful substances. Therefore, it is[3] K. Ooka, R. Liu, and M. Arai: Molecules,

necessary to develop efficient methods to27, 4460 (2022).

improve enzymes. Using the Rosetta software,

231

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

Multiscale Flow Simulations on Complex Fluids
Undergoing Phase Transition
Yohei MORII and Toshihiro KAWAKATSU
Department of Physics, Tohoku University, Sendai, Miyagi 980-8578
The knowledge on the properties of complex

irection to the inflow velocity divided by theflows, such as a viscoelastic polymer flow, aroundinflow speed.

an obstacle plays a crucial role in designing flow

Figure 1 shows a snapshot of thechannels in machines. In the present study, we

molecular orientations

analyze a flow of a dilute polymer solution

stretching of the FENE dumbbells around

passing through a moving wing using thea fixed Joukowski wing. When the wing

multiscale simulation platform named MSSP,

starts to rotate slowly changing its angle

where smoothed particle hydrodynamics (SPH)to the inflow direction as 0  45  0 ,method is coupled with microscopic particle

the value of CL changes as shown in Fig.2,

simulations or phenomenological constitutive

where the data are shown for Newtonian (red),

relations.

Maxwellian(green)

and the bond

and

FENE

dumbbell(blue) models, respectively.

We modelled the wing using Joukowski wingWe

as shown in Fig.1, and the polymer chains arecan confirm that the oscillation in CL ismodelled either by a Maxwellian constitutive

suppressed

relation (infinitely extensible) or FENE dumbbells

polymers(Maxwellian) while such a suppression(finite extensible). Simulations were performed at

is less obvious for FENE dumbbells.

byaddingstretchable

Ohtaka, where 2,560,000 SPH particles with5,000,000,000 dumbbell particles embedded in

References

these SPH particles were used. From the[1] Y.Morii and T.Kawakatsu, Phys. Fluids, 33,

simulation data, we evaluated the rift coefficient

093106 (2021).

CL which is defined as the magnitude of theforce exerted on the wing in the perpendicularFig.1 Flow pattern of a polymer solution
around a Joukowski wing. Shown are thedirection (bars) and the extension (colors) ofpolymers chains.

Fig.2 CL as a function of the angle of the wing.

232

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

Simulations of quantum-classical-hybrid algorithms forsensor materials with considering noise
Wataru Mizukami
Center for Quantum Information and Quantum Biology, Osaka University, Osaka 560-0043,
Japan
In the fiscal year, a domestically produced

classical computers has become crucial. In thisquantum computer was made publicly available

study, we have performed the development andin the cloud, and IBM announced a quantum

validation

computer with 443 quantum bits. Quantum

algorithms using the ISSP’s supercomputer andcomputers with thousands of quantum bits area quantum circuit simulator rather than an actual

expected to appear in the near future, indicating

quantum computer.

ofquantum-classicalhybridsignificant progress in the development ofThis year, following up on last year's work,

quantum computing. A promising application

we focused on accelerating the simulation of thefor quantum computers is the simulation ofvariational quantum eigensolver (VQE), a

quantum many-body systems. In the field ofrepresentative

materials development, it is expected thatalgorithm. We worked on speeding up VQE

simulations will become more accurate thancomputations using GPUs and implemented

current methods. Research into such simulation

VQE with multi-GPU support.

quantum-classicalhybridMoreover, we have been developing atechniques is gaining momentum.
Our aim is to apply quantum computing tomethod that combines the auxiliary-field

the development of materials for use in sensors.

quantum Monte Carlo (AFQMC) technique with

However, performing practical computations on

VQE to compensate for degrees of freedom thatrealistic quantum chemistry problems with acannot be handled by VQE alone. While alarge number of degrees of freedom using only

combination

quantum computers with limit number of qubits

computation was proposed by Google in 2021,

is not feasible. In addition, current quantum

we have designed and implemented an algorithm

computers have short coherence times, making

that results in a more practical quantum circuit

it difficult to perform useful computations using

than their implementation. We tested this

only quantum computers. As a result, thealgorithm on simple systems, such as the OH

development

molecule,ofquantum-classicalhybridandperformance.

algorithms that utilise both quantum and

233

of

AFQMC

confirmed

and

its

quantum

satisfactory

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

Tensor renormalization-group study of random spin
systems
Koji Hukushima
Department of Basic Science, University of Tokyo
3-8-1 Komaba, Meguro-ku, Tokyo 153-8902
Statistical mechanical studies of spin glasses
with both randomness and frustration havebeen conducted for a long time. It is an important issue whether the picture of phase transitions established in the mean-ﬁled spin-glassmodel can be applied to ﬁnite-dimensional
spin-glass models. It is quite diﬃcult to study
the ﬁnite-dimensional spin glass models using analytical methods. Therefore, the main
results for these models have been obtained
mainly by numerical studies, such as Markovchain Monte Carlo (MCMC) methods. However, the long relaxation times of the MCMC
methods for spin-glass systems are a serious
barrier to studying large sizes. Recently, thetensor renormalization group (TRG) method,
numerical calculation methods using tensor
networks have attracted much attention.
In this project, we employ the higherorder tensor renormalization group (HOTRG)
method, one of the TRG-type methods, tostudy the two-dimensional Ising spin-glassmodel. Previous studies have suggested thatthe spin-glass phase transition does not occur
at ﬁnite temperatures in this system when the
distribution of interactions is close to symmetric and highly disordered[1]. We investigate thephase transition of the system near the Nishimori line. Then, the stiﬀness of the system iscalculated for sizes up to the linear dimension
L = 256 using the HOTRG method, which calculates the free energy under diﬀerent boundary conditions. In the calculations at higher
temperatures than the Nishimori line, bothboundary conditions lead to the same transition point and values close to the critical exponent of the pure Ising universality class, thelatter of which is consistent with the picture
dominated by the pure Ising ﬁxed point there,
as pointed out in many previous studies.

In contrast to the phase transitions above
and on the Nishimori line, the transition points
analyzed here have diﬀerent values for thetwist and replica boundary conditions. In
other words, the paramagnetic/ferromagnetic
and paramagnetic/non-paramagnetic transition points are diﬀerent from each other, suggesting the existence of some non-paramagnetic
phase between the paramagnetic and ferromagnetic phases at ﬁnite temperatures in the twodimensional EA Ising model. This may correspond to the spin-glass phase or the previously
proposed random antiphase state. The valueof this transition point at low temperature isroughly consistent with the Nishimori ﬁxed
point, suggesting that the boundary between
the paramagnetic and “spin-glass” phases descends almost vertically from the Nishimori
ﬁxed point toward lower temperatures. Meanwhile, the “spin-glass”/ferromagnetic phaseboundary shows a re-entrant transition.
We also studied the critical phenomena of
the J1 -J2 Ising model by varying the parameter g using HOTRG[2]. This system involves
a competition between the ferromagnetic interaction J1 and antiferromagnetic interaction J2 .
Furthermore, weak ﬁrst-order and second-ordertransitions are observed near the ratio g =J2 /|J1 | = 1/2. Our results, based on HOTRG
calculations for signiﬁcantly larger sizes, indicate that the region of the ﬁrst-order transition
is marginally narrower than that in previous
studies. Moreover, the universality class of thesecond-order transition connected to the transition line is not necessarily fully consistent with
the Ashkin–Teller class considered earlier.
[1] K. Yoshiyama, K. Hukushima: in preparation.
[2] K. Yoshiyama,
K. Hukushima:
arXiv:2303.07733

234

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

Molecular dynamics study of phase separationinduced by binding molecules
Koji Hukushima
Department of Basic Science, University of Tokyo
3-8-1 Komaba, Meguro-ku, Tokyo 153-8902
Chromatin is highly folded into a hierarchicalstructure in the nucleus while interacting with
nuclear molecules. In addition, regions of active and repressed gene expression are formedon chromatin. These epigenetic regions are established and maintained by writer moleculesthat write and reader molecules that recognize chemical modiﬁcations on chromatin. The
relationship between the one-dimensional epigenetic pattern on chromatin and the threedimensional chromatin structure is an interesting question that has not yet been elusidated.
The SBS (Strings and Binders Switch)model[1] is a coarse-grained polymer modelthat takes into account the phase separation
between chromatin and nuclear molecules thatbind to chromatin, and to some extent recapitulates the experimental results such as HiC andFISH. In this model, chromatin forms clusters
due to the eﬀective attraction between chromatin segments via molecules that bind multivalently to chromatin[1,2].
In this project, we consider two types ofbinding molecules in the SBS model that promote or repress transcription as playing theroles of reader and writer of epigenetic modiﬁcations. When polymer segments are boundto these binding molecules, the state stochastically changes depending on the type of thebinding molecule.
The binding molecules
switch between two states, one in which they
can bind chromatin and one in which they cannot, at a constant rate, as in the previous
study[3]. In this model, a one-dimensional pattern of epigenetic modiﬁcation appears on thepolymer, accompanied by a three-dimensional
structure.
In this study, we performed coarse-grained
molecular dynamics simulations on the SBS
model with two types of biding molecules.

From the time-series data of the coordination
dynamics obtained by the molecular dynamics
simulations, we devised a method to extract
the time evolution of cluster formation, ﬁssion,coalescence, and annihilation. With the helpof the method, it is revealed that the clusters
move dynamically along the polymer, repeating coalescence and ﬁssion, and that there are
no independent clusters. The evaluation of the
lifetime of the clusters also suggested the existence of a phase that goes to inﬁnity below
a certain nonzero switching rate, despite a ﬁnite switching rate. This can be regarded as adynamical percolation transition.
Our simulations also showed that the sequence of clusters of alternating types andthe hierarchical structure of these clusters on
the chromatin polymer emerge with the two
types of binding molecules. When the binding molecules do not switch states, the polymer takes on a sausage-like shape with hierarchical structure, and as the switching rate increases, the structure becomes more dynamic
and the clusters become more independent ofeach other. This provides insight into themechanisms underlying the dynamics and stability of one-dimensional epigenetic patterns
and three-dimensional cluster structures onchromatin in the SBS model.
The present work has been done in collaboration with Ryo Nakanishi.
[1] Mariano Barbieri et al., PNAS, 2012 109(40) 16173-16178.
[2] Andrea M. Chiariello et al., Sci. Rep., 2016
29775(6).
[3] Chris A. Brackley et al., Biophys. J., 2017112(6) 1085-1093.

235

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

Numerical simulations combined with machinelearning to search for magnetic skyrmion in
frustrated magnets
Satoru Hayami
Graduate School of Science, Hokkaido University, Sapporo 060-0810, Japan

A magnetic skyrmion, which is characterized by a topologically-nontrivial swirling spin
texture, has attracted considerable attention,
since it leads to unconventional transport phenomena, such as the topological Hall eﬀect. In
order to realize new electronic and spintronic
devices based on the skyrmion, it is highly desired to understand when and how such topological spin textures are stabilized from the microscopic viewpoint. In the project with number 2022-Ca-0002, we have numerically investigated the possibility of magnetic skyrmions byfocusing on magnetic frustration mechanisms
in localized and itinerant systems. We havepresented the main results this year below.
(i) Skyrmion crystals in trilayer systems:
We have theoretically investigated a new stabilization mechanism of the skyrmion crystals ina trilayer system [1]. By focusing on the roleof the layer-dependent Dzyaloshinskii-Moriya
interaction and performing the simulated annealing for the spin model, we obtained various types of skyrmion crystal phases with different skyrmion number in a low-temperaturephase diagram. We found the skyrmion crystal phases characteristics of the layered system:
the twisted surface skyrmion crystal, antiskyrmion crystal, and high-topological-number
skyrmion crystal phases.
(ii) Rectangular and square skyrmion crystals in centrosymmetric tetragonal magnets:
We considered the instability toward the rectangular skyrmion crystal without fourfold
rotational symmetry on a centrosymmetricsquare lattice with magnetic anisotropy [2].
The model includes the eﬀect of frustrated exchange interactions and the easy-axis singleion anisotropy. We have shown that a competition of the interactions in momentum space
triggers the formation of both rectangular andsquare skyrmion crystals. We have also shownthat a phase transition between them occurs
while varying the amplitude of the magnetic
ﬁeld.
(iii) Skyrmion crystals in a nonsymmorphic lattice structure: We have numerically investigated the stabilization mechanism of theskyrmion crystals in a nonsymmorphic latticesystem with screw symmetry [3]. By performing unbiased numerical simulations basedon the simulated annealing for a layered spinmodel to be characterized as a nonsymmorphic system with a threefold screw axis, werevealed that the skyrmion crystal is stabilized in both zero and nonzero magnetic
ﬁelds by the interplay between the ferromagnetic interlayer exchange interaction and thelayer- and momentum-dependent anisotropic
exchange interactions. We also found two new
types of skyrmion crystals characteristics of athreefold-screw system. Our result indicateda possibility to realize the skyrmion crystal in
nonsymmorphic systems without a solely rotational symmetry within the same plane, whichleads to extending the scope of the skyrmion
crystal-hosting materials.
(iv) Skyrmion crystals in a ferromagnetic/antiferromagnetic bilayer: We proposed236

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyoa new stabilization mechanism of the skyrmion
crystal in a bilayer system, which is composedof ferromagnetic and antiferromagnetic layers
but free from the Dzyaloshinskii-Moriya interaction [4]. By performing the variational calculation and Monte Carlo simulation for the frustrated spin model, we found that the skyrmion
crystal is robustly stabilized from zero to ﬁnite temperatures for the strong interlayer coupling. The results provided a new scenario torealize the skyrmion crystal in multi-layer systems including the domain structure and heterostructure by making use of the layer degreeof freedom.
(v) Skyrmion crystal in the Hubbard model:
We investigated a realization of the skyrmion
crystal in the centrosymmetric triangularlattice Hubbard model without localized
spins [5]. On a basis of the self-consistentmean-ﬁeld calculations, we found two types ofskyrmion crystals with spatially nonuniform
charge modulations in the ground state at azero magnetic ﬁeld. We have also shown thatan unconventional triple-Q vortex phase without the net spin scalar chirality emerges. Furthermore, we found that the magnetic ﬁeld induces the topological phase transition from the3Q vortex phase to another skyrmion crystal
phase. We clariﬁed that the microscopic originof these phases is the site-dependent moment
reduction due to the kinetic motion of electrons
in the Hubbard model, which is qualitatively
diﬀerent from the situation in the Kondo lattice model and other spin models.
(vi) Various topological spin textures in an
itinerant honeycomb magnet: We have theoretically revealed new types of topological spin
textures consisting of the magnetic skyrmion
by focusing on the sublattice degree of freedom in the lattice structure [6]. We found ﬁve
types of skyrmion crystals with diﬀerent topological properties from the conventional one in
an itinerant honeycomb magnet based on thesimulated annealing. We have shown that thesublattice degree of freedom in the itinerantmagnet enables us to realize the unobserved
skyrmion crystals, such as the antiferro-typeskyrmion crystal and ferri-type skyrmion crystal.
(vii) Magnetic ﬁeld–temperature phase diagrams for centrosymmetric tetragonal magnets: We systematically constructed the magnetic ﬁeld–temperature phase diagrams of theeﬀective spin model for the itinerant electron
model on a square lattice on the basis of an efﬁcient steepest descent method with a small
computational cost [7]. We found two distinct mechanisms to stabilize the square-typeskyrmion crystal from zero to ﬁnite temperatures: One is the positive biquadratic interaction and the other is the high-harmonic
wave-vector interaction. By showing the phase
diagrams for diﬀerent sets of model parameters, we discussed the stability tendency of theskyrmion crystal and other multiple-Q states.

References[1] S. Hayami, Phys. Rev. B 105, 184426(2022)[2] S. Hayami, Phys. Rev. B 105, 174437(2022).
[3] S. Hayami, Phys. Rev. B 105, 224411(2022).
[4] K. Okigami, R. Yambe, and S. Hayami, J.
Phys. Soc. Jpn. 91, 103701 (2022).
[5] K. Kobayashi and S. Hayami, Phys. Rev.
B 106, L140406 (2022).
[6] R. Yambe and S. Hayami, Phys. Rev. B
107, 014417 (2023).
[7] S. Hayami and Y. Kato, J. Magn. Magn.
Mater. 571, 170547 (2023).

237

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

Universality classes of 3D ±J Ising model
Y. Ozeki and Y. Terasawa
Graduate School of Informatics and Engineering, The University of Electro-Communications

We investigate the universality classes of the ±J
Ising model in three dimensions. Together with theferromagnetic (FM) classes obtained previously [1],the spin glass (SG) transition is analyzed by the useof dynamical scaling analysis for nonequilibrium relaxation (NER) process [2]. the p-dependence ofSG transition temperature TSG(p) and critical exponents γ, ν, z are obtained.
To estimate the dynamical exponent z precisely, an extrapolation scheme is introduced forthe asymptotic value of logarithmic derivative ofthe dynamical Binder parameter at the transition
point. To do this, we performed Monte-Carlo simulation for multi-replica systems. Each replica system has the same bonds Jij , while the initial configuration of relaxation is distinct. The replica overlap q AB (t) between two replicas A and B at timet, is defined by:
1 X (A)(B)
q AB (t) =S (t)Si (t),
(1)
N i i
(A)

where Si is the i-th spin variable in replica A attime t. Te SG Binder parameter is obtained by!
1
q (4) (t)gSG (t, T ) =3 − (2),(2)
2{q (t)}2
where
q

(k)

2
(t) =
n(n − 1)

"

X
A>B

{q

AB

(t)}

k

#

.

(3)

J

When T = TSG , this function is expected to diverge
asN × gSG (t, T ) ∼ td/z
(4)where d is the spacial dimension. This relation allows us to estimate z independently.
The SG region in 0.5 ≤ p ≤ pmc is analyzed. To
estimate TSG and zν, we perform the dynamicalscaling analysis for SG susceptibility χSG (t) withthe linear size L = 61, 512 independent samples
and 64 replicas are used for averaging. To estimatez, we calculate gSG (t) with the linear size L = 31,
460800 independent samples and 64 replicas areused for averaging.
In Figs. 1, the resulting p-dependence of SG critical exponent ν on the PM-SG boundary are shown

238

Figure 1: Result of estimated ν.
together with the previously obtained FM exponents on the PM-FM boundary [1]. Since both the
FM and SG transitions occur at the multi-criticalpoint (MCP), and FM and SG exponents become
identical with each other [3], one may survey an
exponent in whole region including SG and FM
phases. As for the exponent ν, there exist three
FM universality classes associated with the nonrandom FM fixed point (p = 1), the random fixed
point (pmc < p < 1), and the MCP (p = pmc ) [1].
On the other hand, the present result indicates asingle universality class for PM-SG transition in
0.5 ≤ p < pmc , and ν takes a value ∼ 1.55.

References[1] Y. Ozeki and N. Ito, J. Phys.: Math. Theor.
40 R149 (2007).
[2] Y. Echinaka and Y. Ozeki, Phys. Rev. E94
043312 (2016).
[3] Y. Ozeki and N. Ito, J. Phys.A: Math. Gen.
31 5451 (1998).

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

Localization-delocalization transition in
non-Hermitian random system
TOMI OHTSUKI1
TOHRU KAWARABAYASHI2
KEITH SLEVIN3
1) Dept. Phys., Sophia University, Chiyoda-ku, Tokyo 102-8554, Japan
2) Dept. Phys., Toho University, Miyama 2-2-1, Funabashi 274-8510, Japan
3) Dept. Phys., Osaka University, Toyonaka, Osaka 560-0043, Japan

Random quantum systems lack translational/rotational symmetries but can exhibit
time reversal, spin rotation, chiral, andparticle-hole symmetries. These symmetries
classify random quantum systems into 10 distinct symmetry classes. The Wigner-Dyson
classes discovered in the 1950s and 1960s haveundergone signiﬁcant research. Recent experimental ﬁndings regarding topological insulators and Weyl semimetals have inspired researchers to explore random systems outsideof the Wigner-Dyson classes. These newly
discovered 7 classes [3 chiral classes (GadeWegner classes) and 4 Bogoliubov-de Gennes
classes] exhibit discrete chiral or particle-holesymmetries.
We have studied a series of Hamiltonians
that pertain to the newly discovered quantum
classes. We then introduced non-Hermiticity
to the Hamiltonian, and we were able to extend the classiﬁcation of systems to encompass 38 distinct classes. Through our research, we demonstrated that the inclusion of
non-Hermitian elements alters the universality
class [1] and that the widely utilized transfer
matrix method employed in Anderson transitions also applies to non-Hermitian systems [2].
Additionally, our investigations revealed thatthe critical behaviors of the 38 non-Hermitianclasses can be mapped to universality classes
for the 10 Hermitian symmetry classes [3] (seeFig. 1).
Condensed matter physics has seen a signiﬁcant increase in the utilization of machinelearning techniques. Our research involves

using convolutional neural networks (CNN)
to analyze wave functions and study various
quantum phase transitions. Speciﬁcally, we
employ this method to classify wave functionsin quasi-periodic systems[4]. Moreover, we utilize the neural network to interpret the quantum ﬁngerprints in mesoscopic systems [5].
References
1. X. Luo, T. Ohtsuki, R. Shindou, Physical Review Letters 126, 090402 (2021).
2. X. Luo, T. Ohtsuki, and R. Shindou, PhysicalReview B 104, 104203 (2021).
3. X. Luo, Z. Xiao, K. Kawabata, T. Ohtsuki, and
R. Shindou, Physical Review Research 4, L022035(2022).
4. X. Luo, T. Ohtsuki, Physical Review B 106,
104205 (2022).
5. S. Daimon, K. Tsunekawa, S. Kawakami, T.
Kikkawa, R. Ramos, K. Oyanagi, T. Ohtsuki, E.
Saitoh, Nature Communications 13, 3160 (2022).

Figure 1: Correspondence of the Andersontransitions in non-Hermitian and Hermitiansystems.

239

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

Slow spin dynamics near phase boundary between
diﬀerent ordered phases
Hidemaro SUWA
Department of Physics, University of Tokyo, Tokyo 113-0033

The metastable state is ubiquitous in ourdaily lives and in condensed matter physics.
A typical example of the metastable state isdiamond, a solid form of carbon. Graphite is
the most stable phase at room temperature,
but diamond is also stable as a metastablestate, indicating that the time scale of relaxation is much longer than that of experiments.
The metastable state is essential to not only
statistical mechanics but also several functional materials. For example, the use of themetastable state enables higher-density storage and broader applications in phase-change
memory. Recently, a skyrmion-lattice phase of
MnSi has been stabilized by rapid cooling as ametastable state in a wider region of the phase
diagram.
What is the condition for a metastable state
to stay for a long time, or in what situation do
the relaxation dynamics become signiﬁcantly
slow? We can ﬁnd some answers in glass systems. Slow dynamics typically emerge near a
eutectic point, a triple point between a liquid,a crystal, and another crystal phase. In the dynamics of particle systems, the particle number
is usually conserved in time evolution. On theother hand, the dynamics of spin systems are
typically non-conserved: magnetization is notconserved in time evolution. The question weanswer through this project is where slow dynamics appear in the Ising model, a prototype
model of statistical mechanics.
In analogy to the eutectic point, we derived
one of the simplest Ising models showing a“Y-shape” ﬁrst-order phase boundary (Fig. 1):

the J1 -J2 -J3 Ising model with a four-spin interaction. The four-spin interaction could be
induced by integrating out other degrees of
freedom, mimicking an eﬀect of interactions
with charge, phonons, etc. We studied thequench dynamics of the Ising model using theMonte Carlo method on ISSP System B in theclass C project (ID: 2022-Ca-0100). The initialstate is a random spin conﬁguration, and thetime evolution is simulated by Glauber dynamics. The relaxation time becomes signiﬁcantly
longer near the phase boundary. We also studied snapshot conﬁgurations and the structure
factor as a function of time. Near the phaseboundary, clusters of the metastable state remain for a quite long time and prevent vortices
from recombining and annihilating. Our simulation is consistent with experimentally observed slow dynamics in vanadium oxides. The
result of this project sheds new light on thestability of the metastable state and slow dynamics emerging near a phase boundary.

Figure 1: Schematic Y-shape phase diagram of
the Ising model. The phase boundaries are of
the ﬁrst order.

240

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

Polymer Dynamics under Elongational Flow
Takahiro MURASHIMA
Department of Physics, Tohoku University,
Aramaki-aza-Aoba, Aoba-ku, Sendai, Miyagi 980-8578
Fundamental studies of polymer dynamics
under elongational flows are important for abetter understanding of complex flow in themolding process. We have developed a uniform
elongational flow method [1] that can handle
very

large

deformations

underarbitrary(uniaxial, planar, and biaxial) elongationalflows. This code works on LAMMPS [2],which brings the benefits of massive parallelFig. 1: Supplementary Journal Cover ofMacromolecules, Volume 55, Issue 21 [4].

computation at ISSP. Our previous work [3]was the first to discover the stress overshootphenomena of ring/linear blends under biaxial

References

elongational flow using this method. In this[1] T. Murashima, K. Hagita, T. Kawakatsu,

project,

multicyclic/linear

Nihon Reoroji Gakkaishi, 46, 207 (2018);

blend systems under biaxial elongational flows.

https://github.com/t-murash/LAMMPS-UEFEX

We found sharper overshoot phenomena of[2] S. Plimpton, J. Comp. Phys. 117, 1 (1995);

multicyclic/linear

A. P. Thompson et al., Comp. Phys. Comm.,

weinvestigated

blend

than

those

of

monocyclic/linear blend in our previous work271, 10817 (2022); https://www.lammps.org[3].

[3] T. Murashima, K. Hagita, T. Kawakatsu,

The

originofovershootinmulticyclic/linear blend was the open-to-closed

Macromolecules, 54, 7210 (2021).

topological transition. This work [4] was[4] T. Murashima, K. Hagita, T. Kawakatsu,

selected as a supplementary journal cover ofMacromolecules, 55, 9358 (2022).

Macromolecules as shown in Fig. 1.

241

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

Stability of quantum spin liquids in frustrated spin systems
Kota IDO
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
model with other perturbative terms under the

The Kitaev spin liquid has attracted attention

magnetic field along arbitrary direction.

in order to realize topological quantum
computing by Majorana fermions [1, 2].

In this project, we analyze the magnetization

However, because Kitaev candidate materials

process of the Kitaev-Heisenberg model under

have additional magnetic interactions such as

the magnetic field along [001] or [111] direction

the Heisenberg term, the ground state becomes aby using the generalized spinon pair-productmagnetic ordered state. Applying a magnetic

wavefunction, which has been already used forfield is a promising way to suppress theanalysis of quantum spin liquids in strongly

magnetic order. In fact, it is reported that thecorrelated systems [4]. We find that the VMC

half-integer quantization of the thermal Hall

results by using this trial wavefunction can wellconductivity has been measured in a Kitaev

reproduce the exact diagonalization results evencandidate material under a magnetic field, which

when we impose a sublattice structure on thesignals the realization of the Kitaev spin liquid

trial wavefunction. Applications to large size

by applying the magnetic field [3]. Therefore, itsystems will be reported elsewhere.

is desirable to accurately solve the extended

References

Kitaev model under a magnetic field.
The variational Monte Carlo (VMC) method[1] A. Kitaev, Ann. Phys. 321, 2 (2006).

is a numerical method applicable to analysis of[2] H. Takagi et al., Nat. Rev. Phys. 1, 264

two-dimensional systems [4, 5]. Recently, by(2019).

introducing the Jastrow-BCS wavefunction[3] Y. Kasahara et al., Nature 559, 227 (2018).

constructed from the Jordan-Wigner fermions,[4] M. Kurita, Y. Yamaji, S. Morita, and M.

we succeeded in analyzing the magnetization

Imada, Phys. Rev. B 92, 035122 (2015).

process of the Kitaev honeycomb model with

[5] T. Misawa et al., Comput. Phys. Commun.

high accuracy [6]. However, the application of235, 447 (2019).

this approach is limited to only the Kitaev model[6] K. Ido and T. Misawa, Phys. Rev. B 101,with the magnetic field along z-direction. Thus,

045121 (2020).

it is still challenging problem to solve the Kitaev

242

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

Efficient Sampling Simulation of
the Soft Modes Significantly Contribute to Protein Properties
D. P. TRAN, K. TAKEMURA, T.N. WIJAYA, R. OSAWA, M. YOSHIOKA, H.
NAKAYA, S. IKIZAWA, T.Y. CHU, H. T. NGUYEN and A. KITAOSchool of Life Science and Technology, Tokyo Institute of Technology
Ookayama, Meguro, Tokyo 152-8550
In this project, we employed Parallel(Dist1) and Ala146 and Val286 (Dist2) andCascade Selection Moelcular Dynamics

plot in Figure 1. We found that the FEL ofsimulation (PaCS-MD) [1] as the enhanced

CALB/W

sampling method to conduct efficientspace than that of the CALB/I, because ofmolecular simulation. Recently, we used thethe flexibility of the lid region in water. The

conmbination of the PaCS-MD together

minima i1 and i2 in CALB/I FEL are stiffer

with Markov State Modeling (MSM) to

comparing to that of w1-3 in CALB/W FEL.

unveil the energetic and kinetics of the

The free energy level of the closed

interfacial activation of Candida antarctica

conformation is much lower than that of the

Lipase B (CALB) [2].

open conformation in water.

spanslargerconformational

We first constructed the well-relaxed
structures of the CALB/W in water and theinterface of water and tricaprylin (CALB/I).
After that, we performed the PaCSMD/MSM from open-to-close and close-toopen for both CALB/W and CALB/I with30 replicas. Each trajectory length is 0.1ns.

Figure 1. Free energy landscape (FEL) of

We monitored the the α5 helix (residues

CALB/W and /I obtained by PaCS-

142‒146) and half of the α10 helix

MD/MSM(residues278‒287)

inall

PaCS-MD
After that, we performed the PCCA++

simulations.

[3] to explore macroscopic view of theThen we constructed the free energy
landscape (FEL) of both systems by using

open-closetransition

Markov State Modeling from the inter-Cα

employed

thedistances between Arg309 and Leu144

microstates into 6 microstates for both

243

of

coarseCALB.

We

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

cases. The CALB/W conformation can bei) small and hydrophilic substrates bindsplit into three states: closed states (WC1,without interfacial activation, ii) large andWC2, and WC3), one semi-open state (WS),bulky substrates bind via interfacial activation.

and two open states (WO1 and WO2). In

These findings could help expand CALBCALB/W, we found that the transition fromapplications towards a wide variety ofsemi-open state to the closed state is one-

substrates.

order of magnitude faster than that to theopen state. In contrast, the CALB/I exhibits

References

the timescale of semi-open to open states in[1] R. Harada and A. Kitao. J Chem Phys 139,

nanosecond order. CALB prefers the035103 (2013).

transition to the closed state in water but[2] T. N. Wijaya, A. Kitao (2023), ChemRxiv.

moves to the open state in the interface.

doi: 10.26434/chemrxiv-2023-17v42[3] S. Roeblitz, M. Weber, Adv Data Anal
Classif 7, 147 (2013)

Figure 2. Microstates assignment using
PCCA++ for both CALB/W and /I. PCCA++
results and MFPT between macrostates for
CALB/W(left)andCALB/I(right).

Transition time and free energy difference arealso indicated.

Via analysis of PaCS/MSM, we found thatthe closed state is the most stable in the water
system but the stable conformation in theinterface system shifts to the semi-open state.
In the interface system, the transition
probability to the open state is higher than in
the water system. We also suggest double
mechanisms of substrate binding as followed:

244

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

Thermal effects on quantum frustrated magnetisms
Tokuro Shimokawa
Theory of Quantum Matter Unit, Okinawa Institute of Science and Technology Graduate
University, Onna 904-0495, Japan
Investigating the effects of thermal fluctuations

We also report on the investigation of the staticin quantum frustrated magnets is crucial forand dynamic properties of the S=1/2 JA-JB-JC

developing our understanding of quantum spin

honeycomb-lattice Heisenberg model, which is

liquids. Ca10Cr7O28 (CCO) [1] is a quantum

believed to be important for understanding thespin liquid candidate that can be described byquantum spin liquid behavior in different

the S=1/2 bilayer breathing kagome (BBK)honeycomb-lattice material Cu2(pymca)3(ClO4)

Heisenberg model. Here, we examine the[3]. The investigations were carried out using

importance of thermal fluctuations by focusing

quantum Monte Carlo simulations based onon

ALPS application [4] and exact diagonalization.

the

J1-J2

Heisenberg

classical

modelhoneycomb-latticecorrespondingto

Our results demonstrate the realization of athequantum paramagnetic state in this material,classical system of this S=1/2 BBK model.

and suggest that focusing on dynamic properTo do this, we employ classical Monte Carloties is crucial for experimentally identifying

simulations, utilizing a hybrid parallelization

this state [5].

approach combining the heat-bath method,
over-relaxation method, and replica exchange

References

method. We investigate the ordering behavior[1] C. Balz, et al. Nat. Phys. 12 942 (2016).

of the honeycomb lattice model in a finite-[2] T. S., et al, in preparation

temperature magnetic field with a focus on the[3] Z. Honda, et al, J. Phys. Soc. Jpn. 84,

several034601 (2015).

J2/J1

parameter

values.

Physicalquantities, including specific heat, magnetic[4] A. Albuquerque, et al, J. Magn. Magn.

susceptibility, equal-time structure factor, andMater. 310, 1187 (2007).

chirality, are examined to identify the presence[5] T. S., Ken'ichi Takano, Zentaro Honda,

of a Skyrmion lattice state with sublattice

Akira Okutani, and Masayuki Hagiwara, Phys.

structure and other additional multiple-q states

Rev. B 106, 134410 (2022).

[2].

245

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

Matrix Product Renormalization Group
Masahiko G. YAMADAAB
Masahiro O. TAKAHASHIB
Takumi SANNOB
Yutaka AKAGIC
Hidemaro SUWAC
Satoshi FUJIMOTOBD
Masafumi UDAGAWAAA
Department of Physics, Gakushuin University, Mejiro, Tokyo, 171-8588, JapanB
Department of Materials Engineering Science, Osaka University, Toyonaka 560-8531, JapanC
Department of Physics, The University of Tokyo, Hongo, Tokyo, 113-0033, Japan
D
Center for Quantum Information and Quantum Biology, Osaka University,
Toyonaka 560-8531, Japan

In order to support and endorse the results
found in [1], we have reproduced the matrix
product renormalization group (MPRG) calculation done in [1] using the ISSP supercomputer. We newly find that by using the ISSP
supercomputer the calculation with a larger
bond dimension χ is possible. Indeed, the calculations with χ = 128 are possible for bothtransverse-field Ising model and the spin-1/2
antiferromagnetic Heisenberg model.

References[1] M. G. Yamada, T. Sanno, M. O. Takahashi, Y. Akagi, H. Suwa, S. Fujimoto,
M. Udagawa: arXiv:2212.13267 (2022).

The results agree well with an exact solution
or a quantum Monte Carlo simulation. Thismeans that MPRG has a potential to realizea large scale simulation up to χ = 128 anda good scalability. The advantage over othermethods will be discussed in the future, but it
is likely that MPRG is superior to other conventional tensor network methods.
Currently, the algorithm of MPRG is notsuitable for parallel computing, while it may
be straightforward to support a graphics
processing unit (GPU). We will persue thisdirection in the future study.

246

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

Magnetism and superconductivity in ternary
chromium chalcogenides
Harald JESCHKE
Research Institute for Interdisciplinary Science, Okayama University
3-1-1 Tsushima-naka, Kita-ku, Okayama 700-8530

We study CrGeTe3 which is a van-der-Waals
bound layered ferromagnet. Materials like
CrGeTe3 or CrI3 have been intensely investigated in recent years because detached monolayers still order ferromagnetically and are discussed for various nanotechnological applications.
Recently, pressure experiments for CrGeTe3
reported an interesting temperature-pressure
phase diagram; importantly, at high pressures,
an insulator to metal transition roughly coincides with increases of Curie-Weiss and magnetic ordering temperatures. Our aim is to understand the phase diagram and to investigate
the nature of the strongly correlated metal athigh pressures.
We perform calculations using density functional theory (DFT) in combination with dynamical mean field theory (DFMT), as implemented in DCore. Based on experimental crystal structures, we prepare tight binding representations of the DFT bands for pressures
up to P = 10 GPa. Figure 1 shows the spectral function of a ferromagnetic insulating solution at ambient pressure [1]. The difficultyin this material is the fact that Te 6p bands
are very important at the Fermi level. The
insulator to metal transition involves pressure
induced widening of Te 6p bands which thenestablish a Fermi surface at around 5 GPa. By
analyzing the renormalization factor Z whereZ −1 = 1 − ∂ImΣ(iω)/∂ω|ω→0+ , we find that
correlations of some Cr 3d orbitals are enhanced by increasing pressure.

Figure 1: Spectral function for the ferromagnetic state of CrGeTe3 at P = 0 GPa andT = 100 K. Chromium magnetic moments areclose to 3 µB .

References[1] H.-X. Xu, M. Shimizu, J. Otsuki, H. O.
Jeschke, in preparation.

247

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

Fluctuation exchange approximation calculations forthe superconducting transition temperatures andpairing symmetries of organic charge transfer salts
Harald JESCHKE
Research Institute for Interdisciplinary Science, Okayama University
3-1-1 Tsushima-naka, Kita-ku, Okayama 700-8530

BEDT-TTF (bisethylene-dithio-tetrathiafulvalene) is a famous nearly flat moleculethat forms hundreds of charge transfer (CT)salts, among them a number of famous superconductors. We study eight of the known
ambient pressure superconductors in the κphase where dimers of BEDT-TTF molecules
form a triangular lattice which is doped with
one hole per dimer. Our aim is to improveon earlier studies based on random phaseapproximation and two-particle selfconsistent approximation to predict the material
dependence of superconducting transition
temperature and order parameter. Our newly
developed multi-orbital, multi-site fluctuationexchange approximation (flex) code wasparallelized with OpenMP and MPI methods.

Figure 1 shows the leading solutions of theEliashberg equation for one of the materials,
together with the trend of the leading eigenvalues [1]. Tc trends are captured well by thedominant eigenvalue of the Eliashberg equation for five high-Tc materials. However, forthree low-Tc materials, the present approachoverestimates the strength of the s-wave solution. A possible reason for this discrepancy is
the intra-dimer, inter-site interaction V thathas been neglected so far and that is expected
to play an important role in the κ-type CT
salts.

Figure 1: (a), (b) Imaginary part of aninter-dimer element of the anomalous selfenergy for the two leading solutions forκ- (ET)2 Cu[N(CN)2 ]Br. (c) Eigenvalues λ foranisotropic s-wave and for s± + dx2 −y2 -wave.

References[1] M. Shimizu, D. Guterding, J. Otsuki, H.
O. Jeschke, in preparation.

248

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

Numerical studies for science of bulk-edge correspondence
andtopological phases
Yasuhiro HATSUGAI
Department of Physics, University of Tsukuba
Tennodai 1-1-1, Tsukuba, Ibaraki 305-8571
The topological phases are independent ofover the torus spanned by the time and gauge

the symmetry breaking. It implies that there is

twist. In Fig. 1, we have shown numericalno fundamental order parameter to characterize

evaluation of the Chern numbers for theit. Most of them are gapped, and low-energy

topological pump of SU(Q) interacting fermions

gapless excitation is absent. It implies that the[2]. The relation between the edge statesbulk is hidden. However, with edges orspectrumboundaries, there exist low-energy gaplessconsistently discussed by using numericalexcitations that are localized only near thestudies.

and

the

Chern

numbersaresystem boundaries and characterize the phase.
This is the bulk-edge correspondence.
A topological pump (adiabatic pump) is
another kind of topological phenomenon where
the bulk-edge correspondence again governs thephenomena. Here the bulk-edge correspondenceis special, that is, the edge states are hidden sinceFig. 1 Numerical results of the Chern number ofthe gapless nature of the edge states implies that

SU(Q) topological pump (from ref. [2]).

the edge state cannot be observed by an

References

experimental finite speed pump. Experimental[1] Y. Hatsugai, Phys. Rev. Lett. 71, 3697

observables for the topological pump are the(1993).).

center of mass (polarization) of the bulk, which[2] Y. Hatsugai, and Y. Kuno, arXiv:2210.11646

is described by the Chern number integrated

249

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

Eﬀects of ﬂux structures on Majorna excitations in
Kitaev spin liquids
Akihisa Koga
Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551

Quantum spin liquids in the Kitaev model
have been studied intensively in recent years.
It is known that in the Kitaev model, spin
degrees of freedom are split into itinerant
Majorana fermions and localized ﬂuxes dueto the spin fractionalization. The groundstate belongs to the ﬂux-free state and low
energy properties are described by itinerant
Majorana fermions. In fact, the Majoranaedge current has been observed in the halfquantized plateau in the thermal quantum Halleﬀect. Furthermore, the Majorana-mediatedspin transport has been theoretically suggested, which should stimulate further investigation to realize quantum devices using Majorana fermions.
(a)

𝒂!

𝒂"

𝒃!

𝒃"

be described by the toric code. Then a question arises. Are the above gapful states with
distinct origins are adiabatically connected toeach other? To answer this question, we consider several ﬂux conﬁgurations to investigate
how the anisotropy in the exchange couplings
aﬀects the Majorana excitations.
Motivated by this, we deal with theanisotropic Kitaev model with ﬂux conﬁgurations. In particular, we focus on the triangular
ﬂux conﬁgurations, where the ﬂuxes are periodically arranged in the honeycomb sheet, asshown in Fig. 1. Systematic numerical calculations demonstrate how the anisotropy of theexchange couplings and ﬂux conﬁguration induce the excitation gap. The induced gapped
quantum spin liquid states are distinct fromthe gapped one realized in the large anisotropic
limit. We then clarify the nature of gappedstates in the Kitaev model, regarding the ﬂux
conﬁguration as the superlattice potential forthe ﬂux-free state [2, 3].

References
(b)
1

(c)
Figure 1: Flux conﬁguration(q = 2) of theK'
Kitaev model.
X

ky/π

0.5

0

-0.5

-1

0.05 ΔY

In our previous
study [1], we have foundKΓthat the gapped quantum spin liquid is realized in a certain ﬂux conﬁguration. This
indicates that the ﬂux conﬁguration is a keyrole incontrolling Majorana excitation. On
kx/π
the other hand, it is known that the large
anisotropy in the exchange couplings stabilizes
the gapped quantum spin liquid, which should

-1

-0.5

0

0.5

[1] A. Koga, Y. Murakami, and J. Nasu,
Phys. Rev. B 103, 214421 (2021).

0.1

1[2] A. Hashimoto, Y. Murakami, and A.
0
Koga, J. Phys.: Conf. Ser. 2164, 012028(2022).
[3] A. Hashimoto, Y. Murakami, and A.
Koga, arXiv:2303.04295.

250

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

Molecular dynamics simulation of disease-related
biomolecules
Hisashi Okumura
Exploratory Research Center on Life and Living Systems,
Institute for Molecular Science, Okazaki, Aichi 444-8585

Proteins maintain their functions with

the Institute for Solid State Physics, thecorrectly folded structures. However, when

University

of

Tokyo

[1].

The

replica-their concentration increases due to, forpermutation method is an improved version of

example, aging, they aggregate to form

the replica-exchange method, in which theoligomers, spherical aggregates, and amyloidtransition ratio of replicas in the temperature or

fibrils, needle-like aggregates. These protein

parameter space is improved.

aggregates are associated with about 40 humanAs a result, we identified the key residue,

neurodegenerative diseases. For instance, theArg5, for the Aβ42 dimerization. The two C-

toterminal residues in Aβ42 tend to form contact

Alzheimer's disease. We have performed

with Arg5 because of the electrostatic attraction

several molecular dynamics (MD) simulations

between them, stabilizing the β-hairpin. This β-

of Aβ aggregates. In this fiscal year, we

hairpin promotes dimerization through theperformed MD simulations of the dimerizationintermolecular β-sheets.

amyloid-β(Aβ)

peptideisrelatedprocess of Aβ peptides and specified the key

Thus, we asked experimental collaborators

residue for aggregation of Aβ peptides.

toexaminetheeffectsofamino-acid

Aβ has two isoforms: Aβ40 and Aβ42.

substitutions of Arg5 to verify our theoretical

Although the difference between Aβ40 andprediction. We found that the mutations

Aβ42 is only two residues in the C-terminus,

remarkably suppressed the aggregation of Aβ42.

Aβ42 aggregates much faster than Aβ40. It is

In this way, our prediction was approved by the

not clear why the C-terminal two residues

experiments.

accelerate the aggregation. To investigate theperformed

References

molecular[1] S. G. Itoh, M. Yagi-Utsumi, K. Kato, and H.

dynamics simulations of two Aβ40 peptides andOkumura, ACS Chem. Neurosci. 13, 3139

those of two Aβ42 peptides using the(2022).

dimerization
Hamiltonian

process,

wereplica-permutation

supercomputer at the Supercomputer Center,251

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

Nontrivial properties of interface and ﬁnite size
eﬀects of free energies
Naoko NAKAGAWA
Department of Physics, Ibaraki University
2-1-1 Bunkyo, Mito, Ibaraki 310-8512

We investigate liquid-gas coexistence associated with ﬁrst-order transitions via molecular dynamics simulations using LAMMPS. We
managed two kinds of systems, one is a simple
ﬂuid of Lennard-Jones particles and the otheris the two components ﬂuid of argon and krypton.
In the former systems, we notice the conﬁguration of the coexistence states and examine the thermodynamic quantities. Despite theﬁniteness of the system size, we ﬁnd a sharp
interface separating the two phases, whereas
both the shape and the position of the interfaces largely ﬂuctuate.
We set temperature and mean mass density
as almost particles are in liquid, while a small
number of particles be in gas. When the aspect ratio of the box conﬁning the particles isclose to unity, we observe typically a bubble
with a round-shaped interface. Changing the
aspect ratio, the shape of the interface transitioned into a ﬂat shape, where the gas occupies
a stripe region that may situate adjacent to theboundary. To characterize the transition, we
calculate the virial pressure corresponding to
the mean pressure over the system, which may
provide the global thermodynamic properties
in phase coexistence. With the transition forthe shape of the interface, we ﬁnd a significant change in the value of the virial pressure. It drops even to a negative value when
the size of the bubble is too small. This indicates that the system behaves attractively, although the liquid or gas without coexistence,

Figure 1: Phase coexistence of liquid and gas
in the aspect ratio 1 and 4.

or the coexistence with a ﬂat interface, shows
positive pressure as is usual. Such anomalyof the virial pressure suggests that the system
is not explained by thermodynamics in equilibrium, and may possess meta-stable or unstablecharacteristics. We will proceed calculations tocharacterize the local states for identifying themechanism to maintain the metastability.
We also investigated the liquid-gas coexistence for the two components ﬂuid of argon
and krypton using the mixing free energy. The
mixing free energy is numerically determinedby applying the formula derived in [1]. We
ﬁnd that the ﬁnite size eﬀects lead to the violation of the convexity in the free energy. To
examine the stability of the coexistence states
showing the violation of the convexity, we carefully examine the possibility of hysteresis. Up
to now, we conﬁrm that rather smaller systems
do not show hysteresis in our setup. While weneed further investigations with various systemsizes, the result indicates that the ﬁniteness ofthe system size may induce the violation of the252Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

convexity in free energy in a diﬀerent context
from the metastability.
We have used the L4 CPU of System B with
4 nodes in 512 parallel computation using MPI
for the calculation of 2 × 104 particles, the L16
CPU of System B with 16 nodes in 2048 parallel computation using MPI for the calculationof 8 × 104 particles, the L32 CPU of System
B with 24 nodes in 3072 parallel computation
using MPI for the calculations of 5 × 105 and1 × 106 particles.

References[1] A. Yoshida and N. Nakagawa: Phys. Rev.
Res. 4, 023119 (2022)

253

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

Nonequilibrium relaxation temperature scaling:
Applications to quantum systems
Yoshihiko NONOMURA
MANA, National Institute for Materials Science
Namiki 1-1, Tsukuba, Ibaraki 305-0044

χst (t, δ)(δ − δc )γ ∼ fsc [ctσ + ln(δ − δc )γ ]. (1)Scaling results for the L×L periodic-boundary
system with L = 160 and inverse temperature
β = L in the loop algorithm from dimer conﬁgurations are given in Figure 1. Using theestimate γ = 1.400(3) obtained from this ﬁgure and the estimate γ/ν = 1.973(4) based on101χst(t,δ)(δ-δc)γWe numerically revealed that the critical
nonequilibrium relaxation behaviors in cluster algorithms are described by the stretchedexponential simulation-time dependence ofphysical quantities in various classical spin systems [1, 2, 3] and in a quantum phase transition [4], and derived such a relaxation formula
phenomenologically [5]. The nonequilibriumto-equilibrium (NE-EQ) scaling used in thesestudies [1, 3, 4] were generalized to oﬀ-criticalbehaviors, and such a scaling was called as the“temperature scaling” [6].
Here we show that the temperature scaling also holds in quantum phase transitions.
As an example, we consider the Néel-dimer
quantum phase transition in the columnardimerized S = 1/2 antiferromagnetic Heisenberg model on a square lattice [4]. Whenthe ratio of the strength of dimerized bonds
to other bonds are taken as (1 + δ) : 1,
the staggered susceptibility in equilibrium diverges in the vicinity of the critical point δc
as χst (t = ∞, δ) ∼ (δ − δc )−γ with δc =0.90947(6) [4]. Combining this formula with
the initial-time stretched-exponential criticalrelaxation form χst (t; δc ) ∼ exp(ctσ ), we haveχst (t, δ)(δ − δc )γ ∼ exp[ctσ + ln(δ − δc )γ ], orδ = 0.927100

δ = 0.928
δ = 0.930
δ = 0.932

10-1

δ = 0.935

L = 160
γ = 1.400(3)
10-2

-4

-2

0

2

4

δ = 0.940
δ = 0.945
6

8

10

12

cχtσ+ ln(δ-δc)γFigure 1: Temperature scaling of the staggered
susceptibility in the 2D columnar-dimerized
S = 1/2 antiferromagnetic Heisenberg model.
the NE-EQ scaling up to L = 160 [4], we have
ν = 0.710(2) only with a diagonal quantity χst .

References[1] Y. Nonomura, J. Phys. Soc. Jpn. 83,
113001 (2014).
[2] Y. Nonomura and Y. Tomita, Phys. Rev.
E 92, 062121 (2015).
[3] Y. Nonomura and Y. Tomita, Phys. Rev.
E 93, 012101 (2016).
[4] Y. Nonomura and Y. Tomita, Phys. Rev.
E 101, 032105 (2020).
[5] Y. Tomita and Y. Nonomura, Phys. Rev.
E 98, 052110 (2018).
[6] Y. Nonomura and Y. Tomita, Phys. Rev.
E 102, 052118 (2020).

254

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

HPC-based fusion of quantum simulation, experiment
analysis and data-driven science
Takeo Hoshi1,21Department of Mechanical and Physical Engineering, Tottori University,
4-101 Koyama-Minami, Tottori 680-8550, Japan.

2

Slow Positron Facility, Institute of Materials Structure Science, High Energy Accelerator
Research Organization (KEK), Oho 1-1, Tsukuba, Ibaraki, 305-0801, Japan.
We developed several new features on

surface structure as a test of the Bayesian

2DMAT [1-6], an open-source data-analysis

inference by PAMC. [8] We analyzed thesoftwareexperiment diffraction data in the one-beamforadvanced

experimental
wascondition, in which only the 𝑧 -coordinate ofdeveloped by the PASUMS project at FY2020,

atom positions are measured. The Bayesian2021[1]. A major achievement in the present

posterior probability distribution 𝑃(𝑋|𝐷) was

project is the automated optimization of thecalculated for the data set of the atom positions

massivelythe𝑋 with the given experimental diffraction data

population annealing Monte Carlo (PAMC)

set 𝐷 and reproduced the correct atomic

method [7]. As results, we found an almost

position of the 𝑁" = 6 atoms of the top-mostideal parallel efficiency not only on the ISSP

and several subsurface layers. Moreover, thesupercomputer (ohtaka) but also on the Fugaku

correct number of the atoms (𝑁" = 6) was alsosupercomputer with upto 10! nodes [8] . As an

reproduced from the calculated model evidence

application study, we used 2DMAT for the data

𝑃(𝐷). In conclusion, the PAMC analysis, as aanalysis of total-reflection high-energy positron

global search, determined the correct surface

diffraction (TRHEPD)[9]. TRHEPD is an

structure without any initial guess. Other

innovative measurement technology for theapplication studies are ongoing.

measurement

techniques.

parallel

2DMAT

computationofsurface structure analysis and was realized at
Slow

Positron

Facility,

KEK

[9].

ReferencesThe

experimental setup of TRHEPD is similar to

[1]

that

https://www.pasums.issp.utokyo.ac.jp/2DMATofreflection

high-energy

electrondiffraction (RHEED) but TRHEPD reveals the[2] K. Tanaka, et al. (T. Hoshi), Acta. Phys. Pol.

surface sensitivity and was used for the

A 137, 188 (2020).

determination of surface structure.

[3] T. Hoshi, et al. (T. Hoshi), Comp. Phys.

2DMAT

was applied to the analysis of Ge(001)-c4x2

Commun. 271, 108186 (2022).

255

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[4] T. Hanada, et al., Comp. Phys. Commun.

[7] K. Hukushima, Y. Iba, AIP Conf. Proc. 690,

277 108371 (2022).

1604 (2003).

[5] Y. Motoyama, et al.

[8] N. Kinoshita, et al. (T. Hoshi), in JPS(T. Hoshi), Comp.

Phys. Commun. 280, 108465 (2022).

meeting, 22-25, Mar. 2023.

[6] K. Tanaka, et al. (T. Hoshi), JJAP Conf.

[9] https://www2.kek.jp/imss/spf/eng/

Series 9, 011301 (2023).

256

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

Chaotic dynamics in quantum many-body systems
Masaki TEZUKA
Department of Physics, Kyoto University, Kitashirakawa, Sakyo-ku, Kyoto 606-8502
The Sachdev-Ye-Kitaev (SYK) model is amodel of fermions with independently random
all-to-all four-point interactions obeying theGaussian distribution. The model, both for
Majorana and complex fermions, has ananalytic solution in the limit of large number of
fermions. The Lyapunov exponent, defined by
the out-of-time ordered correlator, realizes theuniversal upper bound at low temperatures.
Previously, we studied the effect of arandom hopping term, two-fermion interactions,
added to the model [1]. The chaos boundbehavior is lost and the eigenenergy spectrum
becomes uncorrelated. The localizationtransition can be understood as a Fock-space
localization, a type of many-body localization.
We have obtained physical quantities related tothis transition, such as the location of thetransition point as well as the moments of theeigenstate
wavefunctions
andbipartiteentanglement entropy close to the transition
analytically and confirmed numerically [2].
A “sparse” version of the SYK model, in
which the number of randomly chosen non-zero
interaction terms is only on the order of thenumber of fermions, have been proposed a few
years ago, and for Gaussian random couplings,
the sparse model reproduces essential featuresof the original SYK model [3]. Using the ISSP
Supercomputer, we analyzed the spectral
statistics of a further simplification of thesparse SYK model, in which the magnitude of
the nonzero couplings between Majorana
fermions are set to be a constant [4]. In fact,
this simplification is an improvement, as the
spectral correlation is stronger and thus closer
to the random-matrix universality for the same

number of non-zero couplings, compared to theGaussian random case.
We also studied the quantum error
correction (QEC) capabilities of the unitary
time evolution according to the SYK type
Hamiltonians by the Hayden-Preskill protocol.
This is a process in which one tries to decode
an unknown quantum state thrown into aquantum many-body system with theknowledge of the initial state of the system andthe time evolution, by accessing only a part ofthe system after the evolution. The error
estimate for random unitary evolutions obeying
the circular unitary (Haar) ensemble has been
known. For the sparse SYK model, the QEC
error estimate approaches the Haar value after ashort time if the spectral statistics is randommatrix like. In the presence of two-fermion
interactions, departure from the Haar value isobserved before the Fock space localization [5].

References[1] A. M. Garcia-Garcia, B. Loureiro, A.
Romero-Bermudez, and M. Tezuka, Phys. Rev.
Lett. 120 (2018), 241603.
[2] F. Monteiro, T. Micklitz, M. Tezuka, and A.
Altland, Phys. Rev. Research 3 (2021), 013023;
F. Monteiro, M. Tezuka, D. A. Huse, and T.
Micklitz, Phys. Rev. Lett. 127 (2021), 030601.
[3] S. Su et al., arXiv:2008.02303; A. M.
Garcia-Garcia et al., Phys. Rev. D 103, 106002(2021).
[4] M. Tezuka, O. Oktay, E. Rinaldi, M.
Hanada, and F. Nori, Phys. Rev. B 107,
L081103 (2023).
[5] Y. Nakata and M. Tezuka, arXiv:
2303.02010.

257

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

Molecular Dynamics Simulation of Dissolution of Stretched PEG
Crystal in Aqueous Solution
Koichi Mayumi, Kosuke Aomura, Haru Shinohara

The Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Hydrogel is a material in which

consisting of 24 PEG chains with 42

water is confined in a polymer network.

repeat units were created (Fig. 1). After

Hydrogel has high biocompatibility andenergy minimization calculations, 4200

is expected to be applied to biomaterials,

water molecules are added to the box so

but its mechanical strength has been athat the polymer fraction was 38 wt%. A

problem. In our laboratory, we havestretching force was applied to bothdeveloped

both

ends of each polymer chain for 400 ns.

toughness and resilience by using the

The applied force f was set at 6 points: 0,

self-reinforcing effect [1]. In this self-

10, 30, 70, 100, and 210 kJ/(mol nm).

reinforcing

GROMACSahydrogel

gel,

with

polyethylene

glycol2016.5

was

used

as

(PEG) is used as the main chain of thecomputational software. Simulations

polymer network and water as theweresolvent. When the strain is removed, thesupercomputer system B. The number

PEG chains form hydration with theof atoms in theperformedon

Ohtaka,ISSP

solvent water, and the crystals arequickly melted. Such solvation anddissolution processes are thought toplay an important role in the straininduced crystals that are reversible
against elongation. In this study, we
performed all-atom molecular dynamics
simulations to understand extended
PEG

crystaldissolvesinaqueoussolution.

Fig. 1 Snap-shots of PEG zig-zag crystals
stretched under f = 10 kJ/(mol nm) at time t= 0 and t = 400 ns.

Based on the structures reported from
X-ray diffraction results [2, 3], all-atommodels of zigzag PEG chain crystals

258

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

simulation box was about 24,000, andchain crystal is dissolved and then the athe computation speed was about 200

helical PEG crystal is formed..

ns/day with a parallel number of 512.

References[1] C. Liu, N. Morimoto, L. Jiang, S.
Kawahara, T. Noritomi, H. Yokoyama,
K. Mayumi, K. Ito, Science, 372, 10781081(2021).
[2] Takahashi, Y., Tadokoro, H.
Macromolecules, 6, 672-675 (1973).
[3] Takahashi, Y., Sumita, I., Tadokoro,
H. J. Polym. Sci., B, Polym. Phys., 11,2113-2122 (1973).

For f = 0 kJ/(mol nm), the zigzag
chain crystals are dissolved within 200
ns, which is consistent with the rapid
dissolution of stretched induced PEG
crystals just after the applied stretching
force is unloaded in our experiments.
When f is 10 kJ/(mol nm), the zigzag

259

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

Molecular dynamics simulation for chemical reactions in
multicomponent systems of thermosetting resins
Yutaka Oya
Department of Materials Science and Technology,
Tokyo University of Science, Niijuku, Katsushika-ku, Tokyo 125-8585
In this year, we have studied following two

simulations were performed for reproducing the

topics. The first is to establish the methodology

equilibrium structure of the system including IT-

for reproducing the crosslinked reactions of

4F (donor molecule) and PDTB (accepter

thermosetting polymers based on full-atomistic

molecule) as shown in Figure 1. Further,molecular dynamics simulation. We focused on

frequency distributions for each type of stacking

multicomponent system including base resin

structure in the system were analyzed. As awith epoxy groups, curing agent with amine

result, it was found that J-aggregation between

groups and lignin with hydroxyl groups, which

IT4F and PDTB molecules is dominant in theis one of promising candidates for matrix of bio-

system. In the future study,

based composites. The chemical reactions
among these components were successfully
realized for systems of several hundred thousand
atoms by LAMMPS software.

In the future

study, we systematically investigate the thermal
and mechanical characteristics, such as Young’s
modulus, strength, glass-transition temperature,
for the crosslinked structures formed by these

Figure 1. (a) Equilibrium structure of IT4F andmolecules.

PDTB molecules. (b) Snapshots of typicalThe second is to investigate the stacking

stacking structures obtained by molecular

structure of organic photovoltaic materials. It is

dynamics simulations.

well known that electron transfer occurs throughthe π-π stacking structure of the donor andacceptor molecules. Such photoexcited reactionsof conjugated polymers are expected to be
appliedtonondestructive

inspectionofmechanical damages of a composite. In thisstudy,

full

atomicmolecular

dynamics

260

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

Numerical study on low-energy states of quantum
spin systems
Hiroki NAKANO
Graduate School of Science, University of Hyogo
Kouto, Kamigori, Hyogo 678-1297, Japan

Since researchers are often faced with
many-body problems in the condensed matter physics, it is diﬃcult to estimate physicalquantities precisely. One of the typical issues isa quantum spin system. To study a quantum
spin system, numerical approaches have widelyand eﬀectively been employed. Computational
studies that were carried out have provided us
with useful information of the target systems.
Within the field of quantum spin systems,
three methods have been eﬀectively used. The
first and traditional one is the numerical diagonalization method. The second method is
the quantum Monte Carlo (QMC) simulation.
The third one is the density matrix renormalization group (DMRG) method. Each of thesemethods has advantages; however, it has also
disadvantages at the same time. The QMC
simulations can treat significantly large systems irrespective of the spatial dimensions of
the systems; on the other hand, the negative
sign problem with this approach prevents us
with a diﬃculty in precise evaluation of physical quantities in frustrated systems. When thespatial dimension of a target system is one,however, the DMRG method is very useful irrespective of whether or not the target system includes frustration. Unfortunately, thismethod is still under development for the cases
when the spatial dimension is larger than one.
We can apply the numerical diagonalization
method to various systems irrespective of thepresence of frustrations and the spatial dimension. This method, on the other hand, hasa severe weak point that it can treat only
systems with very small sizes. To overcome
this disadvantage of the numerical diagonalization method, we developed a hybrid-typeparallelized code of Lanczos diagonalization[1].
This code enables us to treat various large systems that have not been treated yet within thismethod. We, therefore, study various quantum spin systems by this method as a primary
approach in this project.
In the project in 2022, we tackled theS = 1 Heisenberg antiferromagnet in one
dimension[2]. We successfully carried out
Lanczos diagonalizations of this system up to
30 sites under the twisted boundary condition. The 30-site system was treated within
Lanczos-diagonalization studies for the firsttime to the best of our knowledge. Note herethat the largest treated matrix dimension is18,252,025,766,941. Our results were analyzed
via a convergence-acceleration technique; wefinaly concluded∆/J = 0.410479239 ± 0.000000041
as the Haldane gap, namely, the spin excitation
above the nondegenerate ground state. We
confirmed that this result agrees with the result from the DMRG method[3]. We also studied the S = 1/2 Heisenberg antiferromagnet on
the floret-pentagonal lattice[4]. Our calculations successfully clarified the behavior aroundthe five-ninth of the saturation magnetizationin its magnetization process. Our studies contribute much to our deeper understandings of261

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

the quantum magnetism. Further investigations would clarify nontrivial quantum eﬀects
in these systems.

References[1] H.
NakanoandA.
Terai:
J. Phys. Soc. Jpn. 78, 014003 (2009).
[2] H. Nakano, H. Tadano, N. Todoroki and
T. Sakai: J. Phys. Soc. Jpn. 91, 074701(2022).
[3] H. Ueda and K. Kusakabe: Phys. Rev. B
84, 054446 (2011).
[4] . R. Furuchi, H. Nakano, and T. Sakai:arXiv:2209.13887 to be published in JPS
Conf. Proc.

262

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

Study of ground states and dynamics in frustrated
quantum magnets using MPS-based simulations
Matthias GOHLKE
Theory of Quantum Matter Unit, Okinawa Institute of Science and Technology
1919-1 Tancha, Onna-son, Okinawa 904-0495

Competing interactions—or frustration—
between the constituents of a many-body system have proven to be an origin of exotic
phases of matter. A famous example in thiscontext are quantum spin liquids (QSL), that
do not exhibit conventional magnetic orderof spin-dipolar moments, but rather are defined by emergent gauge fields, fractionalized
excitations, and topological order. A topicalexample in this context are Kitaev-like frustrated magnets, that exhibit a bond-dependent
spin-exchange resulting in exchange frustration and, for certain choices of parameters,
can result in a QSL ground state. Part ofthe computational resources provided by ISSP
has been used to study the related anisotropic
Kitaev-Γ model [see (II)].

cally been observed in frustrated ferromagnets
with ring-exchange on the square lattice [2].
Experimentally, the nature of such a groundstate is intrinsically difficult to verify, due tothe lack of probes that couple directly to thespin-quadrupole moments. Instead, it is necessary to examine the dynamics of a QSN: A
continuous symmetry for the director of a spinquadrupole remains, that give rise to a gapless
Goldstone mode [3, 4].

On the other hand, ordering of higher moments like spin-quadrupoles may occur as well.
Unlike spin-dipolar order, spin-quadrupolarorder is defined by a director—an arrow without a head—in close analogy to the rods ofnematic order and liquid crystals. Thus, spinquadrupolar order is also named quantum spin
nematic (QSN). Part of the resources havebeen used to study the QSN phase stabilizing in a frustrated ferromagnet with external
magnetic field [see (I)].

Here, we study the square-lattice frustrated
J1 -K model [1], with dominant ferromagnetic Heisenberg exchange J1 between nearest
neighbor spins, cyclic ring-exchange K aroundsquares, and hz the Zeeman coupling to a magnetic field along the z-axis. Using iDMRG andthe matrix product states (MPS) framework as
well as a recently developed exact diagonalization (ED) method near saturation [5], we confirm the existence of the QSN phase in an extended range of K/|J1 |. Overall we find a good
agreement between high-field ED on symmetric clusters and iDMRG on cylindrical geometries with a circumference of Lcirc = 6 sites.
The QSN phase is sandwiched between highfield polarized and 4-sublattice AF phase. Additional phases (Neél, vector-chiral) occur forlarge K.

(I) QSN in a Spin- 12 Frustrated Ferromagnet. Here, we focus on a quantum magnet
with spin- 12 degrees of freedom. Such magnets
can only exhibit a QSN states if two spin- 12
are combined into an effective spin-1 [1]. In
fact, such QSN states on bonds has theoreti-

Given the ground state as an MPS, we
can then obtain dynamical properties via atime-evolution unitary U (dt) represented asa matrix product operator. We confirm thecondensation of two-magnon bound state at
the corresponding wave vector of the QSN

263

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

Figure 1: TM spectrum illustrating dominantcorrelations at kx = 0 as well as incommensurate correlations near kx ≈ 0.35π.
state. Within the QSN phase, we observe (1)the magnon-like excitation that remains fully
gapped, and (2) a gapless mode with its spectral weight in the dipolar structure factor vanishes as q → 0 & ω → 0 as predicted by priormean-field studies. We do, however, observe
relevant qualitative differences: an additionalmode with a gapless node at an incommensurate, field-strength dependent wave vector
occurs. Within this term, we have focusedon identifying the incommensurate mode bymeans of the MPS transfer matrix (TM). The
TM spectrum contains all the information on
correlations. Momentum information is contained through the complex phase of the TM
eigenvalues, and momentum along the circumference being a well defined quantum number
as long as translation symmetry remains intact. Thus, the position in momentum spaceof the gapless modes can be extracted from the
MPS, cf. Fig. 1.. Indeed, we find perfect agreement between the gapless mode in the dynamics and the TM, which will help us in further
elucidating the nature of the incommensurate
gapless mode accompanying the QSN phase.
All MPS-related simulations have been carried out using TenPy [8]. TenPy makes useof Cython and calls heavily optimized Lapack
and Blas routines matrix operations, e.g. singular value decomposition.
(II) Extended QSL in an Anisotropic
Kitaev-Γ Model. The characterization ofquantum spin liquid phases in Kitaev materials has been a subject of intensive studies over
the recent years. Most theoretical studies havefocused on isotropic coupling strength in an attempt to simplify the problem. We, instead,
have focused on an extended spin-1/2 KitaevΓ model on a honeycomb lattice with an additional tuning parameter that controls the coupling strength on one of the bonds [6], allowing us to connect the limit of isolated Kitaev-Γ
chains, which is known to exhibit an emergent
SU (2)1 Tomonaga-Luttinger liquid (TLL) [7],to the two-dimensional model. We find thatthe critical properties of the TLL persist for
finite inter-chain coupling, and a QSL develops in analogy to sliding Luttinger liquids thatdiffers from the Kitaev spin liquid.
We have employed the TenPy package [8] toobtain the ground state in MPS form, study its
scaling behavior via the MPS transfer matrix
spectrum, and obtain dynamical properties.

References[1] A. Andreev and I. Grishchuk, Sov. Phys.
JETP 60 (1984) 267.
[2] N. Shannon, T. Momoi, and P. Sindzingre, Phys. Rev. Lett. 96 (2006) 027213.
[3] R. Shindou, S. Yunoki, and T. Momoi,
Phys. Rev. B 87, (2013) 054429.
[4] A. Smerald, H. T. Ueda, and N. Shannon,
Phys. Rev. B 91 (2015) 174402.
[5] H. Ueda, S. Yunoki, and T. Shimokawa,Comp. Phys. Comm. 277 (2022) 108369.
[6] M. Gohlke, J. C. Pelayo, T. Suzuki,
arXiv:2212.11000 (2022). See also report
by T. Suzuki.
[7] Yang et al. Phys. Rev. Lett. 124, 147205(2020).
[8] J. Hauschild, and F. Pollmann, SciPost
Phys. Lect. Notes 5 (2018).

264

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

Dynamical properties of the extended Kitaev-G model on ahoneycomb lattice
Takafumi SUZUKI
Graduate School of Engineering,
University of Hyogo, Shosha, Himeji, Hyogo 670-2280
dispersion curves for the low-lying excitation
show a linear behavior around gapless wavenumber points, which is similar to thecharacteristics observed in the S=1/2
antiferromagnetic Heisenberg chain. In thisproject, we have investigated the remaining
part; the ground state and the excitations for K≦0 and G≧0. We have found that in the spin
chain limit, an incommensurate phases exists
between the TLL phase and D4-symmetrybreaking phase [5]. In addition, we have foundthat the presence of the extended Kitaev spin
liquid [5] is still controversial [7]. This implies
that the extended Kitaev spin liquid phase isreplaced by a different state from the D3ordered state [5]. These states we newly foundin the chain limit changes two-dimensional
long-range orders when the interchaininteraction becomes finite [7].

We have used ISSP supercomputer resources
for the following two topics: (1) dynamical
properties of the extended honeycomb-lattice
Kitaev-G model and (2) supersolid states in
two-dimensional hard-core bosonic Hubbard
model with dipole interactions.
(1) Dynamical properties of the extended
honeycomb-lattice Kitaev-G model
In this project, we have investigated groundstates and dynamical properties of the S=1/2
Kitaev-G model [1] on a honeycomb lattice by
the numerical exact diagonalization anddensity-matrix-renormalization-group
approaches. The S=1/2 Kitaev-G model on ahoneycomb lattice was proposed as theeffective model of the honeycomb-latticemagnet, a-RuCl3. So far, the ground-statephase diagram of this model with the isotropic
coupling strength has been studied intensively.
In recent experiments on the related compounds,
RuX3 (X=Br, I) [3], it has been reported that
applying pressure can control the anisotropy ofthe coupling strength. The ground-state phase
diagram of the S=1/2 Kitaev-G model withanisotropic interactions has not been fully
understood yet. Note that the anisotropy meansthat by tuning the anisotropy, the modelchanges from the Kitaev-G spin chains to theisolated dimer model via the isotropicallyinteracting model.
In the previous study, we focused on the
ground phase diagram and the low-lying
excitation in the case where the Kitaev
interaction (K) and the G interaction (G) arenegative and positive, respectively [4]. We
clarified that the Tomonaga-Luttinger liquid(TLL) appearing in the spin chain limit (d=0)[5] persist for the finite interchain interactions
and it becomes a proximate TomonagaLuttinger liquid (pTTL), which is analogous tosliding TLL [6]. In the pTTL phase, the(2) Supersolid states in two-dimensional hardcore bosonic Hubbard model with dipole
interactions
In this project, we have studied the phase
diagram of two-dimensional hard-core bosonic
Hubbard model with dipole interactions. Manyefforts have been made to realize supersolidstates so far and ultracold-atomic gases trappedin optical lattices are candidate stages to study
the supersolid states. In lattice Bose-Hubbard
models, frustrated interactions play a key roleto stabilize the supersolid states. In ref. [8],authors investigated the phase diagram of thetwo-dimensional bosonic hardcore Hubbard
model with the dipole interactions by the meanfiled approach and the infinite entangled-pairstate. It was found that several supersolid states
appear in between solid phases with different
commensurate fillings by tilting the angle ofthe dipole axis from the direction perpendicular
to the lattice plane. However, the range of thedipole interactions is restricted within a short265

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

range in ref. [8]. Thus, when the farther
neighbor interactions are included, whether
those supersolid phases are stable or not isunclear.
We have investigated the phase diagram of thetwo-dimensional bosonic hardcore Hubbard
model with the dipole interactions by quantum
Monte Carlo computations and have applied
two methods to discuss the details of eachphase; finite-size-scaling analysis and machinelearning assisted approach. First, we havestudied the phase diagram of the same model
considered in ref. [8]. We have obtained
consistent results with those by the iPEPS
calculations. Next, we have investigated thephase diagram by changing the cut-off range ofthe dipole interactions. It has been found thattwo quarter solid states (r=1/4 and 3/4, where r
is filling.) becomes unstable and diagonal-stripe
solid states (r=1/3 and 2/3) becomes stable.
The supersolid state related with the quarter
solid disappears.

References[1] J. G. Rau, E. Lee, and H. -Y. Kee, Phys. Rev.
Lett. 112, 077204 (2014).
[2] T. Yamada, T. Suzuki, and S. Suga, Phys.
Rev. B 102, 024415 (2020).
[3] Y. Imai, et al., Phys. Rev. B 105, L041112(2022) and Y. Choi, et al., Phys. Rev. B 106,174430 (2022).
[4] M. Gohlke, J. C. Pelayo, and T. Suzuki,
arXiv:2212.11000.
[5] W. Yang, et al., Phys. Rev. Lett. 124,
147205 (2020) and W. Yang, A. Nocera, and I.
Affleck, Phys. Rev. Research 2, 033268 (2020).
[6] R. Mukhopadhyay, C. L. Kane, and T. C.
Lubensky, Phys. Rev. B 64. 045120 (2001).
[7] M. Gohlke, J. C. Pelayo, and T. Suzuki, in
preparation.
[8] H.-K. Wu and W.-L. Tu, Phys. Rev. A 102,053306 (2020).

266

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

Systematic Searches for Grain Boundary Structures
and Structure-based Prediction of Multiple
Properties
Susumu FUJIIDivision of Materials and Manufacturing Science, Osaka University
Yamadaoka, Suita, Osaka, 565-0871
Nanostructures Research Laboratory, Japan Fine Ceramics Center
Mutsuno, Atsuta, Nagoya 456-8587

Grain boundaries (GBs) are ubiquitus in inorganic solid materials except for glasses. At
GBs, electronic and atomic structures thatare diﬀerent from its corresponding crystal
structure are formed, providing unique properties that often determine macroscopic physical responces of materials. The GB structures
vary signiﬁcantly depending on the orientation between adjoined crystal grains, resulting in a variety of GB properties. Revealing
GB structure-property relationship is important to improve the performance of structuraland functional materials based on material microstructures.
To correlate GB structures with material
properties, there are two technological diﬃculties. One is that a large-scale data for GB
structures and their properties, with an accuracy comparable to ﬁrst principles calculations, is necessary to provide a general andsystematic understanding that is applicable toexperimental design of microstructures. Another is that the computational and machinelearning method to quantify GB structures andcorrelate them with properties is not well established. Regarding the latter, the authors
reported that the use of structure descriptor
and hierarchical clustering can clasify the local atomic environments composing GB structures, and multiple linear regression with theclassiﬁed data provide an accurate predictionof GB thermal conductivities for MgO withphysical interpretations [1].
The main drawback of our previous work isthat we use an emprical potential to derive
MgO GB structures and their thermal conductivities. In this study, we applied an accurate machine learing potential for silicon, distributed in MACHINE LEARNING POTENTIAL REPOSITORY [2], to our GB structure searches and property evaluations, to give
more precise GB structure-property relationship [3]. Silicon is known as a typical semiconductor, and considered as a promising candidate for thermoelectric material. To reduce its
high thermal conductivity, nanocrystalline silicon with high GB population has been extensively synthesized, leading to very low thermal
conductivity less than 5 W/mK [4]. Revealing
the impact of Si GBs on thermal conduction
as well as electronic conduction is necessary
for further improvement of its thermoelectric
conversion eﬃciency.
We mainly performed GB structure searcheswith the methods of rigid body translations (RBTs) and structure optimization, andthermal conductivity calculations using perturbed molecular dynamics (MD). The necessary codes for RBTs and perturbed molecular
dynamics were developed by the author, andthe latter was implemented to the Large-scale
Atomic/Molecular Massively Parallel Simula-267

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyotor (LAMMPS) [5].
Using the code for RBTs and LAMMPS forstructure optimization, we derived 45 silicon
GB structures. Then, perturbed MD simulations were performed for calculating thermal
conductivity of each GB structure. The perturbed MD simulations are highly computationally demanding, requiring large computational cells containing from three to thirteen
thousands of atoms and with the time stepsof more than one million. Repeating this simulation ﬁve times with diﬀerent magnitude of
perturbations for each GB yielded a set of statistically accurate thermal conductivities. As
230 separate perturbed MD simulations werenecessary, all of them were performed in a single node rather than multiple nodes.
Compiling the GB structure data and thermal conductivity shows that bond angle strains
at GBs signiﬁcantly decreases their thermal
conductivities. In addition, structural voids,
formed at a few GBs, largely decreases GB
thermal conductivities. These trends resultsfrom the features of covalent silicon, where one
atom bonds strongly with four atoms to forma tetrahedron, and the strong directionality ofthe bonds tends to create voids at GBs.

Figure 1: An example of silicon grain boundary structure derived using a machine learningpotential.
[3] B. Jugdersuren et al., Phys. Rev. B. 96(2017) 014206.
[4] S. Fujii and A. Seko, Comput. Mater. Sci.
204 (2022) 111137.
[5] S. Plimpton, J. Comput. Phys. 117 (1995)
1.
[6] J. Hickman and Y. Mishin, Phys. Rev.
Mater. 4 (2020) 033405.

In contrast to our results, a revious reportusing an empirical force ﬁeld claims that GB
energy is the dominant factor for determining GB thermal conductivity of silicon, rather
than structural features of GBs [6]. A large
scale calculations using supercomputers andthe accurate machine learning potential indicate the importance of precisely estimating GB
structures and their properties to reveal a robust GB structure-property relationship.

References[1] S. Fujii et al., Nature Commun. 11 (2020)
1854.
[2] A. Seko, arXiv:2007.14206.

268

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

Kinetics of phase transition and polyamorphism
Kazuhiro Fuchizaki, Katsumi Nakamura, and Hiroki Naruta
Department of Physics, Ehime University, Matsuyama 790-8577

In FY2022, we continued to undertake thetwo tasks dealt with in FY2021: One was addressing the validity of the two-phase thermodynamic (2PT) model put forth by Goddard
et al. [1], and the other was a new trial application of machine learning (ML) to judge theorder of a phase transition.
One of our major issues is to explain a signof the Clapeyron slope of a liquid-liquid transition (LLT) boundary. Liquid phosphorus [2]and sulfur [3] are known to undergo an LLT.
However, the former (latter) LLT boundary
has a negative (positive) slope, although both
liquids are believed to be brought to highdensity liquids via polymerization. Because a
specific volume must be reduced upon an LLT,
the entropy diﬀerence between high- and lowdensity liquids can determine the slope’s sign.
Therefore, a method to precisely evaluate a liquid entropy is required. The 2PT model has
been widely used for this purpose. However,
because the model was based on many assumptions, we must examine the validity before applying the model.
The second task is a brand-new research
project that has not been tried thus far. We
could finish the project in FY2022, and thesubmitted paper is now under review. Most of
the computation time was exhausted in preparing the learning data.

where D and DHS are the diﬀusion coeﬃcients for the system in question and the corresponding hard-sphere fluid, respectively. ρ
and T are the system’s density and temperature. We explicitly include the dependenceof DHS on the hard-sphere diameter σ. The
neat trick of Goddard et al. in deriving aclosed equation to determine f was to replace
the RHS of Eq. (1) with D(ρ, T )/D0HS (T, ρ; σ),thereby replacing the definition of f as f =D(ρ, T )/D0HS (T, ρ; σ). Enskog’s theory thenallows us to relate DHS to the diﬀusion coeﬃcient D0HS and the compressibility factor
obtainable from the Carnahan–Starling equation of state. Furthermore, D0HS (T, f ρ) =D0HS (T, ρ)/f is valid in the zero-pressure limit
at which the Chapman–Enskog result is valid.
These relations regarding the diﬀusion coeﬃcients resulted in the cubic equation for f .

The self-diﬀusion coeﬃcient of a hard-spherefluid
Our detailed examination ascertained that the
2PT model bears its validity in estimating theliquid fraction f , which is implicitly defined by

On our way to evaluating the RHS of Eq. (1),
we noticed that the empirical formula proposedby Speedy [5] could reproduce the simulated
DHS quite satisfactorily. This confirmation enabled us to interpolate the RHS with other f s
without carrying out simulations. Finally, we
numerically solved Eq. (1), whose RHS was

f = D(T, ρ)/DHS (T, f ρ; σ),(1)Here, we directly solved Eq. (1) instead ofsolving the cubic equation. We employed themodified Lennard-Jones (mLJ) system with
2048 particles for this judgment. We first evaluated the diﬀusion coeﬃcient D, Eq. (1) numerator, for the mLJ system at the thermodynamic conditions reported in [4]. Next, we conducted molecular dynamics simulations for a
hard-sphere fluid to find DHS . We set ρ = 0.05,
0.331, and 0.6 and σ = 1.0, 1.1, 1.2, and 1.3,
where ρ and σ are presented in the mLJ units.

269

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

thus evaluated, to obtain the desired f , which
was found to be relatively close to that obtainable from the cubic equation. Therefore, regarding estimating the molecular fraction, we
concluded that the treatment of Goddard etal. was adequate.
Judging the (dis)continuity of phase transition
using ML [6]Our short-term memory often causes “memory hysteresis” against back-and-forth input
of continuously varying input. For example,
imagine a series composed of gradually varying
pictures from A to B. Suppose that an A-likepicture certainly changes to a B-like when we
see picture C in the series when we see thosepictures from A to B. When we see them in
the reverse order, a picture in which we recognize a definite change from B-like to A-like will
diﬀer from C, thus inducing memory hysteresis. The size of hysteresis may increase withthe smoothness of variation. Because an artificial neural network (ANN) is so constructedthat it can reproduce the functionality of a biological neural network (BNN), an ANN would
behave like a BNN when a set of continuously
changing patterns is input in back-and-forth
directions. We then expected that an ANN
could distinguish a continuous or discontinuous phase transition from its hysteresis behavior when it sequentially learns a series of phasetransition patterns generated as a function of
external parameters such as temperature or external field.
We implemented a convolutional neural network (CNN) to realize an ANN. To test
the phase-transition order, we employed twodimensional Ising and q-state Potts models
with q = 3 through 6. For the Ising model,
we generated spin configurations below andabove the critical temperature Tc as a function of the external field h. For the Potts
model, we obtained them around the phasetransition point as a function of temperature
T . h or T constituted the labels. We let the
CNN learn the diﬀerence between the configu-

rations with adjacent labels in increasing (decreasing) order of the labels to assess the hysteresis size. The learning eﬀects were stored in
the CNN’s weights, and we monitored the sumof the weights, Wsum , connected to the output
layer. For the specifications of the CNN andthe models, see [6].
We have already reported the results of thepreliminary study [4]. For the Ising crossovertransition, hysteresis appeared in Wsum ’s behavior plotted as a label function, whereas nohysteresis was recognizable for the transition
below Tc as expected. However, unlike theIsing case exhibiting an order-to-order transition, the Potts model did not allow the CNN tolearn the configurations in reverse order, i.e.,from the disordered to the ordered state. The
CNN did not seem to capture the disorderedstate uniquely.
We could resolve this issue by “renormalizing” the q states into a single state to
make the CNN facilitate capturing the Potts(dis)ordering aspect [6]. We confirmed that
Wsum ’s hysteresis arose for q = 3 and 4,
whereas Wsum took maximum precisely at thetransition point without hysteresis regardlessof the learning direction. Interestingly, Wsum
behaved marginally for the weak first-ordertransition realized when q = 5.

References[1] S-T. Lin, M. Blanco, and W. A. Goddard:
J. Chem. Phys. 119 (2003) 11792.
[2] Y. Katayama et al.: Nature 403 (2000)
170.
[3] L. Henry et al.: Nature 584 (2020) 382.
[4] K. Fuchizaki et al.: ISSP Supercomputer
Center Activity Report (2021) 256.
[5] R. J. Speedy et al.: Mol. Phys. 66 (1989)
577.
[6] K. Nakamura and K. Fuchizaki: submitted to J. Phys. A.

270

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

Numerical Study of Magnetism in the
Honeycomb-Lattice Spin Systems
Chitoshi YASUDA
Department of Physics and Earth Sciences, Faculty of Science,
University of the Ryukyus, Okinawa 903-0213, Japan
Many low-dimensional antiferromagnets
with excitation modes separated from theground state by a ﬁnite energy gap have beensynthesized, and eﬀects of impurities on thosematerials in the gapped state have attracted
our interest. In the previous work [1], thespin S = 1/2 and S = 1 two-dimensional
quantum Heisenberg antiferromagnets on asquare lattice, which is composed of bondalternating chains and interchain interactions, have been investigated by the quantum
Monte Carlo simulation with the continuousimaginary-time loop algorithm, and groundstate phase diagrams parametrized by thestrength of the bond alternation and that of
the interchain coupling have been determined
precisely. The phase diagrams reveal under what conditions spin-gapped states andmagnetically ordered states are realized. Recently, various honeycomb-lattice systems areinterested in relation to quantum spin liquidstates. For example, the Kitaev model on
the honeycomb lattice has the qunatum spin
liquid state as the ground state.

this technique, the S = 1 system is represented by an S = 1/2 system with special
boundary conditions in the imaginary-time
directions. The phase-transition points of theground state are estimated by the ﬁnite-sizescaling analyses of the correlation lengths
evaluated by the second-moment method.
For valid analyses, it was necessary to calculate up to a system of size Lx × Ly =96 × 48, where the temperature is ﬁxed to
T = 1/Lx . The critical exponent for the correlation length ν is ﬁxed at 0.71 consistent
with the value for the three-dimensional classical Heisenberg model. In the S = 1/2 system, there are the dimerized state and the antiferromagnetic long-range ordered (AFLRO)
state as the ground state. On the other
hand, in the S = 1 system, there is the Haldane state in addition to the dimerized state
and the AFLRO state. The phase diagram
parametrized by the strength of the bond alternation and that of the interchain couplingshow that the region of the dimerized stateof the honeycomb-lattice system is narrower
than that of the square-lattice system.

In this project, we research the S = 1/2
and S = 1 quantum Heisenberg antiferromagnets on the anisotropic dimerized honeycomb lattice using the quantum Monte Carlosimulations with the continuous-imaginarytime loop algorithm [2], and compare the results with those of the square lattice. The
honeycomb lattice is composed of bondalternating chains and interchain interactions. Periodic boundary conditions are imposed in the x and y directions. For the simulation of the S = 1 system, the subspinrepresentation technique is applied [2, 3]. In

References[1] M. Matsumoto, C. Yasuda, S. Todo, and
H. Takayama, Phys. Rev. B 65, 014407(2001).
[2] S. Todo and K. Kato, Phys. Rev. Lett.
87, 047203 (2001).
[3] K. Harada, M. Troyer, and N.
Kawashima, J. Phys. Soc. Jpn. 67,
1130 (1998).

271

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

Spin Dynamics Simulation of the Z2-vortex Fluctuations
Hikaru KAWAMURA
Molecular Photoscience Research Center, Kobe University, Kobe 657-8501
Yo P. MIZUTA
Graduate School of Engineering Science, Osaka University, Toyonaka 560-8531
Kazushi AOYAMA
Graduate School of Science, Osaka University, Toyonaka, 560-0043
In as early as 1984, it was pointed out that theresolution by an order of magnitude.4) The

frustrated isotropic Heisenberg magnets in two

dynamical spin structure factor computed at T

dimensions (2D) could possess a novel vortex

slightly above TV is shown in the figure below,

characterized by the parity-like two-valued

which indeed exhibits a narrow central peak, in

topological number corresponding only to its

addition to the side peaks of spin-wave origin.

presence/absence, a Z2 vortex, which drove a

Quantitative analysis of the data indicates that

topological transition at a finite temperature TV

the width of the central peak indeed becomes asassociated with its binding-unbinding.1) Whilenarrow as ~0.001J in the close vicinity of the K-

the direct experimental observation of the Z2
vortex has remained elusive for years, a very
recent quasi-elastic neutron scattering (QENS)
experiment performed on the powder sample ofthe quasi-two-dimensional (2D) triangularlattice Heisenberg antiferromagnet NaCrO2 has
succeeded in directly probing the signature of
the free Z2 vortex via the observation of a sharp

point slightly above TV (=0.285J), consistently

quasi-elastic (QE) scattering of its energy widthwith the recent QENS experiment. 2)

as narrow as ∼10 μeV corresponding to ~0:001J

References(J the exchange coupling).2) The observed peak[1] H. Kawamura and S. Miyashita, J. Phys. Soc.

width was even narrower than that in the earlaierJpn. 53, 4138 (1984). [2] K. Tomiyasu et al, Phys.

model simulation3) by an order of magnitude.

Rev. B 106, 054407 (2022). [3] T. Okubo and H.

We then have performed further extensive spin-

Kawamura, J. Phys. Soc. Jpn. 79, 084706 (2010).

dynamics simulation on the frustrated 2D[4] Y.P. Mizuta et al, J. Phys. Soc. Jpn. 91,

Heisenberg model by improving the energy

035001 (2022).

272

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

Finite temperature calculation of quantum manybody
system using
random phase product state and neural network
wavefunction
Toshiaki IITAKA
Center for Computational Science,
RIKEN, 2-1 Hirosawa, Wako, Saitama, 351-0198 JAPANCalculation of thermal average of an

Thermal average of observable A atobservable A of a closed quantum many-

inverse temperature β is calculated using

body system in its thermal equilibrium atrandom state Φ as

an inverse temperature is one of the most=A

important problems in computational
physics. It has, however, two major

tr[e − β H A]
=
Z

Φ AΦ
Z

where

computational difficulties. The number ofZ ( β )= tr[e− β H ]=quantum states to be included and theΦΦdimension of each quantum state vector

is the partition function and the double

are both exponentially increases as thebracket indicates statistical average overtemperature and the system size increases,

random states. The thermal state isrespectively. Each difficulty has been

defined by the imaginary time evolution

independently attacked using randomof the random state.

state method [1, 2] and variational

Φ(β ) = e− β H / 2 Φwave function, respectively.

Thermal energy E and specific heat

In this project, we introduced

C are expressed as

random state method to general form ofvariational wave functions. Only few

E (β ) =papers have investigated the applicationof random state method to variationalH

andwave functions [3-7].

2C (β ) kB β 2  E 2 − E  .
=



273

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

Variational wave function is a point
on a manifold embedded in the Hilbert
space,

ψ ( w) = ∑ c( w, σ) σσwhich is parameterized with variationalparameters, w = {wi } , whose number ismuch less than the dimension of theHilbert space.

Fig.1 : Specific heat per site for 14

Thermal state has been introduced insite spin-1/2 antiferromagnetic

some variational wave functions such as

Heisenberg chain with PBC

MPS and PEPS as

calculated with RBM( α = 4 ) and1024 samples.

Ψ ( β ; w) = e − β H /2 Ψ ( w)Referenceswhere Ψ ( w) is random variational[1]

wavefunction.

T. Iitaka, and T. Ebisuzaki, Phys.

Rev. Lett. 90, 047203 (2003).

In this project, the above[2]

formulation was extended to general form

T. Iitaka, and T. Ebisuzaki, Phys.

Rev. E 69, 057701 (2004).

of variational wavefunction by[3]

K. Takai et al., J Phys Soc Jpn 85,

calculating the imaginary time evolution

034601 (2016).

with stochastic reconfiguration[3, 8]. It is[4]

noted that the trace is formally taken over

S. Garnerone, T. R. de Oliveira, andP. Zanardi, Phys. Rev. A 81, 032336 (2010).

the entire Hilbert space in our method[5]

T. Iitaka, arXiv:2006.14459 (2020).

while it was limited to its subspace in the[6]

A. Iwaki, A. Shimizu, and C. Hotta,

previous method[3]. A numerical

Physical Review Research 3, L022015 (2021).

example for Restricted Boltzmann[7]

Machine is computed with NetKet [9]S. Goto, R. Kaneko, and I. Danshita,

Phys. Rev. B 104, 045133 (2021).

and the result is shown in Fig. 1.

[8]

S. Sorella, Phys. Rev. Lett. 80, 4558

(1998).
[9]

G. Carleo et al., SoftwareX 10,100311 (2019).

274

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

Tensor data analysis by tensor network
Kenji Harada
Graduate School of Informatics, Kyoto University
Kyoto 606-8501, Japan

Recently, data in various ﬁelds have been
stored as a multi-dimensional array called atensor. When the number of attributes (indexes) of a tensor increases, the size of thetensor becomes huge because it increases exponentially as the number of indexes increases.
Therefore, we need the high-speed and highquality use of tensor data of huge sizes.
A tensor network (TN) is an eﬀective tool for
representing a huge composite tensor and has
been extensively developed in quantum information and statistical physics research. Since
the tensor network can eﬀectively represent
such high-dimensional data, we may eﬃciently
process the general tensor data.
Previous studies [1, 2] proposed using ofa tree tensor network (TTN) for a multidimensional probability distribution. In thisproject, we try to improve the performance ofa TTN algorithm for generative modeling tomulti-dimensional data distribution.
Here, we consider a generative model basedon a quantum state [1, 2]. Following the Bornrule, we deﬁne p(x) as the square of the amplitude of a wave function:p(x) =|ψ(x)|2,Z(1)

no loop structure. Thus, an MPS is a specialized type of TTN. As shown in [2], the performance depends much on the network topology
of TTNs.
We propose a new algorithm to change
the network structure of a TTN for generative modeling. Several strategies can beused to select a better choice of local network
structure[3]. Our network optimization improves the performance of generative modelingwith a TTN.

References[1] Zhao-Yu Han, Jun Wang, Heng Fan, Lei
Wang, and Pan Zhang. Unsupervised Generative Modeling Using Matrix ProductStates. Phys. Rev. X, 8(3):031012, 2018.
[2] Song Cheng, Lei Wang, T. Xiang, and Pan
Zhang. Tree tensor networks for generativemodeling. Phys. Rev. B, 99(15):155131,
2019.
[3] Kenji Harada, Tsuyoshi Okubo, and Naoki
Kawashima. Network optimization of tree
generative models. in preparation.

where ψ(x) is a wave function and Z is thenormalization factor.
MPSs [1] and balance TTNs [2] have been
proposed to deﬁne the wave function for generative modeling. The number of indexes of atensor in the balance TTN is equal to that in
the MPS. The only diﬀerence is the topology ofthe network; all physical indexes in the TTN,
xi , are connected, and the TTN network has

275

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

Molecular dynamics simulation of BaTiO3 nanostructure II
T. Hashimoto
National Institute of Advanced Industrial Science and Technology (AIST),
Tsukuba Central 2, 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, JapanThe ferroelectric domain wall (FDW) has dielectric, electrical, and mechanical properties
that diﬀer from those of bulk. It can evenbe considered a two-dimensional phase diﬀerent from the bulk. In our previous study, we
found that 71◦ and 109◦ Ising-like FDWs areobserved in BaTiO3 nanoclusters[1]. In thisstudy, we focused on the FDW structures andenergies of rhombohedral BaTiO3 .
We used a shell model proposed by Tinte etal.[2]. The calculations were performed withthe code developed by us. We performed calculations for 71◦ , 109◦ , and 180◦ {110} FDWs in
the rhombohedral phase. We used 10×6×6 5atom unit cells for the 109◦ FDWs, and 8×6×6
and 10×6×6 10-atom building blocks for the71◦ and 180◦ {110} FDWs, respectively. The
initial Ising FDW conﬁgurations were made byshifting the Ti atoms in one direction in halfof the supercell, and in the other direction inanother half. The initial Bloch FDW conﬁgurations were made in the same way as the Isingcase except that the local dielectric polarization vectors were rotated across the FDW. To
mimic the isolated FDW between very thickdomains in the supercell calculations, the inplane lattice constants were set to be the sameas those of the single domain phase, while theout-of-plane lattice constant was relaxed.
Thedomainwidthsandenergies had been calculated with the
Ginzburg–Landau–Devonshire (GLD) modeland the GLD model with reduced gradientterms (GLD-r) [3]. The GLD and the GLD-r

results predict that the 109◦ Bloch wall ismore stable than the 109◦ Ising wall energy,
which is not in accordance with the LDA
results. The GLD-r model predicts narrower
walls and lower wall energies compared to the
GLD model, which improves the energy andwidth of the 71◦ Ising wall and the width of
the 109◦ Ising wall, but overshoots the energyof the 109◦ Ising wall and the energy andwidth of the 109◦ Bloch wall. The GLD-r
model predicts that the 180◦ Bloch wall isunstable compared to the 180◦ Ising wall,
which is not in accordance with the LDA
results.
The shell model results also have inconsistencies compared to the LDA results: The
FDW energy of the 109◦ Ising wall is too low.
The FDW energy of the 180◦ {110} Bloch wall
is too high and the width of it is too narrow.
However, it predicted that the 71◦ Ising wall
has the lowest energy, followed by the 109◦
Ising wall and the 109◦ Bloch wall.

References[1] T. Hashimoto and H. Moriwake, Physica
B, 656, 414768 (2023).
[2] S. Tinte, M. Stachiotti, S. Phillpot,
M. Sepliarsky, D. Wolf, and R. Migoni: J.
Phys.: Condens. Matter 16 (2004) 3495.
[3] M. Taherinejad, D. Vanderbilt, P. Marton, V. Stepkova, and J. Hlinka: Phys.
Rev. B 86 (2012) 155138.

276

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

Numerical study of quantum many-body scars inbosonic systems
Ryui KANEKO
Department of Physics, Kindai University
Higashi-Osaka, Osaka 577-8502, Japan

Ĥ = Ĥ0 + Ĥint ,X †Ĥ0 = −J(âi âi+1 + â†i+1 âi ),

(1)
(2)

i

Ĥint =UXn̂i (n̂i − 1).
2

6
S vN/ ln 2Understanding how isolated quantum manybody systems thermalize is of great interestin physics. The eigenstate thermalization hypothesis (ETH) provides a sufficient condition
for thermalization and is satisfied in many
quantum systems in the presence of interactions [1, 2]. However, recent experimentsof nonintegrable systems suggest ETH breaking. This phenomenon originates from the special eigenstates exhibiting extremely slow thermalization [3, 4]. Such quantum many-bodyscar (QMBS) states exist in various spin models [5–10]. It is desirable to search for the
QMBS states in other quantum systems, which
can be realized in experiments.
Here we focus on the Bose-Hubbard modelsuitable for experiments using ultracold atoms
on optical lattices. The Hamiltonian is given
as

N =8
Scars

4

U/J = 1

2
0
−10

−5

0

5
Energy/J

10

15

Figure 1: Entanglement entropy as a functionof energy. The QMBS states are circled. We
consider the system size L = 8 and the interaction strength U/J = 1 in Eq. (1). The
quantum number sector with the particle number N = L and the even parity I = +1 isgiven by blue dots. Each QMBS state (a red
circle with a dot) with the particle number
N has the energy U N/2 and the expectation
PPL
2
values ⟨ L
i=1 n̂i ⟩ = N and ⟨ i=1 n̂i ⟩ = 2N .
Its bipartite entanglement entropy at N = L,corresponding to a state |Sn=L/2 ⟩, becomesvN → [ln(πL/8) + 1]/2 for L → ∞.
SA

(3)

i

We specifically choose the model in a truncatedHilbert space with the maximum occupation
number nmax = 2. Note that three-body lossesof atoms on optical lattices realize this constraint [11, 12].
We numerically diagonalize the Hamiltonianand find the eigenstates with a small entanglement entropy (see Fig. 1). We also analytically show that they are generated by applying the SU(2) ladder operator consisting ofa linear combination of two-particle annihilation operators. The construction of the QMBS
states is equivalent to the case of the S = 1 XY
model [13, 14].
These findings will stimulate further research on the QMBS states in general BoseHubbard models without any constraints,
which would be more suitable for experimentsof ultracold atoms in optical lattices. The construction of the QMBS states in general spin

277

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

S systems has been discussed recently [15, 16],
which could be extended to the Bose-Hubbard
model by utilizing the improved HolsteinPrimakoff transformation [17, 18].

References[16] L.-H. Tang et al., Phys. Rev. Research 4,043006 (2022).
[17] M. Vogl et al., Phys. Rev. Research 2,043243 (2020).
[18] J. König and A. Hucht, SciPost Phys. 10,
7 (2021).

[1] W. Beugeling et al., Phys. Rev. E 89,
042112 (2014).
[2] H. Kim et al., Phys. Rev. E 90, 052105(2014).
[3] H. Bernien et al., Nature 551, 579 (2017).
[4] D. Bluvstein et al., Science 371, 1355
(2021).
[5] C. J. Turner et al., Nat. Phys. 14, 745(2018).
[6] C. J. Turner et al., Phys. Rev. B 98,
155134 (2018).
[7] A. J. A. James et al., Rev. Lett. 122,
130603 (2019).
[8] N. Shibata et al., Phys. Rev. Lett. 124,
180604 (2020).
[9] M. Serbyn et al., Nat. Phys. 17, 675(2021).
[10] S. Moudgalya et al., Rep. Prog. Phys. 85,
086501 (2022).
[11] A. J. Daley et al., Phys. Rev. Lett. 102,040402 (2009).
[12] M. J. Mark et al., Phys. Rev. Lett. 108,215302 (2012).
[13] M. Schecter and T. Iadecola, Phys. Rev.
Lett. 123, 147201 (2019).
[14] S. Chattopadhyay et al., Phys. Rev. B
101, 174308 (2020).
[15] N. O’Dea et al., Phys. Rev. Research 2,043305 (2020).

278

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

Tensor-network study of the quench dynamics of theKitaev honeycomb model
Ryui KANEKO
Department of Physics, Kindai University
Higashi-Osaka, Osaka 577-8502, Japan⟨i,j⟩γ ,γ∈{x,y,z}energy density

-0.06

D=2
-0.07 D=3
D=4
-0.08 D=5
-0.09 D=6
D=7
-0.1 D=8
-0.11
-0.12
-0.13
0.1

1
time

10

Figure 1: Time evolution of the energy density.
absolute value of z-component magnetization

The Kitaev honeycomb model [1] in themagnetic field displays a rich phase diagramof quantum phases [2, 3]. These findings stimulate further research on the nonequilibrium
dynamics of the model. One of the motivations for studying dynamics is an experimental
interest; the Kitaev candidate materials may
exhibit unconventional transport [4–6]. The
other is an interest in statistical mechanics as
to how the integrable Kitaev model, which israre in quantum 2D systems, thermalizes in thepresence of additional interaction [7, 8]. These
interests motivate us to study the dynamics
of the simple Kitaev model with the magnetic
field. As a benchmark, we specifically focus
on the sudden quench of the ferromagnetic Kitaev model starting from the antiferromagnetic
state without the magnetic field. Previous investigations are limited to mean-field calculations or numerical-diagonalization studies forsmall systems. We use the 2D tensor-network
method using the infinite projected entangledstate (iPEPS) [9, 10], which can directly deal
with inifinite systems accurately, to examine
how far one can go for the quench dynamics ofthe Kitaev model.
We prepare the staggered antiferromagnetic
state, whose spins are oriented along the zaxis, on a honeycomb lattice as the initial state|ψ(0)⟩. The state can be prepared by the bonddimension D = 1, and is time evolved by theS = 1/2 Kitaev HamiltonianX(1)Ĥ = −J γ Ŝiγ Ŝjγ ,1
0.8
0.6 D=2
D=3
0.4 D=4
D=5
D=6
0.2 D=7
D=8
0
0.1

1
time

10

Figure 2: Time evolution of the absolute valueof the z-component magnetization.
where γ corresponds to the bond direction andJ γ is the Kitaev interaction. Here we consider the isotropic and ferromagnetic interaction J γ = 1. The wave function at eachtime |ψ(t)⟩ = exp(−itĤ)|ψ(0)⟩ is calculated by
the simple update algorithm [11, 12] using thetensor-network library TeNeS [13]. The wave
function is optimized up to the bond dimension
D = 8 to examine how physical quantities converge by improving the accuracy of the state279

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

with increasing the bond dimension.
The energy should be conserved for thequench dynamics. We can estimate to what extent the simulation goes well by looking at thedegree of conservation. In the present study,
the energy is nearly conserved for tJ γ ≲ 4when D ≳ 4 (see Fig. 1). Within this reliable
time frame, the magnetization nearly vanishes
at tJ γ ≈ 4. This observation is consistent withthe previous study, where an emergent nonmagnetic steady state has bee proposed [7].
Our study would serve as a reference for understanding the quench dynamics of the Kitaev model and relevant unconventional steadystates, which can be identified by utilizing the
2D tensor-network method. It is left for future
work to investigate the effect of the magnetic
field in the model.

[11] J. Jordan et al., Phys. Rev. Lett. 101,
250602 (2008).
[12] H. C. Jiang et al., Phys. Rev. Lett. 101,
090603 (2008).
[13] Y. Motoyama et al., Comp. Phys. Commun. 279, 108437 (2022).

References[1] A. Kitaev, Ann. Phys. (NY) 321, 2(2006).
[2] Z. Zhu et al., Phys. Rev. B 97, 241110(R)(2018).
[3] C. Hickey and S. Trebst, Nat. Commun.
10, 530 (2019).
[4] J. Nasu and Y. Motome, Phys. Rev. Research 1, 033007 (2019).
[5] T. Minakawa et al., Phys. Rev. Lett. 125,047204 (2020).
[6] M. Udagawa et al., Phys. Rev. Lett. 126,
127201 (2021).
[7] L. Rademaker, SciPost Phys. 7, 071(2019).
[8] G.-Y. Zhu and M. Heyl, Phys. Rev. Research 3, L032069 (2021).
[9] T. Nishino et al., Prog. Theo. Phys. 105,409 (2001).
[10] F. Verstraete and J. I. Cirac, arXiv:condmat/0407066.

280

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

The Study on Superionic Conduction and Phase Transition by Disordering ofComplex Hydride
Ryuhei Sato
Advanced Institute for Materials Research,
Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, JapanThe complex hydrides such as LiCB9H10

order

are promising candidates for electrolytes dueofmagnitudeataroundroom

temperature [2].

to their high ionic conductivity and stability
to

metal

anodes[1].

Althoughthedevelopment of such complex hydrides is
progressing

rapidly,themechanismofsuperionic conductivity induced by the phasetransition due to the disordering (≈ the change
in rotational motion) of the complex anions is
not clear. Therefore, the purpose of this study
is to clarify the effect of the disordering on
the ion conduction and the phase transition by

Figure 1 Arrhenius plot of Li self-diffusion

using ab initio molecular dynamics (AIMD)

coefficient during NNP-MD simulations ofand neural-net potential molecular dynamics

LiCB9H10 with and without CB9H10− rotation(NNP-MD) simulations.

[2].

AIMD simulation at 1000K confirmed

As shown in the discussion above, it isthat LiCB9H10 exhibits concerted motion, in

possible to quantitatively analyze the changes

which

migratein ion transport caused by anion-rotation

simultaneously with the rotation of theusing NNP-MD simulations. However, it iscomplex anion. To clarify the effect of the

difficult to obtain a specific relationship

CB9H10− rotation, we analyzed the change inbetween the anion rotation and Li migration.

diffusion coefficient with and without the

Therefore, in order to visualize the ion

rotation using NNP-MD simulations (Fig. 1).

transport mechanism in such a disordered

As shown in the figure, the activation energy

system, we have developed an ion-migration-

increases when the anion rotation is stopped,mechanism analysis method using topologicaland the conductivity decreases by nearly one

data analysis based on persistent homology.

multiple

Li

ions

281

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

This year, we investigated the applicability ofmolecules improves ionic conductivity is also

topological data analysis on ion migrationunclear. Therefore, we created a workflow tomechanism [3]. Here, we choose α-AgI as aanalyze the effect of such neutral moleculesbenchmark, because this system retains theon complex hydrides. Given a composition

ordered I bcc lattice while Ag ions areincluding neutral molecules, complex anions,disordered in the system and shows complex

and cations, a structure search method based

concerted ion migration. MD simulations at

on genetic algorithms (USPEX) can estimate

various temperatures confirmed specific ring

the stable and reasonable structure including(one-dimensional hole in terms of homology)

neutral molecules after several hundred DFT

structures appear only when the system shows

calculations. The activation energies of ion

high ionic conductivity. Inverse analysis of

migrationthese ring structures revealed that they were

simulations using the obtained structures are

four-membered rings consisting of two Ag

in good agreement with experimental values.

and I ions. By deforming these four-

The obtained trajectories of ion migration

membered rings along the probability density

during the metadynamics simulation confirm

distribution of the ring structure obtained bythat the coordination structure of the neutral

topological data analysis, we succeeded inmolecule

visualizing the representative concerted ion

significant effect on the activation energy. As

migration of Ag ions in α-AgI[3]. In thea result of machine learning using thefuture, we plan to apply this technique toestimatedcomplex hydride systems to clarify theconfirmed that the distance between anions

relationship between anion rotation and ion

and the electronegativity are important ones

migration.

for the activation energy [5].

estimated

duringion

activationfrommetadynamics

migrationenergies,

has

it

a

was

References

Recently, complex hydrides are expected
to be applied as multivalent-ion (ex. Zn2+ and[1] W. Tang et al., Adv. Eng. Mater., 6 (2016)Mg2+) conductors [4]. Although the insertion

1502237.

of neutral molecules into complex hydrides is[2] R. Sato et. al., to be submitted.

known to be effective to improve the ion[3] R. Sato et al., J. Chem. Phys., 158, (2023)

conductivity

144116.

for

these

multivalent-ion

systems, there is no clear guideline to choose[4] K. Kisu et al., J. Mater Chem. A, 10

most effective neutral molecules for these(2022) 24877.

complex hydrides. In addition, the specific[5] S. Egon, R. Sato et al. submitted to NPJ.

mechanism of how the insertion of neutral

Comput. Mater.

282

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

Theoretical studies on dynamics in biomolecular and liquid
systems
Shinji Saito
Graduate University for Advanced Studies (Sokendai) / Institute for Molecular Science,
38 Nishigo-Naka, Myodaiji, Okazaki 444-8585, Japan
Fluctuations are inherent in condensed
phase

systems,such

asliquidattempt to understand the conformational

andorigins of such dynamical disorder.

biomolecular systems, and occur acrossdiversespatio-temporalscales.

In the second half of the 2022 fiscal year,

These

we constructed a Kv.12 model of about 140k

fluctuations play a significant role in structuralatomschanges and reaction dynamics, thereby

concentration of 0.15 M and performed

influencing various properties and functions.

roughly 100 μs molecular dynamics (MD)Consequently,ofsimulations at an external electric field of 500

fluctuations is crucial for a comprehensivemV using GROMACS. During the MD

understanding

relevant

simulation, about 2300 potassium permeation

chemical/biological processes. In the second

events were observed. The time series analysis

half of the fiscal year of 2022, our group

showed that these permeation events roughly

utilized the supercomputer resources at ISSP

followed the pattern of a Poisson process withfor mainly two related projects, 1) the iona slight deviation from it. Further analyses andpermeation process of potassium channel

MD simulations at different ion concentrations

Kv1.2, and 2) the excitation energy transfer ofand external electric field strength are to belight harvesting complex II (LHCII) in higher

performed for a complete picture of theplants.

dynamic disorder in the ion permeation

anofunderstanding
theKv1.2 is a mammalian voltage-controlled

at

the

electrophysiologicalprocess of Kv1.2.

potassium channel from the Shaker K+ channel

LHCII, where the first step of thefamily which regulates the cellular membrane

photosynthetic process occurs, has a complex

potential by conditionally releasing thestructure which makes its excitation energy

intracellular K+ ions. Since its crystal structuretransfer (EET) difficult to be elucidated

was resolved about twenty years ago, it has

experimentally. In order to understand the EET

become a typical and representative model for

dynamics of LHCII, we performed quantum

the study of ion channels. Here, we focus onmechanics (QM) calculation using the CAM-

the dynamical disorder revealed in the ion

B3LYP functional and determined the optimal

permeation processes under various ion

range separate parameters (between CAM andconcentrations and external voltages and

B3LYP) that reproduced the experimental
283

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

excitation energies of Chl a and b in threefunctional theory. The parallelized MMSIC-

solvents [1]. With these optimized parameters,

MD simulation sufficiently reproduced thewe refined the structures for 14 chlorophylls

QM/MM excitation energies on the fly and(Chl) in LHCII by using a modified version of

showed that the excitation energies of Chl a in

GAMESS on the supercomputer at ISSP. In

LHCII varies wider than in solvents.

order to study the dynamical properties, we
introduced

a

model

calledReferences

Molecular

Mechanics Shepard interpolation Correction[1] Zhu, Zhe, Masahiro Higashi, Shinji Saito,(MMSIC) that can give the approximate

Excited states of chlorophyll a and b inexcitation energies during MD simulation

solution by time-dependent density functional

based on the interpolation of optimized

theory, The Journal of Chemical Physics

structures

156(12) (2022) 124111.

using

the

above-mentioned

284

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

Structural formation of inverse patchy particles
Takamichi TERAO
Department of Electrical, Electronic and Computer Engineering,
Gifu University, Yanagido 1-1, Gifu, 501-1132
Inverse patchy colloids are those in which theexhibit a variety of configurations by changing

anisotropy of interparticle interactions is due tothe ratio of positive and negative charges on theinhomogeneous surface charges [1]. Previous

surface.

studies of colloidal particle systems have largely

Figure 1 shows a snapshot of inverse patchy

dealt with the problem of aggregates, which arecolloids when the motion of each sphere wasdominated by attractive interparticle interactions,confined between a pair of two-dimensionaland colloidal crystals, which are dominated byplanes. From these results, we have revealed therepulsive interactions. Conventional patchy

dependence of the orientation order of theparticles are dominated by anisotropic attractive

colloids on the charge coverage ratio, as well asinteractions. Inverse patchy particles, on thethe phase diagram of the system.

other hand, consist of heterogeneously charged
units where the pair potentials between particles
contain both attractive and repulsive properties.
They are expected to exhibit novel properties
and their properties with respect to structureformation are of interest.
In this study, we focused on the structuralformation of inverse charged patchy particles.
The surface of these particles is charged and theFig. 1: Snapshot of inverse patchy colloids

sign of the charge is reversed between theunder geometrical confinement.

bipolar part and the other parts. The structureformation in these systems was studied bymolecular

dynamics

simulations.

ReferencesThe

simulations were performed with the LAMMPS[1] P. D. J. van Oostrum, M. Hejazifar, C.

simulation package in the NVT ensemble. In

Niedermayer, and E. Reimhult: J. Phys.:thermal equilibrium, inverse patchy particles

Condens. Matter 27 (2015) 2341.

285

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

Effect of chain crossing prohibition on phase-separated
structure of block copolymers
Katsumi HAGITA,1 and Takahiro MURASHIMA2
1

Department of Applied Physics, National Defense Academy, Yokosuka 239-8686
2

Department of Physics, Tohoku University, Sendai 980-8578

We have studied phase-separated structures

To confirm fundamental behaviors usingof block copolymers in order to clarify

Gaussian chains, we performed large scale

topological effect of chain crossing prohibition.

DPD simulations to estimate lamellar domain

To achieve this aim, we firstly considered MP-

spacing of symmetric linear, ring, and four-

SRP extensions [1,2] of dissipative particle

arm-star block copolymer blends [7]. We also

dynamics (DPD) simulations [3]. We found that

confirmed practical compatibility of lamellara slight change in the potential (force field)

domain spacing between self-consistent field

between original DPD [3] and MP-SRP

theoryextensions [1,2] affects weakly the lamellar

fundamental understandings of topologicaldomain spacing. The numerical precision iseffects by chain crossing prohibition, we

related to potential validity problems, andperformed CGMD simulations of polymers

comparisonsphysics,

with rings. In [9], ring-filling effect on stress–

Gaussian and excluded volume chains, may notstrain curves of randomly end-linked tetra-arm

be appropriate. Thus, we switched our approach

prepolymers

to the development of force field parameters in

discovered topological transition in multicyclic

order to treat directly with the phase separation

chains with structural symmetry inducing

phenomenon in the coarse-grained molecularstress-overshootdynamics (CGMD) simulations of the Kremer-

multicyclic/linear

Grest

Althoughelongational flow. In addition, we investigatedimplementations of the MP-SRP extensions

topological barriers on the crystallization[1,2] were developed using LAMMPS [5], thebehavior of ring polyethylene melts with trefoil

CGMD approach has the advantage of beingknots via united atom MD simulations [11].

(KG)

betweendifferent

model

[4].

free from the implementation and able to use

andCurrently,

DPD

was

simulations

studied.

[8].

In

[10],

phenomena

weblends

we

in

under

investigatedForbiaxialinteraction

HOOMD-blue [6] on the GPU. Notably,parameters of AB di-block copolymers for

HOOMD-blue has a good performance for the

CGMD simulations. To consider a connection

standard DPD simulations [3].

to Flory-Huggins χ parameter, we are studying

286

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

relationships among interaction parameters,lamellar structure as shown in Fig. 1 (b). We

lamellar domain spacing, and chain length (N).

also confirmed that the lamellar domainFirst, we performed CGMD simulations of AB

spacing depends on the interaction parameters.

di-BCP starting from the lamellar structures asAftershown in Fig. 1 (a). We confirmed theparameters, we will investigate topologicalbehaviors of domain spacing from weak toconstrain effect of polymers such as rings via

strong segregation conditions were consistent

CGMD simulations.

furtherstudiesontheinteraction

with theories, simulations, and experiments.

References[1] R. D. Groot and P. B. Warren: J. Chem. Phys.
107 (1997) 4423.

[2] N. Iwaoka, K. Hagita and H. Takano: J. Chem.
Phys. 149 (2018)114901.
[3] K. Hagita, T. Murashima, H. Shiba, N. Iwaoka,
and T. Kawakatsu: Comput. Mater. Sci., 203(2022) 111104.
[4] K. Kremer and G. S. Grest: J. Chem. Phys. 92(1990) 5057.

[5] S. Plimpton: J. Comp. Phys. 117 (1995) 117.
[6] J. A. Anderson, J. Glaser, and S. C. Glotzer,
Comput. Mat. Sci. 173, 109363 (2020).
[7] K. Hagita, T. Murashima, and T. Kawakatsu,
Macromolecules 55, 8021 (2022).
[8] K. Hagita, and T. Murashima, Polymer 269,
125733 (2023).
[9] K. Hagita, T. Murashima, T. Ohkuma, andH. Jinnai, Macromolecules 55, 6547 (2022).
Fig. 1: Snapshots of CGMD models of AB di-[10] T. Murashima, K. Hagita, and T.

BCP forming lamellar structures.

Kawakatsu, Macromolecules 55, 9358 (2022).
[11] K. Hagita, T. Murashima, N. Sakata, K.

Second, we found that the AB-BCP in

Shimokawa, T. Deguchi, E. Uehara, and S.
Fujiwara, Macromolecules 56, 15 (2023).

CGMD simulations spontaneously forms a287

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

Analysis of Various Crystal Forms of Glycine Nanocrystals
by Large-Scale Metadynamics Simulations
Hiroki NADA
Faculty of Engineering, Tottori University, 4-101 Koyama-Minami, Tottori 680-8552

The control of structure and shape offunctions represented by Gaussian window

glycine nanocrystals is important in connection

functions as collective variables indicated a few

with pharmaceutical engineering [1]. In thislocal

project, a classical molecular dynamics (MD)corresponded to a crystal-like structure with asimulation combined with the metadynamics

different molecular arrangement of glycine(MTD) method (hereafter, MTD simulation)

molecules.

was performed to elucidate stable structuresnanocrystal was also different between thoseand shapes of a glycine nanocrystal.

local minima. These results suggest that notminima.

Each

Moreover,

local

the

shape

minimumofthe

The MTD simulation was performed for a

only the structure of nanocrystals but also the

cubic system in which an α glycine nanocrystal

shape of them can be analyzed using the MTD

consisting of 64 glycine molecules was placed

method.

at its center. The remainder of the system was

In summary, the MTD method is a helpful

filled with water molecules. AMBER force

tool to search for stable structures and shapes offield

nanocrystals consisting of complex molecules.

was

used

for

estimation

of

the

intermolecular interaction energy and force
acting between glycine molecules. The MTD

References

simulation was carried out with DL_POLY 2.20

[1]

[2],

CrystEngComm. 17 (2015) 723.

in

which

PLUMED

1.3

[3]

was

X.

Yang

and

A.

S.

Myerson:implemented to permit combination with the[2] W. Smith and T. R. Forester: J. Mol. Graph.

MTD method, using System C of ISSP

14 (1996) 136.

Supercomputer Center.

[3] M. Bonomi et al.: Comput. Phys. Commun.

A free energy surface (FES) obtained using

180 (2009) 1961.

discrete Cα-Cα and N-N radial distribution

288

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

Structural Analysis of Amorphous Silicon Monoxide by aMetadynamics Method with Unsupervised Machine Learning
Hiroki NADA
Faculty of Engineering, Tottori University, 4-101 Koyama-Minami, Tottori 680-8552

Amorphous SiO has been studied by many

Supercomputer Center. In this project, we

research groups because it can be functional

focused on the structural similarity between thematerials

molecular-scaleexamined amorphous phases. The structural

structure of it remains elusive. The purpose ofsimilarity was analyzed by evaluating thethis project was elucidating molecular-scale

geometrical similarity of Si-Si, Si-O, and Si-O

structure of amorphous SiO using a classicalpair-distribution functions. In this project, themolecular

simulationgeometrical similarity was evaluated using a

combined with the metadynamics (MTD)

method of unsupervised machine learning,

method (hereafter, MTD simulation [2]).

dimensional reduction technique [5].

[1].

However,dynamics(MD)

MTD is an enhanced sampling method that

We confirmed that the present methodology

increases the probability of reaching high free-

is helpful to investigate detailed structures ofenergy states by adding a Gaussian bias

multicomponent amorphous phases.

potential to the Hamiltonian of a state. The
MTD method affords a free energy surface

References(FES) in a space of collective variables (CVs)[1] A. Hirata et al.: Nat. Commun. 7 (2016)

for a system of interest. This FES quantitatively

11591.

represents[2] A. Laio and M. Parrinello: Proc. Natl. Acad.

therelative

thermodynamic

stabilities of many different states of the system.

Sci. U. S. A. 99 (2002) 12562.

The MTD simulation was carried out for[3] W. Smith and T. R. Forester: J. Mol. Graph.

each of amorphous Si, SiO and SiO2 phases

14 (1996) 136.

with DL_POLY 2.20 [3], in which PLUMED[4] M. Bonomi et al.: Comput. Phys. Commun.

1.3 [4] was implemented to permit combination

180 (2009) 1961.

with the MTD method, using System B of ISSP

[5] H. Nada: ACS Omega 3 (2018) 5789.

289

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

Quantum algorithm for combinatorial optimization
Tatsuhiko SHIRAI
Department of Computer Science and Communications Engineering, Waseda University
Ookubo, Shinjuku-ku, Tokyo 169-8555⃗ ⃗γ ⟩ := e−iβp H e−iγp X · · · e−iβ1 H e−iγ1 X |+⟩ ,
|β,
(1)
⃗
where β = (β1 , . . . , βp ) and ⃗γ = (γ1 , . . . , γp ).
p denotes a circuit depth. X is called a∑mixer and adopted as X = i σix where ⃗σi =(σix , σiy , σiz ) is the Pauli operator acting onsite i. |+⟩ is deﬁned as σix |+⟩ = |+⟩ for all⃗ ⃗γ ) toi. The classical computer optimizes (β,⃗ ⃗γ | H |β,⃗ ⃗γ ⟩. We solve
minimize the value of ⟨β,
the quantum dynamics on the super computer
using openMP.
Maximum graph coloring problem colors
each node of a given graph G = (V, E) using ncolors so that the number of edges connecting
the nodes with diﬀerent colors is maximized.
This problem has a constraint that one color
is assigned to each node. When solving a COP
with constraints by QAOA, the Hamiltonian isgiven as H = Hobj + Hcst . Here, Hobj and Hcst

denote the Hamiltonians for objective function and constraints of the COP, respectively.
There are three methods to give the Hamiltonians: penalty method, domain-wall method,and spin-variable reduction method.
Figure 1 shows p dependences of successprobability ps . Here, we set n = 4 and |V | = 6,
and generate 10 random graphs with edge density 0.75. We found that domain-wall method
and spin-variable method outperforms penalty
method.
1

Success probability ps

We have developed a quantum algorithm for
QAOA (Quantum approximate optimization
algorithm), which solves combinatorial optimization problems (COPs) in gate-based quantum computers. This year, we proposed a new
method called spin-variable reduction method
for formulating a constrained COP as an Isingmodel [1] and compared the performance withconventional methods in maximum glaph coloring problems.
QAOA is a hybrid quantum-classical algorithm designed to approximate the groundstate of a Hamiltonian H. QAOA consists ofcalculations in a quantum computer and classical computer. The quantum computer calculates the expectation value of H over a variational state, which is given by

0.5

Penalty
Variable reduction
Domain wall
0

0

2

4

6

8

p

10

12

14

16

Figure 1: Performance comparison among
penalty method (squares), spin-variable reduction method (circles), and domain-wall method(triangles).

References[1] T. Shirai and N. Togawa: IEEE Transactions on Computers (Early Access).

290

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

Topological chirality order and its stability inbreathing-kagome antiferromagnets
Kazushi AOYAMA
Department of Earth and Space Science, Graduate School of Science, Osaka University
Machikaneyama-cho, Toyonaka-shi, Osaka 560-0043

The quest for topological spin textures insolids is one of the central issues in condensed
matter physics because of their possible applications to spin-electronic devices. A typicalexample of such topological spin textures is
the so-called skyrmion crystal (SkX) where thetotal solid angle subtended by all the spin in
the magnetic unit cell takes an integer topological in units of 4π. Noting that the scalarspin chirality χijk = Si · (Sj × Sk ) is related tothe solid angle, the SkX could be regarded asa topological chiral order. Although the SkX
is usually realized in the presence of a magnetic ﬁeld being irrespective of whether thethe Dzaloshinskii-Moriya interaction is active
or not [1], a zero-ﬁeld SkX is also possible in a
speciﬁc metallic system [2]. Previously, we theoretically demonstrated that a breathing bondalternation of the underlying lattice can induce
a zero-ﬁeld topological chiral order which canbe viewed as a miniature version of the SkX[3]. In the work, we investigate eﬀects of an additional single-ion magnetic anisotropy whichmore or less exists in real materials.
We consider the following classical J1 -J3model on the breathing kagome lattice withthe single-ion magnetic anisotropy:
H = J1∑Si · Sj + J1′⟨i,j⟩S+J3∑⟨⟨i,j⟩⟩∑

Si · Sj⟨i,j⟩LSi · Sj + D∑(Siz )2 , (1)iwhere Si is a classical spin, the summation∑⟨⟩S(L) runs over site pairs on small (large)

triangles having the nearest neighbor (NN) interaction J1 (J1′ ), J3 > 0 is the third NN antiferromagnetic interaction along the bond direction, and D < 0 (D > 0) represents an easyaxis (easy-plane) anisotropy. In the isotropic
case of D = 0, a triple-Q state appearing fora relatively large J3 takes the miniature SkX
structure having a nonzero uniform total chirality [3]. By performing Monte Carlo simulations, we investigate the stability of this topological chiral order against the anisotropy D.
In our simulations, 1×106 - 2×106 MC sweeps
are carried out under the periodic boundary
condition, and the ﬁrst half is discarded forthermalization. Si at each site is updated ina conventional random and a successive overrelaxation-like processes, where in the former(latter) process, we try to rotate a spin ina randomly proposed direction (by the angle
π around the local mean ﬁeld) by using theMetropolis algorithm. Observations are doneat every MC sweep, and the statistical average
is taken over 4 independent runs starting from
diﬀerent random initial conﬁgurations.
It is found that the zero-ﬁeld topological chiral order with the miniature SkX structure isrelatively robust against easy-axis (D < 0)and easy-plane (D > 0) anisotropies, although
for strong easy-axis (easy-plane) anisotropy, a
triple-Q collinear (coplanar) state is favored.
Since the SkX and anti-SkX each having positive or negative chirality are energetically degenerate, the topological Hall eﬀect of alternative sign is possible at zero ﬁeld. For moderate291

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of TokyoD < 0 (D > 0), the collinear (coplanar) phasepreempts the SkX phase and in the chiral sector, it shows a random domain structure consisting of positive- and negative-chirality clusters [4]. It is also found that the miniature SkX
is robust against the external magnetic ﬁeld,
but further analysis is necessary to clarify the
in-ﬁeld properties of the model.

References[1] T. Okubo, S. Chung, and H. Kawamura,Phys. Rev. Lett. 108, (2012) 017206.
[2] R. Ozawa, S. Hayami, and Y. Motome,
Phys. Rev. Lett. 118, 147205 (2017).
[3] K. Aoyama and H. Kawamura, Phys. Rev.
B 105, (2022) L100407.
[4] K. Aoyama and H. Kawamura, J. Phys.
Soc. Jpn. 92, (2023) 033701.

292

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

Molecular modeling of ferroelectric nematic phase
Takeaki ARAKI, and Matheus DE MELLO
Department of Physics,
Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo, Kyoto 606-8502
Ferroelectric liquid crystals are partially

The distribution function suggests the molecule

ordered anisotropic fluids with ferroelectricity.

shape can be regarded as a rod without theIts experimental realization has not beenremarkable head-tail asymmetry. Owing to thereported for a long time. Recently, someelectric charges, the local parallel-parallel

experimental

inalignment is favored than the parallel-anti

synthesizing molecules which possibly show

parallel alignment, which is usually observed inferroelectric nematic order [1-2].

the conventional nematic phase.

groupshave

succeeded

Even though the molecular features required
to stabilize the ferroelectric phase weresystematically studied, the origins of the
spontaneous polar ordering remain still unclear.
In this work, we study a molecular liquid DIO,
which were synthesized by Nishikawa et al.,focusing on which of the dipole-dipole
electrostactic interaction or the polar shape ofthe molecules is dominant for the polar ordering.
We performed all-atom molecular dynamics
simulations of 1000-2000 DIO molecules with

Fig. 1: (a) A snapshot of the liquid crystal phase

GROMACS. The partial charge on each atom isof DIO at 330K. The molecule is represented

estimated by using Gaussian. In order to see thewith a rod whose color shows the orientation. (b)effect of the electrostatic interaction, we carried

The distribution function. (c) The correlationout the simulations of the DIO without thefunction.

electric charges. We analyzed the nematic, polar

References

order parameters, the spatial distribution, and the[1] H. Nishikawa, et al. Adv. Mater. 29,

correlation between the molecular axes.

1702354 (2017).

Our

simulation

results

indicate[2] R. J. Mandle, et al., Chem. Eur. J. 23, 14554

the(2017).

spontaneous polar ordering at a low temperature.

293

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

Aggregation of polyelectrolytes in mixed solvents
Takeaki ARAKI, and Dai Matsubara
Department of Physics,
Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo, Kyoto 606-8502
Polymer systems often exhibit phase-

homogeneously. As its composition is decreased

separated structures with diverse patterns from(top), the system is phase-separated and a

microscopic to macroscopic length scales. In

labyrinth pattern is formed. In the polymer-richparticular,

phase, a stripe concentration modulation pattern

diblock

copolymer

melts

and

polyelectrolyte in poor solvents are expected tois observed.

be applied to nanostructures.
In this study, we consider mesoscopic patternformation of weakly charged polyelectrolytes inbinary mixture solvents by means of field
theoretic polymer simulations in lattice space.
The polymer is modeled as a continuous
Gaussian chain, and is described by path
integrals.

The

solvent-polyelectrolyte

andsolvent-solvent interactions are given by Flory
Huggins

parameters.

Theelectrostaticinteraction is obtained by solving PoissonBoltzmann equation. The binary solvent itself ismixed. But the affinities of the solvents to thepolyelectrolyte is different. One of the solvents
prefers the polyelectrolyte more than the other.

Fig. 1: Snapshots of the phase-separated patternFigure 1 shows the concentration field of theof the polyelectrolyte in binary mixture solvents.

polyelectrolyte in the binary solvent mixtures.

The concentration of the more preferred

The composition of the more preferred solvent

component is increased from top to bottom.

is increased from top to bottom with fixing theReferences

polyelectrolyte concentration. When the more

[1] G. H. Fredrickson, The Equilibrium Theoryof Inhomogeneous Polymers, Oxford
University Press, 2013.

preferred component is rich in the solvent(bottom),thepolyelectrolyte

is

mixed

294

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

Rheology Simulation of Dynamically Cross-linked Networks
Yusuke Yasuda
Research Center for Computational Design of Advanced Functional Materials.

National Institute for Advanced Industrial Science and Technology, Ibaraki 305-8568.
Recently a new type of dynamic bondelastomers called entropy-driven has gained
attention for showing superior self-healing andtoughness compared to conventional dynamic
bond elastomers [1], but detailed understanding is
not reached yet. Last year, we succeeded infabricating a new coarse-grained model which
can reproduce typical hysteresis and self-healing
properties of conventional dynamic bondelastomer. This year, we introduced entropydriven properties to the coarse-grained model.
In this work, we made coarse-grained modelfor dynamic bond elastomers based on KremerGrest model [1]. For cross-links between reactive
bead, we introduced harmonic potential withassociating and dissociating condition (Fig. 1).
We tuned harmonic constant K and bond fractureenergy εd. Under this condition we conductedequilibration simulation and mechanical tests.
Calculation was conducted using molecular
dynamics simulation software packages OCTA,
and LAMMPS on ISSP Supercomputer System.

Fig. 2, the relationship among Bond lifetime,
dissociation energy, and harmonic constant K
We also conducted lifetime analysis fordifferent Ks. Then we found that average bondlifetime follows Arrhenius equation depending on
K (Fig. 2), and dynamic bonds with smaller K
shows shorter lifetime in high εd region,
concluding that entropy-driven dynamic bondelastomers with high εd owing to smaller
dissociation activation energy.
We also conducted mechanical and stressrecovery test. Both models show typical large
hysteresis, and stress recovery at hightemperature. Then we analyzed stress-relaxation
behavior by analyzing relationship between stressand healing time, to find that stress-relaxation
times are closely related to bond lifetime. Then
we conclude that faster self-healing in entropydriven dynamic bond elastomers originates fromfaster dissociation / association dynamics, whichcomes from lower activation energy fordissociation [3].
References [1] Kim et al. Macromolecules
2020, 53, 4121. [2] K, Kremer et al., J. Chem.
Phys., 1990, 92, 5057. [3] Y. Yasuda et al.,under review.
Acknowledgement This work is partially supported
by a project JPNP18016, commissioned by the New
Energy and Industrial Technology Development
Organization(NEDO)

Fig.1, Schematic illustration of the coarsegrained model of dynamic bond elastomers.
From equilibration simulation, we found thatequilibrium constant Ka has proportional
dependence to 1/T, consistent to Gibbs equation.
From the slopes and y-intercepts of the 1/T
dependence of Ka, we defined ΔSa and ΔUa. From
this analysis, we found that ΔUa are influenced
primarily by εd, and ΔSa are predominantly
influenced by on K. As K decreases, ΔSa
increases, confirming smaller K represents
entropy-driven systems.

295

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

Structural Analysis and Mechanical Properties Simulations ofPolymer Networks with Various Network Structures
Yusuke Yasuda
Research Center for Computational Design of Advanced Functional Materials.

National Institute for Advanced Industrial Science and Technology, Ibaraki 305-8568.
Rubber elasticity theories proposed so far
qualitatively reproduces the mechanicalproperties of vulcanized rubbers and crosslinked gels, but it is not quantitative. One reasonof this is that the actual polymer networks
contain inhomogeneity, and it is difficult to
detect it by experiments. Then theoretically,
homogeneous network is assumed.
In this work, we aimed to develop amethodology to handle network inhomogeneity,
focusing on detecting structural defects in the
network. We performed coarse-grained
molecular dynamics simulations for networks
with various inhomogeneity. Then we analyzed
the structure during the crosslinking process andthe mechanical properties under uniaxial
deformation. We tried to extract substructuresthat contribute to the mechanical properties.

inhomogeneity, we terminate the reaction
calculation when reaction ratio reaches certain
value in range of 0.24 – 0.96. Then we
conducted uni-axial elongation test (λ=1-4) , andobtain stress-extension ratio curves. Then we
calculated elastic modulus Gsim by fitting theinitial shape of stress-extension ratio curves.
70 x 10-3

Arm length = 5 σ
3-arm
4-arm
Arm length = 10 σ
3-arm
4-arm

3

Gsim [ ε/σ ]

60
50
40
30
20
10
0
0

5

10

15

20

3

Ncycle, eff [ /σ ]25

30 35 x 10-3

Fig. 2 Elastic moduli Gsim against the numberof effective closed loops Ncycle,eff.
Then we tried to calculate the structure which
contribute to elastic moduli. Based on the
Scanlan-Case [2,3] criterion, we analyzed the
network structure. We counted only effective
chain s between branching points with three ormore branches, excluding dangling, loop, superdangling, super-loop structure and eliminate,including bridge and super-bridge structures. In
the network constructed by effective chains, wecounted the number of closed loops in network
based on phantom network theory, finding thatelastic moduli are proportional to Ncycle,eff for all
branch num and arm length.
References 1. 1. Fujiyabu et al., Sci Adv. 2022,
8(14), eabk0010 2. Scanlan, J. Polym. Sci., 1960,
43, 501, 3. Case, J. Polym. Sci., 1960, 45,
Acknowledgement This work is supported by aGrant-in-Aid for Young Scientists (22K14740)

Fig, 1. Schematic illustration of thefabrication method of the system.
In this research, we fabricated coarse-grainedmodel of star-shaped end-crosslinked elastomers
which is known as a homogeneous network, wascreated based of the repulsive Kremer-Grest
model. We prepared 3-arm and 4-arm systems
referring to tetra-PEG, tri-PEG gel [1] tuning thearm length. Two terminal types A, B are
prepared. Total number of beads in system are
set approximately 1.0×105.τ density conditionis set as 0.85σ-3Then system relaxation
calculation and reaction calculations. For
fabrication of networks with various

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

Theory of microwave scattering in a dissipative
superconducting circuit element
TSUYOSHI YAMAMOTO
Faculty of Pure and Applied Sciences, University of Tsukuba
Tsukuba, Ibaraki 305-8571==

Figure 1: (top) A Josephson junction with theJosephson energy EJ and the charging energy
EC shunted by a resistor with the resistanceR. (bottom) Illustration of the phase diagramof the Schmid-Buldagaev transition.
A Josephson junction is a building blockof superconducting circuits. Thanks to its
characteristic nonlinearity, we can develop
a great variety of superconducting devices,such as a superconducting qubit and a parametric amplifier. Although the property of
the isolated Josephson junction is understoodwell, it is little known, especially its manybody eﬀects, when it contacts electromagnetic environments. The most notable phenomenon is a dissipation-induced quantum
phase transition. Forty years ago, Schmidand Buldagaev theoretically found that a single Josephson junction shunted to a resistor, a resistively-shunted Josephson junction,

undergoes a superconductor-insulator transition [1, 2]. When the resistance is larger thana resistance quantum (R > RQ = h/4e2 ),the Josephson junction becomes an insulator,
while when R < RQ , it becomes a superconductor, interestingly not being controlled byjunction parameters (see Figure 1).
Although the theoretical works on the DC
transport checked the Schmid-Bulgadaev transition [3, 4], there is no decisive evidence forthe experimental observation [5], stirring up
controversy. The recent development of superconducting circuits allows experimentalists tomeasure the AC-impedance of the resistivelyshunted Josephson junction using microwave
spectroscopy and then to obtain important
properties, such as the excitation spectrum or
the dynamical property. This endeavor would
be a key to solving the ongoing debate problem.
In this project, we numerically investigatedthe AC-impedance of the resistively-shunted
Josephson junction [6]. To obtain the ACimpedance, we calculated the phase-phase correlation function defined on the imaginarytime axis using the path-integral quantum
Monte Carlo (PIQMC) simulation and thenperformed the analytical continuation for the
PIQMC data using the Padé approximation.
Since the Padè approximation needs highlyaccurate calculations for the applying data,
we performed the parallel computing for theMonte Carlo sampling mainly using “F2cpu”,
succeeding in collecting a large enough Monte297

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[4] S. L. Lukyanov and P. Werner, J. Stat.
Mech. P06002 (2007).
[5] A. Murani et al., Phys. Rev. X 10, 021003(2020).
[6] P. Werner and M. Troyer, Phys. Rev. Lett.
95, 060201 (2005).
[7] T. Yamamoto, T. Kato, L. I. Glazman,
and M. Houzet, in preparation.

Figure 2: Frequency dependence of the realpart of the AC-impedance Z ′ (ω) for various
values of α = RQ /R = 2 (red), 1 (green), and
1/2 (blue), obtained by the PIQMC. At theDC limit (ω → 0), the AC-impedance goes to
0 and R for α > 1 and α < 1, respectively,
which is a signature of the Schmid-Bulgadaev
transition. The dashed lines represents the analytical formulas.
Carlo steps, typically ∼ 5 × 107 , to provide theaccurate calculations for the AC-impedance.
Our numerical results for the AC-impedanceas a function of frequency or Matsubara frequency (the Fourier space of the imaginarytime) confirmed the superconductor-insulatortransition, as originally predicted by Schmidand Bulgadaev and revealed the low- and highenergy physics and their crossover (see Figure 2) [7]. This work would give a deep understanding of the dissipative Josephson junction,
which pushes further development for superconducting devices using Josephson junctions.

References[1] A. Schmid, Phys. Rev. Lett. 51, 1506
(1983).
[2] S. A. Bulgadaev, JETP Lett. 39, 264(1984).
[3] N. Kimura and T. Kato, Phys. Rev. B 69,
012504 (2004).

298

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

Grain boundary Li ion conductivity in NASICON-type Solid
Electrolytes
Masanobu NAKAYAMA
Department of Advanced Ceramics,
Nagoya Institute of Technology, Gokiso, Showa-ku, Nagoya, Aichi 466-8555

LiZr2(PO4)3 (LZP) related materials arecomposition, but (a) with Ca ions randomly

promising solid electrolyte materials with higharranged and (b) with Ca ions segregated at theionic conductivity and chemical stability.[1]grain boundary, as experimentally confirmed.

Many material calculations have been applied

For these models, molecular dynamics was

to lithium ion-conducting NASICON-type solid

performed to evaluate the diffusion coefficientelectrolytes, but most of them are concerned

versus temperature, and the results are shown inwith bulk conduction, and grain boundaries

Figure 2 as an Arrhenius plot. The obtained

have rarely been investigated.[2], [3] In thisresults show that the ionic conductivity (roomstudy, the effect of dopant Ca on the ionictemperature) of the segregated structure (b) is 5

conductivity is discussed.

x 10-4 cm2 s-1, which is about twice theTo calculate lithium ion diffusion in the

diffusion coefficient of the random structure (a).

grain boundary model of LZP, we used theOn the other hand, no significant change in

Bond Valence Force Field (BVFF) proposed by

activation energy was observed. This suggests

Adams et al.[4] plus a Stillinger-Weber (SW)that the grain boundary conductivity can be

type

controlled by controlling the arrangement of

force

field

representing

three-bodyinteractions. Parameters included in the force

dopants.

field were evaluated by a cuckoo search to

A more detailed analysis of the obtained results

match the results of first-principles calculations.

shows that Li ions have a trapping effect near

The same was employed: force-field molecular

Ca ions. In previous studies, the grain boundary

dynamics (FFMD) calculations using the NVT

structure also has the effect of trapping Li ions.

ensemble were performed for 1 ns every 50 K

Therefore, it is suggested that Li ions become

from 400 K to 900 K. The results are presented

excessively

in Table 1. The software package NAP[5], [6]boundary. We believe that such excess of

was used for the above.

lithium ions contributes to the enhancement of

Figure 1 shows the grain boundary model usedionic conduction at the grain boundary.

in this study. (a) and (b) have the same299

concentrated

near

the

grain

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

Figure 1 Grain boundary structure model for Ca-doped LZP. (a) Random substitution of Ca to Zr
sites and (b) Ca segregation model near grain boundaries.

Figure 2 Comparison of diffusion coefficients of Li in two-types of grain boundary models of Cadoped LZP. (Black symbol: random Ca arrangement model; Green symbol: Ca segregation modelnear to grain boundary.

References[4] S. Adams et al., Phys. Stat. Sol. A, 208,

[1] M. Nakayama et al., Chem. Commun., 58,

1745-1753 (2011)

9328–9340, (2022)[5] R. Kobayashi, J. Open Source Softw.,6,[2] K. Nakano et al., J. Phys. Chem. C, 125,

2768, (2020)

23604–23612 (2021)[6] R. Kobayashi et al., APL Mater., 8, 081111,

[3]

R.

Kobayashi

et

al.,

Acta

Mater.,(2020)

226,117596 (2022)

300

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

Surface Tension Stabilizing Cell Motion Order
Katsuyoshi Matsushtia, Naoya Kamamoto, Maki Sudo, Taiko Arakaki
Depertment of Biological Science, Graduate School of Science, Osaka University
Machikaneyama, Toyonaka, Osaka 560-0043

Cells collectively move in particular scaffolds in various biological phenomena, including organ morphogenesis [1]. For example,
neuronal cells move by interacting with surrounding glial tissue[2]. The collective motion
reflects the interaction between the tissue andthe surrounding scaffold. The ordering mechanism in this collective motion needs sufficientexamination through the computer simulation
based on statistical physics. In particular, theessential interaction in statistical physics is thesurface tension between the cells and scaffolds,
which should be incorporated into the simulation. To clarify the effect of the interface
tension, we consider a collective motion of acell cluster and its dependence on the surface
tension.

This result explains the previous result of
the cluster consisting of the two types of cells,
which is reported in the previous project [6].
The cell cluster transitions from the rotational
to straight motion as the differentiation of cells
progresses. The present result implies that thismotion transition originates from the surface
tension decreases associated with the cell typechanges in the differentiation.
We appreciate Prof. H. Kuwayama, Prof.
S. Yabunaka, Dr.
H. Hashimura, Prof.
Sawamoto, Dr. Matsuda, Prof. H. Yoshino,
and Prof. M. Kikuchi their various support on
this work.

We employed the cellular Potts model [3]and extended its cells with the polarity of cellcell adhesion [4]. This model induces the spontaneous order of the cellular motion and successfully explains the collective migration ofcells [5]. We changed the surface tension between the scaffold and the cells and observed
the collective motion. In the case of low surface
tension, the cells have an affinity with the scaffolds and individually exhibit a random walk.
When the surface tension increases above theGraner-Glazier surface tension, the cells forma cluster and collectively move, In this case,the collective motion has a persistence time,in which the cluster almost straightly moves.
As the surface tension increases, the switchingfrom the straight motion to the rotational motion and inverse switching appears. In the caseof large surface tension, the collective motion
becomes the rotation, the center of which corresponds to the center of the cluster. Due tothe rotational motion, the center of the cluster
exhibits only diffusion.

[1] P. Rørth, Ann. Rev. Cell Develop. Biol.
25, 407 (2009).

References[2] C. Nakajima, M. Sawada, K. Sawamoto,Curr. Opi. Neurobio. 66, 1 (2021).
[3] F. Graner and J. Glazier, Phys. Rev. Lett.
69, 2013 (1992).
[4] K Matsushita, Phys. Rev. E 95, 032415(2017).
[5] K Matsushita, Phys. Rev. E 97, 042413(2018).
[6] K. Matsushita, N. Kamamoto, M. Sudo,
K. Fujimoto, Proc. Sympo. Traffic Flow
Self-driven Parti. 27, 23 (2021).

301

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

Calculation of ordered structures, dynamics and opticalproperties of soft materials
Jun-ichi FUKUDA
Department of Physics, Kyushu University, Motooka 744, Nishi-ku, Fukuoka 819-0395
We continued our study that started in theelectromagnetic field of light. We employ

previous year on the structure of a cholestericplane-wave expansions in the directions normal

blue phase liquid crystal (BPLC) under anto the thickness direction, and the Maxwell

electric field. Our study was motivated by aequations were discretized so that the matrix

recent experimental study [1] demonstrating therepresenting the Maxwell equation is of blockformation of various (meta)stable ordering withtridiagonal form.

different periodicity in a reproducible manner

Our calculations [2] showed that structures

by applying an electric field. Particularly, awith different and discretized periodicity are

tetragonal structure with I4122 symmetry was

indeed (meta)stable. The free energy evaluation

shown to be (meta)stable even after theof these structures revealed that the periodicity

cessation of the electric field, which had notofbeen reported.

systematically on the electric field strength. We

the

most

stablestructuredepends

Our numerical study focused on the role ofalso calculated [2] the reflection spectra of theplanar surfaces that confine BPLC on its

structures. We found strong selective reflection

structure and stability, and was based on thedepending on the circular polarization of the

Landau-de Gennes continuum theory thatincident light arising from the chirality ofdescribes the orientational order by a second

BPLC, and the variation of the reflection peakrank tensor. Calculations were carried out on a

position in response to the electric field.

regular

mesh,

which

facilitatedtheparallelization (in our case OpenMP was used).

References

We also calculated the optical properties of the[1] D.-Y. Guo et al., Nature Mater. 19, 94

structures obtained (more specifically, how(2020).

monochromatic incident light with a given[2] J. Fukuda, Front. Soft Matter 2 1011618

wavevector is reflected and transmitted) by(2022).

directly solving the Maxwell equations for the302

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

Measurement-induced phase transitions inU (1)-symmetric monitored circuits
Yohei FUJI
Department of Applied Physics, University of Tokyo, Tokyo 113-8656

Many-body quantum systems evolved under repeated measurements have recently beenshown to reveal a novel type of nonequilibrium phase transition called the measurementinduced phase transition [1, 2]. A paradigmatic example is the entanglement transition from a volume-law to an area-law entangled phase, which is induced by the competition between unitary dynamics and projective measurements. Surprisingly, the entanglement transition exhibits critical phenomena described by conformal field theory as inequilibrium phase transitions. It is then natural to ask whether universality classes due
to symmetry or dimensionality also arise inmeasurement-induced phase transitions.

ground state of one-dimensional quantum systems preserving U (1) symmetry. This result thus highlights the importance of symmetry for classifying and understanding themeasurement-induced phase transitions.

In this study, we considered one-dimensionalrandom quantum circuits consisting of twosite unitary gates and projective measurementsof single qubits, both of which are design topreserve a U (1) charge corresponding to thetotal number of flipped qubits [3]. By employing parallel sampling implemented in Julia on Supercomputer System B at ISSP, we
computed steady-state physical quantities averaged over many circuit realizations and measurement outcomes. In addition to the entanglement transition characterized by the logarithmic scaling of the entanglement entropy,
we found a charge-fluctuation transition characterized by the logarithmic scaling of the bipartite charge fluctuation as shown in Fig. 1.
The latter transition exhibits critical phenomena akin to the Tomonaga-Luttinger liquid
theory, much like critical phenomena for the

Figure 1: Subsystem charge correlation function ⟨nA nC ⟩c between two disjoint subsystems
A and C is plotted as a function of the measurement rate p. It is related to the bipartite
charge fluctuation FA by ⟨nA nC ⟩c = FA +FC −FA∪C .

L=8

0.50

L=12
L=16
L=20

0.20

L=24

0.10 0.8
0.05
0.02

0.6
0.4
0.2
0.0

0.01
0.00

- 0.5

0.0

0.05

0.5

0.10

1.0

0.15

0.20

0.25

0.30

References[1] A. C. Potter and R. Vasseur, Entanglement Dynamics in Hybrid Quantum Circuits in Entanglement in Spin
Chains: From Theory to Quantum Technology Applications, edited by A. Bayat,
S. Bose, and H. Johannesson (Springer
International Publishing, Cham, 2022)
pp. 211–249.

303

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[2] M. P. A. Fisher, V. Khemani, A. Nahum,and S. Vijay, Random Quantum Circuits,
Annu. Rev. Condens. Matter Phys. 14,335 (2023).
[3] H. Oshima and Y. Fuji, Charge fluctuation and charge-resolved entanglement ina monitored quantum circuit with U (1)symmetry, Phys. Rev. B 107, 014308(2023).

304

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

First principles studies of electronic structuresin armchair-edged graphene nanoribbons
Junhuan LI, Kenta ARIMA
Graduate School of Engineering, Osaka University, Yamada-oka, Suita, Osaka 565-0871
Local electronic structures of graphene

The band structure of 39-AGNR (where 39 is

nanosheets have been intensively studied so far

the width measure by the number of rows of[1, 2]. However, the relationship between the

carbon atoms across the AGNR) is shown in theelectronic structure and atomic arrangement isright panel of Figure 1. Based on the simulated

still unclear. In this study, we simulated scanningband structure, the STM images were calculated

tunnelingofusing different sample biases. At a sample bias

nanoribbonsof -1 V, the simulated STM image exhibit

microscopy

armchair-edged(STM)graphene

images(AGNRs) by the first principles calculations.

hexagonal patterns as shown in Figure 2.

Calculations were performed by STATE(Simulation Tool

for AtomThis work has been performed on System B

TEchnology)of the Supercomputer Center, the Institute for

package. The AGNR we simulated was

Solid State Physics, the University of Tokyo.

terminated by hydrogen atoms as shown in
Figure 1. Periodic boundary conditions wereapplied to a supercell setup with length 𝑙.

Fig. 2: The simulated STM image of 39-AGNR.

References[1] T. Wassmann et al., J. Am. Chem. Soc. 132,
Fig. 1: The ribbon configuration and band

3440 (2010).

structure of 39-AGNR.

[2] J. Li et al., Phys. Rev. B 103, 245433 (2021).

305

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

Numerical studies on thermodynamics and excitation
spectra in Heisenberg antiferromagnets on thetriangular and kagome lattices
Yoshiyuki Fukumoto
Department of Physics, Faculty of Science and Technology, Tokyo University of Science
2641 Yamazaki, Noda, Chiba 278-8510

Neutron scattering experiments on twodimensional frustrated spin systems with spontaneous magnetizations, such as a triangularlattice system Ba3 CoSb2 O9 [1] and a kagomelattice system Cs2 Cu3 SnF12 with the q =0 type ordering stabilized by DzyaloshinskiiMoriya (DM) interactions [2,3], have revealed
that the dynamic structure factors consist of
not only collective magnon excitations but also
spinon continuums.
As for the triangular-lattice system, Zhang
and Li used the fermionic spinon operator to
calculate dynamic structure factors within the
RPA approximation [4]. They argued that thestructure of the spinon continuum is closely
related to the roton minimum of the magnon
spectrum, and also found an amplitude mode
composed from two spinons that cannot be
captured by the spin wave approximation.
As for the kagome-lattice system, magnon
dispersions obtained from a neutron scattering experiment on Cs2 Cu3 SnF12 were analyzed
based on the linear spin wave (LSW) approximation, and reported that the spin-wave interactions led to a large negative-renormalization,
where the exchange parameter was found to
be about 40% smaller than a reliable value,
J = 20.7 meV, determined from magnetic susceptibility measurements [2]. Also, a recent
neutron scattering experiment has observed afourth collective excitation that is thought tobe an amplitude mode [3].
In this article, we report our recent results of an Ising expansion (IEP) study on
the large negative-renormalization observed in

306

Cs2 Cu3 SnF12 . In order to stabilize the q = 0
ordering, which is assumed in the xy plane,
we incorporate the DM interaction with the Dvector of D = (0, 0, D∥ ) and write a minimal
Hamiltonian for this material as
H=J

∑
⟨i,j⟩

Si · Sj +∑D · (S i × S j ),

(1)

⟨i,j⟩

where S i represents the S = 1/2 operator atsite i. Introducing a perturbation parameter
λ, we write Hλ = H0 + λH1 , where H0 is
an Ising model and we assume λ → 1 recovers eq. (1). The perturbation calculation was
performed using the connected cluster expansion method, including all clusters up to 9 links
and 10 sites. Multi-block diagonalization was
performed in the calculation of the eﬀective
Hamiltonian [5]. We also incorporated a local field term to improve convergence of theseries [5]. The local field term is the form of
staggered magnetic field stabilizing the q = 0
ordering and its strength is defined as t. We
include the local field term into H0 and thatwith opposite sign into H1 , so as to cancel out
with each other at λ → 1. If we set t = 0,then we have unphysical results such as themagnon energies being negative in most of the
wave vectors. When t > 0, the results change
significantly at first, and the region of negative
magnon energy decreases. When t increases tosome extent, the region of negative magnon energy disappears and the result does not depend
on t very much. As such value of t, we choose
t/J = 1 in this study. As t is further increased,
the resulting t dependence increases slowly, be-Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Figure 1: Renormalization factors rk,α (λ).
The black and red lines, respectively, represent
the ninth and eighth order results.

Figure 2: Comparison of magnon dispersions
between the IEP and LSW. The choice of line
type is the same as in Fig. 1.

cause it becomes harder to correctly handle thelocal field term included in the perturbation.
Let ϵk,α be the magnon excitation energy
to be obtained, where k is the wave vector
and α is the identifier for three modes in the
kagome lattice. We set ϵk,α = rk,α (1) × ϵLSWk,α ,blue lines are, respectively, the ninth-order IEP
and LSW results, where D∥ /J = 0.12. On
the other hand, the red lines are the results of
the LSW with renormalized parameters whereD∥ /J = 0.18 and the energy scale is modified
by a factor 0.62 [2]. Figure 2 tells us that thedispersion relations calculated by the LSW undergo not only renormalization of the interaction parameters but also flattening of the highenergy part and the appearance of roton-likeminimum around M point due to spin-wave interactions. As for the triangular lattice system,
we had modification of high-energy part but norenormalization of interaction parameters [5].
In summary, we have used the IEP to investigate how spin-wave interactions modify the
LSW predictions of magnon dispersions. We
have successfully calculated the renormalization factor rk,α (λ) by the IEP. As a result, inthe kagome lattice system, in addition to therenormalization of the interaction parameter
pointed out by Ono et al. [2], we have foundthat the high-energy part around M point undergoes strong deformation. In the near future, we plan to incorporate more realistic DM
interactions to reproduce the detailed structure of experimental magnon spectrum data.

where ϵLSWk,α is the magnon excitation energy inthe LSW and rk,α (λ) is “renormalization factor”. We calculate rk,α (λ) by the IEP up to
the ninth order, and evaluate value of rk,α (1)by naive summation of the obtained series. In
Fig. 1 we show the eighth and ninth order results of rk,α (λ) at Γ and M points as a functionof λ, and find a good convergence. In particular, for M point, the renormalization factor
represented by the solid line amounts to 0.26 at
λ → 1, which demonstrates that large negative
renormalization occurs in the kagome system.
Experimental results for Cs2 Cu3 SnF12 showthat the two high-energy modes are almost degenerate at the Γ point, and their energies areabout 10 meV [2,3]. By fixing J = 20.7meV,
which has been determined from magnetic susceptibility measurements [2], and searching thevalue of D∥ so that the calculated energy ofthe high-energy modes at Γ point is 10 meV,
we obtain D∥ /J = 0.12.
Figure 2 shows a comparison of the magnon
dispersions between the IEP and the LSW.
The path on the wave-vector plane, which connects Γ and M points, is chosen in the same
way as in the experiment [3]. The black and307[1] S. Ito et al., Nat. commn. 8, 235, 2017.
[2] T. Ono et al., JPSJ 83, 043701, 2014.
[3] M. Saito et al., PRB 105, 064424, 2022.
[4] C. Zhang et al., PRB 102, 075108, 2020.
[5] W. Zheng et al., PRB 74, 224420, 2006.

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

Numerical study of the eﬀect of large bondrandomness on the magnetization process andmagnetic susceptibility in the spherical Kagome
system {W72V30}
Yoshiyuki Fukumoto
Department of Physics, Faculty of Science and Technology, Tokyo University of Science
2641 Yamazaki, Noda, Chiba 278-8510

The spherical kagome system {W72 V30 } is afrustrated cluster magnet consisting of thirty
V4+ with S = 1/2 [1]. The exchange interaction has been estimated as J = 115K. In
the experimental study of the low-temperaturemagnetization process, a step structure uniqueto cluster magnets has not been found [2]. To
understand this result, the eﬀects of bond randomness [2] and Dzyaloshinskii-Moriya (DM)interactions [3] were examined.
Subsequently, specific heat measurements
were performed by Kihara et al. at temperatures T < 11K and magnetic fields up to
15T [4]. Although the magnitude of the specific heat is comparable to that expected in
the Heisenberg model, the shape of the specific heat curve is completely diﬀerent. In particular, in the Heisenberg-model calculation, apeak appeared near 2K, and it was found thateven a magnetic field up to 15T aﬀected theshape of the specific heat curve. However, inthe experiment, the specific heat curve had nopeak and was not aﬀected by the magnetic field
up to 15T. To resolve this problem, the specific heat calculation was performed considering bond randomness or DM interaction. Eachstrength is assumed to be about 0.1J. However, it did not lead to an explanation of theexperimental results.
In a previous study investigating the eﬀect ofbond randomness in a triangular lattice system[5], it was reported that when the randomness
is large compared to the energy scale of the308

magnetic field, the magnetic field dependenceof the specific heat disappears. Therefore, we
decided to investigate the eﬀect of strong bondrandomness also in the spherical kagome system. As a result, it was clarified that the experimental results can be explained by assuming that the strength of the bond randomness
is 40-50% of the exchange interaction. The results were published as a paper [6].
It is also interesting to study if the model
Hamiltonian with the large bond randomness
and a modest DM interaction reproduces experimental magnetic susceptibility. Magnetic
susceptibility calculation in the presence of theboth terms is a challenging problem in termsof computational complexity. If we have DM
interactions, we need to work with vectors ofdimension 230 . In addition, when DM interaction exists, the magnetic susceptibility is obtained by applying a small magnetic field andcalculating the expected value of magnetization by the cTPQ method. Then, it is necessary to average the direction of the magnetic
field, the bond randomness, and the random
vector of the cTPQ method. We will continue our research on eﬃcient implementationof these averaging operations.
[1] A. M. Todae et al., ChemComm, 3351 (2009).
[2] J. Schnack et al., arXiv:1304.2603v1.
[2] Y. Fukumoto et al., JPSJ 87, 124710 (2018).
[4] T. Kihara et al., PRB 99, 064430 (2019).
[5] K. Watanabe et al., JPSJ 83, 034714 (2014).
[6] M. Motohashi et al., PTEP, 113I01 (2022).

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo                            !                                                                             !"        
 #  !  #  " $                       %    &    #  %                            '         (               # !' ' #   '"  ! '  ) '"        '    '  % ' * &                &   +                                            (                     ,                 -    ,  &           .        /    / 0 !//"   1                 //   ,                                 &      ,                                                  //           //        //  
2            -   //    ,3       






















                      &            .         
4  -       !2"             !5"                      !6"   7 .                                  
 !2"   -             8  











  











 9  




 

 



 


9

   




 

 

 
     %


















9









9









    
99#    99#     ' '  '  %  ' '  ' # 
 ' '                          (                                           &        +                    -         :             -    ,                     #          '












 



 



 

 





 







 





309















 

 





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

1.0
ΔAF

F
nTLL

SD

UUDD

0.0

XY1

nTLL
F
−1.0
−0.2

+  ;   ' 


−0.1



J2

0.1

<       ' =*  % ' >





0.0

  3  ,               ! '"  !  " 4        4?  ::??
   ? &      & .                  ,   !A <"  B     
 3  ,          A <0    # ! "! 9 B" # ! 9 B" 
# ! "  ! '"   ! '"    //                       ,      3  ,      ,      '           '          3  ,          # ! "# ! " 
# ! "  !  "   ! '" /      3  ,              +      //   ::??        ! "  ! '"&     .         +             !  # "!  '  ' '=C* ' '"!     "   +  //                       -                     *  .                       -             .          7  $.  7   @  
D 4.      8 .  7 $.  4 6
+     4.  E A  4 E  'B''
!'>" F / 5,  .  4 A  EGA  H'!=C"1 E I 7  ,  ? E    EA  6  > !=C)") 7 @   2 7 , E A  2 C)B !==>"
B I A @    A  2  *H !=CH"* 7 $.     7 @   
4.  A     2   I   A  4   E 4   * ' ' =*  % ' >                          +                !+"    !4?"   
!::??"            /     0  !//"   &      ,      //           !  '"    //       +                         (       !  "   !  "           
 8          
            
 !'      "     , +       !   "              
                ! "                   (           
3  ,        >      -   ,          3   
  ! "   +      ! '"      &          3  ,                                           
-    3  ,      +           4?      // ::??     5,  . 7  ,      &             @  7 , )      










310







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

Phase separation dynamics in hydrodynamic models
Kyosuke ADACHI

Nonequilibrium Physics of Living Matter RIKEN Hakubi Research Team, RIKEN
Center for Biosystems Dynamics Research (BDR), 2-2-3 Minatojima-minamimachi,
Chuo-ku, Kobe 650-0047, JapanRIKEN Interdisciplinary Theoretical and Mathematical Sciences Program (iTHEMS),
2-1 Hirosawa, Wako 351-0198, Japan
Liquid-liquid phase separation (LLPS),

Slyozov-Wagner (LSW) scaling for the growthwhich is widely recognized in condensedof the cluster size. Then, as a simple extension

matter and soft matter physics, has gained

to nonequilibrium conditions that can be

renewed attention in cell biology. Several kinds

relevantof protein condensates have been found in aconsidered a model that breaks the fluctuation-

cell, and LLPS is discussed as a possible

dissipation relation. This model is similar tomechanism of the formation and dynamics ofActive Model B+ [3], which have beencondensates [1]. To understand how the amino

originally proposed as a model of active matter.

acid sequence of proteins and intracellular

Performing numerical simulations with small

nonequilibrium environments can affect theadditional terms to improve numerical stability,phase

we found that the cluster growth is slower thanseparation

simulationsofdynamics,
polymer

numerical

modelstointracellularsituations,wethat in Model B, suggesting the suppression ofandthe phase separation dynamics due to thehydrodynamic models should be useful.

nonequilibrium conditions.

This year, as a preliminary study, we
performed several numerical simulations ofhydrodynamic models, which are PDE models

Referencesof the particle density field. We first tested the[1] S. F. Banani et al., Nat. Rev. Mol. Cell Biol.

Model B dynamics [2] by the Euler-Maruyama

18, 285 (2017).

method with discretization of time and space.

[2] P. C. Hohenberg and B. I. Halperin, Rev.

Using parallel computing to collect many

Mod. Phys. 49, 435 (1977).

samples with different initial conditions and[3] C. Nardini et al., Phys. Rev. X 7, 021007

noise realizations, we confirmed the Lifshitz-(2017).

311

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

Solution of the Eliashberg equation with theelectron-electron Coulomb interaction
Ryosuke AKASHINational Institutes for Quantum Science and Technology
The Eliashberg theory[1] with the Migdal
approximation [2] is the state-of-the-art framework for the first-principles calculations ofphonon-mediated superconductors. In principle, the competition between the phononmediated pairing attraction and screened
Coulomb repulsion is accurately treated there
for weakly correlated metals, which is crucial
for reliable estimation of the superconducting
transition temperature Tc . However, the latter part has been mostly treated as a single
empirical parameter µ∗ [3], with which one cannot apply the theory to undiscovered systems
where the Tc measurement does not exist. It
is only recent that the effect of the screened
Coulomb interaction on the pairing has been
explored non-empirically.
The author has recently studied the effect of
frequency dependent screened Coulomb interaction within the random -phase approximation. He has revealed the cooperative mechanism of the plasmonic and phononic pairing
attractions from first principles, by analysing
the non-diagonal (pairing) part of the self energy [4, 5]. On the other hand, a part of diagonal self energy due to the plasmons were shownto have negative effect, partially cancelling theformer pairing [6]. However, the full treatment of all the components of the screened
Coulomb-induced frequency dependent self energy was missing. In this project, we calculated
the Eliashberg equations with all the screened
Coulomb self energy parts included for uniform
electron gas [7].
The Eliashberg equations with the Migdal

approximation are written as follows:
Zi (iωj ) = 1 − T

X ωj 0 Zi0 (ωj 0 )i0 j 0χi (iωj ) = T

ωj Θj 0 (iωj 0 )X ξi0 + χi0 (ωj 0 )i0 j 0φi (iωj ) = −T

Θj 0 (iωj 0 )X φi0 (ωj 0 )i0 j 0Θj 0 (iωj 0 )Wii0 (ωj − ωj 0 ), (1)Wii0 (ωj − ωj 0 ),Wii0 (ωj − ωj 0 ),

(2)
(3)

Θi (iωj ) = [ωj Zi (iωj )]2 +[ξi + χi (iωj )]2 +[φi (iωj )]2 .
(4)
Indices i (i0 ) and j (j 0 ) denote the basis stateand Matsubara frequency, respectively. ξi isthe energy eigenvalue of the state i measured
from the chemical potential. Wii0 (νj ) is thematrix element of the interaction between the
Cooper pair states labeled by i and i∗ with the
Bosonic Matsubara frequency νj .
For homogeneous electron gas the basis icorresponds to the plane-wave state. We have
implemented a code package that solves theabove equations with an approximate manner. The interaction Wii0 (νj ) was set to thescreened Coulomb interaction Wel within the
random phase approximation plus the attraction mediated by the Einstein phonon with definite frequency and electron-phonon coupling
Wph . The diagonal terms Z and χ were evaluated at T = 0 noniteratively with initial value
Z = 1 and χ = 0 and without Wph , whichcorresponds to the G0 W0 approximation. The
equation for the nondiagonal term φ was laterlinearized with Wph and its largest eigenvalue,
which becomes unity at the transition temperature, was calculated. The iterative method

312

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

was used, which required only no more thana few hours in a serial run at system B at
T ∼ 100K.

with Z
with Z and χ

3.5

1 − ∆ (1/m∗ )Z[1 − ∆ (1/m∗ )Z ] [1 − ∆ (1/m∗ )χ ]

Λ0/Λ

3

References[1] G. M. Eliashberg, Zh. Eksp. Teor. Fiz. 38,(1960) 966.

2.5
2[2] A. B. Migdal, Zh. Eksp. Teor. Fiz. 34,(1958) 1438.

1.5
1using the HSE06 exchange-correlation functional [11] and spin-orbit coupling required
half a day per run with 2 nodes in fat groupin Ohtaka. For the memory requirement thejobs could not run in cpu queue group. Theband structure was interpolated using Wannier90 [12].

0

1

2

3
rs

4

5

6

[3] P. Morel and P. W. Anderson, Phys. Rev.
125, (1962) 1263.
[4] Y. Takada, J. Phys. Soc. Jpn. 45 (1978)
786.

Figure 1: Ratio of the largest eigenvalues of thelinearized Eliashberg equation with and without the self-energy corrections. Lines indicate
approximate formulas. Taken from Ref. [7]We show the results of the corrections to thepairing strength in Fig. 1. The blue open (purple closed) squares denote the inverse ratio of
the largest eigenvalues with Z (Z and χ) to
the non-interacting value. The values larger
than unity indicate suppression of the pairingby the Coulomb self energy. The χ term has
been known to cancel the effect of Z for theeffective mass. In contrast, our result showsthat χ cooperate with Z for the pairing suppression. This result highlight the importanceof χ for superconducting calculations, whichis computationally demanding and has mostly
been ignored in the previous studies.
We also analyzed the band structure of acompound HgCr2Se4 [8]. The calculations
were performed using Quantum Espresso [9]with norm conserving pseudopotentials from
Pseudo dojo [10]. Self consistent calculations[5] R. Akashi and R. Arita, Phys. Rev. Lett.
111 (2013) 057006.
[6] A. Davydov, A. Sanna, C. Pellegrini, J. K.
Dewhurst, S. Sharma, and E. K. U. Gross,
Phys. Rev. B 102, (2020) 214508.
[7] R. Akashi, Phys. Rev. B 105 (2022)
104510.
[8] H. Tanaka et al., Phys. Rev. Lett. accepted.
[9] P. Giannozzi et al., J. Phys. Condens.
Matter 29, (2017) 465901.
[10] G. Pizzi et al., J. Phys. Condens. Matter
32, (2020) 165902.
[11] J. Heyd, G. E. Scuseria, and M. Ernzerhof, J. Chem. Phys. 436 124, (2006)
219906.
[12] N. Marzari and D. Vanderbilt, Phys. Rev.
B 56 (1997) 12847.

313

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

Gapless symmetry-protected topological phase ofquantum antiferromagnets on anisotropic triangular
strip
Yuichiro HIDAKA1 , Shunsuke C. FURUYA2 , Atsushi UEDA1 , and Yasuhiro TADA3
1
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581
2
Department of Basic Science, University of Tokyo
Meguro, Tokyo 153-8902
3
Quantum Matter Program, Graduate School of Advanced Science and Engineering,
Hiroshima University, Higashihiroshima, Hiroshima 739-8530

Classifying phases of matter is one of theparamount issues in condensed matter physics.
One of the important classes of the quantum phase of matter is a symmetry-protected
topological (SPT) phase. An SPT phase is a
gapped phase that preserves some given symmetries but is not adiabatically connected toa trivial phase without breaking the symmetries [1]. Recently many researchers attempt
to extend the concept of the SPT phase into
gapless systems [2, 3].
One of the gapless SPT phases we study isthe anisotropic triangular strip (ATS) model,
which is the three-leg XXZ spin ladder withthe frustrated XXZ rung interaction. The lattice is shown in Fig.1. According to the ﬁeldtheoretic analysis, if the rung coupling is sufﬁciently small, the ground state of the ATSmodel is the gapless SPT state, in that the
ground state has simultaneously the properties of the Haldane state and the TomonagaLuttinger liquid. This property is also supported by the numerical calculation [4].
In the present study, we investigate theground state of the ATS model, where we expand the region of the parameters. We used
the ITensor and TeNPy libraries [5, 6]. We
found that the entanglement spectrum of theATS model is doubly degenerate, which is theimportant feature of the Haldane state if therung coupling J× is suﬃciently small (Fig.2).
Moreover, from the calculation of the entanglement entropy, we found that the ground state
has the central charge c = 1. On the other
hand, if the rung coupling J× is suﬃciently
larger than the ladder coupling J, the groundstate becomes the ferrimagnetic state. This ischecked by the calculation of magnetization.
Finally, if the rung coupling J× is of the sameorder as the ladder coupling, the ground stateis in the trivial phase, in which we cannot see
the even-fold degeneracy of the entanglement
spectrum, and is in a gapless phase. The phase
diagram of the ATS model is determined asshown in Fig.3 [7].

Figure 1: The anisotropic triangular strip(ATS) XXZ model. J denotes the ladder coupling and J× denotes the frustrated rung coupling. All coupling is the XXZ interaction:
Six Sjx + Siy Sjy + ∆Siz Sjz .

314

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[7] Y. Hidaka, S. C. Furuya, A. Ueda, and Y.
Tada, in preparation.

Figure 2: The entanglement spectrum of theATS model with ∆ = 0.8. When J× < 0.54 we
can see the even-fold degeneracy of the entanglement spectrum, which implies the Haldane
state.

Figure 3:model.

The phase diagram of the ATS

References[1] X. Chen, Z.-C. Gu, and X.-G. Wen, Phys.
Rev. B 83, 035107 (2011).
[2] T. Scaﬃdi, D. E. Parker, and R. Vasseur,
Phys. Rev. X 7, 041048 (2017).
[3] S. C. Furuya and M. Oshikawa, Phys.
Rev. Lett. 118, 021601(2017).
[4] Y. Hidaka, S. C. Furuya, A. Ueda, and Y.
Tada, Phys. Rev. B 106, 144436(2022).
[5] M. Fishman, S. R. White, and E. M.
Stoudenmire, arXiv:2007.14 (2020).
[6] J. Hauschild and F. Pollmann, SciPost
Phys. Lect. Notes , 5 (2018).

315

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

Finite temperature calculations of frustrated
quantum spin systems
Katsuhiro MORITA
Department of Physics and Astronomy, Faculty of Science and Technology,
Tokyo University of Science, Noda, Chiba 278-8510

The S = 1/2 kagome antiferromagnets are
known to have exotic properties even at ﬁnite temperatures. The speciﬁc heat has amultiple-peak structure [1], and the 1/3 magnetization plateau is expected to melt asymmetrically at ﬁnite temperatures [2, 3]. These
properties do not appear in the triangular lattice model, for example. Therefore, theseproperties are considered to be unique to the
kagome lattice model.
In this study, we have calculated the magnetic speciﬁc heat, entropy, magnetic susceptibility, and ﬁnite temperature magnetization
curves for the S = 1/2 J1 − J2 kagome lattice
Heisenberg model using the orthogonalized ﬁnite temperature Lanczos method [4, 5] on Supercomputer System B at the ISSP. J1 represents the nearest-neighbor interaction and J2
represents the next nearest-neighbor interaction. Here we set J1 = 1 as the energy unit.
In this study, the calculations were performed
with a 36-site kagome cluster.
We have found that the entropy remains
at low for −0.06 < J2 < 0.02, and outside
this range, the entropy decreases rapidly. At
J2 ≃ 0, the magnetic speciﬁc heat has a triplepeak structure, but as |J2 | increases, the speciﬁc heat tends to have a double-peak structure. In the ﬁnite temperature magnetization
curves, there is asymmetric melting of the 1/3
magnetization plateau around J2 = 0, but the
1/3 plateau becomes clearly visible as |J2 | increases. This is because the q = 0, up-up√√
down state stabilizes for J2 > 0, the 3 × 3

up-up-down state stabilizes for J2 < 0, and theplateau becomes unstable around the boundary.
These results will be compared with experimental results in the future to clarify the magnetic properties of the compounds.

References[1] T. Shimokawa and H. Kawamura, J. Phys.
Soc. Jpn. 85, 113702 (2016).
[2] J. Schnack, J. Schulenburg, and J.
Richter: Phys. Rev. B 98, 094423 (2018).
[3] T. Misawa, Y. Motoyama, and Y. Yamaji:
Phys. Rev. B 102, 094419 (2020).
[4] K. Morita and T. Tohyama: Phys. Rev.
Research 2, 013205 (2020).
[5] K. Morita: Phys. Rev. B 105, 064428(2022).

316

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

Machine learning molecular dynamics studies of complex
structural materials
Masahiko OKUMURA1, Hayato SHIBA2, Masatoshi HANAI3, Yasunobu ANDO4, Akiko
YAMAGUCHI1, Keita KOBAYASHI1, and Hiroki NAKAMURA1
1

Center for Computational Science and e-Systems, Japan Atomic Energy Agency,
Wakashiba, Kashiwa, Chiba 277-0871
2

Information Technology Center, University of Tokyo, Chiba 277-0882
3

Graduate School of Information Science and Technology,
University of Tokyo, Tokyo 133-8658

4

Research Center for Computational Design of Advanced Functional Materials,

National Institute of Advanced Industrial Science and Technology, Tsukuba, Ibaraki 3058568
dynamics

other methods. To overcome this weak point, the(MLMD) is a promising simulation method withmulti-stream EKF method is often used withhigh accuracy and low computational costs [1].

parallel computers [2]. The MD simulations

Artificial neural networks (ANNs) are used towith the ANN potentials for large systems are

reproduce

usually conducted on a supercomputer.

Machinelearning

theresults

molecularoffirst-principles

calculations (FPCs). The ANNs are trained using

Preparing the training dataset and training

datasets that consist of a lot of the results of

the ANNs are the time-consuming steps of

FPCs. For this method, supercomputers are used

MLMD. Notably, the ANNs for complex

for 1) making huge amounts of the FPC results,

materials require a large dataset, because many2) training the ANNs, and 3) conducting

atom configurations are realized, where themolecular dynamics (MD) simulations with thecomplex materials are, in this report, defined as

ANN potentials. Density functional calculations

those with large unit cells with periodic

are popular for making the training data. Several

boundary conditions. For making a dataset of

methods are used to train the ANNs, e.g.,complex materials, a large number of FPCs withstochastic gradient descent, adaptive momentum

many atom configurations in large unit cells are

estimation, extended Kalman filter (EKF), etc.

required. With the large dataset, heavy ANN

The EKF method is known as one of the best

training

optimizers for the ANN potential. However, thecalculations are popular to make the data, whichcalculation speed for an epoch is slower than the

are

317

is

popularneeded.

Densitycalculationsonfunctional

the

ISSP

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

training

Figure 1 shows the relation between thecalculations may not be popular. Therefore, we

number of cores and epoch/time (1/s) for Systemevaluated the scaling of the ANN training on

B and Oakbridge-cs. The blue and orange plots

System B, and the results are compared withshow the results of Oakbridge-cx, where thethose given by Oakbridge-cx operated by theformer and the latter are the results ofInformation Technology Center, The University

calculations using all and half of the cores on aof Tokyo. We prepared 7,007 data of ThO2,

CPU, respectively. We found that the reduction

which shows superionic transition below itsof the active cores improves the performance.

melting point [3]. Therefore, large unit cells areThis may be because the relative memory size

needed to reproduce the superionic transition

was increased by the reduction of the active

and melt, i.e., ThO2 can be considered a complex

cores. The gray and yellow plots represent thematerial although the structure of its solid state

results of System B, where the former and the

is

the

latter are the results of calculations using all andconfigurations of solid, super-ionic, and liquidhalf of the cores on a CPU, respectively. For

states of ThO2. The n2p2 open-source software

System B, the same tendency was obtained, i.e.,

for MLMD [4] was used, which provides thethe reduction of the active cores accelerates themultistream EKF method using MPI.

learning speed.

supercomputer.

simple.

But

The

the

ANN

datasetcontains

epoch/time (1/s)

We investigated the relation between the
0.006

number of cores and the number of epochs per

0.005

second. We found that the reduction of active

0.004

cores improves the calculation speed both for
Oakbridge_1

0.003

System B and Oakbridge-cx.

Oakbridge_2
systemB_1

0.002

systemB_2

References

0.001[1] J. Behler and M. Parrinello, Phys. Rev. Lett.

0
0

200

400

600

800

1000

1200

cores

98, 146401 (2007).

Fig. 1 The relation between the number of cores[2] A. Singraber et al., J. Chem. Theo. Comp. 15,

and the number of steps per second is shown.

3075 (2019).

The blue and orange plots represent full- and[3] K. Kobayashi et al., Sci. Rep. 12, 9808

half-core(2022).

calculations

on

Oakbridge-cs,

respectively. The gray and yellow pots represent[4] A. Singraber et al., J. Chem. Theo. Comp. 15,

full- and half-core calculations on Oakbridge-cs,

1827 (2019).

respectively.

318

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

Multi-impurity method for bond-weighted tensor
renormalization group
Satoshi MORITA
Graduate School of Science and Technology, Keio University
Kohoku-ku, Yokohama, Kanagawa 223-8522

In this study, we propose an algorithm
to calculate higher-order moments of physical quantities based on BWTRG. We introduce a bond-weight impurity which representsa physical quantity. With the appropriate update rules, coarse-grained bond-weights containing multiple impurities can be generated
to calculate the higher-order moments of thephysical quantity. Our proposed method iscompared with conventional methods on thetwo-dimensional classical spin model. The
proposed method achieves higher accuracy atlower computational cost than the higher-order
TRG algorithm [4]. We also show that thefinite-size scaling analysis of the squared magnetization provides critical exponents and distinguishes the weakly first-order and continuous phase transitions.

0.5

TRG
BWTRG
Exact

0.4
0.3
m2

Tensor network (TN) methods are attracting
much attention as powerful tools for computing strongly correlated many-body problems.
The partition function of classical statisticalsystems can be represented by the TN form.
However, the contraction of a large TN stillrequires an exponentially large computational
effort. The concept of the real-space renormalization group resolves this problem. The tensor renormalization group (TRG) method andits variants calculate a coarse-grained tensor
by information compression using the singularvalue decomposition [1, 2]. These methods cancalculate the partition function approximately
in logarithmic time for the system size. Recently, Adachi, et al. have proposed the bondweighted TRG (BWTRG) method, which improves the accuracy of TRG by introducing
a bond weight and distributing it appropriately [3].

0.2
0.1
0.0
2.24

2.25

2.26

2.27
T

2.28

2.29

2.30

Figure 1: The squared magnetization of theIsing model on the square lattice calculated by
TRG and BWTRG with the bond dimension
χ = 128 and the system size L = 225 . The
green solid curve shows the exact solution in
the thermodynamic limit.

References[1] M. Levin and C. P. Nave, Phys. Rev. Lett. 99,
120601 (2007).
[2] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L.
P. Yang, and T. Xiang, Phys. Rev. B 86, 1(2012).
[3] D. Adachi, T. Okubo, and S. Todo, Phys. Rev.
B 105, L060402 (2022).
[4] S. Morita and N. Kawashima, Comput. Phys.
Commun. 236, 65 (2019).

319

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

Majorana Kramers Qubits and Yang-Lee Anyons in
Topological Superconductors
Takeshi MIZUSHIMA
Department of Materials Engineering Science, Osaka University
Toyonaka, Osaka 560-8531

We have investigated non-Abelian statisticsof Majorana Kramers pairs (MKPs) in timereversal invariant topological superconductors(TRITSCs) [1].
Numerically solving thetime-dependent Bogoliubov-de Gennes equation, we have simulated the braiding dynamics of MKPs in a network of one-dimensional
TRITSCs. We have first considered the tolerance of MKPs against a perturbation with broken time-reversal symmetry such as a magnetic
field and demonstrated the robustness of the
non-Abelian braiding of MKPs when the initial and final states of a braiding process maintain the combined symmetry of time reversal
and mirror reflection symmetries. The Majorana Kramers qubit is tolerant even when intermediate states break the combined symmetry. We have also discussed the eﬀect of gateinduced inhomogeneous potentials at junctions
between superconducting nanowires. Such a
potential generally induces a non-Majorana
nearly zero-energy Andreev bound state at thejunctions. Our numerical simulation showsthat the dynamics of MKPs accumulate significant errors when such non-Majorana states exist and interfere with MKPs. Although we concentrated on Majorana particles in a TRITSC
in this work, it is expected that eﬀects of nonMajorana states generated by inhomogeneousgate potentials on the non-Abelian braiding
are ubiquitous for any topological superconductors with broken time-reversal symmetry.
We have also demonstrated that a Yang-Leeanyon system is constructed from collective

Majorana bound states in a network of topological superconducting nanowires [2]. YangLee anyons are nonunitary counterparts of Fibonacci anyons, obeying the same fusion rule.
Because of nonunitarity, the central charge andthe scaling dimension for the one nontrivial
primary field are negative, c = −22/5 and∆ = −2/5, respectively. The nonunitary conformal field theory with c = −22/5 describes
the nonunitary critical phenomenon known as
the Yang-Lee edge singularity. In Ref. [2],we consider a topological superconductor junction system coupled with dissipative electron
baths as a non-hermitian platform for interacting Majorana quasiparticles. Numericallycalculating the central charge and scaling dimension, we find the condition that the nonHermitian Majorana system can simulate theIsing spin model of the Yang-Lee edge singularity and confirmed that, by controlling model
parameters in a feasible way, the Yang-Leeedge criticality is realized. We also discuss thescheme for the measurement and the braidingof Yang-Lee anyons, aiming at the application
to universal quantum computation.

References[1] Y. Tanaka, T. Sanno, T. Mizushima, andS. Fujimoto: Phys. Rev. B 106, 014522(2022).
[2] T. Sanno, M. G. Yamada, T. Mizushima,and S. Fujimoto: Phys. Rev. B 106,174517 (2022).

320

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

Quantum Monte Carlo simulations of
Spin Hall magnetoresistance in quasi-two-dimensional
antiferromagnetic insulator/metal bilayer systems
Takeo Kato
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Spin Hall magnetoresistance (SMR) is a
novel type of magnetoresistance caused by spincurrent. SMR was firstly observed in 2013 by
Nakayama in a bilayer system composed of ametal and a magnetic insulator [1]. The
phenomenological theory for SMR [2] is widelyused for analysis of SMR, but it can’t describe
the temperature dependence of SMR in a bilayer
systemcomposed

of

a

metal

and

Fig. 1:Setup for SMR calculation cited from[5].

anWe used DSQSS for the quantum Monteantiferromagnetic-insulator [3] .

In this study, we discuss the SMR

Carlo method, and we adjustment some code

dependence of temperature, thickness of AFI,

to be able to calculate the local spin

and randomness of exchange interaction. To

susceptibility. The accuracy for the calculation

discuss about these features, we expand theis proportional to the number of Monte Carlomicroscopic theory of SMR based on Keldysh

steps,

Greens function [4]. In this theory the spin

calculations with MPI, which is included incurrent is expressed in terms of the spin

DSQSS, to achieve a high degree of accuracy.

susceptibilityof

antiferromagnetic

We show the result of this calculations

insulator, but analytical calculation based on

based on the paper [5]. The temperature andspin-wave theory are available for only low

thickness dependence of SMR is shown in Fig.

temperature.

widely2. The black points indicate the total SMR

temperature dependence of SMR by thesignal, whereas the red and blue plots arequantum Monte Carlo method, we obtained

contribution from z-component and xy-

thecomponents of spin susceptibility. Green

Due

the

so We performed parallelizationtocalculate

spin conductance formulation by

Matsubara Greens function [5].

arrows indicate the ordering temperature.

321

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

The effect of randomness by disorderedexchange interactions is shown in Fig. 3. It
shows that the ordering temperature issuppressed and the SMR signal becomes small
in the disordered case.

We calculated SMR in AFI by quantum Monte
Carlo simulations, and we discussed about the
SMR dependence of temperature, thickness ofFig. 2: S=1 result for SMR cited from [5].

AFI, and randomness of exchange interaction.

It shows that at high temperature SMR
signals vanish in consistent with the analytical
calculation by the high temperature expansion.
SMR has a peak less than the ordering

References

temperature and the thickness dependence[1] H. Nakayama, et al., Phys. Rev. Lett.

becomes weaker for four or more layers.

110, 206601 (2013).
[2] Y. T. Chen, et al., Phys. Rev. B 87,
144411 (2013).
[3] D. Hou, et al., Phys. Rev. Lett. 118,
147202 (2017).
[4] T. Kato, Y. Ohnuma, and M. Matsuo,
Phys. Rev. B 102, 094437 (2020).
[5] T. Ishikawa, M. Matsuo, and T. Kato, Phys.
Rev. B 107, 054426 (2023).

Fig. 3: Disordered case of result cited from [5].

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

Precise structural analysis of atomic arrangement ofhydrogen-adsorbed ceria surfaces by using 2DMAT
Izumi MOCHIZUKI
Slow Positron Facility, Institute of Materials Structure Science, High Energy Accelerator
Research Organization (KEK), Oho 1-1, Tsukuba, Ibaraki, 305-0801, Japan
We used 2DMAT [1-7], an open-source data-as an innovative feature of TRHEPD. Moreanalysis software that is preinstalled on the ISSP

detailed analysis is ongoing.

supercomputers, developed by the PASUMS

References

project at FY2020 and 2021 for data analyses of[1]

total-reflection high-energy positron diffraction

https://www.pasums.issp.utokyo.ac.jp/2DMAT(TRHEPD). TRHEPD is an innovative method[2] K. Tanaka, et al. (I. Mochizuki), Acta. Phys.

forPol. A 137, 188 (2020).

theprecise determination of atomic

arrangement of crystal surfaces due to its[3] T. Hoshi, et al. (I. Mochizuki), Comp. Phys.

exceeding surface sensitivity. It is realized atCommun. 271, 108186 (2022).

Slow Positron Facility, IMSS, KEK [8].

[4] T. Hanada, et al., Comp. Phys. Commun. 277

We applied 2DMAT to the structural analysis

108371 (2022).

of a Ge(001)-c(4×2) surface as a test of the[5] Y. Motoyama, et al. (I. Mochizuki), Comp.

analysis method by an optimization procedure,Phys. Commun. 280, 108465 (2022).

Nelder-Method method [5, 7].

[6] K. Tanaka, et al. (I. Mochizuki), JJAP Conf.

We also investigated a hydrogen (H)Series 9, 011301 (2023).

adsorbed CeO2(111)-(1×1) surface [9]. It was a[7] N. Kinoshita, et al. (I. Mochizuki), in JPS

challenge for the detection of surface hydrogenmeeting, 22-25, Mar. 2023.

since it is quite difficult for other experimental[8] https://www2.kek.jp/imss/spf/eng/

measurements. In this study, we first calculated[9] I. Mochizuki, in Innovation catalysis science

the structure of the bare CeO2(111)-(1×1)

through the integration of high-performancesurface and that of the H-adsorbed one by using

data-driven science and advanced measurement

Quantum ESPRESSO in the DFT+U method.

technology,

The results were then compared with those

Usage/Research Center for Catalysis (JURCC),

obtained by the TRHEPD experiment. We found30.

that the TRHEPD can detect the surface

https://sites.google.com/view/hpc-catal2023

hydrogen properly with the data-driven analysis,323

onlinemini-symposium,

Mar.

Joint

2023,

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

Equilibration and glass transition in the dense hard
sphere systems
Masaharu ISOBE
Nagoya Institute of Technology
Gokiso-cho, Showa-ku, Nagoya, 466-8555

As the simplest models, the hard
disk/sphere systems have been investigated
via molecular simulation in both equilibriumand non-equilibrium phenomena [1].
Inthis project, we focus on the equilibrationand the non-equilibrium phase transitionin the hard disk/sphere model system withmodern algorithms, i.e., Event-Chain Monte
Carlo (ECMC) [2, 3], Newtonian ECMC
(NEC) [4] and Event-Driven Molecular
Dynamics(EDMD) [5].
Eﬃciency of Newtonian Event-Chain Monte
Carlo in dense hard disk systems
In the statistical mechanics of many bodysystems, the physical properties at the true
equilibrium state are crucial. To invent efficient method of equilibration in the molecular simulation is one of fascinating topics
where the simple model systems such as hardspheres and Ising spins have been investigatedactively and made important contribution. In
2009, Event-Chain Monte Carlo (ECMC) [2]
was invented in hard disk systems that updates the particle position with sequential collisions (i.e., event-chain) and has been applied
to the many typical models studied in statistical physics [3]. ECMC has several variants depending on how the direction of movement of the particles after the collision is determined. Recently, Newtonian-ECMC (NEC) [4]
was invented by applying the velocity collision rule to each collision event to dictate thedirection of particle movement; this is often
used in the Event-Driven Molecular Dynamics(EDMD) of hard-sphere systems. NEC outperformed not only other variants of ECMC but
EDMD in the diﬀusion coeﬃcient and melting
process in hard sphere systems. However, theoptimal performance strongly depends on thephysical properties and chain length (or duration) of ECMC (and its variants). Previously,
we found the high eﬃciency of NEC in diﬀusion coeﬃcient of liquid state in hard disk systems for all chain lengths and various systemsizes, which is caused by the spatially homogeneous sampling of the displacement (i.e., collision) [7]. To figure out the mechanism of equilibration of NEC in dense systems, we focus on
the equilibrium state of dense hard disk systems above the fluid-to-solid transition point,
which show the relaxation via hopping motion
even in the packing fraction at the solid phasein hard disk systems for each methods. We also
found profoundly eﬃcient diﬀusion coeﬃcientby NEC even when the relaxation becomes
only successive hopping motion in highly dense
systems.
Free volume estimator (NELF-A) in dense
poly-disperse hard disk systems — local structure analysis on pressure and inherent structure
In many-particle systems, the free volume
of a tagged particle constructed from excluded
volume by surrounding particles made a crucial contribution to describing the macroscopic
properties in the history of liquid state theory.
Several numerical algorithms for calculating
the free volume (FV) have been invented, such

324

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

as Monte Carlo (MC) sampling and Voronoi
tessellation. However, MC is an approximate
calculation, and the Voronoi tessellation requires high computational costs (especially forpoly-disperse systems). As an alternative algorithm, we propose the simple, eﬃcient, andprecise method, Neighbors for Enclosing Local
Free Area (called NELF-A), which is easily applied to dense poly-disperse hard disk systems
often used in glassy model systems [8].
As an application of the NELF-A, we implement an alternative method to obtain the inherent structure instead of short-time averaged
coarse-graining trajectories. The summary of
the algorithm is as follows: (i) The displacement vectors from a particle position to a geometric center of FV estimated by NELF-A foreach particle are drawn. (ii) The positions in
a whole system are updated with the displacement vectors multiplied by a scaling parameterβ less than unity. (iii) Repeat (ii) until the positions converge toward an inherent structure.
The FV in equilibrium state eventually converging positions after (i)-(iii) procedure andthe isotropic shape of FV located on a triangular lattice could be observed (i.e., inherent
structure). Even in a typical non-equimolar
binary mixture hard disk system, the probability density distribution obtained by (i)-(iii)
procedure is fairly consistent with the conventional method based on the short-time averaged coarse-graining trajectories. We also confirmed the accuracy of this method by estimating pressure based on FV-based as comparedwith the other conventional methods.

tional order (so-called Alder transition) due tothe exclusive eﬀect of hard disk, which causes
anomalous fluctuations around phase transition and transition shifts [9].

References[1] M. Isobe: Mol. Sim. 42 (2016) 1317.
[2] E. P. Bernard, W. Krauth, and D. B. Wilson: Phys. Rev. E 80 (2009) 056704.
[3] W. Krauth: Front. Phys. 9 (2021) 229.
[4] M. Klement and M. Engel: J. Chem.
Phys. 150 (2019) 174108.
[5] M. Isobe: Int. J. Mod. Phys. C 10 (1999)
1281.
[6] M. Klement, S. Lee, J. A. Anderson and
M. Engel: J. Chem. Theor. and Comp. 17(2021) 4686.
[7] H. Banno, D. Mugita, and M. Isobe: J. of
Phys.: Conf. Ser. 2207 (2022) 012011.
[8] K. Souno, H. Koyama, D. Mugita, and M.
Isobe: in preparation.
[9] N. Murase and M. Isobe: Proceedingsof the 28th Symposium on Traﬃc Flowand Self-driven Particles, 28 (2023) 9. (inJapanese)

Anomalous phase transition in self-propelled
hard disk systems
In this study, we investigate the phase behavior of the self-propelled hard disk systems
with the Vecsek-type interaction via EDMD
simulation systematically. In addition to theordinal order-disorder transition of the collective velocity field known in the original point
particle of the Vicsek model, we observed thenovel competition driven by the global posi-

325

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

Molecular dynamics analyses of ion migration at electrodeelectrolyte interfaces and grain boundaries in electrolytes
Ryo KOBAYASHI
Department of Applied Physics, Nagoya Institute of Technology,
Gokiso-cho, Showa-ku, Nagoya, Aichi 466-8555
aresimulation and local ion-flux analyses, weconsidered as a promising candidate for next-

showed that, in poly-crystalline materials, thegeneration lithium-ion batteries. They use aions migrate towards the lower reaches of grainsolid-state electrolyte instead of the liquid-stateand go through spots at GBs, and the high-flux

electrolyte used in conventional lithium-ion

spots contain migration paths as low energy asbatteries.

the bulk ones.

All-solid-stateThebatteries

key

to(ASSBs)achievinghighperformance in ASSBs is to reduce or control

This methodology can be applicable to a

the resistivity at interfaces such as thewide variety of systems for electrolytes andelectrode/electrolyte

interfaces between electrode and electrolyte.

interface

and

grain

boundaries (GBs) within the electrolyte. To
achieve this, we need a deeper understanding of
the atomic structure and precise mechanisms ofion migration at GBs. The goal of this project is

References

to obtain new insights into these processes by[1] R. Kobayashi, Y. Miyaji, Nakano, K. M.

using large-scale molecular dynamics (MD)Nakayama, APL Materials 8, 081111 (2020)

simulations.

[2] Nagoya Atomistic-simulation Package (nap)

In order to obtain atomistic insights on Liion

migrationmechanism

at

GBs,

http://github.com/ryokbys/nap

we[3] Optzer, http://github.com/ryokbys/optzer

performed non-equilibrium MD simulation of[4] R. Kobayashi, K. Nakano, M. Nakayama,

large-scale poly-crystalline Li2Zr(PO4)3 (LZP)

Acta Materialia, 226 (2022) 117596.

system containing half a million atoms using a[5] M. Nakayama, K. Nakano, M. Harada, N.

interatomic potential generated using a meta-

Tanibata, H. Takeda, Y. Noda, R. Kobayashi, M.

heuristic approach [1,2,3] and extracted local

Karasuyama, I. Takeuchi, M. Kotobuki, Chem.

ion flux information around GBs [4]. From theCommun.

326

58,

9328–9340(2022)

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

Theoretical study for firefly bioluminescence
substrate analogs and related molecules
Miyabi HIYAMA
Graduate School of Science and Technology, Gunma University
Tenjin-cho, Kiryu, Gunma 376-8515

In this year, we studied the stable structuresof adenosine triphosphate (ATP)-metal complexes in aqueous solutions and the theoreticalabsorption spectra for the firefly luciferin analog named ”seMpai”.
It is believed that the cofactor ATP and
Mg2+ form a complex and react with thesubstrate luciferin in firefly bioluminescence.
Wang et al. [1] reported that the bioluminescence quantum yields decreased with increasing cadmium or zinc concentration whenthese metal ions were used instead of magnesium one. The effects of metals on luciferin,
on ATP, and on proteins may be responsible
for this feature. We expected that the effect of metals on ATP was most significant
than the others because the ATP is directly
related to the bioluminescence reaction. Thus,
the optimized structures of ATP-Mg2+ , ATPZn2+ , and ATP-Cd2+ in aqueous solutions
were obtained using the density functional theory (DFT). For the initial structures for thesecalculations, the results from Global Reaction
Mapping (GRRM11) calculations [2] were usedbecause there were many local minimum structures of these complexes. We found that there
was no clear difference in the structure of these
metal ion coordination to ATP. It also foundthat ATP-Zn2+ and ATP-Cd2+ were unstable
in aqueous solutions, unlike ATP-Mg2+ .
In these days, one of the firefly luciferin
analogs named ”seMpai”, which is preparedwith the aim of improving water solubility
and produces red luminescence with firefly lu-

ciferase, was reported [3]. To elucidate theabsorption characteristics of this analog, theoptimization structure of its ground state forits conjugate acids and bases expected to be
main component in the aqueous solutions wereobtained from DFT calculations. The time dependent DFT (TDDFT) calculations were carried out to estimate the theoretical absorption
spectra for these structures.
The computational chemistry software
package, Gaussian09 [4], was used for DFT
and TDDFT calculations. All calculations using Gaussian09 and GRRM11 were performedon system B of Super Computer Center in
ISSP.

References[1] Y. Wang, H. Kubota, N. Yamada, T. Irie,
H. Akiyama: Photochem. Photobiol. 87,
846 (2011).
[2] K. Ohno, S. Maeda: Chem. Phys. Lett.
384, 277 (2023). S. Maeda, K. Ohno,:
J. Phys. Chem. A, 109 5742 (2005). K.
Ohno, S. Maeda, J. Phys. Chem. A, 1108933 (2006).
[3] R. Saito, T. Kuchimaru, S. Higashi, S. W.
Lu,M. Kiyama, S. Iwano, R. Obata, T.
Hirano, S. Kizaka-Kondoh, S. A. Maki:
Bull. Chem. Soc. Jpn. 92, 608 (2019).
[4] Gaussian09 Revision D.01, M. J. Frisch et
al.

327

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

Analysis of Ising model in statistical-mechanical informatics
Yuya SEKI
Graduate School of Science and Technology, Keio University, Kanagawa 223-8522
We have developed an optimization method

use multithreading to reduce computational

for discrete black-box optimization problems

time for the Monte Carlo simulation to

by combining machine learning, simulated

minimizes the model equation of the trained

annealing, and integer encoding methods [1].

FM. By using the parallel computing on the

Our approach is based on factorization machinesupercomputer, we obtained numerical results

withefficiently.

annealing

(FMA)

[2]

that

uses

This year we investigated the performancefactorization machines as machine learningmodel. Factorization machine with annealing

of

FMA

approximates an objective function of a black-

optimization problems in detail [1]. Integer

box optimization problem by training avariables can be represented as binary variables

factorization machine. Instead of minimizing

through integer encoding methods such as

the objective function that requires high cost to

binaryevaluate it, FMA minimizes a model equation

domain-wall encoding. Since the encoded bitof FM. Here, simulated annealing can be usedstrings differ depending on the encoding

for the minimization because the model

methods,

equation of FM has a QUBO form that isdetermined by the choice of the encoding

equivalent to Ising model. Since FM is anmethod. Through a numerical simulation, we

important factor in performance of FMA, wefound that the FMA with one-hot encoding

analyze properties of model equations of FM,

method has higher performance than FMA withwhich are Ising models, by using numericalother encoding methods. By using simulation

simulations on the supercomputer.

on

the

for

integer-variable

encoding,

the

one-hot

performancesupercomputer,

weblack-box

encoding,ofFMA

revealed

required to achive a high performance. To thisthe performance of FMA.

end, we performed the calculation in parallel by

Referenceschanging[1] Y. Seki et al.: arXiv:2209.01016 (2022).

parameters

with

is

the

relationship between sparsity of bit string and

Tuning parameters in the FMA method is

the

and

multiple

processes on the supercomputer. In addition to[2] K. Kitai et al.: Phys. Rev. Research 2

the parallel computing mentioned above, we(2020) 013319.

328

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

Quantum Pyrochrole Magnet
H. Kadowaki
Department of Physics, Tokyo Metropolitan University
Hachioji, Tokyo 192-0397, JapanWe are investigating static correlaz σ z ⟩ and dynamic correlations
tions ⟨σQ−Q
z (0)σ z (t)⟩ characterized by the modu⟨σQ−Qlation wave vector k ∼ ( 12 , 12 , 12 ), which wereobserved by neutron scattering experiments,
e.g. Ref. [1, 2], in the frustrated pyrochlore
magnet Tb2+x Ti2−x O7+y (TTO). TTO withx < xc ≃ −0.0025 is very interesting becauseof their quantum spin liquid (QSL) groundstate. A simple question “why TTO shows thestatic short-range order with k ∼ ( 12 , 21 , 12 )?”has been a diﬃcult problem for many years[3]. Recently we have solved this problem
using quantum simulation based on thermal
pure quantum (TPQ) states, and proposedthat small three-spin interaction terms of aform σr± σrz ′ σrz ′′ are a key ingredient of bothk ∼ ( 12 , 12 , 12 ) and the QSL state of TTO at
low temperatures [4].
We are now extending the quantum simuz σ z ⟩ to inelastic neutron speclation for ⟨σQ−Qz (0)σ z (t)⟩ using a
tra, i.e., simulating ⟨σQ−Q
TPQ method [5]. We have written a simulation program, which now works very well,and are performing time-consuming calculations, which we would like to continue in 2023.

References[1] H. Takatsu et al. Phys. Rev. Lett. 116,
(2016) 217201.
[2] H. Kadowaki et al. J. Phys. Soc. Jpn. 87(2018) 064704.

329[3] H. Kadowaki et al. Phys. Rev. B 99 (2019)014406.
[4] H. Kadowaki, M. Wakita, B. Fåk, J.
Ollivier, and S. Ohira-Kawamura, Phys.
Rev. B 105 (2022) 014439.
[5] H. Endo et al. Phys. Rev. Lett. 121 (2018)
220601.

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

Atomistic model study on the coercivity mechanism of
permanent magnets
Hiroshi Hayasaka1, Masamichi Nishino2, 1, and Seiji Miyashita3,4,1
1

Elements Strategy Initiative Center for Magnetic Materials, National Institute for Materials
Science, 1-2-1 Sengen, Tsukuba, Ibaraki 305-0047, Japan
2

Research Center for Advanced Measurement and Characterization,

National Institute for Materials Science, 1-2-1 Sengen, Tsukuba, Ibaraki 305-0047, Japan
3

Institute for Solid State Physics, The University of Tokyo,
5-1-5 Kashiwanoha, Kashiwa 277-8581, Japan

4

The Physical Society of Japan, 2-31-22 Yushima, Tokyo 113-0033, Japan

1. Angular dependence of coercivity

2. Coercivity enhancement by dysprosium
substitution into neodymium permanent

We investigated angular dependences ofnucleation and pinning fields in anisotropic

magnets

Heisenberg model. Unlike a deterministic

We investigated the features and mechanisms ofprocess in magnetization reversal at zero

the enhancement of the threshold field fortemperature, reversal is a stochastic process formagnetization reversal in atomistic Nd-Fe-B

barrier crossing at finite temperatures. We

model with Dy substitution. The results showedstudied the thermal fluctuation effect using thethat a change from surface to bulk nucleation

stochastic LLG equation method [1].

occurs when the number of substituted layers

At

zero

temperature,thenormalized

increases and the anisotropy energy of Dy is

nucleation field shows an asymmetric bowl-likeresistantangular dependence, and the value at large

significantly enhances coercivity, especially atangles rapidly increases. On the other hand, thehigh temperatures [3].

pining field shows almost constant dependence

Referencesat θ < 40◦ and increases at larger θ. The thermal[1] M. Nishino, and S. Miyashita, Phys. Rev. B

fluctuation effect reduces the nucleation and91, 134411 (2015)

pinning fields. At higher temperatures, the[2] H. Hayasaka, M. Nishino, and S. Miyashita,

thermal effect makes both fields closer, in which

Phys. Rev. B 105, 224414 (2022)

the surface nucleation at the hard magnet plays[3] M. Nishino, H. Hayasaka, and S. Miyashita,

an important role [2].

Phys. Rev. B 106, 054422 (2022)

330

to

temperature

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

Large-scale molecular-dynamics simulation of pressure-induced
coordination change in silica glass with ANN potentials
Daisuke WAKABAYASHI
Photon Factory, Institute of Materials Structure Science,High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801
Silica is well known as an archetypal oxide

The test calculations were performed for theand the compression behavior of silica glass andsolidification process of silica melt under highmelt has attracted considerable attention inpressure. The ANN potential was determinedvarious research fields.

Theoretical studies

with the total potential energy, atomic force, andusing molecular-dynamics simulations havepressure data obtained by ab-initio calculations

been actively conducted to obtain information

for a temperature range of 300-5,000 K at 60

on the atomic arrangements and bonding states

GPa with 144 atoms. Using the potential, weof silica glass under high pressure. Since there

made molecular-dynamics calculations on silica

is no long-range order in amorphous materials as

melt at 60 GPa and 4,000 K for a system of about

exists in crystals, calculations on a large number

30,000 atoms.

of atoms are required to reproduce the structural

The use of multiple interatomic potentials

changes. There is a limit to the size of theaveraged together stabilizes the calculations.

system which can be calculated by ab-initio

The initial parameter sets for the fitting were

methods.

randomly generated, and five ANN potentials

In this study, we have fitted the potential with

were determined from the five parameter sets,a machine-learning technique (ANN potential)

respectively. In order to avoid extreme results,

on the basis of the results of ab-initio

molecular-dynamics

calculations, and have applied it for the

performed in an extended system using the

calculations with an extended system. However,average value of the atomic forces calculated byin our previous projects, it has been observed

the five ANN potentials.

simulations

were

that the large-scale calculations of melt wereIn addition, active learning was carried out todifficult to stabilize and quickly diverged,

efficiently compensate for the lack of data. First,

making them impossible to continue. In thismolecular-dynamics

project, potential averaging and active learning

performed at 60 GPa and 4,000 K with theareaveraged ANN potentials in a system of the same

employedincombination

to

avoid

computational instability efficiently [1].

simulations

weresize as in the ab initio calculations. The atomic

331

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

forces calculated from each ANN potential

calculations, but have also been found toshould agree, and therefore it is likely thatsignificantly improve the precision of theatomic configurations at steps having a large

intermediate-range structure and reproduce thevariation in atomic force among the fiveexperimental data in the literature [3]. Thispotentials have not been well learned. Then, thesuggests that the experimentally observed

standard deviation of the atomic forces was

structural transformation of silica glass, in which

estimated for all atoms at each step, and thethe coordination number of Si changes from 4 toatomic configurations for which the maximum

6 with increasing pressure [4], can be simulated

value was larger than 0.25 eV/Å were extracted

in more detail. In fact, the calculations of the[2]. The total potential energies, atomic forces,

depressurization process down to 0 GPa wereand pressure in the extracted configurations

successful with the ANN potentials determined

were calculated by an ab initio method andin the same method as for the cooling process

added to the initial dataset. Active learning was

and reproduced the change in coordination

carried out twice in total.

number from 6 to 4.

The averaged ANN potential obtained after
two active-learning cycles did not cause any

Referencescomputational instability during the simulations[1] D. Wakabayashi, K. Shimamura, A. Koura,at 60 GPa and 4,000 K and in the subsequentand F. Shimojo, submitted to J. Phys. Soc. Jpn.

simulations of the cooling process down to 300[2] Y. Zhang, H. Wang, W. Chen, J. Zeng, L.

K.

Zhang, H. Wang, and W. E, Comput. Phys.

Without active learning, the calculation

becomes unstable before 2,000 steps and cannot

Comm. 253, 107206 (2020).

reproduce the results of the ab initio calculations.

[3] D. Wakabayashi, N. Funamori, T. Sato, and

The time at which the destabilization occurs is

T. Taniguchi, Phys. Rev. B 84, 144103 (2011).

significantly delayed with each active-learning[4] T. Sato and N. Funamori, Phys. Rev. B 82,

cycle.

184102 (2010).

The calculations on silica glass at 60 GPa and300 K have not only reproduced ab initio

332

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

First-principles investigation of ductile deformation
mechanism in silver sulfide
Masaaki MISAWA
Department of Intelligent Mechanical Engineering, Fukuoka Institute of Technology
Wajiro-higashi, Higashi-ku, Fukuoka 811-0295

β-Ag2 S, a low-temperature phase of silver
sulfide, exhibits metal-like ductility at room
temperature even though it is an inorganic
semiconductor[1]. Elucidating the origin of
this ductility will be extremely useful for designing materials for flexible next-generation
semiconductor devices. Although several theoretical investigations have been done to elucidate the mechanism of ductile deformation of
β-Ag2 S from the viewpoint of bonding statesin crystals [1, 2], theoretical verification considering dynamics under finite temperatures has
not yet been performed. Therefore we have investigated the behavior of β-Ag2 S under shear
stress based on first-principles molecular dynamics method [3].
In this study, we have performed firstprinciples molecular dynamics simulations forthe six types of deformations of β-Ag2 S crystals: the (100)[010], (100)[001], (010)[100],(010)[001], (001)[100], and (001)[010] shear,
where (KLM )[klm] shear means that simple shear deformation in which the crystal’s(KLM ) plane slides toward the [klm] direction. The shear rate in the simulations was
equivalent to dγ/dt = 0.05 ps−1 , and shear
deformations up to γ = 0.50 were performed
under room temperature (300 K). The radial
distribution function of the sulfur ion sublattice during the deformation process wasanalyzed to elucidate atomistic behaviors of
β-Ag2 S under shear stress, and it was observed that peaks that once disappeared during deformation were restored in the four sys-

tems: (100)[001], (010)[100], (001)[100], and(001)[010] shear. In addition, we performeda more detailed study for the (100)[010] and(010)[001] shear systems, in which structural
recovery did not occur, by slowing the shear
rate, and observed structural recovery in thesesystems as well. A common feature of thesestructural recovery mechanisms is that thefracture is repaired by the mass movement ofsulfur atoms in a particular crystal plane.
We also pointed out that the BCC-like structure of the sulfur sublattice in β-Ag2 S may be
responsible for this atomistic behavior. Ag2 Se,
which is not ductile, does not have such ahighly symmetric anion-sublattice and thus
cannot undergo the mass movement of atoms
in response to shear stress, resulting in difficulty in structural recovery. Such atomistic
understanding is expected to accelerate research and open a possibility in designing awide variety of inorganic ductile semiconductors for future flexible electronics based on silver chalcogenides.

References[1] X. Shi et al.: Nat. Mater. 17, 421 (2018).
[2] G. Li et al.: npj Comput. Mater. 4, 44(2018).
[3] M. Misawa et al.: Sci. Rep. 12, 19458(2022).

333

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

Development of machine-learning-based phase
diagram construction method for high-throughput
batch experiments
Ryo TAMURAInternational Center for Materials Nanoarchitectonics,
National Institute for Materials Science,
1-1 Namiki, Tsukuba, Ibaraki, 305-0044

Phase diagrams are valuable guidelines forthe development of new materials. Throughexperiments, phase diagrams have been determined in various spaces, such as between processes, between external ﬁelds, and between
compositions. In the ﬁeld of materials informatics, tools focusing on phase diagrams
have been developed using machine learning(ML). Among these tools, we have developed
an ML method based on uncertainty sampling(US)[1] to eﬃciently construct phase diagrams,
which we call the Phase Diagram Construction(PDC) package[2]. This method evaluates theuncertainty based on the ML model learned
under the conditions that the phase domain
is speciﬁed. Using evaluated uncertainty, theexperimental condition with most uncertainty
in the phase diagram is selected as the next
experimental condition. Based on this suggestion, an experiment will be conducted and thecorresponding phase domain is identiﬁed. The
number of known points is then increased andthe next experimental condition is proposedusing the updated ML model. By applying
the PDC to the known phase diagram, it wasdemonstrated that this iterative process canrapidly construct an accurate phase diagramand discover new phase domains[1]. For further validation, PDC was applied to construct
a new diagram for obtaining Zn-Sn-P ﬁlms bymolecular beam epitaxy[3]. The PDC experiments detected new phases that were not ob-

served in the initial stage, and detailed phaseboundaries were quickly determined. Thus,
the PDC proved to be useful in the case where
one condition is proposed in ML and one experiment is performed in each iteration. Further speed-up can be achieved by incorporating
the knowledges of materials science[4].
In addition, we developed PDC for highthroughput batch experiments[6]. To address
the performance of the method, the target wasa Cu-Mg-Zn ternary alloy. Both isothermal
sections of this system (two-dimensional) andits temperature-dependent phase diagram asa whole (three-dimensional) were considered.
We showed that PDC is suitable for highthroughput batch experiments to eﬃciently
construct phase diagrams. Using this version,a temperature-composition phase diagram ofa crosslinked polymer was constructed, and its
usefulness was conﬁrmed[5].

References[1] K. Terayama, R. Tamura, Y. Nose, H. Hiramatsu, H. Hosono, Y. Okuno, and
K. Tsuda, Physical Review Materials 3,033802 (2019).
[2] https://github.com/tsudalab/PDC[3] R. Katsube, K. Terayama, R. Tamura,and Y. Nose, ACS Materials Letters 2,
571-575 (2020).

334

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[4] K. Terayama, K. Han, R. Katsube,
I. Ohnuma, T. Abe, Y. Nose, and
R. Tamura, Scripta Materialia 208,
114335 (2022).
[5] W. Hu, T. Chen, R. Tamura, K. Terayama, S. Wang, I. Watanabe, and
M. Naito, Science and Technology of Advanced Materials 23, 66-75 (2022).
[6] R. Tamura, G. Deﬀrennes, K. Han,
T. Abe, H. Morito, Y. Nakamura,
M. Naito, R. Katsube, Y. Nose, andK. Terayama, Science and Technology ofAdvanced Materials: Methods 2, 153-161(2022).

335

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

Study on Novel Algorithm for Ising Machines
Shu Tanaka
Department of Applied Physics and Physico-Informatics, Keio University
Hiyoshi, Yokohama, Kanagawa, 223-8522
Ising machines are special computers forproperties of the proposed method were not

combinatorial optimization problems. The Ising

clarified. We have used the method of amachine expresses the objective function andprevious study [2,3] to study the slow

constraints of a combinatorial optimization

relaxation due to entropy effects, temperature-

problem in terms of an Ising model, anddependent effective interactions, and local

searches for low-energy states of the Ising

magnetic fields. As a result, it was confirmed

model. Since the Ising machine operates on thethat there is a change in the time scalebasis of stochastic behavior, it is necessary to

depending on the bit-width to be reduced [4].

develop various internal algorithms to improve
the solution accuracy. In this research project,(II) Black-box optimization using Ising

we investigated algorithms to improve themachines

performance of the Ising machine from theAn Ising machine is a piece of hardware thatviewpoint of statistical mechanics.

takes an Ising model or equivalent Quadratic
Unconstrained Binary Optimization (QUBO)(I)

Dynamical

reduction

properties

algorithms

of

for

bit-width

representation as an input format and operates

simulated-

to find its ground state. Ising machines havebeen

annealing-based Ising machine

appliedinvarioussituations

tohardwarecombinatorial optimization problems that canlimitations. We have focused on cases where

be represented by Ising models or QUBOs. On

thethe other hand, a method to use Ising machinesIsingmachinespossiblehavevaluesseveraloftheinteraction

coefficients and local field coefficients in theforIsing model are limited and have proposed a

discrete optimization problems where themethod to make Ising machines work properly

objective function is not explicitly given, waseven in the presence of such limitations [1].

proposed by us in 2020 [5]. We extended it andReference [1] confirms the validity of the

proposed a method to perform integer-valued

proposed method by comparing the groundblack-box optimization with Ising machines [6].

states. On the other hand, the dynamic

336

black-box

discrete

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

References

arXiv:2304.12796 (2023).

[1] D. Oku, M. Tawada, S. Tanaka, and N.

[5] K. Kitai, J. Guo, S. Ju, S. Tanaka, K. Tsuda,

Togawa, IEEE Trans. Comp. 71, 223 (2020).

J. Shiomi, and R. Tamura, Phys. Rev. Res. 2[2] S. Miyashita, S. Tanaka, M. Hirano, J. Phys.

(2020) 013319.

Soc. Jpn. 76, 083001 (2007).

[6] Y. Seki, R. Tamura, and S. Tanaka,[3] S. Tanaka and S. Miyashita, J. Phys. Soc.

arXiv:2209.01016 (2022).

Jpn. 78, 084002 (2009).
[4] S. Kikuchi, N. Togawa, and S. Tanaka,337

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

Charge distribution in metallic and charge order phase inorganic crystals
Masatoshi SAKAI, Hiroki WATANABE, Ryosuke ANDO, Ryo WATANUKI
Department of Electrical and Electronic Engineering,
Chiba University, 1-33 Yayoi-cho Inage, Chiba 263-8522

Ourresearch

group

conductsof the parameter set, which can obtain thedevelopments of novel electronic devices

charge order phase and metallic phase. Next,

using organic charge order materials. In our

we are planning to try calculating real-time

previous work, we observed obvious fieldevolution from the insulator to the metallic

effect transistor (FET) characteristics of -

phase and vice versa..

(BEDT-TTF)2I3[1] with 4-probe electrical

Furthermore, we also proceeded themeasurement and small FET modulation of -

calculation of the charge order phase of -(BEDT-TTF)2PF6 [2] with 2-terminal electrical(BEDT-TTF)2PF6 crystal. This calculation

measurement. To analyze these experimental

started from the first-principle calculation

results, we are trying to calculate the(Quantumelectronic state of -(BEDT-TTF)2I3 and -parameters using RESPACK. By using the(BEDT-TTF)2PF6. This report is a interim

obtained parameters, we tried to represent thereport of this research work in progress.

charge distribution of the experimentally

Espresso),

and

we

obtained

At first, about the electronic state of -observed charge order. At present, we are(BEDT-TTF)2I3 crystal, we adopted new

making a superlattice structure to represent

parameter sets[3] of t, U, V because thethe real charge order. We are also trying toprevious parameter set could not reproduce

calculate the process and probability between

the experimentally observed charge orderthe insulator and metallic phase of -(BEDT-

structure. This time we adopted the parameterTTF)2PF6 to analyze our phase transition

set considering the breaking of inversion

device characteristics.

symmetry.

By

adopting

these

novel

parameters, we obtained calculated charge[1] under preparation

distribution of which the charge distribution[2] M. Sakai et al., Nanosci. and Nanotechnol. 16,

trend agrees with that of the experimentally

3267 (2016).

obtained one. And we could marshal the trend[3] D. Ohki et al., Phys. Rev. B 107, L041108 (2023).

338

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

Magnetic structures on the thin ﬁlms made by the
dipolar interaction and the frictional force stemmingfrom these structures
Hisato Komatsu
Kyushu University, Kasuga, Fukuoka 816-8580 ∗Masamichi Nishino and Yoshihiko Nonomura
National Institute for Materials Science, Tsukuba, Ibaraki 305-0044

We applied for the usage of the ISSP Supercomputer planning to investigate the magnetic
structures of the dipolar system. However, wehave not got notable results on this subject.
Instead, we studied a theoretical model of themagnetic friction.
Magnetic friction, the frictional force generated from the magnetic interaction between
spin variables, has been extensively studied inrecent years. These studies aimed not only toinvestigate the behavior of the magnetic systems, but also to reveal the microscopic mechanisms of the friction. For this purpose, many
types of theoretical models have been proposed
and investigated[1, 2].
In the case of the normal solid surfaces,
there is a well-known empirical law called theDieterich-Ruina law[3, 4]. In our previous
studies, we proposed models of magnetic friction that seemed to obey this law, at least in
the steady state[5, 6]. These studies revealed
that magnetic structures can behave as a kindof potential barrier that prevent the lattice motion. In the study of this academic year, we investigated a simpliﬁed model of the magnetic
friction to consider the ﬁnite-size eﬀect, which∗

Present address : Shiga-university, Hikone, Shiga
522-0069

cannot be discussed in our previous studies.
Speciﬁcally, we introduced a model composedof N Ising spin variables {σi } interacting witheach other by the following Hamiltonian:H=−

J(x) ∑σi σj = −N J(x)m2 .
N i,j

(1)

Here, x is a real value representing the shift ofthe lattice, m is the magnetization per spin:∑m ≡ i σi /N , and J(x) is the periodic function of x. For simplicity, we impose a periodic boundary condition on x. As a concrete
form of J(x), two types, A and B, are considered. Type-A is a piecewise linear function,and type-B is a sinusoidal function.
We introduce the time development of spin
variables by the Glauber dynamics, and deﬁne one Monte Carlo step (MCS) as the unit
time. In addition, the time development of x
under this Hamiltonian and external force Fex
are introduced as the following overdamped
Langevin equation:
γN

dx
∂H √
= Fex −
+ 2γN T R(t),
dt
∂x

(2)

where R(t) is the white Gaussian noise,⟨R(t)R(t′ )⟩ = δ(t − t′ ). Here, we introduce
coeﬃcient N , considering the situation that N
spins move simultaneously.

339

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

References(a) 1.0

10-2
10-3

|m|

We investigated this system by numericalsimulation using ISSP Supercomputer, andfound that this system shows two diﬀerent behaviors depending on the damping constant γ.
When γ is small, the stick and slip states are
separated as the metastable states, and theprobability that each state appears determines
the relation between Fex and the lattice velocity v. When γ is large, on the other hand, suchseparation does not appear and v is determinedby the thermal activation process. We calculated the histogram of m and x like Fig. 1 andthe Fex -v relation, by the simulation. Furthurmore, we also estimated these quantities undersuﬃciently large or small γ theoretically using
approximations, and compared it with the result of the simulation.

10-4

0.5
10-5
10-6

0.0

10-7

-3 -2 -1

0

1

2

3

x

(b) 1.0

10-2[1] D. Kadau, A. Hucht, and D. E. Wolf,
Phys. Rev. Lett. 101, 137205 (2008).

|m|

10-3

[2] M. P. Magiera, L. Brendel, D. E. Wolf,and U. Nowak, Europhys. Lett. 87, 26002(2009).

[5] H. Komatsu, Phys. Rev. E. 100, 052130(2019)[6] H. Komatsu, Phys. Rev. E. 102, 062131(2020)[7] B. N. J. Persson, O. Albohr, F. Mancosu,
V. Peveri, V. N. Samoilov, and I. M. Sivebaek, Wear 254 835 (2003).

0.5
10-5
10-6

[3] A. Ruina, J. Geophys. Res. 88, 10359(1983)[4] J. H. Dieterich, Tectonophysics 144, 127(1987)

10-4

0.0

10-7

-3 -2 -1

0

1

2

3

x
Figure 1: Examples of the histogram of themagnetization m and shift of the lattice x under (a)small γ and (b)large γ. In the case that
γ is small, stick and slip states are separated
as the metastable states like graph (a). Thesegraphs are the results of the type-A model,
type-B model also shows the similar result.

340

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

Numerical Study of One Dimensional Frustrated
Quantum Spin Systems
Kazuo HIDA
Professor Emeritus, Division of Material Science,
Graduate School of Science and Engineering,
Saitama University, Saitama, Saitama 338-8570

1

Model

We investigate the ground-state phases of
mixed diamond chains described by the following Hamiltonian:
H=

L [
∑

(1)

(1 + δ)S l (τ l(1)

+ τ l−1 )l=1(2)+ (1 − δ)S l (τ l
(1)

(2)

(1) (2)]+ τ l−1 ) + λτ l τ l

, (1)

(2)

where S l , τ l and τ l are spin operators with(1)
(2)
magnitudes Sl = τl = 1/2 and τl = 1,
respectively. The number of unit cells is denoted by L, and the total number of sites
is 3L. The lattice structure is depicted inFig. 1. We consider the region λ ≥ 0 and
(1)

τl
1+δ

periodicity of the ground state. Since the spin
operators τ (1) and τ (2) have diﬀerent magnitudes, however, it is not natural to assume
δ = 0 in a realistic situation. Hence, we examine more general cases with δ ̸= 0 in thepresent work.

(1)

1+δ

S=τ =1/2

Sl+1
λ
Sl
1−δ
1−δ
(2)
τl

(2)

τ =1

Figure 1: Structure of the diamond chain investigated in this work.
(1)

(2)

1 ≥ δ ≥ −1. For δ = 0, τ l + τ l commutes with the Hamiltonian (1) for all l. In[1], we made use of this property to determine
the ground-state phase diagram. In addition to
the Lieb-Mattis(LM) type ferrimagnetic phaseand the nonmagnetic Tomonaga-Luttinger liquid (TLL) phase, an inﬁnite series of ferrimagnetic phases with spontaneous magnetizations
msp = 1/p, where p(= 1, 2, ..., ∞) is the spatial

341

2

Analytical Results

Let us start with the analytic consideration ofseveral limiting cases.
(1)

1. λ ≫ 1 + δ, 1 − δ: All pairs of τ l and(2)
τ l form spin-1/2 states. If the eﬀective
interaction between them is antiferromagnetic, the ground state is a TTL state.
If it is ferromagnetic, the ground stateis the LM-type ferrimagnetic state with
msp = 1. Using the perturbation calculation with respect to (1 ± δ)/λ, we ﬁnd thephase boundary between these two phases
is given by δ = 0.6.
2. λ ≤ 0、1 ± δ > 0: The system is unfrustrated and the ground state is the LMtype ferrimagnetic state with msp = 1.
A similar ground state is expected for
0 < λ ≪ 1 ± δ.
3. 1 + δ ≃ 0 and λ ≃ 0 : For 1 + δ = λ = 0,(2)
Sl snd τ l form a ferrimagnetic chain(1)
with msp = 1/2 and τ l are free spinswith magnitude 1/2. For small 1 + δ andλ, we can determine the phase boundary between the LM-type ferrimagnetic

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

phase with msp = 1 and the nonmagnetic phase as λ/(1 + δ) ≃ 0.738 by the
perturbation calculation with respect tothese small parameters 1 + δ and λ combined with the numerical exact diagonalization (NED) calculation of the ferrimagnetic chain.
4. 1 − δ = 0 or λ = 0 : If 1 − δ = λ = 0,
(1)
S l and τ l form a spin-1/2 antiferromag(2)netic Heisenberg chain aand τ l are freespins. If λ > 0 and δ = 1, the ground stateis the LM-type ferrimagnetic phase with
msp = 1. The same is true if δ < 1 and
λ = 0. However, if λ > 0 and δ < 1, these
two interactions frustrate, and a nontrivial ground state is realized.

3

Numerical Results

ferrimagnetic phases with msp = 1 are not asingle phase, but are separated by the partial
ferrimagnetic phase and nonmagnetic phase.
The nonmagnetic phase continues to the TLL
phase for large δ. Hence, the whole nonmagnetic phase is likely to be a TLL phase. However, since the sum of the spin magnitudes in aunit cell is an integer, the Lieb-Schultz-Mattis
theorem does not exclude the gapped phase.
The fate of the inﬁnite series of ferrimagnetic
phases found for δ = 0 in [1] is also investigated. The λ-dependence of the spontaneous
magnetization for δ = 0.1 is shown in Fig. 3.
It seems that the width of the ferrimagnetic
phase with msp = 1/2 (p = 2) remains ﬁnite,while those with p ≥ 3 are smeared out andthe corresponding region turns into the partial
ferrimagnetic phase suggesting the fragility of
the ferrimagnetic phases with large p.
msp

We have carried out the NED calculation for
L = 4, 6 and 8 and the DMRG calculation for
L = 24 and 48. The magnetization curves
and spontaneous magnetization are calculated
to obtain the ground-state phase diagram ofFig. 2. It should be remarked that the two

p=1

0.5

msp

δ=0.1
L=48

1

p=1

1

p=2

0.5

p=2

p=3

p=3

p=4

~0.807λc(∞)−

p=4

p≤38

0
0.6

0.7

0.8

λ 0.9

0
0.6

0.7

λ

0.8

2
λ

1+δ

λ
1/2

1/2

λ
1−δ

nonmag (gapless?)

1

Figure 3: Left panel : λ-dependence of the
spontaneous magnetization in the ground statefor δ = 0.1 calculated by the DMRG method
for open chains with L = 48. Right panel: The
corresponding ﬁgure for δ = 0 taken from ref.
[1].

1+δ

1−δ

msp=1

λ
1−δ

msp=1/2
partial
ferri

msp=1
1+δ
λ
1−δ

0
−1

0

δ

1

References

Figure 2: Ground-state phase diagram. The
open circles are the phase boundaries estimated from the NED data extrapolated to thethermodynamic limit. The ﬁlled circles are thephase boundaries estimated from the DMRG
data for L = 24 and 48. The curves are guides
for the eye. The insets are schematical spin
structures in diﬀerent phases.

342[1] K. Hida, J. Phys. Soc. Jpn. 90, 054701(2021).

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

Novel phenomena in mismatched multilayer systems: Search
for new candidate materials and analysis of electronic
structures
Toshikaze KARIYADO
Research Center for Materials Nanoarchitectonics,
National Institute for Materials Science, Namiki, Tsukuba, Ibaraki 305-0044
a way to valleytronics applications.

Recently, artificially stacked 2D systems
attract much attention as a playground for

In the analysis of monolayer and bilayerelectronic structure manipulations. In particular,

BC3, we make full use of the ab initio method,

systems with nanoscale moiré pattern induced

using mainly Quantum Espresso and sometimes

by

are

OpenMX. A characteristic feature of moiréinteresting, since the nanoscale patterns tend tobilayer systems is that the relation between two

have striking effects on low energy band

layers depends on position, and it is important

structures.

to know the local crystalline and electronicslightmismatchbetweenlayersIn this project, we propose twisted bilayer

structures to build a model to handle theBC3 as a system for realizing interesting valley

nanoscale moiré pattern. This means that we

related physics [1]. BC3 is a graphene

need to work out computations on various local

derivative with a honeycomb network of atoms

structures[2]. The monolayer BC3 is a three-valley

positions in the nanoscale moiré unit cell. Then,

semiconductor with its conduction minima atlarge computational resources are required, andthree M-points in the hexagonal Brillouin zone.

the numerical efficiency is important. For the

Each of the valleys shows a 2D parabolic

numerical efficiency, the Quantum Espresso

dispersion around the conduction minima,

package offers an efficient parallelized version.

corresponding

to

the

variousthough there is anisotropy. It is found that uponmaking twisted bilayers, the 2D dispersion issqueezed into quasi-1D one. Interestingly, theReferences

quasi-1D directionality depends on the valley[1] T. Kariyado, Phys. Rev. B 107, 085127

degrees of freedom. That is, twistronics applied(2023).

on bilayer BC3 leads to valley dependent band[2] H. Tanaka et al., Solid State Commun. 136,

structure manipulation, which potentially paves

22 (2005).

343

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

Structure analysis of borophene by using 2DMAT
Akari TAKAYAMA
Department of Physics,
Waseda University, Ohkubo, Shinjuku-ku, Tokyo 169-8555
TRHEPD were performed by the structure-

Introduction
Atomic-layer sheet composed of group XIV

analysis program, “2DMAT” [4], which iselements such as graphene and silicene haveincluded as standard software in ISSP’s

been actively studied. These atomic-layer sheets

supercomputer.

form a honeycomb structure and are expected tocomputer simulations on supercomputer is very

be applicable to various technologies because of

useful for structure analysis of large periodic

their excellent electron mobility due to theirstructure and surface superstructure with lattice

peculiar band dispersion following the Dirac

reconstruction because they require manyequation. Recently, an atomic-layer sheet

variables for analysis.

Usingmassiveparallelmaterial of boron (the group XIII element),
Results and discussion

called “borophene” is attracting attention,
especially after the success of its fabrication on

The TRHEPD experiment consists ofmetal substrates. Unlike atomic-layer sheets ofmeasuring a series of the diffraction patterns forgroup XIV elements that form a honeycomba fixed incident azimuthal direction at various

lattice, flat layers of boron are theoreticallyglancing angles (). Here, the rocking curve isstable when composed of a triangular lattice anddefined as the diffraction intensity of the (00)hexagonal hollows. Experimental observations

spot plotted as a function of . In the structuralof epitaxial borophene on Ag(111) have been

analysis, the experimental rocking curves are

consistent with the theoretical structural models

compared with those calculated for various[1, 2], but the atomic structure has not been

structural models. Figure 1 shows the rocking

examined directly. More interestingly, scanning

curves of (a) a pristine Ag(111) surface and (b)tunneling microscopy (STM) observations of athe borophene grown on it, measured under thehoneycomb structure on an Al(111) substrate
were reported [3], but the detailed atomic
structure of this honeycomb lattice is unclear.
Thus, there is a strong need to conduct structural
analysisbyanappropriate

diffraction

experiment on the surface.
In our research, we studied the atomic
structure of borophene on Ag(111) and on
Al(111) by means of total-reflection high-energy

Fig. 1: Rocking curves under one-beam

positron diffraction (TRHEPD). For structure

condition with calculated curves for (a)analysis, calculations of the rocking curves forAg(111) and (b) borophene on Ag(111).

344

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

one-beam condition. In the one-beam condition,
the beam is incident along an off-symmetricdirection, thus the rocking curve in the one-beamcondition essentially gives the information on
the atomic positions in the out-of-plane direction.
The experimental profile was reproduced by the
proposed model [2]; layer distance between Ag

Fig. 2: Rocking curves under one-beamand borophene is 2.4 Å and the magnitude of

condition for (a) before and (b) after boron

bucking for borophene 0.04 Å. This indicatesdeposited on Al(111).

that there is no buckling, and we have directly
confirmed that the borophene is a flat sheet [5].
Figure 2 shows the rocking curves of (a) aReferences

pristine Al(111) surface and (b) the borophene[1] R. Wu, et al., Nature Nanotechnology, 14grown on it under the one-beam condition(2019) 44.

obtained from the experiment and from

[2] G. P. Campbell, et al., Nano Letters, 18

calculation of Al(111). After boron deposition,(2018) 2816.

the peak position shifted slightly, but the whole[3] W. Li, et al., Science, 63, (2018) 282.

shape of the rocking curve did not change much.

[4] Y. Motoyama, et al., Computer Physics

According to angle-resolved photoemission

Communications, 280 (2022) 108465.

spectroscopy (ARPES) measurements and band[5] Y. Tsujikawa, AT et al., Molecules. 27

calculations for borophene on Al(111) substrates,(2022) 4219.

it is predicted that the borophene forms AlB2 at[6] D. Geng, et al., Physical Review B 101

the top surface [6]. However, the rocking curve(2020) 161407(R).

of the AlB2 model differs significantly from our[7] A. B. Preobrajenski, et al., ACS Nano, 15,

experimental results, which is inconsistent with

(2021) 15153.

the results of this study. One possible reason forthe difference between our experiment andprevious studies is the effect of ripple structure[3]. Since the ripple structure was not consideredin this analysis, it is possible that the model usedin the analysis did not correctly reproduce theexperimental results. It has also been pointed outthat borophene on Al substrate forms longperiod moiré structures due to lattice mismatchwith the substrate [7]. To determine the structureof honeycomb borophene on Al substrate, wewill analyze the giant lattice system utilizing
parallel computation on a supercomputer.

345

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

Magnetic structures of multiple-Q orders ininversion-symmetric Hubbard Models
Takashi UCHIDA
Hokkaido University of Science
4-1, 7-15, Maeda, Teine-ku, Sapporo 006-8585

In the last decade, multiple-Q orders inmagnetic systems have attracted attention because they sometimes reveal topologically protected magnetic structures such as magnetic
skyrmions. Typically, these structures are realized in the Dzyaloshinskii-Moriya (DM) interaction driven systems under magnetic ﬁeld.
The multiple-Q orders can also be realized ininversion-symmetric systems where the DM interaction is absent [1, 2].
In order to investigate the nature of themultiple-Q orders in inversion-symmetric itinerant systems, we have explored the magnetic
structures of the triangular-lattice single-band
Hubbard model by means of the molecular spin
dynamics (MSD) method [3]. The present formulation of the MSD adopts the static approximation to the functional integral method andthe isothermal molecular dynamics technique,and reduces to the generalized Hartree-Fockapproximation at the ground state.
In the present work, we have performed the
MSD calculation with use of the global charge
neutrality condition instead of using the local charge neutrality condition [3]. This has
the advantage that the charge transfer among
the atoms is allowed in the magnetic structure formation. In the numerical calculations
the most time-consuming process is the magnetic force calculation at each time step, where
the local electronic structures are calculated inreal space by means of the recursion method.
We have adopted the MPI parallel calculation
scheme and found it eﬀective in saving bothcomputing time and CPU resources.
We have performed magnetic structure
calculations on a hexagonal supercell with
972 lattice points, which is embedded in alarge cluster consisting of 6 such supercells,
each of which are connected by the periodic
boundary condition. Under zero magnetic
ﬁeld and the ﬁxed value of the temperature
T /t = 0.0005, we have explored the magnetic structures changing the Coulomb interaction strength U/t and the electron number
n along the antiferromagnetic-ferromagnetic
boundary. We have found that the stablestates for U/t < 4.0 and n = 1.45 ∼ 1.53
are the 3Q magnetic orders with the Q vectors
pointing in the direction of three axes of thetriangular lattice and that they are accompanied by 3Q charge density waves (CDW) withthe same Q vectors as the magnetic ones. The
present results with use of the global charge
neutrality condition reveal the enhancement of
the amplitude of the main 3Q magnetic order as compared to those with use of the local
charge neutrality condition, showing that the
CDW contributes to stabilize the magnetic 3Q
order.

References[1] T. Okubo, S. Chung, and H. Kawamura:Phys. Rev. Lett. 108 (2012) 017206[2] Y. Kakehashi: J. Phys. Soc. Jpn. 89 (2020)
094710.
[3] Y. Kakehashi, S. Akbar, and N. Kimura:Phys. Rev. B57 (1998) 8354.

346

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

Development of integrated interface of eigensolvers
Rokko and application to quantum spin systems
Tatsuya Sakashita
Center for Quantum Information and Quantum Biology, Osaka University
1-2 Machikaneyamachō, Toyonaka, Osaka 560-0043

To establish universal exact diagonalization
package for quantum lattice models including the Heisenberg-Kitaev model, we focusedon developing integrated interfaces for eigensolvers, “Rokko”[1].
In Rokko, we implemented the integrated interfaces for the following types:• Serial solvers for dense matrices (Eigen3,
LAPACK)• MPI parallelized solvers for dense matrices (EigenExa[2], ELPA[3], Elemental[4],ScaLAPACK)• MPI parallelized solvers for sparse matrices (Anasazi in Trilinos[5], SLEPc[6]) tocover matrix representations below:– CRS (Compressed Row Storage)– Matrix-free method (the method togive matrix-vector product routines
to solvers)• Unit and integrated test programs by
GoogleTest• Install scripts of eigensolvers for various
architectures
We prepare a paper to report design policy, software structure, and usage examples ofRokko.

References[1] T. Sakashita, R. Igarashi, Y. Motoyama,
T. Okubo, and S. Todo.
Repository of Rokko.
https://github.com/
t-sakashita/rokko.git, 2012.
[2] T.
Imamura,T.
Hirota,andT. Fukaya.
EigenExa web page.
https://www.r-ccs.riken.jp/labs/
lpnctrt/projects/eigenexa/, 2021.

Rokko has the following features:[3] ELPA Consortium. ELPA (Eigenvalue
soLvers for Petaflop Applications webpage). http://elpa.rzg.mpg.de, 2013.

• Integrated interfaces for eigensolvers andmatrices, independent of individual eigensolver libraries[4] J. Poulson. Distributed-memory dense linear algebra Elemental web page. https://github.com/elemental/Elemental, 2013.

• Rokko’s interfaces are implemented byutilizing factory. It enables the user todynamically select a solver.

[5] M. A. Heroux, R. A. Bartlett, and V. E.
Howle. Trilinos Project web page. https://trilinos.github.io, 2003.

• C, Fortran, and Python bindings of Rokko[6] V. Hernandez, J. E. Roman, and V. Vidal.
SLEPc web page. http://slepc.upv.es,
2002.

• Automatically detecting libraries by using
CMake in building Rokko

347

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

Systematic understanding of the slip behaviors ofslow and fast earthquakes and the eﬀect of the upper
value of porosity
Takehito SUZUKI
Department of Physical Sciences, Aoyama Gakuin University
5-10-1, Fuchinobe, Chuo-ku, Sagamihara, Kanagawa 252-5258

Two qualitatively diﬀerent behaviors ofearthquakes have been known. One is thefast earthquakes, which we feel naturally. The
other one is slow earthquakes, which generate
negligible seismic waves. Although the slow
earthquakes are not disastrous, they are considered to sometimes change to the fast earthquakes after several repetitions. Whether thistransition occurs or not is a fatal problem tothe human society. Systematic understanding
for the slip behaviors of the slow and fast earthquakes is aimed in this study.
The fault rocks can be considered as porous
media, and the pores within the media are usually considered to be ﬁlled with water. If the
dynamic earthquake slip (frictional slip) occurs, we can consider the interaction among
heat, ﬂuid pressure, and porosity [1, 2]. The
slip is interpreted as the deformation aroundthe contact area, called slip zone, which has aﬁnite width wh . During the dynamic slip, thegeneration of pores [frictional heating] reduces[raises] the ﬂuid pressure in the slip zone, inducing the increase [reduction] in the friction
stress and the reduction [increase] in the slip
velocity [1, 2]. The ﬂuid-pressure increase due
to frictional heating will be called FH eﬀect,
whereas the ﬂuid-pressure decrease due to thepore generation is referred to as PG eﬀect. The
slow and ordinary earthquakes are dominatedby the PG and FH eﬀects, respectively.
The spring-block model and the interaction
between three quantities are employed to in-vestigate both the slow and fast earthquakes.
This system repeats the dynamic slip and cessation. The healing of pores in the cessation
time is also taken into account.
With this system, we found an importantfunction determining the system behavior [3].
For the derivation of the function, we employed
three important assumptions: (i) we consider asingle block, (ii) the diﬀusions of heat and ﬂuid
are neglected, and (iii) the change in the bulk
modulus of the medium, M0′ , due to porosity change is neglected; i.e., M0′ (ϕ) = M0′ (ϕ0 ),where ϕ0 is the porosity at the onset of the slip.
With these assumptions, we found a functionof the slip uf ,1
F (uf ) = kp u2f + µstat (σn0 + p0 )uf
2
[ (
)
1
M0′ α00−µslid(1 − e−γuf )σ n + p 0 − α0
−
γ
γ
ϕ∞
]

M ′ ϕ∞+ α0 0(1 − e−α0 uf /ϕ∞ ) , (1)−γϕ∞where kp is the spring constant, µslid and µstat
are the sliding and maximum static friction coeﬃcients, respectively, σn0 is the normal stressacting on the interface, p0 is the ﬂuid pressure at the onset of the slip, γ is the positive
constant, and α0 and ϕ∞ are the positive constants appearing in the porosity evolution law
(

)

∂ϕ
ϕ
wh
wh
= α0 v 1 −[H(y+ )−H(y− )].
∂t
ϕ∞
2
2
(2)

348

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

Two solutions whose magnitudes are different in two orders emerge for the equation
F (uf ) = 0, and the small and large ones correspond to the slow and fast earthquakes, respectively. Detailed analysis shows that thisbehavior of the solution can be interpreted as
the phase transition of the ﬁrst order [3].
We should now remove above-mentioned
three assumptions to treat more realistic situation. We performed numerical calculations
with ISSP supercomputer. The number ofblocks is 40, and the diﬀusions of heat and ﬂuid
are included. Moreover, from Eq. (2), the assumption (iii) should be valid if the value of
ϕ∞ is negligibly small. We consider three values for ϕ∞ , and the condition where we canapply the analytical treatment will be clariﬁed.
If we assume ϕ∞ = 0.01, only the fast earthquakes occur, and the slow earthquakes do notoccur (Fig. 1a). This is because the porosity change is negligible, and the ﬂuid pressurechange due to the PG eﬀect is negligibly small.
If we assume ϕ∞ = 0.1, both the fast and slow
earthquakes occur (Fig. 1b) because both theFH and PG eﬀects work in this case. If we assume ϕ∞ = 0.9, the fast and slow earthquakes
occur (Fig. 1c). The FH and PG eﬀects work
also in this case.
Additionally, the analytical treatment using
F (uf ) was found to be eﬀective for ϕ∞ = 0.01
and 0.1, whereas it was not eﬀective for theslow earthquake observed for ϕ∞ = 0.9. The
value 0.9 may be too large to employ theabove-mentioned assumption (iii). The function F (uf ) is useful for predicting the dynamic
slip behaviors of the natural fault zones if ϕ∞is much smaller than unity.

References

(a)

(b)

(c)

Figure 1: The slip velocity of the 20th blockwith (a) ϕ∞ = 0.01, (b) 0.1, and (c) 0.9.

[1] T. Suzuki, and T. Yamashita: J. Geophys.
Res. 119 (2014)[2] T. Suzuki: Phys. Rev. E 96 (2017)[3] T. Suzuki, and H. Matsukawa, in prep.

349

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

Shape dependence of macroscopic friction between solids
Michio OTSUKI
Graduate School of Engineering Science,
Osaka University, Machikaneyama, Toyonaka, Osaka 560-8531
Amonton's law states that the static friction

static friction coefficient and the behavior of

coefficient does not depend on the external

precursor slip [3].

pressure or the shape of the object. Recent
studies have shown that Amonton's law
systematically breaks down for macroscopic
objects due to precursor slip before the onset ofbulk sliding [1, 2]. However, it is unclear
whether the results apply to more realistic 3D
systems since the analysis in previous works is

Fig. 1 Schematic of a 3D viscoelastic block

restricted to 2D systems.

on a fixed rigid substrate.

First, we have numerically investigated thesliding motion of a 3D viscoelastic block on a

In addition, we have studied the effect of therigid substrate using finite element simulation as

shape of the frictional interface. Frictional

shown in Fig. 1. The maximum computational

interfaces of some industrial products, such as

size requires 500,000 nodes and 10 billion timeshoe soles and tires, are often grooved. However,steps, which takes 100 hours using our code

there is a need to understand how grooves affectcompiled with the Intel Fortran compiler and a

the frictional properties. Therefore, we haveparallel computation of 2,000 MPI processes on

numerically investigated the slip of objects that

the CPU server. The block exhibits stick-slip

are longitudinally grooved as shown in Fig. 2.

motion characterized by the static friction

The static friction coefficient is a decreasing

coefficient. The static friction coefficient

functiondecreases with increasing external pressure,

Furthermore, our theoretical analysis shows that

length, or width of the block, which contradicts

the decrease is caused by a reduction in the

Amonton's law. The precursor slip occurs in theeffective viscosity due to the groove. These

2D interface between the block and substrate

results provide new insight into the design of

before bulk sliding. The instability analysis of agrooves for friction control.

simplified model explains the decrease in the

350

of

groovewidth

anddepth.

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

References[1] M. Otsuki and H. Matsukawa, Sci. Rep. 3,
1586 (2013).
[2] Y. Katano, K. Nakano, M. Otsuki, and H.
Matsukawa, Sci. Rep. 4, 6324 (2014).
[3] W. Iwashita, H. Matsukawa, and M. Otsuki,
Sci. Rep. 13, 2511 (2023).

Fig. 2 Schematic of a 3D viscoelastic
block with grooves block on a fixed rigid
substrate.

351