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

3.4

Cooperative Phenomena
in Complex Macroscopic Systems

205

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

Effects of cavitation on soundwaves
Yuta ASANO, Hiroshi WATANABE, and Hiroshi NOGUCHI
Institute for Materials Research, Tohoku University Katahira, Aoba-ku, Sendai, Miyagi 9808577
Department of Applied Physics and Physico-Informatics,
Keio University, Hiyoshi, Kohoku-ku, Yokohama, Kanagawa 233-8522
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581

Cavitation is a flow phenomenon involving

experiments and numerical calculations, thethe formation of bubbles, caused by local

mechanism remains thoroughly understood. In

pressure-drop in a liquid. Cavitation causes

this project, we investigate the effect of phase

performance degradation, noise, vibration, andtransitions on sound waves by molecularerosion in fluid machinery such as ship

dynamics (MD) simulations [1]. In particular,

propellers and turbomachinery. On the other

we clarify the differences in the effects of first-

hand, the shock waves and localized high-and second-order phase transitions on soundtemperature and -pressure fields generated bywaves.

the bubble collapse are applied to a wide range

The fluids are modeled by a monoatomic

of applications, including cleaning, processing,molecule whose interaction is defined by thechemical reaction acceleration, and medicine.

smoothed-cutoff Lennard-Jones potential. The

Therefore, it is vital to understand thesystem is a rectangular parallelepiped, and thecavitation mechanisms for many fields using

periodic boundary conditions are applied to y-liquid-gas multi-phase flows. In particular,

and z-directions. In x-direction, the oscillating

understanding the effects of cavitation on soundand stationary walls are placed at both ends.

waves is essential for the development of

The total number of particles is about 10

applied technology, such as the realization of

million. We investigate the differences in thehighly efficient chemical reactors. However,waveform change by approaching the phasecavitation

complextransition temperature for the liquid near thephenomenon in which a tremendous number offirst-order transition region and supercriticalbubbles repeatedly generate, grow, split andfluid near the critical point. These waveforms

coalesce, and disappear, making its analysis aare compared with the numerical solution ofchallenging

Burgers' equation.

isanproblem.

extremely

Despitenumerous

206

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

In the first-order transition region, thewaveform varies continuously until just above

References

the phase transition point, in good agreement[1] Y. Asano, H. Watanabe, and H. Noguchi,

with the Burgers' equation. When the phase

arXiv:2112.07139.

transition

occurs,thewaveformchanges

discontinuously due to the appearance of thegas phase, and the Burgers' equation is notapplicable. On the other hand, in the continuoustransition region, no discontinuous changes
occur. Soundwave decay due to the density
fluctuations are significant near the criticalpoint. However, the Burgers' equation isapplicable

even

for

such

large

densityfluctuations.
Finally, we investigated the bubble effects on
the soundwaves by using 60 million particles.
We observed a bubble generation and growthnear the oscillating wall, movement in thesoundwave propagation direction, and thendisappearance. Bubbles repeat this long time(about 30-million-time steps) cycle. We foundthat the bubble motion affects the pressurewave.
The present computational scale could only
treat the interaction between a single bubble
and sound wave. The bubble generation causes
a drastic change in soundwave, known as achange in sound speed, which is caused by theappearance of many bubbles. In order to reveal
the peculiar phenomena, it is necessary to
perform calculations for much larger systems.

207

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

Coarse-Grained Molecular Dynamics Study of Biomolecular
Assemblies and Nanocolloids
Wataru SHINODA
Research Institute for Interdisciplinary Science,
Okayama University, 3-1-1 Tsushima-naka, Kita-ku, Okayama 700-8530
a series of membrane protein assemblies,

We have extended the molecular library
within the SPICA coarse-grained (CG) force

scramblases,

andpoliovirus

capsids

field (FF), which was originally developed fordemonstrated the good performance of SPICA

lipids and surfactants. In this year, we haveFF.

added a new class of lipids, proteins, andA series of MD simulations of Hepatitis B

nucleic acids to the library. Especially, we

virus (HBV) capsid have also been carried out

improved the CG protein model to facilitate

at AA and CG resolutions.[2] In this year, wemolecular dynamics (MD) simulation with ahave constructed a molecular model of the

reproduction of multiple properties from both

entire HBV virion particle including theexperiments and all-atom (AA) simulations.[1]envelope.

The side chain analogs reproduced the transfer

membrane protein within the envelope and thefree energy profiles across the lipid membrane

capsid spike protein plays an important role in aand showed reasonable association free energy

strong anchoring of envelopes around the(potential of mean force) in water compared tocapsid. Due to a significant amount of

those from AA-MD. The CG-MD study of themembrane proteins in the envelope reduces theadsorption

150

lipid diffusion within the envelope, which is

peptides/proteins into the membrane gave

supposed to help increasing the stability ofcorrect predictions of the penetration depth andvirion particle.

orpenetration

of

about

The

interaction

betweenthe

We also analyzed the role of surfactants ininclination angle of peptides/proteins at thetransmembrane

the colloidal stability of nanoparticles in oil byregions, consistent with the OPM (Orientation

CG-MD.[3] Again, SPICA FF was used toof protein in membrane) database. The

observedimerization

nanoparticle aggregation through simulation. A

membrane

periphery

free

and

energies

of

some

the

pathwayofdetailsaggregationoftransmembrane helices in the lipid bilayer werefreefound to be comparable to those obtained fromdetermined for diﬀerent surfactant types andexperimental estimates. Application studies on

system conﬁgurations. The deformation of the208

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

surfactant layer surrounding the nanoparticle

M. Dixit, R. DeVane, C. MacDarmaid, G.

was shown to be the major source of colloid

Fiorin, M. L. Klein, W. Shinoda: J. Chem.

stabilization.

the

Theory Comput. DOI:10.1021/acs.jctc.1c01207.

surprising result that as the packing density of[2] K. Fujimoto, Y. Yamaguchi, R. Urano, W.

the surfactants is reduced, the role of the

Shinoda, T. Ishikawa, K. Omagari, Y. Tanaka, A.

surfactant is reversed from stabilizer to

Nakagawa, S. Okazaki: J. Chem. Phys. 155

nucleation agent.

(2021) 145101.

This

is

evidencedby[3] M. Griffiths and W. Shinoda: J. Phys. Chem.

References

B 125 (2021) 6315. [cover title][1] S. Kawamoto, H. Liu, Y. Miyazaki, S. Seo,

209

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

Finite temperature properties of generalized Kitaev
modelsTsuyoshi OKUBO
Institute for Physics of Intelligence, The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 133-0033

The Kitaev model on the honeycomb lattice
has been widely investigated as a fundamental
model to investigate the nature of the quantum spin liquid in two dimension[1, 2]. In the
Kitaev model, the interactions are Ising type,Siγ Sjγ , and the spin component γ = x, y, z isdetermined from the direction to the neighboring site. The ground state is known to be agapless spin liquid state, in which the nature
is described by Majorana fermions[1]. When
we apply a magnetic ﬁeld to the Kitaev system, the excitation spectrum described by themobile Majorana fermions opens a gap, andthe ground state changes to a topological spin
liquid, where a half-integer thermal Hall conductivity emerges in the low-temperature limit[1]. In the recent experiment on α-RuCl3 under a moderate magnetic ﬁeld, a similar halfinteger thermal Hall conductivity was observed
at a ﬁnite temperature [3], although its groundstate seems to be a magnetically ordered state
at zero magnetic ﬁeld due to oﬀ-diagonal and(long-range) Heisenberg interactions.
In this year’s project, we investigated a ﬁnite
temperature properties of the Kitaev modelwith oﬀ-diagonal interactions, in paticular concerning the possible thermal Hall conductivity
at a ﬁnite temperature. The model Hamiltonian is given asH=∑γ∈x,y,zHγ − ⃗h ·

∑
i

⃗i ,
S

(1)

Figure 1: A typical lattice setup and a mappting to a one-dimensional matrix product operator.
where for γ = z,∑
Hz =[KSiz Sjz + Γ(Six Sjy + Siy Sjx )⟨i,j⟩z+ Γ′ (Siz (Sjx + Sjy ) + (Six + Siy )Sjz )].
(2)
In the cases of γ = x, y, we consider similarinteractions with cyclic rotation of x, y, z components. 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 amatrix product operator (MPO), as shown in
Fig. 1.
The MPO representing a density matrix at agiven temperature was obtained by the exponential tensor renormalization grope (XTRG)
algorith [4], where we calculate the density matrix at an inverse temperature β, ρ(β), throughthe relationship ρ(β) = ρ(β/2)ρ(β/2) with theinitial condition ρ(β0 ) = exp(−β0 H) ≃ 1−β0 H

210

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

for the suﬃciently small β0 . Thus, when we
have an MPO representation of the Hamiltonian, a stadard compression algorithm for thetensor network can be used. When the bonddimension of the MPO is D, both of the computation cost and the memory consumption of
XTRG algorithm scale O(D4 ).
By such an MPO representation of the density matrix, we have successfully calculated thethermal Hall conductivity for various magnetic
ﬁelds and temperatures up to N = 72 spins(see Fig. 1) with a typical bond-dimension
D = 500. Our calculation showed that thethermal Hall conductivity typically has a peak
as a function of the temperature, and it overshoots the half quantized value expected in thezero-temperature limit. Such an overshooting
is consistent with the experimental observation[3].
We also investigated Γ and Γ′ dependenceof the thermal Hall conductivity. We foundthat depending on the sign of Γ and Γ′ , thethermal Hall conductivity is largely modiﬁed
from that of the pure Kitaev model. Some ofthem could be explained by the increase of thegap, previously calculated through the perturbation theory [5], while overall ﬁnite temperature behaviors might not be explained only by
Majorana fermion pictures.
A part of this work has been done in collaboration with J. Nasu, T. Misawa, and Y. Motome.

[5] D.
Takikawa
andS.
Fujimoto,
Phys. Rev. B 102 (2019) 174414.

References[1] A. Kitaev, Ann. Phys. 321 (2006) 2.
[2] Y.
MotomeandJ.
Nasu,
J. Phys. Soc. Jpn. 89 (2020) 012002.
[3] Y. Kasahara, T. Ohnishi, Y. Mizukami, etal., Nature 559 (2018) 227.
[4] H. Li, D.-W. Qu, H.-K. Zhang, et al.
Phys. Rev. Research 2 (2020) 043015.

211

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

Membrane shape deformation by binding ofcurvature-inducing proteins
Hiroshi NOGUCHI
Institute for Solid State Physics, University of Tokyo

In living cells, biomembrane shapes are controlled by curvature-inducing proteins. We
have studied the binding of these proteins using mean-field theories and compared their results with simulation results. Clathrin andcoat protein complexes bend membranes in alaterally isotropic manner and generate spherical buds. In contrast, Bin/Amphiphysin/Rvs(BAR) superfamily proteins bend the membrane anisotropically and generate cylindricalmembrane tubes. The former and latter proteins can be modeled as laterally isotropic objects and anisotropic objects of a banana or
crescent shape, respectively.
For isotropic proteins, we derived the freeenergy including the Gaussian curvature andarea expansion by the protein insertion [1].
Then, we investigated the formation of spherical buds [1] and binding onto a tethered vesicle [2]. A narrow membrane tube (tether) is
elongated from a spherical vesicle by an external force imposed by a micropipette and optical tweezers. A first-order transition occurs between a small number of large buds and a large
number of small buds with increasing chemical
potential of protein binding [1]. For the tethered vesicle, interestingly, a first-order transition occurs twice between low and high protein
densities in the tube. The force-dependence
curves of the protein density in the membrane tube and the tube curvature are reflection symmetric and point symmetric, respectively, from the force point, in which the tube
curvature matches the protein (sensing) curvature [2]. This theory reproduces the meshless
membrane simulation results of the homogeneous phases very well. In addition, beadednecklace-like tubes with microphase separation
are formed in the simulation.

212

For anisotropic proteins, an orientationaldependent excluded volume interaction is included in the free energy. We studied binding onto membrane tubes [3]. The proteins
exhibit a second-order or first-order nematictransition with increasing protein density forintermediate and small radii of the membrane
tube, respectively. The tube curvatures for the
maximum protein binding (sensing) and orientational order are differently varied by theprotein density and rigidity. As the external
force along with the tube axis increases, a firstorder transition from a large tube radius withlow protein density to a small radius with highdensity occurs once, and subsequently, the protein orientation tilts to the tube-axis direction.
The density-dependent sensing curvature andthe number of transitions are characteristics
due to the anisotropic bending energy. Eachisotropic protein has a constant sensing curvature. This density-dependent sensing curvature has been measured in the experimentof BAR proteins. This theory quantitatively
reproduces the results of meshless membrane
simulation for short proteins. For long proteins, the formation of protein clusters generates a quantitative deviation from the theory,
although qualitative agreements still hold.

References[1] H. Noguchi, Phys. Rev. E 104, 014410(2021).
[2] H. Noguchi, Soft Matter 17, 10469 (2021).
[3] H.
Noguchi,C.
Tozzi,
M. Arroyo,
Soft Matter
DOI:10.1039/D2SM00274D.

and(2022)

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

Development of Thermal Functional Materials Based on Firstprinciples Simulations
Junichiro SHIOMI
Department of Mechanical Engineering,
The University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo 113-8656, Japan
We have studied thermal functional

properties and the compositional dependencematerials using massive parallel computerof ternary alloys were studied by first-principle

system. This year, we focused on thecalculation. Fig. 1 shows the predicted energy

thermoelectric performance of a ternary alloy

by cluster expansion and ZT of alloy of asystem and a nanostructured multiphasecomposition. For the multiphase material, thematerial.

grain size dependence of thermoelectric

The main content of these two studies isproperties was studied by first-principle

briefly introduced as below. For the ternary

calculation for two phases to enhancealloys, lattice structures were predicted bythermoelectric performance by nanostructure

Monte Carlo simulation, then thermoelectricengineering. Fig. 2 shows the nanostructured

Fig. 1: (a) Formation energies from cluster

Fig. 2: (a) Thermal conductivity versus grain

expansion; (b) ZT of one of the alloys.

size; (b) ZT of a nanostructured material.

213

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

thermal conductivity at room temperature andPhys.-Condens. Mat. 26 (2014) 225402.

ZT of a nanostructured phase.

[4] G. Kresse and J. Furthmüller: Comp. Mater.
Sci. 6 (1996) 15.

To conduct these studies, the usage ofsupercomputer system was described asfollowing. ICET package [1] was employed to
predict lattice structures through performing
Monte Carlo simulation and cluster expansion
based on element substitution and formationenergy prediction. BoltzTrap2 package [2] wascarried out to calculate electrical transport
properties by solving semiclassical linearized
Boltzmann transport equation. ALAMODE
package [3] was used to fit interatomic force
constants by atomic displacements and forcesfrom first-principle calculation, as well as
calculate phonon thermal conductivity byrelaxation time approximation method. VASP
[4] was carried out to perform first-principle
calculation

forformationenergy,band

structure and atomic forces. By running ourown python script, thermoelectric properties
were

calculated

fromphononthermalconductivity and electrical properties.

References[1] M. Ångqvist, W. A. Muñoz, J. M. Rahm, E.
Fransson, C. Durniak, P. Rozyczko, T. H. Rod,and P. Erhart: Adv. Theor. Simul. 2 (2019)
1900015.
[2] G. K. H. Madsen, J. Carrete, and M. J.
Verstraete: Comput. Phys. Commun. 231(2018) 140.
[3] T. Tadano, Y. Gohda, and S. Tsuneyuki: J.

214

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

Feasibility study of Kitaev quantum spin liquid byfunctional renormalization approach
Kiyu FUKUI and Yusuke KATO
Department of Applied Physics, University of Tokyo, Bunkyo, Tokyo 113-8656
Department of Basic Science, University of Tokyo, Meguro, Tokyo 153-8902

The Kitaev model is an S = 1/2 quantum spin model defined on the honeycomb lattice [1]. In spite of its strong frustration due
to the bond-dependent interaction, this modelis exactly solvable, and its ground state is aquantum spin liquid state. Here we address
the feasibility of the Kitaev quantum spin liquid (KQSL) for two extensions of the model,
by using the pseudofermion functional renormalization group (PFFRG) method [2].
The first one is for ultracold polar molecules
trapped in an optical lattice.
While experimentally accessible implementation of
the Kitaev-type interaction was proposed in2013 [3], the stability of the KQSL has not been
investigated so far. In this proposal, the bonddependent anisotropic interactions are mimicked by angle-dependent dipole interactions
between molecules. However, the interactions
are long ranged with a spatial decay of r−3 ,where r is the distance between the molecules,
it is left as an open question whether the
KQSL can survive against such long-range interactions. We studied the ground state ofa quantum spin model with long-range angledependent Kitaev-type interactions, which was
proposed as an implementation of the Kitaev
model in ultracold polar molecules, by using
the PFFRG method. We clarified that, regardless of the spatial anisotropy of the interactions, the ground state is magnetically orderedin both FM and AFM cases: we found magnetic instabilities toward the FM and zigzag
ordered states in the FM and AFM models, re-215

spectively. By calculation of the anisotropy parameter dependence of the critical cutoﬀ scale,
we concluded that the system is most frustrated and closest to the realization of the
KQSL when the interaction is isotropic in bothcases. Our findings indicate that the QSL
ground state arising from the nearest-neighborbond-dependent anisotropic interactions in the
Kitaev model is destroyed by the long-range interactions. By varying the range of the interactions, we unraveled that the KQSL is unstableeven for the third-neighbor interactions. Thisis the first time that the stability of the KQSL
within the above proposal has been revealed.
The other is the extension to higher-spin systems. Although the Kitaev model was originally introduced for the S = 1/2 moments, its
higher-spin generalization has also attracted
attention. While candidate materials with S =1 [4] and 3/2 [5] have begun to be proposed inrecent years, the stability of the KQSL has notbeen systematically clarified for general spin
S in the presence of non-Kitaev interactions,which are inevitably present in the candidate
materials. In this study, we present our numerical results on the ground state of higher-spin
generalization of the Kitaev model including
the Heisenberg interaction, the spin-S KitaevHeisenberg model, which we consider to be oneof the minimal models for the higher-spin candidate materials, by using an extension of thePFFRG method [6]. Performing the calculations for the models with S = 1, 3/2, 2, 5/2,
and 50 in addition to S = 1/2, we elucidate

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

the ground-state phase diagram by systematically changing the ratio between the Kitaev
and Heisenberg interactions and S. We found
four magnetically ordered phases: Néel AFM,
zigzag AFM, FM, and stripy AFM for all S,
and the phase boundaries between these ordered phases are consistent with the previous
studies for S = 1/2, 1, and ∞. Our results
indicated that the ground state is short-range
QSL for arbitrary S at the ferromagnetic andantiferromagnetic Kitaev points. It is consistent with an analytical result in the previous
study [7]. We clarified that the QSL regions
near these Kitaev points are quickly shrunk
by increasing S, namely, the KQSL is fragile
against the Heisenberg interaction for S ≥ 2.
However, we find the KQSL remains stable forS ≤ 3/2 around the both Kitaev points. Ourfindings will provide the search and design ofthe candidate materials with S > 1/2.

References[1] A. Kitaev, Ann. Phys. 321, 2 (2006).
[2] J. Reuther and P. Wölfle, Phys. Rev. B
81, 144410 (2010).
[3] S. R. Manmana et al., Phys. Rev. B 87,
081106 (2013); A. V. Gorshkov, K. R.
Hazzard, and A. M. Rey, Mol. Phys. 111,
1908 (2013).
[4] P. P. Stavropoulos, D. Pereira, and H.-Y.
Kee, Phys. Rev. Lett. 123, 037203 (2019);
A. M. Samarakoon et al., Phys. Rev. B
104, 184415 (2021).
[5] C. Xu et al., npj Comput. Mater. 4, 57(2018); M. Kim et al., Nat. Electron. 2,457 (2019); I. Lee et al., Phys. Rev. Lett.
124, 017201 (2020); C. Xu et al., ibidem,
087205 (2020).
[6] M. L. Baez and J. Reuther, Phys. Rev. B
96, 045144 (2017).

216[7] G. Baskaran, D. Sen, and R. Shankar,
Phys. Rev. B 78, 115116 (2008).

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

Quantum effects on chiral magnets
SOHEI KODAMA and YUSUKE KATO
Department of Basic Science, the University of Tokyo, Tokyo 153-8902, Japan,

We consider the ground state of quantum

complementary to earlier studies on quantum

spin chain of monoaxial chiral ferromagnet.

effects on Ising domain wall [2] and skyrmions

The Hamiltonian is given by[3] based on semiclassical treatment of̂𝑖 ∙ ̂
̂𝑖 × ̂
𝐻 = ∑𝑖 [−𝐽𝑺
𝑺𝑖+1 − 𝐷(𝑺𝑺𝑖+1 )𝑦 −quantum spin taking account of the Berry phase2𝐻𝑆̂𝑖𝑧 + 𝐾(𝑆̂𝑖𝑦 ) ],term.

where, J, D, H, and K, respectively, denoteexchange interaction, Dzyaloshinskii-Moriya(DM)interaction,themagnetic

field

perpendicular to the chiral axis (y -axis) andsingle ion anisotropy.
The magnetization curves for finite-sized
systems obtained by exact diagonalization are
distinctly different between half-odd-integer
and integer spins; The level crossings occur for
S=1/2,3/2, accompanied by  shift of the crystal
momentum k while the magnetization process
is continuous and k is always zero for S=1,2.
These behaviors are commonly observed for

References

arbitrary values of J/D. We argue that in the[1] S. Kodama, A. Tanaka, and Y. Kato, in

limit J=0 with D finite, single soliton has the

preparation.

minimum energy for k= (k=0) for half-odd[2] H. B. Braun, Adv. Phys. 61, 1 (2012).

integer (integer) S, based on the Perron-[3] R. Takashima, H. Ishizuka, and L. Balents,

Frobenius theorem within the subspace spanned

Phys. Rev. B 94, 134415 (2016).

by a single soliton basis. Our approach [1] is217

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

New scaling analysis for graph percolations
Y. Ozeki and K. Hagiwara v
Graduate School of Informatics and Engineering, The University of Electro-Communications

We investigate the transition point and the critical exponent β for graph percolations. A new sizeindependent scaling analysis is proposed, which enables us to analyze numerical data without finitesize effects. This means that the results are equivalent to those in the infinite system. Such a scaling analysis has been applied as dynamical scaling
analysis in the nonequilibrium relaxation method[1] for phase transitions and critical phenomena.
In graph percolations, a number of vertices areprepared instead of a lattice. Each edge between
two of these vertices is considered to be connected
or not. In the Erdős-Rényi (ER) model, a classical graph percolation, a new edge is connected to agraph at random. The probability of a vertex belonging to the infinite cluster, which is denoted by
P (r), is recognized as the order parameter, wherer is the average number of edges per vertex. The
function P (r) can be written as [2]X
P (r) ≡ 1 −
sns ,
(1)
s

where ns is the number of clusters with size s
divided by N . P (r) is expected to behave algebraically as
P (r) = a|r − rc |β ,(2)where β is referred to as a critical exponent. In thepresent study, we define a new function P (x, r) byextending P (r),and introduce a new parameter x
representing the maximum cluster size:X
P (x, r) ≡ 1 −
sns .
(3)
s<x

This function P (x, r) represents the probability
that a randomly selected vertex belongs to a cluster
greater than or equal to x. We expect the following
asymptotic (x → ∞) behaviors for every r case:(r < rc ) A exp(−x/ξ)P (x, r) ∼x−λP
(r = rc ) (4)
p∞ + A′ exp(−x/ξ) (r > rc )where we introduce a function ξ depending on ras ξ = ξ(r) and exponent λP . The asymptotic
behavior of ξ is expected to beξ(t) = b|r − rc |−ν ,

(5)

218

where we define a new exponent ν satisfying therelation
λP = β/ν.

(6)

In the asymptotic regime in r ∼ rc , P (x, r) is
expected to satisfy the following scaling form:
P (x, r) = ξ −λP Ψ(x/ξ).

(7)

The transition point rc and critical exponent β canbe estimated from the behavior of P (x, t). Appropriate values of rc and β as well as a, b, ν in Eqs.
(2) and (5) should be chosen such that ξ λP P (x, t)
is fitted as a function of x/ξ(r) on a scaling function Ψ. In order to carry out such a scaling plot
efficiently, we apply the method introduced for the
dynamical scaling, in which the Bayesian inference
and the kernel method are used [3].
Here, we demonstrate the above scaling analysis
to the ER model, for which the analytic solution is
known; rc = 1/2, β = 1. For precise analysis, weneed to check the size dependence. We have foundthat, for x ≤ 3000, no size dependence appears in
the data for the system with N = 1 × 108 or larger.
The simulation is performed using these parameters, and 2240 independent samples were used for
averaging. The results for P (x, r) are plotted inFig. 3, and the resultant scaling plot is shownin Fig. 4. The transition points and the critical exponent are evaluated as rc = 0.49978(7) and
β = 1.010(9), respectively, where the error bars are
estimated by the bootstrap method [3]. The deviations from the analytic values rc = 1/2, β = 1
are approximately 0.044% and 1.0%, respectively,
demonstrating that the present analysis is reliable.
We apply the above scaling method to thedCDGM model [4], in which a vertex is selectedfrom randomly chosen m vertices so that the size ofthe connected cluster is smallest and the new edge
connects two vertices selected twice with the aforementioned process; the integer value m is fixed anddistinguishes the model as the dCDGM(m) model.
We investigate the dCDGM(m = 2) model. The
same analysis was performed. It is confirmed thatthe simulation for the system with N = 1 × 108
to x ≤ 3000 is sufficiently large. The simulation is
performed using these parameters, and 4480 inde-

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

0.5
0.4
0.3
0.2
0.1

P(x,r)

P(x,r)

r = 0.480
r = 0.482
r = 0.484
r = 0.486

r = 0.9197

r = 0.488
r = 0.490

r = 0.9202

r = 0.492
r = 0.494

r = 0.9212
r = 0.9217

r = 0.496
r = 0.498

r = 0.9222

r = 0.9207

r = 0.9227

100

1000

1000

x

x

Figure 3: Results of P (x, r) for the dCDGM(2)model.

Figure 1: Results of P (x, r) for the ER model.
100

10

10

Y = ξλPP(x,r)

Y = ξλPP(x,r)

1

1

0.1

0.1

0.01

0.01
0.001

0.01

0.1

0.001

1

1

X = x/ξ

X = x/ξ

Figure 2: Scaling plot of P (x, r) for the ER model.
pendent samples are used for averaging. The transition points and the critical exponent are evaluated as rc = 0.923211(5) and β = 0.056(1), respectively, which are consistent with those of Da Costa(rc = 0.923207509297(2) and β = 0.05557108(1))[4].

References[1] Y. Ozeki and N. Ito, J. Phys.: Math. Theor.
40 R149 (2007).
[2] D. Stauffer and A. Aharony, “Introduction topercolation theory”, (Taylor & Francis,1994)[3] Y. Echinaka and Y. Ozeki, Phys. Rev. E94
043312 (2016).

219

Figure 4: Scaling plot of P (x, r) for the dCDGM(2)model.
[4] R. A. Da Costa, S. N.Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Phys. Rev. E89
042148 (2014); Phys. Rev. E90 022145 (2014).

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

Computational design of proteins useful for medical andindustrial applications
Masataka YOSHIMURA1, Sairi MATSUMOTO1, Mizuki TERANISHI1,
Hirotaro SHIMAMURA2, Rina AOYAMA1, Koji OOKA2, Shino MASUOKA1,
Shunji SUETAKA1, Nao SATO1, Yuuki HAYASHI1, and Munehito ARAI1,2
1

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

Department of Physics, Graduate School of Science,

The University of Tokyo, Komaba, Meguro, Tokyo 153-8902

ProteinsasSevere acute respiratory syndrome (SARS)

therapeutics for various diseases and as enzymes

coronavirus 2 (SARS-CoV-2) is an emerging

for industrial production of useful substances.

virus that causes the pandemic of the

Recently, theoretical design of proteins has

coronavirus disease 2019 (COVID-19). Since

attracted much attention; however, it is still

the viral infection is caused by the binding of the

difficult to theoretically design useful proteins.

receptor binding domain (RBD) of the spike

To solve this problem, we have designed various

proteinuseful proteins using the Rosetta 3.13 suite [1]angiotensin-converting enzyme 2 (ACE2) [3],and have experimentally checked the validity of

we computationally designed a single-chain

the computational design.

variable fragment (scFv) of an antibody against

Magnets

have

are

been

widelyoptogenetic

used

tools

of

SARS-CoV-2

to

the

human

for

the SARS-CoV-2 RBD using the scFv against

manipulating protein-protein interactions (PPIs)

the SARS-CoV RBD as a template. Experiments[2]. They consist of pMag and nMag, whichshowed that the designed antibody has a highheterodimerize upon exposure to blue light.

affinity to SARS-CoV-2 RBD comparable to

However, improvement of their binding affinitythat of the antibody medicines for SARS-CoV-2.
Formation

has been demanded. Here we computationallyofaternary

complexbydesigned the Magnets mutants with higher

interleukin-33 (IL-33), interleukin-1 receptor

pMag/nMag affinity than the original ones byaccessary protein (IL-1RAcP), and the ST2

theoreticalmutagenesis

Rosetta.

receptor of type 2 innate lymphoid cells isExperiments

confirmed

mutants

responsible for allergic asthma [4]. To develop aenhanced the dimer affinity. Thus, these

therapeutic agent for the disease, we are trying

Magnets mutants may be useful for in vivo

to computationally design a protein that tightly

manipulation of PPIs.

binds ST2 but does not bind IL-1RAcP using IL-

using
that

the

220

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

33 as a template. We have finished designing aand product-bound form of ADO. We found that

protein that tightly binds to ST2. Further design

the product has higher flexibility than theof the IL-33 mutant that abolishes the binding

substrate inside the protein molecule, thereby

with IL-1RAcP is under way.

facilitating the product release.

The interaction between the programmed cell

Finally, to predict the folding processes ofdeath protein 1 (PD-1) and the programmed

proteins, we calculated free energy landscapes

death-ligandinof multi-domain proteins using the Wako-Saitô-suppressing T-cell inflammatory activity [5]. To

Muñoz-Eaton model, which is an Ising-like,develop an inhibitor of the T-cell activity, we

structure-based statistical mechanical model [8].

computationally designed the mutants of a PD-

By introducing modifications to the model, we

L1 fragment that tightly bind PD-1. Experiments

succeeded in obtaining the folding pathways andconfirmed that some mutants have improved

intermediates that are fully consistent withaffinity with PD-1. We are trying to further

experimental results. Thus, our method may

improve the affinity by multiple mutations.

pave

1

(PD-L1)isimportant

Alkane biosynthesis has attracted much

the

way

for

solvingthefoldingmechanisms of multi-domain proteins that are

attention as a promising way of producing

abundant in proteomes.

carbon-neutral bioenergy. A key enzyme forcyanobacterial alkane biosynthesis is aldehyde

References

deformylating oxygenase (ADO) [6]. We[1] J. K. Leman, et al.: Nat Methods, 17 (2020)

computationally designed ~10 mutants of ADO

665-680.

that are expected to have high activity. The[2] F. Kawano, et al.: Nat Commun. 6 (2015)

experimental test of the ADO activity is ongoing.

6256.

Understanding protein dynamics is essential[3] J. Shang, et al.: Nature, 581 (2020) 221-224.

for improving catalytic activity of enzymes. We[4] X. Liu, et al.: Proc Natl Acad Sci USA. 110

performed molecular dynamics simulations of(2013) 14918-14923.

ADO using the Gromacs 2021.4 program[5] C. N. Bryan, et al.: Proc Natl Acad Sci USA.

package [7]. An initial structure was constructed

118 (2021) e2102164118.

using the AmberTools21 package. The force[6] M. Arai, et al.: Adv Exp Med Biol. 1080

fields for the substrate (aldehyde) and the(2018) 119-154.

product[7] M. J. Abraham, et al.: SoftwareX, 1-2 (2015)(alkane)were

obtainedusing

Antechamber module. The simulations were

19-25.

performed for ~800 ns each for the apo form,[8] M. Sasai, et al.: Biophys Physicobiol. 13

holo (iron-bound) form, substrate-bound form,(2016) 281-293.

221

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

Topological order and quantum operation inquantum many-body systems
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 methods,
tensor-network algorithms, and other optimization techniques and investigated topological quantum phase transitions in strongly correlated many-body systems and optimizationof quantum operations in quantum computing.
Phase transition in quantum dimer model:
Based on the idea of stochastic series expansion, we implemented a Monte Carlo algorithm for investigating the phase transition of the quantum dimer model in an efﬁcient way to simulate with shorter CPU
time and extended it to simulation on larger
Hilbert space, including monomer-dimer andmonomer-monomer interaction. We obtained
the phase-transition temperature between thecritical and columnar phases by diﬀerent order parameters by simulations at ﬁnite temperatures. We found it vanishes at around
Rohksar-Kivelson (RK) point, while the classical limit merges smoothly to the value obtained by the classical Monte Carlo.
Tensor-network Markov-chain Monte Carlo:
We developed a novel quantum-inspired
method that combines the tensor-network representation of weight function and the sampling scheme based on the Markov-chain
Monte Carlo method. We can remove thesystematic bias due to the ﬁniteness of thebond dimension in tensor network representation from the Markov-chain sampling. Furthermore, the statistical error of the sampling
reduces exponentially as increasing bond di-mension, and we can avoid the negative-sign
problem in quantum lattice models by adopting a high-precision tensor network renormalization scheme.
Many-body localization in random spin chain:
We study the many-body localization of therandom-ﬁeld Heisenberg chain using the nested
shift-invert Lanczos method with an iterative
linear solver. We use the minimum residual
method (MINRES) inside each Lanczos iteration. As a probe of many-body localization transition, we propose a unitary operator
called the twist operator, which has a straightforward interpretation in the real space. We
estimate the transition point of the randomﬁeld Heisenberg chain precisely [1].
Decomposition of multi-qubit gate: We obtaineﬃcient decompositions of CCZ and CCCZ
gates, typical multi-qubit gates, under several qubit connectivities, using the sequential
quantum-circuit optimization algorithm. We
can construct the CCZ gate with only four
CZ-depth when the qubit is square-shaped,
including one auxiliary qubit. In T-shaped
qubit connectivity, which has no closed loop,
we can decompose the CCCZ gate with 17 CZ
gates [2].

References[1] T. Kutsuzawa and S. Todo, preprint:
arXiv:2203.09732.
[2] K. M. Nakanishi, T. Satoh and S. Todo,
preprint: arXiv:2109.13223.

222

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

O(N) algorithm for the Monte Carlo simulation of
the Kitaev model
Masahiko G. YAMADA1
Yutaka AKAGI2
Masafumi UDAGAWA3
1
Department of Materials Engineering Science, Osaka University
Machikaneyama-cho, Toyonaka, Osaka 560-8531
2
Department of Physics, University of Tokyo
Hongo, Bunkyo-ku, Tokyo 113-0033
3
Department of Physics, Gakushuin University
Mejiro, Toshima-ku, Tokyo 171-8588

The thermodynamics of the Kitaev model isof great interest, especially after the characteristic specific heat behavior called a double
peak has been discovered by the Monte Carlosimulation [1].
We extend their simulation using a newO(N ) calculation method for the Kitaev
model Monte Carlo [2]. We newly employed
the infinite product expansion (iPE) method
instead of the previous Trotter decomposition
method. By this iPE method, we succeededin reducing the calculation cost of the Kitaev
model drastically from O(N 4 ) to O(N ). Asa benchmark, we did a simulation for the
Kitaev model on the honeycomb lattice andcomputed the heat capacity and its doublepeak structure. Using the new algorithm anda supermassive parallelization in the ISSP
supercomputer, we succeeded in reproducing
the results in [1] within the errorbar (seeFig. 1). Though currently the calculation
is limited up to the 54-site cluster with anopen boundary, we will eventually extend
this to the large-scale Monte Carlo simulation
up to O(1000)-site simulations. If we also
extend the calculation to three dimensions,
we can potentially discuss the thermodynamic
properties of the vaporization transition in

the three-dimensional Kitaev model, which
was suggested in Ref. [1].

Figure 1: Specific heat data for the N = 54
cluster simulation.

References[1] J. Nasu, M. Udagawa, and Y. Motome:
Phys. Rev. Lett. 113 (2014) 197205.
[2] M. G. Yamada, Y. Akagi, and M. Udagawa: to appear.

223

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of TokyoDMRG algorithm with a full SU(N) implementation
Masahiko G. YAMADA
Department of Materials Engineering Science, Osaka University
Machikaneyama-cho, Toyonaka, Osaka 560-8531

The density matrix renormalization group(DMRG) is one of the most powerful numerical methods for strongly correlated condensed
matter systems. Following the previous work,
we extended DMRG to the case with the
SU(N ) symmetry with N > 2, even for twodimensional systems. As a killer application,
we simulated the ground state of the SU(4)Heisenberg model on the honeycomb lattice,which can potentially be realized in cold
atomic systems and solid state systems like
α-ZrCl3 [2]. We keep up to 12800 SU(4) statesequivalent to more than a million U(1) states.
This supermassive DMRG simulation reveals
the quantum spin-orbital liquid ground state,
which has dubiously been expected for morethan a decade. The methodology developed
here can be extended to any classical Liegroups, paving the way to a next-generation
DMRG with a full symmetry implementation [3].

References[1] P. Nataf, and F. Mila: Phys. Rev. B 97(2018) 134420.
[2] M. G. Yamada, M. Oshikawa, and G.
Jackeli: Phys. Rev. Lett. 121 (2018)
097201.
[3] M. G. Yamada, K. Penc, and F. Pollmann: to appear.

224

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

Efficient Sampling Simulation of
the Soft Modes Significantly Contribute to Protein Properties
K. TAKEMURA, D. TRAN, M.M. SOBEH, Y. MIYAZAWA, R. KIUCHI, T.N.
WIJAYA, T. OGAWA, Y. TAIRA, R. MISU, R. Osawa, M. Yoshioka, H. NAKAYA and
A. KITAO
School of Life Science and Technology, Tokyo Institute of Technology
Ookayama, Meguro, Tokyo 152-8550
Previously, we developed parallel cascadetransitions from the bound to unboundselection molecular dynamics (PaCS-MD)

states. Finally, we repeated the cycles until

as an enhanced sampling strategy that uses

d = 70 Å to complete the dissociation, as

cycleswe can see in the PaCS-MD flowchart inofmultiple

independentMD

Fig. 1

simulations done in parallel to improve thechance of biologically rare events occurring
without external perturbation [1]. PaCSMD

consistsofcyclesofseveralindependent parallel short all-atom MD
simulations followed by selecting initialstructures for the next cycle according to a
specified quantity. We conducted thedissociation simulation of the p53 DNA
binding domain (p53DNA) from DNA [2].
We used the inter-center of mass distance(inter-COM distance, d) between p53-DBD

Fig. 1 The flowchart of PaCS-MD

and DNA as a quantity for ranking thesnapshots generated in each cycle. By

We successfully generated 75 different

repeating a series of cycles for the selecteddissociation pathways of p53-DBD from

top-ranked snapshots (top 10 snapshots),DNA. We used the COM positions of the

dissociation PaCS-MD generates structures

p53-DBD interface residues relative to thosewith larger inter-COM distances than those

of

found

which

pathways, as shown in Fig. 2. Inspection ofsignificantly enhances the probability ofthe dissociation pathways indicates that the

in

the

previous

cycle,

225

DNA tovisualize

thedissociation

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

sampled space formed a cone-like shapebinding free energy (ΔG°) for this system

around the DNA. Around 93% of thepredicted at 11.8 kcal/ mol. This value is very

sampling pathways (70 pathways) dissociatedclose to the binding free energy of p53-DBD

along the +X and −Y directions (namely, the(residues 94−312) with the consensus DNA−Y directions), while the other pathways

sequence of −11.1 kcal/ mol measured by

moved along +X and +Y directions, whichisothermal titration calorimetry, suggesting

only occupy 7% (5 pathways), hereafterthat the calculated free energy profile iscalled the +Y directions.

reasonable.

Fig. 2 Dissociation pathways of 75 PaCS-

Fig. 3 Free energy profile mapped onto the

MD trials of p53-DBD from DNA.

XY-plane

Using the Markov state model (MSM), We
obtained the free energy profile of p53-DBD
dissociation from DNA and showed the
projections onto the XY planes, as shown in
Fig. 3. The global minimum of this profile
represents the bound state. As discussed
above, dissociation mainly (93%) occurredalong the −Y directions on the XY-plane, andFig.4 PMF as a function of X coordinate.

other dissociations (7%) occurred along the+Y directions. As we can see in Figure 4,

References

during the dissociation process, the PMF[1] R. Harada and A. Kitao. J Chem Phys 139,

increased at a higher rate in the bound state to

035103 (2013).

achieve convergence (became flat) in the[2] M.M. Sobeh and A. Kitao J. Chem. Inf.

unbound region (d ≥ 45 Å). We calculated theModel, 62, 1294–1307 (2022).

226

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

Molecular dynamics simulations for assembly anddisassembly of protein aggregates
Hisashi Okumura
Exploratory Research Center on Life and Living Systems,
Institute for Molecular Science, Okazaki, Aichi 444-8585

Proteins are normally folded correctly andthe RPM and RSPM into LAMMPS [1].

perform functions that are necessary to sustain

LAMMPS is an open-source classical MD

life. However, when the concentration ofsimulation package that can be applied fromproteins increases due to aging or other reasons,material to biomolecular systems. We used thethey can aggregate and cause a variety ofsupercomputer at the Supercomputer Center,diseases. For example, Alzheimer’s disease,

the Institute for Solid State Physics, thewhich is one type of dementia, is caused by

University of Tokyo to perform benchmark

amyloid-β (Aβ) peptides, which aggregate into

calculations for this program.

spherical

oligomersoramyloidfibrils.

In addition, continuing from the last year,

Molecular dynamics (MD) is nowadays awe performed non-equilibrium MD simulations

standard tool to investigate these systems.

of the destruction process of amyloid fibrils by

However, once a conformation gets trapped in aan infrared laser. We first revealed the

local-minimum free-energy state, it is hard to

disruption process of Aβ amyloid fibrils [2]. We

escape from there within the time scale of theare now performing MD simulations of poly-

simulation.

alanine amyloid fibrils. We found that poly-

The replica-exchange method (REM), isalanine amyloid fibril is less easily destroyed

widely used to overcome this difficulty because

by infrared laser than Aβ amyloid fibrils. Weof its convenience. Recently, we proposed thewill try to clarify this difference.

replica-permutation method (RPM) and thereplica sub-permutation method (RSPM) to

References

improve the transition ratio of replicas in the[1] M. Yamauchi, G. La Penna, S. G. Itoh, andtemperature space. However, the RPM andH. Okumura, Comput. Phys. Commun. 276,

RSPM are not available for general users

108362 (2022).

because these methods are implemented only in[2] H. Okumura, S. G. Itoh, K. Nakamura, andour own molecular simulation programs and

T. Kawasaki, J. Phys. Chem. B 125, 4964–4976

not in open-source packages. We implemented(2021).

227

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

Critical phenomena in novel Anderson transitions
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 have no translational/rotational symmetries, but can havetime reversal, spin rotation, chiral andparticle-hole symmetries. According to thesesymmetries, the random quantum systems areclassiﬁed to 10 symmetry classes. The threeof them, the so-called Wigner-Dyson classes,
were discovered in 1950’s and 1960’s and havebeen well studied. Recent experimental discoveries of topological insulators and Weyl
semimetal (WSM) have inspired extensive research of random systems other than WignerDyson classes. These 7 new classes containdiscrete chiral or particle-hole symmetries.
Here we have introduced a set of Hamiltonians for these new classes, where the high accuracy numerical studies of the three dimensional Anderson transition is possible [1]. We
have also studied the novel two dimensional
systems, corresponding to a topological superconductor, both analytically [2] and numerically [3]. We especially focused on the criticalbehaviors near the tricritical point, where ordinary insulator, topological insulator and thermal metal phases meet.
When we introduce the non-Hermiticity
to the Hamiltonian, the classiﬁcation is extended and the systems are categorized to 38
classes. We have shown that adding the nonHermiticity changes the universality class [4],and show that the transfer matrix method
widely used in the Anderson transitions also
works for non-Hermitian systems [5]. We havefurther shown that the critical behaviors of 38
classes are mapped to universality classes forthe 10 Hermitian symmetry classes [6].
Machine learning is now widely used in con-densed matter physics. We have studied thevarious quantum phase transitions by analyzing the wave functions via convolutional neural network (CNN). The method is adopted
to classify wave functions obtained via density
functional theory, which is applied to study themetal-insulator transitions in doped semiconductors [7]. In addition, using the long-shortterm memory network (LSTM), we have analyzed the time series of diﬀusion motion ofquantum particles in kicked rotors. The model
is equivalent to the Anderson model, where thedimensions of the Anderson transition is controlled by the modulation of the kick strength.
We have drawn the phase diagram of higher dimensional Anderson transitions by LSTM [8].

References
1. T. Wang, T. Ohtsuki, and R. Shindou, PhysicalReview B 104, 014206 (2021).
2. Z. Pan, T. Wang, T. Ohtsuki, and R. Shindou,,
Physical Review B 104, 174205 (2021).
3. T. Wang, Z. Pan, T. Ohtsuki, I.A. Gruzberg,
and R Shindou, Physical Review B 104, 184201(2021).
4. X. Luo, T. Ohtsuki, R. Shindou, Physical Review Letters 126, 090402 (2021).
5. X. Luo, T. Ohtsuki, and R. Shindou, PhysicalReview B 104, 104203 (2021).
6. X. Luo, Z. Xiao, K. Kawabata, T. Ohtsuki, and
R. Shindou, arXiv:2105.02514
7. Y. Harashima, T. Mano, K. Slevin, and T. Ohtsuki, Journal of the Physical Society of Japan 90,
094001 (2021).
8. T. Mano, and T. Ohtsuki, Annals of Physics
435, 168500 (2021).

228

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

Non-equilibrium relaxation analysis on the J1 − J2
frustrated Ising model
Keisuke Iikura and Hiroshi WATANABE
Department of Applied Physics and Physico-Informatics, Keio University
Yokohama, Kanagawa 223-8522, Japan

We consider frustrated J1 -J2 Ising model on
the square lattice which Hamiltonian is givenby∑∑
H = J1σi σj + J 2
σi σj ,
⟨i,j⟩⟨⟨i,j⟩⟩

where ⟨i, j⟩ and ⟨⟨i, j⟩⟩ denote the summation
over the nearest and the next-nearest neighborpairs, respectively. Spins are Ising-like variable, σi = ±1 and J1 < 0 and J2 > 0 denote ferromagnetic and antiferromagnetic interactions, respectively. The parameter of thesystem is given by g = J2 /|J1 |. This model
exhibits criticality with varying exponents as
the parameter g changes. This behavior isbelieved to be Ashkin-Teller (AT) universality class. The AT-like behavior will end atsome points g ∗ , and the phase transition becomes the ﬁrst-order for 1/2 < g < g ∗ . The
end-point is predicted to be 4-state Potts universality class. However, the localtion of theend-point g ∗ varies from study to study. Furthermore, whether the phase transition in therange 1/2 < g < g ∗ is a ﬁrst-order transition
or a continuous transition is under debate [1].
Since the interactions of this model involves
frustration, the cluster algorithm cannot beapplied. Therefore, we apply non-equilibrium
relaxation (NER) method to analize the criticality of this model.
We ﬁrst observe the non-equilibrium relaxation of the order parameter. The stripe orderof this system is deﬁned as follows,
mx ≡

1 ∑
(−1)xi σi ,
N i

my ≡

where (xi , yi ) are the coordinates of the i-th
spin on the lattice. Then the magnetizaiton is
given by m2 = m2x + m2y . The order parameter
decays exponentially in the disordered phase,converges to a constant in the ordered state,
and exhibits power-law decay at the criticalpoint. Then we can identify the critical pointfrom the non-equilibrium behavior of the parameter. We determine the critical exponents
from the non-equilibrium ﬂuctuation of obervables [2]. The location of Potts point is determined to be g ∗ ∼ 0.58. From the numericalanalysis, the phase transition is continuous for
1/2 < g < g ∗ and we did not ﬁnd any signs of
the ﬁrst order transition. We are currently analyzing the behavior of ordered structure underthe ﬁelds, and we conﬁrm the existence of a nematic phase. The identiﬁcation of the natureof the phase transition is a subject for future
work.

References[1] S. Jin, A. Sen, and A. W. Sandvik, Phys.
Rev. Lett. 108, 045702 (2012).
[2] Y. Ozeki and N. Ito, J. Phys. A: Math.
Theor. 40, R149 (2007).

1 ∑(−1)yi σi ,
N i

229

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

Structure and mechanical properties of crystalline
polymers absorbing water molecules
Yuji HIGUCHI
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Research Institute for Information Technology, Kyushu University
Motooka Nishi-ku, Fukuoka 819-0395

We have studied soft matter by molecularsimulations and revealed the fracture processesof semicrystalline polymers [1] for the improvement of the toughness and dynamics of lipid
bilayers [2] as a model system for biologicalmembranes.
This year, self-assembly processes of amphiphilic molecules were studied by dissipative particle dynamics (DPD) [3] simulations
on systems B and C using the LAMMPS
program [4]. We revealed that two diﬀerent molecules were self-sorted or mixed andthe processes depended on the hydrophobic
length [5].
Furthermore, we studied the adsorption andabsorption processes of water molecules tocrystalline polymers by coarse-grained molecular dynamics simulation. Water molecules are
represented as a sphere in the simulation. The
attractive interaction between water and polymer spheres is suﬃciently strong to treat polymers as hydrophilic. Figure 1 shows a snapshot
at 10 ns and the number of adsorbed watermolecules to polymers along the cell axis. Water molecules selectively attach to and absorb
into amorphous layers. The tendency is thesame as that in the diﬀerent surfaces, where
the water sphere selectively absorbs into amorphous parts. In the hydrophilic semicrystalline
polymers, the absorption of water into amorphous parts leads to the gain of attractive interaction compared to that into crystal parts.

Figure 1: Adsorption and absorption processesof water molecules to polymers. Number ofadsorbed water molecules to polymers along
the cell axis.

References[1] Y. Higuchi: Phys. Rev. E 103 (2021)
042502.
[2] Y. Higuchi, Y. Asano, T. Kuwahara, and
M. Hishida: Langmuir 37 (2021) 5329.
[3] R. D. Groot and P. B. Warren, J. Chem.
Phys. 107 (1997) 4423.
[4] S. Plimpton:(1995) 1.

J. Comput. Phys. 117[5] R. Wakabayashi, R. Imatani, M. Katsuya,
Y. Higuchi, H. Noguchi, N. Kamiya, and
M. Goto: Chem. Commun. 58 (2022) 585.

230

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

Majorana-mediated spin transport in Kitaev modelAkihisa Koga
Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551

Spin transport has recently attracted muchinterest. One of the interesting examples is thespin transport in the Kitaev model [1]. In themodel, quantum spins are fractionalized intoitinerant Majorana fermions and local fluxes
due to quantum many-body eﬀects. Due to theexistence of the local conserved quantities, spin
moments never appear in the ground state.
Nevertheless, the itinerant Majorana fermionscarry the spin excitations and the spin transport is realized without spin oscillations. Since
Majorana and flux excitations have distinct energy scales, it is highly desired to clarify how
stable such a Majorana-mediated spin transport is against thermal fluctuations.

quantum state method to the 28-site cluster,
we examine the dynamics of the system after the magnetic pulse is introduced at theleft edges. When T ∼ TL , where TL corresponds to an energy scale of local flux, larger
oscillations in the spin moments are induced in
the other edge, compared to the results at the
ground state. At higher temperatures, excitedZ2 fluxes disturb the coherent motion of theitinerant Majorana fermions, which suppresses
the spin propagation. Our results demonstratea crucial role of thermal fluctuations in the
Majorana-mediated spin transport [2].
T/J=0.5

T/J=0.2

T/J=0.1
T/J=0∆Sz

T/J=0.01

T/J=0.05
T/J=0.03

tJ

Figure 2: Real-time evolution of the change inthe local magnetization in the right edge.
Figure 1: Kitaev model with armchair edges.
Green, red, and blue lines indicate x-, y-, andz-bonds, respectively. The static magnetic
field hR is applied in the right region, and nomagnetic field is applied in the middle region.
Time-dependent pulsed magnetic field is introduced in the left region.
Motivated by this, we deal with the Kitaev
model with armchair edges, which is shown in
Fig. 1, and consider the spin transport at finite temperatures. Applying the thermal pure

References[1] T. Minakawa, Y. Murakami, A. Koga, and
J. Nasu, Phys. Rev. Lett. 125, 047204(2020).
[2] H. Taguchi, Y. Murakami and A. Koga,
Phys. Rev. B 105, 125137 (2022).

231

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

Development of COMPutation ARchive of Exact
Diagonalization (COMPARED)
Kota IDO
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Development of accurate numerical solvers

highly efficient simulations on supercomputers.

for quantum many-body systems is one of the

Although this parallelization is efficiently done

central issues in the field of computational

for large system sizes such as 36 sites ofmaterials sciences. To achieve this goal, thequantum spin systems, it is not for single

exact diagonalization method has been widely

parameter set for small system sizes. To generateused for benchmarks on the systems smaller thandata for small system sizes, we develop a python

50 sites because this method enables us to obtain

tool

numerically exact results in the systems withparallelly. This tool parallelizes input parameter

arbitrary interactions. Although most numericalsets and performs the simulations independently

data obtained by the exact diagonalizationfor each parameter set. Now we have released

methodthesemore than 20 million data focused on the

representations such as definition of theenergies of the ground state and the first excited

Hamiltonians and the display digits of the results

state on quantum spin systems with small system

are dependent on the articles. Therefore, it issizes [2]. We will upload the results for large

desirable for construction of database to easilysystems and electron systems.

areavailableonarticles,which

performs

HPhi’ssimulations

This project has been done in collaboration

check and obtain the exact diagonalization

with Kazuyoshi Yoshimi and Yuichi Motoyama.

results without definition of notations.

References

In this project, we have constructed an open
database for exact diagonalization results, which[1] K. Yoshimi, K. Ido, Y. Motoyama, R.

is called “COMPutation ARchive of Exact

Tamura, T. Fukushima, Y. Nakanishi-Ohno:

Diagonalization

https://datarepo.mdcl.issp.u-tokyo.ac.jp/repo/3(COMPARED)”[1]. This

database has been available on the ISSP data

[2] https://isspns-gitlab.issp.u-

repository under the CC BY 4.0 license [2]. As a

tokyo.ac.jp/compared/compared

generator of the exact results, we use open-[3] M. Kawamura, K. Yoshimi, T. Misawa, Y.

source software HPhi [3], which supports

Yamaji, S. Todo, and N. Kawashima: Comp.

MPI/OpenMP hybrid parallelization to achieve

Phys. Commun. 217 (2017) 180.

232

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

Ground state and dynamical properties of the
J1-K-Heisenberg model on the square lattice
Matthias GOHLKE, Tokuro SHIMOKAWA
Theory of Quantum Matter Unit, Okinawa Institute of Science and Technology
1919-1 Tancha, Onna-son, Okinawa 904-0495

H = J1
+ K

X
hi,ji

S·S

X 
(i,j,k,l)

(1)


−1
Pijkl + Pijkl
+ hz

X

Siz ,

i

where J1 represents dominant ferromagneticiDMRG Ly = 4

hSiz i0.5

iDMRG Ly = 6

hf-ED N = 6x6
hf-ED N = 8x8
hf-ED N = 10x100.0hSi⊥ i0.2

0.0
1

q = (π, 0)

⊥
SDP
(q)

q = (0, π)01q = (2π, 0)

⊥
SQP
(q)

In the absence of conventional magnetic order of spin-dipolar moments, ordering withhigher-order moments like spin-quadrupoles
may occur. Such spin-quadrupolar—or spinnematic (SN)—order was found theoretically in spin-1 model with biquadratic spinexchange [1].
Magnets with spin- 12 degrees of freedom,
however, can only exhibit a SN states if two
spin- 12 are combined into an effective spin-1 [2].
In fact, such SN states on bonds has theoretically been observed in frustrated ferromagnets
on the square lattice with dominant ferromagnetic Heisenberg exchange, J1 , antiferromagnetic next-nearest neighbor exchange, J2 , andcyclic permutation, K. It was found that thecondensation of a two-magnon bound state, atstrong magnetic fields along the z-axis, can stabilize a phase with bond-nematic order [3].
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 SN: A
continuous symmetry for the director of a spinquadrupole remains, that give rise to a gapless
Goldstone mode [4, 5].
We started by studying the square-latticefrustrated J1 -K model [2],

0
0.0

0.2

0.4

0.6

0.8

1.0

hz

Figure 1: Ground-state of the square-lattice
J1 -K-model in a magnetic field at fixed K/|J1 |.
Different system sizes, geometries, and methods (iDMRG,high-field ED) are compared.
Heisenberg exchange between nearest neighbor spins, K the cyclic ring exchange aroundsquares, and hz the Zeeman coupling to a magnetic field along the z-axis.
Using iDMRG and the matrix product states(MPS) framework as well as a recently developed exact diagonalization method near saturation [6], we confirm the existence of thebond-nematic phase in an extended range of
K/|J1 |. In utilizing different methods combined with different geometries, we can assess
possible finite-size effects, cf. Fig. 1. Overall we find a good agreement between highfield ED on symmetric clusters and iDMRG
on cylindrical geometries with a circumferenceof Lcirc = 6 sites.
A significant amount of computational resources has been spent on mapping out the233

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

non-magnetichz /|J1 |

0.1

ω

0.8
0.6 0.0

0.36

0.38

0.40

101

102
QP
S⊥
(q, ω)S DP (q, ω)

1.2 0.2
1.0

101

0.42

3

3

2

2

ω

1.4

100

hc ED
hc 1-magnon

bond-nematic
4-AF
2-AF Neél

1.6

1

1

0.4
0.2

0

0.0
0.35

0.40

0.45

0.50
K/|J1 |

0.55

0.60

0

Γ

X

Γ0 Γ

X

Γ0

0.65

Figure 2: Phase diagram of the square-lattice
J1 -K-model in a magnetic field obtained withiDMRG using Lcirc = 6 sites circumference.
entire phase diagram using a cylinder with
Lcirc = 6, cf. Fig. 2. The bond-nematic
phase is sandwiched between high-field polarized and 4-sublattice AF phase. Additional
phases (Neél, non-dipolar) occur for large K.
Given the ground state wave function as an
MPS, dynamical properties can then be studied by applying a time-evolution unitary U (dt)represented as a matrix product operator. In
doing so, we observe the condensation of twomagnon bound state at the corresponding wave
vector of the SN state. Within the SN phase,
we observe (1) the magnon-like excitation thatremains fully gapped, and (2) a gapless modewith vanishing spectral weight in the dipolarstructure factor and in the q → 0 & ω → 0
limit as predicted by prior mean-field studies, cf. Fig. 3. We do, however, observe relevant qualitative differences most prominent in
the dynamical quadrupole-structure factor. Its
origin is currently under investigation.

Figure 3: (a) Dynamical spin-structure factorand (b) dynamical quadrupole-structure factor.
[4] R. Shindou, S. Yunoki, and T. Momoi,
Phys. Rev. B 87, (2013) 054429.
[5] A. Smerald, H. T. Ueda, and N. Shannon,
Phys. Rev. B 91 (2015) 174402.
[6] H. Ueda, S. Yunoki, and T. Shimokawa,Comp. Phys. Comm. 277 (2022) 108369.

References[1] K. Harada and N. Kawashima, Phys. Rev.
B 65 (2002) 052403.
[2] A. Andreev and I. Grishchuk, Sov. Phys.
JETP 60 (1984) 267.
[3] N. Shannon, T. Momoi, and P. Sindzingre, Phys. Rev. Lett. 96 (2006) 027213.

234

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

Numerical Diagonalization Study on Quantum Phase
Transitions of Frustrated Spin Systems
Tô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, Japan

1

Magnetization plateau ofthe distorted diamond spin
chain with anisotropic ferromagnetic interaction

The S = 1/2 distorted diamond spin chain
with the anisotropic ferromagnetic interaction
is investigated using the numerical diagonalization and the level spectroscopy analysis. It
is known that the system exhibits a plateau of
the magnetization curve at the 1/3 of the saturation. The present study indicates that as theanisotropy is varied the quantum phase transition occurs between two diﬀerent mechanismsof the 1/3 magnetization plateau. The phase
diagram with respect to the anisotropy and theferromagnetic coupling is also presented[1].

2Magnetization process of
the S=1/2 Heisenberg antiferromagnet on the ﬂoret
pentagon lattice

We study the S = 1/2 Heisenberg antiferromagnet on the ﬂoret pentagonal lattice bynumericaldiagonalization method. We found

235

that magnetization plateaux appear at oneninth height of the saturation magnetization,
at one-third height, and at seven-ninth height.
The magnetization plateaux at one-third andseven-ninth heights come from interactions
linking the sixfold-coordinated spin sites. A
magnetization jump appears from the plateau
at one-ninth height to the plateau at one third
height. Another magnetization jump is observed between the heights corresponding tothe one third and seven-ninth plateaux; however the jump is away from the two plateaux,namely, the jump is not accompanied with any
magnetization plateaux[2].

References[1] T. Sakai, K. Okamoto, H. Nakano and R.
Furuchi, AIP Advances 12 (2022) 035030.
[2] R. Furuchi, H. Nakano, N. Todoroki and
T. Sakai, J. Phys. Commun. 5 (2021)
125008.

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

Numerical Diagonalization Study on Magnetization
Process of Quantum Spin Chain with the
Biquadratic Interaction
Tô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, Japan

1

Quantum spin nematic liquid in the S=1 antiferromagnetic chain with the biquadratic interaction

The magnetization process of the S = 1 antiferromagnetic chain with the biquadratic interaction is invesitgated using the numerical
diagonalization. As a result, it is found thatthe quantum spin nematic liquid phase appears
below the saturation magnetization for suﬃciently large negative biquadratic interaction.
The ground state phase diagram is also presented. [1].

2Field-Induced Quantum
Spin Nematic Liquid Phasein the S=1 Antiferromagnetic Heisenberg Chainwith Additional Interactions

The S = 1/2 distorted diamond spin chain
with the anisotropic ferromagnetic interaction236

is investigated using the numerical diagonalization and the level spectroscopy analysis. It
is known that the system exhibits a plateau of
the magnetization curve at the 1/3 of the saturation. The present study indicates that as theanisotropy is varied the quantum phase transition occurs between two diﬀerent mechanismsof the 1/3 magnetization plateau. The phase
diagram with respect to the anisotropy and theferromagnetic coupling is also presented[2]References[1] T. Sakai, AIP Advances 11 (2021) 015306.
[2] T. Sakai, H. Nakano, R. Furuchi and K.
Okamoto, J. Phys.: Conf. Ser. 2164 (2022)
012030.

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

A proximate Tomonaga-Luttinger state in the S=1/2 Kitaev-
model on a Honeycomb Lattice
Takafumi SUZUKI
Graduate School of Engineering,
University of Hyogo, Shosha, Himeji, Hyogo 670-2280
We have investigated the ground-state phase
diagram and the dynamical properties of the
S=1/2 Kitaev- model [1] on a honeycomb
lattice by the numerical exact diagonalization(ED) method and the density-matrixrenormalization-group (DMRG) method. The
S=1/2 Kitaev-  model on a honeycomb lattice
was proposed as the effective model of thehoneycomb-lattice magnet, -RuCl3. The
ground-state phase diagram of this model with
the isotropic coupling strength has been muchstudied because of the simplicity of the model.
However, the anisotropy of the couplingstrength exists in -RuCl3, which originatesfrom the reduction of the lattice symmetry, andthe effect of the anisotropy has not been fully
understood yet.
In the previous study, we reported the partial
ground phase diagram of the S=1/2 Kitaev-
model on a honeycomb lattice shown in Fig. 1
and indicated that the dimerized phase, whichappears in the isolated dimer limit (d =∞), cansurvives up to the isotropically interacting
model (d=1), when the Kitaev interaction andthe  interaction are negative and positive
respectively [2]. In this study, we investigated
the ground-state phase diagram in 0≤d≤1 [3].
We find that a proximate Tomonaga-Luttinger
liquid (pTTL) is expected in the ground statewith the negative Kitaev interaction and thepositive  interaction (). This pTTL becomes
the Tomonaga-Luttinger liquid (TLL) in thespin chain limit (d=0), which was pointed out
in ref. [4]. When the condition, /|K|=1 (K<0),
is satisfied at d=0, the model is mapped on theantiferromagnetic Heisenberg chains andrecovers the hidden SU(2) symmetry [4]. Thus,
a gapless linear excitation derived from spinon
excitations appears in the low-energy excitation.
These characteristics is also observed in the
low-energy excitation of the pTTL phase. From
the ED calculations for the N=24 site clusterwith the C3 symmetry of the lattice, we findthat the pTTL phase can survive up to d=1.
Next, we calculated the temperaturedependence of the specific heat to investigate
the emergence of the Majorana fermions in thevicinity of d=1 and =0. At d=1 in 0 ≤ /|K|,the presence of the double-peak structure wasindicated in ref. [5]. We find that the highertemperature peak at d=1 corresponds to theprominent single peak in the spin chain limit
and the lower-temperature peak rapidly
develops in the vicinity of d=1. This means thatthe spinon-like excitation is continuously
converted into the Majorana excitation.

Fig. 1: Anisotropically interacting Kitaev-
model. The coupling strength on the X/Y bondis unity and d is the ratio of the couplingstrength on the Z bond against that on the X/Y
bond.

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] M. Gohlke, J. C. Pelayo, and T. Suzuki, in
preparation.
[4] W. Yang, et al., Phys. Rev. Lett. 124,
147205 (2020) and W. Yang, A. Nocera, and I.
Affleck, Phys. Rev. Research 2, 033268 (2020).
[5] A. Catuneanu, et al., npj Quant. Mater. 3, 23(2018).

237

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

Hybrid Simulations on Fluid-Viscoelastic Membrane System
Using Multi-Scale Simulation Platform for Complex Fluids
Yohei MORII and Toshihiro KAWAKATSU
Department of Physics, Tohoku University,
Aza-Aoba, Aramaki, Aoba-ku, Sendai 980-8578
We develop general multiscale simulation

particles with smaller velocity gradient.

Amethods for inhomogeneous complex flows in

preliminary result using linear viscoelastic

viscoelastic fluids and elastoplastic solid. For

constitutive equation was developed during thesuch multiscale simulations, a general-purpose

previous project in 2020 [1]. In the present

platformproject, we extended this hybrid method to

namedMultis-Scale

Simulation

Platform for complex flows (MSSP) was

nonlineardeveloped by the present authors [1]. In our

Figure 1 shows the flow pattern around the

MSSP, the macroscopic flow is simulated byseparation point on a cylindrical obstacle,smoothed(SPH)

where the red SPH particles contains particle

method, to which the microscopic information

MD simulators while blue SPH particles

is introduced through a microscopic molecular

contain nonlinear constitutive equations. We

dynamics (MD) simulator embedded in each of

can confirm that the switching is successfully

the SPH particles.

implemented and 97% of 1,600,000 SPH

particlehydrodynamics

We simulate a viscoelastic flow of a diluteviscoelastic

constitutive

model.

particles are simulated with constitutive model.

polymer solution passing through an obstacle
and clarify the role of added polymers.

To

accurately simulate large scale dynamic flow
patterns behind the obstacle, such as Karmanvortex street, it is essential to use a large-scale
system, where the flow velocity and velocity
Fig.1 MSSP simulation of a polymer
flow around a cylindrical obstacle.
Bars attached to the particles show
the average elongation and its
direction of the polymer chains.

gradient are highly inhomogeneous. To avoid alarge computational cost as well as a large
unbalance in the computational load on thecomputer nodes, we proposed a hybrid method

Referenceswhere we switch the MD simulators with

[1] Y. Morii and T. Kawakatsu: Phys. Fluids 33

simpler constitutive equations for the SPH(2021) 093106.

238

Activity Report 2021 / 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-son Okinawa 904-0412
Quantum frustrated magnets are well-known

them; we found that one of them is ansources of exotic states of matter, but

intriguing six sub-lattice (anti) skyrmion-

investigating their true nature, particularly at

lattice state with topological nature.

finite temperatures, is still a demanding and

2) High-field exact diagonalization code

challenging problem. This fiscal year, wedevelopment and the low-T physics in thehave focused on the following topics with the

S=1/2 J1-K square-lattice model

ISSP supercomputer.

We developed an exact diagonalization code,

1) S=1/2 breathing bilayer kagome (BBK)

quantum spin solver near saturation (QS3), as

magnet for Ca10Cr7O28

an open-source software [2]. Using this high-

Ca10Cr7O28 is a newly quantum spin liquidfield ED code, we investigated the low-T(QSL) material, and the effective

phase diagram of the S=1/2 square-lattice

Hamiltonian was proposed as an S=1/2 BBKHeisenberg model with a ring exchange

Heisenberg model. A recent semi-classicalinteraction (J1-K model).

study [1] supported the scenario of a QSLOur ED calculations revealed the existence

ground state in this material and suggestedof several ordered/disordered states in fields.

possible realizations of exotic low-

Especially in the high-field regime aroundtemperature multiple-q states in fields. But

K/J1~0.5, we clarified the realization of athe true nature of the multiple-q states wasn’t

bond-nematic state without any

uncovered yet because of the difficulty of

approximation to our Hamiltonian. We also

treating large systems.

computed the dynamical spin/quadrupole

We alternatively treated a classical J1-J2

structure factors for future neutron scattering

Heisenberg model on the honeycomb latticemeasurements on candidate materials.

model, which could also be considered aneffective model for Ca10Cr7O28 to investigate

References

the thermodynamic properties of the[1] R. Pohle, et al, PRB 104, 024426 (2021).

multiple-q states. Our large-scale Monte[2] H. Ueda, S. Yunoki, T. S., Comp. Phys.

Carlo simulations succeeded in identifying

Commun. 277, 108369 (2022).

239

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

Steady states realized by a global thermodynamicbalance in nonequilibrium
Naoko NAKAGAWA
Department of Physics, Ibaraki University
2-1-1 Bunkyo, Mito, Ibaraki 310-8512

First-order phase transition and the separation of the phases have attracted scientific interests from various points of view, not only
from physics but also from the biological mechanism. The phase coexistence induces sometimes complicated and singular properties tothe system even though the bulks of respective
phases may behave rather ordinary.
In this project, we work on particle systems
interacting via a Lennard-Jones potential[( )12

V (r) = 4ϵσr−( )6 ]σr+ ϵ,with a diameter σ of each particle and a coupling strength ϵ, which have been a standard
model for numerical experiments. By choosing the total energy and the number density ofthe system appropriately, the system shows aliquid-gas transition. See Fig. 1 for the case of
N = 5 × 104 . We also used Weeks-ChandlerAndersen (WCA) potential when we want toavoid the liquid-gas transition.
Our most interest is the liquid-gas coexistence in heat conduction. Due to the steadyheat current, the coexistence pattern of the liquid and gas becomes much more stable compared with equilibrium systems. The local
states around the liquid-gas interface become
observable due to the induced stability. We
here notice that the rich behaviors related to
the first-order transition are often observed incomplex systems such as mixtures and dense
solutions. Their response to a nonequilibrium
driving will be of great interest. To initiate them, we start with the construction of amethod to know their thermodynamic properties in equilibrium. The details of the results
follow.
We derived a formula for determining themixing free energy for two-component fluids.
The formula is numerically verified with a good
accuracy by the molecular dynamics simulation using the LAMMPS package installed in
the supercomputers at ISSP. We have compared two kinds of protocols creating a mixtureof two pure substances. One is traditional protocol and the other is a protocol with ”alchemical change” of molecules, which may be valid
only in numerical experiments. We calculated
the free energy changes of the two protocols byadopting the Jarzynski equality and information thermodynamics. The comparison of thetwo results in a formula for the mixing Gibbs
free energy ∆mix G asN!
⟨e−β Ŵ(#) ⟩ ∆mix G = −kB T ln ,n!(N − n)!ρX (n) ⟨e−β Ŵ(ii) ⟩where n and N − n are the numbers of particles for the two components and ρX (n) is thenumber density of the first component in theregion x < X. Ŵ(#) is the work in each trajectory according to the alchemical process tocreate a mixture from the first pure substance,
whereas Ŵ(ii) is the work according to the alchemical process to create the second from thefirst pure substance.
By applying the above formula, we examined∆mix G for the mixture of argon and krypton.
By changing the mole fraction of the krypton,

240

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

the system shows liquid-gas transition for a
fixed temperature and pressure. The properties of the transition are well characterized in
the calculated ∆mix G as a double-well-shaped
function. The violation of convexity is likely
due to the finiteness of the system.
We have used the L4 CPU of System B with
4 nodes in 512 parallel computation using MPI
for the calculation of 104 particles, the L16
CPU of System B with 16 nodes in 2048 parallel computation using MPI for the calculationof 5 × 105 particles, the L72 CPU of System
B with 72 nodes in 9216 parallel computation
using MPI for the calculation of 8 × 106 particles.

References[1] A. Yoshida and N. Nakagawa: to appear
in Phys. Rev. Res. (2022)

Figure 1: Snapshot of the particles interactingwith Lennard-Jones potential in heat conduction. The system with N = 5 × 105 exhibits a
sharp liquid-gas interface.

241

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

Scrambling in quantum many-body systems
Masaki TEZUKA
Department of Physics,
Kyoto University, Kitashirakawa, Sakyo-ku, Kyoto 606-8502
Strongly correlated models realizable incondensed matter systems with possible
holographic correspondence to gravitational
systems including black holes have received
much attention in the recent years. Of
particular interest are scrambling and chaotic
dynamics. The scrambling phenomena inquantum systems is delocalization of quantum
information and has been discussed as theformation of quantum error correction (QEC)
code. Thus, we have been trying to study
scrambling and effect of dissipation inquantum many-body systems, toward
experimental demonstration in condensed
matter systems such as ultracold atoms.
The Sachdev-Ye-Kitaev (SYK) model is amodel of Majorana fermions with all-to-all
four-fermion interactions obeying theGaussian distribution. It has an analytic
solution in the limit of large number of
fermions and is maximally chaotic at low
temperatures. We have previously studied theeffect of a random hopping term, two-fermion
interactions, added to the model [1]. The
maximally chaotic behavior is lost and thespectral statistics turns to an uncorrelated one.
Furthermore, the localization transition is a
Fock-space localization, a type of many-body
localization, and quantities such as thelocation of the transition point, the momentsof the eigenstate wavefunctions, and bipartite
entanglement entropy have been obtained
analytically with numerical confirmation [2].
A “sparse” version of the SYK model, in
which the number of interaction terms is on
the order of the number of fermions, insteadof all the choices of four-fermion subsets,

have been proposed a few years ago.
To quantify the scrambling properties ofthese models, we have studied the QEC
capabilities of the unitary time evolution
according to their Hamiltonians according to
the Hayden-Preskill protocol [3]. In thisprotocol, the unknown quantum state throwninto a quantum many-body system is decoded
by a person who, knowing the initial state of
the remainder of the system and the time
evolution, obtains the quantum state of someof the qubits comprising the system, and theerror estimate for the time evolutions obeying
the circular unitary ensemble (CUE) is known.
Using the ISSP Supercomputer, weanalyzed the spectral statistics of a further
simplification of the sparse SYK model [4].
Also, we numerically observed that while forthe sparse SYK models, the QEC error
estimate approaches the CUE value as soon as
the spectral statistics becomes random-matrix
like, for the case of two-fermion interactions,the departure from the CUE value is observed
before the Fock space localization occurs [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] P. Hayden and J. Preskill, JHEP 0709
(2007), 120.
[4][5] M. Tezuka et al., in preparation.

242

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

Investigation of the molecular origins of the mechanical andthermal properties of realistic biopolymers using all-atomistic
molecular dynamics
Okazaki Susumu
Department of Advanced Materials Science,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
Mechanical and thermal properties ofpolymer chains forming network structures.

biopolymers, PBI and its three copolymers

We also investigated the microscopic origins

PBO-co-PBI, PBI-co-PBO-co-pPA, and PBI-of the stress causing the elastic deformation,

co-PBO-co-mPA, have been investigated based

yield point, strain softening, and strain

on all-atomistic molecular dynamics (MD)hardening observed in PC and OMMA, by

calculations to clarify molecular mechanism of

decomposing the stress into contributions fromfracture and glass transition of the polymers. We

various energetic energy terms. Our analysis

also investigated regular polymers, PC andshowed that the stress of the fracture process is

PMMA [1], as a reference system. The

mainly determined by the bond, angle, and vdW

calculations have been conducted on theterms, in addition to the system alignment in the

supercomputers at ISSP using a MD software

pulling direction. In the elastic region, almost

MODYLAS [2] developed by our group.

theModulus of elasticity, yielding stress, andentiresystemexperiences

affinedeformation

despite

the

local

existing

fracture stress were studied in detail. The

heterogeneous deformation in the later elastic

biopolymers showed the high modulus ofregion. As the deformation of the entanglement

elasticity at short strain followed by the highnetwork begins, the stress on the bonds andyielding stress and ductility showing strong

angles rapidly increases.

strain hardening at large strains after theReferences

yielding point.

[1] Z. Tang, K. Fujimoto, S. Okazaki, Polymer,

The biopolymers investigated here all
exhibited

very

high

glass

226(2021)123809.

transition

temperatures at 800 K – 900 K. The[2] Y. Andoh, S. Ichikawa, T. Sakashita, N.

temperatures are higher than their pyrolysis

Yoshii,temperatures. These high glass transition

42(2021)1073.

temperatures are from hydrogen bonds among

243

S.

Okazaki,

J.

Comp,

Chem.,

Activity Report 2021 / 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
The development of quantum computers has

use charges high fees; error rates of today’s

made significant progress in recent years:

NISQs are still too high. Last year, we built aquantum computers with dozens or a hundred

framework to efficiently emulate quantum

qubits have been built. We now expect to see

circuits with (deplolarsing) noises using thequantum devices with hundreds or thousands ofsupercomputer of the ISSP.

qubits near future. Such quantum computers areThis year, we worked on improving theoften called noisy-intermediate scale quantum

efficiency of the quantum circuit optimization

devices (NISQs) since they do not have quantum

part in the VQE. We parallelized the gradienterror correction functionality. There have been

evaluations of the VQE’s cost function withintensive research efforts currently underway torespect to the quantum circuit parameters, using

exploit the computational potential of NISQs.

MPI. We tested our implementation by

One promising algorithm for NISQs is thecomputing the H2O molecule with the 6-31G

quantum-classical hybrid algorithm, which canbasis set. The benchmark calculation used aoperate even in the short coherence time of aminimum of 128 cores and a maximum of 9216

NISQ. Among them, the variational quantum

cores. The observed parallelization efficiency

eigensolver (VQE) can be applied to first-

was 75% with 4608 cores, relative to theprinciples calculations of materials or molecules.

execution time with 128 cores.

We are targeting applying the VQE to practical

Moreover, we developed a method toproblems in material science. In particular, weconsider solvation effects within the framework

aim to use the VQE for analyzing materials usedof the VQE. This was achieved by combining thein sensors. In this study, we use a quantum

VQE and 3D-RISM, an integral equation theory

circuit simulator rather than real quantumof molecular liquids. This method was used tocomputers to develop a new quantum-classicalinvestigate the impacts of solvents on the

hybrid algorithm and validate it. The reasons wecomputational cost of quantum computing. Our

employ such an emulator are the following: thecalculations indicate that solvation effects

number of NISQs is limited at this moment; theirvirtually have no impact on it.

244

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

Molecular Dynamics Simulation of Stress Induced
Crystallization in Robust-Toughening Hydrogel
Koichi Mayumi, Kosuke Aomura, Takeyoshi Masumoto, Haru Shinohara

Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
In

this

project,

westudiedcrystallization of polymer chains inslide-ring gels under stretching by
means

of

coarse-grained

dynamics simulation.

molecular

In addition, we

performed a full-atomistic molecular
dynamics simulation of the melting
process of stretched polymer crystals inwater.
Fig. 1 Stress-extension ratio relation of
fixed crosslink gel and slide-ring gel.

Coarse-Grained Molecular Dynamics
Simulation of Strain-Induced Crystalli-

links.

zation of Slide-Ring Gel
In

this

study,

uniaxial

Further analysis was performed toelongation

simulations of slide-ring gels wereclarify

performed

coarse-grained

stretching. We calculated scattering

molecular dynamics simulation [1]. In

function of SR gel before and afterorder to investigate the strain-induced

stretching. At the extension ratio λ=12,

crystallization of slide-ring gel which isscattering function shows a circularobserved
introduced

using

in

experiments
attractive

[2],

we

aggregatedpeak on 𝑞

interaction

0σ

structureplane at

𝑞

under

𝑞

between polymer chains. Fig. 1 shows

6 σ . The peak suggests the formationstress-extension ratio curve of a fixedofcrosslink (FC) gel and slide-ring (SR) gel

polymer chains are aligned in thewith the same cross-linking density. SR

directionshows smaller stress than FC gelelongation.

because of the slidability of the cross-245

aggregate

structures

perpendicular

in
to

which
the

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

Analysis of Dissolution Process ofExtended Polyethylene Glycol Crystal
in Water
The reversibility of the strain-induced
Fig. 2 Snapshots of a scattering function

crystallization of PEG chains in waterof SR gel at (a) λ=1 (b) λ=12

allows the tough hydrogel to rapidly
recover from mechanical stress. To

A region of neatly aligned polymer

understand the dissolution process of

chains was extracted from the network,

extended PEG chain crystal in water,

and

was

all-atom MD simulation was carried out.

calculated for the region. As a result, we

GROMACS 2016.5 [3] was used as

observed spot peaks corresponding to asimulation software. We investigatedbody-centeredlatticethe dissolution process of helical PEG

crystal structure, as shown in Fig. 3.

chain crystal and zigzag PEG chain

This indicates that the slidability of thecrystal under stretching in water. For

cross-linking points equalizes the chain

helical chain crystal, the crystal partly

deformation and enhances the strain-remain after 200 ns even at no

induced crystallization.

elongation force. On the other hand,thescatteringfunctionorthorhombic

zigzag chain crystal dissolves fasterthan the helix chain crystal.

References[1] Kremer, K.; Grest, G. S.: J. Chem.

Phys. 92 (8), 5057–5086. (1990)[2] C. Liu, N. Morimoto, L. Jiang, S.
orthorhombic

Kawahara, T. Noritomi, H. Yokoyama, K.

lattice crystal structure observed in SR

Mayumi, K. Ito: Science 372 1078-1081.

gel at =12.

(2021)Fig.

3

Body-centered[3] Berendsen, H. J. C.; van der Spoel,D.; van Drunen, R. Comp. Phys. Comm.
91, 43–56. (1995)

246

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

Molecular dynamics study of shear-induced
long-range correlations in simple ﬂuids
Hiroyoshi Nakano
Department of Applied Physics and Physico-Informatics, Keio University
Yokohama, Kanagawa 223-8522, Japan

We studied a nonequilibrium long-range correlation (LRC) in simple ﬂuids using a molecular dynamics (MD) simulation. The nonequilibrium LRC means that same-time spatial correlations of locally conserved quantities exhibit
power-law decays due to a lack of detailedbalance condition [1]. It is known to be oneof the general features of nonequilibrium systems with conserved quantities and anisotropy.
Because in equilibrium the LRCs are speciﬁcto the critical phenomena, the nonequilibrium
LRCs observed far from the critical points
have been extensively studied as the novel
phenomena induced by nonequilibrium conditions. However, the almost all of research
have been performed using phenomenologicalmodels such as ﬂuctuating hydrodynamics andstochastic lattice gases. Indeed, there are few
attempts to directly observe the nonequilibrium LRCs and examine their nature from underlying molecular systems.
In the present study, with the aim of direct observation of the nonequilibrium LRCs,
we focused on the simple ﬂuids under uniform shear ﬂow and performed the MD
simulation [2].
All the MD simulations
were performed by LAMMPS (Large-scale
Atomic/Molecular Massively Parallel Simulator) [3, 4]. The schematic illustration of oursystem is given in Fig. 1. The particles interact via the Weeks–Chandler–Andersen (WCA)
potential. The shear ﬂow is imposed by using
the Lees–Edwards boundary condition along
the z-axis. Along the x and y-axes, the stan-

Figure 1: Schematic illustration of shear ﬂow.
dard periodic boundary conditions are imposed. To maintain a constant temperature T
under shear ﬂow, the dissipative particle dynamics (DPD) thermostat is used. Thus, inthe steady state, the following averaged velocity proﬁle ⟨v(r)⟩ss = (γ̇z, 0, 0) is realized,
where ⟨·⟩ss denotes time average in the steadystate and the ensemble average over diﬀerent
noise realizations, and γ̇ is shear rate.
We observe the same-time correlations of
momentum ﬂuctuations in the steady state.
The momentum ﬂuctuations are given by
δg(r, t) =N∑pi (t)δ(r − ri (t)) − ⟨v(r)⟩ss , (1)i=1where ri (t) and pi (t) are, respectively, the position and momentum of ith particle at timet. Then, the same-time correlations of momentum ﬂuctuations are given by Cαβ (r, r ′ ) =⟨δg α (r, t)δg β (r ′ , t)⟩ss , where α, β = x, y, z.
We plot the Fourier transform of Czz (r, r ′ )in Fig. 2. Because our system has the spe-

247

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

7

7

3
1

5
3
1

0.04

0.08

0.12

0.16

0.08

0.12

0.16

0.04

7

1

5
3
1

0.12

0.16

0.16

Lz = 1024σ
Czz (kx, 0, 0)

3

0.12

7

Lz = 512σ
Czz (kx, 0, 0)

5

0.08

kx

Lz = 256σ

0.08

3

kx

7

0.04

5

1
0.04

kx

Czz (kx, 0, 0)

Lz = 128σ
Czz (kx, 0, 0)

5

Lz = 64σ
Czz (kx, 0, 0)

Czz (kx, 0, 0)

Lz = 32σ

5
3
1

0.04

kx

0.08

0.12

0.16

kx

0.04

0.08

0.12

0.16

kx

Figure 2: Czz (kx , ky = kz = 0) for several system sizes Lz = 32σ, 64σ, 128σ, 256σ, 512σ, 1024σ.
The color circle points represent the MD simulation results, and the black square points theresults calculated from the ﬂuctuating hydrodynamics.
cial Galelian invariance, the Fourier transform of Czz (r, r ′ ) is given as Czz (k) =∫ 3
d rCzz (r, 0)e−ik·r . In Fig. 2, Czz (kx , ky =kz = 0) for diﬀerent Lz is presented. The otherparameters are set to Lx = 1024σ, Ly = 32σ,
−1and γ̇ = 0.02τunit, where σ and τunit are, respectively, the unit of length and time. The
black square points represent Czz (kx , ky =kz = 0) calculated from the phenomenological theory called the ﬂuctuating hydrodynamics [5]. The result obtained from the ﬂuctuating hydrodynamics is independent of Lz .

rium long-range correlations for the suﬃciently
large system that contains over 10 million particles.

From Fig. 2, we ﬁnd that the strong ﬁnitesize eﬀects exist for Czz (kx , ky = kz = 0). For
suﬃciently large Lz , the result calculated fromthe ﬂuctuating hydrodynamics is in quantitative agreement with the MD simulation results.
In contrast, for smaller Lz (e.g. 32σ or 64σ),
there is the clear deviation between their results. In other words, the quite large systemsize is required to observe the shear-induced
long-range correlations.

[3] S. Plimpton, Journal of Computational
Physics 117, 1 (1995).

In summary, we performed the MD simulations of the WCA ﬂuids under uniform shear
ﬂow. We successfully observed the nonequilib-

References[1] J. Dorfman, T. Kirkpatrick, and J. Sengers, Annual Review of Physical Chemistry 45, 213 (1994).
[2] H. Nakano, and Y. Minami, Physical Review Research (in press).

[4] A. P. Thompson, H. M. Aktulga, R.
Berger, D. S. Bolin-tineanu, W. M.
Brown, P. S. Crozier, P. J. in t Veld,A.
Kohlmeyer, S. G. Moore, T. D. Nguyen,
R. Shan, M. J.Stevens, J. Tranchida, C.
Trott, and S. J. Plimpton, Computer
Physics Communications 271, 108171(2022).
[5] J. Lutsko and J. W. Dufty, Phys. Rev. A
32, 3040 (1985).

248

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

Theoretical analysis of absorption and fluorescence
spectra for firefly bioluminescence related molecules
Miyabi HIYAMA
Graduate School of Science and Technology, Gunma University,
Tenjin-cho, Kiryu, Gunma 376-8515

Because firefly bioluminescence is widelyused such as in vivo imaging and reporter assay, many analogs of firefly bioluminescence
substrates (luciferin) are developed [1, 2]. AkaLumine is one of the analogs of luciferin andits bioluminescence with firefly luciferase produces the near-infrared emission, which canpenetrate through biological tissues [3]. For
studying the bioluminescence mechanism, basic knowledge about the chemical structures,electronic states, and absorption properties of
AkaLumine solution at various pH values ofsolution has to be acquired.
In this study, the absorption spectra for
Alakumine in aqueous solutions were studied
theoretically [4]. The optimization structureof ground state for AkaLumine and its conjugate acids and bases, which are expected to be
main component in the aqueous solutions, wereobtained from the density functional theory(DFT). The time dependent DFT (TDDFT)
calculations were carried out to estimate thetheoretical absorption spectra for these structures. The Gibbs free energies for the chemical
structures were obtained using the vibrational
analysis. Form these free energies, pKa values
for protonation/deprotonation of AkaLumine
in aqueous solutions were estimated and therelative concentrations of AkaLumine and its
conjugate acids and bases were obtained. All
calculations were performed using the GAUSSIAN09 program [5] on system B and C of Super Computer Center in ISSP.

for AkaLumine and its conjugate acids adbases appear at 369, 465, 314, 353, 352, 440,
322, and 317 nm and that all of these peaks
correspond to the excitation from S0 to thefirst excited state (S1 ) [4]. These S1 states
can be described dominantly by the one electron excitation configuration of which the coeﬃcients are large (∼0.5). From the theoretical results, the peaks at pH 7–10, pH 4, andpH 2 of AkaLumine experimental absorption
spectra are assigned to a carboxylate anion, acarboxylate anion with an N-protonated thiazoline ring, and a carboxylate anion with an Nprotonated thiazoline ring and N-protonateddimethylamino group, respectively.

References[1] Li et al. , Eur. J. Med. Chem., 211 (2020)
113111.
[2] Iwano et al. , Tetrahedron, 69 (2013)3847.
[3] Kuchimaru et al. , Nat. Comm., 7 (2016)
1.
[4] Ogawa et al. Photochem. Photobiol.
(2021) 97, 1016.
[5] Gaussian 09, Revision D.01, M. J. Frischet al.

It was found that the main absorption peaks

249

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

Skyrmion crystal phase in the RKKY Heisenberg modelin two- and three- dimensions
Kota Mitsumoto and Hikaru Kawamura
Molecular Photoscience Research Center, Kobe University, Kobe 657-8501, Japan
Skyrmion, a swirling noncoplanar spin

three-dimensional stacked-triangular lattices by

texture whose constituent spin directions wrap a

extensive Monte Carlo simulations in order tosphere in spin space, has attracted much recent

get insights into the chiral-degenerate skyrmion

interest and been studied quite extensively. In

crystal of metallic magnets in the weak-coupling

magnets, skyrmion is usually stabilized in theregime [2].

form of the periodic array, i.e., the skyrmion

We find that the symmetric SkX is stabilized

crystal (SkX). While the SkX was first realized

under finite magnetic fields both in two andin chiral ferromagnets with the anti-symmetricthree dimensions, as shown in the temperature

Dzyaloshinskii-Moriya (DM) interaction, avs. magnetic-field phase diagram of the Figure

recent study has revealed that the SkX state canbelow, both the 2D (upper) and 3D (lower) ones.

also be realized in centrosymmetric magnets
without the DM interaction [1]. Such a“symmetric” SkX state is induced by frustration
among exchange interactions, and keeps thechiral degeneracy with respect to the right-left(R-L) symmetry, in sharp contrast to the DMinduced SkX.
Such a frustration-induced symmetric SkX
is possible both in insulating and metallic
magnets. In metallic magnets, the primaryinteraction is the long-range RKKY interaction
References

which oscillates in sign with the distance leading[1] T. Okubo, S. Chung and H. Kawamura, Phys.

to the frustration, in contrast to the insulating

Rev. Letters, 108, 017206-(1-5) (2012).

magnets possessing the short-range interactions.

[2] K. Mitsumoto and H. Kawamura, Phys. Rev.

In this year’s project, we study the ordering

B 104, 184432-(1-17) (2021); Phys. Rev. B 105,

properties of the isotropic RKKY Heisenberg

094427(2022).

model on the two-dimensional triangular and the250

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

Topological phases and science of bulk-edgecorrespondence by numerical methods
Y. Hatsugai
Department of Physics, University of Tsukuba
1-1-1 Tennodai, Tsukuba 305-8571, Ibaraki, JapanThe bulk-edge correspondence[1] is widely used for the studies of topological phenomena
in various systems that include not only microscopic quantum systems but also macroscopic
classical ones such as photonic crystals and meta-materials. It is a scale-free concept. The
abstract mathematical concept of topology is reflected by the physics of boundaries throughgeneric edge states (zero to three-dimensional defects as well). We may further try to discussvarious kinds of localized states as a consequence of non-trivial bulk. A typical example is anequatorial wave on earth where the equator is a boundary for the Coriolis force. Collecting
such hidden reasons why the boundary states exist, we try to establish a universal view for the
localized states in nature and more. This is the science of bulk-edge correspondence. In thecurrent project, we try to apply the viewpoints for various systems using numerical techniques.
Topological pumps focusing on the role of edges are important systems where the bulk-edgecorrespondence governs the phenomena. We have numerically justified the bulk-edge correspondence of topological pumps for quantum spins and bosons with short-range interactions[2, 3]
where a plateau transition analogous to the quantum Hall effects is realized in dimerized high-S
quantum spins. Laughlin’s argument of the quantum Hall effects can be assumed as a variantof the topological pump where bulk-edge correspondence plays a fundamental role. It is numerically considered for a series of anyons with various statistical parameters associated with the
adiabatic heuristic argument[4]. Also, its spinful extension is discussed in details[5].
Fermion quantum Monte Carlo techniques are applied for three dimensional Hubbard model
on a pyrochlore lattice where the bulk-corner correspondence has been justified[6]. We further develop bulk-edge correspondence to non-quantum systems. One is heat conduction in 1dimension[7]. Surprising another example is a bulk-edge correspondence for the evolutional gametheory[8] where non-zero Chern number characterizes a boundary flow in the two-dimensional
game theory (Figure 1.).

References[1] Y. Hatsugai, ”Chern Number and Edge States in the Integer Quantum Hall effect”, Phys. Rev. Lett.
71, 3697 (1993).
[2] Y. Kuno and Y. Hatsugai, “Plateau transitions of a spin pump and bulk-edge correspondence”,
Phys. Rev. B 104, 045113 (2021).
[3] Y. Kuno and Y. Hatsugai, “Topological pump and bulk-edge-correspondence in an extended BoseHubbard model”, Phys. Rev. B 104, 125146 (2021).
[4] K. Kudo, Y. Kuno, and Y. Hatsugai, “Bulk-edge correspondence in the adiabatic heuristic principle”,
Phys. Rev. B 104, L241113(2021)[5] K. Kudo and Y. Hatsugai, “Adiabatic Continuity of the Spinful Quantum Hall States”,
arXiv:2201.07893

251

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

Figure 1: RPS-game on a Kagome lattice (left) and its bundary flow governed by non zero Chern
number of the evolutionary game theory (right). (See ref.[8].)[6] Y. Otsuka, T. Yoshida, K. Kudo, S. Yunoki and Y. Hatsugai , “Higher-order topological Mott
insulator on the pyrochlore lattice”, Scientific Reports 11, 20270 (2021)[7] S. Makino, T. Fukui, T. Yoshida and Y. Hatsugai, “Edge states of a diffusion equation in one
dimension: Rapid heat conduction to the heat bath”, Phys. Rev. E 105, 024137 (2022).
[8] T. Yoshida, T. Mizoguchi, and Y. Hatsugai, “Chiral edge modes in evolutionary game theory: A
kagome network of rock-paper-scissors cycles”, Phys. Rev. E 104, 025003 (2021).

252

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

Molecular dynamics simulation of BaTiO3 nanostructure
T. Hashimoto
National Institute of Advanced Industrial Science and Technology (AIST),
Tsukuba Central 2, 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan

Dielectrics clusters have fascinating newproperties. It is shown that large oﬀ-center displacements of atoms exists even for very small
clusters as small as 5 nm, and vortex structuresof local polarization vectors could be used as
memory devices[1, 2]. However, despite theeﬀorts of preceding studies, the surface termination dependence of the vortex structures
and the toroidal moments are not well understood. In this project, we studied these issues
by molecular dynamics (MD) simulations using a shell model.
The calculations were performed with thecode developed by us. We used a shell model
potential developed by Tinte et al.[3, 4]. The
calculations were performed under constant(N, V, T) conditions. To simulate isolated
clusters, the open circuit boundary condition
was used; i.e., we calculated interactions for
all the atom pairs in the cluster. The temperature was ﬁxed to 10 K to monitor the structures near the ground state while preventing
the system from converging to false structures.
We controlled the temperature by the massive
Nose-Hoover chain method[5]. For each clustermodel, the system was equilibrated for morethan 20000 time steps, and data collection runs
for more than 100000 time steps followed. One
time step was 0.1 fs.
We studied BaO-terminated,
TiO2 terminated,and stoichiometric models.
The former model was rounded apex cubicmodels with eight Ba atoms missing, while thelatter two were cubic models. The number ofthe unit cells on one edge, n, were varied from
3 to 12. The TiO2 terminated surfaces were
unstable, and the TiO2 terminated clusters
were unstable for n ≤ 5. In all terminations,the toroidal moments become non-zero forclusters larger than n = 6. We also studied thetoroidal moment direction dependent energy
for the BaO-terminated clusters, because they
were stable for all n. In these clusters, thestructures with the toroidal moments parallelto [111] which are naturally formed in otherstudies were stable, and structures in which
the toroidal moments were parallel to [110]or [001] were metastable. We also analyzed
the domain structures inside these clusters.
They were found to have domain walls thatare known to be stable in the rhombohedralphase.
[1] H. Fu and L. Bellaiche: Phys. Rev. Lett.
91 (2003) 257601.
[2] I. I. Naumov, L. Bellaiche, and H. Fu: Nature 432 (2004) 737.
[3] S. Tinte, M. Stachiotti, S. Phillpot,
M. Sepliarsky, D. Wolf, and R. Migoni: J.
Phys.: Condens. Matter 16 (2004) 3495.
[4] M. Sepliarsky, A. Asthagiri, S. R.
Phillpot, M. G. Stachiotti, and R. L.
Migoni: Curr. Opin. Solid. St. M. 9 (2005)
107.
[5] G. Martyna, M. Tuckerman, D. Tobias,
and M. Klein: Mol. Phys. 87 (1996) 1117.

253

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

DPD simulations of cross-linked networks to study
topological effect using MP-SRP method
Katsumi HAGITA,1 Takahiro MURASHIMA,2 and Hayato SHIBA3
1

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

Department of Physics, Tohoku University, Sendai 980-8578

Information Technology Center, University of Tokyo, Kashiwa 277-0882

Ring-linear blends have garnered increasing

dynamics (DPD) simulations [6]. Multipointinterest among researchers in recent years forsegmental repulsive potential (MP-SRP) [7]both basic science and industrial applications.

was used to ensure thermodynamic consistency

In ring–linear blends, a deeper understanding ofwhile comparing between the systems thatpenetration behaviors among rings and linear

permitted and prohibited chain crossing in DPD

chainsthepolymer simulations. To extract the topologicalmechanical properties of the blends and theireffect of ring polymers, we observed theapplications. Hagita and Murashima [1–4]distribution of penetration in ring–linear blends

investigated the dependence between ring size

and the radius of gyration of rings.

isimportantforcontrolling

and the penetration of linear chains into the ring

Based on a prototype code as an in-housethrough coarse-grained molecular dynamics

user package of MP-SRP [7] for LAMMPS [8],(CGMD) simulations of the Kremer–Grest

we remade a public open version of the user

(KG) model [5]. They estimated the distribution

packageof the number of chains penetrating a ring bycomputational cost, we improved several

estimating the Gauss linking number (GLN) for

implementations as shown in Figure 1. We also

all pairs of linear chains and rings. The number

confirmed the high scalability in the MPI-basednP of linear chains penetrating a ring wasparallel computing.

of

MP-SRP[7,9].

To

reducepresented as a function of the size of the ring(the number of KG particles in a ring, Nring).
These penetration behaviors are governed by
the spread and concentration of the ring
polymers.
In this study, we investigated the topologicalFig. 1: Schematic of the applied improvement

effect caused by chain crossing prohibition in

to reduce omissible computations by changing

ring-linear blends through dissipative particle

sweep of the SRP points along a bond.

254

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

uncatenated rings repel each other because of

In this study, DPD simulations (MP-SRP

the topological constraints.

and sDPD) of ring–linear blends with a ringfraction of approximately 0.1. Here, sDPD

We performed systematic investigations formeans the standard DPD polymer simulation

other values of (Nring, Nlinear) as shown in thewith the Groot–Warren model [6].

published paper [9]. We concluded that the

Figure 2 shows the results of (Nring, Nlinear)

topological effect prohibits ring–ring crossing,= (80, 80), where Nring and Nlinear denote theresulting in a small radius of gyration of thenumber of DPD particles in a ring and in arings. In the DPD polymer system with MP-

linear chain, respectively. Compared to the

SRP, the minimum size of a ring withsDPD model, the number of penetrations in thepenetration was found to be 30, which is less

MP-SRP model decreased owing to thethan half of the value of 80 observed inprohibition of ring–ring crossing. It should be

previous work using the Kremer–Grest bead-

noted that the penetration state of the initialspring model.

configuration

generatedbythepush-off

procedure was far from equilibrium in both

Referencesmodels.

[1] K. Hagita and T. Murashima: Polymer 218(2021) 123493.
[2] K. Hagita and T. Murashima: Polymer 223(2021) 123705.
[3] K. Hagita and T. Murashima: Macromolecules
54 (2021) 8043-8051.
[4] T. Murashima, K. Hagita, and T. Kawakatsu:
Macromolecules 54 (2021) 7210-7225.
[5] K. Kremer and G. S. Grest: J. Chem. Phys. 92(1990) 5057.
[6] R. D. Groot and P. B. Warren: J. Chem. Phys.

Fig. 2: Probability distributions of linear chains

107 (1997) 4423.

penetrating into a single ring with (Nlinear, Nring)[7] N. Iwaoka, K. Hagita and H. Takano: J. Chem.

= (80,80).

Phys. 149 (2018)114901.
[8] S. Plimpton: J. Comp. Phys. 117 (1995) 117.

We also confirmed the decrease of the[9] K. Hagita, T. Murashima, H. Shiba, N. Iwaoka,squared radius of gyration in comparison to theand T. Kawakatsu: Comput. Mater. Sci., 203

phantom chain of the sDPD model. The(2022) 111104.

255

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

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

The problem of what drives a liquid-liquidtransition motivated us to study entropy estimation for liquids. We first report the results
attained in FY2021. The latter part of the report consists of the outcomes obtained froma study of memory hysteresis using machinelearning (ML).
The two-phase model’s justifiability for entropy
estimation [1]The two-phase model [2] extends the onephase model for describing the thermodynamics of crystalline phases, regarding them as anaggregate of harmonic oscillators with 3N degrees of freedom (DF). Then, the logarithm of
the partition function Z is given by ln Z =∑j ln zj , where zj is the partition function of
the jth oscillator. If the microscopic states
are densely distributed to allow a continuum
description, we may rewrite this relation as∫
ln Z = dϵD(ϵ) ln z(ϵ), where D(ϵ) defines thedensity of states (DOS) in the energy range between ϵ and ϵ + dϵ. Berens et al. [3] prescribedobtaining the DOS from Fourier transforming
the velocity autocorrelation function. Thus,
the entropy for a solid state may be approxi∫∞mated as S = kB 0 dωD(ω)WSHO (ω), where
kB is Boltzmann’s constant. The analytic expression for the entropy spectrum WSHO is easily obtainable.
The one-phase model is not applicable to estimating thermodynamic properties of liquidand amorphous states with non-negligible lowfrequency modes of D(ω). Goddard et al. [2]have constructed the two-phase model basedon their intuitive observation that the DOS ofsuch non-crystalline systems looks like a superposition of DOSs of solid and gas states. They

introduced the notion “fluidicity” f , in termsof which the DOS of non-crystalline systems
could be divided into 3N f DFs of gas component Dg and remaining crystalline component
Ds . Thus, the two-phase model modifies theabove expression for the entropy as∫ ∞S/kB =dω Dg (ω|f )WSg (ω)
0∫ ∞+dω Ds (ω|1 − f )WSHO (ω). (1)0Goddard et al. [2] assumed a hard-sphere gas
to evaluate Dg and WSg . For WSg , they employed the entropy per DF for a HS system
WSHS obtainable from the Carnahan–Starling(CS) equation of state (EOS), whereas they
invoked the Chapman–Enskog (CE) theory tosettle the expression for Dg . The CE theory
also allowed them to determine f . Thus, if theactual DOS D(ω) is available from experiments
or simulations, Eq. (1) gives us an entropy estimate by putting Ds = D − Dg .
We employed a simple liquid whose thermodynamic properties are well studied to test
the two-phase model’s prediction. Thus, wechose the modified Lennard-Jones system [4]with 2048 particles, to which we conducted
isothermal–isochore molecular dynamics simulation at slightly above the critical temperature. We adopted three densities; one was thecritical density ρc , and the other two are below and above ρc . We used the λ- and thermodynamic integrations to estimate the entropies
precisely. We calculated the entropy also fromthe modified Benedict–Webb–Rubin EOS [5].
The entropies resulting from the former
three approaches coincided within the error margin. However, the two-phase model256

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

gave systematically lower entropies than those
three. The deviation increases with density.
However, the most significant deviation was
limited to 6% at the highest density examined. Because both the CS EOS and the CE
theory would not be appropriate to capturesuch a high-density liquid (even though simple), many factors employed in the two-phasemodel would cancel to attain such an accuracy. Unfortunately, because the two-phasemodel was constructed in an ad-hoc way, we
cannot give any suggestions to improve it at
this moment.
Judging the (dis)continuity of phase transition
using ML [7]
Based on the valuable discovery that thesum of weights Wsum in a convolutional neural network (CNN) behaves like an order parameter when it learns spin configurations ofa 2-dimensional Ising system as a functionof (inverse) temperature, we have proposed amethod to detect the phase-transition point at
an early stage of learning processes. The latterwork prompted us to consider whether a CNN
can judge the order of phase transition fromthe viewpoint of memory hysteresis.
When we see a series of continuously changing pictures from, say, A to B, we notice asalient change at C on its way. However, when
we view the same series in reverse order from
B to A, the noticeable point of change does
not necessarily coincide with C. This hysteresis may be related to the mechanism of ourshort-term memory equipped with our brain.
A CNN could already implement such a mechanism. Thus, it will exhibit a hysteretic behavior depending on an impetus applied througha learning process.
Here we only outline what we examined.
We used the same CNN as employed in [6].
However, the crucial diﬀerence was the wayof learning. In [6], learning data has been input randomly, as usually employed in ML. Thistime, the data was input in order of increasing(decreasing) magnitude of external parameterssuch as magnetic field and temperature. Undervarious magnetic fields, we prepared equilibrium spin configurations of the 2-dimensional
64 × 64 Ising system below and above the critical temperature. We also made equilibrium
patterns of 256 q-state Potts model for q = 3
through 7 at various temperatures below andabove the critical temperature. System B tookcare of generating these data sets.
In all of the cases examined, for a discontinuous phase transition, a sharp peak showed up
precisely at the transition point when Wsum
was plotted as a function of the externalparameter, regardless of the direction of itschange. In contrast, Wsum oﬀered a roundedpeak for a continuous transition with the peak
position shifted as if the CNN missed the exacttransition point, leaving a clear-cut hysteresis
against a cyclic change of the external parameter. We thus conclude that our approach will
be applicable to judge the order of a generalphase transition if the learning data are available.

References[1] H. Naruta and K. Fuchizaki: to be presented at The 18th International Conference on Liquid and Amorphous Metals.
[2] S-T. Lin, M. Blanco, and W. A. Goddard:
J. Chem. Phys. 119 (2003) 11792.
[3] P. H. Berens, D. H. J. Mackay, G. M.
White, and K. R. Wilson: J. Chem. 79(1983) 2375.
[4] Y. Asano and K. Fuchizaki: J. Chem.
Phys. 137 (2012) 174502.
[5] Y. Asano and K. Fuchizaki: J. Phys. Soc.
Jpn. 86 (2014) 034601.
[6] K. Fuchizaki, K. Nakamura, and D. Hiroi:
J. Phys. Soc. Jpn. 90 (2021) 055001.
[7] K. Nakamura and K. Fuchizaki: in preparation.

257

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

Three-Dimensional Finite Element Analysis of Friction
between solids
Michio OTSUKI
Graduate School of Engineering Science, Osaka University,
Osaka University, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
According to Amontons’ law, the maximum
static friction force of a solid object is
proportional to the applied loading force andindependent of the apparent contact area, which
Rigid rodindicates that the static friction coefficient does
not depend on the applied pressure and the

Sliding
direction
Viscoelastic
block
Rigid substrate
Friction force

Fig. 1: Schematic of the 3D FEM calculation.

shape of the solid. Previous studies haveclarified that Amontons’ law systematically

The object exhibits macroscopic stick-slip

breaks down for macroscopic objects due tomotions, which are characterized by the

precursor slip before the onset of bulk sliding

macroscopic static friction coefficient 𝜇& . 𝜇&[1], which is consistent with experiments [2].

decreases with the external pressure 𝑃!"# andHowever, the analysis in the previous works is

the width 𝑊 , as shown in Fig. 2, which

restricted to 2D systems. Therefore, it is notcontradicts Amontons' law. We also find thatclear whether the result is applicable to morequasi-static precursor slip occurs in the 2D

realistic 3D systems.

interface between the object and the substrateThis year, we numerically investigate thebefore bulk sliding, which is related to thesliding motion of a 3D viscoelastic object withdecrease of 𝜇& .

length 𝐿 , width 𝑊 , height 𝐻 , and Young’s
modulus 𝐸 on a solid substrate using the finite
element method (FEM) with 500,000 nodes by
MPI of 1000 processors, as shown in Fig. 1. A
Fig. 2: Static friction coefficient 𝜇! against

rigid rod slowly pushes the object with aconstant velocity. A uniform pressure 𝑃!"# isReferences

applied to the top of the object. The friction[1] M. Otsuki and H. Matsukawa: Sci. Rep. 3

force at the bottom locally satisfies Amontons’(2013) 1586.

law with the local static and kinetic friction[2] Y. Katano. et al.: Sci. Rep. 4 (2014) 6324

coefficients 𝜇$ and 𝜇% , respectively.

258

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

Constructing a Dataset of Dielectric and Thermal Properties
of Polymeric Materials
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 dielectric and[3] and MD simulations are performed using

thermal properties of polymeric materials using

the LAMMPS software package [4]. All

quantum chemical (QC) calculations andcalculations are automated using RadonPy [5]molecular

and performed on System B (ohtaka).

dynamics(MD)

simulations.

Seventy-seven polymer systems are constructed
from

the

repeating

units

of

Figure 1(a) compares calculated andpreviously

experimental values of coefficients of linear

synthesized polymers [1,2], where atomic

thermal expansion, exhibiting a weakly linear

charges are derived from QC calculations using

relation

the Gaussian 16 electronic structure package

experimental values have large uncertainties.

betweenthem,

althoughtheCalculated static dielectric constants also show
reasonable agreements with the experimental
values except for several polymers, as shown in
Fig. 1(b). These results validate our simulation
methodology,

althoughmorethoroughinvestigation is required to further improve thepredictive performance.

References[1] L. Chen et al.: npj Comput. Mater. 6 (2020)
61.
[2] S. Otsuka et al.: Int. Conf. on Emerging
Intelligent Data and Web Technologies(2011) 22.
[3] M. J. Frisch et al.: Gaussian 16, Revision

Fig. 1: Comparison of calculated andC.01 (2019).

experimental values for coefficients of[4] A. P. Thompson et al.: Comput. Phys.

–4linear thermal expansion [10 /K] (a) andCommun. 271 (2022) 108171.

static dielectric constants (b).

[5] Y. Hayashi et al.: arXiv:2203.14090 (2022).

259

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

Nanoscale phonon transport across ceramics
interfaces
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

Lattice thermal conductivity is a fundamental property of materials, which is important
for a wide spectrum of applications, e.g., thermoelectrics, thermal barrier coatings, and electronic devices. Interfaces embedded in materials, such as grain boundaries and heterogeneous interfaces, are known to signiﬁcantly
suppress lattice thermal conduction becauseof their distinguished atomic/electronic structures from the corresponding crystal, which induce the scattering of group motion of atoms,
i.e., phonons. Recently, miniaturization ofelectronic devices and development of nanostructuring techniques have led to an increaseof number of cases in which the interfaces are
the rate-limiting factor for the thermal conduction in materials. For more precise control ofthermal conductivity, it is necessary to reveal
the mechanisms of thermal transport across interfaces from the microscopic point of view.
However, there are few studies that quantitatively reveal the correlations between interfacial atomic structures and thermal conduction,and that investigate the scattering mechanismsof phonons at interfaces. The main obstacles
are follows: (1) the lack of method to quantitatively connect the structure and thermal
conduction; and (2) the lack of method to evaluate thermal conductivity of local structure at
atomic-scale.
In this project, we have attempted to elucidate the microscopic mechanisms of inter-

facial thermal transport by using two kindsof molecular dynamics simulations. One is
the perturbed molecular dynamics simulations,in which atomic contributions to overall lattice thermal conduction (atomic thermal conductivity) can be estimated. The other isthe phonon wave packet simulations, in whichtransmission coeﬃcient of individual phonons
can be estimated. Combinations of these two
methods can provide the insights of interfacialphonon scattering from the phenomenologicaland theoretical perspectives. The necessary
codes were developed by the author and implemented to the Large-scale Atomic/Molecular
Massively Parallel Simulator (LAMMPS) [1].
We chose grain boundaries of crystalline silicons and interfaces between (1) crystalline silicon and amorphous silica (SiO2 ) and (2) crystalline and amorphous silicons as model systems, because it has many technologically important applications.
For silicon grain boundaries, we used
a machine learning potential for Si distributed in MACHINE LEARNING POTENTIAL REPOSITORY [2], which accurately reproduce the results of ab initio calculations but
are less computationally demanding. Thermal conductivity calculations using perturbed
molecular dynamics require a large computational cells with at least a few thousands ofatoms, with the time steps of more than one
million. Such large-scale calculations were per-260

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

For heterogenious interfaces between crystalline and amorphous silica, we have performed a large number of perturbed molecular dynamics and phonon wave packet simulations, by using an empirical Tersoﬀ potential,
which is less computationally demanding thanthe machine learning potential. The computational cells typically contains about a million of atoms for phonon wave packet simulations. It was found that the transmission coeﬃcients of acoustic phonons dramatically decrease as a function of the thickness of amorphous SiO2 layer, resulting in very low thermal
conductivity. The signiﬁcant amount of reﬂected phonons at the interface make the spatial depdendence of thermal conductivity evenin the crystalline silicon. Comparisons with
the interfaces between crystalline and amorphous silicons indicates that compositional dif-ferences between two materials enhance thescattering and reﬂection of phonons at the interfaces, because of the mismatch of phonon
frequencies between the two adjoined materials.
All the molecular dynamics calculations
were performed in a single node rather thanmultiple nodes. This is because it is necessary
to perform many calculations with diﬀerent
structures and diﬀerent conditions, and performing calculations on a single node is more
eﬃcient than parallel.

GB thermal resistance (10-9mK-1W-1)

formed with ﬁve diﬀerent magnitudes of perturbations for each grain boundary, to obtain the thermal conductivity with high statistical accuracy. The results indicates that
thermal conductivity of Si grain boundaries
highly depends on their microscopic atomic
structures as well as their crystal misorientations, as shown in Fig. 1 [3]. On the other
hand, phonon wave packet simulations showsthat acoustic phonons mainly conduct heat
across interfaces and the transmission coeﬃcients are not signiﬁcantly diﬀerent between
the grain boundaries with diﬀerent crystal orientation and/or atomic structures. This indicates that anharmonic nature of atomic vibrations at grain boundaries, which is naturally
included in the perturbed molecular dynamics but is intentionally excluded in the phonon
wave packet simulations, determines the interfacial thermal conduction. Additional calculations for 18 silicon grain boundaries suggestthat bond strains induced by the open-core
structures of silicon grain boundaries highly
correlate with the reduction of atomic thermal
conductivities, which should be related withthe magnitude of anhamonicity.

MLP

Tersoff

Σ5(310)/[001]Σ3(112)/[110]

1.5

_

1.0

0.5

0.0

Figure 1: Interfacial thermal resistances ofΣ5(310) and Σ3(112) grain boundaries estimated by perturbed molecular dynamics simulations, with the machine learning potentinal(MLP) and an empirical atomic potential (Tersoﬀ) [3]. The results of two diﬀerent atomic
sturctures are shown for each grain orientation.

References[1] S. Plimpton, J. Comput. Phys. 117 (1995)
1.
[2] A. Seko, arXiv:2007.14206.
[3] S. Fujii, A. Seko, Comput. Mater. Sci. 204(2022) 111137.

261

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

Numerical Study of Spin Systems on the
Honeycomb Lattice
Chitoshi YASUDA
Department of Physics and Earth Sciences, Faculty of Science,
University of the Ryukyus, Okinawa 903-0213, JapanThe quasi-two-dimensional (quasi-2D) antiferromagnetic Heisenberg (AFH) model ona cubic lattice consisting of the nearestneighbor intra-layer (J = 1) and inter-layer(J ′ ) interactions has a ﬁnite-temperaturephase transition at which the critical behavior belongs to the universality class sameas the three-dimensional classical Heisenberg
system. In the previous work [1] on the quasi2D square-lattice AFH model using the quantum Monte Carlo method, it was numericallyfound that the Néel temperature TN obeys anuniversal relation between J ′ and the staggered susceptibility χs (TN ) of the 2D systemon the square lattice, i.e., besides an eﬀective coordination number deﬁned by ζ(J ′ ) ≡1/(J ′ χs (TN (J ′ ))) does not depend on J ′ forJ ′ ≤ 0.02, the value becomes ζ(J ′ ) ≈ 0.65z(z = 2) all for S = 1/2, 1 and ∞. The
universal relation has also been predicted byan analysis of the quantum nonlinear sigma
model [2].

at the Néel temperature. Our results suggestthat ζ/z ≈ 0.65 in the small J ′ region. The
consistent value of ζ/z in the square and honeycomb lattices indicates that ζ is a universal quantity independent of the lattice. However, the accuracy of the result for the smallest J ′ , J ′ = 0.001, is not yet suﬃcient, andfurther calculations for larger sizes might beneeded in the future. Hastings et al. haveproposed that there exists a universal function even in the region of J ′ larger than J ′where ζ is constant [2]. Our results suggestthat there is one lattice-independent function
ζ(J ′ ) at J ′ = 0.05 and 0.1.
In this project, we also research zinc phthalocyanine (ZnC32 N8 H16 (ZnPc)) and theﬂuorinated ZnPc (ZnC32 N8 F16 (F16 ZnPc))
consisting of hexagon as with the honeycomb
lattice as testing calculations for the density functional method using the program
package VASP. The compound ZnPc is anaromatic compound composed of four isoindoline units (C8 H4 N) consisting of a sixmembered ring of benzene bonded to hydrogen and a ﬁve-membered ring containing nitrogen. We conﬁrmed that the energy
levels of F16 ZnPc are deeper than those of
ZnPc. The numerical result corresponds tothe electron-accepting property of F16 ZnPc
and the electron-releasing property of ZnPc.

In this project, we research the renormalized coordination number in the quasi2D honeycomb lattice AFH model with spin
S = 1/2 using the quantum Monte Carlo
method with the continuous imaginary-timeloop algorithm, and compare the results withthose of the square lattice to investigate theuniversality. The Néel temperatures were
evaluated by ﬁnite-size scaling analysis using data up to the system size 80 × 40 × 40,
and it was conﬁrmed that there is a relation
TN ∝ 1/(b − ln(J ′ )). Here, b is a constant.
We have also estimated the renormalized coordination numbers by evaluating the staggered susceptibility on the honeycomb lattice

References[1] C. Yasuda et al, Phys. Rev. Lett. 94(2005) 217201.
[2] M. B. Hastings and C. Mudry, Phys.
Rev. Lett. 96 (2006) 027215.

262

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

Data analysis method using a tensor network
representation
KENJI HARADA
Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan

We study data analysis methods based on
tensor network representations in this study
project.
The first topic is quantum machine learning. Using tensor network schemes, we numerically study the relationship between the performance of a quantum circuit and the entanglement properties. In particular, the quantum circuits on a tree tensor network structure can be efficiently estimated about theentanglement capability[1]. We found that aquantum circuit with a high entanglement capability shows good performance in a classification problem until some hundred qubit
cases. Therefore, the entanglement capability is an important property for estimating
performance[2]. In this study, we use the system B for calculation of a tensor contraction toestimate the entanglement capability throughtensor network representations.
The second topic is tensorized neuralnetworks[3]. A tensor network as MPO represents the weight matrix in a neural network.
We numerically evaluate the profile of an entanglement property in MPO in the learning
process. It much depends on the initial tensors
in MPO[4]. We derive a stable initialization
method for MPO in the tensorized neural network. The deeper tensorized neural network
can also be stably converged by the initialization. In this study, we use the system B for
the systematic calculation of the profile of entanglement in MPO in the learning process.

gling capability of parameterized quantum circuits for hybrid quantum-classicalalgorithms. Advanced Quantum Technologies 2, 1900070 (2019).
[2] H. Manabe and K. Harada, Quantum circuit learning with tensor network structure, 23pA1-1, JPS Autumn meeting(2021).
[3] A. Novikov, D. Podoprikhin, A. Osokin,
D. Vetrov, Tensorizing Neural Networks,
NIPS 2015.
[4] Y. Asoshina and K. Harada, Entanglement Analysis of Neural Networks,
22pL4-9, JPS Autumn meeting (2021).

[1] S. Sim, Peter D. Johnson, and Alán
Aspuru-Guzik, Expressibility and entan-

263

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

Numerical Analysis for Stability and Dynamics
of Magnetic Structures such as SkyrmionsKenji KONDO
Research Institute for Electronic Science, Hokkaido University,
Sapporo, Hokkaido, 001-0020.

Magnetic skyrmions are expected to be promising candidates for representing and transferring one
bit of information. This is because they are topologically protected from external perturbations, which
ensures the stability of information. Moreover,they are very small in size and they are able to be
driven using the current. In this study, we utilize
the spin dynamics code “Spirit” [1] which can incorporate the demagnetization ﬁelds or can calculate three-dimensional magnetic structures in orderto perform the numerical calculations of the dynamics of topological magnetic structures such asskyrmions. In our previous work, we investigatedthe dynamics of skyrmions driven by the spin polarized current [2,3] using a handmade ﬁnite-elementmethod (FEM) program. At that time, we did nottake the demagnetization ﬁelds into consideration
since we considered the inﬂuence to be negligible.
Also, there was an aspect that we avoided considering the demagnetization ﬁelds since it takes
much time to calculate them. Thus, in this study,
we investigate the inﬂuence on the velocities of theskyrmions with the demagnetization ﬁelds or without them using Spirit. Figure 1 shows the current
density dependences of the average velocities of theskyrmions in the ferromagnets (FMs) driven by the
spin-transfer torque (STT) with the demagnetization ﬁelds or without them. In this calculation, we
assume Co thin ﬁlm as the FM and use the material
parameters of the FM from the references [4–6].
The blue and the red lines represent the averagevelocities of the skyrmions in the FMs with andwithout the demagnetization ﬁelds, respectively. It

is found that the existence or the non-existence ofthe demagnetization ﬁelds have only 3.729% inﬂuences on the velocities of the skyrmions, whichis almost negligible. Our previous estimation of
the inﬂuences on the velocities of the skyrmions isfound to be right. On the other hand, we investigate the dynamics of three-dimensional topological magnetic structures, so-called skyrmion strings,
using Spirit. The skyrmions can form string structures by spreading in the out-of-plane direction.
Especially, we investigate the diﬀerences of theskyrmion string structures under the positive andnegative Dzyaloshinskii-Moriya interaction (DMI)
constants. The signs of DMI constants are deﬁned
in Fig. 2. When the DMI vectors point outwards,
we deﬁne the DMI constant as positive. When
the DMI vectors point inwards, we deﬁne the DMI
constant as negative. Figures 2(𝑎) and 2(𝑏) show
the skyrmion string structures in the FMs in thecases of the positive and negative DMI constants,
respectively. It is found that the directions of themagnetization vectors of the skyrmion strings arecounterclockwise and clockwise when seeing fromthe top view in the cases of the positive and negative DMI constants, respectively. Figures 3(𝑎)-3(𝑐)show the contour plots of the 𝑥, 𝑦, and 𝑧 componentsof the eﬀective magnetic ﬁelds 𝐻eﬀ in the case of
the positive DMI constant, respectively. Moreover,ﬁgures 3(𝑑)-3( 𝑓 ) show the contour plots of the 𝑥,𝑦, and 𝑧 components of 𝐻eﬀ in the case of the negative DMI constant, respectively. Comparing theseresults, it is found that the eﬀective magnetic ﬁelds
𝐻eﬀ point to the opposite directions in the 𝑥𝑦-plane.

264

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

Therefore, the positive and negative DMI constantscause the magnetization vectors of the skyrmion
strings to point to the opposite directions in the 𝑥𝑦plane. As a result, there manifest themselves thecounterclockwise skyrmion strings and the clockwise skyrmion strings depending on the signs of
DMI constants. It took approximately one monthto perform all of these calculations as the computer
run time using the ISSP supercomputer system B(Ohtaka). Also, we have calculated the dynamics
of the logic gates consisting of the skyrmions in
the antiferromagnets, which we will not state here
in detail. We are preparing the paper concerning
the logic gates consisting of the skyrmions in theantiferromagnets.
45
Without demagnetization field

40

With demagnetization field

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 3: The contour plots of the 𝑥, 𝑦, and 𝑧 components of 𝐻eﬀ of skyrmion strings are shown inﬁgures (𝑎)-(𝑐) in the case of the positive DMI constant, respectively. The contour plots of the 𝑥, 𝑦,
and 𝑧 components of 𝐻eﬀ of skyrmion strings areshown in ﬁgures (𝑑)-( 𝑓 ) in the case of the negative
DMI constant, respectively.

v (m/s)35

References30
25[1] G. P. Müller, M. Hoﬀmann, C. Dißelkamp,
D. Schürhoﬀ, S. Mavros, M. Sallermann, N.
S. Kiselev. H. Jónsson, and S. Blügel, Phys.
Rev. B. 99, 224414 (2019).

20
15
10
5
0.00

0.25

0.50

0.75

1.00

1.25

2

[2] Y. Ishida and K. Kondo, J. Magn. Magn.
Mater. 493, 165687 (2020).

Currentdensity (TA/m )Fig. 1: The current density dependences of theaverage velocities of the skyrmions in the FMs withthe demagnetization ﬁelds or without them.

(a) Positive DMI constantD(b) Negative DMI constantD[3] Y. Ishida and K. Kondo, Jpn. J. Appl. Phys.
59, SGGI04 (2020).
[4] J. Sampaio, V. Cros, S. Rohart, A. Thiaville,and A. Fert, Nat. Nanotechnology. 8, 839-844(2013).
[5] L. S. E. Alvarez, K.-Y. Wang, S. Lepadatu,
S. Landi, S. J. Bending, and C. H. Marrows,
Phys. Rev. Lett. 104, 137205 (2010).
[6] D. J. Monsma and S. S. P. Parkin, Appl. Phys.
Lett. 77, 883 (2000).

Fig. 2: The 3D vector plots of the skyrmion string
structures in the FMs in the cases of (𝑎) the positiveand (𝑏) the negative DMI constants.

265

Activity Report 2021 / 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

Many-body problems have to been tackled
in the field of the condensed matter physics;
therefore, it is diﬃcult to estimate physicalquantities precisely in many cases. Quantum
spin systems are such a typical one. To investigate these systems, numerical approaches
have widely and eﬀectively been used. A lot ofcomputational studies have been carried out;
useful information of the target systems has
been provided.
Generally speaking, three methods are eﬀectively used within the Field of quantum spin
systems. The first and traditional one is thenumerical diagonalization method; the second
is the quantum Monte Carlo (QMC) method.
The third method is the density matrix renormalization group (DMRG) method. Each ofthem has some advantages; on the other hand,
it has disadvantages, too. The QMC simulations can treat significantly large systems
irrespective of the spatial dimensions of thesystems; however, the negative sign problem
in this approach prevent us with a diﬃculty
in precise evaluation of physical quantities infrustrated systems. The DMRG method is
a very useful method for a one-dimensional
system irrespective of whether the target system includes frustration or not. However,this method is still under development for thecases when the spatial dimension is larger thanone. The numerical diagonalization method
can be applicable irrespective of the presenceof frustrations and the spatial dimension. Thismethod, on the other hand, has a severe weak

point that it can treat only very small system sizes. To overcome this disadvantageof the numerical diagonalization, we developed a hybrid-type parallelized code of Lanczos diagonalization[1]. This code enables us
to treat various large systems that have notbeen previously treated yet within this method
and to evaluate physical quantities precisely.
Therefore, we examine various quantum spin
systems by this method as a primary approachin this project.

Figure 1: Floret-pentagonal lattice.
In the project in 2021, we tackled the S =1/2 Heisenberg antiferromagnet on the floretpentagonal lattice[5]. Among pentagonal lattices, studies concerning the Cairo-pentagonallattice antiferromagnet has been known[2, 3,
4]. However, the floret-pentagonal-lattice an-266

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

tiferromagnet has not been investigated to thebest of our knowledge; this system is the second example as pentagonal lattices that include equivalent two directions in the system.
(See Fig.1.) We have observed the magnetization process of this system and find magnetization plateaux at one-ninth height of the saturation magnetization, at one-third height, andat seven-ninth height. We have also discovered a magnetization jump between the onethird-height and seven-ninth-height plateaux.
The jump is a peculiar phenomenon because
this jump is away from the two plateaux. Our
study contributes much to our deeper understandings of the frustrated quantum antiferromagnets. Further investigations would clarify
nontrivial quantum eﬀects due to frustration.

References[1] H.
NakanoandA.
Terai:
J. Phys. Soc. Jpn. 78, 014003 (2009).
[2] I. Rousochatzakis, A. M. Läuchli and
R. Moessner: Phys. Rev. B 85, 104415(2012).
[3] H. Nakano, M. Isoda, and T. Sakai:
J. Phys. Soc. Jpn. 83, 053702 (2014).
[4] M. Isoda, H. Nakano, and T. Sakai:
J. Phys. Soc. Jpn. 83, 084710 (2014).
[5] R. Furuchi, H. Nakano, N. Todoroki, and
T. Sakai: J. Phys. Commun. 5, 125008(2021).

267

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

Development of the Low-Temperature Oxidation Modelfor Alternative Fuels Using Machine Learning
Minhyeok LEE, Tomotaka MIZUNO, and Yuji SUZUKI
Department of Mechanical Engineering, School of Engineering,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656
To improve thermal efficiencies and reduce10

pollutant emissions of Internal Combustion
simultaneously,itiscrucial4

Number of reactions

engines5

tounderstand the detailed reaction mechanism offuels. Recently, various reaction kinetic models
have been proposed to precisely explain the

10

n-C7-12 (LLNL)
n-C8-16 (LLNL)
n-C7-12 (KUCRS)

4

n-C7 (NUI)
n-C7 (LLNL)
2 AramcoMech 2.0
3 AramcoMech 1.3
4

10

4
2

reaction process of the target fuel. However, as
10

shown in Fig. 1, as increasing the size of the fuelmolecule, the number of reactions involved isn-C1-16 (CRECK)

2n-Ci (KUCRS)USC-Mech II
GRI-Mech 3.0

2

0

2

4

6 8 10 12 14 16
Carbon number

Fig. 1: Reaction size of various kinetic models

increased exponentially. It makes hard to solve

for normal-alkane fuels.

the entire model numerically without reducing
the reaction size due to high computational costs.
1.5 x 10

particularly in the low-temperature region [1].
Therefore, we conducted numerical simulations
with full mechanisms by means of thesupercomputer system at ISSP to evaluate thenC103decrease in accuracy of the numerical estimation,

1.0[kmol/m ]HCHO molar concentration

Note that such model reduction causes a

-4

0.5
nC10

One-dimensional axisymmetric stagnation

KUCRS

nC7
0.0
600

ignition process of cool flames.

nC7

CRECK

650
700
750
800
Wall temperature, Tw [K]850

Fig. 2: Cool flame ignition of n-C7 and n-C10.

flows formed by impinging a mixture of fuel andair on a heated wall were assumed. The nozzle-

K to 850 K under the quasi-steady state. Normal-

to-wall distance was 15 mm. The outlet velocity,

heptane (n-C7H16) and -decane (n-C10H22) weretemperature, and equivalence ratio of the fuel/air

used as the fuel. Two reaction models including

mixture were 50 cm/s, 300 K, and 1, respectively.

both fuels were used, which respectively

The wall temperature was ramped up from 600

proposed by the CRECK modeling group [2]268

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

using

the

KUCRS

estimating not only the low- and intermediate-utilities

forcomplextemperature reactivity of each fuel, but also

reaction systems) [3]. The Cantera code [4] washierarchical difference with respect to theemployed to conduct numerical simulation

molecule (chain) size.

andgeneratedby(Knowledge-basing

Further improvements in the reaction model

under the above condition.

should be made through the experimental

As increasing the wall temperature, coolvalidation.

flames are ignited on the wall due to initiation of
the low-temperature oxidation. Figure 2 showsincreases of the calculated HCHO molar

References

concentration on the wall as a function of the[1] C.B. Reuter, M. Lee, S.H. Won, Y. Ju:

wall

Combust. Flame 179 (2017) 23.

temperature.

While

the

HCHO

concentration estimated by the CRECK model is[2] A. Miyoshi: KUCRS software library,higher, the onset temperature calculated by the

http://akrmys.com/KUCRS/KUCRS model is slightly lower for both fuels.

[3]

Furthermore,

http://creckmodeling.chem.polimi.itadifference

in

the

onset

The

CRECK

Modeling

Group:

temperature between two fuel is negligible for[4] D.G. Goodwin, R.L. Speth, H.K. Moffat,

the KUCRS model, but n-decane’s onset

B.W. Weber: Cantera: An object-oriented

temperature is approximately 12 K lower thansoftwaretoolkitforchemicalkinetics,

that of n-heptane in the CRECK model. It means

thermodynamics,

andtransport

processes,that two models show the discrepancy inhttps://www.cantera.org

269

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

Formation of Lipid Rafts Studied by Molecular Dynamics
Simulation
Satoru G. ITOH
Institute for Molecular Science, Okazaki, Aichi 444-8787
Exploratory Research Center on Life and Living Systems, Okazaki, Aichi 444-8787
The Graduate University for Advanced Studies, Okazaki, Aichi 444-8787

In cell membranes, lipid rafts are formed by

perform MD simulations with various force

aggregation of sphingolipids and cholesterolfields for the all-atom model. We have also

molecules [1]. It has been clarified that lipid

made it possible to perform MD simulations onrafts play important roles in functions such as

molecules other than proteins. In addition, thesignal transduction and protein transport in vivo

program has been modified to allow the[2,3]. Lipid rafts are also related to virus

Coulomb replica-permutation method to scale

infection and protein aggregation [4,5]. It is

onlyessential to elucidate their structures andinteractions

formation

molecules.

processes

tounderstand

their

the

intermolecularelectrostatic

betweenganglioside

GM1

functions. Therefore, we aim to elucidate theseby using molecular dynamics (MD) simulation.

References

In particular, we are interested in lipid rafts[1] K. Simons and E. Ikonen: Nature 387

containing GM1 gangliosides, which are(1997) 569.

abundant[2] K. Simons and D. Toomre: Nat. Rev. Mol.

in

nerve

cell

membranes.

To

investigate the structure of GM1 ganglioside

Cell Biol. 1, (2000) 31.

clusters, we employed the Coulomb replica-[3] D. Lingwood and K. Simons: Sience

permutation molecular dynamics simulation [6].

327 (2010) 46.

The Coulomb replica-permutation method is[4] N. Chazal and D. Gerlier: Microbiol. Mol.

one of the Hamiltonian replica-permutation

Biol. Rev. 67 (2003) 226.

methods, and it is possible to efficiently[5] J. M. Cordy, N. M. Hooper and A. J. Turner:

investigate

Mol. Membr. Biol. 23 (2006) 111.

thecomplexstructureofbiomolecules by using this method.

[6] S. G. Itoh and H. Okumura: J. Comput.

In the last year, we have made significant

Chem. 34 (2013) 2493.

improvements to our own program to be able to270

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

Numerical calculation of spin Hall magnetoresistanceby the quantum Monte Carlo method
Takeo KATO
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

Spin Hall magnetoresistance (SMR), thatis a magnetoresistance caused by the absorption of spin current generated by the spin
Hall eﬀect in a bilayer system composed ofa metal and a magnetic insulator [1], isbeing actively studied for potential applications in spintronics devices such as memory applications. Although SMR is well described quantitatively by a phenomenologicaltheory [2], its temperature dependence for a
metal/antiferromagnetic-insulator bilayer system [3] are diﬃcult to be described by it.
In this numerical study, we evaluate the spin
conductance at the interface by expanding thetheory of Ref. [4]. The spin current is expressed in terms of the spin susceptibility of
the antiferromagnetic insulator and is calculated by the quantum Monte Carlo method. Inorder to obtain dynamic properites, we need toperform the numerical analytic continuation.
However, this is usually an ill-deﬁned probpem
and is diﬃcult to be applied for obtaining thespin conductance. We improved the method
to stably obtain the spin conductance using
only the component of Matsubara frequency ofGreen’s function and succeeded in its accurate
calculation.
As a result of examining the temperaturedependence and the thickness dependence of
the SMR, the sign was reversed as the risingtemperature rose in the negative sign from thelow temperature region, and the result agreed
with the experimental result with the peak of
the positive sign. In addition, it was conﬁrmedthat the behavior of SMR diﬀered dependingon the size of the spin of the antiferromagnetic
insulator, and that SMR at S = 1 was larger
than S = 1/2. We also showed that SMR
vanishes above the magnetic transition temperature unlike the experiment. We are now
preparing for a paper on these results [5].

Figure 1: SMR Measurement Setup.

References[1] H. Nakayama, et al., Phys. Rev. Lett.
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 in
preparation.

271

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

Dynamical property of a resistively shunted
Josephson junction
Tsuyoshi Yamamoto
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

Physical properties of quantum systems are
drastically changed by interaction with anenvironment, because thousand of modes inthe environment are complexly entangled withquantum states of systems. A typical example is the Kondo problem, in which a spin ofa magnetic impurity is screened by surrounding conduction electrons in a normal metal [1].
Such kind of a dissipative system has been acentral topic of condensed matter physics for along time.
One of the most notable problems in dissipative systems is a dissipation-driven quantum phase transition. In 1983, Schmid predicted that a single Josephson junction connected to a resistance, which is called the resistively shunted Josephson junction (RSJJ),
undergoes a superconductor-insulator transition at zero temperature [2]. When the resistance is larger than a resistance quantum(R > RQ = h/4e2 ), the Josephson junction
becomes insulating, while when R < RQ , it becomes superconducting (see Fig. 1). However,there is no clear experimental evidence to haveobserved the transition, and it is still controversial [3]. The recent development of superconducting circuits allows us to access the dynamical property of the RSJJ using microwave
spectroscopy while it has been discussed in thedc limit. This development motivates us to approach the RSJJ from the dynamical response,and it provides a new perspective to solve thecontroversial problem, the Schmid transition.
To this end, we calculated frequency depen-

Figure 1: (left) A Josephson junction with
the Josephson energy EJ and the charging energy EC connected to a resistance R. (right)
Schmid transition diagram. In the phase I (S),the Josephson junction becomes an insulator(a superconductor) at zero temperature.

dence of an impedance of the RSJJ circuit andexamined dynamic properties in the microwave
reﬂection via the RSJJ circuit. We ﬁrst calculated the phase-phase correlation function
deﬁned on the imaginary-time axis using thepath-integral Monte Carlo simulation [4] andthen performed the analytical continuation to
the phase-phase correlation function numerically to obtain the dynamic impedance of the
RSJJ. We checked that our numerical results
agree well with analytical results at some solvable points such as the Toulouse point and thefree-fermion point. This work is expected toprovide a deep understanding of the Josephson junction on environments, which pushes
further development for the superconducting
circuit using Josephson junctions.

272

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

References[1] J. Kondo, Prog. Theor. Phys. 32, 37(1964).
[2] A. Schmid, Phys. Rev. Lett. 51, 1506
(1983).
[3] A. Murani et al., Phys. Rev. X 10, 021003(2020).
[4] P. Werner and M. Troyer, Phys. Rev. Lett.
95, 060201 (2005).

273

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

Elucidation of Pathways for the Crystallization of Ionic
Liquids by Large-Scale Metadynamics Simulation
Hiroki NADA
Environmental Management Research Institute,
National Institute of Advanced Industrial Science and Technology (AIST), Onogawa 16-1,
Tsukuba 305-8569

A classical molecular dynamics (MD)predicted with the FES, did not show any

simulation combined with the metadynamics

energetic barrier.

(MTD) method (hereafter, MTD simulation)

In addition, a crystal-like structure, which

was performed for an ionic liquid, 1-alkyl-3-

was different from the ordered structures for the

methylimidazolium

shallow local minima, also appeared at a highhexafluorophosphates([Cnmin]PF6, n=1 [1]), at 298 K in the NVT

freeensemble. A force field proposed by Canongia

arrangement of ions in this crystal-like structure

Lopes et al. was used for estimation of thewas different from those for the crystalline

interionic interaction [2]. The simulation

forms, which were reported in literature [1].

systemrectangular-parallelepiped

This crystal-like structure might correspond to

consisting of 108 ion pairs. The MTD

a high-pressure phase, which has not yet been

simulation was carried out using DL_POLY

found experimentally.

wasaenergyregion

on

the

FES.

The

2.20 [3], in which PLUMED 1.3 [4] wasimplemented to permit combination with theReferences

MTD method.

[1] T. Endo, H. Masu, K. Fujii, T. Morita, H.

A free energy surface (FES) obtained using

seki, S. Sen and K. Nishikawa: Cryst. Growth

discrete C-C and P-F radial distribution& Des. 13 (2013) 5383.

functions quantified by Gaussian window[2] J. N. Canongia Lopez, J. Deschamps and A.

functions as collective variables indicated a few

A. H. Pádua: J. Phys. Chem. B 108 (2004) 2038.

shallow local minima. The structure of the[3] W. Smith and T. R. Forester: J. Mol. Graph.

system at each shallow local minimum

14 (1996) 136.

represented an ordered arrangement of ions.

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

The pathway for the formation of this ordered

180 (2009) 1961.

structure from the bulk ionic liquid, which was274

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

Creation of a Method for Design of Scale Formation Control
Molecules by a Metadynamics Method
Hiroki NADA
Environmental Management Research Institute,
National Institute of Advanced Industrial Science and Technology (AIST), Onogawa 16-1,
Tsukuba 305-8569

Classical

molecular

dynamics

more pronounced effect than Na+

(MD)

ions.

simulations combined with the metadynamics

Furthermore, the simulation results suggested(MTD) method [1] (hereafter, MTD simulation)that the binding stability at the calcite surface

were conducted to elucidate the stable binding

was higher for the PMA additive than for theconformations of a deprotonated polymaleic

PAA additives, and the PAA additive with aacid (PMA) additive and two deprotonatedpolymerization degree of 10 displayed apolyacrylic acid (PAA) additives with different

binding stability that was similar to or lower

polymerization degrees in the presence ofthan that of the PAA additive with avarious countercations at a hydrated CaCO3

polymerization degree of 5.

calcite (104) surface. The MTD simulations

The present simulation method provides awere performed using DL_POLY 2.20 [2] innewstrategy for

analyzingthebinding

which PLUMED 1.3 [3] was implemented.

conformations of complex additives at material

The simulated free-energy surfaces (FES)

surfaces, developing additives that stably bindsuggested the existence of several slightly

to these surfaces, and designing additives todifferent stable binding conformations for eachcontrol crystal growth.

additive. The appearance of these distinct
binding

conformationsis

speculated

Referencestooriginate from different balances of the[1] A. Laio and M. Parrinello: Proc. Natl. Acad.

interactions between the additive, calcite

Sci. USA 99 (2002) 12562.

surface, and countercations. The binding[2] W. Smith and T. R. Forester: J. Mol. Graph.

conformations and binding stabilities at the

14 (1996) 136.

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

surfacewereaffectedbythecountercations, with Ca2+ ions displaying a

180 (2009) 1961.

275

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

Stability of the hedgehog-lattice topological spin
texture in breathing-pyrochlore 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 magnetic skyrmion in two dimensions and a magnetic hedgehog in three dimensions are typical examples of such topological objects. In
magnetic materials, they often appear in theform of crystal orders, so-called skyrmion crystal (SkX) and hedgehog lattice (HL), where
the Dzaloshinskii-Moriya (DM) interaction is
known to be essential for their stabilities. Concerning the SkX, it has been shown that magnetic frustration serves as another stabilization
mechanism [1], but as for the HL, a counter
ordering mechanism has not been reported sofar. In view of such a situation, we havesearched for a new mechanism other than the
DM interaction for the HL. Previously, we theoretically demonstrated that the HL is realized in breathing-pyrochlore antiferromagnets
in the absence of the DM interaction, where
the nearest-neighbor (NN) exchange interactions for small and large tetrahedra, J1 andJ1′ , are assumed to be antiferromagnetic [2].
In this work, we extend our previous theoretical work, examining the stability of the HL ina wider parameter space.
We consider the J1 -J3 classical Heisenberg
model on the breathing pyrochlore latticewhich is given by

H = J1∑⟨i,j⟩S

Si · Sj + J1′∑⟨i,j⟩LSi · Sj+J3∑

Si · Sj(1)⟨⟨i,j⟩⟩

where J1 and J1′ take diﬀerent positive or
negative values due to the breathing bondalternation and J3 is the third NN antiferromagnetic interaction along the bond direction.
In the model, the breathing bond-alternationof the lattice is characterized by J1′ /J1 , andJ3 is essential for the occurrence of the HL
which in the present case, is a quadruple-Q
state characterized by the ordering vectors of(± 21 , ± 12 , ± 12 ). By performing Monte Carlosimulations, we investigate the ordering properties of the Hamiltonian (1). In our simulations, 2×105 MC sweeps are carried out under the periodic boundary condition and theﬁrst half is discarded for thermalization. Our
1 MC sweep consists of 1 heatbath sweep andsuccessive 10 overrelaxation sweeps, and observations are done at every MC sweep. The
statistical average is taken over 4 independentruns starting from diﬀerent random initial conﬁgurations.
It is found that the quadruple-Q ( 12 , 12 , 12 ) HL
is stable irrespective of the signs of J1 and J1′
and that in an applied magnetic ﬁeld, there appear six quadruple-Q states depending on thevalues of J1 and J1′ , among which three phases
including the in-ﬁeld hedgehog-lattice state exhibit nonzero total chirality [3]. We have also
checked that the HL is robust against relatively
weak single-ion magnetic anisotropies.
In addition to the three-dimensional system,
we investigate an associated two-dimensional

276

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyosystem. We have demonstrated that the J1 J3 Heisenberg model on the breathing kagome
lattice exhibits a chiral transition and that the
low-temperature chiral state can be viewed as aminiature SkX [4]. In contrast to conventional
SkX appearing in an applied magnetic ﬁeld,
the miniature SkX is realized at zero ﬁeld, suggesting that a zero-ﬁeld topological Hall eﬀect
is possible in a metallic system.

References[1] T. Okubo, S. Chung, and H. Kawamura,Phys. Rev. Lett. 108, (2012) 017206.
[2] K. Aoyama and H. Kawamura, Phys. Rev.
B 103 (2021) 014406.
[3] K. Aoyama and H. kawamura, in preparation.
[4] K. Aoyama and H. Kawamura, Phys. Rev.
B 105, (2022) L100407.

277

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

Dimensional reduction eﬀect in the two-dimensional maple-leaf lattice
Heisenberg antiferromagnet
Chisa HOTTA
Department of Basic Science, University of Tokyo
3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan

We aim to clarify the features of uniform
magnetic susceptibilities of a series of frustrated magnets in two dimensions. This year,
we focused on the maple-leaf lattice spin1/2 Heisenberg model, whose lattice structure is regarded as a 1/7-site-depleted triangular lattice. There is recently experimental
evidence that the material called bluebellite
Cu6 I6 O3 (OH)10 Cl shows a particular magnetic
susceptibility that almost perfectly matches
the Bonner-Fisher curve, namely the susceptibility curve of the one-dimensional (1D) spin1/2 antiferromagnetic Heisenberg model[1].
Motivated by this experiment, we analyzed
the ground state and the ﬁnite temperature
properties of the maple leaf lattice Heisenberg
model. By combining the analysis of the several energy eigen states of the model on aﬁnite cluster using the exact diagonalization(ED), we found that the constituent of thelow energy eigenstates have the UUD/DDU
conﬁgurations, i.e. the spins on all the triangular units of the lattice are two-up-onedown or two-down-one-up. This situation issimilar to the antiferromagnetic triangular lattice Ising model, where the UUD/DDU state
contributes macroscopically to the residual entropy. However, in our model, not all the interactions are antiferromagnetic, but the interactions on bonds forming a hexagon aroundthe depleted sites are staggard in their sign.
Even in such a case, there is a substantial
frustration eﬀect, and UUD/DDU are favored.
Due to the quantum ﬂuctuation eﬀect, these
UUD/DDU manifold mixes, and a stripe-typeground state appears. The temperature dependence of the susceptibility turned out not tomatch the Bonner-Fisher curve even though wevary the sign and amplitudes of ﬁve diﬀerent
Heisenberg interaction parameters. Instead,
we found that the susceptibility shows an almost perfect coincidence with that of the 1D
XXZ model having a ﬁnite spin gap. The reason for this dimensional reduction behavior isthat the magnetic excitation from the orderedground state resembles the spinon excitationof the 1D XXZ model above the spin gap. Indeed, ﬂipping one of the spins of the stripge
ground state, the two plaquettes(two+two triangles) with UUUD spins are formed, each carrying spin-1/2. The Heisenberg exchange interactions separate these two plaquettes into
opposite directions, while such spinon propagation takes place only in a particular spatial
direction. A series of calculations on the susceptibility of both the maple leaf lattice andthe 1D XXZ model was done using a TPQ
method combined with SSD[2] The variation ofparameters we needed to search for was large
and the 2D susceptibility is maximally given
for N = 24, 30 sites[3].

References[1] Y. Haraguchi, et.al.: Phys. Rev. B 104,
174439 (2021).
[2] C Hotta and K Asano: Phys. Rev. B 98,
140405 (2018).
[3] R. Makuta and C. Hotta: Phys. Rev. B
104, 224415 (2021).

278

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

The eﬀect of a quantum phase transition on quantum
approximate optimization algorithm
Tatsuhiko Shirai
Department of Computer Science and Communications Engineering, Waseda University
Ookubo, Shinjuku-ku, Tokyo 169-8555−iβp H −iγp X⃗ ⃗γ ⟩ := e|β,−iβ1 H −iγ1 X

|+⟩ ,
(1)
⃗
where β = (β1 , . . . , βp ) and ⃗γ = (γ1 , . . . , γp ).
p denotes a circuit depth. X is called a∑N xmixer and adopted as X =i=1 σi whereyxz⃗σi = (σi , σi , σi ) is the Pauli operator acting
on site i. |+⟩ is deﬁned as σix |+⟩ = |+⟩ for all⃗ ⃗γ ) toi. The classical computer optimizes (β,⃗⃗minimize the value of ⟨β, ⃗γ | H |β, ⃗γ ⟩. Here, we
adopt the transverse ﬁeld Ising model given by
H = −g

e

N
∑
i=1

···e

z
σiz σi+1
−

e

N
∑

σix ,

(2)

10-1

10

6E/N

We have studied the eﬀect of a quantum phase transition on the performance ofquantum approximate optimization algorithm(QAOA) in the transverse-ﬁeld Ising model.
QAOA is a hybrid quantum-classical algorithm
designed to approximate the ground state ofa Hamiltonian H with N spins [1]. It has
a wide application in chemistry, physics, andcomputer science.
QAOA consists of calculations in a quantum
computer and classical computer. The quantum computer calculates the expectation valueof H over a variational state, which is given by

-2

10-3

10-4

10-5 0
10

g=0.8
g=0.9
g=1
g=1.1
g=1.2
slope -1
slope -2101Depth pFigure 1: p-dependences of the residual energies ∆E⃗ ⃗γ ) and the ground state)
timum values of (β,in the ferromagnetic phase, critical point, andparamagnetic phase in the thermodynamic
limit. The numerics is based on the Jordan–
Wigner transformation and time-dependent
Bogoliubov transformation. Figure 1 showsthat the asymptotic behavior discontinuously
changes at the critical point. It indicates theintimate relation between the quantum phasetransition and the performance of QAOA.

i=1

where the periodic boundary condition is imposed (i.e., ⃗σi+N = ⃗σi ). This model shows aquantum phase transition at g = gc = 1. The
ground state is ferromagnetic for g > gc , while
it is paramagnetic for 0 ≤ g < gc .
We numerically study p-dependences of theresidual energies ∆E (i.e., the energy diﬀerence between the variational state at the op-

References[1] E. Fahri, et al., arXiv:1411.4028 (2014).

279

Activity Report 2021 / 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
Our numerical study is based on the

Motivated by a recent experimental study[1] demonstrating the structural transformation

Landau-de

Gennescontinuum

theoryof a cholesteric blue phase liquid crystal

describing the orientational order by a second-(BPLC) by an electric field, we carried out arank tensor. We use the program developed insimulation study on the behavior of BPLC

Ref. [2] where we were particularly interested

under an electric field. Ref. [1] showed that ain the surface behavior, and dealt withstructure with tetragonal (I4122) symmetry,

relatively thin systems. In the present study, towhich had been known to be stable under ansimulate a system whose thickness is close to aelectric field, can be stable even after therealistic experimental value (≃12μm in Ref.

cessation of the electric field when the temporal[1]), 1,351 grid points are taken in the thickness

variation of the applied electric field isdirection. Our numerical system is on a regularappropriately controlled. Possibility of stable

grid and allows a use of OpenMP for efficient(or metastable) tetragonal structures of BPLC

calculations.
Our calculations demonstrate that theunder no electric field has never been discussed

presence of confining boundaries discretize thetheoretically.
Previous numerical studies on BPLCs under

number of periodic structures along thean electric field focused on the bulk behavior

thickness direction, and that the tetragonal(to be precise, in an infinite system withstructure is indeed metastable under no electric

periodic boundaries). We are interested in how

field. We are carrying out systematic studies on

the planar boundaries confining the BPLC will

the free energy versus the number of periodic

affect its structural stability. Elucidation of thestructures, and on the optical properties of the

dynamical behavior under a temporally varying

BPLC structures we find numerically.

electric field turns out to be a daunting task

Referenceseven with the aid of the ISSP supercomputer,[1] D.-Y. Guo et al., Nature Mater. 19 (2020) 94.

and therefore we restrict ourselves to static[2] J. Fukuda and S. Žumer, Phys. Rev.

behavior.

Research 2 (2020) 033407.

280

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

Non-equilibrium response and phase transition in thedense 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]. In thisproject, we focus on the equilibration and thenon-equilibrium response interfered with thephase transition in the hard disk/sphere model
system with modern algorithms, i.e., EventChain Monte Carlo (ECMC) [2, 3], Newtonian
ECMC (NEC) [4] and Event-Driven Molecular
Dynamics(EDMD) [5].
Equilibration and relaxation times in dense
molecular systems
The equilibration and relaxation times of thephysical properties have historically been fundamental issues in molecular simulation. As
practical tasks in the molecular simulation,‘equilibration’ (relaxation toward equilibrium
states) from initial non-equilibrium states isjust preprocessing step before a production
run. Such a step needs to obtain the physicalproperties in the equilibrium state in the specific system. However, the computational costoften becomes high, especially for high-density
and larger scale particle systems, and may
sometimes exceed the human lifetime. In general, the equilibration of particle positions indense systems is much more diﬃcult due to the
excluded volume eﬀect being dominant, which
is why the crystallization and glass/jamming
transition problems remain one of the critical issues and have been actively studied recently. The two standard molecular simulation methods for many-body systems are of-ten used; one is molecular dynamics (MD)using Newton’s equations of motion, and theother is Markov chain Monte Carlo (MC) using a pseudo-random number generator. Various novel algorithms focusing on equilibration
based on MD and MC have been implemented
in the history of molecular simulation, whichenables to reduce the relaxing time. An eﬃcient algorithm for the acceleration of equilibration, termed the ‘event-chain Monte Carlo(ECMC)’ [2] was the breakthrough on this
topic. The ECMC has become one of the milestones, and its variants have been investigatedactively over the past decade [3]. To obtainthe ‘true’ equilibrium states, it is a reasonablechoice to use ECMC for positional relaxation
at first, which is the idea of a hybrid scheme [1].
Recently, an ECMC variant termed the Newtonian ECMC (NEC) was proposed [4]. In NEC,the velocity collision rule was added to dictate
the direction of particle movement. This isoften used in EDMD of hard-sphere systems.
The NEC clearly improves the eﬃciencies ofthe diﬀusion coeﬃcient, nucleation rate, andmelting process. Recently, this eﬃcient algorithm has been extended to anisotropic hard
polygons without any approximations [6].
In these ECMC and its variants (includingNEC), the systems are driven by the sequential
collisions composed of a dozen of particles thatseem to behave like a chain. The eﬃciency depends on chain length, physical properties, andsystem size. However, the microscopic mechanism of equilibration for each parameter re-281

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

mains elusive. To identify key factors aﬀecting equilibration in the hard disk system, we
focus on the event-based algorithms [2, 4, 5]
described above and investigate the comparison between them in cases of (i) the diﬀusional
characteristics per event (collision) and (ii) the
relaxing process toward the liquid states as a
non-equilibrium response after disturbance ofthe homogeneous expansion.
(i) Diﬀusional characteristics of the eventbased algorithms in hard disk systems
To survey the key factors explaining why
the eﬃciency of equilibration diﬀered between
the event-based algorithms, we focus on thediﬀusional characteristics per event (collision)
between event-based algorithms in the liquidstates of the hard disk systems. We defined
the self- and pair-diﬀusion coeﬃcient based on
the displacement number of disks and compared the dependence in terms of the chain
length and the system size systematically. In
ECMC and NEC, as the chain length increases,
the diﬀusion coeﬃcient grows and reaches thehighest values at a specific chain length. We
visualized the spatial distribution of the accumulated number of displacements per particle
and calculated its probability density distribution. We found that NEC has an advantage
over ECMC in the diﬀusion coeﬃcient in thissystem since NEC samples homogeneously inthe system due to the Maxwell Boltzmann distribution of disk velocities for each displacement during the event chain, which can sample
disks uniformly in the liquid state [7].
(ii) Non-equilibrium response and equilibrationin hard disk melting
We then focus on the non-equilibrium response during the equilibration process induced by a disturbance of the homogeneous
expansion. In the case of EDMD, we observed
the anomalous slow equilibration in the harddisk systems around the co-existence phases
toward the liquid states in large-scale simulation via a couple of physical properties. As
the physical mechanism of this anomalous slow

equilibration, we assumed that the spatial inhomogeneity of the initial equilibrated phases
would contribute to the relaxation time significantly. The assumption was confirmed by theprobability distribution of local density andorientational order parameter [8]. We also observed anomalous slow equilibration both in
ECMC and NEC; however, the distinct diﬀerences between methods in the eﬃciency andfunctional forms of relaxation of physical properties are detected. We found that NEC and
NEC are more eﬃcient than EDMD.
As a typical non-equilibrium response, we
also investigated the thermodynamic cycles bymodeling the Stirling engine (SE) at the nanoscale molecular level [9]. We established a numerical model of a β-type SE with a displacer
piston and figured out the minimal necessary
conditions for stable rotation, as well as thelower limits of both the diﬀerence in temperature and the particle number.

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 012011 (2022).
[8] D. Mugita and M. Isobe: EPJ Web Conferences, 249 14004 (2021).
[9] K. Kitaya and M. Isobe: J. of Phys.:Conf. Ser. 2207 012006 (2022).

282

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

Structural Formation of Patchy Particles
Takamichi TERAO
Department of Electrical, Electronic and Computer Engineering,
Gifu University, Yanagido 1-1, Gifu, 501-1132
Patchy particles are spheres in which theeffective

interaction

betweenparticlesisanisotropic. Among the patched particles, there
is interest in those that have anisotropy ininterparticle interactions due to inhomogeneous
surface charges [1]. Charged patchy particles
are expected to have a new structure, different
from previous patchy particle systems, because
Fig. 1: Modeling of a charged patchy particle

the interaction between particles can be both
attractive and repulsive.
In this study, we focused on the structuralformation of charged patchy particles, in which
the particle surface is charged and the sign ofthe surface charge is reversed between the
bipolar part and the other parts. Cluster
formation in these systems was studied bymolecular dynamics simulations. Modeling of
charged Patchy particles is shown in Figure 1.
Figure 2 shows an example of the structure

Fig. 2: Snapshot of charged patchy particles

formation of the system when the motion of thecharged patchy particles was constrained
between a pair of two-dimensional planes. In

References

thermalpatchy[1] P. D. J. van Oostrum, M. Hejazifar, C.

particles exhibit a variety of configurations,

Niedermayer, and E. Reimhult: J. Phys.:including triangular, square, pentagonal, andCondens. Matter 27 (2015) 234105.

equilibrium,thechargedhexagonal lattices.

283

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

Tensor-network study of real-time dynamics in
two-dimensional quantum many-body systems
Ryui KANEKO
Department of Physics, Kindai University
Higashi-Osaka, Osaka 577-8502, JapanThe nonequilibrium physics of quantum
many-body systems have attracted much interest in recent years. Among others, real-time
dynamics of isolated quantum systems has
been investigated intensively in experiments
very recently. Such examples include analog
quantum simulators using ultracold atoms in
optical lattices [1–3] and those using Rydberg
atom arrays [4, 5]. The former is suitable forstudying the Bose-Hubbard models in 1D, 2D,
and 3D, and the latter is often used for studying the quantum Ising systems. The number ofcontrollable Rydberg atoms has now exceeded
200 [6–8], which would be more qubits thanclassical computers can handle. These recent
developments of analog quantum simulators
demand cross-validating the accuracy of experiments and numerical simulations. Althoughnumerical methods based on matrix productstates are powerful in 1D [1, 3, 9], no decisive
methods exist for simulating dynamics in 2D
systems.
In the present study, we examine the applicability of the 2D tensor-network method
based on the infinite projected entangled pairstate (iPEPS) [10, 11] for simulating the realtime dynamics. We first focus on the 2D BoseHubbard model, whose time-dependent singleparticle correlation functions have been measured in the recent experiment [3], and investigate sudden and short-time quench dynamics starting from the Mott insulator [12]. We
calculate the real-time evolved states by thesimple update algorithm. The tensor-network

Figure 1: Comparison of the time-dependent
single-particle correlation functions between
the experiment (red squares) and the iPEPS
simulations (blue line with symbols). We consider a finite-time quench and optimize thewave functions up to the bond dimensions
D = 8.
library TeNeS [13] is adopted.
We show the calculated single-particle correlation functions for the short-time quench to
the Mott insulating region in Fig. 1. They
are in good agreement with the recent experiment. We also predict how group and phasevelocities behave for the sudden quench to theweaker interaction region, which has not been
explored in experiments. Our findings would
serve as a quantitative benchmark for future
experiments.
After checking the applicability of the
iPEPS method to the 2D Bose-Hubbard
model, we apply the method to the 2D
transverse-field Ising model. We focus on asudden quench starting from the fully disordered state (Γ/J → ∞) to the disordered region (Γ > Γc ). Calculated longitudinal spin-

284

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

spin correlation functions agree very well withthe results obtained by the small-size exact diagonalization method and by the spin-wave approximation. We extract the group velocity
from the correlation functions in the stronger
field region. The group velocity is known
to be bounded above by the Lieb-Robinson
bound [14, 15]. The estimated group velocity in the present iPEPS study is much slower
than the recently updated bound [16], suggesting that there is room for improving the analytical bound.

[10] T. Nishino et al., Prog. Theor. Phys. 105,409 (2001).

References[11] F. Verstraete and J. I. Cirac, arXiv:condmat/0407066.

[1] M. Cheneau et al., Nature 481, 484(2012).

[6] S. Ebadi et al., Nature 595, 227 (2021).
[7] P. Scholl et al., Nature 595, 233 (2021).
[8] D. Bluvstein et al., Science 371, 1355
(2021).
[9] P. Barmettler et al., Phys. Rev. A 85,
053625 (2012).

[12] R. Kaneko and I. Danshita, Commun.
Phys. 5, 65 (2022).

[2] P. Jurcevic et al., Nature 511, 202 (2014).
[3] Y. Takasu et al., Sci. Adv. 6, eaba9255(2020).

[13] Y. Motoyama et al., arXiv:2112.13184.
[14] E. H. Lieb and D. W. Robinson, Commun.
Math. Phys. 28, 251 (1972)[4] E. Guardado-Sanche et al., Phys. Rev. X
8, 021069 (2018).

[15] M. B. Hastings, arXiv:1008.5137.

[5] V. Lienhard et al., Phys. Rev. X 8, 021070(2018).

[16] Z. Wang and K. R. Hazzard, PRX Quantum 1, 010303 (2020).

285

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

Tensor-network study of the SU(4) Heisenberg model
on a tetramerized square lattice
Ryui KANEKO
Department of Physics, Kindai University
Higashi-Osaka, Osaka 577-8502, JapanThe SU(N ) Hubbard and Heisenberg models
have attracted much attention recently [1–4].
They can be realized in the ultracold atoms
in optical lattices and in specific antiferromagnets. Since these models possess many internal degrees of freedom and would realize novel
quantum states that are absent in the wellknown SU(2) model, it is desirable to explore
ground states of the models on various lattices.
In the present study, we focus on how multimerization affects the ground states of the
SU(N ) Heisenberg model. To this end, we consider the SU(4) antiferromagnetic Heisenberg
model [5]X

H0 = J

Pij + J 0

hijiintra

X

Pij

X X

+J

X X

Sαβ (i)Sβα (j)Sαβ (i)Sβα (j).

Figure 1: Tetramerized square lattice. WhenJ  J 0 or J  J 0 , the system decouples into
isolated four-site chains. When J = J 0 , thesystem becomes the conventional square lattice, whose ground state is believed to exhibit
simultaneous dimerization and SU(4) symmetry breaking [5].
vors 0 and 1:

hijiinterhijiintra αβ
0

J0

(1)

on a tetramerized square lattice [6] (see Fig. 1),
where Pij is a transposition operator which
exchanges flavors at site i and j, namely,
Pij |αi βj i = |βi αj i (αi , βi ∈ {0, 1, 2, 3}). It canbe written by using the flavor-changing operator Sαβ (i) = |αi ihβj | as
H0 = JJ(2)

hijiinter αβ

We specifically investigate the ground state of
the model in the presence of the population
imbalance (effective field), which favors the fla-

H = H0 + Himb
(3)
X
Himb = Dz[−S00 (i) − S11 (i) + S22 (i) + S33 (i)].
i(4)
In the limit of J 0 /J = 0 or J/J 0 = 0, thesystem decouples into four-site chains. The
ground state is an SU(4) singlet for lower fields,
while it is a resonating-valence-bond (RVB)
state consisting of the 0 and 1 flavors forhigher fields. In the presence of weak couplings
among four-site chains (0 < J 0 /J  1), a
previous cluster-mean-field study with a small
2×1 sublattice structure [6] suggests the emergence of the crystal, which is a mixture of thesinglet and RVB states, in the intermediate
field region.

286

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

0.8
energy density

0.9 J 0/J = 0.2
1.0
1.1
1.2
1.3
1.40.0

RVB
singlet-RVB crystal
SU(4) singlet
Ground state
0.2
0.4
0.6
Dz

0.8

1.0

Figure 2: Energy comparison of selected
ground-state candidate states at J 0 /J = 0.2.
The states are obtained for the bond dimensions D = 7. In the intermediate field region(0.5 . Dz /J . 0.6), we have confirmed thatthe ground state is the singlet-RVB crystal.
To investigate the stability of such a crystal state beyond the mean-field study, weapply the 2D tensor-network method basedon the infinite projected entangled pair state(iPEPS) [7,8]. The ground-state candidate singlet and RVB states can be represented by the
iPEPS with the bond dimensions D = 7 and
D = 3, respectively. The singlet-RVB crystal can also be represented by the iPEPS with
D = 7. We consider a variety of crystal statesthat can be realized within a 2 × 2 sublattice.
We take these states as initial states and optimize these wave functions. We adopt the
TeNeS library [9] and calculate the groundstate candidates at D = 7 using the simple
update algorithm [10].
Figure 2 shows the energy comparison of
selected ground-state candidate states in thepresence of weak couplings J 0 /J = 0.2. The
iPEPS energy is slightly lower than the result obtained by the previous cluster-meanfield study [6]. On the other hand, the phaseboundaries are found to be nearly consistent
with the previous study; the singlet-RVB crystal (0.5 . Dz /J . 0.6) emerges between
the lower-field SU(4) singlet and the higherfield RVB states. The ground-state singletRVB crystal is found to be a mixture of theSU(4) singlet and RVB states in equal proportions, and mixing in different ratios result inhigher energy. As for the ground-state candidate states we have studied, the intermediate
singlet-RVB crystal is stable in the presence of
weak couplings among four-site chains.

References[1] C. Wu et al., Phys. Rev. Lett. 91, 186402(2003).
[2] C. Honerkamp and W. Hofstetter, Phys.
Rev. Lett. 92, 170403 (2004).
[3] M. A. Cazalilla et al., New J. Phys. 11,103033 (2009).
[4] A. V. Gorshkov, Nat. Phys. 6, 289 (2010).
[5] P. Corboz et al., Phys. Rev. Lett. 107,215301 (2011).
[6] Y. Miyazaki et al., AIP Advances 11,025202 (2021).
[7] T. Nishino et al., Prog. Theor. Phys. 105,409 (2001).
[8] F. Verstraete and J. I. Cirac, arXiv:condmat/0407066.
[9] Y. Motoyama et al., arXiv:2112.13184.
[10] H. C. Jiang et al., Phys. Rev. Lett. 101,
090603 (2008).

287

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

Simulated scanning tunneling microscopy images ofgraphene nanoribbons with armchair edges
Junhuan LI, Kenta ARIMA
Graduate School of Engineering, Osaka University, Yamada-oka, Suita, Osaka 565-0871
We are aiming at elucidating local electronic
structures of finite-size graphene sheets and theirribbons. Besides, the bandgaps of ribbons were
investigated based on their band structures. It

impacts on scanning tunneling microscopy

helps to check the validity of our simulations and(STM) images [1]. For this purpose, we

discuss more details in simulated images.

simulated STM images of armchair-edged
graphene nanoribbons (AGNRs), using a(a) W=3a-2

simulation package named STATE (Simulation
Tool for Atom TEchnology).

W=22

W=25

W=28

W=26

W=29

W=27

W=30

We simulated STM images at a 𝑉! of −0.05
V for AGNRs with the width ranging from 𝑊 =22 to 𝑊 = 30 (in terms of ribbon width, 𝑊 is(b) W=3a-1

defined as the number of dimer lines across the

W=23

ribbon width). STM images change periodicallyas the ribbon width increases. Figures 1(a)−1(c)summarize the results for AGNRs with 𝑊 =3𝑎 − 2 , 3𝑎 − 1 , 3𝑎 , respectively. For Figure1(a) and 1(b), the electron density is slightly(c) W=3a
W=24

3.7 Å

high at the armchair edges. The dot patterns forma -√3 × √30𝑅30° superstructure in the interiorof the ribbons. In contrast, for the case of 𝑊 =0.0 Å

3𝑎 in Figure 1(c), electrons are not localized at

10 Å

edges and bright dots in the images appear toform a rectangular lattice, which resembles

Fig. 1: Simulated STM images. A unit structure
is superimposed in each image.

superstructures observed experimentally. This
periodic change of simulated images was also
confirmed to continue until reaching a width ofat least 8 nm, in which the value of 𝑊 was 65.

References

Simulated STM images at higher sample
biases, or larger absolute values than that in[1] J. Li, S. Li, T. Higashi, K. Kawai, K.
Inagaki, K. Yamamura, and K. Arima: Phys.
Rev. B 103 (2021) 245433.

Figure 1， were also calculated in the same288

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

Towards a unified understanding of thermodynamic
properties in S=1/2 spherical kagome systems{W72V30} and {Mo72V30}
Yoshiyuki Fukumoto
Department of Physics, Faculty of Science and Technology, Tokyo University of Science
2641 Yamazaki, Noda, Chiba 278-8510

The spherical kagome systems, or quantumspin icosidodecahedrons, are frustrated magnetic clusters with 30 spins. It has beenargued that {W72 V30 } can be described asa regular, non-distorted, icosidodecahedron,and {Mo72 V30 } contains structural distortion[1]. The experimental magnetization curves for
those increase linearly with magnetic field [2],although a cluster Heisenberg antiferromagnet
should show stepwise behavior. In a theoretical study on {W72 V30 }, Schnack et al. has reported that the discrepancies are attributed tothe distribution of nearest-neighbor exchange
couplings, which is called bond-randomness,where the width of the variance in exchangeinteractions was estimated to be 30% of the average value J = 115K [2]. On the other hand,
our group has reported that the incorporationof about 10% of Dzyaloshinskii-Moriya (DM)interactions also leads to the collapse of thestaircase structure [3].
Subsequently, Kihara et al. experimentallymeasured specific heats of both materials at
temperatures T < 11K under magnetic fields
up to 15T [4], and revealed that magnitudesof the specific heats are comparable to calculated results of the non-distorted and distorted
Heisenberg models [1]. For {Mo72 V30 }, temperature dependence of the specific heats has
a shoulder around 2 ∼ 3K and depends onmagnetic field, which is consistent with calculated results with the distorted Heisenberg
model [1]. However, for {W72 V30 }, the calculated specific heats with the non-distorted

289

Heisenberg model shows a peak around 2K andsignificant magnetic field dependence, but theexperimental result has no peak and does notshow magnetic field dependence. As for thedisappearance of the 2K peak, a possible origin
is blurring of the density of states due to thebond randomness. Because the 2K peak comesfrom low-energy singlet excitations, we expectthat the 2K peak is readily dispersed by about
10% bond randomness, which is much smaller
than Schnack’s value [2]. We actually calculate
the density of state and specific heat to confirm
this expectation [5]. However, quantitative reproduction of experiment is not achieved.
It is interesting to study whether larger bondrandomness leads to quantitative reproductionof experimental specific heat or not. We also
calculate the case of 30% bond randomness
and find that the magnetic field dependence
becomes smaller and the shape of the curve
approaches to the experiment, which suggests
an adequate choice of bond randomness, which
would be 40-50%, gives quantitative reproduction. On the other hand, because the specificheat of {Mo72 V30 } has magnetic field dependence and a shoulder in it’s temperature dependence, we naively expect {Mo72 V30 } has
no bond randomness. The investigation of theorigin of presence or absence of the bond randomness will be an issue in the future.
[1] N. Kunisada et al., PTEP 2014, 41I01 (2014).
[2] J. Schnack et al., arXiv:1304.2603v1.
[3] Y. Fukumoto et al., JPSJ 87, 124710, 2018.
[4] T. Kihara et al., PRB 99, 064430 (2019).
[5] M. Motohashi et al., arXiv:2112.07923v1.

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

Numerical studies on excitation spectra of
Heisenberg antiferromagnets on the triangular andkagome lattices
Yoshiyuki Fukumoto
Department of Physics, Faculty of Science and Technology, Tokyo University of Science
2641 Yamazaki, Noda, Chiba 278-8510

This project consists of two research topics, (i) virational study of the triangular-latticesystem with the fermionic representation and(ii) series expansion study of the kagomelattice system in Cs2 Cu3 SnF12 .
As for (i), a neutron scattering experiment on Ba3 CoSb2 O9 revealed the existenceof spinon continuum in the dynamical structure factor, in addition to magnon excitations[1]. Zhang et al. used the fermionic spinon operators to study the spinon continuum within
random-phase approximation (RPA) [2]. This
RPA treatment, however, does not exclude unphysical states exactly. This year, we havedeveloped a program code using Monte-Carlo
method to obtain the energy of the variational
wave function proposed by Zhang et al. [1],completely considering the constraints. We
also developed a program code to calculate thevariational energy for a 36-site system exactly,
without using Monte-Carlo sampling. The
variational parameter is the strength of thestaggered magnetic field conjugate to the 120
degree structure. We use wavenumbers thatsatisfy periodic or anti-periodic conditions in
two directions of the triangular lattice. There
are four combinations of wave numbers, andwe make four states from each. The optimal
combination of the four states is found by solving the generalized eigenvalue problem. In thenext fiscal year, we will start developing a program code for the particle-hole excitations to
obtain the dynamic structure factor.
As for (ii), Cs2 Cu3 SnF12 has the q=0 order,
and a theoretical analysis, based on the lin-290

ear spin wave theory, of the observed magnon
dispersion resulted in the exchange parameter
which is 40% less than the result by an analysis
of susceptibility. This fact indicates diﬃcultyin the linear spin wave theory for frustrated
systems.
In this study, we calculate the dispersion relation by series expansion from the Ising limit.
We include Dzyaloshinskii-Moriya (DM) interactions into our model, where we set the DM
vector to be perpendicular to the plane on
which ordered moment lies. The direction of
the DM vector is chosen to stabilize the q=0
state. The linked cluster expansion algorithm
is used to make expansion coeﬃcients. We developed a calculation program using the previous research that applied the cluster expansion method to the triangular-lattice system
as a benchmark test [3], and extended it to the
kagome-lattice system. As a result, it is clarified that the bandwidth of the magnon spectrum in the series expansion is about 1/4 ofthe linear spin wave theory, and its shape isalso changed significantly. We also find thatthe correction for the linear spin wave theory
is larger in the kagome lattice than in the caseof the triangular lattice. Also, the shape ofthe magnon spectrum obtained by series expansion is comparable to the experimental results of Cs2 Cu3 SnF12 . It is speculated that theagreement can be improved by making the DM
interaction more realistic in the future.
[1] S. Ito et al., Nat. commn. 8, 235, 2017.
[2] C. Zhang et al., PRB 102, 075108, 2020.
[3] W. Zheng et al., PRB 74, 224420, 2006.

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

Gapless symmetry-protected topological phase ofquantum spin system
Yuichiro HIDAKA1 , Shunsuke C. FURUYA2,3 , Atsushi UEDA1 , and Yasuhiro TADA4
1
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581
2
Department of Physics, Ibaraki University
Mito, Ibaraki 310-8512
3Department of Basic Science, University of Tokyo
Meguro, Tokyo 153-8902
4
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 symmetry-protected topological (SPT) phases. A SPT phase is a gapped
phase which preserves some given symmetries
but not adiabatically connected to a trivial
phase without breaking the symmetries [1].
Recently many researchers attempt to extend
the concept of SPT phase into gapless system[2, 3].
In the present study, we proposed a quantum
spin model that has the properties of gaplessness and nontrivial SPT phase simultaneously.
The model we proposed is the anisotropic triangular strip (ATS) XXZ model, shown in
Fig.1, which consists of three-leg antiferromagnetic ladder with frustrated interchain couplings. We used the density matrix renormalization group (DMRG) method to calculate
the ground state of the ATS XXZ model at
the specific parameters, J = 1.0, ∆ = 0.8, and
J× = 0.5. We used the ITensor and TeNPy
libraries [4, 5].

that the central charge of the ATS XXZ model
is one. The calculation result of central charge
is shown in Fig.2, and the central charge become close to one, which is derived from the extrapolation of the entanglement entropy about
the system size. In addition, Fig.3 shows thecalculation result of entanglement spectrum
and we can see double degeneracy, which is an
important feature of Haldane state in quantum spin system. From these two results, theground state of ATS XXZ model is TomonagaLuttinger liquid which has the property of nontrivial SPT phase. Therefore we can conclude
that the ground state of ATS XXZ model with
J = 1.0, ∆ = 0.8, and J× = 0.5, is in gapless SPT phase. We presented these results atthe 2022 Annual (77th) Meeting of the Physical Society of Japan [6].

We found that the entanglement entropy of
the ATS XXZ model fits with Calabrese-Cardy
formula of open boundary condition. According to the conformal field theory, it is expected

Figure 1: ATS XXZ model

291

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

Figure 2: Entanglement entropy of the ATS
XXZ model. The purple line is the fitting
line of entanglement entropy by the CalabreseCardy formula.

Figure 3: Entanglement spectrum of ATS XXZ
model.

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] M. Fishman, S. R. White, and E. M.
Stoudenmire, arXiv:2007.14 (2020).
[5] J. Hauschild and F. Pollmann, SciPost
Phys. Lect. Notes , 5 (2018).
[6] Y. Hidaka, S. C. Furuya, A. Ueda, and Y.
Tada: in preparation.

292

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

Development of multiscale simulation techniquefor liquids and solids
Takahiro MURASHIMA
Department of Physics, Tohoku University,
6-3, Aramaki-Aza-Aoba, Aoba-Ward, Sendai, Miyagi 980-8578
Predicting the flow and deformation of soft

elongational flow.)matter and composite materials is still a very
difficult problem. Aside from the problem ofsolving macroscale continuum equations, there
are difficulties in determining constitutive lawsthat characterize materials. An approach todetermine the constitutive law (rheology) of amaterial using molecular simulations has been
popular in recent years. Multiscale simulations
that connect these molecular simulations tocontinuum simulations as constitutive laws do

Fig.1: Polymers under uniaxial elongation 𝜅𝜅 =𝜖𝜖̇00
�0 𝜖𝜖̇/2 0 �.
0 0 𝜖𝜖̇/2

not explicitly require constitutive laws, thus
avoiding the difficulty of deriving constitutive
laws [1]. However, there is a problem ofboundary conditions for treating general flow

For general deformations in the molecularfields in molecular dynamics simulations.

dynamics simulation, we proposed a method

Therefore, in this research project, wethat uses QR decomposition to exclude theattempted to extend the boundary conditions torotational component of the deformation tensor

solve this problem.

and impose only pure deformations on the

Recently, we have developed the UEFEX

molecular simulation system [3]. Here, we

method [2], which enables large-elongationconsidered problems in which large

flow simulations. Here, the elongation axis anddeformations and rotations do not occur. This

the unit cell are not parallel to each other totime, we extended this method to include

avoid collapsing the unit cell under

vortex flow in fluids. (Figure 2 shows aelongational deformation. （Figure 1 shows thepolymer chain in a rotating flow.) Before

steady state of polymer chains under uniaxial

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

conducting multiscale simulations using thismethod, we are currently conducting an

References

analysis assuming a corner vortex phenomenon[1] T. Murashima, T. Taniguchi, EPL 96 (2011)that appears in polymer melts. We intend to18002.

summarize the results in a future paper.

[2] T. Murashima, K. Hagita, T. Kawakatsu,
Nihon Reoroji Gakkaishi 46 (2018) 207.
[3] T. Murashima, S. Urata, S. Li, Euro. Phys. J.
B 92 (2019) 211.

Fig. 2: Polymers under rotational flow 𝜅𝜅 =0 𝑎𝑎̇ 0�−𝑏𝑏̇ 0 0�, where 𝑎𝑎̇ > 𝑏𝑏̇.
0 0 0

294

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

Ashkin-Teller-like phase transition in thetwo-dimensional classical dimer-monomer model
Satoshi MORITA
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

An attractive interaction between dimersaligning on a square induces finite-temperaturephase transitions [1, 2]. The phase transition between the liquid and columnar phases
is known as the Kosterlitz-Thouless transition.
When monomers are introduced, the criticalexponents smoothly change and its universality class is the same as the Ashkin-Teller(AT) model. On the other hand, another ordered phase appears between the two phases
when the nearest-neighbor attractive interaction along the dimer direction is sufficiently
large [3]. This phase called the nematic phasebreaks only the lattice rotational symmetry.
Although preceding studies reported these interesting features, an intermediate parameterregion, where two attractive interactions arecomparable, were not investigated. It is notknown yet how the phase transition between
the columnar and disordered phases splits into
two transitions involving the nematic phase.
To address this problem, we investigate aclassical dimer-monomer model on the squarelattice by using the Monte Carlo (MC) method
and the tensor network (TN) method. For the
MC method, we generalize the worm algorithm
to the grand-canonical ensemble of this dimer
model. The TN method for infinite systems iscomplementary to the MC method. From the
TN representation of the partition function,
the row-to-row transfer matrix is constructed
as the matrix product operator (MPO). Itseigenvector with the maximum eigenvalue isapproximated by the uniform matrix prod-

uct state (uMPS). We determine the phasestructure and find the existence of the multicritical point. The universality classes of thephase transitions are discussed by estimating
the scaling dimension and the central charge.
We find that the columnar and nematic orders
correspond to the magnetization and polarization in the Ashkin-Teller model, respectively.
The estimated scaling dimensions clearly indicate the Ashkin-Teller universality class. We
also confirm the 4-state Potts universality at
the multi-critical point.

Figure 1: A typical configuration of the nematic and disordered phases generated by theMC simulations.

References[1] F. Alet, et al., Phys. Rev. Lett. 94, 235702(2005).
[2] F. Alet, et al., Phys. Rev. E 74, 041124 (2006).
[3] S. Papanikolaou, D. Charrier, and E. Fradkin,
Phys. Rev. B 89, 035128 (2014).
[4] S. Morita, H. Lee, K. Damle, N. Kawashima,in preparation

295

Activity Report 2021 / 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

















    99#    ' '   '  %  #  ' '                           &                                   &        (            1       '  '    ' % ' :     ' #  '       +   +    +  )
               !+"    !4;"    ! "           /     0  !//"                ' '      +                         (       !  "   !  "           
 8          
            
 !'      "     , +       !   "              


 












 



 





























296

 

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

1.0

0.6

SD

JAF

ΓF

XY1
nTLL

SD
0.4

0.5

uudd

uudd
F
0.2

nTLL
F

0.0
−1.0

−0.5

0.0

0.0
0.0

0.5 Δ 1.0
AF

XY1
ΓF

0.5

1.0

+  < =           ' +  )< =           ' #  '   ' 








JAF

0.3
XY1

SD

0.2
nTLL
0.1

0.0

F

uudd

       4;     ; &      &.                  ,   !? ="  @     
 3  ,          ? =0    # ! "! 9 @" # ! 9 @" 
# ! "  ! '"   ! '"    //                       ,      3  ,      ,      '           '          3  ,          # ! "# ! " 
# ! "  !  "   ! '" /      3  ,              +      //           ! "  ! '"&          //                     -   3   .          7  $.  7   >  
A 4.      8 .  7 $.  4 6
+     4.  B ?  4 B  '@''
!':" C / 5,  .  4 ?  BD?  E'!FG"1 B H 7  ,  ; B    B?  6  : !FG)") 7 >   2 7 , B ?  2 G)@ !FF:"@ H ? >    ?  2  *E !FGE"


−0.4

−0.2

0.0 ΔAF 0.2+  < =           ' % ' :                ! "                   (           
3  ,        :      -   ,          3   
  ! "   +      ! '"      &          3  ,                                           
-    3  ,      +           4;      //      5,  . 7  ,     &        
     >  
7 , )
      
 3  ,             ! '"  !  " 4 



297



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

First Principle Calculation and Charge Order Structure in
Organic Charge Order Materials
Masatoshi SAKAI, Ryoma ISHII
Department of Electrical and Electronic Engineering,
Chiba University, 1-33 Yayoi-cho Inage, Chiba 263-8522

We tried to calculate the electronic structure

that, Wannier calculation was carried out byof -(BEDT-TTF)2PF6 from the first principle

RESPACK to calculate Wannier orbitals so that

calculation

Espresso,

the band dispersion curve around HOMORESPACK, and mVMC. Starting from theorbitals in scf calculation is reproduced byloading cif file into Quantum Espresso, we

using the Wannier function. Fig. 1 (a) shows

obtained the results of the scf calculation. Afterone of the selected Wannier functions used inusing

Quantum

this calculation, and Fig. 1(b) shows the band(a)dispersion curves of the Wannier function
superimposed into the result of the scf
calculation. Wannier function represents themolecular orbital of the BEDT-TTF, and band
dispersion curves are also well reproduced.
Therefore, we obtained reasonable results in(b)

RESPACK calculation in the 3x3x3 K-meshcondition. Then we carried out the cRPA
calculation

toobtaindielectric

function,

Coulomb integrals, and exchange integrals. We(c)spent a long time converging these calculations;however,weobtainedtheinteraction

parameters in -(BEDT-TTF)2PF6. Now we aretrying to execute mVMC calculation both in Fig. 1 (a) One of the calculated Wannierorbital around the Fermi level of (BEDT-TTF)2PF6. (b) Band dispersion curvesof Wannier function superimposed in the scf
calculation. (c) Coulomb interaction as afunction of inter center distance.

(BEDT-TTF)2PF6

and

-(BEDT-TTF)2I3

organic charge order system to elucidate ourexperiment.

298

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

Structural analysis of surface superstructures and atomic layermaterials using by 2DMAT on supercomputer system at ISSP
Akari TAKAYAMA
Department of Physics,
Waseda University, Ohkubo, Shinjuku-ku, Tokyo 169-8555
In our research, we use supercomputers forsupercomputer is very useful for structure

structural analysis in terms of determining theanalysis of large periodic structure and surface

structure of materials. Since structure andsuperstructure

physical properties are closely related, it isbecause they require many variables for analysis.

essential to know the exact structure in order toFigure 1 shows the rocking curve under thewithlatticereconstruction

correctly understand the physical properties.

one-beam

Generally, surface and interface structures areexperiment (circles) and from calculation (solid

complicated than bulk crystal due to structuralline) for 0-BL graphene (buffer layer) of

relaxation and structural reconstruction at

SiC(0001). From the fitting calculation for thissurface and interface. For instance, structure ofrocking curve, we determined the buckling

epitaxial graphene growth on SiC is still

structure of buffer layer and found that ourcontroversial due to the underlying-buffer layer

determined structure corresponds well with thewith a long-period (6√3×6√3-R30°). Also thestructural model proposed by Lima et al. [2].

√3×√3-Pb/Si(111) surface superstructure, which

Moreover, the ratio of bounding and free-

is known to undergo a structural phase transition

standing regions with/from substrate wereat T = 86 K, forms a “mosaic” structure, thus its

determined to be 22 % and 78 %, respectively,

detailed structure is fully understood. In thisin this sample. The ratio is consistent withstudy, we investigated the structure of graphene

previous studies [3, 4].

condition

obtained

from

the

on SiC including a buffer layer and Pb/Si(111)
surface superstructure by using total reflection
high-energy positron diffraction (TRHEPD)
method. For structure analysis, we used
“2DMAT” [1], an open source framework for 2D
material structure analysis, which is included as

Fig. 1: Structure analysis results of buffer

standard software in ISSP’s supercomputer.

layer on SiC.

Using massive parallel computer simulations on

299

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

References

Figure 2 shows the rocking curve under the
one-beam

condition

obtainedfrom[1] Y. Motoyama, et al., arXiv:2204.04484

theexperiment (circles) and from calculation (solid[2] L. H. de Lima, et al, Chem. Commun., 50line) for mosaic phase of √3×√3-Pb/Si(111). We(2014) 13571.

have succeeded in determining the structure of[3] M. Conrad, et al., Phys. Rev. B 96 (2017)

the mosaic phase by analyzing the Pb occupancy

195304.

as a variable. In addition, no noticeable changes[4] J. D. Emery, et al., Phys. Rev. Lett. 111

in the rocking curves were observed at low (T =(2013) 215501.

13 K) and room temperatures. This suggests that[5] I. Brihuega et al., Phys. Rev. Lett. 94 (2005)

the phase transition in this structure is caused by

046101.

soft phonons [5].

Fig. 2: Structure analysis results of mosaic
phase of Pb/Si(111).

300

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

Minority Driving Transition in Collective Movement
Katsuyoshi Matsushita
Department of Biological Science, Osaka University
Machikaneyama, Toyonaka 560-0043

In biological organisms, cells aggregate andform tissues to function. Cell aggregation
exhibits various dynamics, including rotationand translational movements. The transition
from the rotation to the translational movements is well known in the developmental process in Dictyostelium Discoideum (dicty) [1].
The physical mechanism of this transition has
been theoretically investigated on the basis of
the behavior of cells [2, 3]. In particular, theworks consider the response of cells to chemical gradients, which is observed in a earlier
stage of development. These works successfully explained the transition and apperently
unified the cell movements throughout the developmental process of dicty combined withthe works for later stages [4].
These works assume homogeneous cell systems. However, as well known, the transition istightly associated with cell differentiation. Cell
differentiation means the change in the material properties of the cells. Therefore, the response property to the chemical gradient may
not be conserved in the transition. Recent
advances in the investigations of this transition showed the lack of the response of dicty
cells [5]. Furthermore, both the rotational andtranslational movements do not depend on the
response [6, 7]. These results imply the existence of another physical mechanism associated with cell differentiation.
A possible mechanism is the material parameter change in the cells. In particular, therate of differentiation in population graduallychanges during the transition process with several hours [1]. Therefore, we can expect thefollowing situation: The collective movement
is rotational when the differentiated cells are aminority in the cell population. Otherwise, thecollective movement is a transnational movement. Namely, the minority change in cell population in the differentiation process drives thetransition in collective movements.
To examine this scenario of the transition,
we theoretically considered the process of cell
population change [8]. We use the cellular
Potts model [9, 10]. The cell population is assumed to consist of leader and follower cells[11]. The leader cells correspond to the differentiated cells and can lead the following cells
using the cell-cell adhesion [12]. We also consider cell polarities by memorizing the cell trajectories [13, 14]. These polarities correlatewith the adhesion [15] and stabilize the ordering of the cell movements [16, 17].
From the simulation results based on theabove model, we examine the population dependence of the collective cell movements. We
showed that the minority change in cell population in the differentiation process drives thetransition of collective movements. A motilityinducing phase separation (MIPS) inner thecell aggregation appears for intermediate populations. The MIPS divides the cell aggregation into two aggregations: a homogeneous
un-differentiated population and a highly differentiated population. Therefore, this minority change cannot explain the transition fromrotation to translational movements for stable
cell aggregations. The improvement of the scenario based on the minority change remains in
the future.

301

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

The clarification of the origin of this aggregation division is expected to give us hints toimprove the scenario. In our recent work, weinvestigate a homogeneous cell system to identify the parameters of collective movements[18]. We find that the occurrence of the collective movement, including the rotation andtransnational movements, depends on the interface tension of the cell aggregations. This
result implies the transition due to minoritychange originates from a change in the interface tension of the aggregation. The minority change is speculated to reduce the interface
tension, and aggregation division is associatedwith the transition. This speculation implies
the fine-tuning of the interface tension in the
differentiation. However, it is not realizable ina realistic system and the necessity for anotheridea.

References[1] J. T. Bonner, The Social Amoebae: The
Biology of Cellular Slime Molds (Princeton University Press, Princeton, 2009).
[2] T. Umeda and K. Inouye, J. theor. Biol.
219, 301 (2002).
[3] B. Vasiev C. J. Weijer. J. Theor. Biol.
223, 347 (2003)

Fujimoto, Proc. Sympo. Traffic Flow
Self-Driven Particles 26 39-42 (2020).
[9] F. Graner and J. A. Glazier, Phys. Rev.
Lett. 69, 2013 (1992).
[10] F. Graner, J. Theor. Biol. 164, 455(1993).
[11] A. J. Kabla, J. R. Soc. Interface 9, 3268
(2012).
[12] K. Matsushita, Phys. Rev. E 97, 042413(2018).
[13] P. Szabó, M. Nagy, and T. Vicsek, Phys.
Rev. E 79,021908 (2007).
[14] K. Matsushita, Phys. Rev. E 97, 042413(2018).
[15] J. C. Coates and A. J. Wood. J. Cell Sci.
114, 4349 (2001)[16] K. Matsushita, K. Horibe, N. Kamamoto,and K. Fujimoto, J. Phys. Soc. Jpn. 88,
103801 (2019).
[17] K. Matsushita, S. Yabunaka, and K. Fujimoto, J. Phys. Soc. Jpn. 90, 054801(2021).
[18] K. Matsushita, N. Kamamoto, M. Sudo,
and K. Fujimoto, Proc. Sympo. Traffic
Flow and Self-Driven Part. 27, 23 (2021).

[4] A. Maree et al. J Their. Biol. 199, 297(1997); Proc. Natl. Acad. Sci. USA 98,
3879 (2001) (2000).
[5] H. Hashimura, Y. V. Morimoto, M. Yasui,
and M. Ueda, Communications Biology 2,
34 (2019).
[6] W.-J. Rappel, A. Nicol, A. Sarkissian, H.
Levine, and W. F. Loomis, Phys. Rev.
Lett. 83, 1247 (1999)[7] Y. Kida, K. Pan, and H. Kuwayama, Differentiation 105, 71 (2019).
[8] K. Matsushita, S. Yabunaka, H.
Hashimura, H. Kuwayama, and K.

302

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

Non-uniform thermal transport properties in proteins
Takahisa YAMATO
Graduate School of Science,
Nagoya University, Furo-cho, Chikusa-ku, Nagoya 464-8602

Protein molecules are thermally fluctuating andautocorrelation function is shown in blue. Its

tightly packed amino acid residues strongly

short time behavior is shown in the inset. Theinteract with each other. Such interactions arevalue of λ is shown in orange as a function ofcharacterized in terms of heat current at the

integration time, τ, for τ > 3.0 ps. The red line

atomic level. We calculated the thermal

indicates the running mean of λ with a window

conductivity of a small globular protein, villin

size of 100 fs.

headpiece subdomain, based on the linear
response theory using equilibrium molecular
dynamics simulation. The value of its thermal
conductivity was 0.3 ± 0.01 [W m –1 K–1],
which is in good agreement with experimental
and computational studies on the other proteinsin the literature (Figure 1). Heat current along
the main chain was dominated by local
vibrations in the polypeptide bonds, with amide

Figure 2: Fourier transform of the backbone

I, II, III, and A bands on the Fourier transformheat current autocorrelation function. Theof the heat current autocorrelation function [1]Fourier transform of the average heat current(Figure 2).

autocorrelation function is shown for the
backbone heat current.

References[1] T. Yamato, T. Wang, W. Sugiura, O.
Laprévote, T. Katagiri: J. Phys. Chem. B (2022)
https://doi.org/10.1021/acs.jpcb.2c00958

Figure 1: Thermal conductivity. The time
evolution of the average heat current

303

Activity Report 2021 / 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]. We have also demonstrated
that a Yang-Lee anyon system is constructedfrom Majorana bound states in topological superconductors [2].
In Ref. [1], we numerically simulated braiding dynamics of MKPs in a network of onedimensional TRITSCs, and examined the tolerance against various perturbations which
may cause decoherence of MKPs. We first considered eﬀects of a magnetic field which breaks
time-reversal symmetry. In contrast to a naive
expectation, we find the non-Abelian braidingof MKPs is robust against applied magnetic
fields provided that the initial and final statesof a braiding process are invariant under the
combination of a time reversal and a mirror reflection, even when intermediate states breakthe combined symmetry. In addition, we investigated the tolerance of non-Abelian braidings against gate-induced inhomogeneous potentials at junctions between superconducting
nanowires. The potential generally generates
a non-Majorana nearly zero-energy Andreeev
bound state at the junctions. We demonstrated that the non-Majorana states interferewith MKPs, resulting in the failure of nonAbelian braidings.
Yang-Lee anyons, described by the nonunitary conformal field theory with the central
charge c = −22/5, are non-unitary counterparts of Fibonacci anyons, obeying the sameron bath)

Metal (electFigure 1: Schematics of Yang-Lee anyon system with Majorana quasiparticles.
fusion rule. Therefore, the Yang-Lee anyon
system is expected to be a platform for realizing universal topological quantum computation. In Ref. [2], we considered a topologicalsuperconductor junction system coupled with
dissipative electron baths (see Fig. 1), which
realizes a non-Hermitian interacting Majorana
system. Numerically estimating the central
charge, we examined the condition that the
non-Hermitian Majorana system can simulate
the Ising spin model of the Yang-Lee edge singularity, and confirmed that, by controllingmodel parameters in a feasible way, the YangLee edge criticality is realized. We also discuss
the scheme for the measurement and the braiding of Yang-Lee anyons, aiming at the application to the universal quantum computation.

References[1] Y. Tanaka, T. Sanno, T. Mizushima, andS. Fujimoto: submitted.
[2] T. Sanno, M. G. Yamada, T. Mizushima,and S. Fujimoto: submitted.

304

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

Extension of susceptibilities, screened exchange andspin-fluctuation integrals into ultrasoft
pseudopotentials
Mitsuaki KAWAMURA
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

In the Eliashberg theory or density
functional
theoryforsuperconductors(SCDFT) [1], where the electron-phonon,
screened Coulomb, spin-fluctuation mediatedinteractions are included non-empilically, weneed to compute the folowing product of two
Kohn-Sham orbitals
ρnkn0 k0 (r) = ϕ∗nk (r)ϕn0 k0 (r),

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0

15

20

25

30

(b) Nb : Spin fluctuation

0.2
Averaged interaction

Bloch wavenumber, respectively. To perform
this product together with the ultrasoft pseudopotentials (USPP) or projector augumented
waves (PAW) that are widely used because of
the good accuracy and reasonable numericalcosts [2], we need a correction term for the
norm-conservation. Such a correction is originally proposed for the calculation of the susceptibility [3] asX∆ρnkn0 k0 (r) =hϕnk |βτ i0 ihβτ i |ϕn0 k0 iQτ ii0 (r),

NC :
US :
US wo Q :0.10q-point index
0.8(c) V : Screened exchange

0.6
0.4
0.2
0
0
1.4
1.2
1
0.8
0.6
0.4
0.2
0

τ ii0

(2)

To utilise this correction to the calculationof the spin-fluctuation mediated interaction,

10

0.3

where n (n0 ) and k (k0 ) are the band index and

P S∗
PS
AE
Qτ ii0 (r) ≡ ψτAE∗i (r)ψτ i0 (r) − ψτ i (r)ψτ i0 (r).
(3)

5

0.4

(1)

where βτ i is the projector dual to the atomic
pseudo orbital ψτPiS at atom τ and orbital i,and the augumentation charge Qτ ii0 (r) is computed from the pseudo (PS) and all-electron(AE) atomic orbitals as follows:(a) Nb : Screened exchange

5

10

15

20

25

30

35

(d) V : Spin fluctuation

Figure 1: Spin-fluctuation mediated- (a, c) andscreened exchange- (b, d) interactions averaged over the Fermi surface of Nb (a, b) and V(c, d). “+”, “×”, and “◦” idicate results with
norm-conserving (NC) ultrasoft (US) with andwithout augumentation charge Q(r), respectivelly.

305

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

References[1] M. Lüders, M. A. L. Marques, N. N. Lathiotakis, A. Floris, G. Profeta, L. Fast, A.
Continenza, S. Massidda, and E. K. U.
Gross Phys. Rev. B 72, 024545 (2005).
[2] G. Prandini, A. Marrazzo, I. E. Castelli,
N. Mounet and N. Marzari:npjComputational Materials 4, 72 (2018).
https://materialscloud.org/sssp[3] M. Gajdoš, K. Hummer, G. Kresse, J.
Furthmüller and F. Bechstedt: Phys. Rev.
B 73, 045112 (2006).

100
Computed Tc [K]we implemented this formalism into our firstprinciples program package SuperconductingToolkit [4] which is based on SCDFT.
Figure 1 shows the screened exchange andspin-fluctuation interactions averaged over the
Fermi surface of Nb and V. In Nb, the effect of the augmentation charges Q is small
because the USPP for this atom is almost
norm-conserving. At the same time, V has
a significant contribution from Q due to thenodeless 3d orbitals. In Fig. 1(d), we cansee a deviation between the result by NC and
US pseudopotential even if we include thecorrection. This deviation may be because
the exchange-correlation kernel included in the
spin-fluctuation is sensitive to the charge density in the vicinity of atoms.
We also performed the benchmark of the calculation of Tc for 15 materials, namely Al, V,
Ta, In, Zn, Cd, Sn, ZrN, TaC, MgB2 , H3 S
(at a pressure of 200 GPa), CaC6 , YNi2 B2 C,
and V3 Si. Figure 2 shows the experimental,
and calculated Tc ; we performed four kindof calculations by changing superconducting
density functional, namely the conventional
plasmon-assisted [5], Sanna’s Eliashbergmimic (Sanna) [6], Sanna+Coulomb renormalization (ZC ) [7], and Sanna + ZC + Spinfluctuation [9, 8] functional.

10
Conventional :Sanna :
Sanna/ZC :Sanna/ZC/SF :

1
1

10

100

Experimental Tc [K]Figure 2:
Computed- and experimental
Tc . “+”, “×”, “4”, and “5” indicate result with conventional plasmon-assisted [5],Sanna’s Eliashberg-mimic (Sanna) [6],Sanna+Coulomb renormalization (ZC ) [7],and Sanna + ZC + Spin-fluctuation [9, 8]
functional, respectively.
[4] M. Kawamura, Y. Hizume, and T.
Ozaki: Phys. Rev. B 101, 134511 (2020).
https://sctk.osdn.jp/index.html[5] R. Akashi and R. Arita Phys. Rev. Lett.
111, 057006 (2013).
[6] A. Sanna, C. Pellegrini, and E. K. U.
Gross Phys. Rev. Lett. 125, 057001(2020).
[7] A. Davydov, A. Sanna, C. Pellegrini, J. K.
Dewhurst, S. Sharma, and E. K. U. Gross
Phys. Rev. B 102, 214508 (2020).
[8] K. Tsutsumi, Y. Hizume, M. Kawamura,
R. Akashi, and S. Tsuneyuki Phys. Rev.
B 102, 214515 (2020).
[9] F. Essenberger, A. Sanna, A. Linscheid,
F. Tandetzky, G. Profeta, P. Cudazzo,
and E. K. U. Gross Phys. Rev. B 90,
214504 (2014).

306

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

Development of Algorithms for Ising Machines Basedon Statistical Mechanics
Shu TANAKA
Department of Applied Physics and Physico-Informatics, Keio University
Hiyoshi, Yokohama, Kanagawa, 223-8522

Ising machines are a promising computational technique for highly eﬃcient combinatorial optimization problems. Simulated annealing and quantum annealing are known as
the internal algorithms of Ising machines, anda thorough investigation of the statistical dynamic properties of the Ising model is expected
to improve the performance of Ising machines.
In addition, investigating new ways of utilizing Ising machines is also an essential topic in
the ﬁeld of Ising machines. In this research
project, we concentrated on the following two
points.
(I) Dynamical properties of onedimensional random Ising model withrandom transverse ﬁeld [1]Quantum annealing is a method of introducing an Ising Hamiltonian and Hamiltonian
expressing noncommutative quantum ﬂuctuations. We wish to search the ground state, andby varying the coeﬃcients of each Hamiltonian with time, we eventually ﬁnd the groundstate of the Ising Hamiltonian. One of the factors that deteriorate the computational performance of quantum annealing is the reduction of the energy gap. The study of quantum phase transitions, where the energy gap
systematically decreases with system size, is acommon strategy for understanding the performance of quantum annealing from the perspective of statistical mechanics.
We considered an Ising model with random interactions deﬁned as a one-dimensional
chain. We investigate the system size depen-

dence of the energy gap and the behavior ofthe dynamical critical exponent when a transverse magnetic ﬁeld with randomness uncorrelated with the randomness of the interaction isapplied and when a transverse magnetic ﬁeld
with randomness correlated with the randomness of the interaction is applied. As a result,it was found that there is a diﬀerence between
the two cases.
(II) Black-box continuous optimization
using Ising machines [2]
An Ising machine is a hardware that takes
an Ising model or equivalent QUBO representation as input format and operates to search
its ground state. Ising machines have been applied in various situations to combinatorial optimization problems that Ising models or QUBOs can represent. On the other hand, amethod for utilizing Ising machines for blackbox discrete optimization, i.e., discrete optimization problems for which the objectivefunction is not explicitly given, was proposedby us in 2020 [3]. We extended it and proposeda method to perform black-box continuous optimization with Ising machines.

References[1] T. Shirai and S. Tanaka: Ann. Phys. 435(2021) 168483[2] S. Izawa, K. Kitai, S. Tanaka, R. Tamura,and K. Tsuda: Phys. Rev. Res. 4 (2022)
023062.

307

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[3] K. Kitai, J. Guo, S. Ju, S. Tanaka, K.
Tsuda, J. Shiomi, and R. Tamura: Phys.
Rev. Res. 2 (2020) 013319.

308

Activity Report 2021 / 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

We investigate the ground-state phases of
mixed diamond chains with bond alternation
described by the following Hamiltonian : [1]
H=

L [
∑

(1)

(1 + δ)S l (τ l
(1)

(2)

(2)

(1) (2)+ τ l )S l+1 + λτ l τ l

(1)

]
, (1)

(2)

where S l , τ l and τ l are spin operators with(1)
(2)
magnitudes Sl = τl = 1/2 and τl = 1. The
number of unit cells is denoted by L, and thetotal number of sites is 3L. Here, the parameters λ and δ control the frustration and bondalternation, respectively, as depicted in Fig. 1.
(1)

τl
1+δ
Sl
1+δ

δ=0
δ=0.004
δ=0.01
δ=0.02
δ=0.03
δ=0.04
δ=0.05+ τl )

l=1+ (1 − δ)(τ lmsp (a)1~0.807λc0(0)−
0
0.65

0.7

(b)

1−δ
λ

0.5

1−δ
(2)

τl

Figure 1: Structure of the diamond chain investigated in this work.
Deﬁning the composite spin operators T l as
(1)
(2)
T l ≡ τ l + τ l , it is evident that ∀l [T 2l , H] =0. Thus, we have L good quantum numbers
T 2l ≡ (Tl + 1)Tl where Tl = 1/2 and 3/2. The
total Hilbert space of the Hamiltonian (1) consists of separated subspaces, each of which is
speciﬁed by a deﬁnite set {Tl }.
For large λ, ∀l Tl = 1/2. Hence, for δ ̸= 0,
the ground state is equivalent to that of thespin 1/2 chain with bond alternation δ. Therefore, the ground state is a gapped spin liquid.

309

λ

0.8

δ=0
δ=0.004
δ=0.006
δ=0.008δ=0.01

msp
0.15

Sl+1

0.75

~0.807λc0(0)−

0.1

0.05

0

0.79

0.8

λ

0.81

Figure 2: (a) λ-dependence of msp calculatedby the inﬁnite-size DMRG method. (b) Magniﬁed ﬁgure around λ ∼ λc0 (0).

Activity Report 2021 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo0.81
λc0(δ)

This energy gap vanishes for δ = 0. For small
λ, ∀l Tl = 3/2. Hence, the ground state isequivalent to that of the spin 1/2-3/2 alternating chain. Therefore, the ground state is aferrimagnetic phase with a spontaneous magnetization msp = 1 per unit cell. Between
these two phases, we carry out the inﬁnite-size
DMRG calculations for various conﬁgurationsof {Tl } to obtain the ground-state energy perunit cell. As a result, we ﬁnd that the groundstate undergoes a transition at λ = λc0 (δ) to aseries of ferrimagnetic phases with msp = 1/p
where p takes positive integer values. Thistransition is a ﬁrst-order transition for δ ̸= 0
with a discontinuous change mcsp in msp , while
no discontinuity is found for δ = 0.[2]The critical behaviors of mcsp and λc0 (δ)
around the critical point (δ, λ) = (0, λc0 (0))
turned out to be

(a)

0.8

0.79

0.78
0

0.15
c
msp

0.05

0.1

2/3

δ

(b)

mcsp ∝ δ 2/3 ,(2)λc0 (δ) − λc0 (0) ∝ δ 2/3 .

(3)

These
relationsanalytically.[1]

are

also

explained

0.1

References[1] K. Hida, J. Phys. Soc. Jpn. 91, 024706(2022).
[2] K. Hida, J. Phys. Soc. Jpn. 90, 054701(2021).

0.05

0
0

0.05

2/3

δ

0.1

Figure 3: δ-dependence of (a) λc0 (δ) and (b)mcsp .plotted against δ 2/3 .

310

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

Quantum simulation using Ising models
Yuya SEKI
Graduate School of Science and Technology, Keio University, Kanagawa 223-8522
consider the problem of finding the groundWe have studied properties of factorization
machine with annealing (FMA) [1] for quantum

state

simulations. The FMA method approximates

Hamiltonian consists of fermions. In thisthe Hamiltonian of the target system with Ising

project, we used three encoding methods:

model through a machine learning model calledbinaryfactorization machine (FM) [2]. Instead of

domain-wallfinding the low-energy states of the target

performance for each encoding method, we

system that is hard to analysis, we investigate

calculate the distribution of the accuracy of

low-energy states of the approximated Ising

solution of the developed method by runningspin system. In FMA procedure, we use Montethe developed method repeatedly. As a result,

Carlo methods to generate low-energy samples

we obtained the optimal region of parametersfrom the trained FM. In addition, the training

where the developed method provides aprocess of FM includes the backpropagation

solution that is close to the true ground statemethod for updating the model parameters inwith a high probability for each encoding

FM, which requires a linear algebra library formethod. Adjusting the parameters, we obtaincalculation. These calculations were performed

the ground-state of the hydrogen moleculeusing multithreading to reduce the execution

where relative energy error is less than 10!" .

time. FMA includes various hyperparameters

Also, we determined the proper encodingthat should be optimized to achieve highmethod from the obtained data. Among theaccuracy solutions. To search for optimal

encoding methods, one-hot encoding shows

hyperparameters efficiently, we performed thebetter performance than the others. This result

calculation

indicates that the sparsity of data expression

inparallelbychangingtheparameters.

of

the

hydrogenencoding,

one-hot

encoding.

molecule,

whose

encoding,Toandinvestigate

determines the performance of FMA.

References

This year we developed a new method thatcombines FMA with integer-encoding methods.

[1] K. Kitai et al.: Phys. Rev. Research 2

The developed method enables us to calculate(2020) 013319.

the[2] S. Rendle: 2010 IEEE Int. Conf. Data Min.

groundstateofquantum

systems

approximately. As a benchmark problem, we(2010) 995.

311

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

Large-scale molecular-dynamics simulation of silica melt and glass
under high pressure with ANN potentials by active learning
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

obtaining required training data through active

and the compression behavior of silica glass andlearning.

melt has attracted considerable attention in

An active-learning method can be roughly

various research fields of physical sciences. Due

divided into three stages. First, multiple ANN

to the strong covalent bonding, the melting point

potentials are decided for the same set of training

is extremely high.

data and MD simulations are made with each

It is difficult to perform

experiments on silica melt under high pressure,potential.

and theoretical calculations are a powerful

made for the atomic arrangement corresponding

research tool.

to the steps for which the calculation results in

In contrast to crystalline materials, non-

Second, ab-initio calculations are

the first stage differ significantly among the

crystalline materials do not have long-range

potentials.

order, and therefore, in order to obtain reliable

added as new training data for the next potential

information,

fitting.

itisnecessary

toperform

Finally, the ab-initio results aretheoretical calculations on a system as large as

The test calculations were performed for the

possible. There is a limit to the size of thesolidification process of silica melt under highsystem which can be calculated by ab-initio

pressure. The ANN potential was determined bymethods.

In this study, we have fitted thethe machine-learning method with the potential,

potential with a machine-learning techniqueforce, and stress data obtained by ab-initio(ANN potential) on the basis of the results of ab-

calculations for a temperature range of 300-

initio calculations, and have applied it for the

5,000 K at 60 GPa with 144 atoms. Active

calculations with an extended system. However,learning was performed three times to search forin our previous projects, it has been observed

the conditions under which calculations on silica

that the large-scale calculations of melt weremelt do not diverge at 60 GPa and 4,000 K for adifficult to stabilize and quickly diverged,

system expanded to about 30,000 atoms.

making them impossible to continue. In thisThe number of data was originally 15,000

project, we tested a new method of efficiently

and about 2,000 were added for each active

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

learning. Although there was improvement inlearning. This may be due to the fact that thethe differences in calculation results between

active learning process increases the likelihood

ANN potentials and in the number of steps atof mixing training data of unrealistic atomic

which divergence occurred in the expanded

arrangements. It was found that active learning

system up to the second active learning, there

enables us to delay, if not completely prevent,

was almost no improvement after the third

divergence.

313

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

Theoretical study of thermoelectric properties
in doped Fe2VAl: A weak-coupling approach
Kazutaka NISHIGUCHI
Graduate School of System Informatics, Kobe University
Rokkodai-cho, Nada-ku, Kobe 657-8501

1 Heusler compound Fe2 VAl

Fe (d)
V (d)

4.0

2.0

Energy (eV)

Heusler compound Fe2 VAl is one of the promising thermoelectric materials, which exhibits large
power factor P = σS 2 with σ and S being theelectrical conductivity and Seebeck coeﬃcient, respectively, although the thermal the dimensionlessfigure of merit ZT = σS 2 T/κ is not so large yetreflecting its large thermal conductivity κ.
Recently, it is experimentally observed in
doped Fe2 VAl (i.e., Fe2 V0.9 Cr0.1 Al0.9 Si0.1 andFe2.2 V0.8 Al1−y Siy ) as a weekly ferromagnetic material that ferromagnetic fluctuations enhance thethermoelectric properties such as S and P aroundthe Curie temperature (T c = 285 K). [1] It is not
only a desirable property for practical use in realistic (room) temperatures, but also a fundamental
and intriguing quantum phenomena where quantum fluctuations and thermoelectric eﬀects are entangled by electron correlations.
Motivated above, we have studied electronic
properties of Fe2 VAl using first-principles calculations based on the density functional theory (DFT) within the generalized gradient approximation (GGA). We performed DFT calculations with the Perdew–Burke–Ernzerhof (PBE)
exchange-correlation functional and the projector
augmented wave (PAW) method using The Vienna Ab Initio Simulation Package (VASP) [2] andQUANTUM ESPRESSO package. [3] One could
observe from the band structure and density ofstates (DOS) are shown in Figure 1 that the Fe andV d-orbitals around the Fermi level mainly contribute to the low-energy electronic structures.
We have also considered the dopant eﬀects on

6.0

0.0

−2.0

−4.0

−6.0

−8.0

W

L

Γ

X

W

K

Γ

a.u.

Figure 1: (Left) Band structure and (Right) DOS
of Fe2 VAl. The energy is measured from the Fermi
level. The components of Fe and V d-orbitals areprojected (by red and blue), respectively.

the doped Fe2 VAl to investigate its magnetic properties: the Cr and Fe substitution for V and the Si
substitution for Si were taken into account by the(2 × 2 × 2) supercell calculations. We could observe that the Cr and Fe substitution for V induce
the spin-polarized states and the magnetic momentarises at the Cr and Fe, respectively, while the electronic structure for the Si substitution for Si is stillparamagnetic state.
Furthermore, we have constructed the firstprinciples eﬀective model for Fe2 VAl using the
maximally localized Wannier functions with Wannier90. [4] We here found that the eﬀective modelhas to involve not only Fe and V d-orbitals mentioned above but also Al s- and p-orbitals in or-314

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

der to reproduce the original band structure aroundthe Fermi level since these orbitals highly entangle with each other and have large hybridization.
Therefore, we adopted a 19-orbital eﬀective model
for Fe2 VAl: 5(Fe − d) × 2 + 5(V − d) + 4(Al − sp).

+1/0

0/-1

Sc
+2/+1

+1/0 0/-1
Y+2/+1+1/0 0/-1La

2 122 Zintl Phase Compound
CaZn2 X2 (X = As, P)
In order to search for high-performance thermoelectric materials, the possibility of n-type dopingof CaAl2 Si2 -type Zintl phase compound CaZn2 X2(X = As, P) is explored using first-principles calculations based on DFT. [5] We consider n-type(electron carrier) doping of CaZn2 X2 with the following two situations: interstitial-site doping of alkaline earth metals AE (= Mg, Ca, Sr, Ba) andgroup 3 elements G3 (= Sc, Y, La), and G3 substitutional doping for the Ca site.
To see this, the defect formation energy of thesecharged states is evaluated within GGA using thesupercell approach. Among the considered possibilities, the interstitial-site doping of AE = Ca,
Mg or G3 = Sc, Y, and the G3 = La, Y substitutional doping for the Ca have been found to haverelatively small formation energies. In particular,
the formation energy of the La substitutional doping is found to be negative for both CaZn2 As2 and
CaZn2 P2 (Figure 2). This suggests that La can substitute the Ca site spontaneously and hence provide electron carriers, as far as these calculation
results are concerned. We have also found thatthe formation energies of the defects are smaller
for CaZn2 As2 than for CaZn2 P2 , which suggeststhat n-type doping is relatively easier for the former than for the latter.

+2/+1

Figure 2: The formation energy of the G3 substitutional doping for the Ca site Eform against the
Fermi energy ∆EF = EF − EV for CaZn2 As2 . The
group 3 elements G3 (= Sc, Y, La) are chosen as
dopant D. The dotted line denotes the position of
the conduction band minimum (EF = 0.408 eV).
tational Materials Science 6 (1996) 15; Phys.
Rev. B 54 (1996) 11169.
[3] P. Giannozzi et. al.: J. Phys.: Condens. Matter 21 (2009) 395502; J. Phys.: Condens.
Matter 29 (2017) 465901; J. Chem. Phys. 152(2020) 154105.
[4] I. Souza, N. Marzari, and D. Vanderbilt:
Phys. Rev. B 65 (2001) 035109; A. A.
Mostofi et. al.: Comput. Phys. Commun. 178(2008) 685; A. A. Mostofi et. al.: Comput.
Phys. Commun. 185 (2014) 2309; G. Pizzi
et. al.: J. Phys.: Condens. Matter 32 (2020)
165902.
[5] K. Nishiguchi, M. Ochi, C. H. Lee, and K.
Kuroki: J. Phys. Soc. Jpn. (to be published).

References[1] N. Tsujii, A. Nishide, J. Hayakawa, and T.
Mori: Sci. Adv. 5 (2019) eaat5935.
[2] G. Kresse et. al.: Phys. Rev. B 47 (1993)
558; Phys. Rev. B 49 (1994) 14251; Compu-

315

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

Spontaneous formation 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
in inversion-symmetric systems where the DM
interaction is absent [1, 2]. In the present
work, we have attempted to explore theoretically the possibility of multiple-Q orders in
the inversion-symmetric systems in the absence of external magnetic ﬁeld. For this purpose, we have applied the molecular spin dynamics (MSD) method [3] to the triangularlattice single-band Hubbard model.
The MSD method is based on the functional
integral method for the spin ﬂuctuation theories and the isothermal molecular dynamics
method. The method allows us to ﬁnd automatically the magnetic structure of a large system with thousands of atoms in a unit cell at
ﬁnite temperatures. Starting from the Hamiltonian expressed in terms of the locally rotated coordinates and by adopting the static
approximation to the functional integral technique, the MSD method reduces to the generalized Hartree-Fock approximation at the
ground state.
In the numerical calculations the most timeconsuming process is the magnetic force calculation at each time step, where the local electronic structures are calculated in real

space by means of the recursion method. We
have adopted the automatic parallel calculation scheme and found it eﬀective in saving
both computing time and CPU resources.
We have performed magnetic structure calculations on a hexagonal supercell with 972
lattice points, which is embedded in a large
cluster consisting of 6 such supercells, each ofwhich 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 structureschanging the Coulomb interaction strength
U/t (U/t = 3.0 ∼ 8.0) and the electron number n (n = 1.10 ∼ 1.60) along
the antiferromagnetic-ferromagnetic boundary. We have found that for U/t ∼ 3.5 andn = 1.48 ∼ 1.53 the dominant state is 3Qstructure with Q1 = (0.3889, −0.2245)2π/a,Q2 = (0.3889, 0.2245)2π/a, and Q3 =(0.0, 0.4491)2π/a (with a being the lattice constant), which is superimposed with satellite 3Q
structures with the Q vectors pointing in thesame direction ( in the direction of three axesof the triangular lattice).

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.

316

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

Molecular Dynamics Simulations of Reversibly Cross-linked Rubbers
Yusuke Yasuda
Research Center for Computational Design of Advanced Functional Materials.

National Institute for Advanced Industrial Science and Technology, Ibaraki 305-8568.
that

Reversibly cross-linked rubbers, whose

dissociationsynchronizes

tostresscross-links are made of reversible bonds, suchretardation. Self-healing test was conducted via

as hydrogen bonds, ionic bonds, and dynamic

equilibrium calculation at high temperaturecovalent bonds, have self-healing property andafter loading and unloading. In the elongation

toughness. Though many researchers haveafter equilibrium, stress recovered completely,

studiedmolecular

and our model can reproduce the feature ofmechanisms of these properties are still unclear

self-healing. (Fig.(b)) From these results, we

due to the difficulty of the observations ofhave succeeded in constructing coarse-grained

molecular dynamics. To understand thosemodel of reversibly cross-linked rubbers.

thosematerials,

thematerials, it is needed to study the relationshipbetween structure, dynamics, and properties ofthese rubbers from several research methods.
In this work, we constructed coarse-grainedmodel for reversibly cross-linked rubbers basedon Kremer-Grest model [1], by introducing
associating and dissociating condition (Fig.(a))
changing bond dissociation energy and bondpotential shape. Calculation was conductedusing molecular dynamics simulation software
packages OCTA and LAMMPS on ISSP
Supercomputer

System.

Fromequilibrium

simulation, we found that association constants
and bond lifetime are controlled by bondFig. (a) The schematic illustration of coarse-

dissociation energy and bond potential shape.

grained dynamic bonds, (b) loading-unloading

Mechanical simulation was conducted via

curves of the model.

loading-unloading test conducted in the rangeof elongation ratio λ= 1 – 4. From loading-

References [1] K, Kremer, G. S. Grest; J. Chem.

unloading simulation, stress-extension rate

Phys., 92, (1990) 5057.

curve shows stress retardation and specifichysteresis (Mullins effect, Fig.(b)) which is

Acknowledgement

observed experimentally in reversibly cross-

This

work

is

partially

supported by a project JPNP18016, commissioned

linked elastomer. We also found that theby the New Energy and Industrial Technology

dynamic bonds dissociate by elongation, andDevelopment

317

Organization(NEDO).

Activity Report 2021 / 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

318

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

Molecular dynamics analyses of ion migration at grainboundaries in solid state electrolyte
Ryo KOBAYASHI
Department of Applied Physics, Nagoya Institute of Technology,
Gokiso-cho, Showa-ku, Nagoya, Aichi 466-8555
All-solid-state battery (ASSB) is considered

GBs. From the simulation and local ion-flux

as a candidate of next generation Li-ion

analyses[1], we showed that, in poly-crystalline

batteries because it uses a solid-state electrolytematerials, the ions migrate towards the lower

instead of liquid-state electrolytes that are usedreaches of grain and go through spots at GBs,

in conventional Li-ion batteries. The key toand the high-flux spots contain migration paths

achieve high-performance in ASSB is toas low energy as the bulk ones.

reduce/control the resistivities at interfaces

In this study, we have developed a newsuch as electrode/electrolyte interface and graintechnique to obtain local ion-flux information

boundary (GB) inside electrolyte. For thisfrom non-equilibrium MD implemented in our

purpose, we need deeper understandings about

own code, nap[2]. And several MD simulation

the atomic structure and precise mechanisms ofruns and the local-ion flux analyses were

ion migration of the GBs, and the goal of thisperformed on the ISSP supercomputers.

project is to obtain new insights about them byusing large-scale molecular dynamics (MD)simulation.

References

In order to obtain atomistic insights on Liion

migrationmechanism

at

GBs,

we

[1] R. Kobayashi, K. Nakano, M. Nakayama,

performed non-equilibrium MD simulation of

Acta Materialia, 226 (2022) 117596.

large-scale poly-crystalline Li2Zr(PO4)3 (LZP)[2] Nagoya Atomistic-simulation Package (nap)system containing half a million atoms andhttp://github.com/ryokbys/nap

extracted local ion flux information around.

319

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

Electronic structure in mismatched multilayer systems:
evaluation and application of the modeling scheme
Toshikaze KARIYADOInternational Center for Materials Nanoarchitectonics,
National Institute for Materials Science, Tsukuba, Ibaraki 305-0044
After the report of superconductivity inwhich gives us an intuitive understanding of thetwisted bilayer graphene, artificially stacked

physics behind the anisotropic band flattening.

two-dimensional (2D) materials with moiré

The first-principles calculations have been

pattern induced by some mismatches between

done using Quantum Espresso package, which

the layers attracts much attention. Generically,is adapted to the parallel computation andthe moiré pattern affects the electronic bandsuitable to use the supercomputing resources

structure, and can be used to realize uniqueefficiently.

band structures. Then, designing interesting

Not

only

calculating

electronic

band

band structures in artificially stacked systems is

structures itself, we also tried to extract somea current urgent topic to explore.

characteristic response tied to the unique bandIn this project, we have analyzed electronic

structure, focusing on the orbital magnetic

structures of twisted bilayer GeSe, where thesusceptibility. For that, we first derived a new

anisotropic band flattening is expected. Here,formula to calculate the orbital magnetic

the anisotropic band flattening means thatsusceptibility,

originally 2D band dispersion is squeezed into

validity using several types of toy models by

quasi 1D one in twisted bilayers. To analyze theour homemade program [2]. In order to accessband structure, we derived an effective model

the low temperature regime, we need fine meshfor twisted bilayer GeSe. There are multiple

in 𝑘-point sampling, where the parallelization

important quantities in the effective model suchby MPI is applied in our homemade program.

andnumericallytestedits

as the effective mass for each layer or theeffective interlayer tunneling. We obtain these

References

quantities making use of the first-principles[1] M. Fujimoto and T. Kariyado: Phys. Rev. B

calculations. Our analysis clarifies the relation

104 (2021) 125427.

between the crystal structure of GeSe and the[2] T. Kariyado, H. Matsuura, and M. Ogata: J.

anisotropy in the effective interlayer tunneling,

Phys. Soc. Jpn. 90 (2021) 124708.

320

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

Compact spin qubits using the common gate
structure of ﬁn ﬁeld-eﬀect transistors
Tetsufumi Tanamoto
Department of Information and Electronic Engineering, Teikyo University
Utsunomiya 320-8511, Japan.

By using the ISSP supercomputer, we numerically investigate a transport properies of
the compact spin-qubit system embedded into
the common multi-gate FinFET transistors [1],where all gates are electrically tied togther as
the common gate. The quantum dots (QDs)
as qubits are coupled with their nearest Fin
conducting channels. The Hamiltonian of the
QDs and the channel is given by the tunneling
Hamiltonian:tance gyy (ω) is calculated from the Kubo for1
R
mula given by gyy (ω) = − iω[ΦRyy (ω) − Φyy (0)],iRΦyy [t] = − h̄V θ(t)⟨Jy (t)Jy (0) − Jy (0)Jy (t)⟩.
The current operator Jyi of the ith channel∑
is given by Jyi = (eh̄/(m∗ L)) ki ki c†ki s cki s ,where L is the channel length and the summation of ki is carried out over the channel.
(3)
The expression for gyy is

2e2
4[∆2 + δ 2 ]2,kdh[(∆2 − 2s33 ∆ − δ 2 )2 + 4∆2 Γ23 ]2
∑ ∑
(2)
†
†
†
c
H0 = E2 d2s d2s + E4 d4s d4s +Ek ci ki s ki s
2where k1 = 1, k2 = πne2 W , ∆ = (2EkF −i=1,3,5 ki ,s∑†ESL − ESR − s11 − s55 )/2, and δ = ESL −∗ †+[Vk1 ck1 s d2s + Vk1 d2s ck1 s ]ESR . ne2 is the number of the carriers per nm2 ,
k1
∫
∑
2
+
[Vk3 c†k3 s (d2s + d4s ) + Vk∗3 (d†2s + d†4s )ck3 s ] and sij ≡ |Vtun (ki )| /(Eki − Ekj ) is the selfenergy. Γi ≈ 2π|Vtun (ki )|2 ρF (ρF is the density
k3∑+[Vk5 c†k5 s d4s + Vk∗5 d†4s ck5 s ],(1) of state at Fermi energy EF , and Vtun is theoverlap of wave functions between the channel
k5
and the QDs in the tunneling Hamiltonian).
where the channels are numbered as 1, 3, andFigure 1(b) shows the conductance gyy of
5, and the two QDs are numbered as 2 and(i)the summation of the three current lines gyy
4. dis and ck,s are the annihilation operators(i ∈ 1, 3, 5) as a function of the energy levels ofof the QD i and the conducting electrons in
the two QDs, where ESL and ESR are either
the channel, respectively. The qubit states areES↓ or ES↑ . We can observe a double-peakdetected by the channel currents. The conducstructure around the Fermi energy where ESL
is close to ESR but ESL ̸= ESR . The symmetDEΓ Γ Γ Γric case δ = 0 gives the conventional resonant*DWH
tunneling form g = 4/[(∆ − 2s33 )2 + 4Γ23 ]2 .
1

ESL

3

3

(3)
gyy
=

5

ESR

VXEVWUDWH

References

ESL /EF
ESR /EF

Figure 1: (a) Three current lines. (b)gyy as
functions of energy levels of two quantum dots.

[1] T. Tanamoto and K. Ono, AIP Advances
11 (4), 045004 (2021).

321

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

Change in period of ordinary earthquakes due to theinteraction with slow earthquakes
Takehito SUZUKI
Dept. Phys. Sci., Aoyama Gakuin University
5-10-1, Fuchinobe, Chuo-ku, Sagamihara, Kanagawa 252-5258

Two qualitatively diﬀerent behaviors of
earhtquakes have been known. One is the ordinary earthquake, 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 ordinary
earthquakes after several repetitions. Whether
this transition occurs or not is a fatal problem
to the human society, and we should clarify thecondition for the transition to occur.

constant is kp . As widely known, this systemrepeats the dynamic slip and cessation. The
healing of pores in the cessation time is also
taken into account. However, the choice of thehealing law does not induce qualitative changes
in the result below.
We neglect the heat and ﬂuid diﬀusions during the dynamic slips. The governing equations for the ﬂuid pressure, porosity, and temperature in the slip zone during the dynamic
slip are given by

The fault rocks can be considered as poroelastic media, including pores inside it. The
pores are usually considered to be ﬁlled withwater. 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 around the contact area, called slip zone,
which has a ﬁnite width. During the dynamic
slip, the generation of pores [frictional heating] reduces [raises] the ﬂuid pressure in theslip zone, inducing the increase [reduction] in
the friction stress and the reduction [increase]in the slip velocity [1, 2]. Therefore, the slow
and ordinary earthquakes are dominated by
the generation of pores and the frictional heating, respectively.

ṗ = C1 Ṫ − M0 ϕ̇,The spring-block model and the interaction
between three quantities are employed to analyze the slow-ordinary transition. The single
block is assumed here. The upper substrateis drived with the constant speed Vp , and it
also loads the block via a spring whose spring

(

(1)

)

ϕ
,
ϕ∞

(2)

Ṫ = −C2 (σn0 + p)u̇,(3)ϕ̇ = α0 u̇ 1 −respectively, where C1 and C2 are positive constants. From here and below, quantities not
explained are summarized in Ref. [2]. From
Eq. (2), we obtainϕ = ϕ∞ − (ϕ∞ − ϕ0 )e−α0 u/ϕ∞ ,(4)

where ϕ0 is the porosity at the intant of theslip onset. Note that this value need not equal
to the value at time t = 0. Additionally, from
Eqs. (1)–(4), we obtain
(

(

)

)

ϕ0
ṗ =
1−
α0 e−α0 u/ϕ∞ u̇,ϕ∞(5)
where γ = C1 C2 > 0. From this representation, we can write p in terms of u to obtain−C1 C2 (σn0 + p) − M0

(

p=

322

M ′ α0σn0 + p0 − α0 0
ϕ∞ − γ)M ′ α0 −α0 u/ϕ∞ 0e−γu + α0 0e−σn ,ϕ∞ − γ(6)

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

where M0′ ≡ M0 (1 − ϕ0 /ϕ∞ ) and p0 is the ﬂuid
pressure at the instant of slip onset.
We now consider the energy conservation
law. First, note that the loading stress fromthe upper substrate must be equal to themacroscopic static friction stress at the instantof the slip onset, which leads tokp Vp ti = −µstat (σn0 + p0 ),(7)

where ti is the time when the slip initiated.
Using Eq. (7), the energy stored in the upper
spring just before the dynamic slip is given by1
1 µ2stat 0(σn + p0 )2 .
E1 = kp (Vp ti )2 =
2
2 kp

(8)

Moreover, the energy stored in the upper
spring after the cessation is written as (we assume that the upper spring does not move during the dynamic slip)1
kp (Vp ti − uf )22()2
1µstat 0(σn + p0 ) − uf , (9)
kp −
2
kp

E2 =
=

where uf is the ﬁnal slip distance for the single
slip. Furthermore, the work done to the blockby the friction stress during the dynamic slip
is obtained using Eq. (6):∫ uf (

Efric =
= −µslid[= −µslid1γ∫ uf [(
0(0)−µslid (σn0 + p) du

M 0 α0σn0 + p0 − α0
ϕ∞ − γ)e−γu]

M ′ ϕ∞+ α0 0(1 − e−α0 uf /ϕ∞ ) = 0.(11)−γϕ∞We deﬁne F (uf ) as the left hand side of Eq.
(11). The condition F (uf ) = 0 corresponds tothe condition where the elastic energy change
is equal to the energy loss due to the friction
during the single slip. It is also important that
F (uf ) includes p0 and ϕ0 . As mentioned, these
initial values are not values at t = 0, but thevalues at the slip onset of each slip.
Notably, F (uf ) is the double-well form, andthe ﬁrst well is much shallower than the second one. If F (uf ) = 0 has a solution between the wells, the slow earthquakes occur
because we can consider that the slip is inhibited. On the other hand, if F (uf ) = 0 does
not have a solution there, the solution jumps
to almost two orders larger one. This behavior
is interpreted as the ﬁrst-order phase change.
The slow-ordinary earthquake transition canbe understood as the ﬁrst-order phase change.
Additionally, we can conclude that this phasechange is governed by p0 and ϕ0 .
The analytical result described here was conﬁrmed by ISSP supercomputer. The numerical study can treat the BK model with several
blocks, and we conﬁrmed that the results here
can give useful implications to such a severalblock case. Detailed and comprehensive numerical treatments will be a potential future
work.

]

M0 α0 −α0 u/ϕ∞
+ α0
e
du
ϕ∞ − γM0 α0
σn0 + p0 − α0
ϕ∞ − γ)(1 − e−γuf )References][1] T. Suzuki, and T. Yamashita: J. Geophys.
M0 ϕ ∞+ α0(1 − e−α0 uf /ϕ∞ ) .(10)Res. 119 (2014)
ϕ∞ − γTherefore, from the energy conservation law,
we found an equation which uf must satisfy:[2] T. Suzuki: Phys. Rev. E 96 (2017)1
kp u2f + µstat (σn0 + p0 )uf
2
[ (
)
1
M0′ α00
σn + p0 − α0
(1 − e−γuf )−µslidγ−γ
ϕ∞323