On the algorithm to perform Monte Carlo simulations in cells with constant volume and variable shape

2025-05-02 0 0 1.39MB 18 页 10玖币

侵权投诉

Condensed Matter Physics, 2022, Vol. 25, No. 3, 33201: 1–18

DOI: 10.5488/CMP.25.33201

http://www.icmp.lviv.ua/journal

On the algorithm to perform Monte Carlo simulations

in cells with constant volume and variable shape

A. Baumketner ∗

Institute for Condensed Matter Physics of the National Academy of Sciences of Ukraine, 1 Svientsitskii Str.,

79011, Lviv, Ukraine

Received May 02, 2022, in ﬁnal form May 02, 2022

In simulations of crystals, unlike liquids or gases, it may happen that the properties of the studied system depend

not only on the volume of the simulation cell but also on its shape. For such cases it is desirable to change

the shape of the box on the ﬂy in the course of the simulation as it may not be known ahead of time which

geometry ﬁts the studied system best. In this work we derive an algorithm for this task based on the condition

that the distribution of speciﬁc geometrical parameter observed in simulations at a constant volume matches

that observed in the constant-pressure ensemble. The proposed algorithm is tested for the system of hard-

core ellipses which makes lattices of different types depending on the asphericity parameter of the particle.

It is shown that the performance of the algorithm critically depends on the range of the sampled geometrical

parameter. If the range is narrow, the impact of the sampling method is minimal. If the range is large, inadequate

sampling can lead to signiﬁcant distortions of the relevant distribution functions and, as a consequence, errors

in the estimates of free energy.

Key words: hard-ellipse ﬂuid, Monte Carlo simulation, constant volume, varying shape, umbrella sampling

1. Introduction

Today computer simulations play a key role in fundamental research across multiple disciplines,

including physics, chemistry and materials science [1]. A large share of computational studies employ

simulation boxes with ﬁxed volume. This choice is mainly motivated by convenience as constant-

volume/constant-temperature ensemble is easier to program than the equivalent ensemble with constant

pressure. But this is also due to the involvement of the constant-volume ensemble in other, specialized

simulation techniques such as free energy calculations [2], Gibbs ensemble [3] or replica-exchange

method [4–6]. Regardless of the particular context, it is always understood that the eﬀect of volume

vanishes in the thermodynamic limit where the results are thought to be independent of the employed

ensemble. This claim is certainly true for liquids or gases, whose properties are independent of the

geometry of the box.

In the case of crystals, however, the situation could be quite diﬀerent [7, 8]. In crystalline materials

there could be properties that depend explicitly on the volume as well as on the shape of the box. Take

for instance the example of a rectangular lattice with lattice constants 𝑎and 𝑏, as shown in ﬁgure 1. The

dimensions of the box that accommodates 𝑛columns and 𝑚rows are 𝐿𝑥=𝑎𝑛 along 𝑥axis and 𝐿𝑦=𝑏𝑚

along 𝑦axis as shown in ﬁgure 1 (a). The corresponding aspect ratio is 𝜏=𝐿𝑥/𝐿𝑦=𝑎𝑛/𝑏𝑚. Now, let

us assume that the lattice constants are not known ahead of time but are meant to be determined in the

course of the simulations. If we initially choose a box with the wrong aspect ratio, 𝜏0> 𝜏 for instance,

see ﬁgure 1 (b) for appropriate illustration, the lattice constant determined in simulations will also be

incorrect. Since the geometry of the cell drives the structure of the lattice, wrong geometry translates

into wrong structure. Importantly, lattice distortions will not go away easily even when the size of the

simulation box is increased.

∗Corresponding author: andrĳ@icmp.lviv.ua.

This work is licensed under a Creative Commons Attribution 4.0 International License. Further distribution

of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.

33201-1

arXiv:2210.00485v1 [physics.comp-ph] 2 Oct 2022

A. Baumketner

Figure 1. Aspect ratio of simulation boxes is governed by the geometry of the modelled lattice. Incorrect

initial guess will result in lattice distortions in constant volume simulations, for instance leading to 𝑎0> 𝑎.

To circumvent this problem, the box should be capable of changing its shape.

One way to deal with this issue is to estimate the free energy of the studied system 𝐹as a function

of some geometrical parameter of the cell, for instance 𝜏, and then choose the geometry with the

lowest 𝐹. This can be done by a variety of tools, including the Einstein crystal method for the free

energy computation [2, 9]. Although formally correct, this approach is cumbersome and carries a large

computational cost. An alternative is to allow the shape of the box to change in the course of the

simulation. The proper geometry then corresponds to the free energy minimum, so it will be seen as the

most frequently visited structure. An additional beneﬁt will be for systems that can populate multiple

geometries at the same time as their relative free energy in this case can be determined from a single

Monte Carlo (MC) trajectory.

The question then is how does one allow the shape of the box to change? How is that accomplished

in practice? Would, for instance, generating randomly, from time to time, a new aspect ratio 𝜏in the

course of the simulation constitute a good method? If not, what is the good method? These questions had

been addressed before as multiple studies report using simulation cells with constant volume but variable

shape (MCVS) [7, 8, 10, 11]. Unfortunately, the details of the performed simulations are scarce and to

establish the speciﬁcs of the used algorithms appears diﬃcult. Yet, we ﬁnd evidence that the method

one employs to sample the trial geometries may have measurable consequences for physical properties

extracted from simulations. Thus, which geometry parameters are best to choose in MCVS simulations,

and how to choose them remains unclear. These are the questions that we answer in the present article.

We show that in order for the constant-volume ensemble with variable shape to be consistent with the

constant-pressure ensemble, the aspect ratio should be sampled from the 1/𝜏distribution. Any other

sampling law will lead to erroneous results. We illustrate this point for the system of impenetrable

ellipses in two-dimensional space, for which we evaluate the performance of the method that relies on

uniformly sampled 𝜏, or on the so-called 𝜏-sampling and show that it produces wrong free energy for the

relevant states of the system.

2. Theory

2.1. Sampling law for 𝝉

Let us assume that, in addition to volume, the partition function in the canonical, or NVT for short,

ensemble 𝑄(𝑁, 𝑇, 𝑉 ;𝜏)=𝑉;𝜏exp[−𝛽U(Γ)]dΓexplicitly depends on some geometrical parameter 𝜏,

for instance the aspect ratio of the box sides. Here, 𝛽is the inverse temperature, 𝑉is the volume, 𝑈(Γ)is

the potential energy and Γis the abbreviation for the point in the conﬁguration space. The integration is

carried out over volume 𝑉with the set parameter 𝜏. The full partition function then should be constructed

as a weighted sum (or integral) over all possible realizations of the additional degree of freedom [12]. In

33201-2

MC simulations in boxes with variable shape

the most general case, one ﬁnds:

𝑄(𝑁, 𝑇, 𝑉 )=𝑄(𝑁, 𝑇, 𝑉;𝜏)𝑓(𝜏)d𝜏=𝑃(𝜏)d𝜏, (2.1)

where 𝑓(𝜏)is the weighting function of the extended ensemble deﬁned by both volume and the shape

of the box and 𝑃(𝜏)is the probability distribution function of 𝜏, characteristic of the NVT ensemble

with variable shape. In principle, one is free to choose 𝑓(𝜏)arbitrarily, provided that it satisﬁes certain

conditions typically imposed on distribution functions, such as positive deﬁniteness or integrability. Here,

we will select 𝑓(𝜏)on the condition that the extended ensemble with ﬁxed volume satisﬁes distribution

of 𝜏speciﬁc for the constant pressure, or NPT, ensemble. In this way, modelling in the two ensembles,

constant-volume and constant-pressure, will be consistent, hence minimizing ﬁnite size eﬀects.

The distribution function in the constant-pressure ensemble reads:

𝑃(Γ, 𝑁, 𝑇, P) ∼ e−𝛽P𝑉e−𝛽𝑈 (Γ),(2.2)

where Pis the pressure. Let us focus now on 2D space and obtain formulas for this simpler case ﬁrst.

The relevant volume is 𝑉=𝐿𝑥𝐿𝑦sin 𝛼, where 𝐿𝑥and 𝐿𝑦are the lengths of the simulation box and sin 𝛼

is the sine of the angle between them. In NPT simulations, volume is sampled randomly from a uniform

distribution. This can be achieved either by uniformly sampling one of the variables involved in volume,

𝐿𝑥,𝐿𝑦or sin 𝛼, or by non-uniformly sampling some combination of these variables which leads to a

uniformly distributed volume [2]. Let us assume that the ﬁrst scenario takes place, as it is more general

and can be applied to both liquids and crystals, and each concerned variable is sampled uniformly, one

at a time and in random order. The joint distribution function measured in such simulations for variables

𝐿𝑥,𝐿𝑦and sin 𝛼will be given by the following expression:

𝑃𝐿𝑥, 𝐿𝑦,sin 𝛼∼e−𝛽P𝐿𝑥𝐿𝑦sin 𝛼𝑄(𝑁, 𝑇, 𝑉;𝜏),(2.3)

which can also be obtained by integrating distribution function (2.2) over all conﬁgurations Γ. The

dependence on 𝜏in the partition function arises because the integration over Γis carried out for the

box with speciﬁc dimensions 𝐿𝑥and 𝐿𝑦, which in addition to the volume 𝑉=𝐿𝑥𝐿𝑦also deﬁne other

geometrical parameters including 𝜏=𝐿𝑥/𝐿𝑦.

Figure 2. Cartoon illustrating hypothetical joint distribution function of variables 𝐿𝑥and 𝐿𝑦for the

constant-pressure ensemble. Constant-volume conﬁgurations correspond to the sub-ensemble bound to

the line 𝐿𝑦=𝑉/𝐿𝑥, where 𝐿𝑥is treated as the independent variable. Sampling along this line uniquely

deﬁnes the distribution function 𝑃𝑉(𝜏)for geometrical parameter 𝜏.

As an illustration, consider a schematic distribution 𝑃(𝐿𝑥, 𝐿𝑦)shown in ﬁgure 2, speciﬁc for rectan-

gular boxes with sin 𝛼=1. Among all possible conﬁgurations, those that correspond to volume 𝑉satisfy

the constraint 𝑉=𝐿𝑥𝐿𝑦, creating a one-dimensional sub-ensemble characterized by a single degree of

freedom. If 𝐿𝑥is chosen as the independent degree of freedom, conﬁgurations with given 𝑉and 𝐿𝑥will

appear in simulations with probability 𝑃(𝐿𝑥, 𝑉/𝐿𝑥). Any other geometric parameter of the system will

33201-3

A. Baumketner

also be characterized by a unique distribution function. This includes the aspect ratio 𝜏=𝐿𝑥/𝐿𝑦, for

which the distribution function 𝑃𝑉(𝜏)represents the relative frequency of seeing conformations with

diﬀerent 𝜏. This function is directly accessible in simulations.

Our goal is to evaluate 𝑃𝑉(𝜏)and then try to reproduce it in simulations at a constant volume. Using

𝑃(𝐿𝑥, 𝐿𝑦,sin 𝛼)as the starting point, 𝑃𝑉(𝜏)can be evaluated as 𝑃𝑉(𝜏)=𝛿𝜏−𝐿𝑥/𝐿𝑦𝑉, where

the brackets h···i𝑉indicate a statistical average over sub-ensemble with ﬁxed 𝑉. The constant-volume

constraint can be imposed through another delta function, leading to the following expression in terms

of the unbiased distribution function:

𝑃𝑉(𝜏) ∼ d𝐿𝑥d𝐿𝑦𝛿𝑉−𝐿𝑥𝐿𝑦sin 𝛼𝛿(𝜏−𝐿𝑥/𝐿𝑦)e−𝛽P𝐿𝑥𝐿𝑦sin 𝛼𝑄(𝑁, 𝑇, 𝐿𝑥𝐿𝑦sin 𝛼;𝐿𝑥/𝐿𝑦).

(2.4)

This integral can be evaluated in two steps. First, a change of variables 𝑦=𝐿𝑥𝐿𝑦sin 𝛼helps to eliminate

the integration over 𝐿𝑦while keeping the other variable constant. The result is:

𝑃𝑉(𝜏) ∼ d𝐿𝑥𝛿𝜏−𝐿2

𝑥sin 𝛼

𝑉1

𝐿𝑥sin 𝛼e−𝛽P𝑉𝑄(𝑁, 𝑇, 𝑉 ;𝐿2

𝑥sin 𝛼/𝑉).(2.5)

The second integration can be carried out by making a change of variables 𝑦=𝐿2

𝑥sin 𝛼/𝑉, which after

dropping the terms that are independent of 𝜏yields:

𝑃𝑉(𝜏)=1

𝜏𝑄(𝑁, 𝑇, 𝑉 ;𝜏).(2.6)

Let us now require that 𝑃𝑉(𝜏)=𝑃(𝜏), i.e., the distribution functions in the NPT and NVT ensembles are

equal. It is easy to see then from equation (2.1) that 𝑄(𝑁, 𝑇 , 𝑉;𝜏)𝑓(𝜏)d𝜏=𝑃(𝜏)d𝜏=𝑃𝑉(𝜏)d𝜏=

𝑄(𝑁, 𝑇, 𝑉 ;𝜏)(1/𝜏)d𝜏, indicating that 𝑓(𝜏)=1/𝜏. In other words, we arrive at the conclusion that

new aspect ratios in constant-volume simulations should be drawn from the 1/𝜏distribution in order

to reproduce the result of the NPT simulations. Furthermore, it is seen from equation (2.6) that free

energy 𝛽𝐹 (𝜏)=−log [𝑄(𝑁, 𝑇, 𝑉 ;𝜏)]associated with the degree of freedom 𝜏can be computed as

𝛽𝐹 (𝜏)=−log [𝜏𝑃𝑉(𝜏)]. Therefore, the aspect ratio 𝜏is not suitable for the computation of free

energy diﬀerences directly from the distribution function. In other words, 𝛽Δ𝐹=𝛽[𝐹(𝜏1)−𝐹(𝜏2)] ≠

log [𝑃𝑉(𝜏2)/𝑃𝑉(𝜏1)], where 𝜏1and 𝜏2are some values deﬁning two macroscopic states. It is easy to

show, however, that 𝑃𝑉(𝑧)=𝜏(𝑧)𝑃𝑉(𝜏(𝑧)) is the distribution function for a new variable 𝑧=log(𝜏).

Thus, in terms of this variable 𝛽𝐹 [𝜏(𝑧)]=𝛽𝐹 (𝑧)=−log[𝜏(𝑧)𝑃𝑉(𝜏(𝑧))]=−log[𝑃𝑉(𝑧)], making 𝑧

the proper order parameter associated with the geometry parameter 𝜏. In the Appendix we show that 1/𝜏

sampling law also applies in the three-dimensional space.

2.2. Simulation algorithm

Given that the sampling law is known, how does one conduct a constant-volume simulation with

the variable shape? Let us ﬁrst point out the following auxiliary results. The lengths of the box can be

expressed in terms of 𝜏and volume when the cell angle is considered constant: 𝐿𝑥=𝑉 𝜏/sin 𝛼and

𝐿𝑦=𝑉/(𝜏sin 𝛼). Given that the volume is ﬁxed, one can ﬁnd the appropriate distributions for the

lengths using the standard identity: 𝑃𝑒𝑞 (𝑥)d𝑥=𝑃𝑒𝑞 (𝑦)d𝑦, where 𝑦(𝑥)is some function of 𝑥and 𝑃𝑒𝑞 (𝑥)

is the distribution function of this variable. It can be shown that they are given by the same expression

as for 𝜏:𝑓(𝐿𝜈)∼1/𝐿𝜈, where 𝜈=𝑥, 𝑦. In other words, the sampling distributions of the size in 𝑥

and 𝑦directions obey the same law. This result is of fundamental importance. The invariance with respect

to the swap of 𝑥and 𝑦coordinates is central to simulations of condensed matter (with the exception

of cases where external ﬁelds are imposed that break the symmetry). Any algorithm that violates this

condition should be considered ﬂawed. For instance, uniform sampling of 𝜏assumes that 𝑓(𝜏) ∼ const,

and so one can ﬁnd that 𝑓(𝐿𝑥)∼𝐿𝑥while 𝑓𝐿𝑦∼1/𝐿3

𝑦. Both distributions are incorrect and will lead

to a bias in the sampled ensemble. Another point that should be made regarding the 1/𝜏law concerns

33201-4

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

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

10 玖币 0人已下载

立即下载

摘要：

CondensedMatterPhysics,2022,Vol.25,No.3,33201:118DOI:10.5488/CMP.25.33201http://www.icmp.lviv.ua/journalOnthealgorithmtoperformMonteCarlosimulationsincellswithconstantvolumeandvariableshapeA.Baumketner*InstituteforCondensedMatterPhysicsoftheNationalAcademyofSciencesofUkraine,1SvientsitskiiStr.,7901...

展开>> 收起<<

On the algorithm to perform Monte Carlo simulations in cells with constant volume and variable shape.pdf

共18页,预览4页

还剩页未读，继续阅读

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

On the algorithm to perform Monte Carlo simulations in cells with constant volume and variable shape

相关推荐

开通VIP享超值会员特权

作者详情

相关内容

热门标签

举报选择: