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

3.4

Cooperative Phenomena
in Complex Macroscopic Systems

220

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

Tensor-Network Renormalization-Group Studyof Critical Phenomena
Naoki KAWASHIMA
Institute for Solid State Physics,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
The objective of the project is to develop
newmethods

forstudyingtheIn [2], we studied the J-K- model, a model

criticalthat may represent some aspect of thephenomena, both classical and quantum, andcompound -RuCl3, which is expected to be

apply them to strongly correlated systems.

close to the Kitaev spin liquid and presents

Especially, we aim at improving the real-spacesome features hinting a non-trivial topologicalrenormalization group (RSRG) method based

structure of its state. We discovered that theon the tensor network representation with thechiral spin liquid phase proposed by ahope for fuller understanding of the structure ofpreceding study is realized in a narrower region

the scaling dimensions of 2+1 or 3 dimensional

than predicted. In addition, we discovered two

fixed-point theories.

more non-magnetic phases that our calculation

In [1], we proposed a method for obtaining

suggested possesses also non-trivial topologicalan efficient low-rank approximation of a given

structure.

tensor. This was motivated by the conventional

In [3], we studied the spin-1 Kitaev model.

RSRG methods such as the loop-TNR andThis was motivated by the preceding work byentanglement filtering. In those methods, theLee, Kaneko, Okubo and I on the conventional

crucial part is removal of redundancy in the

spin-1/2

tensor-network representation arising from

relationship between the gapless Kitaev spinloops. They can be regarded as a special form

liquid and the classical loop gas model. In [3],of low-rank approximation of a given tensor. It

in contrast to the S=1/2 case, we discovered themight then be useful also for more general class

gapful spin liquid with Z2 structure of theof tensors, e.g., the image data in the form ofquasi-particles.

case,

in

which

we

discovered

tensors. The proposed method successfully

In [4], we improved the high-order tensor

identified the redundant structure of therenormalization group (HOTRG), one of theentanglement naturally formed in the process of

conventional RSRG besed on the tensor-

thenetwork representation. In the HOTRG, wetensorrenormalization

group

(TRG)

calculation.

bind two parallel bonds together to form one

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

renormalized bond in much the same spirit as

characterization

the Migdal-Kadanoff RG. In doing so, we needphenomena has been a tough problem even fora projection operator for the rank reduction. In

DMRG.

its simplest implementation of the HOTRG,

characterize various phases including symmetry

only the local environment is taken into

protected topological phases by the 2nd

account in computing the projection operator,

cohomology group represented by the resulting

while it would obviously be better to optimize

fixed-point tensors. Critical behaviors between

the projector in more global environment.

the phases are also accurately characterized.

In

of

our

the
study,

quantum
wecritical

successfully

Several methods for better optimization has
been

proposed

at

the

cost

of

References

highercomputational complexity. The method we[1] Hyun-Yong Lee, and Naoki Kawashima: J.

proposed is based on the corner transfer-matrix

Phys. Soc. Jpn. 89, 054003 (2020)

technique and reduces the computational[2] Hyun-Yong Lee, Ryui Kaneko, Li Erncomplexity without serious degrading of theChern, Tsuyoshi Okubo, Youhei Yamaji, Naoki

accuracy.

Kawashima and Yong Baek Kim: Nature

In [5], we applied the method of variational

Communications 11, 1639 (2020).

uniform matrix-product-state (VUMPS), an[3] Hyun-Yong Lee, Naoki Kawashima and

MPS-based

Yong Baek Kim: Phys. Rev. Res. 2, 033318

methodforone-dimensional

quantum systems to the 2-leg ladder system(2020).

with[4] Satoshi Morita and Naoki Kawashima: Phys.

four-bodyinteraction,

i.e.,

a

one-

dimensional analogue of the JQ model that

Rev. B103, 045131（1-6） (2020)shows a deconfined critical phenomena. While[5] Takuhiro Ogino, Ryui Kaneko, Satoshi

DMRG is known to be a very powerful

Morita,numerical

Kawashima: Phys. Rev. B 103, 085117 (2021).

dimensional

tool

in investigating thequantum

systems,one
direct222

Shunsuke

Furukawa

and

Naoki

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

Molecular Dynamics Simulations of Bio-Inspired
Materials
Wataru Shinoda
Department of Materials Chemistry, Nagoya University
Nagoya 464-8603

We have studied the surface hydrophobicity
of self-assembled monolayers (SAM) [1] using
massive parallel computer simulations. The
surface hydrophobicity was investigated for arange chemically functionalized (–CF3 , –CH3 ,–OCH3 , –OH) SAMs as illustrated in Fig. 1.
These SAMs have many industrial and biological applications such as for making well-defined
non-fouling surfaces, cell supports, receptors,
sensors, microarrays, separation devices, andfor surface reactions [1, 2]. A good understanding of surface properties of these SAMs can allow us to further modify and control them forthe specific uses.

Figure 1: Surface hydrophobicity of chemically functionalized (–CF3 , –CH3 , –OCH3 , –OH) self-assembled monolayers. (a) Structureof a water droplet on the fluoroalkyl SAM surface (b) The formation of a cavity near the
SAM-water interface.
The study was conducted with all-atom(AA) simulation model using the LAMMPS

223

molecular dynamics (MD) simulation package[3]. The systems were modeled using the optimized potentials for liquid simulation (OPLS)
force field [4]. SPC/E rigid-body water model
was used for the water molecules [5]. A wide
range of packing densities and chain lengths
were considered in the study to understandtheir eﬀects on the surface hydrophobicity.
Analyses of the MD trajectories were done using the in-house codes. The complete detailsof the computational method can be found in
Ref. [6, 7].
The surface hydrophobicity was investigatedby means of (1) the structure, (2) the contact
angle of a water droplet on the surface, and(3) the cavity formation free energy. In the
structure, we mainly focused on density distribution, hydrogen-bonding, chain flexibility via
root mean square displacement (RMSD), androughness of the SAM surfaces. The free energy of cavity formation was calculated fromindirect umbrella sampling methods [6, 7],
in which used massive parallel computing resource. The results were found quite consistent with experiments. Through these studies,
we demonstrated a systematic path to examine the hydrophobicity of such synthetic surfaces. For the hydrophobic surfaces such as forthe fluorinated chains, it was found that thesurface hydrophobicity is mainly governed bypacking density rather than the chain length.
It increases on reducing the packing densityof SAMs. The SAM flexibility causes an enhancement in the surface hydrophobicity inde-

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

pendent of chemical functional groups. These
information would be certainly useful for thematerial design and modifying the surfaces.
Additionally, we developed a coarse-grained(CG) model for the semi-fluorocarbon diblocks
and studied the hemimicelle formation of thesechains at the air-water interface using ISSP
supercomputer. These systems contain morethan 100,000 particles even in the CG model.
We showed a mechanism of micelle formation
for these diblocks and successfully explained
the non-coalescence behavior of micelles, which
was a puzzle for the experimentalists [8]. It
was found that an unique alternating P-phaseforms in between the micelles, which has lateral and vertical dimensions of the nanoscaleorder (see Fig. 2). We believe that this Pphase could be reason for the non-coalescencebehavior of micelles under strong compression.

Figure 3: Adsorption of protein on synthetic
SAM surfaces.
time scale. In this study, we will also consider
the adsorption of more larger systems such as
virus capsids and organelles.

References[1] J. C. Love, L. A. Estroﬀ,J. K. Kriebel, R.
G. Nuzzo, and G. M. Whitesides: Chem.
Rev. 2005, 105 (2005) 1103.
[2] J. E. Raynor, J. R. Capadona, D. M.
Collard, T. A. Petrie, and A. J. Garcia:
Biointerphases 4 (2009) FA3.
[3] S. Plimpton:(1995) 1-19.

J. Comput. Phys. 117[4] W. L. Jorgensen, D. S. Maxwell, and J.
Tirado-Rives: J. Am. Chem. Soc. 118(1996) 11225.
Figure 2: The hemimicelle formation of semifluorocarbons at the air-water interface andthe structure of P-phase existing in between
the hemimicelles.
More recently, we studied the amino acids
and protein adsorption on the SAM surfaces(see Fig 3). We are in the process of summarizing results of these simulations for the publication. We found ISSP supercomputer really
useful for conducting simulation of such largescale complex systems. We are further looking
to simulate protein adsorption on the polymers
using a coarse-grained model to cover extensive

224[5] H. J. C. Berendsen, J. R. Grigera, andT. P. Straatsma: J. Phys. Chem. A. 91(1987) 6269.
[6] H. O. S. Yadav, A. -T. Kuo, S. Urata andW. Shinoda: Langmuir 35 (2019) 14316.
[7] H. O. S. Yadav, A. -T. Kuo, S. Urataand W. Shinoda: J. Phys. Chem. C 124(2020) 14237.
[8] H. O. S. Yadav and S. Harada, A. -T. Kuo,
S. Urata, and W. Shinoda: Mol. Phys.
(2021) e1910355.

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

Construction of data assimilation in materials science
Munehisa MATSUMOTO
Institute of Materials Structure Science, High Energy Accelerator Research Organization
Oho 1-1, Tsukuba 305-0801

We have been developing methods to use
both of experimental data and theoretical data
to better characterize and predict good magnetic materials [1, 2]. Experimental data takenfrom actual materials and corresponding computational data constitute an integrated data
model. The desired properties of materials areoften well defined and it now looks as though it
is getting more and more feasible to implementa scheme for rational materials design based on
the integrated data among various experimental probes and ab initio materials data. However, during the course of our developments weare rather finding the relevance of stochastic
aspects of the materials design.
The problem is caused by multiple requirements which are often simultaneously imposed on the target materials. For permanent
magnets, a candidate compound should havestrong magnetization, magnetic anisotropy,
high Curie temperature, and sufficient structure stability. Unfortunately it is hardly thecase that all of these can be maximized witha particular chemical composition. Trade-off’s
are encountered and we have to work with agood compromise. Besides, another problem
is imposed by the multi-scale nature of materials. The material is put into practical use
in a macroscopic world in non-equilibrium andat the moment it is hardly possible to keep all
those extrinsic parameters under control in themodeling.
A starting point is to restrict the working
parameter space and implement an optimization scheme therein referring to the given setof requirements for the target material [2]. An-other solution which potentially has the flexibility to deal with the effect of the extrinsic parameters is a stochastic sampling over the landscape of the parameter space. A set of equivalently preferable parameter points can be sampled out and proposed as a set of candidate
materials to meet the given needs. We notethat the external needs can vary depending onsocial and economical trends. Even some unexpected requirements can pop up depending on
how the material is used in practical applications. These are not predictable, which points
to an unavoidable stochasticity in the materials design. Thus it seems to be a good idea toconsider the materials design as a stochastic
process working on the integrated data spaceconstructed over various experimental probes
and theoretical calculations. We adapt thedata assimilation techniques in weather forecasts into the materials science. Developments
are in progress and will be reported in detailelsewhere.
The author benefited from the online schoolfor data assimilation in the autumn 2020 [3].

References[1] MM, T. Hawai, K. Ono, Phys. Rev. Applied 13, 064028 (2020).
[2] Y. Harashima, K. Tamai et al.,Phys. Rev. Materials 5, 013806 (2021).
[3] http://www.data-assimilation.
riken.jp/jp/events/riken_da_tr_
2020/index.html

225

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

Data integration between experiments andcalculations: the new stages
Munehisa MATSUMOTO
Institute of Materials Structure Science, High Energy Accelerator Research Organization
Oho 1-1, Tsukuba 305-0801

We have been working on magnetic materials made of transition metals and rareearth elements. Since magnetism is a relatively low-energy phenomenon coexisting withhigh-energy properties of conducting electrons
in the metallic parts, fundamental problems
caused by the dichotomy between localized
magnetic moments and delocalized electrons
is already implied. On top of this, we encounter the multiple-scale problem where theintrinsic properties determined on the basis of
atomic-scale structure is combined with extrinsic properties controlled by the microstructure
to yield the macroscopic properties of magnetic
materials. They both involve different layersof physics and we have accordingly used combined approaches such as LDA+DMFT [1] toaddress the magnetism of localized f -electrons
embedded in the metallic environments. Herea quantum impurity problem is solved withquantum Monte Carlo method within the dynamical mean field theory (DMFT), and therealistic electronic structure is addressed withdensity functional theory based on local density approximation (LDA). In an analogous
spirit, we have worked on a computational
scheme where the electronic structure calculations for metallic ferromagnets and the input structure information taken from the Rietveld analysis of neutron diffraction data proceed via feedback to each other, which is tentatively termed “LDA+Rietveld” [2]. Here theexperimental data play the role of an effective
self-energy in the language of LDA+DMFT.

226

Both approaches can be considered as developments for data integration of different types.
The LDA+DMFT was originally developed forthe simulation of strongly correlated electron
systems which may be considered as a prototype to bridge over multiple scales in thetarget systems, which is often encountered inphysics and materials science. While the pastfiscal year 2020 was dominated by the serious
pandemic problem, remote working and online learning tools have significantly improved
which helped us to import several techniquesfrom experts outside of our conventional discipline. This included the data assimilation
techniques in weather forecast where unpredictable phenomenon coming from nonlinear
dynamics can be put under numerical control
by combining a limited number of observatory
data and relatively abundant simulation data.
This methodology adapted into the materials
science might help the upcoming materials design and discovery to be done in an accelerated
and predictable way, prospectively including
a more comprehensive data integration withsample-fabrication processes. Further developments are in progress.

References[1] MM, Phys. Rev. Materials 4, 054401(2020).
[2] MM, T. Hawai, K. Ono, Phys. Rev. Applied 13, 064028 (2020).

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

Co-equivalent of Ce in 4f -3d intermetallics
Munehisa MATSUMOTO

Rare-earth permanent magnets based on
Nd-Fe-B alloys have been widely used and expected to be even more of the key material in
the upcoming decades. The weak anisotropy
field at high temperatures and the relatively
low Curie temperature of the main-phase compound, Nd2 Fe14 B, sometimes require supplementary rare metals in order to put the material in practical use. Since such rare metals
are expensive and not always straightforwardly
available, we have been in urgent quest forgood ferromagnets on which permanent magnets may be developed in a cost-effective, robust, and sustainable manner. In this regard,
clarification of the potential utility of Ce canbe of a help. Indeed from the viewpoint ofsolid state physics, many interesting phenomena associated with Ce, such as valence fluctuations, heavy fermions, and unconventional
superconductivity, have been known. Some of
the fundamental ideas can help in answering
the problems posed from material science.
We have looked into the potential utility of
Ce for most of the representative propertiesof a ferromagnet to qualify as a constituent
of a permanent magnet, namely, magnetization, magnetic anisotropy [1], Curie temperature, and structure stability. Among all theseproperties, we have observed that the largest
merits actually seem to originate in the valence fluctuations, in particular the contribution from Ce4+ states to the bulk properties.
One of such aspects can be illustrated with
an effect of Ce substitution in SmFe12 [2]: wecan computationally get an analogue of the
Slater-Pauling curve with the Ce-substituted
ferromagnets. The magnetization of the Febased ferromagnet shows a peak as a function

227of the concentration of dopants, typically Co.
Remarkably, the same thing can happen evenwhen the dopant is Ce.
Obviously it is hard to observe this phenomenon in real experiments. There seems
to be a strict constraint on the inter-atomic
distance between the rare-earth sublattice andthe transition-metal sublattice in order for the
delocalized 4f -electron to contribute to thebulk magnetization that is mostly made of 3delectrons. Especially for Nd2 Fe14 B with therelatively inflated lattice, extracting the meritof Ce with the Slater-Pauling analogue is notstraightforward. Nevertheless there has beensome actual observation [3] confirmed bothexperimentally and computationally. Several
findings as exemplified by these cases point tosome hope for promising ferromagnets with reduced costs. So far we identified and understood several singular cases. Systematic quest
to pop up further candidates on the basis ofdata assimilation between theory and experiments are in progress.
Helpful discussions with K. Saito, T. Ueno,
and H. Shishido through the JSPS KAKENHI
15K13525 project are gratefully acknowledged.

References[1] H. Shishido et al., arXiv:2103.10202.
[2] MM, T. Hawai, K. Ono, Phys. Rev. Applied 13, 064028 (2020).
[3] X. Tang et al., Acta Materialia 175, 1(2019).

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

Finite temperature properties of frustrated systems
Tsuyoshi OKUBO
Department of Physics, The University of Tokyo
7-3-1 Hongo, Bunkyo-ku, Tokyo 133-0033

A variety of exciting phenomena, including non-collinear magnetic order, magnetization plateaus, and spin liquids, occur in thefrustrated spin systems. Frustrated interactions often appear as geometrical frustrations
based on the antiferromagnetic interaction on
triangular units. Another origin of the frustration might be the competition of several interactions, such as the J1 -J2 model on the squarelattice. The honeycomb lattice Kitaev model
could be considered as an example of the latter frustration [1]. In the Kitaev model, theinteractions are Ising type, Siγ Sjγ , and the spincomponent γ = x, y, z is determined from thedirection to the neighboring site. The Kitaev
model is exactly solvable, and its ground stateis known to be a spin liquid without long-range
magnetic order. Such Kitaev interaction could
appear in the realistic situation through thestrong spin-orbit interaction [2]. α-RuCl3 is
an example of such a Kitaev compound. In
this compound, the ground state is not the
spin liquid: it is a magnetically ordered state.
However, when we apply a magnetic ﬁeld, thecompound looks show the spin-liquid-like behaviors. In particular, a half-integer thermal
Hall conductivity was observed at ﬁnite temperature [3], indicating the appearance of Kitaev spin liquid.
In this year’s project, we investigated ﬁnite
temperature properties of the Kitaev modelwith oﬀ-diagonal interactions motivated by the
recent experiments on a Kitaev compound.
The model Hamiltonian is given as∑h ∑ x(S + S y + S z ), (1)H=
Hγ − √3 iγ∈x,y,zwhere for γ = z,∑
[KSiz Sjz + Γ(Six Sjy + Siy Sjx )Hz =⟨i,j⟩z+ Γ′ (Siz (Sjx + Sjy ) + (Six + Siy )Sjz )].
(2)
In the cases of γ = x, y, we consider similarinteraction with cyclic rotation of x, y, z components. To numerically calculate the physicalquantities at a ﬁnite temperature, we represent the system’s density matrix as a tensor
network and optimize its elements.
Firstly, we considered pure ferromagnetic
Kitaev model (K = −1, Γ = Γ′ = 0 and
h = 0) in the thermodynamic limit. In ourprevious calculation using the inﬁnite tensor
product operator (iTPO) tensor network, wehave found that the accuracy at the low temperature is not suﬃcient; we could not quantitatively reproduce the low-temperature speciﬁc heat peak previously calculated by quantum Monte Carlo [4]. In this year, we considered a cluster optimization of the tensor to improve the accuracy. In the standard optimization based on the imaginary time evolution,
we perform Suzuki-Trotter decomposition andonly optimize two sites simultaneously (simple update). In the new cluster optimization,
we consider the imaginary time evolution ofsix sites and optimize them simultaneously. In
Fig. 1, we plot thus obtained speciﬁc heat for
the Kitaev model in the thermodynamic limit.
Although both the simple update and the cluster optimization reproduce two expected peakstructures of the speciﬁc heat, the lower peakheight is largely increased by using the cluster

228

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

optimization. Because the speciﬁc heat calculated by unbiased QMC shows a higher peakthan the present data obtained by the cluster
optimization, the accuracy is still not enough.
However, we hope it will be improved when we
consider larger bond dimensions in the future.

0.2

Simple update
Cluster update

imental observation [3]. We also investigatedthe eﬀect of the weak Γ and Γ′ terms, and we
found that depending on the sign of Γ and Γ′ ,the thermal Hall conductivity is largely modiﬁed from that of the pure Kitaev model.
A part of this work has been done in collaboration with J. Nasu, T. Misawa and Y. Motome.

References

0.15[1] A. Kitaev: Ann. Phys. 321 (2006) 2.

0.1[2] G.
JackeliandG.
Khaliullin,
Phys. Rev. Lett. 102 (2009) 017205.

0.05[3] Y. Kasahara, T. Ohnishi, Y. Mizukami, etal., Nature 559 (2018) 227.

0[4] J. Nasu, M. Udagawa, and Y. Motome:
Phys. Rev. B 92 (2015) 115122.

0.001

0.01

0.1

1

10

100

[5] H. Li, D.-W. Qu, H.-K. Zhang, et al.
Phys. Rev. ResearchBw 2 (2020) 043015.
Figure 1: Speciﬁc heat of the inﬁnite Kitaev
model calculated by iTPO method. (open circle) Speciﬁc haet calculated by the standard
simple update. (ﬁlled circle) Speciﬁc heat calculated by the new cluster optimization.
Secondary, we investigated the thermal Hall
conductivity under a magnetic ﬁeld. As wediscussed, the accuracy of the inﬁnite size simulation is still insuﬃcient. Thus, we calculated
the thermal Hall conductivity for the ﬁnite size
cluster, and the density matrices are represented by matrix product operators (MPO) instead of iTPO. It has been demonstrated thatby using such MPO representation, indeed, one
can calculate ﬁnite temperature properties of
the Kitaev model very accurately [5]. By using
this MPO representation, we successfully calculated the thermal Hall conductivity undera magnetic ﬁeld. For the pure Kitaev model,
the thermal Hall conductivity has a peak at aﬁnite temperature, and it overshoots the halfquantized value; it is consistent with the exper-

229

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

Conformation of ultra-long-chain fatty acid in lipid
bilayer
Kazutomo KAWAGUCHIa , Koh M. NAKAGAWAb , Satoshi NAKAGAWAa ,Hideo SHINDOUc , Hidemi NAGAOa , and Hiroshi NOGUCHIba
Institute of Science and Engineering, Kanazawa University,
Kanazawa, Ishikawa 920-1192, Japanb
Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581
cDepartment of Lipid Signaling, National Center for Global Health and Medicine,
Shinjuku-ku, Tokyo 162-8655, Japan; Department of Lipid Science, Graduate Schoolof Medicine, University of Tokyo, Bunkyo-ku, Tokyo, 113-0033, Japan

In living cells, the most abundant lipids
are phospholipids, which have a polar headgroup and two hydrocarbon tails. Each tailtypically contains between 14 and 22 carbon
atoms. Around C22 fatty acids such as docosahexaenoic acid (C22:6, DHA) are calledvery long-chain fatty acids (VLCFAs). Moreover, much longer chains, called ultra-longchain fatty acids (ULCFAs), with 32 to 36
carbons with 6 double bonds were found at
the sn-1 position of phosphatidylcholine (PC)in photoreceptors, spermatocytes, fibroblasts,
and keratinocytes. However, the biosynthestic mechanisms and biological roles of ULCFAs
have not understood yet. In this study, we examined the conformation of one of ULCFAs,
dotriacontahexaenoic acid (C32:6) containing
phosphatidylcholine (dTSPC, C32:6-C18:0), in
phospholipid bilayers using all-atom molecular
dynamics simulations [1].
One dTSPC is embedded in the lipid bilayer
consisting of one of the following three lipids:
Distearoyl PC (DSPC, C18:0-C18:0), StearoylDHA PC (SDPC, C18:0-C22:6), and stearoyloleoyl PC (SOPC, C18:0-C18:1). Moreover, we
examined the effects of the difference in lipid
density between the two leaflets.
We found that the ultra-long tail of the ULCFA flips between two leaflets and fluctuates
among an elongation into the opposite leaflet,lying between two leaflets, and turning back.
The time scale of the conformational change
between the elongated and turned shapes is ∼230

10 ns, so that it is quite fast. The ratio of thesethree states depends on the lipid-density difference between the two leaflets. The sn-1 chain is
located at the opposite leaflet more frequently,
as the lipid density of the opposite leaflet relatively decreases. We have clarified the linear
relationships between the position of the sn-1
terminal of dTSPC and the lipid-density difference in all three types of membranes.
The main difference among the three types
of host lipids is the number of the double
bonds, and clear effects appear in the orderprofiles. Nevertheless, the conformation of dTSPC exhibits no qualitative differences. A minor influence is found in the distribution of C32
of dTSPC. Different shapes are obtained between SDPC and the others, whereas those of
C18 and C3 are not. In the SDPC membrane,
C32 has a rounded triangular distribution. In
contrast, a small peak or shoulder shape appears for the DSPC and SOPC membranes.
As described above, we revealed that ULCFA can sense and rapidly respond to thelipid-density difference to reduce such a difference. This behavior may be essential for thefunctions of ULCFAs in living cells.

References[1] K. Kawaguchi, K. M. Nakagawa, S. Nakagawa, H. Shindou, H. Nagao, and H.
Noguchi: J. Chem. Phys. 15, 165101(2020).

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

Giant magnetic response of hidden SU(2) symmetricantiferromagnets induced by impurity disorder
Hidemaro SUWA
Department of Physics, University of Tokyo
Hongo, Bunkyo, Tokyo 113-0033

In the research area of antiferromagnetic(AFM) spintronics, one of the most fundamental and challenging subjects is the efficient
manipulation of the AFM order parameter.
The key to extreme susceptibility lies in twodimensional (2D) systems with spin isotropy.
Although the ground state typically has a longrange order, any continuous symmetry is notbroken at finite temperature in a 2D systemwith short-range interactions, according to theMermin-Wagner theorem. As a result, themagnetic susceptibility exhibits an exponential divergence in the low-temperature region.
Here the strongest divergence is obtained in
the presence of the spin isotropy.
To achieve a giant magnetic response, wehave proposed exploiting the hidden SU(2)symmetry built in a spin-orbit system [1]. The
large spin-orbit coupling of iridates gives rise toJeff = 1/2 quantum states and pins the Jeff =1/2 pseudo-spins on the IrO6 octahedral network. The octahedral rotation around the axis
perpendicular to the plane can lead to a canted
in-plane magnetic order and a linear coupling
to uniform magnetic fields without breaking
the spin isotropy. This hidden SU(2) symmetry is revealed by a staggered rotation of thespin reference frame under the frustration-free
condition of the Dzyaloshinskii-Moriya interaction. Nevertheless, real materials have otherperturbations to the pure 2D spin isotropic system: in particular, spin anisotropy introducedby the mixing with the Jeff = 3/2 orbitals. For
efficient manipulation of the AFM order, it is

crucial to make the perturbations smaller thanthe external fields we control.
We have shown that the AFM tunability
can be systematically and significantly enhanced by introducing magnetic dilution anddiminishing the spin anisotropy effect in the
2D antiferromagnets. We used ISSP System B
and C in class C projects (ID: 2020-Ca-0097,
2020-Cb-0091) and calculated the large-scalequantum spin system consisting of morethan one million spins using ALPS/looper(http://github.com/wistaria/alps-looper).
Hybrid parallel computation was performed
using up to 3,456 cores. We numericallyshowed that the non-magnetic substitution
significantly suppresses the spin stiffness.
Therefore, the perturbation of the controllable
magnetic field, quantified by the ratio of
the Zeeman energy to the intralayer spin
stiffness, is significantly enhanced and dominant over the other perturbations. We haveaccomplished an extreme response of SU(2)symmetric AFM fluctuations in a diluted
system close to the percolation threshold. The
quasi-2D system of iridates thus represents a
prototype of the magnetic dilution problem in
2D pseudo-spin-half antiferromagnets, leading
to a tunability of effective spin anisotropy.

References[1] L. Hao, D. Meyers, H. Suwa, et al., Nat.
Phys. 14, 806 (2018)

231

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

Development of Thermal Functional Materials Using
Materials Informatics
Junichiro SHIOMI
Department of Mechanical Engineering,
The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo 113-8656
principles lattice dynamics, the Monte Carlo

The ability to effectively reduce a material’s
its

ray-tracing method, and effective medium

electrical properties unchanged is the ultimate

theory was used to understand the mechanism

goal for thermoelectric materials. Due to theofrelative difference in the mean free paths

nanostructured silicon as well as the weak

(MFPs) of electrons and phonons, introducing

temperaturenanostructures can effectively scatter phonon

conductivity values as well as its temperature

while affecting less to the electron transport.

dependence

For example, the thermal conductivities ofexperimental measurements.

thermalconductivity

nanostructured

whilematerials

likekeeping

phonontransportindependence.
are

also

The

verified

multiphasethermalby

thesuperlattice,nanowire, and nano-porous material are about
one or two orders of magnitude lower than thebulk counterpart.
On the other hand, Silicon also has many
allotropes other than the commonly seen
diamond phase (referred to as Si-I phase) andhas nature low thermal conductivity. Underhigh pressure, it can also form the bodycentered cubic (BC8, referred to as Si-III)phase and rhombohedral (R8, referred to as SiXII) phase as well as many other intermediate

Figure 1(a) Bulk thermal conductivity of Si-I, Si-III, andphases. How those multiphase nanostructure

Si-XII phases of silicon ure (b) The thermal conductivity

will affect the thermal conductivity remains

spectrum as a function of effective mean free (c) The effective

unknown.

thermal conductivity of nanograin Si-I, Si-III, and Si-XII at

In this work, we utilized a multiscalemodeling

techniquethat

includingdifferent temperature. (d) The effective thermal conductivity

first-of multiphase nanograined silicon as function of temperature.

232

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

The bulk phonon properties of the three

a temperature-independent behavior. This weakphases are shown in Fig. 1(a), which aretemperature dependence suggests that thecalculated from anharmonic lattice dynamics

phonon-boundarywith interatomic force constants from first-

temperature-independent process in the elastic

principle calculations. The bulk phonon thermal

regime, is the dominate process that limits the

conductivities of these three phases of silicon

thermal transport. To further verify theare quite high, which is unfavored for thesimulation, we also compared the simulation

thermoelectric energy conversion. Even for the

results with the experiment measurements, and

Si-III phase, which has the lowest thermal

the results are shown in Fig. 1(d). The inputs

conductivity among the three phases, theparameters like the volume fractions and thethermal conductivity is still above 10 W/mK at

averaged grain sizes are extracted from theroom temperature. Base on the Monte Carlosamples, and the good agreement between the

ray-tracingthat

simulations and the experimental measurements

introducing nanograin is an efficient way tosuggests the prediction power of our multiscale

reduce its lattice thermal conductivity. Taken

modeling scheme.

simulations,

we

found

scattering,

which

is

a

the Si-phase as an example [Fig. 1(b)], wefound that with the decreasing in the grain sizes,

References

both the effective mean free path and the[1] C. Shao, K. Matsuda, S. Ju, Y. Ikoma, M.

thermal conductivity will be shifted to low

Kohno, & J. Shiomi, Journal of Appliedvalues. The effective thermal conductivity the

Physics, 129(8), 085101 (2021)

nanograined

siliconasafunctionoftemperature is shown in Fig. 1(c), which shows

233Activity Report 2020 / 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, K. TAKABA, M.M. SOBEH, Y. MIYAZAWA, R. KIUCHI,
T.N. WIJAYA, T. OGAWA, Y. TAIRA, R. MISU, and A. KITAOSchool of Life Science and Technology, Tokyo Institute of Technology
Ookayama, Meguro, Tokyo 152-8550
As an efficient conformation sampling

sampled during previous cycles of eePaCS-

method, we previously developed parallel

MD and solving the “convex hull problem”

cascadeprovide the edges and vertices.

selection

molecular

dynamics(PaCS-MD) [1]. PaCS-MD consists ofcycles of multiple independent short MD
simulations. The initial structures of eachcycle were selected as the closest structures
to a target structure based on suitableproperties

such

as

root-mean-square-

deviation (RMSD), center-of-mass distance.
These selections increase the probability ofrare event occurrences, and PaCS-MD
enhances conformation sampling without
external perturbations.

Fig. 1 Schematic illustration of the evolution

One problem in the original PaCS-MDof the vertices and edges in eePaCS-MD.

is a requirement prior knowledge of thetarget structure. To address this issue, we
proposededge

expansion

The sampling efficiency of the method

PaCS-MD

can be examined by conducting eePaCS-MD(eePaCS-MD) [2]. As illustrated in Fig. 1,

from one state and examinimg whether theeePaCS-MD

initialsampled conformations contain the target

structure (black cross). The initial structures

structures or not. For this purpose, we appliedof each cycle are randomly taken from theeePaCS-MD to open-close transitions (Fig.2)

vertices (stars) of a multi-dimensionalof

principal

componentThemaltose/maltodextrin binding protein (MBP),

subspaceconsists

conformations

adenylate kinase (ADK). We conducted

requires

only

the

subspace.
of234

glutaminebinding

protein(QBP),

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

eePaCS-MD from both open and close states.

the close state when simulations started fromthe open state). For QBP and MBP, closeopen transitions were also successfully
sampled.

We

also

showedthatthecombination of eePaCS-MD and accelerated
MD(eePaCS-aMD) achieved

improvementofsampling

a furtherefficiency.

Minimum RMSD (RMSDmin), time to reach
RMSDmin (tmin), and the total computational
cost used to reach the opposite state for the
first time (T1st) of eePaCS-aMD are shown in
Table 1. OC and CO represent open to closeand close to open transition, respectively. As
shown in Table 1, open to close transitions
were observed within 10 ns, which is several
orders of magnitude shorter simulation timethan conventional MD. Thus, eePaCS-MD isa very efficient sampling method that does
not require prior knowledge.
Table 1 Summary of results (eePaCS-aMD)
TargetQBP (OC)QBP (CO)MBP (OC)MBP (CO)ADK (OC)ADK (CO)Fig. 2 Structures of target complexes. (A)Glutamine
Maltosebindingbinding

protein(QBP).

(B)

protein(MBP).

(C)

Adenylate kinase (ADK). Structures colored

RMSDmin (Å)
1.3
1.3
1.2
1.3
1.9
3.2

tmin (ns)
4.9
6.4
6.7
7.2
8.9
9.3

T1st (ns)
40.4
52.8
41.6
49.4
114
-

Referencesin black and white correspond to close and[1] R. Harada and A. Kitao. J Chem Phys 139,

open states, respectively.

035103 (2013).
[2] K. Takaba, D.P. Tran and A. Kitao. J

For all targets, eePaCS-MD efficiently

Chem Phys 152, 225101 (2020).

sampled open-close transitions (generated
conformations containing structures similar to235

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

Anderson-Kitaev spin liquid
Masahiko G. YAMADA and Satoshi FUJIMOTO
Department of Materials Engineering Science, Osaka University
Toyonaka, 560-8531, Japan

More than half a century, the Anderson
localization has been investegated intensively
on the low-dimensional electron systems or
fermionic systems. Normal magnets with along-range order does not have a fermionic excitation, and thus the randomness effect has
been investigated in a different manner. On
the other hand, in Kitaev spin liquids excitations are described by Majorana fermions,and the Anderson localization of Majorana
fermions can be expected. We name this new
localized state Anderson-Kitaev spin liquid.
This state potentially explains the experimental observation in α-RuCl3 or A3 LiIr2 O6 with
A = H, D, Ag.
We discovered Anderson-Kitaev spin liquidin the bond-disordererd Kitaev model witha numerical simulation. In the presense ofa magnetic field, we need to recompute theflux gap of the disordered Kitaev model inthe third-order perturbation. We compute
this quantity in two methods. One is the exact diagonalization of the Majorana one-body
Hamiltonian, and the other is the kernel polynomial method. Especially, the kernel polynomial method is very much suitable for theparallel computing, and we implement a hybrid parallelization code for this algorithm fora large-scale simulation. We note in the small
scale two methods agree well.

a rapid decay as the disorder strength grows.
After the extrapolation, the region where thethermal Hall conducitivity is quantized disappears when only 5% disorder is added. These
results are published in npj Quantum Materials [1].
Overall we discover the new Anderson transition into an Anderson-Kitaev spin liquid
thanks to the supercomputer in the Institutefor Solid State Physics. The calculation has
been done in the manner suitable for the supercomputer based on our newly establishedmethod.

References[1] M. G. Yamada: npj Quantum Mater. 5(2020) 82.

With a massive parallel computation, we
succeed in computing the flux gap, andhence the thermal Hall conductivity up to aO(10000)-site system for the bond-disordered
Kitaev model. In the 100 × 100 honeycomb
lattice, the thermal Hall conductivity shows236

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

Analyses on complex fluids using multiscale simulation
platform
Ｙohei MORII and Toshihiro KAWAKATSU
Department of Physics, Tohoku University, Sendai, Miyagi 980-8578
We have been developing a simulation platform

The obstacle is an elastic cylinder that rotates

named Multiscale Simulation Platform forwith a rotation rate (i.e. the speed at the surface

complex fluids (MSSP), with which one canof the obstacle/flow speed) 1.0, and thesimulate macroscopic complex flows coupled

Reynolds

with microscopic molecular simulators. MSSP

number 4.0, respectively.

uses smoothed particle hydrodynamics (SPH)

SPH particles, each of which contains 1000

method, where each hydrodynamic particle

FENE dumbbells. As a result, total number ofcontains microscopic molecular simulator or

dumbbells are 230,000,000.

microscopic stress calculator.

We can confirm that the rotation of thenumber

40.0,

andWeisenberg

We used 230,000

In this project, we implemented moving

obstacle induces an asymmetric flow patternboundary conditions into MSSP so that we canaround the obstacle, and that the dumbbells are

simulate moving and deforming objects such as

strongly stretched at the obstacle surface, which

floating filler particles and deformable vesicles

is the origin of the nonlinear behavior of the[1].

viscoelastic fluid.

Figure 1 shows several characteristic flow

We also develop a simulation technique withproperties of a viscoelastic fluid passing by awhich one can simulate an elastic membrane,

rotating cylindrical obstacle. The viscoelastic

i.e. a vesicle, in a channel flow.

fluid is modelled as a Newtonian fluid where

preliminary work for simulating red blood cells

many dumbbells, i.e. each consists of two

in a blood Bessel.

particles connected by a finite extensible

References

nonlinear elastic (FENE) spring, are immersed.

[1] Y.Morii and T.Kawakatsu, in preparation.

This is a

Figure 1: A FENE dumbbell flow around a rotating cylindrical obstacle. (a) Flow field,(b) pressure distribution, and (c) orientation distribution of FENE dumbbells.

237

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

Topological Order and Quantum Dynamics 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 have developed various novel numerical methods for quantum/classical many-bodysystems with strong correlations and performed large-scale simulations on ISSP supercomputer systems.
Rokhsar and Kivelson proposed the quantum dimer model in 1988 as a low-energy eﬀective model for frustrated magnets. Although
there is no negative sign problem in the Hamiltonian of the quantum dimer model, Monte
Carlo simulations are challenging because of
the strong geometrical restriction on the conﬁguration of the dimers. With the help of eﬃcient cluster update based on stochastic series
expansion, we now extend our simulation to amore general class of models, including dimers
and monomers and multiple interactions between them. It can include a larger Hilbert
space based on the quantum dimer model, under the same framework of Monte Carlo simulation and update techniques, by which we
can detect the deconﬁnement of monomers
under the ﬁnite temperature phase transition(Fig. 1). We also implemented the exchange
Monte Carlo method on the quantum dimer
model to make a winding number-free simulation, with the kinetic terms to break a dimer
into two monomers and still obtain the measurements for only dimer conﬁgurations.
In the meantime, the real-space renormalization group method using tensor networks
has recently been widely used as a numerical method for many-body spin systems. The
tensor network method can eﬃciently com-

Figure 1: Deconﬁnement of monomers in quantum dimer model.
pute the physical quantities of large-scale classical and quantum systems. However, existing methods, such as TRG and HOTRG, havea problem that the computational complexity becomes enormous as the spatial dimension increases. Our proposed Anisotropic Tensor Renormalization Group can dramaticallyreduce the computational complexity of highdimensional systems such as three-dimensional
quantum systems. Furthermore, by extending
the tensor network to include tensors not only
on the vertices but also on the lines connecting the vertices, the accuracy of the real-space
renormalization method can be increased bya factor of 100 compared to existing numerical methods that require a similar amount ofcomputation time.

238

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

Confirmation of stress-overshoot phenomena under biaxial
elongational flow of ring-linear mixtures
Katsumi HAGITA1 and Takahiro MURASHIMA2
1

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

Department of Physics, Tohoku University, Sendai 980-8578

Recently, we have studied (linear-rich) ring-

complexes such as catenanes and bonded-rings

linear mixtures [1,2] using Kremer-Grest (KG)

consist of two or three rings per complex as[3] bases coarse grained molecular dynamics

shown in Fig. 1. We studied the cases with

(CGMD) simulations. To study linear-chain(Nlinear, Nring) = (160,40), (160,80), (160,120),

penetrations into rings, we performed topology(160,160), (40,160), (10,160) and ring fraction

analyses based on the Gauss linking number

was fixed to be about 0.1. Here, Nlinear and Nring(GLN) for all pairs of linear and ring polymers.

denotes number of beads per a linear chain and

The probability distributions of the number ofa ring, respectively. To keep enough statisticallinear chain penetrations per ring (nP) wereprecisions, we used large system sizes withevaluated. We have investigated the possibility

approximate 0.7 M beads. As the initialof increase of nP by some conditions in order torelaxation run, we performed a tun with 109

make rings to work as a movable cross-linkers.

MD steps for each system. As the MD solvers,

In the preliminary simulations at the last

we used LAMMPS [5] for preparations at theyear, we discovered the increase of nP under

ISSP supercomputer. At later we also usedbiaxial elongational flows [4]. Here, we used

HooMD-bule [6] on the GPUs for continuous

our extended uniform extensional flow (UEF)

product runs.

method [5]. We also discovered the stressovershoot under the biaxial elongational flows.
Although the stress overshoot under uniaxial
elongational flows were reported, that under
biaxial elongational flows had been not. Thus,
we proposed the D-Class project to publish apaper as quickly as possible.
In the proposed D-Class project, we started
systematic preparation of ring-linear blends.

Fig. 1: Schematics of single ring, bonded-rings,

Here, we treated uncatenated rings and ring

poly-catenanes and ring-linear mixture.

239

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

In the present project, we developed amethod to estimate nP by using GLN. Here, theends of linear chains are virtually connected toeach other, but we prepared an extra linear
chain and connected it to the original linear
chain to form a cyclic chain as explained in ourwork [1,2]. For computation of GLN among
cyclic chain and ring polymer, we used theTopoly Python package [7]. Figure 2 shows theprobability

distributionsoflinearchainspenetrating into a single ring with (Nlinear, Nring)

Fig.

3:

Stressovershootunderbiaxial= (160,40), (160,80), (160,120), (160,160).

elongational flows. (Nlinear, Nring) = (160,40)

References[1] K. Hagita and T. Murashima: Polymer 218(2021) 123493.
[2] K. Hagita and T. Murashima: Polymer 223(2021) 123705.
[3] K. Kremer and G. S. Grest: J. Chem. Phys.
92 (1990) 5057.
[4]T.

Murashima,

K.

Hagita,andT.

Kawakatsu: Under review in Macromolecules.
[5]T.

Murashima,

K.

Hagita,andT.

Fig. 2: Probability distributions of linear chains

Kawakatsu: Nihon Reoroji Gakkaishi (J. Soc.

penetrating into a single ring with (Nlinear, Nring)

Rheol. Jpn.), 46 (2018) 207.

= (160,40), (160,80), (160,120), (160,160).

[6] S. Plimpton: J. Comp. Phys. 117 (1995) 117.
[7] J. A. Anderson, et. al.: Comput. Mat. Sci.
173 (2020) 173, 109363

In this project, we also performed deformed
simulations with the deformation rate of 0.001[8] P. Dąbrowski-Tumański, et. al.: Brief.

in order to grub preliminary results. As results,

Bionformatics (2020) bbaa196.

we obtained the key result for biaxial
elongational flows as shown in Fig. 3.

240

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

Frustration-induced symmetric skyrmion lattices
in three dimensions
K. MITSUMOTO1, R. OSAMURA2, K. AOYAMA2 and H. KAWAMURA1
1

Molecular Photoscience Research Center, Kobe University, Kobe 657-8501
2

Graduate School of Science, Osaka University, Toyonaka 560-0043

Considerable attention has recently been paid

dimensional (3D) with certain amount of

to “skyrmion”, a topologically stable nano-scaleinterlayer coupling. In this year’s project, we

spin texture, both from fundamental interest ofinvestigate the properties of the possible

topology and also from possible application toskyrmion and skyrmion lattice in 3D, by

spintronics. In the usual situation, skyrmion andperforming extensive Monte Carlo simulations“skyrmion lattice” are stabilized in the non-on the frustrated classical Heisenberg model on

centrosymmetric magnets without the inversion

the 3D stacked-triangular lattice. We investigate

center via the anti-symmetric Dzyaloshinskii-

by means of extensive Monte Carlo simulations

Moriya (DM) interaction, which energeticallyboth (i) the short-range interaction model

discriminates

spirals.

mimicking insulators, i.e., J1-J3-J1c Heisenberg

Meanwhile, it was proposed in Ref.[1] that themodel, and (ii) the long-range RKKY interaction

symmetric skyrmion and skyrmion lattice might(falling with 1/r3) model mimicking metals.

“right”

and“left”

be stabilized in centrosymmetric magnets withIn the short-range model (i), the transition

the inversion center by the frustrated exchange

behavior of the 3D model turns out to be similarinteraction. Examples might be the J1-J3 (J1-J2)

to that of the 2D model for the ferromagnetic

Heisenberg magnet on the two-dimensional

interlayer coupling, exhibiting the skyrmion(2D) triangular lattice under magnetic fields. In(anti-skyrmion)-lattice state, while, for thesuch frustrated magnets, in contrast to the DM

antiferromagnetic interlayer coupling, the naturecase, “right” and “left” spirals are energeticallyof the skyrmion-lattice state changes somewhat.

equivalent, and the skyrmion with both signs ofIn the RKKY model (ii), we also find thethe scalar chirality (or the topological charge) isskyrmion-lattice state, which, however, seems topossible, giving rise to both the skyrmion andaccompany an exotic replica-symmetry breaking

anti-skyrmion lattices. Interestingly, skyrmion(RSB) unlike the short-range model.

and anti-skyrmion lead to the topological Hall
effect of mutually opposite signs.
Realmagnets

areofcourse[1] T. Okubo, S. Chung and H. Kawamura, Phys.
Rev. Letters 108, 017206 (2012).

three-

241

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

Rotational dynamics of water molecules on the lipid
membranes
Yuji HIGUCHI
Institute for Solid State Physics, The University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

We calculate the rotational relaxation ofwater molecules on the lipid membranes
with phosphatidylcholine (PC) and phosphatidylethanolamine (PE) head groups. The
relaxation is faster in PE than that in PC(Fig. 1), which is consistent with the previous
experiment [4]. To reveal the mechanism of
the diﬀerent relaxation, hydrogen bond (HB)
networks among water molecules are analyzed.
The number of HBs in PE is less than that in
PC. The PE head group is strongly hydrated,which breaks the HB network. In contrast, theweaker hydration around the PC head groupdoes not break the HB network. In DFTB-MD
simulation, the head groups of PC and PE aremodeled as N(CH3)4+ and NH4+, respectively.
The rotational relaxation of water molecules
and the number of HBs around ions are consistent with those on the lipid membranes, conﬁrming the mechanism. Thus, we successfully
reveal the atomistic mechanism of the diﬀerent

water dynamics on the lipid membranes [5].

Rotational relaxation

We have studied the structures and properties of soft matter. We revealed the fracture processes of semicrystalline polymers [1]by coarse-grained molecular dynamics simulations. Understanding the mechanical properties is essential to improve the toughness ofpolymeric materials. This year, we study thewater dynamics on the lipid membranes, which
are considered as a model system for biological membranes. All atomic molecular dynamics (MD) simulations are performed on System
B and C using the LAMMPS program [2] andDFTB+ software [3].

1.0
0.8

PC
PE

0.6
0.4
0.2
0.0
0

10

20
30
Time (ps)40

Figure 1: The rotational relaxation of watermolecules just above the lipid membranes. Ten
trajectories are shown, which conﬁrms the signiﬁcant diﬀerence.

References[1] Y. Higuchi: Phys. Rev. E 103 (2021)
042502.
[2] S. Plimpton:(1995) 1.

J. Comput. Phys. 117[3] B. Hourahine, et al.: J. Chem. Phys. 152(2020) 124101.
[4] M. Hishida, K. Tanaka, Y. Yamamura,and K. Saito, J. Phys. Soc. Jpn. 83 (2014)
044801.
[5] Y. Higuchi, Y. Asano, T. Kuwahara,
and M. Hishida: Langmuir (2021). DOI:
10.1021/acs.langmuir.1c00417

242

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

FRG study on realization of the Kitaev quantum
spin liquid in condensed matter physics
Kiyu FUKUI, Yukitoshi MOTOME, and Yusuke KATO
Department of Physics, The University of Tokyo, Bunkyo, Tokyo 113-0033
Department of Applied Physics, The University of Tokyo, Bunkyo, Tokyo 113-8656
Department of Basic Science, The University of Tokyo, Meguro, Tokyo 153-8902

Frustrated quantum spin systems have beenstudied intensively for very long time. In
general, however, it is impossible to solve
them exactly mathematically. Therefore, numerical approaches to investigate them areindispensable. Numerical approaches usedwidely to study them are, for example, exact diagonalization (ED), quantum Monte
Carlo (QMC), density matrix renormalizationgroup (DMRG). Each method has pros andcons, of course. Since 2010s, a new method
called pseudo-fermion functional renormalization group (PFFRG or pf-FRG) has been usedgradually [1]. Although PFFRG overestimatesordering tendencies to magnetic orders [2],it can treat large-size systems and stronglyfrustrated systems and detect quantum paramagnetism due to strong frustration. We done
two studies on the feasibility of the Kitaev
quantum spin liquid by PFFRG [3].
The first one is the study on the feasibilityof it in ultracold molecular systems trappedin the optical lattice rather than in solids. We
define dipolar Kitaev model on the honeycomb
lattice as
1∑ 1
H=−
3
2
3rij
i̸=j

(1)

i,j

4π x x)]Si Sj
3
2π y y+ Jy [1 − 2 cos(2Φij −)]Si Sj
3+ Jz [1 − 2 cos(2Φij )]Siz Sjz },(2)× {Jx [1 − 2 cos(2Φij −based on the realization of Kitaev-type inter-243

actions proposed in 2013 [4]. In the dipolar
Kitaev model, the interactions between thenearest-neighboring sites are consistent with
the Kitaev model, but those between sites farther than the nearest-neighbor are more complicated. We apply PFFRG to this model toinvestigate the ground state at each parameter. The results show that FM order andzigzag AFM order are realized in the FM andAFM dipolar Kitaev model, respectively, forall anisotropy parameters. Furthermore, in order to investigate the connection between the(nearest-neighbor) Kitaev model and the dipolar Kitaev model, we introduce an artificial
range of interactions and calculate the susceptibility when approaching the dipolar Kitaev
model with long-range interactions from the
Kitaev model with only nearest-neighbor interactions. From this calculation, it is clarified
that the spin liquid state realized in the Kitaev model is quickly collapses, as the range
of interactions is extended. After the proposals of Kitaev-type interaction in ultracold polar molecular systems by microwave irradiationin 2013, the calculation based on these proposals has not been performed, and whether
Kitaev quantum spin liquid state is actuallyrealized has remained an open question. We
address this issue with PFFRG and elucidated
the above results for the first time.
The other is the study on the realization of
the Kitaev quantum spin liquid in the highspin materials. We calculate the phase dia-

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

gram of the spin-S Kitaev-Heisenberg model[4] S. R. Manmana et al., Phys. Rev. B 87,
with S = 1/2-5/2 and S = 50. Its Hamilto081106 (2013); A. V. Gorshkov, K. R.
nian isHazzard, and A. M. Rey, Mol. Phys. 111,
1908 (2013).
∑∑[
]
H=A
2 sin(2πξ)Siµ Sjµ + cos(2πξ)Si · Sj .
[5] M. L. Baez and J. Reuther, Phys. Rev. Bµ ⟨i,j⟩µ96, 045144 (2017).
(3)
We apply the extension of PFFRG to treat
spin-S systems proposed in 2017 [5]. The
obtained phase diagrams of the KitaevHeisenberg model for S = 1/2 and S = 1
are in general good agreement with the previous studies by other numerical methods. The
phase diagram for S = 50 is also in good agreement with the previous study on the classicalKitaev-Heisenberg model by Monte Carlo simulation, except for some special points. As aresult of systematic calculations with diﬀerent
S, for S ≤ 3/2, both the AFM and FM Kitaev spin liquid regions have a finite extent.
For S ≥ 2, no region showing Kitaev spin liquid state was found. Therefore, we believethat S = 3/2 gives an upper bound on thespins possessed by the candidate materials inwhich Kitaev quantum spin liquid is realized.
The phase diagram calculation of the KitaevHeisenberg model with a systematic change ofspin S, as we perform here, has not been done
before. This is the first study of the applicationof spin-S PFFRG to the Kitaev-Heisenberg
systems. The results we obtain here provide
a guideline for the recent intensive search forcandidate materials of S > 1/2 Kitaev quantum spin liquid.

References[1] J. Reuther and P. Wölfle, Phys. Rev. B
81, 144410 (2010).
[2] S. Göttel, S. Andergassen, C. Honerkamp,
D. Schuricht, and S. Wessel, Phys. Rev. B
85, 214406 (2012).
[3] K. Fukui, Ph.D Thesis, The University ofTokyo (2021).

244

Activity Report 2020 / 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 andbonds formed between C=O and N–H in theperform functions that are necessary to sustain

amyloid fibril are broken with each pulse oflife. However, when the concentration oflaser irradiation. In most cases, these bonds areproteins increases due to aging or other reasons,reformed spontaneously after the irradiation.

they can aggregate and cause a variety ofHowever, if a water molecule happens to enter

diseases. Alzheimer’s disease, which is one

the gap between C=O and N–H made by thetype of dementia, is caused by amyloid-β (Aβ)

laser irradiation, the reformation of thepeptides,sphericalhydrogen bonds is inhibited. This role of wateroligomers or amyloid fibrils. We are conducting

molecules is quite different from other known

theoretical studies on the aggregation andmechanisms. This new mechanism can explain

disaggregation of Aβ peptides [1]. In this fiscal

recent experiments showing that amyloid fibrils

year, we investigated the destruction process ofare not destroyed by laser irradiation under dry

an Aβ amyloid fibril by infrared free electron

conditions [2].

which

aggregate

intolaser irradiation using molecular dynamics

Furthermore, we found that more α-helix

simulation.

structures are formed after the laser-irradiatedIn recent years, Aβ aggregates have been

amyloid-fibril destruction. This is because theirradiated with infrared free electron lasers andα-helix structure has a different resonancedestroyed.

non-equilibrium

frequency from the β-sheet structure. Our

molecular dynamics simulations of Aβ amyloid

findings provide a theoretical basis for the

fibrils under a time-varying electric field

application of the laser to the future treatment

mimicking that of an infrared free electron laserof amyloidosis.

to clarify the destruction process and theReferences

structure of amyloid fibrils after destruction.

[1] H. Okumura, S. G. Itoh, J. Am. Chem. Soc.

We

performed

136 (2014) 10549-10552.

As a result, we have discovered a newmechanism by which water molecules destroy[2] T. Kawasaki, et al. Cellular and Molecular

the Aβ amyloid fibril. Intermolecular hydrogen

Neurobiology 38 (2018) 1039-1049.

245

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

Study on Complex Systems by Generalized-Ensemble Algorithms
Takuya HAYASHI and Yuko OKAMOTO
Department of Physics, Graduate School of Science, Nagoya University
Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602

We had two projects regarding generalized-

This estimate is in good agreement with theensemble algorithms applied to complex

results of several other research groups (seesystems.

Fig. 1).

The density of states (DOS) is one of themostimportantphysicalquantities

instatistical mechanics. Recently, we proposed asimulation protocol, REWL-MUCAREM [1],in order to obtain the DOS with high accuracy
in complex systems by combining the
Replica-Exchange

Wang-Landau(REWL)

Figure 1: The estimates of residual entropy ofmethod [2] and the Multicanonical Replica-

ice Ih by several research groups.

Exchange Method (MUCAREM) [3]. The
However, our previous results disagreed witheffectiveness of REWL-MUCAREM wasthe above results [5]. We found that thedemonstrated by using the 2-dimensional

discrepancy resulted from the choice of theIsing model in [1].

random number generators [4]. The present

We also applied the REWL-MUCAREM
protocol to the ice Ih system in order toestimate the residual entropy with highaccuracy. The residual entropy of ice has
become one of good examples to test theeffectiveness of sampling algorithms.
Our latest estimate of residual entropy
per one water molecule [4] is:

W0 = 1.507472 ± 0.000047
S 0 = kB lnW0

Figure 2: Generated random numbers by the= 0.815615 ± 0.000063 [cal/(mol K)].

Marsaglia generator (green) and Mersenne
Twister generator (purple).

246

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

work used Mersenne Twister generator [6],We see that the helix-coil transition

while our previous work used Marsaglia

temperature is higher for N=20 than for N=10.

generator [7]. In Fig. 2 we compare the twoItrandom number generators. Although thetemperature is higher for in vacuum than in

Mersenne Twister generator gives uniform

water.

is

By

distribution, Marsaglia takes on only nine

also

shown

calculating

that

the

this

transition

microcanonicalvalues in this narrow range. We haveentropy and its derivative with respect toreminded the reader that the choice of random

energy, we found that the latter has an

number generators is very important.

inflection point in vacuum but no inflection
point in water. This suggests that the helix-

In the second project, we studied thehelix-coiltransitionsofcoil transition is first-order-like in vacuum,

homo-alanine

but it is second-order-like in water [8].

polymers [8] based on the MicrocanonicalInflection-Point Analysis (MIPA) Method [9].
Two lengths of homo-polymers were taken,
N=10

and

20,

and

two

References

environment[1] T. Hayashi and Y. Okamoto, Phys. Rev. E.

conditions (in vacuum and in water) were100, 043304 (2019).

considered.

REWL-[2] T. Vogel, Y. W. Li, T. Wüst and D. P.

MUCAREM simulations. In Fig. 3, we

Landau, Phys. Rev. Lett. 110, 210603 (2013).

compare the specific heat.

[3] Y. Sugita and Y. Okamoto, Chem. Phys.

We

performed

theLett. 329, 261 (2000).
[4] T. Hayashi, C. Muguruma, and Y.
Okamoto, J. Chem. Phys. 154, 044503(2021).
[5] B. A. Berg, C. Muguruma, and Y.
Okamoto, Mol. Sim. 38, 856 (2012).
[6] M. Matsumoto and T. Nishiura, ACM
Trans. Modle. Comput. Simul. 8, 3 (1998).
[7] G. Marsaglia, A. Zaman, and W. W.
Tsang, Stat. Probab. Lett. 9, 35 (1990).
[8] T. Hayashi, Y. Sakae, and Y. Okamoto,submitted for publication.
Figure 3: Specific heat as functions of[9] K. Qi and M. Bachmann, Phys. Rev. Lett.

temperature for N=10 and 20.

120, 180601 (2018).

247

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

Numerical Diagonalization Study on the
Field-Induced Spin Nematic LiquidTôru SAKAI1,2 , Hiroki NAKANO1 , Rito FURUCHI1 , and Tokuro SHIMOKAWA3 ,1
Graduate School of Material Science, University of Hyogo,
Kouto, Kamigori, Hyogo 678-1297, Japan
2
Synchrotron Radiation Research Center, Kansai Photon Science Institute,
Quantum Beam Science Research Directorate,
National Institute for Quantum and Radiological Science and Technology (QST)
SPring-8, Kouto, Sayo, Hyogo 679-5148, Japan
3
Okinawa Institute of Science and Technology Graduate University, Onna,
Okinawa 904-0945, JapanThe quantum spin nematic state has attracted a lot of interest in the field of themagnetism. It is a kind of state between theconventional long-range antiferromagnetic order and the quantum spin liquid. This state
has been proposed as one of mechanisms of
the hidden order which was observed in somestrongly correlated electron systems. The spin
nematic state was theoretically predicted tobe realized in several frustrated systems; thesquare lattice model with the ferromagnetic
nearest-neighbor and antiferromagnetic nextnearest-neighbor interactions[1], the triangular lattice antiferromagnet with multi spin exchange interactions[2], and the ferromagnetic
and antiferromagnetic zigzag chain[3, 4]. The
previous density matrix renormalization group(DMRG) analysis[5] indicated that the fieldinduced spin nematic liquid phase appears in
the S = 1 antiferromagnetic chain with the biquadratic interaction. Thus it would be important to investigate this system using anothermethod, in order to confirm the spin nematic
liquid phase. We investigated the magnetization process of the S = 1 antiferromagnetic
chain with the biquadratic interaction, using
the numerical exact diagonalization of finitesize clusters[6]. The present analysis confirmed

248

the appearance of the field-induced spin nematic Tomonaga-Luttinger liquid phase andgives the ground state phase diagram whichis qualitatively consistence with the previous
DMRG study[5]References[1] N. Shannon, T. Momoi and P. Sindzingre,
Phys. Rev. Lett. 96 (2006) 027213.
[2] T. Momoi, P. Sindzingre and N. Shannon,
Phys. Rev. Lett. 97 (2006) 257204.
[3] A. V. Chubukov, Phys. Rev. B 44 (1991)
R4693.
[4] T. Hikihara, I. Kecke, T. Momoi and A.
Furusakai, Phys. Rev. B 78 (2008) 144404.
[5] T. S. R. Manmana, A. M. Lauchli, F. H.
L. Essler and F. Mila, Phys. Rev. B 83
(2011) 184433.
[6] T. Sakai, AIP Advances 11 (2021) 015306.

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

Symmetry Protected Topological Phase of the S=2
Antiferromagnetic Chain
Tôru SAKAI1,2 , Hiroki NAKANO1 , Rito FURUCHI1 , and Tokuro SHIMOKAWA3 ,1
Graduate School of Material Science, University of Hyogo,
Kouto, Kamigori, Hyogo 678-1297, Japan
2
Synchrotron Radiation Research Center, Kansai Photon Science Institute,
Quantum Beam Science Research Directorate,
National Institute for Quantum and Radiological Science and Technology (QST)
SPring-8, Kouto, Sayo, Hyogo 679-5148, Japan
3
Okinawa Institute of Science and Technology Graduate University, Onna,
Okinawa 904-0945, JapanThe spin gap of the integer-spin antiferromagnetic Heisenberg chain was proposed by
Haldane[1, 2] many years ago. Recently it has
been generalized as the symmetry protected
topological phase[3, 4], which is the gapped
phase of the odd-integer-spin antiferromagnetic chain. On the other hand, the numerical diagonalization and the level spetroscopy
analyses on the S = 2 antiferromagnetic chain
with the coupling anisotropy ∆ and the singleion anisotropy D, indicated that the symmetry
protected topological phase appears at a small
region in the ∆ − D phase diagram[5, 6]. This
phase corresponds to the intermediate D phasewhich had been predicted by Oshikawa[7]. In
the present study, the S = 2 antiferromagnetic
Heisenbergy chain with the biquadratic interaction JBQ and the single-ion anisotropy D isinvestigated using the numerical diagonalization and the level spectroscopy analyses. As
a result, it is found that the intermediate D
phase appears at a wide region in the JBQ − Dphase diagram.

References[1] F. D. M. Haldane, Phys. Lett. 93A (1983)
464.

249

[2] Phys. Rev. Lett. 50 (1983) 1153.
[3] F. Pollmann, A. M. Tuner, E. Berg and M.
Oshikawa, Phys. Rev. B 81 (2010) 064439.
[4] F. Pollmann, E. Berg, A. M. Turner and
M. Oshikawa, Phys. Rev. B 85 (2012)
075125.
[5] T. Tonegawa, K. Okamoto, H. Nakano, T.
Sakai, K. Nomura and M. Kaburagi, J.
Phys. Soc. Jpn. 80 (2011) 043001.
[6] K. Okamoto, T. Tonegawa and T. Sakai,
J. Phys. Soc. Jpn. 85 (2016) 063704.
[7] M. Oshikawa, J. Phys.: Condens. Matter
4 (1992) 7469.

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

Ground-State Phase Diagram of the S=1/2 Heisenberg-G
model on a Honeycomb Lattice
Takafumi SUZUKI
Graduate School of Engineering,
University of Hyogo, Shosha, Himeji, Hyogo 670-2280
from the isolated dimer limit, we found that a

In this study, we investigated the ground-

triplet dimer phase can survive up to thestate phase diagram of the S=1/2 Heisenberg-G

isotropically interacting system in a large partmodel on a honeycomb lattice [1] by theof the phase diagram, where the Heisenberg

numerical exact diagonalization method andcluster-series

expansion

method[2].

and G interactions are ferromagnetic andWe

antiferromagnetic, respectively. Otherwise, afocused on the effects of the anisotropic

phase transition to a magnetically orderedinteraction; by tuning the coupling constants,

phase takes place before the interaction

the system changes between the spin-chain andbecomes isotropic. This means that theisolated dimer models. We found that, in thequantum spin liquid discussed in the G modelspin-chain limit, there are three kinds of states,[3] is unstable against the anisotropy of thenamely a Tomonaga-Luttinger liquid and twointeractions. The obtained results, namely themagnetically long-range-ordered states. All

unitability of this spin liquid and the stable

three states change two-dimensional long-range

triplet dimer phase, are the same as those of theordered states, when the interchain interaction

previous results for the Kitaev-G model on ais included infinitesimally except for the casehoneycomb lattice [4].

where the Heisenberg interaction is much
weakerthan

theoff-diagonal

symmetric

Referencesinteraction (G). When the antiferromagnetic G[1] T. Suzuki, T. Yamada, and S. Suga, ininteraction is large enough and the system

preparation.

locates near the spin chain limit, we observe[2] J. Oitmaa, et al., Series Expansion Methods

that there are no prominent peaks of the staticforstructure factor and the feature of the low-

Cambridge University Press, 2006.

energy excitation is quite similar to that in the[3] A. Catuneanu, et al., npj. Quantum

spin chain limit. We consider that two-Materials, 3 23 (2018).

dimensional true-long-range-ordered state is[4] T. Yamada, T. Suzuki, and S. Suga, Phys.

suppressed in such parameter region. Starting

Rev. B 102, 024415 (2020).

250

Strongly

Interacting

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

Numerical analysis for nonequilibrium dynamics inelectronic systems on quasicrystals & Study of newordered phase and nonequilibrium phenomena inquasicrystals

4

(a)

(b)

Akihisa Koga∆
Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo
152-8551∆!"!
!"! "#!"

Recently, quasicrystals, which have no translational symmetry but have ordered lattice
structures, are attracting interests in the condensed matter physics. So far, their static
properties such as resistivity and specific heats
have been mainly focused on. Thus, it is important to reveal nonequilibrium features ofquasicrystal systems and understand how special lattice structures affect the dynamics.
Motivated by this, we studied nonequilibrium dynamics of the excitonic insulating (EI)
phase described by the two-band Hubbard
model on the Penrose lattice. The EI phase is
the condensation of the electron-hole pairs (excitons) in semiconductors or semimetals. The
EI phase is analogous to superconductors andis attracting much interests in the context of
nonequilbrium condensed matter physics. Bymeans of the real-space time-dependent meanfield theory for large systems, we revealed qualitatively different behaviors in the BCS and
BEC regimes, see Fig.1. Namely, in the BCS
regime, the order parameter is suppressed after
a single cycle pulse, while, in the BEC regime,
it is increased by the pulse. Such a behavior
has been reported in a previous work on thesquare lattice, and its origin was discussed in
the momentum space. However, in quasicrystals, one cannot define the momentum spaceand such argument does not apply. Our results indicate existence of a more fundamental
origin regardless of the momentum picture. In
addition, we revealed a characteristic dynamics in the BCS regime using the perpendicularspace analysis.

!
(a)

(b)
(c)

BCS

∆
∆ !"

!
(c)

4

BEC

!"∆
!

!"!

!"∆
! "#
!"

!"! "#!!

culatedin the area without edges [see appendix. A]. We define!
!NNB
those averaged MF parameters as |∆(τ)|NB = i=1|∆i (τ)|/NNB
!NNB
! "#
and 2nc (τ)NB = 2 i=1 nci (τ)/NNB , where NNB is the nubmerof sites away from the edge. We set NNB = 7936 in this study.
In Fig. 4, we plot |∆(τ)|NB and 2nc (τ)NB with dashed lines foreach F0 . Note that, for τ > 2π/ω, 2nc (τ) is conserved, but
2nc (τ)NB is not conserved. We find that the behavior of |∆(τ)|
and |∆(τ)|NB are qualitatively similar. The same can be said
about the relation between 2nc (τ) and 2nc (τ)NB . Therefore,
we can consider that the boundary eﬀect (finite size eﬀect) on
FIG. 4. (Color online) Time evolution of |∆| and 2nc for V = 1.95 andtheV ≃4.28 withistheconditionfor(N,qualitative
U, D, ω) = (11006,
4, 4, 0.4).
dynamics
negligible
discussions.
On(a)(b)theand 2nhand,≃ 4.28.
We plotthat|∆(τ)|
and 2nsize
with thedashedTime evolution of |∆| and 2nc for V = 1.95. (c)(d) Time evolution of |∆|otherc for Vwec (τ)NB
confirmed
theNBfinite
eﬀectsappears(|∆(0)|
line for each F0 . Horizontal black (dashed) lines represent |∆(0)| and 2ninc (0)the
NB and 2nc (0)NB ). in N = 11006 system used here
quantitative
aspects[see appendix B]. Therefore, it is hard to estimate the accuratequantitative values such as the frequency of the Higgs (ampli[see mode.[
Fig. 4(a)].
2nc (τ) slightly increases from 2nc (0) except
culated in the area without edges [see appendix. A]. We define tude)Θ͟Θ͟‫ݴ‬Θͳͯ͘΋͍͍‫͕͢ؾ‬Δɻͦ΋ͦ
!NNB−0.15. Note that the change in nc (τ) is] not
monoF0 =ϞʔυΒ͖͠΋ͷ͸‫͍ͯ͑ݟ‬Δͷ͔ʁthose averaged MF parameters as |∆(τ)|NB = i=1|∆i (τ)|/NNB ΋for
Higgs
Thus,
in
!NNB
tonic
againstthe fielddiscuss
strength
F0 . For example,
and 2nc (τ)NB = 2 i=1
nci (τ)/NNB , where NNB is the nubmer thispaper,
we mainly
qualitative
features. nc (τ) − nc (0)for F0 = −0.10 is smaller than that for F0 = −0.07 and nc (τ)
of sites away from the edge. We set NNB = 7936 in this study.
ForFV= 1.95 (the BCS region), |∆(τ)| is smaller than thefor
0 = −0.15 is smaller than the initial value. Therefore,
In Fig. 4, we plot |∆(τ)|NB and 2nc (τ)NB with dashed lines forequilibrium|∆(0)|isfor
the values
here
it is found valuethat there
thealloptimalvalueofofF0Fstudied
0 for increaseach F0 . Note that, for τ > 2π/ω, 2nc (τ) is conserved, but
ing 2nc (τ). For V ≃ 4.28 (the BEC region), both |∆(τ)| and2nc (τ)NB is not conserved. We find that the behavior of |∆(τ)|
2nc (τ) are increased from |∆(0)| and 2nc (0) for all the valand |∆(τ)|NB are qualitatively similar. The same can be saidues of F0 studied here [see Figs. 4(c) and (d)]. Again, one
about the relation between 2nc (τ) and 2nc (τ)NB . Therefore,
can see that the changes in nc (τ) and |∆(τ)| are not monotonic
we can consider that the boundary eﬀect (finite size eﬀect) on
against the field strength F0 . |∆(τ)|−|∆(0)| and nc (τ)−nc (0) for
the dynamics is negligible for qualitative discussions. On theF0 = −0.10 and −0.15 are smaller than those for F0 = −0.07.
other hand, we confirmed that the finite size eﬀects appears
Therefore, as in the case of V = 1.95, there is the optimalin the quantitative aspects in N = 11006 system used herevalueof F0 for increasing the MF parameters. It is found that[see appendix B]. Therefore, it is hard to estimate the accuratequantitative values such as the frequency of the Higgs (amplithe enhancements of |∆(τ)| and 2nc (τ) are prominent aroundtude) mode.[Θ͟Θ͟‫ݴ‬Θͳͯ͘΋͍͍‫͕͢ؾ‬Δɻͦ΋ͦ
F0 = −0.07.
΋ Higgs ϞʔυΒ͖͠΋ͷ͸‫͍ͯ͑ݟ‬Δͷ͔ʁ] Thus, in
These nonequilibrium features for the Penrose tiling arethis paper, we mainly discuss qualitative features.
qualitatively similar to those of the previous study on the
For V = 1.95 (the BCS region), |∆(τ)| is smaller than thesquare lattice19 . In particular, we observe the enhancement
equilibrium value |∆(0)| for all the values of F0 studied hereof the order in the BEC regime and the suppression of the or-

!"

!

!

phases in quasicrystals. We have clarified thenature of magnetic orders in several types ofquasi-periodic lattices [2].
Furthermore, we have also studied realspace dynamics and correlations in the Kitaev spin liquid states. We revealed peculiar
spin transport mediated by itinerant Majorana
fermions [3] and discuss how to modify the Majorana correlations using intrinsic degrees of
freedom of the Kitaev model [4].

References[1] K. Inayoshi, Master thesis (2021); K. Inayoshi, Y. Murakami, A. Koga, in prep.
[2] A. Koga, PRB 102, 115125 (2020):
Mater. Trans. 62, 360-366 (2021).
[3] T. Minakawa, Y. Murakami, A. Koga and
J. Nasu, PRL 125, 047204 (2020).
[4] A. Koga, Y. Murakami and J. Nasu, submitted to PRB.

251!

FIG. 4. (Color(d)online) Time evolution of |∆| and 2nc for V = 1.95 and V ≃ 4.28 with the condition (N, U, D, ω) = (11006, 4, 4, 0.4). (a)(b)Time evolution of |∆| and 2nc for V = 1.95. (c)(d) Time evolution of |∆| and 2nc for V ≃ 4.28. We plot |∆(τ)|NB and 2nc (τ)NB with the dashedline for each F0 . Horizontal black (dashed) lines represent |∆(0)| and 2nc (0) (|∆(0)|NB and 2nc (0)NB ).

Figure 1: Dynamics of the EI order parame∆!"ter∆ for the two-band Hubbard
model on the!"Penrose lattice. Horizontal lines indicate equilibrium values.

Another interesting topic is possible ordered!
(d)[see Fig. 4(a)]. 2nc (τ) slightly increases from 2nc (0) except
for F0 = −0.15. Note that the change in nc (τ) is not monotonic against the field strength F0 . For example, nc (τ) − nc (0)for F0 = −0.10 is smaller than that for F0 = −0.07 and nc (τ)
for F0 = −0.15 is smaller than the initial value. Therefore,
it is found that there is the optimal value of F0 for increasing 2nc (τ). For V ≃ 4.28 (the BEC region), both |∆(τ)| and2nc (τ) are increased from |∆(0)| and 2nc (0) for all the values of F0 studied here [see Figs. 4(c) and (d)]. Again, one
can see that the changes in nc (τ) and |∆(τ)| are not monotonic
against the field strength F0 . |∆(τ)|−|∆(0)| and nc (τ)−nc (0) for
F0 = −0.10 and −0.15 are smaller than those for F0 = −0.07.
Therefore, as in the case of V = 1.95, there is the optimalvalue of F0 for increasing the MF parameters. It is found thatthe enhancements of |∆(τ)| and 2nc (τ) are prominent aroundF0 = −0.07.
These nonequilibrium features for the Penrose tiling are
qualitatively similar to those of the previous study on thesquare lattice19 . In particular, we observe the enhancementof the order in the BEC regime and the suppression of the or-Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Rational computational design of proteins toward
drug discovery and industrial applications
Nao SATO1, Sairi MATSUMOTO1, Mizuki TERANISHI1, Katsuhisa IWAYA2,
Koji OOKA2, Masataka YOSHIMURA1, Shunji SUETAKA1, Yoshiki OKA1,
Yuuki HAYASHI1, and Munehito ARAI1,2
Department of Life Sciences, Graduate School of Arts and Sciences, and1

Department of Physics, Graduate School of Science,

2

The University of Tokyo, Komaba, Meguro, Tokyo 153-8902
Proteins have been extensively used intranscriptional

activator

mixed

lineageindustrial and pharmaceutical applications;leukemia protein (MLL) was used as a template.

many protein-based drugs are among the top

Theoretical design of the mutants of the MLL

selling drugs. However, these proteins have

TAD fragment that may bind KIX more tightly

been developed by experimental approaches

than the wild type was performed using Rosetta.

that require a lot of time and cost. Therefore,

Among the designed mutants, we could obtaintheoretical approaches to efficiently design

the mutant that binds KIX two-fold more

useful proteins are demanded. To solve thistightly than the wild type, suggesting that ourissue, we have been rationally designing

strategy is useful for designing PPI inhibitors.

variousproteins,

potentially

Severeacute

applicable to industry and medicine, using thecoronavirus

2

Rosetta 3.8 or 3.12 suite [1-3]. We have also

responsible for coronavirus disease 2019

performed experimental verification of the(COVID-19), infects human cells through thecomputational protein design.

PPI between the receptor-binding domain

Protein-protein

which

are

which

is

(RBD) of the viral-surface Spike protein andrelated to many diseases, and thus designing

the cell-surface receptor angiotensin-converting

PPI inhibitors is a promising way for drug

enzyme 2 (ACE2) [5]. To inhibit this PPI, we

discovery. We rationally designed inhibitors forhave computationally designed the antibodies

the interactions between the KIX domain of athat can bind the ACE2-binding region of thetranscriptionaland

RBD of SARS-CoV-2. We are currently

transcriptional activators, which are involved in

producing the antibodies and the RBD proteins,

many diseases including leukemia [4]. The

to experimentally verify the predictions.

coactivator

domain(PPIs)(SARS-CoV-2),

syndromearetransactivationinteractions

respiratory

CBP

(TAD)

of

the

Allergic asthma is known to occur through252

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

the PPI between interleukin-33 (IL-33) released

protein folding problem [9,10]. By extending

from damaged epithelial cells and the ST2

this model, we have previously developed a

receptor of type 2 innate lymphoid cells [6]. To

statistical mechanical model of protein folding

develop the IL-33–ST2 inhibitors, we rationallythat can explain the folding pathways of multi-

designed the IL-33 mutants that tightly bind to

domain proteins having multiple disulfide

one of the two IL-33 binding sites on ST2 but

bonds. To further extend this model for any

does not to another site. The interaction of ST2

types of proteins, here we applied this model towith one of the designed mutants of IL-33 wasa protein that does not have disulfide bonds andmeasured by fluorescence anisotropy. However,succeeded in predicting its folding mechanism.

the experiment did not support the theoreticalThus, our model may be useful for a unified

prediction,

theoretical description of protein folding

indicatingthenecessityofmechanisms.

improving protein design methods.
Alkane biosynthesis

hasgainedgreat

attention as an alternative to fossil fuels and is

References

expected as one of the promising ways to[1] B.J. Bender, et al.: Biochemistry, 55 (2016)

produce4748-4763.

carbon-neuralrenewableenergy.

Cyanobacterial alkane biosynthesis involves[2] R. F. Alford, et al.: J Chem Theory Comput.

two enzymes, an acyl-ACP reductase (AAR)

13 (2017) 3031-3048.

and an aldehyde deformylating oxygenase[3] J. K. Leman, et al.: Nat Methods, 17 (2020)(ADO) [7]. Interaction between these proteins665-680.

allows an efficient delivery of an aldehyde from

[4] P. E. Wright, et al.: Nat Rev Mol Cell Biol.

AAR to ADO [8]. To improve bioalkane

16 (2015) 18-29.

production, we rationally designed the mutants[5] J. Shang, et al.: Nature, 581 (2020) 221-224.

of ADO that have higher affinity with AAR.

[6] X. Liu, et al.: Proc Natl Acad Sci USA. 110

Our preliminary experimental results show that(2013) 14918-14923.

some of the designed mutants enhanced alkane[7] M. Arai, et al.: Adv Exp Med Biol. 1080

production as predicted.

(2018) 119-154.
[8] M. Chang, et al.: Biosci Biotech Biochem.

Finally, solving the protein folding problem
is important in improving rational design of

84 (2020) 228-237.

proteins. The WSME model can explain the[9] M. Sasai, et al.: Biophys Physicobiol. 13

folding pathways of small proteins and is(2016) 281-293.

considered to be a promising model to solve the[10] M. Arai: Biophys Rev. 10 (2018) 163-181.

253

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

Development of Giant Thermal-property Databases for
Amorphous Polymer Materials
Yuxuan Liao, Yoshihiro Hayashi
The Department of Mechanical Engineering
The University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo, JAPANWe intended to develop a giant thermal-

compute the long-range Coulombic interaction.

property database for amorphous polymer

In equilibration MD, NPT calculations usingmaterials. This year, we have established an

Nose-Hoover thermostat and the barostat wereautomation framework for developing the giant

run for 5 ns annealing and 8 ns equilibration.

databases,

TheAll bonds and angles, including those of theautomation framework allows us to generatehydrogen atoms, were constrained by theamorphous polymer structures and calculate its

SHAKE algorithm, and the time step was set tothermal conductivity (TC) automatically. The

1 fs. The TC is calculated by performing thesteps are: 1) creation of polymer chain byreverse NEMD (RNEMD) simulation proposed

random walk algorithm, 2) assignment of theby Müller-Plathe, as shown in Figure 2.

as

shown

in

Figure1.

force field and atomic charge into each atom, 3)

We have calculated 500 amorphous polymer

generation of the amorphous polymer cell, 4)systems, with an experimental validation shown

running

in Figure 2(B). The results are organized into atheequilibration

simulationbypaper, which will be coming online soon.

molecular dynamics (MD), 5) checking theconvergence of total energy and density, etc. forMoreover, we also calculate the thermal

the equilibration MD, 6) if not converged,transport properties of amorphous materials,

restarting the equilibration MD, 7) running thewhich is published in Nano Energy 2021 [2].

non-equilibrium MD (NEMD) to compute theReferences

TC value, 8) collection of TC values for all[1] Yoshihiro Hayashi, Ruimin Ma, Yuxuan

polymers.

Liao, et al, Automated calculation system of

All MD calculations were performed using

physical properties with molecular dynamics

the corrected heat flux branch of LAMMPS

simulation for polymer informatics (Preparing

package. We used an AMBER type GAFF2

for submitting)forcefield. Atomic charges were derived by the[2] Yuxuan Liao et al, Heat conduction below

Gasteiger

diffusive limit in amorphous superlatticemethod.

The

particle-particle

particle-mesh (pppm) method was employed tostructures,

254

Nano

Energy84,105903.

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

Figure 1 Automation framework for developing the giant thermal-property polymer
databases

Figure 2. (a) Schematic representation of the simulation box for reverse nonequilibrium
molecular dynamics (RNEMD). (b) Comparing the calculated thermal conductivity (κ)
and experimental κ for 13 polymers. A dots line indicates calculated κ equal toexperimental κ.

255

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

Development of Giant Thermal-property Databases for
Amorphous Polymer Materials
Yuxuan Liao
The Department of Mechanical Engineering
The University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo, JAPANin terms of the chain length.

We intended to develop a giant thermal-

Figure. 1 shows the thermal conductivity as

property database for amorphous polymer
materials. This year, we have 1) optimized theafunction

of

chain

lengths,

which

is

thermal conductivity calculation method and 2)

characterized by the number of atoms per chain.

investigated the influence of the classical size

The results show that 20000 atoms per chain

effect and the chain length of polymers on theare needed to obtain a converged thermal

thermal conductivity of amorphous polymers.

conductivity. The results indicate that heat

The motivation is that 1) Since thecarriers with long mean free paths (propagons)

calculation is computationally heavy, we needplay an important role in the heat conduction of

to find out the best NEMD method for ouramorphous polymers.

automation system. The Muller-Plathe (MP)To give a deeper insight into the underlying

method has the advantage of a rapidly

physics, we have explored the relaxation

converging quantity in temperature gradient

process and mean free path of propagons inrather than a slowly converging heat flux, thus,

amorphous materials. Here, we have used

it can speed up the calculation speed. The

amorphous silicon, amorphous silica, andresults from the MP method are validated byamorphous silicon nitride as examples. We

comparing them with the results from thefound that the Akhiezer mechanism dominates

Langevin thermostat, as shown in Fig. 1. Here,the relaxation process of propagons [1]. Thea polymer with the ID of P010080 is used as an

next step is how to generalize the analysis into

example. 2) Since the real polymer has

amorphous polymers, which is still under

extremely long chains consisting of over tens orinvestigation.

hundreds

Referencesofthousandatoms,

whichischallenging to be calculated with NEMD. The[1] Yuxuan Liao et al, Akhiezer Mechanism

thermal conductivity of polymers typically

Dominates

increases with the chain length. Thus, we need

Amorphous at Room Temperature, arXiv

to find out the converged thermal conductivity

preprint

256

RelaxationofPropagons

inarXiv:2006.07773.

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

257

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

Tensor renormalization-group study of spin glasses
Koji Hukushima
Department of Basic Science, University of Tokyo
3-8-1 Komaba, Meguro-ku, Tokyo 153-8902
One of the most important remaining problems in statistical physics is to clarify the
low-temperature thermodynamic properties ofIsing spin glasses in three dimensions. The
search for the ground states is computationallyintractable in the sense that it is one of the NPhard problems in the theory of computation,and the Markov-chain Monte Carlo (MCMC)
calculations at ﬁnite temperature have extremely slow relaxation phenomena, both ofwhich make large size calculations diﬃcult. Although many computational physical methods
derived from spin-glass research, such as the extended ensemble method, have aﬀected otherresearch ﬁelds, they have not yet solved thespin-glass problem.

We addressed two speciﬁc issues in Ising spinglasses.
(1) Vertically of ferromagnetic phase transition
and ferromagnetic limit on the Nishimori line
It is rigorously known that the ferromagnetic transition on the Nishimori line is theferromagnetic limit on disorder. MCMC calculations suggested that the phase boundary
is somewhat reentrant, but this is not deﬁnite.
The negative partition function problem is considered to appear at the point where the spinglassiness is strong. Therefore, it is expected tobe reduced by getting closer to the ferromagnetic side. Here, we determined the ferromagnetic transition temperature on the Nishimoriline with high accuracy by calculating the freeenergy diﬀerence under the twisted boundary
condition using TRG.
(2)Development of an algorithm to evaluate thespin-glass order parameter by TRG
Here, we developed a method to calculate thespin-glass order parameter by TRG, and performed the calculation for the Ising spin glassmodel. The spin-glass order parameter is an
overlap function, which is deﬁned using two independent systems with the same interactions.
There are two main ways to evaluate the expectation values of physical quantities using
the TRG method. One is to use the numerical derivative of the free energy, and the otheris to use the impurity tensor to construct themeasurement tensor of the observables. In any
case, it is necessary to devise a way to calculate
the spin-glass order parameter using TRG. We
constructed a tensor network for a replica coupled system that introduces an interaction between two systems, and succeeded in calculating the order parameter on the Nishimori line
as its derivative with respect to the coupling.

In this study, we aim to break the limitof computational size by applying the tensor
renormalization group (TRG) method, which
has attracted much attention in recent years,
to the spin glass system instead of the MCMC
method. We have previously applied TRG toa lattice glass model in two dimensions thatexhibits a glass transition in the mean-ﬁeld
limit[1]. There have been a few previous studies that have applied the TRG to spin-glass
systems, but they have not been studied extensively. One reason seems to be the problem of numerical accuracy. It has been pointed
out that, as a result of frustration, a “negative
partition function” appears even though it isoriginally positive deﬁnite, and there is anotherstudy showing that this can be solved by very
high precision ﬂoating-point arithmetic. However, it has been suggested that the numerical accuracy is not realistic for such calculations at low temperature. We proposed a numerical method to detect the ordered phase in
Ising spin glasses and explore the possibility of
performing the calculation while avoiding the [1] K. Yoshiyama and K. Hukushima:
negative partition function problem in TRG. J. Phys. Soc. Jpn. 89, (2020) 104003.

258

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

Parallel Bayesian computation in material science
Koji Hukushima
Department of Basic Science, University of Tokyo
3-8-1 Komaba, Meguro-ku, Tokyo 153-8902
Recent development of various experimental and measurement technologies has allowed
us to obtain a large amount of high-precision
data, but the methodologies for data analysis have not yet caught up with the improvement of experimental techniques. On the other
hand, the development of machine learning,such as deep learning, has also been remarkable in recent years, and the practical application of machine learning to materials sciencehas attracted much attention. However, theblack-boxed application of deep learning to materials science is powerful for “prediction” but
not necessarily suitable for “understanding”.
In this project, we attempted to apply anddevelop Bayesian inference based on statistical models that incorporate physical laws asmuch as possible, using data obtained from experiments as input. Since Bayesian statisticalmodels using physical processes are nonlinear
in general, they cannot be computed analytically. In the 1980s, the Markov-chain Monte
Carlo (MCMC) method was found to be usefulfor Bayesian inference, and it started to be applied in large scale. In recent years, the amountof data to be analyzed has increased and theamount of computation per step has becomelarge, and Markov chains which are bound by
causality, are not necessarily easy to parallelize.
Therefor, we applied population annealing(PA), which is a population type Monte Carlo
method, not based on Markov chains, tothe problem of Bayesian statistics of materials science. PA, proposed by HukushimaIba(2003), has been introduced as a computational method for statistical mechanical models, and recently it has been working on largescale parallel computations including GPGPU.
This method has a large number of nearly independent computations and is compatible withparallel computation, which signiﬁcantly re-duces the computational time by parallel computation. In addition, an advantage of PA over
MCMC methods is that the partition function
in statistical mechanics, which is the normalization constant for the probability distribution,
can be obtained as a by-product of the calculation. The partition function, which is called
evidence in the context of Bayesian inference, isan important quantity that gives an indicatorof model selection.
Furthermore, since the leave-one-out crossvalidation (LOOCV), which is one of the crossvalidation methods, can be formulated by thereweighting method often used in statisticalphysics, its implementation is also carried out
at the same time in PA. By systematically comparing the evidence and LOOCV, we developed
a methodology for data analysis that integrates
hyper-parameter estimation and model selection.
The proposed method has been implemented
in linear regression models with noise. It has
been pointed out that L1 regularization-based
methods such as LASSO lead to bias, and in order to avoid this, a direct method of analyzing
sparsity as an L0 regularization problem has
been considered. However, the L0 regularization problem is known to be NP hard, and noeﬃcient algorithm is known. We successfully
solved this problem in a realistic time using PA,
and also succeeded in evaluating the magnitudeof the noise by computing the evidence. Thismeans that we have established a new dataanalysis method for sparse modeling, including
uncertainty evaluation. We believe that it candevelop materials science through our sampling
methods such as noise magnitude estimationand model selection, which are a step further
than conventional Bayesian estimation usually
attributed to optimization calculations.

259

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

Study on relaxations of fluctuation with theevent-chain algorithm
Y. Ozeki and Y. Osada
Graduate School of Informatics and Engineering, The University of Electro-Communications

X
i

Si

!2 +1/2

T = 2.214
T = 2.212
T = 2.210

0.5

T = 2.208
T = 2.206
T = 2.204
T = 2.200
T = 2.198
T = 2.196
T = 2.194
T = 2.192
T = 2.190
T = 2.188
T = 2.186
T = 2.184
T = 2.182

0.4
0.3

0.2

0.1

1

10

t
Figure 1: Relaxation of |m(t) for the 3D XY modelin 2.182 ≤ T ≤ 2.222 with the interval ∆T = 0.002.
10(1)

and confirmed the efficiency of it.
Let us show the result of dynamical scaling analysis. The relaxation of the order parameter |m(t|)
are calculated for 2.182 ≤ T ≤ 2.222 on a 3003 simple cubic lattice up to an observation time of 100
MCSs. About 864 samples are taken for statistical
averaging. The estimations are plotted in Fig. 1.
Using the improved dynamical scaling analysis for
the NER data [4], we analyze the dynamical scalingform,|m(t, T )| = τ −λ Y [t/τ ],(2)
where τ (T ) is the relaxation time and λ is a dynamical critical exponent. The result is shown in
Fig. 2 with Tc = 2.202 and z = 1.47.
The exponents estimated by dynamical scaling
is used to deviate slightly because of corrections toscaling. To estimate the exponent more precisely,

260

l

N

−1

T = 2.220
T = 2.218
T = 2.216

Y = t m(t,J)|m(t)| =*T = 2.222|m(t)|

We investigate the applicability of the event
chain algorithm [1, 2] to the relaxation of fluctuation for the nonequilibrium relaxation (NER)
method, [3,4] which has been an efficient numericalmethod to evaluate critical exponents. The eventchain algorithm has introduced for the Monte Carlo(MC) simulation with a multi-spin-flip algorithm,
and shows estimation of the dynamical critical exponent z ∼ 1 [1], which indicates a faster dynamics
as compared with those with a single-spin-flip algorithm such as the the Metropolis one.
Previously, we investigated it for its applicability to the NER method and estimated the dynamical exponent z by means of the dynamical scaling
analysis. The initial state of relaxation is preparedas an all-aligned state. As an example, we analyze the Pclassical XY model in three dimensions;
H = −J hiji Si · Sj , where summation for hiji istaken over all nearest-neighboring sites on a simple
cubic lattice. Since the relaxation is much fast in
the simulation with the event chain algorithm, we
proposed the relaxation of the absolute value of themagnetization as the dynamical order parameter,

1

0.1
0.1

1

X = t/tFigure 2: Scaling plot for the data in2.182 ≤ T ≤2.222 in Fig. 1

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

2.60

fmm(T=2.202)zz(1/t0.704)

2.40
103

2.20

fmm

2.00

z

102

1.80
1.60

101

1.40
1.20

100

1.00
1

10

100

0

0.03 0.06 0.09 0.12 0.15 0.18 0.21

1/t0.704tFigure 3: Relaxation of fluctuation fmm calculated
at T = 2.202

Figure 4: Critical exponent z is estimated as t →∞.

we analyze the relaxation of fluctuation,h|m(t)|2 i
fmm (t) ≡ N−
1,
h|m(t)i2[2] M. Michel et. al., Europhys. Lett. 112
20003 (2015)(3)

at the transition temperature Tc = 2.202 estimatedabove. This function is expected to diverge as
fmm (t) ∼ tλmm ,(4)

where λmm = d/z. Calculations are carried out on
a 2003 simple cubic lattice up to the observation
time of 100 MCSs. About 645120 samples are taken
for statistical averaging. The result is plotted inFig. 3. We evaluate the numerical derivative
λmm (t) ≡d log fmm (t),d log t(5)and plot d/λmm (t) in fig. 4, where the horizontal axis is modulated as 1/t0.704 to show a precise
asymptotic value for t → ∞. The result is z = 1.67,
which is a little greater than the value obtained by
the dynamical scaling. Consequently, the eventchain algorithm can be used in the NER analysis
with an efficient relaxation performance, and would
be applicable to various slowly relaxing problems
such as frustrated and/or random systems.

References[1] Y. Nishikawa et. al., Phys. Rev. E92 063306(2015)

261[3] Y. Ozeki and N. Ito, J. Phys.: Math.
Theor. 40 R149 (2007).
[4] Y. Echinaka and Y. Ozeki, Phys. Rev. E94
043312 (2016).
[5] A. P. Aloysius and M. Harsenbusch,
Physica A201 593 (1993)

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

Numerical study of bulk-edge correspondence andtopological phases:
From quantum to classical mechanics
Y. Hatsugai
Department of Physics, University of Tsukuba
1-1-1 Tennodai, Tsukuba 305-8571, Ibaraki, JapanAs for standard topological phenomena, non trivial topology of the bulk is reflected in theboundary physics which can be, in principle, easily accessible by experiments. This is the bulkedge correspondence[1] that is widely used for various systems from quantum to classical ones.
In this project, we have tried to extend the universal feature of the bulk-edge correspondence byusing numerical methods. In this context, the adiabatic pump proposed by Thouless more than
30 years ago is quite special since the center of mass of the bulk, that has topological origin, isthe physical observable although the edge physics is not directly observed by a realistic finite
speed. Experimental realization the adiabatic limit including the gapless edge states is difficult
although it is quite useful for clear understanding of topological character of the pumping[2].
We have extended the bulk-edge correspondence of the topological pumping for various systems
with strong correlation[3, 4]. Another extension of the bulk-edge correspondence for the higher
order topological phases is the bulk-corner correspondence. We demonstrated its validity forclassical mechanical systems on a Kagome lattice[5]. We further applied numerical methods tovarious topological phenomena among quantum and classical ones. [6, 7, 8].

References[1] Y. Hatsugai, ”Chern Number and Edge States in the Integer Quantum Hall effect”, Phys. Rev. Lett.
71, 3697 (1993).
[2] Y. Hatsugai and T. Fukui, “ Bulk-edge correspondence in topological pumping”, Phys. Rev. B 94,
041102(R) (2016).
[3] Y. Kuno and Y. Hatsugai, “Interaction-induced topological charge pump”, Phys. Rev. Research 2,
042024 (2020)[4] K. Kudo and Y. Hatsugai, “Adiabatic heuristic principle on a torus and generalized Streda formula”,
Phys. Rev. B 102, 125108 (2020). Also K. Kudo, Y. Kuno and Y. Hatsugai in preparation.
[5] H. Wakao, T. Yoshida, H. Araki, T. Mizoguchi, and Y. Hatsugai, “Higher-order topological phases
in a spring-mass model on a breathing kagome lattice”, Phys. Rev. B 101, 094107 (2020).
[6] T. Mizoguchi, Y. Kuno, and Y. Hatsugai, “Detecting Bulk Topology of Quadrupolar Phase from
Quench Dynamics”, Phys. Rev. Lett. 126, 016802 (2021).
[7] T. Yoshida, K. Kudo, H. Katsura, and Y. Hatsugai, “Fate of fractional quantum Hall states in open
quantum systems: Characterization of correlated topological states for the full Liouvillian”, Phys.
Rev. Research 2, 033428 (2020).
[8] T. Yoshida, T. Mizoguchi, and Y. Hatsugai, “Mirror skin effect and its electric circuit simulation”,
Phys. Rev. Research 2, 022062(R) (2020).

262

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

Numerical study on low-energy statesof quantum spin systems
Hiroki NAKANO
Graduate School of Material Science, University of Hyogo

In condensed matter physics, we often haveto tackle many-body problems, in which it isdiﬃcult to estimate physical quantities precisely. Quantum spin systems are such typicalcases. To examine the systems, under circumstances, numerical approaches have widely andeﬀectively been employed in various studies. A
lot of computational investigations have been
carried out and gave us useful information of
the target systems.

we developed a hybrid-type parallelized codeof Lanczos diagonalization[1]. This Lanczosdiagonalization code enables us to treat various
large systems that have not been previously
treated yet within this method. We, thus,
study various quantum spin systems by thismethod as a primary approach in this project.

Within the field of quantum spin systems,
three methods are eﬀectively used. One is
the numerical diagonalization method. The
second is the quantum Monte Carlo (QMC)
method. The third one is the density matrix renormalization group (DMRG) method.
Each of the methods has advantages; at thesame time, it also disadvantages. In the QMC
simulations, large systems can be treated irrespective of the spatial dimensions of the systems although it is diﬃcult to precisely evaluate physical quantities in frustrated systems
due to the negative sign problem. When
spatial dimension of a target system is one,on the other hand, the DMRG method isvery useful irrespective of whether the target system includes frustration or not. For
the cases when the spatial dimension is larger
than one, however, this method is still under development. The numerical diagonalization method can be applied irrespective of thepresence of frustrations and the spatial dimension. However, this method has a serious weak
point that only very small system sizes canbe treated. To overcome this disadvantage,

In the project in 2020, we tackled theS = 1/2 Heisenberg antiferromagnet on theorthogonal-dimer lattice. This system was
examined based on the calculations of 36and 40-site systems of the same Lanczos
diagonalization[2], which pointed out that thethird characteristic ratio defined as J2 /J1 appears between the edge of the exact-dimer
phase and the edge of the Neel-ordered phase;
here, J1 is the amplitude of the exchange coupling on orthogonal-dimer bonds and J2 is theamplitude of the exchange coupling forming
the simple square lattice. In order to clarify
whether or not the behavior showing the presence of the third characteristic ratio of J2 /J1
survives when the system becomes larger, weadditionally calculated the 44-site cluster ofthis system, which has not been treated before.
As a result, our calculations for 44-site system show the behavior that is consistent withthose for 40-site system. Our examination contributes to our deeper understandings of thequantum antiferromagnets with frustrations.
Further investigations would clarify nontrivial
eﬀects of quantum nature and frustration.

263

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

References[1] H. Nakano and A. Terai: J. Phys. Soc.
Jpn. 78 (2009) 014003.
[2] H. Nakano and T. Sakai, J. Phys. Soc.
Jpn. 87 (2018) 123702.

264

Activity Report 2020 / 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 magnetism is wellknown as a source of exotic states of matter, but investigating the true nature, particularly at finite temperatures, is still a toughand challenging problem. In this fiscal year,
we investigated both the ground state andfinite-temperature properties of several twodimensional frustrated magnetism by means ofexact diagonalization (ED) and quantum typicality (QT) methods. We will here briefly report our numerical activities with the ISSP supercomputer.
1) S=1/2 breathing bilayer kagome(BBK) magnets for Ca10 Cr7 O28 magnet
Motivated by recent experimental results
on a new quantum spin liquid candidate,
Ca10 Cr7 O28 , we investigated the ground stateand finite-T properties of the S=1/2 BBK
model. We developed a new ED code specialized for high fields [1], which enables us totreat much larger system sizes at most 2000
spin clusters and compute spin dynamics without any bias in our quantum Hamiltonian. We
succeeded in reproducing bow-tie and ring-likefeatures in quantum spin dynamics, which inelastic neutron scattering (INS) measurements
reported as signatures of the quantum spin liquid state of this material. We also computed
the temperature dependences of specific heatand equal-time spin structure factors by meansof the QT method, and we found a realizationof a spiral spin liquid state at moderate temperatures having a ring-like structure in equaltime spin structure factors. [2]Recent experimental results under high pressures revealed the intermediate plaquettesinglet phase in SrCu2 (BO3 )2 compound. In
this intermediate phase, we naturally expect
to have a spontaneous mirror symmetry breaking at a finite temperature; however, detecting
the possible finite-T transition has been a challenging problem, such as in the traditional ED
method because of the smallness of accessible
system sizes. We used the quantum typicality(QT) method for this problem and succeededin capturing some clear signatures of the finiteT phase transition. For example, the temperature dependence of the real-space spin-spincorrelation is one of the valuable quantities to
see the signatures: we could see two possible
mirror symmetry broken correlation patterns
at low temperatures depending on the random
initial spin configurations used in QT method.
We also discussed the relationship to the material SrCu2 (BO3 )2 . These results are publishedin Phys. Rev. B. [3]References[1] H. Ueda, S. Yunoki, and T. Shimokawa,in preparation.
[2] T. Shimokawa, R. Pohle, and N. Shannon,
in preparation.
[3] T. Shimokawa: Phys. Rev. B 103 (2021)
134419.

2) Signatures of finite-T phase transition in S=1/2 Shastry-Sutherland model265

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

Global thermodynamic functions extended tononequilibrium steady states
Naoko NAKAGAWA
Department of Physics, Ibaraki University
2-1-1 Bunkyo, Mito, Ibaraki 310-8512

We study spatially inhomogeneous systems,
such as liquid-gas coexistence and multicomponent systems, from the viewpoint of global
thermodynamics [1]. We conducted numericalresearch on the following two topics by utilizing the supercomputer at ISSP. In both topics,
we performed molecular dynamics simulation
using LAMMPS.
(1) Numerical experimental protocol to determine the standard entropy and free energyof two-component fluids:
The “standard” entropy or free energy plays
a fundamental role to describe the thermodynamic properties of pure substances. This isalso the case in mixtures to distinguish theirvariety and industrial applicability. We havedeveloped a numerical protocol for determining the standard free energy of mixtures fromthermodynamic measurements by using knowledge of fluctuating thermodynamics and information thermodynamics. The protocol is valid
to numerical experiments, especially to molecular dynamics. In order to verify the validityof our protocol, we organized the molecular dynamics simulation taking a simple example onwhich we can make a theoretical prediction.
We concluded the validity and the eﬃciencyof our protocol as a method to determine thestandard free energy, mixing entropy, solubility, and so on.
(2) Osmotic pressure in a binary mixture
separated by a semipermeable membrane:
We have applied global thermodynamics to abinary mixture in thermal conduction, which

is separated by a semipermeable membrane.
It predicts how the osmotic pressure changes
from equilibrium. One of the point to be examined is whether the osmotic pressure is deviated from a conventional expectation, whichcorresponds to the balance of local chemical
potential. To clarify this point, we start largescale numerical experiments. The resources
used were mainly the L9 CPU of System C(because LAMMPS has not been installed in
the new System B). The simulation was run
on 9 nodes in 360 parallel computation using
MPI. This allowed us to perform long calculations on a system with up to 105 particles tomeasure the local and global thermodynamic
quantities in the steady state with a suﬃcientrelaxation time and a suﬃciently large numberof samples. We are now checking the eﬀect of
the finite system size.

References[1] N. Nakagawa and S.-i. Sasa: J. Stat. Phys.
177 (2019) 825.
(NAL )eq(NAR )eq

NAL =?

Tm

Tm

等温

NAR =?

TL

TR

熱伝導

Figure 1: Configuration of the investigatedsystem with two component fluid and a semipermeable membrane. The distribution of particles is modified by heat current.

266

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

Multiscale simulation of polymeric fluids and solids
Takahiro MURASHIMA
Department of Physics, Graduate School of Science, Tohoku University,
6-3, Aramaki-aza-Aoba, Aoba-ku, Sendai, Miyagi 980-8578
To develop multiscale simulation directly

problems. In this algorithm, the flow frame andconnecting macroscopic continuum system andthe simulation frame were separately treated. In

microscopic molecular particle system, we need

the similar way, we have developed a QR-

a suitable boundary condition for moleculardecomposition based deformation algorithm [2]dynamics

years,

and applied to a tensile deformation ofextensional flows were difficult problems inpolymeric solid. In this fiscal year, we extended

molecular dynamics simulation. Recently, we

this QR-decomposition based algorithm to be

have

[1]

applicable for rotational flow, rotational shear

applicable for the Langevin thermostat. Thisflow, pure shear flow, and so on. Fig.1 shows

new algorithm can solve the extensional flow

the rotational shear flow.

simulation.

developed

For

UEFEX

long

algorithm

Fig. 1: Snapshots of rotational shear flow. Non-symmetric deformation rate tensor applied.
Simulation box (the flow frame) is rotating clockwise. Black frame represents “the MD frame”.
The rotational shear flow is typical in cavity

References

flow and contraction and expansion flow. I[1] T. Murashima, K. Hagita, T. Kawakatsu, J.

would like to apply this technique to the ring-

Soc. Rheol. Jpn. 46 (2018) 207.

linear polymer blend [3], where the topological[2] T. Murashima, S. Urata, S. Li, Euro. Phys. J.

constraint of ring polymer causes uniqueB, 92 (2019) 211.

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

forexample,theviscosity

overshoot in a biaxial elongational flow.

arXiv:2102.02452.

267

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

Ground state and dynamical properties of the
J1J2K-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

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- 21 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],H = J1

X

S·S

+ K

(i,j,k,l)References

(1)

hi,ji

X 

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, we confirm the existenceof the bond-nematic phase in an extended
range of K/|J1 | while being sandwiched between the fully-polarized phase at high fields
and a long-range magnetically ordered phaseat low fields. A recently developed exact diagonalization method near saturation [6], showsa good agreement with iDMRG on cylindricalgeometries with circumference Lcirc ≥ 6 sites.
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 a gapless mode with vanishing spectral weight in the q → 0 & ω → 0 limit aspredicted by prior mean-field studies. We do,
however, observe some relevant qualitative differences, whose origin we are currently investigating.

−1
Pijkl + Pijkl
+ hz

X

Siz ,

i

where J1 represents dominant ferromagnetic[1] K. Harada and N. Kawashima, Phys. Rev.
B 65 (2002) 052403.
[2] A. Andreev and I. Grishchuk, Sov. Phys.
JETP 60 (1984) 267.

268

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[3] N. Shannon, T. Momoi, and P. Sindzingre, Phys. Rev. Lett. 96 (2006) 027213.
[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,in preparation.

269

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

Development of entanglement optimization method
Kenji HARADA
Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan

There are various interesting algorithms on
tensor networks to calculate the free energy
and critical properties of classical systems andthe ground state of quantum systems. In thesealgorithms, we often need to control the entanglement structure in a tensor network. For
example, we need to remove the short entanglement loop structure in grid tensor networks
to calculate a critical fixed point tensor.
In this project, we try to optimize an entanglement structure in a local tensor network
directly. In particular, we consider the development of a new optimization method fora branching operator[1]. To make a branching operator, we usually use the truncation
optimization for a local tensor network. Instead, we directly reduce the entanglement after applying a branching operator. Combining automatic differentiation and manifold optimization techniques, we can construct an iteration method to reduce the entanglement by
a branching operator.
Applying the new method to a tensor network representation of a Born machine[2, 3],
we can directly transform a MPS into a MERA
tensor network[4, 5]. In the case of the Bornmachine of the Ising model, we can constructa correct MERA tensor network from an exact
MPS of the Born machine. In two-dimensional
tensor networks, it is hard to transform the
PEPS into a tree tensor network. We needto transform the PEPS into the MERA tensor
network directly. It is also necessary to develop
a similar optimization method of a branching
operator for the two-dimensional network.

[2] Z.-Y. Han, J. Wang, H. Fan, L. Wang,and P. Zhang, “Unsupervised Generative
Modeling Using Matrix Product States,”
Phys. Rev. X, 8, 031012 (2018).
[3] S. Cheng, L. Wang, T. Xiang, and P.
Zhang, “Tree tensor networks for generative modeling,” Phys. Rev. B, 99, 155131(2019).
[4] G. Vidal, “Entanglement Renormalization,” Phys. Rev. Lett., 99, 220405(2007).
[5] G. Vidal, “Algorithms for entanglement
renormalization,”
arXiv:0707.1454v1(2007).

[1] K. Harada, “Entanglement branching operator,” Phys. Rev. B 97,045124 (2018).

270

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

In FY2020, we mainly treated the structuralvariation during the amorphous–amorphoustransition in SnI4 . We conducted in situ
synchrotron x-ray measurements along theisotherm at room temperature from 30 GPa
down to 1.1 GPa, using a diamond anvil cellto obtain the structure factor required for the
analysis. We also succeeded in applying anidea of nonequilibrium relaxation to machinelearning (ML) for predicting the order-disordertransition temperature in a 2-dimensional (2D)
Ising system. The last project was unscheduled at the outset. We exploited System B in
the following scene.
Amorphous–amorphous transition in SnI4 [1]We decided to analyze the structure employing reverse Monte Carlo (RMC) simulation because it is much advantageous than a usual
Fourier inversion in the following respect: we
can extract 3D structural information, fromwhich partial radial distribution functions are
obtainable. Before conducting the RMC ﬁtof the experimental structure factor S(k), we
had to determine the system’s density and initial conﬁguration transferred to the ﬁt as input. For this purpose, we prepared a system,
which “mimics” the actual system at the measured thermodynamic conditions, using classical isothermal and isobaric molecular dynamics simulation for a system composed of 2744
rigid regular tetrahedral molecules conﬁned ina cubic box with periodic boundary conditions.
System B took care of this process.
We call the system, whose density is to be
determined, a target system. Our model system [2] was known to capture the real systembelow 1 GPa. Hence, we ﬁrst searched in a271

trial-and-error manner thermodynamic conditions at which the model system “looks like”−1 ,the real system, at least in the length scale km
where km denotes the principal-peak positionof S(k). The conditions thus found were listed
in Table 1 of Ref. [1]. The resultant model
systems after equilibration were all metastableliquids.
Our method for density estimation from
S(k) presumes the existence of a reference
state that has a similar structure on the length−1 . Let n be a trial number density.
scale of kmt
We redeﬁne the system’s length by multiply−1/3
ing ntand calculate the structure factor,
St (k), using the scaled length. If we denotethe principal-peak position of St (k) by kt , wecan deﬁne the diﬀerence, δk = kt − km , asa function of nt . The target’s number density, n, is then estimated by the condition,
δk/δnt |n = 0. After ﬁnalizing the density, we
resized the whole system uniformly so that not
only the length of sides of the system but also
the interatomic distances were updated. We
transferred the conﬁguration along with thedensity to the RMC program.
The RMC simulation (see Fig. 1 for an
example) could reveal that the high-density
amorphous (HDA) state is divided into twostates: the high-pressure HDA state beyond 14
GPa containing isolated Sn atoms and the lowpressure HDA state below 14 GPa consistingof deformed molecules connected by metallic I2
bonds. In the latter state, the molecular shape
becomes C3v -like just before the transition to
the low-density amorphous (LDA) state, inwhich molecules recover the original Td symmetry. This local symmetry change has been

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyodetected on the liquid–liquid transition of SnI4 ,suggesting the strong coupling between the local symmetry and the global order parameterof density.

Figure 1: The RMC-simulation result obtained
from S(k) at 3.3 GPa, just before the transition to the LDA state on decompression, isvisualized as an example. SnI4 units (not necessarily molecules) are depicted in red, green,and purple when they are regularly tetrahedral
shaped, severely deformed, and even dissociated, respectively.
Applying nonequilibrium-relaxation scheme toML for detecting a phase transition [3]There seems no need to explain the signiﬁcance of applying the ML technique toproblems across the various disciplines today.
One example in physics problems is detecting the phase transition without any information about the “answers” a priori, as intelligibly illustrated for a 2D Ising model [4].
Tanaka and Tomiya identiﬁed that the weights,
more strictly, their sum Wsum , in a convolutional neural network (CNN) can play a relevant order parameter [4]. The identiﬁcation
then readily prompted us to apply the idea of
nonequilibrium relaxation [5] to the detection.
We employed the same CNN with one hidden layer besides the input and output layers.
To investigate the learning processes in detail,we realized the CNN in-house. We used thesame 2D Ising system on an L × L lattice withperiodic boundary conditions imposed for the272

examination, where L = 32, 64, 128, and 256.
We used as training data the systems well equilibrated at the speciﬁed temperature. Here,the temperatures constituted the labels, whichwere one-hot encoded into the binary vector.
Initially, we set Wsum to zero, which corresponds to the completely “unlearned” state.
Wsum soon evolved into one of the two stateswith learning, depending on the pattern’s temperature, leaving Wsum at the transition temperature Tc unlearned. We could thus identify Tc as the ﬁxed point in the label space.
Because Wsum with labels on both sides of
the ﬁxed point soon bifurcates with learning,a CNN with proper ﬁlters implemented canjudge the location of a Tc at an early stage oflearning. In this context, the existence of hidden layers is practically relevant to eﬃciency.
However, such a Tc irregularly depends on
L; Tc does not follow the relationship, Tc (L) ∼L−1 , expected from the ﬁnite-size scaling. The
unusual L-dependence happened to make suchan L = 32 system predict the transition temperature remarkably close to the exact one,and large systems used as training data do notnecessarily bring about a better result for predicting Tc .

References[1] K. Fuchizaki, A. Ohmura, H. Naruta, andT. Nishioka, submitted to J. Phys.: Condens. Matter.
[2] K. Fuchizaki and Y. Asano, J. Phys. Soc.
Jpn. 84, 064601 (2015).
[3] K. Fuchizaki, K. Nakamura, and D. Hiroi,
J. Phys. Soc. Jpn. 90 (2021), in press.
[4] A. Tanaka and A. Tomiya, J. Phys. Soc.
Jpn. 86, 063001 (2017).
[5] N. Ito, Physica A 192, 604 (1993).

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

Investigation of the molecular origins of the mechanical andthermal properties of realistic bio-polymers using allatomistic molecular dynamics
Okazaki Susumu
Department of Advanced Materials Science,
The University of Tokyo, Kashiwa-no-ha, Kashiwa, Chiba 277-8581
We have studied the impact fractures ofat different temperatures, within the observation

amorphous glassy polymer materials using all-time of current MD simulation time scale. Thisatomistic molecular dynamics (MD) simulations

causes the volume, the order parameter usuallyfor the past few years [1, 2]. In the second halfused to define glass transition, to show a broadof the past year, we continued the pursuit of theglass transition over decades of kelvins. The

molecular scale understanding of the yielding

same phenomenon is believed to occur inand the glass transition phenomena using theexperimental conditions, although the impact ofprotocols employing large-scale parallel MD

the faster degrees of freedom are less.

simulations and analyses established earlier, by

References

utilizing the supercomputer resources at ISSP.
In the first half, we investigated the stress[1] K. Fujimoto, Z. Tang, W. Shinoda, S.

origins of the yielding process of the brittle

Okazaki, All-atom molecular dynamics study ofmaterial poly-(methyl methacrylate) (PMMA)

impact fracture of glassy polymers. I: Molecularusing the stress decomposition method wemechanism of brittleness of pmma and ductility

developed previously, and found that 1) theof pc, Polymer 178(12) (2019) 121570.

brittle materials share the same stress origins as[2] Z. Tang, K. Fujimoto, S. Okazaki, All-atomthe ductile materials; 2) the dense and bulky

molecular dynamics study of impact fracture ofsidechains of PMMA makes it a brittle material

glassy polymers. II: Microscopic origins of

by lowering its flexibility and increasing thestresses in elasticity, yielding, and strainstress strength. [3]hardening, Polymer 207(20) (2020) 122908.

In the second half, we thoroughly studied the[3] Z. Tang, K. Fujimoto, S. Okazaki, A

glass transition phenomenon by analyzing thecomparison of the brittle PMMA with the ductile

MD trajectories, and realized that various

PC on the elasticity and yielding from adegrees of freedom, including the translation

molecular

and rotation of polymer chains, the dihedral(2021), in press.

angle rotation, and some angle bending, defreeze

273

dynamics

perspective,

Polymer

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

Global Optimization of Tensor Renormalization
Group using the Corner Transfer Matrix
Satoshi MORITA
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

To resolve this problem, we propose anotherapproximation for global optimization [5]. We
replace the environment tensor with the cornertransfer matrices (CTM) and the edge tensors,which can be updated by using the CTM renormalization group (CTMRG) [4]. Moreover, we
introduce additional decomposition, which reduces computational cost of tensor contraction for the coarse-grained tensor. The computational time of our algorithm (CTM-TRG)
in two dimensions scales as O(χ6 ) against
the bond dimension χ while the HOTRG and
HOSRG have O(χ7 ) scaling.
We perform benchmark calculations in the
Ising model on the square lattice and showthat the time-to-solution of the proposed algorithm is faster than that of other methods. As
shown in Fig. 1, our proposed method shows
better accuracy than HOTRG and compatible
to HOSRG.

Relative error in the free energy

Tensor network methods attract much attention as powerful tools for computing
strongly correlated many-body problems. The
tensor renormalization group method (TRG)[1] and the higher-order TRG (HOTRG) [2]provide a way to calculate coarse-grained tensor based on the singular value decomposition.
Both the methods do information compression
by solving local optimization problems. Approximations in these methods are locally optimal but not so for contraction of the wholetensor network. The second renormalizationgroup method (SRG) [3] and the higher-order
SRG method (HOSRG) [2] introduce global
optimization of tensor renormalization group.
Although these methods drastically improveaccuracy, calculation of the environment tensor requires performing the forward-backward
iterations.

10−7
10−8
10−9
10−10
10−11

HOTRG
χ=24
HOSRG
χ=24
HOTRG
χ=64

CTM-TRG
χ=24 FBC
χ=24 OBC
χ=64 FBC

10−12
2.24

2.25

2.26

2.27

2.28

2.29

2.3

Temperature

Figure 1: Relative errors in the free energy
measured in the Ising model on the square lattice.

References[1] M. Levin and C. P. Nave, Phys. Rev. Lett. 99,
120601 (2007).
[2] Z. Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P.
Yang, and T. Xiang, Phys. Rev. B 86, 045139(2012).
[3] Z. Y. Xie, H. C. Jiang, Q. N. Chen, Z. Y.
Weng, and T. Xiang, Phys. Rev. Lett. 103,
160601 (2009).
[4] T. Nishino and K. Okunishi, J. Phys. Soc. Jpn.
65, 891 (1996).
[5] S. Morita and N. Kawashima: Phys. Rev. B
103, 045131 (2021).

274

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

Dynamics of Phonon Entanglement Creation
between Remote Electron-phonon Systems
Kunio ISHIDA
School of Engineering and Center for Optics Research and Education
Utsunomiya University, Yoto, Utsunomiya, Tochigi 321-8585
Coherent control of quantum mechanical stateselectrons and photons. In this case, theis important to realize various states of quantum

wavefunction|Φ⟩is described bymaterials. In particular, quantum entanglement|Φ⟩ = ∑ √𝜆𝑛 |𝜂𝑛 ⟩|𝜃𝑛 ⟩,is a measure for quantum mechanical features,nand thus we have focused on its control method

where |𝜂𝑛 ⟩ and |𝜃𝑛 ⟩ denote the phonon statesby photoirradiation[1-3]. In the present study,

and the electron-photon states, respectively.

we study the dynamics of entanglement
generation

betweenremotesystems

We calculated the dynamical behavior of 𝜆𝑛byand the entanglement entropy between phonons

irradiation of a quantized light pulse byin |𝜂𝑛 ⟩ as functions of time, and found that theemploying a model of coupled electron-phonon-

entanglement in the singular vector |𝜂0 ⟩ for thephoton systems described by[2]largest singular value 𝜆0 is weak, while |𝜂1 ⟩ and|𝜂2 ⟩ for the second and the third largest singularℋ = ∑3𝑖=1 𝛺𝑖 𝑐𝑖 † 𝑐𝑖 + ∑2𝑗=1 [𝜔𝑎𝑗† 𝑎𝑗 + {𝜇(𝑎𝑗† +𝑗

𝜎 +1

values 𝜆1 and 𝜆2 are strongly entangled. The
phonon entanglement is generated as 𝜆1 and 𝜆2

𝑗𝑎𝑗 ) + 𝜀} 𝑧2 + {∑3𝑖=1 𝜈𝑖 (𝑐𝑖 † + 𝑐𝑖 ) + 𝜆}𝜎𝑥 ].
Solvingthetime-dependent

increase, i.e., the entanglement generation

Schrödinger

process corresponds to the conversion of theequation by numerical calculation on the Systemphonon states from a pure decomposable state to

B at ISSP, we found that the quantum mutuala mixed state of entangled states.

information for phonons reveals the dynamics ofReferences

phonon entanglement generation[4].
In order to investigate the dynamical behavior[1] K. C. Lee, et al., Science 334, 1253(2011).

of phonon states in both material systems, we[2] F. Borjans, et al., Nature 577, 195 (2020).

performed a Schmidt decomposition on the[3] K. Ishida, Eur. Phys. J. D73, 117 (2019).

electron-phonon-photon states in which the[4]K.

entire system is divided into the phonons and the2005.14615.

275

IshidaandH.

Matsueda,

arXiv:

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

A Theoretical Study for Thermal Unfolding of Proteins with
Quite Similar Native Structure
Takashi YOSHIDOME
Department of Applied Physics, Graduate School of Engineering, Tohoku University,
6-6-05, Aoba, Aramaki, Aoba-ku, Sendai 980-8579

from difference in the length of the secondary

Recently, drug design for inhibiting thestructures among the two proteins.

folding of proteins has begun. This indicatesthe importance of the elucidation of the

The calculations were performed using themolecular mechanism of folding process of

F4cpu and L4cpu in the system B, and L4cpu inproteins for drug design.

the system C. We used Gromacs program suite

In the present study, we focused on the two[2] for the MD simulations. CHARMM22 force

proteins that have quite similar native structures

field [3] and TIP3P model [4] were used for thebut different structures at the transition state [1].

proteins and ions and for the water molecule,

They are goat alpha-lactalbumin and humanrespectively.

lysozyme,

respectively.

To

elucidatethemolecular mechanism of the difference in the

Referencestransition state structures among goat alpha-[1] T. Oroguchi, M. Ikeguchi, M. Ota, K.

lactalbumin

andKuwajima, and A. Kidera: J. Mol. Biol. 371

performed

molecularhumanlysozyme,
dynamics

we(MD)(2007) 1354.

simulations for their unfolding processes at 400[2] H. J. C. Berendsen, D. Van der Spoel, R.

K.

van Drunen: Comput. Phys. Commu. 91 (1995)

We

successfully

reproduced

theexperimental results for the transition state

43 (1995).

structures.

[3] A. D. MacKerell, Jr. et al.: J. Phys. Chem. B,

By

carefully

analyzingthesimulation results, we elucidated that the

102 (1998) 3586.

difference in the transition state structures arose[4] W. L. Jorgensen, J. Chandrasekhar, and J. D.

from the existence of the Ca2+ ion in the native

Madura: J. Chem. Phys. 79 (1983) 926.

structure of the goat alpha-lactalbumin, and276

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

A Method for Analyzing Protein Dynamics: A Hybrid of
Cryo-Electron Microscopy Experiment and Molecular
Simulation
Takashi YOSHIDOME
Department of Applied Physics, Graduate School of Engineering, Tohoku University,
6-6-05, Aoba, Aramaki, Aoba-ku, Sendai 980-8579
cryo-EM

Cryo-electron (Cryo-EM) microscopy is one

microscopy

experiments

wereof the methods for analyzing the three-performed. In computing the maps, the protein

dimensional structures of biomolecules like

was rotated by assigning the polar coordinates

proteins. In this method, the three-dimensionalz=(θ, ϕ) as follows: the variable θ was fixed toelectronis0°; and the variable ϕ was randomly assigned inreconstructed using a lot of two-dimensional

the range of 0°≤ϕ≤360°. It was found thatelectron density maps of the protein projected

estimating the projection angle for eachfrom several projection angles. To reconstructprojection image was successful using thethe three-dimensional electron density map, theGTM method. In future, our method will be

projection angle for each map is necessary to

improved in order to apply to actual cryo-

estimate. A method is the “common-line

electron microscopy experimental data.

density

map

of

a

proteinmethod” [1], but verification method for the

The calculations were performed using theestimation using the common-line method has

F4cpu and L4cpu in the system B, and L4cpu inbeen lacking.

the system C. We used our custom-madeprograms for the calculations. The program of

In the present study, we investigatedwhether

generative

topographic

the GTM method was implemented using the

mapping

Python3 program language.

(GTM) method [2] can be used for an
estimating method. The method enables us toestimate the function describing the two-

References

dimensional electron density maps as the[1] S. J. Ludtke, P. R. Baldwin, and W. Chiu: J.

function of the projection angle. Thus, theStruct. Biol. 128 (1999) 82.

estimation of the projection angle can be[2] C. M. Bishop, M. Svensén and C. K. I.

possible with the GTM method. To demonstrate

Williams: Neural Comput. 10 (1998) 215.

the possibility of the estimation, simulations for277

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

Large-Scale Metadynamics Simulations for Water andAqueous Solutions
Hiroki NADA
National Institute of Advanced Industrial Science and Technology (AIST),
16-1 Onogawa, Tsukuba, Ibaraki 305-8569

Phase transitions in water and aqueousconsisting of 2880 water molecules using

solutions, such as crystallization and

ISSP system C with high efficiency.

amorphous formation, have been studied

Following a previous study [1], two

in fields of physics and chemistry fromdiscretelong ago. Nevertheless, there still exist a

distribution functions represented by

lot of unresolved issues concerning phase

Gaussian window functions were used astransitions of water and aqueous solutions.

collective variables.

For example, some of metastable phaseoxygen-oxygenradialA free energy landscape obtained bytransitions that are postulated to occur inthe

MTD-MD

simulation

for

the

water at a large supercooling, such as

TIP4P/Ice model at 233 K indicated two

liquid-liquid phase separation (LLPS),

local minima. One of the minima

have not yet been demonstrated.

corresponded to a low-density water

Molecular dynamics (MD) simulations

phase, and the other corresponded to ahave often been used to study thehigh-density water phase. The highstructure of water at a large supercooling.

density water phase had a structure

However, evidence for the occurrence of

resembling the structure of high-density

LLPS can hardly be obtained by standard

ice VII. At present, it is difficult to judge

MD simulations.

whether the observed two different water

In this study, a metadynamics (MTD)

phases corresponded to water phases

method was introduced to search forformed by the postulated LLPS ofmetastable phases of water at a large

supercooled water. More detailed studies

supercooling. An MD simulation inare needed to elucidate it in the future.

which

the

MTD

methodwas

implemented (MTD-MD simulation) was

References

performed for a liquid water phase[1] H. Nada: Sci. Rep. 10 (2020) 4708

278

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

Metadynamics Simulation Analysis of Various Cluster
Structures Appearing in Calcium Carbonate Supersaturated
Solution
Hiroki NADA
National Institute of Advanced Industrial Science and Technology (AIST),
16-1 Onogawa, Tsukuba, Ibaraki 305-8569

Terrestrial water on the Earth contains

structures that may stably appear in amineral components, such as calcium
carbonate(CaCO3),withasolution.

high

In this study, a metadynamics (MTD)

concentration and, hence, the formation

method was introduced to overcome theof CaCO3 crystals occurs ubiquitously.

timescale problem of MD simulation [1].

Industrially, the control of CaCO3 crystal

An MD simulation in which the MTD

nucleation is crucial in connection withmethod was implemented (MTD-MD

issues concerning the formation of scales,

simulation)

causing blockage of drainage pipes andsupersaturated CaCO3 aqueous solution.

reduction of thermal efficiency of boilers.

The simulation was performed efficiently

However, the control of it is quite

by parallel computing with multiple

difficult. This is because the mechanism

nodes of ISSP system B.

of CaCO3 crystal formation from a

A

was

free

energyperformed

landscape

for

a

(FEL)

supersaturated solution still has remained

obtained with the MTD-MD simulationunclear.

indicated

Atomistic

simulations,suchstableandmetastableasaggregates formed in the solution, each ofmolecular dynamics (MD) simulations,which correspond to a local minimum on

are helpful to elucidate the structure andthe FEL). The results are helpful tothermodynamic stability of precursors.

consider the mechanism of CaCO3 crystal

However, the timescale of phenomenon

nucleation via precursors.

that can be analyzed by an MD
simulation is normally on the order ofReferences

microsecond or shorter. This timescale is[1] A. Laio and M. Parrinello: Proc. Natl. Acad.

too short to pursue all possible precursor

Sci. USA 99 (2002) 12562.

279

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

Transport properties of multiple-q states infrustrated magnets
Kazushi AOYAMA
Department of Earth and Space Science, Graduate School of Science, Osaka University
Machikaneyama-cho, Toyonaka-shi, Osaka 560-0043

Transport phenomena in magnetic systems
reflect dynamical properties of interacting
spins, such as magnetic excitations and fluctuations. In this work, we theoretically investigate transport properties of two-dimensional
antiferromagnetic insulators, putting emphasis
on how nature of spin textures characterizedby more than one ordering wave vectors Q’s(multiple-Q states) is reflected in spin transport. In the J1 -J3 classical Heisenberg modelon the triangular lattice in a magnetic field,
the Hamiltonian is given by
X
1X
H=
Jij Si · Sj − H
Siz
(1)
2
i,j

i

with Jij = J1 δj∈N 1(i) +J3 δj∈N 3(i) , where N n(i)denotes the nth nearest neighbor sites of site
i. For ferromagnetic J1 and antiferromagnetic
J3 , it was shown by Okubo, Chung, and Kawamura [1] that single-Q, double-Q, and tripleQ states are stabilized in the low, middle, andhigh field regions, respectively. In particular,
the triple-Q state is known to be the skyrmion
lattice phase, and its higher-temperature region is another phase, the so-called Z phase, inwhich there are randomly distributed domainsof the skyrmion and antiskyrmion lattices.
The spin current and the associated spins are respectively given by
conductivity σµνJzs (t)
s
σµν

1X
= −
Jij rij (Si × Sj )z ,
2
i,j
Z ∞
1
z
z
=
dt hJs,ν(0) Js,µ(t)i, (2)
T L2 0where rij and L denote a vector connecting two

280

sites i and j and the linear system size, respectively. We numerically integrate the semiclassical equation of motion with initial equilibrium spin configurations generated by Monte
Carlo simulations and calculate the time correz (0) J z (t)i at each time step [2, 3].
lations hJs,ν
s,µIt is found that the longitudinal spin-currents is significantly enhanced at
conductivity σxx
the transition temperature between the tripleQ and Z phases, whereas at the highertemperature transition between the Z andparamagnetic phases, such an anomalous behavior is not obtained. This enhancement ofs is quite similar to the divergent behavior ofσxx
s at the Kosterlitz-Thouless-type topological
σxxtransition [2, 3] in which long-lifetime vortex
excitations play a crucial role. In the present
case, the topological object of the skyrmion
may be relevant to the spin transport, but further analysis is necessary to understand theorigin of the significant enhancement of thespin-current conductivity at the transition into
the skyrmion lattice phase.

References[1] T. Okubo, S. Chung, and H. Kawamura,Phys. Rev. Lett. 108, 017206 (2012).
[2] K. Aoyama and H. Kawamura, Phys. Rev.
B 100, 144416 (2019).
[3] K. Aoyama and H. Kawamura, Phys. Rev.
Lett. 124, 047202 (2020).

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

1/N bias of the hump of Binder parameter for Potts
Model
Hiroshi WATANABE
Department of Applied Physics and Physico-Informatics, Keio University
Yokohama, Kanagawa 223-8522, Japan

from the strong 1/N bias of the Binder parameter as it is a function of the expected value.
Therefore, this bias can be removed by the
jackknife resampling method. However, sincethe 1/N bias is removed and the N -infinity
limit still exhibit system-size-dependence, thefinite size scaling function appears to be system size dependent. The size-dependence ofthe peak position is stronger for larger valuesof Q, and the peak value is expected to be
infinite at N infinity for systems involving thefirst-order transition. The existence of the universal finite-size scaling function in the limit of
infinite system size for systems involving thecontinuous transitions is still non-trivial, andtherefore, the further study is required.
1.6

L= 64
L=128
L=256
L=512

1.58
1.56
1.54
1.52
Ulow

The Binder parameter has been widely used
to analyze critical phenomena [1]. Since theparameter is the dimensionless, the parameters for different system sizes are collapsed
into a single curve by scaling only horizontal
axis. The Binder parameter is usually a monotonic function, but it sometimes peculiar behavior. A typical example of this is the Potts
model. The Binder parameter of the Potts
model is non-monotonic function and exhibits
a hump near the criticality. Since the hump
has the system-size dependence, the Biner parameter of Potts model cannot be collapsed
into a single curve by scaling only horizontal
axis. The humps become larger as the system size increases and these humps interferewith the finite-size scaling analysis especiallyfor the high-temperature side of the criticalpoint. We also found that the peak value ofthe hump strongly depends on the MCs as
well as the system size. In order to identify
the origin of the hump, we study the Binder
parameter of Q = 3 Potts model. We foundthat the Binder parameter can be decomposed
into a high- and a low-temperature components, and that the humps originate from thelow-temperature one. THe system size andMCs dependence of the peak value of the lowtemperature component of the Binder parameter Ulow is shown in Fig. 1. One can see that
the MCs dependence becomes large as the system size increases. These MCs dependence arefound to be almost independent of the relaxation time. The MCs dependence originate

1.5
1.48
1.46
1.44
1.42
1.4

0

0.0005

0.001
1/MCs

0.00150.002

Figure 1: The system size and MCs dependence of the peak value.

References[1] K. Binder, Z. Phys. B Cond. Matt. , vol.
43, pp. 119-140 (1981).

281

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

Magnetism in the Multiple-Spin Exchange Modelon the Honeycomb Lattice
Chitoshi YASUDA
Department of Physics and Earth Sciences, Faculty of Science,
University of the Ryukyus, Okinawa 903-0213, Japan
3 He atoms adsorbed on graphite form theto a metastable state. When the external
magnetic ﬁeld is zero, the ground state in
the region where the contribution of the sixbody exchange interaction is small becomesan antiferromagnetic state. This is the sameground state as the antiferromagnetic Heisenberg model on the honeycomb lattice. On theother hand, in the region where the contribution of the six-body exchange interaction
is large, the eight-sublattice structure with anearest-neighbor correlation function and aﬁnite sublattice scalar chirality is stabilized
as the ground state. When an external magnetic ﬁeld is applied, the antiferromagnetic
state becomes the canted state at very smallﬁelds, while the eight-sublattice structure canbe maintained up to a certain ﬁnite ﬁeld.
When the applied magnetic ﬁeld is increased,
the eight-sublattice structure changes to thelarge sublattice structure, the four-sublattice
structure, and the canted state, sequentially.
In this project, we also research thermodynamic properties of the classical MSE modelusing the Monte Carlo simulation with theexchange Monte Carlo method. We conﬁrm
its usefulness in the classical MSE model on
the honeycomb lattice.

commensurate solid layer regarded as the
spin-1/2 quantum spin system on a triangular lattice at a certain density of 3 He.
There are not only an exchange of 3 He atoms
between two neighbor atoms but also thatamong three or more neighbor atoms. In relation to the solid 3 He, the multiple-spin exchange (MSE) model on the triangular lattice has been extensively studied. Groundstates and thermodynamic properties of themodel have been investigated using various
methods [1,2,3]. In particular, the model has
attracted attention owing to the existence ofnovel states such as the quantum spin liquid and the spin-nematic state, and owing to
the eﬀects of ﬂuctuations on ordered phaseswith chirality. Recently, a new quantum spin
liquid is experimentally observed in the solid
3 He layer. This layer is a monolayer 3 He adsorbed on graphite preplated with atomic layers of deuterium hydride. At low density of
3 He, it was suggested that there exists thequantum spin liquid with novel dependencesof the heat capacity and the magnetic susceptibility on temperature. The details of thisexperimental work are not yet clear, however,
it is suggested that the solid 3 He layer forms
a honeycomb lattice.

References[1] C. Yasuda, D. Kinouchi, K. Kubo: J.
Phys. Soc. Jpn. 75 (2006) 104705.

In this project, we research the classical
MSE model with the two-spin and six-bodyring exchange interactions on the honeycomb
lattice using numerical methods. The groundstate is investigated using the conjugate gradient method. In this method, we need toprepare many random spin states as the initial states because of avoidance of the trap[2] C. Yasuda, Y. Uchihira, S. Taira, K.
Kubo: J. Phys. Soc. Jpn. 87 (2018)
104704.
[3] S. Taira, C. Yasuda, T. Momoi, K.
Kubo: J. Phys. Soc. Jpn. 88 (2019)014701.

282

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

Microwave transmission through a Josephson
junction array
Tsuyoshi YAMAMOTO
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581

Recent technological development enables us
to perform quantum simulation using artiﬁcial
quantum systems . Its use for solving quantum many-body problems is one of hot topics in condensed matter physics. The superconducting circuits are one of promising platforms for such an attempt. However, background random charge in superconducting circuits hinders a “good” simulation. Our study
facilitates the opposite view. Namely, randomness in a circuit provides to simulate quantum
disordered systems [1].
In one-dimensional systems, localization occurs no matter how the disorder is weak.
While disorder in the one-dimension systemwith long-range interactions induces a transition into an insulating phase, its properties
is known little. In this study [2], we show
microwave transmission through a Josephson
junction array which represents one of the simplest one-dimension disordered systems with
the Coulomb interaction.
A one-dimensional nature of the system allows for the most comprehensive analytical andnumerical approaches. We ﬁnd signatures of
the insulating modes in the microwave transmission for a wide frequency range. At highfrequencies, the high accurate numerics reveals
that the localization property aﬀects the microwave transmission, which has been missedin the previous work [3], by the detailed comparison with the analytics (Fig. 1). At low
frequencies, by performing the parallel computation with respect to conﬁgurations to take a

Figure 1: The frequency dependence of thetransmission. The green line is the numerics
and the black dashed line is the analytics.
disorder average, we gather several non-trivial
pieces of information about statistics of thetransmission.
The versatility of the design of Josephsonjunction arrays motivates their further development for quantum simulation.

References[1] R. Kuzmin et al., Nat. Phys. 15, 930(2019).
[2] T. Yamamoto, L. I. Glazman, and M.
Houzet, Phys. Rev. B 103, 224211 (2021).
[3] M. Houzet and L. I. Glazman, Phys. Rev.
Lett. 122, 237701 (2019).

283

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

Molecular dynamics simulation of liquid BaTiO3
T. Hashimoto
National Institute of Advanced Industrial Science and Technology (AIST),
Tsukuba Central 2, 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, JapanThe structure of oxides in a liquid state
at high temperature is being studied by thecontainerless method using a floating furnace. The field strength is expressed as F =ZM /(rMO − rO )2 , where ZM , rMO , and rO are
the cation charge, the bond length of the cationand oxygen, and the ionic radius of oxygen. In
SiO2 , which has a large F , the coordination
number of oxygen around Si, whether solid or
liquid, is 4, but in ZrO2 , Y2 O3 , etc., where F
is small, the coordination number of oxygenaround the cation drops significantly in thecase of liquids[1]. Ti has an F value in the
middle of them, but in the case of BaTi2 O5 , it
has been reported that nTiO drops significantly
from the crystalline case to the liquid case[2].
Furthermore, it has been reported that thetemperature dependence of rTiO and nTiO is∂rTiO /∂T < 0 and ∂rTiO /∂T < 0. In the caseof BaTiO3 , it was reported by X-ray scatteringand neutron scattering by 46 Ti-rich and 48 Tirich samples that nTiO is reduced from 6 forcrystals to about 4.4 for liquids at 2073 K[4].
However, the temperature dependences of rTiO
and nTiO are not known. In addition, the validity of the partial radial distribution function
obtained by the empirical potential structurerefinement (EPSR) method using experimental
data has not been fully investigated. We calculated the structure of liquid BaTiO3 by firstprinciples molecular dynamics simulations.
The calculations were performed using
the plane wave pseudopotential method(QUANTUM-ESPRESSO), ultrasoft pseudopotential (GBRV), and GGA (PBE). The
density data from Paradis et al.[4] was ex-

trapolated to the high temperature side andlow temperature side. We calculated at 3500,
3000, 2500, 2000, 1500, 900 K using 135
atoms, using the (N, V, T ) ensemble and thevelocity scaling method.
The partial radial distribution functions of
the liquid BaTiO3 at 2000 K obtained bythis study[5] were compared with those by the
EPSR method at 2073 K[3]. The Ti-O results were in good agreement with each other,but Ba-Ba and Ti-Ti results were inconsistent.
The reason for this is that the pair weighting
factor contained in S(Q) is small for both neutron scattering and X-ray scattering, especiallyfor Ti-Ti, so there is a possibility that the interatomic interaction cannot be obtained correctly by the EPSR method. However, rTiO
and nTiO were in good agreement with the experiment. Liquid BaTiO3 had the same temperature dependence of rTiO and nTiO as liquid
BaTi2 O5 .
[1] L. Skinner et al., Phys. Rev. Lett. 112,157801 (2014).
[2] O. L. G. Alderman et al., J. Phys. Chem.
C 120, 26974 (2016).
[3] O. L. G. Alderman et al., J. Phys.: Condens. Matter 31, 20LT01 (2019).
[4] P.-F. Paradis et al., Appl. Phys. A 79,
1965 (2004).
[5] T. Hashimoto et al., J. Phys. Soc. Jpn 90,
044604 (2021).

284

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

Tensor-Network Renormalization Study ofFinite-Size Conformal Spectrum
Atsushi Ueda and Masaki Oshikawa
Institute for Solid State Physics, University of Tokyo
Kashiwa-no-ha, Kashiwa, Chiba 277-8581
April 30, 2021

We have studied the finite-size effects of classical 2D systems, using Tensor Network Renormalization(TNR). Based on the finite-size scaling theory of conformal field theory(CFT) [1],we demonstrated that it is possible to reconstruct an effective quantum Hamiltonian
from the spectrum of the renormalized tensors. Tracking the effective Hamiltonian ateach scale allows to numerically construct therenormalization group (RG) flow. We performed it for the classical Ising and XY model
as concrete examples, as shown in Figs. 1
and 2.
Figure 1: The numerically computed RG flowof the Ising model. t and h denote the deviations from the fixed point, corresponding to
the reduced temperature T − Tc and the magnetic field. The RG flows away from the fixed
point induced by the relevant perturbations t
and h are visualized.

285

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

of Physics C: Solid State Physics 7, 1046
(1974).
[3] Frank Pollmann, Subroto Mukerjee, Ari
M Turner, and Joel E Moore, Theory of finite-entanglement scaling at
one-dimensional quantum critical points,
Physical Review Letters, 102(25):255701,(2009).
[4] A. Ueda and M. Oshikawa, in preparation.

Figure 2: The numerically computed RG flowof the XY model. yK and yV are the perturbations corresponding to the spin-wave stiffness
and the vortex fugacity, respectively. Our numerical construction gives a quantitative verification of the celebrated Kosterlitz RG flow [2].
Throughout the study above, we also investigated the effect of the finite bond cutoff,which had been often overlooked previously.
We showed that there emerges an effective correlation length ξD = Dκ that is observed for
Matrix Product State(MPS) [3], by mapping
renormalized tensors to a reduced density matrix. As the exponent κ becomes smaller when
the central charge c is large, the focus of ourstudy (c = 2) turns out to call for a massive
calculation. Thus, we are in attempt to implement a parallel computation of TNR with thesupercomputer. More details of our study will
be reported elsewhere [4].

References[1] J. L. Cardy, Conformal invariance anduniversality in finite-size scaling, Journalof Physics A: Mathematical and General
17, L385 (1984).
[2] J. M. Kosterlitz, The critical propertiesof the two-dimensional xy model, Journal286

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

Magnetic field dependence of the thermal Hall eﬀect
based on the augmented quasiclassical equations
Takafumi KITA
Department of Physics, Hokkaido University
Sapporo, 060-0810

Charging of vortices in the equilibrium superconducting state was studied before calculating linear responses such as the thermal Halleﬀect. The results are reported below.
Vortices in type-II superconductors have not
only a single magnetic flux quantum but also
accumulated charge. It is known that thevortex-core charging is caused by three forces:(i) the Lorentz force [1], (ii) pair-potentialgradient (PPG) force [2], and (iii) the pressure diﬀerence arising from the slope in thedensity of states (SDOS) [3]. Recently, augmented quasiclassical (AQC) equations withthese forces were derived by incorporating thenext-to-leading-order contributions in the expansion of the Gor’kov equations in terms of
the quasiclassical parameter δ ≡ 1/kF ξ0 [4],where kF is the Fermi wavenumber and ξ0 is
the coherence length. Numerous studies on
the charging in an isolated vortex were performed [4], but the magnetic field dependenceof the vortex-charging has not been fully calculated. To this end, we calculated the vortexcore charging in two-dimensional s-wave superconductors with the Abrikosov vortex lattice
due to the Lorentz and PPG forces by using the
AQC equations. The SDOS pressure is now
absent for this case with the cylindrical Fermi
surface. The fixed parameters are the magnetic penetration depth as λ0 = 5ξ0 , Thomas–
Fermi screening length as λTF = 0.03ξ0 , andquasiclassical parameter as δ = 0.03.
Figure 1 plots the spatial variation of thecharge density due to the Lorentz and PPG

forces at temperature T = 0.2Tc and the average flux densities B̄ = 0.15Bc2 , 0.42Bc2 , and
0.88Bc2 , respectively, where Tc is the transition temperature and Bc2 is the upper critical field. Figure 2 shows the spatial variationof the charge density due to the PPG force attemperatures T = 0.2Tc and 0.5Tc , respectively.
The charge caused by the Lorentz force has astrong field dependence with a peak and canbe enhanced substantially from the value of an
isolated vortex as shown in our previous work[5]. Moreover, we show in this report that the
charge caused by the PPG force monotonicallydecreases as the magnetic field increases.
These calculations require a fine mesh ofspace and momenta when diﬀerentiating theGreen’s functions and self-energies numerically. Therefore, we performed parallel calculations using the ISSP supercomputer.

References[1] F. London, Superfluids (Dover, New York,1961), Vol. 1, p. 56.
[2] N. B. Kopnin, J. Low Temp. Phys. 97(1994) 157.
[3] D. I. Khomskii and A. Freimuth, Phys.
Rev. Lett. 75 (1995) 1384.
[4] H. Ueki, M. Ohuchi, and T. Kita, J. Phys.
Soc. Jpn. 87 (2018) 044704.
[5] W. Kohno, H. Ueki, and T. Kita, J. Phys.
Soc. Jpn. 85 (2016) 083705.

287

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(a) Lorentz, B̄ = 0.15Bc2(b) Lorentz, B̄ = 0.42Bc2(c) Lorentz, B̄ = 0.88Bc2(d) PPG, B̄ = 0.15Bc2(e) PPG, B̄ = 0.42Bc2(f) PPG, B̄ = 0.88Bc2

Figure 1: Spatial dependence of the charge density ρ(r) due to the Lorentz force ((a), (b), and(c)), and the PPG force ((d), (e), and (f)) at temperature T = 0.2Tc in units of ρ0 ≡ ∆0 ϵ0 d/|e|ξ02
on a square grid with x and y ranging from [−2ξ0 , +2ξ0 ] for the average flux densities B̄ =0.15Bc2 , 0.42Bc2 , and 0.88Bc2 from left to right, respectively. ∆0 is the energy gap at zerotemperature and zero fields, ϵ0 is the the vacuum permittivity, d is the thickness, and e < 0 is
the electron charge.
(a) T = 0.2Tc(b) T = 0.5Tc

70

2

Total

Total
1.5

50

ρ(0)/ρ0ρ(0)/ρ060

Lorentz
40

Lorentz
1

30
20

PPG

0.5

PPG

10
00

0.2

0.4

0.6

0.8

00

1

B̄/Bc2

0.2

0.4

0.6

0.8

1

B̄/Bc2

Figure 2: Magnetic field dependence of the charge density at the vortex center ρ(0) due to theLorentz force (green circular points), the PPG force (blue square points), and the total force(red triangular points), in units of ρ0 ≡ ∆0 ϵ0 d/|e|ξ02 as a function of the magnetic field, attemperatures (a) T = 0.2Tc and (b) 0.5Tc from left to right.

288

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

Calculation of ordered structures, dynamics and opticalproperties of soft materials
Jun-ichi FUKUDA
Department of Physics, Kyushu University, Motooka 744, Nishi-ku, Fukuoka 819-0395
We studied the structural properties offree energy with the variation of the direction

cholesteric blue phases, three dimensionalof the alignment imposed by the surface. We

ordered structures exhibited by a chiral liquidconsidered two cases commonly observed

crystal. We particularly focused on how theexperimentally: BPI with its (110) planes

cubic lattice of blue phases (BPI with O8

parallel to the surface(s) and BPII with its (100)symmetry and BPII with O2 symmetry) is

planes parallel to the surface(s). In both cases,

oriented when in contact with confining

the blue phase lattice adopts one specific latticesurface(s) that impose unidirectional surface

orientation that minimizes the free energy. The

alignment of the liquid crystal (unidirectionalsurface anchoring strength of 10-5 Jm-2, easily

surface anchoring). Our study was motivatedachievable experimentally, was shown to be

by a recent experimental study [1] thatsufficient to lock the lattice orientation. Ourdemonstrated specific lattice orientation of blue

finding is consistent with some of thephases in contact with such surfaces.

experimental observations, but unfortunately

Our study [2] was based on a continuum

not all. Our study provides a starting point fortheory describing the orientational order of athe understanding of how the blue phase latticeliquid crystal by a second-rank tensor (Landau-

is oriented in response to the surface anchoring,de Gennes theory). The free energy of thewhich is strongly desired for practical opticalliquid crystal was given as a functional of theapplication of blue phases that requires thetensor order parameter, and the profile ofpreparation of large-scale monodomain lattices.

orientational order was obtained by minimizing
the free energy numerically. The discretization

Referencesof the liquid crystal system by a regular grid[1] M. Takahashi et al.: J. Phys. D: Appl. Phys.

allows an efficient use of OpenMP (Our system

51, 104003 (2018).

was not large enough to require the use of MPI).

[2] J. Fukuda and S. Zumer: Phys. Rev. Res. 2,

We carried out systematic evaluation of the033407 (2020).

289

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

Molecular simulation of colloidal particles
Takamichi TERAO
Department of Electrical, Electronic and Computer Engineering,
Gifu University, Yanagido 1-1, Gifu, 501-1132
Structural formation of colloidal particles

in aggregates were also quantitatively clarified.

have been studied both theoretically andFrom these studies, we have been able toexperimentally for many years. In previous

clarify the behavior of patchy particle systems

studies, spherical colloidal particles with

with respect to the formation of various

isotropic

structures.

interparticle

interaction

wereinvestigated. Recently, there has been a great
deal of research on colloidal systems withanisotropic interactions. A typical example isthe so-called patchy particle system, in whichanisotropic interparticle interaction exists.
The aggregation phenomena of colloidal
particles have been a long-standing interest [1].
In

this

study,

we

performed

computerFig. 1: Snapshot of patchy particles.

simulations of the aggregation phenomenon inpatchy particles and clarified their structureformation. In Fig. 1, an example of thesnapshot of patchy particles with four patches
is shown. In particular, we have quantitatively
clarified the sol-gel transition and the phase
diagram of the system (Figure 2). In this figure,the abscissa and the ordinate denote therescaled dimensionless temperature and theFig. 2: Phase diagram of patchy particles.

volume fraction of the system, respectively. In
this figure, the red squares show that the patchy

References

particles form a gel network with these sets

[1] N. E. Valadez-Pérez, Y. Liu, and R.

volume fraction and temperatures. The cluster

Castañeda-Priego: Phys. Rev. Lett. 120 (2018)

size distribution and the static structure factor

248004.

290

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo                                                           












    ! 



Ð





    Ö         #



Ð




 !  !"   











! 
 !    $   %       
                &  '          &  ' (    
                   
( "  #     
            (             ) $   *            +        ),                "         +    + # +-) $          ./ /% &.//'        &  +'        &  +') .  .//    *                
     & ' $              )    & '   
   & '     

.//     
& '  & '    .//  ) .   .//          .//        .//  )0        
      .//   
Ö



¼

¼

¼



 

¼



··

··

  *1            
         $
       2    )3      &0'        &4'         
     
 &5'     6 2        
   ).         &4'      
        ) 7            
 &'         " &"  +' # &++  #  +'     +  +  &'         "&+8  "  +9'   &++    +:'   +  # +8  &'        # &++  #  +'   &++    +'   +  "  + ).          &' &'  &'    ;) ;)   ;) : ) $            & '  <=
&<'  &>? '   &;'       .//  ) $         + "  +   ++  #  +  .//        
           & ;) :');                ) 0                   .//          .//  <     
4*26*.    ) 0       6*    :             " ++) .      
    ) 7     &  '   &  '
                

291

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

nTLL

ΓF
0.5

F

0.3

H

XY1

JAF

1.0

0.2
sN

−0.5

0.0

0.0

0.5 ∆ 1.0
AF

−0.5

sN

0.0

∆AF 0.5

0.20
JAF

ΓF

0.9
H

XY1

H

XY1

0.15

0.8
F

0.7

nTLL

F

0.1
0.0
−1.0

H

XY1

−0.5

sN

nTLL

0.0

0.10

∆AF 0.5

sN

nTLL

F

−0.4

−0.2

0.0

0.2 ∆

AF

; @ A       + ; @ A       +  + ) 0      1    # +8) 0      1  
      )     )1.0
JAF                 * & +     ') 7        1*             *  ).       1     & '    ;     & +'     ) $          1*           ).                        1*  ) ;          <      B $ ) .     2        
        
* &C A'   D ) .     1*    
  C A %  " &'& ! ' " & ! '
 " &'  & +'   & +') 3  .//     
       ) .     *     1*   *    +           + 
      ) .    1* 
     " &'" &'
 " &'  &  '   & +') /       1*   

XY1
0.5
nTLL
0.0
0.0

F
0.5

ΓF 1.0; :@ A       + "  + )         ;       .//        & '  & +')$   6) E2 6) <  .) 32   2 ).) .  .) 2  6) E2 3) 5);  .) 32 F) C ) 3) F)  +D++
&+ G') H) /) 4*2 3) C ) FI.C  9 +& J8 '( F) K) 6*  B) F) .  F)
C ) 5  G & J8:'): 0) 6*  6) L  6) < F) C )
0@ K ) A)  /:- &++:')D K) C) <  C  0  -9 & J89')

292

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

Coarsening mechanism in mass-conserved
reaction-diﬀusion systems
Michio Tateno and Shuji Ishihara
Graduate School of Arts and Sciences, The University of Tokyo
Universal Biology Institute, The University of Tokyo
Komaba 3-8-1, Meguro-ku, Tokyo 153-8902

Reaction-diﬀusion systems (RDs) are central
mathematical scheme that comprehensively
describe molecular assemblies driven by chemical reactions, and have tremendously contributed to our understanding of pattern formation observed in nature, particularly thosein biological systems. Nearly 14 years ago, it
was reported that in RDs satisfying the massconservation law (MCRDs), patterns exhibit
coarsening and result in formation of single isolated domain, unlike usual behavior of RDs [1].
MCRDs were originally introduced and havebeen discussed as models for molecular localization such as membrane-bounded GTPases
which are responsible for the formation of cell
polarity [1]. Recently the relevance of MCRDs
has been recognized in a broad range of biological phenomena such as oscillatory motion
in min-Protein and chemical turbulence (see,e.g,[2]. Despite its uniquity and fundamental importance, the physical mechanism underlying coarsening behavior in MCRDs remains
elusive. This is partially because most of earlier studies only investigated one-dimensional
systems. Considering the fact that morphology of the patterns is strongly inﬂuenced byspatial dimensions, an intensive study in highdimension is highly desirable.
In this project, we realize a large-scale
MCRD simulation in both two and three dimensions using multiple GPUs, and successfully capture coarsening dynamics of MCRDs
in high dimension with statistical signiﬁcanceFigure 1: A spatial pattern obtained from athree-dimensional simulation. A contour surface indicates the interface between bistablestates. The color is labeled according to
the value of chemical-potential-like quantity in
MCRDs. The numerical simulation was performed in a cubic lattice of 10243 .

for the ﬁrst time (see Fig. 1). We havean eye on similarity between coarsening dynamics in MCRDs and classical phase separation systems, and apply analysis methods
established in thermodynamic phase transition (dynamic scaling law, nucleation-growth
process, etc.)
to the simulation results.
We reveal that droplets forming in MCRDs
obey the Young-Laplace law and coarsen following the evaporation-condensation (LifshitzSlyosov-Wagner) mechanism (see Fig. 2).
These outcomes indicate that in the presenceof conserved variable, a physical quantity similar to surface tension is relevant to RDs,

293

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

Figure 2: a Droplet-radius distribution ρ(R)for various time t. b. Dynamic scaling of
ρ based on Lifshitz–Slyozov–Wagner theory.
The symbol shows the same simulation results
as in panel a. The red dashed lines represent
theoretical predictions for the scaled distribution N (R/Rc ) [4], where Rc is critical radius.
which provides new insight into molecular selfassembly driven by chemical reactions. The
result is now under review in an international
peer reviewed journal, and the preprint is uploaded in an online paper repository [3].

References[1] M. Otsuji, S. Ishihara, et al.: PLOS Comput. Biol. 3, (2007) 1.; S. Ishihara, M. Otsuji and A. Mochizuki: Phys. Rev. E 75,(2007) 015203.
[2] E. Frey, J. Halatek, et al.: Physics of Biological Membranes, Springer International
Publishing, pp.229-260 ; J. Halatek and E.
Frey: Nat. Phys. 14, (2018) 507.
[3] http://arxiv.org/abs/2010.03900[4] J. H. Yao, K. R. Elder, et al.: Phys. Rev.
B 47, (1993) 14110.

294

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

Multicellular simulation by phase ﬁeld method
Michio Tateno and Shuji Ishihara
Graduate School of Arts and Sciences, The University of Tokyo
Universal Biology Institute, The University of Tokyo
Komaba 3-8-1, Meguro-ku, Tokyo 153-8902

The central aim of this project is to establish a multicellular simulation method thatcan capture multiscale dynamics, ranging fromintracellular chemical concentration, shapechanges of cells, to their collective motions in
a tissue scale. To realize this, we adopted thephase ﬁeld (PF) model, which can express cellswith arbitrary shape. In this method, an auxiliary ﬁeld (phase ﬁeld) which obeys a bistableenergy functional is introduced to distinguishthe inside and outside of a cell. This procedure allows us to describe the time evolutionof cell shape via interfacial dynamics without
suﬀering from moving boundary problems.
One of the authors (SI) and coworkers combined PF model with a reaction-diﬀusion system, in which a bistable and excitable chemical network is incorporated to model cell migration driven by F-actin. They demonstrated
that this model can precisely reproduce complex motion of amoebic cells (Dictyostelium atearly starvation stage) as well as intracellularchemical waves related to actin polymerization[1]. In this term, we have made an attempt toextend this PF model to multicellular systems.
In PF model for a single cell, a square moving frame with lattice size L2 is provided tocompute the phase ﬁeld. Here L should be
suﬃciently larger than the cell size (typically L
being hundreds of order) to capture the complex morphology of a cell and avoid the cellfrom overhanging the frame. When increase
the cells, of course, the computational cost is
not simply scaled by the number of cells N ;one need to deal with overlap regions among

Figure 1: a. Adaptive moving frame (gray areas) is used in our PF model. This largely reduces the computational cost for intercellularinteractions (the region surrounded by a dotted line). The white and blue regions represent
the distributions of PIP3 and PTEN, respectively. b. A snapshot of multicelluar simulation with 100 cells.
the frames for diﬀerent cells, to compute intercelluar interactions such as steric eﬀect or celladhesion, which is a computation bottleneck.
The resulting computational cost was roughly
the order of 10 × N L2 .
To overcome this problem, we improve thecomputational eﬃciency by updating the moving frames to adapt to the shape of cells. Figure 1a show the computational regions for twocells (gray rectangle areas). The area surrounded by a dotted line represent the region
required to compute the intercelluar interaction between the cells. Using this algorithm
and GPU, we have speeded up the computation by four times and realized simulation withorder of 100 cells as shown in Fig. 1b. In thefurther, we plan to apply our program to speciﬁc problems such ﬂuidity transition or cell295

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

sorting of a epithelial tissue.

References[1] D. Taniguchi, Ishihara, et al.: Proc. Natl.
Acad. Sci. USA 110, (2013) 5016.

296

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

Universal properties of dissipated many-bodyquantum systems
Tatsuhiko Shirai
Department of Computer Science and Communications Engineering, Waseda University
Ookubo, Shinjuku-ku, Tokyo 169-8555

We have studied two diﬀerent types of nonequilibrium quantum systems: bulk-dissipatedsystems and boundary-dissipated systems.

1

Bulk-dissipated systems

In the bulk-dissipated systems, the dissipativeenvironments are uniformly coupled to the system of interest. We consider that the systemenvrionment coupling is weak and the environment is nonequilibrium (i.e., the detailed balance condition is broken). The reference [1]has argued that the Gibbs state at an eﬀective temperature gives a good description of
the nonequilibrium steady state provided thatthe system Hamiltonian obeys the eigenstate
thermalization hypothesis and the perturbation theory in the weak system-environment
coupling is valid in the thermodynamic limit.
The numerics are based on the exact diagonalization method, which restricts the system sizestudied. We have developed a tensor-network
approach that directly studies not only thesteady states but also the relaxation dynamics in the thermodynamic limit [2]. The study
reveals that when an initial state is given by athermal Gibbs state, the states during the relaxation dynamics are well described by Gibbs
states with a time-dependent eﬀective temperature.

2

Boundary-dissipated
tems

sys-

We discussed the relaxation time of theboundary-dissipated systems, where the environments with diﬀerent thermodynamic potentials (e.g., temperature and chemical potential) are coupled to the edges of the nonintegrable systems. The presence of the diffusive transport indicates that the relaxation
time τ is proportional to the square of thesystem size L (i.e., τ ∼ L2 ). The relaxation time is naively related to the gap of thetime-evolution operator called the Liouvillion,
but the previous studies have shown that thegap closing is slower than L−2 [3]. We have
solved this discrepancy by numerically showingthat the relaxation time is not determined by
the lower-lying eigenvalues of the Liouvillion
but by superexponentially large expansion coeﬃcients for Liouvillion eigenvector with nonsmall eigenvalues at an initial state [4].

References[1] T. Shirai, and T. Mori: Phys. Rev. E 101(2020) 042116.
[2] H. Nakano, T. Shirai, and T. Mori:arXiv:2012.12274.
[3] M. Žnidarič, Phys. Rev. E 92 (2015)
042143.
[4] T. Mori, and T. Shirai: Phys. Rev. Lett.
125 (2020) 230604.

297

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

Quantum annealing in transverse-ﬁeld Ising chainswith a correlated disorder
Tatsuhiko Shirai
Department of Computer Science and Communications Engineering, Waseda University
Ookubo, Shinjuku-ku, Tokyo 169-8555

We have studied disordered Ising chains under transverse ﬁelds to examine the eﬀects ofa correlated disorder on quantum phase transitions. The Hamiltonian is given by
H=−

N
−1
∑
i=1

z
Ji σiz σi+1−N∑

Γi σix ,(1)

i=1

where ⃗σi = (σix , σiy , σiz ) are Pauli spin matrices, and Ji and Γi denote nearest-neighbor coupling strength and transverse-ﬁeld strength,
respectively. Ji is chosen from two diﬀerent
types of distribution: weak disorder and strong
disorder. The distribution in the weak disorder case is uniform over (J (0) , 1] and gapped(i.e., J (0) > 0), and thus given as{(1 − J (0) )−1 for J (0) < Ji ≤ 1,
πw (Ji ) =0 otherwise.
(2)
On the other hand, the distribution in thestrong disorder case is a gapless power-law distribution over (0, 1] given as
1 1 J −1+ D for 0 < J ≤ 1,iiπs (Ji ) = D(3)0 otherwise,
where D > 0 denotes the disorder strength.
In the correlated disorder case, the transverseﬁelds are given byΓ1= (1 − s) ln J1 ,
ln


Γ




 ln Γi = s ln J
i−1 + (1 − s) ln Ji
Γ
(4)


for
2
≤
i
≤
N
−
1,






 ln ΓN = s ln JN −1 ,Γwhere s ∈ [0, 1] is a parameter for tuning thetransverse ﬁelds.
In. [1], we analytically show that the dynamical critical exponents z in the system with thecorrelated disorder are ﬁnite: z = 1 in the
weak disorder case and max(D(1/2+|s−1/2|+
1/2), 1) ≤ z ≤ D + 1 in the strong-disordercase. We numerically estimate z in the strong
disorder case. The numerics is based on the
Jordan–Wigner transformation and an exact
diagonalization method [2]. The ﬁnite z is incontrast to inﬁnite z obtained in the uncorrelated disordered Ising chains (i.e., no correlation between the transverse ﬁelds and nearestneighbor couplings) [2, 3]. The suppression ofz is useful to enhance the performance of adiabatic quantum computations including quantum annealing with an argument based on the
Kibble–Zurek mechanism.

References[1] T.
Shirai,andS.
arXiv:2007.07439 (2020).

Tanaka:[2] A. P. Young and H. Rieger: Phys. Rev. B
53 (1996) 8486.
[3] D. S. Fisher: Phys. Rev. Lett. 69 (1992)
534.

298

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

Study of superstructure induced novel phenomena: stacking
of atomically thin materials
Toshikaze KARIYADOInternational Center for Materials Nanoarchitectonics,
National Institute for Materials Science, Namiki, Tsukuba, Ibaraki 305-0044
The purpose of this study is to propose ashift dependence of layer-layer distance, which

new design for van der Waals heterostructures

is derived by computing the binding energy as aof atomically thin materials to realize novel

function of layer-layer distance with rev-vdW-

phases or phenomena. Particularly, we focus on

DF2 type functional to take account of the van

twisted bilayers, where two layers are stacked

der Waal force. Once the layer-layer distance iswith angle mismatch. When the relative angle

fixed, the electronic structure is computed

is small, the system shows a long-range

using PBE-GGA type functional.

structure called moiré pattern that can induce

Overall, we have to scan over the possible

novel physics. With the small relative angle, a

relative in-plane shift, and therefore, it is

twisted bilayer system locally resembles to aimportant to make a calculation for eachbilayer system without twist, but globally, theparameter as light as possible. In that

relative in-plane shift between two layers

perspective, MPI parallelization implemented

depends on position due to the twist. Then, it isin Quantum Espresso package helps us a lot if

a powerful scheme to analyze crystalline andit is used on a massively parallel computational

electronic structures for untwisted bilayers andsystem.
In addition to the above project, we havescan over the possible relative in-plane shift.

analyzed

In this project, we have applied the above

diamagnetic

responses

inanprocedure on several candidate materials,

antiperovskite type Dirac electron system [1].

including graphene relatives and group-IV

The

monocalcogenide monolayers. For the analysis

compares well with the experimental result.

numericallyobtainedsusceptibilityof crystalline and electronic structures, we haveused the first-principles density functional

References

theory method implemented in Quantum[1] S. Suetsugu, K. Kitagawa, T. Kariyado, A.

Espresso package. One of the key quantities in

W. Rost, J. Nuss, C. Mühle, M. Ogata, and H.

crystalline structure analysis is relative in-plane

Takagi: Phys. Rev. B 103 (2021) 115117.

299

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

τl

(1)

1

S=τ =1/2

λ

Sl

Sl+1

1

(2)

τ =1

1
τl

p=1

1
msp

(1)

1

0.5

(2)

p=2
~0.807λc(∞)−

Figure 1: Structure of the diamond chain investigated in this work.
We investigate the ground-state phases of
mixed diamond chains described by the following Hamiltonian [1]:
H=

L [
∑

(1)

S l (τ l(2)+ τl )

l=1
]
(1)
(2)
(1) (2)+ (τ l + τ 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 the unit cells is denoted by L. Here,the parameter λ controls the frustration as depicted in Fig. 1. Deﬁning the composite spin
(1)(2)
operators T l as 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 whereTl = 1/2 and 3/2. The total Hilbert spaceof the Hamiltonian (1) consists of separated
subspaces, each of which is speciﬁed by a definite set of {Tl }. For large λ, ∀l Tl = 1/2.
Hence, the ground state is equivalent to thatof the uniform spin 1/2 chain, namely a gapless spin liquid. For λ ≤ 0, ∀l Tl = 3/2.
Hence, the ground state is a Lieb-Mattis ferrimagnetic phase with spontaneous magnetiza-

300

p≤38
0
0.6

0.7

p=3

p=4
λ

0.8

Figure 2: λ-dependence of msp calculated by
the inﬁnite size DMRG method.
tion msp = 1 per unit cell. For intermediate
λ, we ﬁnd a series of ferrimagnetic phases with
msp = 1/p where p takes positive integer values as shown in Fig. 2 based on the inﬁnitesize DMRG calculations for various conﬁgurations of {Tl }. The calculations for many different {Tl } are independent from each other.
Hence, they are suitable for the massively parallel system of the ISSP supercomputer center. The phases with p ≥ 2 are accompaniedby the spontaneous breakdown of the p-foldtranslational symmetry. It is suggested thatthe phase with arbitrarily large p, namely inﬁnitesimal spontaneous magnetization, is allowed as λ approaches the transition point to
the gapless spin liquid phase.

References[1] K. Hida, arXiv:2102.02116, to be published in J. Phys. Soc. Jpn.

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

Development of observation-noise estimation methodby machine learning
Ryo TAMURAInternational Center for Materials Nanoarchitectonics,
National Institute for Materials Science,
1-1 Namiki, Tsukuba, Ibaraki, 305-0044

To know the Hamiltonian of target materials, a data-driven approach is useful, in whichmodel parameters of Hamiltonian are determined so as to ﬁt the experimentally measured
data[1, 2, 3, 4]. In a data-driven approach, theposterior distribution is deﬁned, and the appropriate parameters are searched by a maximum a posteriori (MAP) estimation. Thus,
the Hamiltonian, which can well explain theexperimental results, is obtained by the framework of MAP estimation. However, by MAP
estimation, observation noise cannot be evaluated. The Hamiltonian estimation method
based on Bayesian inference is formulated under the assumption of ﬁnite observation noise.
Then, there must be errors in the estimatedparameters. Therefore, we proposed a method
to evaluate the error of the parameters by estimating the observation noise[5]．
In our framework, the noise amplitude σ isconsidered to be a hyperparameter, and a plausible value is determined by minimizing theBayes free energy F (σ), deﬁned as
F (σ) = − log Z(σ),(1)

where Z(σ) is the normalization of the posterior distribution given by(Z(σ) =

1
2πσ 2

)L ∫
2

dx exp [−E(x, σ)] . (2)ΩxHere, Ωx is the support of the posterior distribution determined by the prior distribution.
In addition, x is the model parameters in the

Hamiltonian and E(x, σ) is the energy function to be minimized in our model estimation.
To evaluate F (σ) , it is convenient to extend
the normalization factor with a “ﬁnite temperature”, which is deﬁned as(

Zβ (σ) =

1
2πσ 2

)L ∫
2

dx exp [−βE(x, σ)] , (3)Ωxwhere β is the inverse temperature. By usingZβ (σ), the Bayesian free-energy is deﬁned as
∫ 1

(

)

d
dβ
F (σ) = −log Zβ (σ) − log Z0 (σ)
dβ
0
∫ 1
)
(
L
dβ⟨E(x, σ)⟩β + log 2πσ 2
=
2
0
∫
− log

dx,

(4)

Ωx

where the ensemble average ⟨E(x, σ)⟩β can beobtained by the MCMC method. Here, thethird term on the right-hand side does not depend on σ and can be omitted. The noise amplitude of the experimental data σ ∗ is evaluated as the value where F (σ) is minimized.
This noise estimation method is applied forthe eﬀective model estimation of KCu4 P3 O12 .
As a result, we obtained the noise amplitude,
and the magnetic interactions with error bars
can be estimated as J1 = −8.54 ± 0.51 meV,J2 = −2.67 ± 1.13 meV, J3 = −3.90 ±0.15 meV, and J4 = 6.24 ± 0.95 meV, by using
our noise estimation method.
This work is the collaboration work withKoji Hukushima, Akira Matsuo, Koichi Kindo,and Masashi Hase.

301

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

References[1] R. Tamura and K. Hukushima, PLoS
ONE 13, e0193785 (2018).
[2] H. Takenaka, K. Nagata, T. Mizokawa,and M. Okada, J. Phys. Soc. Jpn. 83,
124706 (2014).
[3] R. Tamura and K. Hukushima, Phys. Rev.
B 95, 064407 (2017).
[4] H. Fujita, Y. O. Nakagawa, S Sugiura, and M. Oshikawa, Phys. Rev. B 97,
075114 (2018).
[5] R. Tamura, K. Hukushima, A. Matsuo,
K. Kindo, and M. Hase, Phys. Rev. B
101, 224435 (2020).

302

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

Systematic Investigation on Phonon Transport at
Nanoscale Interfaces
Susumu FUJII
Nanostructures Research Laboratory, Japan Fine Ceramics Center
Mutsuno, Atsuta, Nagoya 456-8587

Grain boundaries and hetero interfaces are
known to suppress lattice thermal conductivity of a material dramatically, because of theirdisordered structures diﬀerent from the corresponding crystal structure(s). Recently, there
is an increasing demand for revealing the mechanisms of interfacial thermal conduction at thenanoscale, as electronic devices have become
miniaturized and nanostructuring techniques
have been extensively developed. However,there is a few studies that investigate the impact of interface atomic structures or interfacial scattering mechanisms of phonons that aredeﬁned in the reciprocal space, making it diﬃcult to obtain a guideline for controlling thermal conductivity on the basis of interfaces.
In the present study, we have mainly investigated the interfacial thermal conduction
across silicon grain boundaries (Fig. 1). We
chose silicon as a model material because it has
many technologically important applications
such as electronic devices and thermoelectrics.
We performed two types of molecular dynamics simulations: (1) perturbed molecular dynamics and (2) phonon wave packet. The former adds a small magnitude of perturbation toatoms and calculate thermal conductivity fromthe average of heat ﬂux of the system. The
latter introduces an phonon wave packet generated from a single phonon mode from phonon
dispersions and simulate its scattering process
at the interface. The necessary codes weredeveloped by the author and implemented to
the Large-scale Atomic/Molecular Massively

Parallel Simulator (LAMMPS) [1]. In addition, we used a machine learning potential
for Si distributed in MACHINE LEARNINGPOTENTIAL REPOSITORY [2], which aremuch less computationally demanding compared with density functional theory calculations but more demanding than the empiricalinteratomic potentials such as Stillinger-Weberand Tersoﬀ potentials.

Figure 1: An example of silicon grain boundary structure derived from the machine learning potential.
The perturbed molecular dynamics simulations require a large number of timesteps
more than one million and large computational
cells with a few thousands of atoms. On theother hand, phonon wave packet simulations
require a smaller number of timesteps less than
100000 but the computational cells must bevery long for the direction perpendicular to thegrain boundary plane and thus contain tensof thousands of atoms. The former were per-303

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

formed with ﬁve diﬀerent magnitudes of perturbations, and the latter were performed for
40 - 60 diﬀerent phonon modes for each grainboundary. Each of these calculations was performed in a single node.
The results perturbed molecular dynamics simulations show that thermal conductivity of Si grain boundaries signiﬁcantly correlates with their microscopic atomic structures rather than their energies. On the other
hand, our phonon wave packet simulations reveal that the phonon transmission (or reﬂaction) at the grain boundaries does not depends
on the grain boundary structures much, especially for the acoustic phonons that transfer most of heat. This contradicting trends
obtained from two kind of molecular dynamics simulations suggest that anharmonic eﬀect
of atomic vibrations at the grain boundaries,
which is naturally included in the perturbed
molecular dynamics while is intentionally neglected in the phonon wave packet simulations,is one of the dominant factors for determinig
interfacial thermal conduction.
In additon, comparisons of grain boundary structures, phonon properties and thermal conductivity between the machine learning potential and a classical interatomic potential (Tersoﬀ potential) demonstrate the highpredictive power of the machine learning potential. By combining machine learning potentials with supercomputers and advanced computational programs, it is becoming possible toquantitatively analyze interfacial thermal conduction.

References[1] S. Plimpton, J. Comput. Phys. 117 (1995)
1.
[2] A. Seko, arXiv:2007.14206.

304

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

Developing numerical tool for open quantum dynamics basedon neural networks
Nobuyuki YOSHIOKA
Department of Applied Physics,
The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, Japan, 113-8656.
In this project, we have planned to use neural

have shown that, the thermodynamic limit of the

networks to study the property of open quantum

strongly correlated hydrogen chain can be

many-body systems which are subject tocomputed within chemical accuracy.

dissipations that give rise to competition similar

In the second work, we have developed two

to “frustration” in closed quantum many-bodyalgorithms

that

simulates

systems. While the proposed project is currently

equilibrium of an isolated quantum many-bodyon-going, we have achieved some other works

system using deep neural networks. In one

regarding scalable simulations using neural

method, we exploit the quantum-to-classicalnetworks. In the following, we report the resultscorrespondencein such two works [1, 2].

representation of Gibbs state using the deep

toconstruct

the

the

thermalexact

In the first work, we have successfully

Boltzmann machine (DBM). This demonstrates

demonstrated a novel simulation algorithm thatthe expressibility of the DBM. In the otherperform first-principles calculation of solid-statemethod, we approximate the imaginary-timesystems. Namely, we have shown that theevolution to generate finite-temperature density

expressive power of neural networks arematrix starting from the infinite-temperaturesufficiently high so that the essential properties

Gibbs state. We have numerically demonstratedof solids, the ground state and the quasiparticle

that the second approach opens a door toband spectra, can be simulated both efficiently

scalable simulation even in frustrated systems

and accurately. The ground states are computed

where the thermal property is still under debate.

by the standard variational Monte Carlo manner,

Referencessuch that the imaginary-time evolution is[1] N. Yoshioka, W. Mizukami, and F. Nori,

realized approximately. This enables us toarXiv:2010.01358 (2020).

simulate the potential curve of various real[2] Y. Nomura, N. Yoshioka, and F. Nori,

solids including polymers, non-organic crystals,

arXiv:2103.04791 (2021).

or strongly correlated systems. In particular, we

305

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

Quantum Monte Carlo study of theantiferromagnetic Ising model on a square latticewith longitudinal and transverse magnetic fields
Ryui KANEKO
Department of Physics, Kindai University
Higashi-Osaka, Osaka 577-8502, Japan
Analog quantum simulations have attracted
much interest recently. Ground-state properties and nonequilibrium dynamics of quantum
spin, SU(N) Heisenberg, and SU(N) Hubbard
systems have been extensively studied. It is
important to verify how far the simulation of
the experiment is correct by comparing numerical simulations on classical computers. It is aninteresting challenge to investigate the properties of complex models in static systems andthose in nonequilibrium systems using numerical methods such as tensor network algorithms.
On the other hand, it is also indispensable tostudy static properties in simple systems because it is much easier to compare the experiments and numerical simulations.
The Rydberg-atom systems are suitable forrealizing the fundamental spin models. The
quantum Ising models are frequently experimented with. The Hamiltonian is given as
XXXH=
Jij Siz Sjz − h
Siz − ΓSix . (1)
i,j

i

i

The longitudinal (transverse) field h (Γ) can
be controlled by the frequency detuning (theRabi frequency) of the laser [1]. The Ising interaction is given as van der Waals interaction,
whose long-range part is often ignored for simplicity [2].
Although the ground-state phase diagram of
the quantum Ising models on complex lattices
have been extensively investigated, the simple
model on a square lattice has yet to be ex-plored. The ground-state phase diagram was
obtained only for few dozen sites using the exact diagonalization method [3]. We investigate
the ground-state phase diagram of the modelon the square lattice using the quantum Monte
Carlo method with the DSQSS library [4]. We
typically choose 105 Monte Carlo steps for calculation of physical quantities and for thermalization. We use the system sizes Ns = L2
with L ≤ 32 and consider the periodic-periodic
boundary condition.
To determine the phase boundary, we calculate the staggered susceptibility defined ashM̂ z (Q)2 i, Q = (π, π),
(2)
βLd
Z β X
dτ
Ŝjz (τ )e−iq·rj , (3)and M̂ z (q) =χzzstag =

0j

where rj is the real space coordinate at site j,d(= 2) is the spatial dimension, and β is theinverse temperature.
We perform the finite-size scaling analysis
to determine the phase boundary between theantiferromagnetic and disordered phases. The
scaling form of χzz
stag is given as2−ηχzz
F(δL1/ν ),stag ∼ L(4)

where F is a scaling function, η is the anomalous dimension, and ν is the correlation length
exponent. The difference between the field andthe critical point is given as δ. Since the dynamical exponent would satisfy z = 1 for the306Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Disorder

h/J

2

1

0

AF order
0

0.5

1

1.5

Γ/J
Figure 1: Phase diagram of the antiferromagnetic Ising model on a square lattice with longitudinal and transverse magnetic fields.
Ising universality class, we choose the inverse
temperature β to be proportional to L during
the simulation.
Figure 1 shows the ground-state phase diagram of the mixed-field Ising model on thesquare lattice. When the transverse field Γ ∼0, the transition longitudinal field hc is nearly
intact as a function of Γ. On the other hand,when h ∼ 0, hc drops nearly vertically as afunction of Γ. The latter behavior is consis-

tent with the phase boundary obtained by themean-field approximation [5].
By scrutinizing the phase diagram near Γ ∼0, we have also uncovered the narrow region where the disordered phase exhibits reentrance. Our phase diagram [6] would be usefulfor examining the accuracy of the analog quantum simulation in the Rydberg-atom systems.

References[1] F. Robicheaux and J. V. Hernández,
Phys. Rev. A 72, 063403 (2005).
[2] A. Browaeys and T. Lahaye, Nat. Phys.
16, 132 (2020).
[3] V. Lienhard et al., Phys. Rev. X 8, 021070(2018).
[4] Y. Motoyama et al., Comput. Phys. Commun. 264, 107944 (2021).
[5] A. A. Ovchinnikov et al., Phys. Rev. B
68, 214406 (2003).
[6] R. Kaneko et al., arXiv:2103.12364

307

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

Theoretical study for caged compound and its
conjugate acid/bases
Miyabi HIYAMA
Graduate School of Science and Technology, Gunma University,
Tenjin-cho, Kiryu, Gunma 376-8515

Caged luciferins, which can generate luciferin by UV photolytic reaction, would be
useful to understand the firefly bioluminescence. We synthesized the new caged luciferin named coumarin caged luciferin and obtained its absorption spectra [1]. When we use
this caged-luciferin for spectroscopic studies,
we need the detail of electronic states of thismolecule in aqueous solutions.
The optimization structure of ground state
for coumarin caged-luciferin anion expected
to be main component in the pH 8 aqueous solutions was obtained from the density
functional theory. The time dependent DFT(TDDFT) calculations were carried out to estimate the theoretical absorption spectra forthis structure. These theoretical spectra wereused to assign the experimental spectra measured at pH 8 [2]. However, these theoreticalresults cannot explain the diﬀerence between
the shape of absorption spectra of coumarin
caged luciferin at pH 3 and those at pH 8. We
expect that the main chemical species in theacidic experimental condition would be neutral
coumarin caged luciferin molecule of which total charge is 0.

of water solute molecules. We found thatthe absorption energy for neutral coumarin
caged luciferin is larger than that for coumarin
caged luciferin anion. It was also found thatthe shape of absorption spectra for neutral
coumarin caged luciferin has a large main
peak, while that for coumarin caged luciferin
anion has a double peak.
All calculations were performed using the
GAUSSIAN09 program [3] on system B and
C of Super Computer Center in ISSP.

References[1] Kurata et al. , J. Photochem. Photobiol.
B, 189 (2018) 81.
[2] Usukura et al. Photochem. Photobiol.
(2020) Photochem. Photobiol. (2020) 96,
805.
[3] Gaussian 09, Revision D.01, M. J. Frischet al.

In this year, we performed the DFT calculations for the neutral coumarin caged luciferin
molecule. The excitation energies and the oscillator strengths at the ground state were obtained using TDDFT calculations for the optimization structure of the neutral coumarin
caged luciferin molecule. The polarized continuum model was used for the description

308

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

Exact diagonalization calculations of density of statesand dynamical structure factor for the
spherical-kagome spin-system {W72V30}
Yoshiyuki Fukumoto
Department of Physics, Faculty of Science and Technology, Tokyo University of Science
2641 Yamazaki, Noda, Chiba 278-8510

For the spherical-kagome cluster with thirty
spin-1/2s, as well as for the 2D kagome lattice,many low-energy singlet excitations have been
expected to exist in the energy region below
spin gap [1], which has been actually confirmed
by Kihara et al. in their specific heat measurements at low temperature in {W72 V30 } [2].
However, the experimental curve of the specific heat can not be reproduced by the theoretical curve in the Heisenberg model H =∑J ⟨i,j⟩ S i · S j , where ⟨i, j⟩ denote nearest
neighbors and the exchange parameter was estimated as J = 115K [1]. The experimental
curve does not have the peak around 2K, although the theoretical one has the 2K peak[1,2].
In last year, incorporating DzyaloshinskiiMoriya (DM) interactions and bond randomness into our model Hamiltonian, both ofwhich were originally proposed to explain theluck of the step-wise behaviour in the lowtemperature magnetization curve of {W72 V30 }[3,4], we used the method of thermal purequantum (TPQ) state [5] to calculate the specific heat. Then, we found that 10% of randomness collapse the 2K peak.
In this year, we calculate density of states,
entropy, and specific heat at low temperaturesby using the Lanczos method [6]. In particular,
we aim to understand low-temperature specificheat qualitatively via observation of the density of states. We find that DM interactions do
not significantly aﬀect the energy distribution309

of a dozen or so singlet states above the groundstate, which are involved in the peak structureof the specific heat around 2K. On the other
hand, we find that 10% randomness disperses
this distribution to collapse the 2K peak, whichis consistent with the above-mentioned TPQ
result.
In conclusion, we have clarified the similarities and diﬀerences between the eﬀects of
the DM interaction and bond randomness: theDM interaction breaks the conservation of thetotal spin and cause magnetic components tomix dense singlet states at low energies, while
the bond randomness disperses the energy distribution of singlets and triplets without breaking the conservation of the total spin, whichleads to similar eﬀects on the magnetization
process but does to diﬀerent eﬀects on the specific heat.

References[1] N. Kunisada et al., PTEP 2014, 41I01(2014).
[2] T. Kihara et al., PRB 99, 064430 (2019).
[3] Y. Fukumoto et al., J. Phys. Soc. Jpn. 87,
124710, 2018.
[4] J.
Schnacketmat/1304.2603v1.

al.,

arXiv:cond-[5] S. Sugiura and A.Shimizu, Phys. Rev.
Lett. 108, 240401(2012).
[6] H. Otsuka, Phys. Rev. B 51, 305 (1995).

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

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

We have constructed an algorithm for the
Ising machines. Ising machines are expected tobe highly eﬃcient solution machines for combinatorial optimization problems, but a signiﬁcant challenge is to develop an algorithm
that can exploit its potential and expand its
range of applications. With this background,
we have studied the following three topics.
(I) Minor-embedding algorithm forsimulated-annealing-based Ising machines
To input the Ising model, which represents
the objective function and constraints, into theactual Ising machine to solve the combinatorial
optimization problem using the Ising machine,
an operation called embedding is required.
The accuracy of the solution depends on thevalue of the hyperparameters in embedding algorithms. Therefore, appropriate hyperparameter tuning is necessary to obtain a highly accurate solution using the Ising machine. We
have proposed an appropriate method for adjusting this hyperparameter based on the viewpoint of statistical mechanics. We also veriﬁed
the properties of the proposed method by simulation. The results suggest that the proposed
method is superior to the commonly used hyperparameter adjustment method for embedding algorithms [1]. The work was done incollaboration with Tatsuhiko Shirai (Waseda
Univ.) and Nozomu Togawa (Waseda Univ.).
(II) Quantum annealing with correlateddisorder transverse ﬁelds
When solving a combinatorial optimization
problem in conventional quantum annealing, a
uniform transverse magnetic ﬁelds are applied

to all spins, and the transverse ﬁelds are gradually weakened on the same schedule. It wasrecently reported that the solution accuracy isimproved by introducing inhomogeneity in thetransverse ﬁelds. We have succeeded in obtaining rigorous inequalities for the upper andlower bounds of the dynamical critical exponents for the disordered one-dimensional transverse ﬁeld Ising model, which can be analyzed
rigorously through the free-fermion representation [2]. In addition, the dynamic criticalexponents are conﬁrmed by numerical calculations to see the ﬁnite size eﬀect. The work
was done in collaboration with Tatsuhiko Shirai (Waseda Univ.).
(III) Development of black-box optimization using Ising machines
We have previously proposed a method forblack-box optimization using Ising machinesand demonstrated it for automated metamaterial design [3]. We are currently working on
expanding the scope of this method and arepreparing a paper on it. The work was done incollaboration with Ryo Tamura (NIMS/Univ.
of Tokyo).

References[1] T. Shirai, S. Tanaka, and N. Togawa: IEEE
Access 8 (2020) 210490.
[2] T. Shirai and S. Tanaka: Annals of Physics,
Accepted.
[3] K. Kitai, J. Guo, S. Ju, S. Tanaka, K.
Tsuda, J. Shiomi, and R. Tamura: Physical Review Research 2 (2020) 013319.

310

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

Diversity in memory effects of flow in paste
Akio NAKAHARA
College of Science and Technology, Nihon University
7-24-1 Narashino-dai, Funabashi, Chiba 274-8501
A dense packed colloidal suspension, called

When NaCl is added to CaCO3 paste, thea paste, remembers the direction of its motion,

paste gets the ability of remembering flow

and its memory is visualized as a morphology ofdirection due to the screening effect by Cl - ions.

desiccation crack patterns. When a paste

Recently it is experimentally found that byremembers the direction of its vibrational

adding filtered starch solution into CaCO3 paste,motion, desiccation cracks propagate in thethe paste also gets the ability of remembering

direction perpendicular to the direction of theflow direction. However, if we add more filteredvibration, while, when a paste remembers its

starch solution into the paste, the paste loses theflow motion, cracks propagate along the flow

ability of remembering flow direction. To

direction, as are shown in Fig. 1 [1].

explain this phenomenon, the adsorption ofHere, CaCO3 paste remembers the direction

starch polymer onto colloidal particles should beof vibrational motion but cannot remember flow

included into the model, and numericaldirection. Note that CaCO3 particles are charged

simulations based on the model are inin water. Numerical simulations of shear motion

preparation.

of colloidal suspension using LAMMPS show
that, for the emergence of memory of flow, theattractive interaction should be dominant inshort

range interaction between colloidal

particles. In numerical simulations, colloidal
particles attract each other via Lennard-Jones
potential and receives Stokes’s drag force from

Fig. 1: Desiccation crack patterns induced bysurrounding fluid. Elongated clusters are formed

memory effect of vibration in (a) and flow in (b).

along flow direction, but when attractive
interaction is eliminated in the model, colloidal

References

particles could not form any elongated clusters[1] A. Nakahara, T. Hiraoka, R. Hayashi, Y.

along flow direction and that is why CaCO3

Matsuo, S. Kitsunezaki, Phil. Trans. R. Soc.

paste cannot remember flow direction.

A377:20170395 (2018).

311

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

Magnetic structures under the non-equilibrium stateof the magnetic thin films with the dipolarinteraction
Hisato KOMATSU, Yoshihiko NONOMURA, and Masamichi NISHINONational Institute for Materials Science, Namiki, Tsukujba, Ibaraki 305-0044

Magnetic structure of spin systems is an important subject of solid-state physics and statistical mechanics, which aﬀects the various
behaviors of systems. In particular, magnetic
friction, the frictional force caused by the magnetic interaction between spins, is an example
of non-equilibrium phenomena that the magnetic structure have an important role[1, 2].
We already introduced the model of magnetic friction which shows the crossover ortransition from the Dieterich-Ruina law to the
Stokes law in our previous study [3]. Spinsof this model interact with each other by theshort-range interaction, and we have not considered whether the behavior of the frictionchanges in the case of the long-range interaction system. However, consideration on general long-range interaction systems, including
the dipolar interaction system, is diﬃcult, because they make complicated magnetic structures depending on the condition[4] and needO(N 2 ) computational complexity. Hence, weshould first consider a relatively simple example of the long-range interaction system.
In this study, we introduce a model where
the spins interact with each other by theinfinite-range interaction. This model resembles that of our previous study, except for theinteraction range. The reason why we choose
the infinite-range interaction systems is thepoints that the computational complexity of
the numerical simulation is O(N ), and the behavior at the thermodynamic limit can be investigated by mean field analysis.

312

According to the numerical and theoretical study, this model always obeys the Stokes
law when the temperature is higher than thecritical value, Tc , whereas the short-range
model shows a crossover or transition from theDieterich-Ruina law to the Stokes law evenwhen the temperature is higher than the equilibrium transition temperature. This is because the frictional force of the present model
depends directly on the long-range order anddisappears when the order becomes zero. It is
the large diﬀerence between two models. In thecase that the temperature is lower than Tc , themodel shows a crossover or transition similar
to that of the short-range model when N < ∞,whereas it shows hysteresis dependence whenN → ∞. This seemingly contradictory behavior is caused by the divergence of the relaxation
time in the large-N limit.

Result of this study is published as Ref. [5].
Systems with realistic long-range interaction
such as the dipolar interaction should be studied in future works.

The numerical calculations of this study
were mainly performed on the Numerical Materials Simulator at National Institute for Materials Science. We applied for the usage of the
ISSP Supercomputer in preparation for transfer of aﬃliation. However, we have ended without any transfer, and have not used the ISSP
Supercmputer.

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

References[1] D. Kadau, A. Hucht, and D. E. Wolf:
Phys. Rev. Lett. 101 (2008) 137205.
[2] M. P. Magiera, L. Brendel, D. E. Wolf,and U. Nowak: Europhys. Lett. 87 (2009)
26002.
[3] H. Komatsu: Phys. Rev. E. 100 (2019)052130.
[4] K. De’Bell, A. B. MacIsaac, and J. P.
Whitehead: Rev. Mod. Phys 72 (2000)
225.
[5] H. Komatsu: Phys. Rev. E. 102 (2020)
062131.

313

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

Establishment of new analysis method
for extend X-ray absorption fine structurewith sparse modeling
Hiroyuki KUMAZOE
Institute of Industrial Nanomaterials, Kumamoto University
Kurokami, Chuo-ku, Kumamoto, 860-8555

Measurement of extended X-ray absorption
fine structure (EXAFS) is one of the primarymethods to obtain local structure information
around a specific element with atomic scale.
We are building a new method for the analysis
of EXAFS by the sparse modeling [1]. This
year, we focused on using new basis function
based on the multiple scattering theory.
To perform the sparse modeling on the EXAFS signals, we employed a new model of the
EXAFS formalism based on the multiple scattering theory [3]. we consider two-body terms
which are the main contribution of the EXAFS:Z ∞χ(k) =γ 2 (r, k)n(r) dr,(1)0

where, χ(k) is the EXAFS signal as a function of the photoelectron wavenumber; n(r) is
the coordination number of the distance r andγ 2 (r, k) contains the two-body scattering term
corresponding to the same geometrical configuration. We use a complex Hedin–Lundqvist
potential for inelastic loss effects in the framework of the multiple scattering theory. To obtain the local structure information, we solve
Eq. (1) with L1 regularization.
12ŵ = arg min ky − Xwk2 + λkwk1 , (2)2wwhere the regression coefficient is w. The response vector, y comes from the target data
and the predictor matrix, X incorporates other314

terms in Eq. (1). In this method, some elements of the coefficient, w is are suppressed toexactly zero with the moderate value of λ.
Let us consider the sparse regression
problem for the LASSO estimations within
Bayesian inference to choose the regression parameter, λ. In Bayesian inference, this linear
regression problem is optimized by maximizing
the posterior probability. Here, we introduce
the Bayesian free energy [2] as an information
criterion to optimize a regularization parameter and to extract the physical model appropriately. As a result, we estimated the radial
distribution function with fewer bases than inprevious studies [1].
This work was supported by JST CREST
Grant Number JPMJCR 1861, Japan.

References[1] I. Akai, K. Iwamitsu, Y. Igarashi,
M. Okada, H. Setoyama, T. Okajima, and
Y. Hirai, J. Phys. Soc. Jpn. 87, 074003(2018).
[2] Y. Igarashi, H. Takenaka, Y. NakanishiOhno, M. Uemura, S. Ikeda, andM. Okada, J. Phys. Soc. Jpn. 87, 044802(2018).
[3] F. Iesari, H. Setoyama, Y. Igarashi,
M. Okada, H. Kumazoe, K. Iwamitsu,
I. Akai, Y. Seno, and T. Okajima,
arXiv:2104.01805 [physics.data-an].

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

Nonequilibrium phase transition and slow dynamics
in the dense hard sphere systems
Masaharu ISOBE
Nagoya Institute of Technology
Gokiso-cho, Showa-ku, Nagoya, 466-8555

As one of the simplest models, the hard
disk/sphere systems have been investigatedin the field of both equilibrium and nonequilibrium statistical physics. In this project,
we investigated non-equilibrium phenomenain the hard disk/sphere model system withmodern algorithms, especially for Event-Chain
Monte Carlo(ECMC) [1] and Event-Driven
Molecular Dynamics(EDMD) [2], where we
propose the “Hybrid Scheme”, namely, ECMC
for equilibration and EDMD for calculation of
dynamical properties [3].
Direct Evidence of Void-Induced Structural
Relaxations in Colloidal Glass Formers:
The mystery of the microscopic origin of kinetic arrest and slow relaxation mechanismsof glass-forming materials has puzzled scientists for decades in condensed and statisticalphysics. In general, when a liquid is rapidly
cooled (compressed), its molecular arrangement remains disordered, and its motion becomes frozen, resulting in a glassy state. This
causes “slow relaxation” due to the increase inviscosity. Most of the particles are frozen anddo not move. Still, in rare cases, it is known
that “active particles” that cause large displacements are generated in a spatially heterogeneous manner and move cooperatively. The
cooperative motions of these “active particles”
are roughly classified into two types: cooperative motions as a core-like domain region andstring-like hopping chain motions, which havebeen observed in simulations at first. In thelatter type of hopping chain motion, a voidwith a size of diameter is required for the initial hopping. However, in high densities, sucha void does not exist, which is a paradox. Oneof the crucial issues is understanding whether
the essential properties of glass forming materials are fundamentally thermodynamic or dynamic in origin. One perspective that favors
a dynamic origin is called Dynamic Facilitation (DF) theory [4, 5]. However, the microscopic origin of structural (so-called α) slow
relaxation in deeply supercooled liquids at an
atomic scale has remained elusive due to thelimitation of electron microscopy experiments
and computer simulation.
Recently, we devoted to investigating theapplicability of DF theory to athermal systems, i.e., systems of hard particles where therelevant control parameter is pressure, under“super-compressed” conditions, using “Hybrid
Scheme” [6]. Besides, the optical microscopy
experiments on colloidal glass former system
play an important role because detailed motions of individual particles can be accessible.
International collaboration project with Hong
Kong, we precisely traced the movement of
individual colloidal particles for a long timein experiments As a result, (A) the colloidal
particles undergo a transition from collective“creeping” to the string-like “hopping” motion
on approaching the glass transition, and thestructural relaxation of the supercooled glass
is just caused by string-like escape hopping
motion. (B) In the supercooled glassy states,
small fragmented voids distributed over a few

315

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

particle distances, called quasi-voids, cooperatively promote the formation of large voids.
(C) The voids are transported like particles,
which drive a string-like hopping chain motion. In addition, we found (a) the cooperative
motion of the core-like domain regions of the“active particles” observed in previous studies
can be decomposed into a hopping chain motion with a higher temporal resolution, and (b)the hopping chain motion is reversed with highprobability and causes a string-repetition motion, which is proved to be the leading contribution to the kinetic arrest for slow relaxation.
These new findings provide fundamental insights into the structural relaxation of glassforming materials and may pave the way for a
full understanding of the microscopic mechanism [7].
Non-equilibrium response and slow equilibration in hard disk systems
The “equilibration” (relaxation toward equilibrium states) from initial non-equilibrium
states in the molecular simulations is one of the
crucial preliminary tasks to obtain the physical properties in “true” equilibrium. The issue such as crystallization and glass/jamming
transition have been actively studied recently.
It often requires a large-scale simulation indense molecular systems, which needs a longcomputational cost for equilibration. In general, the equilibration of particle positions indense systems is much diﬃcult due to the excluded volume eﬀect being dominant. To obtain the equilibrium states, it is a reasonablechoice to use ECMC for positional relaxation
at first, which is the central idea of a hybrid
scheme [3]. To elucidate the microscopic originof equilibration in a hard disk system, we investigate the relaxing process toward the liquidstates as a non-equilibrium response inducedby the disturbance of the homogeneous expansion. After such disturbance around the Aldertransition is induced, we performed EDMD
and estimated the relaxation time of four physical properties. As the preliminary results, we

do not expect, an anomalous slow equilibration
toward the liquid states to emerge when starting from the co-existence phase in large-scale
EDMD. To identify the physical mechanism of
anomalous slow equilibration, we found thatthe spatial inhomogeneity of the initial equilibrated phases would contribute to the relaxation time, which was confirmed by the probability distribution of local density and orientational order parameter[8]. We plan to investigate further by changing methodologies systematically.

References[1] E. P. Bernard, W. Krauth, and D. B. Wilson: Phys. Rev. E 80 (2009) 056704.
[2] M. Isobe: Int. J. Mod. Phys. C 10 (1999)
1281.
[3] M. Isobe: Mol. Sim. 42 (2016) 1317.
[4] D. Chandler and J. P. Garrahan: Annu.
Rev. Phys. Chem. 61 (2010) 191.
[5] A. S. Keys, L. O. Hedges, J. P. Garrahan,
S. C. Glotzer, and D. Chandler, Phys.
Rev. X, 1 (2011) 021013.
[6] M. Isobe, A. S. Keys, D. Chandler, and J.
P. Garrahan: Phys. Rev. Lett. 117 (2016)
145701.
[7] C.-T. Yip, M. Isobe, C.-H. Chan, S. Ren,
K.-P. Wong, Q. Huo, C.-S. Lee, Y.-H.
Tsang, Y. Han, and C.-H. Lam: Phys.
Rev. Lett. 125 (2020) 258001.
[8] D. Mugita and M. Isobe: EPJ Web Conferences, in press.

316

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

Quantum liquid crystals in the finite-field KΓ modelfor α-RuCl3
Masahiko G. YAMADA and Satoshi FUJIMOTO
Department of Materials Engineering Science, Osaka University
Toyonaka, 560-8531, Japan

We study the extended Kitaev model called
the KΓ model, using a perturbative expansion. From the Kitaev limit K = 0, thethird-order perturbation in Γ leads to interacting Majorana models, which are solved by awell-controlled mean-field approximation anda cutting-edge exact diagonalization. In theexact diagonalization we use a thread parallelization up to 112 threads. By this large-scalecalculation, the exact diagnolization up to 54
sites of Majorana fermions in the symmetrichoneycomb geometry is possible, while usually
the exact diagonalization has been done in thesmall system size up to 24 or 32 sites.
From the exact diagonalization combinedwith the mean-field theory, we found the following phases [1].

high-field state with zero Chern number in αRuCl3 [3]. This observation leads to a potential control of quantum states by a domain wallmotion in quantum liquid crystals.

References[1] M. G. Yamada, and S. Fujimoto:preparation.

in[2] S. Vijay, T. H. Hsieh, and L. Fu: Phys.
Rev. X 5 (2020) 041038.
[3] O. Tanaka et al.:(2020).

1. a Kekulé Kitaev spin liquid.
2. a non-Abelian chiral spin liquid.
3. an Abelian chiral spin liquid.
4. a nematic Kitaev spin liquid.
5. a Vijay-Hsieh-Fu surface code [2].
Especially, a Kekulé Kitaev spin liquid with a
Chern number 0 and a Vijay-Hsieh-Fu surface
code are new phases discovered in this largescale calculation.
The results imply the existence of a vast region of quantum liquid crystal states (a nematic Kitaev spin liquid and a Kekulé Kitaev
spin liquid). Both phases break the threefold rotation symmetry of the system, and potentially explains the experimentally observed

317

arXiv:2007.06757

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

Crowding Movement of Cells
Katsuyoshi Matsushtia
Depertment of Biological Science, Graduate School of Science, Osaka University
Machikaneyama, Toyonaka, Osaka 560-0043

Cells are indispensable constituents of theorganism’s body. Cells are distributed toproper positions to function in the various organs in the developmental processes of organisms [1]. In this distribution, cells sometimes
take a crowding state and move spontaneously.
This crowding movement is expected to finallylead to a jamming state of cells through theirsteric interactions [2, 3]. In contrast to thisexpectation, cells smoothly move, similarly to
fluid, and, further, enhances their motility as
shown in the observations of Dictyostelium discodeum, neuron, and blood cells [4]. The understanding of this crowding movement is a basic common issue in the fields of soft matterphysics and developmental biology.
Some of the eukaryote cells including Dictyostelium discodeum use the persistent memory of cell trajectory to stabilize their movements [5]. We additionally showed that the
memory induces the collective movement ofcells [6, 7]. The formation of collective movements is expected to have the same origin of
the self-propelled disks [8]. To get hints to clarify the mechanism of the crowding movement,
the comparison to the self-propelled disks is effective. The self-propelled disks have been wellconfirmed to stabilize collective movement in
the sparse condition through the inelastic multiple collisions. In this case, the memory of trajectories aligns movement directions of disks
in the collisions. In contrast, the investigation of the crowding self-propelled disks so far
has clarified the destabilization of the crowding movements and take a disordered state [9].
This result is contradictory to the crowding

movement of cells and implies that the crowding movements have an uncovered necessary
condition.
To clarify the condition of the cellular
crowding movement, we investigate the crowding movement on the basis of the cellular Potts
model [10]. At the first step of this investigation, we attempt to confirm that the crowding
model cells can stabilize their collective movements. To this end, we develop the model extension for dealing with the general polarity ofcells [11, 12] and by using this extension remodels the memory of cell trajectory [13, 14].
By using Monte Carlo simulation based on ourmodel, we successfully reproduced the collective movements even in crowding cells [15].
Further, the models show that the crowding
state rather stabilizes the collective movementof cells. More concretely, we show that thethreshold propulsion of the collective movement is reduced as the concentration of cells
increases. This result at least shows that thecrowding movement is possible for cells andclearly indicates that the model cells satisfies
the aforementioned necessary condition of thecrowding movements.
The key to clarify the condition is the originof ordering, namely the memory of trajectories.
To deeply examine effects of the memory, we
try to evaluate the relation between the memory and the stability of collective movement in
the crowding cells [15]. In a case, the crowding cells with movement show a solid to fluid
transition, and then their configurations highly
fluctuate. By using this transition, we confirm
the stability in the order of movement in the318Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

crowding model cells. We can show that evenat the transition point, the cells exhibit theordering of the movement under the high fluctuation of cell configurations. To further approach the origin of this stability [16], we consider memory time dependence of this movement order. We show that the stability depends highly on the memory time. Namely,
the long memory time is the origin of the stability of movements in the crowding cells.
Based on the comparison of this result andthat of the self-propelled disks, we attempt toshed light on the uncovered necessary condition. For the case of the self-propelled disks,
the movement has an intrinsic fluctuation. In
contrast to these disks, the model cells reduce the intrinsic fluctuation in their setting
to resemble the cells with the memory of celltrajectory. In the disks, the intrinsic fluctuation shortens the memory time of the disk
trajectory. Thus, the disks have a possibility
that the disks lose the stabilization propertyof the collective movements. This differencein the intrinsic fluctuation of movements reflects the difference in the stability in crowding
movements between the model cells and selfpropelled disks. Thus, we speculate that the
reduction of intrinsic fluctuation is the necessary condition of crowding movements.
These successes of the works are attributed
to kind cooperation with Prof. H. Kuwayama,
Prof. S. Yabunaka, Prof. K. Fujimoto, Dr.
H. Hashimura, Prof. H. Yoshino, and Prof.
M. Kikuchi. The author deeply thank them.

[5] H. Takagi, PLoS ONE 3 e2648 (2008).
[6] K. Matsushita, et al., Proc. Sympo. Traffic Flow Self-driven Parti. 24, 5 (2018a).
[7] K. Matsushita, et al., J. Phys. Soc. Jpn
88, 103801 (2019).
[8] C. A. Weber, et al., Phys. Rev. Lett. 110,
208001 (2013).
[9] T. Hanke, C. A. Weber, and E. Frey Phys.
Rev. E 88, 052309 (2013).
[10] F. Graner and J. Glazier, Phys. Rev. Lett.
69, 2013 (1992).
[11] K Matsushita, Phys. Rev. E 95, 032415(2017).
[12] K Matsushita, Phys. Rev. E 97, 042413(2018).
[13] B. Szabó, et al., Phys. Rev. E, 74, 061908(2006).
[14] A. J. Kabla J. Royal Soc. Interface 9, 3268
(2012).
[15] K. Matsushita, et al., , Proc. Sympo. Traffic Flow Self-driven Parti. 25 5 (2019).
[16] K. Matsushita, S. Yabunaka, K. Fujimoto
J. Phys. Soc. Jpn 90, 054801 (2021).

References[1] C. J. Weijer, J. Cell Sci. 122, 3215 (2009).
[2] S. Henkes, Y. Fily and M. C. Marchetti,
Phys. Rev. E 84, 040301(R) (2011).
[3] L. Berthier Phys. Rev. Lett. 112, 220602(2014)[4] P. Rørth, Ann. Rev. Cell Develop. Biol.
25, 407 (2009).

319

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

Multiple helical spin density waves ininversion-symmetric itinerant magnets
Takashi UCHIDA
Hokkaido University of Science
4-1, 7-15, Maeda, Teine-ku, Sapporo 006-8585

Since the discovery of magnetic skyrmions in
B20 transition metal compounds, multiple helical spin density waves have attracted attention because the magnetic skyrmions emerge
as 3Q multiple helical spin density waves
in these materials. Typically, the magnetic
skyrmions are realized in the DzyaloshinskiiMoriya (DM) interaction driven systems under magnetic field. Recently, frustrated systems without the DM interaction have beenfound to reveal magnetic skyrmions [1, 2]. In
order to explore theoretically the possibility ofvortex-type multiple helical spin density waves
such as magnetic skyrmions in the inversionsymmetric itinerant systems, we have applied
the molecular spin dynamics (MSD) method[3] to the triangular-lattice single-band Hubbard model.
The MSD method is based on the functional
integral method for the spin fluctuation theories and the isothermal molecular dynamics
method. The method allows us to find automatically the magnetic structure of a large system with thousands of atoms in a unit cell at
finite 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 lo-cal electronic structures are calculated in real
space by means of the recursion method. We
have adopted the automatic parallel calculation scheme and found it to be eﬀective in saving both computing time and CPU resources.
We have performed the magnetic structure
calculations on the supercell with 20×20 triangular lattice, which is embedded in a large cluster consisting of 3 × 3 supercells, each of which
are connected by the periodic boundary condition. Under zero magnetic field and the fixedvalue of the temperature T /t = 0.0005, wehave explored the magnetic structures changing the Coulomb interaction strength U/t andthe electron number n along the antiferromagnetic (AF)-ferromagnetic (F) boundary (n =1.3 ∼ 1.4) . We have performed rather
long-step MSD calculations (more than 15000
steps ) to find the thermodynamically stablestates among the various metastable multipleQ states and have found that for U/t = 5.0 ∼8.0 the 1Q helical spin density waves are stabilized in the vicinity of the AF-F boundary.
To explore the possibility of 3Q multiple helical spin density waves, the MSD calculations
for U/t < 5.0 are now in progress.

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.

320

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

Robustness of cluster states and stabilizer states forsurface codes against random local ﬁelds

By using the ISSP supercomputer, we numerically solve a time-dependent Schrödinger
equation of a many-qubit system and theoretically describe how random local ﬁelds affect cluster states and stabilizer states for surface codes [1]. While there are already several
known sources of decoherence in gate operationand measurement processes, we here focus ona new source of decoherence due to eﬀects ofﬁeld ﬂuctuations on idling qubits which should
become increasingly more serious as quantum
computers scale up. We ﬁnd similar temporal ﬁdelity degradation for both cluster statesand stabilizer states of surface codes for up
to ten qubits. We also ﬁnd that the eﬀect of
local-ﬁeld ﬂuctuations is greatly mitigated if
the magnitude of the ﬂuctuations can be suppressed to below 10% of the energy gap ∆ forboth cluster states and stabilizer states of thesurface codes. If the magnitude of the ﬂuctuations exceeds ∆/2, the state ﬁdelities forboth states deteriorate dramatically. A simple estimate based on the average ﬁdelity up
to time t ∼ 2h̄/∆ shows that the maximum
number of qubits that can suppress the system inﬁdelity below the 1% is less than 28 forthe surface code stabilizer states and 21 for thecluster states when the initial states are groundstates (Fig.1).

We have also used the ISSP supercomputer
to calculate the sucess probability of our proposed pulse sequence in the quantum annealingmechanism [2].

Average δ/∆ at fidelity=0.99

Tetsufumi Tanamoto
Department of Information and Electronic Engineering, Teikyo University
Utsunomiya 320-8511, Japan.

0.25

0.25

(a)

0.2

(b)

0.2
Cluster

Cluster

0.15

0.15

0.1

0.1

0.05

Surface

Surface

0.05
0

0
4

6
8
10
Number of qubits N

4

6
8
10
Number of qubits N

Figure 1: Average variation δ/∆ at which theﬁdelity is 0.99 for the cluster states and thesurface-code states. The error bars are calculated from the standard deviation of δ/∆.
The extrapolation equations for the cluster
states and surface-code states in (a) are givenby δ/∆ = −0.0053N + 0.1079 and δ/∆ =−0.0033N + 0.0924, respectively. Those in (b)are given by δ/∆ = −0.0004N + 0.0567 andδ/∆ = −0.0014N + 0.0732, respectively. The
extrapolation of these data points to the horizontal axis gives a rough estimate for the maxcsand Nmax, forimum number of qubits Nmax
which the errors can be corrected.

References[1] T.
Tanamoto
ArXiv:1910.05649

and

M

Ueda:

[2] T. Tanamoto: J. Appl. Phys. 84 (2015)
31415.

321

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

322

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

Comparison of rocking curves of 4 2 2 Bragg-reflected X-rays
from thin silicon crystals calculated based on the conventional
method and by solving the eigenvalue problem
Kouhei OKITSU
Institute of Engineering Innovation, School of Engineering,
The University of Tokyo, 2-11-16 Yayoi Bunkyo-ku Tokyo 113-8656

Fig. 1 Rocking curves of 4 2 2 Bragg-reflected X-rays from thin (14, 21 and 30 mm thick) silicon
crystals. (a), (b) and (c) were calculated based on the conventional method. (d), (e) and (f) wereobtained by solving the eigenvalue problem. (g), (h) and (i) were experimentally obtained [1].
The experiments were performed with the temperature carefully stabilized with the incidence ofhighly collimated and monochromatized X-rays generated from a Cu K source.

The present author has worked on the n-beam

Professor H. Hashizume kindly taught theEwald-Laue (E-L) and Takagi (T-T) X-ray

present author the existence of rocking curves

dynamical diffraction theories and numericalexperimentally obtained in 1970 [1]. He is one of

method to solve them for over twenty years [2,3].

theHowever, this report describes two-beam rocking

diffraction theory of X-rays. He also taught the

curves calculated based on the E-L theory andpresent author that the description of theexperimentally obtained over fifty years ago.

dynamical theory as the eigenvalue problem is not

323

rankingauthorities

of

the

dynamical

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

general and is not widely recognized.

dynamical diffraction theory.

The E-L theory was developed by Ewald (1917)

The conventional method to solve the X-ray two-and Laue (1931). The fundamental equation given

beam dynamical theory has two steps. At first, theby Laue is described as follows:

dispersion surfaces are calculated to give the𝐤 2𝑖 − 𝐊 2
𝐃𝑖 = ∑ 𝜒ℎ𝑖 −ℎ𝑗 [𝐃𝑗 ]⊥𝐤𝑖
𝐤 2𝑖

condition that the amplitudes of Bloch wave has(1)𝑗nonzero values. Next, for cross points (tie points)

Here, 𝐊 is the wave vector of the incident X-rays,of the dispersion surfaces and the downward

𝐤 𝑖 and 𝐤𝑗 are the wave vectors of the Bloch wave,

surface normal of the crystal, the amplitudes of𝐃𝑖 and 𝐃𝑗 are amplitude vectors of the 𝑖th and 𝑗th

Bloch wave are calculated to obtain the rocking

numbered Bloch wave, 𝜒ℎ𝑖 −ℎ𝑗 is the 𝐡𝒊 − 𝐡𝒋 order

curves.

Fourier coefficient of the electric susceptibility.

However, eq.(2) for n=2 (two-beam case) can beis the vector component ofsolved directly without calculating the dispersion[𝐃𝑗 ]⊥𝐤𝑖

𝐃𝑗perpendicular to 𝐤 𝑖 .

surfaces. The eigenvalue problem can be solved

By applying an approximation that 𝑘𝑖 + 𝐾 ≈e.g. by using the LAPACK to calculate the rocking

2𝑘𝑖 , the following equation is obtained [2,3]:
𝑛−1
1
𝐾
(𝑙)
(𝑙,𝑚) (𝑚)𝜉𝑖 𝐷𝑖 = ∑ 𝜒ℎ𝑖 −ℎ𝑗 ∑ 𝐶𝑖,𝑗 𝐷𝑗 .
2
𝑗=0
𝑚=0

curves of Bragg-reflected X-rays as found in Figs.
1 (d), (e) and (f). They are in excellent agreement(2)with those experimentally obtained [Figs. 1 (g), (h)(𝑙)is theand (i)] and calculated by the conventional way(𝑙,𝑚)is the[Figs. 1 (a), (b) and (c)].

Here, 𝜉𝑖 = 𝑘𝑖 − 𝐾 (the eigenvalue), 𝐷𝑖
scalar amplitude of Bloch wave and 𝐶𝑖,𝑗The present author is indebted to Mr. T. Sasaki(𝑙)polarization factor. The scalar amplitudes 𝐷𝑖(𝑚)

and 𝐷𝑗of Bloch wave are defined as follows:

who graduated from Tokyo Gakugei University in(0) (0)(1) (1)(3𝑎)2021 for his assistance when coding the program.

(0) (0)+ 𝐷𝑗 𝐞𝑗 .

(1) (1)(3𝑏)References

𝐃𝑖 = 𝐷𝑖 𝐞𝑖 + 𝐷𝑖 𝐞𝑖 ,𝐃𝑗 = 𝐷𝑗 𝐞𝑗

are unit vectors defined such that[1] H. Hashizume, K. Nakayama, T. Matsushita

they are perpendicular to 𝐤 𝑖 and mutually

and K. Kohra: J. Phys. Soc. Jpn. (1970) 29 806-

perpendicular to each other.

807.

(0)

𝐞𝑖

(1)

and 𝐞𝑖

Joko and Fukuhara pointed out [4], for the first[2] K. Okitsu, Y. Imai, and Y. Yoda: Recent

time, that eq. (2) can be described as an eigenvalue

Advances in Crystallography, (2012) pp. 67-86.

problem of a matrix for an n-beam dynamical

InTech Open. http://dx.doi.org/10.5772/47846.

theory. When n is 2, eq.(2) is nothing but the E-L[3] K. Okitsu: J. Jpn. Soc. Synchrotron Rad. Res.
33 (2020) 61-80 [Invited].

two-beam dynamical theory that is the most widespread X-ray dynamical theory. This fact has been

[4] T. Joko and A. Fukuhara: J. Phys.Soc. Jpn.
(1967) 22 597-604.

almost completely overlooked even by Ewald,
Laue and many other authorities of the X-ray324

Activity Report 2020 / 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
We have studied thermal transport properties(https://youtu.be/zCXmIXskBFE)

in proteins using molecular simulations [1-5].

is

under

development.

Within protein molecules, tightly packed amino

We implemented an accelerated CURP codeacid residues interact with each other throughvia GPU computing with the OpenACC library.

heat and energy exchanges. We illustrate non-

As a test calculation, this program was applied

uniform heat flow with our own model based on

to a sensory domain of an oxygen sensor protein,“local heat/energy conductivity” between each

FixL. As a result, the computation time of theresidue pair.

auto-correlation function of heat currents in

Native contacts in proteins are classified into

FixL was accelerated by 98.4%.

nonpolar, polar, and charged types. Harmonic

Human superoxide dismutase (SOD1) isspring picture applies to energy transfer (ET)known to bind metal ions for its function. To

through nonpolar and polar contacts, and ET

study the selective binding metal binding

rates vary inversely with the variance of themechanism of SOD1, QM/MM molecular

contact length: whereas diffusion picture is

dynamics simulation using Gromacs/DFTB is

relevant to ET through charged contacts and ET

under progress. We plan to analyze therates correlate inversely with the mean-square-fluctuations the metal binding sites using the

distance between charged atoms of a residue pair.

CURP program.

The CURP program permits to computeinter-residue flow or energy/heat and atomic

References

stress tensor in a protein, given atomic[1] P. Humanath, K. M. Reid, T. Yamato, D. M.

coordinates and velocity trajectories obtained

Leitner: J. Phys. Chem. B 124 (2020) 9852.

through molecular dynamics (MD). Energy flow[2] K. M. Reid, T. Yamato, D. M. Leitner: J. Phys.

data permit to picture an inter-reidue Energy

B. 124 (2020) 1148.

Exchange Network (EEN) as a graph. For[3] T. Yamato: Seibutsu Butsuri 60 (2020) 94.

interactive analysis of EEN graphs using[4] D. M. Leitner, T. Yamato: Biophys. Rev. 12

pointing devices, a new visualization tool, EEN(2020) 317.

VIEWER, with JavaScript and Python codes[5] T. Yamato: Biophys. Rev. 12 (2020) 291

325