Application of Local Information Entropy in Cluster Monte Carlo Algorithms

Abstract

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.

1. Introduction

The problem of simulation magnetization processes of multiphase magnetic materials seems to be important, regarding the tendency for applications of high-efficient 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₂Fe₁₄B grains [1, 2]. In this system, Tb and Fe magnetic moments are coupled antiferromagnetically, which is responsible for relatively low magnetic remanence (μ₀M_r ≈ 0.3 T) and in consequence |BH|_max (about 13 kJ/m³). 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 cluster-building 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.

2. Local disorder-based CMC method

Let the considered physical system be characterized by a discrete spectrum of microscopic energy states (microstates) labeled by $E_{\alpha }$. 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: $\mathit{exp}\left(-\beta E_{\alpha }\right)$, where $\beta =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]:

$$\overline{F}=\frac{1}{Z}\sum _{\alpha }{F}_{\alpha }\mathit{exp}\left(-\beta E_{\alpha }\right),$$E1

where $Z={\sum }_{\alpha }\mathit{exp}\left(-\beta E_{\alpha }\right)$ 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¹⁰⁰ states—the number that makes the summation occurring in Eq. (1) impossible.

In order to estimate the average $\acute{F}$, a nonuniform sampling of the system states can be applied. If $\left\{\alpha \right\}$ denotes a set of indices of $M$ system states selected with the probability $p_{\alpha }$, the equilibrium value of the observable $F$ can be modeled [13] by the estimator:

$$F_S=\frac{{\sum }_{\left\{\alpha \right\}}{F}_{\alpha }\mathit{exp}\left(-\beta E_{\alpha }\right)p_{\alpha }^{-1}}{{\sum }_{\left\{\alpha \right\}}\mathit{exp}\left(-\beta E_{\alpha }\right)p_{\alpha }^{-1}}=\frac{1}{N}\sum _{\left\{\alpha \right\}}{F}_{\alpha },$$E2

where probabilities $p_{\alpha }$ were selected to be equal $\mathit{exp}\left(-\beta E_{\alpha }\right)/Z$. All we need is a method that generates a set of system states with the Boltzmann probabilities $p_{\alpha }$. 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_{\mathit{\alpha \beta }}$ (from α to β state) is independent of the states preceding α. Moreover, It is also assumed that the detailed balance condition, $p_{\beta }{W}_{\mathit{\beta \alpha }}=p_{\alpha }{W}_{\mathit{\alpha \beta }}$, is satisfied when the system is in a state of equilibrium [12, 13, 14, 15].

The transition probability $W_{\mathit{\alpha \beta }}$ can be considered as a product of the selection probability $g_{\mathit{\alpha \beta }}$ and the acceptance ratio (probability)$A_{\mathit{\alpha \beta }}$. In general, the selection probabilities can be chosen to a large extent freely, e.g., they can be symmetrical $g_{\mathit{\alpha \beta }}=g_{\mathit{\beta \alpha }}$ [12]. In that case the acceptance probabilities satisfy the equation:

$$A_{\mathit{\alpha \beta }}/A_{\mathit{\beta \alpha }}=\mathit{exp}\left(-\beta \left(E_{\beta }-E_{\alpha }\right)\right)$$E3

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

$$E=-J\sum _{i\ne j=1}^N{s}_i{s}_j$$E4

where $s_i=\pm 1$ describes the spin state at the $i$th 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_{\mathit{\alpha \beta }}=1/N$. Then, the new spin configuration can be accepted or rejected with an acceptance ratio $A_{\mathit{\alpha \beta }}=\mathit{exp}\left(-\beta \left(E_{\beta }-E_{\alpha }\right)\right)$ (for $E_{\beta }-E_{\alpha }>0$) and $A_{\mathit{\alpha \beta }}=1$ (for all other cases). 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_{\mathrm{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_{\mathrm{add}}$ is defined so that the detailed balance condition (with acceptance ratio equal to one) is met:

$$P_{\mathit{add}}=1-\mathit{exp}\left(-2\mathit{\beta J}\right)$$E5

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,\dots ,K_{N_K}$. The local distribution of the $K$ system property is defined by the numbers $\left\{n_i\right\}=\left\{n_1,n_2,\dots ,n_{N_K}\right\}$, 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 $\left\{n_i\right\}$. The number of possible realizations of the $K$-state denoted by ${\Omega }_{\left\{n_i\right\}}$ is given by the expression

$${\Omega }_{\left\{n_i\right\}}=\frac{N!}{n_1!n_2!\dots n_{N_K}!}$$E6

Using the Stirling formula, the previous equation takes the form

$$\mathit{ln}\left({\Omega }_{\left\{n_i\right\}}\right)\approx -N\sum _{i=1}^{N_K}{p}_i\mathit{ln}\left(p_i\right)$$E7

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]:

$$S_{\left\{n_i\right\}}^{\mathit{loc}}=-\sum _{i=1}^{N_K}{p}_i\mathit{ln}\left(p_i\right)$$E8

So the number of “microstates” corresponding to the $K$-state is

$${\Omega }_{\left\{n_i\right\}}\approx \mathit{exp}\left(N{S}_{\left\{n_i\right\}}^{\mathit{loc}}\right)$$E9

The probability of particular local distribution $\left\{n_i\right\}$ of the system property $K$ is proportional to ${\Omega }_{\left\{n_i\right\}}$. If the only one value of $K$ dominates its distribution inside the sphere $V$, e.g., $\left\{\mathrm{N},0,\dots 0\right\}$, the local information entropy achieves the lowest value $S^{\mathit{loc}}=0$. Other local distribution $\left\{n_i\right\}$ that involves more than one value of $K$ results in $S^{\mathit{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:

$$p_K=\mathit{exp}\left(-S^{\mathit{loc}}\right)$$E10

which gives $p_K=1$ if the case of fully ordered system property ($S^{\mathit{loc}}=0$). One can see that $p_K$ decreases when local disorder of $K$ system property increases ($S^{\mathit{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_{\mathit{add}}^{\prime }$ among the system nodes. In order to derive a reasonable formula for $P_{\mathit{add}}^{\prime }$, 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_{\mathit{add}}^{\prime }$ by $p_k^{-1}$. Then, assuming that this product satisfies the Wolff’s formula $P_{\mathit{add}}^{\prime }\mathit{exp}\left(S^{\mathit{loc}}\right)=1-\mathit{exp}\left(-2\mathit{\beta J}\right)$, the adding probability that takes into account the impact of local disorder of some system property can be easily obtained:

$$P_{\mathit{add}}^{\prime }=\left[1-\mathit{exp}\left(-2\mathit{\beta J}\right)\right]\mathit{exp}\left(-{\mathit{\alpha S}}^{\mathit{loc}}\right)$$E11

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).

3. 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.

Figure 1.

Schematic diagram of the used algorithm.

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 ±θ, 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:

$$E=-\sum _{i,j}{J}_{\mathit{ij}}{S}_i·S_j-\sum _i{K}_i{\left(S_i·n_i\right)}^2-g{\mu }_B{\mu }_0\sum _i{H}_i·S_i+D\sum _{i,j}\frac{S_i·S_j-3\left(S_i·e_{\mathit{ij}}\right)\left(S_j·e_{\mathit{ij}}\right)}{r_{\mathit{ij}}^3}$$E12

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, μ₀ 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_{\mathit{add}}^{\prime }$ expressed by the formula

$$P_{\mathit{add}}^{\prime }=\left[1-\mathit{exp}\left(-{\mathit{\beta E}}_{\mathit{\text{coupling}}}\right)\right]\mathit{exp}\left(-{\mathit{\alpha S}}^{\mathit{loc}}\right)$$E13

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). 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.

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
Pcl	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
Srange	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
Niter	The number of iterations contained in one MC step	Typically, we suggest Niter = 3 N where N is the number of all spins in the system
Nstep	The number of MC steps used to calculate the average spin of the system <SZ>	We propose Nstep = 400
α	Modifier of the impact of local information entropy	For α < 1 and for α > 1, the Sloc impact will be depressed and strengthened, respectively. The best approach is to draw an α value before each cluster search
PWolff	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 PWolff = 0.2

Table 1.

A guide of the parameters used in the algorithm.

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_{\mathit{add}}^{\prime }$ 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 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).

Figure 2.

Comparison of computation time (400 MC steps) for different P_cl, system size, and its magnetic properties.

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.

4. 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.

Figure 3.

Example of the system with one-bridge coupling between the hard and soft spheres. The graphs on the right show cross sections (z-x plane) for different y values equal to 10 (a), 25 (b) and 40 (c).

Parameter	Soft phase	Hard phase
Exchange coupling J	1.5e–2 eV	1.5e–2 eV
Anisotropy constant K	0	1e–3 eV
	Simulation procedure
Dipolar constant D	1.8e–7 eV
Thermal energy kBT	1e–4, 1 e–5, 1e–6 eV
α parameter	Random value from 0 to 25
Cluster analysis probability Pcl	10^−3
Range of the entropy Srange	±3 nodes
Spin change angle θ	π/100
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.

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 cluster-building successes per 100 trials.

Figure 4.

Comparison of the adding probabilities for the one-bridge system.

Figure 5.

Comparison of the adding probabilities for the seven-bridge system.

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 cluster-building procedure can distinguish the hard and soft phase with relatively high 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.

Figure 6.

Reverse magnetization curves simulated using the classical P_add and modified $P_{\mathit{add}}^{\prime }$ adding probability.

The difference between P_add and $P_{\mathit{add}}^{\prime }$ 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_BT 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_BT = 1e–4 eV (see Figure 7).

Figure 7.

Spin configurations depicting the reverse magnetization process for the seven-bridge system at k_B T = 1e‒4 eV. Red and blue colors indicate the soft and hard phase, respectively.

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_{\mathit{add}}^{\prime }$ 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_{\mathit{add}}^{\prime }$, 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_{\mathit{add}}^{\prime }$ results in “faster” (in lower fields) minimization of the system energy than the procedure based on the classical adding probability.

Figure 8.

Reverse magnetization curve and energy of the system in the fields related to the spin flip of the soft phase.

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.

5. 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.

Acknowledgments

This work was supported by the National Science Centre in Poland by grant no. 2015/19/B/ST8/02636 (A. Chrobak, G. Ziółkowski). D. Chrobak acknowledges support from the National Science Centre of Poland, grant no. 2016/21/B/ST8/02737.