=Paper=
{{Paper
|id=Vol-2426/paper25
|storemode=property
|title=Core Method of Numerical Calculation of Vector Models Density of States
|pdfUrl=https://ceur-ws.org/Vol-2426/paper25.pdf
|volume=Vol-2426
|authors=Anton A. Kuzin,Vladislav S. Strongin,Sergey A. Anisimov,Maxim V. Lebedev,Konstantin V. Nefedev,Valery I. Belokon
}}
==Core Method of Numerical Calculation of Vector Models Density of States==
https://ceur-ws.org/Vol-2426/paper25.pdf
Core Method of Numerical Calculation of Vector Models Density of
States
Anton A. Kuzin¹, Vladislav S. Strongin¹, Sergey A. Anisimov¹, Maxim V. Lebedev¹,
Konstantin V. Nefedev², Valery I. Belokon¹
¹Far Eastern Federal University, Vladivostok, Russia, kuzin_aa@dvfu.ru
²Institute of Applied Mathematics, Far Eastern Branch, Russian Academy of Science,
Vladivostok, Russia
Abstract
This article presents a new method for the numerical calculation of the
density of states in vector models. The method is based on quasi-Markov
random walks in the state space of vector models and is multi-spin. An
exhaustive enumeration of all cluster states for each of the 2L states of the
boundary of length L makes it possible to obtain a local distribution of
configurations, and such a piecewise method allows one to approximately
calculate the entire state space. The method is designed to calculate the
degeneracy rate G, the interaction energy E, and the magnetization M for
any given lattice geometry, knowing which Gibbs approach can be used
in order to calculate the partition function and thermodynamic averages.
1 Introduction
The models under consideration should be studied from different points of view [3].
To study the behavior of spin systems, you need to know magnetic susceptibility, heat capacity and entropy.
According to the definition, magnetic susceptibility characterizes the connection between the magnetization of the
magnetic system and the external field:
𝜕〈𝑀〉(𝑇)
𝜒(𝑇)|ℎ→0 = (1)
𝜕ℎ
In this case, the average magnetization at a given temperature
1 𝐸 +ℎ𝑀𝑖
〈𝑀〉(𝑇) |ℎ→0 = ∑𝑖 𝑀𝑖 𝑒𝑥𝑝 [− 𝑖 ] (2)
𝑍 𝑘𝐵 𝑇
a statistical amount
𝐸 +ℎ𝑀𝑖
𝑍 = ∑𝑖 𝑒𝑥𝑝 [− 𝑖 ] (3)
𝑘𝐵 𝑇
Combining equations (2) and (3) we get
𝑖 𝐸 +ℎ𝑀
𝑖
𝜕〈𝑀〉 𝜕 ∑ 𝑀𝑖 𝑒𝑥𝑝[− 𝑘𝐵 𝑇 ]
𝜒(𝑇)|ℎ→0 = | = | (4)
𝜕ℎ ℎ→0 𝜕ℎ ∑ 𝑒𝑥𝑝[−𝐸𝑖 +ℎ𝑀𝑖 ]
𝑘𝐵 𝑇
ℎ→0
For convenience, you can enter an abbreviation:
1
𝛽= .
𝑘𝐵 𝑇
Because
𝓊 / 𝓊 / 𝜐−𝜐/ 𝓊
( ) = (5)
𝜐 𝜐2
Imagine the numerator and denominator as:
𝓊 = ∑ 𝑀𝑖 𝑒𝑥𝑝 [−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )] (6)
𝜐 = ∑ 𝑒𝑥𝑝 [−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )]
and their derivatives:
Copyright © 2019 for the individual papers by the papers' authors. Copyright © 2019 for the volume as a collection by its
editors. This volume and its papers are published under the Creative Commons License Attribution 4.0 International (CC BY 4.0).
In: Sergey I. Smagin, Alexander A. Zatsarinnyy (eds.): V International Conference Information Technologies and High-
Performance Computing (ITHPC-2019), Khabarovsk, Russia, 16-19 Sep, 2019, published at http://ceur-ws.org
167
Mathematical Modeling in Physics and Technology
______________________________________________________________________________________
𝜕𝓊
𝓊/ = = 𝛽 ∑ 𝑀𝑖 2 𝑒𝑥𝑝 [−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )] (7)
𝜕ℎ
𝜕𝜐
𝜐/ = = 𝛽 ∑ 𝑀𝑖 𝑒𝑥𝑝 [−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )],
𝜕ℎ
this implies
𝛽(𝑀𝑖2 exp ∑ [−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )])exp ∑ [−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )] 𝛽(𝑀𝑖 exp[−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )])2
𝜒(𝑇)|ℎ→0 = − =
(∑ exp[−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )])2 (∑ exp[−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )])2
2
(∑𝑀𝑖2 exp[−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )]) (𝑀𝑖 exp[−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )])
𝛽( −( ) )
(∑ exp[−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )]) (∑ exp[−𝛽(𝐸𝑖 + ℎ𝑀𝑖 )])
⟨𝑀2 ⟩−⟨𝑀⟩2
𝜒(𝑇)|ℎ→0 = (8)
𝑘𝐵 𝑇
In turn, the heat capacity
𝜕⟨𝐸⟩(𝑇) 𝜕𝛽 𝜕⟨𝐸⟩(𝑇)
𝐶(𝑇) = = , (9)
𝜕𝑇 𝜕𝑇 𝜕𝛽
where the average energy at a given temperature is
∑𝐸 exp[−𝛽𝐸𝑖 ]
⟨𝐸⟩(𝑇) = 𝑖 (10)
∑ exp[−𝛽𝐸𝑖 ]
Similar to magnetic susceptibility, the abbreviation is introduced:
1
𝛽=
𝑘𝐵 𝑇
𝐸
𝑢 = ∑𝐸𝑖 exp [− 𝑖 ] = ∑𝐸𝑖 exp[−𝛽𝐸𝑖 ] (11)
𝑘𝐵𝑇
𝐸𝑖
𝑣 = ∑𝐸𝑖 exp [− ] = ∑ exp[−𝛽𝐸𝑖 ]
𝑘𝐵𝑇
and their derivatives:
−𝐸𝑖2 exp[−𝛽𝐸𝑖 ]∑ exp[−𝛽𝐸𝑖 ] − 𝐸𝑖 exp[−𝛽𝐸𝑖 ]∑ 𝐸𝑖 exp[−𝛽𝐸𝑖 ]
𝐶(𝑇) = (−𝛽 2 − − )
(∑ exp[−𝛽𝐸𝑖 ])2 (∑ exp[−𝛽𝐸𝑖 ])2
2
2
𝛴𝐸𝑖2 exp[−𝛽𝐸𝑖 ] ∑𝐸𝑖 exp[−𝛽𝐸𝑖 ]
=𝛽 ( −( ))
∑ exp[−𝛽𝐸𝑖 ] ∑ exp[−𝛽𝐸𝑖 ]
⟨𝐸 2 ⟩−⟨𝐸⟩2
𝐶(𝑇) = . (12)
𝑘𝐵 𝑇 2
If we represent the multiplicity of degeneracy of states as a function of the energy 𝐺(𝐸), then the entropy for a
given energy level, according to the Boltzmann definition, is given by 𝑆(𝐸) = 𝑘𝐵 ln 𝐺(𝐸). Entropy can be considered
as a thermodynamic quantity that determines the measure of energy dissipation in the system and is given by the
relation
𝑈−𝐹 〈𝐸〉(𝑇)
𝑆(𝑇) = = 𝑙𝑛𝑍 + (13)
𝑘𝐵 𝑇 𝑘𝐵 𝑇
where U is the internal energy of the system, 𝐹 = 𝑈 − 𝑘𝐵 𝑙𝑛𝑍 is the Helmholtz free energy. In equation (10), 〈𝐸〉(𝑇)
is the Gibbs thermodynamic averaging for the energy level. Similarly, the distribution of the energy level of entropy
can be obtained:
−𝐸𝑖
∑𝑖 𝐺(𝐸𝑖 )𝑙𝑛𝐺(𝐸𝑖 ) exp[ ]
𝑘𝐵 𝑇
〈𝑆〉 = −𝐸𝑖 (14)
∑ 𝐺(𝐸𝑖 )exp[ ]
𝑘𝐵 𝑇
2 A Brief Presentation of the Method and the Results of the Software Product.
Visualization of Results
2.1 Theoretical Component of the Method
For vector models of Ising spins, the interaction energy is usually given by the Hamiltonian of the form
Е = − ∑<𝑖,𝑗> 𝐽𝑖𝑗 𝑆𝑖 𝑆𝑗 , (15)
where 𝑆𝑖 is «spins», or one-component vectors whose components take the values either «+1» or «-1», and 𝐽𝑖𝑗 is a
constant of exchange interaction, usually takes the value «+1» for ferromagnets, t. e. if the vectors 𝑆𝑖 and 𝑆𝑗 are co-
directed, then this corresponds to a minimum of energy, «-1» is an antiferromagnet, where the antiparallel orientation
of vectors leads to minimization of energy. The magnetization of the N spin Ising system is calculated by the formula:
М = ∑𝑁 𝑖=1 𝑆𝑖 (16)
Ising spins are placed in space, for example, in a lattice - a metal model, or randomly for modeling amorphous
structures.
To solve the problem of thermodynamics of the system of interacting Ising spins (calculate the heat capacity,
magnetization, magnetic susceptibility, and other statistical characteristics as a function of temperature), you need to
168
Mathematical Modeling in Physics and Technology
______________________________________________________________________________________
know how many states G (configurations) will have given values of E, M. Thus, we are interested in Any given
system geometry that has only three global numbers G, E, M. In total there will be 2 N states.
Even in the interaction model of the nearest environment, there are serious computational problems associated
with an exponentially rapid increase in the number of possible configurations. For example, for a system of Ising
spins on a three-dimensional pyrochlore lattice with 16 nodes there will be 2 16 configurations.
Figure 1: Structure of the pyrochlore lattice (left) [1]. The state space of the Ising spin system with the interaction of
the nearest environment on the pyrochlore lattice obtained by the method of exhaustive enumeration (right) [2]
histogram G, E, M.
2.2 Brief Introduction of the Method
The essence of the proposed method is to get one random configuration from the 2 N space; this is quite easy to do;
we simply set random values of the vector component Si. One configuration corresponds to one point in the histogram
of fig. 1 (right). Select the nucleus (region) of closely spaced spins. We carry out an exhaustive enumeration for each
of the 2L states of the boundary, where L is the number of spins on the boundary and each configuration of the
nucleus from 2N and we obtain local numbers g, e, m taking into account the interaction with the boundaries.
Figure 2: Schematic representation of the Core (cluster) method.
If the core contains 4x4 spins, then there will be a total of 65,536 kernel configurations for one configuration.
However, not all of them will have the same values of E i, Mk, since in systems with the interaction of the nearest
environment, degeneracy in E and M is essential. It is possible to use any geometry of the kernel.
We find on the histogram of fig. 1 (on the right), the point with these coordinates corresponding to each value of Е i
and Mk is added to the value of gi that exists at this point g. Thus, we must go through each point in this projection of
the histogram, along the horizontal and vertical axes. You may need to make several passes.
The absence of new points Ei Mk in the process of full coverage of the projection of the histogram can be selected
as a criterion for stopping the numerical calculation. After stopping, we will need to normalize G. There are several
ways of normalization. For example:
The sum of all histogram columns is 2N
(1−𝑀𝑘 )
The sum of the histogram columns with the given value Ei, Mk is the binomial coefficient 𝐵 = (𝑁, ).
2
If we know the degeneracy of one or several points, for example, when all the spins are «up» or all «down» in an
antiferromagnet, we can control the normalization of these points.
169
Mathematical Modeling in Physics and Technology
______________________________________________________________________________________
Knowing the numbers G, E, M, we can use the Gibbs approach to calculate the partition function and
thermodynamic averages.
To test our method, a software product was implemented, the algorithm of which is presented below.
2.3 Algorithm for the Implementation of a Software Product
We fill the double array with all kernel variants (each spin value takes the value -1 or +1), where L is the linear
size of the core. Then we display an array of a kernel of size LхL, write it into a file. Next, we find and return the
numbers of configurations with minimal energy and magnetization. We determine the nearest neighbors of all spins
and return the energy of the system. We define the point for building the kernel and derive it. We define the
periodicity condition and output the array of the kernel.
Set the boundaries of the nucleus and determine the boundary conditions. Calculate and return the magnetization
of the entire system. Find the interaction energy of spins inside the nucleus without taking into account the
boundaries, where n is the linear size of the system, N is the number of spins in the system, L is the linear size of the
nucleus, Nbound is the number of spins on the boundary, Ncore is the number of spins in the core, allconfig is the
number of all configurations - maximum number of neighbors of each particle, Emax - maximum number of system
energy, Mmax - maximum number of system magnetization, dE and dM - system energy without core energy /
without core magnetization, Ecore and Mcore - core energy / magnetization, Etot and Mtot - total energy /
magnetization of the whole system, spins - array in which the values of spins are stored (±1), bound is the array in
which the boundaries are stored, core is the array in which the core is stored, e is the array in which the core energy is
stored, m-array in which the core magnetization is stored, b is the double array (allconfig size on Ncore) with all
kernel variants (± 1), gem is the array in which EM and g are stored (E and M are array indices, and g is the value of
these indices).
In the next stage, creating an array with all the configurations and display the indexes of the system spins. We
prepare the starting configuration, that is, we scatter +1 or -1 values throughout the system. Calculate and derive the
energy and magnetization of the system. We determine the interaction of the spins with the boundary, then sum up
with what was calculated without considering the boundary and obtaining the energy.
The next step is to calculate the energy and magnetization for all configurations by iterating through all
configurations. Next, fill the core with spins, which correspond to the configuration. We consider the magnetization
and energy of the system without taking into account the interaction with the boundaries. We consider the energy of
interaction of spins with the boundary and record the energy of the nucleus. We write in the gem energy and
magnetization of the nucleus and then we find the configuration with the minimum energy and magnetization and
write the core with the minimum energy and magnetization into a common system. Output to file: energy and
magnetization of the nucleus; energy and magnetization of the system with a new core; E M and g.
2.4. Calculated Density of States
According to the results of calculations, the following histograms were constructed with a grid size of 4x4 (Fig. 3)
and 5x5 (Fig. 4)
Figure 3: Two-dimensional representation of DOS for systems of 4x4. a) one of two possible configurations with
minimum energy (E=-32, M=0), b) configuration with maximum energy (E=32, M=16)
170
Mathematical Modeling in Physics and Technology
______________________________________________________________________________________
Figure 4: DOS 3D view for 5x5 systems
3 Conclusion
The obtained approximate solutions by the cluster method of numerical calculation of the density of states of
vector models exhaustively repeat the exact solution in relation to E and M. The accuracy of the calculations is
controlled by comparing the values of the multiplicity of degenerations G obtained approximately and the exact
values.
References
1. Andriushchenko, P.: Large peaks in the entropy of the diluted nearest-neighbor spin-ice model on the pyrochlore
lattice in a [111] magnetic field / P. Andriushchenko, K. Soldatov, A. Peretyatko, Y. Shevchenko, K.
Nefedev [and etc.] // PHYSICAL REVIEW. - № 022138. - P. 1-12. (2019)
2. Kuzin, A. A.: Diluted spin ice in an external magnetic field / A. A. Kuzin, A. A. Peretyatko, K. S. Soldatov, K. V.
Nefediev // Far Eastern Mathematical Journal. - № 1. - P. 59–81. (2017)
3. Kim, S.: Exact Solution of the Ising Model on the Square Lattice with Free Boundary Conditions/ Kim Seung-
Yeon // World Academy of Science, Engineering and Technology International Journal of Mathematical,
Computational, Physical, Electrical and Computer Engineering Vol:5, No:11 (2011)
171