3.1. Classical spin Hamiltonian

The magnetic systems simulated using Vegas are described by the classical Heisenberg model. Interactions at the atomic scale between spin moments in a magnetic material are modeled by the Heisenberg spin Hamiltonian, which includes the exchange, anisotropy and applied field interactions. The Heisenberg spin Hamiltonian is given by

where J _ij is the exchange interaction constant between sites i and j, B and n _i are the magnetic field intensity and direction of site i, respectively, μs is the atomic magnetic moment, S _i and S _j are the spin moment directions of sites i and j, respectively, and Hani is the anisotropy term. The anisotropy term can take the form of uniaxial anisotropy

where kiuni and e _i are the uniaxial anisotropy constant and direction of site i, respectively, or cubic anisotropy ¹⁵

where kicub is the cubic anisotropy constant, and S_x, S_y and S_z are the x, y and z components of the spin moment S, respectively.

3.2. Monte Carlo Metropolis algorithm

The time evolution of a classical spin magnetic system can be given by the Monte Carlo Metropolis algorithm, where new configurations of a system are generated from a previous state using a transition probability. The algorithm is developed as follows: an initial state is chosen for the magnetic system. Then, one single spin moment is randomly chosen and its direction is changed to a new trial direction according to an update policy (trial move). The energy of the system before and after the trial move is then evaluated according to Eq. (1) and the energy change ΔE is calculated. Finally, the trial move is accepted with a probability given by

where kB is the Boltzmann constant and T is the absolute temperature of the system. Repeating this process once per spin moment is called a MCS. The update policy has influence on the efficiency of the algorithm and its physical interpretation ³. Common update policies, including the spin-flip, random, small step and Gaussian moves ¹⁶, are implemented in Vegas. Furthermore, an adaptive move which keeps the acceptance rate of the new states near 50%, enhancing the efficiency of the Monte Carlo Metropolis algorithm is implemented.

For a comprehensive review of the atomistic spin model and the integration methods employed in the atomistic simulation of magnetic materials see ¹³.

3.3. File Formats

In order to facilitate the exchange and archiving of data, Vegas uses the standardized file formats, Json and Hdf5, for the input and output files of the simulations, respectively.

JSON uses universal data structures through which the simulation parameters can be defined using a name/value pair structure. An example of a JSON input file is shown in Fig. 1. In this file, the simulation parameters, including the temperature, magnetic field range and step, number of MCS and random number seed, for a hysteresis loop simulation are defined. Also, the material, anisotropy and initial state (if required) files, and the results output file are given. In order to take into account the atomic magnetic moment, the field input parameters corresponds to quantities of μsB (see Eq. (1)). In this case, the simulation parameters are normalized such that kB=1. Then, the exchange interaction, anisotropy and field parameters must be given in units of temperature. However, if the user intends to use other kind of normalization, it is necessary to use coherent units for all the parameters. For instance, if the user intends to use real units in the MKS system, the Boltzmann constant must be given in J/K and the exchange interaction, anisotropy and field parameters in J.

FIGURE 1 Example of a JSON input file defining the simulation parameters. 

On the other hand, Hdf5 is a data model and file format for storing and managing extremely large data, which is useful for storing and analyzing the simulation results in Vegas, specially the history of the simulation. Other input data, including the anisotropy and the material files, are given in plain text format. For these files, we designed simple data schemes to specify the anisotropy type and the magnetic and structural parameters of the material.

3.4. Vegas tools

Based on different python and C++ libraries, Vegas provides tools to build, simulate and analyze magnetic systems.

4.1.Materials with different shapes and structures

Materials with simple and complex shapes, such as core/shell nanoparticles, nanotubes, nanowires, nanorings, bit patterned media, thin films, multilayers and bulk systems, and crystal structures, such as simple, body-centered and face-centered cubic, can be simulated in Vegas. Samples of some materials are shown in Fig. 2.

FIGURE 2 Samples of a a) spherical core/shell nanoparticle, b) core/shell nanowire and c) multilayer. 

4.2. Thermal averages

Thermal averages can be easily extracted and visualized from the simulation results output file. Figure 3 shows the thermal dependence of the magnetization (M) and the susceptibility (X), and a hysteresis loop of a bulk material simulated using the structure and magnetic parameters of iron, namely, J=44.01 meV/link and k=3.53×10-3 meV/atom ¹³. From these results, important magnetic properties of the material, such as the critical temperature (T _c) and the coercivity (H_c), can be determined.

FIGURE 3 a) Thermal dependence of the magnetization and susceptibility, and b) hysteresis loop at 10 K of an iron bulk material. Inset in Fig. 3a) shows the BCC crystal structure of iron. 

The hard disk space required for a simulation depends on parameters as the size of the system and the number of Monte Carlo, temperature, and field steps of the simulation. For instance, the simulation of the thermal dependence of M and X (see Fig. 3a) was made considering a system with N=2000 sites with 8 nearest neighbors each one, BCC structure and uniaxial anisotropy. The temperature was varied from a high to a low value with a total of 200 temperature steps, while the magnetic field was set to a single value of 0 T. 20000 MCS were performed at each temperature step. For this simulation, the size of the Hdf5 file, which stores the history of the simulation, is approximately 209 MB. The running time, in a Intel(R) Xeon(R) CPU E5-2680 v2 @ 2.80 GHz with an x86_64 architecture, was 1.31 hours ± 24.48 seconds. All the information, including the thermal averages, magnetization states and dynamic data of the system, is extracted from this single file.

4.3. Dynamic behavior

Besides the static magnetic properties, the dynamic properties of a magnetic system can be extracted from a single simulation ¹⁹. Figure 4 shows the equilibrium correlation times (τ) for a ferromagnetic bulk material with generic parameters. The correlation times were computed for some of the different spin update policies implemented in Vegas. Overall, the adaptive policy produces the lowest correlation times at all temperatures, indicating that it is a very efficient update policy for this kind of system.

FIGURE 4 Equilibrium correlation times for a bulk material using adaptive, Hinzke-Nowak combinational ¹⁶ (for a small-step opening angle of 15 and 30±) and random update policies. 

4.4. Critical exponents

Phase transitions can be characterized by simulating the critical behavior of the magnetic systems. Using the finite size scaling method, it is possible to extract values for critical exponents by observing the variation of the thermal averages with the system size ^19,20. Figure 5 shows log-log plots of the maxima of the susceptibility (χmax), and the first ( V1max) and second order (V2max) cumulants of the magnetization as a function of the system size (L) for a ferromagnetic 2D Ising system with generic parameters. The critical exponents γ and v are obtained from the slope of the plots in Figs. 5a and 5b, respectively.

FIGURE 5 Log-log plots of (a) the maxima of the susceptibility, and (b) the first and second order cumulants of the magnetization as a function of system size for a 2D Ising system. 

4.5. Magnetization switching processes

Simulation of magnetization switching processes are of great importance for the development of recording devices based on magnetic nanoparticles. Figures 6a and 6b show cross-sectional views of magnetization states in a ferromagnetic torus nanoring. Two magnetization states, vortex (see Fig. 6a) and reverse vortex (see Fig. 6) states, can be formed in the torus nanoring and convenientbly employed to represent the bits 0 and 1. Figure 6c shows the temperature dependence of the magnetization in the θ^ direction (Mθ) of two independent simulations where a torus nanoring with generic parameters was cooled down from a high temperature. During the cooling down process and in a short temperature range (shaded region), an external circular magnetic field was applied. In one of the simulations, the external magnetic field switches the magnetization of the torus nanoring producing a reverse vortex state (solid line), while in the other simulation, it does not switch the magnetization producing a vortex state (dashed line). This kind of behavior indicates that there is a switching probability which depends on the temperature range at which the external magnetic field is applied ⁸.

FIGURE 6 a) Vortex and b) reverse vortex magnetization states in a torus nanoring, and c) temperature dependence of the magnetization in the θ^ direction of a torus nanoring cooled down from a high temperature. 

4.6. Histogram and multiple histogram methods

From results obtained using Vegas, it is also possible to apply data analysis techniques such as the single and multiple histogram method. These methods allow to take a single simulation performed at some temperature and extrapolate or interpolate results to give predictions of observable quantities at other temperatures ²⁰. Figure 7 show histograms of the energy per ion site (E/N) obtained at different temperatures in a ferromagnetic 2D Ising system with generic parameters. Using the multiple histogram method, it is possible to interpolate observable quantities, such as the energy and magnetization, at other temperatures with high accuracy.

FIGURE 7 Histograms of the energy per ion site obtained at different temperatures in a ferromagnetic 2D Ising system. 

4.7. Identification of magnetization states

The feature of Vegas that allows to store the direction of every spin moment at the last MCS is very useful to identify different magnetization states, which is difficult from the thermal average quantities alone. Figure 8 shows a cross-sectional view of different magnetization states produced in a spherical ferromagnetic nanoparticle with core kc and Néel surface ks anisotropy. The magnetization state depends on the ratio ks/kc. When ks/kc→0, 100, ∞ and -100, “collinear”, “throttled”, “hedgehog” and “artichoke” states are produced, respectively ⁶.

FIGURE 8 a) “Collinear”, b) “throttled”, c) “hedgehog” and d) “artichoke” magnetization states produced in a ferromagnetic spherical nanoparticle with core and Néel surface anisotropy. 

Figures 2, 6a, 6b, 8 and the inset in Fig. 3a were generated using POV-Ray ²¹.

Worked examples, tutorials and a detailed description of the use of Vegas can be found at its web page https://pcm-ca.github.io/vegas/.