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

Some Recent Developments in ab initio Thermodynamics ofIon Disorder in Solids
Shusuke KASAMATSU
Academic Assembly (Faculty of Science), Yamagata University
Kojirakawa, Yamagata-shi, Yamagata, JAPAN 990-8560

Abstract
Configurational disorder of ions in solids often
play key roles in determining the properties offunctional materials. A natural way to obtain thermodynamically relevant structures from the combinatorially increasing number of configurations in
multicomponent systems is to combine computational statistical physics methods such as Metropolis sampling with ab initio methods. Since ab initio methods incur high computational costs, theusual approach is to fit ab inito datasets to lowcost models such as cluster expansion, then use
those models for sampling. However, obtaining reliable low-cost models for complex oxide systems
have turned out to be challenging. To circumvent
this problem, we have recently been examining
the feasibility of direct sampling on ab initio energies, as well as turning to more flexible neural network models as an alternative to cluster expansion.
Here, we review some of our recent works in thisdirection; we also introduce our open source code
abICS, which is designed to harness the computational capacity of modern supercomputers to perform configurational sampling on complex multicomponent systems.

1

perature and chemical potentials. Such disorderoften plays significant roles in determining materials properties, and many materials in e.g., battery or magnet applications utilize impurities tocause disorder on purpose so as to obtain desirable
properties. Thus, the ability to simulate such disorder as a function of temperature should be immensely helpful for prediction and design of properties in materials with realistic levels of disorder.
To achieve this, it is necessary to perform sufficient thermodynamic sampling based on statistical mechanics. First-principles molecular dynamics is the method of choice for sampling quicklyrelaxing (e.g., liquid) systems, but it is nearly useless for sampling solid state systems with very
high energy barriers between different configurations. On the other hand, Monte Carlo methods
are not limited by realistic relaxation time scales
since they can employ “unphysical” trial steps suchas swapping of atom positions. However, standard Metropolis Monte Carlo sampling can still get
trapped in local minima, and it is usually not efficient enough to allow for a direct combination with
first-principles calculations. Thus, many workers
have opted to first derive lightweight models fromfirst-principles calculations, then use those models
for Monte Carlo sampling.

Introduction

Functional materials for various applications including energy conversion, nanoelectronics devices, and structural materials exhibit varying levels of configurational disorder which depend on
the details of materials processing including tem-

For modeling of systems that can be mapped
to a lattice (i.e., crystalline systems with configurational disorder), the cluster expansion approach
has seen much success [1–8]. This method expands the total energy in terms of “clusters” of16

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

Figure 1: Clusters within the conventional BCC
cell. From Chang et al. [8] (Creative Commons
Attribution License 3.0, https://doi.org/10.
1088/1361-648X/ab1bbc).
atoms (Fig. 1):
E(~
σ) =

X

mα Vα Φα (~
σ),

(1)

α

where α denotes empty, single, and multicomponent clusters, mα is the cluster multiplicity,~ is the occupation vector. Φα is the clusterand σbasis function, which for a binary system takes thesimple form of
Y
Φα (~
σ) =
σp,
(2)
p∈α

where p are site indices and σ p = ±1 dependingon which component occupies site p. Vα is the effective cluster interaction (ECI), which represents
the contribution of each cluster to the total energy.
This model reduces to the classical Ising model
when considering only nearest neighbor pair clusters. The ECIs are usually obtained by fitting to energies of relaxed structures from many small-scale
first-principles density functional theory (DFT)
calculations, then used to perform Monte Carlosampling on a larger cell.
The efficiency of the cluster expansion approach
depends on the ability to quantitatively reproduce
configuration energies using as few clusters (andcorresponding ECIs) as possible. This is rarely
an issue for simple two-component metallic alloys
where considering only near-neighbor interactions

often leads to sufficient accuracy. However, complex systems such as multicomponent oxides often
require many clusters due to longer range coulombinteractions as well as complicated strain patterns.
In addition to this, increasing the number of components (i.e., ion types) leads to a combinatorial increase in the number of clusters which makes cluster expansion computationally unfeasible. Anotherissue that has been recognized when using cluster
expansion is that surface calculations require many
more clusters than bulk calculations for sufficientaccuracy [9].
There are two possible ways to tackle the above
issue. One approach is to bypass the use of models and sample directly on first-principles energies. Although this approach is computationallyexpensive as mentioned above, the use of parallel sampling algorithms combined with modernday supercomputers have made it just barely feasible in recent years. Another approach is to utilize recently-developed machine-learning models
which is expected to be more flexible than the cluster expansion model. We have been working on
both fronts; we will review some of the literature aswell as our recent developments on the former approach in Sec. 2 and the latter approach in Sec. 3.
A software framework that incorporates these developments, which we have named ab Initio Configuration Sampling Toolkit, or abICS for short, isintroduced in Sec. 5.

2

Directcouplingoffirstprinciples calculations withreplica exchange sampling [10]There have already been some reports of firstprinciples Monte Carlo simulation in the literature, but they were limited to relatively simple
and quickly relaxing systems where first-principles
molecular dynamics may be just as efficient [11,12]. The issue, as noted above, is the tendency
to get stuck in local minima, especially at low
temperatures. One may think that with the increase in computer power available to researchers
today, simply performing many Metropolis sam-17

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

7HPSHUDWXUHH[FKDQJH

0&VDPSOLQJFRQWURO')7 9$63

5HSOLFD

')7 9$63

5HSOLFD

')7 9$63

5HSOLFD(QHUJ\
6WUXFWXUH(QHUJ\
6WUXFWXUH(QHUJ\

0&VDPSOLQJFRQWUROʞOn the other hand, we have been working extensively on combining first-principles calculation
with another popular extended ensemble method,i.e., the replica exchange Monte Carlo (RXMC)
method [17]. This approach considers an ensemble of Nrep replicas of a given system at Nrep different temperatures. Each replica is simulated in
the canonical ensemble, so the low temperature
replicas would still have the minima trapping issue. However, the problem is alleviated by allowing swaps of temperatures between replicas according to the following Metropolis criterion for

0&VDPSOLQJFRQWURO

6WUXFWXUHʞThis is where the so-called extended ensemble
sampling comes to the rescue. Using histogram
reweighting [13], it is possible, in theory, to convert between arbitrary ensembles. Thus, we are notlimited to sampling in a single canonical ensemble
at a given temperature, and this allows for choosing ensembles that can avoid the minima trapping
issue. For example, one can choose to sample
an ensemble where each configuration is realizedwith probabilities proportional to the inverse of thedensity of states. This is the idea behind multicanonical [14] and Wang-Landau [15] methods. In
fact, the Wang-Landau method was recently usedsuccessfully in combination with order-N firstprinciples Korringa-Kohn-Rostoker approach (albeit without structural relaxation) to calculate theorder-disorder transition in a 250-atom CuZn supercell [16].

0&VDPSOLQJFRQWUROʞpling calculations in parallel would lead to a sufficient number of samples for obtaining statistically
relevant results. However, such embarrassingly
parallel sampling is useless if each parallel image
cannot achieve equilibrium. The reason for getting trapped in local minima is that the Metropolis algorithm samples the canonical ensemble; agiven configuration is sampled with a probability
proportional to exp(−βE), where β is the inverse
temperature and E is the energy of the configuration. Thus, configurations residing at or near energy barriers on the potential energy surface will
rarely be accepted and the difficulty to cross suchenergy barriers increases exponentially with βE.

6WUXFWXUH')7 9$63

5HSOLFD1UHS(QHUJ\

Figure 2: Schematic of our scheme for coupling
the RXMC method with DFT codes. Reproduced
from Ref. [10]the ensemble of Nrep replicas:
P = min{1, exp[(βi − β j )(Ei − E j )]}.

(3)

Here, i and j are the indices of the replicas attempting the swap, and βi( j) and Ei( j) are the inverse temperatures and energies of each replica. According to this criterion, a replica with lower energy
will be assigned progressively lower energies, andvice versa. Essentially, a high-temperature replica
is given the task of a global scan of configuration
space; when it finds a new local minimum, it will
be assigned a lower temperature and will perform amore careful search within the newly found energy
basin. In this manner, the RXMC method provides
a good balance between global search and local optimization. At the same time, the canonical ensembles are easily recovered by collecting the replicas
at each temperature.
Figure 2 shows our computational scheme.
There are Nrep Metropolis samplers running at different temperatures; the samplers are coded using python and parallelized using mpi4py [18–20].
At each Metropolis step, each sampler performs aswap between atoms of different types (possibly
including vacancies), prepares input files for the
specified DFT code, and spawns parallel DFT processes using MPI-2 dynamic process management,i.e., the MPI COMM SPAWN function. The DFT

18

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

Figure 3: Comparison of the degree of inversion
as a function of temperature calculated by CESPCM (Seko et al. [21]; circles), the present approach (randomly initialized: squares; initializedwith an ordered spinel structure: triangles andexperimental neutron diffraction data [22]). The
spinel structure is shown in the inset. Reproduced
from Ref. [10] with addition of neutron diffraction
data.
processes perform local structural relaxation andenergy calculation, then the results are read automatically by the sampler to determine whether
to accept or reject the new configuration. The
samplers also attempt temperature swaps between
replicas at preset intervals according to Eq. 3. This
multilayered parallelism over replicas and DFT
processes makes it ideal for harnessing the powerof modern-day parallel supercomputers.
As a first benchmark for testing the feasibilityof this approach, we chose to calculate the degree of cation disorder in MgAl2 O4 spinel oxide[10]. Spinel oxides have a general formula ofA(II)B(III)2 O3 with A (B) representing a divalent(trivalent) cation, and are of interest in magnetic(spintronic) and electronic applications as well as
mineralogy. The structure is shown in Fig. 3. Usually, the divalent cation occupies the tetrahedral
sites and the trivalent cation occupies the octahedral sites. However, some inversion between thesites can occur, and some A(II)–B(III) combinations are known to exhibit completely inverse ordering, where all trivalent cations occupy the oc-

tahedral sites. The degree of inversion (DOI) is
quantified by the ratio of divalent cations at octahedral sites. Seko and coworkers have performed
calculations of the DOI in several spinels using
cluster expansion [4, 21], and we took their resultfor MgAl2 O4 as reference. We employed a calculation supercell containing 16 Mg, 32 Al, and64 O atoms, and sampled the configuration of Mg
and Al on the cation sublattice. No disorder on
the O sublattice was considered. We performed
two separate RXMC runs to check for convergenceand initial configuration dependence; one run was
initialized with random configurations and anotherwas initialized from the ordered spinel configuration. We used 16 replicas spaced evenly between
600 and 1500 K and performed 16750 steps forthe randomly initialized run and 13000 steps forthe run initialized from the ordered configuration.
This amounts to roughly 4000 node hours for eachrun on System C (enaga) at ISSP which has two
Xeon Gold 6148 processors per node. Figure 3
compares the calculated DOI vs. temperature forthe two runs as well as that obtained by a cluster
expansion model augmented by a screened point
charge model. We find the results to be virtually
identical, showing that sufficient sampling steps
were performed and it is indeed feasible to perform such calculations on meaningfully large supercells without resorting to fitted models. Moreover, the RXMC calculation initialized in the random configuration discovered the ground state ordered spinel configuration within 300 steps. We
think this is rather impressive considering that
there are 48 C16 = 2254848913647 possible configurations, and attests to the efficiency of the RXMC
method not only as a sampler but also as an optimizer. We also note that the deviation from experiment at lower temperature is seen in both our calculation and Seko’s. This is most likely due to thelow mobility of cations at low temperature which
inhibits the equilibration in the experiment.

19

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

Figure 4: The structure of the HDNNP. The
atomic coordinates are converted to a rotationallyand translationally invariant environment descriptor before feeding in to the neural network.

3

Accelerated sampling utilizing aneural network model combinedwith active learning [23]Although direct sampling using DFT-relaxed energies was shown to be feasible as detailed above,
the amount of required computer resources is quite
large, taking several days to a week to sample justone composition on a modern supercomputer. A
much higher throughput is necessary to enable materials design based on such simulations. Also,
because of the computational cost of DFT calculations, it is not feasible to go beyond a systemsize of about one hundred atoms and perform sufficient numbers of sampling steps. This is a problem when treating systems with, e.g., dilute doping
or long correlation lengths. To provide the necessary acceleration, we have no choice at this point
but to go back to lightweight models. However,as mentioned above, the de facto standard cluster
expansion approach has difficulty handling complex multicomponent systems that are of interestfor various applications.
To overcome this issue, we turned to the highdimensional neural network potential (HDNNP)
approach pioneered by Behler and Parinello [24,25]. This approach relies on the ansatz that the total energy can be expressed as a sum of atomic en-ergies determined by the environment around each
atom and uses a a feed-forward neural network to
represent those atomic energies (Fig. 4). The network weights are trained to reproduce DFT energies and forces (i.e., energy derivatives). Acceleration by 2–4 orders of magnitude compared to
DFT has been achieved while maintaining similarlevels of accuracy.
The HDNNP approach was developed for accelerating molecular dynamics simulations. There,
it is necessary to accurately reproduce the energy
and forces in continuous coordinate space. Thisis overkill when considering the lattice configuration problem; we only need to map the ideal lattice
structure with configurational disorder to the energy after local relaxation. Thus, we proposed totrain the HDNNP model to predict the relaxed energies from the set of atom coordinates and atom
species on the ideal lattice [23]. The new ansatz
here is that the total energy of a relaxed configuration can be expressed as a sum of atomic energies determined from atom environment descriptors calculated on the ideal lattice before relaxation. This may be expressed as follows:
Erel (
σ) =

atoms

i

 ∈ {
NNPrelσRi c ]) for σσlattice },ti ( f [(4) represents the coordinates of all atoms in
where σthe system, which is restricted, in our scheme, tothose on the ideal lattice. i represents atom in Ri cdices, ti is the corresponding atom type, and σrepresents the configuration of atoms within a cutoff radius from atom i. A separate neural net NNPti
is trained for each atom type; it takes as input
a rotationally and translationally invariant fingerprint of the atomic environment f [σRi c ] and outputs the contribution of the atom i to the relaxedtotal energy. In this work, we employed the Chebychev fingerprint proposed by Artrith and coworkers [26]. In this fingerprinting scheme, Chebychev expansion coefficients of the radial and angular distribution functions are used as the fingerprint. The atom type is encoded by multiplying with type-dependent factors when calculating
the radial/angular distributions. We use ænet code20

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo(http://ann.atomistic.net/) for the training
and evaluation of the neural network model.
An issue with neural network approaches ingeneral is that they are good at interpolating between training data but not at extrapolating. To
overcome this problem, we adopted the so-called“active learning” approach, where the neural network is retrained on-the-fly when the simulation
wanders into regions of structure space that were
not included in the original training data set [27–35]. Our approach is summarized as follows:
1. Prepare a training set consisting of DFTrelaxed energies of randomly generated configurations.
2. Perform neural network training.
3. Perform RXMC sampling using the trainedneural network model.
4. Perform DFT relaxations on a subset of configurations that appeared during RXMC sampling and check the prediction accuracy. If
sufficient accuracy is achieved, stop here.
5. If the accuracy is found to be insufficient, add
those configurations to the training set and repeat from step 2.
We benchmarked this approach, again, on the temperature dependence of the DOI in MgAl2 O4 . We
used a supercell with 192 cation sites compared to
48 in the previous section. The supercell size combined with the explosion of the possible number ofconfigurations (192 C64 ∼ 1051 ) makes this calculation completely out of reach for direct sampling on
DFT energies.
Figure 5 shows the correlation plot between theneural network predictions vs. reference DFT energies. The model trained on randomly generated
configurations performs well in the high energy region but shows a sizable deviation at lower energies. This is because the lower energy structureswith more ordering were obviously not included in
the original training data set. In our active learning scheme, these structures are then added to the

Figure 5: The correlation between reference DFT
energies and neural network predictions whentraining was performed on randomly generated
configurations (closed circles) and after performing active learning (open squares).
training set and the model is retrained. The resulting correlation after active learning is nearly perfect with deviations of less than a few meV/atom.
The active-learned neural network model was used
to sample the DOI vs. temperature (Fig. 6) from
480,000 RXMC steps with 15 replicas, and good
agreement was obtained with cluster expansion
and direct sampling on DFT energies explained in
the previous section.

4

Applications

We have only just started to apply the thermodynamic sampling methods outlined in previous sections to complex oxide systems that were difficult
to treat using previously available methods. Here,
we briefly present some results obtained for a couple of systems of interest in solid state ionics andelectrochemistry.

Acceptor-doped BaZrO3
Ceramic-based electrolytes are under intense research for application in battery, fuel cell, andelectrolyzer applications. Recently, we have been
focusing on acceptor-doped BaZrO3 , which is
known as a promising material to be used as the21Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

Figure 6: Comparison of the DOI calculated using
cluster expansion (blue circles), direct samplingon DFT energies (orange triangles), and the neural network model (open green squares) compared
to experiment (red triangles).
proton conducting electrolyte in solid oxide fuelcells. The material is an insulator without doping,but oxygen vacancies with an effective charge of+2 can be induced by substituting Zr4+ by trivalent
acceptor cations such as Y3+ . This is expressed inKröger-Vink notation as follows:
2Zr×Zr + O×O + Y2 O3 → 2Y0Zr + V••
O + 2ZrO2 . (5)These oxygen vacancies are immobile in the temperature range of interest, but highly mobile protons can be induced by a hydration reaction withthese vacancies:×•
H2 O + V••
O + OO → 2(OHO ).

(6)

The maximum dopant concentration can reach up
to 60% depending on the dopant species [36] (although it may be questionable whether it is appropriate to refer to such high concentrations as “doping”). Conventionally, dopants were assumed to
distribute randomly on the Zr sites. On the other
hand, some recent works have pointed out thatthe configuration of dopants can impact the proton
mobility in both positive and negative ways [37–39]. For example, an isolated dopant or a cluster of
dopants can act as a trapping site for protons due
to coulomb attraction. On the other hand, dopants
may form a long-range percolating pathway along

which protons can diffuse relatively freely. Thus,
we performed RXMC sampling combined directlywith DFT for revealing what the realistic dopant
configurations will be under processing conditions
in Y-doped BaZrO3 with dopant concentration up
to 30%. We note that many works have discussed
the interaction between few dopants and O vacancies or protons, but none have treated such highly
doped systems rigorously within a statistical thermodynamics framework like ours. Our approachrevealed that the dopant configurations are far fromrandom even at the usual sintering temperature of∼ 1800 K. This originates from complex manybody interactions that cannot be explained simply
in terms of coulomb interactions between charged
defects. Detailed analysis based on the statistics
obtained from RXMC samples can be found in
Ref. [40]. We also fed the obtained configurations into a master equation model for proton diffusion [41] and found that the realistic dopant configurations obtained by RXMC sampling predicts
slightly higher conductivities compared to random
configurations [42]. Another insight obtained thereis that in general, the formation of deep trappingsites due to dopant clustering has much more profound effect than the formation of a percolating
pathway. Thus, avoiding dopant clustering would
be the route to higher conductivity in this system.
We are now performing the active learning NNRXMC approach of Sec. 3 to understand the hydration behavior and the recently reported highconductivity of 0.01 S/cm at 400 ◦ C in 60% Scdoped BaZrO3 [36]. Prior to this report, Y dopingof ∼ 20% was reported to result in the highest proton conductivity using BaZrO3 as the parent material. An issue was that further doping results inthe decline of proton conductivity, which was also
reproduced by our work mentioned above [42]. In
the case of Sc, however, the conductivity continues to increase up to 60% doping, which is the solubility limit for this system. It will be interestingand insightful to see how the character of different
dopants with the same formal charge can lead tosuch different behavior. We are working to contribute in this regard using ab initio thermodynam-

22

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

(b)

Figure 7: (a) Evaluation of NN accuracy after active learning on The Pt/YSZ slab model and (b) the
obtained Y and O vacancy concentration profiles.
ics.

The ionic space charge layer
Another issue we are currently working on is theionic space charge effect at solid-solid interfaces.
The space charge concept is a central idea of semiconductor device physics, where a depletion or accumulation of holes or electrons occurs to alignthe Fermi level across the interface. The concept
should naturally be extended to mobile ions, which
are abundant in solid electrolytes, but the situation
is now much more complicated because various
ionic species can contribute to the effect, and ionscan even be introduced spontaneously (e.g., hydrogen in the atmosphere can spontaneously dope oxides [43]). A related (or maybe equivalent) concept in electrochemistry of liquid electrolytes isthe electrical double layer. The nanoscale atomistic details of the space charge layer is expected
to modify the ion dynamics in various ways. It
should impact the performance of solid state electrochemical devices such as solid-state batteries,
fuel cells, and electrolyzers, or may even lead tonovel functionalities. Because of this, the topic
is under intense study from both experimental
and theoretical standpoints (see, e.g., https://interface-ionics.jp/en/index.html).

In the past, we have combined ion defect formation energetics calculated using DFT with aPoisson-Boltzmann type model to calculate the ion
distribution near metal/solid electrolyte interfaces[44–46]. However, the energetics were calculated
in the bulk and do not include effects such as thechemical characteristics of the interface or the concentration dependencies of the dielectric constant.
Also, the interactions between defects were modeled considering only ideal point charge electrostatics. In other words, it is a mean field modelwith many assumptions and may not be reliable
enough for certain cases.
ab initio statistical thermodynamics approaches
introduced above are ideal for tackling this issue
and going beyond the simple mean-field picture.
As a first step, we calculated the ion distribution in
the Pt (111)/Yttria-stabilized zirconia (YSZ) (111)slab model (Fig. 7). Note that YSZ is a wellstudied oxide ion conductor for application in solidoxide fuel cells. The parent lattice is cubic fluorite
ZrO2 with 68 Zr sites and 126 O sites, and configurations of 12 Y ions on Zr sites and 6 O vacancies
residing on the O sites were sampled using our active learning NN-RXMC approach (Sec. 3). The
correlation between the DFT energies and NN predictions of the training, test, and verification sets

23

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

seems to be good enough, although not as perfectas the bulk MgAl2 O4 case (Fig. 5). The calculated
concentration profiles show a segregation of Y on
the outermost Zr layer and depletion on the second
and third Zr layers. There is also an increase in
the vacancy concentration on the second O layer.
The next challenge would be to consider variationsin the number of oxygen vacancies and perform
grand canonical sampling, so that the results can becompared directly with experiments which control
the oxygen partial pressure. We expect these calculations to be immensely useful for understanding the thermodynamics of space charge layer formation, and also to serve as a basis for understanding the modification of various materials properties
due to formation of interfaces.

5

ab Initio Configuration Sampling Toolkit (abICS)The software framework for direct couplingof RXMC sampling with DFT (Sec. 2) was
selected as an ISSP Project for Advancementof Software Usability in Materials Science for
SY2019. Parts of the prototype code developedby the author (S.K.) were rewritten in modular
fashion so as to enable easy extension of theframework. The original prototype only supported
VASP as the energy calculator, but interfaces
for OpenMX and Quantum Espresso are now
supported. In addition, a user interface employing
the human-readable and computer-friendly TOML
format (https://toml.io/en/) was developed
for setting up the calculations. The code has
been named ab Initio Configuration Sampling
Toolkit, or abICS for short, and released tothe public under an open source license (GPLver. 3), along with a user manual. It is under continous development (https://www.pasums.
issp.u-tokyo.ac.jp/abics/,
https://github.com/issp-center-dev/abICS),
and we are in the process of implementing theactive learning approach explained in Sec. 3.

6

Summary and outlook

In this report, we presented our recent efforts in enabling configurational sampling in complex oxidesystems from first principles. By harnessing thepower of modern day supercomputers in combination with extended ensemble sampling, it has now
become just barely feasible to perform sampling
directly on DFT energies on unit cells of ∼ 100
atoms. We also demonstrated that similar accuracy can be achieved with a speedup by a factor of∼ 104 by the use of a neural network model that
reproduces the configuration energetics. This has
made possible routine sampling on unit cells of afew hundred atoms with ab initio accuracy.
We have been implementing these approaches
in an open source framework abICS. We believethat this framework will be immensely useful formaking efficient use of the ever-increasing capacity of modern supercomputers. We are also aiming to bring together information science, materials simulation, and statistical physics using abICS
as a hub. Because of modular coding practices employed in our project, it should be relatively easy
to implement interfaces for other solvers, or to implement new sampling schemes. Please and do notrefrain from contacting us if you are interested inusing or extending this software.

Acknowledgements
abICS is developed by Yuichi Motoyama,
Kazuyoshi Yoshimi, Yoshiyuki Yamamoto, Osamu Sugino, Taisuke Ozaki, and the author (S.K.)
as part of the ISSP Project for Advancement ofSoftware Usability in Materials Science. The calculations presented in this work were performedon the ISSP Joint-Use Supercomputer Systems“sekirei”, “enaga”, and “ohtaka” as well as on the
supercomputer Fugaku provided by RIKEN for
the Fugaku Battery and Fuel Cell Project. This
report reviewed works supported over the years
by KAKENHI (Grant numbers JP18H05519,
19K15287, 20H05284), MEXT LEADER program, and JST CREST (JPMJCR15Q3) and24

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

FOREST (JPMJFR2037) programs.
Atomicstructure figures were drawn using VESTA [47].

[17] K. Hukushima and K. Nemoto, J. Phys. Soc.
Jpn. 65, 1604 (1996).

References[18] L. Dalcı́n, R. Paz, and M. Storti, J. Parallel
Distrib. Comput. 65, 1108 (2005).

[1] J. M. Sanchez, F. Ducastelle, and D. Gratias,
Physica A 128A, 334 (1984).

[19] L. Dalcı́n, R. Paz, M. Storti, and J. D’Elı́a, J.
Parallel Distrib. Comput. 68, 655 (2008).

[2] G. Ceder, Comput. Mater. Sci. 1, 144 (1993).

[20] L. D. Dalcin, R. R. Paz, P. A. Kler, and
A. Cosimo, Adv. Water Resour. 34, 1124
(2011).

[3] A. van de Walle, M. Asta, and G. Ceder, Calphad 26, 539 (2002).
[4] A. Seko, Y. Koyama, and I. Tanaka, Phys.
Rev. B 80, 165122 (2009).
[5] J. M. Sanchez, Phys. Rev. B 81, 224202(2010).
[6] J. M. Sanchez, J. Phase Equilib. Diffus. 38,
238 (2017).
[7] Q. Wu, B. He, T. Song, J. Gao, and S. Shi,
Comput. Mater. Sci. 125, 243 (2016).
[8] J. H. Chang, D. Kleiven, M. Melander,
J. Akola, J. M. Garcia-Lastra, and T. Vegge,
J. Phys.: Condens. Matter 31, 325901 (2019).
[9] K. Yuge, A. Seko, A. Kuwabara, F. Oba, and
I. Tanaka, Phys. Rev. B 76, 045407 (2007).
[10] S. Kasamatsu and O. Sugino, J. Phys.: Condens. Matter 31, 085901 (2019).
[11] J. Jellinek, S. Srinivas, and P. Fantucci,
Chem. Phys. Lett. 288, 705 (1998).
[12] S. Wang, S. J. Mitchell, and P. A. Rikvold,
Comput. Mater. Sci. 29, 145 (2004).
[13] A. M. Ferrenberg and R. H. Swendsen, Phys.
Rev. Lett. 61, 2635 (1988).
[14] B. A. Berg and T. Neuhaus, Phys. Rev. Lett.
68, 9 (1992).

[21] A. Seko and I. Tanaka, J. Phys.: Condens.
Matter 26, 115403 (2014).
[22] S. A. T. Redfern, R. J. Harrison, H. S. C.
O’Neill, and D. R. R. Wood, Am. Mineral.
84, 12 (1999).
[23] S. Kasamatsu, Y. Motoyama, K. Yoshimi,
U. Matsumoto, A. Kuwabara, and T. Ogawa,(2020), arXiv preprint: 2008.02572.
[24] J. Behler, Int. J. Quantum Chem. 115, 1032
(2015).
[25] S. Watanabe, W. Li, W. Jeong, D. Lee,
K. Shimizu, E. Mimanitani, Y. Ando, and
S. Han, J. Phys.: Energy 3, 012003 (2020).
[26] N. Artrith, A. Urban, and G. Ceder, Phys.
Rev. B 96, 014112 (2017).
[27] E. V. Podryabinkin and A. V. Shapeev, Comput. Mater. Sci. 140, 171 (2017).
[28] E. V. Podryabinkin, E. V. Tikhonov, A. V.
Shapeev, and A. R. Oganov, Phys. Rev. B 99,
064114 (2019).
[29] K. Gubaev, E. V. Podryabinkin, G. L. Hart,
and A. V. Shapeev, Comput. Mater. Sci. 156,
148 (2019).

[15] F. Wang and D. P. Landau, Phys. Rev. Lett.
86, 2050 (2001).

[30] J. S. Smith, B. Nebgen, N. Lubbers,
O. Isayev, and A. E. Roitberg, J. Chem. Phys.
148, 241733 (2018).

[16] S. N. Khan and M. Eisenbach, Phys. Rev. B
93, 024203 (2016).

[31] L. Zhang, D.-Y. Lin, H. Wang, R. Car, andW. E, Phys. Rev. Mater. 3, 023804 (2019).

25

Activity Report 2020 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo[32] T. Kostiuchenko, F. Körmann, J. Neugebauer,and A. Shapeev, npj Comput. Mater. 5, 1
(2019).

[46] K. Shimizu, W. Liu, W. Li, S. Kasamatsu,
Y. Ando, E. Minamitani, and S. Watanabe,
Phys. Rev. Mater. 4, 015402 (2020).

[33] T. D. Loeffler, T. K. Patra, H. Chan,
M. Cherukara, and S. K. R. S. Sankaranarayanan, J. Phys. Chem. C 124, 4907
(2020).

[47] K. Momma and F. Izumi, J. Appl. Crystallogr. 41, 653 (2008).

[34] J. Vandermause, S. B. Torrisi, S. Batzner,
Y. Xie, L. Sun, A. M. Kolpak, and B. Kozinsky, npj Comput. Mater. 6, 20 (2020).
[35] R. Jinnouchi, K. Miwa, F. Karsai, G. Kresse,and R. Asahi, J. Phys. Chem. Lett. 11, 6946
(2020).
[36] J. Hyodo, K. Kitabayashi, K. Hoshino,
Y. Okuyama, and Y. Yamazaki, Adv. Energy
Mater. n/a, 2000213 (2020).
[37] K. Toyoura, W. Meng, D. Han, and T. Uda, J.
Mater. Chem. A 6, 22721 (2018).
[38] K. Toyoura, T. Fujii, N. Hatada, D. Han, andT. Uda, J. Phys. Chem. C 123, 26823 (2019).
[39] F. M. Draber, C. Ader, J. P. Arnold, S. Eisele,
S. Grieshammer, S. Yamaguchi, and M. Martin, Nature Mater. (2019).
[40] S. Kasamatsu, O. Sugino, T. Ogawa, andA. Kuwabara, J. Mater. Chem. A 8, 12674(2020).
[41] K. Toyoura, T. Fujii, N. Hatada, D. Han, andT. Uda, J. Phys. Chem. C 123, 26823 (2019).
[42] T. Fujii, K. Toyoura, T. Uda, and
S. Kasamatsu, Phys. Chem. Chem. Phys. 23,5908 (2021).
[43] C. G. Van de Walle, Phys. Rev. Lett. 85, 1012
(2000).
[44] S. Kasamatsu, T. Tada, and S. Watanabe,
Solid State Ionics 183, 20 (2011).
[45] S. Kasamatsu, T. Tada, and S. Watanabe,
Solid State Ionics 226, 62 (2012).

26