Analyzing and validating simulated tempering implementations at phase transition regimes

doi:10.1016/j.cpc.2020.107256

Computer Physics Communications

Volume 260, March 2021, 107256

https://doi.org/10.1016/j.cpc.2020.107256 Get rights and content

Abstract

This work presents a comprehensive analysis of key elements determining the efficiency of Simulated Tempering implementations for lattice models. Important aspects like the proper choice of the replicas temperature set $ϒ_{R}$ , their quantity $R$ and adequate estimation of the weight factors (establishing the probabilities of jumps between the replicas) are considered. A number of tests to validate distinct ST implementations – as well as to characterize convergence and fluctuations (this latter, after reaching the steady condition) of the thermodynamic estimators – is proposed. Also, by combining the ST method with the so called fixed exchange frequency protocol (to determine $ϒ_{R}$ ), a new procedure to calculate critical temperatures is developed. As cases studies, two commonly addressed lattice models are discussed, the Blume–Emery–Griffiths and Bell–Lavis. Their detailed investigation at the coexistence condition or under first order phase transition regimes allows to properly inspect factors like, pace of convergence, tunneling between replicas and minimal values of the highest replica temperature $T_{R}$ , affecting a good overall performance of the ST method.

Introduction

Exchange Monte Carlo (MC) approaches, such as Parallel Tempering [1] and Simulated Tempering (ST) [2] have been used to study several systems with rough energy landscapes. They are important alternatives to standard numerical simulations – implemented through one-flip, e.g. Metropolis, algorithms – which usually get trapped into local free-energy minima. The central idea in tempering methods is to generate ergodic paths at low temperatures based on the states obtained from higher temperatures. Conceivably, it is easier to visit distinct regions of the phase space (specially if full of entropic barriers, see [3] and the refs. therein) when the temperature is increased. Also, these random walks in distinct energy regions (at large $T$ ’s) should allow the system (at low $T$ ’s) to properly inspect the phase space, yielding appropriate statistical sampling for the calculations of relevant thermodynamic quantities.

In the particular case of the ST, it may be considered an optimized version of the Simulated Annealing, where the temperature is dynamically changed during the simulations run. Remarkably, the generality of the ST relies on the fact that these changes ( $T_{r^{'}} \to T_{r^{''}}$ , for $r^{'}, r^{''} = 1, 2, \dots, R$ ) are performed without a detailed survey of the system potential landscape. But in practice, the process of moving in the phase space (and among the replicas $T_{r}$ ) is not so simple. Exactly because of this, the ST efficiency is strongly dependent on two elements: (a) how to obtain the weight factors $g_{r} = β ψ_{r}$ ’s (with $β = {(k_{B} T)}^{- 1}$ and $ψ_{r}$ the free energy per volume at $T_{r}$ ), fundamental for defining the probabilities of the temperature jumps (refer to Eq. (1) in Section 2); and (b) how to specify the set of temperatures $ϒ_{R} = {T_{1}, T_{2}, \dots, T_{R}}$ of the distinct replicas.

Concerning (b), the correct final steady state of the system can be reached only by means of appropriate tunneling between the replicas, which are mediated by the transition rates $T_{1} \leftrightarrow T_{2} \leftrightarrow \dots \leftrightarrow T_{R}$ . This will be achieved (more rapidly and accurately) through good, ideally optimal, choices for the set $ϒ_{R} = {T_{1}, T_{2}, \dots, T_{R}}$ . Since there is no universal rule of thumb for such task, different methodologies for $ϒ_{R}$ have been proposed in recent works [4], [5], [6]. So, here we discuss two of the most readily protocols for establishing $ϒ_{R}$ . The first is the constant entropy (CE) [5], in which the same value for the entropy difference $Δ s$ between two successive replicas is imposed. The second is the fixed exchange frequency (FEF) — a purely algorithmic procedure [6] — in which given any two replicas $r$ and $r + 1$ , the exchange frequency between them has the same value $f$ .

As for the number of replicas, large $R$ ’s ( $R \geq 5$ ) often guarantee convergence, but at the expense of long computational times. On the other hand, $R = 3$ and $4$ may impose numerical problems (e.g., of weak convergence or strong fluctuations). We examine certain conditions in which $R = 4$ or even $R = 3$ already suffice for adequate simulations. We also address how the estimators statistical distributions, not only at the temperature of interest $T_{1}$ , but also at the other replicas temperatures $T_{r}$ (and specially at $T_{R}$ ) are extremely useful to evaluate the quality of the set $ϒ_{R}$ and of the overall ST implementation performance. Lastly, we unveil a complete new application for the ST plus FEF framework. We show how to construct the frequency map of a lattice model, which in a rather straightforward way, helps to identify the system critical temperatures.

Regarding (a), a simple technique to calculate the $g_{r}$ ’s is the well known Transfer Matrix (TM) method [7]. However, a key issue is how accurately the TM does work for a given specific system. We establish a set of direct tests to estimate if the obtained ST weights are appropriate for simulations. They are constituted by thermodynamic validations of quantities relatively easy to compute (e.g., from basic Simulated Annealing protocols) like entropy and pressure, also not demanding heavy computation.

As specific case studies, we consider two important lattice models: the Blume–Emery–Griffiths (BEG) and the Bell–Lavis (BL) in 2D (but we observe all the present analysis and procedures would also hold in higher dimensions). We assume these systems either at a phase coexistence regime or undergoing a discontinuous phase transition (whose control parameters are either the chemical potential or a magnetic field, this latter just for the BEG). The BEG and BL – already studied under distinct points of view and methods (see next Sections) – display phases which are separated by large entropic barriers, hence demanding more sophisticated simulation protocols. Also, the BEG is a relatively easy to treat system through tempering methods. On the other hand, the BL is harder to converge to the equilibrium. Such important contrast between the models provides a nice way to benchmark our results, allowing to establish a parallel between distinct ST implementations at phase transition conditions.

The work is organized as the following. We present a short overview about the key factors for ST implementations in Section 2, including discussions about the TM approach and the CE and FEF protocols. For the latter, we explain the construction of the corresponding frequency maps. In Section 3 we summarize the models to be studied, BEG and BL. In Section 4 we propose a series of thermodynamic validations to test the accuracy of the ST weights (calculated from the TM). A detailed collection of analysis, describing convergence, fluctuations, critical temperatures and first order phase transition characterization, are given in Section 5. Finally, we draw our final remarks and conclusions in Section 6.

Section snippets

The ST method and some of its distinct implementations

Hereafter, $H$ is the system Hamiltonian with $S$ representing the space of all its possible microstates $σ$ ’s. Also, $σ_{i}$ generically represents the value of the microstate variable $σ$ at the lattice site $i$ (specific proper labeling for $i$ will be set according to each context convenience). Thus, $H_{σ}$ denotes the energy of a specific configuration $σ \in S$ . In very general terms, when simulating $H$ , one is implementing a certain set of evolution rules making the system to pass through a succession of

The test models: BEG and BL

To discuss the simulation methods analyzed in this work, we consider two important discrete lattice models (studied in the literature with different purposes). One is the Blume–Emery–Griffiths (BEG) [25], which considerably extends the Ising model. The other is the Bell–Lavis (BL) [26], a simplified version of a two-dimensional bonded fluid, whose main focus is the orientational properties of hydrogen bonds on triangular lattices.

In both systems, the quantities of interest are the particle

Proper estimation of the ST weights: a thermodynamic validation for the TM method

As previously stated, our aim is to explore distinct ST implementations, so to identify relevant ways to improve the whole protocol. In special, the $g_{r}$ ’s weights (Section 2) constitute one of its key elements. Thence, it is compelling to determine the $g_{r}$ ’s by means of the computationally simple and fast TM approach.

In many lattice models, such as the BEG and BL, the partition function $Z$ can be calculated through a $n$ -layer summation over a (hopefully) reasonable small sampling of microstates

Exchange frequency maps and the identification of critical points

To obtain the replicas temperature set $ϒ_{R}$ from the FEF protocol, we have built a frequency map (see Section 2.2.2), for which few contour plots for both the BEG and BL models are depicted in Fig. 5. These plots are generated by: (1) setting values for $T_{r^{'}}$ and $T_{r^{''}}$ ; (2) running $N_{S T}$ times the ST method (using the Metropolis algorithm) with $ϒ_{2} = {T_{r^{'}}, T_{r^{''}}}$ ; (3) recording the number of times, $N_{r^{'} \to r^{''}}$ , in which a temperature change takes place; and (4) computing the exchange frequency as [6], [40] $f_{r \to r}$

Final remarks and conclusion

In this contribution we have developed an in-depth simulational analysis of distinct aspects and implementations of the simulated tempering (ST) method for systems under phase coexistence conditions. Taking as examples two very well known lattice models in the literature, Blume–Emery–Griffiths (BEG) and Bell–Lavis (BL), we have addressed key issues for the ST efficiency: (a) the calculation of the weight factors $g_{r} = β ψ_{r}$ ’s (used in the expressions for the probabilities of jumps between replicas);

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

C.E.F. and M.G.E. da Luz acknowledge research grants from CNPq (304532/2019-3) and C.E.F from FAPESP (2018/02405-1). A.E.M.M. acknowledges the Organization of American States-CAPES for a Ph.D. scholarship. Finally, we thank the anonymous referees for very constructive suggestions.

References (55)

ValentimA. et al.
Comput. Phys. Comm.
(2014)
FabianV.
Comput. Math. Appl.
(1997)
AnD. et al.
Rel. Eng. Sys. Saf.
(2015)
WosniackM.E. et al.
Comput. Phys. Comm.
(2015)
WosniackM.E. et al.
J. Theoret. Biol.
(2017)
KeskinM. et al.
Physica A
(2006)
GrigelionisG. et al.
Physica A
(1996)
ValentimA. et al.
Comput. Phys. Comm.
(2015)
BergB.A.
Comput. Phys. Comm.
(2002)
IngberL. et al.
Mathl. Compat. Model.
(1992)

HendersonD. et al.

WosniackM.E. et al.

Europhys. Lett.

(2019)

Cited by (1)

Speedup of the Metropolis protocol via algorithmic optimization
2023, Journal of Computational Science
The Metropolis algorithm is widely used in Monte Carlo (MC) simulations in diverse areas of science and technology, especially for problems formulated in terms of lattice models. A common situation is the necessity to perform long sequential processing, e.g., when looking for equilibrium states of distinct physical systems. Hence, even marginal increases in efficiency of the algorithm individual steps can lead to significant reductions in absolute execution runtimes. Usual speedup procedures include hardware updates, parallelization (when possible) and sampling methods. Here we follow a different direction in trying to decrease the full execution times of MC approaches: algorithmic optimization. We show that the algorithms can be improved by implementing relatively few and simple changes in their organization and structure. First, we discuss some refinements for the pseudo-random number generator, addressing the broadly employed Mersenne-Twister algorithm (MT19937-64). Second, we develop a protocol to precalculate the Boltzmann factor, thereby avoiding the high cost of repeatedly calls to this exponential function (indeed, a very recurring step in the standard Metropolis method). To benchmark our proposals we choose the Ising model since it is one of the best known and more extensively studied problems in statistical physics. We consider the mentioned optimizations and different computational elements, like compilers and Hamiltonian variables ranges, testing the efficiency to obtain the system solutions. Our results suggest that the present set of improvement schemes—namely; decreasing the processing time for both, to generate a random number and to implement the one-flip Metropolis step; systematically enforcing optimization for the maximum quantity of algorithm structures accessing random numbers in a code; and considerably reducing the amount of required computations of the MC probabilistic actualization term—might constitute a relevant addition to the existing collection of expediting techniques in MC computational routines.

^☆: The review of this paper was arranged by Prof. D.P. Landau.

View full text

Analyzing and validating simulated tempering implementations at phase transition regimes☆

Abstract

Introduction

Section snippets

The ST method and some of its distinct implementations

The test models: BEG and BL

Proper estimation of the ST weights: a thermodynamic validation for the TM method

Exchange frequency maps and the identification of critical points

Final remarks and conclusion

Declaration of Competing Interest

Acknowledgments

Comput. Phys. Comm.

Comput. Math. Appl.

Rel. Eng. Sys. Saf.

Comput. Phys. Comm.

J. Theoret. Biol.

Physica A

Physica A

Comput. Phys. Comm.

Comput. Phys. Comm.

Mathl. Compat. Model.

Phys. Chem. Chem. Phys.

Europhys. Lett.

Phys. Rev. E

J. Chem. Phys.

J. Chem. Phys.

J. Chem. Phys.

J. Chem. Phys.

J. Comput. Chem.

J. Phys. Soc. Japan

ACM Trans. Model. Comput. Simul.

Phys. Rev. B

Phys. Rev. Lett.

Phys. Rev. E

J. Chem. Phys.

Science

Europhys. Lett.