Application of Local Information Entropy in Cluster Monte Carlo Algorithms

The chapter refers to a modification of the so-called adding probability used in cluster Monte Carlo algorithms. The modification is based on the fact that in real systems, different properties can influence its clusterization. Finally, an additional factor related to property disorder was introduced into the adding probability, which leads to more effective free energy minimization during MC iteration. As a measure of the disorder, we proposed to use a local information entropy. The proposed approach was tested and compared with the classical methods, showing its high efficiency in simulations of multiphase magnetic systems where magnetic anisotropy was used as the property influencing the system clusterization.


Introduction
The problem of simulation magnetization processes of multiphase magnetic materials seems to be important, regarding the tendency for applications of highefficient permanent magnets with reduced or without rare earth elements. Recently, we reported unique hard magnetic properties of Tb-Fe-B-Nb bulk alloys, i.e., coercive field over 7 T at room temperature, attributed to a specific microstructure of dendrite-like Tb 2 Fe 14 B grains [1,2]. In this system, Tb and Fe magnetic moments are coupled antiferromagnetically, which is responsible for relatively low magnetic remanence (μ 0 M r ≈ 0.3 T) and in consequence |BH| max (about 13 kJ/m 3 ). However, the Fe-Nb-B-Tb bulk alloys can be considered as a material with extremely high resistance to the external magnetic field and can be a source of magnetic anisotropy in powders as well as bulk spring-exchange composites containing magnetically soft and ultra-high coercive phases. For this reason, the ability of simulating of such systems is very useful in the process of examining and designing of spring-exchange composites in the pre-lab phase.
Monte Carlo (MC) algorithms are now widely used to clarify various physical phenomena, as well as to investigate their potential application in modern technology. Among many simulation methods, the Metropolis MC (MMC) approach [3,4] is especially attractive in statistical physics for the determination of system physical quantities in thermodynamic equilibrium. The MMC algorithm realizes an ergodic stochastic process, ensuring the fulfillment of the detailed balance condition.
One of the bright examples of the application of the MMC algorithm is the Ising model of spins located on the nodes of some lattice. Indeed, using this method one can study a course of magnetic ordering and its dependence of temperature and details of interactions between the spins.
The MMC method utilizes the single-spin-flip procedure to change the spin configuration; however, in many cases (e.g., simulations of magnetization processes) a more effective algorithm is needed. The simplest approach relies on the generation of a cluster of uniformly oriented spins and their subsequent flip to reach new state of the system. The main question is how to determine the cluster and how to establish a rule of its acceptance, simultaneously satisfying the detailed balance condition. Classical approaches, based on the Kasteleyn-Fortuin theorem [5,6], were proposed by Swendsen and Wang (SW) [7] as well as by Wolff [8] who assumed a specific cluster-building procedure controlled by the so-called adding probability. It is known that the cluster Monte Carlo methods (CMC) are very efficient in the analysis of critical phenomena, e.g., transformation from ferromagnetic to paramagnetic phase [9,10]. In contrast, their application for studying magnetization processes of systems far below the Curie point produces artificial results, which can be demonstrated for the systems containing magnetically different phases (e.g., hard and soft) as well as geometrical irregularities. The clusterbuilding algorithms implemented within SW and Wolf approaches are steered by the exchange interactions and the system temperature, but they are not sensitive to other features, potentially affecting the clusterization of spins.
In order to broaden the applications of the CMC methods to simulations of real magnetic composites, we proposed a new method based on some modification of the SW/Wolff adding probability and a particular Metropolis-like algorithm, ensuring the principle of detailed balance [11]. The idea is based on the fact that some kind of regions of the system, characterized by a local disorder of selected system property, constitutes natural barriers for the extension of clusters. In the case of magnetic multiphase composites, spatial distribution of the magnetic anisotropy can be considered as the property affecting the cluster formation.
In the chapter, the disorder-based CMC algorithm is introduced and discussed in a context of classical CMC methods. We show that the new simulation procedure is efficient leading to physically reliable results, especially for multiphase magnetic composites.

Local disorder-based CMC method
Let the considered physical system be characterized by a discrete spectrum of microscopic energy states (microstates) labeled by E α . Furthermore, let the system be in equilibrium with a thermostat having a temperature T. According to basic principles of statistical mechanics, the probability that the system occupies the state α is proportional to the Boltzmann factor: exp ÀβE α ð Þ, where β ¼ 1=k B T, k B refers to the Boltzmann constants. Then, the equilibrium value of some system quantity (observable) F can be calculated using the following formula [12,13]: where Z ¼ P α exp ÀβE α ð Þis the partition function. The direct use of Eq. (1) is impractical due to a very large number of states that should be taken into account. Indeed, even for small systems as, for example, a two-dimensional 10 Â 10 lattice of spins, we get a total of 2 100 states-the number that makes the summation occurring in Eq. (1) impossible.
In order to estimate the average´F, a nonuniform sampling of the system states can be applied. If α f g denotes a set of indices of M system states selected with the probability p α , the equilibrium value of the observable F can be modeled [13] by the estimator: where probabilities p α were selected to be equal exp ÀβE α ð Þ=Z. All we need is a method that generates a set of system states with the Boltzmann probabilities p α . Because the exact value of the partition function Z is unknown, the generation of the states is usually carried out by the ergodic Markov process. This process produces a proper chain of states under the assumption that transition probability W αβ (from α to β state) is independent of the states preceding α. Moreover, It is also assumed that the detailed balance condition, p β W βα ¼ p α W αβ , is satisfied when the system is in a state of equilibrium [12][13][14][15].
The transition probability W αβ can be considered as a product of the selection probability g αβ and the acceptance ratio (probability)A αβ . In general, the selection probabilities can be chosen to a large extent freely, e.g., they can be symmetrical g αβ ¼ g βα [12]. In that case the acceptance probabilities satisfy the equation: As an example, let us consider Ising model of N interacting spins placed at the nodes of a two-dimensional regular. The energy of the system is, then, given by the formula where s i ¼ AE1 describes the spin state at the ith lattice node and J refers to the exchange integral. In order to determine the physical properties of our magnetic system, the Metropolis algorithm can be employed. It relies on the particular choice of both the selection probability and the acceptance ratio. Having some configuration of spins, the next one is obtained by the flip of a single spin (single-spin-flip algorithm) [12,13]. This procedure results in uniform distribution of the selection probabilities, i.e., each new state participates in simulations with probability g αβ ¼ 1=N. Then, the new spin configuration can be accepted or rejected with an acceptance ratio . Although the Metropolis algorithm can be applied to a variety of physical problems, when applied to magnetic systems, it has disadvantage that relies on a very rapid increase of the correlation time as well as correlation length near the critical point. As a result the system contains domains of the same oriented spins and therefore becomes configurationally frozen. This unexpected behavior (critical slowing down) of the Metropolis algorithm is the reason for the difficulties in the generation of statistically independent spin configurations that are needed for the calculation of the estimator F S .
The solution of the critical slowing down problem was proposed by Swendsen-Wang [7] and later by Wolff [8]. The approach developed by Wolf (cluster-flipping algorithm) based on the generation of the uniformly oriented spin cluster and its subsequent flipping [12]. In contrast to single-spin-flip algorithm, this procedure easily destroys domains of correlated spins and allows the system to walk through the configuration space. The Wolff algorithm is recognized to be more effective than the Swendsen-Wang one [16][17][18]. The selection of a spin cluster starts from randomly chosen spin to which the neighbors occupying the same spin state are added with the probability P add . The cluster grows up until no spin is added to it. It is a great advantage that the Wolff algorithm is a rejection-free one. Indeed, adding probability P add is defined so that the detailed balance condition (with acceptance ratio equal to one) is met: Despite the great achievements, the Wolff model is not able to correctly simulate magnetic phenomena that occur far below the critical point in real magnetic materials, including those that are composed of various magnetic phases as well as those containing geometrical irregularities. To be more precise, one can consider, as an example, remagnetization of the system build of two magnetically hard and soft ferromagnetic spheres coupled by a narrow bridge. Let us assume that in the initial state, the magnetization of the system is collinear with the direction of the external magnetic field and then the magnetic field is switched in the opposite direction. What happens to the system is that the magnetization of the soft sphere will follow the change of the magnetic, and then the similar behavior of the hard sphere is expected. Unfortunately, for systems with strong spin-spin coupling, the two-step behavior of the considered system cannot be modeled using Wolff clusterization algorithm. Indeed, every attempt to build a cluster within the two-sphere magnetic system results in that the hard and soft spheres belong to the same cluster independently on magnetic anisotropy and geometry of the system. So, even if the two-step remagnetization process is energetically preferred, the simulated magnetization curve consists of one step related to the common spin rotation. Taking into account the problems encountered during modeling the remagnetization of magnetically inhomogeneous system, we propose a modification of the Wolff algorithm that relies on an assumption that some regions of the system, characterized by a disorder of selected system property, can serve barriers for an extension of magnetic clusters.
We began from the introduction a distribution of system property K that can refer, for example, magnetic anisotropy or some other properties potentially affecting the clusterization of the system. Let us define a sphere V around a node of the spin lattice-the sphere containing N nodes in total. Furthermore, let the system property be characterized by a discrete and finite set of values: K 1 , K 2 , …, K N K . The local distribution of the K system property is defined by the numbers n i f g ¼ n 1 ; n 2 ; …; n N K f g , where n i stands for the number of nodes having the value K i . The sum overall is n i equals N. Thus, the particular K-state of the sphere V is defined by the set of numbers n i f g. The number of possible realizations of the K-state denoted by Ω n i f g is given by the expression Using the Stirling formula, the previous equation takes the form where p i ¼ N i =N stands for the probability of finding K i value inside the sphere V. The expression on the right side of Eq. (7) contains the term called local information entropy [19]: So the number of "microstates" corresponding to the K-state is The probability of particular local distribution n i f g of the system property K is proportional to Ω n i f g . If the only one value of K dominates its distribution inside the sphere V, e.g., N; 0; …0 f g , the local information entropy achieves the lowest value S loc ¼ 0. Other local distribution n i f g that involves more than one value of K results in S loc > 0. Now we can introduce the perception of local disorder of the K system property and its measure p K . We will say that the K system property is fully ordered inside the sphere V if all the nodes inside the sphere have the same value of K. Consequently, one can say that the sphere V exhibits some local disorder of K if more values of K are distributed inside the sphere. The measure of the local disorder can be defined by as follows: which gives p K ¼ 1 if the case of fully ordered system property (S loc ¼ 0). One can see that p K decreases when local disorder of K system property increases (S loc > 0).
We assume that local disorder of K system property decreases the probability of cluster growth and constitutes barriers for cluster expansion. Thus, one could expect nonuniform distribution of adding probability P 0 add among the system nodes. In order to derive a reasonable formula for P 0 add , we follow an analogy with derivation of the van der Waals equation. Our "ideal gas" corresponds to the uniform distribution of adding probability given by Wolff algorithm, while "real gas" refers to the case when the distribution of adding probability is affected by local disorder of some system property. In order to still play with Wolf algorithm, an influence of local disorder should be compensated in some way, e.g., by multiplication P 0 add by p À1 k . Then, assuming that this product satisfies the Wolff's formula P 0 add exp S loc À Á ¼ 1 À exp À2βJ ð Þ, the adding probability that takes into account the impact of local disorder of some system property can be easily obtained: The above equation is the clue of our method which expresses the decrease of the classical adding probability by the disorder-based factor. A set of property, for which the disorder is determined (by the local information entropy), should be chosen dependently on a specific application. The α parameter expresses some enhancement or weakness of the influence of the entropy on the system clusterization. It should be defined regarding the specific problem (see next paragraphs).

Simulation procedure
The main simulation procedure consists of a series of Monte Carlo steps, and it is based on the classical Metropolis algorithm applied to the so-called spin continuous approach. However, the main difference lies in the cluster-building procedure which is executed with a small probability P cl , as shown in Figure 1.
First of all, a random node i of the system is chosen, and then it is decided, regarding the probability P cl , whether to analyze the selected node or build a cluster, starting from the node as a seed. In the first case, the algorithm goes to the typical Metropolis procedure and spin continuous method, i.e., the spin direction is randomly modified by angle AEθ, and then the energy difference ΔE between the new and the old configuration is calculated. The energy of the system is computed in the frame of the 3D Heisenberg model: where J ij is the exchange parameter, S i is the spin vector on site i, K i is the anisotropy constant (per site), n i is the easy magnetization axis, g is the Lande factor, μ B is the Bohr magneton, μ 0 is the vacuum permeability, H i is the magnetic field on site i, D is the dipolar constant, e ij is the directional versor between the ith and jth nodes, and r ij is the distance between the ith and jth nodes. The change is accepted with the Metropolis acceptance probability. In the second case, the algorithm is similar; however, it is necessary to find a cluster around the selected node and, if it exists, carry out a coherent rotation of all spins belonging to the cluster. The procedure of the cluster building is the key and important point. We used the Wolff algorithm but with the modified adding probability P 0 add expressed by the formula where E coupling is the exchange interaction energy between the neighboring spins. The local information entropy is computed in the defined sphere around each node. In general, the choice of properties (used in the entropy calculation) will depend on the problem being considered; however, for magnetic systems, the natural limit of cluster growth is the change in the value and direction of magnetic anisotropy K. Therefore, an optimal feature is a set of three components [K x , K y , K z ], whereby the nonmagnetic nodes do not participate in the entropy calculations. In addition, at the beginning of each cluster-building procedure, the α coefficient is drawn that will weaken (α < 1) or strengthen (α > 1) the influence of a local property disorder on the system clusterization. Moreover, it is recommended that from time to time (typically about 20% cases), it completely ignores the impact of disorder and builds a cluster based on the standard Wolff method.
The presented algorithm gives some freedom of the cluster building for which the thermodynamic balance is fulfilled. Indeed, even if the cluster rotation slightly disturbs the balance, the remaining single-spin MCM iterations restore it again. The only condition is that the P cl and θ values are relatively low (it can be determined experimentally for cases when the results do not depend on the parameters).

Parameter Description
Recommendation and notes θ The angle at which the direction of the spin or cluster can be changed in a single iteration Note, it affects the energy changes in a single step. A large θ value will freeze the system, while the small will make it very loose as well as more iterations will be needed to obtain the same change P cl Probability of cluster analysis instead of spin in a single iteration If many independent clusters in the system are expected, then the value of the parameter should be raised to give everyone a chance to be analyzed. However, high values may destroy the thermodynamic equilibrium S range The range of local information entropy The range of local entropy should be selected depending on the geometry of the system. The larger range increases the chance to separation between the magnetic grains and a small common area For α < 1 and for α > 1, the S loc impact will be depressed and strengthened, respectively. The best approach is to draw an α value before each cluster search The probability of searching for a cluster based only on the Wolff method In some cases, you can build a cluster ignoring the impact of local entropy. Typically, we suggest P Wolff = 0.2 Table 1.
A guide of the parameters used in the algorithm.
Finally, one MC step consists of N iter iterations presented in Figure 1, and the N step steps are taken to obtain the magnetization of the system <S z > (the average spin in the external field direction). The most important parameters of the algorithm, their possible values as well as our recommendations from a practical point of view are summarized in Table 1.
The algorithm can be, in some aspects, paralleled [20][21][22][23]. The main limitations of the parallelization process are due to three reasons. First of all, each step in the Monte Carlo procedure should be based on the system modified in the previous step. In particular, the decision to accept the new spin direction depends on the current direction of the neighboring spins due to the exchange energy. Consequently, a situation in which different threads are testing a new configuration for the two neighboring spins at the same time should be refused. Similarly, simultaneous analysis of the whole cluster and individual spins inside it (as well as spins interacting with it) is not allowed. The second thing to keep in mind is to ensure that each state of the system can be selected with the same probability. This means that none of the actions can interfere with the probability of choosing a spin for analysis. For example, the spin cannot be temporarily omitted in the draw as well as the analyzed cluster must always be found in the real time and cannot be taken from a fixed database. Finally, it should be taken into account that a change of the spin direction inside a cluster, which is constructed based on the P 0 add probability, may disrupt the thermodynamic balance of the system. Therefore, from the statistical point of view, after each change of the cluster spins, there should be many iterations analyzing single spins in the cluster and whole system. Despite all the above restrictions, there are three time-consuming problems in the algorithm which can be parallelized: the main MC loop, the calculation of cluster energy as well as the cluster-building procedure around the chosen spin. Our experience shows that the parallelization of the algorithm accelerates the computation time Comparison of computation time (400 MC steps) for different P cl , system size, and its magnetic properties. more than 10 times (depending on the probability of cluster analysis and the system size) and maintains the correctness of results. Figure 2 shows a comparison of computation time (needed for 400 MC steps) for different P cl , system size, and its magnetic properties. On this graph, the independent quantity is a number of threads applied for the simulations (HP ProLiant DL580 G7, 4x Intel Xeon 8C X7560 2.27 GHz).
Generally, this aspect of calculations is very complex. Nevertheless, one can conclude that the parallelization is more effective for the systems with small number of the cluster-building trials.

Application of the method to multiphase magnetic systems
In order to show an efficiency of the disorder-based cluster MC algorithm and its comparison with the classical Wolff method, the two systems containing magnetically hard and soft phases were analyzed. The 3D system space consists of 50 Â 50 Â 50 (1,25,000) nodes; the spins are arranged in shapes of joined spheres with hard and soft magnetic properties. The bonding between the spheres, further called "bridges," is different for the two cases, and it equals 1 or 7 nodes, receptively. Figure 3 depicts an example of the system with one-bridge coupling between the hard and soft spheres. The parameters of the magnetic phases and simulation procedure are listed in Table 2.
The difference between classical and disorder-based approaches lies in the definition of adding probability. It can be demonstrated by the cluster-building successes when the procedure starts from the same cluster "seed." Figures 4 and 5 show the comparison for the one-and seven-bridge systems, respectively. The cluster seed is placed in the center, and the white color is attributed to 100 clusterbuilding successes per 100 trials.
As shown, using the classical Wolff algorithm, all spins in the system belong to the same cluster independently on their magnetic properties. This means that during simulations, a possibility of closing the clusters inside the hard or soft phase is not available. The modification of adding probability causes that the clusterbuilding procedure can distinguish the hard and soft phase with relatively high System size (in one direction) n 50 Number of iteration in one MC step 3n 3 = 3,75,000 Number of MC steps for magnetization averaging 400 Table 2.
Parameters of the analyzed systems and simulation procedure.   Reverse magnetization curves simulated using the classical P add and modified P 0 add adding probability.
probability. However, with lower probability, the cluster contains whole spins in the system (like in the Wolff algorithm) which is related to a random value of the α parameter. The modification of adding probability causes that the configurations necessary for modeling of real magnetization processes are available during MC iterations. In other words, if magnetization processes require separate behavior of the hard and soft phase, the simulation procedure will test such a possibility. As a final test, the so-called reverse magnetization curves were simulated using the classical as well as modified adding probability. Initially, the system was saturated in the field direction (all spins are directed up), and then the field was switched off. During calculations the field was increased in the opposite direction. Magnetization, determined as average spin value in the field direction, as a function of the external magnetic field for all examined cases is shown in Figure 6.
The difference between P add and P 0 add appeared in all studied examples. Note that for the one-bridge system (relatively low interactions between the hard and soft phases), the two-step reverse magnetization curve is expected. The first step is attributed to the spin flip of the soft phase, while the second one related to spin flip of the hard phase. For the seven-bridge system, the reverse magnetization is different due to stronger interphase coupling. However, at temperatures k B T above 1e-6 eV, this process consists of a subsequent change of the "soft" spins and next, coherent flip of the "hard" spins. Such a scenario is fully confirmed by the spin configurations depicting the reverse magnetization process for the seven-bridge system at k B T = 1e-4 eV (see Figure 7).
Regarding the fact that in some cases also for P add the two-step behavior occurs, the main question is which curve is physically correct. For this reason, it is worth to analyze the curves in fields when the soft phase changed spin direction using P 0 add in the cluster-building procedure. It is known that thermal equilibrium is related to a minimum of the free energy of the system. Comparing the spin configurations for the Wolff and our modified algorithm, one can state that entropy (as a thermodynamic function) is higher for the second one, i.e., applying P 0 add , which contributes to a decreasing of the free energy. Therefore, this approach is correct for which the energy of the system is lower. Figure 8 shows a focus of the interesting region as well as the energy of the system for the two tested algorithms. In all cases the application of P 0 add results in "faster" (in lower fields) minimization of the system energy than the procedure based on the classical adding probability.
It is clear that the introduced modification of adding probability results in more effective finding of the free energy minimum and, therefore, produces physically reliable results which allows modeling multiphase magnetic systems.

Conclusions
As it was shown, in some conditions (low temperature and/or strong exchange coupling) the classical CMC algorithm can produce incorrect results in the application of multiphase magnetic systems. The problem lies in the adding probability used in the cluster-building procedures. This problem can be overcome by the proposed modification of P add accompanied by the specific Metropolis-like algorithm. The key point is to recognize a property that influences the system clusterization and, next, introduction of the local information entropy of this property into the adding probability.
In summary, one can conclude that the method proposed for modeling of magnetic systems allows sampling spin configurations practically inaccessible by original Wolff clusterization procedure. An interesting future of our method is the acceleration of the energy minimization process.