Chapter 8

Computational Materials Science
Research Team
8.1

Members
Seiji Yunoki (Team Leader)
Yuichi Otsuka (Research Scientist)
Shigetoshi Sota (Research Scientist)
Shixun Zhang (Postdoctoral Researcher)Ahmad Ranjbar (Postdoctoral Researcher)
James S. M. Anderson (Foreign Postdoctoral Researcher)Kazuhiro Seki (Special Postdoctoral Researcher, joint)
Sandro Sorella (Senior Visiting Scientist)
Takami Tohyama (Senior Visiting Scientist)
Yutaka Imamura (Visiting Scientist)
Susumu Yamada (Visiting Scientist)
Michele Casula (Visiting Scientist)
Keiko Matsuoka (Assistant)

8.2

Research Activities

The main purpose of our team is to understand and predict quantum states of condensed mattersystems and solid state materials for the next generation highly functional materials and for thebasis of future revolutionary technologies. More emphasis is put specially on strong many-bodycorrelations that require the treatment going beyond the single-particle approximation. For thispurpose, we develop many-body numerical methods, which include numerically exact diagonalization
method, classical Monte Carlo and molecular dynamics methods, quantum Monte Carlo method,and density matrix renormalization group method. We also develop highly parallelized application
software based on these numerical methods to efficiently perform large-scale simulations.

67

Chapter 8. Computational Materials Science Research Team

Figure 8.1: Peak performance of benchmark calculations for the Hubbard model with various systemsizes.

8.3

Research Results and Achievements

8.3.1

Development of large-scale quantum Monte Carlo (QMC) simulations

The QMC method is one of the most powerful tools to investigate quantum many-body systems.
There are several variants of QMC techniques depending on target systems. Among them, we havedeveloped the determinant Monte Carlo method to study the interacting fermions on lattice, i.e.,the Hubbard model, with a specific aim to calculate physical observables with a high degree of
accuracy on quite large systems. The accuracy and the system size are inseparable in the latticesimulations as the energy resolution is roughly determined by the inverse of the linear size L of thesystem. Therefore, it is inevitable to perform simulations on very large clusters in order to accurately
calculate physical observables, which demands huge computational resources. Since the numericalcalculations are based mostly on linear algebraic operations such as matrix-matrix product andnumerical orthogonalization, we can take advantage of the highly optimized numerical library on K
computer. We have also optimized the delayed update algorithm to efficiently use the cache memory.
We have confirmed that our simulation code achieves 80% of the peak performance in the single nodecalculation with increasing the system size (see Fig. 8.1). Since in this case one Slater determinant,representing one ensemble in the Monte Carlo sampling, can be stored in the single node, we areable to use up to 24,576 nodes with quite high efficiency of about 50% of the peak performance intotal. We have utilized this code in the practical study of the two-dimensional (2D) Hubbard modelon honeycomb lattice with the system size up to N = 2, 596 sites, which is about 4 times larger thanthe previous study (N = 648) [Z. Y. Meng et al., Nature 464, 847 (2010)] and represents currently
the largest system size ever. The computational complexity scales as N 3 and thus it could havenot been possible without K computer. We have also implemented more non-trivial parallelization
for the case where one Slater determinant is distributed to a group of multiple nodes, in which theparallelized numerical library ScaLAPACK is employed to update the Slater determinant in eachMonte Carlo sweep. With this implementation, it has become possible to perform the simulations
with N up to 8,100 within reasonable CPU time.

8.3.2

Absence of spin liquid in the half-filled Hubbard model on the honeycomb latticeFirst, we have applied our improved QMC code to elucidate the ground state phase diagram of thehalf-filled Hubbard model on the honeycomb lattice model (honeycomb lattice model), in which agapped spin liquid (SL) phase was proposed previously [Z. Y. Meng et al., Nature 464, 847 (2010)].
Since it is widely believed that not only strong quantum fluctuations but also geometrical frustrations
are responsible for stabilizing a SL phase, the possibility of SL phase in the unfrustrated honeycomb
lattice is rather surprising, and thus has been one of the most debated issues in recent years. The

68

Chapter 8. Computational Materials Science Research Team

Figure 8.2: Finite size scaling of AF order parameter squared at U/t = 4.

Figure 8.3: Data collapse fits for (a) honeycomb lattice and (b) π-flux models.
proposed SL phase in the previous report was claimed to exist for 3.4 < U/t < 4.3 (U : Hubbard
interaction, t: nearest neighbor hopping) as a spin-gapped insulating phase without any broken
symmetry. We thus have first tried to clarify the existence of SL at U/t = 4, which corresponds tothe middle of the proposed SL region. Taking a full advantage of K computer, we have performed
the QMC simulations on the lattice with N = 2L2 up to 2,596 sites. Figure 8.2 shows our resultsof the antiferromagnetic (AF) spin structure factors, SAF , and spin-spin correlation functions at
the maximum distance, Cs (Lmax ), at U/t = 4. The extrapolated values of both quantities are
confirmed to be finite within statistical errors, indicating the AF long-range order. Complemented
with simulations performed at other U/t, our results strongly support the conventional scenario thata single and direct phase transition occurs between semi-metal and antiferromagnetic Mott insulator
with increasing U/t.

8.3.3

Quantum criticality of metal-insulator transition in 2D interacting
Dirac electrons

Next, having established the continuous character of the transition, we have performed careful finitesize scaling analysis of our QMC data to obtain critical exponents of this metal-insulator transition.
As shown in Fig. 8.3, the data of the staggered magnetization calculated on finite-size clusters are
excellently collapsed into a universal function in the honeycomb lattice model and also the Hubbard
model on square lattice with a magnetic flux π per plaquette (π-flux model). The latter model isknown to have massless Dirac dispersions as in the honeycomb lattice model. It is turned out thatthe critical exponents are exactly the same within the statistical errors for the two lattice models,
which strongly suggests that the metal-insulator transitions in these models belong to the same
universality class. This class should be coincident with that in the Gross-Neveu model, a model
extensively studied in the particle physics, since it has recently been recognized that the effectivemodel for the 2D interacting Dirac fermions in the continuous limit is described by the Gross-Neveu
model. Thus, we expect that our finding based on the unbiased large-scale QMC simulations havean impact also in an interdisciplinary field.

69

Chapter 8. Computational Materials Science Research Team

Figure 8.4: Performance of our massively parallel DMRG code on K computer.

8.3.4

Development of massively parallel density matrix renormalizationgroup (DMRG) algorithms

The DMRG method was originally proposed by S. White in 1992. Although the DMRG method
was soon recognized as one of the best numerically exact methods for strongly correlated many-bodyquantum systems, the application is severely limited mostly to one-dimensional (1D) systems becauseof the exponential increase of degrees of freedom in higher spatial dimensions. However, owning tohigh performance computers that are recently availably, it has become realistic to apply the DMRG
method even to 2D systems. We have been developing massively parallel DMRG algorithms toinvestigate strongly correlated quantum systems in two dimensions. We have already successfully
developed our massively parallel DMRG code that employs many efficient parallelization techniques,
and the first version of our massively parallel DMRG code has been opened in public for generalusers of K computer. A typical performance of the present version of our massively parallel DMRG
code is summarized in Fig. 8.4. For this example, we have calculated the optical conductivity for the
1D Hubbard model, using K computer with up to 82,488 nodes. We have performed the dynamical
DMRG calculations with the DMRG truncation number m = 8, 064 (the number of bases kept during
the calculation). The energy parallelization number is varied in such a way to correspond to theweek scaling. As shown in Fig. 8.4, we have achieved the extremely high peak performance ratio of
73.6% when 82,488 nodes are used on K computer, which is about 7.8 PFLOPS.

8.3.5

Magnetic excitations of the doped Hubbard model in two dimensions: Dynamical DMRG study

In order to disentangle magnetic excitations in the resonant inelastic X-ray scattering (RIXS) spectrum for electron-doped high-Tc cuprate Nd2−x Cex CuO4 (x: electron doping rate from half filling),
we have performed the dynamical DMRG calculations using our massively parallel DMRG code on K
computer. It has been reported experimentally that the peak position corresponding to the magnetic
excitation in the RIXS spectrum sifts toward higher energy with increasing doping rate x [K. Ishii,
et. al., Nat. Commun. 5, 3714 (2014)]. To understand this unexpected doping dependence, we havecalculated dynamical spin structure factor S(q, ω) for one of the simplest models for the high-Tc
cuprate , i.e., the 2D Hubbard model with long range hoppings on the square lattice. Using K computer, we were able to do the calculations for clusters up to 6x6 sites, which represents currently thelargest dynamical calculations in the world for the Hubbard model away from the half filling. We
have found that the lowest peak in the dynamical spin structure factor S(q, ω) shift toward higher
energy with electron doping, thus in good qualitative agreement with the RIXS observation.

8.3.6

Time-dependent DMRG method for real time quantum dynamics
in two dimensions

Up to now, no reliable numerical method has existed for the real time dynamics of 2D strongly
correlated quantum systems. We have challenged this issue by extending our massively parallel
DMRG method combined with the kernel polynomial expansion. The previously proposed time70

Chapter 8. Computational Materials Science Research Teamdependent DMRG method employs the Suzuki-Trotter decomposition for the time evolution operatoreiδtĤ . Although this scheme costs computationally quite small since the time-evolution operator isgiven a product of local operators, it is limited only to 1D systems. Therefore, we have proposed a new
time-dependent DMRG method combined with the kernel polynomial expansion to calculate directly
the time evolution operator in any spatial dimension. Employing P
Legendre polynomials Pl (x), the∞time-evolved state at T = t + δt is given as |t + δti = limL→∞ l=0 (−i)l (2l + 1)jl (δt)Pl (Ĥ)|ti,
where jl (x) is a spherical Bessel function. Since the special functions are satisfied with the threeterm recursive formula, the time-evolved state is calculated recursively. In the practical calculation,
we can truncate the number L of kernel polynomials to a finite value because the convergence of
the kernel polynomial expansion is quite fast when δt is small. With our new algorithm, we weresuccessfully able to perform the reliable real time dynamics simulations in two dimensions for the
first time. We have investigated, for example, the quantum annealing for the 2D Ising model withthe time-dependent transverse magnetic field, a relevant system for quantum information.

8.3.7

Development of ab-initio DMRG method for strongly correlated
molecular systems

Ab-initio DMRG method is developed in the quantum chemistry to perform the full configurationinteraction (full-CI) calculations. Traditionally, the exact diagonalization method is used in the
full-CI calculations. However, the degrees of freedom for strongly correlated molecular systems arequite large and increase exponentially with the number of orbitals considered. Thus, the exact diagonalization method is severely limited to small molecules. To overcome this difficulty, we have beendeveloping massively parallel ab-initio DMRG method for the full-CI calculations. The computational complexity depends strongly on the connectivity of the system, which relates to the area-low
of the entanglement entropy of the quantum state, and thus the order of the orbitals is crucial forthe efficient ab-initio DMRG calculations. We have developed a general scheme to give the mostefficient order of the orbitals in the DMRG calculation by using the graph theory. This allows us toperform the ab-initio DMRG calculations without explicitly considering the orbital order. Our abinitio DMRG method has been extended for the dynamical calculations using the dynamical DMRG
method.

8.3.8

Magnetic Skyrmions in chiral magnets: Electrons coupled with classical spins

Spins in chiral magnetic materials such as MnSi with a non-centrosymmetric structure are modulated significantly by Dzyaloshinskii–Moriya (DM) interaction and thus exhibit complex spin texture.
Notably, there exists a vortex-like spin configuration named Skyrmion that shows non-trivial topological properties such as topological Hall effect. Due to the potential application of Skyrmion in
the future information storage, Skyrmion has attracted intensive attention in recent years. We haveexamined chiral magnetic structure including Skyrmion using Monte Carlo and molecular dynamics
simulations.
The first class of systems studied is modeled by the double exchange model where the classicalspins are coupled to the non-interacting electrons described by a simple tight binding model. To
simulate this model, we have developed a semi-quantum Monte Carlo method, in which the electron
degrees of freedom are numerically solved exactly for a give spin configurations and the spin degrees
of freedom are treated using the importance sampling. Because the system size N has to be large
enough to accommodate many Skyrmions with typically tens of nanometers, the exact diagonalization of electron degrees of freedom requires huge computation cost of O(N 3 ). Therefore, we haveadopted a kernel polynomial method (KPM) to expand the electron Green’s function using Chebyshev polynomials, which reduces the computational complexity down to O(N ). Applying the KPM,
we were able to simulate up to 10,000 electrons coupled with 10,000 classical spins. The Skyrmion
phase digram with temerature vs magnetic field is shown in Fig. 8.5, which is in e xcellent agreementwith experimental observations.

71

Chapter 8. Computational Materials Science Research Team

Figure 8.5: Normalized Skyrmion number distribution in terms of magnetic field (B/t) and Temperature (T/t).

Figure 8.6: Spin structure of the 2D FM Heisenberg model with the easy axis spin anisotropy (K1)but no external magnetic field.

8.3.9

Classical Monte Carlo and molecular dynamics simulations for chiral spin texture

In addition, we have developed classical Monte Carlo and molecular dynamics simulations for pureclassical spin models. Because the spin interactions are only between the nearest neighbors, we canparallelize our codes very efficiently by separating spins into different sublattice groups, for which
spins in each group can be modified or updated simultaneously. We were capable to simulate over2 million spins for millions of MC sweeps or molecular dynamics time steps. We have studied thechiral spin texture of a 2D ferromagnetic (FM) Heisenberg model with the DM interaction and theeasy axis spin anisotropy. As shown in Fig. 8.6 (a), when the magnetic field and the anisotropy
are both absent, the stable spin structure is helical, where spins are aligned in strip configuration(the blue strips and red strips are separated by Bloch domain walls). The anisotropy breaks the
stripe structure to form larger FM domains separated by Bloch wall, as shown in Figs. 8.6 (b) and(c). When the anisotropy becomes very strong, isolated Skyrmions emerge in a FM background, as
shown in Fig. 8.6 (d). Our results explain why Skyrmions have been observed experimentally even
without magnetic field in some materials with strong easy axis anisotropy.

72Chapter 8. Computational Materials Science Research Team

8.4

Schedule and Future Plan

8.4.1

Large-scale quantum Monte Carlo (QMC) simulations for strongly
correlated electrons

QMC is one of the most powerful methods for quantum simulations in condensed matter physics,
especially when there is no notorious negative sign problem. We have developed a highly efficient
determinant QMC code to perform the ground state calculation for the two-dimensional (2D) Hubbard model on the largest lattice sites ever using K computer. We plan to continue to develop thedeterminant based QMC method for strongly correlated electrons at zero temperature as well as at
finite temperatures, to treat even larger lattice sites on exascale computers. One of the challenges
for the large-scale simulations is to study strongly correlated electrons in three dimensions. We alsoplan to develop the first principles QMC for strongly correlated materials and molecules throughinternational collaboration.

8.4.2

Density matrix renormalization group (DMRG) method for quantum dynamics in two dimensions

The DMRG method was originally proposed and is best performed for the ground state calculationsof one-dimensional (1D) quantum systems. Owning to high performance computers that are available
recently, it has become possible to apply the DMRG method not only to the ground state calculationsin 2D quantum systems but also to the excited state calculations, thermodynamics, as well as realtime dynamics. However, these dynamical calculations are limited mostly to 1D systems. This ispartially because of the exponential increase of degrees of freedom to be kept for 2D systems but
mainly because of the lack of reliable algorithms for thermodynamics and real time dynamics in two
dimensions. By combining the kernel polynomial method, we plan to develop a new scheme basedon the DMRG method for thermodynamics and real time dynamics in two dimensions. This will
provide the first reliable quantum dynamics calculation for 2D strongly correlated quantum systems.

8.4.3

First-principles calculations for strongly correlated materials

Strongly correlated materials are promising as highly functional materials for future technological applications. However, the standard first-principles calculations based on the density functional theory(DFT) is not reliable because the strong many-body interactions require the treatment that goes beyond the single-particle approximation. We plan to develop new first-principles calculation schemes
for strongly correlated materials. For strongly correlated solid state materials such as transition
metal oxides, we apply the many-body cluster approximation scheme such as dynamical mean field
theory (DMFT), variational cluster approximation (VCA), and cluster perturbation theory (CPT),
combined with the first-principles calculations based on the DFT. These cluster approximation cantreat the local electron correlation exactly within a cluster considered and therefore is apparently
superior to the single-particle approximation. We also plan to apply directly the DMRG method
as well as the tensor network states to develop an interdisciplinary research on strongly correlated
molecules such as metalloprotein, in collaboration with Computational Molecular Science Research
Team in AICS.

8.4.4

Tensor network scheme for strongly correlated quantum systems

Although the DMRG method is very powerful for 1D quantum systems, it meets the severe difficulty
due to the exponential increase of degrees of freedom that have to be kept to guarantee the numericalaccuracy when it is applied to quantum systems in higher spatial dimensions. This difficulty can bepartially overcome by aggressively using high performance computers. However, the fundamental
resolution has to rely on a radical revolution of the algorithm itself. Such a revolution has begun
to happen when condensed matter physics meets quantum information. Based on the idea of information compression, a tensor network scheme has been recently proposed to represent a quantum
many-body states with a product of tensors. There still exists a big challenge on how to treat thebasic tensor operations highly efficiently in large scale calculations. We will address this issue todevelop the renormalization group method based on tensor network states for strongly correlated
quantum systems in higher dimensions.

73

Chapter 8. Computational Materials Science Research Team

8.5

Publications

Journal Articles[1]

S. Hikino and S. Yunoki. “0-π transition driven by magnetic proximity effect in a Josephson
junction”. In: J. Phys. Soc. Jpn. 84, 024712/1-5 (2015).

[2]

S. Hikino and S. Yunoki. “Magnetization induced by odd-frequency spin-triplet Cooper pairs
in a Josephson junction with metallic trilayers”. In: Phys. Rev. B 92, 024512/1-13 (2015).

[3]

M. Khazaei et al. “OH-terminated two-dimensional transition metal carbides and nitrides as
ultralow work function materials”. In: Phys. Rev. B 92, 075411/1-10 (2015).

[4]

H. Matsuzaki et al. “Ultrafast charge and lattice dynamics in one-dimensional Mott insulatorof CuO-chain compound Ca2 CuO3 investigated by femtosecond absorption spectroscopy”. In:
Phys. Rev. B 91, 081114/1-5 (2015).

[5]

T. Sato, T. Shirakawa, and S. Yunoki. “Density-matrix renormalization group study of extended Kitaev-Heisenberg model”. In: Phys. Rev B 91, 054401/1-10 (2015).

[6]

K. Shinjo, S. Sota, and T. Tohyama. “Spin-orbit-induced exotic insulators in a three-orbital
Hubbard model with (t2g )5 electrons”. In: Phys. Rev B 91, 054401/1-10 (2015).

[7]

S. Sota, S. Yunoki, and T. Tohyama. “Density-matrix renormalization group study of third
harmonic generation in one-dimensional Mott insulator coupled with phonon”. In: J. Phys.
Soc. Jpn. 84, 054403/1-4 (2015).

[8]

T. Tohyama et al. “Enhanced charge excitations in electron-doped cuprates by resonant inelastic x-ray scattering”. In: Phys. Rev. B 92, 014515/1-8 (2015).

[9]

H. Watanabe, K. Seki, and S. Yunoki. “Charge-density wave induced by combined electronelectron and electron-phonon interactions in 1T-TiSe2 : A variational Monte Carlo study”. In:
Phys. Rev. B 91, 205135/1-6 (2015).

[10]

H. Weng et al. “Large-gap two-dimensional topological insulator in oxygen functionalized MXene”. In: Phys. Rev. B 92, 075436/1-7 (2015).

[11]Z. Zhong et al. “Giant switchable Rashba effect in oxide heterostructures”. In: Adv. Mater.
Interfaces 2, 1400445/1-5 (2015).

Conference Papers[12]

W. Fan and S. yunokii. “Electronic and Magnetic Structure under Lattice Distortion in SrIrO3/
SrTiO3 Superlattice: A First-Principles Study”. In: J. Phys.: Conf. Ser. 592, 012139/1-6. 2015.

[13]

H. Watanabe, K. Seki, and S. yunoki. “A Variational Monte Carlo Study of Exciton Condensation”. In: J. Phys.: Conf. Ser. 592, 012097/1-7. 2015.

Invited Talks[14]

J. Anderson. Relativistic formulation of the quantum theory of atoms and molecules. Mount
St. Vincent Theoretical Chemistry Meeting, Halifax (Canada). 2015.

[15]

A. Ranjbar. GGA and GGA+U modeling for atomic, electronic and magnetic structures of
Cr2 AlC. The 10th Anniversary General Meeting of Asian Consortium on Computational Materials Science (ACCMS), Sendai (Japan). 2015.

[16]

S. Yunoki. Block Lanczos density-matrix renormalization group method for Anderson impurity models: Application to pseudogap Kondo problems. Workshop and Symposium on DMRG
Technique for Strongly Correlated Systems in Physics and Chemistry, Natal (Brazil). 2015.

[17]

S. Yunoki. Carrier doped Mott insulator in iridates. The 1st International Conference on Computational Design and Simulation of Materials (CDSM 2015), Shenyang (China). 2015.

[18]

S. Yunoki. Carrier doped Mott insulators: iridates vs. cuprates. EMN Qingdao Meeting, Energy
Materials Nanotechnology, Qingdao (China). 2015.

74

Chapter 8. Computational Materials Science Research Team
[19]

S. Yunoki. Superconductivity and metal-insulator transition in 5d iridium oxides. The 10th
Anniversary General Meeting of Asian Consortium on Computational Materials Science (ACCMS), Sendai (Japan). 2015.

[20]

S. Yunoki. Universal quantum criticality in the metal-insulator transition of two-dimensional
interacting Dirac electrons. Novel Computational Methods for Quantitative Electronic Structure Calculations, Kobe (Japan). 2015.

[21]

S. Zhang. A Monte Carlo Study of Complex Magnetic Structures. The 5th International Symposium on Frontier Technology, Kunming (China). 2015.

[22]

S. Zhang. A Monte Carlo study of surface spin texture in a three-dimensional Skyrmion lattice.
Workshop on Topological Magnets, Wako (Japan). 2015.

75