Ab initio Canonical Sampling based on Variational Inference Alo ıs Castellano1Fran cois Bottin1Johann Bouchet2Antoine Levitt3and Gabriel Stoltz3 1CEA DAM DIF F-91297 Arpajon France and Universit e Paris-Saclay CEA

2025-04-30 0 0 4.44MB 20 页 10玖币

侵权投诉

Ab initio Canonical Sampling based on Variational Inference

Alo¨ıs Castellano,1Fran¸cois Bottin,1Johann Bouchet,2Antoine Levitt,3and Gabriel Stoltz3

1CEA, DAM, DIF, F-91297 Arpajon, France, and Universit´e Paris-Saclay, CEA,

Laboratoires des Mat´eriaux en Conditions Extrˆemes, 91680 Bruy`eres-le-Chˆatel, France.

2CEA, DES, IRESNE, DEC, Cadarache, F-13018 St Paul Les Durance, France.

3CERMICS, Ecole des Ponts, Marne-la-Vall´ee, France

MATHERIALS team-project, Inria Paris, France

Finite temperature calculations, based on ab initio molecular dynamics (AIMD) simulations, are

a powerful tool able to predict material properties that cannot be deduced from ground state calcu-

lations. However, the high computational cost of AIMD limits its applicability for large or complex

systems. To circumvent this limitation we introduce a new method named Machine Learning As-

sisted Canonical Sampling (MLACS), which accelerates the sampling of the Born–Oppenheimer

potential surface in the canonical ensemble. Based on a self-consistent variational procedure, the

method iteratively trains a Machine Learning Interatomic Potential to generate conﬁgurations that

approximate the canonical distribution of positions associated with the ab initio potential energy.

By proving the reliability of the method on anharmonic systems, we show that the method is able

to reproduce the results of AIMD with an ab initio accuracy at a fraction of its computational cost.

Molecular dynamics (MD) simulations, beside Monte-

Carlo calculations [1], are nowadays a popular way to ob-

tain ﬁnite temperature properties and explore the phase

diagram of materials. If the seminal works of Alder [2–

4], Rahman [5,6] and coworkers only treated classical

potentials, all the powerful tools introduced in these pa-

pers were reused in ab initio molecular dynamics (AIMD)

codes thereafter. In particular, by using the density func-

tional theory (DFT) [7,8] and the Born–Oppenheimer

(BO) [9] approximation, but also by assuming that the

ground state electronic density is obtained at each MD

time step [10,11], all the ideas proposed previously apply.

However, the high computational cost, due to both the

evaluation of the Hellmann–Feynman forces using DFT

and the high number of MD time steps required to sam-

ple the BO surface potential in the canonical ensemble,

limits the range of systems that can be studied. To accel-

erate the computation of ﬁnite temperature properties,

two main strategies have been proposed. In the ﬁrst one,

AIMD simulations are replaced by MD with ab initio-

based numerical potentials. In the second one, the sam-

pling of the canonical distribution is performed through

a direct generation of atomic conﬁgurations, bypassing

AIMD simulations.

Recent progress in the ﬁeld of Machine Learning In-

teratomic Potentials (MLIP) [12–16] promises an accel-

eration of ﬁnite-temperature studies with a near-DFT

accuracy. This high accuracy is ensured by the ﬂexibility

of MLIPs, which enables to reproduce a large variety of

BO surfaces. However, the construction of a MLIP is a

tedious task as MLIP show poor extrapolative capabili-

ties. Consequently, the set of conﬁgurations used for the

MLIP requires a careful construction, which led to the de-

velopment of learn-on-the-ﬂy molecular dynamics [17–19]

and other active learning based dataset selection meth-

ods [20–22]. Moreover, usual ﬁnite-temperature works

using MLIP shift the studied system from the DFT de-

scription to the MLIP one, with atomic positions (and

computed properties) distributed according to the Boltz-

mann weights associated with the MLIP potential. Con-

sequently, one may need to ensure that simulations using

MLIP are not in the extrapolative regime. This veriﬁca-

tion can be done by computing corrections from free en-

ergy perturbation methods [23], which however requires

additional DFT single-point calculations.

Several groups have recently proposed methods to by-

pass AIMD and generate conﬁgurations by adapting the

Self-Consistent Harmonic Approximation (SCHA) [24–

27] to modern ab initio methods, using a varational in-

ference strategy. Among those, we can cite the stochas-

tic Temperature Dependent Eﬀective Potential [28],

the Stochastic SCHA [29–31] and the Quantum Self-

Consistent Ab Initio Lattice Dynamics [32]. Within

those methods (named EHCS in the following, for Eﬀec-

tive Harmonic Canonical Sampling), conﬁgurations are

generated with displacements around equilibrium posi-

tions according to a distribution corresponding to an ef-

fective harmonic Hamiltonian. The self-consistent (SC)

construction of this Hamiltonian, based on a variational

procedure, allows to include explicit temperature eﬀects.

However, the harmonic form of the eﬀective potential

means that the atomic displacements follow a Gaussian

distribution. This can entice diﬀerences for actual dis-

tributions of displacements in highly anharmonic solids

or close to the melting temperature, and makes this ap-

proach completely inapplicable on liquids.

In this Letter, we propose a new method named Ma-

chine Learning Assisted Canonical Sampling (MLACS),

which can be thought of as a generalization of the vari-

ational inference strategy used in the EHCS methods to

linear MLIP. MLACS consists in a SC variational pro-

cedure to generate conﬁgurations in order to best ap-

proximate the DFT canonical distribution and obtain

a near-DFT accuracy in the properties computed at a

arXiv:2210.11531v1 [cond-mat.mtrl-sci] 20 Oct 2022

fraction of the cost of AIMD. In the present work, the

MLIP is built to produce an eﬀective canonical distri-

bution which reproduces the DFT one at a single ther-

modynamical point, rather than sketching an eﬀective

BO potential surface for all thermodynamic conditions.

Therefore, MLACS is a sampling method which accel-

erates the computation of ﬁnite-temperature properties

compared to AIMD simulations.

Let us consider an arbitrary system of Nat atoms (a

crystal or a liquid, elemental or alloyed) at a tempera-

ture T, described by the coordinates R∈R3Nat . The

potential energy V(R) induces the canonical equilibrium

distribution p(R) = e−βV (R)/Z, with β= 1/kBT and

Z=RdRe−βV (R)the partition function. For this sys-

tem, the average of an observable O(R), depending only

on positions, is given by

hOi=ZdRO(R)p(R)≈X

O(Rn)wn, (1)

with wnthe weight of a conﬁguration n, which follows

the normalization Pnwn= 1. Rather than performing

AIMD simulations to generate a set of Rnand wn, the

objective of MLACS is to generate a reduced set of con-

ﬁgurations and weights using a MLIP and to compute

them using DFT.

Let us deﬁne to this purpose a MLIP potential

Vγ

γ(R) = PK

k=1 e

Dk(R)γkwith partition function e

Zγ

γ=

RdRe−β

Vγ

γ(R)and canonical equilibrium distribution

eqγ

γ(R) = e−β

Vγ

γ(R)/e

Zγ

γ. The column vector γ

γ∈RKcon-

tains adjustable parameters and the functions included

in the row vector e

D:R3Nat →RKare named descrip-

tors. While a linear dependence is assumed between the

potential energy and the descriptor space, no assumption

is made on the form of e

D(R), enabling the use of a wide

variety of descriptors, from a harmonic form to more so-

phisticated atom centered descriptors such as Atom Cen-

tered Symmetry Functions [33] or the Smooth Overlap of

Atomic Positions [14,34]. For the sake of simplicity and

robustness, we use the Spectral Neighbor Analysis Po-

tential [15,35,36] in this work.

A key idea in MLACS is to use the distribution eqγ

γ(R)

instead of p(R) in Eq. (1). With this approach, the ac-

celeration is provided by the reduced computational cost

to generate conﬁgurations using e

Vγ

γ(R) instead of V(R).

To ensure accurate results, one needs to adjust the pa-

rameters γ

γof the eﬀective potential so that eqγ

γ(R) is the

best approximation to the true distribution p(R). As a

measure of the similarity between the two distributions

pand eqγ

γ, we use the Kullback–Leibler divergence (KLD)

DKL(eqγ

γkp) = ZdReqγ

γ(R) ln eqγ

γ(R)

p(R)≥0 . (2)

The KLD is a non-negative and asymmetric measure of

the discrepancy between two distributions. The smaller

the KLD, the closer the distribution pand eqγ

γare, with

DKL(eqγ

γkp) = 0 meaning that the two distributions are

identical. So, by minimizing this quantity with respect to

the parameters γ

γ, we obtain the distribution eqγ

γ(R) that

best reproduces p(R). A more tractable formulation can

be obtain by transforming Eq. (2) to an equivalent free

energy minimisation (see SM I):

F ≤ min

γe

Fγ

γwith e

Fγ

γ=e

γ+hV(R)−e

Vγ

γ(R)ie

Vγ

γ, (3)

where hOie

Vγ

γ=RdRO(R)eqγ

γ(R) is the canonical aver-

age for the potential e

Vγ

γ(R), F=−kBT ln(Z) the free

energy associated to V(R), and e

γ=−kBT ln( e

Zγ

γ) the

free energy associated to e

Vγ

γ(R). The Gibbs–Bogoliubov

(GB) free energy e

Fγ

γwritten in Eq. (3) is the starting

point for various variational procedures, in particular the

SCHA [24,26,30]. By minimizing Eq. (3) with respect

to γ

γ, we obtain (see SM II)

γ=De

D(R)Te

D(R)E−1

Vγ

γDe

D(R)TV(R)Ee

Vγ

, (4)

a non-trivial least-squares solution showing a circular de-

pendency over γ

γ, solved using a SC variational procedure

described below. Thus, the initial problem turns into an

optimisation problem, solved with a ﬂexible and simple

MLIP potential, using variational inference [37,38]. So

far, only the DFT data V(R) are used, and not their

derivatives. In order to lower the number of DFT calcu-

lations needed we also use the Fisher divergence [38,39]

DF(eqγ

γkp) = RdReqγ

γ(R)|∇Rln[eqγ

γ(R)/p(R)]|2≥0, which

measures the diﬀerence between two distributions us-

ing information on the gradients. The Fisher diver-

gence can be rewritten (see SM III) as DF(eqγ

γkp) =

β2h|F−e

Fγ

γ|2ie

Vγ

γ≥0, with F=−∇RV(R) and e

Fγ

γ=

−∇Re

Vγ

γ(R), the DFT and MLIP forces, respectively. Fi-

nally, the optimal parameters of the MLIP are obtained

by minimizing the cost function ∆ = αKLDKL +αFDF

with respect to γ

γ, with αKL and αFthe weights of each

contribution.

In order to perform this minimization, we adopt the

following SC procedure (see Fig. 1): start with N

atomic conﬁgurations, compute their energies, forces and

stresses at the DFT level, check the convergence of

speciﬁc observables (phonon spectrum, pair distribution

function...), compute γ

γ, build a MLIP potential e

Vγ

γ, per-

form a MD run using the MLIP (i.e. sample using e

Vγ

γ)

and extract new Natomic conﬁgurations from the tra-

jectory. This SC procedure is repeated until convergence.

Formally, the least-squares ﬁtting has to be made using

data distributed according to the last e

Vγ

γ. To still reduce

the number of DFT calculations, all the conﬁgurations

N conﬁgurations

Density Functional

Theory

Weights & Properties

Computation

Machine-Learning

Molecular Dynamics

Machine-Learning

Interatomic Potential

Training

MLACS

Temperature dependent properties

Convergence

Start here

FIG. 1. Workﬂow of MLACS. All the computational details

are given in SM IV.

from previous iterations are reused and reweighted (see

SM IV) using the Multistate Bennett Acceptance Ratio

(MBAR) [40,41]. The weights wnobtained by this pro-

cedure are used to ﬁt the MLIP potential and to evaluate

physical properties (see Eq. (1)). Therefore, all observ-

ables computed using DFT (forces, energies and stresses,

but also electronic properties) can be averaged at no ex-

tra cost, in contrast with usual MLIPs. Moreover, the

free energy of the DFT system can be computed using

both thermodynamic perturbation [42,43] and cumulant

expansions [44], without extra calculations (see SM V).

Concerning the computational cost, the acceleration

enabled by MLACS with respect to AIMD calculations

is twofold. First, MLACS strongly reduces the time re-

quired to generate independent conﬁgurations (around

a hundred DFT calculations are needed and no longer

thousands as in AIMD). Secondly, the calculation of the

Nconﬁgurations can be parallelized at each iteration.

Thus, an acceleration factor of 50 in computation time

and 1000 in human time can be reached with N= 20.

Moreover, an additional acceleration can be obtained by

using, as a starting point of MLACS, a MLIP potential

coming from a previous calculation.

To demonstrate the accuracy and versatility of the

method, we compare AIMD, MLACS and EHCS results

on eight diﬀerent systems (crystals, liquid and alloy) by

performing three classical simulations (see SM VII) and

ﬁve DFT calculations (see SM VIII). The ab initio cal-

culations are performed using the ABINIT code [45] over

thousands of processors [46]. As a stringent test, we fo-

cus on the phonon spectrum, because it is particularly

sensitive to the quality of the sampling so that its con-

vergence guarantees the convergence of thermodynamical

and elastic properties. These data are extracted using

the Temperature Dependent Eﬀective Potential (TDEP)

method [47] as implemented in ABINIT [48].

As a ﬁrst DFT example, we consider the diamond

phase of Si with a (3 ×3×3) supercell including 216

atoms at T = 900 K. It has been shown using EHCS that

anharmonicity is manifest at this temperature [49,50],

FIG. 2. Results for Si at T = 900 K. Main: TDEP

phonon spectra for AIMD (blue), MLACS (red) and EHCS

(green) simulations. Inset: weight of each conﬁguration for

the MLACS simulation.

whereas Si is more harmonic at lower temperatures [51].

Even though MLACS only requires Nconfs = 60 conﬁgu-

rations (20 for the initialization and 40 for the produc-

tion), the phonon spectrum of Si (see Fig. 2) is in excel-

lent agreement with the AIMD one (Nconfs = 3546 during

9 ps). In contrast, the EHCS simulations (Nconfs = 160)

overestimate the optical branches by several meV. This

fast and good convergence comes from both the ability

of MLACS to generate uncorrelated conﬁgurations and

the eﬃciency of the reweighting. Consequently, a near-

DFT accuracy associated with a strong acceleration of

the computational time, from one week with AIMD to 3

hours using MLACS, can be achieved.

We next consider Zr in its bcc βphase, with a (4 ×

4×4) supercell containing 128 atoms, at T = 1000 K.

The β-Zr phase being unstable at low temperatures, the

stabilisation at higher temperatures comes from strong

anharmonic eﬀects [47,48,52]. As previously, EHCS

overestimates the whole phonon spectrum compared to

AIMD (see Fig. 3). Conversely, the phonon spectrum

computed using MLACS is in excellent agreement with

the AIMD one. Here, the acceleration in computational

time is from two months with AIMD (with Nconfs = 7758

during 19 ps) to two days using MLACS (Nconfs = 160).

The third DFT example is the bcc γphase of Ura-

FIG. 3. Same caption as Fig. 2for β-Zr at T = 1000 K.

FIG. 4. Correlations between AIMD and MLACS forces: a) Si at T = 900 K, b) β-Zr at T = 1000 K and c) γ-U at T = 1200 K.

nium with a (4 ×4×4) supercell containing 128 atoms,

at T = 1200 K. As for β-Zr, γ-U is unstable at low tem-

perature and is stabilized by anharmonic eﬀects [53–55].

The EHCS calculations never converged, despite several

attempts and starting points (especially using the MLIP

ﬁtted using AIMD simulations). The phonon spectrum

obtained using MLACS (see Fig. 5) reproduces correctly

the AIMD one, even if some small discrepancies remain

(lower than 1 meV). Here, the reweighting is crucial,

since the ﬁrst 40 atomic conﬁgurations depart from the

bcc phase (they look like a glass) and are discarded by

MBAR from the statistical averages (see Fig. 4c and the

inset in Fig. 5). After them, the next 80 ones get closer

to the ﬁnal structure and have a non-zero weight in the

equilibrium canonical distribution. This diﬀers strongly

from Si (see Fig. 4a and the inset in Fig. 2) and β-Zr

(see Fig. 4b and the inset in Fig. 3), for which all the

conﬁgurations almost equally contribute. Despite this

loss of data, MLACS still leads to a signiﬁcant reduction

in the computational cost from two months with AIMD

(Nconfs = 5981 during 22 ps) to two days using MLACS

(Nconfs = 140).

We also present three systems simulated with classi-

cal potentials and with suﬃciently long MD trajectories

FIG. 5. Same caption as Fig. 2for γ-U at T = 1200 K.

FIG. 6. Phonon density of states of Al0.5Cu0.5at T = 600 K

obtained using MLACS and MD simulations.

FIG. 7. Pair distribution function of U liquid at T = 2500 K

obtained using MLACS and MD simulations.

(including tens of thousands time steps), which would

be impossible to simulate using DFT. The ﬁrst one is

Silicon at T = 1500 K using a Tersoﬀ potential, the sec-

ond one is the Al0.5Cu0.5alloy at T = 600 K using an

Angular-Dependent Potential (ADP), and the third one

is the liquid phase of Uranium at T = 2500 K, using a

Modiﬁed Embedded-Atom Method (MEAM) potential.

In Fig. 6and 7, we show that both the phonon density

of states (DOS) of Al0.5Cu0.5and the pair distribution

文档加载中……请稍候！
如果长时间未打开，您也可以点击刷新试试。

下载文档到电脑，查找使用更方便

10 玖币 0人已下载

立即下载

摘要：

AbinitioCanonicalSamplingbasedonVariationalInferenceAlosCastellano,1FrancoisBottin,1JohannBouchet,2AntoineLevitt,3andGabrielStoltz31CEA,DAM,DIF,F-91297Arpajon,France,andUniversiteParis-Saclay,CEA,LaboratoiresdesMateriauxenConditionsExtr^emes,91680Bruyeres-le-Ch^atel,France.2CEA,DES,IRESNE,DEC,...

展开>> 收起<<

Ab initio Canonical Sampling based on Variational Inference Alo ıs Castellano1Fran cois Bottin1Johann Bouchet2Antoine Levitt3and Gabriel Stoltz3 1CEA DAM DIF F-91297 Arpajon France and Universit e Paris-Saclay CEA.pdf

共20页,预览4页

还剩页未读，继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

Ab initio Canonical Sampling based on Variational Inference Alo ıs Castellano1Fran cois Bottin1Johann Bouchet2Antoine Levitt3and Gabriel Stoltz3 1CEA DAM DIF F-91297 Arpajon France and Universit e Paris-Saclay CEA

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: