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

Training method for reﬁnement of machine-learninginteratomic potential and its applications
Kohei Shimamura, Akihide Koura, and Fuyuki Shimojo
Department of Physics, Kumamoto University,
Kumamoto 860-8555, Japan

Abstract
Machine-learning interatomic potential (MLIP) isbased on low-cost ML models trained with ﬁrstprinciples data to achieve high accuracy in molecular dynamics simulations, expanding the range ofapplicable time-space domains. However, the current widely used training method has limitationsin accurately estimating stress and fundamental
physical quantities such as free energy and thermal conductivity. We have therefore emphasized
the importance of stress training for MLIPs andinvestigated other issues behind the current training method. The report discusses the role of stresstraining in reproducing solid-liquid phase transitions and estimating thermal conductivity, as wellas the effect of training data with different k points
on free energy.

ing on the reﬁning the training method of MLIPs.
In particular, we have emphasized the importanceof stress training to MLIP [4, 5]. Recent mainstream training methods set total potential energy
and atomic force as the targets in training [6]. Despite the fact that the accuracy of stress is not guaranteed, some previous studies to investigate highpressure phenomena have utilized such MLIPs. It
is moreover found that the stress training is necessary to accurately estimate fundamental physical quantities such as thermal conductivity (TC),and becomes one of the indispensable training
conditions [5]. This report discusses the importance of stress training of MLIPs for reproducing aphase transition and estimating TC. In addition, we
also discuss the effect of training data with different k points on thermodynamic quantities obtained
through MLIPs, such as free energy. As the MLIP,
we adopted an interatomic potential using an artiﬁcial neural network (ANN potential).

1 Introduction
The machine-learning interatomic potential(MLIP) is used in a wide range of researchﬁelds. By training low-cost machine learningmodels with ﬁrst-principles molecular dynamics(FPMD) data, the MLIP achieves high accuracy,
signiﬁcantly expanding the applicable space-time
domain of MD simulations. MLIP is used to
calculate physical quantities that require large
statistical amounts, to investigate phenomena
occurring in long-term dynamic simulations
such as shock waves, and to explore reaction
pathways [1, 2, 3].
Under these circumstances, we have been work-

1.1

Stress training

Originally MLIPs were trained with only total potential energy. Afterward, atomic force, which is
crucial for MD simulations, was integrated into
the training process, leading to a signiﬁcant improvement in accuracy compared to MLIPs that
only consider the total potential energy. Recently,
some MLIPs have been obtained by training with
FP data of the total potential energy, atomic force,and stress. However, such MLIPs that target stressin training may still be in a minority. In this report,
we show our results that the accuracy of stress is
not guaranteed by training with the total potential

15

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

energy and atomic force alone, and that the solidliquid phase transition of Na is accurately reproduced by training with stress [4].

1.2Free energy calculation

Thermodynamic integration (TI) is a popularmethod for computing free energy, but when combined with MLIPs such as ANN potentials, someconsiderations need to be taken for liquids. We
discovered a method that involved two steps, using an ideal gas and a soft potential as reference
systems [7]. We show how to calculate the free energies of solid and liquid states by the TI method
using ANN potential for Na [4].
Furthermore, in numerous previous studies,
MLIPs were trained using FP data with multiple
k points instead of using of large-scale systems.
We however show that as precision requirements
increase, larger systems with a greater number of kpoints will be necessary for the free energy calculation [4].

1.3

TC calculation

The TC calculation method based on the GreenKubo (GK) formula has been widely used because
it can be applied to even disordered systems. However, MLIPs including the ANN potential belong
to many-body potentials, and the many-body effects have to be taken into account in the formulaof the heat ﬂux [8]. We found that using the rigorous heat ﬂux formula considering the effects together with applying stress training to the ANN
potential led to highly accurate TC estimates [9].
Here, we show the result of the TC for the superionic conducting Ag2 Se (α -Ag2 Se) by homogeneous nonequilibrium MD (HNEMD) based on the
GK formula [9].

2Computational Details

In this study, FPMD and MD simulations using MLIPs were conducted using the QXMD
code [10].

2.1

ANN potential

The ANN potentials were constructed usingænet [11], which is a training code developed byArtrith et al. The ANN potential comprises feed
forward neural networks (FFNNs) created for eachatomic species. The total potential energy U ANN
is deﬁned by the sum of the atomic potential energies {εiANN } output from FFNNs for all atoms inthe system [12].
Natom

U ANN = ∑ εiANN ,(1)iwhere Natom denotes the total number of atoms.
The atomic potential energy depends on the
structural descriptor used. Chebyshev descriptors [13] were used in this study. The atomic potential energy of the Chebyshev descriptor is only
a function of the relative coordinates ri j [5]()εiANN = εiANN {ri j } j̸=i .
(2)

2.2

Training methods

The following cost function C was used to train the
ANN potentials. It consists of three loss functions,
i.e., the total potential energy, atomic force, andvirial.
()2
pU 1 NI UIANN −UIRef
C =2 NI ∑
Natom
I
+

+

)2
pF 1 NI1 Natom ( ANN
FI,i − FRef∑∑
I,i
2 NI I 3Natom i()2
WI,ANN−WI,RefpW 1 NI 1 6
j
j
, (3)
∑
2 NI ∑NatomI 6 j

where NI denotes the number of training data. The
symbol “Ref” denotes the training data. Factor 6
for the virial term reﬂects the number of independent degrees of freedom of the virial tensor WANN.
I
Because the three loss functions differ in dimension and size, pU , pF , and pW are introduced as
adjustable parameters. With appropriate settingsof these parameters, highly accurate ANN potentials can be constructed [14].
Note that the training for virial tensor WANN
isI
ANN
equivalent to that for the stress tensor PI
calcuANN
lated by dividing WI
by supercell volume.

16

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

The ANN potential trained with total potential
energy and atomic force is called EF-ﬁt, while theone trained including virial is called EFP-ﬁt.

2.3

Generation of training data on Na

We obtained the training data of Na by running
FPMD simulations. We calculated the electronic
states using the projector-augmented-wave (PAW)
method [15] within the framework of density functional theory (DFT) [16], in which the generalized
gradient approximation (GGA) [17] was used forthe exchange-correlation energy. The plane-wave
cutoff energies were 10 and 50 Ry for the electronic pseudo-wave functions and pseudo-charge
density, respectively. Projector functions weregenerated for the 3s, 3p, and 3d state of Na atom.
The FPMD simulations were mainly performed
using 128 Na atoms with periodic boundary conditions. For the Brillouin-zone sampling, 4 and32 special k points were utilized in addition to theΓ point. These are equivalent to sampling at theΓ point in systems consisting of 1024 and 8192
atoms, respectively. For comparison, FPMD simulations were also conducted on 1024 Na atoms
using the Γ point. The equations of motion weresolved under the isothermal-isobaric (NPT ) ensemble through an explicit reversible integrator[18]. The pressure was set to be ambient. The MD
time step ∆t was 2.9 fs. The simulation temperature was mainly set to 350 K. In addition, FPMD
data using 128 atoms with Γ point were preparedfor heating and cooling processes in the range between 450 and 200 K to investigate whether the
ANN potential can reproduce the melting and solidiﬁcation processes of Na.

2.4

TI method

In the TI method [4, 7], the interaction energyU(λ ) with parameter λ and the Helmholtz free energy F(λ ) are connected by the following relation:∂ F(λ )=∂λ⟨∂ U(λ )
∂λwhere ⟨· · · ⟩NV T denotes the canonical (NV T ) ensemble average. U(λ ) is given asU(λ ) = λ UNa + (1 − λ )Uref (0 ≤ λ ≤ 1),

where UNa and Uref are the total potential energiesof the target (Na) and reference systems, respectively. By substituting Eq. (5) into Eq. (4), andintegrating both sides with respect to λ from 0 to
1, we obtain the following formula:F(λ = 1) − F(λ = 0) =∫ 1
0⟨UNa −URef ⟩NV T d λ .

(6)F(λ = 1) and F(λ = 0) are the free energies of
the target and reference systems, respectively. As
reference systems for the solid and liquid phases,
we employed an Einstein solid and an ideal gas,
respectively.
In the liquid phase, the soft core (SC) potential USC , deﬁned as follows, was employed to
avoid shortening the interatomic distance as λ approaches zero [7].
N

USC ({rr i }) = ∑ µ
i< j

(

σ
ri j

)n
.

(7)

The parameters were used to be µ = 7.92 ×
10−4 hartree, σ = 6.752 bohr, and n = 12. The
Helmholtz free energy in the liquid phase was calculated by two TIs. The ﬁrst TI is calculated along
the ideal gas to SC liquid and the second one isfrom SC to the target liquid Na. In addition, to
integrate the ﬁrst while preventing divergence atsmaller λ , we adopted the coordinate transformation from λ to x proposed in the previous study
[19]
(
) 1
x + 1 1−kλ (x) =,
(8)2
where the interval of integration is changed from0 ≤ λ ≤ 1 to −1 ≤ x ≤ 1. We employed the parameter k = 0.85. Applying this transformation to
the integration with respect to λ , we have the following formula:
∫ 1

⟩
,

(5)

0

(4)

NV T

17

f (λ )d λ =12(1 − k)∫ 1
−1

f (λ (x))λ (x)k dx .
(9)Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

2.5

Generation of training data on α Ag2 Se

We generated training data using MD simulations
with an empirical interatomic potential (EIP) instead of FPMD. One reason for this is that we
can obtain a reference value of TC. The other reason is that the EIP of Ag2 Se proposed by Rino etal. [20] is a physically excellent EIP that can describe phase transitions among nonsuperionic, superionic conducting, and liquid phases [21]. The
training data were obtained by MD simulation with
the NV T ensemble at 500 K for the α -Ag2 Se system that comprises 256 Ag + 128 Se atoms in experimental density [22].

the heat ﬂux formula and showed that the atomic
virial is expressed as follows.
WANN-I= ∑ ri j ⊗ij̸=i

2.6.1 GK formula and heat ﬂux
According to the GK formula [23, 24], the αβcomponent of the TC tensor καβ can be estimatedas (α , β ∈ {x, y, z})

καβ

=

Ω
kB T 2∫ ∞0

WANN-II= ∑ ri j ⊗ij̸=iwhere T , kB , and Ω are temperature, the Boltzmann constant, and volume of the supercell, respectively. ⟨JQ,α (t)JQ,β (0)⟩ denotes the autocorrelation function of the heat ﬂux JQ , and thefollowing expression is employed for both the EIP
and ANN potentials:
JQ =,∂ εiANN,∂ r ji[

1 Natom
1 Natom
1 Natom
ε
t
v
+
v
+
Wi vi
i i
i i
Ω ∑
Ω ∑
Ω ∑
i
i
i
(11)

where ti , vi , εi , and Wi are atomic kinetic energy,
velocity, potential energy, and virial tensor for the
ith atom, respectively.
2.6.2 Atomic virials in heat ﬂux formula
To derive a rigorous formula of the heat ﬂux forthe many-body potentials, an expression for theatomic potential energy is required [8], which corresponds to Eq. (2) for the ANN potentials. In ourprevious study [5], we have succeeded in deriving(13)∂ εiANN ∂ ε j
ri j ⊗+∑
2 j̸=i∂ r ji∂ r ji1
WANN-III
=
i

ANN

]
. (14)

WANN-II
appears naturally in the stress-tensor
iderivation for the ANN potentials [25].
In
Ref. [26], WANN-III
wasusedbecause
the
pairi
[
]
wise force deﬁned as FANN=−ij

⟨JQ,α (t)JQ,β (0)⟩dt, (10)(12)

Incidentally, for the many-body potentials suchas the ANN potential, there are countless deﬁnitions of atomic virial which give identical total
virials [9, 5]. We deﬁne two more apparently reasonable atomic virials used for heat ﬂuxes in theprevious studies.

and

2.6 TC calculation∂ ε ANNj.
∂ r ji∂εj∂ εiANN∂ r ji + ∂ r ji

ANNsatisﬁes Newton’s third law [8]. The summation
Nis equal to FANN
. While the mathe∑ j atom FANN
ij
imatical properties of the three atomic virials differ, their sums over Natom yield the same total virial
WANN . Hereinafter, the heat ﬂux formulae deﬁned
by substituting WANN-I, WANN-II, and WANN-IIIiiiI
into Eq. (11) are referred to as JQ , JIIQ , and JIII
Q , respectively. Some studies used ANN potentials that
estimated TCs with JIIQ [27] and JIII
Q [26]. Recent
studies highlighted that the atomic virial can signiﬁcantly impact the TC when dealing with manybody potentials [8, 28, 29]. It is crucial to compare
the TCs obtained from these heat ﬂux formulae.
2.6.3

HNEMD method

The HNEMD method [30], which is based on the
GK formula, has successfully reduced the computational cost by enabling the calculation of TC using the time average of the heat ﬂux instead of integrating its auto-correlation function.

κ=

18

Ω

lim ⟨JQ,x ⟩t ,Fext T t→∞(15)

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

where Fext denotes a magnitude of perturbation
along the x direction of the system. Fext = 0.01
bohr−1 was selected for the HNEMD simulationsin this study [5].

43.1

Solid-liquid phase transition of Na

We ﬁrst investigated whether the EF-ﬁt and EFPﬁt models reproduced the hysteresis of melting andsolidiﬁcation in Na [4]. These two types of ANN
potentials were constructed from the FPMD data inthe temperature range of 200 to 450 K, performed
using 128 atoms with Γ point as mentioned in section 2.3. The EFP-ﬁt was able to accurately reproduce the heating and cooling processes of FPMD.
The temperature dependence of radial distributionfunction g(r) obtained from the MD simulations
is shown in Fig. 1. g(r)’s by the EFP-ﬁt (green
dashed lines) are quite consistent with the FPMD
results (black solid lines). In contrast, the g(r) obtained from the MD simulations using EF-ﬁt (reddashed lines) showed a signiﬁcant deviation fromthose obtained from FPMD . The above results
demonstrate the risk by using MLIPs such as EFﬁt, which do not guarantee the accuracy of stress.

3.2

Size dependence of the free energy ofNa

In many previous studies, the training data of
MLIPs were created by FP calculation with manyk points to avoid using large-scale systems. However, it is not obvious whether the respective
MLIPs constructed from the training data obtainedwith a speciﬁc number of k points and corresponding system sizes show equivalent accuracy. Therefore, we conducted the following analysis to clarify the inﬂuence of training conditions such as system size and the number of k points on thermodynamic quantities such as free energy [4].
We created the four types of EFP-ﬁt models.
These were trained with the FPMD data obtained
from the 128-atom system with Γ point, 128-

FPMD
EF-fitEFP-fit350 K

cooling2

g(r)Results

400 K

2
0
4

3

heating

0
4
300 K
2
0
4

250 K

2
0

2

4

6

2

4

6

8

r [Å]

Figure 1: Temperature dependence of the radial
distribution function g(r) in the heating and cooling processes. The black solid lines correspond to
the results from FPMD simulations. The red andgreen dashed lines correspond to those from MD
simulations with EF-ﬁt and EFP-ﬁt, respectively.
Taken from Ref. [4].

atom system with 4k points, 128-atom systemwith 32k points, and 1024-atom system with Γpoint, which are hereinafter referred to as 128ΓANN, 128-4kANN, 128-32kANN, and 1024ΓANN potentials, respectively. To investigate thesize-dependent behavior of thermodynamic quantities, we calculated the Helmholtz free energy
for 128 and 1024 atoms using 128-ΓANN and
128-4kANN potentials. Additionally, for comparison, we computed the same quantities for 1024
atoms using 1024-ΓANN and 128-32kANN potentials. Figure 2 shows the (a) Helmholtz free energy F, (b) total potential energy U, and (c) product of entropy and temperature T S as a functionof the number of atoms N. The triangles and circles represent the solid and liquid phases, respectively. The black, red, green, and blue symbols de-19

Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of TokyoF [eV/atom]

note the values calculated using the 1024-ΓANN,
128-4kANN, 128-ΓANN, and 128-32kANN potentials, respectively. As observed from this ﬁgure,the outcomes obtained using the128-ΓANN potential exhibit signiﬁcant deviation from the others.
Moreover, the values computed using the 1024ΓANN, 128-4kANN, and 128-32kANN potentials
are nearly identical with regard to the vertical axis.

-391.86(a)-391.88(b)

U [eV/atom]-391.68-391.70-391.72∆F between 128 and 1024 atoms is merely 0.13
and 0.11 meV/atom when calculated with the 128ΓANN and 128-4kANN, respectively, suggestingthat for systems with more than 128 atoms, the size
dependence on ∆F is almost negligible in comparison with the difference of Helmholtz free energy;(iii) the difference between the ∆F values computed using the128-4kANN and 1024-ΓANN potentials is relatively small, with a difference of only
0.50 meV/atom for 1024 atoms; (iv) the value of∆F obtained using the 128-32kANN is more similar to that obtained using the 1024-ΓANN than tothat obtained using the 128-4kANN. Despite being
trained with FPMD data at a single temperature of
350 K, the 128-4kANN potential is capable of accurately calculating ∆F with an error of only a few
tenths of meV/atom. However, if higher precision
is required, larger systems with more k points will
be necessary.

TS [eV/atom](c)1024-GANN
128-4kANN
128-GANN
128-32kANN

0.18

0.16
0

500

Table 1: Differences of Helmholtz free energy∆F = Fliq − Fsol between the solid and liquid statesat 350 K in meV/atom. Taken from Ref. [4]Atoms 128-Γ 128-4k 1024-Γ 128-32k
128
1024

1000

N

Figure 2: (a) Helmholtz free energy F, (b) total
potential energy U, and (c) product of entropy andtemperature T S obtained using the TI method at
350 K as a function of the number of atoms N.
Triangles and circles correspond to the solid andliquid states, respectively. The black, red, green,
and blue symbols indicate the values calculatedusing the 1024-ΓANN, 128-4kANN, 128-ΓANN,
and 128-32kANN potentials, respectively. Taken
from Ref. [4].
To further examine the size dependence, we
computed the differences in Helmholtz free energy between the solid and liquid states ∆F =Fliq − Fsol , as shown in Table 1. From the table,
we obtained the following ﬁndings: (i) the values obtained from the 128-ΓANN potential differs largely from the others; (ii) the difference in

3.3

-2.47
-2.34

-0.63
-0.74

-0.24

-0.33

TC calculation of α -Ag2 Se

The stress training is also associated with the accuracy of TC. Here, we present the results of aninvestigation of the inﬂuence of the heat ﬂux formulae (JIQ , JIIQ , and JIII
Q ) and stress training on TC,
using α -Ag2 Se as the test system [5].
To examine the effect of stress training, we used
the EF-ﬁt and EFP-ﬁt models, and also took into
account the dependence of the initial weight parameters in ANN. We constructed ﬁve each of
the EF-ﬁt and EFP-ﬁt with different initial weightparameters. Figure 3 plots ﬁve TCs (κ ANN ) obtained from the respective EF-ﬁt and EFP-ﬁt models as a function of root mean square errors ofstress for the validation data (∆Pvalid ). The vertical axis shows the absolute value |κ ANN − κ EIP |,20

4

Summary

In this report, we have investigated and discussed
issues that are unveriﬁed with the current training
method of MLIPs. Speciﬁcally, we reported thatthe stress training is necessary to accurately reproduce the solid-liquid phase transition of Na and theTC of α -Ag2 Se. In addition, it was suggested thatin order to construct more accurate MLIPs fromthe viewpoint of free energy calculations, increasing the number of k points is not the only solution,
and a large-scale system is also necessary.
We would like here to comment on other unestablished aspects of MLIP training methods. It
is important to recognize that adjusting the coefﬁcients of the cost function can dramatically alter
the accuracy of MLIPs [14]. The heat ﬂux regu-

-1
-1
EIP

0.1

EF-fit withJIEF-fit with JIIIII
EF-fit with J
I
EFP-fit with J
II
EFP-fit with J
EFP-fit with JIII


ANN

-k

0.2

|k

where κ EIP is the reference TC obtained from theEIP (0.274 Wm−1 K−1 ). For the EF-ﬁt, while three
TCs calculated by JIQ are sufﬁciently close to 0,
all of the ﬁve ∆Pvalid are larger than 0.47 GPa.
Furthermore, |κ ANN − κ EIP | exhibits a large valueof ∼ 0.2 Wm−1 K−1 with higher values of ∆Pvalid
when using JIQ . JIIQ gives large errors of the TCs∼0.25 Wm−1 K−1 with higher ∆Pvalid . In the caseof the TC calculated with JIII
Q , the deviation in the
TC is larger than 0.02 Wm−1 K−1 . The plots in
Fig. 3 indicate that when ∆Pvalid is larger, errors in
TC become larger and sometimes result in signiﬁcant deviations.
The good agreement with the reference values
suggests that stress training plays a crucial role as
anticipated. For EFP-ﬁt, all |κ ANN − κ EIP | calculated with JIQ is plotted near the origin as shownin Fig. 3. There are no cases showing larger TC
errors as seen in the EF-ﬁt. The maximum valueof |κ ANN − κ EIP | is ∼0.005 W m−1 K−1 , which
is much smaller than even the minimum valuesof |κ ANN − κ EIP | calculated with JIIQ and JIII
Q , i.e.
−1−10.055 and 0.041 W m K , respectively.
From the above results, only JIQ gives a correct
estimation of TC. Further, the correct TC may not
be obtained even with JIQ unless ∆Pvalid is sufﬁciently minimized through stress training.

| (Wm K )Activity Report 2022 / Supercomputer Center, Institute for Solid State Physics, The University of Tokyo

0.0
0.0

0.5

1.0

1.5

2.0

ΔPvalid (GPa)Figure 3: Absolute difference between the TC calculated from the ANN potentials (κ ANN ) and theEIP (κ EIP ) as a function of root mean square errors
for the validation data (∆Pvalid ) of ﬁve ANN potentials belonging EF-ﬁt (ﬁlled circles) and EFPﬁt (cross marks) for the validation data. κ ANN
and κ EIP are calculated from HNEMD simulations.
κ ANN are computed using heat ﬂux formulae JIQ(black), JIIQ (red), and JIII
Q (green). Taken from
Ref. [5].

larization method [31] introduced into the costfunction to estimate partial TCs may also lead to
improved robustness of MLIPs. It should be emphasized that more attention should be paid to theimportance of the modiﬁcation of the cost function
for the improvement of MLIPs as well as the developments in their architectures and descriptors.

Ackwnowledgements
This study was supported by MEXT/JSPS
KAKENHI Grant Numbers Nos. 22K03454,
21H01766, and 19K14676, and JST CREST Grant
Number JPMJCR18I2, Japan. The authors thankthe Supercomputer Center, the Institute for Solid
State Physics, University of Tokyo for the use of
the facilities. The computations were also carriedout using the facilities of the Research Institute forInformation Technology, Kyushu University.

21

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

References[1] Morawietz T and Artrith N 2021 J. Comput.
Aided Mol. Des. 35 557–586[2] Misawa M, Fukushima S, Koura A, Shimamura K, Shimojo F, Tiwari S, Nomura K i,
Kalia R K, Nakano A and Vashishta P 2020
J. Phys. Chem. Lett. 11 4536–4541[3] Yang M, Bonati L, Polino D and Parrinello M
2022 Catal. Today 387 143–149[4] Irie A, Fukushima S, Koura A, Shimamura
K and Shimojo F 2021 J. Phys. Soc. Jpn 90
094603[5] Shimamura K, Takeshita Y, Fukushima S,
Koura A and Shimojo F 2021 Chem. Phys.
Lett. 778 138748[15] Kresse G and Joubert D 1999 Phys. Rev. B 59
1758–1775[16] Hohenberg P and Kohn W 1964 Phys. Rev.
136 B864[17] Perdew J P, Burke K and Ernzerhof M 1996
Phys. Rev. Lett. 77 3865–3868[18] Tuckerman M, Berne B J and Martyna G J
1992 J. Chem. Phys. 97 1990–2001[19] Rang M and Kresse G 2019 Phys. Rev. B
99(18) 184103[20] Rino J P, Hornos Y M M, Antonio G A,
Ebbsjö I, Kalia R K and Vashishta P 1988 J.
Chem. Phys. 89 7542–7555[21] Shimojo F and Okazaki H 1991 J. Phys. Soc.
Jpn. 60 3745–3753[6] Behler J 2021 Chem. Rev. 121 10037–10072[7] Fukushima S, Ushijima E, Kumazoe H,
Koura A, Shimojo F, Shimamura K, Misawa
M, Kalia R K, Nakano A and Vashishta P2019 Phys. Rev. B 100(21) 214108[8] Fan Z, Pereira L F C, Wang H Q, Zheng J C,
Donadio D and Harju A 2015 Phys. Rev. B
92(9) 094301[22] Okazaki H 1967 J. Phys. Soc. Jpn. 23 355–
360[23] Kubo R 1957 J. Phys. Soc. Jpn. 12 570–586[24] Green M S 1952 J. Chem. Phys. 20 1281–
1295[25] Behler J 2011 J. Chem. Phys. 134 074106[9] Shimamura K, Takeshita Y, Fukushima S,
Koura A and Shimojo F 2020 J. Chem. Phys.
153 234301[26] Mangold C, Chen S, Barbalinardo G, Behler
J, Pochet P, Termentzidis K, Han Y, Chaput
L, Lacroix D and Donadio D 2020 J. Appl.
Phys. 127 244901[10] Shimojo F, Fukushima S, Kumazoe H, Misawa M, Ohmura S, Rajak P, Shimamura K,
Bassman L, Tiwari S, Kalia R K, Nakano A
and Vashishta P 2019 SoftwareX 10 100307[27] Sosso G C, Donadio D, Caravati S, Behler J
and Bernasconi M 2012 Phys. Rev. B 86(10)104301[11] Artrith N and Urban A 2016 Comput. Mater.
Sci. 114 135 – 150[12] Behler J and Parrinello M 2007 Phys. Rev.
Lett. 98(14) 146401[13] Artrith N, Urban A and Ceder G 2017 Phys.
Rev. B 96(1) 014112[14] Irie A, Shimamura K, Koura A and Shimojo
F 2022 J. Phys. Soc. Jpn. 91 045002[28] Surblys D, Matsubara H, Kikugawa G andOhara T 2019 Phys. Rev. E 99(5) 051301[29] Boone P, Babaei H and Wilmer C E 2019 J.
Chem. Theory Comput. 15 5579–5587[30] Evans D J 1982 Phys. Lett. A 91 457 – 460[31] Shimamura K, Koura A and Shimojo F 2022
arXiv: 2101.10468

22