Global Thermodynamic Properties of Complex Spin Systems Calculated from Density of States and Indirectly by Thermodynamic Integration Method

Evaluation of global thermodynamic properties, such as the entropy or the free energy, of complex systems featuring a high degree of frustration or disorder is often desirable. Nevertheless, they cannot be measured directly in standard Monte Carlo simulation. Therefore, they are either evaluated indirectly from the directly measured quantities, for example by the thermodynamic integration method (TIM), or by applying more sophisticated simulation methods, such as the Wang-Landau (WL) algorithm, which can directly sample density of states. In the present investigation we compare the performance of the WL and TIM methods in terms of calculation of the entropy of an Ising antiferromagnetic system on a Kagome lattice - a typical example of a complex spin system with high geometrical frustration resulting in a non-zero residual entropy the value of which is exactly known. It is found that in terms of accuracy the implementationally simpler TIM can deliver results comparable with the more involved WL method.


Introduction
Calculation of global thermodynamic properties, which cannot be measured in Monte Carlo (MC) simulation directly, such as the free energy and the entropy, is generally a difficult task. A brute force approach of scanning the entire configuration space to obtain density of states is feasible only for sufficiently small systems. However, the exponential increase of the configurational space with the system size N makes this approach intractable even for moderate sizes and small number of degrees of freedom, such as the Ising model with the configuration space increasing as 2 N .
In statistical physics commonly used standard MC methods, such as the Metropolis algorithm (MA) [1], allow direct evaluation of several thermodynamic quantities, such as the internal energy or magnetization, but not global ones, such as the free energy and the entropy. One possible approach that allows a direct evaluation of density of states (DOS) and consequently the entire thermodynamics is the Wang-Landau (WL) algorithm [2], which has been successfully applied to a variety of problems, e.g., an efficient study of first-order and second-order phase transitions. It is a powerful tool for the investigation of systems with rough energy landscapes with large energy barriers separating local minima, which make the use of other standard methods infeasible. An alternative indirect approach to calculation of the global thermodynamic quantities, that avoids the calculation of DOS, is the so-called thermodynamic integration method (TIM) [3]. Since its introduction in 1977 it has been sparingly used, even though several studies pointed to its competitiveness [4,5], for example in calculation of the ground-state entropy of some typical disordered/frustrated spin systems, such as the ±J Ising model and the spin-s triangular lattice Ising antiferromagnet.
In the present study, we compare the performance of the WL and TIM methods in terms of calculation of global thermodynamic quantities of a highly frustrated Ising antiferromagnet on a Kagome lattice (IAKL) [7] with the focus on the entropy, the ground-state value of which is exactly known [6].

IAKL model
The Hamiltonian of the studied spin s = 1/2 IAKL system is given by where the first summation goes over the nearest neighbors and σ i = ±1 is the spin at the ith site. In order to introduce frustration, interactions between neighboring spins were chosen to be antiferromagnetic (J < 0). A schematic illustration of the Kagome lattice is shown in the inset of Fig. 1(a). IAKL is a typical example of a complex spin system with high geometrical frustration resulting in a massive ground-state degeneracy with a finite residual entropy and no long-range ordering at any temperature.

MA and TIM
MA is a well-known, general, easy to implement and therefore widely used MC method [1]. The algorithm performs a random walk in the energy space. In every MC step a new state with the energy H is proposed and accepted with the probability p(s old → s new ) = min(1, e −β ∆H ), where ∆H is the energy difference between the new state and the old state, β = 1/(k B T ) is the inverse temperature, and k B is the Boltzmann constant (hereafter set to k B = 1). MA can be used to directly calculate and investigate several quantities, such as the internal energy e = H /N and magnetization m = M /N, N is the number of spins and . . . denotes a thermal average. The entropy of the magnetic system with a discrete spin number s can be obtained as a function of the inverse temperature by TIM [3] as: where E = Ne. Assuming equilibrium conditions, thermal averages calculated based on a fixed number of MC sweeps N sweeps , for a given temperature range [T 1 , T N T ] and a fixed lattice size L the only relevant parameter that can influence accuracy of the entropy estimation is the temperature mesh density, characterized by the number of temperature points N T and their distribution. Generally, a denser mesh (large N T ) leads to a smaller quadrature error and thus a more accurate estimation.

WL algorithm
WL method is a relatively new MC method producing accurate results, including the global thermodynamic functions [2]. A random walk is performed in the energy space to extract an estimate of DOS, g(E), from which one can calculate the partition function at any temperature and consequently all other thermodynamic quantities. In particular, the partition function can be obtained as: where the summation goes over all possible energy values E. Consequently, a mean value of any thermodynamic quantity, including the entropy, can be evaluated by using the standard statistical physics relations. For a given lattice size, the user-defined parameters in the WL method are the flatness criterion F C < 1 and the modification factor f final > 1 [2]. Generally, the closer are the values of F C and f final to 1.0 the more precise results can be expected. Typically chosen values include F C = 0.8 or 0.9 and f final = 1 + 10 −k , for k = 8, 9, and 10.

Results and discussion
We performed several simulations using both MA and WL to calculate the entropy per spin (entropy density) and the free energy. Below, we only present the former quantity, as the latter one is just a simple function of the former. The presented results were calculated as averages obtained from 10 independent runs. Fist, we present results obtained by MA, with the following parameters: N sweeps = 5 × 10 5 and [T 1 , T N T ] = [0.0064, ∞] (or [β 1 , β N T ] = [156.3731, 0]). The temperature mesh was chosen nonuniformly with the largest density of points in the region of the largest variation of the energy E(β ) in order to increase the precision of the numerical quadrature (Eq. (2)). In Figure 1(a) we examine the lattice size dependence of TIM with N T = 301, for L = 16, . . . , 84. All the curves are found to collapse on a single curve within the error bars and, hence, we can conclude that the finite-size effects are very small. In Figure 1(b) we study the influence of the parameter N T , for a fixed value of L = 32. N T = 151 and 76 case were obtained from the initial N T = 301 cases by repeatedly removing every second nod from the previous denser grid. Again, all the curve appear to coincide within statistical errors. Nevertheless, as shown in the inset, by comparing the ground-state value estimate S GS MA /N with the exact value S GS exact /N = 0.5018 [6] one can notice a gradual improvement with the increasing N T , even though accuracy is fairly high for all the values of N T .
In Figure 2(a) we present the WL results for different L = 16, . . . , 84 and the simulation parameters set to F C = 0.8 and f final = 1 + 10 −8 . Like for the TIM results in Figure 1(a), all the curves coincide within the error bars. There is also a good coincidence between the WL and TIM results, as evidenced in the inset of Figure 2(a) that shows the deviation of the residual entropies obtained by the WL method and TIM (for N T = 301) from the exact value. Finally, in Figure 2(b) we demonstrate the effect of the choice of the WL method parameters F C and f final . Again, for the standard values the differences are within the error bars.
To conclude, when the standard values of parameters are chosen, both TIM and WL methods deliver sufficiently accurate results that mutually indistinguishable. The advantage of the TIM is its implementational simplicity but the simulation has to be run at all desired temperatures. On the other hand, having obtained DOS from the WL simulation enables straightforward calculation of any thermodynamic quantity at all temperatures.