A guide of the parameters used in the algorithm.
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.
- Monte Carlo simulations
- cluster Monte Carlo methods
- magnetic simulations
- entropy methods
- ultra-high coercive alloys
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 Tb2Fe14B grains [1, 2]. In this system, Tb and Fe magnetic moments are coupled antiferromagnetically, which is responsible for relatively low magnetic remanence (μ0
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)  as well as by Wolff  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 . 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 . Furthermore, let the system be in equilibrium with a thermostat having a temperature . According to basic principles of statistical mechanics, the probability that the system occupies the state
where 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 2100 states—the number that makes the summation occurring in Eq. (1) impossible.
In order to estimate the average , a nonuniform sampling of the system states can be applied. If denotes a set of indices of system states selected with the probability , the equilibrium value of the observable can be modeled  by the estimator:
where probabilities were selected to be equal . All we need is a method that generates a set of system states with the Boltzmann probabilities . Because the exact value of the partition function 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 (from
The transition probability can be considered as a product of the selection probability and the acceptance ratio (probability). In general, the selection probabilities can be chosen to a large extent freely, e.g., they can be symmetrical . In that case the acceptance probabilities satisfy the equation:
As an example, let us consider Ising model of interacting spins placed at the nodes of a two-dimensional regular. The energy of the system is, then, given by the formula
where describes the spin state at the th lattice node and 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 (
The solution of the critical slowing down problem was proposed by Swendsen-Wang  and later by Wolff . The approach developed by Wolf (
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 that can refer, for example, magnetic anisotropy or some other properties potentially affecting the clusterization of the system. Let us define a sphere around a node of the spin lattice—the sphere containing nodes in total. Furthermore, let the system property be characterized by a discrete and finite set of values: . The local distribution of the system property is defined by the numbers , where stands for the number of nodes having the value . The sum overall is equals . Thus, the particular -state of the sphere is defined by the set of numbers . The number of possible realizations of the -state denoted by is given by the expression
Using the Stirling formula, the previous equation takes the form
So the number of “microstates” corresponding to the -state is
The probability of particular local distribution of the system property is proportional to . If the only one value of dominates its distribution inside the sphere , e.g., , the local information entropy achieves the lowest value . Other local distribution that involves more than one value of results in .
Now we can introduce the perception of
which gives if the case of fully ordered system property (). One can see that decreases when local disorder of system property increases ().
We assume that local disorder of system property decreases the probability of cluster growth and constitutes barriers for cluster expansion. Thus, one could expect nonuniform distribution of adding probability among the system nodes. In order to derive a reasonable formula for , 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 by . Then, assuming that this product satisfies the Wolff’s formula , 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
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
First of all, a random node
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
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
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 |
|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|
|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|
|The number of iterations contained in one MC step||Typically, we suggest |
|The number of MC steps used to calculate the average spin of the system <SZ>||We propose |
|Modifier of the impact of local information entropy||For α < 1 and for α > 1, the |
|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 |
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 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
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.
|Parameter||Soft phase||Hard phase|
|Exchange coupling ||1.5e–2 eV||1.5e–2 eV|
|Anisotropy constant ||0||1e–3 eV|
|Dipolar constant ||1.8e–7 eV|
|Thermal energy kB||1e–4, 1 e–5, 1e–6 eV|
|α parameter||Random value from 0 to 25|
|Cluster analysis probability ||10−3|
|Range of the entropy ||±3 nodes|
|Spin change angle ||π/100|
|System size (in one direction) ||50|
|Number of iteration in one MC step||3n3 = 3,75,000|
|Number of MC steps for magnetization averaging||400|
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.
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.
The difference between
Regarding the fact that in some cases also for
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.
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
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.
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.
Chrobak A, Ziółkowski G, Randrianantoandro N. Magnetic hardening of Fe–Nb–B–Tb type of bulk nanocrystalline alloys. Journal of Alloys and Compounds. 2014; 583:48-54
Chrobak A, Zółkowski G, Randrianantoandro N, Klimontko J, Chrobak D, Prusik K, et al. Ultra-high coercivity of (Fe86xNbxB14)0.88Tb0.12 bulk nanocrystalline magnets. Acta Materialia. 2015; 98:318-326
Binder K, editor. The Monte Carlo Method in Condensed Matter Physics. 2nd ed. Berlin: Springer; 1995
Metropolis N, Rosenbluth AW, Teller AH, Teller E. Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics. 1953; 21:1087
Fortuin M, Kasteleyn PW. On the random-cluster model: I. Introduction and relation to other models. Physica. 1972; 57:536
Coniglio A, Klein W. Clusters and Ising critical droplets: A renormalisation group approach. Journal of Physics A. 1980; 13:2775
Swendsen RH, Wang JS. Nonuniversal critical dynamics in Monte Carlo simulations. Physical Review Letters. 1987; 58:8688
Wolff U. Collective Monte Carlo Updating for Spin Systems. Physical Review Letters. 1989; 62:361
Plascak JA, Ferrenberg AM, Landau DP. Cluster hybrid Monte Carlo simulation algorithms. Physical Review E. 2002; 65:066702-1-066702-8
Ferrenberg AM, Xu J, Landau DP. Pushing the limits of Monte Carlo simulations for the three-dimensional Ising model. Physical Review E. 2018; 97:043301-1-043301-12
Chrobak A, Ziółkowski G, Granek K, Chrobak D. Computer Physics Communucation. 2019; 238:157-164
Binder K. Applications of Monte Carlo methods to statistical physics. Reports on Progress in Physics. 1997; 60:487
Newman MEJ, Barkema GT. Monte Carlo Methods in Statistical Physics. New York: Oxford University Press Inc; 2001
Kalos MH, Whitlock P. Monte Carlo Methods I: Basics. New York: Wiley; 1986
Wood WW. Physics of Simple Liquids. Amsterdam: North-Holland; 1968
Coddington PD, Baillie CF. Empirical relations between static and dynamic exponents for Ising model cluster algorithms. Physical Review Letters. 1992; 68:962
Matz R, Hunter DL, Jan N. The dynamic critical exponent of the three-dimensional Ising model. Journal of Statistical Physics. 1994; 74:903
Nightingale MP, Blote HWJ. Dynamic Exponent of the Two-Dimensional Ising Model and Monte Carlo Computation of the Subdominant Eigenvalue of the Stochastic Matrix. Physical Review Letters. 1996; 76:4548
Jaynes ET. Information Theory and Statistical Mechanics. Physics Review. 1957; 106:620
Rosenthal JS. Parallel computing and Monte Carlo algorithms. Far East Journal of Theoretical Statistics. 2000; 4:207-236
Hallekalek P. Proceedings of the 1998 Workshop on Parallel and Distributed Simulation; Vol. 1: pp. 82-89
Khan MO, Kennedy G, Derek YC, Chan J. Computers & Chemistry. 2005; 26:72-77
Campillo F, Rakotozafy R, Rossi V. Mathematics and Computers in Simulation. 2009; 79(12):3424-3433